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