! $Id: ras.F90,v 1.24.2.8.2.1.4.1 2007/11/09 20:46:45 bacmj Exp $
MODULE RAS 1
use GEOS_Mod
use ESMF_Mod
use GEOS_UtilsMod, only : DQSAT=>GEOS_DQsat
IMPLICIT NONE
PRIVATE
PUBLIC RASE
CONTAINS
SUBROUTINE RASE(IDIM, IRUN, K0, ICMIN, DT , & 1,8
CPO,ALHLO,ALHL1,TICE,GRAVO, &
SEEDRAS,IRAS,JRAS,SIGE, &
KCBL,WGT0,WGT1,ZCBL,TPERT,QPERT, &
THO, QHO, UHO, VHO, &
QSS, DQS, &
PLE, PKE, CLW, FLX, FLXD, FLXC, &
CNV_PRC3, &
CNV_UPDFRC, &
CNV_QC, &
RASPARAMS, &
ITRCR, XHO, IRC, &
TRIEDLEV_DIAG, &
FSCAV , DISSKE )
!*********************************************************************
!*********************************************************************
!******************** Relaxed Arakawa-Schubert ***********************
!************************ Parameterization ***************************
!********************** SCALAR RAS-1 VERSION ************************
!************************* 31 DECEMBER 1999 **************************
!*********************************************************************
!************************** Developed By *****************************
!*********************************************************************
!************************ Shrinivas Moorthi **************************
!******************************* and *********************************
!************************** Max J. Suarez ****************************
!*********************************************************************
!******************** Laboratory for Atmospheres *********************
!****************** NASA/GSFC, Greenbelt, MD 20771 *******************
!*********************************************************************
!*********************************************************************
! Input:
! ------
!
! K0 : Number of vertical levels (increasing downwards)
!
! DT : Time step in seconds
!
! RASAL : Array of dimension K-1 containing relaxation parameters
! for cloud-types detraining at those levels
!
! CPO : Specific heat at constant pressure (J/kg/K)
!
! ALHLO : Latent Heat of condensation (J/kg)
!
! ALHL1 : Latent Heat of condensation + fusion (J/kg)
!
! GRAVO : Acceleration due to gravity (m/s^2)
!
! PLE : 2D array of dimension (IDIM,K0+1) containing pressure
! in hPa at the interfaces of K-layers from top of the
! atmosphere to the bottom (mb)
!
! PKE : 2D array of dimension (IDIM,K0+1) containing (PRS/P00) **
! RKAP. i.e. Exner function at layer edges.
!
! PKL : 2D array of dimension (IDIM,K0) ) containing the
! Exner function at the layers.
!
! QSS : 2D array of dimension (IDIM,K0 ) containing the
! saturation specific humidity at the layers. (kg/kg)
!
! DQS : 2D array of dimension (IDIM,K0 ) containing
! d(qss)/dt at the layers. (1/K)
!
! Update:
! -------
!
! THO : 2D array of dimension (IDIM,K0) containing potential
! temperature (K)
!
! QHO : 2D array of dimension (IDIM,K0) containing specific
! humidity (kg/kg)
!
! UHO : 2D array of dimension (IDIM,K0) containing u-wind (m/s)
!
! VHO : 2D array of dimension (IDIM,K0) containing v-wind (m/s)
!
! Output:
! -------
!!
! CLW : 2D array of dimension (IDIM,K0) containing the
! detrained cloud liquid water. (kg/m^2/s)
!
! FLX : 2D array of dimension (IDIM,K0) containing the
! cloud-base mass flux for each cloud type ordered by
! detrainment level. (kg/m^2/s)
!
! FLXD : 2D array of dimension (IDIM,K0) containing the
! detrained mass flux for each cloud type ordered by
! detrainment level. (kg/m^2/s)
!
! FLXC : 2D array of dimension (IDIM,K0) containing the
! total cloud mass flux for all cloud types through
! the top of each level. (e.g., FLXC(K)=SUM(FLX(ICMIN:K))
! and FLXD(L) = FLXC(L+1)-FLSD(L) )
! (kg/m^2/s)
!
!************************************************************************
! ARGUMENTS
INTEGER, INTENT(IN ) :: IDIM, IRUN, K0, ICMIN
REAL, DIMENSION (IDIM,K0 ), INTENT(INOUT) :: THO, QHO, UHO, VHO
REAL, DIMENSION (IDIM,K0+1), INTENT(IN ) :: PLE, PKE
REAL, DIMENSION (IDIM,K0 ), INTENT(IN ) :: QSS, DQS
REAL, DIMENSION ( K0+1), INTENT(IN ) :: SIGE
REAL, DIMENSION (IDIM,K0 ), INTENT( OUT) :: FLX, CLW , FLXD, FLXC
REAL, DIMENSION (IDIM,K0 ), INTENT( OUT) :: CNV_PRC3
REAL, DIMENSION (IDIM,K0 ), INTENT( OUT) :: CNV_UPDFRC, CNV_QC
REAL, INTENT(IN ) :: DT, CPO, ALHLO, GRAVO
REAL, INTENT(IN ) :: ALHL1, TICE
INTEGER, INTENT( IN) :: ITRCR
INTEGER, DIMENSION ( 2 ), INTENT(IN ) :: SEEDRAS
INTEGER, DIMENSION (IDIM ), INTENT(IN ) :: IRAS,JRAS,KCBL
REAL, DIMENSION (IDIM ), INTENT(IN ) :: ZCBL,TPERT,QPERT
REAL, DIMENSION (K0 ), INTENT(IN ) :: WGT0,WGT1
REAL, DIMENSION(:), INTENT(IN ) :: RASPARAMS
REAL, OPTIONAL, INTENT(INOUT) :: XHO(IDIM,K0,ITRCR)
INTEGER, OPTIONAL, INTENT( OUT) :: IRC(IDIM,K0)
REAL, OPTIONAL, INTENT( OUT) :: TRIEDLEV_DIAG(IDIM,K0)
REAL, OPTIONAL, INTENT( OUT) :: DISSKE(IDIM,K0)
REAL, OPTIONAL, INTENT(IN) :: FSCAV(ITRCR) ! Fraction scavenged per km
! = 0 (no scav), = 1 (full scav)
! LOCALS
REAL, DIMENSION(K0) :: POI_SV, QOI_SV, UOI_SV, VOI_SV
REAL, DIMENSION(K0) :: POI, QOI, UOI, VOI, DQQ, BET, GAM, CLL
REAL, DIMENSION(K0) :: POI_c, QOI_c
REAL, DIMENSION(K0) :: PRH, PRI, GHT, DPT, DPB, PKI, DISSK0,DISSK1
REAL, DIMENSION(K0) :: TCU, QCU, UCU, VCU, CLN, RNS, POL
REAL, DIMENSION(K0) :: QST, SSL, RMF, RNN, RN1, RMFC, RMFP
REAL, DIMENSION(K0) :: GMS, ETA, GMH, EHT, GM1, HCC, RMFD
REAL, DIMENSION(K0) :: HOL, HST, QOL, ZOL, HCLD, CLL0,CLLX,CLLI,CLLB
REAL, DIMENSION(K0) :: WSP, LAMBDSV, BKE , CVW, UPDFRC
REAL, DIMENSION(K0) :: RASAL, MTKWI, UPDFRP,BK2,BK3,DLL0,DLLX
REAL, DIMENSION(ITRCR) :: XHT
REAL, DIMENSION(K0,ITRCR) :: XOI, XCU, XOI_SV
REAL, DIMENSION(K0+1) :: PRJ, PRS, QHT, SHT ,ZET, XYD, XYD0
INTEGER, DIMENSION(K0-1) :: RC
INTEGER :: K,MY_PE
REAL, DIMENSION(IDIM,K0) :: LAMBDSV2
REAL TX2, TX3, UHT, VHT, AKM, ACR, ALM, TTH, QQH, SHTRG, WSPBL, DQX
REAL WFN, TEM, TRG, TRGEXP, EVP, PCU, WLQ, QCC,MTKW_MAX !, BKE
REAL SHTRG_FAC, SIGE_MINHOL, WFNOG
INTEGER I, IC, L, ICL , ITR , ICL_C, N_DTL
INTEGER NDTLEXPON
INTEGER, ALLOCATABLE, DIMENSION(:) :: ICL_V
! RASE GLOBAL CONSTANTS
REAL GRAV, CP, ALHL, CPBG, ALHI, CPI, GRAVI, DDT, LBCP, OBG, AFC
!!!!!!!!!
REAL CO_AUTO, FRICFAC,DPTH_BL,WUPDRFT,PBLFRAC,AUTORAMPB,CO_ZDEP
REAL RASAL1, RASAL2, CO_T, RASNCL,FRICLAMBDA,SDQVT1,SDQV2
REAL LAMBDA_FAC,STRAPPING,ACRITFAC,HMINTRIGGER,LLDISAGGXP
REAL LAMBMX_FAC, DIAMMN_MIN,RDTLEXPON, CLI_CRIT,SDQV3
INTEGER KSTRAP
real cld_radius, areal_frac, spect_mflx, cvw_cbase
!!!!!!!!!
REAL, PARAMETER :: ONEPKAP = 1.+ 2./7., DAYLEN = 86400.0
! REAL, PARAMETER :: PBLFRAC = 0.5
REAL, PARAMETER :: RHMAX = 0.9999
! LAMBDA LIMITS
REAL :: LAMBDA_MIN
REAL :: LAMBDA_MAX
! TRIGGER PARAMETERS
REAL, PARAMETER :: RHMN = .5, RHMX = .65
real, parameter :: RHO_W = 1.0e3 ! Density of liquid water in kg/m^3
LOGICAL DYNA_STRAPPING, DO_TRACERS, SMOOTH_HST
character(len=ESMF_MAXSTR) :: CBL_STYLE
! SCAVANGING RELATED PARAMETERS
REAL :: DELZKM ! layer thickness in km
REAL :: FNOSCAV ! fraction of tracer *not* scavenged
REAL :: FSCAV_(ITRCR) ! Fraction scavenged per km
! *********************************************************************
IF ( PRESENT(FSCAV) ) then
FSCAV_ = FSCAV
ELSE
FSCAV_ = 0.0 ! NO SCAVENGING BY DEFAULT
ENDIF
IF(IRUN <= 0) RETURN
IF (PRESENT(TRIEDLEV_DIAG)) TRIEDLEV_DIAG=0.
CNV_PRC3 =0.
CNV_UPDFRC=0.
CNV_QC =0.
if ( PRESENT( DISSKE )) DISSKE =0.
!!LAMBDSV2 = 0.
IF(PRESENT(IRC)) IRC = -2
! SMOOTH_HST = .TRUE.
SMOOTH_HST = .FALSE.
FRICFAC = RASPARAMS(1) ! --- 1
SHTRG_FAC = RASPARAMS(2) ! --- 2
CO_AUTO = RASPARAMS(3) ! --- 3
CLI_CRIT = RASPARAMS(4) ! --- 4
RASAL1 = RASPARAMS(5) ! --- 5
RASAL2 = RASPARAMS(6) ! --- 6
RASNCL = RASPARAMS(7) ! --- 7
LAMBDA_FAC = RASPARAMS(8) ! --- 8
LAMBMX_FAC = RASPARAMS(9) ! --- 9
DIAMMN_MIN = RASPARAMS(10) ! --- 10
FRICLAMBDA = RASPARAMS(11) ! --- 11
RDTLEXPON = RASPARAMS(12) ! --- 12
STRAPPING = RASPARAMS(13) ! --- 13
SDQV2 = RASPARAMS(14) ! --- 14
SDQV3 = RASPARAMS(15) ! --- 15
SDQVT1 = RASPARAMS(16) ! --- 16
ACRITFAC = RASPARAMS(17) ! --- 17
HMINTRIGGER = RASPARAMS(18) ! --- 18
LLDISAGGXP = RASPARAMS(19) ! --- 19
PBLFRAC = RASPARAMS(20) ! --- 20
AUTORAMPB = RASPARAMS(21) ! --- 21
CO_ZDEP = RASPARAMS(22) ! --- 22
IF ( STRAPPING <= 0.0 ) THEN
DYNA_STRAPPING = .TRUE.
ELSE
DYNA_STRAPPING = .FALSE.
KSTRAP = INT( STRAPPING )
ENDIF
DO_TRACERS = present(XHO) .and. ITRCR>0
WUPDRFT = 2.500
GRAV = GRAVO
ALHL = ALHLO
CP = CPO
CPI = 1.0/CP
ALHI = 1.0/ALHL
GRAVI = 1.0/GRAV
CPBG = CP*GRAVI
DDT = DAYLEN/DT
AFC = -1.04E-4*SQRT(DT*113.84)
LBCP = ALHL*CPI
OBG = 100.*GRAVI
DO I=1,IRUN
!!CALL FINDBASE
K = KCBL(I)
CALL FINDDTLS
IF ( K > 0 ) THEN
CALL STRAP( FINAL=0 )
call HTEST
DO ICL_C = 1,N_DTL
ICL = ICL_V( ICL_C )
IF (DO_TRACERS) XCU(ICMIN:,:) = 0.
IF (PRESENT(TRIEDLEV_DIAG)) TRIEDLEV_DIAG(I,ICL) = 1.
UCU(ICMIN:) = 0. ! This change makes cumulus friction
VCU(ICMIN:) = 0. ! correct.
IF( ICL > ICMIN ) THEN
CALL CLOUDE(ICL)
ENDIF
ENDDO
IF ( SUM( RMF(ICMIN:K) ) > 0.0 ) THEN
CALL RNEVP
CALL STRAP( FINAL=1 )
ELSE
CALL STRAP( FINAL=2 )
ENDIF
ELSE
CALL STRAP( FINAL=2 )
ENDIF
ENDDO
IF(ALLOCATED(ICL_V)) DEALLOCATE(ICL_V)
RETURN
CONTAINS
!*********************************************************************
SUBROUTINE CLOUDE(IC) 1,2
INTEGER, INTENT(IN ) :: IC
REAL :: DEEP_FACT,CU_DIAM,WSCALE
REAL :: CLI , TE_A, C00_X, CLI_CRIT_X, PETE, TOKI, GMHx, HSTx !, dQx
REAL :: DT_LYR, RATE, CVW_X, CLOSS, F2, F3, F4, F5
INTEGER :: K700
TRG = AMIN1(1.,(QOI(K)/QST(K)-RHMN)/(RHMX-RHMN))
F4 = MIN( 1.0, MAX( 0.0 , (AUTORAMPB-SIGE(IC))/0.2 ) ) ! F4 should ramp from 0 at SIG=AUTORAMPB
! to 1 at SIG=AUTORAMPB-0.2
if ( SIGE(IC) >= 0.5 ) then
F5 = 1.0
else
F5 = 1.0 - 2.*CO_ZDEP *( 0.5 - SIGE(IC) )
F5 = MAX( F5 , 0.0 )
endif
IF(TRG <= 1.0E-5) THEN ! TRIGGER =========>>
RC(IC) = 7
RETURN
ENDIF
! RECOMPUTE SOUNDING UP TO DETRAINMENT LEVEL
POI_c = POI
QOI_c = QOI
POI_c(K) = POI_c(K) + TPERT(I)
QOI_c(K) = QOI_c(K) + QPERT(I)
ZET(K+1) = 0.
SHT(K+1) = CP*POI_c(K)*PRJ(K+1)
DO L=K,IC,-1
QOL(L) = AMIN1(QST(L)*RHMAX,QOI_c(L))
QOL(L) = AMAX1( 0.000, QOL(L) ) ! GUARDIAN/NEG.PREC.
SSL(L) = CP*PRJ(L+1)*POI_c(L) + GRAV*ZET(L+1)
HOL(L) = SSL(L) + QOL(L)*ALHL
HST(L) = SSL(L) + QST(L)*ALHL
TEM = POI_c(L)*(PRJ(L+1)-PRJ(L))*CPBG
ZET(L) = ZET(L+1) + TEM
ZOL(L) = ZET(L+1) + (PRJ(L+1)-PRH(L))*POI_c(L)*CPBG
ENDDO
DO L=IC+1,K
TEM = (PRJ(L)-PRH(L-1))/(PRH(L)-PRH(L-1))
SHT(L) = SSL(L-1) + TEM*(SSL(L)-SSL(L-1))
QHT(L) = .5*(QOL(L)+QOL(L-1))
ENDDO
! SMOOTH HSTAR W/ 1-2-1 Filter
if ( SMOOTH_HST ) then
HSTx = HST(IC) ! save for later
DO L=K-1,IC+1,-1
HST(L) = 0.25*(HST(L+1)+HST(L-1))+0.5*HST(L)
END DO
DO L=IC,IC
HST(L) = 0.5*HST(L+1)+0.5*HST(L)
END DO
endif
! CALCULATE LAMBDA, ETA, AND WORKFUNCTION
LAMBDA_MIN = .2/(LAMBDA_FAC*DPTH_BL)
LAMBDA_MAX = .2/( MAX( LAMBMX_FAC*DPTH_BL , DIAMMN_MIN ) )
IF (HOL(K) <= HST(IC)) THEN ! CANNOT REACH IC LEVEL ======>>
RC(IC) = 1
RETURN
ENDIF
! LAMBDA CALCULATION: MS-A18
TEM = (HST(IC)-HOL(IC))*(ZOL(IC)-ZET(IC+1))
DO L=IC+1,K-1
TEM = TEM + (HST(IC)-HOL(L ))*(ZET(L )-ZET(L +1))
ENDDO
IF(TEM <= 0.0) THEN ! NO VALID LAMBDA ============>>
RC(IC) = 2
RETURN
ENDIF
ALM = (HOL(K)-HST(IC)) / TEM
IF(ALM > LAMBDA_MAX) THEN
RC(IC) = 3
RETURN
ENDIF
TOKI=1.0
IF(ALM < LAMBDA_MIN) THEN
TOKI = ( ALM/LAMBDA_MIN )**2
!RC(IC) = 6
!RETURN
ENDIF
!LAMBDSV(IC) = ALM
! ETA CALCULATION: MS-A2
DO L=IC+1,K
ETA(L) = 1.0 + ALM * (ZET(L )-ZET(K))
ENDDO
ETA(IC) = 1.0 + ALM * (ZOL(IC)-ZET(K))
! WORKFUNCTION CALCULATION: MS-A22
WFN = 0.0
HCC(K) = HOL(K)
DO L=K-1,IC+1,-1
HCC(L) = HCC(L+1) + (ETA(L) - ETA(L+1))*HOL(L)
TEM = HCC(L+1)*DPB(L) + HCC(L)*DPT(L)
EHT(L) = ETA(L+1)*DPB(L) + ETA(L)*DPT(L)
WFN = WFN + (TEM - EHT(L)*HST(L))*GAM(L)
ENDDO
HCC(IC) = HST(IC)*ETA(IC)
WFN = WFN + (HCC(IC+1)-HST(IC)*ETA(IC+1))*GAM(IC)*DPB(IC)
! VERTICAL VELOCITY/KE CALCULATION (ADDED 12/2001 JTB)
BK3(K) = 0.0
BK2(K) = 0.0
BKE(K) = 0.0
HCLD(K) = HOL(K)
DO L=K-1,IC,-1
HCLD(L) = ( ETA(L+1)*HCLD(L+1) + &
(ETA(L) - ETA(L+1))*HOL(L) ) / ETA(L)
TEM = (HCLD(L)-HST(L) )*(ZET(L)-ZET(L+1))/ (1.0+LBCP*DQQ(L))
BKE(L) = BKE(L+1) + GRAV * TEM / ( CP*PRJ(L+1)*POI(L) )
BK2(L) = BK2(L+1) + GRAV * MAX(TEM,0.0) / ( CP*PRJ(L+1)*POI(L) )
BK3(L) = BK3(L+1) + GRAV * MIN(TEM,0.0) / ( CP*PRJ(L+1)*POI(L) )
CVW(L) = SQRT( 2.0* MAX( BK2(L) , 0.0 ) )
ENDDO
CU_DIAM = 1000. ! 1.0 / ( 5.0*ALM )
! ALPHA CALCULATION
IF ( ZET(IC) < 2000. ) RASAL(IC) = RASAL1
IF ( ZET(IC) >= 2000. ) THEN
RASAL(IC) = RASAL1 + (RASAL2-RASAL1)*(ZET(IC) - 2000.)/8000.
ENDIF
RASAL(IC) = MIN( RASAL(IC) , 1.0e5 )
RASAL(IC) = DT / RASAL(IC)
! HEAL/REPAIR/RATIONALIZE STUPID CALCULATION OF CVW
! IT IS STUPID BECAUSE IT IGNORES A NUMBER OF IMPORTANT EFFECTS
! 1) CONVECTIVE PRESSURE PERT. COULD GO EITHER WAY + OR -
! 2) CONDENSATE LOADING/MECHANICAL DRAG
! 3) ENVIRONMENTAL CONDENSATE COULD GO +
! 4) SUBGRID SCALE CLOUD BASE ENHANCEMENT OF Q AND T
!
! SO, HERE AFTER NAIVELY USING ITS VALUE TO TEST FOR CLOUD TOP
! W, WE FLOOR IT AT 1 M/S. THIS ENSURES THAT IF A PLUME
! MAKES IT THROUGH ALL TEST CONDITIONS, I.E., GETS A MASS-FLUX
! IT WILL ALSO HAVE SOME DIAGNOSED "UPWARD MOTION"
CVW(IC:K) = MAX( CVW(IC:K) , 1.00 )
! NOTE THIS "CENTRALIZES" A KLUGE PRESENT IN OTHER LOCATIONS.
! CLEAN UP SOME TIME. -JTB 12/04/03
! TEST FOR CRITICAL WORK FUNCTION
CALL ACRITN(POL(IC), PRS(K), ACR)
IF(WFN <= ACR) THEN ! SUB-CRITICAL WORK FUNCTION ======>>
RC(IC) = 4
RETURN
ENDIF
! CLOUD TOP WATER AND MOMENTUM (TIMES ETA(IC)) MS-A16
IF (DO_TRACERS) THEN
DO ITR=1,ITRCR
XHT(ITR) = XOI(K,ITR)
END DO
END IF
WLQ = QOL(K)
UHT = UOI(K)
VHT = VOI(K)
RNN(K) = 0.
CLL0(K) = 0.
DO L=K-1,IC,-1
TEM = ETA(L) - ETA(L+1)
WLQ = WLQ + TEM * QOL(L)
UHT = UHT + TEM * UOI(L)
VHT = VHT + TEM * VOI(L)
IF (DO_TRACERS) THEN
DO ITR=1,ITRCR
! Scavenging -- fscav is the fraction scavenged per km
! delz is the layer thickness in km
! fnoscav is the fraction not scavenged
!
DELZKM = ( ZET(L) - ZET(L+1) ) / 1000.
FNOSCAV = 1. - FSCAV_(ITR) * DELZKM
FNOSCAV = max( min(FNOSCAV,1.), 0.) ! must be in [0,1]
! Before scanveging: XHT(ITR) = XHT(ITR) + TEM * XOI(L,ITR)
XHT(ITR) = XHT(ITR) + TEM * XOI(L,ITR) * FNOSCAV
END DO
END IF
!!!! How much condensate (CLI) is present here?
IF(L>IC) THEN
TX2 = 0.5*(QST(L)+QST(L-1))*ETA(L)
TX3 = 0.5*(HST(L)+HST(L-1))*ETA(L)
QCC = TX2 + GM1(L)*(HCC(L)-TX3)
CLL0(L) = (WLQ-QCC)
ELSE
CLL0(L) = (WLQ-QST(IC)*ETA(IC))
ENDIF
CLL0(L) = MAX( CLL0(L) , 0.00 )
CLI = CLL0(L) / ETA(L) ! condensate (kg/kg)
TE_A = POI(L)*PRH(L) ! Temperature (K)
CALL SUNDQ3_ICE( TE_A, SDQV2, SDQV3, SDQVT1, F2 , F3 )
C00_X = CO_AUTO * F2 * F3 * F4 ! * F5 ! F4 reduces AUTO for shallow clouds, F5 modifies auto for deep clouds
CLI_CRIT_X = CLI_CRIT / ( F2 * F3 )
RATE = C00_X * ( 1.0 - EXP( -(CLI)**2 / CLI_CRIT_X**2 ) )
CVW_X = MAX( CVW(L) , 1.00 ) ! Floor conv. vel. since we don't
! really trust it at low values
DT_LYR = ( ZET(L)-ZET(L+1) )/CVW_X ! l.h.s. DT_LYR => time in layer (L,L+1)
CLOSS = CLL0(L) * RATE * DT_LYR
CLOSS = MIN( CLOSS , CLL0(L) )
CLL0(L) = CLL0(L) - CLOSS
DLL0(L) = CLOSS
IF(CLOSS>0.) then
WLQ = WLQ - CLOSS
RNN(L) = CLOSS
else
RNN(L) = 0.
ENDIF
ENDDO
WLQ = WLQ - QST(IC)*ETA(IC)
! CALCULATE GAMMAS AND KERNEL
GMS(K) = (SHT(K)-SSL(K))*PRI(K) ! MS-A30 (W/O GRAV)
GMH(K) = GMS(K) + (QHT(K)-QOL(K))*PRI(K)*ALHL ! MS-A31 (W/O GRAV)
AKM = GMH(K)*GAM(K-1)*DPB(K-1) ! MS-A37 (W/O GRAV)
TX2 = GMH(K)
DO L=K-1,IC+1,-1
GMS(L) = ( ETA(L )*(SHT(L)-SSL(L )) &
+ ETA(L+1)*(SSL(L)-SHT(L+1)) ) *PRI(L)
GMH(L) = GMS(L) &
+ ( ETA(L )*(QHT(L)-QOL(L )) &
+ ETA(L+1)*(QOL(L)-QHT(L+1)) )*ALHL*PRI(L)
TX2 = TX2 + (ETA(L) - ETA(L+1)) * GMH(L)
AKM = AKM - GMS(L)*EHT(L)*PKI(L) + TX2*GHT(L)
ENDDO
GMS(IC) = ETA(IC+1)*(SSL(IC)-SHT(IC+1))*PRI(IC)
AKM = AKM - GMS(IC)*ETA(IC+1)*DPB(IC)*PKI(IC)
GMH(IC) = GMS(IC) + ( ETA(IC+1)*(QOL(IC)-QHT(IC+1))*ALHL &
+ ETA(IC )*(HST(IC)-HOL(IC )) )*PRI(IC)
if ( SMOOTH_HST ) then
GMHx = GMS(IC) + ( ETA(IC+1)*(QOL(IC)-QHT(IC+1))*ALHL &
+ ETA(IC )*(HSTx -HOL(IC )) )*PRI(IC)
endif
! CLOUD BASE MASS FLUX
IF (AKM >= 0.0 .OR. WLQ < 0.0) THEN ! =========>
RC(IC) = 5
RETURN
ENDIF
WFN = - (WFN-ACR)/AKM ! MS-A39 MASS-FLUX IN Pa/step
WFN = MIN( ( RASAL(IC)*TRG*TOKI )*WFN , (PRS(K+1)-PRS(K) )*(100.*PBLFRAC))
! CUMULATIVE PRECIP AND CLOUD-BASE MASS FLUX FOR OUTPUT
WFNOG = WFN*GRAVI
TEM = WFN*GRAVI
CLL (IC) = CLL (IC) + WLQ*TEM ! (kg/m^2/step)
RMF (IC) = RMF (IC) + TEM ! (kg/m^2/step)
RMFD(IC) = RMFD(IC) + TEM * ETA(IC) ! (kg/m^2/step)
DO L=IC+1,K
RMFP(L) = TEM * ETA(L) ! (kg/m^2/step)
RMFC(L) = RMFC(L) + RMFP(L) ! (kg/m^2/step)
DLLX(L) = DLLX(L) + TEM*DLL0(L)
IF ( CVW(L) > 0.0 ) THEN
UPDFRP(L) = rmfp(L) * (DDT/DAYLEN) * 1000. / ( CVW(L) * PRS(L) )
ELSE
UPDFRP(L) = 0.0
ENDIF
CLLI(L) = CLL0(L)/ETA(L) ! current cloud; incloud condensate
CLLB(L) = CLLB(L) + UPDFRP(L) * CLLI(L) ! cumulative grid mean convective condensate
UPDFRC(L) = UPDFRC(L) + UPDFRP(L)
ENDDO
! THETA AND Q CHANGE DUE TO CLOUD TYPE IC
DO L=IC,K
RNS(L) = RNS(L) + RNN(L)*TEM ! (kg/m^2/step)
GMH(L) = GMH(L) * WFN
GMS(L) = GMS(L) * WFN
QOI(L) = QOI(L) + (GMH(L) - GMS(L)) * ALHI
POI(L) = POI(L) + GMS(L)*PKI(L)*CPI
QST(L) = QST(L) + GMS(L)*BET(L)*CPI
ENDDO
if ( SMOOTH_HST ) then
GMHx = GMHx * WFN
dQx = (GMHx - GMH(IC)) * ALHI
RNS(IC) = RNS(IC) + dQx / ( PRI(IC)*GRAV )
endif
IF (DO_TRACERS) THEN
WFN = WFN*0.5 *1.0 !*FRICFAC*0.5
TEM = WFN*PRI(K)
DO ITR=1,ITRCR
XCU(K,ITR) = XCU(K,ITR) + TEM * (XOI(K-1,ITR) - XOI(K,ITR))
END DO
DO ITR=1,ITRCR
DO L=K-1,IC+1,-1
TEM = WFN*PRI(L)
XCU(L,ITR) = XCU(L,ITR) + TEM * &
( (XOI(L-1,ITR) - XOI(L,ITR )) * ETA(L) &
+ (XOI(L,ITR ) - XOI(L+1,ITR)) * ETA(L+1) )
ENDDO
ENDDO
TEM = WFN*PRI(IC)
DO ITR=1,ITRCR
XCU(IC,ITR) = XCU(IC,ITR) + &
(2.*(XHT(ITR) - XOI(IC,ITR)*(ETA(IC)-ETA(IC+1))) &
- (XOI(IC,ITR)+XOI(IC+1,ITR))*ETA(IC+1))*TEM
ENDDO
DO ITR=1,ITRCR
DO L=IC,K
XOI(L,ITR) = XOI(L,ITR) + XCU(L,ITR)
ENDDO
ENDDO
else
WFN = WFN*0.5 *1.0 !*FRICFAC*0.5
endif
LAMBDSV(IC) = 1.000
! CUMULUS FRICTION
IF(FRICFAC <= 0.0) THEN
RC(IC) = 0
RETURN ! NO CUMULUS FRICTION =========>>
ENDIF
WFN = WFN*FRICFAC*EXP( -ALM / FRICLAMBDA )
TEM = WFN*PRI(K)
UCU(K) = UCU(K) + TEM * (UOI(K-1) - UOI(K))
VCU(K) = VCU(K) + TEM * (VOI(K-1) - VOI(K))
DO L=K-1,IC+1,-1
TEM = WFN*PRI(L)
UCU(L) = UCU(L) + TEM * &
( (UOI(L-1) - UOI(L )) * ETA(L ) &
+ (UOI(L ) - UOI(L+1)) * ETA(L+1) )
VCU(L) = VCU(L) + TEM * &
( (VOI(L-1) - VOI(L )) * ETA(L) &
+ (VOI(L ) - VOI(L+1)) * ETA(L+1) )
ENDDO
TEM = WFN*PRI(IC)
UCU(IC) = UCU(IC) + (2.*(UHT - UOI(IC)*(ETA(IC)-ETA(IC+1))) &
- (UOI(IC)+UOI(IC+1))*ETA(IC+1))*TEM
VCU(IC) = VCU(IC) + (2.*(VHT - VOI(IC)*(ETA(IC)-ETA(IC+1))) &
- (VOI(IC)+VOI(IC+1))*ETA(IC+1))*TEM
DISSK0(IC) = ETA(IC)* GRAV * WFNOG * PRI(IC) * 0.5 &
*( (UHT/ETA(IC) - UOI(IC) )**2 &
+ (VHT/ETA(IC) - VOI(IC) )**2 )
DO L=IC,K
UOI(L) = UOI(L) + UCU(L)
VOI(L) = VOI(L) + VCU(L)
ENDDO
RC(IC) = 0
RETURN
END SUBROUTINE CLOUDE
SUBROUTINE ACRITN( PL, PLB, ACR) 1
REAL, INTENT(IN ) :: PL, PLB
REAL, INTENT(OUT) :: ACR
INTEGER IWK
!!REAL, PARAMETER :: FACM=0.5
REAL, PARAMETER :: &
PH(15)=(/150.0, 200.0, 250.0, 300.0, 350.0, 400.0, 450.0, 500.0, &
550.0, 600.0, 650.0, 700.0, 750.0, 800.0, 850.0/)
REAL, PARAMETER :: &
A(15)=(/ 1.6851, 1.1686, 0.7663, 0.5255, 0.4100, 0.3677, &
0.3151, 0.2216, 0.1521, 0.1082, 0.0750, 0.0664, &
0.0553, 0.0445, 0.0633 /) !!*FACM
IWK = PL * 0.02 - 0.999999999
IF (IWK .GT. 1 .AND. IWK .LE. 15) THEN
ACR = A(IWK-1) + (PL-PH(IWK-1))*.02*(A(IWK)-A(IWK-1))
ELSEIF(IWK > 15) THEN
ACR = A(15)
ELSE
ACR = A(1)
ENDIF
ACR = ACRITFAC * ACR * (PLB - PL)
RETURN
END SUBROUTINE ACRITN
SUBROUTINE RNEVP 1
ZET(K+1) = 0
DO L=K,ICMIN,-1
TEM = POI(L)*(PRJ(L+1)-PRJ(L))*CPBG
ZET(L) = ZET(L+1) + TEM
ENDDO
DO L=ICMIN,K
TEM = PRI(L)*GRAV
CNV_PRC3(I,L) = RNS(L)*TEM
ENDDO
!! If hst is smoothed then adjusted precips may be negative
if ( SMOOTH_HST ) then
DO L=ICMIN,K
if ( CNV_PRC3(I,L) < 0. ) then
QOI(L) = QOI(L) + CNV_PRC3(I,L)
POI(L) = POI(L) - CNV_PRC3(I,L) * (ALHL/CP) / PRJ(L+1)
CNV_PRC3(I,L) = 0.
endif
END DO
endif
RETURN
END SUBROUTINE RNEVP
SUBROUTINE HTEST 1
REAL, DIMENSION(K0) :: HOL1
integer :: lminhol
real :: minhol
hol=0. ! HOL initialized here in order not to confuse Valgrind debugger
lminhol = K+1
MINHOL = -999999.
ZET(K+1) = 0
SHT(K+1) = CP*POI(K)*PRJ(K+1)
DO L=K,ICMIN,-1
QOL(L) = AMIN1(QST(L)*RHMAX,QOI(L))
QOL(L) = AMAX1( 0.000, QOL(L) ) ! GUARDIAN/NEG.PREC.
SSL(L) = CP*PRJ(L+1)*POI(L) + GRAV*ZET(L+1)
HOL(L) = SSL(L) + QOL(L)*ALHL
HST(L) = SSL(L) + QST(L)*ALHL
TEM = POI(L)*(PRJ(L+1)-PRJ(L))*CPBG
ZET(L) = ZET(L+1) + TEM
ZOL(L) = ZET(L+1) + (PRJ(L+1)-PRH(L))*POI(L)*CPBG
ENDDO
HOL1=HOL
DO L=K-1,ICMIN+1,-1
HOL1(L) = (0.25*HOL(L+1)+0.50*HOL(L)+0.25*HOL(L-1))
if ( ( MINHOL>=HOL1(L) ) .OR. (MINHOL<0.) ) THEN
MINHOL = HOL1(L)
LMINHOL = L
end if
ENDDO
SIGE_MINHOL = SIGE(LMINHOL)
end subroutine HTEST
subroutine FINDDTLS 1
real :: SIGDT0,sigmax,sigmin
integer :: LL
integer :: THE_SEED(2)
THE_SEED(1)=SEEDRAS(1)*IRAS(I) + SEEDRAS(2)*JRAS(I)
THE_SEED(2)=SEEDRAS(1)*JRAS(I) + SEEDRAS(2)*IRAS(I)
THE_SEED(1)=THE_SEED(1)*SEEDRAS(1)/( SEEDRAS(2) + 10)
THE_SEED(2)=THE_SEED(2)*SEEDRAS(1)/( SEEDRAS(2) + 10)
call random_seed(PUT=THE_SEED)
SIGMAX=SIGE(K)
SIGMIN=SIGE(ICMIN)
if ( RASNCL < 0.0 ) then
N_DTL = K - ICMIN
else
N_DTL = min( int( RASNCL ) , K-ICMIN )
endif
if(allocated(ICL_V)) deallocate(ICL_V)
allocate(ICL_V(N_DTL))
if ( ( RASNCL < 0.0 ) .and. ( RASNCL >=-100.) ) then
do L=1,N_DTL
ICL_V(L) = ICMIN + L - 1
enddo
else if ( RASNCL < -100.0 ) then
do L=1,N_DTL
ICL_V(L) = K - L
enddo
else
do L=1,N_DTL
call random_number ( SIGDT0 )
SIGDT0 = 1.00 - ( SIGDT0**RDTLEXPON )
SIGDT0 = SIGMIN+SIGDT0*(SIGMAX-SIGMIN)
do LL=ICMIN,K
if ( (SIGE(LL+1)>=SIGDT0) .and. (SIGE(LL)<SIGDT0 ) ) ICL_V(L)=LL
enddo
end do
endif
end subroutine FINDDTLS
SUBROUTINE STRAP(FINAL) 4
INTEGER :: FINAL
REAL , DIMENSION(K0) :: WGHT, MASSF
REAL :: WGHT0, PRCBL
! !DESCRIPTION:
! {\tt STRAP} is called: FINAL=0, to compute cloud base layer CBL properties
! given a value K for the index of the upper {\em EDGE} of the CBL; FINAL=1
! to redistribute convective tendencies within CBL
integer :: KK
! LOCAL VARIABLES FOR USE IN CLOUDE
!!IF (.NOT. PRESENT(FINAL)) THEN
IF (FINAL==0) THEN
!!PRJ(ICMIN:K+1) = PKE(I,ICMIN:K+1)
do kk=icmin,k+1
PRJ(kk) = PKE(I,kk)
enddo
poi=0. ! These initialized here in order not to confuse Valgrind debugger
qoi=0. ! Do not believe it actually makes any difference.
uoi=0.
voi=0.
PRS(ICMIN:K0+1) = PLE(I,ICMIN:K0+1)
POI(ICMIN:K) = THO(I,ICMIN:K)
QOI(ICMIN:K) = QHO(I,ICMIN:K)
UOI(ICMIN:K) = UHO(I,ICMIN:K)
VOI(ICMIN:K) = VHO(I,ICMIN:K)
WSP(ICMIN:K) = SQRT( ( UOI(ICMIN:K) - UOI(K) )**2 &
+ ( VOI(ICMIN:K) - VOI(K) )**2 )
QST(ICMIN:K) = QSS(I,ICMIN:K)
DQQ(ICMIN:K) = DQS(I,ICMIN:K)
IF (DO_TRACERS) THEN
DO ITR=1,ITRCR
XOI(ICMIN:K,ITR) = XHO(I,ICMIN:K,ITR)
END DO
END IF
!!! Mass fraction of each layer below cloud base
!!! contributed to aggregate cloudbase layer (CBL)
MASSF(:) = WGT0(:)
!!! RESET PRESSURE at bottom edge of CBL
PRCBL = PRS(K)
do l= K,K0
PRCBL = PRCBL + MASSF(l)*( PRS(l+1)-PRS(l) )
end do
PRS(K+1) = PRCBL
PRJ(K+1) = (PRS(K+1)/1000.)**(MAPL_RGAS/MAPL_CP)
DO L=K,ICMIN,-1
POL(L) = 0.5*(PRS(L)+PRS(L+1))
PRH(L) = (PRS(L+1)*PRJ(L+1)-PRS(L)*PRJ(L)) &
/ (ONEPKAP*(PRS(L+1)-PRS(L)))
PKI(L) = 1.0 / PRH(L)
DPT(L) = PRH(L ) - PRJ(L)
DPB(L) = PRJ(L+1) - PRH(L)
PRI(L) = .01 / (PRS(L+1)-PRS(L))
ENDDO
!!!!! RECALCULATE PROFILE QUAN. IN LOWEST STRAPPED LAYER
if( K<=K0) then
POI(K) = 0.
QOI(K) = 0.
UOI(K) = 0.
VOI(K) = 0.
!! SPECIFY WEIGHTS GIVEN TO EACH LAYER WITHIN SUBCLOUD "SUPERLAYER"
WGHT = 0.
DO L=K,K0
WGHT(L) = MASSF(L) * &
( PLE(I,L+1) - PLE(I,L) ) &
/( PRS(K+1) - PRS(K) )
END DO
DO L=K,K0
POI(K) = POI(K) + WGHT(L)*THO(I,L)
QOI(K) = QOI(K) + WGHT(L)*QHO(I,L)
UOI(K) = UOI(K) + WGHT(L)*UHO(I,L)
VOI(K) = VOI(K) + WGHT(L)*VHO(I,L)
ENDDO
IF (DO_TRACERS) THEN
XOI(K,:)=0.
DO ITR=1,ITRCR
DO L=K,K0
XOI(K,ITR) = XOI(K,ITR) + WGHT(L)*XHO(I,L,ITR)
END DO
END DO
END IF
DQQ(K) = DQSAT( POI(K)*PRH(K) , POL(K), qsat=QST(K) )
endif
!!DPTH_BL = CPBG*POI(K)*( PRJ(K+1)-PRJ(K) )
DPTH_BL = ZCBL(I)
DO L=K,ICMIN,-1
BET(L) = DQQ(L)*PKI(L) !*
GAM(L) = PKI(L)/(1.0+LBCP*DQQ(L)) !*
IF(L<K) THEN
GHT(L+1) = GAM(L)*DPB(L) + GAM(L+1)*DPT(L+1)
GM1(L+1) = 0.5*LBCP*(DQQ(L )/(ALHL*(1.0+LBCP*DQQ(L ))) + &
DQQ(L+1)/(ALHL*(1.0+LBCP*DQQ(L+1))) )
ENDIF
ENDDO
TCU(ICMIN:K) = -POI(ICMIN:K)*PRH(ICMIN:K)
QCU(ICMIN:K) = -QOI(ICMIN:K)
RNS = 0.
CLL = 0.
RMF = 0.
RMFD = 0.
RMFC = 0.
RMFP = 0.
CLL0 = 0.
DLL0 = 0.
CLLX = 0.
DLLX = 0.
CLLI = 0.
CLLB = 0.
POI_SV = POI
QOI_SV = QOI
UOI_SV = UOI
VOI_SV = VOI
IF (DO_TRACERS) XOI_SV = XOI
LAMBDSV = 0.0
CVW = 0.0
UPDFRC = 0.0
UPDFRP = 0.0
DISSK0 = 0.0
END IF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! IF (PRESENT(FINAL)) THEN
IF (FINAL==1) THEN
THO(I,ICMIN:K-1) = POI(ICMIN:K-1)
QHO(I,ICMIN:K-1) = QOI(ICMIN:K-1)
UHO(I,ICMIN:K-1) = UOI(ICMIN:K-1)
VHO(I,ICMIN:K-1) = VOI(ICMIN:K-1)
CNV_UPDFRC(I,ICMIN:K-1) = UPDFRC(ICMIN:K-1)
CNV_QC(I,ICMIN:K-1) = CLLB(ICMIN:K-1)
!! De-strap tendencies from RAS
!! specify weighting "SHAPE"
WGHT = WGT1
!! Scale properly by layer masses
wght0 = 0.
DO L=K,K0
wght0 = wght0 + WGHT(L)* ( PLE(I,L+1) - PLE(I,L) )
END DO
wght0 = ( PRS(K+1) - PRS(K) )/wght0
WGHT = wght0 * WGHT
DO L=K,K0
THO(I,L) = THO(I,L) + WGHT(L)*(POI(K) - POI_SV(K))
QHO(I,L) = QHO(I,L) + WGHT(L)*(QOI(K) - QOI_SV(K))
UHO(I,L) = UHO(I,L) + WGHT(L)*(UOI(K) - UOI_SV(K))
VHO(I,L) = VHO(I,L) + WGHT(L)*(VOI(K) - VOI_SV(K))
END DO
IF (DO_TRACERS) THEN
XHO(I,ICMIN:K-1,:) = XOI(ICMIN:K-1,:)
DO ITR=1,ITRCR
DO L=K,K0
XHO(I,L,ITR) = XHO(I,L,ITR) + WGHT(L)*(XOI(K,ITR) - XOI_SV(K,ITR))
END DO
END DO
END IF
FLX (I,ICMIN:K) = RMF (ICMIN:K) * DDT/DAYLEN ! (KG/m^2/s @ CLOUD BASE)
FLXD(I,ICMIN:K) = RMFD(ICMIN:K) * DDT/DAYLEN ! (KG/m^2/s @ CLOUD TOP)
FLXC(I,ICMIN:K) = RMFC(ICMIN:K) * DDT/DAYLEN ! (KG/m^2/s @ CLOUD TOP)
CLW (I,ICMIN:K) = CLL (ICMIN:K) * DDT/DAYLEN ! (KG/m^2/s )
if ( PRESENT( DISSKE )) DISSKE(I,ICMIN:K-1) = DISSK0(ICMIN:K-1) * DDT/DAYLEN
FLX (I,1:ICMIN-1) = 0.
FLXD(I,1:ICMIN-1) = 0.
FLXC(I,1:ICMIN-1) = 0.
CLW (I,1:ICMIN-1) = 0.
IF ( K < K0 ) THEN
FLX (I,K:K0) = 0.
FLXD(I,K:K0) = 0.
FLXC(I,K:K0) = 0.
CLW (I,K:K0) = 0.
END IF
IF(PRESENT(IRC)) THEN
! IRC (I,1:ICMIN-1) = -2
IRC (I,ICMIN:K-1) = RC(ICMIN:K-1)
! IRC (I,K:K0 ) = -2
ENDIF
IF(ALLOCATED(ICL_V)) DEALLOCATE( ICL_V )
ENDIF
IF (FINAL==2) THEN
FLX (I,:) = 0.
FLXD(I,:) = 0.
FLXC(I,:) = 0.
CLW (I,:) = 0.
IF(PRESENT(IRC)) THEN
! IRC (I,1:ICMIN-1) = -2
IRC (I,ICMIN:K-1) = RC(ICMIN:K-1)
! IRC (I,K:K0 ) = -2
ENDIF
ENDIF
RETURN
END SUBROUTINE STRAP
!*********************************************************************
END SUBROUTINE RASE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE SUNDQ3_ICE( TEMP,RATE2,RATE3,TE1, F2, F3) 1
REAL, INTENT( IN) :: TEMP,RATE2,RATE3,TE1
REAL, INTENT(OUT) :: F2, F3
REAL :: XX, YY,TE0,TE2,JUMP1 !,RATE2,RATE3,TE1
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Ice - phase treatment totally invented
!! Sharp increase in autoconversion in range
!! ~~TE1 K ~< T < TE0 K .
!! (JTB, 3/25/2003)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
TE0=273.
TE2=200.
JUMP1= (RATE2-1.0) / ( ( TE0-TE1 )**0.333 )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Ice - phase treatment !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF ( TEMP .GE. TE0 ) THEN
F2 = 1.0
F3 = 1.0
ENDIF
IF ( ( TEMP .GE. TE1 ) .AND. ( TEMP .LT. TE0 ) )THEN
F2 = 1.0 + JUMP1 * (( TE0 - TEMP )**0.3333)
F3 = 1.0
ENDIF
IF ( TEMP .LT. TE1 ) THEN
F2 = RATE2 + (RATE3-RATE2)*(TE1-TEMP)/(TE1-TE2)
F3 = 1.0
ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF ( F2 .GT. 27.0 ) F2=27.0
end subroutine sundq3_ice
END MODULE RAS