C $Header: /u/gcmpack/MITgcm/model/src/ini_p_ground.F,v 1.11 2016/04/04 21:29:00 jmc Exp $
C $Name: $
#include "CPP_OPTIONS.h"
#undef CHECK_ANALYLIC_THETA
CBOP
C !ROUTINE: INI_P_GROUND
C !INTERFACE:
SUBROUTINE INI_P_GROUND(selectMode,
& Hfld, Pfld,
I myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE INI_P_GROUND
C | o Convert Topography [m] to (reference) Surface Pressure
C | according to tRef profile,
C | using same discretisation as in calc_phi_hyd
C |
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ==
#include "SIZE.h"
#include "GRID.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "SURFACE.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C selectMode :: > 0 = find Pfld from Hfld ; < 0 = compute Hfld from Pfld
C selectFindRoSurf = 0 : use Theta_Ref profile
C selectFindRoSurf = 1 : use analytic fct Theta(Lat,P)
C Hfld (input/outp) :: Ground elevation [m]
C Pfld (outp/input) :: reference Pressure at the ground [Pa]
INTEGER selectMode
_RS Hfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RS Pfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
INTEGER myThid
C !LOCAL VARIABLES:
C == Local variables ==
C msgBuf :: Informational/error message buffer
C-
C For an accurate definition of the reference surface pressure,
C define a High vertical resolution (in P):
C nLevHvR :: Number of P-level used for High vertical Resolution (HvR)
C plowHvR :: Lower bound of the High vertical Resolution
C dpHvR :: normalized pressure increment (HvR)
C pLevHvR :: normalized P-level of the High vertical Resolution
C pMidHvR :: normalized mid point level (HvR)
C thetaHvR :: potential temperature at mid point level (HvR)
C PiHvR :: Exner function at P-level
C dPiHvR :: Exner function difference between 2 P-levels
C recip_kappa :: 1/kappa = Cp/R
C PiLoc, zLoc, dzLoc, yLatLoc, phiLoc :: hold on temporary values
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CHARACTER*(MAX_LEN_MBUF) msgBuf
INTEGER bi,bj,i,j,k, ks
_RL Po_surf
_RL hRef(2*Nr+1), rHalf(2*Nr+1)
LOGICAL findPoSurf
INTEGER nLevHvR
PARAMETER ( nLevHvR = 60 )
_RL plowHvR, dpHvR, pLevHvR(nLevHvR+1), pMidHvR(nLevHvR)
_RL thetaHvR(nLevHvR), PiHvR(nLevHvR+1), dPiHvR(nLevHvR)
_RL recip_kappa, PiLoc, zLoc, dzLoc, yLatLoc, phiLoc
_RL psNorm, rMidKp1
_RL ratioRm(Nr), ratioRp(Nr)
INTEGER kLev
#ifdef CHECK_ANALYLIC_THETA
_RL tmpVar(nLevHvR,61)
#endif
CEOP
IF ( selectFindRoSurf.LT.0 .OR. selectFindRoSurf.GT.1 ) THEN
WRITE(msgBuf,'(A,I2,A)')
& 'INI_P_GROUND: selectFindRoSurf =', selectFindRoSurf,
& ' <== bad value !'
CALL PRINT_ERROR( msgBuf, myThid )
STOP 'INI_P_GROUND'
ENDIF
DO k=1,Nr
rHalf(2*k-1) = rF(k)
rHalf(2*k) = rC(k)
ENDDO
rHalf(2*Nr+1) = rF(Nr+1)
C- Reference Geopotential at Half levels :
C Tracer level: hRef(2k) ; Interface_W level: hRef(2k+1)
C- Convert phiRef to Z unit :
DO k=1,2*Nr+1
hRef(k) = phiRef(k)*recip_gravity
ENDDO
IF (selectFindRoSurf.EQ.0 .AND. selectMode .GT. 0 ) THEN
C- Find Po_surf : Linear between 2 half levels :
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
DO j=1,sNy
DO i=1,sNx
ks = 1
DO k=1,2*Nr
IF (Hfld(i,j,bi,bj).GE.hRef(k)) ks = k
ENDDO
Po_surf = rHalf(ks) + (rHalf(ks+1)-rHalf(ks))*
& (Hfld(i,j,bi,bj)-hRef(ks))/(hRef(ks+1)-hRef(ks))
c IF (Hfld(i,j,bi,bj).LT.hRef(1)) Po_surf= rHalf(1)
c IF (Hfld(i,j,bi,bj).GT.hRef(2*Nr+1)) Po_surf=rHalf(2*Nr+1)
Pfld(i,j,bi,bj) = Po_surf
ENDDO
ENDDO
ENDDO
ENDDO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- endif selectFindRoSurf=0 & selectMode > 0
ENDIF
IF ( selectFindRoSurf.EQ.1 ) THEN
C-- define high resolution Pressure discretization:
recip_kappa = 1. _d 0 / atm_kappa
plowHvR = 0.4 _d 0
dpHvR = nLevHvR
dpHvR = (1. - plowHvR) / dpHvR
pLevHvR(1)= rF(1)/atm_Po
PiHvR(1) = atm_Cp*(pLevHvR(1)**atm_kappa)
DO k=1,nLevHvR
pLevHvR(k+1)= pLevHvR(1) - k*dpHvR
PiHvR(k+1) = atm_Cp*(pLevHvR(k+1)**atm_kappa)
pMidHvR(k)= (pLevHvR(k)+pLevHvR(k+1))*0.5 _d 0
dPiHvR(k) = PiHvR(k) - PiHvR(k+1)
ENDDO
C-- to modify pressure when using non fully linear relation between Phi & p
C (Integr_GeoPot=2 & Tracer Point at middle between 2 interfaces)
DO k=1,Nr
ratioRm(k) = 1. _d 0
ratioRp(k) = 1. _d 0
IF (k.GT.1 ) ratioRm(k) = 0.5 _d 0*drC(k)/(rF(k)-rC(k))
IF (k.LT.Nr) ratioRp(k) = 0.5 _d 0*drC(k+1)/(rC(k)-rF(k+1))
ENDDO
#ifdef CHECK_ANALYLIC_THETA
_BEGIN_MASTER( myThid )
DO j=1,61
yLatLoc =-90.+(j-1)*3.
CALL ANALYLIC_THETA( yLatLoc, pMidHvR,
& tmpVar(1,j), nLevHvR, myThid )
ENDDO
OPEN(88,FILE='analytic_theta',
& STATUS='unknown',FORM='unformatted')
WRITE(88) tmpVar
CLOSE(88)
_END_MASTER( myThid )
STOP 'CHECK_ANALYLIC_THETA'
#endif /* CHECK_ANALYLIC_THETA */
C-- endif selectFindRoSurf=1
ENDIF
IF ( selectFindRoSurf*selectMode .GT. 0) THEN
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C- Find Po_surf such as g*Hfld = Phi[Po_surf,theta(yLat,p)]:
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
C- start bi,bj loop:
DO j=1,sNy
DO i=1,sNx
phiLoc = Hfld(i,j,bi,bj) - hRef(1)
IF ( phiLoc .LE. 0. _d 0 ) THEN
Pfld(i,j,bi,bj) = rF(1)
ELSE
yLatLoc = yC(i,j,bi,bj)
CALL ANALYLIC_THETA( yLatLoc, pMidHvR,
& thetaHvR, nLevHvR, myThid )
zLoc = 0.
DO k=1,nLevHvR
IF (zLoc.GE.0.) THEN
C- compute phi/g corresponding to next p_level:
dzLoc = dPiHvR(k)*thetaHvR(k)*recip_gravity
IF ( phiLoc .LE. zLoc+dzLoc ) THEN
C- compute the normalized surf. Pressure psNorm
PiLoc = PiHvR(k)
& - gravity*( phiLoc - zLoc )/thetaHvR(k)
psNorm = (PiLoc/atm_Cp)**recip_kappa
C- use linear interpolation:
c psNorm = pLevHvR(k)
c & - dpHvR*( phiLoc - zLoc )/dzLoc
zLoc = -1.
ELSE
zLoc = zLoc + dzLoc
ENDIF
ENDIF
ENDDO
IF (zLoc.GE.0.) THEN
WRITE(msgBuf,'(2A)')
& 'INI_P_GROUND: FAIL in trying to find Pfld:',
& ' selectMode,i,j,bi,bj=',selectMode,i,j,bi,bj
CALL PRINT_ERROR( msgBuf, myThid )
WRITE(msgBuf,'(A,F7.1,2A,F6.4,A,F8.0)')
& 'INI_P_GROUND: Hfld=', Hfld(i,j,bi,bj), ' exceeds',
& ' Zloc(lowest P=', pLevHvR(1+nLevHvR),' )=',zLoc
CALL PRINT_ERROR( msgBuf, myThid )
STOP 'ABNORMAL END: S/R INI_P_GROUND'
ELSE
Pfld(i,j,bi,bj) = psNorm*atm_Po
ENDIF
ENDIF
ENDDO
ENDDO
IF (selectMode.EQ.2 .AND. integr_GeoPot.NE.1) THEN
C---------
C Modify pressure to account for non fully linear relation between
C Phi & p (due to numerical truncation of the Finite Difference scheme)
C---------
DO j=1,sNy
DO i=1,sNx
Po_surf = Pfld(i,j,bi,bj)
IF ( Po_surf.LT.rC(1) .AND. Po_surf.GT.rC(Nr) ) THEN
findPoSurf = .TRUE.
DO k=1,Nr
IF ( findPoSurf .AND. Po_surf.GE.rC(k) ) THEN
Po_surf = rC(k) + (Po_surf-rC(k))/ratioRm(k)
findPoSurf = .FALSE.
ENDIF
rMidKp1 = rF(k+1)
IF (k.LT.Nr) rMidKp1 = (rC(k)+rC(k+1))*0.5 _d 0
IF ( findPoSurf .AND. Po_surf.GE.rMidKp1 ) THEN
Po_surf = rC(k) + (Po_surf-rC(k))/ratioRp(k)
findPoSurf = .FALSE.
ENDIF
ENDDO
IF ( findPoSurf )
& STOP 'S/R INI_P_GROUND: Pb with selectMode=2'
ENDIF
Pfld(i,j,bi,bj) = Po_surf
ENDDO
ENDDO
C---------
ENDIF
C- end bi,bj loop.
ENDDO
ENDDO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- endif selectFindRoSurf*selectMode > 0
ENDIF
IF (selectMode .LT. 0) THEN
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C--- Compute Hfld = Phi(Pfld)/g.
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
C- start bi,bj loop:
C-- Compute Hfld from g*Hfld = hRef(Po_surf)
DO j=1,sNy
DO i=1,sNx
C- compute phiLoc = hRef(Ro_surf):
ks = kSurfC(i,j,bi,bj)
IF (ks.LE.Nr) THEN
C- sigma coord. case (=> uniform kSurfC), find true ks:
IF ( selectSigmaCoord.NE.0 ) THEN
DO k=2,Nr
IF ( Pfld(i,j,bi,bj).LT.rF(k) ) ks = k
ENDDO
ENDIF
IF ( Pfld(i,j,bi,bj).GE.rC(ks) ) THEN
phiLoc = hRef(2*ks)
& + (hRef(2*ks-1)-hRef(2*ks))
& *(Pfld(i,j,bi,bj)-rC(ks))/(rHalf(2*ks-1)-rHalf(2*ks))
ELSE
phiLoc = hRef(2*ks)
& + (hRef(2*ks+1)-hRef(2*ks))
& *(Pfld(i,j,bi,bj)-rC(ks))/(rHalf(2*ks+1)-rHalf(2*ks))
ENDIF
Hfld(i,j,bi,bj) = phiLoc
ELSE
Hfld(i,j,bi,bj) = 0.
ENDIF
ENDDO
ENDDO
IF (selectFindRoSurf.EQ.1) THEN
C-----
C goal: evaluate phi0surf (used for computing geopotential_prime = Phi - PhiRef):
C phi0surf = g*Ztopo-PhiRef(Ro_surf) if no truncation was applied to Ro_surf;
C but because of hFacMin, surf.ref.pressure (=Ro_surf) is truncated and
C phi0surf = Phi(Theta-Analytic,P=Ro_surf) - PhiRef(P=Ro_surf)
C-----
C-- Compute Hfld from g*Hfld = Phi[Po_surf,theta(yLat,p)]:
DO j=1,sNy
DO i=1,sNx
zLoc = hRef(1)
IF ( Pfld(i,j,bi,bj) .LT. rF(1) ) THEN
Po_surf = Pfld(i,j,bi,bj)
C---------
C Modify pressure to account for non fully linear relation between
C Phi & p (due to numerical truncation of the Finite Difference scheme)
IF (selectMode.EQ.-2 .AND. integr_GeoPot.NE.1) THEN
IF ( Po_surf.LT.rC(1) .AND. Po_surf.GT.rC(Nr) ) THEN
findPoSurf = .TRUE.
DO k=1,Nr
IF ( findPoSurf .AND. Po_surf.GE.rC(k) ) THEN
Po_surf = rC(k) + (Po_surf-rC(k))*ratioRm(k)
findPoSurf = .FALSE.
ENDIF
IF ( findPoSurf .AND. Po_surf.GE.rF(k+1) ) THEN
Po_surf = rC(k) + (Po_surf-rC(k))*ratioRp(k)
findPoSurf = .FALSE.
ENDIF
ENDDO
ENDIF
ENDIF
C---------
psNorm = Po_surf/atm_Po
kLev = 1 + INT( (pLevHvR(1)-psNorm)/dpHvR )
yLatLoc = yC(i,j,bi,bj)
CALL ANALYLIC_THETA( yLatLoc, pMidHvR,
& thetaHvR, kLev, myThid )
C- compute height at level pLev(kLev) just below Pfld:
DO k=1,kLev-1
dzLoc = dPiHvR(k)*thetaHvR(k)*recip_gravity
zLoc = zLoc + dzLoc
ENDDO
dzLoc = ( PiHvR(kLev)-atm_Cp*(psNorm**atm_kappa) )
& * thetaHvR(kLev)*recip_gravity
zLoc = zLoc + dzLoc
ENDIF
C- compute phi0surf = Phi[Po_surf,theta(yLat,p)] - PhiRef(Po_surf)
phi0surf(i,j,bi,bj) = gravity*(zLoc - Hfld(i,j,bi,bj))
C- save Phi[Po_surf,theta(yLat,p)] in Hfld (replacing PhiRef(Po_surf)):
Hfld(i,j,bi,bj) = zLoc
ENDDO
ENDDO
C- endif selectFindRoSurf=1
ENDIF
C- end bi,bj loop.
ENDDO
ENDDO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C-- endif selectMode < 0
ENDIF
RETURN
END
CBOP
C !SUBROUTINE: ANALYLIC_THETA
C !INTERFACE:
SUBROUTINE ANALYLIC_THETA( yLat, pNlev,
O thetaLev,
I kSize, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE ANALYLIC_THETA
C | o Conpute analyticaly the potential temperature Theta
C | as a function of Latitude (yLat) and normalized
C | pressure pNlev.
C | approximatively match the N-S symetric, zonal-mean and
C | annual-mean NCEP climatology in the troposphere.
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C yLat :: latitude (degre)
C pNlev :: normalized pressure levels
C kSize :: dimension of pNlev & ANALYLIC_THETA
C myThid :: Thread number for this instance of the routine
INTEGER kSize
_RL yLat
_RL pNlev (kSize)
_RL thetaLev(kSize)
INTEGER myThid
C !LOCAL VARIABLES:
C == Local variables ==
INTEGER k
_RL yyA, yyB, yyC, yyAd, yyBd, yyCd
_RL cAtmp, cBtmp, ttdC
_RL ppN0, ppN1, ppN2, ppN3a, ppN3b, ppN4
_RL ttp1, ttp2, ttp3, ttp4, ttp5
_RL yAtmp, yBtmp, yCtmp, yDtmp
_RL ttp2y, ttp4y, a1tmp
_RL ppl, ppm, pph, ppr
CEOP
DATA yyA , yyB , yyC , yyAd , yyBd , yyCd
& / 45. _d 0, 65. _d 0, 65. _d 0, .9 _d 0, .9 _d 0, 10. _d 0 /
DATA cAtmp , cBtmp , ttdC
& / 2.6 _d 0, 1.5 _d 0, 3.3 _d 0 /
DATA ppN0 , ppN1 , ppN2 , ppN3a , ppN3b , ppN4
& / .1 _d 0, .19 _d 0, .3 _d 0, .9 _d 0, .7 _d 0, .925 _d 0 /
DATA ttp1 , ttp2 , ttp3 , ttp4 , ttp5
& / 350. _d 0, 342. _d 0, 307. _d 0, 301. _d 0, 257. _d 0 /
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
yAtmp = ABS(yLat) - yyA
yAtmp = yyA + MIN(0. _d 0,yAtmp/yyAd) + MAX(yAtmp, 0. _d 0)
yAtmp = COS( deg2rad*MAX(yAtmp, 0. _d 0) )
yBtmp = ABS(yLat) - yyB
yBtmp = yyB + yBtmp/yyBd
yBtmp = COS( deg2rad*MAX( 0. _d 0, MIN(yBtmp,90. _d 0) ) )
yCtmp = ABS(yLat) - yyC
yCtmp = MAX( 0. _d 0, 1. _d 0 - (yCtmp/yyCd)**2 )
yDtmp = ppN3a +(ppN3b - ppN3a)*yCtmp
ttp2y = ttp3 + (ttp2-ttp3)*yAtmp**cAtmp
ttp4y = ttp5 + (ttp4-ttp5)*yBtmp**cBtmp
a1tmp = (ttp1-ttp2y)*ppN1*ppN2/(ppN2-ppN1)
DO k=1,kSize
ppl = MIN(pNlev(k),ppN1)
ppm = MIN(MAX(pNlev(k),ppN1),ppN2)
pph = MAX(pNlev(k),ppN2)
ppr =( ppN0 + ABS(ppl-ppN0) - ppN1 )/(ppN2-ppN1)
thetaLev(k) =
& ( (1. _d 0 -ppr)*ttp1*ppN1**atm_kappa
& + ppr*ttp2y*ppN2**atm_kappa
& )*ppl**(-atm_kappa)
& + a1tmp*(1. _d 0 /ppm - 1. _d 0/ppN1)
& + (ttp4y-ttp2y)*(pph-ppN2)/(ppN4-ppN2)
& + (ttdC+yCtmp)*MAX(0. _d 0,pNlev(k)-yDtmp)/(1-yDtmp)
ENDDO
C---------------------------------------------------
C matlab script, input: pN, yp=abs(yLat)
C pN0=.1; pN1=.19 ; pN2=.3; pN4=.925;
C pm=min(max(pN,pN1),pN2); pp=max(pN,pN2);
C pl=min(pN,pN1); pr=(pN0+abs(pl-pN0)-pN1)/(pN2-pN1);
C
C yA=yp-45; yA=45+min(0,yA/.9)+max(0,yA); yA=max(0,yA); cosyA=cos(yA*rad) ;
C yB=yp-65; yB=65+yB/.9; yB=min(max(0,yB),90); cosyB=cos(yB*rad) ;
C tp1=350*ones(nyA,1);
C tp2=307+(342-307)*cosyA.^2.6;
C tp4=257+(301-257)*cosyB.^1.5;
C a1=(tp1-tp2)*pN1*pN2/(pN2-pN1);
C pF=max(0,1.-((yp-65)/10).^2); pT=.9+(0.7-.9)*pF;
C
C tA0=((1-pr(k))*tp1(j)*pN1^kappa+pr(k)*tp2(j)*pN2^kappa)*pl(k)^-kappa;
C tA1=a1(j)*(1./pm(k)-1./pN1);
C tA2=(tp4(j)-tp2(j))*(pp(k)-pN2)/(pN4-pN2);
C tA3=(3.3+pF(j))*max(0,pN(k)-pT(j))/(1-pT(j));
C tAn(j,k)=tA0+tA1+tA2+tA3;
C---------------------------------------------------
RETURN
END