GCC Code Coverage Report

Directory: ./
File: rad/srtm_reftra.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 62 0.0%
Branches: 0 14 0.0%
Line Branch Exec Source
1 ✗ SUBROUTINE SRTM_REFTRA &
2 & ( KLEV , KMODTS, &
3 & LDRTCHK, &
4 & PGG , PRMUZ, PTAU , PW, &
5 & PREF , PREFD, PTRA , PTRAD &
6 & )
7
8 !**** *SRTM_REFTRA* - REFLECTIVITY AND TRANSMISSIVITY
9
10 ! PURPOSE.
11 ! --------
12 ! COMPUTES THE REFLECTIVITY AND TRANSMISSIVITY OF A CLEAR OR
13 ! CLOUDY LAYER USING A CHOICE OF VARIOUS APPROXIMATIONS.
14
15 !** INTERFACE.
16 ! ----------
17 ! *SRTM_REFTRA* IS CALLED BY *SRTM_SPCVRT*
18
19 ! EXPLICIT ARGUMENTS :
20 ! --------------------
21 ! INPUTS
22 ! ------
23 ! KMODTS = 1 EDDINGTON (JOSEPH ET AL., 1976)
24 ! = 2 PIFM (ZDUNKOWSKI ET AL., 1980)
25 ! = 3 DISCRETE ORDINATES (LIOU, 1973)
26 ! LDRTCHK = .T. IF CLOUDY
27 ! = .F. IF CLEAR-SKY
28 ! PGG = ASSYMETRY FACTOR
29 ! PRMUZ = COSINE SOLAR ZENITH ANGLE
30 ! PTAU = OPTICAL THICKNESS
31 ! PW = SINGLE SCATTERING ALBEDO
32
33 ! OUTPUTS
34 ! -------
35 ! PREF : COLLIMATED BEAM REFLECTIVITY
36 ! PREFD : DIFFUSE BEAM REFLECTIVITY
37 ! PTRA : COLLIMATED BEAM TRANSMISSIVITY
38 ! PTRAD : DIFFUSE BEAM TRANSMISSIVITY
39
40 ! METHOD.
41 ! -------
42 ! STANDARD DELTA-EDDINGTON, P.I.F.M., OR D.O.M. LAYER CALCULATIONS.
43
44 ! EXTERNALS.
45 ! ----------
46 ! NONE
47
48 ! REFERENCE.
49 ! ----------
50
51 ! AUTHOR.
52 ! -------
53 ! JEAN-JACQUES MORCRETTE *ECMWF*
54
55 ! MODIFICATIONS.
56 ! --------------
57 ! ORIGINAL : 03-02-27
58 ! M.Hamrud 01-Oct-2003 CY28 Cleaning
59 ! Mike Iacono, AER, Mar 2004: bug fix
60
61 ! ------------------------------------------------------------------
62
63 !* 0.1 ARGUMENTS
64 ! ---------
65
66 USE PARKIND1 ,ONLY : JPIM ,JPRB
67 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
68
69 USE PARSRTM , ONLY : JPLAY
70 !USE YOESW , ONLY : NDBUG
71 USE YOERDU , ONLY : REPLOG
72
73 IMPLICIT NONE
74
75 INTEGER(KIND=JPIM),INTENT(IN) :: KLEV
76 INTEGER(KIND=JPIM),INTENT(OUT) :: KMODTS
77 LOGICAL ,INTENT(IN) :: LDRTCHK(JPLAY)
78 REAL(KIND=JPRB) ,INTENT(IN) :: PGG(JPLAY)
79 REAL(KIND=JPRB) ,INTENT(IN) :: PRMUZ
80 REAL(KIND=JPRB) ,INTENT(IN) :: PTAU(JPLAY)
81 REAL(KIND=JPRB) ,INTENT(IN) :: PW(JPLAY)
82 REAL(KIND=JPRB) ,INTENT(INOUT) :: PREF(JPLAY)
83 REAL(KIND=JPRB) ,INTENT(INOUT) :: PREFD(JPLAY)
84 REAL(KIND=JPRB) ,INTENT(INOUT) :: PTRA(JPLAY)
85 REAL(KIND=JPRB) ,INTENT(INOUT) :: PTRAD(JPLAY)
86 ! ------------------------------------------------------------------
87
88 INTEGER(KIND=JPIM) :: JK
89
90 REAL(KIND=JPRB) :: ZA, ZA1, ZA2
91 REAL(KIND=JPRB) :: ZBETA, ZDEND, ZDENR, ZDENT
92 REAL(KIND=JPRB) :: ZE1, ZE2, ZEM1, ZEM2, ZEMM, ZEP1, ZEP2
93 REAL(KIND=JPRB) :: ZG, ZG3, ZGAMMA1, ZGAMMA2, ZGAMMA3, ZGAMMA4, ZGT
94 REAL(KIND=JPRB) :: ZR1, ZR2, ZR3, ZR4, ZR5, ZRK, ZRK2, ZRKG, ZRM1, ZRP, ZRP1, ZRPP
95 REAL(KIND=JPRB) :: ZSR3, ZT1, ZT2, ZT3, ZT4, ZT5, ZTO1
96 REAL(KIND=JPRB) :: ZW, ZWCRIT, ZWO
97 REAL(KIND=JPRB) :: ZHOOK_HANDLE,EXP500,ZTEMP
98
99 ! ------------------------------------------------------------------
100
101 ✗ IF (LHOOK) CALL DR_HOOK('SRTM_REFTRA',0,ZHOOK_HANDLE)
102 EXP500=EXP(500.0_JPRB)
103 ZSR3=SQRT(3._JPRB)
104 ZWCRIT=0.9995_JPRB
105 ✗ KMODTS=2
106
107 !NDBUG=3
108
109 ✗ DO JK=1,KLEV
110 ! if (NDBUG < 2) then
111 ! print 9000,JK,LRTCHK(JK),PTAU(JK),PW(JK),PGG(JK),PRMUZ
112 9000 format(1x,'SRTM_REFTRA:inputs:',I3,L8,4E13.6)
113 ! end if
114 ✗ IF (.NOT.LDRTCHK(JK)) THEN
115 ✗ PREF(JK) =0.0_JPRB
116 ✗ PTRA(JK) =1.0_JPRB
117 ✗ PREFD(JK)=0.0_JPRB
118 ✗ PTRAD(JK)=1.0_JPRB
119 ! if (NDBUG < 2) then
120 ! print 9001,JL,JK,PREF(JK),PTRA(JK),PREFD(JK),PTRAD(JK)
121 9001 format(1x,'SRTM_REFTRA:not.LRTCKH:',2I3,4F10.6)
122 ! end if
123 ELSE
124 ✗ ZTO1=PTAU(JK)
125 ✗ ZW =PW(JK)
126 ✗ ZG =PGG(JK)
127
128 !-- GENERAL TWO-STREAM EXPRESSIONS
129
130 ✗ ZG3= 3._JPRB * ZG
131 IF (KMODTS == 1) THEN
132 ZGAMMA1= (7._JPRB - ZW * (4._JPRB + ZG3)) * 0.25_JPRB
133 ZGAMMA2=-(1._JPRB - ZW * (4._JPRB - ZG3)) * 0.25_JPRB
134 ZGAMMA3= (2._JPRB - ZG3 * PRMUZ ) * 0.25_JPRB
135 ELSEIF (KMODTS == 2) THEN
136 ✗ ZGAMMA1= (8._JPRB - ZW * (5._JPRB + ZG3)) * 0.25_JPRB
137 ✗ ZGAMMA2= 3._JPRB *(ZW * (1._JPRB - ZG )) * 0.25_JPRB
138 ✗ ZGAMMA3= (2._JPRB - ZG3 * PRMUZ ) * 0.25_JPRB
139 ELSEIF (KMODTS == 3) THEN
140 ZGAMMA1= ZSR3 * (2._JPRB - ZW * (1._JPRB + ZG)) * 0.5_JPRB
141 ZGAMMA2= ZSR3 * ZW * (1._JPRB - ZG ) * 0.5_JPRB
142 ZGAMMA3= (1._JPRB - ZSR3 * ZG * PRMUZ ) * 0.5_JPRB
143 ENDIF
144 ✗ ZGAMMA4= 1._JPRB - ZGAMMA3
145
146 !-- RECOMPUTE ORIGINAL S.S.A. TO TEST FOR CONSERVATIVE SOLUTION
147 ✗ ZWO= ZW / (1._JPRB - (1._JPRB - ZW) * (ZG / (1._JPRB - ZG))**2)
148 ! ZTEMP=(1._JPRB - ZG)**2
149 ! ZWO= ZW*ZTEMP/ (ZTEMP - (1._JPRB - ZW)*(ZG **2))
150
151 ✗ IF (ZWO >= ZWCRIT) THEN
152 !!-- conservative scattering
153
154 ✗ ZA = ZGAMMA1 * PRMUZ
155 ✗ ZA1 = ZA - ZGAMMA3
156 ✗ ZGT = ZGAMMA1 * ZTO1
157
158 !-- Homogeneous reflectance and transmittance
159
160 ! collimated beam
161
162 ✗ ZE1 = MIN ( ZTO1 / PRMUZ , 500._JPRB)
163 ✗ ZE2 = EXP ( - ZE1 )
164 ✗ ZTEMP=1.0_JPRB/(1._JPRB + ZGT)
165 ✗ PREF(JK) = (ZGT - ZA1 * (1._JPRB - ZE2)) *ZTEMP
166 ✗ PTRA(JK) = 1._JPRB - PREF(JK)
167
168 ! isotropic incidence
169
170 ✗ PREFD(JK) = ZGT *ZTEMP
171 ✗ PTRAD(JK) = 1._JPRB - PREFD(JK)
172
173 ! if (NDBUG < 2) then
174 ! print 9002,JL,JK,PREF(JK),PTRA(JK),PREFD(JK),PTRAD(JK)
175 9002 format(1x,'SRTM_REFTRA: consrv: LDRTCHK:',2I3,4F10.6)
176 ! end if
177
178 ELSE
179
180 !-- non-conservative scattering
181
182 ✗ ZA1 = ZGAMMA1 * ZGAMMA4 + ZGAMMA2 * ZGAMMA3
183 ✗ ZA2 = ZGAMMA1 * ZGAMMA3 + ZGAMMA2 * ZGAMMA4
184 ! ZRK = SQRT ( ZGAMMA1**2 - ZGAMMA2**2)
185 ✗ ZRK = SQRT ( MAX ( REPLOG, ZGAMMA1**2 - ZGAMMA2**2) )
186 ✗ ZRP = ZRK * PRMUZ
187 ✗ ZRP1 = 1._JPRB + ZRP
188 ✗ ZRM1 = 1._JPRB - ZRP
189 ✗ ZRK2 = 2._JPRB * ZRK
190 ✗ ZRPP = 1._JPRB - ZRP*ZRP
191 ✗ ZRKG = ZRK + ZGAMMA1
192 ✗ ZR1 = ZRM1 * (ZA2 + ZRK * ZGAMMA3)
193 ✗ ZR2 = ZRP1 * (ZA2 - ZRK * ZGAMMA3)
194 ✗ ZR3 = ZRK2 * (ZGAMMA3 - ZA2 * PRMUZ )
195 ✗ ZR4 = ZRPP * ZRKG
196 ✗ ZR5 = ZRPP * (ZRK - ZGAMMA1)
197 ✗ ZT1 = ZRP1 * (ZA1 + ZRK * ZGAMMA4)
198 ✗ ZT2 = ZRM1 * (ZA1 - ZRK * ZGAMMA4)
199 ✗ ZT3 = ZRK2 * (ZGAMMA4 + ZA1 * PRMUZ )
200 ZT4 = ZR4
201 ZT5 = ZR5
202 ✗ ZBETA = - ZR5 / ZR4
203
204 !-- Homogeneous reflectance and transmittance
205
206 ✗ IF(ZRK * ZTO1 > 500._JPRB)THEN
207 ZEP1=EXP500
208 ELSE
209 ✗ ZEP1=EXP(ZRK * ZTO1)
210 ENDIF
211 ✗ ZEM1=1.0_JPRB/ZEP1
212 ✗ IF(ZTO1 > 500._JPRB*PRMUZ)THEN
213 ZEP2=EXP500
214 ELSE
215 ✗ ZEP2=EXP(ZTO1 / PRMUZ)
216 ENDIF
217 ✗ ZEM2=1.0_JPRB/ZEP2
218
219 ! ZE1 = MIN ( ZRK * ZTO1, 500._JPRB)
220 ! ZE2 = MIN ( ZTO1 / PRMUZ , 500._JPRB)
221
222 ! ZEP1 = EXP( ZE1 )
223 ! ZEM1 = EXP(-ZE1 )
224 ! ZEM1=1.0_JPRB/ZEP1
225
226 ! ZEP2 = EXP( ZE2 )
227 ! ZEM2 = EXP(-ZE2 )
228 ! ZEM2=1.0_JPRB/ZEP2
229
230 ! collimated beam
231
232 ✗ ZDENR = ZR4*ZEP1 + ZR5*ZEM1
233 !- bug noticed by Mike Iacono
234 ! PREF(JK) = ZWO * (ZR1*ZEP1 - ZR2*ZEM1 - ZR3*ZEM2) / ZDENR
235 ✗ PREF(JK) = ZW * (ZR1*ZEP1 - ZR2*ZEM1 - ZR3*ZEM2) / ZDENR
236
237 ZDENT = ZT4*ZEP1 + ZT5*ZEM1
238 !- bug noticed by Mike Iacono
239 ! PTRA(JK) = ZEM2 * (1._JPRB - ZWO * (ZT1*ZEP1 - ZT2*ZEM1 - ZT3*ZEP2) / ZDENT)
240 ✗ PTRA(JK) = ZEM2 * (1._JPRB - ZW * (ZT1*ZEP1 - ZT2*ZEM1 - ZT3*ZEP2) / ZDENT)
241
242 ! diffuse beam
243
244 ✗ ZEMM = ZEM1*ZEM1
245 ✗ ZDEND = 1._JPRB / ( (1._JPRB - ZBETA*ZEMM ) * ZRKG)
246 ✗ PREFD(JK) = ZGAMMA2 * (1._JPRB - ZEMM) * ZDEND
247 ✗ PTRAD(JK) = ZRK2*ZEM1*ZDEND
248
249 ! if (NDBUG < 2) then
250 ! print 9003,JL,JK,PREF(JK),PTRA(JK),PREFD(JK),PTRAD(JK)
251 9003 format(1x,'SRTM_REFTRA: OMG<1: LDRTCHK:',2I3,4F10.6)
252 ! end if
253 ENDIF
254
255 ENDIF
256
257 ENDDO
258
259 ! ------------------------------------------------------------------
260 ✗ IF (LHOOK) CALL DR_HOOK('SRTM_REFTRA',1,ZHOOK_HANDLE)
261 ✗ END SUBROUTINE SRTM_REFTRA
262