*IF DEF,A01_1A,OR,DEF,A01_1B,OR,DEF,A01_2A AWI1F403.428
C ******************************COPYRIGHT****************************** GTS2F400.9937
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.9938
C GTS2F400.9939
C Use, duplication or disclosure of this code is subject to the GTS2F400.9940
C restrictions as set forth in the contract. GTS2F400.9941
C GTS2F400.9942
C Meteorological Office GTS2F400.9943
C London Road GTS2F400.9944
C BRACKNELL GTS2F400.9945
C Berkshire UK GTS2F400.9946
C RG12 2SZ GTS2F400.9947
C GTS2F400.9948
C If no contract has been raised with this copy of the code, the use, GTS2F400.9949
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.9950
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.9951
C Modelling at the above address. GTS2F400.9952
C ******************************COPYRIGHT****************************** GTS2F400.9953
C GTS2F400.9954
CLL Subroutine SWCLOP ---------------------------------------------- SWCLOP1A.3
CLL SWCLOP1A.4
CLL Purpose : SWCLOP1A.5
CLL It calculates cloud optical properties for use by SWMAST. SWCLOP1A.6
CLL It is suitable for single-column use. SWCLOP1A.7
CLL SWCLOP1A.8
CLL Author: William Ingram SWCLOP1A.9
CLL SWCLOP1A.10
CLL Model Modification history from model version 3.0: SWCLOP1A.11
CLL version Date SWCLOP1A.12
CLL 4.0 4/7/95 Code added for safety & accuracy near the singular AWI3F400.1
CLL case when the e-folding depths for direct & diffuse light are equal AWI3F400.2
CLL 4.2 22/10/96 The code for the singular case when the e-folding ADR1F402.1
CLL depths for direct & diffuse light are nearly equal ADR1F402.2
CLL rearranged to a form algebraically the same but ADR1F402.3
CLL avoiding numerical problems for very thick cloud ADR1F402.4
CLL which gave failure on the T3E (also slightly more ADR1F402.5
CLL efficient). W Ingram ADR1F402.6
CLL ADR1F402.7
CLL SWCLOP1A.16
CLL Programming standard : SWCLOP1A.17
CLL Contains ! comments, but otherwise conforms to the FORTRAN 77 SWCLOP1A.18
CLL It conforms to programming standard A of UMDP 4 (version 2, SWCLOP1A.19
CLL 18/1/90). SWCLOP1A.20
CLL standard with no features deprecated by 8X. SWCLOP1A.21
CLL SWCLOP1A.22
CLL Logical components covered : P234 SWCLOP1A.23
CLL (interaction of shortwave radiation with the atmosphere) SWCLOP1A.24
CLL SWCLOP1A.25
CLL Project task : P23 (radiation) SWCLOP1A.26
CLL SWCLOP1A.27
CLL External documentation: UMDP 23, sub-section "cloud optical SWCLOP1A.28
CLL properties" of shortwave section. SWCLOP1A.29
CLL SWCLOP1A.30
CLLEND ----------------------------------------------------------------- SWCLOP1A.31
C*L SWCLOP1A.32
SUBROUTINE SWCLOP (CW, RE, COSZ, 6SWCLOP1A.33
& L1, L2, NCLDS, REF, CTR) SWCLOP1A.34
*CALL SWNBANDS 
SWCLOP1A.35
INTEGER!, INTENT (IN) SWCLOP1A.36
& NCLDS, ! Number of clouds SWCLOP1A.37
& L1, ! First dimension of input arrays SWCLOP1A.38
& L2 ! Number of points to be treated SWCLOP1A.39
REAL!, INTENT (IN) SWCLOP1A.40
& CW(L1,NCLDS), ! Cloud condensed water path SWCLOP1A.41
& RE(L1,NCLDS), ! Effective cloud droplet radius SWCLOP1A.42
& COSZ(L1) ! Cos(solar zenith angle) SWCLOP1A.43
REAL!, INTENT (OUT) SWCLOP1A.44
& REF(L2,NBANDS,NCLDS,2), ! Cloud reflectivity and SWCLOP1A.45
& CTR(L2,NBANDS,NCLDS,2) ! transmissivity for direct and SWCLOP1A.46
C ! diffuse radiation respectively SWCLOP1A.47
C SWCLOP1A.48
C ! SWCLOP has no dynamically allocated workspace, no EXTERNAL SWCLOP1A.49
C ! calls and no significant structure - just nested loops. SWCLOP1A.50
C* SWCLOP1A.51
REAL B0CON, U1 ! Used in the quadrature SWCLOP1A.52
PARAMETER (B0CON=.375, U1=2.) SWCLOP1A.53
C ! These values are set for PIFM as defined by Zdunkowski & al SWCLOP1A.54
C ! (1980, Beitraege, 53, 147-166). SWCLOP1A.55
REAL SWCTOL ! Minimum value of D2 for which AWI3F400.3
C ! the standard formulae for direct-beam optical properties will AWI3F400.4
C ! be used - if D2 is too small, the singular solution for equal AWI3F400.5
C ! direct & diffuse e-folding depths is used. AWI3F400.6
C ! ("Too small" rather than "zero", because near the singular case AWI3F400.7
C ! the general formulae become ill-conditioned and only the AWI3F400.8
C ! asymptotic ones give sensible answers.) AWI3F400.9
PARAMETER ( SWCTOL = .005 ) AWI3F400.10
C ! Surprisingly, the errors are not machine-dependent. AWI3F400.11
REAL AFIT(NBANDS), BFIT(NBANDS), ! Used in (4.8.23)-(4.8.25) SWCLOP1A.56
& CFIT(NBANDS), DFIT(NBANDS), ! but called just a-f SWCLOP1A.57
& EFIT(NBANDS), FFIT(NBANDS) ! in UMDP 23. SWCLOP1A.58
REAL TAU, OMEGA, G, BETA0, ! Intermediate quantities AWI3F400.12
& BETAMU, F, U2, ALPHA1, ! as defined in UMDP 23 AWI3F400.13
& ALPHA2, ALPHA3, ALPHA4, EPS,! AWI3F400.14
& M, E, H, A1, D1, D2, GAMMA1,! AWI3F400.15
& GAMMA2, RT, RM, C1N, P1, P2,! AWI3F400.16
& OM1MF, E2 ! OMEGA*(1-F) & E**2 AWI3F400.17
INTEGER CLOUD, J, BAND ! Loopers over clouds, points & SWCLOP1A.65
C bands SWCLOP1A.66
*CALL SWCLOPNU SWCLOP1A.67
C SWCLOP1A.68
DO 100 CLOUD=1, NCLDS SWCLOP1A.69
DO 100 BAND=1, NBANDS SWCLOP1A.70
Cfpp$ Select(CONCUR) SWCLOP1A.71
DO 100 J=1, L2 SWCLOP1A.72
C ! First calculate basic optical parameters. SWCLOP1A.73
C ! N.B. The code does not check for OMEGA0 < 0 or G > 1 (which SWCLOP1A.74
C ! will occur if Re is too large) or OMEGA0 > 1 (which will SWCLOP1A.75
C ! occur if Re is too small), and these unphysical cases will SWCLOP1A.76
C ! produce a negative argument for the square root. SWCLOP1A.77
IF ( RE(J,CLOUD) .GT. 0. ) THEN ! (4.8.23) & (4.8.24) SWCLOP1A.78
TAU = CW(J,CLOUD) * ( AFIT(BAND) + BFIT(BAND)/RE(J,CLOUD) ) SWCLOP1A.79
OMEGA = 1. - CFIT(BAND) - DFIT(BAND) * RE(J,CLOUD) SWCLOP1A.80
ELSE ! RE may be unset where no cloud SWCLOP1A.81
TAU = 1. SWCLOP1A.82
OMEGA = .5 SWCLOP1A.83
ENDIF SWCLOP1A.84
C ! The special-case code ends up with 0/0 if tau v.small AWI3F400.18
C ! N.B. 32bit machines may need larger small value AWI3F400.19
*IF DEF,T3E,AND,-DEF,SCMA AJC0F405.203
IF ( TAU .LE. 1.0E-10 ) THEN ! AWI3F400.20
*ELSE AJC0F405.204
IF ( TAU .LE. 1.0E-4 ) THEN ! AJC0F405.205
*ENDIF AJC0F405.206
CTR(J,BAND,CLOUD,1) = 1. AWI3F400.21
CTR(J,BAND,CLOUD,2) = 1. AWI3F400.22
REF(J,BAND,CLOUD,1) = 0. AWI3F400.23
REF(J,BAND,CLOUD,2) = 0. AWI3F400.24
ELSE AWI3F400.25
G = EFIT(BAND) + FFIT(BAND) * RE(J,CLOUD) ! (4.8.25) SWCLOP1A.85
C ! Now apply the chosen quadrature - Zdunkovski's PIFM SWCLOP1A.86
BETA0 = B0CON * (1.-G) ! (4.8.18) SWCLOP1A.87
BETAMU = 0.5 - 0.75 * COSZ(J) * G/(1.+G) ! (4.8.19) SWCLOP1A.88
F = G * G ! (4.8.20) SWCLOP1A.89
U2 = 2. ! (4.8.22) SWCLOP1A.90
C ! Now have quadrature-dept quantities. (U1 is constant) SWCLOP1A.91
C ! Next find alphas etc: SWCLOP1A.92
ALPHA1 = U1 * ( 1. - OMEGA*(1.-BETA0) ) ! (4.8.14) SWCLOP1A.93
ALPHA2 = U2 * OMEGA * BETA0 ! (4.8.15) SWCLOP1A.94
OM1MF = OMEGA * (1.-F) SWCLOP1A.95
ALPHA3 = OM1MF * BETAMU ! (4.8.16) SWCLOP1A.96
ALPHA4 = OM1MF - ALPHA3 ! (4.8.17) SWCLOP1A.97
EPS = SQRT ( ALPHA1**2 - ALPHA2**2 ) ! (4.8.13) SWCLOP1A.98
M = ALPHA2 / (ALPHA1+EPS) ! (4.8.8) SWCLOP1A.99
E = EXP ( - EPS * TAU ) ! (4.8.7) SWCLOP1A.100
E2 = E * E AWI3F400.26
H = 1. - OMEGA * F ! (4.8.12) SWCLOP1A.102
C ! Now can find a1, the contribution to the transmissivity for a SWCLOP1A.103
C ! direct incident beam due to transmission as a direct beam. SWCLOP1A.104
A1 = EXP ( - H * TAU / COSZ(J) ) ! (4.8.1) SWCLOP1A.105
C ! Now find transmissivities & reflectivities for diffuse SWCLOP1A.106
C ! incident beam, a4 & a5: SWCLOP1A.107
D1 = 1. - E2 * M * M ! (4.8.6) AWI3F400.27
CTR(J,BAND,CLOUD,2) = E * (1.-M*M) / D1 ! (4.8.4) SWCLOP1A.109
REF(J,BAND,CLOUD,2) = M * (1.-E2) / D1 ! (4.8.5) AWI3F400.28
C ! Next find the gammas: SWCLOP1A.111
D2 = ( COSZ(J)*EPS )**2 - H*H ! (4.8.11) SWCLOP1A.112
P1 = COSZ(J) * ( ALPHA1*ALPHA3 + ALPHA2*ALPHA4 ) - H*ALPHA3 AWI3F400.29
P2 = COSZ(J) * ( ALPHA1*ALPHA4 + ALPHA2*ALPHA3 ) + H*ALPHA4 AWI3F400.30
CTR(J,BAND,CLOUD,1) = .5 ! Ensure valid REALs set for 2nd IF AWI3F400.31
REF(J,BAND,CLOUD,1) = .5 ! block to test, even if 1st skipped. AWI3F400.32
IF ( D2 .NE. 0. ) THEN AWI3F400.33
C ! For direct incident beam, (a1+a2) & a3 (eqns 4.8.2-4.8.5): AWI3F400.34
GAMMA1 = P1 / D2 AWI3F400.35
GAMMA2 = P2 / D2 AWI3F400.36
CTR(J,BAND,CLOUD,1) = A1 * ( 1.-REF(J,BAND,CLOUD,2)*GAMMA1 ) AWI3F400.37
& - GAMMA2 * ( CTR(J,BAND,CLOUD,2) - A1 ) AWI3F400.38
REF(J,BAND,CLOUD,1) = GAMMA1 * ( 1.-CTR(J,BAND,CLOUD,2)*A1 ) AWI3F400.39
& - GAMMA2 * REF(J,BAND,CLOUD,2) AWI3F400.40
ENDIF AWI3F400.41
IF ( ABS(D2) .LT. SWCTOL .OR. AWI3F400.42
& CTR(J,BAND,CLOUD,1) .GT. 1. .OR. CTR(J,BAND,CLOUD,1) .LT. 0. .OR. AWI3F400.43
& REF(J,BAND,CLOUD,1) .GT. 1. .OR. REF(J,BAND,CLOUD,1) .LT. 0.)THEN AWI3F400.44
C ! With the standard 4 bands and rE, this can only be AWI3F400.45
C ! executed in band 4 (only a 30th of the total insolation.) AWI3F400.46
GAMMA1 = .5 * P1 / EPS AWI3F400.47
GAMMA2 = .5 * P2 / EPS AWI3F400.48
RT = ALPHA2 * TAU / ( 1. + ALPHA3 / GAMMA1 ) ADR1F402.8
RM = M * ( 1. + M * RT ) / ( RT + M ) ADR1F402.9
C1N = GAMMA1 * TAU / ( RM - 1. ) ADR1F402.10
CTR(J,BAND,CLOUD,1) = A1 + ADR1F402.11
& E * ( C1N * M * (1.-E2) + TAU * GAMMA2 ) / COSZ(J) ADR1F402.12
REF(J,BAND,CLOUD,1) = C1N * ( E2 - RM ) / COSZ(J) ADR1F402.13
ENDIF AWI3F400.55
ENDIF AWI3F400.56
100 CONTINUE SWCLOP1A.131
RETURN SWCLOP1A.132
END SWCLOP1A.133
*ENDIF DEF,A01_1A,OR,DEF,A01_1B,OR,DEF,A01_2A SWCLOP1A.134