GCC Code Coverage Report


Directory: ./
File: rad/srtm_cldprop.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 96 0.0%
Branches: 0 76 0.0%

Line Branch Exec Source
1 SUBROUTINE SRTM_CLDPROP &
2 & ( KLEV, K_ICLDATM, K_INFLAG, K_ICEFLAG, K_LIQFLAG, K_NSTR,&
3 & P_CLDFRAC, P_CLDDAT1, P_CLDDAT2, P_CLDDAT3, P_CLDDAT4, P_CLDDATMOM,&
4 & P_TAUCLDORIG, P_TAUCLOUD, P_SSACLOUD, P_XMOM &
5 & )
6
7 ! path: $Source: /storm/rc1/cvsroot/rc/rrtm_sw/src/cldprop_sw.f,v $
8 ! author: $Author: jdelamer $
9 ! revision: $Revision: 2.6 $
10 ! created: $Date: 2002/04/04 18:29:47 $
11
12 ! PURPOSE: COMPUTE THE CLOUD OPTICAL DEPTH(S) FOR EACH CLOUDY
13 ! LAYER. NOTE: ONLY INFLAG = 0 AND INFLAG=2/LIQFLAG=1,
14 ! ICEFLAG=3
15 ! (HU & STAMNES, Q. FU) ARE IMPLEMENTED.
16
17 USE PARKIND1 ,ONLY : JPIM ,JPRB
18 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK
19
20 USE PARSRTM , ONLY : JPLAY, JPBAND, JPB1, JPB2
21 USE YOESRTOP, ONLY : EXTLIQ1, SSALIQ1, ASYLIQ1 &
22 & , EXTICE3, SSAICE3, ASYICE3, FDLICE3 &
23 & , FDELTA , EXTCOICE, SSACOICE, GICE, FORWICE &
24 & , EXTCOLIQ, SSACOLIQ, GLIQ, FORWLIQ
25
26 ! ------------------------------------------------------------------
27
28 IMPLICIT NONE
29
30 !-- real arguments
31
32 INTEGER(KIND=JPIM),INTENT(IN) :: KLEV
33 INTEGER(KIND=JPIM),INTENT(OUT) :: K_ICLDATM
34 INTEGER(KIND=JPIM),INTENT(IN) :: K_INFLAG
35 INTEGER(KIND=JPIM),INTENT(IN) :: K_ICEFLAG
36 INTEGER(KIND=JPIM),INTENT(IN) :: K_LIQFLAG
37 INTEGER(KIND=JPIM),INTENT(IN) :: K_NSTR
38 REAL(KIND=JPRB) ,INTENT(IN) :: P_CLDFRAC(JPLAY)
39 REAL(KIND=JPRB) ,INTENT(IN) :: P_CLDDAT1(JPLAY)
40 REAL(KIND=JPRB) ,INTENT(IN) :: P_CLDDAT2(JPLAY)
41 REAL(KIND=JPRB) ,INTENT(IN) :: P_CLDDAT3(JPLAY)
42 REAL(KIND=JPRB) ,INTENT(IN) :: P_CLDDAT4(JPLAY)
43 REAL(KIND=JPRB) ,INTENT(IN) :: P_CLDDATMOM(0:16,JPLAY)
44 REAL(KIND=JPRB) ,INTENT(INOUT) :: P_TAUCLDORIG(JPLAY,JPBAND)
45 REAL(KIND=JPRB) ,INTENT(INOUT) :: P_TAUCLOUD(JPLAY,JPBAND)
46 REAL(KIND=JPRB) ,INTENT(OUT) :: P_SSACLOUD(JPLAY,JPBAND)
47 REAL(KIND=JPRB) ,INTENT(OUT) :: P_XMOM(0:16,JPLAY,JPBAND)
48 !-- integer arguments
49
50 !-- real locals
51 REAL(KIND=JPRB) :: Z_EPS
52 REAL(KIND=JPRB) :: Z_TAUCLDORIG_A, Z_FFP, Z_FFP1, Z_FFPSSA, Z_SSACLOUD_A, Z_TAUCLOUD_A
53 REAL(KIND=JPRB) :: Z_CWP, Z_FICE, Z_RADICE, Z_FACTOR, Z_FINT, Z_FLIQ, Z_RADLIQ
54 REAL(KIND=JPRB) :: Z_TAUICEORIG, Z_SCATICE, Z_SSAICE, Z_TAUICE &
55 & , Z_TAULIQORIG, Z_SCATLIQ, Z_SSALIQ, Z_TAULIQ
56
57 !-- integer locals
58 INTEGER(KIND=JPIM) :: I_NCBANDS, I_NLAYERS
59 INTEGER(KIND=JPIM) :: IB, IB1, IB2, I_LAY, ISTR , INDEX
60 INTEGER(KIND=JPIM) :: I_NDBUG
61 REAL(KIND=JPRB) :: ZHOOK_HANDLE
62
63 ! INCLUDE 'param.f'
64
65 ! COMMON /CONTROL/ IAER, NSTR, IOUT, ISTART, IEND, ICLD
66
67 ! COMMON /PROFILE/ NLAYERS,PAVEL(MXLAY),TAVEL(MXLAY),
68 ! & PZ(0:MXLAY),TZ(0:MXLAY)
69
70 ! COMMON /CLOUDIN/ INFLAG,CLDDAT1(MXLAY),CLDDAT2(MXLAY),
71 ! & ICEFLAG,LIQFLAG,CLDDAT3(MXLAY),CLDDAT4(MXLAY),
72 ! & CLDDATMOM(0:16,MXLAY)
73
74 ! COMMON /CLOUDDAT/ NCBANDS,CLDFRAC(MXLAY),
75 ! & TAUCLOUD(MXLAY,NBANDS),SSACLOUD(MXLAY,NBANDS),
76 ! & XMOM(0:16,MXLAY,NBANDS),TAUCLDORIG(MXLAY,NBANDS)
77
78 ! COMMON /HVERSN/ HVRRTM,HVRRTR,HVRATM,HVRSET,HVRTAU,
79 ! * HVDUM1(4),HVRUTL,HVREXT,
80 ! * HVRD1M,HVRR1M,HVREPK,HVRLPK,HVRAER,HVRBKA,
81 ! * HVRBKB,HVRCLD,HVRDIS,HVRLAM,HVRPAR
82
83 ! CHARACTER*15 HVRRTM,HVRRTR,HVRATM,HVRSET,HVRTAU,
84 ! * HVDUM1,HVRUTL,HVREXT,
85 ! * HVRD1M,HVRR1M,HVREPK,HVRLPK,HVRAER,HVRBKA,
86 ! * HVRBKB,HVRCLD,HVRDIS,HVRLAM,HVRPAR
87
88 ! DIMENSION ABSICE0(2), ABSICE1(2,5), ABSICE2(40,16)
89 ! DIMENSION ABSLIQ1(58,16)
90
91 ! DIMENSION ABSCOICE(NBANDS), ABSCOLIQ(NBANDS)
92 ! DIMENSION EXTCOLIQ(NBANDS),SSACOLIQ(NBANDS),GLIQ(NBANDS)
93 ! DIMENSION EXTCOICE(NBANDS),SSACOICE(NBANDS),GICE(NBANDS)
94 ! DIMENSION FORWLIQ(NBANDS),FORWICE(NBANDS)
95
96 ! Added for SW
97 ! DIMENSION EXTLIQ1(58,NBANDS), SSALIQ1(58,NBANDS),
98 ! & ASYLIQ1(58,NBANDS)
99 ! DIMENSION EXTICE3(27,NBANDS), SSAICE3(27,NBANDS),
100 ! & ASYICE3(27,NBANDS),FDLICE3(28,NBANDS)
101 ! DIMENSION FDELTA(NBANDS)
102
103 !DATA EPS /1.E-6/
104
105 ! HVRCLD = '$Revision: 2.6 $'
106
107 IF (LHOOK) CALL DR_HOOK('SRTM_CLDPROP',0,ZHOOK_HANDLE)
108 Z_EPS = 1.E-06_JPRB
109 I_NDBUG = 3
110
111 K_ICLDATM = 0
112 I_NCBANDS = 29
113 I_NLAYERS = KLEV
114 IB1 = JPB1
115 IB2 = JPB2
116
117 IF (I_NDBUG <= 2) THEN
118 print *,'cldprop before loop K_INFLAG, K_ICEFLAG, K_LIQFLAG:',K_INFLAG,K_ICEFLAG,K_LIQFLAG,IB1,IB2
119 ENDIF
120
121 DO I_LAY = 1, I_NLAYERS
122
123 IF (P_CLDFRAC(I_LAY) >= Z_EPS) THEN
124 K_ICLDATM = 1
125
126 IF (I_NDBUG <= 2) THEN
127 print 9101,I_LAY,K_ICLDATM,P_CLDFRAC(I_LAY),P_CLDDAT1(I_LAY),P_CLDDAT2(I_LAY),P_CLDDAT3(I_LAY)&
128 & ,P_CLDDAT4(I_LAY),(P_CLDDATMOM(ISTR,I_LAY),ISTR=0,K_NSTR)
129 9101 format(1x,'Cld :',2I3,f7.4,7E12.5)
130 ENDIF
131
132 ! Ice clouds and water clouds combined.
133 IF (K_INFLAG == 0) THEN
134 Z_TAUCLDORIG_A = P_CLDDAT1(I_LAY)
135 Z_FFP = P_CLDDATMOM(K_NSTR,I_LAY)
136 Z_FFP1 = 1.0 - Z_FFP
137 Z_FFPSSA = 1.0 - Z_FFP * P_CLDDAT2(I_LAY)
138 Z_SSACLOUD_A = Z_FFP1*P_CLDDAT2(I_LAY)/Z_FFPSSA
139 Z_TAUCLOUD_A = Z_FFPSSA*Z_TAUCLDORIG_A
140
141 ! DO IB = 16,NBANDS
142 DO IB = IB1 , IB2
143 P_TAUCLDORIG(I_LAY,IB) = Z_TAUCLDORIG_A
144 P_SSACLOUD(I_LAY,IB) = Z_SSACLOUD_A
145 P_TAUCLOUD(I_LAY,IB) = Z_TAUCLOUD_A
146
147 DO ISTR = 0,K_NSTR
148 P_XMOM(ISTR,I_LAY,IB) = (P_CLDDATMOM(ISTR,I_LAY) - Z_FFP)/ &
149 & (Z_FFP1)
150 ENDDO
151 ENDDO
152
153 ! Separate treatement of ice clouds and water clouds.
154 ELSEIF(K_INFLAG == 2) THEN
155 Z_CWP = P_CLDDAT1(I_LAY)
156 Z_FICE = P_CLDDAT2(I_LAY)
157 Z_RADICE = P_CLDDAT3(I_LAY)
158
159 IF (I_NDBUG <= 1) THEN
160 print 9102,I_LAY,Z_CWP,Z_FICE,Z_RADICE
161 9102 format(1x,'A',I3,3E13.6)
162 ENDIF
163
164 ! Calculation of absorption coefficients due to ice clouds.
165 IF (Z_FICE == 0.0) THEN
166 ! NCBANDS = 29
167 ! DO IB = 16, NCBANDS
168 DO IB = IB1 , IB2
169 EXTCOICE(IB) = 0.0_JPRB
170 SSACOICE(IB) = 1.0_JPRB
171 GICE(IB) = 1.0_JPRB
172 FORWICE(IB) = 0.0_JPRB
173
174 IF (I_NDBUG <= 1) THEN
175 print 9103,I_LAY,Z_FICE,Z_CWP,Z_RADICE,IB,EXTCOICE(IB),SSACOICE(IB),GICE(IB),FORWICE(IB)
176 9103 format(1x,'B',I3,F6.3,2E13.6,I3,4E12.5)
177 ENDIF
178
179 ENDDO
180
181 ELSEIF (K_ICEFLAG == 3) THEN
182 IF (Z_RADICE < 10.0 .OR. Z_RADICE > 140.0) STOP 'ICE EFFECTIVE SIZE OUT OF BOUNDS'
183 ! NCBANDS = 29
184 Z_FACTOR = (Z_RADICE - 5._JPRB)/5._JPRB
185 INDEX = INT(Z_FACTOR)
186 IF (INDEX == 27) INDEX = 26
187 Z_FINT = Z_FACTOR - REAL(INDEX)
188
189 ! DO IB=16,NCBANDS
190 DO IB = IB1 , IB2
191 EXTCOICE(IB) = Z_FICE * (EXTICE3(INDEX,IB) + Z_FINT * &
192 & (EXTICE3(INDEX+1,IB) - EXTICE3(INDEX,IB)))
193 SSACOICE(IB) = SSAICE3(INDEX,IB) + Z_FINT * &
194 & (SSAICE3(INDEX+1,IB) - SSAICE3(INDEX,IB))
195 GICE(IB) = ASYICE3(INDEX,IB) + Z_FINT * &
196 & (ASYICE3(INDEX+1,IB) - ASYICE3(INDEX,IB))
197 FDELTA(IB) = FDLICE3(INDEX,IB) + Z_FINT * &
198 & (FDLICE3(INDEX+1,IB) - FDLICE3(INDEX,IB))
199 if (fdelta(ib) < 0.0) STOP 'FDELTA LESS THAN 0.0'
200 if (fdelta(ib) > 1.0) STOP 'FDELTA GT THAN 1.0'
201 FORWICE(IB) = FDELTA(IB) + 0.5 / SSACOICE(IB)
202 ! See Fu 1996 p. 2067
203 IF (FORWICE(IB) > GICE(IB)) FORWICE(IB) = GICE(IB)
204 ! Check to ensure all calculated quantities are within physical limits.
205 if (extcoice(ib) < 0.0_JPRB) STOP 'ICE EXTINCTION LESS THAN 0.0'
206 if (ssacoice(ib) > 1.0_JPRB) STOP 'ICE SSA GRTR THAN 1.0'
207 if (ssacoice(ib) < 0.0_JPRB) STOP 'ICE SSA LESS THAN 0.0'
208 if (gice(ib) > 1.0_JPRB) STOP 'ICE ASYM GRTR THAN 1.0'
209 if (gice(ib) < 0.0_JPRB) STOP 'ICE ASYM LESS THAN 0.0'
210
211 IF (I_NDBUG <= 1) THEN
212 print 9104,I_LAY,Z_FICE,Z_CWP,Z_RADICE,IB,EXTCOICE(IB),SSACOICE(IB),GICE(IB),FORWICE(IB),FDELTA(IB)
213 9104 format(1x,'C',I3,F5.3,2E13.6,I3,5E12.5)
214 ENDIF
215
216 ENDDO
217 ENDIF
218 print *,'end of ice computations for I_LAY=',I_LAY
219
220 ! Calculation of absorption coefficients due to water clouds.
221 Z_FLIQ = 1. - Z_FICE
222 IF (Z_FLIQ == 0.0) THEN
223 ! NCBANDS = 29
224 ! DO IB = 16, NCBANDS
225 DO IB = IB1 , IB2
226 EXTCOLIQ(IB) = 0.0
227 SSACOLIQ(IB) = 1.0
228 GLIQ(IB) = 1.0
229 FORWLIQ(IB) = 0.0
230
231 IF (I_NDBUG <= 1) THEN
232 print 9105,I_LAY,Z_FLIQ,Z_CWP,IB,EXTCOLIQ(IB),SSACOLIQ(IB),GLIQ(IB),FORWLIQ(IB)
233 9105 format(1x,'D',I3,F5.3,1E13.6,I3,4E12.5)
234 ENDIF
235
236 ENDDO
237
238 ELSEIF (K_LIQFLAG == 1) THEN
239 Z_RADLIQ = P_CLDDAT4(I_LAY)
240 IF (Z_RADLIQ < 1.5 .OR. Z_RADLIQ > 60.) STOP 'LIQUID EFFECTIVE RADIUS OUT OF BOUNDS'
241 INDEX = INT(Z_RADLIQ - 1.5)
242 IF (INDEX == 0) INDEX = 1
243 IF (INDEX == 58) INDEX = 57
244 Z_FINT = Z_RADLIQ - 1.5 - REAL(INDEX)
245 ! NCBANDS = 29
246
247 ! DO IB = 16, NCBANDS
248 DO IB = IB1 , IB2
249 EXTCOLIQ(IB) = Z_FLIQ * (EXTLIQ1(INDEX,IB) + Z_FINT * &
250 & (EXTLIQ1(INDEX+1,IB) - (EXTLIQ1(INDEX,IB))))
251 SSACOLIQ(IB) = SSALIQ1(INDEX,IB) + Z_FINT * &
252 & (SSALIQ1(INDEX+1,IB) - SSALIQ1(INDEX,IB))
253 GLIQ(IB) = ASYLIQ1(INDEX,IB) + Z_FINT * &
254 & (ASYLIQ1(INDEX+1,IB) - ASYLIQ1(INDEX,IB))
255 FORWLIQ(IB) = GLIQ(IB)**K_NSTR
256 ! Check to ensure all calculated quantities are within physical limits.
257 if (extcoliq(ib) < 0.0_JPRB) STOP 'LIQUID EXTINCTION LESS THAN 0.0'
258 if (ssacoliq(ib) > 1.0_JPRB) STOP 'LIQUID SSA GRTR THAN 1.0'
259 if (ssacoliq(ib) < 0.0_JPRB) STOP 'LIQUID SSA LESS THAN 0.0'
260 if (gliq(ib) > 1.0_JPRB) STOP 'LIQUID ASYM GRTR THAN 1.0'
261 if (gliq(ib) < 0.0_JPRB) STOP 'LIQUID ASYM LESS THAN 0.0'
262
263 IF (I_NDBUG <= 1) THEN
264 print 9106,I_LAY,Z_FLIQ,Z_CWP,Z_RADLIQ,IB,EXTCOLIQ(IB),SSACOLIQ(IB),GLIQ(IB),FORWLIQ(IB)
265 9106 format(1x,'E',I3,F5.3,2E13.6,I3,5E12.5)
266 ENDIF
267
268 ENDDO
269 ENDIF
270
271 IF (I_NDBUG <= 1) THEN
272 print *,'end of liquid water computations for I_LAY=',I_LAY
273 ENDIF
274
275 ! DO IB = 16, NCBANDS
276 DO IB = IB1 , IB2
277 Z_TAULIQORIG = Z_CWP * EXTCOLIQ(IB)
278 Z_TAUICEORIG = Z_CWP * EXTCOICE(IB)
279 P_TAUCLDORIG(I_LAY,IB) = Z_TAULIQORIG + Z_TAUICEORIG
280
281 IF (I_NDBUG <= 1) THEN
282 print 9107,IB,Z_TAULIQORIG,Z_TAUICEORIG,P_TAUCLDORIG(I_LAY,IB),Z_CWP &
283 & ,EXTCOLIQ(IB),EXTCOICE(IB),SSACOLIQ(IB),SSACOICE(IB) &
284 & ,FORWLIQ(IB),FORWICE(IB)
285 9107 format(1x,'F',I3,10E12.5)
286 ENDIF
287
288 Z_SSALIQ = SSACOLIQ(IB) * (1. - FORWLIQ(IB)) / &
289 & (1. - FORWLIQ(IB) * SSACOLIQ(IB))
290 Z_TAULIQ = (1. - FORWLIQ(IB) * SSACOLIQ(IB)) * &
291 & Z_TAULIQORIG
292 Z_SSAICE = SSACOICE(IB) * (1. - FORWICE(IB)) / &
293 & (1. - FORWICE(IB) * SSACOICE(IB))
294 Z_TAUICE = (1. - FORWICE(IB) * SSACOICE(IB)) * &
295 & Z_TAUICEORIG
296 Z_SCATLIQ = Z_SSALIQ * Z_TAULIQ
297 Z_SCATICE = Z_SSAICE * Z_TAUICE
298 P_TAUCLOUD(I_LAY,IB) = Z_TAULIQ + Z_TAUICE
299 P_SSACLOUD(I_LAY,IB) = (Z_SCATLIQ + Z_SCATICE) / &
300 & P_TAUCLOUD(I_LAY,IB)
301 P_XMOM(0,I_LAY,IB) = 1.0
302
303 IF (I_NDBUG <= 1) THEN
304 print 9108,IB,Z_TAULIQORIG,Z_TAUICEORIG,Z_SSALIQ,Z_TAULIQ,Z_SCATLIQ,Z_SSAICE,Z_TAUICE,Z_SCATICE
305 9108 format(1x,'G',I3,8E13.6)
306 ENDIF
307
308 DO ISTR = 1, K_NSTR
309 !This commented code is the standard method for delta-m scaling. In accordance
310 ! with the 1996 Fu paper, equation A.3, the moments for ice were calculated
311 ! as in the uncommented code.
312 ! XMOM(ISTR,LAY,IB) = (SCATLIQ * &
313 ! & (GLIQ(IB)**ISTR - FORWLIQ(IB)) / &
314 ! & (1. - FORWLIQ(IB)) &
315 ! & + SCATICE * &
316 ! & (GICE(IB)**ISTR - FORWICE(IB)) / &
317 ! & (1. - FORWICE(IB)))/(SCATLIQ + SCATICE)
318
319 P_XMOM(ISTR,I_LAY,IB) = (1.0/(Z_SCATLIQ+Z_SCATICE))* &
320 & (Z_SCATLIQ*(GLIQ(IB)**ISTR - FORWLIQ(IB)) / &
321 & (1. - FORWLIQ(IB)) &
322 & + Z_SCATICE * &
323 & ((gice(ib)-forwice(ib))/(1.0-forwice(ib)))**ISTR)
324 ENDDO
325 ENDDO
326
327 ENDIF
328
329 ENDIF
330
331 ENDDO
332
333 IF (I_NDBUG <= 1) THEN
334 print *,'about to leave SRTM_CLDPROP'
335 ENDIF
336
337 !-----------------------------------------------------------------------
338 IF (LHOOK) CALL DR_HOOK('SRTM_CLDPROP',1,ZHOOK_HANDLE)
339 END SUBROUTINE SRTM_CLDPROP
340
341