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