pkg/gridalt/make_phys_grid.F

32f2418b0c Jean*0001 #include "GRIDALT_OPTIONS.h"
                0002
                0003        SUBROUTINE MAKE_PHYS_GRID(drF,hfacC,im1,im2,jm1,jm2,Nr,
                0004      & nSx,nSy,i1,i2,j1,j2,bi,bj,Nrphys,Lbot,dpphys,numlevphys,nlperdyn)
                0005 C***********************************************************************
                0006 C SUBROUTINE MAKE_PHYS_GRID
                0007 C
                0008 C Purpose: Define the grid that the will be used to run the high-end
                0009 C          atmospheric physics.
                0010 C
                0011 C Algorithm: Fit additional levels of some (~) known thickness in
                0012 C          between existing levels of the grid used for the dynamics
                0013 C
                0014 C Need:    Information about the dynamics grid vertical spacing
                0015 C
                0016 C Input:   drF         - delta r (p*) edge-to-edge
                0017 C          hfacC       - fraction of grid box above topography
                0018 C          im1, im2    - beginning and ending i - dimensions
                0019 C          jm1, jm2    - beginning and ending j - dimensions
                0020 C          Nr          - number of levels in dynamics grid
                0021 C          nSx,nSy     - number of processes in x and y direction
                0022 C          i1, i2      - beginning and ending i - index to fill
                0023 C          j1, j2      - beginning and ending j - index to fill
                0024 C          bi, bj      - x-dir and y-dir index of process
                0025 C          Nrphys      - number of levels in physics grid
                0026 C
                0027 C Output:  dpphys      - delta r (p*) edge-to-edge of physics grid
                0028 C          numlevphys  - number of levels used in the physics
                0029 C          nlperdyn    - physics level number atop each dynamics layer
                0030 C
                0031 C NOTES: 1) Pressure levs are built up from bottom, using p0, ps and dp:
                0032 C              p(i,j,k)=p(i,j,k-1) + dp(k)*ps(i,j)/p0(i,j)
                0033 C        2) Output dp(s) are aligned to fit EXACTLY between existing
                0034 C           levels of the dynamics vertical grid
                0035 C        3) IMPORTANT! This routine assumes the levels are numbered
                0036 C           from the bottom up, ie, level 1 is the surface.
                0037 C           IT WILL NOT WORK OTHERWISE!!!
                0038 C        4) This routine does NOT work for surface pressures less
                0039 C           (ie, above in the atmosphere) than about 350 mb
                0040 C***********************************************************************
                0041
                0042        IMPLICIT NONE
                0043
                0044        integer im1,im2,jm1,jm2,Nr,nSx,nSy,Nrphys
8aa94f3b31 Andr*0045        integer i1,i2,j1,j2,bi,bj
                0046        integer numlevphys
32f2418b0c Jean*0047        _RS hfacC(im1:im2,jm1:jm2,Nr,nSx,nSy)
                0048        _RL dpphys(im1:im2,jm1:jm2,Nrphys,nSx,nSy)
8b4e6631e3 Jean*0049        _RS drF(Nr)
32f2418b0c Jean*0050        integer Lbot(im1:im2,jm1:jm2,nSx,nSy)
                0051        integer nlperdyn(im1:im2,jm1:jm2,Nr,nSx,nSy)
                0052
8aa94f3b31 Andr*0053        integer i,j,L,Lbotij,Lnew
32f2418b0c Jean*0054 C Require 12 (or 15) levels near the surface (300 mb worth) for fizhi.
                0055 C   the dp(s) are in the dptry arrays:
13837d01f6 Andr*0056        integer ntry,ntry10,ntry40
3768cb558d Andr*0057        parameter (ntry10 = 12)
                0058        parameter (ntry40 = 12)
                0059        _RL dptry(15),dptry10(ntry10),dptry40(ntry40)
                0060        _RL bot_thick,bot_thick40
13837d01f6 Andr*0061        _RL dptry_accum(15)
                0062        data dptry10/300.000,600.000,1000.000,1400.000,1700.000,2500.000,
32f2418b0c Jean*0063      &            2500.000,2500.000,2500.000,5000.000,5000.000,5000.000/
3768cb558d Andr*0064        data dptry40/300.000,600.000, 800.000, 800.000,1250.000,
32f2418b0c Jean*0065      &            1250.000,2500.000,2500.000,2500.000,2500.000,2500.000,
                0066      &            2500.000/
3768cb558d Andr*0067        data bot_thick40/20000.000/
8aa94f3b31 Andr*0068        _RL deltap, dpstar_accum
                0069        integer nlbotmax, nstart, nlevs, nlphys, ndone
b8ce59ac21 Andr*0070        _RL thindp
32f2418b0c Jean*0071
13837d01f6 Andr*0072        if( (Nr.eq.10) .or. (Nr.eq.20) ) then
                0073         ntry = ntry10
3768cb558d Andr*0074         bot_thick = bot_thick40
13837d01f6 Andr*0075         do L = 1,ntry
                0076          dptry(L) = dptry10(L)
                0077         enddo
b8ce59ac21 Andr*0078        elseif((Nr.eq.40).or.(Nr.eq.46).or.(Nr.eq.70)) then
13837d01f6 Andr*0079         ntry = ntry40
3768cb558d Andr*0080         bot_thick = bot_thick40
13837d01f6 Andr*0081         do L = 1,ntry
                0082          dptry(L) = dptry40(L)
                0083         enddo
                0084        else
                0085         print *,' Dont know how to make fizhi grid '
                0086         stop
                0087        endif
b8ce59ac21 Andr*0088
                0089        thindp=100.
                0090        if(Nr.eq.70)thindp=0.02
32f2418b0c Jean*0091
8aa94f3b31 Andr*0092        do L = 1,Nr
                0093         do j = j1,j2
                0094         do i = i1,i2+1
                0095          nlperdyn(i,j,L,bi,bj) = 0
                0096         enddo
                0097         enddo
                0098        enddo
32f2418b0c Jean*0099
                0100 C Figure out how many physics levels there will be
                0101 C   (need 12 between sfc and 300 mb above it - see how many
                0102 C    there are in the dynamics if the surface pressure is at
                0103 C    the sum of drF, ie, the maximum dynamics grid layers possible)
8aa94f3b31 Andr*0104        nlevs = 0
                0105        dpstar_accum = 0.
                0106        do L = 1,Nr
                0107         dpstar_accum = dpstar_accum + drF(L)
3768cb558d Andr*0108         if(dpstar_accum.le.bot_thick) nlevs = nlevs+1
8aa94f3b31 Andr*0109        enddo
                0110        numlevphys = Nr - nlevs + ntry + 1
32f2418b0c Jean*0111
8aa94f3b31 Andr*0112        dptry_accum(1) = dptry(1)
                0113        do Lnew = 2,ntry
                0114         dptry_accum(Lnew) = dptry_accum(Lnew-1) + dptry(Lnew)
                0115        enddo
32f2418b0c Jean*0116
                0117 C      do for each grid point:
8aa94f3b31 Andr*0118        do j = j1,j2
be5ec7d59d Andr*0119        do i = i1,i2
8aa94f3b31 Andr*0120         Lbotij = Lbot(i,j,bi,bj)
32f2418b0c Jean*0121
                0122 C Find the maximum number of physics levels to fit in the bottom level
                0123
8aa94f3b31 Andr*0124         nlbotmax = 0
                0125         do Lnew = 1,ntry
32f2418b0c Jean*0126          if ( (nlbotmax.eq.0) .and.
                0127      & (dptry_accum(Lnew).gt.(hfacC(i,j,Lbotij,bi,bj)*drF(Lbotij))))then
8aa94f3b31 Andr*0128           nlbotmax = Lnew
                0129          endif
                0130         enddo
                0131         if(nlbotmax.eq.0)then
                0132          nlbotmax = ntry
                0133         endif
32f2418b0c Jean*0134
                0135 C See how many of the physics levs can fit in the bottom level
                0136
8aa94f3b31 Andr*0137         nlphys = 0
                0138         deltap = 0.
                0139         do Lnew = 1,nlbotmax
32f2418b0c Jean*0140 C Test to see if the next physics level fits, if yes, add it
8aa94f3b31 Andr*0141          if((hfacC(i,j,Lbotij,bi,bj)*drF(Lbotij)).ge.
32f2418b0c Jean*0142      &                                           deltap+dptry(Lnew))then
8aa94f3b31 Andr*0143           nlphys = nlphys + 1
                0144           dpphys(i,j,nlphys,bi,bj) = dptry(Lnew)
                0145           deltap = deltap + dptry(Lnew)
                0146          else
32f2418b0c Jean*0147 C If the level does not fit, decide whether to make a new thinner
                0148 C  one or make the one below a bit thicker
8aa94f3b31 Andr*0149           if((dptry(Lnew-1)+(hfacC(i,j,Lbotij,bi,bj)*
32f2418b0c Jean*0150      &             drF(Lbotij)-deltap)) .gt. (dptry(Lnew-1)*1.5) ) then
                0151 C Add a new thin layer
8aa94f3b31 Andr*0152            nlphys = nlphys + 1
32f2418b0c Jean*0153            dpphys(i,j,nlphys,bi,bj) =
                0154      &                      (hfacC(i,j,Lbotij,bi,bj)*drF(Lbotij))-deltap
8aa94f3b31 Andr*0155           else
32f2418b0c Jean*0156 C Make the one below thicker
                0157            dpphys(i,j,nlphys,bi,bj) = dpphys(i,j,nlphys,bi,bj) +
                0158      &                      (hfacC(i,j,Lbotij,bi,bj)*drF(Lbotij)-deltap)
8aa94f3b31 Andr*0159           endif
                0160           deltap = deltap+(hfacC(i,j,Lbotij,bi,bj)*drF(Lbotij)-deltap)
                0161          endif
                0162         enddo
32f2418b0c Jean*0163
8aa94f3b31 Andr*0164         nlperdyn(i,j,Lbotij,bi,bj) = nlphys
32f2418b0c Jean*0165
                0166 C Now proceed upwards - see how many physics levels fit in each
                0167 C   subsequent dynamics level - go through all 12 required levels
                0168
8aa94f3b31 Andr*0169         do L = Lbotij+1,Nr
                0170          ndone = 0
                0171          if(nlphys.lt.ntry)then
                0172           deltap = 0.
                0173           nstart = nlphys + 1
                0174           do Lnew = nstart,ntry
                0175            if((hfacC(i,j,L,bi,bj)*drF(L)).ge.deltap+dptry(Lnew))then
                0176             nlphys = nlphys + 1
                0177             dpphys(i,j,nlphys,bi,bj) = dptry(Lnew)
                0178             deltap = deltap + dptry(Lnew)
                0179             ndone = 0
                0180            elseif (ndone.eq.0) then
32f2418b0c Jean*0181 C If the level does not fit, decide whether to make a new thinner
                0182 C  one or make the one below a bit thicker
8aa94f3b31 Andr*0183             ndone = 1
                0184             if( (dptry(Lnew-1)+(hfacC(i,j,L,bi,bj)*drF(L)-deltap))
32f2418b0c Jean*0185      &                    .gt. (dptry(Lnew-1)*1.5) ) then
                0186 C Add a new thin layer
8aa94f3b31 Andr*0187              nlphys = nlphys + 1
32f2418b0c Jean*0188              dpphys(i,j,nlphys,bi,bj) =
                0189      &                             (hfacC(i,j,L,bi,bj)*drF(L))-deltap
8aa94f3b31 Andr*0190              deltap = hfacC(i,j,L,bi,bj)*drF(L)
                0191             else
32f2418b0c Jean*0192 C Make the one below thicker
8aa94f3b31 Andr*0193              dpphys(i,j,nlphys,bi,bj) = dpphys(i,j,nlphys,bi,bj) +
32f2418b0c Jean*0194      &                      (hfacC(i,j,L,bi,bj)*drF(L)-deltap)
8aa94f3b31 Andr*0195              deltap = hfacC(i,j,L,bi,bj)*drF(L)
                0196             endif
                0197            endif
                0198           enddo
be5ec7d59d Andr*0199 C Need one more peice of logic - if we finished Lnew loop and
                0200 C  now we are done adding new physics layers, we need to be sure
                0201 C  that we are at the edge of a dynamics layer. if not, we need
                0202 C  to add one more layer.
                0203           if(nlphys.ge.ntry)then
                0204            if(abs(deltap-hfacC(i,j,L-1,bi,bj)*drF(L-1)).gt.0.001)then
                0205             nlphys = nlphys + 1
                0206             dpphys(i,j,nlphys,bi,bj) = hfacC(i,j,L-1,bi,bj)*drF(L-1)
32f2418b0c Jean*0207      &                                                          - deltap
be5ec7d59d Andr*0208            endif
                0209           endif
                0210
8aa94f3b31 Andr*0211          elseif(nlphys.eq.ntry)then
32f2418b0c Jean*0212 C Mostly done with new layers - make sure we end at dynamics edge,
                0213 C      if not, make one more thinner (thinner than dyn grid) layer
8aa94f3b31 Andr*0214           if(abs(deltap-hfacC(i,j,L-1,bi,bj)*drF(L-1)).gt.0.001)then
                0215            nlphys = nlphys + 1
32f2418b0c Jean*0216            dpphys(i,j,nlphys,bi,bj) = hfacC(i,j,L-1,bi,bj)*drF(L-1)
                0217      &                                                          - deltap
8aa94f3b31 Andr*0218            nlphys = nlphys + 1
                0219            dpphys(i,j,nlphys,bi,bj) = hfacC(i,j,L,bi,bj)*drF(L)
                0220           else
                0221            nlphys = nlphys + 1
                0222            dpphys(i,j,nlphys,bi,bj) = hfacC(i,j,L,bi,bj)*drF(L)
                0223           endif
                0224          else
32f2418b0c Jean*0225 C we are done adding new physics layers, just copy the rest
                0226 C    of the dynamics grid onto the physics grid
8aa94f3b31 Andr*0227           nlphys = nlphys + 1
                0228           dpphys(i,j,nlphys,bi,bj) = hfacC(i,j,L,bi,bj)*drF(L)
                0229          endif
                0230          nlperdyn(i,j,L,bi,bj) = nlphys
                0231         enddo
32f2418b0c Jean*0232
                0233 C  All done adding layers - if we need more to make numlevphys, put
                0234 C     them as thin (1 mb) layers near the top
8aa94f3b31 Andr*0235        if(nlphys.lt.numlevphys)then
                0236         nlevs = numlevphys-nlphys
b8ce59ac21 Andr*0237         dpphys(i,j,nlphys,bi,bj)=dpphys(i,j,nlphys,bi,bj)-thindp*nlevs
8aa94f3b31 Andr*0238         do Lnew = nlphys+1,numlevphys
b8ce59ac21 Andr*0239          dpphys(i,j,Lnew,bi,bj) = thindp
8aa94f3b31 Andr*0240         enddo
                0241         nlperdyn(i,j,Nr,bi,bj) = numlevphys
                0242        endif
32f2418b0c Jean*0243 C END OF LOOP OVER GRID POINTS
13837d01f6 Andr*0244
8aa94f3b31 Andr*0245        enddo
                0246        enddo
                0247
                0248        return
                0249        end