1 |
molod |
1.11 |
C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_swrad.F,v 1.10 2004/07/14 17:31:57 molod Exp $ |
2 |
molod |
1.1 |
C $Name: $ |
3 |
molod |
1.2 |
|
4 |
|
|
#include "CPP_OPTIONS.h" |
5 |
molod |
1.6 |
#include "PACKAGES_CONFIG.h" |
6 |
molod |
1.4 |
subroutine swrio (nymd,nhms,bi,bj,ndswr,myid,istrip,npcs, |
7 |
molod |
1.5 |
. low_level,mid_level, |
8 |
|
|
. pz,plz,plze,dpres,pkht,pkz,tz,qz,oz,co2, |
9 |
molod |
1.2 |
. albvisdr,albvisdf,albnirdr,albnirdf, |
10 |
molod |
1.3 |
. dtradsw,dtswclr,radswg,swgclr, |
11 |
molod |
1.2 |
. fdifpar,fdirpar,osr,osrclr, |
12 |
molod |
1.5 |
. im,jm,lm,ptop, |
13 |
molod |
1.2 |
. nswcld,cldsw,cswmo,nswlz,swlz, |
14 |
|
|
. lpnt,imstturb,qliqave,fccave,landtype,xlats,xlons) |
15 |
molod |
1.1 |
|
16 |
|
|
implicit none |
17 |
molod |
1.2 |
#ifdef ALLOW_DIAGNOSTICS |
18 |
molod |
1.7 |
#include "SIZE.h" |
19 |
|
|
#include "diagnostics_SIZE.h" |
20 |
molod |
1.2 |
#include "diagnostics.h" |
21 |
|
|
#endif |
22 |
molod |
1.1 |
|
23 |
|
|
c Input Variables |
24 |
|
|
c --------------- |
25 |
molod |
1.5 |
integer nymd,nhms,bi,bj,ndswr,myid,istrip,npcs |
26 |
|
|
integer mid_level,low_level |
27 |
molod |
1.2 |
integer im,jm,lm |
28 |
molod |
1.5 |
real ptop |
29 |
|
|
real pz(im,jm),plz(im,jm,lm),plze(im,jm,lm+1),dpres(im,jm,lm) |
30 |
|
|
real pkht(im,jm,lm+1),pkz(im,jm,lm) |
31 |
|
|
real tz(im,jm,lm),qz(im,jm,lm) |
32 |
|
|
real oz(im,jm,lm) |
33 |
|
|
real co2 |
34 |
|
|
real albvisdr(im,jm),albvisdf(im,jm),albnirdr(im,jm) |
35 |
|
|
real albnirdf(im,jm) |
36 |
|
|
real radswg(im,jm),swgclr(im,jm),fdifpar(im,jm),fdirpar(im,jm) |
37 |
|
|
real osr(im,jm),osrclr(im,jm),dtradsw(im,jm,lm),dtswclr(im,jm,lm) |
38 |
molod |
1.2 |
integer nswcld,nswlz |
39 |
molod |
1.5 |
real cldsw(im,jm,lm),cswmo(im,jm,lm),swlz(im,jm,lm) |
40 |
molod |
1.2 |
logical lpnt |
41 |
|
|
integer imstturb |
42 |
molod |
1.5 |
real qliqave(im,jm,lm),fccave(im,jm,lm) |
43 |
molod |
1.2 |
integer landtype(im,jm) |
44 |
molod |
1.5 |
real xlats(im,jm),xlons(im,jm) |
45 |
molod |
1.1 |
|
46 |
|
|
c Local Variables |
47 |
|
|
c --------------- |
48 |
molod |
1.5 |
integer i,j,L,nn,nsecf |
49 |
molod |
1.8 |
integer ntmstp,nymd2,nhms2 |
50 |
|
|
real getcon,grav,cp,undef |
51 |
molod |
1.1 |
real ra,alf,reffw,reffi,tminv |
52 |
|
|
|
53 |
molod |
1.2 |
parameter ( reffw = 10.0 ) |
54 |
|
|
parameter ( reffi = 65.0 ) |
55 |
molod |
1.1 |
|
56 |
molod |
1.5 |
real tdry(im,jm,lm) |
57 |
|
|
real TEMP1(im,jm) |
58 |
|
|
real TEMP2(im,jm) |
59 |
|
|
real zenith (im,jm) |
60 |
|
|
real cldtot (im,jm,lm) |
61 |
|
|
real cldmxo (im,jm,lm) |
62 |
|
|
real totcld (im,jm) |
63 |
|
|
real cldlow (im,jm) |
64 |
|
|
real cldmid (im,jm) |
65 |
|
|
real cldhi (im,jm) |
66 |
|
|
real taulow (im,jm) |
67 |
|
|
real taumid (im,jm) |
68 |
|
|
real tauhi (im,jm) |
69 |
|
|
real tautype(im,jm,lm,3) |
70 |
|
|
real tau(im,jm,lm) |
71 |
|
|
real albedo(im,jm) |
72 |
|
|
|
73 |
|
|
real PK(ISTRIP,lm) |
74 |
|
|
real qzl(ISTRIP,lm),CLRO(ISTRIP,lm) |
75 |
|
|
real TZL(ISTRIP,lm) |
76 |
|
|
real OZL(ISTRIP,lm) |
77 |
|
|
real PLE(ISTRIP,lm+1) |
78 |
molod |
1.8 |
real COSZ(ISTRIP) |
79 |
molod |
1.5 |
real dpstrip(ISTRIP,lm) |
80 |
|
|
|
81 |
|
|
real albuvdr(ISTRIP),albuvdf(ISTRIP) |
82 |
|
|
real albirdr(ISTRIP),albirdf(ISTRIP) |
83 |
|
|
real difpar (ISTRIP),dirpar (ISTRIP) |
84 |
|
|
|
85 |
|
|
real fdirir(istrip),fdifir(istrip) |
86 |
|
|
real fdiruv(istrip),fdifuv(istrip) |
87 |
|
|
|
88 |
|
|
real flux(istrip,lm+1) |
89 |
|
|
real fluxclr(istrip,lm+1) |
90 |
|
|
real dtsw(istrip,lm) |
91 |
|
|
real dtswc(istrip,lm) |
92 |
|
|
|
93 |
|
|
real taul (istrip,lm) |
94 |
|
|
real reff (istrip,lm,2) |
95 |
|
|
real tauc (istrip,lm,2) |
96 |
|
|
real taua (istrip,lm) |
97 |
|
|
real tstrip (istrip) |
98 |
molod |
1.1 |
|
99 |
molod |
1.5 |
logical first |
100 |
|
|
data first /.true./ |
101 |
molod |
1.1 |
|
102 |
|
|
C ********************************************************************** |
103 |
|
|
C **** INITIALIZATION **** |
104 |
|
|
C ********************************************************************** |
105 |
|
|
|
106 |
|
|
grav = getcon('GRAVITY') |
107 |
|
|
cp = getcon('CP') |
108 |
|
|
undef = getcon('UNDEF') |
109 |
|
|
|
110 |
|
|
NTMSTP = nsecf(NDSWR) |
111 |
|
|
TMINV = 1./float(ntmstp) |
112 |
|
|
|
113 |
|
|
C Compute Temperature from Theta |
114 |
|
|
C ------------------------------ |
115 |
|
|
do L=1,lm |
116 |
|
|
do j=1,jm |
117 |
|
|
do i=1,im |
118 |
|
|
tdry(i,j,L) = tz(i,j,L)*pkz(i,j,L) |
119 |
|
|
enddo |
120 |
|
|
enddo |
121 |
|
|
enddo |
122 |
|
|
|
123 |
|
|
if (first .and. myid.eq.0 ) then |
124 |
|
|
print * |
125 |
|
|
print *,'Low-Level Clouds are Grouped between levels: ', |
126 |
|
|
. lm,' and ',low_level |
127 |
|
|
print *,'Mid-Level Clouds are Grouped between levels: ', |
128 |
|
|
. low_level-1,' and ',mid_level |
129 |
|
|
print * |
130 |
|
|
first = .false. |
131 |
|
|
endif |
132 |
|
|
|
133 |
|
|
C ********************************************************************** |
134 |
|
|
C **** CALCULATE COSINE OF THE ZENITH ANGLE **** |
135 |
|
|
C ********************************************************************** |
136 |
|
|
|
137 |
molod |
1.5 |
CALL ASTRO ( NYMD, NHMS, XLATS,XLONS, im*jm, TEMP1,RA ) |
138 |
molod |
1.1 |
NYMD2 = NYMD |
139 |
|
|
NHMS2 = NHMS |
140 |
|
|
CALL TICK ( NYMD2, NHMS2, NTMSTP ) |
141 |
molod |
1.5 |
CALL ASTRO ( NYMD2, NHMS2, XLATS,XLONS, im*jm, TEMP2,RA ) |
142 |
molod |
1.1 |
|
143 |
|
|
do j = 1,jm |
144 |
|
|
do i = 1,im |
145 |
|
|
zenith(I,j) = TEMP1(I,j) + TEMP2(I,j) |
146 |
|
|
IF( zenith(I,j) .GT. 0.0 ) THEN |
147 |
|
|
zenith(I,j) = (2./3.)*( zenith(I,j)-TEMP2(I,j)*TEMP1(I,j) |
148 |
|
|
. / zenith(I,j) ) |
149 |
|
|
ENDIF |
150 |
|
|
ENDDO |
151 |
|
|
ENDDO |
152 |
|
|
|
153 |
|
|
|
154 |
|
|
C ********************************************************************** |
155 |
|
|
c **** Compute Two-Dimension Total Cloud Fraction (0-1) **** |
156 |
|
|
C ********************************************************************** |
157 |
|
|
|
158 |
|
|
c Initialize Clear Sky Probability for Low, Mid, and High Clouds |
159 |
|
|
c -------------------------------------------------------------- |
160 |
|
|
do j =1,jm |
161 |
|
|
do i =1,im |
162 |
|
|
cldlow(i,j) = 0.0 |
163 |
|
|
cldmid(i,j) = 0.0 |
164 |
|
|
cldhi (i,j) = 0.0 |
165 |
|
|
enddo |
166 |
|
|
enddo |
167 |
|
|
|
168 |
|
|
c Adjust cloud fractions and cloud liquid water due to moist turbulence |
169 |
|
|
c --------------------------------------------------------------------- |
170 |
|
|
if(imstturb.ne.0) then |
171 |
|
|
do L =1,lm |
172 |
|
|
do j =1,jm |
173 |
|
|
do i =1,im |
174 |
molod |
1.2 |
cldtot(i,j,L)=min(1.0,max(cldsw(i,j,L),fccave(i,j,L)/imstturb)) |
175 |
|
|
cldmxo(i,j,L)=min(1.0,cswmo(i,j,L)) |
176 |
|
|
swlz(i,j,L)=swlz(i,j,L)+qliqave(i,j,L)/imstturb |
177 |
molod |
1.1 |
enddo |
178 |
|
|
enddo |
179 |
|
|
enddo |
180 |
|
|
else |
181 |
|
|
do L =1,lm |
182 |
|
|
do j =1,jm |
183 |
|
|
do i =1,im |
184 |
|
|
cldtot(i,j,L) = min( 1.0,cldsw(i,j,L) ) |
185 |
|
|
cldmxo(i,j,L) = min( 1.0,cswmo(i,j,L) ) |
186 |
|
|
enddo |
187 |
|
|
enddo |
188 |
|
|
enddo |
189 |
|
|
endif |
190 |
|
|
|
191 |
|
|
c Compute 3-D Cloud Fractions |
192 |
|
|
c --------------------------- |
193 |
|
|
if( nswcld.ne.0 ) then |
194 |
|
|
do L = 1,lm |
195 |
|
|
do j = 1,jm |
196 |
|
|
do i = 1,im |
197 |
|
|
c Compute Low-Mid-High Maximum Overlap Cloud Fractions |
198 |
|
|
c ---------------------------------------------------- |
199 |
|
|
if( L.lt.mid_level ) then |
200 |
|
|
cldhi (i,j) = max( cldtot(i,j,L),cldhi (i,j) ) |
201 |
|
|
elseif( L.ge.mid_level .and. |
202 |
|
|
. L.lt.low_level ) then |
203 |
|
|
cldmid(i,j) = max( cldtot(i,j,L),cldmid(i,j) ) |
204 |
|
|
elseif( L.ge.low_level ) then |
205 |
|
|
cldlow(i,j) = max( cldtot(i,j,L),cldlow(i,j) ) |
206 |
|
|
endif |
207 |
|
|
|
208 |
|
|
enddo |
209 |
|
|
enddo |
210 |
|
|
enddo |
211 |
|
|
endif |
212 |
|
|
|
213 |
|
|
c Totcld => Product of Clear Sky Probabilities for Low, Mid, and High Clouds |
214 |
|
|
c -------------------------------------------------------------------------- |
215 |
|
|
do j = 1,jm |
216 |
|
|
do i = 1,im |
217 |
|
|
totcld(i,j) = 1.0 - (1.-cldhi (i,j)) |
218 |
|
|
. * (1.-cldmid(i,j)) |
219 |
|
|
. * (1.-cldlow(i,j)) |
220 |
|
|
enddo |
221 |
|
|
enddo |
222 |
|
|
|
223 |
|
|
c Compute Cloud Diagnostics |
224 |
|
|
c ------------------------- |
225 |
|
|
if(icldfrc.gt.0) then |
226 |
|
|
do j=1,jm |
227 |
|
|
do i=1,im |
228 |
molod |
1.4 |
qdiag(i,j,icldfrc,bi,bj) = qdiag(i,j,icldfrc,bi,bj) + totcld(i,j) |
229 |
molod |
1.1 |
enddo |
230 |
|
|
enddo |
231 |
|
|
ncldfrc = ncldfrc + 1 |
232 |
|
|
endif |
233 |
|
|
|
234 |
|
|
if( icldras.gt.0 ) then |
235 |
|
|
do L=1,lm |
236 |
|
|
do j=1,jm |
237 |
|
|
do i=1,im |
238 |
molod |
1.4 |
qdiag(i,j,icldras+L-1,bi,bj) = qdiag(i,j,icldras+L-1,bi,bj) + |
239 |
|
|
. cswmo(i,j,L) |
240 |
molod |
1.1 |
enddo |
241 |
|
|
enddo |
242 |
|
|
enddo |
243 |
|
|
ncldras = ncldras + 1 |
244 |
|
|
endif |
245 |
|
|
|
246 |
|
|
if( icldtot.gt.0 ) then |
247 |
|
|
do L=1,lm |
248 |
|
|
do j=1,jm |
249 |
|
|
do i=1,im |
250 |
molod |
1.4 |
qdiag(i,j,icldtot+L-1,bi,bj) = qdiag(i,j,icldtot+L-1,bi,bj) + |
251 |
|
|
. cldtot(i,j,L) |
252 |
molod |
1.1 |
enddo |
253 |
|
|
enddo |
254 |
|
|
enddo |
255 |
|
|
ncldtot = ncldtot + 1 |
256 |
|
|
endif |
257 |
|
|
|
258 |
|
|
if( icldlow.gt.0 ) then |
259 |
|
|
do j=1,jm |
260 |
|
|
do i=1,im |
261 |
molod |
1.4 |
qdiag(i,j,icldlow,bi,bj) = qdiag(i,j,icldlow,bi,bj) + cldlow(i,j) |
262 |
molod |
1.1 |
enddo |
263 |
|
|
enddo |
264 |
|
|
ncldlow = ncldlow + 1 |
265 |
|
|
endif |
266 |
|
|
|
267 |
|
|
if( icldmid.gt.0 ) then |
268 |
|
|
do j=1,jm |
269 |
|
|
do i=1,im |
270 |
molod |
1.4 |
qdiag(i,j,icldmid,bi,bj) = qdiag(i,j,icldmid,bi,bj) + cldmid(i,j) |
271 |
molod |
1.1 |
enddo |
272 |
|
|
enddo |
273 |
|
|
ncldmid = ncldmid + 1 |
274 |
|
|
endif |
275 |
|
|
|
276 |
|
|
if( icldhi.gt.0 ) then |
277 |
|
|
do j=1,jm |
278 |
|
|
do i=1,im |
279 |
molod |
1.4 |
qdiag(i,j,icldhi,bi,bj) = qdiag(i,j,icldhi,bi,bj) + cldhi(i,j) |
280 |
molod |
1.1 |
enddo |
281 |
|
|
enddo |
282 |
|
|
ncldhi = ncldhi + 1 |
283 |
|
|
endif |
284 |
|
|
|
285 |
|
|
if( ilzrad.gt.0 ) then |
286 |
|
|
do L=1,lm |
287 |
|
|
do j=1,jm |
288 |
|
|
do i=1,im |
289 |
molod |
1.4 |
qdiag(i,j,ilzrad+L-1,bi,bj) = qdiag(i,j,ilzrad+L-1,bi,bj) + |
290 |
|
|
. swlz(i,j,L)*1.0e6 |
291 |
molod |
1.1 |
enddo |
292 |
|
|
enddo |
293 |
|
|
enddo |
294 |
|
|
nlzrad = nlzrad + 1 |
295 |
|
|
endif |
296 |
|
|
|
297 |
|
|
c Albedo Diagnostics |
298 |
|
|
c ------------------ |
299 |
|
|
if( ialbvisdr.gt.0 ) then |
300 |
|
|
do j=1,jm |
301 |
|
|
do i=1,im |
302 |
molod |
1.4 |
qdiag(i,j,ialbvisdr,bi,bj) = qdiag(i,j,ialbvisdr,bi,bj) + |
303 |
|
|
. albvisdr(i,j) |
304 |
molod |
1.1 |
enddo |
305 |
|
|
enddo |
306 |
|
|
nalbvisdr = nalbvisdr + 1 |
307 |
|
|
endif |
308 |
|
|
|
309 |
|
|
if( ialbvisdf.gt.0 ) then |
310 |
|
|
do j=1,jm |
311 |
|
|
do i=1,im |
312 |
molod |
1.4 |
qdiag(i,j,ialbvisdf,bi,bj) = qdiag(i,j,ialbvisdf,bi,bj) + |
313 |
|
|
. albvisdf(i,j) |
314 |
molod |
1.1 |
enddo |
315 |
|
|
enddo |
316 |
|
|
nalbvisdf = nalbvisdf + 1 |
317 |
|
|
endif |
318 |
|
|
|
319 |
|
|
if( ialbnirdr.gt.0 ) then |
320 |
|
|
do j=1,jm |
321 |
|
|
do i=1,im |
322 |
molod |
1.4 |
qdiag(i,j,ialbnirdr,bi,bj) = qdiag(i,j,ialbnirdr,bi,bj) + |
323 |
|
|
. albnirdr(i,j) |
324 |
molod |
1.1 |
enddo |
325 |
|
|
enddo |
326 |
|
|
nalbnirdr = nalbnirdr + 1 |
327 |
|
|
endif |
328 |
|
|
|
329 |
|
|
if( ialbnirdf.gt.0 ) then |
330 |
|
|
do j=1,jm |
331 |
|
|
do i=1,im |
332 |
molod |
1.4 |
qdiag(i,j,ialbnirdf,bi,bj) = qdiag(i,j,ialbnirdf,bi,bj) + |
333 |
|
|
. albnirdf(i,j) |
334 |
molod |
1.1 |
enddo |
335 |
|
|
enddo |
336 |
|
|
nalbnirdf = nalbnirdf + 1 |
337 |
|
|
endif |
338 |
|
|
|
339 |
|
|
C Compute Optical Thicknesses and Diagnostics |
340 |
|
|
C ------------------------------------------- |
341 |
molod |
1.2 |
call opthk(tdry,plz,plze,swlz,cldtot,cldmxo,landtype,im,jm,lm, |
342 |
|
|
. tautype) |
343 |
molod |
1.1 |
|
344 |
|
|
do L = 1,lm |
345 |
|
|
do j = 1,jm |
346 |
|
|
do i = 1,im |
347 |
molod |
1.2 |
tau(i,j,L) = tautype(i,j,L,1)+tautype(i,j,L,2)+tautype(i,j,L,3) |
348 |
molod |
1.1 |
enddo |
349 |
|
|
enddo |
350 |
|
|
enddo |
351 |
|
|
|
352 |
|
|
if( itauave.gt.0 ) then |
353 |
|
|
do L=1,lm |
354 |
|
|
do j=1,jm |
355 |
|
|
do i=1,im |
356 |
molod |
1.4 |
qdiag(i,j,itauave+L-1,bi,bj) = qdiag(i,j,itauave+L-1,bi,bj) + |
357 |
molod |
1.2 |
. tau(i,j,L)*100/(plze(i,j,L+1)-plze(i,j,L)) |
358 |
molod |
1.1 |
enddo |
359 |
|
|
enddo |
360 |
|
|
enddo |
361 |
|
|
ntauave = ntauave + 1 |
362 |
|
|
endif |
363 |
|
|
|
364 |
|
|
if( itaucld.gt.0 ) then |
365 |
|
|
do L=1,lm |
366 |
|
|
do j=1,jm |
367 |
|
|
do i=1,im |
368 |
molod |
1.2 |
if( cldtot(i,j,L).ne.0.0 ) then |
369 |
molod |
1.4 |
qdiag(i,j,itaucld +L-1,bi,bj) = qdiag(i,j,itaucld +L-1,bi,bj) + |
370 |
molod |
1.2 |
. tau(i,j,L)*100/(plze(i,j,L+1)-plze(i,j,L)) |
371 |
molod |
1.4 |
qdiag(i,j,itaucldc+L-1,bi,bj) = |
372 |
|
|
. qdiag(i,j,itaucldc+L-1,bi,bj) + 1.0 |
373 |
molod |
1.2 |
endif |
374 |
molod |
1.1 |
enddo |
375 |
|
|
enddo |
376 |
|
|
enddo |
377 |
|
|
endif |
378 |
|
|
|
379 |
|
|
c Compute Low, Mid, and High Cloud Optical Depth Diagnostics |
380 |
|
|
c ---------------------------------------------------------- |
381 |
|
|
if( itaulow.ne.0 ) then |
382 |
molod |
1.4 |
do j = 1,jm |
383 |
|
|
do i = 1,im |
384 |
|
|
if( cldlow(i,j).ne.0.0 ) then |
385 |
|
|
taulow(i,j) = 0.0 |
386 |
|
|
do L = low_level,lm |
387 |
|
|
taulow(i,j) = taulow(i,j) + tau(i,j,L) |
388 |
|
|
enddo |
389 |
|
|
qdiag(i,j,itaulow,bi,bj ) = qdiag(i,j,itaulow,bi,bj ) + |
390 |
|
|
. taulow(i,j) |
391 |
|
|
qdiag(i,j,itaulowc,bi,bj) = qdiag(i,j,itaulowc,bi,bj) + 1.0 |
392 |
|
|
endif |
393 |
|
|
enddo |
394 |
|
|
enddo |
395 |
molod |
1.1 |
endif |
396 |
|
|
|
397 |
|
|
if( itaumid.ne.0 ) then |
398 |
molod |
1.4 |
do j = 1,jm |
399 |
|
|
do i = 1,im |
400 |
|
|
if( cldmid(i,j).ne.0.0 ) then |
401 |
|
|
taumid(i,j) = 0.0 |
402 |
|
|
do L = mid_level,low_level+1 |
403 |
|
|
taumid(i,j) = taumid(i,j) + tau(i,j,L) |
404 |
|
|
enddo |
405 |
|
|
qdiag(i,j,itaumid,bi,bj ) = qdiag(i,j,itaumid,bi,bj ) + |
406 |
|
|
. taumid(i,j) |
407 |
|
|
qdiag(i,j,itaumidc,bi,bj) = qdiag(i,j,itaumidc,bi,bj) + 1.0 |
408 |
|
|
endif |
409 |
|
|
enddo |
410 |
|
|
enddo |
411 |
molod |
1.1 |
endif |
412 |
|
|
|
413 |
|
|
if( itauhi.ne.0 ) then |
414 |
molod |
1.4 |
do j = 1,jm |
415 |
|
|
do i = 1,im |
416 |
|
|
if( cldhi(i,j).ne.0.0 ) then |
417 |
|
|
tauhi(i,j) = 0.0 |
418 |
|
|
do L = 1,mid_level+1 |
419 |
|
|
tauhi(i,j) = tauhi(i,j) + tau(i,j,L) |
420 |
|
|
enddo |
421 |
|
|
qdiag(i,j,itauhi,bi,bj ) = qdiag(i,j,itauhi,bi,bj ) + |
422 |
|
|
. tauhi(i,j) |
423 |
|
|
qdiag(i,j,itauhic,bi,bj) = qdiag(i,j,itauhic,bi,bj) + 1.0 |
424 |
|
|
endif |
425 |
|
|
enddo |
426 |
|
|
enddo |
427 |
molod |
1.1 |
endif |
428 |
|
|
|
429 |
|
|
C*********************************************************************** |
430 |
|
|
C **** LOOP OVER GLOBAL REGIONS **** |
431 |
|
|
C ********************************************************************** |
432 |
|
|
|
433 |
|
|
do 1000 nn = 1,npcs |
434 |
|
|
|
435 |
|
|
C ********************************************************************** |
436 |
|
|
C **** VARIABLE INITIALIZATION **** |
437 |
|
|
C ********************************************************************** |
438 |
|
|
|
439 |
|
|
CALL STRIP ( zenith,COSZ,im*jm,ISTRIP,1,NN ) |
440 |
|
|
|
441 |
molod |
1.5 |
CALL STRIP ( plze, ple ,im*jm,ISTRIP,lm+1,NN) |
442 |
|
|
CALL STRIP ( pkz , pk ,im*jm,ISTRIP,lm ,NN) |
443 |
|
|
CALL STRIP ( dpres,dpstrip,im*jm,ISTRIP,lm ,NN) |
444 |
|
|
CALL STRIP ( tdry, tzl ,im*jm,ISTRIP,lm ,NN) |
445 |
|
|
CALL STRIP ( qz , qzl ,im*jm,ISTRIP,lm ,NN) |
446 |
|
|
CALL STRIP ( oz , ozl ,im*jm,ISTRIP,lm ,NN) |
447 |
|
|
CALL STRIP ( tau , taul ,im*jm,ISTRIP,lm ,NN) |
448 |
molod |
1.1 |
|
449 |
|
|
CALL STRIP ( albvisdr,albuvdr,im*jm,ISTRIP,1,NN ) |
450 |
|
|
CALL STRIP ( albvisdf,albuvdf,im*jm,ISTRIP,1,NN ) |
451 |
|
|
CALL STRIP ( albnirdr,albirdr,im*jm,ISTRIP,1,NN ) |
452 |
|
|
CALL STRIP ( albnirdf,albirdf,im*jm,ISTRIP,1,NN ) |
453 |
|
|
|
454 |
|
|
call strip ( cldtot,clro,im*jm,istrip,lm,nn ) |
455 |
|
|
|
456 |
|
|
c Partition Tau between Water and Ice Particles |
457 |
|
|
c --------------------------------------------- |
458 |
|
|
do L= 1,lm |
459 |
|
|
do i= 1,istrip |
460 |
|
|
alf = min( max((tzl(i,l)-253.15)/20.,0.) ,1.) |
461 |
|
|
taua(i,L) = 0. |
462 |
|
|
|
463 |
|
|
if( alf.ne.0.0 .and. alf.ne.1.0 ) then |
464 |
|
|
tauc(i,L,1) = taul(i,L)/(1.+alf/(1-alf) * (reffi/reffw*0.80) ) |
465 |
|
|
tauc(i,L,2) = taul(i,L)/(1.+(1-alf)/alf * (reffw/reffi*1.25) ) |
466 |
|
|
endif |
467 |
|
|
|
468 |
|
|
if( alf.eq.0.0 ) then |
469 |
|
|
tauc(i,L,1) = taul(i,L) |
470 |
|
|
tauc(i,L,2) = 0.0 |
471 |
|
|
endif |
472 |
|
|
|
473 |
|
|
if( alf.eq.1.0 ) then |
474 |
|
|
tauc(i,L,1) = 0.0 |
475 |
|
|
tauc(i,L,2) = taul(i,L) |
476 |
|
|
endif |
477 |
|
|
|
478 |
|
|
reff(i,L,1) = reffi |
479 |
|
|
reff(i,L,2) = reffw |
480 |
|
|
enddo |
481 |
|
|
enddo |
482 |
|
|
|
483 |
|
|
call sorad ( istrip,1,1,lm,ple,tzl,qzl,ozl,co2, |
484 |
|
|
. tauc,reff,clro,mid_level,low_level,taua, |
485 |
|
|
. albirdr,albirdf,albuvdr,albuvdf,cosz, |
486 |
|
|
. flux,fluxclr,fdirir,fdifir,dirpar,difpar, |
487 |
|
|
. fdiruv,fdifuv ) |
488 |
|
|
|
489 |
|
|
C ********************************************************************** |
490 |
|
|
C **** Compute Mass-Weighted Theta Tendencies from Fluxes **** |
491 |
|
|
C ********************************************************************** |
492 |
|
|
|
493 |
|
|
do l=1,lm |
494 |
|
|
do i=1,istrip |
495 |
molod |
1.5 |
alf = grav*(ple(i,L+1)-ptop)/(cp*dpstrip(i,L)*100) |
496 |
molod |
1.1 |
dtsw (i,L) = alf*( flux (i,L)-flux (i,L+1) )/pk(i,L) |
497 |
|
|
dtswc(i,L) = alf*( fluxclr(i,L)-fluxclr(i,L+1) )/pk(i,L) |
498 |
|
|
enddo |
499 |
|
|
enddo |
500 |
|
|
|
501 |
|
|
call paste ( dtsw , dtradsw ,istrip,im*jm,lm,nn ) |
502 |
|
|
call paste ( dtswc, dtswclr ,istrip,im*jm,lm,nn ) |
503 |
|
|
|
504 |
|
|
call paste ( flux (1,1),osr ,istrip,im*jm,1,nn ) |
505 |
|
|
call paste ( fluxclr(1,1),osrclr,istrip,im*jm,1,nn ) |
506 |
|
|
|
507 |
|
|
call paste ( flux (1,lm+1),radswg,istrip,im*jm,1,nn ) |
508 |
|
|
call paste ( fluxclr(1,lm+1),swgclr,istrip,im*jm,1,nn ) |
509 |
|
|
|
510 |
|
|
call paste ( dirpar,fdirpar,istrip,im*jm,1,nn ) |
511 |
|
|
call paste ( difpar,fdifpar,istrip,im*jm,1,nn ) |
512 |
|
|
|
513 |
|
|
c Calculate Mean Albedo |
514 |
|
|
c --------------------- |
515 |
|
|
do i=1,istrip |
516 |
molod |
1.2 |
if( cosz(i).gt.0.0 ) then |
517 |
|
|
tstrip(i) = 1.0 - flux(i,lm+1)/ |
518 |
|
|
. ( fdirir(i)+fdifir(i)+dirpar(i)+difpar(i) + fdiruv(i)+fdifuv(i) ) |
519 |
molod |
1.1 |
if( tstrip(i).lt.0.0 ) tstrip(i) = undef |
520 |
molod |
1.2 |
else |
521 |
|
|
tstrip(i) = undef |
522 |
|
|
endif |
523 |
molod |
1.1 |
enddo |
524 |
|
|
call paste ( tstrip,albedo,istrip,im*jm,1,nn ) |
525 |
|
|
|
526 |
|
|
1000 CONTINUE |
527 |
|
|
|
528 |
|
|
c Mean Albedo Diagnostic |
529 |
|
|
c ---------------------- |
530 |
|
|
if( ialbedo.gt.0 ) then |
531 |
|
|
do j=1,jm |
532 |
|
|
do i=1,im |
533 |
|
|
if( albedo(i,j).ne.undef ) then |
534 |
molod |
1.4 |
qdiag(i,j,ialbedo,bi,bj ) = qdiag(i,j,ialbedo,bi,bj )+albedo(i,j) |
535 |
|
|
qdiag(i,j,ialbedoc,bi,bj) = qdiag(i,j,ialbedoc,bi,bj) + 1.0 |
536 |
molod |
1.1 |
endif |
537 |
|
|
enddo |
538 |
|
|
enddo |
539 |
|
|
endif |
540 |
|
|
|
541 |
|
|
C ********************************************************************** |
542 |
|
|
C **** ZERO-OUT OR BUMP COUNTERS **** |
543 |
|
|
C ********************************************************************** |
544 |
|
|
|
545 |
|
|
nswlz = 0 |
546 |
|
|
nswcld = 0 |
547 |
|
|
imstturb = 0 |
548 |
|
|
|
549 |
|
|
do L = 1,lm |
550 |
|
|
do j = 1,jm |
551 |
|
|
do i = 1,im |
552 |
|
|
fccave(i,j,L) = 0. |
553 |
|
|
qliqave(i,j,L) = 0. |
554 |
|
|
enddo |
555 |
|
|
enddo |
556 |
|
|
enddo |
557 |
|
|
|
558 |
|
|
return |
559 |
|
|
end |
560 |
|
|
subroutine opthk ( tl,pl,ple,lz,cf,cfm,lwi,im,jm,lm,tau ) |
561 |
|
|
C*********************************************************************** |
562 |
|
|
C |
563 |
|
|
C PURPOSE: |
564 |
|
|
C ======== |
565 |
|
|
C Compute cloud optical thickness using temperature and |
566 |
|
|
C cloud fractions. |
567 |
|
|
C |
568 |
|
|
C INPUT: |
569 |
|
|
C ====== |
570 |
|
|
C tl ......... Temperature at Model Layers (K) |
571 |
|
|
C pl ......... Model Layer Pressures (mb) |
572 |
|
|
C ple ........ Model Edge Pressures (mb) |
573 |
|
|
C lz ......... Diagnosed Convective and Large-Scale Cloud Water (g/g) |
574 |
|
|
C cf ......... Total Cloud Fraction (Random + Maximum Overlap) |
575 |
|
|
C cfm ........ Maximum Overlap Cloud Fraction |
576 |
|
|
C lwi ........ Surface Land Type |
577 |
|
|
C im ......... Longitudinal dimension |
578 |
|
|
C jm ......... Latitudinal dimension |
579 |
|
|
C lm ......... Vertical dimension |
580 |
|
|
C |
581 |
|
|
C OUTPUT: |
582 |
|
|
C ======= |
583 |
|
|
C tau ........ Cloud Optical Thickness (non-dimensional) |
584 |
|
|
C tau(im,jm,lm,1): Suspended Ice |
585 |
|
|
C tau(im,jm,lm,2): Suspended Water |
586 |
|
|
C tau(im,jm,lm,3): Raindrops |
587 |
|
|
C |
588 |
|
|
C*********************************************************************** |
589 |
|
|
C* GODDARD LABORATORY FOR ATMOSPHERES * |
590 |
|
|
C*********************************************************************** |
591 |
|
|
|
592 |
|
|
implicit none |
593 |
|
|
|
594 |
|
|
integer im,jm,lm,i,j,L |
595 |
|
|
|
596 |
|
|
real tl(im,jm,lm) |
597 |
|
|
real pl(im,jm,lm) |
598 |
|
|
real ple(im,jm,lm+1) |
599 |
|
|
real lz(im,jm,lm) |
600 |
|
|
real cf(im,jm,lm) |
601 |
|
|
real cfm(im,jm,lm) |
602 |
|
|
real tau(im,jm,lm,3) |
603 |
|
|
integer lwi(im,jm) |
604 |
|
|
|
605 |
|
|
real dp, alf, fracls, fraccu |
606 |
|
|
real tauice, tauh2o, tauras |
607 |
|
|
|
608 |
|
|
c Compute Cloud Optical Depths |
609 |
|
|
c ---------------------------- |
610 |
|
|
do L=1,lm |
611 |
|
|
do j=1,jm |
612 |
|
|
do i=1,im |
613 |
|
|
alf = min( max(( tl(i,j,L)-233.15)/20.,0.) ,1.) |
614 |
|
|
dp = ple(i,j,L+1)-ple(i,j,L) |
615 |
|
|
tau(i,j,L,1) = 0.0 |
616 |
|
|
tau(i,j,L,2) = 0.0 |
617 |
|
|
tau(i,j,L,3) = 0.0 |
618 |
|
|
if( cf(i,j,L).ne.0.0 ) then |
619 |
|
|
|
620 |
|
|
c Determine fraction of large-scale and cumulus clouds |
621 |
|
|
c ---------------------------------------------------- |
622 |
|
|
fracls = ( cf(i,j,L)-cfm(i,j,L) )/cf(i,j,L) |
623 |
|
|
fraccu = 1.0-fracls |
624 |
|
|
|
625 |
|
|
c Define tau for large-scale ice and water clouds |
626 |
|
|
c Note: tauice is scaled between (0.02 & .2) for: 0 < lz < 2 mg/kg over Land |
627 |
|
|
c Note: tauh2o is scaled between (0.20 & 20) for: 0 < lz < 5 mg/kg over Land |
628 |
|
|
c Note: tauh2o is scaled between (0.20 & 12) for: 0 < lz < 50 mg/kg over Ocean |
629 |
|
|
c ------------------------------------------------------------------------------- |
630 |
|
|
|
631 |
|
|
c Large-Scale Ice |
632 |
|
|
c --------------- |
633 |
|
|
tauice = max( 0.0002, 0.002*min(500*lz(i,j,L)*1000,1.0) ) |
634 |
|
|
tau(i,j,L,1) = fracls*(1-alf)*tauice*dp |
635 |
|
|
|
636 |
|
|
c Large-Scale Water |
637 |
|
|
c ----------------- |
638 |
molod |
1.2 |
C Over Land |
639 |
molod |
1.1 |
if( lwi(i,j).le.10 ) then |
640 |
molod |
1.2 |
tauh2o = max( 0.0020, 0.200*min(200*lz(i,j,L)*1000,1.0) ) |
641 |
|
|
tau(i,j,L,3) = fracls*alf*tauh2o*dp |
642 |
molod |
1.1 |
else |
643 |
molod |
1.2 |
C Non-Precipitation Clouds Over Ocean |
644 |
|
|
if( lz(i,j,L).eq.0.0 ) then |
645 |
|
|
tauh2o = .12 |
646 |
|
|
tau(i,j,L,2) = fracls*alf*tauh2o*dp |
647 |
|
|
else |
648 |
|
|
C Over Ocean |
649 |
|
|
tauh2o = max( 0.0020, 0.120*min( 20*lz(i,j,L)*1000,1.0) ) |
650 |
|
|
tau(i,j,L,3) = fracls*alf*tauh2o*dp |
651 |
|
|
endif |
652 |
molod |
1.1 |
endif |
653 |
|
|
|
654 |
|
|
c Sub-Grid Convective |
655 |
|
|
c ------------------- |
656 |
|
|
if( tl(i,j,L).gt.210.0 ) then |
657 |
|
|
tauras = .16 |
658 |
|
|
tau(i,j,L,3) = tau(i,j,L,3) + fraccu*tauras*dp |
659 |
|
|
else |
660 |
|
|
tauras = tauice |
661 |
|
|
tau(i,j,L,1) = tau(i,j,L,1) + fraccu*tauras*dp |
662 |
|
|
endif |
663 |
|
|
|
664 |
|
|
c Define tau for large-scale and cumulus clouds |
665 |
|
|
c --------------------------------------------- |
666 |
|
|
ccc tau(i,j,L) = ( fracls*( alf*tauh2o + (1-alf)*tauice ) |
667 |
|
|
ccc . + fraccu*tauras )*dp |
668 |
|
|
|
669 |
|
|
endif |
670 |
|
|
|
671 |
|
|
enddo |
672 |
|
|
enddo |
673 |
|
|
enddo |
674 |
|
|
|
675 |
|
|
return |
676 |
|
|
end |
677 |
|
|
subroutine sorad(m,n,ndim,np,pl,ta,wa,oa,co2, |
678 |
|
|
* taucld,reff,fcld,ict,icb, |
679 |
|
|
* taual,rsirbm,rsirdf,rsuvbm,rsuvdf,cosz, |
680 |
|
|
* flx,flc,fdirir,fdifir,fdirpar,fdifpar, |
681 |
|
|
* fdiruv,fdifuv) |
682 |
|
|
c******************************************************************** |
683 |
|
|
c |
684 |
|
|
c This routine computes solar fluxes due to the absoption by water |
685 |
|
|
c vapor, ozone, co2, o2, clouds, and aerosols and due to the |
686 |
|
|
c scattering by clouds, aerosols, and gases. |
687 |
|
|
c |
688 |
|
|
c This is a vectorized code. It computes the fluxes simultaneous for |
689 |
|
|
c (m x n) soundings, which is a subset of the (m x ndim) soundings. |
690 |
|
|
c In a global climate model, m and ndim correspond to the numbers of |
691 |
|
|
c grid boxes in the zonal and meridional directions, respectively. |
692 |
|
|
c |
693 |
|
|
c Ice and liquid cloud particles are allowed to co-exist in any of the |
694 |
|
|
c np layers. Two sets of cloud parameters are required as inputs, one |
695 |
|
|
c for ice paticles and the other for liquid particles. These parameters |
696 |
|
|
c are optical thickness (taucld) and effective particle size (reff). |
697 |
|
|
c |
698 |
|
|
c If no information is available for reff, a default value of |
699 |
|
|
c 10 micron for liquid water and 75 micron for ice can be used. |
700 |
|
|
c |
701 |
|
|
c Clouds are grouped into high, middle, and low clouds separated by the |
702 |
|
|
c level indices ict and icb. For detail, see the subroutine cldscale. |
703 |
|
|
c |
704 |
|
|
c----- Input parameters: |
705 |
|
|
c units size |
706 |
|
|
c number of soundings in zonal direction (m) n/d 1 |
707 |
|
|
c number of soundings in meridional direction (n) n/d 1 |
708 |
|
|
c maximum number of soundings in n/d 1 |
709 |
|
|
c meridional direction (ndim) |
710 |
|
|
c number of atmospheric layers (np) n/d 1 |
711 |
|
|
c level pressure (pl) mb m*ndim*(np+1) |
712 |
|
|
c layer temperature (ta) k m*ndim*np |
713 |
|
|
c layer specific humidity (wa) gm/gm m*ndim*np |
714 |
|
|
c layer ozone concentration (oa) gm/gm m*ndim*np |
715 |
|
|
c co2 mixing ratio by volumn (co2) parts/part 1 |
716 |
|
|
c cloud optical thickness (taucld) n/d m*ndim*np*2 |
717 |
|
|
c index 1 for ice particles |
718 |
|
|
c index 2 for liquid drops |
719 |
|
|
c effective cloud-particle size (reff) micrometer m*ndim*np*2 |
720 |
|
|
c index 1 for ice particles |
721 |
|
|
c index 2 for liquid drops |
722 |
|
|
c cloud amount (fcld) fraction m*ndim*np |
723 |
|
|
c level index separating high and middle n/d 1 |
724 |
|
|
c clouds (ict) |
725 |
|
|
c level index separating middle and low clouds n/d 1 |
726 |
|
|
c clouds (icb) |
727 |
|
|
c aerosol optical thickness (taual) n/d m*ndim*np |
728 |
|
|
c solar ir surface albedo for beam fraction m*ndim |
729 |
|
|
c radiation (rsirbm) |
730 |
|
|
c solar ir surface albedo for diffuse fraction m*ndim |
731 |
|
|
c radiation (rsirdf) |
732 |
|
|
c uv + par surface albedo for beam fraction m*ndim |
733 |
|
|
c radiation (rsuvbm) |
734 |
|
|
c uv + par surface albedo for diffuse fraction m*ndim |
735 |
|
|
c radiation (rsuvdf) |
736 |
|
|
c cosine of solar zenith angle (cosz) n/d m*ndim |
737 |
|
|
c |
738 |
|
|
c----- Output parameters |
739 |
|
|
c |
740 |
|
|
c all-sky flux (downward minus upward) (flx) fraction m*ndim*(np+1) |
741 |
|
|
c clear-sky flux (downward minus upward) (flc) fraction m*ndim*(np+1) |
742 |
|
|
c all-sky direct downward ir (0.7-10 micron) |
743 |
|
|
c flux at the surface (fdirir) fraction m*ndim |
744 |
|
|
c all-sky diffuse downward ir flux at |
745 |
|
|
c the surface (fdifir) fraction m*ndim |
746 |
|
|
c all-sky direct downward par (0.4-0.7 micron) |
747 |
|
|
c flux at the surface (fdirpar) fraction m*ndim |
748 |
|
|
c all-sky diffuse downward par flux at |
749 |
|
|
c the surface (fdifpar) fraction m*ndim |
750 |
|
|
* |
751 |
|
|
c----- Notes: |
752 |
|
|
c |
753 |
|
|
c (1) The unit of flux is fraction of the incoming solar radiation |
754 |
|
|
c at the top of the atmosphere. Therefore, fluxes should |
755 |
|
|
c be equal to flux multiplied by the extra-terrestrial solar |
756 |
|
|
c flux and the cosine of solar zenith angle. |
757 |
|
|
c (2) Clouds and aerosols can be included in any layers by specifying |
758 |
|
|
c fcld(i,j,k), taucld(i,j,k,*) and taual(i,j,k), k=1,np. |
759 |
|
|
c For an atmosphere without clouds and aerosols, |
760 |
|
|
c set fcld(i,j,k)=taucld(i,j,k,*)=taual(i,j,k)=0.0. |
761 |
|
|
c (3) Aerosol single scattering albedos and asymmetry |
762 |
|
|
c factors are specified in the subroutines solir and soluv. |
763 |
|
|
c (4) pl(i,j,1) is the pressure at the top of the model, and |
764 |
|
|
c pl(i,j,np+1) is the surface pressure. |
765 |
|
|
c (5) the pressure levels ict and icb correspond approximately |
766 |
|
|
c to 400 and 700 mb. |
767 |
|
|
c |
768 |
|
|
c************************************************************************** |
769 |
|
|
|
770 |
|
|
implicit none |
771 |
|
|
|
772 |
|
|
c-----Explicit Inline Directives |
773 |
|
|
|
774 |
molod |
1.11 |
#ifdef CRAY |
775 |
|
|
#ifdef f77 |
776 |
molod |
1.1 |
cfpp$ expand (expmn) |
777 |
|
|
#endif |
778 |
|
|
#endif |
779 |
|
|
real expmn |
780 |
|
|
|
781 |
|
|
c-----input parameters |
782 |
|
|
|
783 |
|
|
integer m,n,ndim,np,ict,icb |
784 |
|
|
real pl(m,ndim,np+1),ta(m,ndim,np),wa(m,ndim,np),oa(m,ndim,np) |
785 |
|
|
real taucld(m,ndim,np,2),reff(m,ndim,np,2) |
786 |
|
|
real fcld(m,ndim,np),taual(m,ndim,np) |
787 |
|
|
real rsirbm(m,ndim),rsirdf(m,ndim), |
788 |
|
|
* rsuvbm(m,ndim),rsuvdf(m,ndim),cosz(m,ndim),co2 |
789 |
|
|
|
790 |
|
|
c-----output parameters |
791 |
|
|
|
792 |
|
|
real flx(m,ndim,np+1),flc(m,ndim,np+1) |
793 |
|
|
real fdirir(m,ndim),fdifir(m,ndim) |
794 |
|
|
real fdirpar(m,ndim),fdifpar(m,ndim) |
795 |
|
|
real fdiruv(m,ndim),fdifuv(m,ndim) |
796 |
|
|
|
797 |
|
|
c-----temporary array |
798 |
|
|
|
799 |
molod |
1.10 |
integer i,j,k |
800 |
molod |
1.1 |
real cc(m,n,3),tauclb(m,n,np),tauclf(m,n,np) |
801 |
|
|
real dp(m,n,np),wh(m,n,np),oh(m,n,np),scal(m,n,np) |
802 |
|
|
real swh(m,n,np+1),so2(m,n,np+1),df(m,n,np+1) |
803 |
molod |
1.10 |
real sdf(m,n),sclr(m,n),csm(m,n),x |
804 |
molod |
1.1 |
|
805 |
|
|
c----------------------------------------------------------------- |
806 |
|
|
|
807 |
|
|
do j= 1, n |
808 |
|
|
do i= 1, m |
809 |
|
|
|
810 |
|
|
swh(i,j,1)=0. |
811 |
|
|
so2(i,j,1)=0. |
812 |
|
|
|
813 |
|
|
c-----csm is the effective secant of the solar zenith angle |
814 |
|
|
c see equation (12) of Lacis and Hansen (1974, JAS) |
815 |
|
|
|
816 |
|
|
csm(i,j)=35./sqrt(1224.*cosz(i,j)*cosz(i,j)+1.) |
817 |
|
|
|
818 |
|
|
enddo |
819 |
|
|
enddo |
820 |
|
|
|
821 |
|
|
|
822 |
|
|
do k= 1, np |
823 |
|
|
do j= 1, n |
824 |
|
|
do i= 1, m |
825 |
|
|
|
826 |
|
|
c-----compute layer thickness and pressure-scaling function. |
827 |
|
|
c indices for the surface level and surface layer |
828 |
|
|
c are np+1 and np, respectively. |
829 |
|
|
|
830 |
|
|
dp(i,j,k)=pl(i,j,k+1)-pl(i,j,k) |
831 |
|
|
scal(i,j,k)=dp(i,j,k)*(.5*(pl(i,j,k)+pl(i,j,k+1))/300.)**.8 |
832 |
|
|
|
833 |
|
|
c-----compute scaled water vapor amount, unit is g/cm**2 |
834 |
|
|
|
835 |
|
|
wh(i,j,k)=1.02*wa(i,j,k)*scal(i,j,k)* |
836 |
|
|
* (1.+0.00135*(ta(i,j,k)-240.)) |
837 |
|
|
swh(i,j,k+1)=swh(i,j,k)+wh(i,j,k) |
838 |
|
|
|
839 |
|
|
c-----compute ozone amount, unit is (cm-atm)stp. |
840 |
|
|
|
841 |
|
|
oh(i,j,k)=1.02*oa(i,j,k)*dp(i,j,k)*466.7 |
842 |
|
|
|
843 |
|
|
enddo |
844 |
|
|
enddo |
845 |
|
|
enddo |
846 |
|
|
|
847 |
|
|
|
848 |
|
|
c-----scale cloud optical thickness in each layer from taucld (with |
849 |
|
|
c cloud amount fcld) to tauclb and tauclf (with cloud amount cc). |
850 |
|
|
c tauclb is the scaled optical thickness for beam radiation and |
851 |
|
|
c tauclf is for diffuse radiation. |
852 |
|
|
|
853 |
|
|
call cldscale(m,n,ndim,np,cosz,fcld,taucld,ict,icb, |
854 |
|
|
* cc,tauclb,tauclf) |
855 |
|
|
|
856 |
|
|
c-----initialize fluxes for all-sky (flx), clear-sky (flc), and |
857 |
|
|
c flux reduction (df) |
858 |
|
|
|
859 |
|
|
do k=1, np+1 |
860 |
|
|
do j=1, n |
861 |
|
|
do i=1, m |
862 |
|
|
flx(i,j,k)=0. |
863 |
|
|
flc(i,j,k)=0. |
864 |
|
|
df(i,j,k)=0. |
865 |
|
|
enddo |
866 |
|
|
enddo |
867 |
|
|
enddo |
868 |
|
|
|
869 |
|
|
c-----compute solar ir fluxes |
870 |
|
|
|
871 |
|
|
call solir (m,n,ndim,np,wh,taucld,tauclb,tauclf,reff,ict,icb |
872 |
|
|
* ,fcld,cc,taual,csm,rsirbm,rsirdf,flx,flc,fdirir,fdifir) |
873 |
|
|
|
874 |
|
|
c-----compute and update uv and par fluxes |
875 |
|
|
|
876 |
|
|
call soluv (m,n,ndim,np,oh,dp,taucld,tauclb,tauclf,reff,ict,icb |
877 |
|
|
* ,fcld,cc,taual,csm,rsuvbm,rsuvdf,flx,flc |
878 |
|
|
* ,fdirpar,fdifpar,fdiruv,fdifuv) |
879 |
|
|
|
880 |
|
|
c-----compute scaled amount of o2 (so2), unit is (cm-atm)stp. |
881 |
|
|
|
882 |
|
|
do k= 1, np |
883 |
|
|
do j= 1, n |
884 |
|
|
do i= 1, m |
885 |
|
|
so2(i,j,k+1)=so2(i,j,k)+165.22*scal(i,j,k) |
886 |
|
|
enddo |
887 |
|
|
enddo |
888 |
|
|
enddo |
889 |
|
|
|
890 |
|
|
c-----compute flux reduction due to oxygen following |
891 |
|
|
c chou (J. climate, 1990). The fraction 0.0287 is the |
892 |
|
|
c extraterrestrial solar flux in the o2 bands. |
893 |
|
|
|
894 |
|
|
do k= 2, np+1 |
895 |
|
|
do j= 1, n |
896 |
|
|
do i= 1, m |
897 |
|
|
x=so2(i,j,k)*csm(i,j) |
898 |
|
|
df(i,j,k)=df(i,j,k)+0.0287*(1.-expmn(-0.00027*sqrt(x))) |
899 |
|
|
enddo |
900 |
|
|
enddo |
901 |
|
|
enddo |
902 |
|
|
|
903 |
|
|
c-----compute scaled amounts for co2 (so2). unit is (cm-atm)stp. |
904 |
|
|
|
905 |
|
|
do k= 1, np |
906 |
|
|
do j= 1, n |
907 |
|
|
do i= 1, m |
908 |
|
|
so2(i,j,k+1)=so2(i,j,k)+co2*789.*scal(i,j,k) |
909 |
|
|
enddo |
910 |
|
|
enddo |
911 |
|
|
enddo |
912 |
|
|
|
913 |
|
|
c-----compute and update flux reduction due to co2 following |
914 |
|
|
c chou (J. Climate, 1990) |
915 |
|
|
|
916 |
|
|
call flxco2(m,n,np,so2,swh,csm,df) |
917 |
|
|
|
918 |
|
|
c-----adjust for the effect of o2 cnd co2 on clear-sky fluxes. |
919 |
|
|
|
920 |
|
|
do k= 2, np+1 |
921 |
|
|
do j= 1, n |
922 |
|
|
do i= 1, m |
923 |
|
|
flc(i,j,k)=flc(i,j,k)-df(i,j,k) |
924 |
|
|
enddo |
925 |
|
|
enddo |
926 |
|
|
enddo |
927 |
|
|
|
928 |
|
|
c-----adjust for the all-sky fluxes due to o2 and co2. It is |
929 |
|
|
c assumed that o2 and co2 have no effects on solar radiation |
930 |
|
|
c below clouds. clouds are assumed randomly overlapped. |
931 |
|
|
|
932 |
|
|
do j=1,n |
933 |
|
|
do i=1,m |
934 |
|
|
sdf(i,j)=0.0 |
935 |
|
|
sclr(i,j)=1.0 |
936 |
|
|
enddo |
937 |
|
|
enddo |
938 |
|
|
|
939 |
|
|
do k=1,np |
940 |
|
|
do j=1,n |
941 |
|
|
do i=1,m |
942 |
|
|
|
943 |
|
|
c-----sclr is the fraction of clear sky. |
944 |
|
|
c sdf is the flux reduction below clouds. |
945 |
|
|
|
946 |
|
|
if(fcld(i,j,k).gt.0.01) then |
947 |
|
|
sdf(i,j)=sdf(i,j)+df(i,j,k)*sclr(i,j)*fcld(i,j,k) |
948 |
|
|
sclr(i,j)=sclr(i,j)*(1.-fcld(i,j,k)) |
949 |
|
|
endif |
950 |
|
|
flx(i,j,k+1)=flx(i,j,k+1)-sdf(i,j)-df(i,j,k+1)*sclr(i,j) |
951 |
|
|
|
952 |
|
|
enddo |
953 |
|
|
enddo |
954 |
|
|
enddo |
955 |
|
|
|
956 |
|
|
c-----adjust for the direct downward ir flux. |
957 |
|
|
do j= 1, n |
958 |
|
|
do i= 1, m |
959 |
|
|
x = (1.-rsirdf(i,j))*fdifir(i,j) + (1.-rsirbm(i,j))*fdirir(i,j) |
960 |
|
|
x = (-sdf(i,j)-df(i,j,np+1)*sclr(i,j))/x |
961 |
|
|
fdirir(i,j)=fdirir(i,j)*(1.+x) |
962 |
|
|
fdifir(i,j)=fdifir(i,j)*(1.+x) |
963 |
|
|
enddo |
964 |
|
|
enddo |
965 |
|
|
|
966 |
|
|
return |
967 |
|
|
end |
968 |
|
|
|
969 |
|
|
c******************************************************************** |
970 |
|
|
|
971 |
|
|
subroutine cldscale (m,n,ndim,np,cosz,fcld,taucld,ict,icb, |
972 |
|
|
* cc,tauclb,tauclf) |
973 |
|
|
|
974 |
|
|
c******************************************************************** |
975 |
|
|
c |
976 |
|
|
c This subroutine computes the covers of high, middle, and |
977 |
|
|
c low clouds and scales the cloud optical thickness. |
978 |
|
|
c |
979 |
|
|
c To simplify calculations in a cloudy atmosphere, clouds are |
980 |
|
|
c grouped into high, middle and low clouds separated by the levels |
981 |
|
|
c ict and icb (level 1 is the top of the atmosphere). |
982 |
|
|
c |
983 |
|
|
c Within each of the three groups, clouds are assumed maximally |
984 |
|
|
c overlapped, and the cloud cover (cc) of a group is the maximum |
985 |
|
|
c cloud cover of all the layers in the group. The optical thickness |
986 |
|
|
c (taucld) of a given layer is then scaled to new values (tauclb and |
987 |
|
|
c tauclf) so that the layer reflectance corresponding to the cloud |
988 |
|
|
c cover cc is the same as the original reflectance with optical |
989 |
|
|
c thickness taucld and cloud cover fcld. |
990 |
|
|
c |
991 |
|
|
c---input parameters |
992 |
|
|
c |
993 |
|
|
c number of grid intervals in zonal direction (m) |
994 |
|
|
c number of grid intervals in meridional direction (n) |
995 |
|
|
c maximum number of grid intervals in meridional direction (ndim) |
996 |
|
|
c number of atmospheric layers (np) |
997 |
|
|
c cosine of the solar zenith angle (cosz) |
998 |
|
|
c fractional cloud cover (fcld) |
999 |
|
|
c cloud optical thickness (taucld) |
1000 |
|
|
c index separating high and middle clouds (ict) |
1001 |
|
|
c index separating middle and low clouds (icb) |
1002 |
|
|
c |
1003 |
|
|
c---output parameters |
1004 |
|
|
c |
1005 |
|
|
c fractional cover of high, middle, and low clouds (cc) |
1006 |
|
|
c scaled cloud optical thickness for beam radiation (tauclb) |
1007 |
|
|
c scaled cloud optical thickness for diffuse radiation (tauclf) |
1008 |
|
|
c |
1009 |
|
|
c******************************************************************** |
1010 |
|
|
|
1011 |
|
|
implicit none |
1012 |
|
|
|
1013 |
|
|
c-----input parameters |
1014 |
|
|
|
1015 |
|
|
integer m,n,ndim,np,ict,icb |
1016 |
|
|
real cosz(m,ndim),fcld(m,ndim,np),taucld(m,ndim,np,2) |
1017 |
|
|
|
1018 |
|
|
c-----output parameters |
1019 |
|
|
|
1020 |
|
|
real cc(m,n,3),tauclb(m,n,np),tauclf(m,n,np) |
1021 |
|
|
|
1022 |
|
|
c-----temporary variables |
1023 |
|
|
|
1024 |
|
|
integer i,j,k,im,it,ia,kk |
1025 |
molod |
1.10 |
real fm,ft,fa,xai,taux |
1026 |
molod |
1.1 |
|
1027 |
|
|
c-----pre-computed table |
1028 |
|
|
|
1029 |
|
|
integer nm,nt,na |
1030 |
|
|
parameter (nm=11,nt=9,na=11) |
1031 |
|
|
real dm,dt,da,t1,caib(nm,nt,na),caif(nt,na) |
1032 |
|
|
parameter (dm=0.1,dt=0.30103,da=0.1,t1=-0.9031) |
1033 |
|
|
|
1034 |
|
|
c-----include the pre-computed table for cai |
1035 |
|
|
|
1036 |
molod |
1.8 |
#include "cai-dat.h" |
1037 |
molod |
1.9 |
c save caib,caif |
1038 |
molod |
1.1 |
|
1039 |
|
|
|
1040 |
|
|
c-----clouds within each of the high, middle, and low clouds are |
1041 |
|
|
c assumed maximally overlapped, and the cloud cover (cc) |
1042 |
|
|
c for a group is the maximum cloud cover of all the layers |
1043 |
|
|
c in the group |
1044 |
|
|
|
1045 |
|
|
do j=1,n |
1046 |
|
|
do i=1,m |
1047 |
|
|
cc(i,j,1)=0.0 |
1048 |
|
|
cc(i,j,2)=0.0 |
1049 |
|
|
cc(i,j,3)=0.0 |
1050 |
|
|
enddo |
1051 |
|
|
enddo |
1052 |
|
|
|
1053 |
|
|
do k=1,ict-1 |
1054 |
|
|
do j=1,n |
1055 |
|
|
do i=1,m |
1056 |
|
|
cc(i,j,1)=max(cc(i,j,1),fcld(i,j,k)) |
1057 |
|
|
enddo |
1058 |
|
|
enddo |
1059 |
|
|
enddo |
1060 |
|
|
|
1061 |
|
|
do k=ict,icb-1 |
1062 |
|
|
do j=1,n |
1063 |
|
|
do i=1,m |
1064 |
|
|
cc(i,j,2)=max(cc(i,j,2),fcld(i,j,k)) |
1065 |
|
|
enddo |
1066 |
|
|
enddo |
1067 |
|
|
enddo |
1068 |
|
|
|
1069 |
|
|
do k=icb,np |
1070 |
|
|
do j=1,n |
1071 |
|
|
do i=1,m |
1072 |
|
|
cc(i,j,3)=max(cc(i,j,3),fcld(i,j,k)) |
1073 |
|
|
enddo |
1074 |
|
|
enddo |
1075 |
|
|
enddo |
1076 |
|
|
|
1077 |
|
|
c-----scale the cloud optical thickness. |
1078 |
|
|
c taucld(i,j,k,1) is the optical thickness for ice particles, and |
1079 |
|
|
c taucld(i,j,k,2) is the optical thickness for liquid particles. |
1080 |
|
|
|
1081 |
|
|
do k=1,np |
1082 |
|
|
|
1083 |
|
|
if(k.lt.ict) then |
1084 |
|
|
kk=1 |
1085 |
|
|
elseif(k.ge.ict .and. k.lt.icb) then |
1086 |
|
|
kk=2 |
1087 |
|
|
else |
1088 |
|
|
kk=3 |
1089 |
|
|
endif |
1090 |
|
|
|
1091 |
|
|
do j=1,n |
1092 |
|
|
do i=1,m |
1093 |
|
|
|
1094 |
|
|
tauclb(i,j,k) = 0.0 |
1095 |
|
|
tauclf(i,j,k) = 0.0 |
1096 |
|
|
|
1097 |
|
|
taux=taucld(i,j,k,1)+taucld(i,j,k,2) |
1098 |
|
|
if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then |
1099 |
|
|
|
1100 |
|
|
c-----normalize cloud cover |
1101 |
|
|
|
1102 |
|
|
fa=fcld(i,j,k)/cc(i,j,kk) |
1103 |
|
|
|
1104 |
|
|
c-----table look-up |
1105 |
|
|
|
1106 |
|
|
taux=min(taux,32.) |
1107 |
|
|
|
1108 |
|
|
fm=cosz(i,j)/dm |
1109 |
|
|
ft=(log10(taux)-t1)/dt |
1110 |
|
|
fa=fa/da |
1111 |
|
|
|
1112 |
|
|
im=int(fm+1.5) |
1113 |
|
|
it=int(ft+1.5) |
1114 |
|
|
ia=int(fa+1.5) |
1115 |
|
|
|
1116 |
|
|
im=max(im,2) |
1117 |
|
|
it=max(it,2) |
1118 |
|
|
ia=max(ia,2) |
1119 |
|
|
|
1120 |
|
|
im=min(im,nm-1) |
1121 |
|
|
it=min(it,nt-1) |
1122 |
|
|
ia=min(ia,na-1) |
1123 |
|
|
|
1124 |
|
|
fm=fm-float(im-1) |
1125 |
|
|
ft=ft-float(it-1) |
1126 |
|
|
fa=fa-float(ia-1) |
1127 |
|
|
|
1128 |
|
|
c-----scale cloud optical thickness for beam radiation. |
1129 |
|
|
c the scaling factor, xai, is a function of the solar zenith |
1130 |
|
|
c angle, optical thickness, and cloud cover. |
1131 |
|
|
|
1132 |
|
|
xai= (-caib(im-1,it,ia)*(1.-fm)+ |
1133 |
|
|
* caib(im+1,it,ia)*(1.+fm))*fm*.5+caib(im,it,ia)*(1.-fm*fm) |
1134 |
|
|
|
1135 |
|
|
xai=xai+(-caib(im,it-1,ia)*(1.-ft)+ |
1136 |
|
|
* caib(im,it+1,ia)*(1.+ft))*ft*.5+caib(im,it,ia)*(1.-ft*ft) |
1137 |
|
|
|
1138 |
|
|
xai=xai+(-caib(im,it,ia-1)*(1.-fa)+ |
1139 |
|
|
* caib(im,it,ia+1)*(1.+fa))*fa*.5+caib(im,it,ia)*(1.-fa*fa) |
1140 |
|
|
|
1141 |
|
|
xai= xai-2.*caib(im,it,ia) |
1142 |
|
|
xai=max(xai,0.0) |
1143 |
|
|
|
1144 |
|
|
tauclb(i,j,k) = taux*xai |
1145 |
|
|
|
1146 |
|
|
c-----scale cloud optical thickness for diffuse radiation. |
1147 |
|
|
c the scaling factor, xai, is a function of the cloud optical |
1148 |
|
|
c thickness and cover but not the solar zenith angle. |
1149 |
|
|
|
1150 |
|
|
xai= (-caif(it-1,ia)*(1.-ft)+ |
1151 |
|
|
* caif(it+1,ia)*(1.+ft))*ft*.5+caif(it,ia)*(1.-ft*ft) |
1152 |
|
|
|
1153 |
|
|
xai=xai+(-caif(it,ia-1)*(1.-fa)+ |
1154 |
|
|
* caif(it,ia+1)*(1.+fa))*fa*.5+caif(it,ia)*(1.-fa*fa) |
1155 |
|
|
|
1156 |
|
|
xai= xai-caif(it,ia) |
1157 |
|
|
xai=max(xai,0.0) |
1158 |
|
|
|
1159 |
|
|
tauclf(i,j,k) = taux*xai |
1160 |
|
|
|
1161 |
|
|
endif |
1162 |
|
|
|
1163 |
|
|
enddo |
1164 |
|
|
enddo |
1165 |
|
|
enddo |
1166 |
|
|
|
1167 |
|
|
return |
1168 |
|
|
end |
1169 |
|
|
c*********************************************************************** |
1170 |
|
|
|
1171 |
|
|
subroutine solir (m,n,ndim,np,wh,taucld,tauclb,tauclf,reff, |
1172 |
|
|
* ict,icb,fcld,cc,taual,csm,rsirbm,rsirdf,flx,flc,fdirir,fdifir) |
1173 |
|
|
|
1174 |
|
|
c************************************************************************ |
1175 |
|
|
c compute solar flux in the infrared region. The spectrum is divided |
1176 |
|
|
c into three bands: |
1177 |
|
|
c |
1178 |
|
|
c band wavenumber(/cm) wavelength (micron) |
1179 |
|
|
c 1 1000-4400 2.27-10.0 |
1180 |
|
|
c 2 4400-8200 1.22-2.27 |
1181 |
|
|
c 3 8200-14300 0.70-1.22 |
1182 |
|
|
c |
1183 |
|
|
c----- Input parameters: units size |
1184 |
|
|
c |
1185 |
|
|
c number of soundings in zonal direction (m) n/d 1 |
1186 |
|
|
c number of soundings in meridional direction (n) n/d 1 |
1187 |
|
|
c maximum number of soundings in n/d 1 |
1188 |
|
|
c meridional direction (ndim) |
1189 |
|
|
c number of atmospheric layers (np) n/d 1 |
1190 |
|
|
c layer water vapor content (wh) gm/cm^2 m*n*np |
1191 |
|
|
c cloud optical thickness (taucld) n/d m*ndim*np*2 |
1192 |
|
|
c index 1 for ice paticles |
1193 |
|
|
c index 2 for liquid particles |
1194 |
|
|
c scaled cloud optical thickness n/d m*n*np |
1195 |
|
|
c for beam radiation (tauclb) |
1196 |
|
|
c scaled cloud optical thickness n/d m*n*np |
1197 |
|
|
c for diffuse radiation (tauclf) |
1198 |
|
|
c effective cloud-particle size (reff) micrometer m*ndim*np*2 |
1199 |
|
|
c index 1 for ice paticles |
1200 |
|
|
c index 2 for liquid particles |
1201 |
|
|
c level index separating high and n/d m*n |
1202 |
|
|
c middle clouds (ict) |
1203 |
|
|
c level index separating middle and n/d m*n |
1204 |
|
|
c low clouds (icb) |
1205 |
|
|
c cloud amount (fcld) fraction m*ndim*np |
1206 |
|
|
c cloud amount of high, middle, and n/d m*n*3 |
1207 |
|
|
c low clouds (cc) |
1208 |
|
|
c aerosol optical thickness (taual) n/d m*ndim*np |
1209 |
|
|
c cosecant of the solar zenith angle (csm) n/d m*n |
1210 |
|
|
c near ir surface albedo for beam fraction m*ndim |
1211 |
|
|
c radiation (rsirbm) |
1212 |
|
|
c near ir surface albedo for diffuse fraction m*ndim |
1213 |
|
|
c radiation (rsirdf) |
1214 |
|
|
c |
1215 |
|
|
c----- output (updated) parameters: |
1216 |
|
|
c |
1217 |
|
|
c all-sky flux (downward-upward) (flx) fraction m*ndim*(np+1) |
1218 |
|
|
c clear-sky flux (downward-upward) (flc) fraction m*ndim*(np+1) |
1219 |
|
|
c all-sky direct downward ir flux at |
1220 |
|
|
c the surface (fdirir) fraction m*ndim |
1221 |
|
|
c all-sky diffuse downward ir flux at |
1222 |
|
|
c the surface (fdifir) fraction m*ndim |
1223 |
|
|
c |
1224 |
|
|
c----- note: the following parameters must be specified by users: |
1225 |
|
|
c aerosol single scattering albedo (ssaal) n/d nband |
1226 |
|
|
c aerosol asymmetry factor (asyal) n/d nband |
1227 |
|
|
c |
1228 |
|
|
c************************************************************************* |
1229 |
|
|
|
1230 |
|
|
implicit none |
1231 |
|
|
|
1232 |
|
|
c-----Explicit Inline Directives |
1233 |
|
|
|
1234 |
molod |
1.11 |
#ifdef CRAY |
1235 |
|
|
#ifdef f77 |
1236 |
molod |
1.1 |
cfpp$ expand (deledd) |
1237 |
|
|
cfpp$ expand (sagpol) |
1238 |
|
|
cfpp$ expand (expmn) |
1239 |
|
|
#endif |
1240 |
|
|
#endif |
1241 |
|
|
real expmn |
1242 |
|
|
|
1243 |
|
|
c-----input parameters |
1244 |
|
|
|
1245 |
|
|
integer m,n,ndim,np,ict,icb |
1246 |
|
|
real taucld(m,ndim,np,2),reff(m,ndim,np,2),fcld(m,ndim,np) |
1247 |
|
|
real tauclb(m,n,np),tauclf(m,n,np),cc(m,n,3) |
1248 |
|
|
real rsirbm(m,ndim),rsirdf(m,ndim) |
1249 |
|
|
real wh(m,n,np),taual(m,ndim,np),csm(m,n) |
1250 |
|
|
|
1251 |
|
|
c-----output (updated) parameters |
1252 |
|
|
|
1253 |
|
|
real flx(m,ndim,np+1),flc(m,ndim,np+1) |
1254 |
|
|
real fdirir(m,ndim),fdifir(m,ndim) |
1255 |
|
|
|
1256 |
|
|
c-----static parameters |
1257 |
|
|
|
1258 |
|
|
integer nk,nband |
1259 |
|
|
parameter (nk=10,nband=3) |
1260 |
|
|
real xk(nk),hk(nband,nk),ssaal(nband),asyal(nband) |
1261 |
|
|
real aia(nband,3),awa(nband,3),aig(nband,3),awg(nband,3) |
1262 |
|
|
|
1263 |
|
|
c-----temporary array |
1264 |
|
|
|
1265 |
|
|
integer ib,ik,i,j,k |
1266 |
|
|
real ssacl(m,n,np),asycl(m,n,np) |
1267 |
|
|
real rr(m,n,np+1,2),tt(m,n,np+1,2),td(m,n,np+1,2), |
1268 |
|
|
* rs(m,n,np+1,2),ts(m,n,np+1,2) |
1269 |
|
|
real fall(m,n,np+1),fclr(m,n,np+1) |
1270 |
|
|
real fsdir(m,n),fsdif(m,n) |
1271 |
|
|
|
1272 |
|
|
real tauwv,tausto,ssatau,asysto,tauto,ssato,asyto |
1273 |
|
|
real taux,reff1,reff2,w1,w2,g1,g2 |
1274 |
|
|
real ssaclt(m,n),asyclt(m,n) |
1275 |
|
|
real rr1t(m,n),tt1t(m,n),td1t(m,n),rs1t(m,n),ts1t(m,n) |
1276 |
|
|
real rr2t(m,n),tt2t(m,n),td2t(m,n),rs2t(m,n),ts2t(m,n) |
1277 |
|
|
|
1278 |
|
|
c-----water vapor absorption coefficient for 10 k-intervals. |
1279 |
|
|
c unit: cm^2/gm |
1280 |
|
|
|
1281 |
|
|
data xk/ |
1282 |
|
|
1 0.0010, 0.0133, 0.0422, 0.1334, 0.4217, |
1283 |
|
|
2 1.334, 5.623, 31.62, 177.8, 1000.0/ |
1284 |
|
|
|
1285 |
|
|
c-----water vapor k-distribution function, |
1286 |
|
|
c the sum of hk is 0.52926. unit: fraction |
1287 |
|
|
|
1288 |
|
|
data hk/ |
1289 |
|
|
1 .01074,.08236,.20673, .00360,.01157,.03497, |
1290 |
|
|
2 .00411,.01133,.03011, .00421,.01143,.02260, |
1291 |
|
|
3 .00389,.01240,.01336, .00326,.01258,.00696, |
1292 |
|
|
4 .00499,.01381,.00441, .00465,.00650,.00115, |
1293 |
|
|
5 .00245,.00244,.00026, .00145,.00094,.00000/ |
1294 |
|
|
|
1295 |
|
|
c-----aerosol single-scattering albedo and asymmetry factor |
1296 |
|
|
|
1297 |
|
|
data ssaal,asyal/nband*0.999,nband*0.850/ |
1298 |
|
|
|
1299 |
|
|
c-----coefficients for computing the single scattering albedo of |
1300 |
|
|
c ice clouds from ssa=1-(aia(*,1)+aia(*,2)*reff+aia(*,3)*reff**2) |
1301 |
|
|
|
1302 |
|
|
data aia/ |
1303 |
|
|
1 .08938331, .00215346,-.00000260, |
1304 |
|
|
2 .00299387, .00073709, .00000746, |
1305 |
|
|
3 -.00001038,-.00000134, .00000000/ |
1306 |
|
|
|
1307 |
|
|
c-----coefficients for computing the single scattering albedo of |
1308 |
|
|
c liquid clouds from ssa=1-(awa(*,1)+awa(*,2)*reff+awa(*,3)*reff**2) |
1309 |
|
|
|
1310 |
|
|
data awa/ |
1311 |
|
|
1 .01209318,-.00019934, .00000007, |
1312 |
|
|
2 .01784739, .00088757, .00000845, |
1313 |
|
|
3 -.00036910,-.00000650,-.00000004/ |
1314 |
|
|
|
1315 |
|
|
c-----coefficients for computing the asymmetry factor of ice clouds |
1316 |
|
|
c from asycl=aig(*,1)+aig(*,2)*reff+aig(*,3)*reff**2 |
1317 |
|
|
|
1318 |
|
|
data aig/ |
1319 |
|
|
1 .84090400, .76098937, .74935228, |
1320 |
|
|
2 .00126222, .00141864, .00119715, |
1321 |
|
|
3 -.00000385,-.00000396,-.00000367/ |
1322 |
|
|
|
1323 |
|
|
c-----coefficients for computing the asymmetry factor of liquid clouds |
1324 |
|
|
c from asycl=awg(*,1)+awg(*,2)*reff+awg(*,3)*reff**2 |
1325 |
|
|
|
1326 |
|
|
data awg/ |
1327 |
|
|
1 .83530748, .74513197, .79375035, |
1328 |
|
|
2 .00257181, .01370071, .00832441, |
1329 |
|
|
3 .00005519,-.00038203,-.00023263/ |
1330 |
|
|
|
1331 |
|
|
c-----initialize surface fluxes, reflectances, and transmittances |
1332 |
|
|
|
1333 |
|
|
do j= 1, n |
1334 |
|
|
do i= 1, m |
1335 |
|
|
fdirir(i,j)=0.0 |
1336 |
|
|
fdifir(i,j)=0.0 |
1337 |
|
|
rr(i,j,np+1,1)=rsirbm(i,j) |
1338 |
|
|
rr(i,j,np+1,2)=rsirbm(i,j) |
1339 |
|
|
rs(i,j,np+1,1)=rsirdf(i,j) |
1340 |
|
|
rs(i,j,np+1,2)=rsirdf(i,j) |
1341 |
|
|
td(i,j,np+1,1)=0.0 |
1342 |
|
|
td(i,j,np+1,2)=0.0 |
1343 |
|
|
tt(i,j,np+1,1)=0.0 |
1344 |
|
|
tt(i,j,np+1,2)=0.0 |
1345 |
|
|
ts(i,j,np+1,1)=0.0 |
1346 |
|
|
ts(i,j,np+1,2)=0.0 |
1347 |
|
|
enddo |
1348 |
|
|
enddo |
1349 |
|
|
|
1350 |
|
|
c-----integration over spectral bands |
1351 |
|
|
|
1352 |
|
|
do 100 ib=1,nband |
1353 |
|
|
|
1354 |
|
|
c-----compute cloud single scattering albedo and asymmetry factor |
1355 |
|
|
c for a mixture of ice and liquid particles. |
1356 |
|
|
|
1357 |
|
|
do k= 1, np |
1358 |
|
|
|
1359 |
|
|
do j= 1, n |
1360 |
|
|
do i= 1, m |
1361 |
|
|
|
1362 |
|
|
ssaclt(i,j)=1.0 |
1363 |
|
|
asyclt(i,j)=1.0 |
1364 |
|
|
|
1365 |
|
|
taux=taucld(i,j,k,1)+taucld(i,j,k,2) |
1366 |
|
|
if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then |
1367 |
|
|
|
1368 |
|
|
reff1=min(reff(i,j,k,1),130.) |
1369 |
|
|
reff2=min(reff(i,j,k,2),20.0) |
1370 |
|
|
|
1371 |
|
|
w1=(1.-(aia(ib,1)+(aia(ib,2)+ |
1372 |
|
|
* aia(ib,3)*reff1)*reff1))*taucld(i,j,k,1) |
1373 |
|
|
w2=(1.-(awa(ib,1)+(awa(ib,2)+ |
1374 |
|
|
* awa(ib,3)*reff2)*reff2))*taucld(i,j,k,2) |
1375 |
|
|
ssaclt(i,j)=(w1+w2)/taux |
1376 |
|
|
|
1377 |
|
|
g1=(aig(ib,1)+(aig(ib,2)+aig(ib,3)*reff1)*reff1)*w1 |
1378 |
|
|
g2=(awg(ib,1)+(awg(ib,2)+awg(ib,3)*reff2)*reff2)*w2 |
1379 |
|
|
asyclt(i,j)=(g1+g2)/(w1+w2) |
1380 |
|
|
|
1381 |
|
|
endif |
1382 |
|
|
|
1383 |
|
|
enddo |
1384 |
|
|
enddo |
1385 |
|
|
|
1386 |
|
|
do j=1,n |
1387 |
|
|
do i=1,m |
1388 |
|
|
ssacl(i,j,k)=ssaclt(i,j) |
1389 |
|
|
enddo |
1390 |
|
|
enddo |
1391 |
|
|
do j=1,n |
1392 |
|
|
do i=1,m |
1393 |
|
|
asycl(i,j,k)=asyclt(i,j) |
1394 |
|
|
enddo |
1395 |
|
|
enddo |
1396 |
|
|
|
1397 |
|
|
enddo |
1398 |
|
|
|
1399 |
|
|
c-----integration over the k-distribution function |
1400 |
|
|
|
1401 |
|
|
do 200 ik=1,nk |
1402 |
|
|
|
1403 |
|
|
do 300 k= 1, np |
1404 |
|
|
|
1405 |
|
|
do j= 1, n |
1406 |
|
|
do i= 1, m |
1407 |
|
|
|
1408 |
|
|
tauwv=xk(ik)*wh(i,j,k) |
1409 |
|
|
|
1410 |
|
|
c-----compute total optical thickness, single scattering albedo, |
1411 |
|
|
c and asymmetry factor. |
1412 |
|
|
|
1413 |
|
|
tausto=tauwv+taual(i,j,k)+1.0e-8 |
1414 |
|
|
ssatau=ssaal(ib)*taual(i,j,k) |
1415 |
|
|
asysto=asyal(ib)*ssaal(ib)*taual(i,j,k) |
1416 |
|
|
|
1417 |
|
|
c-----compute reflectance and transmittance for cloudless layers |
1418 |
|
|
|
1419 |
|
|
tauto=tausto |
1420 |
|
|
ssato=ssatau/tauto+1.0e-8 |
1421 |
|
|
|
1422 |
|
|
c if (ssato .gt. 0.001) then |
1423 |
|
|
|
1424 |
|
|
c ssato=min(ssato,0.999999) |
1425 |
|
|
c asyto=asysto/(ssato*tauto) |
1426 |
|
|
|
1427 |
|
|
c call deledd(tauto,ssato,asyto,csm(i,j), |
1428 |
|
|
c * rr1t(i,j),tt1t(i,j),td1t(i,j)) |
1429 |
|
|
|
1430 |
|
|
c call sagpol (tauto,ssato,asyto,rs1t(i,j),ts1t(i,j)) |
1431 |
|
|
|
1432 |
|
|
c else |
1433 |
|
|
|
1434 |
|
|
td1t(i,j)=expmn(-tauto*csm(i,j)) |
1435 |
|
|
ts1t(i,j)=expmn(-1.66*tauto) |
1436 |
|
|
tt1t(i,j)=0.0 |
1437 |
|
|
rr1t(i,j)=0.0 |
1438 |
|
|
rs1t(i,j)=0.0 |
1439 |
|
|
|
1440 |
|
|
c endif |
1441 |
|
|
|
1442 |
|
|
c-----compute reflectance and transmittance for cloud layers |
1443 |
|
|
|
1444 |
|
|
if (tauclb(i,j,k) .lt. 0.01) then |
1445 |
|
|
|
1446 |
|
|
rr2t(i,j)=rr1t(i,j) |
1447 |
|
|
tt2t(i,j)=tt1t(i,j) |
1448 |
|
|
td2t(i,j)=td1t(i,j) |
1449 |
|
|
rs2t(i,j)=rs1t(i,j) |
1450 |
|
|
ts2t(i,j)=ts1t(i,j) |
1451 |
|
|
|
1452 |
|
|
else |
1453 |
|
|
|
1454 |
|
|
tauto=tausto+tauclb(i,j,k) |
1455 |
|
|
ssato=(ssatau+ssacl(i,j,k)*tauclb(i,j,k))/tauto+1.0e-8 |
1456 |
|
|
ssato=min(ssato,0.999999) |
1457 |
|
|
asyto=(asysto+asycl(i,j,k)*ssacl(i,j,k)*tauclb(i,j,k))/ |
1458 |
|
|
* (ssato*tauto) |
1459 |
|
|
|
1460 |
|
|
call deledd(tauto,ssato,asyto,csm(i,j), |
1461 |
|
|
* rr2t(i,j),tt2t(i,j),td2t(i,j)) |
1462 |
|
|
|
1463 |
|
|
tauto=tausto+tauclf(i,j,k) |
1464 |
|
|
ssato=(ssatau+ssacl(i,j,k)*tauclf(i,j,k))/tauto+1.0e-8 |
1465 |
|
|
ssato=min(ssato,0.999999) |
1466 |
|
|
asyto=(asysto+asycl(i,j,k)*ssacl(i,j,k)*tauclf(i,j,k))/ |
1467 |
|
|
* (ssato*tauto) |
1468 |
|
|
|
1469 |
|
|
call sagpol (tauto,ssato,asyto,rs2t(i,j),ts2t(i,j)) |
1470 |
|
|
|
1471 |
|
|
endif |
1472 |
|
|
|
1473 |
|
|
enddo |
1474 |
|
|
enddo |
1475 |
|
|
|
1476 |
|
|
c-----the following code appears to be awkward, but it is efficient |
1477 |
|
|
c in certain parallel processors. |
1478 |
|
|
|
1479 |
|
|
do j=1,n |
1480 |
|
|
do i=1,m |
1481 |
|
|
rr(i,j,k,1)=rr1t(i,j) |
1482 |
|
|
enddo |
1483 |
|
|
enddo |
1484 |
|
|
do j=1,n |
1485 |
|
|
do i=1,m |
1486 |
|
|
tt(i,j,k,1)=tt1t(i,j) |
1487 |
|
|
enddo |
1488 |
|
|
enddo |
1489 |
|
|
do j=1,n |
1490 |
|
|
do i=1,m |
1491 |
|
|
td(i,j,k,1)=td1t(i,j) |
1492 |
|
|
enddo |
1493 |
|
|
enddo |
1494 |
|
|
do j=1,n |
1495 |
|
|
do i=1,m |
1496 |
|
|
rs(i,j,k,1)=rs1t(i,j) |
1497 |
|
|
enddo |
1498 |
|
|
enddo |
1499 |
|
|
do j=1,n |
1500 |
|
|
do i=1,m |
1501 |
|
|
ts(i,j,k,1)=ts1t(i,j) |
1502 |
|
|
enddo |
1503 |
|
|
enddo |
1504 |
|
|
|
1505 |
|
|
do j=1,n |
1506 |
|
|
do i=1,m |
1507 |
|
|
rr(i,j,k,2)=rr2t(i,j) |
1508 |
|
|
enddo |
1509 |
|
|
enddo |
1510 |
|
|
do j=1,n |
1511 |
|
|
do i=1,m |
1512 |
|
|
tt(i,j,k,2)=tt2t(i,j) |
1513 |
|
|
enddo |
1514 |
|
|
enddo |
1515 |
|
|
do j=1,n |
1516 |
|
|
do i=1,m |
1517 |
|
|
td(i,j,k,2)=td2t(i,j) |
1518 |
|
|
enddo |
1519 |
|
|
enddo |
1520 |
|
|
do j=1,n |
1521 |
|
|
do i=1,m |
1522 |
|
|
rs(i,j,k,2)=rs2t(i,j) |
1523 |
|
|
enddo |
1524 |
|
|
enddo |
1525 |
|
|
do j=1,n |
1526 |
|
|
do i=1,m |
1527 |
|
|
ts(i,j,k,2)=ts2t(i,j) |
1528 |
|
|
enddo |
1529 |
|
|
enddo |
1530 |
|
|
|
1531 |
|
|
300 continue |
1532 |
|
|
|
1533 |
|
|
c-----flux calculations |
1534 |
|
|
|
1535 |
|
|
call cldflx (m,n,np,ict,icb,cc,rr,tt,td,rs,ts, |
1536 |
|
|
* fclr,fall,fsdir,fsdif) |
1537 |
|
|
|
1538 |
|
|
do k= 1, np+1 |
1539 |
|
|
do j= 1, n |
1540 |
|
|
do i= 1, m |
1541 |
|
|
flx(i,j,k) = flx(i,j,k)+fall(i,j,k)*hk(ib,ik) |
1542 |
|
|
enddo |
1543 |
|
|
enddo |
1544 |
|
|
do j= 1, n |
1545 |
|
|
do i= 1, m |
1546 |
|
|
flc(i,j,k) = flc(i,j,k)+fclr(i,j,k)*hk(ib,ik) |
1547 |
|
|
enddo |
1548 |
|
|
enddo |
1549 |
|
|
enddo |
1550 |
|
|
|
1551 |
|
|
do j= 1, n |
1552 |
|
|
do i= 1, m |
1553 |
|
|
fdirir(i,j) = fdirir(i,j)+fsdir(i,j)*hk(ib,ik) |
1554 |
|
|
fdifir(i,j) = fdifir(i,j)+fsdif(i,j)*hk(ib,ik) |
1555 |
|
|
enddo |
1556 |
|
|
enddo |
1557 |
|
|
|
1558 |
|
|
200 continue |
1559 |
|
|
100 continue |
1560 |
|
|
|
1561 |
|
|
return |
1562 |
|
|
end |
1563 |
|
|
|
1564 |
|
|
c************************************************************************ |
1565 |
|
|
|
1566 |
|
|
subroutine soluv (m,n,ndim,np,oh,dp,taucld,tauclb,tauclf,reff, |
1567 |
|
|
* ict,icb,fcld,cc,taual,csm,rsuvbm,rsuvdf,flx,flc |
1568 |
|
|
* ,fdirpar,fdifpar,fdiruv,fdifuv) |
1569 |
|
|
|
1570 |
|
|
c************************************************************************ |
1571 |
|
|
c compute solar fluxes in the uv+visible region. the spectrum is |
1572 |
|
|
c grouped into 8 bands: |
1573 |
|
|
c |
1574 |
|
|
c Band Micrometer |
1575 |
|
|
c |
1576 |
|
|
c UV-C 1. .175 - .225 |
1577 |
|
|
c 2. .225 - .245 |
1578 |
|
|
c .260 - .280 |
1579 |
|
|
c 3. .245 - .260 |
1580 |
|
|
c |
1581 |
|
|
c UV-B 4. .280 - .295 |
1582 |
|
|
c 5. .295 - .310 |
1583 |
|
|
c 6. .310 - .320 |
1584 |
|
|
c |
1585 |
|
|
c UV-A 7. .320 - .400 |
1586 |
|
|
c |
1587 |
|
|
c PAR 8. .400 - .700 |
1588 |
|
|
c |
1589 |
|
|
c----- Input parameters: units size |
1590 |
|
|
c |
1591 |
|
|
c number of soundings in zonal direction (m) n/d 1 |
1592 |
|
|
c number of soundings in meridional direction (n) n/d 1 |
1593 |
|
|
c maximum number of soundings in n/d 1 |
1594 |
|
|
c meridional direction (ndim) |
1595 |
|
|
c number of atmospheric layers (np) n/d 1 |
1596 |
|
|
c layer ozone content (oh) (cm-atm)stp m*n*np |
1597 |
|
|
c layer pressure thickness (dp) mb m*n*np |
1598 |
|
|
c cloud optical thickness (taucld) n/d m*ndim*np*2 |
1599 |
|
|
c index 1 for ice paticles |
1600 |
|
|
c index 2 for liquid particles |
1601 |
|
|
c scaled cloud optical thickness n/d m*n*np |
1602 |
|
|
c for beam radiation (tauclb) |
1603 |
|
|
c scaled cloud optical thickness n/d m*n*np |
1604 |
|
|
c for diffuse radiation (tauclf) |
1605 |
|
|
c effective cloud-particle size (reff) micrometer m*ndim*np*2 |
1606 |
|
|
c index 1 for ice paticles |
1607 |
|
|
c index 2 for liquid particles |
1608 |
|
|
c level indiex separating high and n/d m*n |
1609 |
|
|
c middle clouds (ict) |
1610 |
|
|
c level indiex separating middle and n/d m*n |
1611 |
|
|
c low clouds (icb) |
1612 |
|
|
c cloud amount (fcld) fraction m*ndim*np |
1613 |
|
|
c cloud amounts of high, middle, and n/d m*n*3 |
1614 |
|
|
c low clouds (cc) |
1615 |
|
|
c aerosol optical thickness (taual) n/d m*ndim*np |
1616 |
|
|
c cosecant of the solar zenith angle (csm) n/d m*n |
1617 |
|
|
c uv+par surface albedo for beam fraction m*ndim |
1618 |
|
|
c radiation (rsuvbm) |
1619 |
|
|
c uv+par surface albedo for diffuse fraction m*ndim |
1620 |
|
|
c radiation (rsuvdf) |
1621 |
|
|
c |
1622 |
|
|
c----- output (updated) parameters: |
1623 |
|
|
c |
1624 |
|
|
c all-sky net downward flux (flx) fraction m*ndim*(np+1) |
1625 |
|
|
c clear-sky net downward flux (flc) fraction m*ndim*(np+1) |
1626 |
|
|
c all-sky direct downward par flux at |
1627 |
|
|
c the surface (fdirpar) fraction m*ndim |
1628 |
|
|
c all-sky diffuse downward par flux at |
1629 |
|
|
c the surface (fdifpar) fraction m*ndim |
1630 |
|
|
c |
1631 |
|
|
c----- note: the following parameters must be specified by users: |
1632 |
|
|
c |
1633 |
|
|
c aerosol single scattering albedo (ssaal) n/d 1 |
1634 |
|
|
c aerosol asymmetry factor (asyal) n/d 1 |
1635 |
|
|
c |
1636 |
|
|
* |
1637 |
|
|
c*********************************************************************** |
1638 |
|
|
|
1639 |
|
|
implicit none |
1640 |
|
|
|
1641 |
|
|
c-----Explicit Inline Directives |
1642 |
|
|
|
1643 |
molod |
1.11 |
#ifdef CRAY |
1644 |
|
|
#ifdef f77 |
1645 |
molod |
1.1 |
cfpp$ expand (deledd) |
1646 |
|
|
cfpp$ expand (sagpol) |
1647 |
|
|
#endif |
1648 |
|
|
#endif |
1649 |
|
|
|
1650 |
|
|
c-----input parameters |
1651 |
|
|
|
1652 |
|
|
integer m,n,ndim,np,ict,icb |
1653 |
|
|
real taucld(m,ndim,np,2),reff(m,ndim,np,2),fcld(m,ndim,np) |
1654 |
|
|
real tauclb(m,n,np),tauclf(m,n,np),cc(m,n,3) |
1655 |
|
|
real oh(m,n,np),dp(m,n,np),taual(m,ndim,np) |
1656 |
|
|
real rsuvbm(m,ndim),rsuvdf(m,ndim),csm(m,n) |
1657 |
|
|
|
1658 |
|
|
c-----output (updated) parameter |
1659 |
|
|
|
1660 |
|
|
real flx(m,ndim,np+1),flc(m,ndim,np+1) |
1661 |
|
|
real fdirpar(m,ndim),fdifpar(m,ndim) |
1662 |
|
|
real fdiruv(m,ndim),fdifuv(m,ndim) |
1663 |
|
|
|
1664 |
|
|
c-----static parameters |
1665 |
|
|
|
1666 |
|
|
integer nband |
1667 |
|
|
parameter (nband=8) |
1668 |
|
|
real hk(nband),xk(nband),ry(nband) |
1669 |
|
|
real asyal(nband),ssaal(nband),aig(3),awg(3) |
1670 |
|
|
|
1671 |
|
|
c-----temporary array |
1672 |
|
|
|
1673 |
|
|
integer i,j,k,ib |
1674 |
|
|
real taurs,tauoz,tausto,ssatau,asysto,tauto,ssato,asyto |
1675 |
|
|
real taux,reff1,reff2,g1,g2,asycl(m,n,np) |
1676 |
|
|
real td(m,n,np+1,2),rr(m,n,np+1,2),tt(m,n,np+1,2), |
1677 |
|
|
* rs(m,n,np+1,2),ts(m,n,np+1,2) |
1678 |
|
|
real fall(m,n,np+1),fclr(m,n,np+1),fsdir(m,n),fsdif(m,n) |
1679 |
|
|
real asyclt(m,n) |
1680 |
|
|
real rr1t(m,n),tt1t(m,n),td1t(m,n),rs1t(m,n),ts1t(m,n) |
1681 |
|
|
real rr2t(m,n),tt2t(m,n),td2t(m,n),rs2t(m,n),ts2t(m,n) |
1682 |
|
|
|
1683 |
|
|
c-----hk is the fractional extra-terrestrial solar flux. |
1684 |
|
|
c the sum of hk is 0.47074. |
1685 |
|
|
|
1686 |
|
|
data hk/.00057, .00367, .00083, .00417, |
1687 |
|
|
* .00600, .00556, .05913, .39081/ |
1688 |
|
|
|
1689 |
|
|
c-----xk is the ozone absorption coefficient. unit: /(cm-atm)stp |
1690 |
|
|
|
1691 |
|
|
data xk /30.47, 187.2, 301.9, 42.83, |
1692 |
|
|
* 7.09, 1.25, 0.0345, 0.0539/ |
1693 |
|
|
|
1694 |
|
|
c-----ry is the extinction coefficient for Rayleigh scattering. |
1695 |
|
|
c unit: /mb. |
1696 |
|
|
|
1697 |
|
|
data ry /.00604, .00170, .00222, .00132, |
1698 |
|
|
* .00107, .00091, .00055, .00012/ |
1699 |
|
|
|
1700 |
|
|
c-----aerosol single-scattering albedo and asymmetry factor |
1701 |
|
|
|
1702 |
|
|
data ssaal,asyal/nband*0.999,nband*0.850/ |
1703 |
|
|
|
1704 |
|
|
c-----coefficients for computing the asymmetry factor of ice clouds |
1705 |
|
|
c from asycl=aig(*,1)+aig(*,2)*reff+aig(*,3)*reff**2 |
1706 |
|
|
|
1707 |
|
|
data aig/.74625000,.00105410,-.00000264/ |
1708 |
|
|
|
1709 |
|
|
c-----coefficients for computing the asymmetry factor of liquid |
1710 |
|
|
c clouds from asycl=awg(*,1)+awg(*,2)*reff+awg(*,3)*reff**2 |
1711 |
|
|
|
1712 |
|
|
data awg/.82562000,.00529000,-.00014866/ |
1713 |
|
|
|
1714 |
|
|
c-----initialize surface reflectances and transmittances |
1715 |
|
|
|
1716 |
|
|
do j= 1, n |
1717 |
|
|
do i= 1, m |
1718 |
|
|
rr(i,j,np+1,1)=rsuvbm(i,j) |
1719 |
|
|
rr(i,j,np+1,2)=rsuvbm(i,j) |
1720 |
|
|
rs(i,j,np+1,1)=rsuvdf(i,j) |
1721 |
|
|
rs(i,j,np+1,2)=rsuvdf(i,j) |
1722 |
|
|
td(i,j,np+1,1)=0.0 |
1723 |
|
|
td(i,j,np+1,2)=0.0 |
1724 |
|
|
tt(i,j,np+1,1)=0.0 |
1725 |
|
|
tt(i,j,np+1,2)=0.0 |
1726 |
|
|
ts(i,j,np+1,1)=0.0 |
1727 |
|
|
ts(i,j,np+1,2)=0.0 |
1728 |
|
|
enddo |
1729 |
|
|
enddo |
1730 |
|
|
|
1731 |
|
|
c-----compute cloud asymmetry factor for a mixture of |
1732 |
|
|
c liquid and ice particles. unit of reff is micrometers. |
1733 |
|
|
|
1734 |
|
|
do k= 1, np |
1735 |
|
|
|
1736 |
|
|
do j= 1, n |
1737 |
|
|
do i= 1, m |
1738 |
|
|
|
1739 |
|
|
asyclt(i,j)=1.0 |
1740 |
|
|
|
1741 |
|
|
taux=taucld(i,j,k,1)+taucld(i,j,k,2) |
1742 |
|
|
if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then |
1743 |
|
|
|
1744 |
|
|
reff1=min(reff(i,j,k,1),130.) |
1745 |
|
|
reff2=min(reff(i,j,k,2),20.0) |
1746 |
|
|
|
1747 |
|
|
g1=(aig(1)+(aig(2)+aig(3)*reff1)*reff1)*taucld(i,j,k,1) |
1748 |
|
|
g2=(awg(1)+(awg(2)+awg(3)*reff2)*reff2)*taucld(i,j,k,2) |
1749 |
|
|
asyclt(i,j)=(g1+g2)/taux |
1750 |
|
|
|
1751 |
|
|
endif |
1752 |
|
|
|
1753 |
|
|
enddo |
1754 |
|
|
enddo |
1755 |
|
|
|
1756 |
|
|
do j=1,n |
1757 |
|
|
do i=1,m |
1758 |
|
|
asycl(i,j,k)=asyclt(i,j) |
1759 |
|
|
enddo |
1760 |
|
|
enddo |
1761 |
|
|
|
1762 |
|
|
enddo |
1763 |
|
|
|
1764 |
|
|
do j=1,n |
1765 |
|
|
do i=1,m |
1766 |
|
|
fdiruv(i,j)=0.0 |
1767 |
|
|
fdifuv(i,j)=0.0 |
1768 |
|
|
enddo |
1769 |
|
|
enddo |
1770 |
|
|
|
1771 |
|
|
c-----integration over spectral bands |
1772 |
|
|
|
1773 |
|
|
do 100 ib=1,nband |
1774 |
|
|
|
1775 |
|
|
do 300 k= 1, np |
1776 |
|
|
|
1777 |
|
|
do j= 1, n |
1778 |
|
|
do i= 1, m |
1779 |
|
|
|
1780 |
|
|
c-----compute ozone and rayleigh optical thicknesses |
1781 |
|
|
|
1782 |
|
|
taurs=ry(ib)*dp(i,j,k) |
1783 |
|
|
tauoz=xk(ib)*oh(i,j,k) |
1784 |
|
|
|
1785 |
|
|
c-----compute total optical thickness, single scattering albedo, |
1786 |
|
|
c and asymmetry factor |
1787 |
|
|
|
1788 |
|
|
tausto=taurs+tauoz+taual(i,j,k)+1.0e-8 |
1789 |
|
|
ssatau=ssaal(ib)*taual(i,j,k)+taurs |
1790 |
|
|
asysto=asyal(ib)*ssaal(ib)*taual(i,j,k) |
1791 |
|
|
|
1792 |
|
|
c-----compute reflectance and transmittance for cloudless layers |
1793 |
|
|
|
1794 |
|
|
tauto=tausto |
1795 |
|
|
ssato=ssatau/tauto+1.0e-8 |
1796 |
|
|
ssato=min(ssato,0.999999) |
1797 |
|
|
asyto=asysto/(ssato*tauto) |
1798 |
|
|
|
1799 |
|
|
call deledd(tauto,ssato,asyto,csm(i,j), |
1800 |
|
|
* rr1t(i,j),tt1t(i,j),td1t(i,j)) |
1801 |
|
|
|
1802 |
|
|
call sagpol (tauto,ssato,asyto,rs1t(i,j),ts1t(i,j)) |
1803 |
|
|
|
1804 |
|
|
c-----compute reflectance and transmittance for cloud layers |
1805 |
|
|
|
1806 |
|
|
if (tauclb(i,j,k) .lt. 0.01) then |
1807 |
|
|
|
1808 |
|
|
rr2t(i,j)=rr1t(i,j) |
1809 |
|
|
tt2t(i,j)=tt1t(i,j) |
1810 |
|
|
td2t(i,j)=td1t(i,j) |
1811 |
|
|
rs2t(i,j)=rs1t(i,j) |
1812 |
|
|
ts2t(i,j)=ts1t(i,j) |
1813 |
|
|
|
1814 |
|
|
else |
1815 |
|
|
|
1816 |
|
|
tauto=tausto+tauclb(i,j,k) |
1817 |
|
|
ssato=(ssatau+tauclb(i,j,k))/tauto+1.0e-8 |
1818 |
|
|
ssato=min(ssato,0.999999) |
1819 |
|
|
asyto=(asysto+asycl(i,j,k)*tauclb(i,j,k))/(ssato*tauto) |
1820 |
|
|
|
1821 |
|
|
call deledd(tauto,ssato,asyto,csm(i,j), |
1822 |
|
|
* rr2t(i,j),tt2t(i,j),td2t(i,j)) |
1823 |
|
|
|
1824 |
|
|
tauto=tausto+tauclf(i,j,k) |
1825 |
|
|
ssato=(ssatau+tauclf(i,j,k))/tauto+1.0e-8 |
1826 |
|
|
ssato=min(ssato,0.999999) |
1827 |
|
|
asyto=(asysto+asycl(i,j,k)*tauclf(i,j,k))/(ssato*tauto) |
1828 |
|
|
|
1829 |
|
|
call sagpol (tauto,ssato,asyto,rs2t(i,j),ts2t(i,j)) |
1830 |
|
|
|
1831 |
|
|
endif |
1832 |
|
|
|
1833 |
|
|
enddo |
1834 |
|
|
enddo |
1835 |
|
|
|
1836 |
|
|
do j=1,n |
1837 |
|
|
do i=1,m |
1838 |
|
|
rr(i,j,k,1)=rr1t(i,j) |
1839 |
|
|
enddo |
1840 |
|
|
enddo |
1841 |
|
|
do j=1,n |
1842 |
|
|
do i=1,m |
1843 |
|
|
tt(i,j,k,1)=tt1t(i,j) |
1844 |
|
|
enddo |
1845 |
|
|
enddo |
1846 |
|
|
do j=1,n |
1847 |
|
|
do i=1,m |
1848 |
|
|
td(i,j,k,1)=td1t(i,j) |
1849 |
|
|
enddo |
1850 |
|
|
enddo |
1851 |
|
|
do j=1,n |
1852 |
|
|
do i=1,m |
1853 |
|
|
rs(i,j,k,1)=rs1t(i,j) |
1854 |
|
|
enddo |
1855 |
|
|
enddo |
1856 |
|
|
do j=1,n |
1857 |
|
|
do i=1,m |
1858 |
|
|
ts(i,j,k,1)=ts1t(i,j) |
1859 |
|
|
enddo |
1860 |
|
|
enddo |
1861 |
|
|
|
1862 |
|
|
do j=1,n |
1863 |
|
|
do i=1,m |
1864 |
|
|
rr(i,j,k,2)=rr2t(i,j) |
1865 |
|
|
enddo |
1866 |
|
|
enddo |
1867 |
|
|
do j=1,n |
1868 |
|
|
do i=1,m |
1869 |
|
|
tt(i,j,k,2)=tt2t(i,j) |
1870 |
|
|
enddo |
1871 |
|
|
enddo |
1872 |
|
|
do j=1,n |
1873 |
|
|
do i=1,m |
1874 |
|
|
td(i,j,k,2)=td2t(i,j) |
1875 |
|
|
enddo |
1876 |
|
|
enddo |
1877 |
|
|
do j=1,n |
1878 |
|
|
do i=1,m |
1879 |
|
|
rs(i,j,k,2)=rs2t(i,j) |
1880 |
|
|
enddo |
1881 |
|
|
enddo |
1882 |
|
|
do j=1,n |
1883 |
|
|
do i=1,m |
1884 |
|
|
ts(i,j,k,2)=ts2t(i,j) |
1885 |
|
|
enddo |
1886 |
|
|
enddo |
1887 |
|
|
|
1888 |
|
|
300 continue |
1889 |
|
|
|
1890 |
|
|
c-----flux calculations |
1891 |
|
|
|
1892 |
|
|
call cldflx (m,n,np,ict,icb,cc,rr,tt,td,rs,ts, |
1893 |
|
|
* fclr,fall,fsdir,fsdif) |
1894 |
|
|
|
1895 |
|
|
do k= 1, np+1 |
1896 |
|
|
do j= 1, n |
1897 |
|
|
do i= 1, m |
1898 |
|
|
flx(i,j,k)=flx(i,j,k)+fall(i,j,k)*hk(ib) |
1899 |
|
|
enddo |
1900 |
|
|
enddo |
1901 |
|
|
do j= 1, n |
1902 |
|
|
do i= 1, m |
1903 |
|
|
flc(i,j,k)=flc(i,j,k)+fclr(i,j,k)*hk(ib) |
1904 |
|
|
enddo |
1905 |
|
|
enddo |
1906 |
|
|
enddo |
1907 |
|
|
|
1908 |
|
|
if(ib.eq.nband) then |
1909 |
|
|
do j=1,n |
1910 |
|
|
do i=1,m |
1911 |
|
|
fdirpar(i,j)=fsdir(i,j)*hk(ib) |
1912 |
|
|
fdifpar(i,j)=fsdif(i,j)*hk(ib) |
1913 |
|
|
enddo |
1914 |
|
|
enddo |
1915 |
|
|
else |
1916 |
|
|
do j=1,n |
1917 |
|
|
do i=1,m |
1918 |
|
|
fdiruv(i,j)=fdiruv(i,j)+fsdir(i,j)*hk(ib) |
1919 |
|
|
fdifuv(i,j)=fdifuv(i,j)+fsdif(i,j)*hk(ib) |
1920 |
|
|
enddo |
1921 |
|
|
enddo |
1922 |
|
|
endif |
1923 |
|
|
|
1924 |
|
|
100 continue |
1925 |
|
|
|
1926 |
|
|
return |
1927 |
|
|
end |
1928 |
|
|
|
1929 |
|
|
c********************************************************************* |
1930 |
|
|
|
1931 |
|
|
subroutine deledd(tau,ssc,g0,csm,rr,tt,td) |
1932 |
|
|
|
1933 |
|
|
c********************************************************************* |
1934 |
|
|
c |
1935 |
|
|
c-----uses the delta-eddington approximation to compute the |
1936 |
|
|
c bulk scattering properties of a single layer |
1937 |
|
|
c coded following King and Harshvardhan (JAS, 1986) |
1938 |
|
|
c |
1939 |
|
|
c inputs: |
1940 |
|
|
c |
1941 |
|
|
c tau: the effective optical thickness |
1942 |
|
|
c ssc: the effective single scattering albedo |
1943 |
|
|
c g0: the effective asymmetry factor |
1944 |
|
|
c csm: the effective secant of the zenith angle |
1945 |
|
|
c |
1946 |
|
|
c outputs: |
1947 |
|
|
c |
1948 |
|
|
c rr: the layer reflection of the direct beam |
1949 |
|
|
c tt: the layer diffuse transmission of the direct beam |
1950 |
|
|
c td: the layer direct transmission of the direct beam |
1951 |
|
|
|
1952 |
|
|
c********************************************************************* |
1953 |
|
|
|
1954 |
|
|
implicit none |
1955 |
|
|
|
1956 |
|
|
c-----Explicit Inline Directives |
1957 |
|
|
|
1958 |
molod |
1.11 |
#ifdef CRAY |
1959 |
|
|
#ifdef f77 |
1960 |
molod |
1.1 |
cfpp$ expand (expmn) |
1961 |
|
|
#endif |
1962 |
|
|
#endif |
1963 |
|
|
real expmn |
1964 |
|
|
|
1965 |
|
|
real zero,one,two,three,four,fourth,seven,tumin |
1966 |
|
|
parameter (one=1., three=3.) |
1967 |
|
|
parameter (seven=7., two=2.) |
1968 |
|
|
parameter (four=4., fourth=.25) |
1969 |
|
|
parameter (zero=0., tumin=1.e-20) |
1970 |
|
|
|
1971 |
|
|
c-----input parameters |
1972 |
|
|
real tau,ssc,g0,csm |
1973 |
|
|
|
1974 |
|
|
c-----output parameters |
1975 |
|
|
real rr,tt,td |
1976 |
|
|
|
1977 |
|
|
c-----temporary parameters |
1978 |
|
|
|
1979 |
|
|
real zth,ff,xx,taup,sscp,gp,gm1,gm2,gm3,akk,alf1,alf2, |
1980 |
|
|
* all,bll,st7,st8,cll,dll,fll,ell,st1,st2,st3,st4 |
1981 |
|
|
c |
1982 |
|
|
zth = one / csm |
1983 |
|
|
|
1984 |
|
|
c delta-eddington scaling of single scattering albedo, |
1985 |
|
|
c optical thickness, and asymmetry factor, |
1986 |
|
|
c K & H eqs(27-29) |
1987 |
|
|
|
1988 |
|
|
ff = g0*g0 |
1989 |
|
|
xx = one-ff*ssc |
1990 |
|
|
taup= tau*xx |
1991 |
|
|
sscp= ssc*(one-ff)/xx |
1992 |
|
|
gp = g0/(one+g0) |
1993 |
|
|
|
1994 |
|
|
c extinction of the direct beam transmission |
1995 |
|
|
|
1996 |
|
|
td = expmn(-taup*csm) |
1997 |
|
|
|
1998 |
|
|
c gamma1, gamma2, gamma3 and gamma4. see table 2 and eq(26) K & H |
1999 |
|
|
c ssc and gp are the d-s single scattering |
2000 |
|
|
c albedo and asymmetry factor. |
2001 |
|
|
|
2002 |
|
|
xx = three*gp |
2003 |
|
|
gm1 = (seven - sscp*(four+xx))*fourth |
2004 |
|
|
gm2 = -(one - sscp*(four-xx))*fourth |
2005 |
|
|
gm3 = (two - zth*xx )*fourth |
2006 |
|
|
|
2007 |
|
|
c akk is k as defined in eq(25) of K & H |
2008 |
|
|
|
2009 |
|
|
xx = gm1-gm2 |
2010 |
|
|
akk = sqrt((gm1+gm2)*xx) |
2011 |
|
|
|
2012 |
|
|
c alf1 and alf2 are alpha1 and alpha2 from eqs (23) & (24) of K & H |
2013 |
|
|
|
2014 |
|
|
alf1 = gm1 - gm3 * xx |
2015 |
|
|
alf2 = gm2 + gm3 * xx |
2016 |
|
|
|
2017 |
|
|
c all is last term in eq(21) of K & H |
2018 |
|
|
c bll is last term in eq(22) of K & H |
2019 |
|
|
|
2020 |
|
|
xx = akk * two |
2021 |
|
|
all = (gm3 - alf2 * zth )*xx*td |
2022 |
|
|
bll = (one - gm3 + alf1*zth)*xx |
2023 |
|
|
|
2024 |
|
|
xx = akk * zth |
2025 |
|
|
st7 = one - xx |
2026 |
|
|
st8 = one + xx |
2027 |
|
|
|
2028 |
|
|
xx = akk * gm3 |
2029 |
|
|
cll = (alf2 + xx) * st7 |
2030 |
|
|
dll = (alf2 - xx) * st8 |
2031 |
|
|
|
2032 |
|
|
xx = akk * (one-gm3) |
2033 |
|
|
fll = (alf1 + xx) * st8 |
2034 |
|
|
ell = (alf1 - xx) * st7 |
2035 |
|
|
|
2036 |
|
|
st3 = max(abs(st7*st8),tumin) |
2037 |
|
|
st3 = sign (st3,st7) |
2038 |
|
|
|
2039 |
|
|
st2 = expmn(-akk*taup) |
2040 |
|
|
st4 = st2 * st2 |
2041 |
|
|
|
2042 |
|
|
st1 = sscp / ((akk+gm1 + (akk-gm1)*st4) * st3) |
2043 |
|
|
|
2044 |
|
|
c rr is r-hat of eq(21) of K & H |
2045 |
|
|
c tt is diffuse part of t-hat of eq(22) of K & H |
2046 |
|
|
|
2047 |
|
|
rr = ( cll-dll*st4 -all*st2)*st1 |
2048 |
|
|
tt = - ((fll-ell*st4)*td-bll*st2)*st1 |
2049 |
|
|
|
2050 |
|
|
rr = max(rr,zero) |
2051 |
|
|
tt = max(tt,zero) |
2052 |
|
|
|
2053 |
|
|
return |
2054 |
|
|
end |
2055 |
|
|
|
2056 |
|
|
c********************************************************************* |
2057 |
|
|
|
2058 |
|
|
subroutine sagpol(tau,ssc,g0,rll,tll) |
2059 |
|
|
|
2060 |
|
|
c********************************************************************* |
2061 |
|
|
c-----transmittance (tll) and reflectance (rll) of diffuse radiation |
2062 |
|
|
c follows Sagan and Pollock (JGR, 1967). |
2063 |
|
|
c also, eq.(31) of Lacis and Hansen (JAS, 1974). |
2064 |
|
|
c |
2065 |
|
|
c-----input parameters: |
2066 |
|
|
c |
2067 |
|
|
c tau: the effective optical thickness |
2068 |
|
|
c ssc: the effective single scattering albedo |
2069 |
|
|
c g0: the effective asymmetry factor |
2070 |
|
|
c |
2071 |
|
|
c-----output parameters: |
2072 |
|
|
c |
2073 |
|
|
c rll: the layer reflection of diffuse radiation |
2074 |
|
|
c tll: the layer transmission of diffuse radiation |
2075 |
|
|
c |
2076 |
|
|
c********************************************************************* |
2077 |
|
|
|
2078 |
|
|
implicit none |
2079 |
|
|
|
2080 |
|
|
c-----Explicit Inline Directives |
2081 |
|
|
|
2082 |
molod |
1.11 |
#ifdef CRAY |
2083 |
|
|
#ifdef f77 |
2084 |
molod |
1.1 |
cfpp$ expand (expmn) |
2085 |
|
|
#endif |
2086 |
|
|
#endif |
2087 |
|
|
real expmn |
2088 |
|
|
|
2089 |
|
|
real one,three,four |
2090 |
|
|
parameter (one=1., three=3., four=4.) |
2091 |
|
|
|
2092 |
|
|
c-----output parameters: |
2093 |
|
|
|
2094 |
|
|
real tau,ssc,g0 |
2095 |
|
|
|
2096 |
|
|
c-----output parameters: |
2097 |
|
|
|
2098 |
|
|
real rll,tll |
2099 |
|
|
|
2100 |
|
|
c-----temporary arrays |
2101 |
|
|
|
2102 |
|
|
real xx,uuu,ttt,emt,up1,um1,st1 |
2103 |
|
|
|
2104 |
|
|
xx = one-ssc*g0 |
2105 |
|
|
uuu = sqrt( xx/(one-ssc)) |
2106 |
|
|
ttt = sqrt( xx*(one-ssc)*three )*tau |
2107 |
|
|
emt = expmn(-ttt) |
2108 |
|
|
up1 = uuu + one |
2109 |
|
|
um1 = uuu - one |
2110 |
|
|
xx = um1*emt |
2111 |
|
|
st1 = one / ((up1+xx) * (up1-xx)) |
2112 |
|
|
rll = up1*um1*(one-emt*emt)*st1 |
2113 |
|
|
tll = uuu*four*emt *st1 |
2114 |
|
|
|
2115 |
|
|
return |
2116 |
|
|
end |
2117 |
|
|
|
2118 |
|
|
c******************************************************************* |
2119 |
|
|
|
2120 |
|
|
function expmn(fin) |
2121 |
|
|
|
2122 |
|
|
c******************************************************************* |
2123 |
|
|
c compute exponential for arguments in the range 0> fin > -10. |
2124 |
molod |
1.10 |
c******************************************************************* |
2125 |
|
|
implicit none |
2126 |
|
|
real fin,expmn |
2127 |
molod |
1.1 |
|
2128 |
molod |
1.10 |
real one,expmin,e1,e2,e3,e4 |
2129 |
molod |
1.1 |
parameter (one=1.0, expmin=-10.0) |
2130 |
|
|
parameter (e1=1.0, e2=-2.507213e-1) |
2131 |
|
|
parameter (e3=2.92732e-2, e4=-3.827800e-3) |
2132 |
|
|
|
2133 |
|
|
if (fin .lt. expmin) fin = expmin |
2134 |
|
|
expmn = ((e4*fin + e3)*fin+e2)*fin+e1 |
2135 |
|
|
expmn = expmn * expmn |
2136 |
|
|
expmn = one / (expmn * expmn) |
2137 |
|
|
|
2138 |
|
|
return |
2139 |
|
|
end |
2140 |
|
|
|
2141 |
|
|
c******************************************************************* |
2142 |
|
|
|
2143 |
|
|
subroutine cldflx (m,n,np,ict,icb,cc,rr,tt,td,rs,ts, |
2144 |
|
|
* fclr,fall,fsdir,fsdif) |
2145 |
|
|
|
2146 |
|
|
c******************************************************************* |
2147 |
|
|
c compute upward and downward fluxes using a two-stream adding method |
2148 |
|
|
c following equations (3)-(5) of Chou (1992, JAS). |
2149 |
|
|
c |
2150 |
|
|
c clouds are grouped into high, middle, and low clouds which are |
2151 |
|
|
c assumed randomly overlapped. It involves eight sets of calculations. |
2152 |
|
|
c In each set of calculations, each atmospheric layer is homogeneous, |
2153 |
|
|
c either with clouds or without clouds. |
2154 |
|
|
|
2155 |
|
|
c input parameters: |
2156 |
|
|
c m: number of soundings in zonal direction |
2157 |
|
|
c n: number of soundings in meridional direction |
2158 |
|
|
c np: number of atmospheric layers |
2159 |
|
|
c ict: the level separating high and middle clouds |
2160 |
|
|
c icb: the level separating middle and low clouds |
2161 |
|
|
c cc: effective cloud covers for high, middle and low clouds |
2162 |
|
|
c tt: diffuse transmission of a layer illuminated by beam radiation |
2163 |
|
|
c td: direct beam tranmssion |
2164 |
|
|
c ts: transmission of a layer illuminated by diffuse radiation |
2165 |
|
|
c rr: reflection of a layer illuminated by beam radiation |
2166 |
|
|
c rs: reflection of a layer illuminated by diffuse radiation |
2167 |
|
|
c |
2168 |
|
|
c output parameters: |
2169 |
|
|
c fclr: clear-sky flux (downward minus upward) |
2170 |
|
|
c fall: all-sky flux (downward minus upward) |
2171 |
|
|
c fdndir: surface direct downward flux |
2172 |
|
|
c fdndif: surface diffuse downward flux |
2173 |
|
|
c*********************************************************************c |
2174 |
|
|
|
2175 |
|
|
implicit none |
2176 |
|
|
|
2177 |
|
|
c-----input parameters |
2178 |
|
|
|
2179 |
|
|
integer m,n,np,ict,icb |
2180 |
|
|
|
2181 |
|
|
real rr(m,n,np+1,2),tt(m,n,np+1,2),td(m,n,np+1,2) |
2182 |
|
|
real rs(m,n,np+1,2),ts(m,n,np+1,2) |
2183 |
|
|
real cc(m,n,3) |
2184 |
|
|
|
2185 |
|
|
c-----temporary array |
2186 |
|
|
|
2187 |
|
|
integer i,j,k,ih,im,is |
2188 |
|
|
real rra(m,n,np+1,2,2),tta(m,n,np+1,2,2),tda(m,n,np+1,2,2) |
2189 |
|
|
real rsa(m,n,np+1,2,2),rxa(m,n,np+1,2,2) |
2190 |
|
|
real ch(m,n),cm(m,n),ct(m,n),flxdn(m,n,np+1) |
2191 |
|
|
real fdndir(m,n),fdndif(m,n),fupdif |
2192 |
|
|
real denm,xx |
2193 |
|
|
|
2194 |
|
|
c-----output parameters |
2195 |
|
|
|
2196 |
|
|
real fclr(m,n,np+1),fall(m,n,np+1) |
2197 |
|
|
real fsdir(m,n),fsdif(m,n) |
2198 |
|
|
|
2199 |
|
|
c-----initialize all-sky flux (fall) and surface downward fluxes |
2200 |
|
|
|
2201 |
|
|
do k=1,np+1 |
2202 |
|
|
do j=1,n |
2203 |
|
|
do i=1,m |
2204 |
|
|
fall(i,j,k)=0.0 |
2205 |
|
|
enddo |
2206 |
|
|
enddo |
2207 |
|
|
enddo |
2208 |
|
|
|
2209 |
|
|
do j=1,n |
2210 |
|
|
do i=1,m |
2211 |
|
|
fsdir(i,j)=0.0 |
2212 |
|
|
fsdif(i,j)=0.0 |
2213 |
|
|
enddo |
2214 |
|
|
enddo |
2215 |
|
|
|
2216 |
|
|
c-----compute transmittances and reflectances for a composite of |
2217 |
|
|
c layers. layers are added one at a time, going down from the top. |
2218 |
|
|
c tda is the composite transmittance illuminated by beam radiation |
2219 |
|
|
c tta is the composite diffuse transmittance illuminated by |
2220 |
|
|
c beam radiation |
2221 |
|
|
c rsa is the composite reflectance illuminated from below |
2222 |
|
|
c by diffuse radiation |
2223 |
|
|
c tta and rsa are computed from eqs. (4b) and (3b) of Chou |
2224 |
|
|
|
2225 |
|
|
c-----for high clouds. indices 1 and 2 denote clear and cloudy |
2226 |
|
|
c situations, respectively. |
2227 |
|
|
|
2228 |
|
|
do 10 ih=1,2 |
2229 |
|
|
|
2230 |
|
|
do j= 1, n |
2231 |
|
|
do i= 1, m |
2232 |
|
|
tda(i,j,1,ih,1)=td(i,j,1,ih) |
2233 |
|
|
tta(i,j,1,ih,1)=tt(i,j,1,ih) |
2234 |
|
|
rsa(i,j,1,ih,1)=rs(i,j,1,ih) |
2235 |
|
|
tda(i,j,1,ih,2)=td(i,j,1,ih) |
2236 |
|
|
tta(i,j,1,ih,2)=tt(i,j,1,ih) |
2237 |
|
|
rsa(i,j,1,ih,2)=rs(i,j,1,ih) |
2238 |
|
|
enddo |
2239 |
|
|
enddo |
2240 |
|
|
|
2241 |
|
|
do k= 2, ict-1 |
2242 |
|
|
do j= 1, n |
2243 |
|
|
do i= 1, m |
2244 |
|
|
denm = ts(i,j,k,ih)/( 1.-rsa(i,j,k-1,ih,1)*rs(i,j,k,ih)) |
2245 |
|
|
tda(i,j,k,ih,1)= tda(i,j,k-1,ih,1)*td(i,j,k,ih) |
2246 |
|
|
tta(i,j,k,ih,1)= tda(i,j,k-1,ih,1)*tt(i,j,k,ih) |
2247 |
|
|
* +(tda(i,j,k-1,ih,1)*rr(i,j,k,ih) |
2248 |
|
|
* *rsa(i,j,k-1,ih,1)+tta(i,j,k-1,ih,1))*denm |
2249 |
|
|
rsa(i,j,k,ih,1)= rs(i,j,k,ih)+ts(i,j,k,ih) |
2250 |
|
|
* *rsa(i,j,k-1,ih,1)*denm |
2251 |
|
|
tda(i,j,k,ih,2)= tda(i,j,k,ih,1) |
2252 |
|
|
tta(i,j,k,ih,2)= tta(i,j,k,ih,1) |
2253 |
|
|
rsa(i,j,k,ih,2)= rsa(i,j,k,ih,1) |
2254 |
|
|
enddo |
2255 |
|
|
enddo |
2256 |
|
|
enddo |
2257 |
|
|
|
2258 |
|
|
c-----for middle clouds |
2259 |
|
|
|
2260 |
|
|
do 10 im=1,2 |
2261 |
|
|
|
2262 |
|
|
do k= ict, icb-1 |
2263 |
|
|
do j= 1, n |
2264 |
|
|
do i= 1, m |
2265 |
|
|
denm = ts(i,j,k,im)/( 1.-rsa(i,j,k-1,ih,im)*rs(i,j,k,im)) |
2266 |
|
|
tda(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*td(i,j,k,im) |
2267 |
|
|
tta(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*tt(i,j,k,im) |
2268 |
|
|
* +(tda(i,j,k-1,ih,im)*rr(i,j,k,im) |
2269 |
|
|
* *rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm |
2270 |
|
|
rsa(i,j,k,ih,im)= rs(i,j,k,im)+ts(i,j,k,im) |
2271 |
|
|
* *rsa(i,j,k-1,ih,im)*denm |
2272 |
|
|
enddo |
2273 |
|
|
enddo |
2274 |
|
|
enddo |
2275 |
|
|
|
2276 |
|
|
10 continue |
2277 |
|
|
|
2278 |
|
|
c-----layers are added one at a time, going up from the surface. |
2279 |
|
|
c rra is the composite reflectance illuminated by beam radiation |
2280 |
|
|
c rxa is the composite reflectance illuminated from above |
2281 |
|
|
c by diffuse radiation |
2282 |
|
|
c rra and rxa are computed from eqs. (4a) and (3a) of Chou |
2283 |
|
|
|
2284 |
|
|
c-----for the low clouds |
2285 |
|
|
|
2286 |
|
|
do 20 is=1,2 |
2287 |
|
|
|
2288 |
|
|
do j= 1, n |
2289 |
|
|
do i= 1, m |
2290 |
|
|
rra(i,j,np+1,1,is)=rr(i,j,np+1,is) |
2291 |
|
|
rxa(i,j,np+1,1,is)=rs(i,j,np+1,is) |
2292 |
|
|
rra(i,j,np+1,2,is)=rr(i,j,np+1,is) |
2293 |
|
|
rxa(i,j,np+1,2,is)=rs(i,j,np+1,is) |
2294 |
|
|
enddo |
2295 |
|
|
enddo |
2296 |
|
|
|
2297 |
|
|
do k=np,icb,-1 |
2298 |
|
|
do j= 1, n |
2299 |
|
|
do i= 1, m |
2300 |
|
|
denm=ts(i,j,k,is)/( 1.-rs(i,j,k,is)*rxa(i,j,k+1,1,is) ) |
2301 |
|
|
rra(i,j,k,1,is)=rr(i,j,k,is)+(td(i,j,k,is) |
2302 |
|
|
* *rra(i,j,k+1,1,is)+tt(i,j,k,is)*rxa(i,j,k+1,1,is))*denm |
2303 |
|
|
rxa(i,j,k,1,is)= rs(i,j,k,is)+ts(i,j,k,is) |
2304 |
|
|
* *rxa(i,j,k+1,1,is)*denm |
2305 |
|
|
rra(i,j,k,2,is)=rra(i,j,k,1,is) |
2306 |
|
|
rxa(i,j,k,2,is)=rxa(i,j,k,1,is) |
2307 |
|
|
enddo |
2308 |
|
|
enddo |
2309 |
|
|
enddo |
2310 |
|
|
|
2311 |
|
|
c-----for middle clouds |
2312 |
|
|
|
2313 |
|
|
do 20 im=1,2 |
2314 |
|
|
|
2315 |
|
|
do k= icb-1,ict,-1 |
2316 |
|
|
do j= 1, n |
2317 |
|
|
do i= 1, m |
2318 |
|
|
denm=ts(i,j,k,im)/( 1.-rs(i,j,k,im)*rxa(i,j,k+1,im,is) ) |
2319 |
|
|
rra(i,j,k,im,is)= rr(i,j,k,im)+(td(i,j,k,im) |
2320 |
|
|
* *rra(i,j,k+1,im,is)+tt(i,j,k,im)*rxa(i,j,k+1,im,is))*denm |
2321 |
|
|
rxa(i,j,k,im,is)= rs(i,j,k,im)+ts(i,j,k,im) |
2322 |
|
|
* *rxa(i,j,k+1,im,is)*denm |
2323 |
|
|
enddo |
2324 |
|
|
enddo |
2325 |
|
|
enddo |
2326 |
|
|
|
2327 |
|
|
20 continue |
2328 |
|
|
|
2329 |
|
|
c-----integration over eight sky situations. |
2330 |
|
|
c ih, im, is denotes high, middle and low cloud groups. |
2331 |
|
|
|
2332 |
|
|
do 100 ih=1,2 |
2333 |
|
|
|
2334 |
|
|
c-----clear portion |
2335 |
|
|
|
2336 |
|
|
if(ih.eq.1) then |
2337 |
|
|
do j=1,n |
2338 |
|
|
do i=1,m |
2339 |
|
|
ch(i,j)=1.0-cc(i,j,1) |
2340 |
|
|
enddo |
2341 |
|
|
enddo |
2342 |
|
|
|
2343 |
|
|
else |
2344 |
|
|
|
2345 |
|
|
c-----cloudy portion |
2346 |
|
|
|
2347 |
|
|
do j=1,n |
2348 |
|
|
do i=1,m |
2349 |
|
|
ch(i,j)=cc(i,j,1) |
2350 |
|
|
enddo |
2351 |
|
|
enddo |
2352 |
|
|
|
2353 |
|
|
endif |
2354 |
|
|
|
2355 |
|
|
do 100 im=1,2 |
2356 |
|
|
|
2357 |
|
|
c-----clear portion |
2358 |
|
|
|
2359 |
|
|
if(im.eq.1) then |
2360 |
|
|
|
2361 |
|
|
do j=1,n |
2362 |
|
|
do i=1,m |
2363 |
|
|
cm(i,j)=ch(i,j)*(1.0-cc(i,j,2)) |
2364 |
|
|
enddo |
2365 |
|
|
enddo |
2366 |
|
|
|
2367 |
|
|
else |
2368 |
|
|
|
2369 |
|
|
c-----cloudy portion |
2370 |
|
|
|
2371 |
|
|
do j=1,n |
2372 |
|
|
do i=1,m |
2373 |
|
|
cm(i,j)=ch(i,j)*cc(i,j,2) |
2374 |
|
|
enddo |
2375 |
|
|
enddo |
2376 |
|
|
|
2377 |
|
|
endif |
2378 |
|
|
|
2379 |
|
|
do 100 is=1,2 |
2380 |
|
|
|
2381 |
|
|
c-----clear portion |
2382 |
|
|
|
2383 |
|
|
if(is.eq.1) then |
2384 |
|
|
|
2385 |
|
|
do j=1,n |
2386 |
|
|
do i=1,m |
2387 |
|
|
ct(i,j)=cm(i,j)*(1.0-cc(i,j,3)) |
2388 |
|
|
enddo |
2389 |
|
|
enddo |
2390 |
|
|
|
2391 |
|
|
else |
2392 |
|
|
|
2393 |
|
|
c-----cloudy portion |
2394 |
|
|
|
2395 |
|
|
do j=1,n |
2396 |
|
|
do i=1,m |
2397 |
|
|
ct(i,j)=cm(i,j)*cc(i,j,3) |
2398 |
|
|
enddo |
2399 |
|
|
enddo |
2400 |
|
|
|
2401 |
|
|
endif |
2402 |
|
|
|
2403 |
|
|
c-----add one layer at a time, going down. |
2404 |
|
|
|
2405 |
|
|
do k= icb, np |
2406 |
|
|
do j= 1, n |
2407 |
|
|
do i= 1, m |
2408 |
|
|
denm = ts(i,j,k,is)/( 1.-rsa(i,j,k-1,ih,im)*rs(i,j,k,is) ) |
2409 |
|
|
tda(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*td(i,j,k,is) |
2410 |
|
|
tta(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*tt(i,j,k,is) |
2411 |
|
|
* +(tda(i,j,k-1,ih,im)*rr(i,j,k,is) |
2412 |
|
|
* *rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm |
2413 |
|
|
rsa(i,j,k,ih,im)= rs(i,j,k,is)+ts(i,j,k,is) |
2414 |
|
|
* *rsa(i,j,k-1,ih,im)*denm |
2415 |
|
|
enddo |
2416 |
|
|
enddo |
2417 |
|
|
enddo |
2418 |
|
|
|
2419 |
|
|
c-----add one layer at a time, going up. |
2420 |
|
|
|
2421 |
|
|
do k= ict-1,1,-1 |
2422 |
|
|
do j= 1, n |
2423 |
|
|
do i= 1, m |
2424 |
|
|
denm =ts(i,j,k,ih)/(1.-rs(i,j,k,ih)*rxa(i,j,k+1,im,is)) |
2425 |
|
|
rra(i,j,k,im,is)= rr(i,j,k,ih)+(td(i,j,k,ih) |
2426 |
|
|
* *rra(i,j,k+1,im,is)+tt(i,j,k,ih)*rxa(i,j,k+1,im,is))*denm |
2427 |
|
|
rxa(i,j,k,im,is)= rs(i,j,k,ih)+ts(i,j,k,ih) |
2428 |
|
|
* *rxa(i,j,k+1,im,is)*denm |
2429 |
|
|
enddo |
2430 |
|
|
enddo |
2431 |
|
|
enddo |
2432 |
|
|
|
2433 |
|
|
c-----compute fluxes following eq (5) of Chou (1992) |
2434 |
|
|
|
2435 |
|
|
c fdndir is the direct downward flux |
2436 |
|
|
c fdndif is the diffuse downward flux |
2437 |
|
|
c fupdif is the diffuse upward flux |
2438 |
|
|
|
2439 |
|
|
do k=2,np+1 |
2440 |
|
|
do j=1, n |
2441 |
|
|
do i=1, m |
2442 |
|
|
denm= 1./(1.- rxa(i,j,k,im,is)*rsa(i,j,k-1,ih,im)) |
2443 |
|
|
fdndir(i,j)= tda(i,j,k-1,ih,im) |
2444 |
|
|
xx = tda(i,j,k-1,ih,im)*rra(i,j,k,im,is) |
2445 |
|
|
fdndif(i,j)= (xx*rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm |
2446 |
|
|
fupdif= (xx+tta(i,j,k-1,ih,im)*rxa(i,j,k,im,is))*denm |
2447 |
|
|
flxdn(i,j,k)=fdndir(i,j)+fdndif(i,j)-fupdif |
2448 |
|
|
enddo |
2449 |
|
|
enddo |
2450 |
|
|
enddo |
2451 |
|
|
|
2452 |
|
|
do j=1, n |
2453 |
|
|
do i=1, m |
2454 |
|
|
flxdn(i,j,1)=1.0-rra(i,j,1,im,is) |
2455 |
|
|
enddo |
2456 |
|
|
enddo |
2457 |
|
|
|
2458 |
|
|
c-----summation of fluxes over all (eight) sky situations. |
2459 |
|
|
|
2460 |
|
|
do k=1,np+1 |
2461 |
|
|
do j=1,n |
2462 |
|
|
do i=1,m |
2463 |
|
|
if(ih.eq.1 .and. im.eq.1 .and. is.eq.1) then |
2464 |
|
|
fclr(i,j,k)=flxdn(i,j,k) |
2465 |
|
|
endif |
2466 |
|
|
fall(i,j,k)=fall(i,j,k)+flxdn(i,j,k)*ct(i,j) |
2467 |
|
|
enddo |
2468 |
|
|
enddo |
2469 |
|
|
enddo |
2470 |
|
|
|
2471 |
|
|
do j=1,n |
2472 |
|
|
do i=1,m |
2473 |
|
|
fsdir(i,j)=fsdir(i,j)+fdndir(i,j)*ct(i,j) |
2474 |
|
|
fsdif(i,j)=fsdif(i,j)+fdndif(i,j)*ct(i,j) |
2475 |
|
|
enddo |
2476 |
|
|
enddo |
2477 |
|
|
|
2478 |
|
|
100 continue |
2479 |
|
|
|
2480 |
|
|
return |
2481 |
|
|
end |
2482 |
|
|
|
2483 |
|
|
c***************************************************************** |
2484 |
|
|
|
2485 |
|
|
subroutine flxco2(m,n,np,swc,swh,csm,df) |
2486 |
|
|
|
2487 |
|
|
c***************************************************************** |
2488 |
|
|
|
2489 |
|
|
c-----compute the reduction of clear-sky downward solar flux |
2490 |
|
|
c due to co2 absorption. |
2491 |
|
|
|
2492 |
|
|
implicit none |
2493 |
|
|
|
2494 |
|
|
c-----input parameters |
2495 |
|
|
|
2496 |
|
|
integer m,n,np |
2497 |
|
|
real csm(m,n),swc(m,n,np+1),swh(m,n,np+1),cah(22,19) |
2498 |
|
|
|
2499 |
|
|
c-----output (undated) parameter |
2500 |
|
|
|
2501 |
|
|
real df(m,n,np+1) |
2502 |
|
|
|
2503 |
|
|
c-----temporary array |
2504 |
|
|
|
2505 |
|
|
integer i,j,k,ic,iw |
2506 |
|
|
real xx,clog,wlog,dc,dw,x1,x2,y2 |
2507 |
|
|
|
2508 |
|
|
c******************************************************************** |
2509 |
|
|
c-----include co2 look-up table |
2510 |
|
|
|
2511 |
molod |
1.8 |
#include "cah-dat.h" |
2512 |
molod |
1.9 |
c save cah |
2513 |
molod |
1.1 |
|
2514 |
|
|
c******************************************************************** |
2515 |
|
|
c-----table look-up for the reduction of clear-sky solar |
2516 |
|
|
c radiation due to co2. The fraction 0.0343 is the |
2517 |
|
|
c extraterrestrial solar flux in the co2 bands. |
2518 |
|
|
|
2519 |
|
|
do k= 2, np+1 |
2520 |
|
|
do j= 1, n |
2521 |
|
|
do i= 1, m |
2522 |
|
|
xx=1./.3 |
2523 |
|
|
clog=log10(swc(i,j,k)*csm(i,j)) |
2524 |
|
|
wlog=log10(swh(i,j,k)*csm(i,j)) |
2525 |
|
|
ic=int( (clog+3.15)*xx+1.) |
2526 |
|
|
iw=int( (wlog+4.15)*xx+1.) |
2527 |
|
|
if(ic.lt.2)ic=2 |
2528 |
|
|
if(iw.lt.2)iw=2 |
2529 |
|
|
if(ic.gt.22)ic=22 |
2530 |
|
|
if(iw.gt.19)iw=19 |
2531 |
|
|
dc=clog-float(ic-2)*.3+3. |
2532 |
|
|
dw=wlog-float(iw-2)*.3+4. |
2533 |
|
|
x1=cah(1,iw-1)+(cah(1,iw)-cah(1,iw-1))*xx*dw |
2534 |
|
|
x2=cah(ic-1,iw-1)+(cah(ic-1,iw)-cah(ic-1,iw-1))*xx*dw |
2535 |
|
|
y2=x2+(cah(ic,iw-1)-cah(ic-1,iw-1))*xx*dc |
2536 |
|
|
if (x1.lt.y2) x1=y2 |
2537 |
|
|
df(i,j,k)=df(i,j,k)+0.0343*(x1-y2) |
2538 |
|
|
enddo |
2539 |
|
|
enddo |
2540 |
|
|
enddo |
2541 |
|
|
|
2542 |
|
|
return |
2543 |
|
|
end |