pkg/gridalt/make_phys_grid.F

C $Header: /u/gcmpack/MITgcm/pkg/gridalt/make_phys_grid.F,v 1.7 2006/07/06 13:47:28 molod Exp $
C $Name:  $

       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
c
#include "CPP_OPTIONS.h"

       integer im1,im2,jm1,jm2,Nr,Nsx,Nsy,Nrphys
       integer i1,i2,j1,j2,bi,bj
       integer numlevphys
       _RL hfacC(im1:im2,jm1:jm2,Nr,Nsx,Nsy)
       _RL dpphys(im1:im2,jm1:jm2,Nrphys,Nsx,Nsy)
       _RL drF(Nr)
       integer Lbot(im1:im2,jm1:jm2,Nsx,Nsy)
       integer nlperdyn(im1:im2,jm1:jm2,Nr,Nsx,Nsy)
c
       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
c
       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)) 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
c
       do L = 1,Nr
        do j = j1,j2
        do i = i1,i2+1
         nlperdyn(i,j,L,bi,bj) = 0
        enddo
        enddo
       enddo
c
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
c
       dptry_accum(1) = dptry(1)
       do Lnew = 2,ntry
        dptry_accum(Lnew) = dptry_accum(Lnew-1) + dptry(Lnew)
       enddo
c
c      do for each grid point:
       do j = j1,j2
       do i = i1,i2
        Lbotij = Lbot(i,j,bi,bj)
c
c Find the maximum number of physics levels to fit in the bottom level
c
        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
c See how many of the physics levs can fit in the bottom level
c
        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
c
        nlperdyn(i,j,Lbotij,bi,bj) = nlphys
c
c Now proceed upwards - see how many physics levels fit in each
c   subsequent dynamics level - go through all 12 required levels
c                
        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
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)-100. * nlevs
        do Lnew = nlphys+1,numlevphys
         dpphys(i,j,Lnew,bi,bj) = 100.
        enddo
        nlperdyn(i,j,Nr,bi,bj) = numlevphys
       endif
c END OF LOOP OVER GRID POINTS

       enddo
       enddo

       return
       end
1	molod	1.8	C $Header: /u/gcmpack/MITgcm/pkg/gridalt/make_phys_grid.F,v 1.7 2006/07/06 13:47:28 molod Exp $
2	edhill	1.3	C $Name: $
3
4	molod	1.1	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			_RL hfacC(im1:im2,jm1:jm2,Nr,Nsx,Nsy)
50			_RL dpphys(im1:im2,jm1:jm2,Nrphys,Nsx,Nsy)
51			_RL 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	molod	1.4	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	molod	1.8	parameter (ntry10 = 12)
60			parameter (ntry40 = 12)
61			_RL dptry(15),dptry10(ntry10),dptry40(ntry40)
62			_RL bot_thick,bot_thick40
63	molod	1.4	_RL dptry_accum(15)
64			data dptry10/300.000,600.000,1000.000,1400.000,1700.000,2500.000,
65	molod	1.1	. 2500.000,2500.000,2500.000,5000.000,5000.000,5000.000/
66	molod	1.8	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	molod	1.1	_RL deltap, dpstar_accum
71			integer nlbotmax, nstart, nlevs, nlphys, ndone
72	molod	1.2	c
73	molod	1.4	if( (Nr.eq.10) .or. (Nr.eq.20) ) then
74			ntry = ntry10
75	molod	1.8	bot_thick = bot_thick40
76	molod	1.4	do L = 1,ntry
77			dptry(L) = dptry10(L)
78			enddo
79	molod	1.8	elseif((Nr.eq.40).or.(Nr.eq.46)) then
80	molod	1.4	ntry = ntry40
81	molod	1.8	bot_thick = bot_thick40
82	molod	1.4	do L = 1,ntry
83			dptry(L) = dptry40(L)
84			enddo
85			else
86			print *,' Dont know how to make fizhi grid '
87			stop
88			endif
89			c
90	molod	1.1	do L = 1,Nr
91			do j = j1,j2
92			do i = i1,i2+1
93			nlperdyn(i,j,L,bi,bj) = 0
94			enddo
95			enddo
96			enddo
97			c
98			c Figure out how many physics levels there will be
99			c (need 12 between sfc and 300 mb above it - see how many
100			c there are in the dynamics if the surface pressure is at
101			c the sum of drF, ie, the maximum dynamics grid layers possible)
102			nlevs = 0
103			dpstar_accum = 0.
104			do L = 1,Nr
105			dpstar_accum = dpstar_accum + drF(L)
106	molod	1.8	if(dpstar_accum.le.bot_thick) nlevs = nlevs+1
107	molod	1.1	enddo
108			numlevphys = Nr - nlevs + ntry + 1
109			c
110			dptry_accum(1) = dptry(1)
111			do Lnew = 2,ntry
112			dptry_accum(Lnew) = dptry_accum(Lnew-1) + dptry(Lnew)
113			enddo
114			c
115			c do for each grid point:
116			do j = j1,j2
117	molod	1.2	do i = i1,i2
118	molod	1.1	Lbotij = Lbot(i,j,bi,bj)
119			c
120			c Find the maximum number of physics levels to fit in the bottom level
121			c
122			nlbotmax = 0
123			do Lnew = 1,ntry
124			if ( (nlbotmax.eq.0) .and.
125			. (dptry_accum(Lnew).gt.(hfacC(i,j,Lbotij,bi,bj)*drF(Lbotij))))then
126			nlbotmax = Lnew
127			endif
128			enddo
129			if(nlbotmax.eq.0)then
130			nlbotmax = ntry
131			endif
132			c
133			c See how many of the physics levs can fit in the bottom level
134			c
135			nlphys = 0
136			deltap = 0.
137			do Lnew = 1,nlbotmax
138			c Test to see if the next physics level fits, if yes, add it
139			if((hfacC(i,j,Lbotij,bi,bj)*drF(Lbotij)).ge.
140			. deltap+dptry(Lnew))then
141			nlphys = nlphys + 1
142			dpphys(i,j,nlphys,bi,bj) = dptry(Lnew)
143			deltap = deltap + dptry(Lnew)
144			else
145			c If the level does not fit, decide whether to make a new thinner
146			c one or make the one below a bit thicker
147			if((dptry(Lnew-1)+(hfacC(i,j,Lbotij,bi,bj)*
148			. drF(Lbotij)-deltap)) .gt. (dptry(Lnew-1)*1.5) ) then
149			c Add a new thin layer
150			nlphys = nlphys + 1
151			dpphys(i,j,nlphys,bi,bj) =
152			. (hfacC(i,j,Lbotij,bi,bj)*drF(Lbotij))-deltap
153			else
154			c Make the one below thicker
155			dpphys(i,j,nlphys,bi,bj) = dpphys(i,j,nlphys,bi,bj) +
156			. (hfacC(i,j,Lbotij,bi,bj)*drF(Lbotij)-deltap)
157			endif
158			deltap = deltap+(hfacC(i,j,Lbotij,bi,bj)*drF(Lbotij)-deltap)
159			endif
160			enddo
161			c
162			nlperdyn(i,j,Lbotij,bi,bj) = nlphys
163			c
164			c Now proceed upwards - see how many physics levels fit in each
165			c subsequent dynamics level - go through all 12 required levels
166			c
167			do L = Lbotij+1,Nr
168			ndone = 0
169			if(nlphys.lt.ntry)then
170			deltap = 0.
171			nstart = nlphys + 1
172			do Lnew = nstart,ntry
173			if((hfacC(i,j,L,bi,bj)*drF(L)).ge.deltap+dptry(Lnew))then
174			nlphys = nlphys + 1
175			dpphys(i,j,nlphys,bi,bj) = dptry(Lnew)
176			deltap = deltap + dptry(Lnew)
177			ndone = 0
178			elseif (ndone.eq.0) then
179			c If the level does not fit, decide whether to make a new thinner
180			c one or make the one below a bit thicker
181			ndone = 1
182			if( (dptry(Lnew-1)+(hfacC(i,j,L,bi,bj)*drF(L)-deltap))
183			. .gt. (dptry(Lnew-1)*1.5) ) then
184			c Add a new thin layer
185			nlphys = nlphys + 1
186			dpphys(i,j,nlphys,bi,bj) =
187			. (hfacC(i,j,L,bi,bj)*drF(L))-deltap
188			deltap = hfacC(i,j,L,bi,bj)*drF(L)
189			else
190			c Make the one below thicker
191			dpphys(i,j,nlphys,bi,bj) = 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			endif
195			endif
196			enddo
197	molod	1.2	C Need one more peice of logic - if we finished Lnew loop and
198			C now we are done adding new physics layers, we need to be sure
199			C that we are at the edge of a dynamics layer. if not, we need
200			C to add one more layer.
201			if(nlphys.ge.ntry)then
202			if(abs(deltap-hfacC(i,j,L-1,bi,bj)*drF(L-1)).gt.0.001)then
203			nlphys = nlphys + 1
204			dpphys(i,j,nlphys,bi,bj) = hfacC(i,j,L-1,bi,bj)*drF(L-1)
205			. - deltap
206			endif
207			endif
208
209	molod	1.1	elseif(nlphys.eq.ntry)then
210			c Mostly done with new layers - make sure we end at dynamics edge,
211			c if not, make one more thinner (thinner than dyn grid) layer
212			if(abs(deltap-hfacC(i,j,L-1,bi,bj)*drF(L-1)).gt.0.001)then
213			nlphys = nlphys + 1
214			dpphys(i,j,nlphys,bi,bj) = hfacC(i,j,L-1,bi,bj)*drF(L-1)
215			. - deltap
216			nlphys = nlphys + 1
217			dpphys(i,j,nlphys,bi,bj) = hfacC(i,j,L,bi,bj)*drF(L)
218			else
219			nlphys = nlphys + 1
220			dpphys(i,j,nlphys,bi,bj) = hfacC(i,j,L,bi,bj)*drF(L)
221			endif
222			else
223			c we are done adding new physics layers, just copy the rest
224			c of the dynamics grid onto the physics grid
225			nlphys = nlphys + 1
226			dpphys(i,j,nlphys,bi,bj) = hfacC(i,j,L,bi,bj)*drF(L)
227			endif
228			nlperdyn(i,j,L,bi,bj) = nlphys
229			enddo
230			c
231			c All done adding layers - if we need more to make numlevphys, put
232	molod	1.7	c them as thin (1 mb) layers near the top
233	molod	1.1	if(nlphys.lt.numlevphys)then
234			nlevs = numlevphys-nlphys
235	molod	1.7	dpphys(i,j,nlphys,bi,bj) = dpphys(i,j,nlphys,bi,bj)-100. * nlevs
236	molod	1.1	do Lnew = nlphys+1,numlevphys
237	molod	1.7	dpphys(i,j,Lnew,bi,bj) = 100.
238	molod	1.1	enddo
239			nlperdyn(i,j,Nr,bi,bj) = numlevphys
240			endif
241			c END OF LOOP OVER GRID POINTS
242	molod	1.4
243	molod	1.1	enddo
244			enddo
245
246			return
247			end