! $Id: cloudnew.F90,v 1.40.2.7.2.2.10.1.4.1.4.8.2.1.2.3.6.10 2008/04/25 19:15:33 ltakacs Exp $
! $Name: GEOSagcm-Eros_7_24 $
module cloudnew 1,2
use GEOS_UtilsMod, only : QSAT=>GEOS_Qsat, DQSAT=>GEOS_DQsat
use MAPL_ConstantsMod
use DDF
implicit none
private
PUBLIC PROGNO_CLOUD
PUBLIC ICE_FRACTION
PUBLIC T_CLOUD_CTL
type T_CLOUD_CTL
real :: SCLMFDFR
real :: RSUB_RADIUS
end type T_CLOUD_CTL
real, parameter :: T_ICE_MAX= MAPL_TICE ! -7.0+TICE
real, parameter :: RHO_W = 1.0e3 ! Density of liquid water in kg/m^3
real, parameter :: MIN_CLD_FRAC =1.0e-8
!! Some parameters set by PHYSPARAMS
real :: T_ICE_ALL ! =-40.+MAPL_TICE ! -40.+TICE
real :: RH00 ! = PHYSPARAMS(11) ! Critical relative humidity
real :: C_00 ! = PHYSPARAMS(12)
real :: LWCRIT ! = PHYSPARAMS(13)
real :: C_ACC ! = PHYSPARAMS(14) ! m^2/kg
real :: C_EV ! = PHYSPARAMS(19) ! 1/s
!!real :: ICEFRPWR
integer :: ICEFRPWR
integer :: NSMAX ! = PHYSPARAMS(46)
real :: MIN_CLD_WATER ! = PHYSPARAMS(18)
real :: CLD_EVP_EFF ! = PHYSPARAMS(28) ! 1/s
real :: PRC_MAX_FALL ! = PHYSPARAMS(34) / DAYLEN
real :: CNV_BETA ! = PHYSPARAMS(1) ! Area factor for convective rain showers (non-dim)
real :: ANV_BETA
real :: LS_BETA
real :: SHR_EVAP_FAC
real :: CLDVOL2FRC, RHSUP_ICE
real :: LS_SDQV2 , LS_SDQV3 , LS_SDQVT1
real :: ANV_SDQV2 , ANV_SDQV3 , ANV_SDQVT1
real :: ANV_TO_LS, ANV_ICEFALL_C, LS_ICEFALL_C
real :: REVAP_OFF_P
REAL :: N_ICE,N_WARM,N_ANVIL,N_PBL
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
contains
subroutine progno_cloud ( & 1,18
!!! first vars are (in) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
DT , &
ZLO , &
ZLE , &
PP , &
PPE , &
EXNP , &
EXNPE , &
FRLAND , &
KH , &
RMFDTR , &
QLWDTR , &
QRN_CU , &
CNV_MFC , &
CNV_UPDFRC , &
U , &
V , &
TH , &
Q , &
QLW_LS , &
QLW_AN , &
QIW_LS , &
QIW_AN , &
ANVFRC , &
CLDFRC , &
RAD_CLDFRC , &
RAD_QL , &
RAD_QI , &
CLDREFFL , &
CLDREFFI , &
CLDNCCN , &
PRELS , &
PRECU , &
PREAN , &
LSARF , &
CUARF , &
ANARF , &
SNRLS , &
SNRCU , &
SNRAN , &
PHYSPARAMS , &
CTLV , &
RHX , &
REV_LS , &
REV_AN , &
REV_CN , &
RSU_LS , &
RSU_AN , &
RSU_CN , &
ACLL_CN,ACIL_CN, &
ACLL_AN,ACIL_AN, &
ACLL_LS,ACIL_LS, &
PFL_CN,PFI_CN, &
PFL_AN,PFI_AN, &
PFL_LS,PFI_LS, &
PDFL,PDFI,FIXL,FIXI, &
AUT, EVAPC, SDM, SUBLC, &
FRZ_TT, DCNVL, DCNVi, &
ALPHT, ALPH1, ALPH2, &
CFPDF, CFPDFX, RHXX, &
DQRL,FRZ_PP )
type (T_CLOUD_CTL) , intent(in) :: CTLV
real , intent(in) :: DT
real , intent(in) , dimension(:,:,:) :: PP, EXNP, ZLO
real , intent(in) , dimension(:,:,0:) :: PPE, EXNPE, KH, CNV_MFC, ZLE
real , intent(in) , dimension(:,: ) :: FRLAND
real , intent(inout) , dimension(:,:,:) :: TH !,TEMP
real , intent(inout) , dimension(:,:,:) :: Q
real , intent(in) , dimension(:,:,:) :: RMFDTR
real , intent(in) , dimension(:,:,:) :: QLWDTR
real , intent(inout) , dimension(:,:,:) :: QRN_CU
real , intent(inout) , dimension(:,:,:) :: CNV_UPDFRC ! on edges, but dims=1:LM
real , intent(in) , dimension(:,:,:) :: U
real , intent(in) , dimension(:,:,:) :: V
real , intent(inout) , dimension(:,:,:) :: QLW_LS
real , intent(inout) , dimension(:,:,:) :: QLW_AN
real , intent(inout) , dimension(:,:,:) :: QIW_LS
real , intent(inout) , dimension(:,:,:) :: QIW_AN
real , intent(inout) , dimension(:,:,:) :: ANVFRC
real , intent(inout) , dimension(:,:,:) :: CLDFRC
real , intent(out) , dimension(:,:,:) :: RAD_CLDFRC
real , intent(out) , dimension(:,:,:) :: RAD_QL
real , intent(out) , dimension(:,:,:) :: RAD_QI
real , intent(out) , dimension(:,:,:) :: CLDREFFL
real , intent(out) , dimension(:,:,:) :: CLDREFFI
real , intent(out) , dimension(:,:,:) :: CLDNCCN
real , intent(out) , dimension(:,: ) :: PRELS
real , intent(out) , dimension(:,: ) :: PRECU
real , intent(out) , dimension(:,: ) :: PREAN
real , intent(out) , dimension(:,: ) :: LSARF
real , intent(out) , dimension(:,: ) :: CUARF
real , intent(out) , dimension(:,: ) :: ANARF
real , intent(out) , dimension(:,: ) :: SNRLS
real , intent(out) , dimension(:,: ) :: SNRCU
real , intent(out) , dimension(:,: ) :: SNRAN
real , intent(in) , dimension(:) :: PHYSPARAMS
real , pointer , dimension(:,:,:) :: REV_CN
real , pointer , dimension(:,:,:) :: REV_AN
real , pointer , dimension(:,:,:) :: REV_LS
real , pointer , dimension(:,:,:) :: RSU_CN
real , pointer , dimension(:,:,:) :: RSU_AN
real , pointer , dimension(:,:,:) :: RSU_LS
real , pointer , dimension(:,:,:) :: ACIL_CN,ACLL_CN
real , pointer , dimension(:,:,:) :: ACIL_AN,ACLL_AN
real , pointer , dimension(:,:,:) :: ACIL_LS,ACLL_LS
real , pointer , dimension(:,:,:) :: PFI_CN,PFL_CN
real , pointer , dimension(:,:,:) :: PFI_AN,PFL_AN
real , pointer , dimension(:,:,:) :: PFI_LS,PFL_LS
real , pointer , dimension(:,:,:) :: PDFL,PDFI,FIXL,FIXI
real , pointer , dimension(:,:,:) :: AUT, EVAPC, SDM, SUBLC
real , pointer , dimension(:,:,:) :: FRZ_TT, DCNVL, DCNVi,FRZ_PP
real , pointer , dimension(:,:,:) :: DQRL, VF_LS, VF_AN
real , pointer , dimension(:,:,:) :: RHX, ALPHT, ALPH1, ALPH2, CFPDF
real , pointer , dimension(:,:,:) :: RHXX, CFPDFX
integer :: IM , JM , LM , UNIT , I , J , K , ITEST, L
integer :: DISABLE_RAD, FRACTION_REMOVAL
real , dimension( 1:size(Q,1) , 1:size(Q,2) , 1:size(Q,3) ) :: MASS, iMASS
real , dimension( 1:size(Q,1) , 1:size(Q,2) , 1:size(Q,3) ) :: TOTAL_WATER
real , dimension( 1:size(Q,1) , 1:size(Q,2) , 1:size(Q,3) ) :: DQST3,QST3
real , dimension( 1:size(Q,1) , 1:size(Q,2) , 1:size(Q,3) ) :: TOTFRC,OTHFRC ! these are needed, but should not be
real , dimension( 1:size(Q,1) , 1:size(Q,2) , 1:size(Q,3) ) :: QRN_LS, QRN_AN ,QV,TE
real , dimension( 1:size(Q,1) , 1:size(Q,2) , 1:size(Q,3) ) :: QSN_LS, QSN_AN, QSN_CU
real , dimension( 1:size(Q,1) , 1:size(Q,2) , 1:size(Q,3) ) :: QTMP1,QTMP2,QTMP3
real , dimension( 1:size(Q,1) , 1:size(Q,2) ) :: TOT_WATER_VMI_0
real , dimension( 1:size(Q,1) , 1:size(Q,2) ) :: TOT_WATER_VMI_2
real , dimension( 1:size(Q,1) , 1:size(Q,2) ) :: TOT_WATER_VMI_1
real , dimension( 1:size(Q,1) , 1:size(Q,2) ) :: TOT_WATER_VMI_3
real , dimension( 1:size(Q,1) , 1:size(Q,2) ) :: TOT_WATER_VMI_4
real , dimension( 1:size(Q,1) , 1:size(Q,2) ) :: TOT_WATER_VMI_5
real , dimension( 1:size(Q,1) , 1:size(Q,2) ) :: TOT_PREC_VMI
real , dimension( 1:size(Q,1) , 1:size(Q,2) ) :: TOT_ENERG_VMI_0
real , dimension( 1:size(Q,1) , 1:size(Q,2) ) :: TOT_ENERG_VMI_1
real , dimension( 1:size(Q,1) , 1:size(Q,2) ) :: TOT_ENERG_VMI_2
real , dimension( 1:size(Q,1) , 1:size(Q,2) ) :: TOT_ENERG_VMI_3,TOT_ENERG_VMI_4,TOT_ENERG_VMI_5
real , dimension( 1:size(Q,1) , 1:size(Q,2) ) :: TOT_ENERG_PREC
real , dimension( 1:size(Q,1) , 1:size(Q,2) ) :: AREA_UPD_PRC
real , dimension( 1:size(Q,1) , 1:size(Q,2) ) :: AREA_ANV_PRC
real , dimension( 1:size(Q,1) , 1:size(Q,2) ) :: AREA_LS_PRC
real , dimension( 1:size(Q,1) , 1:size(Q,2) ) :: SHEAR_ENHANCE
real , dimension( 1:size(Q,1) , 1:size(Q,2) , 1:size(Q,3) ) :: TEMP
real , dimension( 1:size(Q,1) , 1:size(Q,2) , 1:size(Q,3) ) :: THETV,DZET,SHEAR !! EXNP
real , dimension( 1:size(Q,1) , 1:size(Q,2) , 1:size(Q,3) ) :: ICEFRCT
real , dimension( 1:size(Q,1) , 1:size(Q,2) , 1:size(Q,3) ) :: RHCRIT, ALPHA
real , dimension( 1:size(Q,1) , 1:size(Q,2) , 1:size(Q,3) ) :: AA3,BB3,VFALL
real , dimension( 1:size(Q,1) , 1:size(Q,2) ) :: TOT_PREC_EVAP
real :: LW ,QO ,TEMPR ,Q_SAT ,ACF ,PREC ,PRESS_0 ,TAU, AREA, SHOWER_INTEN0,NN_ANVIL,NN_ICE,NN_WARM,CNVENVFC
real :: WRHODEP,CNVICEPARAM,CNVDDRFC,ANVDDRFC,LSDDRFC
CNV_BETA = PHYSPARAMS(1) ! 1 Area factor for convective rain showers (non-dim)
ANV_BETA = PHYSPARAMS(2) ! 2 Area factor for anvil rain showers (non-dim)
LS_BETA = PHYSPARAMS(3) ! 3 Area factor for Large Scale rain showers (non-dim)
RH00 = PHYSPARAMS(4) ! 4 ! Critical relative humidity
C_00 = PHYSPARAMS(5) ! 5
LWCRIT = PHYSPARAMS(6) ! 6
C_ACC = PHYSPARAMS(7) ! 7 ! m^2/kg
C_EV = PHYSPARAMS(8) ! 8 ! 1/s
CLDVOL2FRC = PHYSPARAMS(9) ! 9
RHSUP_ICE = PHYSPARAMS(10) ! 10
SHR_EVAP_FAC = PHYSPARAMS(11) ! 11
MIN_CLD_WATER = PHYSPARAMS(12) ! 12
CLD_EVP_EFF = PHYSPARAMS(13) ! 13
NSMAX = INT( PHYSPARAMS(14) ) ! 14
LS_SDQV2 = PHYSPARAMS(15) ! 15
LS_SDQV3 = PHYSPARAMS(16) ! 16
LS_SDQVT1 = PHYSPARAMS(17) ! 17
ANV_SDQV2 = PHYSPARAMS(18) ! 18
ANV_SDQV3 = PHYSPARAMS(19) ! 19
ANV_SDQVT1 = PHYSPARAMS(20) ! 20
ANV_TO_LS = PHYSPARAMS(21) ! 21
N_WARM = PHYSPARAMS(22)
N_ICE = PHYSPARAMS(23)
N_ANVIL = PHYSPARAMS(24)
N_PBL = PHYSPARAMS(25)
DISABLE_RAD = INT( PHYSPARAMS(26) )
ANV_ICEFALL_C = PHYSPARAMS(28)
LS_ICEFALL_C = PHYSPARAMS(29)
REVAP_OFF_P = PHYSPARAMS(30)
CNVENVFC = PHYSPARAMS(31)
WRHODEP = PHYSPARAMS(32)
T_ICE_ALL = PHYSPARAMS(33) + MAPL_TICE
CNVICEPARAM = PHYSPARAMS(34)
ICEFRPWR = INT( PHYSPARAMS(35) + .001 )
CNVDDRFC = PHYSPARAMS(36)
ANVDDRFC = PHYSPARAMS(37)
LSDDRFC = PHYSPARAMS(38)
IM = size( Q , 1)
JM = size( Q , 2)
LM = size( Q , 3)
!Zero out/initialize precips, except QRN_CU which comes from RAS
QRN_LS =0.
QRN_AN =0.
QSN_LS =0.
QSN_AN =0.
QSN_CU =0.
MASS(:,:,1:LM) = ( PPE(:,:,1:LM) - PPE(:,:,0:LM-1) )*100./MAPL_GRAV ! layer-mass (kg/m**2)
iMASS = 1.0 / MASS
TEMP = EXNP * TH
DZET(:,:,1:LM) = TH(:,:,1:LM)*( EXNPE(:,:,1:LM) - EXNPE(:,:,0:LM-1) )*MAPL_CP/MAPL_GRAV
DQST3 = DQSAT ( &
TEMP , &
PP , QSAT = QST3 )
! Energy defined relative to liquid phase
TOT_ENERG_VMI_0 = SUM( ( MAPL_CP*TEMP &
+ MAPL_ALHL*Q &
- MAPL_ALHF*(QIW_LS+QIW_AN) )*MASS, 3 )
! Intitial total water (all phases, temporary reservoirs, etc.)
!! TOT_WATER_VMI_0 = SUM( ( Q + QLW_LS + QLW_AN + QIW_LS + QIW_AN + QRN_CU)*MASS + QLWDTR*DT , 3 )
IF( ASSOCIATED(FRZ_PP) ) FRZ_PP = 0.00
!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Total Condensate Source
!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!
IF( ASSOCIATED(FIXL) ) FIXL = QLW_AN + QLW_LS
IF( ASSOCIATED(FIXI) ) FIXI = QIW_AN + QIW_LS
CALL fix_up_clouds( &
IM,JM,LM , &
Q , &
TEMP , &
QLW_LS , &
QIW_LS , &
CLDFRC , &
QLW_AN , &
QIW_AN , &
ANVFRC )
IF( ASSOCIATED(FIXL) ) FIXL = -( QLW_AN + QLW_LS - FIXL ) / DT
IF( ASSOCIATED(FIXI) ) FIXI = -( QIW_AN + QIW_LS - FIXI ) / DT
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF( ASSOCIATED(FRZ_TT) ) FRZ_TT = QIW_AN + QIW_LS
CALL meltfrz ( &
IM,JM,LM , &
DT , &
TEMP , &
QLW_LS , &
QIW_LS )
CALL meltfrz ( &
IM,JM,LM , &
DT , &
TEMP , &
QLW_AN , &
QIW_AN )
IF( ASSOCIATED(FRZ_TT) ) FRZ_TT = ( QIW_AN + QIW_LS - FRZ_TT ) / DT
!! TOT_WATER_VMI_1 = SUM( ( Q + QLW_LS + QLW_AN + QIW_LS + QIW_AN + QRN_CU)*MASS + QLWDTR*DT , 3 )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF( ASSOCIATED(DCNVi) ) DCNVi = QIW_AN
IF( ASSOCIATED(DCNVL) ) DCNVL = QLW_AN
CALL cnvsrc ( &
IM,JM,LM , &
DT , &
CNVICEPARAM , &
CTLV%SCLMFDFR , &
MASS , &
iMASS , &
PP , &
TEMP , &
Q , &
QLWDTR , &
RMFDTR , &
QLW_AN , &
QIW_AN , &
CLDFRC , &
ANVFRC , &
QST3 )
IF( ASSOCIATED(DCNVi) ) DCNVi = ( QIW_AN - DCNVi ) / DT
IF( ASSOCIATED(DCNVL) ) DCNVL = ( QLW_AN - DCNVL ) / DT
!! TOT_WATER_VMI_1 = SUM( ( Q + QLW_LS + QLW_AN + QIW_LS + QIW_AN + QRN_CU)*MASS , 3 )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF( ASSOCIATED(PDFL) ) PDFL = QLW_LS+QLW_AN
IF( ASSOCIATED(PDFI) ) PDFI = QIW_LS+QIW_AN
call pdf_spread (U,V,KH,DZET,CNV_UPDFRC,ALPHA,ALPHT,ALPH1,ALPH2 )
! impose a minimum amount of variability
ALPHA = MAX( ALPHA , 1.0 - RH00 )
RHCRIT = 1.0 - ALPHA
call hystpdf( &
IM,JM,LM , &
DT,ALPHA , &
PP , &
Q , &
QLW_LS , &
QLW_AN , &
QIW_LS , &
QIW_AN , &
TEMP , &
CLDFRC , &
ANVFRC , &
RHXX , &
CFPDFX )
IF( ASSOCIATED(RHX) ) RHX = Q/QSAT( TEMP, PP )
IF( ASSOCIATED(CFPDF) ) CFPDF = CLDFRC
IF( ASSOCIATED(PDFL) ) PDFL = ( QLW_LS + QLW_AN - PDFL ) / DT
IF( ASSOCIATED(PDFI) ) PDFI = ( QIW_LS + QIW_AN - PDFI ) / DT
!! TOT_WATER_VMI_1 = SUM( ( Q + QLW_LS + QLW_AN + QIW_LS + QIW_AN + QRN_CU)*MASS , 3 )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
TOTFRC = CLDFRC + ANVFRC
WHERE ( TOTFRC > 1.00 )
CLDFRC = CLDFRC*(1.00 / TOTFRC )
ANVFRC = ANVFRC*(1.00 / TOTFRC )
ENDWHERE
TOTFRC = CLDFRC + ANVFRC
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Total Condensate Source
!!!!!!!!!!!!!!!!!!!!!!!!!!!
QTMP1 = QLW_LS + QLW_AN
!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! CONDENSATE/FRACTION SOURCES FINISHED. NOW LOSE CLOUD CONDENSATE !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Total Condensate Sink
!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! E V A P O R A T I O N
!! A N D
!! S U B L I M A T I O N
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF( ASSOCIATED(EVAPC) ) EVAPC = QLW_LS+QLW_AN
IF( ASSOCIATED(SUBLC) ) SUBLC = QIW_LS+QIW_AN
call evap3( &
IM,JM,LM , &
DT , &
RHCRIT , &
PP , &
TEMP , &
Q , &
QLW_AN , &
QIW_AN , &
ANVFRC , &
CLDFRC , &
QST3 )
call subl3( &
IM,JM,LM , &
DT , &
RHCRIT , &
PP , &
TEMP , &
Q , &
QLW_AN , &
QIW_AN , &
ANVFRC , &
CLDFRC , &
QST3 )
IF( ASSOCIATED(EVAPC) ) EVAPC = ( EVAPC - (QLW_LS+QLW_AN) ) / DT
IF( ASSOCIATED(SUBLC) ) SUBLC = ( SUBLC - (QIW_LS+QIW_AN) ) / DT
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
QTMP1 = QLW_LS + QLW_AN
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! TOT_WATER_VMI_1 = SUM( ( Q + QLW_LS + QLW_AN + QIW_LS + QIW_AN + QRN_CU)*MASS , 3 )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! A U T O C O N V E R S I O N
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! FRACTION_REMOVAL = 0 -> none
! 1 -> constant in-cloud QC
! 2 -> trim high edge of PDF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF( ASSOCIATED(AUT) ) AUT = QLW_LS+QLW_AN
FRACTION_REMOVAL = 1
call autocon3( &
IM,JM,LM , &
DT , &
QLW_LS , &
QRN_LS , &
TEMP , &
PP , &
KH , &
CLDFRC , &
LS_SDQV2 , &
LS_SDQV3 , &
LS_SDQVT1 , &
FRACTION_REMOVAL )
FRACTION_REMOVAL = 0
call autocon3( &
IM,JM,LM , &
DT , &
QLW_AN , &
QRN_AN , &
TEMP , &
PP , &
KH , &
ANVFRC , &
ANV_SDQV2 , &
ANV_SDQV3 , &
ANV_SDQVT1 , &
FRACTION_REMOVAL )
IF( ASSOCIATED(AUT) ) AUT = ( AUT - ( QLW_AN + QLW_LS ) )/DT
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Ice cloud settling
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF( ASSOCIATED(SDM) ) SDM = QIW_AN+QIW_LS
FRACTION_REMOVAL = 0
CALL SETTLE_VEL( &
IM, JM, LM , &
WRHODEP , &
QIW_AN , &
PP , &
TEMP , &
ANVFRC , &
KH , &
VFALL , &
ANVIL = 1 )
VFALL = ANV_ICEFALL_C * VFALL
CALL ICEFALL( &
QIW_AN , &
DZET , &
QSN_AN , &
VFALL , &
ANVFRC , &
DT , &
FRACTION_REMOVAL )
FRACTION_REMOVAL = 1
CALL SETTLE_VEL( &
IM, JM, LM , &
WRHODEP , &
QIW_LS , &
PP , &
TEMP , &
CLDFRC , &
KH , &
VFALL , &
LARGESCALE = 1 )
VFALL = LS_ICEFALL_C * VFALL
CALL ICEFALL( &
QIW_LS , &
DZET , &
QSN_LS , &
VFALL , &
CLDFRC , &
DT , &
FRACTION_REMOVAL )
IF( ASSOCIATED(SDM) ) SDM = ( SDM - (QIW_LS + QIW_AN) )/DT
IF( ASSOCIATED(DQRL) ) DQRL = ( QRN_LS + QRN_AN + QSN_LS + QSN_AN ) / DT
!! TOT_WATER_VMI_2 = SUM( ( Q + QLW_LS + QLW_AN + QIW_LS + QIW_AN + QRN_CU &
!! + QRN_LS + QRN_AN + QSN_LS + QSN_AN)*MASS , 3 )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Add in convective rain
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! CU-FREEZE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Also "freeze" out any conv. precip that needs
! to be since this isn't done in RAS. This is
! precip w/ large particles, so freezing is
! strict. Check up on this!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
QTMP2 = 0.
where( TEMP < MAPL_TICE )
QTMP2 = QRN_CU
QSN_CU = QRN_CU
QRN_CU = 0.
TEMP = TEMP + QSN_CU*(MAPL_ALHS-MAPL_ALHL) / MAPL_CP
endwhere
IF( ASSOCIATED(FRZ_PP) ) FRZ_PP = FRZ_PP + QTMP2/DT
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Water after all local conversions
!! TOT_WATER_VMI_2 = SUM( ( Q + QLW_LS + QLW_AN + QIW_LS + QIW_AN + QRN_CU + QSN_CU &
!! + QRN_LS + QRN_AN + QSN_LS + QSN_AN)*MASS , 3 )
!----------------------------------------------------------------------------------------------
! Column will now be swept from top-down for precip accumulation/accretion/re-evaporation
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! ZERO DIAGNOSTIC OUTPUTS BEFORE SHOWERS !!
PRELS(:,:) = 0.
PRECU(:,:) = 0.
PREAN(:,:) = 0.
SNRCU(:,:) = 0.
SNRLS(:,:) = 0.
SNRAN(:,:) = 0.
!!FRAC_SHOWERS(:) = 1./NSHOWERS
TOT_PREC_VMI = SUM( ( QRN_CU + QSN_CU ) * MASS , 3 )
AREA_UPD_PRC = SUM( CNV_UPDFRC* ( QRN_CU + QSN_CU )* MASS , 3 )
WHERE ( TOT_PREC_VMI > 0.0 )
AREA_UPD_PRC = MAX( AREA_UPD_PRC/TOT_PREC_VMI, 1.E-6 )
ENDWHERE
TOT_PREC_VMI = SUM( ( QRN_AN + QSN_AN ) * MASS , 3 )
AREA_ANV_PRC = SUM( ANVFRC*( QRN_AN + QSN_AN ) * MASS , 3 )
WHERE ( TOT_PREC_VMI > 0.0 )
AREA_ANV_PRC = MAX( AREA_ANV_PRC/TOT_PREC_VMI, 1.E-6 )
ENDWHERE
TOT_PREC_VMI = SUM( ( QRN_LS + QSN_LS ) * MASS , 3 )
AREA_LS_PRC = SUM( CLDFRC* ( QRN_LS + QSN_LS ) * MASS , 3 )
WHERE ( TOT_PREC_VMI > 0.0 )
AREA_LS_PRC = MAX( AREA_LS_PRC/TOT_PREC_VMI , 1.E-6 )
ENDWHERE
AREA_LS_PRC = LS_BETA * AREA_LS_PRC
AREA_UPD_PRC= CNV_BETA * AREA_UPD_PRC
AREA_ANV_PRC= ANV_BETA * AREA_ANV_PRC
!! "couple" to diagnostic areal fraction output
!! Intensity factor in PRECIP3 is floored at
!! 1.0. So this is fair.
LSARF = MIN( AREA_LS_PRC, 1.0 )
CUARF = MIN( AREA_UPD_PRC, 1.0 )
ANARF = MIN( AREA_ANV_PRC, 1.0 )
! First a downdraft involving anvil precip
!call DNDRFTDRVR ( DT, TEMP, Q, QLW_AN, QIW_AN, QRN_AN, QSN_AN , &
! ZLO, ZLE, PP , PPE, EXNP, CNV_MFC )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! GET SOME MICROPHYSICAL QUANTITIES
CALL MICRO_AA_BB_3( TEMP,PP,QST3,AA3,BB3 )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
QTMP1 = QLW_LS + QLW_AN
QTMP2 = QIW_LS + QIW_AN
call PRECIP3 ( &
IM,JM,LM , &
DT , &
FRLAND , &
RHCRIT , &
QRN_CU , &
QSN_CU , &
QTMP1 , &
QTMP2 , &
TEMP , &
Q , &
mass , &
imass , &
PPE , &
PP , &
DZET , &
AA3 , &
BB3 , &
AREA_UPD_PRC , &
PRECU , &
SNRCU , &
REV_CN , &
RSU_CN , &
ACLL_CN , &
ACIL_CN , &
PFL_CN , &
PFI_CN , &
FRZ_PP , &
ENVFC = CNVENVFC , DDRFC = CNVDDRFC )
call PRECIP3 ( &
IM,JM,LM , &
DT , &
FRLAND , &
RHCRIT , &
QRN_AN , &
QSN_AN , &
QTMP1 , &
QTMP2 , &
TEMP , &
Q , &
mass , &
imass , &
PPE , &
PP , &
DZET , &
AA3 , &
BB3 , &
AREA_ANV_PRC , &
PREAN , &
SNRAN , &
REV_AN , &
RSU_AN , &
ACLL_AN , &
ACIL_AN , &
PFL_AN , &
PFI_AN , &
FRZ_PP , &
DDRFC = ANVDDRFC )
call PRECIP3 ( &
IM,JM,LM , &
DT , &
FRLAND , &
RHCRIT , &
QRN_LS , &
QSN_LS , &
QTMP1 , &
QTMP2 , &
TEMP , &
Q , &
mass , &
imass , &
PPE , &
PP , &
DZET , &
AA3 , &
BB3 , &
AREA_LS_PRC , &
PRELS , &
SNRLS , &
REV_LS , &
RSU_LS , &
ACLL_LS , &
ACIL_LS , &
PFL_LS , &
PFI_LS , &
FRZ_PP , &
DDRFC = LSDDRFC )
WHERE ( (QLW_LS+QLW_AN) > 0.00 )
QTMP3 = 1./(QLW_LS+QLW_AN)
ELSEWHERE
QTMP3 = 0.0
ENDWHERE
QLW_LS = QLW_LS * QTMP1 * QTMP3
QLW_AN = QLW_AN * QTMP1 * QTMP3
WHERE ( (QIW_LS+QIW_AN) > 0.00 )
QTMP3 = 1./(QIW_LS+QIW_AN)
ELSEWHERE
QTMP3 = 0.0
ENDWHERE
QIW_LS = QIW_LS * QTMP2 * QTMP3
QIW_AN = QIW_AN * QTMP2 * QTMP3
!! TOT_WATER_VMI_3 = SUM( ( Q + QLW_LS + QLW_AN + QIW_LS + QIW_AN)*MASS , 3 ) &
!! + (PRECU + SNRCU + PRELS + SNRLS + PREAN + SNRAN )*DT
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
TH = TEMP / EXNP
NN_ANVIL = N_ANVIL * 1.e6
NN_ICE = N_ICE * 1.e6
NN_WARM = N_WARM * 1.e6
!!#if 0
call RADCOUPLE ( &
IM, JM, LM, &
TEMP, &
PP, &
CLDFRC, &
ANVFRC, &
QLW_LS, &
QIW_LS, &
QLW_AN, &
QIW_AN, &
RAD_QL, &
RAD_QI, &
RAD_CLDFRC, &
CLDREFFL, &
CLDREFFI, &
CLDVOL2FRC, &
NN_ANVIL,NN_ICE,NN_WARM, &
DISABLE_RAD )
END SUBROUTINE PROGNO_CLOUD
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! P R O C E S S S U B R O U T I N E S !!
!! * * * * * !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! P R O C E S S S U B R O U T I N E S !!
!! * * * * * !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!! !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! P R O C E S S S U B R O U T I N E S !!
!! * * * * * !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine pdf_spread (U,V,KH,DZ,UPDF,ALPHA, ALPHT, ALPH1, ALPH2 )
real , intent(in), dimension(:,:,:) :: U,V,DZ,UPDF
real , intent(in), dimension(:,:,0:) :: KH
real , intent(out), dimension(:,:,:) :: ALPHA
real , pointer, dimension(:,:,:) :: ALPH1, ALPH2, ALPHT
integer :: i,j,l,im,jm,lm
real, dimension( size(u,1) , size(u,2) , size(u,3) ) :: A1,A2,A3
! why SHOULDNT RH-crit = 100% ????
! alpha is the 1/2*width so RH_crit=1.0-alpha
IM = size(u,1)
JM = size(u,2)
LM = size(u,3)
alpha = 0.001 ! 0.1% RH SLOP
!! DIRECTIONAL SHEAR == ABS( e_normal dot [U_z,V_z] )
A1=0.
A3 = 1./SQRT( U**2 + V**2 + 0.01 ) ! inverse of wind mag
A2 = V * A3 ! x-component of unit normal to (U,V)
A1(:,:,2:lm-1) = ( A2(:,:,2:lm-1) * ( U(:,:,1:LM-2) - U(:,:,3:LM) ) &
/ ( DZ(:,:,1:LM-2)+DZ(:,:,3:LM) ) )
A2 = -U * A3 ! y-component of unit normal to (U,V)
A1(:,:,2:lm-1) = ( A2(:,:,2:lm-1) * ( V(:,:,1:LM-2) - V(:,:,3:LM) ) &
/ ( DZ(:,:,1:LM-2)+DZ(:,:,3:LM) ) ) + A1(:,:,2:lm-1)
A1 = ABS( A1 ) ! A1 is now magnitude of veering shear at layers in (m/s) /m. Thus, A1=.001 ==> 1 m/s/km
ALPHA = ALPHA + 10.*A1
if(associated(ALPH1)) ALPH1 = 10.*A1
!! Total shear = SQRT( [U_z,V_z] dot [U_z,V_z] )
A1 = 0.
A1(:,:,2:lm-1) = ( ( U(:,:,1:LM-2) - U(:,:,3:LM) )/ ( DZ(:,:,1:LM-2)+DZ(:,:,3:LM) ) )**2 &
+ ( ( V(:,:,1:LM-2) - V(:,:,3:LM) )/ ( DZ(:,:,1:LM-2)+DZ(:,:,3:LM) ) )**2
A1 = SQRT ( A1 ) ! A1 is now magnitude of TOTAL shear at layers in (m/s) /m. Thus, A1=.001 ==> 1 m/s/km
ALPHA = ALPHA + 3.33*A1
!! KH values ~100 m+2 s-1 typical of strong PBLs
ALPHA = ALPHA + 0.002*KH(:,:,0:LM-1)
if(associated(ALPH2)) ALPH2 = 0.002*KH(:,:,0:LM-1)
if(associated(ALPHT)) ALPHT = ALPHA
ALPHA = MIN( ALPHA , 0.25 ) ! restrict RHcrit to > 75%
end subroutine pdf_spread
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine fix_up_clouds ( & 1
IM,JM,LM , &
QV , &
TE , &
QLC , &
QIC , &
CF , &
QLA , &
QIA , &
AF )
integer, intent(in) :: IM,JM,LM
real , intent(inout), dimension(IM,JM,LM) :: TE,QV,QLC,CF,QLA,AF,QIC,QIA
! Fix if Anvil cloud fraction too small
WHERE ( AF < 1.E-5 )
QV = QV + QLA + QIA
TE = TE - (MAPL_ALHL/MAPL_CP)*QLA - (MAPL_ALHS/MAPL_CP)*QIA
AF = 0.
QLA = 0.
QIA = 0.
ENDWHERE
! Fix if LS cloud fraction too small
WHERE ( CF < 1.E-5 )
QV = QV + QLC + QIC
TE = TE - (MAPL_ALHL/MAPL_CP)*QLC - (MAPL_ALHS/MAPL_CP)*QIC
CF = 0.
QLC = 0.
QIC = 0.
ENDWHERE
! LS LIQUID too small
WHERE ( QLC < 1.E-8 )
QV = QV + QLC
TE = TE - (MAPL_ALHL/MAPL_CP)*QLC
QLC = 0.
ENDWHERE
! LS ICE too small
WHERE ( QIC < 1.E-8 )
QV = QV + QIC
TE = TE - (MAPL_ALHS/MAPL_CP)*QIC
QIC = 0.
ENDWHERE
! Anvil LIQUID too small
WHERE ( QLA < 1.E-8 )
QV = QV + QLA
TE = TE - (MAPL_ALHL/MAPL_CP)*QLA
QLA = 0.
ENDWHERE
! Anvil ICE too small
WHERE ( QIA < 1.E-8 )
QV = QV + QIA
TE = TE - (MAPL_ALHS/MAPL_CP)*QIA
QIA = 0.
ENDWHERE
! Fix ALL cloud quants if Anvil cloud LIQUID+ICE too small
WHERE ( ( QLA + QIA ) < 1.E-8 )
QV = QV + QLA + QIA
TE = TE - (MAPL_ALHL/MAPL_CP)*QLA - (MAPL_ALHS/MAPL_CP)*QIA
AF = 0.
QLA = 0.
QIA = 0.
ENDWHERE
! Ditto if LS cloud LIQUID+ICE too small
WHERE ( ( QLC + QIC ) < 1.E-8 )
QV = QV + QLC + QIC
TE = TE - (MAPL_ALHL/MAPL_CP)*QLC - (MAPL_ALHS/MAPL_CP)*QIC
CF = 0.
QLC = 0.
QIC = 0.
ENDWHERE
end subroutine fix_up_clouds
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine meltfrz ( & 2
IM,JM,LM , &
DT , &
TE , &
QL , &
QI )
real , intent(in) :: DT
integer, intent(in) :: IM,JM,LM
real , intent(inout), dimension(:,:,:) :: TE,QL,QI
real , dimension(im,jm,lm) :: fQi,dQil
real :: taufrz
taufrz = 1000.
fQi = ice_fraction( TE )
dQil = 0.
! freeze liquid
where( TE <= T_ICE_MAX )
dQil = Ql *(1.0 - EXP( -Dt * fQi / taufrz ) )
endwhere
dQil = max( 0., dQil )
Qi = Qi + dQil
Ql = Ql - dQil
TE = TE + (MAPL_ALHS-MAPL_ALHL)*dQil/MAPL_CP
dQil = 0.
! melt ice instantly above 0^C
where( TE > T_ICE_MAX )
dQil = -Qi
endwhere
dQil = max( 0., dQil )
Qi = Qi + dQil
Ql = Ql - dQil
TE = TE + (MAPL_ALHS-MAPL_ALHL)*dQil/MAPL_CP
end subroutine meltfrz
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine hystpdf ( & 1
IM,JM,LM , &
DT,ALPHA , &
PL , &
QV , &
QCl , &
QAl , &
QCi , &
QAi , &
TE , &
CF , &
AF , &
RH_DIAG , &
CF_DIAG )
real , intent(in) :: DT ! ,ALPHA
integer, intent(in) :: IM,JM,LM
real , intent(in) , dimension(IM,JM,LM) :: PL,ALPHA
real , intent(inout), dimension(IM,JM,LM) :: TE,QV,QCl,QCi,CF,QAl,QAi,AF
! optional diagnostics
real , pointer, dimension(:,:,:) :: RH_DIAG,CF_DIAG
! internal arrays
real , dimension(IM,JM,LM) :: QCO, QVO, CFO, QAO, TAU
real , dimension(IM,JM,LM) :: QT, QMX, QMN, DQ, QVtop
real , dimension(IM,JM,LM) :: TEO,QSx,DQsx,QS,DQs
real , dimension(IM,JM,LM) :: QCx, QVx, CFx, QAx, QC, QA, fQi, fQi_A
real , dimension(IM,JM,LM) :: dQAi, dQAl, dQCi, dQCl
real , dimension(IM,JM,LM) :: tmpARR
! internal scalars
integer :: N
QC = QCl + QCi
QA = QAl + QAi
where( QA > 0.0 )
fQi_A = QAi / QA
elsewhere
fQi_A = 0.0
endwhere
TEo = TE
DQSx = DQSAT( &
TEo , &
PL , QSAT=QSx )
where( AF < 1.0 )
tmpARR = 1./(1.-AF)
elsewhere
tmpARR = 0.0
endwhere
CFx = CF*tmpARR
QCx = QC*tmpARR
QVx = ( QV - QSx*AF )*tmpARR
where( AF >= 1.0 )
QVx = QSx*1.e-4
endwhere
where( AF > 0. )
QAx = QA/AF
elsewhere
QAx = 0.
endwhere
QT = QCx + QVx
fQi = ice_fraction( TE )
do n=1,3
! DQS = DQSAT( &
! TEo , &
! PL ,QSAT = QS )
DQS = DQSx
QS = QSx + DQS*(TEo - TE )
QS = MAX( QS , 1.0e-7 )
DQ = 2.0*ALPHA*QS
DQ = MIN( DQ , QT ) ! Disallow PDF from going into negative QTs
QMX = QT + 0.5*DQ
QMN = QT - 0.5*DQ
WHERE( QMX < QS )
CFo = 0.
QCo = 0.
ENDWHERE
WHERE( (QMN < QS) .AND. (QS <= QMX) )
CFo = (QMX - QS)/DQ
QCo = 0.5*CFo*(QMX -QS)
ENDWHERE
WHERE( QMN >= QS )
CFo = 1.0
QCo = QT - QS
ENDWHERE
TAU = DT
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! These lines represent adjustments
! to anvil condensate due to the
! assumption of a stationary TOTAL
! water PDF subject to a varying
! QSAT value during the iteration
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
where( AF > 0. )
QAo = QAx ! + QSx - QS
elsewhere
QAo = 0.
endwhere
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
CFo = CFx + (DT/TAU)*( CFo - CFx )
QCo = QCx + (DT/TAU)*( QCo - QCx ) / ( 1.0 + DQS*MAPL_ALHL/MAPL_CP )
QVo = QVx - (QCo - QCx)
TEo = TE + (1.0-fQi)*(MAPL_ALHL/MAPL_CP)*( (QCo - QCx)*(1.-AF) + (QAo-QAx)*AF ) &
+ fQi* (MAPL_ALHS/MAPL_CP)*( (QCo - QCx)*(1.-AF) + (QAo-QAx)*AF )
if(associated(RH_DIAG)) RH_DIAG = QVo / QS
if(associated(CF_DIAG)) CF_DIAG = CFo
enddo ! qsat iteration
! Update prognostic variables. Deal with special case of AF=1
! Temporary variables QCo, QAo become updated grid means.
where( AF < 1.0 )
CF = CFo * ( 1.-AF)
QCo = QCo * ( 1.-AF)
QAo = QAo * AF
elsewhere
! Special case AF=1, i.e., box filled with anvil.
! - Note: no guarantee QV_box > QS_box
CF = 0. ! Remove any other cloud
QAo = QA + QC ! Add any LS condensate to anvil type
QCo = 0. ! Remove same from LS
QT = QAo + QV ! Total water
! Now set anvil condensate to any excess of total water
! over QSx (saturation value at top)
QAo = MAX( QT - QSx, 0. )
endwhere
! Now take {\em New} condensate and partition into ice and liquid
! taking care to keep both >=0 separately. New condensate can be
! less than old, so $\Delta$ can be < 0.
QCx = QCo - QC
dQCl = (1.0-fQi)*QCx
dQCi = fQi * QCx
where((QCl+dQCl)<0.)
dQCi = dQCi + (QCl+dQCl)
dQCl = -QCl !== dQCl - (QCl+dQCl)
endwhere
where((QCi+dQCi)<0.)
dQCl = dQCl + (QCi+dQCi)
dQCi = -QCi !== dQCi - (QCi+dQCi)
endwhere
QAx = QAo - QA
dQAl = QAx ! (1.0-fQi)*QAx
dQAi = 0. ! fQi * QAx
where((QAl+dQAl)<0.)
dQAi = dQAi + (QAl+dQAl)
dQAl = -QAl
endwhere
where((QAi+dQAi)<0.)
dQAl = dQAl + (QAi+dQAi)
dQAi = -QAi
endwhere
! Clean-up cloud if fractions are too small
where( AF < 1.e-5 )
dQAi = -QAi
dQAl = -QAl
endwhere
where( CF < 1.e-5 )
dQCi = -QCi
dQCl = -QCl
endwhere
QAi = QAi + dQAi
QAl = QAl + dQAl
QCi = QCi + dQCi
QCl = QCl + dQCl
QV = QV - ( dQAi+dQCi+dQAl+dQCl)
!!TE = TE + (MAPL_ALHS/MAPL_CP)*(dQAi+dQCi) + (MAPL_ALHL/MAPL_CP)*(dQAl+dQCl)
TE = TE + (MAPL_ALHL*( dQAi+dQCi+dQAl+dQCl)+MAPL_ALHF*(dQAi+dQCi))/ MAPL_CP
! We need to take care of situations where QS moves past QA
! during QSAT iteration. This should be only when QA/AF is small
! to begin with. Effect is to make QAo negative. So, we
! "evaporate" offending QA's
!
! We get rid of anvil fraction also, although strictly
! speaking, PDF-wise, we should not do this.
where( QAo <= 0. )
QV = QV + QAi + QAl
TE = TE - (MAPL_ALHS/MAPL_CP)*QAi - (MAPL_ALHL/MAPL_CP)*QAl
QAi = 0.
QAl = 0.
AF = 0.
endwhere
end subroutine hystpdf
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine cnvsrc ( & 1
IM,JM,LM , &
DT , &
ICEPARAM , &
SCLMFDFR , &
MASS , &
iMASS , &
PL , &
TE , &
QV , &
DCF , &
DMF , &
QLA , &
QIA , &
CF , &
AF , &
QS )
!INPUTS:
!
! ICEPARAM: 0-1 controls how strongly new conv condensate is partitioned in ice-liquid
! 1 means partitioning follows ice_fraction(TE). 0 means all new condensate is
! liquid
!
! SCLMFDFR: Scales detraining mass flux to a cloud fraction source - kludge. Thinly justified
! by fuzziness of cloud boundaries and existence of PDF of condensates (for choices
! 0.-1.0) or by subgrid layering (for choices >1.0)
integer , intent(in) :: IM,JM,LM
real , intent(in) :: DT,ICEPARAM,SCLMFDFR
real , dimension(IM,JM,LM), intent(in) :: MASS,iMASS
real , dimension(IM,JM,LM), intent(in) :: DMF,PL
real , dimension(IM,JM,LM), intent(in) :: DCF,CF,QS
real , dimension(IM,JM,LM), intent(inout) :: AF,TE,QV
real , dimension(IM,JM,LM), intent(inout) :: QLA, QIA
real , dimension(IM,JM,LM) :: TEND,QVx,QCA,fQi
integer :: STRATEGY
real :: minrhx
STRATEGY = 1
!Minimum allowed env RH
minrhx = 0.001
!Addition of condensate from RAS
TEND = DCF*iMASS
fQi = 0.0 + ICEPARAM*ice_fraction( TE )
QLA = QLA + (1.0-fQi)* TEND*DT
QIA = QIA + fQi * TEND*DT
! dont forget that conv cond has never frozen !!!!
TE = TE + (MAPL_ALHS-MAPL_ALHL) * fQi * TEND*DT/MAPL_CP
QCA = QLA + QIA
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Tiedtke-style anvil fraction !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
TEND=DMF*iMASS * SCLMFDFR
AF = AF + TEND*DT
AF = MIN( AF , 0.99 )
! where ( (AF+CF) > 1.00 )
! AF=1.00-CF
! endwhere
!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Check for funny (tiny, negative)
! external QV s, resulting from assumed
! QV=QSAT within anvil.
!
! Two strategies to fix
! 1) Simply constrain AF assume condensate
! just gets "packed" in
! 2) Evaporate QCA to bring up QVx leave AF alone
! Should include QSAT iteration, but ...
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! QS = QSAT( &
! TE , &
! PL )
where ( AF < 1.0 )
QVx = ( QV - QS * AF )/(1.-AF)
elsewhere
QVx = QS
endwhere
if (STRATEGY==1) then
where ( (( QVx - minrhx*QS ) < 0.0 ) .and. (AF > 0.) )
AF = (QV - minrhx*QS )/( QS*(1.0-minrhx) )
endwhere
where ( AF < 0. ) ! If still cant make suitable env RH then destroy anvil
AF = 0.
QV = QV + QLA + QIA
TE = TE - (MAPL_ALHL*QLA + MAPL_ALHS*QIA)/MAPL_CP
QLA = 0.
QIA = 0.
endwhere
endif
if (STRATEGY==2) then
where ( (( QVx - minrhx*QS ) < 0.0 ) .and. (AF > 0.) )
QV = QV + (1.-AF)*( minrhx*QS - QVx )
QCA = QCA - (1.-AF)*( minrhx*QS - QVx )
TE = TE - (1.-AF)*( minrhx*QS - QVx )*MAPL_ALHL/MAPL_CP
endwhere
endif
end subroutine cnvsrc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine evap3 ( & 1
IM,JM,LM , &
DT, RHCR , &
PL , &
TE , &
QV , &
QL , &
QI , &
F , &
XF , &
QS ) ! , &
!! NN )
integer, intent(in) :: IM,JM,LM
real , intent(in ) :: DT ! , RHCR
real , dimension(IM,JM,LM) , intent(in ) :: PL, RHCR
real , dimension(IM,JM,LM) , intent(inout) :: TE
real , dimension(IM,JM,LM) , intent(inout) :: QV
real , dimension(IM,JM,LM) , intent(inout) :: QL,QI
real , dimension(IM,JM,LM) , intent(inout) :: F
real , dimension(IM,JM,LM) , intent(in ) :: XF
real , dimension(IM,JM,LM) , intent(in ) :: QS
!real , dimension(IM,JM,LM) , intent(in ) :: NN
real , dimension(IM,JM,LM) :: ES,NN,RADIUS,K1,K2,TEFF,QCm,EVAP,RHx,QC !,QS
real, parameter :: epsilon = MAPL_H2OMW/MAPL_AIRMW
real, parameter :: K_COND = 2.4e-2 ! J m**-1 s**-1 K**-1
real, parameter :: DIFFU = 2.2e-5 ! m**2 s**-1
real :: A_eff
A_EFF = CLD_EVP_EFF
!NN = N_CCN*1.0E6 ! convert to particles per cubic METER from N/(cm**3)
NN = 50.*1.0e6
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! EVAPORATION OF CLOUD WATER. !!
!! !!
!! DelGenio et al (1996, J. Clim., 9, 270-303) !!
!! formulation (Eq.'s 15-17) !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! QS = QSAT( &
! TE , &
! PL )
ES = 100.* PL * QS / ( EPSILON + (1.0-EPSILON)*QS ) ! (100's <-^ convert from mbar to Pa)
RHx = MIN( QV/QS , 1.00 )
K1 = (MAPL_ALHL**2) * RHO_W / ( K_COND*MAPL_RVAP*(TE**2))
K2 = MAPL_RVAP * TE * RHO_W / ( DIFFU * (1000./PL) * ES )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Here DIFFU is given for 1000 mb !!
!! so 1000./PR accounts for inc- !!
!! reased diffusivity at lower !!
!! pressure. !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
where ( ( F > 0.) .and. ( QL > 0. ) )
QCm=QL/F
elsewhere
QCm=0.
endwhere
RADIUS = LDRADIUS3(PL,TE,QCm,NN)
where( (RHx < RHCR ) .and.(RADIUS > 0.0) )
TEFF = (RHCR - RHx) / ((K1+K2)*RADIUS**2) ! / (1.00 - RHx)
elsewhere
TEFF = 0.0 ! -999.
endwhere
EVAP = a_eff*QL*DT*TEFF
EVAP = MIN( EVAP , QL )
QC=QL+QI
where (QC > 0.)
F = F * ( QC - EVAP ) / QC
endwhere
QV = QV + EVAP
QL = QL - EVAP
TE = TE - (MAPL_ALHL/MAPL_CP)*EVAP
end subroutine evap3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine subl3 ( & 1
IM,JM,LM , &
DT, RHCR , &
PL , &
TE , &
QV , &
QL , &
QI , &
F , &
XF , &
QS ) ! , &
!! NN )
integer, intent(in) :: IM,JM,LM
real , intent(in ) :: DT !, RHCR
real , dimension(IM,JM,LM) , intent(in ) :: PL, RHCR
real , dimension(IM,JM,LM) , intent(inout) :: TE
real , dimension(IM,JM,LM) , intent(inout) :: QV
real , dimension(IM,JM,LM) , intent(inout) :: QL,QI
real , dimension(IM,JM,LM) , intent(inout) :: F
real , dimension(IM,JM,LM) , intent(in ) :: XF
real , dimension(IM,JM,LM) , intent(in ) :: QS
!real , dimension(IM,JM,LM) , intent(in ) :: NN
real , dimension(IM,JM,LM) :: ES,NN,RADIUS,K1,K2,TEFF,QCm,SUBL,RHx,QC !, QS
real, parameter :: epsilon = MAPL_H2OMW/MAPL_AIRMW
real, parameter :: K_COND = 2.4e-2 ! J m**-1 s**-1 K**-1
real, parameter :: DIFFU = 2.2e-5 ! m**2 s**-1
real :: A_eff
A_EFF = CLD_EVP_EFF
!NN = N_CCN*1.0E6 ! convert to particles per cubic METER from N/(cm**3)
NN = 5.*1.0e6
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! SUBLORATION OF CLOUD WATER. !!
!! !!
!! DelGenio et al (1996, J. Clim., 9, 270-303) !!
!! formulation (Eq.'s 15-17) !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! QS = QSAT( &
! TE , &
! PL )
ES = 100.* PL * QS / ( EPSILON + (1.0-EPSILON)*QS ) ! (100's <-^ convert from mbar to Pa)
RHx = MIN( QV/QS , 1.00 )
K1 = (MAPL_ALHL**2) * RHO_W / ( K_COND*MAPL_RVAP*(TE**2))
K2 = MAPL_RVAP * TE * RHO_W / ( DIFFU * (1000./PL) * ES )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Here DIFFU is given for 1000 mb !!
!! so 1000./PR accounts for inc- !!
!! reased diffusivity at lower !!
!! pressure. !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
where ( ( F > 0.) .and. ( QI > 0. ) )
QCm=QI/F
elsewhere
QCm=0.
endwhere
RADIUS = LDRADIUS3(PL,TE,QCm,NN)
where( (RHx < RHCR ) .and.(RADIUS > 0.0) )
TEFF = ( RHCR - RHx) / ((K1+K2)*RADIUS**2) ! / (1.00 - RHx)
elsewhere
TEFF = 0.0 ! -999.
endwhere
SUBL = a_eff*QI*DT*TEFF
SUBL = MIN( SUBL , QI )
QC=QL+QI
where (QC > 0.)
F = F * ( QC - SUBL ) / QC
endwhere
QV = QV + SUBL
QI = QI - SUBL
TE = TE - (MAPL_ALHS/MAPL_CP)*SUBL
end subroutine subl3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine autocon3 ( & 2,1
IM,JM,LM, &
DT , &
QC , &
QP , &
TE , &
PL , &
KH , &
F , &
SUNDQV2 , &
SUNDQV3 , &
SUNDQT1 , &
FRACTION_REMOVAL )
integer, intent(in) :: IM,JM,LM,FRACTION_REMOVAL
real , intent(in ) :: DT
real , dimension(IM,JM,LM) , intent(in ) :: TE
real , dimension(IM,JM,LM) , intent(in ) :: PL
real, intent(in ) , dimension(IM,JM,0:LM) :: KH
real , dimension(IM,JM,LM) , intent(inout) :: QC
real , dimension(IM,JM,LM) , intent(inout) :: QP
real , dimension(IM,JM,LM) , intent(inout) :: F
real , intent(in) :: SUNDQV2, SUNDQV3, SUNDQT1
integer :: NSMX
real :: ACF0, ACF, DTX
real , dimension(IM,JM,LM) :: C00x, iQCcrx, F2, F3,RATE,dQP,QCm
integer :: NS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Precip. conversion from Smith (1990, !
! QJRMS, 116, 435, Eq. 2.29). Similar !
! to Del Genio's Eq.(10). !
! !
! Coalesence term needs to be determined !
! through entire column and is done in !
! subroutine "ACCRETE_EVAP_PRECIP" !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
NSMX = NSMAX
DTX = DT/NSMX
CALL SUNDQ3_ICE3( IM, JM, LM, TE, SUNDQV2, SUNDQV3, SUNDQT1, F2, F3 )
C00x = C_00 * F2 * F3
!QCcrx = LWCRIT / ( F2 * F3 )
iQCcrx = F2 * F3 / LWCRIT
where ( ( F > 0.) .and. ( QC > 0. ) )
QCm = QC/F
elsewhere
QCm = 0.
endwhere
RATE = C00x * ( 1.0 - EXP( - ( QCm * iQCcrx )**2 ) )
!! temporary kluge until we can figure a better to make
!! thicker low clouds ( reuse arrays F2 and F3 )
F2 = 1.0
!where( KH(:,:,0:LM-1) > 10. )
! F2 = 0.1
!endwhere
F3 = 1.0
where( PL > 700. )
F3 = MAX( 1.0 - (PL - 700.)/300. , 0.1 )
endwhere
F3 = MAX( F3*F2 , 0.1 )
RATE = F3 * RATE
dQP = QC*( 1.0 - EXP( -RATE * DT ) )
dQP = MAX( dQP , 0.0 ) ! Protects against floating point problems for tiny RATE
QC = QC - dQP
QP = QP + dQP
SELECT CASE( FRACTION_REMOVAL )
CASE( 0 )
! do NOTHING
CASE( 1 )
where( (QC + dQP) > 0. )
F = QC * F / (QC + dQP )
endwhere
CASE( 2 )
where( (QC + dQP) > 0. .and. QC > 0. )
F = F * SQRT( QC / (QC + dQP ) )
endwhere
CASE( 3 )
where( (QC + dQP) > 0. .and. QC > 0. )
F = F * ( QC / (QC + dQP ) )**0.333
endwhere
END SELECT
RETURN
END SUBROUTINE AUTOCON3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine DNDRFTDRVR ( DT, TEMP, QV, QL, QI, QPL, QPI, &
ZLO, ZLE, PLO, PLE, PK, UMFC )
real, intent(in ) :: DT
real, dimension(:,:,:), intent(inout) :: TEMP,QV,QL,QI,QPL,QPI
real, dimension(:,:,:), intent(in ) :: ZLO,PLO,PK
real, dimension(:,:,0:), intent(in ) :: UMFC,ZLE,PLE
type (T_DDF_CTL) :: CTL
integer :: I,J,L,IM,JM,LM
!!real, dimension( 1 : size(TEMP,3) ) :: DDF_DQDTz, DDF_DTDTz
!!real, dimension( 0 : size(TEMP,3) ) :: DDF_MFCz
!!real :: DDF_ZSCALEz
real, dimension(:), pointer :: DDF_DQDTz, DDF_DTDTz, DDF_MUPHz, DDF_TCz,DDF_QCz
real, dimension(:), pointer :: DDF_MFCz, DDF_RH1z, DDF_RH2z, DDF_BYNCz
real :: DDF_ZSCALEz
IM=size(TEMP,1)
JM=size(TEMP,2)
LM=size(TEMP,3)
allocate( DDF_MFCz(0:LM) )
allocate( DDF_MUPHz(1:LM) )
allocate( DDF_DQDTz(1:LM) )
allocate( DDF_DTDTz(1:LM) )
allocate( DDF_BYNCz(1:LM) )
allocate( DDF_RH1z(1:LM) )
allocate( DDF_RH2z(1:LM) )
allocate( DDF_QCz(1:LM) )
allocate( DDF_TCz(1:LM) )
CTL%ALPHA_DDF_UDF = 0.1
CTL%AREAL_FRACTION = 0.1
! CTL%PARTITION%TYPE = "UNIFORM"
! CTL%PARTITION%RATIO = 0.5
CTL%PARTITION%TYPE = "HEIGHT_DEP"
CTL%PARTITION%MAX_RATIO = 0.9
CTL%PARTITION%MIN_RATIO = 0.1
CTL%THETA_IS_THVAR = .FALSE.
! Comment out to avoid conflicts at compilation if ddf changes
!--------------------------------------------------------------
! do I=1,IM
! do J=1,JM
!
! call DDF1( CTL , DT , TEMP(I,J,:) , &
! QV(I,J,:), QL(I,J,:), QI(I,J,:), &
! QPL(I,j,:), QPI(I,j,:), &
! ZLO(I,J,:), ZLE(I,J,:), PLO(I,J,:), &
! PLE(I,J,:), PK(I,J,:), UMFC(I,J,:) , &
! DDF_ZSCALEz, DDF_DQDTz, DDF_DTDTz, &
! DDF_MUPHz, DDF_RH1z, DDF_RH2z, DDF_BYNCz, DDF_TCz, DDF_QCz, DDF_MFCz )
!
!end do
!end do
! clean up
!--------------------
deallocate( DDF_MFCz )
deallocate( DDF_BYNCz )
deallocate( DDF_DQDTz )
deallocate( DDF_DTDTz )
deallocate( DDF_MUPHz )
deallocate( DDF_RH1z )
deallocate( DDF_RH2z )
deallocate( DDF_TCz )
deallocate( DDF_QCz )
end subroutine DNDRFTDRVR
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine PRECIP3 ( & 3,2
IM,JM,LM , &
DT, FRLAND , &
RHCR3 , &
QPl , &
QPi , &
QCl , &
QCi , &
TE , &
QV , &
mass , &
imass , &
PE , &
PL , &
dZE , &
AA , &
BB , &
AREA , &
RAIN , &
SNOW , &
REVAP , &
RSUBL , &
ACRLL , &
ACRIL , &
PFL_DIAG , &
PFI_DIAG , &
FRZ_DIAG , &
ENVFC,DDRFC )
integer, intent(in) :: IM,LM,JM
real, intent(in) :: DT !, RHCR
real, intent(inout), dimension(IM,JM,LM) :: QPl,QPi,QCl,QCi,TE,QV
real, intent(in), dimension(IM,JM,LM) :: mass,imass,dZE,PL,AA,BB,RHCR3
real, intent(in), dimension(IM,JM,0:LM) :: PE
real, intent(inout), dimension(IM,JM) :: RAIN,SNOW
real, intent(in ), dimension(IM,JM) :: AREA,FRLAND
real , pointer , dimension(:,:,:) :: REVAP
real , pointer , dimension(:,:,:) :: RSUBL
real , pointer , dimension(:,:,:) :: ACRLL,ACRIL
real , pointer , dimension(:,:,: ) :: PFL_DIAG, PFI_DIAG, FRZ_DIAG
real, intent(in), optional :: ENVFC,DDRFC
!!real, intent(in), optional, dimension(IM,JM) :: shearf
real, dimension(IM,JM,LM+1) :: ZE
real, dimension(IM,JM,LM) :: PFi,PFl,QS,dQS,QDDF3,ENVFRAC
real, dimension(IM,JM) :: TKo,QKo,QSTKo,DQSTKo,RH_BOX,T_ED,RHCR,QPlKo,QPiKo
real, dimension(IM,JM) :: Ifactor,AREAx,RAINRAT0,SNOWRAT0
real, dimension(IM,JM) :: FALLRN,FALLSN,VEsn,VErn,NRAIN,NSNOW,Efactor
real, dimension(IM,JM) :: TinLAYERrn,DIAMrn,DROPRAD
real, dimension(IM,JM) :: TinLAYERsn,DIAMsn,FLAKRAD
real :: B_SUB
real, dimension(IM,JM) :: EVAP,SUBL,ACCR,MLTFRZ,EVAPx,SUBLx
real, dimension(IM,JM) :: EVAP_DD,SUBL_DD,DDFRACT
real, dimension(IM,JM) :: PSCL_QP,VMIP,LANDSEAF
real, parameter :: TRMV_L = 1.0 ! m/s
real :: TAU_FRZ, TAU_MLT
integer :: K, NS, NSMX, itr,L
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! fraction of precip falling through "environment" vs
! through cloud
IF(PRESENT(ENVFC)) THEN
ENVFRAC = ENVFC
ELSE
ENVFRAC = 1.00
ENDIF
AREAx = AREA
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
WHERE ( AREAx > 0. )
Ifactor = 1./ ( AREAx )
ELSEwhere
Ifactor = 1.00
ENDwhere
Ifactor = MAX( Ifactor, 1.) !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Start at top of precip column:
!
! a) Accrete
! b) Evaporate/Sublimate
! c) Rain/Snow-out to next level down
! d) return to (a)
!
! ....................................................................
!
! Accretion formulated according to Smith (1990, Q.J.R.M.S., 116, 435
! Eq. 2.29)
!
! Evaporation (ibid. Eq. 2.32)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
PFl=0.
PFi=0.
!!! INITIALIZE DIAGNOSTIC ARRAYS !!!!!!!!!!!!!!!!!!!!!
IF (ASSOCIATED(PFL_DIAG)) PFL_DIAG = 0.
IF (ASSOCIATED(PFI_DIAG)) PFI_DIAG = 0.
IF (ASSOCIATED(ACRIL)) ACRIL = 0.
IF (ASSOCIATED(ACRLL)) ACRLL = 0.
IF (ASSOCIATED(REVAP)) REVAP = 0.
IF (ASSOCIATED(RSUBL)) RSUBL = 0.
!!!!!!!!!!!!!! UPDATE SATURATED HUMIDITY !!!!!!!!!!!!!
dQS = DQSAT ( &
TE , &
PL, QSAT = QS )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Calculate relative edge heights
!________________________________
ZE(:,:,LM+1)=0.
DO K=LM,1,-1
ZE(:,:,K)=ZE(:,:,K+1)+dZE(:,:,K)
end do
! Determine a pressure depth scale
! for mass wgtd prec (==PSCL_QP)
!----------------------------------
VMIP = SUM( (QPi+QPl)*MASS, 3 )
PSCL_QP = SUM( (QPi+QPl)*MASS*PL, 3 )
where(VMIP > 0.)
PSCL_QP = PSCL_QP/VMIP
elsewhere
PSCL_QP = 999999.
endwhere
B_SUB = 1.00
PFl(:,:,1)=QPl(:,:,1)*MASS(:,:,1)
PFi(:,:,1)=QPi(:,:,1)*MASS(:,:,1)
EVAP_DD(:,:) = 0.
SUBL_DD(:,:) = 0.
IF(PRESENT(DDRFC)) THEN
DDFRACT(:,:) = DDRFC
ELSE
DDFRACT(:,:) = 0.00
ENDIF
DO K=2,LM
RHCR = RHCR3(:,:,K)
QPl(:,:,K) = QPl(:,:,K) + PFl(:,:,K-1) * iMASS(:,:,K)
PFl(:,:,K-1) = 0.00
QPl(:,:,K-1) = 0.00
QPi(:,:,K) = QPi(:,:,K) + PFi(:,:,K-1) * iMASS(:,:,K)
PFi(:,:,K-1) = 0.00
QPi(:,:,K-1) = 0.00
ACCR(:,:) = B_SUB * C_ACC* &
( QPl(:,:,K)*MASS(:,:,K) &
) *QCl(:,:,K)
ACCR = MIN( ACCR , QCl(:,:,K) )
QPl(:,:,K) = QPl(:,:,K) + ACCR
QCl(:,:,K) = QCl(:,:,K) - ACCR
IF (ASSOCIATED(ACRLL)) THEN
ACRLL(:,:,K) = ACCR / DT
END IF
!! Accretion of liquid condensate by falling ice/snow
ACCR(:,:) = B_SUB * C_ACC* &
( QPi(:,:,K)*MASS(:,:,K) &
) *QCl(:,:,K)
ACCR = MIN( ACCR , QCl(:,:,K) )
QPi(:,:,K) = QPi(:,:,K) + ACCR
QCl(:,:,K) = QCl(:,:,K) - ACCR
!! Liquid freezes when accreted by snow
TE(:,:,K) = TE(:,:,K) + MAPL_ALHF*ACCR/MAPL_CP
IF (ASSOCIATED(ACRIL)) THEN
ACRIL(:,:,K) = ACCR / DT
END IF
RAINRAT0 = Ifactor*QPl(:,:,K)*MASS(:,:,K)/DT
SNOWRAT0 = Ifactor*QPi(:,:,K)*MASS(:,:,K)/DT
call MARSHPALMQ2(IM, JM, RAINRAT0,PL(:,:,K),DIAMrn,NRAIN,FALLrn,VErn)
call MARSHPALMQ2(IM, JM, SNOWRAT0,PL(:,:,K),DIAMsn,NSNOW,FALLsn,VEsn)
where( FRLAND < 0.1 )
DIAMsn = MAX( DIAMsn, 1.0e-3 ) ! Over Ocean
endwhere
TinLAYERrn = dZE(:,:,K) / ( FALLrn+0.01 )
TinLAYERsn = dZE(:,:,K) / ( FALLsn+0.01 )
!*****************************************
! Melting of Frozen precipitation
!*****************************************
TAU_FRZ = 5000. ! time scale for freezing (s).
MLTFRZ = 0.0
where ( (TE(:,:,K) > MAPL_TICE ) .and.(TE(:,:,K) <= MAPL_TICE+5. ) )
MLTFRZ = TinLAYERsn * QPi(:,:,K) *( TE(:,:,K) - MAPL_TICE ) / TAU_FRZ
MLTFRZ = MIN( QPi(:,:,K) , MLTFRZ )
TE(:,:,K) = TE(:,:,K) - MAPL_ALHF*MLTFRZ/MAPL_CP
QPl(:,:,K) = QPl(:,:,K) + MLTFRZ
QPi(:,:,K) = QPi(:,:,K) - MLTFRZ
endwhere
IF (ASSOCIATED(FRZ_DIAG)) FRZ_DIAG(:,:,K) = FRZ_DIAG(:,:,K) - MLTFRZ / DT
MLTFRZ = 0.0
where ( TE(:,:,K) > MAPL_TICE+5. ) ! Go Ahead and melt any snow/hail left above 5 C
MLTFRZ = QPi(:,:,K)
TE(:,:,K) = TE(:,:,K) - MAPL_ALHF*MLTFRZ/MAPL_CP
QPl(:,:,K) = QPl(:,:,K) + MLTFRZ
QPi(:,:,K) = QPi(:,:,K) - MLTFRZ
endwhere
IF (ASSOCIATED(FRZ_DIAG)) FRZ_DIAG(:,:,K) = FRZ_DIAG(:,:,K) - MLTFRZ / DT
MLTFRZ = 0.0
if ( K >= LM-1 ) THEN
where ( TE(:,:,K) > MAPL_TICE+0. ) ! Go Ahead and melt any snow/hail left above 0 C in lowest layers
MLTFRZ = QPi(:,:,K)
TE(:,:,K) = TE(:,:,K) - MAPL_ALHF*MLTFRZ/MAPL_CP
QPl(:,:,K) = QPl(:,:,K) + MLTFRZ
QPi(:,:,K) = QPi(:,:,K) - MLTFRZ
endwhere
endif
IF (ASSOCIATED(FRZ_DIAG)) FRZ_DIAG(:,:,K) = FRZ_DIAG(:,:,K) - MLTFRZ / DT
!*****************************************
! Freezing of liquid precipitation
!*****************************************
MLTFRZ = 0.0
where ( TE(:,:,K) <= MAPL_TICE )
TE(:,:,K) = TE(:,:,K) + MAPL_ALHF*QPl(:,:,K)/MAPL_CP
QPi(:,:,K) = QPl(:,:,K) + QPi(:,:,K)
MLTFRZ = QPl(:,:,K)
QPl(:,:,K) = 0.
endwhere
IF (ASSOCIATED(FRZ_DIAG)) FRZ_DIAG(:,:,K) = FRZ_DIAG(:,:,K) + MLTFRZ / DT
! ******************************************
! In the exp below, evaporation time
! scale is determined "microphysically"
! from temp, press, and drop size. In this
! context C_EV becomes a dimensionless
! fudge-fraction.
! Also remember that these microphysics
! are still only for liquid.
! ******************************************
QKo = QV(:,:,K)
TKo = TE(:,:,K)
QPlKo = QPl(:,:,K)
QPiKo = QPi(:,:,K)
do itr = 1,3
!!! do itr = 1,1
! DQSTKo = DQSAT ( &
! TKo , &
! PL(:,:,K),QSAT=QSTko )
DQSTKo = dQS(:,:, K)
QSTKo = QS(:,:,K) + DQSTKo * ( TKo - TE(:,:,K) )
QSTKo = MAX( QSTKo , 1.0e-7 )
RH_BOX = QKo/QSTKo
RH_BOX = MAX( RH_BOX , 0.55 ) ! attempt to increase precip over arid land. JTB- 11/09/2007
where( TE(:,:,K) > 294. )
RH_BOX = MAX( RH_BOX , 0.70 ) !
endwhere
QKo = QV(:,:,K)
TKo = TE(:,:,K)
WHERE ( RH_BOX < RHCR )
Efactor = RHO_W * ( AA(:,:,K) + BB(:,:,K) ) / (RHCR - RH_BOX )
elsewhere
Efactor = 9.99e9
endwhere
where( FRLAND < 0.1 )
LANDSEAF = 0.5 ! Over Ocean
elsewhere
LANDSEAF = 0.5 ! Over Land
endwhere
!!!!! RAin falling !!!!!!!!!!!!!!!!!!!!!!!
WHERE ( ( RH_BOX < RHCR ) .AND. ( DIAMrn > 0.00 ) .AND. ( PL(:,:,K) > 100. ) .AND. ( PL(:,:,K) < REVAP_OFF_P ) )
DROPRAD=0.5*DIAMrn
T_ED = Efactor * DROPRAD**2
T_ED = T_ED * ( 1.0 + DQSTKo*MAPL_ALHL/MAPL_CP )
EVAP = QPl(:,:,K)*(1.0 - EXP( -C_EV * VErn * LANDSEAF * ENVFRAC(:,:,K) * TinLAYERrn / T_ED ) )
ELSEwhere
EVAP = 0.0
ENDwhere
!!!!! Snow falling !!!!!!!!!!!!!!!!!!!!!!!
WHERE ( ( RH_BOX < RHCR ) .AND. ( DIAMsn > 0.00 ) .AND. ( PL(:,:,K) > 100. ) .AND. ( PL(:,:,K) < REVAP_OFF_P ) )
FLAKRAD=0.5*DIAMsn
T_ED = Efactor * FLAKRAD**2
T_ED = T_ED * ( 1.0 + DQSTKo*MAPL_ALHS/MAPL_CP )
SUBL = QPi(:,:,K)*(1.0 - EXP( -C_EV * VEsn * LANDSEAF * ENVFRAC(:,:,K) * TinLAYERsn / T_ED ) )
ELSEwhere
SUBL = 0.0
ENDwhere
if (itr == 1) then
EVAPx = EVAP
SUBLx = SUBL
else
EVAP = (EVAP+EVAPx) /2.0
SUBL = (SUBL+SUBLx) /2.0
endif
QKo=QV(:,:,K) + EVAP + SUBL
TKo=TE(:,:,K) - EVAP * MAPL_ALHL / MAPL_CP - SUBL * MAPL_ALHS / MAPL_CP
enddo
! Knock out weak drizzle near surface if box is dry enough (keep for MERRA??)
! if ( K >= LM-3 ) then
! where ( ( QPl(:,:,K)<1.0e-5 ) .AND. ( RH_BOX < RHCR ) )
! EVAP = QPl(:,:,K)
! endwhere
! end if
QPi(:,:,K) = QPi(:,:,K) - SUBL
QPl(:,:,K) = QPl(:,:,K) - EVAP
!! Put some re-evap/re-subl precip in to a \quote{downdraft} to be applied later
EVAP_DD = EVAP_DD + DDFRACT*EVAP*MASS(:,:,K)
EVAP = EVAP - DDFRACT*EVAP
SUBL_DD = SUBL_DD + DDFRACT*SUBL*MASS(:,:,K)
SUBL = SUBL - DDFRACT*SUBL
! -----
QV(:,:,K) = QV(:,:,K) + EVAP + SUBL
TE(:,:,K) = TE(:,:,K) - EVAP * MAPL_ALHL / MAPL_CP - SUBL * MAPL_ALHS / MAPL_CP
IF (ASSOCIATED(REVAP)) THEN
REVAP(:,:,K) = EVAP / DT
END IF
IF (ASSOCIATED(RSUBL)) THEN
RSUBL(:,:,K) = SUBL / DT
END IF
PFl(:,:,K) = QPl(:,:,K)*MASS(:,:,K)
PFi(:,:,K) = QPi(:,:,K)*MASS(:,:,K)
IF (ASSOCIATED(PFL_DIAG)) PFL_DIAG(:,:,K) = PFl(:,:,K)/DT
IF (ASSOCIATED(PFI_DIAG)) PFI_DIAG(:,:,K) = PFi(:,:,K)/DT
END DO
! Now re-add evap and subl from \quote(downdraft) to QV and TE field
! QDDF3 is fraction of total added in a given layer. Vertical sum
! should add up to 1.00
!---------------------------------------------------------------------------
where( ZE(:,:,1:LM) < 3000. )
QDDF3 = -( ZE(:,:,1:LM)-3000. ) * ZE(:,:,1:LM) * MASS
elsewhere
QDDF3 = 0.
endwhere
VMIP = SUM(QDDF3, 3 )
DO K=1,LM
QDDF3(:,:,K) = QDDF3(:,:,K) / VMIP
end do
vmip = sum( qddf3 , 3) ! check
DO K=1,LM
EVAP = QDDF3(:,:,K)*EVAP_DD(:,:)/MASS(:,:,K)
SUBL = QDDF3(:,:,K)*SUBL_DD(:,:)/MASS(:,:,K)
QV(:,:,K) = QV(:,:,K) + EVAP + SUBL
TE(:,:,K) = TE(:,:,K) - EVAP * MAPL_ALHL / MAPL_CP - SUBL * MAPL_ALHS / MAPL_CP
IF (ASSOCIATED(REVAP)) THEN
REVAP(:,:,K) = REVAP(:,:,K) + EVAP / DT
END IF
IF (ASSOCIATED(RSUBL)) THEN
RSUBL(:,:,K) = RSUBL(:,:,K) + SUBL / DT
END IF
END DO
RAIN = PFl(:,:,LM)/DT
SNOW = PFi(:,:,LM)/DT
! For budgets etc. clear precipitating cond out layer LM
QPi(:,:,LM) = 0.
QPl(:,:,LM) = 0.
end subroutine precip3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine ICEFALL( QI, DZ, QP, VF, F, DT, FRACTION_REMOVAL ) 2
integer, intent(in) :: FRACTION_REMOVAL
real, dimension(:,:,:) :: QI , QP, DZ, F, VF
real :: DT
real, dimension(size(QI,1),size(QI,2),size(QI,3)) :: QIxP
QIxP = QI * ( VF * DT / DZ )
QIxP = MIN( QIxP , QI )
QIxP = MAX( QIxP, 0.0 ) ! protects against precision problem
QP = QP + QIxP
QI = QI - QIxP
SELECT CASE( FRACTION_REMOVAL )
CASE( 0 )
! do NOTHING
CASE( 1 )
where( (QI + QIxP) > 0. )
F = QI * F / (QI + QIxP )
endwhere
CASE( 2 )
where( (QI + QIxP) > 0. .and. QI > 0. )
F = F * SQRT( QI / (QI + QIxP ) )
endwhere
CASE( 3 )
where( (QI + QIxP) > 0. .and. QI > 0. )
F = F * ( QI / (QI + QIxP ) )**0.333
endwhere
END SELECT
end subroutine ICEFALL
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine SETTLE_VEL( IM, JM, LM, WXR, QI , PL, TE , F , KH , VF , LARGESCALE, ANVIL ) 2
integer, intent(in) :: IM,LM,JM
real, intent(in ) :: WXR
real, intent(in ) , dimension(IM,JM,LM) :: QI , F, TE, PL
real, intent(in ) , dimension(IM,JM,0:LM) :: KH
real, intent(out ) , dimension(IM,JM,LM) :: VF
integer, intent(in), optional :: ANVIL, LARGESCALE
real, dimension(IM,JM,LM) :: RHO, XIm,LXIm
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Uses Eq. 1 Lawrence and Crutzen (1998, Tellus 50B, 263-289)
! Except midlat form is taken to be for LS cloud, and tropical
! form is taken to be for ANVIL cloud
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
RHO = 1000.*100.*PL/(MAPL_RGAS*TE) ! 1000 TAKES TO g m^-3 ; 100 takes mb TO Pa
where ( ( F > 0.) .and. ( QI > 0. ) )
XIm = (QI/F)*RHO
elsewhere
XIm = 0.
endwhere
WHERE( XIm > 0.)
LXIm = ALOG10(XIm)
elsewhere
LXIm = 0.0
endwhere
IF (PRESENT(ANVIL)) then
VF = 128.6 + 53.2*LXIm + 5.5*LXIm**2
ENDIF
IF (PRESENT(LARGESCALE)) then
VF = 109.0*(XIm**0.16)
ENDIF
!VF = VF*100./MAX(PL,10.) ! Reduce/increase fall speeds for high/low pressure (NOT in LC98!!! )
! Assume unmodified they represent situation at 100 mb
if (WXR > 0.) then
VF = VF * ( 100./MAX(PL,10.) )**WXR
endif
VF = VF/100.
where( KH(:,:,0:LM-1) > 2.0 )
VF = 0.01 * VF
endwhere
!where(PL > 700.)
! VF = 0.1*VF
!endwhere
!where( (TE >= 250.) .and. (TE<260.))
! VF = ( ( -0.75/10.)*(TE-250.) + 1.0 ) * VF
!endwhere
!where(TE >= 260.)
! VF = 0.25*VF
!endwhere
end SUBROUTINE SETTLE_VEL
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine VENTILATION(PRESS,TEMPR,DIAM,W,VENT)
real, intent(in ) :: PRESS,TEMPR,DIAM,W ! in kg m**-2 s**-1
real, intent(out) :: VENT
real :: RHO,SCHM,ALPH,BET,RE,Wc,DIAMc,PR,DIFFx
real, parameter :: mu = 2.0e-4 ! "g cm^-1 s^-1 "
real, parameter :: diff = 0.25 ! "cm^2 s^-1 "
PR = PRESS*100.
RHO=(1.0e-3)*PR/(MAPL_RGAS*TEMPR)
DIFFx = DIFF*1000./PRESS
SCHM = mu/(rho*diffx)
SCHM = SCHM**0.3333
Wc = W*100.
DIAMc = DIAM*100.
re = rho*Wc*DIAMc/mu
alph=1.00
bet =0.22
VENT = alph + bet * SCHM * sqrt(Re)
end subroutine VENTILATION
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine MARSHPALM(RAIN,DIAM1,DIAM3,NTOTAL,W)
real, intent(in ) :: RAIN ! in kg m**-2 s**-1
real, intent(out) :: DIAM1,DIAM3,NTOTAL,W
real :: RAIN_HR,LAMBDA,A,B
real, parameter :: N0 = 0.08 ! # cm**-3
RAIN_HR = RAIN * 3600.
IF ( RAIN_HR <= 0.00 ) THEN
DIAM1 = 0.00
DIAM3 = 0.00
NTOTAL= 0.00
W = 0.00
RETURN
END IF
LAMBDA=41.*(RAIN_HR**(-0.21))
ntotal = 0.08 / lambda
DIAM1 = 1.0 / ( lambda**2 ) /ntotal
DIAM3 = 6.0 / ( lambda**4 ) /ntotal
DIAM3 = DIAM3**0.3333
A=2115.
B=0.8
W=A*(DIAM3**B)
! Change back to MKS units
DIAM1 = DIAM1/100.
DIAM3 = DIAM3/100.
W = W/100.
NTOTAL = NTOTAL*1.0e6
end subroutine MARSHPALM
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine MARSHPALMX(RAIN,DIAM1,DIAM3,NTOTAL,W)
real, intent(in ) :: RAIN ! in kg m**-2 s**-1
real, intent(out) :: DIAM1,DIAM3,NTOTAL,W
real :: RAIN_HR,LAMBDA,A,B
real, parameter :: N0 = 0.08 ! # cm**-3
RAIN_HR = RAIN * 3600.
IF ( RAIN_HR <= 0.00 ) THEN
DIAM1 = 0.00
DIAM3 = 0.00
NTOTAL= 0.00
W = 0.00
RETURN
END IF
DIAM1 = 2.0e-3
DIAM3 = 2.0e-3
W = 5.00
NTOTAL = 0.1*1.0e6
end subroutine MARSHPALMX
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine MARSHPALMQ2(IM,JM,RAIN,PR,DIAM3,NTOTAL,W,VE) 2
integer, intent(in) :: IM, JM
real, dimension(IM,JM), intent(in ) :: RAIN,PR ! in kg m**-2 s**-1, mbar
real, dimension(IM,JM), intent(out) :: DIAM3,NTOTAL,W,VE
real, dimension(IM,JM) :: RAIN_DAY,LAMBDA,A,B,SLOPR,DIAM1
real, parameter :: N0 = 0.08 ! # cm**-3
real :: RX(8) , D3X(8)
INTEGER :: IQD
! Marshall-Palmer sizes at different rain-rates: avg(D^3)
RX = (/ 0. , 5. , 20. , 80. , 320. , 1280., 4*1280., 16*1280. /) ! rain per in mm/day
D3X= (/ 0.019, 0.032, 0.043, 0.057, 0.076, 0.102, 0.137 , 0.183 /)
RAIN_DAY = RAIN * 3600. *24.
where( RAIN_DAY <= 0.00 )
DIAM1 = 0.00
DIAM3 = 0.00
NTOTAL= 0.00
W = 0.00
ENDwhere
DO IQD = 1,7
where( (RAIN_DAY <= RX(IQD+1)) .AND. (RAIN_DAY > RX(IQD) ) )
SLOPR =( D3X(IQD+1)-D3X(IQD) ) / ( RX(IQD+1)-RX(IQD) )
DIAM3 = D3X(IQD) + (RAIN_DAY-RX(IQD))*SLOPR
endwhere
ENDDO
where ( RAIN_DAY >= RX(8) )
DIAM3=D3X(8)
endwhere
NTOTAL = 0.019*DIAM3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
DIAM3 = 0.664 * DIAM3 !! DRYING/EVAP SHOULD PROBABLY GO AS !!
!! D_1.5 == <<D^(3/2)>>^(2/3) NOT AS !!
!! D_3 == <<D^3>>^(1/3) !!
!! RATIO D_1.5/D_3 =~ 0.66 (JTB 10/17/2002) !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
W = (2483.8 * DIAM3 + 80.)*SQRT(1000./PR)
!VE = 1.0 + 28.0*DIAM3
VE = MAX( 0.99*W/100. , 1.000 )
DIAM1 = 3.0*DIAM3
! Change back to MKS units
DIAM1 = DIAM1/100.
DIAM3 = DIAM3/100.
W = W/100.
NTOTAL = NTOTAL*1.0e6
end subroutine MARSHPALMQ2
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine MARSHPALMQ(RAIN,PR,DIAM3,NTOTAL,W,VE)
real, intent(in ) :: RAIN,PR ! in kg m**-2 s**-1, mbar
real, intent(out) :: DIAM3,NTOTAL,W,VE
real :: RAIN_DAY,LAMBDA,A,B,SLOPR,DIAM1
real, parameter :: N0 = 0.08 ! # cm**-3
real :: RX(8) , D3X(8)
INTEGER :: IQD
! Marshall-Palmer sizes at different rain-rates: avg(D^3)
RX = (/ 0. , 5. , 20. , 80. , 320. , 1280., 4*1280., 16*1280. /) ! rain per in mm/day
D3X= (/ 0.019, 0.032, 0.043, 0.057, 0.076, 0.102, 0.137 , 0.183 /)
! D3X= (/ 0.045, 0.074, 0.098, 0.132, 0.178, 0.238, 0.318 , 0.423 /) ! Old incorrectly scaled values
RAIN_DAY = RAIN * 3600. *24.
IF ( RAIN_DAY <= 0.00 ) THEN
DIAM1 = 0.00
DIAM3 = 0.00
NTOTAL= 0.00
W = 0.00
RETURN
END IF
DO IQD = 1,7
IF( (RAIN_DAY <= RX(IQD+1)) .AND. (RAIN_DAY > RX(IQD) ) ) THEN
SLOPR =( D3X(IQD+1)-D3X(IQD) ) / ( RX(IQD+1)-RX(IQD) )
DIAM3 = D3X(IQD) + (RAIN_DAY-RX(IQD))*SLOPR
ENDIF
ENDDO
IF ( RAIN_DAY >= RX(8) ) DIAM3=D3X(8)
NTOTAL = 0.019*DIAM3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
DIAM3 = 0.664 * DIAM3 !! DRYING/EVAP SHOULD PROBABLY GO AS !!
!! D_1.5 == <<D^(3/2)>>^(2/3) NOT AS !!
!! D_3 == <<D^3>>^(1/3) !!
!! RATIO D_1.5/D_3 =~ 0.66 (JTB 10/17/2002) !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
W = (2483.8 * DIAM3 + 80.)*SQRT(1000./PR)
!VE = 1.0 + 28.0*DIAM3
VE = MAX( 0.99*W/100. , 1.000 )
DIAM1 = 3.0*DIAM3
! Change back to MKS units
DIAM1 = DIAM1/100.
DIAM3 = DIAM3/100.
W = W/100.
NTOTAL = NTOTAL*1.0e6
end subroutine MARSHPALMQ
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine MICRO_AA_BB(TEMP,PR,Q_SAT,AA,BB) 1
real, intent(in ) :: TEMP,PR,Q_SAT
real, intent(out) :: AA,BB
real :: E_SAT
real, parameter :: K_COND = 2.4e-2 ! J m**-1 s**-1 K**-1
real, parameter :: DIFFU = 2.2e-5 ! m**2 s**-1
real, parameter :: epsilon = MAPL_H2OMW/MAPL_AIRMW
E_SAT = 100.* PR * Q_SAT /( EPSILON + (1.0-EPSILON)*Q_SAT ) ! (100 converts
! from mbar to Pa)
AA = ( ALHX(Temp)**2 ) / ( K_COND*MAPL_RVAP*(TEMP**2) )
BB = MAPL_RVAP*TEMP / ( DIFFU*(1000./PR)*E_SAT )
end subroutine MICRO_AA_BB
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine MICRO_AA_BB_3(TEMP,PR,Q_SAT,AA,BB) 1
real, dimension(:,:,:), intent(in ) :: TEMP,PR,Q_SAT
real, dimension(:,:,:), intent(out) :: AA,BB
real, dimension(size(TEMP,1),size(TEMP,2),size(TEMP,3)) :: E_SAT
real, parameter :: K_COND = 2.4e-2 ! J m**-1 s**-1 K**-1
real, parameter :: DIFFU = 2.2e-5 ! m**2 s**-1
real, parameter :: epsilon = MAPL_H2OMW/MAPL_AIRMW
E_SAT = 100.* PR * Q_SAT /( EPSILON + (1.0-EPSILON)*Q_SAT ) ! (100 converts
! from mbar to Pa)
AA = ( GET_ALHX3(Temp)**2 ) / ( K_COND*MAPL_RVAP*(TEMP**2) )
! AA = ( MAPL_ALHL**2 ) / ( K_COND*MAPL_RVAP*(TEMP**2) )
BB = MAPL_RVAP*TEMP / ( DIFFU*(1000./PR)*E_SAT )
end subroutine MICRO_AA_BB_3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real function OPTDPTH(PL,PE0,PE1,T,LW,FRCVOL,RADIUS,N_CCN)
real, intent(in) :: PL,PE0,PE1,T,LW,FRCVOL
integer, intent(in) :: N_CCN
real, intent(inout) :: RADIUS
real :: MUU , DELZ , RHO , LWC
IF (FRCVOL < MIN_CLD_FRAC) THEN
OPTDPTH=0.0
RADIUS =0.0
RETURN
ENDIF
IF (LW < MIN_CLD_WATER) THEN
OPTDPTH=0.0
RADIUS =0.0
RETURN
ENDIF
RADIUS=PRADIUS(PL,T,LW,FRCVOL,N_CCN,RHO_o=RHO,MUU_o=MUU) ! Cloud Drop Radius in m
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Optical thickness !!
!! formulation from DelGenio !!
!! (J. Clim., 9, 270, [Eq.26]) !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
DELZ = 100.*( PE1 -PE0 ) / (RHO*MAPL_GRAV)
OPTDPTH = 3*MUU*DELZ / (2*RHO_W*RADIUS)
end function OPTDPTH
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function LDRADIUS3(PL,TE,QCL,NN) RESULT(RADIUS)
real, dimension(:,:,:), intent(in) :: TE,PL,NN,QCL
real, dimension( size(TE,1), size(TE,2), size(TE,3) ) :: RADIUS
real, dimension( size(TE,1), size(TE,2), size(TE,3) ) :: MUU,RHO
RHO = 100.*PL / (MAPL_RGAS*TE )
MUU = QCL * RHO
RADIUS = MUU/(NN*RHO_W*(4./3.)*MAPL_PI)
RADIUS = RADIUS**(1./3.) ! Equiv. Spherical Cloud Particle Radius in m
end function LDRADIUS3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real function PRADIUS(PL,T,LW,FRCVOL,N_CCN,RHO_o,MUU_o)
real, intent(in) :: PL,T,LW,FRCVOL
integer, intent(in) :: N_CCN
real, intent(out) , optional :: MUU_o, RHO_o
real :: DELZ ,LWC, RHO, MUU
IF (FRCVOL < MIN_CLD_FRAC ) THEN
PRADIUS =0.0
RETURN
ENDIF
IF (LW < MIN_CLD_WATER ) THEN
PRADIUS =0.0
RETURN
ENDIF
LWC=LW/FRCVOL
RHO = 100.*PL / (MAPL_RGAS*T )
MUU = LWC * RHO
PRADIUS = MUU/(N_CCN*RHO_W*(4./3.)*MAPL_PI)
PRADIUS = PRADIUS**(1./3.) ! Equiv. Spherical Cloud Particle Radius in m
IF (PRESENT(MUU_O)) MUU_O=MUU
IF (PRESENT(RHO_O)) RHO_O=MUU
end function PRADIUS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real function ALHX(T,T_TRANS)
real, intent(in) :: T
real, intent(in),optional :: T_TRANS
real :: T_X
if (present( T_TRANS )) then
T_X = T_TRANS
else
T_X = T_ICE_MAX
endif
if ( T < T_ICE_ALL ) ALHX=MAPL_ALHS
if ( T > T_X ) ALHX=MAPL_ALHL
if ( T <= T_X .and. T >= T_ICE_ALL ) then
ALHX = MAPL_ALHS + (MAPL_ALHL-MAPL_ALHS)*( T - T_ICE_ALL ) /( T_X - T_ICE_ALL )
endif
end function ALHX
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function ICE_FRACTION (TEMP) RESULT(ICEFRCT)
real, dimension(:,:,:), intent(in) :: TEMP
real, dimension( size(TEMP,1), size(TEMP,2), size(TEMP,3) ) :: ICEFRCT
ICEFRCT(:,:,:) = 0.00
WHERE( TEMP <= T_ICE_ALL )
ICEFRCT = 1.000
ENDWHERE
WHERE( (TEMP > T_ICE_ALL) .AND. (TEMP <= T_ICE_MAX) )
ICEFRCT = 1.00 - ( TEMP - T_ICE_ALL ) / ( T_ICE_MAX - T_ICE_ALL )
ENDWHERE
ICEFRCT = MIN(ICEFRCT,1.00)
ICEFRCT = MAX(ICEFRCT,0.00)
!!ICEFRCT = ICEFRCT**4
ICEFRCT = ICEFRCT**ICEFRPWR
end function ICE_FRACTION
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function GET_ALHX3(T,T_TRANS) RESULT(ALHX3)
real, dimension(:,:,:), intent(in) :: T
real, dimension( size(T,1), size(T,2), size(T,3) ) :: ALHX3
real ,intent(in),optional :: T_TRANS
real :: T_X
if (present( T_TRANS )) then
T_X = T_TRANS
else
T_X = T_ICE_MAX
endif
where ( T < T_ICE_ALL )
ALHX3=MAPL_ALHS
endwhere
where ( T > T_X )
ALHX3=MAPL_ALHL
endwhere
where( (T <= T_X) .and. (T >= T_ICE_ALL) )
ALHX3 = MAPL_ALHS + (MAPL_ALHL-MAPL_ALHS)*( T - T_ICE_ALL ) /( T_X - T_ICE_ALL )
endwhere
end function GET_ALHX3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real function ICEFRAC(T,T_TRANS,T_FREEZ)
real, intent(in) :: T
real, intent(in),optional :: T_TRANS
real, intent(in),optional :: T_FREEZ
real :: T_X,T_F
if (present( T_TRANS )) then
T_X = T_TRANS
else
T_X = T_ICE_MAX
endif
if (present( T_FREEZ )) then
T_F = T_FREEZ
else
T_F = T_ICE_ALL
endif
if ( T < T_F ) ICEFRAC=1.000
if ( T > T_X ) ICEFRAC=0.000
if ( T <= T_X .and. T >= T_F ) then
ICEFRAC = 1.00 - ( T - T_F ) /( T_X - T_F )
endif
end function ICEFRAC
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE SUNDQ_ICE( TEMP, F2, F3)
REAL, INTENT( IN) :: TEMP
REAL, INTENT(OUT) :: F2, F3
REAL :: XX, YY
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Ice - phase treatment from Sundquist 1988
!! (also~ Sud and Walker 1998)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF ( TEMP >= 268. ) THEN
F2 = 1.0
F3 = 1.0
ENDIF
IF ( ( TEMP >= 253. ) .AND. ( TEMP < 268. ) )THEN
F2 = 1.0 + 0.65 * SQRT( 268. - TEMP )
F3 = 1.0
ENDIF
IF ( ( TEMP >= 233. ) .AND. ( TEMP < 253. ) )THEN
F2 = 1.0
XX = ABS( TEMP - 233.) / 17.
YY = XX*( 1.0 + XX + (4./3.)*XX**2 )
F3 = 0.15 * ( 1.07 + YY/(1.0+YY) )
F3 = 1./F3
ENDIF
IF ( ( TEMP >= 213. ) .AND. ( TEMP < 233. ) )THEN
F2 = 1.0
XX = ABS( TEMP - 233.) / 17.
YY = XX*( 1.0 + XX + (4./3.)*XX**2 )
F3 = 0.15 * ( 1.07 - YY/(1.0+YY) )
F3 = 1./F3
ENDIF
IF ( ( TEMP >= 0.00 ) .AND. ( TEMP < 213. ) )THEN
F2 = 1.0
XX = ABS( 213. - 233.) / 17.
YY = XX*( 1.0 + XX + (4./3.)*XX**2 )
F3 = 0.15 * ( 1.07 - YY/(1.0+YY) )
F3 = 1./F3
ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end subroutine sundq_ice
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE SUNDQ2_ICE( TEMP, F2, F3)
REAL, INTENT( IN) :: TEMP
REAL, INTENT(OUT) :: F2, F3
REAL :: XX, YY
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Ice - phase treatment ADAPTED
!! from Sundquist 1988
!! (also~ Sud and Walker 1998).
!! Modified to give earlier and more rapid
!! increase in autoconversion in range
!! ~~260K ~< T < 273K .
!! (JTB, 1/16/2002)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
IF ( TEMP >= 273. ) THEN
F2 = 1.0
F3 = 1.0
ENDIF
IF ( ( TEMP >= 253. ) .AND. ( TEMP < 273. ) )THEN
F2 = 1.0 + 0.93 * ( ( 273. - TEMP )*(0.333) )
F3 = 1.0
ENDIF
IF ( ( TEMP >= 233. ) .AND. ( TEMP < 253. ) )THEN
F2 = 1.0
XX = ABS( TEMP - 233.) / 17.
YY = XX*( 1.0 + XX + (4./3.)*XX**2 )
F3 = 0.15 * ( 1.07 + YY/(1.0+YY) )
F3 = 1./F3
ENDIF
IF ( ( TEMP >= 213. ) .AND. ( TEMP < 233. ) )THEN
F2 = 1.0
XX = ABS( TEMP - 233.) / 17.
YY = XX*( 1.0 + XX + (4./3.)*XX**2 )
F3 = 0.15 * ( 1.07 - YY/(1.0+YY) )
F3 = 1./F3
ENDIF
IF ( ( TEMP >= 0.00 ) .AND. ( TEMP < 213. ) )THEN
F2 = 1.0
XX = ABS( 213. - 233.) / 17.
YY = XX*( 1.0 + XX + (4./3.)*XX**2 )
F3 = 0.15 * ( 1.07 - YY/(1.0+YY) )
F3 = 1./F3
ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end subroutine sundq2_ice
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
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.
!TE1=263.
TE2=200.
!RATE2= 10.
JUMP1= (RATE2-1.0) / ( ( TE0-TE1 )**0.333 )
!RATE3= 25.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE SUNDQ3_ICE3( IM,JM,LM,TEMP,RATE2,RATE3,TE1, F2, F3) 1
INTEGER, INTENT(IN) :: IM,JM,LM
REAL, INTENT( IN) :: RATE2,RATE3,TE1
REAL, DIMENSION(IM,JM,LM), INTENT( IN) :: TEMP
REAL, DIMENSION(IM,JM,LM), 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.
!TE1=263.
TE2=200.
!RATE2= 10.
JUMP1= (RATE2-1.0) / ( ( TE0-TE1 )**0.333 )
!RATE3= 25.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Ice - phase treatment !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
where ( TEMP .GE. TE0 )
F2 = 1.0
F3 = 1.0
ENDWHERE
WHERE ( ( TEMP .GE. TE1 ) .AND. ( TEMP .LT. TE0 ) )
F2 = 1.0 + JUMP1 * (( TE0 - TEMP )**0.3333)
F3 = 1.0
ENDWHERE
WHERE ( TEMP .LT. TE1 )
F2 = RATE2 + (RATE3-RATE2)*(TE1-TEMP)/(TE1-TE2)
F3 = 1.0
ENDWHERE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
F2 = MIN(F2,27.0)
end subroutine sundq3_ice3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine RADCOUPLE ( & 1
IM, JM, LM, &
TE, &
PL, &
CF, &
AF, &
QClLS, &
QCiLS, &
QClAN, &
QCiAN, &
RAD_QL, &
RAD_QI, &
RAD_CF, &
RAD_RL, &
RAD_RI, &
CLDVOL2FRC, &
NN_ANVIL,NN_ICE,NN_WARM, &
DISABLE_RAD )
integer, intent(in ) :: IM,JM,LM
real, intent(in ) :: NN_ANVIL,NN_ICE,NN_WARM
real, intent(in ) :: CLDVOL2FRC
real, dimension(IM,JM,LM), intent(in ) :: TE,PL
real, dimension(IM,JM,LM), intent(in ) :: AF,CF, QClAN, QCiAN, QClLS, QCiLS
real, dimension(IM,JM,LM), intent(out) :: RAD_QL,RAD_QI,RAD_CF,RAD_RL,RAD_RI
real, dimension(IM,JM,LM) :: RElAN, REiAN, RElLS, REiLS, QCm, NN, ss, RAD_RI_AN
real, dimension(IM,JM,LM) :: QClANm, QCiANm, QClLSm, QCiLSm, QCtot, AFx
integer, optional, intent(in) :: DISABLE_RAD
integer :: L
real :: MIN_RI, MAX_RI, MIN_RL, MAX_RL, ALPH, POLAR_RL
! Limits on Radii needed to ensure
! correct behavior of cloud optical
! properties currently calculated in
! sorad and irrad (1e-6 m = micron)
MIN_RI = 20.0e-6
MAX_RI = 75.0e-6 ! 03/23/2007 JTB - to bring SWCF/LWCF UP IN MID-LEV/MID-LAT ICE CLOUDS
!!MIN_RL = 10.0e-6 ! 03/15/2007 JTB - TO BRING MERRA SWCF UP IN LOW CLOUD
!!MIN_RL = 15.0e-6 ! 06/20/2007 JTB - SWCF WAS TOO HIGH IN LOW CLOUDS. THIS SHOULD BRING IT DOWN.
! THIS IS BAD WAY TO GET BETTER SWCF. LWP IN LOW CLOUDS
! IS TOO HIGH. WOULD BE BETTER TO TUNE AUTOCON, BUT NO
! TIME BEFORE MERRA
MIN_RL = 10.0e-6 ! 09/12/2007 JTB - BACK!! SWCF NEED TO ADDRESS LWP.
POLAR_RL= 5.0e-6 ! 11/09/2007 JTB - COLD low level clouds
MAX_RL = 21.0e-6
IF (PRESENT(DISABLE_RAD )) THEN
IF (DISABLE_RAD==1 ) THEN
RAD_QL = 0.
RAD_QI = 0.
RAD_CF = 0.
RAD_RL = 0.
RAD_RI = 0.
RETURN
ENDIF
ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Adjust Anvil fractions for
! warm clouds
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ALPH = 0.1
SS = (280.-TE)/20.
SS = MIN( 1.0 , SS )
SS = MAX( 0.0 , SS )
SS = ALPH + (SS**3) * ( 1.0 - ALPH )
AFx = AF * SS
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Total cloud fraction
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
RAD_CF = MIN( CF + AFx, 1.00 )
! Total In-cloud liquid
where( RAD_CF > 0. )
RAD_QL = ( QClLS + QClAN ) / RAD_CF
elsewhere
RAD_QL = 0.0
endwhere
RAD_QL = MIN( RAD_QL, 0.01 )
! Total In-cloud ice
where( RAD_CF >0. )
RAD_QI = ( QCiLS + QCiAN ) / RAD_CF
elsewhere
RAD_QI = 0.0
endwhere
RAD_QI = MIN( RAD_QI, 0.01 )
where (PL < 150. )
RAD_RI = MAX_RI
endwhere
where (PL >= 150. )
RAD_RI = MAX_RI*150./PL
endwhere
!! weigh in a separate R_ice for Anvil Ice according to
!
! R_net_eff = (q_anv + q_ls) / ( q_anv/R_ice_anv + q_ls/R_ice_ls )
!-------------------------------------------------------------------------
RAD_RI_AN = RAD_RI ! 40.0e-6 ! MIN_RI
where ( ( QCiLS + QCiAN ) > 0.0 )
RAD_RI_AN = ( QCiLS + QCiAN ) / ( (QCiLS/RAD_RI) + (QCiAN/20.0e-6) )
endwhere
RAD_RI = MIN( RAD_RI, RAD_RI_AN )
RAD_RI = MAX( RAD_RI, MIN_RI )
where (PL < 300. )
RAD_RL = MAX_RL
endwhere
where (PL >= 300. )
RAD_RL = MAX_RL*300./PL
endwhere
RAD_RL = MAX( RAD_RL, MIN_RL )
where ( ( PL >= 900. ) .AND. (TE < 280. ) )
RAD_RL = POLAR_RL
endwhere
where ( RAD_CF < 1.e-5 )
RAD_QL = 0.
RAD_QI = 0.
RAD_CF = 0.
endwhere
end subroutine RADCOUPLE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE GET_REFF(TL,PL,QCL,QCI,ANV,REFFL,REFFI,NCCN,BKH)
REAL, INTENT(IN) :: TL(:,:,:)
REAL, INTENT(IN) :: PL(:,:,:)
REAL, INTENT(IN) :: ANV(:,:,:)
REAL, INTENT(IN) :: QCL(:,:,:)
REAL, INTENT(IN) :: QCI(:,:,:)
REAL, INTENT(OUT) :: REFFL(:,:,:)
REAL, INTENT(OUT) :: REFFI(:,:,:)
REAL, INTENT(OUT) :: NCCN(:,:,:)
REAL, OPTIONAL, INTENT(IN) :: BKH(:,:,:)
REAL :: QC(SIZE(TL,1),SIZE(TL,2),SIZE(TL,3))
!!REAL :: NCCN(SIZE(TL,1),SIZE(TL,2),SIZE(TL,3))
REAL :: RHOD(SIZE(TL,1),SIZE(TL,2),SIZE(TL,3))
REAL :: BKHF(SIZE(TL,1),SIZE(TL,2),SIZE(TL,3))
INTEGER :: LM
REAL :: RHO_STD
REAL, PARAMETER :: RHOI=1.0 * 1000. ! DENSITY OF ICE H2O (kg/m**3)
REAL, PARAMETER :: RHOL=1.0 * 1000. ! DENSITY OF LIQ. H2O (kg/m**3)
LM=SIZE(TL,3)
RHO_STD = 100000./(MAPL_RGAS*300.)
RHOD = 100.*PL/(MAPL_RGAS*TL)
RHOD = ( RHOD/RHO_STD )
WHERE( TL > MAPL_TICE )
NCCN = N_WARM*RHOD
ENDWHERE
WHERE( (TL > MAPL_TICE-20.) .AND. (TL <=MAPL_TICE ) )
NCCN = (N_ICE + (N_WARM-N_ICE )*( TL - (MAPL_TICE-20.) )/20.)*RHOD
ENDWHERE
WHERE( TL <= MAPL_TICE-20. )
NCCN = N_ICE*RHOD
ENDWHERE
IF(PRESENT(BKH)) THEN
BKHF = BKH
BKHF(:,:,LM) = BKH(:,:,LM-1)
IF (N_PBL > 0) THEN
WHERE( BKHF >= 1. )
NCCN = N_PBL*RHOD
ENDWHERE
ENDIF
ENDIF
IF (N_ANVIL > 0) THEN
WHERE( ( ANV >= 0.001 ) .AND. (TL <=MAPL_TICE ) )
NCCN = ANV * N_ANVIL * RHOD + &
(1.0-ANV)*NCCN
ENDWHERE
ENDIF
NCCN = NCCN * 1.0e6 ! go to N/m^3 from N/cm^3
RHOD = ( RHOD*RHO_STD )
QC = QCI
QC = QC*RHOD
REFFI= 3.0 * QC / ( RHOI * 4.*MAPL_PI*NCCN )
QC = QCL
QC = QC*RHOD
REFFL= 3.0 * QC / ( RHOL * 4.*MAPL_PI*NCCN )
REFFI = REFFI**(0.3333)
REFFL = REFFL**(0.3333)
RETURN
END SUBROUTINE GET_REFF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function ARE_WE_INSIDE(I,J, I1,IN,J1,JN) RESULT(IS_INSIDE)
INTEGER, INTENT(IN) :: I1,I,IN,J1,J,JN
LOGICAL :: IS_INSIDE
if ( (I <= IN) .and. (I >= I1 ) .and. (J <= JN) .and. (J >= J1 ) ) then
IS_INSIDE = .TRUE.
ELSE
IS_INSIDE = .FALSE.
ENDIF
END FUNCTION ARE_WE_INSIDE
end module cloudnew