pkg/aim_v23/phy_convmf.F

d676f916b2 Jean*0001 #include "AIM_OPTIONS.h"
                0002
                0003       SUBROUTINE CONVMF (PSA,dpFac,SE,QA,QSAT,
                0004      O                   IDEPTH,CBMF,PRECNV,DFSE,DFQA,
                0005      I                   kGrd,bi,bj,myThid)
                0006 C--
                0007 C--   SUBROUTINE CONVMF (PSA,SE,QA,QSAT,
                0008 C--  *                   IDEPTH,CBMF,PRECNV,DFSE,DFQA)
                0009 C--
                0010 C--   Purpose: Compute convective fluxes of dry static energy and moisture
                0011 C--            using a simplified mass-flux scheme
                0012 C--   Input:   PSA    = norm. surface pressure [p/p0]            (2-dim)
                0013 C              dpFac  = cell delta_P fraction                    (3-dim)
                0014 C--            SE     = dry static energy                        (3-dim)
                0015 C--            QA     = specific humidity [g/kg]                 (3-dim)
                0016 C--            QSAT   = saturation spec. hum. [g/kg]             (3-dim)
                0017 C--   Output:  IDEPTH = convection depth in layers               (2-dim)
                0018 C--            CBMF   = cloud-base mass flux                     (2-dim)
                0019 C--            PRECNV = convective precipitation [g/(m^2 s)]     (2-dim)
                0020 C--            DFSE   = net flux of d.s.en. into each atm. layer (3-dim)
                0021 C--            DFQA   = net flux of sp.hum. into each atm. layer (3-dim)
                0022 C    Input:    kGrd   = Ground level index                       (2-dim)
                0023 C              bi,bj  = tile index
                0024 C              myThid = Thread number for this instance of the routine
                0025 C-------
                0026 C  Note: dry static energy has been replaced by Pot.Temp.
                0027 C-------
                0028
                0029       IMPLICIT NONE
                0030
                0031 C     Resolution parameters
                0032
                0033 C-- size for MITgcm & Physics package :
                0034 #include "AIM_SIZE.h"
                0035
                0036 #include "EEPARAMS.h"
                0037
                0038 C     Physical constants + functions of sigma and latitude
                0039
                0040 #include "com_physcon.h"
                0041
                0042 C     Convection constants
                0043
                0044 #include "com_cnvcon.h"
                0045
                0046 C-- Routine arguments:
                0047       _RL PSA(NGP), SE(NGP,NLEV), QA(NGP,NLEV), QSAT(NGP,NLEV)
                0048       _RL dpFac(NGP,NLEV)
                0049       INTEGER IDEPTH(NGP)
                0050       _RL CBMF(NGP), PRECNV(NGP), DFSE(NGP,NLEV), DFQA(NGP,NLEV)
                0051       INTEGER  kGrd(NGP)
                0052       INTEGER  bi,bj,myThid
                0053
                0054 #ifdef ALLOW_AIM
                0055
                0056 C-- Local variables:
                0057       INTEGER ITOP(NGP)
                0058 c_FM  REAL SM(NGP,NLEV), QATHR(NGP), ENTR(2:NLEV-1)
                0059       _RL  QATHR(NGP), ENTR(2:NLEV-1)
                0060       _RL  ENTR_PS(NGP,2:NLEV-1), FM0(NGP)
                0061
                0062       INTEGER J, K, K1, Ktmp
                0063       _RL  dSEdp(NGP,NLEV), factP, PSA_1
                0064       _RL  dSEdpTot, stab_crit, FDMUS
                0065 C- jmc: declare all local variables:
                0066       _RL  QMAX, DELQ, QB, QSATB, FMASS, ENMASS, SENTR
                0067       _RL  FPSA, FQMAX, RDPS, FUQ, FDQ, FSQ
                0068 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0069
                0070 C--   1. Initialization of output and workspace arrays
                0071
                0072       PSA_1 = 1.
                0073       FQMAX=  5.
                0074
                0075       RDPS = 2. _d 0 /(1. _d 0 - PSMIN)
                0076 c_FM  FM0=P0*DSIG(NLEV)/(GG*TRCNV*3600)
                0077 c_FM    FPSA=PSA(J)*MIN(1.,(PSA(J)-PSMIN)*RDPS)
                0078 C-    compute FM0(J) = FM0*FPSA
                0079       DO J=1,NGP
                0080        FM0(J)=0.
                0081        Ktmp = kGrd(J)
                0082        IF ( Ktmp .NE. 0 ) THEN
                0083         FPSA = MIN(1. _d 0 ,(PSA(J)-PSMIN)*RDPS)
                0084         FM0(J)=P0*DSIG(Ktmp)*dpFac(J,Ktmp)/(GG*TRCNV*3600. _d 0)
                0085        ENDIF
                0086       ENDDO
                0087
                0088       DO K=1,NLEV
                0089         DO J=1,NGP
                0090           DFSE(J,K)=0.0
                0091           DFQA(J,K)=0.0
                0092         ENDDO
                0093       ENDDO
                0094       DO K=2,NLEV-1
                0095         DO J=1,NGP
                0096           ENTR_PS(J,K)=0.
                0097         ENDDO
                0098       ENDDO
                0099
                0100       DO J=1,NGP
                0101         ITOP(J)  =kGrd(J)
                0102         CBMF(J)  =0.0
                0103         PRECNV(J)=0.0
                0104       ENDDO
                0105
                0106 C     Saturation moist static energy
                0107 c_FM  DO K=1,NLEV
                0108 c_FM    DO J=1,NGP
                0109 c_FM      SM(J,K)=SE(J,K)+ALHC*QSAT(J,K)
                0110 c_FM    ENDDO
                0111 c_FM  ENDDO
                0112
                0113 C ---------------------------------------------
                0114 C    Write Conditional stability based on Pot.Temp :
                0115 C    dSEdp(K) = Delta[Static-Energy] between 2 Plevels(k,k+1)
                0116 C    and corresponds to SE(K+1)-SE(K) in the original code
                0117 C -------
                0118       DO K=1,NLEV-1
                0119        factP = CP*SIGH(K)**(RD/CP)
                0120        DO J=1,NGP
                0121          dSEdp(J,K)=(SE(J,K+1)-SE(J,K))*factP
                0122        ENDDO
                0123       ENDDO
                0124
                0125 C     Entrainment profile (up to sigma = 0.5)
                0126
                0127 c_FM  SENTR=0.
                0128 c_FM  DO K=2,NLEV-1
                0129 c_FM    ENTR(K)=( MAX( 0. _d 0, SIG(K)-0.5 _d 0) )**2
                0130 c_FM    SENTR=SENTR+ENTR(K)
                0131 c_FM  ENDDO
                0132
                0133 c_FM  SENTR=ENTMAX/SENTR
                0134 c_FM  DO K=2,NLEV-1
                0135 c_FM    ENTR(K)=ENTR(K)*SENTR
                0136 c_FM  ENDDO
                0137
                0138       DO J=1,NGP
                0139        DO K=2,NLEV-1
                0140         ENTR_PS(J,K)=0.
                0141        ENDDO
                0142        Ktmp = kGrd(J)
                0143        IF (Ktmp.GT.2) THEN
                0144          SENTR=0.
                0145          DO K=2,Ktmp-1
                0146            ENTR(K)= ( MAX( 0. _d 0, SIG(K)/PSA(J) - 0.5 _d 0) )**2
                0147            SENTR=SENTR+ENTR(K)
                0148          ENDDO
                0149          IF (SENTR.GT.0.) THEN
                0150           SENTR=ENTMAX/SENTR
                0151           DO K=2,Ktmp-1
                0152            ENTR_PS(J,K) = ENTR(K)*SENTR*PSA(J)
                0153           ENDDO
                0154          ENDIF
                0155        ENDIF
                0156       ENDDO
                0157
                0158 C--   2. Check of conditions for convection
                0159
                0160 C     2.1 Conditional instability
                0161
                0162 c_FM  DO K=NLEV-2,2,-1
                0163 c_FM    DO J=1,NGP
                0164 c_FM      SMB=SM(J,K)+WVI(K,2)*(SM(J,K+1)-SM(J,K))
                0165 c_FM      IF (SM(J,NLEV).GT.SMB) ITOP(J)=K
                0166 c_FM    ENDDO
                0167 c_FM  ENDDO
                0168
                0169       DO J=1,NGP
                0170        Ktmp = kGrd(J)
                0171        IF ( Ktmp .GE. 2 ) THEN
                0172         dSEdpTot = dSEdp(J,Ktmp-1)
                0173         DO k=Ktmp-2,2,-1
                0174           dSEdpTot = dSEdpTot + dSEdp(J,K)
                0175           stab_crit = dSEdpTot + ALHC*(QSAT(J,Ktmp)-QSAT(J,K))
                0176      &     -WVI(K,2)*(dSEdp(J,K) + ALHC*(QSAT(J,K+1)-QSAT(J,K)) )
                0177           IF (stab_crit.GT.0.) ITOP(J) = K
                0178         ENDDO
                0179        ENDIF
                0180       ENDDO
                0181
                0182
                0183 C     2.2 Humidity exceeding prescribed threshold
                0184
                0185       DO J=1,NGP
                0186        Ktmp = kGrd(J)
                0187        IF ( Ktmp .NE. 0 ) THEN
                0188         QATHR(J)=MIN(QBL,RHBL*QSAT(J,Ktmp))
                0189         IF (QA(J,Ktmp).LT.QATHR(J).OR.PSA(J).LT.PSMIN)
                0190      &      ITOP(J)=Ktmp
                0191        ENDIF
                0192         IDEPTH(J)=Ktmp-ITOP(J)
                0193       ENDDO
                0194
                0195 C--   3. Convection over selected grid-points
                0196
                0197       DO 300 J=1,NGP
                0198        Ktmp = kGrd(J)
                0199       IF (ITOP(J).EQ.Ktmp) GO TO 300
                0200
                0201 C-      3.1 Boundary layer (cloud base)
                0202
                0203         K = Ktmp
                0204         K1=K-1
                0205
                0206 C       Maximum specific humidity in the PBL
                0207         QMAX=MAX(1.01 _d 0 *QA(J,K),QSAT(J,K))
                0208
                0209 C       Dry static energy and moisture at upper boundary
                0210 c_FM    SB=SE(J,K1)+WVI(K1,2)*(SE(J,K)-SE(J,K1))
                0211         QB=QA(J,K1)+WVI(K1,2)*(QA(J,K)-QA(J,K1))
                0212         QB=MIN(QB,QA(J,K))
                0213
                0214 C       Cloud-base mass flux, computed to satisfy:
                0215 C       fmass*(qmax-qb)*(g/dp)=(q-qthr)/trcnv
                0216
                0217 c_FM    FPSA=PSA(J)*MIN(1.,(PSA(J)-PSMIN)*RDPS)
                0218 c_FM    FMASS=FM0*FPSA*MIN(FQMAX,(QA(J,K)-QATHR(J))/(QMAX-QB))
                0219         FMASS = FM0(J)*MIN(FQMAX,(QA(J,K)-QATHR(J))/(QMAX-QB))
                0220         CBMF(J)=FMASS
                0221
                0222 C       Upward fluxes at upper boundary
                0223 c_FM    FUS=FMASS*SE(J,K)
                0224         FUQ=FMASS*QMAX
                0225
                0226 C       Downward fluxes at upper boundary
                0227 c_FM    FDS=FMASS*SB
                0228         FDQ=FMASS*QB
                0229
                0230 C       Net flux of dry static energy and moisture
                0231         FDMUS = FMASS*dSEdp(J,K1)*(WVI(K1,2)-1.)
                0232         DFSE(J,K)=FDMUS
                0233 c_FM    DFSE(J,K)=FDS-FUS
                0234         DFQA(J,K)=FDQ-FUQ
                0235
                0236 C-      3.2 Intermediate layers (entrainment)
                0237
                0238         DO K=Ktmp-1,ITOP(J)+1,-1
                0239         K1=K-1
                0240
                0241 C         Fluxes at lower boundary
                0242 c_FM      DFSE(J,K)=FUS-FDS
                0243           DFQA(J,K)=FUQ-FDQ
                0244
                0245 C         Mass entrainment
                0246 c_FM      ENMASS=ENTR(K)*PSA(J)*CBMF(J)
                0247           ENMASS=ENTR_PS(J,K) * CBMF(J)
                0248           FMASS=FMASS+ENMASS
                0249
                0250 C         Upward fluxes at upper boundary
                0251 c_FM      FUS=FUS+ENMASS*SE(J,K)
                0252           FUQ=FUQ+ENMASS*QA(J,K)
                0253
                0254 C         Downward fluxes at upper boundary
                0255 c_FM      SB=SE(J,K1)+WVI(K1,2)*(SE(J,K)-SE(J,K1))
                0256           QB=QA(J,K1)+WVI(K1,2)*(QA(J,K)-QA(J,K1))
                0257 c_FM      FDS=FMASS*SB
                0258           FDQ=FMASS*QB
                0259
                0260 C         Net flux of dry static energy and moisture
                0261           DFSE(J,K) = FMASS*(WVI(K1,2)-1.)*dSEdp(J,K1)
                0262      &             -(FMASS-ENMASS)*WVI(K,2)*dSEdp(J,K)
                0263           FDMUS = FDMUS + DFSE(J,K)
                0264 c_FM      DFSE(J,K)=DFSE(J,K)+FDS-FUS
                0265           DFQA(J,K)=DFQA(J,K)+FDQ-FUQ
                0266
                0267 C         Secondary moisture flux
                0268           DELQ=RHIL*QSAT(J,K)-QA(J,K)
                0269           IF (DELQ.GT.0.0) THEN
                0270             FSQ=SMF*CBMF(J)*DELQ
                0271             DFQA(J,K)   =DFQA(J,K)   +FSQ
                0272             DFQA(J,Ktmp)=DFQA(J,Ktmp)-FSQ
                0273           ENDIF
                0274
                0275         ENDDO
                0276
                0277 C-      3.3 Top layer (condensation and detrainment)
                0278
                0279         K=ITOP(J)
                0280
                0281 C       Flux of convective precipitation
                0282         QSATB=QSAT(J,K)+WVI(K,2)*(QSAT(J,K+1)-QSAT(J,K))
                0283         PRECNV(J)=MAX(FUQ-FMASS*QSATB, 0. _d 0)
                0284
                0285 C       Net flux of dry static energy and moisture
                0286         DFSE(J,K)= -FDMUS+ALHC*PRECNV(J)
                0287 c_FM    DFSE(J,K)=FUS-FDS+ALHC*PRECNV(J)
                0288         DFQA(J,K)=FUQ-FDQ-PRECNV(J)
                0289
                0290  300  CONTINUE
                0291
                0292 #endif /* ALLOW_AIM */
                0293
                0294       RETURN
                0295       END