*IF DEF,A01_2B SWCLOP2B.2
C ******************************COPYRIGHT****************************** SWCLOP2B.3
C (c) CROWN COPYRIGHT 1997, METEOROLOGICAL OFFICE, All Rights Reserved. SWCLOP2B.4
C SWCLOP2B.5
C Use, duplication or disclosure of this code is subject to the SWCLOP2B.6
C restrictions as set forth in the contract. SWCLOP2B.7
C SWCLOP2B.8
C Meteorological Office SWCLOP2B.9
C London Road SWCLOP2B.10
C BRACKNELL SWCLOP2B.11
C Berkshire UK SWCLOP2B.12
C RG12 2SZ SWCLOP2B.13
C SWCLOP2B.14
C If no contract has been raised with this copy of the code, the use, SWCLOP2B.15
C duplication or disclosure of it is strictly prohibited. Permission SWCLOP2B.16
C to do so must first be obtained in writing from the Head of Numerical SWCLOP2B.17
C Modelling at the above address. SWCLOP2B.18
C ******************************COPYRIGHT****************************** SWCLOP2B.19
C SWCLOP2B.20
CLL Subroutine SWCLOP ---------------------------------------------- SWCLOP2B.21
CLL SWCLOP2B.22
CLL Purpose : SWCLOP2B.23
CLL It calculates cloud optical properties for use by SWMAST. SWCLOP2B.24
CLL It is suitable for single-column use. SWCLOP2B.25
CLL SWCLOP2B.26
CLL Author: William Ingram SWCLOP2B.27
CLL SWCLOP2B.28
CLL Model Modification history: SWCLOP2B.29
CLL version Date SWCLOP2B.30
CLL 4.3 18/4/97 This deck, SWCLOP2B, added to restrict the SWCLOP2B.31
CLL treatment of the singular case when the e-folding depths for direct SWCLOP2B.32
CLL & diffuse light are nearly equal to when there is a danger of SWCLOP2B.33
CLL failure on the T3E - for use in HadCM2, which for true backwards SWCLOP2B.34
CLL compatibility with the original HadCM2 on the C90 would never use SWCLOP2B.35
CLL this form, but which, thanks to the T3E's far lower range of SWCLOP2B.36
CLL numbers, is liable to a failure that never occurred on the C90. SWCLOP2B.37
CLL William Ingram SWCLOP2B.38
CLL SWCLOP2B.39
CLL SWCLOP2B.40
CLL Programming standard : SWCLOP2B.41
CLL Contains ! comments, but otherwise conforms to the FORTRAN 77 SWCLOP2B.42
CLL It conforms to programming standard A of UMDP 4 (version 2, SWCLOP2B.43
CLL 18/1/90). SWCLOP2B.44
CLL standard with no features deprecated by 8X. SWCLOP2B.45
CLL SWCLOP2B.46
CLL Logical components covered : P234 SWCLOP2B.47
CLL (interaction of shortwave radiation with the atmosphere) SWCLOP2B.48
CLL SWCLOP2B.49
CLL Project task : P23 (radiation) SWCLOP2B.50
CLL SWCLOP2B.51
CLL External documentation: UMDP 23, sub-section "cloud optical SWCLOP2B.52
CLL properties" of shortwave section. SWCLOP2B.53
CLL SWCLOP2B.54
CLLEND ----------------------------------------------------------------- SWCLOP2B.55
C*L SWCLOP2B.56
SUBROUTINE SWCLOP (CW, RE, COSZ, 6SWCLOP2B.57
& L1, L2, NCLDS, REF, CTR) SWCLOP2B.58
INTEGER NBANDS ! Number of spectral bands in the shortwave SWCLOP2B.59
PARAMETER (NBANDS=4) ! This run uses the standard set of SWCLOP2B.60
C shortwave bands as described by Slingo (May 1989, J.Atmos.Sci., 46, SWCLOP2B.61
C 10, 1419-1427) or UMDP 23. SWCLOP2B.62
INTEGER!, INTENT (IN) SWCLOP2B.63
& NCLDS, ! Number of clouds SWCLOP2B.64
& L1, ! First dimension of input arrays SWCLOP2B.65
& L2 ! Number of points to be treated SWCLOP2B.66
REAL!, INTENT (IN) SWCLOP2B.67
& CW(L1,NCLDS), ! Cloud condensed water path SWCLOP2B.68
& RE(L1,NCLDS), ! Effective cloud droplet radius SWCLOP2B.69
& COSZ(L1) ! Cos(solar zenith angle) SWCLOP2B.70
REAL!, INTENT (OUT) SWCLOP2B.71
& REF(L2,NBANDS,NCLDS,2), ! Cloud reflectivity and SWCLOP2B.72
& CTR(L2,NBANDS,NCLDS,2) ! transmissivity for direct and SWCLOP2B.73
C ! diffuse radiation respectively SWCLOP2B.74
C SWCLOP2B.75
C ! SWCLOP has no dynamically allocated workspace, no EXTERNAL SWCLOP2B.76
C ! calls and no significant structure - just nested loops. SWCLOP2B.77
C* SWCLOP2B.78
REAL B0CON, U1 ! Used in the quadrature SWCLOP2B.79
PARAMETER (B0CON=.375, U1=2.) SWCLOP2B.80
C ! These values are set for PIFM as defined by Zdunkowski & al SWCLOP2B.81
C ! (1980, Beitraege, 53, 147-166). SWCLOP2B.82
REAL SWCTOL ! Minimum value of D2 for which SWCLOP2B.83
C ! the standard formulae for direct-beam optical properties will SWCLOP2B.84
C ! be used - if D2 is too small, the singular solution for equal SWCLOP2B.85
C ! direct & diffuse e-folding depths is used. SWCLOP2B.86
C ! ("Too small" rather than "zero", now because of the risk of SWCLOP2B.87
C ! failure on dividing by it.) SWCLOP2B.88
C ! Note that this value is intended for the T3E only ! SWCLOP2B.89
C ! 1.7977E+308 is the largest 64-bit REAL the T3E can cope with. SWCLOP2B.90
*IF DEF,T3E PXSWCLOP.1
PARAMETER ( SWCTOL = 1.E-305 ) PXSWCLOP.2
*ELSE PXSWCLOP.3
PARAMETER ( SWCTOL = TINY(1.0) ) PXSWCLOP.4
*ENDIF PXSWCLOP.5
REAL AFIT(NBANDS), BFIT(NBANDS), ! Used in (4.8.23)-(4.8.25) SWCLOP2B.92
& CFIT(NBANDS), DFIT(NBANDS), ! but called just a-f SWCLOP2B.93
& EFIT(NBANDS), FFIT(NBANDS) ! in UMDP 23. SWCLOP2B.94
REAL TAU, OMEGA, G, BETA0, ! Intermediate quantities SWCLOP2B.95
& BETAMU, F, U2, ALPHA1, ! as defined in UMDP 23 SWCLOP2B.96
& ALPHA2, ALPHA3, ALPHA4, EPS,! SWCLOP2B.97
& M, E, H, A1, D1, D2, GAMMA1,! SWCLOP2B.98
& GAMMA2, RT, RM, C1N, P1, P2,! SWCLOP2B.99
& OM1MF, E2 ! OMEGA*(1-F) & E**2 SWCLOP2B.100
INTEGER CLOUD, J, BAND ! Loopers over clouds, points & SWCLOP2B.101
C bands SWCLOP2B.102
C ! These are the values for the standard set of 4 bands. SWCLOP2B.103
C ! The limits on Re they imply are that it must be between SWCLOP2B.104
C ! 0.344 microns (where OMEGA0 in band 1 becomes greater than 1) SWCLOP2B.105
C ! and 37.6 microns (where G in band 3 becomes greater than 1). SWCLOP2B.106
DATA AFIT / 28.17, 26.82, 22.64, 12.81 / SWCLOP2B.107
DATA BFIT / 1.305E-3, 1.346E-3, 1.454E-3, 1.641E-3 / SWCLOP2B.108
DATA CFIT / -56.17E-9, -6.938E-6, 463.6E-6, 0.2011 / SWCLOP2B.109
DATA DFIT / 0.1632, 23.49, 1238., 7556. / SWCLOP2B.110
DATA EFIT / 0.8287, 0.7935, 0.7535, 0.8255 / SWCLOP2B.111
DATA FFIT / 2482., 4226., 6560., 4353. / SWCLOP2B.112
C SWCLOP2B.113
DO 100 CLOUD=1, NCLDS SWCLOP2B.114
DO 100 BAND=1, NBANDS SWCLOP2B.115
Cfpp$ Select(CONCUR) SWCLOP2B.116
DO 100 J=1, L2 SWCLOP2B.117
C ! First calculate basic optical parameters. SWCLOP2B.118
C ! N.B. The code does not check for OMEGA0 < 0 or G > 1 (which SWCLOP2B.119
C ! will occur if Re is too large) or OMEGA0 > 1 (which will SWCLOP2B.120
C ! occur if Re is too small), and these unphysical cases will SWCLOP2B.121
C ! produce a negative argument for the square root. SWCLOP2B.122
IF ( RE(J,CLOUD) .GT. 0. ) THEN ! (4.8.23) & (4.8.24) SWCLOP2B.123
TAU = CW(J,CLOUD) * ( AFIT(BAND) + BFIT(BAND)/RE(J,CLOUD) ) SWCLOP2B.124
OMEGA = 1. - CFIT(BAND) - DFIT(BAND) * RE(J,CLOUD) SWCLOP2B.125
ELSE ! RE may be unset where no cloud SWCLOP2B.126
TAU = 1. SWCLOP2B.127
OMEGA = .5 SWCLOP2B.128
ENDIF SWCLOP2B.129
C ! The special-case code ends up with 0/0 if tau v.small SWCLOP2B.130
C ! N.B. 32bit machines may need larger small value SWCLOP2B.131
IF ( TAU .LE. 1.0E-10 ) THEN ! SWCLOP2B.132
CTR(J,BAND,CLOUD,1) = 1. SWCLOP2B.133
CTR(J,BAND,CLOUD,2) = 1. SWCLOP2B.134
REF(J,BAND,CLOUD,1) = 0. SWCLOP2B.135
REF(J,BAND,CLOUD,2) = 0. SWCLOP2B.136
ELSE SWCLOP2B.137
G = EFIT(BAND) + FFIT(BAND) * RE(J,CLOUD) ! (4.8.25) SWCLOP2B.138
C ! Now apply the chosen quadrature - Zdunkovski's PIFM SWCLOP2B.139
BETA0 = B0CON * (1.-G) ! (4.8.18) SWCLOP2B.140
BETAMU = 0.5 - 0.75 * COSZ(J) * G/(1.+G) ! (4.8.19) SWCLOP2B.141
F = G * G ! (4.8.20) SWCLOP2B.142
U2 = 2. ! (4.8.22) SWCLOP2B.143
C ! Now have quadrature-dept quantities. (U1 is constant) SWCLOP2B.144
C ! Next find alphas etc: SWCLOP2B.145
ALPHA1 = U1 * ( 1. - OMEGA*(1.-BETA0) ) ! (4.8.14) SWCLOP2B.146
ALPHA2 = U2 * OMEGA * BETA0 ! (4.8.15) SWCLOP2B.147
OM1MF = OMEGA * (1.-F) SWCLOP2B.148
ALPHA3 = OM1MF * BETAMU ! (4.8.16) SWCLOP2B.149
ALPHA4 = OM1MF - ALPHA3 ! (4.8.17) SWCLOP2B.150
EPS = SQRT ( ALPHA1**2 - ALPHA2**2 ) ! (4.8.13) SWCLOP2B.151
M = ALPHA2 / (ALPHA1+EPS) ! (4.8.8) SWCLOP2B.152
E = EXP ( - EPS * TAU ) ! (4.8.7) SWCLOP2B.153
E2 = E * E SWCLOP2B.154
H = 1. - OMEGA * F ! (4.8.12) SWCLOP2B.155
C ! Now can find a1, the contribution to the transmissivity for a SWCLOP2B.156
C ! direct incident beam due to transmission as a direct beam. SWCLOP2B.157
A1 = EXP ( - H * TAU / COSZ(J) ) ! (4.8.1) SWCLOP2B.158
C ! Now find transmissivities & reflectivities for diffuse SWCLOP2B.159
C ! incident beam, a4 & a5: SWCLOP2B.160
D1 = 1. - E2 * M * M ! (4.8.6) SWCLOP2B.161
CTR(J,BAND,CLOUD,2) = E * (1.-M*M) / D1 ! (4.8.4) SWCLOP2B.162
REF(J,BAND,CLOUD,2) = M * (1.-E2) / D1 ! (4.8.5) SWCLOP2B.163
C ! Next find the gammas: SWCLOP2B.164
D2 = ( COSZ(J)*EPS )**2 - H*H ! (4.8.11) SWCLOP2B.165
P1 = COSZ(J) * ( ALPHA1*ALPHA3 + ALPHA2*ALPHA4 ) - H*ALPHA3 SWCLOP2B.166
P2 = COSZ(J) * ( ALPHA1*ALPHA4 + ALPHA2*ALPHA3 ) + H*ALPHA4 SWCLOP2B.167
CTR(J,BAND,CLOUD,1) = .5 ! Ensure valid REALs set for 2nd IF SWCLOP2B.168
REF(J,BAND,CLOUD,1) = .5 ! block to test, even if 1st skipped. SWCLOP2B.169
IF ( ABS(D2) .GT. SWCTOL ) THEN SWCLOP2B.170
C ! For direct incident beam, (a1+a2) & a3 (eqns 4.8.2-4.8.5): SWCLOP2B.171
GAMMA1 = P1 / D2 SWCLOP2B.172
GAMMA2 = P2 / D2 SWCLOP2B.173
CTR(J,BAND,CLOUD,1) = A1 * ( 1.-REF(J,BAND,CLOUD,2)*GAMMA1 ) SWCLOP2B.174
& - GAMMA2 * ( CTR(J,BAND,CLOUD,2) - A1 ) SWCLOP2B.175
REF(J,BAND,CLOUD,1) = GAMMA1 * ( 1.-CTR(J,BAND,CLOUD,2)*A1 ) SWCLOP2B.176
& - GAMMA2 * REF(J,BAND,CLOUD,2) SWCLOP2B.177
ELSE SWCLOP2B.178
GAMMA1 = .5 * P1 / EPS SWCLOP2B.179
GAMMA2 = .5 * P2 / EPS SWCLOP2B.180
RT = ALPHA2 * TAU / ( 1. + ALPHA3 / GAMMA1 ) SWCLOP2B.181
RM = M * ( 1. + M * RT ) / ( RT + M ) SWCLOP2B.182
C1N = GAMMA1 * TAU / ( RM - 1. ) SWCLOP2B.183
CTR(J,BAND,CLOUD,1) = A1 + SWCLOP2B.184
& E * ( C1N * M * (1.-E2) + TAU * GAMMA2 ) / COSZ(J) SWCLOP2B.185
REF(J,BAND,CLOUD,1) = C1N * ( E2 - RM ) / COSZ(J) SWCLOP2B.186
ENDIF SWCLOP2B.187
ENDIF SWCLOP2B.188
100 CONTINUE SWCLOP2B.189
RETURN SWCLOP2B.190
END SWCLOP2B.191
*ENDIF SWCLOP2B.192