C $Header: /u/gcmpack/MITgcm/pkg/fizhi/getpwhere.F,v 1.11 2010/03/16 00:19:33 jmc Exp $
C $Name: $
#include "FIZHI_OPTIONS.h"
#undef CHECK_GETPWHERE
subroutine GETPWHERE(myThid,numpress,pressures,levpressures)
C***********************************************************************
C subroutine getpwhere
C
C Purpose: Approximate (!) the level at which the mid-level pressure
C is less than (ie, above in the atmosphere) a given value.
C
C Algorithm: Assume surface pressure of 1000 mb, and the pressure
c thicknesses set in make_phys_grid, with drF thicknesses above
C
C Need: Information about the dynamics grid vertical spacing
C
C Input: myThid - process(or) number
C numpres - Number of pressures to process
C pressures - Pressure values to find levels for
C
C Output: levpressures - Array of levels at which pressures are found
C These pressure levels correspond to the fizhi
C levels, and assume that the levels are counted
C from top to bottom (CRITICAL!)
C
C NOTE: The new physics levels specified here MUST correspond to the
C physics levels specified in make_phys_grid from gridalt package.
C***********************************************************************
implicit none
c
#include "SIZE.h"
#include "fizhi_SIZE.h"
#include "GRID.h"
#ifdef CHECK_GETPWHERE
#include "EEPARAMS.h"
#include "PARAMS.h"
#endif /* CHECK_GETPWHERE */
integer myThid,numpress
_RL pressures(numpress)
integer levpressures(numpress)
c
#ifdef CHECK_GETPWHERE
INTEGER ioUnit, upperBnd
c CHARACTER*(MAX_LEN_MBUF) msgBuf
#endif /* CHECK_GETPWHERE */
integer n,L,dynlev
C Code that MUST correspond to make_phys_grid in the gridalt package!
C Require 12 bottom levels (300 mb worth) for the physics,
C Counting from bottom to top, the dp(s) are:
integer ntry,ntry10,ntry40
parameter (ntry40=15)
parameter (ntry10=12)
_RL dptry(ntry40), dptry10(ntry10), dptry40(ntry40)
_RL dptry_pedge(ntry40+1)
#ifdef TRY_NEW_GETPWHERE
_RL rC_dyn(Nr), dpTop
#else
_RL rF_pmid(Nr),rF_edge(Nr+1)
#endif
_RL pref(Nrphys)
integer tmplev
data dptry10 /3.00, 6.00,10.00,14.00,17.00,25.00,
. 25.00,25.00,25.00,50.00,50.00,50.00/
data dptry40 /3.00, 6.00,10.00,14.00,17.00,25.00,
. 25.00,25.00,25.00,25.00,25.00,25.00,
. 25.00,25.00,25.00/
if( (Nr.eq.10).or.(Nr.eq.20) ) then
ntry = ntry10
do L = 1,ntry
dptry(L) = dptry10(L)
enddo
elseif ((Nr.eq.40).or.(Nr.eq.70)) then
ntry = ntry40
do L = 1,ntry
dptry(L) = dptry40(L)
enddo
else
print *,' Dont know how to set levels for given pressures '
stop
endif
C define the mid pressure for the levels that are specified - bottom 300 mb.
#ifdef TRY_NEW_GETPWHERE
dptry_pedge(1) = rF(1)*1. _d -2
#else /* TRY_NEW_GETPWHERE */
dptry_pedge(1) = 1000.
#endif /* TRY_NEW_GETPWHERE */
do L = 1,ntry
dptry_pedge(L+1) = dptry_pedge(L) - dptry(L)
enddo
do L = 1,ntry
pref(L) = (dptry_pedge(L) + dptry_pedge(L+1))/2.
enddo
#ifdef CHECK_GETPWHERE
ioUnit = errorMessageUnit
WRITE(ioUnit,'(A)') '===== GETPWHERE: CHECK start ====='
WRITE(ioUnit,'(4(A,I6),A)') ' Nr =', Nr,
& ' , Nrphys=', Nrphys,
& ' , ntry =', ntry,
& ' , pref(1:ntry):'
WRITE(ioUnit,'(10F10.4)') (pref(L),L=1,ntry)
#endif /* CHECK_GETPWHERE */
#ifdef TRY_NEW_GETPWHERE
C top levels DP of 1 mb (or 0.01 mb for strat version)
dpTop = 1. _d 0
IF (Nr.EQ.70) dpTop = 1. _d -2
DO L = ntry+1,Nrphys
pref(L) = (Nrphys-L+0.5)*dpTop
ENDDO
C define the rest of the mid pressures from the dynamics levels
DO L = 1,Nr
rC_dyn(L) = rC(L)*1. _d -2
ENDDO
dynlev = 0
DO L = 1,Nr
IF ( rC_dyn(L).GE.dptry_pedge(ntry+1) ) dynlev = L
ENDDO
DO L = ntry+1,MIN(Nrphys,ntry+Nr-dynlev)
pref(L) = rC_dyn(dynlev+L-ntry)
ENDDO
#else /* TRY_NEW_GETPWHERE */
C define the rest of the mid pressures from the dynamics levels
rF_edge(1) = 1000.
do L = 2,Nr+1
rF_edge(L) = rF_edge(L-1) - (drF(L-1)/100.)
enddo
do L = 1,Nr
rF_pmid(L) = (rF_edge(L) + rF_edge(L+1))/2.
enddo
dynlev = 0
do L = 1,Nr
if(rF_pmid(L).ge.pref(ntry)) dynlev = L
enddo
#ifdef CHECK_GETPWHERE
IF ( rF_pmid(dynlev).ge.pref(ntry)-25. ) THEN
upperBnd = ntry+Nr-dynlev
ELSE
upperBnd = ntry+Nr-dynlev-1
ENDIF
WRITE(ioUnit,'(1(A,I5),A)') ' Up-Bnd=', upperBnd,
& ' , rF_pmid:'
WRITE(ioUnit,'(10F10.4)') rF_pmid
IF ( upperBnd.LT.Nrphys-1 ) THEN
WRITE(ioUnit,'(A)')
& 'ERROR: exeeding "rF_pmid" array bounds => pref ill defined'
c STOP 'ABNORMAL END: S/R GETPWHERE'
ENDIF
#endif /* CHECK_GETPWHERE */
if(rF_pmid(dynlev).ge.pref(ntry)-25.) then
do L = ntry+1,Nrphys-1
pref(L) = rF_pmid(dynlev+L-ntry)
enddo
else
pref(ntry) = rF_pmid(dynlev)
do L = ntry+1,Nrphys-1
pref(L) = rF_pmid(dynlev+L-ntry+1)
enddo
endif
C Add top level DP of 1 mb - p mid is at 0.5 mb (or 0.05 mb for strat version)
pref(Nrphys) = 0.5
if(Nr.eq.70)pref(Nrphys) = 0.005
#endif /* TRY_NEW_GETPWHERE */
DO n = 1,numpress
DO L = 1,Nrphys
IF ( pref(L).GE.pressures(n) ) tmplev = L
ENDDO
C and now flip the level numbers for the top down counting in fizhi
levpressures(n) = Nrphys+1-tmplev
ENDDO
#ifdef CHECK_GETPWHERE
WRITE(ioUnit,'(3(A,I5),A)') ' dynlev=', dynlev,
& ' , numpress=', numpress,
& ' , pressures:'
WRITE(ioUnit,'(10F10.4)') pressures
WRITE(ioUnit,'(A)') 'pref(ntry:Nrphys):'
WRITE(ioUnit,'(10F10.4)') (pref(L),L=ntry,Nrphys)
WRITE(ioUnit,'(A)') 'levpressures:'
WRITE(ioUnit,'(20I5)') levpressures
WRITE(ioUnit,'(A)') '===== GETPWHERE: CHECK end ====='
#endif /* CHECK_GETPWHERE */
RETURN
END