| 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 |
Parent Directory
| 
Revision Log
| 
Revision Graph