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