Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:36:52 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
6c9bcec054 Jean*0001 #include "CPP_OPTIONS.h"
463053c692 Jean*0002 #undef CHECK_ANALYLIC_THETA
6c9bcec054 Jean*0003 
9366854e02 Chri*0004 CBOP
                0005 C     !ROUTINE: INI_P_GROUND
                0006 C     !INTERFACE:
f15994caab Jean*0007       SUBROUTINE INI_P_GROUND(selectMode,
                0008      &                        Hfld, Pfld,
6c9bcec054 Jean*0009      I                        myThid )
9366854e02 Chri*0010 C     !DESCRIPTION: \bv
                0011 C     *==========================================================*
f15994caab Jean*0012 C     | SUBROUTINE INI_P_GROUND
                0013 C     | o Convert Topography [m] to (reference) Surface Pressure
                0014 C     |   according to tRef profile,
                0015 C     |   using same discretisation as in calc_phi_hyd
                0016 C     |
9366854e02 Chri*0017 C     *==========================================================*
                0018 C     \ev
                0019 
                0020 C     !USES:
6c9bcec054 Jean*0021       IMPLICIT NONE
                0022 C     == Global variables ==
                0023 #include "SIZE.h"
                0024 #include "GRID.h"
                0025 #include "EEPARAMS.h"
                0026 #include "PARAMS.h"
463053c692 Jean*0027 #include "SURFACE.h"
6c9bcec054 Jean*0028 
9366854e02 Chri*0029 C     !INPUT/OUTPUT PARAMETERS:
6c9bcec054 Jean*0030 C     == Routine arguments ==
46ac67bdf3 Jean*0031 C     selectMode ::  > 0 = find Pfld from Hfld ; < 0 = compute Hfld from Pfld
                0032 C                   selectFindRoSurf = 0 : use Theta_Ref profile
                0033 C                   selectFindRoSurf = 1 : use analytic fct Theta(Lat,P)
f15994caab Jean*0034 C     Hfld (input/outp) :: Ground elevation [m]
                0035 C     Pfld (outp/input) :: reference Pressure at the ground [Pa]
463053c692 Jean*0036       INTEGER selectMode
6c9bcec054 Jean*0037       _RS Hfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0038       _RS Pfld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0039       INTEGER myThid
                0040 
9366854e02 Chri*0041 C     !LOCAL VARIABLES:
6c9bcec054 Jean*0042 C     == Local variables ==
f15994caab Jean*0043 C     msgBuf :: Informational/error message buffer
463053c692 Jean*0044 C-
f15994caab Jean*0045 C       For an accurate definition of the reference surface pressure,
463053c692 Jean*0046 C       define a High vertical resolution (in P):
                0047 C     nLevHvR :: Number of P-level used for High vertical Resolution (HvR)
                0048 C     plowHvR :: Lower bound of the High vertical Resolution
                0049 C     dpHvR   :: normalized pressure increment (HvR)
                0050 C     pLevHvR :: normalized P-level of the High vertical Resolution
                0051 C     pMidHvR :: normalized mid point level (HvR)
                0052 C     thetaHvR :: potential temperature at mid point level (HvR)
                0053 C     PiHvR  :: Exner function at P-level
                0054 C     dPiHvR :: Exner function difference between 2 P-levels
                0055 C     recip_kappa :: 1/kappa = Cp/R
                0056 C     PiLoc, zLoc, dzLoc, yLatLoc, phiLoc :: hold on temporary values
                0057 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0058       CHARACTER*(MAX_LEN_MBUF) msgBuf
893c662779 Jean*0059       INTEGER bi,bj,i,j,k, ks
797ab70dc1 Jean*0060       _RL Po_surf
769c28d1ef Jean*0061       _RL hRef(2*Nr+1), rHalf(2*Nr+1)
46ac67bdf3 Jean*0062       LOGICAL findPoSurf
463053c692 Jean*0063 
                0064       INTEGER nLevHvR
                0065       PARAMETER ( nLevHvR = 60 )
                0066       _RL plowHvR, dpHvR, pLevHvR(nLevHvR+1), pMidHvR(nLevHvR)
                0067       _RL thetaHvR(nLevHvR), PiHvR(nLevHvR+1), dPiHvR(nLevHvR)
                0068       _RL recip_kappa, PiLoc, zLoc, dzLoc, yLatLoc, phiLoc
f15994caab Jean*0069       _RL  psNorm, rMidKp1
46ac67bdf3 Jean*0070       _RL ratioRm(Nr), ratioRp(Nr)
463053c692 Jean*0071       INTEGER kLev
                0072 #ifdef CHECK_ANALYLIC_THETA
46ac67bdf3 Jean*0073       _RL tmpVar(nLevHvR,61)
463053c692 Jean*0074 #endif
9366854e02 Chri*0075 CEOP
6c9bcec054 Jean*0076 
463053c692 Jean*0077       IF ( selectFindRoSurf.LT.0 .OR. selectFindRoSurf.GT.1 ) THEN
                0078         WRITE(msgBuf,'(A,I2,A)')
                0079      &   'INI_P_GROUND: selectFindRoSurf =', selectFindRoSurf,
                0080      &        ' <== bad value !'
893c662779 Jean*0081         CALL PRINT_ERROR( msgBuf, myThid )
f15994caab Jean*0082         STOP 'INI_P_GROUND'
463053c692 Jean*0083       ENDIF
                0084 
893c662779 Jean*0085       DO k=1,Nr
                0086         rHalf(2*k-1) = rF(k)
                0087         rHalf(2*k)   = rC(k)
6c9bcec054 Jean*0088       ENDDO
                0089        rHalf(2*Nr+1) = rF(Nr+1)
                0090 
769c28d1ef Jean*0091 C- Reference Geopotential at Half levels :
893c662779 Jean*0092 C    Tracer level: hRef(2k)  ;  Interface_W level: hRef(2k+1)
6c9bcec054 Jean*0093 C- Convert phiRef to Z unit :
893c662779 Jean*0094       DO k=1,2*Nr+1
                0095         hRef(k) = phiRef(k)*recip_gravity
6c9bcec054 Jean*0096       ENDDO
463053c692 Jean*0097 
46ac67bdf3 Jean*0098       IF (selectFindRoSurf.EQ.0 .AND. selectMode .GT. 0 ) THEN
6c9bcec054 Jean*0099 C- Find Po_surf : Linear between 2 half levels :
                0100       DO bj = myByLo(myThid), myByHi(myThid)
                0101        DO bi = myBxLo(myThid), myBxHi(myThid)
                0102         DO j=1,sNy
                0103          DO i=1,sNx
893c662779 Jean*0104            ks = 1
                0105            DO k=1,2*Nr
                0106              IF (Hfld(i,j,bi,bj).GE.hRef(k)) ks = k
6c9bcec054 Jean*0107            ENDDO
893c662779 Jean*0108            Po_surf = rHalf(ks) + (rHalf(ks+1)-rHalf(ks))*
                0109      &       (Hfld(i,j,bi,bj)-hRef(ks))/(hRef(ks+1)-hRef(ks))
6c9bcec054 Jean*0110 
769c28d1ef Jean*0111 c          IF (Hfld(i,j,bi,bj).LT.hRef(1)) Po_surf= rHalf(1)
                0112 c          IF (Hfld(i,j,bi,bj).GT.hRef(2*Nr+1)) Po_surf=rHalf(2*Nr+1)
6c9bcec054 Jean*0113            Pfld(i,j,bi,bj) = Po_surf
                0114          ENDDO
                0115         ENDDO
                0116        ENDDO
                0117       ENDDO
                0118 
                0119 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
46ac67bdf3 Jean*0120 C-- endif selectFindRoSurf=0 & selectMode > 0
463053c692 Jean*0121       ENDIF
                0122 
46ac67bdf3 Jean*0123       IF ( selectFindRoSurf.EQ.1 ) THEN
463053c692 Jean*0124 C-- define high resolution Pressure discretization:
                0125 
                0126       recip_kappa = 1. _d 0 / atm_kappa
                0127       plowHvR = 0.4 _d 0
                0128       dpHvR = nLevHvR
                0129       dpHvR = (1. - plowHvR) / dpHvR
2ad79bdf32 Jean*0130         pLevHvR(1)= rF(1)/atm_Po
463053c692 Jean*0131         PiHvR(1) = atm_Cp*(pLevHvR(1)**atm_kappa)
                0132       DO k=1,nLevHvR
893c662779 Jean*0133         pLevHvR(k+1)= pLevHvR(1) - k*dpHvR
463053c692 Jean*0134         PiHvR(k+1)  = atm_Cp*(pLevHvR(k+1)**atm_kappa)
f15994caab Jean*0135         pMidHvR(k)= (pLevHvR(k)+pLevHvR(k+1))*0.5 _d 0
463053c692 Jean*0136         dPiHvR(k) = PiHvR(k) - PiHvR(k+1)
                0137       ENDDO
                0138 
46ac67bdf3 Jean*0139 C--   to modify pressure when using non fully linear relation between Phi & p
                0140 C       (Integr_GeoPot=2 & Tracer Point at middle between 2 interfaces)
                0141       DO k=1,Nr
                0142          ratioRm(k) = 1. _d 0
                0143          ratioRp(k) = 1. _d 0
                0144          IF (k.GT.1 ) ratioRm(k) = 0.5 _d 0*drC(k)/(rF(k)-rC(k))
f15994caab Jean*0145          IF (k.LT.Nr) ratioRp(k) = 0.5 _d 0*drC(k+1)/(rC(k)-rF(k+1))
46ac67bdf3 Jean*0146       ENDDO
                0147 
463053c692 Jean*0148 #ifdef CHECK_ANALYLIC_THETA
893c662779 Jean*0149       _BEGIN_MASTER( myThid )
463053c692 Jean*0150       DO j=1,61
                0151         yLatLoc =-90.+(j-1)*3.
893c662779 Jean*0152         CALL ANALYLIC_THETA( yLatLoc, pMidHvR,
                0153      &                       tmpVar(1,j), nLevHvR, myThid )
f15994caab Jean*0154       ENDDO
463053c692 Jean*0155       OPEN(88,FILE='analytic_theta',
                0156      &      STATUS='unknown',FORM='unformatted')
                0157       WRITE(88) tmpVar
                0158       CLOSE(88)
893c662779 Jean*0159       _END_MASTER( myThid )
463053c692 Jean*0160       STOP 'CHECK_ANALYLIC_THETA'
                0161 #endif /* CHECK_ANALYLIC_THETA */
                0162 
46ac67bdf3 Jean*0163 C-- endif selectFindRoSurf=1
463053c692 Jean*0164       ENDIF
                0165 
46ac67bdf3 Jean*0166       IF ( selectFindRoSurf*selectMode .GT. 0) THEN
463053c692 Jean*0167 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0168 C- Find Po_surf such as g*Hfld = Phi[Po_surf,theta(yLat,p)]:
                0169 
                0170       DO bj = myByLo(myThid), myByHi(myThid)
                0171        DO bi = myBxLo(myThid), myBxHi(myThid)
                0172 C- start bi,bj loop:
46ac67bdf3 Jean*0173 
463053c692 Jean*0174         DO j=1,sNy
                0175          DO i=1,sNx
893c662779 Jean*0176           phiLoc = Hfld(i,j,bi,bj) - hRef(1)
                0177           IF ( phiLoc .LE. 0. _d 0 ) THEN
2ad79bdf32 Jean*0178            Pfld(i,j,bi,bj) = rF(1)
463053c692 Jean*0179           ELSE
                0180            yLatLoc  = yC(i,j,bi,bj)
893c662779 Jean*0181            CALL ANALYLIC_THETA( yLatLoc, pMidHvR,
                0182      &                          thetaHvR, nLevHvR, myThid )
463053c692 Jean*0183            zLoc = 0.
f15994caab Jean*0184            DO k=1,nLevHvR
463053c692 Jean*0185             IF (zLoc.GE.0.) THEN
                0186 C-    compute  phi/g corresponding to next p_level:
f15994caab Jean*0187              dzLoc = dPiHvR(k)*thetaHvR(k)*recip_gravity
893c662779 Jean*0188              IF ( phiLoc .LE. zLoc+dzLoc ) THEN
46ac67bdf3 Jean*0189 C-    compute the normalized surf. Pressure psNorm
f15994caab Jean*0190                PiLoc = PiHvR(k)
893c662779 Jean*0191      &               - gravity*( phiLoc - zLoc )/thetaHvR(k)
f15994caab Jean*0192                psNorm = (PiLoc/atm_Cp)**recip_kappa
463053c692 Jean*0193 C- use linear interpolation:
f15994caab Jean*0194 c              psNorm = pLevHvR(k)
893c662779 Jean*0195 c    &                - dpHvR*( phiLoc - zLoc )/dzLoc
463053c692 Jean*0196                zLoc = -1.
                0197              ELSE
f15994caab Jean*0198                zLoc = zLoc + dzLoc
463053c692 Jean*0199              ENDIF
                0200             ENDIF
                0201            ENDDO
                0202            IF (zLoc.GE.0.) THEN
                0203              WRITE(msgBuf,'(2A)')
                0204      &        'INI_P_GROUND: FAIL in trying to find Pfld:',
                0205      &        ' selectMode,i,j,bi,bj=',selectMode,i,j,bi,bj
893c662779 Jean*0206              CALL PRINT_ERROR( msgBuf, myThid )
463053c692 Jean*0207              WRITE(msgBuf,'(A,F7.1,2A,F6.4,A,F8.0)')
                0208      &        'INI_P_GROUND: Hfld=', Hfld(i,j,bi,bj), ' exceeds',
                0209      &        ' Zloc(lowest P=', pLevHvR(1+nLevHvR),' )=',zLoc
893c662779 Jean*0210              CALL PRINT_ERROR( msgBuf, myThid )
f15994caab Jean*0211              STOP 'ABNORMAL END: S/R INI_P_GROUND'
463053c692 Jean*0212            ELSE
46ac67bdf3 Jean*0213              Pfld(i,j,bi,bj) = psNorm*atm_Po
463053c692 Jean*0214            ENDIF
                0215           ENDIF
                0216          ENDDO
                0217         ENDDO
46ac67bdf3 Jean*0218 
                0219         IF (selectMode.EQ.2 .AND. integr_GeoPot.NE.1) THEN
                0220 C---------
f15994caab Jean*0221 C     Modify pressure to account for non fully linear relation between
46ac67bdf3 Jean*0222 C      Phi & p (due to numerical truncation of the Finite Difference scheme)
                0223 C---------
                0224           DO j=1,sNy
                0225            DO i=1,sNx
                0226              Po_surf = Pfld(i,j,bi,bj)
                0227               IF ( Po_surf.LT.rC(1) .AND. Po_surf.GT.rC(Nr) ) THEN
                0228                 findPoSurf = .TRUE.
                0229                 DO k=1,Nr
                0230                   IF ( findPoSurf .AND. Po_surf.GE.rC(k) ) THEN
                0231                     Po_surf = rC(k) + (Po_surf-rC(k))/ratioRm(k)
                0232                     findPoSurf = .FALSE.
                0233                   ENDIF
                0234                   rMidKp1 = rF(k+1)
f15994caab Jean*0235                   IF (k.LT.Nr) rMidKp1 = (rC(k)+rC(k+1))*0.5 _d 0
46ac67bdf3 Jean*0236                   IF ( findPoSurf .AND. Po_surf.GE.rMidKp1 ) THEN
                0237                     Po_surf = rC(k) + (Po_surf-rC(k))/ratioRp(k)
                0238                     findPoSurf = .FALSE.
                0239                   ENDIF
                0240                 ENDDO
f15994caab Jean*0241                 IF ( findPoSurf )
46ac67bdf3 Jean*0242      &               STOP 'S/R INI_P_GROUND: Pb with selectMode=2'
f15994caab Jean*0243               ENDIF
46ac67bdf3 Jean*0244              Pfld(i,j,bi,bj) = Po_surf
                0245            ENDDO
                0246           ENDDO
                0247 C---------
                0248         ENDIF
                0249 
463053c692 Jean*0250 C- end bi,bj loop.
                0251        ENDDO
                0252       ENDDO
                0253 
                0254 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
46ac67bdf3 Jean*0255 C-- endif selectFindRoSurf*selectMode > 0
463053c692 Jean*0256       ENDIF
                0257 
46ac67bdf3 Jean*0258       IF (selectMode .LT. 0) THEN
463053c692 Jean*0259 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
46ac67bdf3 Jean*0260 C--- Compute Hfld = Phi(Pfld)/g.
463053c692 Jean*0261 
                0262       DO bj = myByLo(myThid), myByHi(myThid)
                0263        DO bi = myBxLo(myThid), myBxHi(myThid)
                0264 C- start bi,bj loop:
                0265 
769c28d1ef Jean*0266 C--  Compute Hfld from g*Hfld = hRef(Po_surf)
463053c692 Jean*0267         DO j=1,sNy
                0268          DO i=1,sNx
769c28d1ef Jean*0269 C-   compute phiLoc = hRef(Ro_surf):
f15994caab Jean*0270           ks = kSurfC(i,j,bi,bj)
463053c692 Jean*0271           IF (ks.LE.Nr) THEN
f15994caab Jean*0272 C-   sigma coord. case (=> uniform kSurfC), find true ks:
                0273            IF ( selectSigmaCoord.NE.0 ) THEN
                0274              DO k=2,Nr
                0275                IF ( Pfld(i,j,bi,bj).LT.rF(k) ) ks = k
                0276              ENDDO
                0277            ENDIF
463053c692 Jean*0278            IF ( Pfld(i,j,bi,bj).GE.rC(ks) ) THEN
769c28d1ef Jean*0279             phiLoc = hRef(2*ks)
                0280      &       + (hRef(2*ks-1)-hRef(2*ks))
463053c692 Jean*0281      &        *(Pfld(i,j,bi,bj)-rC(ks))/(rHalf(2*ks-1)-rHalf(2*ks))
                0282            ELSE
769c28d1ef Jean*0283             phiLoc = hRef(2*ks)
                0284      &       + (hRef(2*ks+1)-hRef(2*ks))
463053c692 Jean*0285      &        *(Pfld(i,j,bi,bj)-rC(ks))/(rHalf(2*ks+1)-rHalf(2*ks))
                0286            ENDIF
46ac67bdf3 Jean*0287            Hfld(i,j,bi,bj) = phiLoc
463053c692 Jean*0288           ELSE
                0289            Hfld(i,j,bi,bj) = 0.
                0290           ENDIF
                0291          ENDDO
                0292         ENDDO
46ac67bdf3 Jean*0293 
                0294         IF (selectFindRoSurf.EQ.1) THEN
                0295 C-----
ff02675122 Jean*0296 C  goal: evaluate phi0surf (used for computing geopotential_prime = Phi - PhiRef):
46ac67bdf3 Jean*0297 C   phi0surf = g*Ztopo-PhiRef(Ro_surf) if no truncation was applied to Ro_surf;
                0298 C  but because of hFacMin, surf.ref.pressure (=Ro_surf) is truncated and
                0299 C   phi0surf = Phi(Theta-Analytic,P=Ro_surf) - PhiRef(P=Ro_surf)
                0300 C-----
                0301 C--   Compute Hfld from g*Hfld = Phi[Po_surf,theta(yLat,p)]:
                0302          DO j=1,sNy
                0303           DO i=1,sNx
893c662779 Jean*0304            zLoc = hRef(1)
2ad79bdf32 Jean*0305            IF ( Pfld(i,j,bi,bj) .LT. rF(1) ) THEN
46ac67bdf3 Jean*0306             Po_surf = Pfld(i,j,bi,bj)
                0307 C---------
f15994caab Jean*0308 C     Modify pressure to account for non fully linear relation between
46ac67bdf3 Jean*0309 C      Phi & p (due to numerical truncation of the Finite Difference scheme)
                0310              IF (selectMode.EQ.-2 .AND. integr_GeoPot.NE.1) THEN
                0311               IF ( Po_surf.LT.rC(1) .AND. Po_surf.GT.rC(Nr) ) THEN
                0312                 findPoSurf = .TRUE.
                0313                 DO k=1,Nr
                0314                   IF ( findPoSurf .AND. Po_surf.GE.rC(k) ) THEN
                0315                     Po_surf = rC(k) + (Po_surf-rC(k))*ratioRm(k)
                0316                     findPoSurf = .FALSE.
                0317                   ENDIF
                0318                   IF ( findPoSurf .AND. Po_surf.GE.rF(k+1) ) THEN
                0319                     Po_surf = rC(k) + (Po_surf-rC(k))*ratioRp(k)
                0320                     findPoSurf = .FALSE.
                0321                   ENDIF
                0322                 ENDDO
f15994caab Jean*0323               ENDIF
                0324              ENDIF
46ac67bdf3 Jean*0325 C---------
                0326             psNorm = Po_surf/atm_Po
                0327             kLev = 1 + INT( (pLevHvR(1)-psNorm)/dpHvR )
                0328             yLatLoc  = yC(i,j,bi,bj)
893c662779 Jean*0329             CALL ANALYLIC_THETA( yLatLoc, pMidHvR,
                0330      &                           thetaHvR, kLev, myThid )
46ac67bdf3 Jean*0331 C-    compute height at level pLev(kLev) just below Pfld:
                0332             DO k=1,kLev-1
f15994caab Jean*0333               dzLoc = dPiHvR(k)*thetaHvR(k)*recip_gravity
46ac67bdf3 Jean*0334               zLoc = zLoc + dzLoc
                0335             ENDDO
                0336             dzLoc = ( PiHvR(kLev)-atm_Cp*(psNorm**atm_kappa) )
f15994caab Jean*0337      &            * thetaHvR(kLev)*recip_gravity
46ac67bdf3 Jean*0338             zLoc = zLoc + dzLoc
                0339            ENDIF
                0340 C-    compute phi0surf = Phi[Po_surf,theta(yLat,p)] - PhiRef(Po_surf)
                0341            phi0surf(i,j,bi,bj) = gravity*(zLoc - Hfld(i,j,bi,bj))
                0342 C-    save Phi[Po_surf,theta(yLat,p)] in Hfld (replacing PhiRef(Po_surf)):
                0343            Hfld(i,j,bi,bj) = zLoc
                0344           ENDDO
                0345          ENDDO
                0346 C- endif selectFindRoSurf=1
                0347         ENDIF
                0348 
463053c692 Jean*0349 C- end bi,bj loop.
                0350        ENDDO
                0351       ENDDO
                0352 
                0353 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
46ac67bdf3 Jean*0354 C-- endif selectMode < 0
463053c692 Jean*0355       ENDIF
                0356 
                0357       RETURN
                0358       END
                0359 
                0360 CBOP
                0361 C     !SUBROUTINE: ANALYLIC_THETA
                0362 C     !INTERFACE:
                0363       SUBROUTINE ANALYLIC_THETA( yLat, pNlev,
                0364      O                           thetaLev,
893c662779 Jean*0365      I                           kSize, myThid )
463053c692 Jean*0366 
                0367 C     !DESCRIPTION: \bv
                0368 C     *==========================================================*
f15994caab Jean*0369 C     | SUBROUTINE ANALYLIC_THETA
                0370 C     | o Conpute analyticaly the potential temperature Theta
                0371 C     |   as a function of Latitude (yLat) and normalized
                0372 C     |   pressure pNlev.
463053c692 Jean*0373 C     |   approximatively match the N-S symetric, zonal-mean and
                0374 C     |   annual-mean NCEP climatology in the troposphere.
                0375 C     *==========================================================*
                0376 C     \ev
                0377 
                0378 C     !USES:
                0379       IMPLICIT NONE
                0380 
                0381 C     == Global variables ==
                0382 #include "SIZE.h"
                0383 #include "EEPARAMS.h"
                0384 #include "PARAMS.h"
                0385 
                0386 C     !INPUT/OUTPUT PARAMETERS:
                0387 C     == Routine arguments ==
                0388 C     yLat   :: latitude (degre)
                0389 C     pNlev  :: normalized pressure levels
                0390 C     kSize  :: dimension of pNlev & ANALYLIC_THETA
                0391 C     myThid :: Thread number for this instance of the routine
                0392       INTEGER kSize
                0393       _RL  yLat
                0394       _RL  pNlev  (kSize)
                0395       _RL  thetaLev(kSize)
                0396       INTEGER myThid
                0397 
                0398 C     !LOCAL VARIABLES:
                0399 C     == Local variables ==
                0400       INTEGER k
                0401       _RL  yyA, yyB, yyC, yyAd, yyBd, yyCd
                0402       _RL  cAtmp, cBtmp, ttdC
f15994caab Jean*0403       _RL  ppN0, ppN1, ppN2, ppN3a, ppN3b, ppN4
463053c692 Jean*0404       _RL  ttp1, ttp2, ttp3, ttp4, ttp5
                0405       _RL  yAtmp, yBtmp, yCtmp, yDtmp
f15994caab Jean*0406       _RL  ttp2y, ttp4y, a1tmp
463053c692 Jean*0407       _RL  ppl, ppm, pph, ppr
                0408 CEOP
                0409 
                0410       DATA yyA ,    yyB ,     yyC ,     yyAd ,   yyBd ,   yyCd
                0411      &  / 45. _d 0, 65. _d 0, 65. _d 0, .9 _d 0, .9 _d 0, 10. _d 0 /
                0412       DATA  cAtmp ,   cBtmp ,   ttdC
                0413      &   /  2.6 _d 0, 1.5 _d 0, 3.3 _d 0 /
f15994caab Jean*0414       DATA  ppN0  ,   ppN1  ,  ppN2  ,  ppN3a ,  ppN3b ,  ppN4
463053c692 Jean*0415      &   / .1 _d 0, .19 _d 0, .3 _d 0, .9 _d 0, .7 _d 0, .925 _d 0 /
                0416       DATA ttp1 ,     ttp2 ,     ttp3 ,     ttp4 ,     ttp5
                0417      &   / 350. _d 0, 342. _d 0, 307. _d 0, 301. _d 0, 257. _d 0 /
                0418 
                0419 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0420 
                0421        yAtmp = ABS(yLat) - yyA
                0422        yAtmp = yyA + MIN(0. _d 0,yAtmp/yyAd) + MAX(yAtmp, 0. _d 0)
                0423        yAtmp = COS( deg2rad*MAX(yAtmp, 0. _d 0) )
                0424        yBtmp = ABS(yLat) - yyB
                0425        yBtmp = yyB + yBtmp/yyBd
                0426        yBtmp = COS( deg2rad*MAX( 0. _d 0, MIN(yBtmp,90. _d 0) ) )
                0427        yCtmp = ABS(yLat) - yyC
                0428        yCtmp = MAX( 0. _d 0, 1. _d 0 - (yCtmp/yyCd)**2 )
                0429        yDtmp = ppN3a +(ppN3b - ppN3a)*yCtmp
                0430        ttp2y = ttp3 + (ttp2-ttp3)*yAtmp**cAtmp
                0431        ttp4y = ttp5 + (ttp4-ttp5)*yBtmp**cBtmp
                0432        a1tmp = (ttp1-ttp2y)*ppN1*ppN2/(ppN2-ppN1)
                0433       DO k=1,kSize
                0434        ppl = MIN(pNlev(k),ppN1)
                0435        ppm = MIN(MAX(pNlev(k),ppN1),ppN2)
                0436        pph = MAX(pNlev(k),ppN2)
                0437        ppr =( ppN0 + ABS(ppl-ppN0) - ppN1 )/(ppN2-ppN1)
f15994caab Jean*0438        thetaLev(k) =
463053c692 Jean*0439      &       ( (1. _d 0 -ppr)*ttp1*ppN1**atm_kappa
                0440      &        + ppr*ttp2y*ppN2**atm_kappa
                0441      &       )*ppl**(-atm_kappa)
                0442      &     + a1tmp*(1. _d 0 /ppm - 1. _d 0/ppN1)
f15994caab Jean*0443      &     + (ttp4y-ttp2y)*(pph-ppN2)/(ppN4-ppN2)
463053c692 Jean*0444      &     + (ttdC+yCtmp)*MAX(0. _d 0,pNlev(k)-yDtmp)/(1-yDtmp)
                0445       ENDDO
6c9bcec054 Jean*0446 
f15994caab Jean*0447 C---------------------------------------------------
463053c692 Jean*0448 C matlab script, input: pN, yp=abs(yLat)
                0449 C pN0=.1; pN1=.19 ; pN2=.3; pN4=.925;
                0450 C pm=min(max(pN,pN1),pN2); pp=max(pN,pN2);
                0451 C pl=min(pN,pN1); pr=(pN0+abs(pl-pN0)-pN1)/(pN2-pN1);
f15994caab Jean*0452 C
463053c692 Jean*0453 C  yA=yp-45; yA=45+min(0,yA/.9)+max(0,yA); yA=max(0,yA); cosyA=cos(yA*rad) ;
                0454 C  yB=yp-65; yB=65+yB/.9; yB=min(max(0,yB),90); cosyB=cos(yB*rad) ;
                0455 C  tp1=350*ones(nyA,1);
                0456 C  tp2=307+(342-307)*cosyA.^2.6;
                0457 C  tp4=257+(301-257)*cosyB.^1.5;
                0458 C  a1=(tp1-tp2)*pN1*pN2/(pN2-pN1);
                0459 C  pF=max(0,1.-((yp-65)/10).^2); pT=.9+(0.7-.9)*pF;
f15994caab Jean*0460 C
463053c692 Jean*0461 C  tA0=((1-pr(k))*tp1(j)*pN1^kappa+pr(k)*tp2(j)*pN2^kappa)*pl(k)^-kappa;
                0462 C  tA1=a1(j)*(1./pm(k)-1./pN1);
                0463 C  tA2=(tp4(j)-tp2(j))*(pp(k)-pN2)/(pN4-pN2);
                0464 C  tA3=(3.3+pF(j))*max(0,pN(k)-pT(j))/(1-pT(j));
                0465 C  tAn(j,k)=tA0+tA1+tA2+tA3;
f15994caab Jean*0466 C---------------------------------------------------
6c9bcec054 Jean*0467 
                0468       RETURN
                0469       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!