--- trunk/phylmd/Radlwsw/swde.f 2013/11/15 18:45:49 76 +++ trunk/phylmd/Radlwsw/swde.f90 2014/03/05 14:38:41 81 @@ -1,133 +1,128 @@ - SUBROUTINE SWDE (PGG,PREF,PRMUZ,PTO1,PW, - S PRE1,PRE2,PTR1,PTR2) - use dimens_m - use dimphy - use raddim - IMPLICIT none -C -C ------------------------------------------------------------------ -C PURPOSE. -C -------- -C COMPUTES THE REFLECTIVITY AND TRANSMISSIVITY OF A CLOUDY -C LAYER USING THE DELTA-EDDINGTON'S APPROXIMATION. -C -C METHOD. -C ------- -C -C STANDARD DELTA-EDDINGTON LAYER CALCULATIONS. -C -C REFERENCE. -C ---------- -C -C SEE RADIATION'S PART OF THE MODEL'S DOCUMENTATION AND -C ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE IFS -C -C AUTHOR. -C ------- -C JEAN-JACQUES MORCRETTE *ECMWF* -C -C MODIFICATIONS. -C -------------- -C ORIGINAL : 88-12-15 -C ------------------------------------------------------------------ -C* ARGUMENTS: -C - DOUBLE PRECISION PGG(KDLON) ! ASSYMETRY FACTOR - DOUBLE PRECISION PREF(KDLON) ! REFLECTIVITY OF THE UNDERLYING LAYER - DOUBLE PRECISION PRMUZ(KDLON) ! COSINE OF SOLAR ZENITH ANGLE - DOUBLE PRECISION PTO1(KDLON) ! OPTICAL THICKNESS - DOUBLE PRECISION PW(KDLON) ! SINGLE SCATTERING ALBEDO - DOUBLE PRECISION PRE1(KDLON) ! LAYER REFLECTIVITY (NO UNDERLYING-LAYER REFLECTION) - DOUBLE PRECISION PRE2(KDLON) ! LAYER REFLECTIVITY - DOUBLE PRECISION PTR1(KDLON) ! LAYER TRANSMISSIVITY (NO UNDERLYING-LAYER REFLECTION) - DOUBLE PRECISION PTR2(KDLON) ! LAYER TRANSMISSIVITY -C -C* LOCAL VARIABLES: -C - INTEGER jl - DOUBLE PRECISION ZFF, ZGP, ZTOP, ZWCP, ZDT, ZX1, ZWM - DOUBLE PRECISION ZRM2, ZRK, ZX2, ZRP, ZALPHA, ZBETA, ZARG - DOUBLE PRECISION ZEXMU0, ZARG2, ZEXKP, ZEXKM, ZXP2P, ZXM2P, ZAP2B - DOUBLE PRECISION ZAM2B - DOUBLE PRECISION ZA11, ZA12, ZA13, ZA21, ZA22, ZA23 - DOUBLE PRECISION ZDENA, ZC1A, ZC2A, ZRI0A, ZRI1A - DOUBLE PRECISION ZRI0B, ZRI1B - DOUBLE PRECISION ZB21, ZB22, ZB23, ZDENB, ZC1B, ZC2B - DOUBLE PRECISION ZRI0C, ZRI1C, ZRI0D, ZRI1D -C ------------------------------------------------------------------ -C -C* 1. DELTA-EDDINGTON CALCULATIONS -C - 100 CONTINUE -C - DO 131 JL = 1, KDLON -C -C* 1.1 SET UP THE DELTA-MODIFIED PARAMETERS -C - 110 CONTINUE -C - ZFF = PGG(JL)*PGG(JL) - ZGP = PGG(JL)/(1.+PGG(JL)) - ZTOP = (1.- PW(JL) * ZFF) * PTO1(JL) - ZWCP = (1-ZFF)* PW(JL) /(1.- PW(JL) * ZFF) - ZDT = 2./3. - ZX1 = 1.-ZWCP*ZGP - ZWM = 1.-ZWCP - ZRM2 = PRMUZ(JL) * PRMUZ(JL) - ZRK = SQRT(3.*ZWM*ZX1) - ZX2 = 4.*(1.-ZRK*ZRK*ZRM2) - ZRP=ZRK/ZX1 - ZALPHA = 3.*ZWCP*ZRM2*(1.+ZGP*ZWM)/ZX2 - ZBETA = 3.*ZWCP* PRMUZ(JL) *(1.+3.*ZGP*ZRM2*ZWM)/ZX2 -CMAF ZARG=MIN(ZTOP/PRMUZ(JL),200.) - ZARG=MIN(ZTOP/PRMUZ(JL),2.0d+2) - ZEXMU0=EXP(-ZARG) -CMAF ZARG2=MIN(ZRK*ZTOP,200.) - ZARG2=MIN(ZRK*ZTOP,2.0d+2) - ZEXKP=EXP(ZARG2) - ZEXKM = 1./ZEXKP - ZXP2P = 1.+ZDT*ZRP - ZXM2P = 1.-ZDT*ZRP - ZAP2B = ZALPHA+ZDT*ZBETA - ZAM2B = ZALPHA-ZDT*ZBETA -C -C* 1.2 WITHOUT REFLECTION FROM THE UNDERLYING LAYER -C - 120 CONTINUE -C - ZA11 = ZXP2P - ZA12 = ZXM2P - ZA13 = ZAP2B - ZA22 = ZXP2P*ZEXKP - ZA21 = ZXM2P*ZEXKM - ZA23 = ZAM2B*ZEXMU0 - ZDENA = ZA11 * ZA22 - ZA21 * ZA12 - ZC1A = (ZA22*ZA13-ZA12*ZA23)/ZDENA - ZC2A = (ZA11*ZA23-ZA21*ZA13)/ZDENA - ZRI0A = ZC1A+ZC2A-ZALPHA - ZRI1A = ZRP*(ZC1A-ZC2A)-ZBETA - PRE1(JL) = (ZRI0A-ZDT*ZRI1A)/ PRMUZ(JL) - ZRI0B = ZC1A*ZEXKM+ZC2A*ZEXKP-ZALPHA*ZEXMU0 - ZRI1B = ZRP*(ZC1A*ZEXKM-ZC2A*ZEXKP)-ZBETA*ZEXMU0 - PTR1(JL) = ZEXMU0+(ZRI0B+ZDT*ZRI1B)/ PRMUZ(JL) -C -C* 1.3 WITH REFLECTION FROM THE UNDERLYING LAYER -C - 130 CONTINUE -C - ZB21 = ZA21- PREF(JL) *ZXP2P*ZEXKM - ZB22 = ZA22- PREF(JL) *ZXM2P*ZEXKP - ZB23 = ZA23- PREF(JL) *ZEXMU0*(ZAP2B - PRMUZ(JL) ) - ZDENB = ZA11 * ZB22 - ZB21 * ZA12 - ZC1B = (ZB22*ZA13-ZA12*ZB23)/ZDENB - ZC2B = (ZA11*ZB23-ZB21*ZA13)/ZDENB - ZRI0C = ZC1B+ZC2B-ZALPHA - ZRI1C = ZRP*(ZC1B-ZC2B)-ZBETA - PRE2(JL) = (ZRI0C-ZDT*ZRI1C) / PRMUZ(JL) - ZRI0D = ZC1B*ZEXKM + ZC2B*ZEXKP - ZALPHA*ZEXMU0 - ZRI1D = ZRP * (ZC1B*ZEXKM - ZC2B*ZEXKP) - ZBETA*ZEXMU0 - PTR2(JL) = ZEXMU0 + (ZRI0D + ZDT*ZRI1D) / PRMUZ(JL) -C - 131 CONTINUE - RETURN - END +SUBROUTINE swde(pgg, pref, prmuz, pto1, pw, pre1, pre2, ptr1, ptr2) + USE dimens_m + USE dimphy + USE raddim + IMPLICIT NONE + + ! ------------------------------------------------------------------ + ! PURPOSE. + ! -------- + ! COMPUTES THE REFLECTIVITY AND TRANSMISSIVITY OF A CLOUDY + ! LAYER USING THE DELTA-EDDINGTON'S APPROXIMATION. + + ! METHOD. + ! ------- + + ! STANDARD DELTA-EDDINGTON LAYER CALCULATIONS. + + ! REFERENCE. + ! ---------- + + ! SEE RADIATION'S PART OF THE MODEL'S DOCUMENTATION AND + ! ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE IFS + + ! AUTHOR. + ! ------- + ! JEAN-JACQUES MORCRETTE *ECMWF* + + ! MODIFICATIONS. + ! -------------- + ! ORIGINAL : 88-12-15 + ! ------------------------------------------------------------------ + ! * ARGUMENTS: + + DOUBLE PRECISION pgg(kdlon) ! ASSYMETRY FACTOR + DOUBLE PRECISION pref(kdlon) ! REFLECTIVITY OF THE UNDERLYING LAYER + DOUBLE PRECISION prmuz(kdlon) ! COSINE OF SOLAR ZENITH ANGLE + DOUBLE PRECISION pto1(kdlon) ! OPTICAL THICKNESS + DOUBLE PRECISION pw(kdlon) ! SINGLE SCATTERING ALBEDO + DOUBLE PRECISION pre1(kdlon) ! LAYER REFLECTIVITY (NO UNDERLYING-LAYER REFLECTION) + DOUBLE PRECISION pre2(kdlon) ! LAYER REFLECTIVITY + DOUBLE PRECISION ptr1(kdlon) ! LAYER TRANSMISSIVITY (NO UNDERLYING-LAYER REFLECTION) + DOUBLE PRECISION ptr2(kdlon) ! LAYER TRANSMISSIVITY + + ! * LOCAL VARIABLES: + + INTEGER jl + DOUBLE PRECISION zff, zgp, ztop, zwcp, zdt, zx1, zwm + DOUBLE PRECISION zrm2, zrk, zx2, zrp, zalpha, zbeta, zarg + DOUBLE PRECISION zexmu0, zarg2, zexkp, zexkm, zxp2p, zxm2p, zap2b + DOUBLE PRECISION zam2b + DOUBLE PRECISION za11, za12, za13, za21, za22, za23 + DOUBLE PRECISION zdena, zc1a, zc2a, zri0a, zri1a + DOUBLE PRECISION zri0b, zri1b + DOUBLE PRECISION zb21, zb22, zb23, zdenb, zc1b, zc2b + DOUBLE PRECISION zri0c, zri1c, zri0d, zri1d + ! ------------------------------------------------------------------ + + ! * 1. DELTA-EDDINGTON CALCULATIONS + + + DO jl = 1, kdlon + + ! * 1.1 SET UP THE DELTA-MODIFIED PARAMETERS + + + zff = pgg(jl)*pgg(jl) + zgp = pgg(jl)/(1.+pgg(jl)) + ztop = (1.-pw(jl)*zff)*pto1(jl) + zwcp = (1-zff)*pw(jl)/(1.-pw(jl)*zff) + zdt = 2./3. + zx1 = 1. - zwcp*zgp + zwm = 1. - zwcp + zrm2 = prmuz(jl)*prmuz(jl) + zrk = sqrt(3.*zwm*zx1) + zx2 = 4.*(1.-zrk*zrk*zrm2) + zrp = zrk/zx1 + zalpha = 3.*zwcp*zrm2*(1.+zgp*zwm)/zx2 + zbeta = 3.*zwcp*prmuz(jl)*(1.+3.*zgp*zrm2*zwm)/zx2 + ! MAF ZARG=MIN(ZTOP/PRMUZ(JL),200.) + zarg = min(ztop/prmuz(jl), 2.0D+2) + zexmu0 = exp(-zarg) + ! MAF ZARG2=MIN(ZRK*ZTOP,200.) + zarg2 = min(zrk*ztop, 2.0D+2) + zexkp = exp(zarg2) + zexkm = 1./zexkp + zxp2p = 1. + zdt*zrp + zxm2p = 1. - zdt*zrp + zap2b = zalpha + zdt*zbeta + zam2b = zalpha - zdt*zbeta + + ! * 1.2 WITHOUT REFLECTION FROM THE UNDERLYING LAYER + + + za11 = zxp2p + za12 = zxm2p + za13 = zap2b + za22 = zxp2p*zexkp + za21 = zxm2p*zexkm + za23 = zam2b*zexmu0 + zdena = za11*za22 - za21*za12 + zc1a = (za22*za13-za12*za23)/zdena + zc2a = (za11*za23-za21*za13)/zdena + zri0a = zc1a + zc2a - zalpha + zri1a = zrp*(zc1a-zc2a) - zbeta + pre1(jl) = (zri0a-zdt*zri1a)/prmuz(jl) + zri0b = zc1a*zexkm + zc2a*zexkp - zalpha*zexmu0 + zri1b = zrp*(zc1a*zexkm-zc2a*zexkp) - zbeta*zexmu0 + ptr1(jl) = zexmu0 + (zri0b+zdt*zri1b)/prmuz(jl) + + ! * 1.3 WITH REFLECTION FROM THE UNDERLYING LAYER + + + zb21 = za21 - pref(jl)*zxp2p*zexkm + zb22 = za22 - pref(jl)*zxm2p*zexkp + zb23 = za23 - pref(jl)*zexmu0*(zap2b-prmuz(jl)) + zdenb = za11*zb22 - zb21*za12 + zc1b = (zb22*za13-za12*zb23)/zdenb + zc2b = (za11*zb23-zb21*za13)/zdenb + zri0c = zc1b + zc2b - zalpha + zri1c = zrp*(zc1b-zc2b) - zbeta + pre2(jl) = (zri0c-zdt*zri1c)/prmuz(jl) + zri0d = zc1b*zexkm + zc2b*zexkp - zalpha*zexmu0 + zri1d = zrp*(zc1b*zexkm-zc2b*zexkp) - zbeta*zexmu0 + ptr2(jl) = zexmu0 + (zri0d+zdt*zri1d)/prmuz(jl) + + END DO + RETURN +END SUBROUTINE swde