SUBROUTINE GETCLDS(IM, IX, levs, TGRS, QGRS, 1,3
& VVEL, RH, QST,
& slmsk, RLON, RLAT,
& CVR, CVTR, CVBR, RHCL,
& ISTRAT, PRSI, PRSL, CLDTOT, CLDCNV,
& CLDSA, MBOTA, MTOPA, CLW, NCLD)
!
USE MACHINE
, ONLY : kind_phys,kind_rad
USE PHYSCONS, PI => con_PI
implicit none
! include 'constant.h'
!
integer MCLD,NSEAL,NBIN,NLON,NLAT
!bl ... 4 tuned cld/rh curves..bl,l,m,h ... KAC APR96
PARAMETER (MCLD=4,NSEAL=2,NBIN=100,NLON=2,NLAT=4)
!
integer IM, IX, I
integer K,IREGNH,ISLA,KC,LO,NBI,KLA,ILSEA,ILONGT,KLO,IKN,ILFT,
1 IRGT,INCV,ITOP,IBTC,levs,ncld,istrat
integer MTOPA(IX,3),MBOTA(IX,3)
integer INVR,IVVA
!
!
real (kind=kind_rad) PRSI(IX,LEVS+1), PRSL(IX,LEVS)
real (kind=kind_rad) TGRS(IX,LEVS), QGRS(IX,LEVS),
& VVEL(IM,LEVS), RH(IM,LEVS)
real (kind=kind_rad) CLW(IM,LEVS), CLDSA(IM,5),
& QST(IM,LEVS),
& CLDTOT(IM,LEVS), CLDCNV(IM,LEVS)
! real (kind=kind_rad) PGR, CVR, CVTR, CVBR
real (kind=kind_rad) CVR(IM), CVTR(IM), CVBR(IM)
&, slmsk(im), rlon(im), rlat(im)
!....
!TUNE
!... ARRAY ADDED FOR RH-CL CALCULATION
! INDICES FOR LON,LAT,CLD TYPE(L,M,H), LAND/SEA RESPECTIVELY
! NLON=1-2, FOR EASTERN AND WESTERN HEMISPHERES
! NLAT=1-4, FOR 60N-30N,30N-EQU,EQU-30S,30S-60S
! LAND/SEA=1-2 FOR LAND(AND SEAICE),SEA
!
real (kind=kind_rad) RHCL (NBIN,NLON,NLAT,MCLD,NSEAL)
real (kind=kind_rad) RHCLA(NBIN,NLON, MCLD,NSEAL)
real (kind=kind_rad) RHCLD(IM,NBIN,MCLD)
real (kind=kind_rad) XLABDY(3),XLOBDY(3)
!... XLABDY = LAT BNDRY BETWEEN TUNING REGIONS,+/- XLIM FOR TRANSITION
!. XLOBDY = LON BNDRY BETWEEN TUNING REGIONS
!
real (kind=kind_rad) tem, xlatdg, xlondg
&, xlnn, xlss,rhmax,xlim,xrgt,xlft,diflo
!
c$$$ SAVE INVR,IVVA,RHMAX,XLABDY,XLOBDY,XLIM
!
DATA XLABDY / 30.E0 , 0.E0 , -30.E0 /
DATA XLOBDY / 0.E0 , 180.E0 , 360.E0 /
DATA XLIM / 5.E0 /
!TUNE
!..... INITIAL RH CRIT. SET 1 FOR OCEAN, SET 2 FOR LAND.
! INVR=0 NO LAPSE RATE INVERSION TYPE CLD, =1 WIHT IT
! IVVA=0 NO VERTICAL VELOCITY ADJ. FOR LOW CLD, =1 WITH ADJ.
DATA RHMAX/1.00E0/, INVR/1/, IVVA/1/
! DATA ISTRAT,IBL,ICONV,IEMIS,ITHK/1,1,1,0,1/
!
integer ipnGlobal,its,mype
logical :: DiagPrint = .false.
! call PhysicsGetIpnItsMype(ipnGlobal,its,mype,DiagPrint)
!TUNE
if (ncld .eq. 0) then ! Do the following only for RH clouds!
TEM = 180.0 / PI
DO I=1,IM
XLATDG = RLAT(I) * TEM
XLONDG = RLON(I) * TEM
IREGNH = 4
!.... GET RH-CLD RELATION FOR THIS LAT
DO K=1,3
IF (XLATDG .GT. XLABDY(K)) THEN
IREGNH = K
GO TO 215
END IF
ENDDO
215 CONTINUE
DO 230 ISLA=1,NSEAL
DO 230 KC=1,MCLD
DO 230 LO=1,NLON
DO 230 NBI=1,NBIN
RHCLA(NBI,LO,KC,ISLA) = RHCL(NBI,LO,IREGNH,KC,ISLA)
230 CONTINUE
!..... LINEAR TRANSITION BETWEEN LATITUDINAL REGIONS...
DO 240 KLA=1,3
XLNN = XLABDY(KLA)+XLIM
XLSS = XLABDY(KLA)-XLIM
IF (XLATDG.LT.XLNN.AND.XLATDG.GT.XLSS) THEN
DO 235 ISLA=1,NSEAL
DO 235 KC=1,MCLD
DO 235 LO=1,NLON
DO 235 NBI=1,NBIN
RHCLA (NBI,LO,KC,ISLA) =
& (RHCL(NBI,LO,KLA,KC,ISLA)-RHCL(NBI,LO,KLA+1,KC,ISLA))
& * (XLATDG-XLSS)/(XLNN-XLSS) + RHCL(NBI,LO,KLA+1,KC,ISLA)
235 CONTINUE
END IF
240 CONTINUE
!... GET RH-CLD RELATIONSHIP FOR EACH GRID POINT, INTERPOLATING
! LONGITUDINALLY BETWEEN REGIONS IF NECESSARY..
ILSEA = 1
IF (slmsk(I).LT.1.E0) THEN
ILSEA = 2 !... OPEN OCEAN POINT....
END IF
!... WHICH HEMISPHERE (E,W)
ILONGT = 1
IF (XLONDG.GT.180.E0) ILONGT = 2
DO 246 K=1,MCLD
DO 241 NBI=1,NBIN
RHCLD(I,NBI,K) = RHCLA(NBI,ILONGT,K,ILSEA)
241 CONTINUE
IKN = 0
DO 243 KLO=1,3
DIFLO = ABS(XLONDG-XLOBDY(KLO))
IF (DIFLO.LT.XLIM) THEN
IKN = KLO
GO TO 244
END IF
243 CONTINUE
GO TO 246
244 CONTINUE
ILFT = ILONGT
IRGT = ILFT + 1
IF (IRGT.GT.NLON) IRGT = 1
XLFT = XLOBDY(IKN) - XLIM
XRGT = XLOBDY(IKN) + XLIM
DO 245 NBI=1,NBIN
RHCLD (I,NBI,K) =
& (RHCLA(NBI,ILFT,K,ILSEA)-RHCLA(NBI,IRGT,K,ILSEA))
& * (XLONDG-XRGT)/(XLFT-XRGT)+RHCLA(NBI,IRGT,K,ILSEA)
245 CONTINUE
246 CONTINUE
ENDDO
endif
!TUNE
!... VERTICAL MOTION (CB/SEC) IN VVEL
!... GET MODEL DIAGNOSED CLDS
!
CALL CLDJMS(IX, IM, LEVS, NBIN, MCLD,
& QGRS, TGRS, VVEL, RH, QST,
& CVR, CVTR, CVBR, PRSI, PRSL, SLMSK,
& CLDSA, MTOPA, MBOTA, CLDTOT, CLDCNV,
& IVVA,INVR,RLAT,RHCLD,ISTRAT,CLW,NCLD)
!
RETURN
END
SUBROUTINE CLDJMS(IDIMT,IDIMS,KDIM,NBIN,MCLD, 1,2
& Q,T,VVEL,RHRH,QST,
& CV,CVT,CVB,PRSI,PRSL,SLMSK,
& CLD,MTOP,MBOT,CLDTOT,CLDCNV,
& IVVA,INVR,XLATRD,RHCLD,ISTRAT,CLW,NCLD)
!TUNE
!FPP$ NOCONCUR R
!.... FROM YH.RAD.MDL93(CLDNEW28).......
!.... LATER UPDATED FROM YH.RAD.MDL94(CLDMUL28)...22JAN94
!.... LATER UPDATED FROM YH.RAD.MDL94(CLDML28A)... 1FEB94
!.... LATER UPDATED FROM YH.RAD.MDL94(CLDML28B)... 5FEB94
!. SUBR CLDPRP REPLACED
!. ADDED VERTICAL INTERP OF CLD-RH RELATIONS(ISTRAT GT 1)
!.... LATER UPDATED FROM YH.RAD.MDL94(CLDML28E)... 11MAR94
!. SUBR CLDPRP REPLACED,GCL ADJUSTED
!.... LATER UPDATED FROM CLOUD6................... 24MAR94
!. SUBR CLDPRP , LOW ENHANCED TO OLD VALUE..0.14..
!. SUBR GCLNEW , LLYR CALCULATION ADJ TO OLD VALU(KL-1)
!. LLYRL WAS OK.. IVE REMOVED IT AND
!. REPLACED IT BY ITS EQUIVALENT, KLOWB
!.... LATER UPDATED FROM CLOUD6................... 30MAR94
!. SUBR CLDPRP , LOW AND MIDDLE (NOT CV) ENHANCED=0.10
!---------------------------------------------------------------------
! NOV., 1992 - Y.H., K.A.C., AND A.K.
! CLOUD PARAMETERIZATION PATTERNED AFTER SLINGO AND SLINGO'S
! WORK (JGR, 1991).
! STRATIFORM CLOUDS ARE ALLOWED IN ANY LAYER EXCEPT THE SURFACE
! AND UPPER STRATOSPHERE. THE RELATIVE HUMIDITY CRITERION MAY
! VARY IN DIFFERENT MODEL LAYERS.
!YH94
! OUTPUT CLOUD AMOUNTS ARE IN CLDARY(I,K), K=1 IS THE LOWEST
! MODEL LAYER, STRATIFORM (STR) AND CONVECTIVE (CNV) TYPES OF
! CLOUD ARE COMPRESSED INTO ONE WORD: CAMT = STR + 1.0E4*CNV
! LOW MARINE STRATUS AMT'S ARE FLAGED BY ADDING 2.
!YH94
!TUNE
!.. FOR ISTRAT = 0, THERE IS RH-CLD RELATION FOR EACH LAYER..
! CRIT RH COMPUTED WITHIN..
!.. FOR ISTRAT = 1, RH-CLD RELATION FROM TABLES CREATED USING
! MITCHELL-HAHN TUNING TECHNIQUE (A.F. RTNEPH OBS)
! ...STRATUS COMPUTED SIMILAR TO OLD OPNL CLDJMS.....
! EXCEPT NO CLOUD BELOW LAYER=KLOWB..APPROX 955MB
!TUNE
! CONVECTIVE CLOUDS ARE FROM MODEL CONVECTIVE SCHEME AND ARE
! NO LONGER BROKEN INTO .75,.25,.25..RATHER CC ITSELF IS USED..
! CONVECTIVE STILL TAKES PRECEDENCE OVER STRATIFORM IN RADFS
! .(IN RADIATION USE OF
! CC GIVES IMPROVEMENT TO TROPICAL MIDDLE CLD (AS DID ST+CV))
!
! CLOUDS ARE ALSO DIVIDED INTO 3 ATMOSPHERIC DOMAINS (L,M,H)
! plus a boundary layer cloud for
! DIAGNOSTIC PURPOSES. THEY ARE COMPUTED FROM RANDOM OVERLAP
! ASSUMPTION FOR SEPARATED CLOUD LAYERS AND MAXIMUM OVERLAP
! FOR ADJACENT CLOUD LAYERS. A TOTAL CLOUD FRACTION IS ALSO
! COMPUTED.
!
! H,M,L DOMAIN PRESSURE TOPS 'PTOP1(K)' VARY LINEARLY FROM
! 'PTOPC(K,1)' AT 45DEG TO 'PTOPC(K,2)' AT THE POLE
!
!bl
! 4 cloud relationships used
! so that 1-4 is for BL,L,M,H now... clouds down to layer 2
! everywhere, with no vvel filter for these lower-10percent-atm clds..
!..... the BL relationships used below LLYR, except in marine stratus regions
! where the old method was used (down to layer 3)..2Nov95 .. kac
!....... old method of marine stratus not used...3Mar96
!..... save bl cldamt (max overlap) in cld( ,5)...11apr96
!bl
!YH0898
! AUG 1998 Y-T HOU RECODED THE DEFAULT PART OF SLINGO'S
! CLOUD SCHEME PATTERNED AFTER CCM3 (KIEHL ET AL. 1998,J.CLIM;
! AND 1994,JGR). OUTPUT RH AND CONV CLOUDS USE SEPERATED
! ARRAYS CLDTOT AND CLDCNV TO AVOID PACKING-UNPACKING PROCESS.
!
!--------------------------------------------------------------------
! INPUT VARIABLES:
! PS (CB) - SURFACE PRESSURE
! Q (KG/KG) - SPECIFIC HUMIDITY
! T (DEG K) - ABSOLUTE TEMPERATURE
! VVEL(CB/SEC) - VERTICAL VELOCITY
! CV,CVT,CVB - CONV CLD FRACTION, TOP, BOTTOM PRESSURE as obtained
! - from the convection scheme layers
! SI,SL - MDL SIGMA INTERFACE AND LAYER MEAN
! SLMSK - SEA/LAND MASK ARRAY(SEA:0.,LAND:1.,SNOW:2.)
! IVVA - FLAG TO CONTROL VERTICAL VELOCITY ADJ.
! =1: WITH, =0: WITHOUT
! INVR - FLAG TO CONTROL LAPSE RATE INVERSION CLD
! =1: WITH, =0: WITHOUT
! RHMAX - UPPER LIMIT OF RELATIVE HUMIDITY TO
! FORM OVERCAST CLOUD (CLD FRACTN = 1.)
!TUNE
! --------------- MODIFY TO AS AN ARRAY (H.-M. H. JUANG)
!********XLATRD - CURRENT LATITUDE IN RADIANS (1ST DATA PT)
!******** FOR MODELS WITH DIFF LAT AT EACH PT, NEED TO
!******** USE THE LAT OF ALL POINTS....CAREFUL.....
! RHCLD - CLOUD-RH RELATIONS FROM MITCHELL+HAHN,
! USING A.F. RTNEPH ANALYSES
! ISTRAT - 0 OR 1:FOR DEFAULT OR 'RHCLD' TABLES
! IN THE STRATIFORM CLOUD CALCULATION
!TUNE
! OUTPUT VARIABLES:
! CLDTOT - VERTICAL COLUMN ARRAY OF STRATIFORM CLOUD
! CLDCNV - VERTICAL COLUMN ARRAY OF CONVECTIVE CLOUD
! CLD - CLD FRACTION IN 3 TYPES OF DOMAINS (L,M,H)
!bl AND TOTAL IN 4TH LAYER and bl takes the 5th
! MTOP,MBOT - TOP, BOTTOM LAYERS OF CLOUDS (L,M,H)
!
!--------------------------------------------------------------------
!
cmy USE MACHINE , ONLY : kind_phys,kind_rad
USE PHYSCONS, ROCP => con_ROCP, grav => con_g, RD => con_RD
&, PI => con_PI
USE COMCD1
implicit none
! include 'constant.h'
!
! real (kind=kind_rad) EPS, EPSM1, GRV
! PARAMETER (EPS=RD/RV, EPSM1=RD/RV-1.0, GRV=GRAV)
!
integer idimt, idims, kdim,ncld
integer NBIN,MCLD,IVVA,INVR,ISTRAT
! --- INPUT VARIABLES
!
real (kind=kind_rad) CV(IDIMT), CVT(IDIMT),
& CVB(IDIMT), SLMSK(IDIMT),
& PRSL(IDIMT,KDIM), PRSI(IDIMT,KDIM+1),
& T(IDIMS,KDIM), Q(IDIMS,KDIM),
& CLW(IDIMT,KDIM), RHRH(IDIMT,KDIM),
& VVEL(IDIMT,KDIM), QST(IDIMT,KDIM),
& XLATRD(IDIMT)
! --- OUTPUT VARIABLES
real (kind=kind_rad) CLD(IDIMT,5),
& CLDTOT(IDIMT,KDIM), CLDCNV(IDIMT,KDIM)
integer MTOP(IDIMT,3), MBOT(IDIMT,3)
!TUNE
!... RH-CLD RELATIONSHIPS FOR EACH POINT
real (kind=kind_rad) RHCLD(IDIMT,NBIN,MCLD)
!TUNE
!
! --- PTOPC(K,L): TOP PRESURE OF EACH CLD DOMAIN (K=1-4 ARE SFC,L,M,H;
! L=1,2 ARE LOW-LAT (<45 DEGREE) AND POLE REGIONS)
! --- WORKSPACE ---
!
integer KDIMP,LEVM1,LEVM2,KLOW,K
integer KLEV,KC,NX,NHALF,KCVB,KCVT,KK,L,LLYR11,kstr,i
integer kbase, kinver(IDIMT)
logical BITX(IDIMT), BITY(IDIMT), BITM(IDIMT), BIT1, BIT2
&, inversn(IDIMT)
real (kind=kind_rad) PRSLY(IDIMT,KDIM),
& DTHDP(IDIMT,KDIM), THETA(IDIMT,KDIM),
& DTDPM(IDIMT), CL1 (IDIMT),
& OMEG(IDIMT), CL2 (IDIMT),
& PTOP1(IDIMT,4)
real (kind=kind_rad) RHH, RHM, RHL, FAC, EXNR
&, ES, QS, TMT0, TMT15, AI, BI, QI
&, DPTOP, QINT, AA, BB, RHCR, CRH
&, CFILTR, ONEMRH, VALUE, TEM, tem1
&, CLWMIN, CLWM, CR1, CR2, CR3
&, VVC1, VVC2, DVVCLD
&, CRK, XX
&, clwt
integer KCUT(IDIMT), KBT1 (IDIMT), KTH1 (IDIMT), KBT2 (IDIMT),
& KTH2(IDIMT), ITYP1(IDIMT), ITYP2(IDIMT), KSAVE(IDIMT),
& kdomn(idimt)
cc
cc--------------------------------------------------------------------
cc
real(kind=kind_rad) cons_0 !constant
real(kind=kind_rad) cons_1pdm10 !constant
real(kind=kind_rad) cons_p0001 !constant
real(kind=kind_rad) cons_p01 !constant
real(kind=kind_rad) cons_p3 !constant
real(kind=kind_rad) cons_p95 !constant
real(kind=kind_rad) cons_1 !constant
real(kind=kind_rad) cons_50 !constant
integer ipnGlobal,its,mype
logical :: DiagPrint = .false.
! call PhysicsGetIpnItsMype(ipnGlobal,its,mype,DiagPrint)
cc
cons_0 = 0.d0 !constant
cons_1pdm10 = 1.d-10 !constant
cons_p0001 = .0001d0 !constant
cons_p01 = .01d0 !constant
cons_p3 = .3d0 !constant
cons_p95 = .95d0 !constant
cons_1 = 1.d0 !constant
cons_50 = 50.d0 !constant
cc
cc--------------------------------------------------------------------
cc
!bl
!===> BEGIN HERE ................................................
! LLYR = 1 !JFM llyr is set in CGLJMS
! if(DiagPrint) then
! print"('JFMclds341',4i5,2L5,1p3e15.7)",its,KLOWB,KLOWT,LLYR
! endif
KDIMP = KDIM + 1
LEVM1 = KDIM - 1
LEVM2 = KDIM - 2
RHH = 0.95E0 ! CRITICAL VALUE FOR HI CLD
RHM = 0.70E0 ! CRITICAL VALUE FOR MID CLD
RHL = 0.70E0 ! CRITICAL VALUE FOR LOW CLD
!
!... FIND TOP PRESSURE FOR EACH CLOUD DOMAIN
DO K=1,4
DPTOP = PTOPC(K,2) - PTOPC(K,1)
DO I=1,IDIMT
FAC = MAX(cons_0, 4.0*ABS(XLATRD(I))/PI-1.0) !constant
PTOP1(I,K) = PTOPC(K,1) + DPTOP * FAC
END DO
END DO
! --- LOW CLOUD TOP SIGMA LEVEL, COMPUTED FOR EACH LAT CAUSE
! DOMAIN DEFINITION CHANGES WITH LATITUDE...
KLOW=KDIM
DO K=KDIM,1,-1
DO I=1,IDIMT
IF (PRSI(I,K) .LT. PTOP1(I,2) * 1.0E-3) KLOW = MIN(KLOW,K)
ENDDO
ENDDO
! --- POTENTIAL TEMP AND LAYER RELATIVE HUMIDITY
!
DO K=1,KDIM
DO I=1,IDIMT
CLDTOT(I,K) = 0.0
CLDCNV(I,K) = 0.0
PRSLY(I,K) = PRSL(I,K) * 10.0
EXNR = (PRSLY(I,K)*0.001) ** (-ROCP)
THETA(I,K) = EXNR * T(I,K)
ENDDO
ENDDO
! --- POTENTIAL TEMP LAPSE RATE
DO K=1,LEVM1
DO I=1,IDIMT
DTHDP(I,K) = (THETA(I,K+1) - THETA(I,K)) /
& (PRSLY(I,K+1) - PRSLY(I,K))
ENDDO
ENDDO
! ------------------------------------------------------------------
! FIND THE STRATOSPHERE CUT OFF LAYER FOR HIGH CLOUD. IT
! IS ASSUMED TO BE ABOVE THE LAYER WITH DTHDP LESS THAN
! -0.25 IN THE HIGH CLOUD DOMAIN (FROM LOOKING AT 1 CASE).
! ------------------------------------------------------------------
DO I=1,IDIMT
KCUT(I) = LEVM2
END DO
DO K=KLOW+1,LEVM2
BIT1 = .FALSE.
DO I=1,IDIMT
IF (KCUT(I).EQ.LEVM2 .AND. PRSLY(I,K).LE.PTOP1(I,3) .AND.
& DTHDP(I,K).LT.-0.25E0) THEN
KCUT(I) = K
END IF
BIT1 = BIT1 .OR. KCUT(I).EQ.LEVM2
END DO
IF (.NOT. BIT1) GO TO 85
END DO
85 CONTINUE
!
! --- GET CONVECTIVE CLOUD TOP AND BASE LEVEL, SET LOGICAL ARRAY BITX
!
BIT1 = .FALSE.
DO I=1,IDIMT
! BITX(I) = CV(I).GT.0.0E0 .AND. CVT(I).GE.CVB(I)
!... pressure decreases upward
BITX(I) = CV(I).GT.0.0E0 .AND. CVT(I).LT.CVB(I)
BIT1 = BIT1 .OR. BITX(I)
END DO
IF (BIT1 .and. ncld .eq. 0) THEN
! IF (BIT1) THEN
DO I=1,IDIMT
KBT1(I) = 1
KTH1(I) = 1
IF (BITX(I)) THEN
! KBT1(I) = NINT(CVB(I))
! KTH1(I) = MIN(LEVM2, NINT(CVT(I)))
!....use layer pressure in mb, as converted to cb
do k=2,kdim
tem = prsly(i,k) * 0.1
if (cvt(i) .le. tem) KTH1(I) = k - 1
if (cvb(i) .le. tem) KBT1(I) = k - 1
end do
END IF
END DO
! --- PUT CONVECTIVE CLOUD INTO 'CLDCNV', NO MERGE AT THIS POINT..
DO K=2,LEVM2
DO I=1,IDIMT
IF (BITX(I) .AND. KBT1(I).LE.K .AND. KTH1(I).GE.K) THEN
CLDCNV(I,K) = CV(I)
! RHRH(I,K) = MAX(0.0E0,
! & RHRH(I,K)-CLDCNV(I,K))/(1.0E0-CLDCNV(I,K))
END IF
END DO
END DO
! The following loop commented by Moorthi on 12/29/98
!
! --- IF MEAN CVT LAYER HIGHER THAN 400MB ADD ANVIL CIRRUS
! DO I=1,IDIMT
! IF (BITX(I) .AND. PRSLY(I,KTH1(I)).LE.CVTOP) THEN
! KK = KTH1(I) + 1
! KK = KTH1(I)
! CLDCNV(I,KK) = MAX(0.0E0, MIN(1.0E0,
! & 2.0E0 * (CV(I) - 0.3E0) ))
! END IF
! END DO
!
END IF
! ------------------------------------------------------------------
IF (ISTRAT.LE.0) THEN
IF (NCLD .GT. 0) THEN
DO K=KLOWB,LEVM2
DO I=1,IDIMT
ONEMRH = MAX(cons_1pdm10,1.-RHRH(I,K)) !constant
VALUE = MAX(MIN(1000.*CLW(I,K)/ONEMRH,cons_50),cons_0) !constant
CLDTOT(I,K) = MAX(RHRH(I,K)*(1.-EXP(-VALUE)),cons_0) !constant
ENDDO
ENDDO
ELSE
! ------------------------------------------------------------------
! ....DEFAULT SCHEME .... PATTERNED AFTER SLINGO'S SCHEME
! ------------------------------------------------------------------
! CALCULATE LARGE SCALE STRATIFORM CLOUD FROM RELATIVE HUMIDITY
! AND PUT INTO ARRAY 'CLDTOT'
! ------------------------------------------------------------------
XX = -grav * grav / (0.00035E0*RD)
DO K=KLOWB,LEVM2
DO I=1,IDIMT
IF (PRSLY(I,K) .GE. PTOP1(I,2)) THEN
! --- LOW CLOUD DOMAIN
RHCR = RHL
! IF (NINT(SLMSK(I)).EQ.1)
! & RHCR = RHL - 0.10E0
ELSE
! --- MID-HIGH CLOUD DOMAIN
IF (PRSLY(I,K) .GE. PTOP1(I,3)) THEN
CRH = RHM
IF (NINT(SLMSK(I)) .EQ. 1) CRH = RHM - 0.10E0
ELSE
CRH = RHH
IF (NINT(SLMSK(I)) .EQ. 1) CRH = RHH - 0.10E0
END IF
AA = XX*PRSLY(I,K)*DTHDP(I,K) / (T(I,K)*THETA(I,K))
RHCR = 0.999E0 - (1.0E0 - CRH)
& * (1.0E0 - MIN(cons_1, MAX(cons_0, AA)) ) !constant
END IF
BB = MAX(cons_0, (RHRH(I,K)-RHCR)/(1.0E0-RHCR)) !constant
CLDTOT(I,K) = MAX(CLDTOT(I,K), MIN(cons_1, BB**2)) !constant
END DO
END DO
! ------------------------------------------------------------------
! SPECIAL TREATMENT ON LOW CLOUDS
! ------------------------------------------------------------------
DO I=1,IDIMT
BITX(I) = .FALSE.
KBT1(I) = KLOWB-1
END DO
DO K=KLOWB-1,LLYR
DO I=1,IDIMT
IF (.NOT.BITX(I)) THEN
BITX(I) = NINT(SLMSK(I)).EQ.0
& .AND. PRSLY(I,K).GE.PSTRT .AND. DTHDP(I,K).LE.CLAPKC
& .AND. RHRH(I,K+1).LE.0.6E0 .AND. RHRH(I,K+2).LE.0.6E0
KBT1(I) = K
END IF
END DO
END DO
IF (IVVA .LE. 0) GO TO 250
! --- VERTICAL VELOCITY ADJUSTMENT ON LOW CLOUDS
VVC1 = 0.0000E0
VVC2 = -0.0006E0 ! mb/sec
DVVCLD = VVC1 - VVC2
! DVVCLD = VVCLD(1) - VVCLD(2)
DO K=KLOWB,KLOW
DO I=1,IDIMT
IF (.NOT.BITX(I) .OR. K.GT.LLYR) THEN
CR1 = MIN(cons_1, MAX(cons_0, !constant
& (VVCLD(1) - 10.0E0*VVEL(I,K)) / DVVCLD ) )
CLDTOT(I,K) = CLDTOT(I,K) * SQRT(CR1)
END IF
END DO
END DO
250 IF (INVR .LE. 0) GO TO 350
! --- T INVERSION RELATED STRATUS CLOUDS
BIT1 = .FALSE.
DO I=1,IDIMT
BIT1 = BIT1 .OR. BITX(I)
END DO
IF (.NOT. BIT1) GO TO 350
! --- SMOOTH TRANSITION FOR CLOUD WHEN DTHDP BETWEEN
! CLAPSE AND CLAPSE+DCLPS (-0.05 AND -0.06)
DO K=KLOWB,KLOW+2
DO I=1,IDIMT
IF (BITX(I)) THEN
CFILTR = MAX(cons_p3, !constant
& 1.0E0 - ((CLPSE-DTHDP(I,KBT1(I))) / DCLPS) )
CLDTOT(I,K) = MIN(cons_p95, CLDTOT(I,K) * CFILTR) !constant
END IF
END DO
END DO
! DO K=KLOWB-1,KLOW+3
! DO I=1,IDIMT
! K1 = MAX(KLOWB-1,KBT1(I)-1)
! K2 = MIN(KLOW+3, KBT1(I)+3)
! IF (K.GE.K1 .AND. K.LE.K2 .AND. VVEL(I,K).GE.0.0000E0) THEN
! AA = MIN(1.0E0, MAX(0.0E0,
! & 1.0E0 - (RHL - RHRH(I,K-1)) / (RHL - 0.6E0) ))
! BB = MIN(1.0E0, MAX(0.0E0,
! & (PRSLY(I,K)-PTOP1(I,2)) / (900.0E0-PTOP1(I,2)) ))
! CST = AA * BB * CL2(I)
! CLDTOT(I,K) = MAX(CLDTOT(I,K), CST)
! END IF
! END DO
! END DO
350 CONTINUE
END IF !END of NCLD if
! ------------------------------------------------------------------
END IF !END DEFAULT SCHEME
! ------------------------------------------------------------------
IF (ISTRAT.GT.0) THEN
IF (NCLD .GT. 0) THEN ! Liquid Water based Clouds!
inversn = .false.
kinver = kdim
DO K=2,KDIM
DO I=1,IDIMT
if (prsly(i,k) .gt. 600.0 .and. (.not. inversn(I))) then
tem = T(I,K+1) - T(I,K)
if (tem .gt. 0.1 .and. T(I,K) .GT. 278.0) then
inversn(I) = .true.
kinver(I) = k
endif
endif
ENDDO
ENDDO
CLWMIN = 0.0E-6
DO K=1,KDIM
! clwt = 1.0e-6 * sl(k)
! clwt = 2.0e-6 * sl(k)
DO I=1,IDIMT
clwt = 1.0e-6 * (prsly(I,k)*0.001)
! if (clw(i,k) .gt. clwt) then
if (clw(i,k) .gt. clwt .or.
& (inversn(I) .and. k .le. kinver(I)) ) then
ONEMRH = MAX(cons_1pdm10, 1.-RHRH(I,K)) !constant
!
TEM = 1.0 / MAX(PRSLY(I,K)*0.001, cons_p01) !constant
CLWM = CLWMIN * TEM
TEM = MIN(MAX(sqrt(sqrt(onemrh*qst(i,k))),cons_p0001), !constant
& cons_1) !constant
! TEM = 2500.0 / TEM
TEM = 2000.0 / TEM
! TEM = 1000.0 / TEM
! TEM = 100.0 / TEM
! if (inversn(I) .and. k .le. kinver(I)) tem = tem * 2.0
if (inversn(I) .and. k .le. kinver(I)) tem = tem * 5.0
VALUE = MAX(MIN(TEM*(CLW(I,K)-CLWM), cons_50),cons_0) !constant
TEM = SQRT(SQRT(RHRH(I,K)))
CLDTOT(I,K) = MAX(TEM*(1.0-EXP(-VALUE)), cons_0) !constant
endif
ENDDO
ENDDO
! ------------------------------------------------------------------
! SEPARATE CLOUDS INTO 3 PRESSURE DOMAINS (L,M,H).
!bl PLUS BL cloud amount only(5th location). Within each
! OF THE DOMAINS, ASSUME CLOUD LAYERS ARE RANDOMLY OVERLAPPED.
! ...Random/Max til Autumn 2001
! VERTICAL LOCATION OF EACH TYPE OF CLOUD IS DETERMINED BY
! THE THICKEST CONTINUING CLOUD LAYERS IN THE DOMAIN.
! Convective cloud no-longer used, remove it from diagnostic calc
! -------------------------------------------------------------------
!
! --- LOOP OVER 3 CLOUD DOMAINS (L,M,H)
!bl grabbing bl cloud fraction during L=1
!..
IF (ISTRAT.GT.0) THEN
copnl KSTR = LLYR
KSTR = LLYR+1
ELSE
KSTR = 2
END IF
!
DO L=1,3
!
!bl for 1st pass extract the bl cloud..just take random cloud
! was MAX til Autumn 2001
if (l.eq.1) then
DO I=1,IDIMT
CLD (I,5) = 0.0E0
END DO
copnl llyr11=llyr-1
copnl DO k=2,llyr11
DO k=1,llyr
DO i=1,idimt
cld(I,5) = cld(I,5)+ CLDTOT(I,k)-cld(I,5)*CLDTOT(I,k)
END DO
END DO
end if
!
DO I=1,IDIMT
CLD (I,L) = 0.0E0
MTOP(I,L) = 1
MBOT(I,L) = 1
CL1 (I) = 0.0E0
CL2 (I) = 0.0E0
KBT1(I) = 1
KBT2(I) = 1
KTH1(I) = 0
KTH2(I) = 0
copnl BITM(I) = CLDTOT(I,KSTR).GT.0.0E0 .OR.
copnl & CLDCNV(I,KSTR).GT.0.0E0
BITM(I) = CLDTOT(I,KSTR).GT.0.0E0
END DO
!
DO 700 K=KSTR,KDIM
!
BIT1 = .FALSE.
DO I=1,IDIMT
! --- BITX LOGICAL ARRAY FOR LAYER K
BITX(I) = (PRSLY(I,K).GE.PTOP1(I,L+1)) .AND.
& (PRSLY(I,K).LT.PTOP1(I,L)) .AND. BITM(I)
BIT1 = BIT1 .OR. BITX(I)
! --- BITM LOGICAL ARRAY FOR LAYER K+1
copnl BITM(I) = CLDTOT(I,K+1).GT.0.0E0 .OR.
copnl & CLDCNV(I,K+1).GT.0.0E0
if (k .lt. kdim) BITM(I) = CLDTOT(I,K+1).GT.0.0E0
END DO
IF (.NOT. BIT1) GO TO 700
DO I=1,IDIMT
copnl CR1 = CLDTOT(I,K)
copnl CR2 = CLDCNV(I,K)
IF (BITX(I)) THEN
IF(KTH2(I).LE.0) THEN
! --- KTH2 LE 0 : 1ST CLD LAYER.
KBT2(I) = K
KTH2(I) = 1
ELSE
! --- KTH2 GT 0 : CONSECUTIVE CLD LAYER.
KTH2(I) = KTH2(I) + 1
ENDIF
! --- PHYSICAL CLOUD AS SEEN BY RADIATION (random)
CL2(I)=CL2(I) + CLDTOT(I,K)-CL2(I)*CLDTOT(I,K)
copnl IF (CLDCNV(I,K).GT.0.0E0) THEN
copnl CL2 (I) = MAX(CL2(I), CLDCNV(I,K))
copnl ELSE
copnl CL2 (I) = MAX(CL2(I), CLDTOT(I,K))
copnl END IF
ENDIF
END DO
BIT2 = .FALSE.
!.... BITY=TRUE IF NEXT LYR=CLEAR OR WE CHANGE CLOUD DOMAINS..
if (k .ge. kdim) then
do i=1,idimt
bity(i) = .true.
BIT2 = .true.
enddo
else
DO 650 I=1,IDIMT
BITY(I) = BITX(I) .AND. (.NOT. BITM(I)
& .OR. PRSLY(I,K+1).LT.PTOP1(I,L+1) )
BIT2 = BIT2 .OR. BITY(I)
650 CONTINUE
endif
IF (.NOT. BIT2) GO TO 700
! --- AT THE DOMAIN BOUNDARY, RANDOM OVERLAP.
! CHOOSE THE THICKEST OR THE LARGEST FRACTION AMT AS THE CLD
! LAYER IN THAT DOMAIN
DO I=1,IDIMT
IF (BITY(I)) THEN
IF (CL1(I).GT.0.0E0) THEN
KBT1(I) = INT( (CL1(I)*KBT1(I) + CL2(I)*KBT2(I))
& / (CL1(I) + CL2(I)) )
KTH1(I) = NINT( (CL1(I)*KTH1(I) + CL2(I)*KTH2(I))
& / (CL1(I) + CL2(I)) ) + 1
CL1 (I) = CL1(I) + CL2(I) - CL1(I)*CL2(I)
ELSE
KBT1(I) = KBT2(I)
KTH1(I) = KTH2(I)
CL1 (I) = CL2 (I)
ENDIF
KBT2(I) = 1
KTH2(I) = 0
CL2 (I) = 0.0E0
ENDIF
END DO
700 CONTINUE
! --- FINISH ONE DOMAIN, SAVE EFFECTIVE CLOUDS
DO I=1,IDIMT
CLD(I,L) = CL1(I)
MTOP(I,L) = MAX(KBT1(I), KBT1(I)+KTH1(I)-1)
MBOT(I,L) = KBT1(I)
END DO
END DO ! end of domain loop (L=1,3)
! -------------------------------------------------------------------
! CALCULATE TOTAL CLOUD FROM THE MULTI-LYR CLOUD ARRAY, IN A MANNER
! AS SEEN BY THE RADIATION CODE...Random OVERLAP IS USED FOR all lyrs.
! -------------------------------------------------------------------
DO I=1,IDIMT
CLD(I,4) = 0.E0
END DO
!
DO K=1,KDIM
!
DO I=1,IDIMT
CR1 = CLDTOT(I,K)
CLD(I,4) = CLD(I,4) + CR1-CLD(I,4)*CR1
END DO
END DO
!
DO I=1,IDIMT
DO K=1,5
IF (CLD(I,K).GT.1.) CLD(I,K) = 1.
END DO
END DO
ELSE ! Old Campana Clouds!
!
!TUNE
! ------------------------------------------------------------------
! CALCULATE STRATIFORM CLOUD AND PUT INTO ARRAY 'CLDTOT' USING
! THE CLOUD-REL.HUMIDITY RELATIONSHIP FROM TABLE LOOK-UP..WHERE
! TABLES OBTAINED USING K.MITCHELL FREQUENCY DISTRIBUTION TUNING
!bl (OBSERVATIONS ARE DAILY MEANS FROM US AF RTNEPH).....K.A.C.
!bl TABLES CREATED WITHOUT LOWEST 10 PERCENT OF ATMOS.....K.A.C.
! (observations are synoptic using -6,+3 window from rtneph)
! tables are created with lowest 10-percent-of-atmos, and are
! now used.. 25 october 1995 ... kac.
! ------------------------------------------------------------------
!bl.... ... determine marine stratus regions ...NOT.....
DO 831 I=1,IDIMT
BITM(I) = .TRUE.
831 CONTINUE
! THIS LOOP TO RETRIEVE CLOUD FROM RH REWRITTEN 950113 -MI
!bl DO KLEV=KLOWB,LEVM2
do klev=2,levm2
DO I=1,IDIMT
kdomn(i)=0
BITX(I)=.FALSE.
ENDDO
!mcl3 DO KC=MCLD,1,-1
if (klev.gt.llyr) then
do kc=mcld,2,-1
DO I=1,IDIMT
!mcl3 IF(PRSLY(I,KLEV).GE.PTOP1(I,KC+1)) KBASE(I)=KC
if(prsly(i,klev).ge.ptop1(i,kc)) kdomn(i)=kc
ENDDO
ENDDO
!mcl4
else
do i=1,idimt
if (bitm(i)) then
!...if true not a suspected marine stratus region, use BL curves
kdomn(i)=1
else
!...if a marine stratus region, reinstate old method (cause very
! little marine stratus probably went into tuned curves,
! since RTNEPH has trouble here)..so use L curves
!....THIS IS BY-PASSED SINCE BITM is set TRUE everywhere..kac
kdomn(i)=2
end if
enddo
endif
!mcl4
NX=0
NHALF=(NBIN+1)/2
DO I=1,IDIMT
IF(kdomn(I).LE.0.OR.KLEV.GT.KCUT(I)) THEN
CLDTOT(I,KLEV)=0.
ELSEIF(RHRH(I,KLEV).LE.RHCLD(I,1,kdomn(I))) THEN
CLDTOT(I,KLEV)=0.
ELSEIF(RHRH(I,KLEV).GE.RHCLD(I,NBIN,kdomn(I))) THEN
CLDTOT(I,KLEV)=1.
ELSE
BITX(I)=.TRUE.
KSAVE(I)=NHALF
NX=NX+1
ENDIF
ENDDO
DOWHILE(NX.GT.0)
NHALF=(NHALF+1)/2
DO I=1,IDIMT
IF(BITX(I)) THEN
CRK=RHRH(I,KLEV)
CR1=RHCLD(I,KSAVE(I),kdomn(I))
CR2=RHCLD(I,KSAVE(I)+1,kdomn(I))
IF(CRK.LE.CR1) THEN
KSAVE(I)=MAX(KSAVE(I)-NHALF,1)
ELSEIF(CRK.GT.CR2) THEN
KSAVE(I)=MIN(KSAVE(I)+NHALF,NBIN-1)
ELSE
CLDTOT(I,KLEV)=0.01*(KSAVE(I)+(CRK-CR1)/(CR2-CR1))
BITX(I)=.FALSE.
NX=NX-1
ENDIF
ENDIF
ENDDO
ENDDO
ENDDO
! ------------------------------------------------------------------
! SPECIAL TREATMENT ON LOW CLOUDS
! ------------------------------------------------------------------
DVVCLD = VVCLD(1) - VVCLD(2)
!
do 950 k=2,klow
!
DO 904 I=1,IDIMT
OMEG(I) = 10.0E0 * VVEL(I,K)
CL1 (I) = 0.0E0
904 CONTINUE
IF (IVVA .LE. 0) GO TO 920
! --- VERTICAL VELOCITY ADJUSTMENT ON LOW CLOUDS
BIT1 = .FALSE.
DO 906 I=1,IDIMT
BITX(I) = PRSLY(I,K).GE.PTOP1(I,2) .AND. CLDTOT(I,K).GT.0.0E0
BIT1 = BIT1 .OR. BITX(I)
906 CONTINUE
IF (.NOT. BIT1) GO TO 920
IF(K.GT.LLYR) THEN
DO 910 I=1,IDIMT
IF (BITX(I)) THEN
IF(OMEG(I).GE.VVCLD(1)) THEN
CLDTOT(I,K) = 0.0E0
ELSE IF(OMEG(I).GT.VVCLD(2)) THEN
CR1 = (VVCLD(1) - OMEG(I)) / DVVCLD
! CLDTOT(I,K) = CLDTOT(I,K) * CR1
CLDTOT(I,K) = CLDTOT(I,K) * SQRT(CR1)
ENDIF
ENDIF
910 CONTINUE
ENDIF
!blC --- T INVERSION RELATED STRATUS CLOUDS .. not used cause bitm=T
920 IF (INVR .LT. 1) GO TO 950
IF (K.GT.LLYR) GO TO 950
BIT1 = .TRUE.
DO 930 I=1,IDIMT
BIT1 = BIT1 .AND. BITM(I)
930 CONTINUE
IF (BIT1) GO TO 950
DO 940 I=1,IDIMT
IF (.NOT.BITM(I)) THEN
!ERROR IF (DTHDP(I,KBASE(I)).GT.CLPSE) THEN
IF (DTHDP(I,kdomn(I)).GT.CLPSE) THEN
!--- SMOOTH TRANSITION FOR CLOUD WHEN DTHDP BETWEEN
! CLAPSE AND CLAPSE+DCLPS (-0.05 AND -0.06)
!ERROR CFILTR = 1.0E0 - ((CLPSE - DTHDP(I,KBASE(I))) / DCLPS)
CFILTR = 1.0E0 - ((CLPSE - DTHDP(I,kdomn(I))) / DCLPS)
CLDTOT(I,K) = CLDTOT(I,K)*CFILTR
END IF
! --- FOR T INVERSION TYPE CLOUD, ADD FLAG VALUE OF 2.0
!1098 CLDARY(I,K) = CLDARY(I,K)+2.0E0
END IF
940 CONTINUE
950 CONTINUE
! ------------------------------------------------------------------
! SEPARATE CLOUDS INTO 3 PRESSURE DOMAINS (L,M,H).
!bl PLUS BL cloud amount only(5th location). Within each
! OF THE DOMAINS, ASSUME SEPARATED CLOUD LAYERS ARE RANDOMLY
! OVERLAPPED AND ADJACENT CLOUD LAYERS ARE MAXIMUM OVERLAPPED.
! VERTICAL LOCATION OF EACH TYPE OF CLOUD IS DETERMINED BY
! THE THICKEST CONTINUING CLOUD LAYERS IN THE DOMAIN.
! -------------------------------------------------------------------
!
! --- LOOP OVER 3 CLOUD DOMAINS (L,M,H)
!bl grabbing bl cloud fraction during L=1
!..
IF (ISTRAT.GT.0) THEN
KSTR = LLYR
ELSE
KSTR = 2
END IF
!
DO L=1,3
!
!bl for 1st pass extract the bl cloud..just take max cloud
if (l.eq.1) then
DO I=1,IDIMT
CLD (I,5) = 0.0E0
END DO
llyr11=llyr-1
DO k=2,llyr11
DO i=1,idimt
CR1 = CLDTOT(I,K)
CR2 = CLDCNV(I,K)
cr3 = cr1
if (cr2.gt.0.E0) cr3 = cr2
cld(i,5) = max(cld(i,5) , cr3)
END DO
END DO
end if
!
DO I=1,IDIMT
CLD (I,L) = 0.0E0
MTOP(I,L) = 1
MBOT(I,L) = 1
CL1 (I) = 0.0E0
CL2 (I) = 0.0E0
KBT1(I) = 1
KBT2(I) = 1
KTH1(I) = 0
KTH2(I) = 0
BITM(I) = CLDTOT(I,KSTR).GT.0.0E0 .OR.
& CLDCNV(I,KSTR).GT.0.0E0
END DO
!
DO K=KSTR,LEVM2
!
BIT1 = .FALSE.
DO I=1,IDIMT
! --- BITX LOGICAL ARRAY FOR LAYER K
BITX(I) = (PRSLY(I,K).GE.PTOP1(I,L+1)) .AND.
& (PRSLY(I,K).LT.PTOP1(I,L)) .AND. BITM(I)
BIT1 = BIT1 .OR. BITX(I)
! --- BITM LOGICAL ARRAY FOR LAYER K+1
BITM(I) = CLDTOT(I,K+1).GT.0.0E0 .OR.
& CLDCNV(I,K+1).GT.0.0E0
END DO
IF (.NOT. BIT1) GO TO 710
DO I=1,IDIMT
CR1 = CLDTOT(I,K)
CR2 = CLDCNV(I,K)
IF (BITX(I)) THEN
IF(KTH2(I).LE.0) THEN
! --- KTH2 LE 0 : 1ST CLD LAYER.
KBT2(I) = K
KTH2(I) = 1
ELSE
! --- KTH2 GT 0 : CONSECUTIVE CLD LAYER.
KTH2(I) = KTH2(I) + 1
ENDIF
! --- PHYSICAL CLOUD AS SEEN BY RADIATION..CONV TAKES PRECEDENCE
! --- EXCEPT ANVIL CIRRUS NOT RANDOM OVERLAPPED WITH CV TOWER AS
! --- IN RADIATION CODE(SO HI MAY BE SLIGHT UNDERESTIMATE)....
IF (CLDCNV(I,K).GT.0.0E0) THEN
CL2 (I) = MAX(CL2(I), CLDCNV(I,K))
ELSE
CL2 (I) = MAX(CL2(I), CLDTOT(I,K))
END IF
ENDIF
END DO
BIT2 = .FALSE.
!.... BITY=TRUE IF NEXT LYR=CLEAR OR WE CHANGE CLOUD DOMAINS..
DO 640 I=1,IDIMT
BITY(I) = BITX(I) .AND. (.NOT. BITM(I)
& .OR. PRSLY(I,K+1).LT.PTOP1(I,L+1) )
BIT2 = BIT2 .OR. BITY(I)
640 CONTINUE
IF (.NOT. BIT2) GO TO 710
! --- AT THE DOMAIN BOUNDARY OR SEPARATED CLD LYRS, RANDOM OVERLAP.
! CHOOSE THE THICKEST OR THE LARGEST FRACTION AMT AS THE CLD
! LAYER IN THAT DOMAIN
DO I=1,IDIMT
IF (BITY(I)) THEN
IF (CL1(I).GT.0.0E0) THEN
KBT1(I) = INT( (CL1(I)*KBT1(I) + CL2(I)*KBT2(I))
& / (CL1(I) + CL2(I)) )
KTH1(I) = NINT( (CL1(I)*KTH1(I) + CL2(I)*KTH2(I))
& / (CL1(I) + CL2(I)) ) + 1
CL1 (I) = CL1(I) + CL2(I) - CL1(I)*CL2(I)
ELSE
KBT1(I) = KBT2(I)
KTH1(I) = KTH2(I)
CL1 (I) = CL2 (I)
ENDIF
KBT2(I) = 1
KTH2(I) = 0
CL2 (I) = 0.0E0
ENDIF
END DO
710 CONTINUE
END DO ! end of k-loop
! --- FINISH ONE DOMAIN, SAVE EFFECTIVE CLOUDS
DO I=1,IDIMT
CLD(I,L) = CL1(I)
MTOP(I,L) = MAX(KBT1(I), KBT1(I)+KTH1(I)-1)
MBOT(I,L) = KBT1(I)
END DO
END DO ! end of domain loop, L=1,3
! -------------------------------------------------------------------
! CALCULATE TOTAL CLOUD FROM THE MULTI-LYR CLOUD ARRAY, IN A MANNER
! AS SEEN BY THE RADIATION CODE. WHERE, MAX OVERLAP IS USED FOR
! VERTICALLY ADJACENT CLOUD LAYERS (BOTH STRATIFORM AND CONVECTIVE).
! FOR CONVECTION, ANY ANVIL IS CONSIDERED A SEPARATE RANDOMLY
! OVERLAPPED CLOUD..
!
! CL1, ITYP1 STORE PREVIOUS LAYER CLOUD AMOUNT AND TYPE,
! CL2, ITYP2 STORE CURRENT LAYER CLOUD AMOUNT AND TYPE.
! TYPE=0,1,2 FOR CLEAR LAYER, STRATIFORM AND CONVECTIVE CLOUD.
! -------------------------------------------------------------------
DO I=1,IDIMT
CLD(I,4) = 0.E0
ITYP1(I) = 0
CL1(I) = 0.0E0
END DO
!
DO K=1,KDIM
!
DO I=1,IDIMT
CL2(I) = 0.0E0
ITYP2(I) = 0
IF (CLDTOT(I,K) .GT. 0.0E0) THEN
CL2(I) = CLDTOT(I,K)
ITYP2(I) = 1
END IF
IF (CLDCNV(I,K) .GT. 0.0E0) THEN
CL2(I) = CLDCNV(I,K)
ITYP2(I) = 2
END IF
END DO
!
DO I=1,IDIMT
IF (ITYP1(I) .EQ. ITYP2(I)) THEN
IF (ITYP1(I).EQ.2 .AND. CL1(I).NE.CL2(I)) THEN
! --- CNV ANVIL CLOUDS ARE TREATED AS SEPARATED AND RADOMLY
! OVERLAP
CLD(I,4) = CLD(I,4) + CL1(I) - CLD(I,4)*CL1(I)
CL1(I) = CL2(I)
ELSE IF (ITYP1(I).EQ.1) THEN
! --- MAX OVERLAP FOR NON CONVECTIVE ADJACENT CLD LAYERS
CL1(I) = MAX(CL1(I), CL2(I))
END IF
ELSE
! --- DIFF TYPES CLOUD, RANDOM OVERLAP
IF (CL1(I) .GT. 0.0E0)
& CLD(I,4) = CLD(I,4) + CL1(I) - CLD(I,4)*CL1(I)
CL1(I) = CL2(I)
ITYP1(I) = ITYP2(I)
END IF
END DO
END DO
!
DO I=1,IDIMT
IF (CL1(I) .GT. 0.0E0) THEN
CR1 = CLD(I,4)
CLD(I,4) = CR1 + CL1(I) - CR1*CL1(I)
END IF
END DO
! ------------------------------------------------------------------
END IF ! If for Liquid water or campana
END IF ! ISTRAT GT 0
RETURN
END
SUBROUTINE CLDPRP(IDIMT,PRSL,PRSI,T,Q 1,2
&, CLDTOT,CLDCNV
&, ICWP,CLWP,SLMSK,KTOP,KBTM,NCLDS
&, CLDLW,TAUCL,CFAC,CLDSW,CWP,CIP,REW,REI,FICE
&, tau_rrtm, L, LP1, IMAX,lprnt)
!FPP$ NOCONCUR R
!---------------------------------------------------------------------
! FEB., 1993 - Y.H.
! CLOUD RADIATIVE PROPERTIES CALCULATIONS AFTER DAVIS (1982)
! AND HARSHVARDHAN ET AL. (1987).
! NOV., 1995 - Y.H.
! MODIFIED TO PROVIDE MIXED CLOUD OVERLAPPING SCHEME FOR
! CHOU'S SW RADIATION SCHEME (1992)_.
! AUG., 1998 - Y.H.
! MODIFIED TO ESTIMATE EFFECTIVE RADIUS OF WATER/ICE CLOUD
! DROP, FRACTION OF ICE WATER CONTENT, AND CLOUD EMISSIVITY
! FOR LW RADIATION CALCULATION (KIEHL ET AL. 1998,J.CLIM).
!--------------------------------------------------------------------
! INPUT VARIABLES:
! PRSL(I,K) - Layer Pressure K=1 IS TOA
! PRSI(I,K) - Interface PRESSURE (CB) K=1 IS TOA
! T (I,K) - ABSOLUTE TEMPERATURE, K=1 IS TOP LAYER (K)
! CLDTOT(I,K) - STRATIFORM CLOUD K=1 IS MDL SFC LAYER
! CLDCNV(I,K) - CONVECTIVE CLOUD K=1 IS MDL SFC LAYER
! --- MODIFY XLATRD TO GENERAL FOR REGIONAL AND GLOBAL (H.-M.H. JUANGK
! *** XLATRD - CURRENT LATITUDE IN RADIANS (1ST DATA PT)
! FOR MODELS WITH DIFF LAT AT EACH PT, NEED TO
! USE THE LAT OF EACH POINT....CAREFUL.....
! ICWP - FLAG INDICATES THE METHOD USED FOR CLOUD
! PROPERTIES, =0 USE T-P; =1 USE CLWP.
! CLWP - LAYER CLOUD WATER+ICE PATH (G/M**2), K=1 IS TOA
! SLMSK - LAND/SEA/ICE MASK (0:SEA.1:LAND,2:ICE)
! OUTPUT VARIABLES:
! --- CLOUDS FOR LW RAD. K=1 IS THE SFC, K=2 IS THE
! LOWEST CLOUD LAYER, AND SO ON
! KTOP,KBTM(I,K)- CLOUD TOP AND BOTTOM INDECES, KTOP AND
! KBTM VALUES FROM 1 TO L MODEL LAYERS,
! WITH VALUE OF 1 BEING THE TOP MDL LAYER
! NCLDS(I) - NO. OF SEPARATED CLOUD LAYERS IN A COLUMN
! CLDLW(I,K) - CLOUD FRACTIONS FOR LW, EMISSIVITY ADJUSTED
! *** ITYP(I,K) - TYPE OF CLOUDS, ITYP=1, AND 2 ARE FOR
! THE RH, AND CONV TYPES
! --- CLOUDS FOR SW RAD. K=1 IS THE TOP LAYER OR LEVEL
! TAUCL(I,K) - CLOUD OPTICAL DEPTH IN EVERY MODEL LAYER
! K=1 IS THE TOP LAYER (USED FOR ICWP=0 ONLY)
! CFAC(I,K) - FRACTION OF CLEAR SKY VIEW AT THE LAYER INTERFA
! CLDSW(I,K) - LAYER CLOUD FRACTIONS FOR SW
! - - FOLLOWING ARE OUTPUT VARIABLES FOR ICWP=1
! CWP - LAYER CLOUD WATER PATH (G/M**2)
! CIP - LAYER CLOUD ICE PATH (G/M**2)
! REW (I,K) - EFFECTIVE LAYER CLOUD WATER DROP SIZE (MICRON)
! REI (I,K) - EFFECTIVE LAYER CLOUD ICE DROP SIZE (MICRON)
! FICE(I,K) - FRACTION OF CLOUD ICE CONTENT
!
!--------------------------------------------------------------------
!
!- INPUT ARRAY
USE MACHINE , ONLY : kind_rad, kind_phys
USE PHYSCONS, grav => con_g, RD => con_RD, fv => con_fvirt
implicit none
! include 'constant.h'
!
integer idimt, imax, l, lp1, icwp
!
real (kind=kind_rad) CLDTOT(IDIMT,L), CLDCNV(IDIMT,L)
&, CLWP(IDIMT,L), PRSL(IMAX,L)
&, SLMSK(IDIMT), PRSI(IMAX,LP1)
&, T (IMAX,LP1), Q(IMAX,L)
! - OUTPUT ARRAY
real (kind=kind_rad) TAUCL(IMAX,L), CLDLW(IMAX,LP1)
&, CLDSW(IMAX,L), CFAC(IMAX,LP1)
&, REW (IMAX,L), REI (IMAX,L)
&, FICE(IMAX,L), CWP (IMAX,L)
&, CIP (IMAX,L)
integer KTOP(IMAX,LP1), KBTM(IMAX,LP1)
&, ITYP(IMAX,LP1), NCLDS(IMAX)
! --- WORKSPACE ---
real (kind=kind_rad) XAMT(IMAX), TAUC(IMAX), EMIS(IMAX)
&, CL1 (IMAX), CL2 (IMAX), ALFA (IMAX)
integer MTYP(IMAX), KCLD(IMAX), MBTM (IMAX)
logical BITX(IMAX), BITY(IMAX), BITW(IMAX)
& , BIT1, BIT2
real (kind=kind_rad) DELT, WGT, DELP, TCLD, TAU0, AWLW, AILW
real (kind=kind_rad) reimax, reimin, tem, icec
real (kind=kind_rad) tau_rrtm(IMAX,L)
integer i,k,kkcl,mclds,nc,minktp,maxkbt,mkbtp1,kk,nncld
&, ipts1
logical lprnt
cc
cc--------------------------------------------------------------------
cc
real(kind=kind_rad) cons_0 !constant
real(kind=kind_rad) cons_p1dm3 !constant
real(kind=kind_rad) cons_p08 !constant
real(kind=kind_rad) cons_1 !constant
integer ipnGlobal,its,mype
logical :: DiagPrint = .false.
! call PhysicsGetIpnItsMype(ipnGlobal,its,mype,DiagPrint)
cc
cons_0 = 0.d0 !constant
cons_p1dm3 = .1d-3 !constant
cons_p08 = .08d0 !constant
cons_1 = 1.d0 !constant
cc
cc--------------------------------------------------------------------
cc
!===> BEGIN HERE ................................................
DO I=1,IMAX
KCLD(I) = 2
MBTM(I) = 1
XAMT(I) = 0.0E0
ITYP(I,1) = 0
KTOP(I,1) = LP1
KBTM(I,1) = LP1
CLDLW(I,1) = 1.0E0
CFAC(I,1) = 1.0E0
END DO
DO K=2,LP1
DO I=1,IMAX
ITYP(I,K) = 0
KTOP(I,K) = 1
KBTM(I,K) = 1
CLDLW(I,K) = 0.0E0
CFAC(I,K) = 1.0E0
END DO
END DO
DO K=1,L
DO I=1,IMAX
CLDSW(I,K) = 0.0E0
TAUCL(I,K) = 0.0E0
CWP (I,K) = 0.0E0
CIP (I,K) = 0.0E0
REW (I,K) = 10.0E0
REI (I,K) = 10.0E0
FICE(I,K) = 0.0E0
tau_rrtm(I,K)= 0.0E0
END DO
END DO
!0499 IF (ICWP .EQ. 1) THEN ! 0898 - USE CLOUD WATER/ICE CONTENT
! reimax = 80.0
! reimin = 15.0
! reimax = 100.0
! reimin = 50.0
reimax = 0.0 ! JFM reimax was undefined
reimin = 0.0 ! JFM reimin was undefined
tem = (reimax-reimin)/(60.0-20.0)
DO K=1,L
DO I=1,IMAX
! DELT = 263.16 - T(I,K)
DELT = 273.16 - T(I,K)
IF (NINT(SLMSK(I)) .EQ. 1) THEN
REW(I,K) = 5.0E0 + 5.0E0*MIN(cons_1, MAX(cons_0, !constant
& DELT*0.05E0 )) ! - EFFECTIVE RADIUS FOR WATER
END IF
! REI(I,K) = 60.75 + (2.47 - (0.11 - 0.001 * DELT) * DELT)
! & * DELT
! REI(I,K) = reimin + (60.0 - DELT) * TEM
! REI(I,K) = MAX(reimin, MIN(reimax,REI(I,K)))
! WGT = MIN(1.0E0, MAX(0.0E0,
! & (SIGL(LP1-K) - 0.1E0) / 0.7E0 ))
! & (SIGL(LP1-K) - 0.3E0) / 0.5E0 ))
! & (SIGL(LP1-K) - 0.4E0) / 0.4E0 ))
! REI(I,K) = 30.0 - 20.0*WGT ! - EFFECTIVE RADIUS FOR ICE
! REI(I,K) = 50.0 - 40.0*WGT ! - EFFECTIVE RADIUS FOR ICE !02/16/2000
! REI(I,K) = 90.0 - 70.0*WGT ! - EFFECTIVE RADIUS FOR ICE
FICE(I,K) = MIN(cons_1, MAX(cons_0, DELT/20.0E0 )) !constant
! - FRACTION OF ICE IN CLOUD
CIP (I,K) = CLWP(I,K) * FICE(I,K)
CWP (I,K) = CLWP(I,K) - CIP(I,K)
ENDDO
ENDDO
! if (lprnt) then
! print *,' clwp=',clwp
! print *,' cwp1=',cwp
! print *,' cip1=',cip
! endif
!
DO K=1,L
DO I=1,IMAX
tem = (grav*prsl(i,k))/ (RD*(prsi(i,k+1)-prsi(i,k)))
DELT = T(I,K) - 273.16
if (cip(i,k) .gt. 0.0) then
icec = tem * cip(i,k) / (T(i,k)*(1.0 + fv*q(i,k)))
if (DELT .lt. -50.0) then
rei(i,k) = (1250.0/9.917) * icec ** 0.109
elseif (DELT .lt. -40.0) then
rei(i,k) = (1250.0/9.337) * icec ** 0.08
elseif (DELT .lt. -30.0) then
rei(i,k) = (1250.0/9.208) * icec ** 0.055
else
rei(i,k) = (1250.0/9.387) * icec ** 0.031
endif
! if(lprnt .and. k .eq. l) print *,' reiL=',rei(i,k),' icec=',icec
! *,' cip=',cip(i,k),' tem=',tem,' delt=',delt
! rei(i,k) = max(10.0, min(rei(i,k), 300.0))
endif
ENDDO
ENDDO
! if(lprnt) print *,' rei1=',rei
!
!0499 END IF
!
! --- LOOP OVER MDL LAYERS (BOTTOM UP)
!
! DO K=2,L
DO K=1,L
!
DO I=1,IMAX
CL1(I) = CLDTOT(I,K)
CL2(I) = CLDCNV(I,K)
END DO
! if(DiagPrint.and.its>14.and.kcld(1)==19) then
! print"('JFMclds1316',3i5,1p3e15.7)",its,k,imax,cl2(1)
! endif
!
BIT1 = .FALSE.
DO I=1,IMAX
BITX(I) = CL1(I).GT.0.001E0 .OR. CL2(I).GT.0.001E0
BIT1 = BIT1 .OR. BITX(I)
END DO
IF (BIT1) THEN
DO I=1,IMAX
! --- MTYP=1,2 FOR RH+STRATUS, AND CONV CLOUD TYPES
IF (CL2(I) .GT. 0.001E0) THEN
MTYP(I) = 2
ELSE IF (CL1(I) .GT. 0.001E0) THEN
MTYP(I) = 1
ELSE
MTYP(I) = 0
END IF
END DO
!
IF (K .LT. L) THEN
DO I=1,IMAX
! --- SET BITW FOR CLEAR GAP ABOVE CLOUD
! BITW(I) = BITX(I) .AND. CLDTOT(I,K+1).LE.0.001E0
! & .AND. CLDCNV(I,K+1).LE.0.001E0
BITW(I) = .TRUE.
END DO
END IF
BIT2 = .FALSE.
DO I=1,IMAX
BIT2 = BIT2 .OR. BITW(I)
IF (BITX(I)) THEN
KKCL = KCLD(I)
IF(ITYP(I,KKCL) .EQ. 0) THEN
ITYP(I,KKCL) = MTYP(I)
XAMT(I) = CL1(I)
IF (MTYP(I) .EQ. 2) XAMT(I) = CL2(I)
MBTM(I) = K
ELSE IF(ITYP(I,KKCL).NE.MTYP(I) .OR.
& (MTYP(I).EQ.2 .AND. XAMT(I).NE.CL2(I)) ) THEN
CLDLW(I,KKCL) = XAMT(I)
KTOP(I,KKCL) = LP1 - (K - 1)
KBTM(I,KKCL) = LP1 - MBTM(I)
ITYP(I,KKCL+1) = MTYP(I)
MBTM(I) = K
XAMT(I) = CL1(I)
IF (MTYP(I).EQ.2) XAMT(I) = CL2(I)
KCLD(I) = KKCL + 1
ELSE IF(MTYP(I).EQ.1) THEN
XAMT(I) = MAX(XAMT(I), CL1(I))
ENDIF
END IF
END DO
IF (BIT2) THEN
DO I=1,IMAX
IF (BITW(I)) THEN
KKCL = KCLD(I)
CLDLW(I,KKCL) = XAMT(I)
KTOP(I,KKCL) = LP1 - K
KBTM(I,KKCL) = LP1 - MBTM(I)
KCLD(I) = KKCL + 1
MTYP(I) = 0
MBTM(I) = 1
XAMT(I) = 0.0E0
END IF
END DO
ENDIF ! BIT2 ENDIF
ENDIF ! BIT1 ENDIF
ENDDO ! The K Loop ends here!
! --- RECORD NUM OF CLD LYRS AND FIND MAX NUM OF CLD LYRS
MCLDS = 0
DO I=1,IMAX
NCLDS(I) = KCLD(I) - 2
MCLDS = MAX(MCLDS, NCLDS(I))
END DO
! WRITE(6,221) MCLDS, IBEG
!221 FORMAT(' IN CLDPRP: MAXCLDS =',I4,' IBEG=',I4)
! IF (MCLDS .EQ. 0) RETURN
!
! --- ESTIMATE CLOUD OPTICAL PROPERTIES FROM T AND Q -- (TOP DOWN)
!
DO NNCLD=1,MCLDS
NC = MCLDS - NNCLD + 2
!
DO I=1,IMAX
BITX(I) = CLDLW(I,NC) .GE. 0.001E0
! BITY(I) = BITX(I) .AND. ITYP(I,NC).EQ.2
BITW(I) = BITX(I)
END DO
! --- FIND TOP PRESSURE FOR MID CLOUD (3) DOMAIN=FUNCTION OF LATITUDE
MINKTP=LP1
MAXKBT=1
DO I=1,IMAX
IF (BITX(I)) THEN
MINKTP = MIN(MINKTP,KTOP(I,NC))
MAXKBT = MAX(MAXKBT,KBTM(I,NC))
END IF
END DO
! write(6,241) NC,MINKTP, MAXKBT
!241 format(3x,'NC, MINKTP, MAXKBT =',3I6/3X,'BITX, BITW :')
! --- FIND CLEAR SKY VIEW AT EACH LEVELS
DO KK=MINKTP,MAXKBT
DO I=1,IMAX
IF (KK.GE.KTOP(I,NC) .AND.
& KK.LE.KBTM(I,NC) .AND. BITX(I)) THEN
CLDSW(I,KK) = CLDLW(I,NC)
IF (BITW(I)) THEN
CFAC(I,KK+1) = CFAC(I,KK) * (1.0E0 - CLDSW(I,KK))
BITW(I) = .FALSE.
ELSE
CFAC(I,KK+1) = CFAC(I,KK)
END IF
ELSEIF (KK.GT.KBTM(I,NC) .AND. BITX(I)) THEN
CFAC(I,KK+1) = CFAC(I,KK)
END IF
END DO
END DO
!
MKBTP1 = MAXKBT + 1
DO K=MKBTP1,L
DO I=1,IMAX
IF (BITX(I)) CFAC(I,K+1) = CFAC(I,MKBTP1)
END DO
END DO
!
DO I=1,IMAX
TAUC(I) = 0.0E0
END DO
IF (ICWP .NE. 1) THEN
!CONV - REDUCE CONV CLOUD AMOUNT FOR SW RAD
DO I=1,IMAX
IF (ITYP(I,NC) .EQ. 2) THEN
!0799 ALFA(I) = CLDLW(I,NC)
ALFA(I) =
& MAX(0.25E0, 1.0E0-0.125E0*(KBTM(I,NC)-KTOP(I,NC)))
ELSE
! ALFA(I) = SQRT(CLDLW(I,NC))
ALFA(I) = 1.0
END IF
END DO
!
! --- CALC CLD THICKNESS DELP AND MEAN TEMP (CELSIUS)
DO KK=MINKTP,MAXKBT
DO I=1,IMAX
IF (KK.GE.KTOP(I,NC) .AND.
& KK.LE.KBTM(I,NC) .AND. BITX(I)) THEN
DELP = 10.0 * (PRSI(I,KK+1) - PRSI(I,KK))
TCLD = T(I,KK) - 273.16E0
! --- CONVECTIVE CLOUD
IF (ITYP(I,NC) .EQ. 2) THEN
! TAU0 = DELP * 0.05E0 ! Yu Tai's new value
TAU0 = DELP * 0.06E0 ! Operational value
!0499 - IF CONV CLD, SET TO WATER CLOUD ONLY
FICE(I,KK) = 0.0E0
! --- RH CLOUDS
ELSE
IF (TCLD .LE. -10.0E0) THEN
TAU0 = DELP
& * MAX(cons_p1dm3, 2.00E-6*(TCLD+82.5E0)**2) !constant
ELSE
TAU0 = DELP
& * MIN(cons_p08, 6.949E-3*TCLD+0.08E0) !constant
END IF
END IF
TAUC(I) = TAUC(I) + TAU0
TAUCL(I,KK) = TAU0 * ALFA(I)
tau_rrtm(i,kk) = tau0
END IF
END DO
END DO
! --- CALC CLD EMIS AND EFFECTIVE CLOUD COVER FOR LW
DO I=1,IMAX
IF (BITX(I))
& CLDLW(I,NC) = CLDLW(I,NC)*(1.0E0-EXP(-0.75E0*TAUC(I)))
ENDDO
!
ELSE ! prognostic liquid water
!
! --- CALC TAUC AND EMIS
DO KK=MINKTP,MAXKBT
DO I=1,IMAX
IF (KK .GE. KTOP(I,NC) .AND.
& KK .LE. KBTM(I,NC) .AND. BITX(I)) THEN
! AWLW = 0.070000E0 * CWP(I,KK)
AWLW = 0.078000E0 * CWP(I,KK)
! AWLW = 0.090361E0 * CWP(I,KK)
! AILW = (0.005E0 + 1.0E0/REI(I,KK)) * CIP(I,KK)
AILW = (0.0029E0 + 1.0E0/REI(I,KK)) * CIP(I,KK)
TAUC(I) = TAUC(I) + AWLW + AILW
tau_rrtm(i,kk) = awlw + ailw
ENDIF
END DO
END DO
DO I=1,IMAX
IF (BITX(I))
& CLDLW(I,NC) = CLDLW(I,NC) * max(cons_0,MIN(cons_1, !constant
& 1.0E0 - EXP(-1.66E0*TAUC(I)) ) )
ENDDO
END IF
!
END DO ! End of NNCLD loop!
!
! if(lprnt) print *,' cip2=',cip
! if(lprnt) print *,' cldsw=',cldsw
IF (ICWP .EQ. 1) THEN
DO K=1,L
DO I=1,IMAX
IF (CLDSW(I,K) .LT. 0.001E0) THEN
TAUCL(I,K) = 0.0E0
CLDSW(I,K) = 0.0E0
CWP (I,K) = 0.0E0
CIP (I,K) = 0.0E0
END IF
END DO
END DO
END IF
! if(lprnt) print *,' cip3=',cip
!
! write(6,572)
!572 format(/3x,'=== CLOUD LIQUID/ICE FRACTION ===')
! do k=3,L,2
! write(6,574) k
!574 format(5x,'K =',i4,3x,'I=1,64,12')
! write(6,576) (fice(i,k),i=1,64,12)
!576 format(10f10.6)
! enddo
! WRITE(6,581) IBEG
!581 FORMAT(' IN CLDPRP: PRINT LW CLOUDS, IBEG=',I5)
! DO 586 K=2,MCLDS+1
! WRITE(6,582) K
!582 FORMAT(' FROM SFC AND UP, CLD NUM =',I5,' CLDLW,EMIS,TOP,BOT')
! WRITE(6,583) (CLDLW(I,K),EMIS(I,K),KTOP(I,K),
! 1 KBTM(I,K),I=1,IMAX)
!583 FORMAT(6(2F6.3,2I4))
!586 CONTINUE
!
! WRITE(6,591)
!591 FORMAT(' IN CLDPRP: PRINT SW CLOUDS')
! DO 596 K=KSTRT,L
! WRITE(6,592) K
!592 FORMAT(' FROM TOP AND DN, LEV NUM =',I5,' CLDSW,CFAC,TAU')
! WRITE(6,593) (CLDSW(I,K),CFAC(I,K+1),TAUCL(I,K),I=1,IMAX)
!593 FORMAT(6(3F6.3,2X))
!596 CONTINUE
!
RETURN
END
SUBROUTINE ALBAER(SLMSK,SNOWF,ZORLF,COSZF,TSEAF,HPRIF,TGDF, 1,1
& ALVSF,ALNSF,ALVWF,ALNWF,FACSF,FACWF,PI,DLTG,
& KPRFG,IDXCG,CMIXG,DENNG,XLAT,XLON,
& FICE, ! FOR SEA-ICE - XW Nov04
& ALVBR,ALNBR,ALVDR,ALNDR,KAER,KPRF,IDXC,
& CMIX,DENN,NXC,NDN,IMXAE,JMXAE,LEN,T0C)
!FPP$ NOCONCUR R
!*******************************************************************
! THIS PROGRAM COMPUTES FOUR COMPONENTS OF SURFACE ALBEDOS (I.E.
! VIS-NIR, DIRECT-DIFFUSED), from vis/nir strong/weak zenith angle
! dependency, BASED ON BRIEGLEB'S SCHEME. AND
! BILINEARLY INTERPOLATES ALBEDO AND AEROSOL DISTRIBUTION TO
! RADIATION GRID....SEASONAL INTERPOLATION DONE IN CYCLE..HLP FEB98
!
! INPUT VARIABLES:
! SLMSK - SEA(0),LAND(1),ICE(2) MASK ON FCST MODEL GRID
! SNOWF - SNOW DEPTH WATER EQUIVALENT IN MM
! ZORLF - SURFACE ROUGHNESS IN CM
! COSZF - COSIN OF SOLAR ZENITH ANGLE
! TSEAF - SEA SURFACE TEMPERATURE IN K
! HPRIF - TOPOGRAPHIC SDV IN M
! TGDF - GROUND SOIL TEMPERATURE IN K
! ALVSF - MEAN VIS ALBEDO WITH STRONG COSZ DEPENDENCY
! ALNSF - MEAN NIR ALBEDO WITH STRONG COSZ DEPENDENCY
! ALVWF - MEAN VIS ALBEDO WITH WEAK COSZ DEPENDENCY
! ALNWF - MEAN NIR ALBEDO WITH WEAK COSZ DEPENDENCY
! FACSF - FRACTIONAL COVERAGE WITH STRONG COSZ DEPENDENCY
! FACWF - FRACTIONAL COVERAGE WITH WEAK COSZ DEPENDENCY
! KPRFG, IDXCG, CMIXG, DENNG
! - GLOBAL DIST OF TROPOSPHERIC AEROSOLS DATA:
! PROFILE TYPE, COMPONENT INDEX, COMPONENT MIXING RATIO,
! AND NUMBER DENSITY (FOR FIRST AND SECOND LAYERS)
! XLON, XLAT
! - LONGITUDE AND LATITUDE OF GIVEN POINTS IN RADIANCE
!
! OUTPUT VARIABLES: (ALL ON RADIATION GRID)
! ALVBR - VIS BEAM SURFACE ALBEDO
! ALNBR - NIR BEAM SURFACE ALBEDO
! ALVDR - VIS DIFF SURFACE ALBEDO
! ALNDR - NIR DIFF SURFACE ALBEDO
! KPRF, IDXC, CMIX, DENN
! - AEROSOL DATA FOR THE SELECTED GRID POINTS
!
! Nov. 97 - MODIFIED SNOW ALBEDO (USE HPRIF, JSNO) - YH
! JAN 98 - TO INCLUDE GRUMBINE'S SEAICE SCHEME
! AND USE SLMSKR TO SETUP TWO BASIC
! TYPES OF AEROSOLS (CONT-I, MAR-I) - YH
! MAR 00 - MODIFIED TO USE OPAC AEROSOL DATA(1998) - YH
!
!******************************************************************
! --- INPUT
!
USE MACHINE , ONLY : kind_rad
implicit none
!
integer len, imxae, jmxae, nxc, ndn, kaer
! --- INPUT
real (kind=kind_rad) SLMSK(LEN), SNOWF(LEN), ZORLF(LEN)
&, TSEAF(LEN), COSZF(LEN), HPRIF(LEN)
&, ALVSF(LEN), ALNSF(LEN), ALVWF(LEN)
&, ALNWF(LEN), FACSF(LEN), FACWF(LEN)
&, TGDF(LEN), PI, DLTG
&, CMIXG(NXC,IMXAE,JMXAE), XLON(LEN)
&, DENNG(NDN,IMXAE,JMXAE), XLAT(LEN)
!c-- XW: FOR SEA-ICE Nov04
real (kind=kind_rad) FICE(LEN)
!c-- XW: END SEA-ICE
!
integer IDXCG(NXC,IMXAE,JMXAE),KPRFG(IMXAE,JMXAE)
! --- OUTPUT
real (kind=kind_rad) ALVBR(LEN), ALNBR(LEN), ALVDR(LEN)
&, ALNDR(LEN), CMIX(NXC,LEN), DENN(NDN,LEN)
&, TEMP1, TEMP2
integer IDXC(NXC,LEN), KPRF(LEN)
! --- Local VARIABLES
!c-- XW: FOR SEA-ICE Nov04
real (kind=kind_rad) FFW, TICE, TFW
SAVE TFW
DATA TFW / 271.2 /
!c-- XW: END SEA-ICE
real (kind=kind_rad) ASNVB, ASNNB, ASNVD, ASNND, ASEVB
&, ASENB, ASEVD, ASEND, FSNO, FSEA
&, RFCS, RFCW, FLND
!
&, ASNOW, ARGH, HRGH, FSNO0, FSNO1
&, FLND0, FSEA0, DTSQ, DTGD, CSNOW
&, A1, A2, B1, B2, T0C
&, tmp1, tmp2, hdlt, rdg
integer i,i1, i2, j1, j2, j3, k
cc
cc--------------------------------------------------------------------
cc
real(kind=kind_rad) cons_0 !constant
real(kind=kind_rad) cons_p025 !constant
real(kind=kind_rad) cons_p2 !constant
real(kind=kind_rad) cons_p5 !constant
real(kind=kind_rad) cons_p98 !constant
real(kind=kind_rad) cons_1 !constant
real(kind=kind_rad) cons_5 !constant
cc
cons_0 = 0.d0 !constant
cons_p025 = .025d0 !constant
cons_p2 = .2d0 !constant
cons_p5 = .5d0 !constant
cons_p98 = .98d0 !constant
cons_1 = 1.d0 !constant
cons_5 = 5.d0 !constant
cc
cc--------------------------------------------------------------------
cc
DO I=1,LEN
!
! --- MODIFIED SNOW ALBEDO SCHEME - UNITS CONVERT TO M
! (ORIGINALLY SNOWF IN MM; ZORLF IN CM)
ASNOW = 0.02*SNOWF(I)
ARGH = MIN(cons_p5, MAX(cons_p025, 0.01*ZORLF(I))) !constant
HRGH = MIN(cons_1, MAX(cons_p2, 1.0577-1.1538E-3*HPRIF(I) ) ) !constant
FSNO0 = ASNOW / (ARGH + ASNOW) * HRGH
IF (SLMSK(I).EQ.0.0 .AND. TSEAF(I).GT.271.2) FSNO0 = 0.0
FSNO1 = 1.0 - FSNO0
FLND0 = MIN(cons_1, FACSF(I)+FACWF(I)) !constant
FSEA0 = MAX(cons_0, 1.0 - FLND0) !constant
FSNO = FSNO0
FSEA = FSEA0 * FSNO1
FLND = FLND0 * FSNO1
! --- DIFFUSED SEA SURFACE ALBEDO
IF (TSEAF(I) .GE. 271.5) THEN
ASEVD = 0.06
ASEND = 0.06
ELSE IF (TSEAF(I) .LT. 271.1) THEN
ASEVD = 0.70
ASEND = 0.65
ELSE
DTSQ = (TSEAF(I) - 271.1)**2
ASEVD = 0.7 - 4.0*DTSQ
ASEND = 0.65 - 3.6875*DTSQ
END IF
! --- DIFFUSED SNOW ALBEDO
!c-- XW: FOR SEA-ICE Nov04
! IF (SLMSK(I) .EQ. 2.0) THEN
!! (Bob Grumbine's method for snow covered sea ice)
! DTGD = MAX(0.0, MIN(5.0, 273.16-TGDF(I)) )
! ASNVD = 0.70 + 0.03 * DTGD
! ASNND = 0.60 + 0.03 * DTGD
! ASEVD = 0.70
! ASEND = 0.65
IF (SLMSK(I) .GE. 2.0) THEN
FFW = 1.0 - FICE(I)
IF (FFW .LT. 1.0) THEN
TICE = (TSEAF(I)-FFW*TFW)/FICE(I)
!mr DTGD = MAX( 0.0, MIN( 5.0, 273.16-TICE) )
DTGD = MAX(cons_0, MIN(cons_5, 273.16-TICE) ) !constant
ELSE
DTGD = 0.0
ENDIF
ASNVD = (0.70 + 0.03*DTGD)*FICE(I) + 0.06*FFW
ASNND = (0.60 + 0.03*DTGD)*FICE(I) + 0.06*FFW
ASEVD = 0.70*FICE(I) + 0.06*FFW
ASEND = 0.65*FICE(I) + 0.06*FFW
!c-- XW: END SEA-ICE
ELSE
ASNVD = 0.90
ASNND = 0.75
END IF
!
! --- DIRECT SNOW ALBEDO
IF (COSZF(I) .LT. 0.5) THEN
CSNOW = 0.5 * (3.0 / (1.0+4.0*COSZF(I)) - 1.0)
ASNVB = MIN( cons_p98, ASNVD+(1.0-ASNVD)*CSNOW ) !constant
ASNNB = MIN( cons_p98, ASNND+(1.0-ASNND)*CSNOW ) !constant
ELSE
ASNVB = ASNVD
ASNNB = ASNND
END IF
! --- DIRECT SEA SURFACE ALBEDO
IF (COSZF(I) .GT.0.0) THEN
RFCS = 1.4 / (1.0 + 0.8*COSZF(I))
!ERROR RFCS = 1.4 / (1.0 + 0.4*COSZF(I))
RFCW = 1.1 / (1.0 + 0.2*COSZF(I))
IF (TSEAF(I) .GE. T0C) THEN
ASEVB = MAX(ASEVD, 0.026/(COSZF(I)**1.7+0.065)
& + 0.15 * (COSZF(I)-0.1) * (COSZF(I)-0.5)
& * (COSZF(I)-1.0))
ASENB = ASEVB
ELSE
ASEVB = ASEVD
ASENB = ASEND
END IF
ELSE
RFCS = 1.0
RFCW = 1.0
ASEVB = ASEVD
ASENB = ASEND
END IF
!
A1 = ALVSF(I) * FACSF(I)
B1 = ALVWF(I) * FACWF(I)
A2 = ALNSF(I) * FACSF(I)
B2 = ALNWF(I) * FACWF(I)
ALVBR(I) = (A1*RFCS + B1*RFCW) *FLND
& + ASEVB*FSEA + ASNVB*FSNO
ALVDR(I) = (A1 + B1 )*FLND
& + ASEVD*FSEA + ASNVD*FSNO
ALNBR(I) = (A2*RFCS + B2*RFCW) *FLND
& + ASENB*FSEA + ASNNB*FSNO
ALNDR(I) = (A2 + B2 )*FLND
& + ASEND*FSEA + ASNND*FSNO
!
ENDDO
!
if (kaer .gt. 0) then
!
!.... AEROSOL DISTRIBUTIONS
!YH0400 USE OPAC AEROSOL DATA
RDG = 180.0 / PI
HDLT = 0.5 * DLTG
!
DO I=1,LEN
I2 = 1
J2 = 1
TMP1 = XLON(I) * RDG
DO I1=1,IMXAE
TMP2 = DLTG * (I1 - 1)
IF (TMP2 .GT. 360.0-HDLT) THEN
TMP2 = TMP2 - 360.0
END IF
IF (ABS(TMP1-TMP2) .LE. HDLT) THEN
I2 = I1
GO TO 40
END IF
END DO
40 TMP1 = XLAT(I) * RDG
DO J1=1,JMXAE
TMP2 = 90.0 - DLTG * (J1 - 1)
IF (ABS(TMP1-TMP2) .LE. HDLT) THEN
J2 = J1
GO TO 50
END IF
END DO
!
50 KPRF(I) = KPRFG(I2,J2)
DO K=1,NXC
IDXC(K,I) = IDXCG(K,I2,J2)
CMIX(K,I) = CMIXG(K,I2,J2)
END DO
DO K=1,NDN
DENN(K,I) = DENNG(K,I2,J2)
END DO
ENDDO
endif
!
RETURN
END
SUBROUTINE GCLJMS (SI,KDIM) 2,1
cmy USE MACHINE , ONLY : kind_phys,kind_rad
! USE PHYSCONS
USE COMCD1
implicit none
! include 'constant.h'
integer i,j,kl,k,kdim,kk
real (kind=kind_rad) xthk,silowxthk,silow
real (kind=kind_rad) SI(KDIM+1), PPPTOP(4,2)
! --- PRESSURE LIMITS FOR SFC AND TOP OF EACH CLOUD DOMAIN (L,M,H)
! IN MB, MODEL LAYERS FOR CLD TOPS ARE L=7,M=11,H=15 AT LOW
! LATITUDES AND L= ,M= ,H= , AT POLE REGION.
!.... PTOP ABOVE H CHANGED FROM 150 TO 100, CAUSE
!C DATA PPPTOP /1050.,642.,350.,150., 1050.,750.,500.,150./
! CODE WAS TRUNCATING TOPS OF CONVECTIVE CLOUDS
!C DATA PPPTOP /1050.,642.,350.,100., 1050.,750.,500.,100./
! use 0. for high domain, as CLW parameterization uses
! all the upper model layers for clouds
DATA PPPTOP /1050.,642.,350., 0., 1050.,750.,500., 0./
integer ipnGlobal,its,mype
logical :: DiagPrint = .false.
! call PhysicsGetIpnItsMype(ipnGlobal,its,mype,DiagPrint)
!
! ROCP = RD / CP
! --- INVERSON TYPE CLD CRITICAL VALUE-ISTRAT=0
!YH94 CLAPSE = -0.055E0
CLAPSE = -0.06E0
! --- INVERSON TYPE CLD CRITICAL VALUE-ISTRAT=1
CLAPKC = -0.05E0
!....CRITICAL DTHETA/DP FOR OCEAN STRATUS(WGT VARIES 0 TO 1
! LINEARLY FROM CLAPSE TO CLPSE)
DCLPS = -0.01E0
CLPSE = CLAPKC + DCLPS
CVTOP = 400.0E0
PSTRT = 800.0E0
! --- LOW CLD BOTTOM (AT SIGMA=0.95) AND TOP SIGMA LEVEL
!JFM DO 5 K=1,KDIM
DO 5 K=2,KDIM ! K-1 must be greater than 0
KK=K
! if(DiagPrint) then
! print"('JFMclds1872',2i5,1p3e15.7)",its,kk,si(kk)
! endif
IF (SI(KK) .LE. 0.95E0) GO TO 10
5 CONTINUE
10 KLOWB = KK - 1
SILOW = PPPTOP(2,1) * 1.0E-3
DO 20 K=1,KDIM
KK=K
IF (SI(KK) .LT. SILOW) GO TO 30
! if(DiagPrint) then
! print"('JFMclds1882',2i5,1p3e15.7)",its,kk,si(kk),silow
! endif
20 CONTINUE
30 KLOWT = KK ! KLOWT not used 3/19/07
! --- PRESURE LIMIT AT SFC AND AT TOP OF CLOUD DOMAINS (L,M,H) IN MB
DO 40 J=1,2
DO 40 I=1,4
PTOPC(I,J) = PPPTOP(I,J)
40 CONTINUE
! --- L CLD VERTICAL VEL ADJ BOUNDARIES
VVCLD(1) = 0.0003E0
VVCLD(2) = -0.0005E0
CRHRH = 0.60E0
!--- COMPUTE LLYR--WHICH IS TOPMOST NON CLD(LOW) LAYER, FOR STRATIFORM
XTHK = 0.E0
!.... DEFAULT LLYR
KL = KDIM + 1
!.... TOPMOST NONCLOUD LAYER WILL BE THE ONE AT OR ABOVE LOWEST
! 0.1 OF THE ATMOSPHERE..
!JFM DO 202 K=1,KDIM
DO 202 K=2,KDIM ! k-1 must be greater than 0
! XTHK = XTHK + SI(K) - SI(K+1)
! IF (XTHK.LT.0.1E0) GO TO 202
KL = K
! GO TO 204
IF (SI(K).LT.0.9E0) GO TO 204
202 CONTINUE
204 LLYR = KL-1
!yt PRINT 205,LLYR,KLOWB,KLOWT
205 FORMAT(1H ,'-------LLYR,KLOWB =',2I5)
RETURN
END