Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:37:19 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
d676f916b2 Jean*0001 #include "AIM_OPTIONS.h"
                0002 
54052ec14b Jean*0003 CBOP
                0004 C     !ROUTINE: AIM_DYN2AIM
                0005 C     !INTERFACE:
d676f916b2 Jean*0006       SUBROUTINE AIM_DYN2AIM(
                0007      O           TA, QA, ThA, Vsurf2, PSA, dpFac,
                0008      O           kGrd,
                0009      I           bi,bj, myTime, myIter, myThid)
54052ec14b Jean*0010 C     !DESCRIPTION: \bv
d676f916b2 Jean*0011 C     *==========================================================*
                0012 C     | S/R AIM_DYN2AIM
                0013 C     | o Map dynamics conforming arrays to AIM internal arrays.
                0014 C     *==========================================================*
                0015 C     | this routine transfers grid information
                0016 C     | and all dynamical variables (T,Q, ...) to AIM physics
                0017 C     *==========================================================*
54052ec14b Jean*0018 C     \ev
d676f916b2 Jean*0019 
54052ec14b Jean*0020 C     !USES:
                0021       IMPLICIT NONE
                0022 C     === Global variables ===
d676f916b2 Jean*0023 C-- size for MITgcm & Physics package :
54052ec14b Jean*0024 #include "AIM_SIZE.h"
d676f916b2 Jean*0025 
                0026 #include "EEPARAMS.h"
                0027 #include "PARAMS.h"
                0028 #include "GRID.h"
                0029 #include "SURFACE.h"
                0030 #include "DYNVARS.h"
                0031 
                0032 #include "AIM_GRID.h"
                0033 #include "com_physcon.h"
                0034 
54052ec14b Jean*0035 C     !INPUT/OUTPUT PARAMETERS:
d676f916b2 Jean*0036 C--   Input:
54052ec14b Jean*0037 C     bi,bj  :: Tile index
                0038 C     myTime :: Current time of simulation ( s )
                0039 C     myIter :: Current iteration number in simulation
                0040 C     myThid :: Number of this instance of the routine
                0041 C--   Output:  TA     :: temperature  [K}                        (3-dim)
                0042 C              QA     :: specific humidity [g/kg]                (3-dim)
                0043 C              ThA    :: Pot.Temp. [K] (replace dry stat. energy)(3-dim)
                0044 C              Vsurf2 :: square of surface wind speed            (2-dim)
                0045 C              PSA    :: norm. surface pressure [p/p0]           (2-dim)
                0046 C              dpFac  :: cell delta_P fraction                   (3-dim)
                0047 C              kGrd   :: Ground level index                      (2-dim)
                0048 C--  Updated common blocks: AIM_GRID_R
                0049 C             WVSurf  :: weights for near surf interpolation     (2-dim)
                0050 C             fOrogr  :: orographic factor (for surface drag)    (2-dim)
                0051 C         snLat,csLat :: sin(Lat) & cos(Lat)                     (2-dim)
d676f916b2 Jean*0052       _RL TA(NGP,NLEV), QA(NGP,NLEV), ThA(NGP,NLEV)
                0053       _RL Vsurf2(NGP), PSA(NGP), dpFac(NGP,NLEV)
                0054       INTEGER kGrd(NGP)
54052ec14b Jean*0055       INTEGER bi, bj, myIter, myThid
                0056       _RL myTime
                0057 CEOP
d676f916b2 Jean*0058 
                0059 #ifdef ALLOW_AIM
54052ec14b Jean*0060 C     !LOCAL VARIABLES:
                0061 C     i, j, I2 :: Loop counters
                0062 C     k, Katm  :: Loop counters
                0063 C     msgBuf   :: Informational/error message buffer
d676f916b2 Jean*0064       CHARACTER*(MAX_LEN_MBUF) msgBuf
54052ec14b Jean*0065       INTEGER i, j, I2, k, Katm
d676f916b2 Jean*0066       _RL conv_theta2T, temp1, temp2
                0067 
                0068 c     _RL hInitC(5), hInitF(5)
54052ec14b Jean*0069 c     DATA    hInitC / 17338.0,10090.02,5296.88,2038.54,418.038/
                0070 c     DATA    hInitF / 15090.4, 8050.96, 4087.75, 1657.54, 0. /
d676f916b2 Jean*0071 
                0072 C-    Compute Sin & Cos (Latitude) :
54052ec14b Jean*0073       DO j = 1,sNy
                0074        DO i = 1,sNx
                0075         I2 = i+(j-1)*sNx
                0076         snLat(I2,myThid) = SIN(yC(i,j,bi,bj)*deg2rad)
                0077         csLat(I2,myThid) = COS(yC(i,j,bi,bj)*deg2rad)
d676f916b2 Jean*0078        ENDDO
                0079       ENDDO
                0080 
                0081 C-    Set surface level index :
54052ec14b Jean*0082       DO j = 1,sNy
                0083        DO i = 1,sNx
                0084         I2 = i+(j-1)*sNx
                0085         kGrd(I2) = (Nr+1) - kSurfC(i,j,bi,bj)
d676f916b2 Jean*0086        ENDDO
                0087       ENDDO
                0088 
                0089 C-    Set (normalized) surface pressure :
54052ec14b Jean*0090       DO j=1,sNy
                0091        DO i=1,sNx
                0092         I2 = i+(j-1)*sNx
                0093         k = kSurfC(i,j,bi,bj)
                0094         IF ( k.LE.Nr ) THEN
                0095 c        PSA(I2) = rF(k)/atm_Po
                0096          PSA(I2) = Ro_surf(i,j,bi,bj)/atm_Po
d676f916b2 Jean*0097         ELSE
                0098          PSA(I2) = 1.
                0099         ENDIF
                0100        ENDDO
                0101       ENDDO
                0102 
                0103 C-    Set cell delta_P fraction (of the full delta.P = drF_k):
26494fa735 Jean*0104 #ifdef NONLIN_FRSURF
                0105       IF ( staggerTimeStep .AND. nonlinFreeSurf.GT.0 ) THEN
                0106        IF ( select_rStar.GT.0 ) THEN
                0107         DO k = 1,Nr
                0108          Katm = _KD2KA( k )
                0109          DO j = 1,sNy
                0110           DO i = 1,sNx
                0111            I2 = i+(j-1)*sNx
                0112            dpFac(I2,Katm) = h0FacC(i,j,k,bi,bj)*rStarFacC(i,j,bi,bj)
d676f916b2 Jean*0113 c          dpFac(I2,Katm) = 1. _d 0
26494fa735 Jean*0114           ENDDO
                0115          ENDDO
d676f916b2 Jean*0116         ENDDO
26494fa735 Jean*0117        ELSE
                0118         DO k = 1,Nr
                0119          Katm = _KD2KA( k )
                0120          DO j = 1,sNy
                0121           DO i = 1,sNx
                0122            I2 = i+(j-1)*sNx
54052ec14b Jean*0123            IF ( k.EQ.kSurfC(i,j,bi,bj) ) THEN
26494fa735 Jean*0124             dpFac(I2,Katm) = hFac_surfC(i,j,bi,bj)
                0125            ELSE
                0126             dpFac(I2,Katm) = hFacC(i,j,k,bi,bj)
                0127            ENDIF
                0128 c          dpFac(I2,Katm) = 1. _d 0
                0129           ENDDO
                0130          ENDDO
                0131         ENDDO
                0132        ENDIF
                0133       ELSE
                0134 #else /* ndef NONLIN_FRSURF */
                0135       IF (.TRUE.) THEN
                0136 #endif /* NONLIN_FRSURF */
                0137         DO k = 1,Nr
                0138          Katm = _KD2KA( k )
                0139          DO j = 1,sNy
                0140           DO i = 1,sNx
                0141            I2 = i+(j-1)*sNx
                0142            dpFac(I2,Katm) = hFacC(i,j,k,bi,bj)
                0143 c          dpFac(I2,Katm) = 1. _d 0
                0144           ENDDO
                0145          ENDDO
                0146         ENDDO
                0147       ENDIF
d676f916b2 Jean*0148 
                0149 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0150 
                0151 C     Physics package works with sub-domains 1:sNx,1:sNy,1:Nr.
                0152 C     Internal index mapping is linear in X and Y with a second
                0153 C     dimension for the vertical.
                0154 
                0155 C-    Dynamical var --> AIM var :
                0156 C       note: UA & VA are not used  => removed
                0157       temp1 = lwTemp1
                0158       temp2 = lwTemp2
54052ec14b Jean*0159       DO k = 1,Nr
                0160        conv_theta2T = (rC(k)/atm_Po)**atm_kappa
                0161        Katm = _KD2KA( k )
                0162        DO j = 1,sNy
                0163         DO i = 1,sNx
                0164          I2 = i+(j-1)*sNx
                0165          IF ( maskC(i,j,k,bi,bj).EQ.oneRS ) THEN
                0166 c         UA(I2,Katm)  = uVel(i,j,k,bi,bj)
                0167 c         VA(I2,Katm)  = vVel(i,j,k,bi,bj)
d676f916b2 Jean*0168 C     Physics works with temperature - not potential temp.
54052ec14b Jean*0169           TA(I2,Katm)  = theta(i,j,k,bi,bj)*conv_theta2T
d676f916b2 Jean*0170 c         TA(I2,Katm)  = max(temp1,min(temp2,
54052ec14b Jean*0171 c    &                       theta(i,j,k,bi,bj)*conv_theta2T ))
                0172 C     In atm.Phys, water vapor must be > 0 :
                0173           QA(I2,Katm)  = MAX( salt(i,j,k,bi,bj), zeroRL )
d676f916b2 Jean*0174 C     Dry static energy replaced by Pot.Temp:
54052ec14b Jean*0175           ThA(I2,Katm) = theta(i,j,k,bi,bj)
d676f916b2 Jean*0176          ELSE
                0177           TA(I2,Katm)  = 300. _d 0
                0178           QA(I2,Katm)  =   0. _d 0
                0179           ThA(I2,Katm) = 300. _d 0
                0180          ENDIF
                0181         ENDDO
                0182        ENDDO
54052ec14b Jean*0183 #ifdef NONLIN_FRSURF
                0184        IF ( select_rStar.GE.1 ) THEN
                0185         DO j = 1,sNy
                0186          DO i = 1,sNx
                0187           I2 = i+(j-1)*sNx
                0188           TA(I2,Katm)  = TA(I2,Katm)*pStarFacK(i,j,bi,bj)
                0189          ENDDO
                0190         ENDDO
                0191        ENDIF
                0192 #endif /* NONLIN_FRSURF */
d676f916b2 Jean*0193       ENDDO
                0194 
                0195 C_jmc: add square of surface wind speed (center of C grid) = 2 * KE_surf
54052ec14b Jean*0196       DO j = 1,sNy
                0197        DO i = 1,sNx
                0198         I2 = i+(j-1)*sNx
                0199         k = kSurfC(i,j,bi,bj)
                0200         IF (k.LE.Nr) THEN
d676f916b2 Jean*0201          Vsurf2(I2) = 0.5 * (
54052ec14b Jean*0202      &                uVel(i,j,k,bi,bj)*uVel(i,j,k,bi,bj)
                0203      &              + uVel(i+1,j,k,bi,bj)*uVel(i+1,j,k,bi,bj)
                0204      &              + vVel(i,j,k,bi,bj)*vVel(i,j,k,bi,bj)
                0205      &              + vVel(i,j+1,k,bi,bj)*vVel(i,j+1,k,bi,bj)
d676f916b2 Jean*0206      &                        )
                0207         ELSE
                0208          Vsurf2(I2) = 0.
                0209         ENDIF
                0210        ENDDO
                0211       ENDDO
                0212 
                0213 C-    Check that Temp is OK for LW Radiation scheme :
54052ec14b Jean*0214       DO k = 1,Nr
                0215        Katm = _KD2KA( k )
d676f916b2 Jean*0216        DO I2=1,NGP
                0217         IF (  TA(I2,Katm).LT.lwTemp1 .OR.
                0218      &        TA(I2,Katm).GT.lwTemp2 ) THEN
                0219          i = 1 + mod((I2-1),sNx)
                0220          j = 1 + int((I2-1)/sNx)
                0221          WRITE(msgBuf,'(A,1PE20.13,A,2I4)')
                0222      &    'AIM_DYN2AIM: Temp=', TA(I2,Katm),
                0223      &    ' out of range ',lwTemp1,lwTemp2
                0224          CALL PRINT_ERROR( msgBuf , myThid)
                0225          WRITE(msgBuf,'(A,3I4,3I3,I6,2F9.3)')
                0226      &    'AIM_DYN2AIM: Pb in i,j,k,bi,bj,myThid,I2,X,Y=',
                0227      &        i,j,k,bi,bj,myThid,I2,xC(i,j,bi,bj),yC(i,j,bi,bj)
                0228          CALL PRINT_ERROR( msgBuf , myThid)
54052ec14b Jean*0229          STOP 'ABNORMAL END: S/R AIM_DYN2AIM'
d676f916b2 Jean*0230         ENDIF
                0231        ENDDO
                0232       ENDDO
                0233 
                0234 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0235 
                0236 C-    Set geopotential surfaces
                0237 c     DO Katm=1,NLEV
                0238 c      DO I2=1,NGP
                0239 c       PHIG1(I2,Katm) = gravity*HinitC(Katm)
                0240 c      ENDDO
                0241 c     ENDDO
                0242 
                0243 C-    Weights for vertical interpolation down to the surface
                0244 C     Fsurf = Ffull(nlev)+WVS*(Ffull(nlev)-Ffull(nlev-1))
54052ec14b Jean*0245       DO j = 1,sNy
                0246        DO i = 1,sNx
                0247         I2 = i+(j-1)*sNx
d676f916b2 Jean*0248         WVSurf(I2,myThid) = 0.
54052ec14b Jean*0249         k = kGrd(I2)
                0250         IF (k.GT.1) THEN
cf526db02d Jean*0251 C- full cell version of Franco Molteni formula:
54052ec14b Jean*0252 c         WVSurf(I2,myThid) = (LOG(SIGH(k))-SIGL(k))*WVI(k-1,2)
cf526db02d Jean*0253 C- partial cell version using true log-P extrapolation:
54052ec14b Jean*0254           WVSurf(I2,myThid) = (LOG(PSA(I2))-SIGL(k))*WVI(k-1,1)
d676f916b2 Jean*0255 C- like in the old code:
54052ec14b Jean*0256 c         WVSurf(I2,myThid) = WVI(k,2)
d676f916b2 Jean*0257         ENDIF
                0258        ENDDO
                0259       ENDDO
54052ec14b Jean*0260       IF (myIter.EQ.nIter0)
7f9e3ec2e7 Jean*0261      &  CALL AIM_WRITE_PHYS( 'aim_WeightSurf', '', 1, WVSurf,
                0262      &                       1, bi, bj, 1, myIter, myThid )
d676f916b2 Jean*0263 
                0264 #endif /* ALLOW_AIM */
                0265 
                0266       RETURN
                0267       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!