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