forked from darwinproject/darwin3
-
Notifications
You must be signed in to change notification settings - Fork 1
/
make_phys_grid.F
249 lines (228 loc) · 9.03 KB
/
make_phys_grid.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
#include "GRIDALT_OPTIONS.h"
SUBROUTINE MAKE_PHYS_GRID(drF,hfacC,im1,im2,jm1,jm2,Nr,
& nSx,nSy,i1,i2,j1,j2,bi,bj,Nrphys,Lbot,dpphys,numlevphys,nlperdyn)
C***********************************************************************
C SUBROUTINE MAKE_PHYS_GRID
C
C Purpose: Define the grid that the will be used to run the high-end
C atmospheric physics.
C
C Algorithm: Fit additional levels of some (~) known thickness in
C between existing levels of the grid used for the dynamics
C
C Need: Information about the dynamics grid vertical spacing
C
C Input: drF - delta r (p*) edge-to-edge
C hfacC - fraction of grid box above topography
C im1, im2 - beginning and ending i - dimensions
C jm1, jm2 - beginning and ending j - dimensions
C Nr - number of levels in dynamics grid
C nSx,nSy - number of processes in x and y direction
C i1, i2 - beginning and ending i - index to fill
C j1, j2 - beginning and ending j - index to fill
C bi, bj - x-dir and y-dir index of process
C Nrphys - number of levels in physics grid
C
C Output: dpphys - delta r (p*) edge-to-edge of physics grid
C numlevphys - number of levels used in the physics
C nlperdyn - physics level number atop each dynamics layer
C
C NOTES: 1) Pressure levs are built up from bottom, using p0, ps and dp:
C p(i,j,k)=p(i,j,k-1) + dp(k)*ps(i,j)/p0(i,j)
C 2) Output dp(s) are aligned to fit EXACTLY between existing
C levels of the dynamics vertical grid
C 3) IMPORTANT! This routine assumes the levels are numbered
C from the bottom up, ie, level 1 is the surface.
C IT WILL NOT WORK OTHERWISE!!!
C 4) This routine does NOT work for surface pressures less
C (ie, above in the atmosphere) than about 350 mb
C***********************************************************************
IMPLICIT NONE
integer im1,im2,jm1,jm2,Nr,nSx,nSy,Nrphys
integer i1,i2,j1,j2,bi,bj
integer numlevphys
_RS hfacC(im1:im2,jm1:jm2,Nr,nSx,nSy)
_RL dpphys(im1:im2,jm1:jm2,Nrphys,nSx,nSy)
_RS drF(Nr)
integer Lbot(im1:im2,jm1:jm2,nSx,nSy)
integer nlperdyn(im1:im2,jm1:jm2,Nr,nSx,nSy)
integer i,j,L,Lbotij,Lnew
C Require 12 (or 15) levels near the surface (300 mb worth) for fizhi.
C the dp(s) are in the dptry arrays:
integer ntry,ntry10,ntry40
parameter (ntry10 = 12)
parameter (ntry40 = 12)
_RL dptry(15),dptry10(ntry10),dptry40(ntry40)
_RL bot_thick,bot_thick40
_RL dptry_accum(15)
data dptry10/300.000,600.000,1000.000,1400.000,1700.000,2500.000,
& 2500.000,2500.000,2500.000,5000.000,5000.000,5000.000/
data dptry40/300.000,600.000, 800.000, 800.000,1250.000,
& 1250.000,2500.000,2500.000,2500.000,2500.000,2500.000,
& 2500.000/
data bot_thick40/20000.000/
_RL deltap, dpstar_accum
integer nlbotmax, nstart, nlevs, nlphys, ndone
_RL thindp
if( (Nr.eq.10) .or. (Nr.eq.20) ) then
ntry = ntry10
bot_thick = bot_thick40
do L = 1,ntry
dptry(L) = dptry10(L)
enddo
elseif((Nr.eq.40).or.(Nr.eq.46).or.(Nr.eq.70)) then
ntry = ntry40
bot_thick = bot_thick40
do L = 1,ntry
dptry(L) = dptry40(L)
enddo
else
print *,' Dont know how to make fizhi grid '
stop
endif
thindp=100.
if(Nr.eq.70)thindp=0.02
do L = 1,Nr
do j = j1,j2
do i = i1,i2+1
nlperdyn(i,j,L,bi,bj) = 0
enddo
enddo
enddo
C Figure out how many physics levels there will be
C (need 12 between sfc and 300 mb above it - see how many
C there are in the dynamics if the surface pressure is at
C the sum of drF, ie, the maximum dynamics grid layers possible)
nlevs = 0
dpstar_accum = 0.
do L = 1,Nr
dpstar_accum = dpstar_accum + drF(L)
if(dpstar_accum.le.bot_thick) nlevs = nlevs+1
enddo
numlevphys = Nr - nlevs + ntry + 1
dptry_accum(1) = dptry(1)
do Lnew = 2,ntry
dptry_accum(Lnew) = dptry_accum(Lnew-1) + dptry(Lnew)
enddo
C do for each grid point:
do j = j1,j2
do i = i1,i2
Lbotij = Lbot(i,j,bi,bj)
C Find the maximum number of physics levels to fit in the bottom level
nlbotmax = 0
do Lnew = 1,ntry
if ( (nlbotmax.eq.0) .and.
& (dptry_accum(Lnew).gt.(hfacC(i,j,Lbotij,bi,bj)*drF(Lbotij))))then
nlbotmax = Lnew
endif
enddo
if(nlbotmax.eq.0)then
nlbotmax = ntry
endif
C See how many of the physics levs can fit in the bottom level
nlphys = 0
deltap = 0.
do Lnew = 1,nlbotmax
C Test to see if the next physics level fits, if yes, add it
if((hfacC(i,j,Lbotij,bi,bj)*drF(Lbotij)).ge.
& deltap+dptry(Lnew))then
nlphys = nlphys + 1
dpphys(i,j,nlphys,bi,bj) = dptry(Lnew)
deltap = deltap + dptry(Lnew)
else
C If the level does not fit, decide whether to make a new thinner
C one or make the one below a bit thicker
if((dptry(Lnew-1)+(hfacC(i,j,Lbotij,bi,bj)*
& drF(Lbotij)-deltap)) .gt. (dptry(Lnew-1)*1.5) ) then
C Add a new thin layer
nlphys = nlphys + 1
dpphys(i,j,nlphys,bi,bj) =
& (hfacC(i,j,Lbotij,bi,bj)*drF(Lbotij))-deltap
else
C Make the one below thicker
dpphys(i,j,nlphys,bi,bj) = dpphys(i,j,nlphys,bi,bj) +
& (hfacC(i,j,Lbotij,bi,bj)*drF(Lbotij)-deltap)
endif
deltap = deltap+(hfacC(i,j,Lbotij,bi,bj)*drF(Lbotij)-deltap)
endif
enddo
nlperdyn(i,j,Lbotij,bi,bj) = nlphys
C Now proceed upwards - see how many physics levels fit in each
C subsequent dynamics level - go through all 12 required levels
do L = Lbotij+1,Nr
ndone = 0
if(nlphys.lt.ntry)then
deltap = 0.
nstart = nlphys + 1
do Lnew = nstart,ntry
if((hfacC(i,j,L,bi,bj)*drF(L)).ge.deltap+dptry(Lnew))then
nlphys = nlphys + 1
dpphys(i,j,nlphys,bi,bj) = dptry(Lnew)
deltap = deltap + dptry(Lnew)
ndone = 0
elseif (ndone.eq.0) then
C If the level does not fit, decide whether to make a new thinner
C one or make the one below a bit thicker
ndone = 1
if( (dptry(Lnew-1)+(hfacC(i,j,L,bi,bj)*drF(L)-deltap))
& .gt. (dptry(Lnew-1)*1.5) ) then
C Add a new thin layer
nlphys = nlphys + 1
dpphys(i,j,nlphys,bi,bj) =
& (hfacC(i,j,L,bi,bj)*drF(L))-deltap
deltap = hfacC(i,j,L,bi,bj)*drF(L)
else
C Make the one below thicker
dpphys(i,j,nlphys,bi,bj) = dpphys(i,j,nlphys,bi,bj) +
& (hfacC(i,j,L,bi,bj)*drF(L)-deltap)
deltap = hfacC(i,j,L,bi,bj)*drF(L)
endif
endif
enddo
C Need one more peice of logic - if we finished Lnew loop and
C now we are done adding new physics layers, we need to be sure
C that we are at the edge of a dynamics layer. if not, we need
C to add one more layer.
if(nlphys.ge.ntry)then
if(abs(deltap-hfacC(i,j,L-1,bi,bj)*drF(L-1)).gt.0.001)then
nlphys = nlphys + 1
dpphys(i,j,nlphys,bi,bj) = hfacC(i,j,L-1,bi,bj)*drF(L-1)
& - deltap
endif
endif
elseif(nlphys.eq.ntry)then
C Mostly done with new layers - make sure we end at dynamics edge,
C if not, make one more thinner (thinner than dyn grid) layer
if(abs(deltap-hfacC(i,j,L-1,bi,bj)*drF(L-1)).gt.0.001)then
nlphys = nlphys + 1
dpphys(i,j,nlphys,bi,bj) = hfacC(i,j,L-1,bi,bj)*drF(L-1)
& - deltap
nlphys = nlphys + 1
dpphys(i,j,nlphys,bi,bj) = hfacC(i,j,L,bi,bj)*drF(L)
else
nlphys = nlphys + 1
dpphys(i,j,nlphys,bi,bj) = hfacC(i,j,L,bi,bj)*drF(L)
endif
else
C we are done adding new physics layers, just copy the rest
C of the dynamics grid onto the physics grid
nlphys = nlphys + 1
dpphys(i,j,nlphys,bi,bj) = hfacC(i,j,L,bi,bj)*drF(L)
endif
nlperdyn(i,j,L,bi,bj) = nlphys
enddo
C All done adding layers - if we need more to make numlevphys, put
C them as thin (1 mb) layers near the top
if(nlphys.lt.numlevphys)then
nlevs = numlevphys-nlphys
dpphys(i,j,nlphys,bi,bj)=dpphys(i,j,nlphys,bi,bj)-thindp*nlevs
do Lnew = nlphys+1,numlevphys
dpphys(i,j,Lnew,bi,bj) = thindp
enddo
nlperdyn(i,j,Nr,bi,bj) = numlevphys
endif
C END OF LOOP OVER GRID POINTS
enddo
enddo
return
end