1 | ! |
---|
2 | ! $Header$ |
---|
3 | ! |
---|
4 | SUBROUTINE ISCCP_CLOUD_TYPES( |
---|
5 | & debug, |
---|
6 | & debugcol, |
---|
7 | & npoints, |
---|
8 | & sunlit, |
---|
9 | & nlev, |
---|
10 | & ncol, |
---|
11 | & seed, |
---|
12 | & pfull, |
---|
13 | & phalf, |
---|
14 | & qv, |
---|
15 | & cc, |
---|
16 | & conv, |
---|
17 | & dtau_s, |
---|
18 | & dtau_c, |
---|
19 | & top_height, |
---|
20 | & overlap, |
---|
21 | & tautab, |
---|
22 | & invtau, |
---|
23 | & skt, |
---|
24 | & emsfc_lw, |
---|
25 | & at,dem_s,dem_c, |
---|
26 | & fq_isccp, |
---|
27 | & totalcldarea, |
---|
28 | & meanptop, |
---|
29 | & meantaucld, |
---|
30 | & boxtau, |
---|
31 | & boxptop |
---|
32 | &) |
---|
33 | |
---|
34 | !$Id: isccp_cloud_types.F 730 2006-11-08 16:05:38Z aborella $ |
---|
35 | |
---|
36 | ! Copyright Steve Klein and Mark Webb 2002 - all rights reserved. |
---|
37 | ! |
---|
38 | ! This code is available without charge with the following conditions: |
---|
39 | ! |
---|
40 | ! 1. The code is available for scientific purposes and is not for |
---|
41 | ! commercial use. |
---|
42 | ! 2. Any improvements you make to the code should be made available |
---|
43 | ! to the to the authors for incorporation into a future release. |
---|
44 | ! 3. The code should not be used in any way that brings the authors |
---|
45 | ! or their employers into disrepute. |
---|
46 | |
---|
47 | implicit none |
---|
48 | |
---|
49 | ! NOTE: the maximum number of levels and columns is set by |
---|
50 | ! the following parameter statement |
---|
51 | |
---|
52 | INTEGER ncolprint |
---|
53 | |
---|
54 | ! ----- |
---|
55 | ! Input |
---|
56 | ! ----- |
---|
57 | |
---|
58 | INTEGER npoints ! number of model points in the horizontal |
---|
59 | c PARAMETER(npoints=6722) |
---|
60 | INTEGER nlev ! number of model levels in column |
---|
61 | INTEGER ncol ! number of subcolumns |
---|
62 | |
---|
63 | INTEGER sunlit(npoints) ! 1 for day points, 0 for night time |
---|
64 | |
---|
65 | INTEGER seed(npoints) ! seed value for random number generator |
---|
66 | c ! ( see Numerical Recipes Chapter 7) |
---|
67 | c ! It is recommended that the seed is set |
---|
68 | c ! to a different value for each model |
---|
69 | c ! gridbox it is called on, as it is |
---|
70 | c ! possible that the choice of the samec |
---|
71 | c ! seed value every time may introduce some |
---|
72 | c ! statistical bias in the results, particularly |
---|
73 | c ! for low values of NCOL. |
---|
74 | c |
---|
75 | REAL pfull(npoints,nlev) ! pressure of full model levels (Pascals) |
---|
76 | c ! pfull(npoints,1) is top level of model |
---|
77 | c ! pfull(npoints,nlev) is bottom level of model |
---|
78 | |
---|
79 | REAL phalf(npoints,nlev+1) ! pressure of half model levels (Pascals) |
---|
80 | c ! phalf(npoints,1) is top of model |
---|
81 | c ! phalf(npoints,nlev+1) is the surface pressure |
---|
82 | |
---|
83 | REAL qv(npoints,nlev) ! water vapor specific humidity (kg vapor/ kg air) |
---|
84 | c ! on full model levels |
---|
85 | |
---|
86 | REAL cc(npoints,nlev) ! input cloud cover in each model level (fraction) |
---|
87 | c ! NOTE: This is the HORIZONTAL area of each |
---|
88 | c ! grid box covered by clouds |
---|
89 | |
---|
90 | REAL conv(npoints,nlev) ! input convective cloud cover in each model level (fraction) |
---|
91 | c ! NOTE: This is the HORIZONTAL area of each |
---|
92 | c ! grid box covered by convective clouds |
---|
93 | |
---|
94 | REAL dtau_s(npoints,nlev) ! mean 0.67 micron optical depth of stratiform |
---|
95 | c ! clouds in each model level |
---|
96 | c ! NOTE: this the cloud optical depth of only the |
---|
97 | c ! cloudy part of the grid box, it is not weighted |
---|
98 | c ! with the 0 cloud optical depth of the clear |
---|
99 | c ! part of the grid box |
---|
100 | |
---|
101 | REAL dtau_c(npoints,nlev) ! mean 0.67 micron optical depth of convective |
---|
102 | c ! clouds in each |
---|
103 | c ! model level. Same note applies as in dtau_s. |
---|
104 | |
---|
105 | INTEGER overlap ! overlap type |
---|
106 | |
---|
107 | ! 1=max |
---|
108 | |
---|
109 | ! 2=rand |
---|
110 | ! 3=max/rand |
---|
111 | |
---|
112 | INTEGER top_height ! 1 = adjust top height using both a computed |
---|
113 | c ! infrared brightness temperature and the visible |
---|
114 | c ! optical depth to adjust cloud top pressure. Note |
---|
115 | c ! that this calculation is most appropriate to compare |
---|
116 | c ! to ISCCP data during sunlit hours. |
---|
117 | c ! 2 = do not adjust top height, that is cloud top |
---|
118 | c ! pressure is the actual cloud top pressure |
---|
119 | c ! in the model |
---|
120 | c ! 3 = adjust top height using only the computed |
---|
121 | c ! infrared brightness temperature. Note that this |
---|
122 | c ! calculation is most appropriate to compare to ISCCP |
---|
123 | c ! IR only algortihm (i.e. you can compare to nighttime |
---|
124 | c ! ISCCP data with this option) |
---|
125 | |
---|
126 | REAL tautab(0:255) ! ISCCP table for converting count value to |
---|
127 | c ! optical thickness |
---|
128 | |
---|
129 | INTEGER invtau(-20:45000) ! ISCCP table for converting optical thickness |
---|
130 | c ! to count value |
---|
131 | ! |
---|
132 | ! The following input variables are used only if top_height = 1 or top_height = 3 |
---|
133 | ! |
---|
134 | REAL skt(npoints) ! skin Temperature (K) |
---|
135 | REAL emsfc_lw ! 10.5 micron emissivity of surface (fraction) |
---|
136 | REAL at(npoints,nlev) ! temperature in each model level (K) |
---|
137 | REAL dem_s(npoints,nlev) ! 10.5 micron longwave emissivity of stratiform |
---|
138 | c ! clouds in each |
---|
139 | c ! model level. Same note applies as in dtau_s. |
---|
140 | REAL dem_c(npoints,nlev) ! 10.5 micron longwave emissivity of convective |
---|
141 | c ! clouds in each |
---|
142 | c ! model level. Same note applies as in dtau_s. |
---|
143 | cIM reg.dyn BEG |
---|
144 | REAL t1, t2 |
---|
145 | ! REAL w(npoints) !vertical wind at 500 hPa |
---|
146 | ! LOGICAL pct_ocean(npoints) !TRUE if oceanic point, FALSE otherway |
---|
147 | ! INTEGER iw(npoints) , nw |
---|
148 | ! REAL wmin, pas_w |
---|
149 | ! INTEGER k, l, iwmx |
---|
150 | ! PARAMETER(wmin=-100.,pas_w=10.,iwmx=30) |
---|
151 | ! REAL fq_dynreg(7,7,iwmx) |
---|
152 | ! REAL nfq_dynreg(7,7,iwmx) |
---|
153 | ! LOGICAL pctj(7,7,iwmx) |
---|
154 | cIM reg.dyn END |
---|
155 | ! ------ |
---|
156 | ! Output |
---|
157 | ! ------ |
---|
158 | |
---|
159 | REAL fq_isccp(npoints,7,7) ! the fraction of the model grid box covered by |
---|
160 | c ! each of the 49 ISCCP D level cloud types |
---|
161 | |
---|
162 | REAL totalcldarea(npoints) ! the fraction of model grid box columns |
---|
163 | c ! with cloud somewhere in them. This should |
---|
164 | c ! equal the sum over all entries of fq_isccp |
---|
165 | |
---|
166 | |
---|
167 | c ! The following three means are averages over the cloudy areas only. If no |
---|
168 | c ! clouds are in grid box all three quantities should equal zero. |
---|
169 | |
---|
170 | REAL meanptop(npoints) ! mean cloud top pressure (mb) - linear averaging |
---|
171 | c ! in cloud top pressure. |
---|
172 | |
---|
173 | REAL meantaucld(npoints) ! mean optical thickness |
---|
174 | c ! linear averaging in albedo performed. |
---|
175 | |
---|
176 | REAL boxtau(npoints,ncol) ! optical thickness in each column |
---|
177 | |
---|
178 | REAL boxptop(npoints,ncol) ! cloud top pressure (mb) in each column |
---|
179 | |
---|
180 | |
---|
181 | ! |
---|
182 | ! ------ |
---|
183 | ! Working variables added when program updated to mimic Mark Webb's PV-Wave code |
---|
184 | ! ------ |
---|
185 | |
---|
186 | REAL frac_out(npoints,ncol,nlev) ! boxes gridbox divided up into |
---|
187 | c ! Equivalent of BOX in original version, but |
---|
188 | c ! indexed by column then row, rather than |
---|
189 | c ! by row then column |
---|
190 | |
---|
191 | REAL tca(npoints,0:nlev) ! total cloud cover in each model level (fraction) |
---|
192 | c ! with extra layer of zeroes on top |
---|
193 | c ! in this version this just contains the values input |
---|
194 | c ! from cc but with an extra level |
---|
195 | REAL cca(npoints,nlev) ! convective cloud cover in each model level (fraction) |
---|
196 | c ! from conv |
---|
197 | |
---|
198 | REAL threshold(npoints,ncol) ! pointer to position in gridbox |
---|
199 | REAL maxocc(npoints,ncol) ! Flag for max overlapped conv cld |
---|
200 | REAL maxosc(npoints,ncol) ! Flag for max overlapped strat cld |
---|
201 | |
---|
202 | REAL boxpos(npoints,ncol) ! ordered pointer to position in gridbox |
---|
203 | |
---|
204 | REAL threshold_min(npoints,ncol) ! minimum value to define range in with new threshold |
---|
205 | c ! is chosen |
---|
206 | |
---|
207 | REAL dem(npoints,ncol),bb(npoints) ! working variables for 10.5 micron longwave |
---|
208 | c ! emissivity in part of |
---|
209 | c ! gridbox under consideration |
---|
210 | |
---|
211 | REAL ran(npoints) ! vector of random numbers |
---|
212 | REAL ptrop(npoints) |
---|
213 | REAL attrop(npoints) |
---|
214 | REAL attropmin (npoints) |
---|
215 | REAL atmax(npoints) |
---|
216 | REAL atmin(npoints) |
---|
217 | REAL btcmin(npoints) |
---|
218 | REAL transmax(npoints) |
---|
219 | |
---|
220 | INTEGER i,j,ilev,ibox,itrop(npoints) |
---|
221 | INTEGER ipres(npoints) |
---|
222 | INTEGER itau(npoints),ilev2 |
---|
223 | INTEGER acc(nlev,ncol) |
---|
224 | INTEGER match(npoints,nlev-1) |
---|
225 | INTEGER nmatch(npoints) |
---|
226 | INTEGER levmatch(npoints,ncol) |
---|
227 | |
---|
228 | c !variables needed for water vapor continuum absorption |
---|
229 | real fluxtop_clrsky(npoints),trans_layers_above_clrsky(npoints) |
---|
230 | real taumin(npoints) |
---|
231 | real dem_wv(npoints,nlev), wtmair, wtmh20, Navo, grav, pstd, t0 |
---|
232 | real press(npoints), dpress(npoints), atmden(npoints) |
---|
233 | real rvh20(npoints), wk(npoints), rhoave(npoints) |
---|
234 | real rh20s(npoints), rfrgn(npoints) |
---|
235 | real tmpexp(npoints),tauwv(npoints) |
---|
236 | |
---|
237 | character*1 cchar(6),cchar_realtops(6) |
---|
238 | integer icycle |
---|
239 | REAL tau(npoints,ncol) |
---|
240 | LOGICAL box_cloudy(npoints,ncol) |
---|
241 | REAL tb(npoints,ncol) |
---|
242 | REAL ptop(npoints,ncol) |
---|
243 | REAL emcld(npoints,ncol) |
---|
244 | REAL fluxtop(npoints,ncol) |
---|
245 | REAL trans_layers_above(npoints,ncol) |
---|
246 | real isccp_taumin,fluxtopinit(npoints),tauir(npoints) |
---|
247 | real meanalbedocld(npoints) |
---|
248 | REAL albedocld(npoints,ncol) |
---|
249 | real boxarea |
---|
250 | integer debug ! set to non-zero value to print out inputs |
---|
251 | c ! with step debug |
---|
252 | integer debugcol ! set to non-zero value to print out column |
---|
253 | c ! decomposition with step debugcol |
---|
254 | |
---|
255 | integer index1(npoints),num1,jj |
---|
256 | real rec2p13,tauchk |
---|
257 | |
---|
258 | character*10 ftn09 |
---|
259 | |
---|
260 | DATA isccp_taumin / 0.3 / |
---|
261 | DATA cchar / ' ','-','1','+','I','+'/ |
---|
262 | DATA cchar_realtops / ' ',' ','1','1','I','I'/ |
---|
263 | |
---|
264 | tauchk = -1.*log(0.9999999) |
---|
265 | rec2p13=1./2.13 |
---|
266 | |
---|
267 | ncolprint=0 |
---|
268 | |
---|
269 | cIM |
---|
270 | c PRINT*,' isccp_cloud_types npoints=',npoints |
---|
271 | c |
---|
272 | c if ( debug.ne.0 ) then |
---|
273 | c j=1 |
---|
274 | c write(6,'(a10)') 'j=' |
---|
275 | c write(6,'(8I10)') j |
---|
276 | c write(6,'(a10)') 'debug=' |
---|
277 | c write(6,'(8I10)') debug |
---|
278 | c write(6,'(a10)') 'debugcol=' |
---|
279 | c write(6,'(8I10)') debugcol |
---|
280 | c write(6,'(a10)') 'npoints=' |
---|
281 | c write(6,'(8I10)') npoints |
---|
282 | c write(6,'(a10)') 'nlev=' |
---|
283 | c write(6,'(8I10)') nlev |
---|
284 | c write(6,'(a10)') 'ncol=' |
---|
285 | c write(6,'(8I10)') ncol |
---|
286 | c write(6,'(a10)') 'top_height=' |
---|
287 | c write(6,'(8I10)') top_height |
---|
288 | c write(6,'(a10)') 'overlap=' |
---|
289 | c write(6,'(8I10)') overlap |
---|
290 | c write(6,'(a10)') 'emsfc_lw=' |
---|
291 | c write(6,'(8f10.2)') emsfc_lw |
---|
292 | c write(6,'(a10)') 'tautab=' |
---|
293 | c write(6,'(8f10.2)') tautab |
---|
294 | c write(6,'(a10)') 'invtau(1:100)=' |
---|
295 | c write(6,'(8i10)') (invtau(i),i=1,100) |
---|
296 | c do j=1,npoints,debug |
---|
297 | c write(6,'(a10)') 'j=' |
---|
298 | c write(6,'(8I10)') j |
---|
299 | c write(6,'(a10)') 'sunlit=' |
---|
300 | c write(6,'(8I10)') sunlit(j) |
---|
301 | c write(6,'(a10)') 'seed=' |
---|
302 | c write(6,'(8I10)') seed(j) |
---|
303 | c write(6,'(a10)') 'pfull=' |
---|
304 | c write(6,'(8f10.2)') (pfull(j,i),i=1,nlev) |
---|
305 | c write(6,'(a10)') 'phalf=' |
---|
306 | c write(6,'(8f10.2)') (phalf(j,i),i=1,nlev+1) |
---|
307 | c write(6,'(a10)') 'qv=' |
---|
308 | c write(6,'(8f10.3)') (qv(j,i),i=1,nlev) |
---|
309 | c write(6,'(a10)') 'cc=' |
---|
310 | c write(6,'(8f10.3)') (cc(j,i),i=1,nlev) |
---|
311 | c write(6,'(a10)') 'conv=' |
---|
312 | c write(6,'(8f10.2)') (conv(j,i),i=1,nlev) |
---|
313 | c write(6,'(a10)') 'dtau_s=' |
---|
314 | c write(6,'(8g12.5)') (dtau_s(j,i),i=1,nlev) |
---|
315 | c write(6,'(a10)') 'dtau_c=' |
---|
316 | c write(6,'(8f10.2)') (dtau_c(j,i),i=1,nlev) |
---|
317 | c write(6,'(a10)') 'skt=' |
---|
318 | c write(6,'(8f10.2)') skt(j) |
---|
319 | c write(6,'(a10)') 'at=' |
---|
320 | c write(6,'(8f10.2)') (at(j,i),i=1,nlev) |
---|
321 | c write(6,'(a10)') 'dem_s=' |
---|
322 | c write(6,'(8f10.3)') (dem_s(j,i),i=1,nlev) |
---|
323 | c write(6,'(a10)') 'dem_c=' |
---|
324 | c write(6,'(8f10.2)') (dem_c(j,i),i=1,nlev) |
---|
325 | c enddo |
---|
326 | c endif |
---|
327 | |
---|
328 | ! ---------------------------------------------------! |
---|
329 | |
---|
330 | ! assign 2d tca array using 1d input array cc |
---|
331 | |
---|
332 | do j=1,npoints |
---|
333 | tca(j,0)=0 |
---|
334 | enddo |
---|
335 | |
---|
336 | do ilev=1,nlev |
---|
337 | do j=1,npoints |
---|
338 | tca(j,ilev)=cc(j,ilev) |
---|
339 | enddo |
---|
340 | enddo |
---|
341 | |
---|
342 | ! assign 2d cca array using 1d input array conv |
---|
343 | |
---|
344 | do ilev=1,nlev |
---|
345 | cIM pas besoin do ibox=1,ncol |
---|
346 | do j=1,npoints |
---|
347 | cca(j,ilev)=conv(j,ilev) |
---|
348 | enddo |
---|
349 | cIM enddo |
---|
350 | enddo |
---|
351 | |
---|
352 | cIM |
---|
353 | ! do j=1, iwmx |
---|
354 | ! do l=1, 7 |
---|
355 | ! do k=1, 7 |
---|
356 | ! fq_dynreg(k,l,j) =0. |
---|
357 | ! nfq_dynreg(k,l,j) =0. |
---|
358 | ! enddo !k |
---|
359 | ! enddo !l |
---|
360 | ! enddo !j |
---|
361 | cIM |
---|
362 | cIM |
---|
363 | c if (ncolprint.ne.0) then |
---|
364 | c do j=1,npoints,1000 |
---|
365 | c write(6,'(a10)') 'j=' |
---|
366 | c write(6,'(8I10)') j |
---|
367 | c write (6,'(a)') 'seed:' |
---|
368 | c write (6,'(I3.2)') seed(j) |
---|
369 | |
---|
370 | c write (6,'(a)') 'tca_pp_rev:' |
---|
371 | c write (6,'(8f5.2)') |
---|
372 | c & ((tca(j,ilev)), |
---|
373 | c & ilev=1,nlev) |
---|
374 | |
---|
375 | c write (6,'(a)') 'cca_pp_rev:' |
---|
376 | c write (6,'(8f5.2)') |
---|
377 | c & ((cca(j,ilev),ibox=1,ncolprint),ilev=1,nlev) |
---|
378 | c enddo |
---|
379 | c endif |
---|
380 | |
---|
381 | if (top_height .eq. 1 .or. top_height .eq. 3) then |
---|
382 | |
---|
383 | do j=1,npoints |
---|
384 | ptrop(j)=5000. |
---|
385 | atmin(j) = 400. |
---|
386 | attropmin(j) = 400. |
---|
387 | atmax(j) = 0. |
---|
388 | attrop(j) = 120. |
---|
389 | itrop(j) = 1 |
---|
390 | enddo |
---|
391 | |
---|
392 | do 12 ilev=1,nlev |
---|
393 | do j=1,npoints |
---|
394 | if (pfull(j,ilev) .lt. 40000. .and. |
---|
395 | & pfull(j,ilev) .gt. 5000. .and. |
---|
396 | & at(j,ilev) .lt. attropmin(j)) then |
---|
397 | ptrop(j) = pfull(j,ilev) |
---|
398 | attropmin(j) = at(j,ilev) |
---|
399 | attrop(j) = attropmin(j) |
---|
400 | itrop(j)=ilev |
---|
401 | end if |
---|
402 | if (at(j,ilev) .gt. atmax(j)) atmax(j)=at(j,ilev) |
---|
403 | if (at(j,ilev) .lt. atmin(j)) atmin(j)=at(j,ilev) |
---|
404 | enddo |
---|
405 | 12 continue |
---|
406 | |
---|
407 | end if |
---|
408 | |
---|
409 | ! -----------------------------------------------------! |
---|
410 | |
---|
411 | ! ---------------------------------------------------! |
---|
412 | |
---|
413 | cIM |
---|
414 | c do 13 ilev=1,nlev |
---|
415 | cnum1=0 |
---|
416 | c do j=1,npoints |
---|
417 | c if (cc(j,ilev) .lt. 0. .or. cc(j,ilev) .gt. 1.) then |
---|
418 | c num1=num1+1 |
---|
419 | c index1(num1)=j |
---|
420 | c end if |
---|
421 | c enddo |
---|
422 | c do jj=1,num1 |
---|
423 | c j=index1(jj) |
---|
424 | c write(6,*) ' error = cloud fraction less than zero' |
---|
425 | c write(6,*) ' or ' |
---|
426 | c write(6,*) ' error = cloud fraction greater than 1' |
---|
427 | c write(6,*) 'value at point ',j,' is ',cc(j,ilev) |
---|
428 | c write(6,*) 'level ',ilev |
---|
429 | c STOP |
---|
430 | c enddo |
---|
431 | cnum1=0 |
---|
432 | c do j=1,npoints |
---|
433 | c if (conv(j,ilev) .lt. 0. .or. conv(j,ilev) .gt. 1.) then |
---|
434 | c num1=num1+1 |
---|
435 | c index1(num1)=j |
---|
436 | c end if |
---|
437 | c enddo |
---|
438 | c do jj=1,num1 |
---|
439 | c j=index1(jj) |
---|
440 | c write(6,*) |
---|
441 | c & ' error = convective cloud fraction less than zero' |
---|
442 | c write(6,*) ' or ' |
---|
443 | c write(6,*) |
---|
444 | c & ' error = convective cloud fraction greater than 1' |
---|
445 | c write(6,*) 'value at point ',j,' is ',conv(j,ilev) |
---|
446 | c write(6,*) 'level ',ilev |
---|
447 | c STOP |
---|
448 | c enddo |
---|
449 | |
---|
450 | cnum1=0 |
---|
451 | c do j=1,npoints |
---|
452 | c if (dtau_s(j,ilev) .lt. 0.) then |
---|
453 | c num1=num1+1 |
---|
454 | c index1(num1)=j |
---|
455 | c end if |
---|
456 | c enddo |
---|
457 | c do jj=1,num1 |
---|
458 | c j=index1(jj) |
---|
459 | c write(6,*) |
---|
460 | c & ' error = stratiform cloud opt. depth less than zero' |
---|
461 | c write(6,*) 'value at point ',j,' is ',dtau_s(j,ilev) |
---|
462 | c write(6,*) 'level ',ilev |
---|
463 | c STOP |
---|
464 | c enddo |
---|
465 | cnum1=0 |
---|
466 | c do j=1,npoints |
---|
467 | c if (dtau_c(j,ilev) .lt. 0.) then |
---|
468 | c num1=num1+1 |
---|
469 | c index1(num1)=j |
---|
470 | c end if |
---|
471 | c enddo |
---|
472 | c do jj=1,num1 |
---|
473 | c j=index1(jj) |
---|
474 | c write(6,*) |
---|
475 | c & ' error = convective cloud opt. depth less than zero' |
---|
476 | c write(6,*) 'value at point ',j,' is ',dtau_c(j,ilev) |
---|
477 | c write(6,*) 'level ',ilev |
---|
478 | c STOP |
---|
479 | c enddo |
---|
480 | |
---|
481 | cnum1=0 |
---|
482 | c do j=1,npoints |
---|
483 | c if (dem_s(j,ilev) .lt. 0. .or. dem_s(j,ilev) .gt. 1.) then |
---|
484 | c num1=num1+1 |
---|
485 | c index1(num1)=j |
---|
486 | c end if |
---|
487 | c enddo |
---|
488 | c do jj=1,num1 |
---|
489 | c j=index1(jj) |
---|
490 | c write(6,*) |
---|
491 | c & ' error = stratiform cloud emissivity less than zero' |
---|
492 | c write(6,*)'or' |
---|
493 | c write(6,*) |
---|
494 | c & ' error = stratiform cloud emissivity greater than 1' |
---|
495 | c write(6,*) 'value at point ',j,' is ',dem_s(j,ilev) |
---|
496 | c write(6,*) 'level ',ilev |
---|
497 | c STOP |
---|
498 | c enddo |
---|
499 | |
---|
500 | cnum1=0 |
---|
501 | c do j=1,npoints |
---|
502 | c if (dem_c(j,ilev) .lt. 0. .or. dem_c(j,ilev) .gt. 1.) then |
---|
503 | c num1=num1+1 |
---|
504 | c index1(num1)=j |
---|
505 | c end if |
---|
506 | c enddo |
---|
507 | c do jj=1,num1 |
---|
508 | c j=index1(jj) |
---|
509 | c write(6,*) |
---|
510 | c & ' error = convective cloud emissivity less than zero' |
---|
511 | c write(6,*)'or' |
---|
512 | c write(6,*) |
---|
513 | c & ' error = convective cloud emissivity greater than 1' |
---|
514 | c write (6,*) |
---|
515 | c & 'j=',j,'ilev=',ilev,'dem_c(j,ilev) =',dem_c(j,ilev) |
---|
516 | c STOP |
---|
517 | c enddo |
---|
518 | c13 continue |
---|
519 | |
---|
520 | |
---|
521 | do ibox=1,ncol |
---|
522 | do j=1,npoints |
---|
523 | boxpos(j,ibox)=(ibox-.5)/ncol |
---|
524 | enddo |
---|
525 | enddo |
---|
526 | |
---|
527 | ! ---------------------------------------------------! |
---|
528 | ! Initialise working variables |
---|
529 | ! ---------------------------------------------------! |
---|
530 | |
---|
531 | ! Initialised frac_out to zero |
---|
532 | |
---|
533 | do ibox=1,ncol |
---|
534 | do ilev=1,nlev |
---|
535 | do j=1,npoints |
---|
536 | frac_out(j,ibox,ilev)=0.0 |
---|
537 | enddo |
---|
538 | enddo |
---|
539 | enddo |
---|
540 | |
---|
541 | cIM |
---|
542 | c if (ncolprint.ne.0) then |
---|
543 | c write (6,'(a)') 'frac_out_pp_rev:' |
---|
544 | c do j=1,npoints,1000 |
---|
545 | c write(6,'(a10)') 'j=' |
---|
546 | c write(6,'(8I10)') j |
---|
547 | c write (6,'(8f5.2)') |
---|
548 | c & ((frac_out(j,ibox,ilev),ibox=1,ncolprint),ilev=1,nlev) |
---|
549 | |
---|
550 | c enddo |
---|
551 | c write (6,'(a)') 'ncol:' |
---|
552 | c write (6,'(I3)') ncol |
---|
553 | c endif |
---|
554 | c if (ncolprint.ne.0) then |
---|
555 | c write (6,'(a)') 'last_frac_pp:' |
---|
556 | c do j=1,npoints,1000 |
---|
557 | c write(6,'(a10)') 'j=' |
---|
558 | c write(6,'(8I10)') j |
---|
559 | c write (6,'(8f5.2)') (tca(j,0)) |
---|
560 | c enddo |
---|
561 | c endif |
---|
562 | |
---|
563 | ! ---------------------------------------------------! |
---|
564 | ! ALLOCATE CLOUD INTO BOXES, FOR NCOLUMNS, NLEVELS |
---|
565 | ! frac_out is the array that contains the information |
---|
566 | ! where 0 is no cloud, 1 is a stratiform cloud and 2 is a |
---|
567 | ! convective cloud |
---|
568 | |
---|
569 | !loop over vertical levels |
---|
570 | DO 200 ilev = 1,nlev |
---|
571 | |
---|
572 | ! Initialise threshold |
---|
573 | |
---|
574 | IF (ilev.eq.1) then |
---|
575 | ! If max overlap |
---|
576 | IF (overlap.eq.1) then |
---|
577 | ! select pixels spread evenly |
---|
578 | ! across the gridbox |
---|
579 | DO ibox=1,ncol |
---|
580 | do j=1,npoints |
---|
581 | threshold(j,ibox)=boxpos(j,ibox) |
---|
582 | enddo |
---|
583 | enddo |
---|
584 | ELSE |
---|
585 | DO ibox=1,ncol |
---|
586 | call ran0_vec(npoints,seed,ran) |
---|
587 | ! select random pixels from the non-convective |
---|
588 | ! part the gridbox ( some will be converted into |
---|
589 | ! convective pixels below ) |
---|
590 | do j=1,npoints |
---|
591 | threshold(j,ibox)= |
---|
592 | & cca(j,ilev)+(1-cca(j,ilev))*ran(j) |
---|
593 | enddo |
---|
594 | enddo |
---|
595 | ENDIF |
---|
596 | cIM |
---|
597 | c IF (ncolprint.ne.0) then |
---|
598 | c write (6,'(a)') 'threshold_nsf2:' |
---|
599 | c do j=1,npoints,1000 |
---|
600 | c write(6,'(a10)') 'j=' |
---|
601 | c write(6,'(8I10)') j |
---|
602 | c write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint) |
---|
603 | c enddo |
---|
604 | c ENDIF |
---|
605 | ENDIF |
---|
606 | |
---|
607 | c IF (ncolprint.ne.0) then |
---|
608 | c write (6,'(a)') 'ilev:' |
---|
609 | c write (6,'(I2)') ilev |
---|
610 | c ENDIF |
---|
611 | |
---|
612 | DO ibox=1,ncol |
---|
613 | |
---|
614 | ! All versions |
---|
615 | do j=1,npoints |
---|
616 | if (boxpos(j,ibox).le.cca(j,ilev)) then |
---|
617 | cIM REAL maxocc(j,ibox) = 1 |
---|
618 | maxocc(j,ibox) = 1.0 |
---|
619 | else |
---|
620 | cIM REAL maxocc(j,ibox) = 0 |
---|
621 | maxocc(j,ibox) = 0.0 |
---|
622 | end if |
---|
623 | enddo |
---|
624 | |
---|
625 | ! Max overlap |
---|
626 | if (overlap.eq.1) then |
---|
627 | do j=1,npoints |
---|
628 | threshold_min(j,ibox)=cca(j,ilev) |
---|
629 | cIM REAL maxosc(j,ibox)=1 |
---|
630 | maxosc(j,ibox)=1.0 |
---|
631 | enddo |
---|
632 | endif |
---|
633 | |
---|
634 | ! Random overlap |
---|
635 | if (overlap.eq.2) then |
---|
636 | do j=1,npoints |
---|
637 | threshold_min(j,ibox)=cca(j,ilev) |
---|
638 | cIM REAL maxosc(j,ibox)=0 |
---|
639 | maxosc(j,ibox)=0.0 |
---|
640 | enddo |
---|
641 | endif |
---|
642 | |
---|
643 | ! Max/Random overlap |
---|
644 | if (overlap.eq.3) then |
---|
645 | do j=1,npoints |
---|
646 | threshold_min(j,ibox)=max(cca(j,ilev), |
---|
647 | & min(tca(j,ilev-1),tca(j,ilev))) |
---|
648 | if (threshold(j,ibox) |
---|
649 | & .lt.min(tca(j,ilev-1),tca(j,ilev)) |
---|
650 | & .and.(threshold(j,ibox).gt.cca(j,ilev))) then |
---|
651 | cIM REAL maxosc(j,ibox)= 1 |
---|
652 | maxosc(j,ibox)= 1.0 |
---|
653 | else |
---|
654 | cIM REAL maxosc(j,ibox)= 0 |
---|
655 | maxosc(j,ibox)= 0.0 |
---|
656 | end if |
---|
657 | enddo |
---|
658 | endif |
---|
659 | |
---|
660 | ! Reset threshold |
---|
661 | call ran0_vec(npoints,seed,ran) |
---|
662 | |
---|
663 | do j=1,npoints |
---|
664 | threshold(j,ibox)= |
---|
665 | !if max overlapped conv cloud |
---|
666 | & maxocc(j,ibox) * ( |
---|
667 | & boxpos(j,ibox) |
---|
668 | & ) + |
---|
669 | !else |
---|
670 | & (1-maxocc(j,ibox)) * ( |
---|
671 | !if max overlapped strat cloud |
---|
672 | & (maxosc(j,ibox)) * ( |
---|
673 | !threshold=boxpos |
---|
674 | & threshold(j,ibox) |
---|
675 | & ) + |
---|
676 | !else |
---|
677 | & (1-maxosc(j,ibox)) * ( |
---|
678 | !threshold_min=random[thrmin,1] |
---|
679 | & threshold_min(j,ibox)+ |
---|
680 | & (1-threshold_min(j,ibox))*ran(j) |
---|
681 | & ) |
---|
682 | & ) |
---|
683 | enddo |
---|
684 | |
---|
685 | ENDDO ! ibox |
---|
686 | |
---|
687 | ! Fill frac_out with 1's where tca is greater than the threshold |
---|
688 | |
---|
689 | DO ibox=1,ncol |
---|
690 | do j=1,npoints |
---|
691 | if (tca(j,ilev).gt.threshold(j,ibox)) then |
---|
692 | cIM REAL frac_out(j,ibox,ilev)=1 |
---|
693 | frac_out(j,ibox,ilev)=1.0 |
---|
694 | else |
---|
695 | cIM REAL frac_out(j,ibox,ilev)=0 |
---|
696 | frac_out(j,ibox,ilev)=0.0 |
---|
697 | end if |
---|
698 | enddo |
---|
699 | ENDDO |
---|
700 | |
---|
701 | ! Code to partition boxes into startiform and convective parts |
---|
702 | ! goes here |
---|
703 | |
---|
704 | DO ibox=1,ncol |
---|
705 | do j=1,npoints |
---|
706 | if (threshold(j,ibox).le.cca(j,ilev)) then |
---|
707 | ! = 2 IF threshold le cca(j) |
---|
708 | cIM REAL frac_out(j,ibox,ilev) = 2 |
---|
709 | frac_out(j,ibox,ilev) = 2.0 |
---|
710 | else |
---|
711 | ! = the same IF NOT threshold le cca(j) |
---|
712 | frac_out(j,ibox,ilev) = frac_out(j,ibox,ilev) |
---|
713 | end if |
---|
714 | enddo |
---|
715 | ENDDO |
---|
716 | |
---|
717 | ! Set last_frac to tca at this level, so as to be tca |
---|
718 | ! from last level next time round |
---|
719 | |
---|
720 | cIM |
---|
721 | c if (ncolprint.ne.0) then |
---|
722 | |
---|
723 | c do j=1,npoints ,1000 |
---|
724 | c write(6,'(a10)') 'j=' |
---|
725 | c write(6,'(8I10)') j |
---|
726 | c write (6,'(a)') 'last_frac:' |
---|
727 | c write (6,'(8f5.2)') (tca(j,ilev-1)) |
---|
728 | |
---|
729 | c write (6,'(a)') 'cca:' |
---|
730 | c write (6,'(8f5.2)') (cca(j,ilev),ibox=1,ncolprint) |
---|
731 | |
---|
732 | c write (6,'(a)') 'max_overlap_cc:' |
---|
733 | c write (6,'(8f5.2)') (maxocc(j,ibox),ibox=1,ncolprint) |
---|
734 | |
---|
735 | c write (6,'(a)') 'max_overlap_sc:' |
---|
736 | c write (6,'(8f5.2)') (maxosc(j,ibox),ibox=1,ncolprint) |
---|
737 | |
---|
738 | c write (6,'(a)') 'threshold_min_nsf2:' |
---|
739 | c write (6,'(8f5.2)') (threshold_min(j,ibox),ibox=1,ncolprint) |
---|
740 | |
---|
741 | c write (6,'(a)') 'threshold_nsf2:' |
---|
742 | c write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint) |
---|
743 | |
---|
744 | c write (6,'(a)') 'frac_out_pp_rev:' |
---|
745 | c write (6,'(8f5.2)') |
---|
746 | c & ((frac_out(j,ibox,ilev2),ibox=1,ncolprint),ilev2=1,nlev) |
---|
747 | c enddo |
---|
748 | c endif |
---|
749 | |
---|
750 | 200 CONTINUE !loop over nlev |
---|
751 | |
---|
752 | ! |
---|
753 | ! ---------------------------------------------------! |
---|
754 | |
---|
755 | |
---|
756 | ! |
---|
757 | ! ---------------------------------------------------! |
---|
758 | ! COMPUTE CLOUD OPTICAL DEPTH FOR EACH COLUMN and |
---|
759 | ! put into vector tau |
---|
760 | |
---|
761 | !initialize tau and albedocld to zero |
---|
762 | do 15 ibox=1,ncol |
---|
763 | do j=1,npoints |
---|
764 | tau(j,ibox)=0. |
---|
765 | albedocld(j,ibox)=0. |
---|
766 | boxtau(j,ibox)=0. |
---|
767 | boxptop(j,ibox)=0. |
---|
768 | box_cloudy(j,ibox)=.false. |
---|
769 | enddo |
---|
770 | 15 continue |
---|
771 | |
---|
772 | !compute total cloud optical depth for each column |
---|
773 | do ilev=1,nlev |
---|
774 | !increment tau for each of the boxes |
---|
775 | do ibox=1,ncol |
---|
776 | do j=1,npoints |
---|
777 | cIM REAL if (frac_out(j,ibox,ilev).eq.1) then |
---|
778 | if (frac_out(j,ibox,ilev).eq.1.0) then |
---|
779 | tau(j,ibox)=tau(j,ibox) |
---|
780 | & + dtau_s(j,ilev) |
---|
781 | endif |
---|
782 | cIM REAL if (frac_out(j,ibox,ilev).eq.2) then |
---|
783 | if (frac_out(j,ibox,ilev).eq.2.0) then |
---|
784 | tau(j,ibox)=tau(j,ibox) |
---|
785 | & + dtau_c(j,ilev) |
---|
786 | end if |
---|
787 | enddo |
---|
788 | enddo ! ibox |
---|
789 | enddo ! ilev |
---|
790 | cIM |
---|
791 | c if (ncolprint.ne.0) then |
---|
792 | |
---|
793 | c do j=1,npoints ,1000 |
---|
794 | c write(6,'(a10)') 'j=' |
---|
795 | c write(6,'(8I10)') j |
---|
796 | c write(6,'(i2,1X,8(f7.2,1X))') |
---|
797 | c & ilev, |
---|
798 | c & (tau(j,ibox),ibox=1,ncolprint) |
---|
799 | c enddo |
---|
800 | c endif |
---|
801 | ! |
---|
802 | ! ---------------------------------------------------! |
---|
803 | |
---|
804 | |
---|
805 | |
---|
806 | ! |
---|
807 | ! ---------------------------------------------------! |
---|
808 | ! COMPUTE INFRARED BRIGHTNESS TEMPERUATRES |
---|
809 | ! AND CLOUD TOP TEMPERATURE SATELLITE SHOULD SEE |
---|
810 | ! |
---|
811 | ! again this is only done if top_height = 1 or 3 |
---|
812 | ! |
---|
813 | ! fluxtop is the 10.5 micron radiance at the top of the |
---|
814 | ! atmosphere |
---|
815 | ! trans_layers_above is the total transmissivity in the layers |
---|
816 | ! above the current layer |
---|
817 | ! fluxtop_clrsky(j) and trans_layers_above_clrsky(j) are the clear |
---|
818 | ! sky versions of these quantities. |
---|
819 | |
---|
820 | if (top_height .eq. 1 .or. top_height .eq. 3) then |
---|
821 | |
---|
822 | |
---|
823 | !---------------------------------------------------------------------- |
---|
824 | ! |
---|
825 | ! DO CLEAR SKY RADIANCE CALCULATION FIRST |
---|
826 | ! |
---|
827 | !compute water vapor continuum emissivity |
---|
828 | !this treatment follows Schwarkzopf and Ramasamy |
---|
829 | !JGR 1999,vol 104, pages 9467-9499. |
---|
830 | !the emissivity is calculated at a wavenumber of 955 cm-1, |
---|
831 | !or 10.47 microns |
---|
832 | wtmair = 28.9644 |
---|
833 | wtmh20 = 18.01534 |
---|
834 | Navo = 6.023E+23 |
---|
835 | grav = 9.806650E+02 |
---|
836 | pstd = 1.013250E+06 |
---|
837 | t0 = 296. |
---|
838 | cIM |
---|
839 | c if (ncolprint .ne. 0) |
---|
840 | c & write(6,*) 'ilev pw (kg/m2) tauwv(j) dem_wv' |
---|
841 | do 125 ilev=1,nlev |
---|
842 | do j=1,npoints |
---|
843 | !press and dpress are dyne/cm2 = Pascals *10 |
---|
844 | press(j) = pfull(j,ilev)*10. |
---|
845 | dpress(j) = (phalf(j,ilev+1)-phalf(j,ilev))*10 |
---|
846 | !atmden = g/cm2 = kg/m2 / 10 |
---|
847 | atmden(j) = dpress(j)/grav |
---|
848 | rvh20(j) = qv(j,ilev)*wtmair/wtmh20 |
---|
849 | wk(j) = rvh20(j)*Navo*atmden(j)/wtmair |
---|
850 | rhoave(j) = (press(j)/pstd)*(t0/at(j,ilev)) |
---|
851 | rh20s(j) = rvh20(j)*rhoave(j) |
---|
852 | rfrgn(j) = rhoave(j)-rh20s(j) |
---|
853 | tmpexp(j) = exp(-0.02*(at(j,ilev)-t0)) |
---|
854 | tauwv(j) = wk(j)*1.e-20*( |
---|
855 | & (0.0224697*rh20s(j)*tmpexp(j)) + |
---|
856 | & (3.41817e-7*rfrgn(j)) )*0.98 |
---|
857 | dem_wv(j,ilev) = 1. - exp( -1. * tauwv(j)) |
---|
858 | enddo |
---|
859 | cIM |
---|
860 | c if (ncolprint .ne. 0) then |
---|
861 | c do j=1,npoints ,1000 |
---|
862 | c write(6,'(a10)') 'j=' |
---|
863 | c write(6,'(8I10)') j |
---|
864 | c write(6,'(i2,1X,3(f8.3,3X))') ilev, |
---|
865 | c & qv(j,ilev)*(phalf(j,ilev+1)-phalf(j,ilev))/(grav/100.), |
---|
866 | c & tauwv(j),dem_wv(j,ilev) |
---|
867 | c enddo |
---|
868 | c endif |
---|
869 | 125 continue |
---|
870 | |
---|
871 | !initialize variables |
---|
872 | do j=1,npoints |
---|
873 | fluxtop_clrsky(j) = 0. |
---|
874 | trans_layers_above_clrsky(j)=1. |
---|
875 | enddo |
---|
876 | |
---|
877 | do ilev=1,nlev |
---|
878 | do j=1,npoints |
---|
879 | |
---|
880 | ! Black body emission at temperature of the layer |
---|
881 | |
---|
882 | bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. ) |
---|
883 | !bb(j)= 5.67e-8*at(j,ilev)**4 |
---|
884 | |
---|
885 | ! increase TOA flux by flux emitted from layer |
---|
886 | ! times total transmittance in layers above |
---|
887 | |
---|
888 | fluxtop_clrsky(j) = fluxtop_clrsky(j) |
---|
889 | & + dem_wv(j,ilev)*bb(j)*trans_layers_above_clrsky(j) |
---|
890 | |
---|
891 | ! update trans_layers_above with transmissivity |
---|
892 | ! from this layer for next time around loop |
---|
893 | |
---|
894 | trans_layers_above_clrsky(j)= |
---|
895 | & trans_layers_above_clrsky(j)*(1.-dem_wv(j,ilev)) |
---|
896 | |
---|
897 | |
---|
898 | enddo |
---|
899 | cIM |
---|
900 | c if (ncolprint.ne.0) then |
---|
901 | c do j=1,npoints ,1000 |
---|
902 | c write(6,'(a10)') 'j=' |
---|
903 | c write(6,'(8I10)') j |
---|
904 | c write (6,'(a)') 'ilev:' |
---|
905 | c write (6,'(I2)') ilev |
---|
906 | |
---|
907 | c write (6,'(a)') |
---|
908 | c & 'emiss_layer,100.*bb(j),100.*f,total_trans:' |
---|
909 | c write (6,'(4(f7.2,1X))') dem_wv(j,ilev),100.*bb(j), |
---|
910 | c & 100.*fluxtop_clrsky(j),trans_layers_above_clrsky(j) |
---|
911 | c enddo |
---|
912 | c endif |
---|
913 | |
---|
914 | enddo !loop over level |
---|
915 | |
---|
916 | do j=1,npoints |
---|
917 | !add in surface emission |
---|
918 | bb(j)=1/( exp(1307.27/skt(j)) - 1. ) |
---|
919 | !bb(j)=5.67e-8*skt(j)**4 |
---|
920 | |
---|
921 | fluxtop_clrsky(j) = fluxtop_clrsky(j) + emsfc_lw * bb(j) |
---|
922 | & * trans_layers_above_clrsky(j) |
---|
923 | enddo |
---|
924 | |
---|
925 | cIM |
---|
926 | c if (ncolprint.ne.0) then |
---|
927 | c do j=1,npoints ,1000 |
---|
928 | c write(6,'(a10)') 'j=' |
---|
929 | c write(6,'(8I10)') j |
---|
930 | c write (6,'(a)') 'id:' |
---|
931 | c write (6,'(a)') 'surface' |
---|
932 | |
---|
933 | c write (6,'(a)') 'emsfc,100.*bb(j),100.*f,total_trans:' |
---|
934 | c write (6,'(4(f7.2,1X))') emsfc_lw,100.*bb(j), |
---|
935 | c & 100.*fluxtop_clrsky(j), |
---|
936 | c & trans_layers_above_clrsky(j) |
---|
937 | c enddo |
---|
938 | c endif |
---|
939 | |
---|
940 | |
---|
941 | ! |
---|
942 | ! END OF CLEAR SKY CALCULATION |
---|
943 | ! |
---|
944 | !---------------------------------------------------------------- |
---|
945 | |
---|
946 | |
---|
947 | cIM |
---|
948 | c if (ncolprint.ne.0) then |
---|
949 | |
---|
950 | c do j=1,npoints ,1000 |
---|
951 | c write(6,'(a10)') 'j=' |
---|
952 | c write(6,'(8I10)') j |
---|
953 | c write (6,'(a)') 'ts:' |
---|
954 | c write (6,'(8f7.2)') (skt(j),ibox=1,ncolprint) |
---|
955 | |
---|
956 | c write (6,'(a)') 'ta_rev:' |
---|
957 | c write (6,'(8f7.2)') |
---|
958 | c & ((at(j,ilev2),ibox=1,ncolprint),ilev2=1,nlev) |
---|
959 | |
---|
960 | c enddo |
---|
961 | c endif |
---|
962 | !loop over columns |
---|
963 | do ibox=1,ncol |
---|
964 | do j=1,npoints |
---|
965 | fluxtop(j,ibox)=0. |
---|
966 | trans_layers_above(j,ibox)=1. |
---|
967 | enddo |
---|
968 | enddo |
---|
969 | |
---|
970 | do ilev=1,nlev |
---|
971 | do j=1,npoints |
---|
972 | ! Black body emission at temperature of the layer |
---|
973 | |
---|
974 | bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. ) |
---|
975 | !bb(j)= 5.67e-8*at(j,ilev)**4 |
---|
976 | enddo |
---|
977 | |
---|
978 | do ibox=1,ncol |
---|
979 | do j=1,npoints |
---|
980 | |
---|
981 | ! emissivity for point in this layer |
---|
982 | cIM REAL if (frac_out(j,ibox,ilev).eq.1) then |
---|
983 | if (frac_out(j,ibox,ilev).eq.1.0) then |
---|
984 | dem(j,ibox)= 1. - |
---|
985 | & ( (1. - dem_wv(j,ilev)) * (1. - dem_s(j,ilev)) ) |
---|
986 | cIM REAL else if (frac_out(j,ibox,ilev).eq.2) then |
---|
987 | else if (frac_out(j,ibox,ilev).eq.2.0) then |
---|
988 | dem(j,ibox)= 1. - |
---|
989 | & ( (1. - dem_wv(j,ilev)) * (1. - dem_c(j,ilev)) ) |
---|
990 | else |
---|
991 | dem(j,ibox)= dem_wv(j,ilev) |
---|
992 | end if |
---|
993 | |
---|
994 | |
---|
995 | ! increase TOA flux by flux emitted from layer |
---|
996 | ! times total transmittance in layers above |
---|
997 | |
---|
998 | fluxtop(j,ibox) = fluxtop(j,ibox) |
---|
999 | & + dem(j,ibox) * bb(j) |
---|
1000 | & * trans_layers_above(j,ibox) |
---|
1001 | |
---|
1002 | ! update trans_layers_above with transmissivity |
---|
1003 | ! from this layer for next time around loop |
---|
1004 | |
---|
1005 | trans_layers_above(j,ibox)= |
---|
1006 | & trans_layers_above(j,ibox)*(1.-dem(j,ibox)) |
---|
1007 | |
---|
1008 | enddo ! j |
---|
1009 | enddo ! ibox |
---|
1010 | |
---|
1011 | cIM |
---|
1012 | c if (ncolprint.ne.0) then |
---|
1013 | c do j=1,npoints,1000 |
---|
1014 | c write (6,'(a)') 'ilev:' |
---|
1015 | c write (6,'(I2)') ilev |
---|
1016 | |
---|
1017 | c write(6,'(a10)') 'j=' |
---|
1018 | c write(6,'(8I10)') j |
---|
1019 | c write (6,'(a)') 'emiss_layer:' |
---|
1020 | c write (6,'(8f7.2)') (dem(j,ibox),ibox=1,ncolprint) |
---|
1021 | |
---|
1022 | c write (6,'(a)') '100.*bb(j):' |
---|
1023 | c write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint) |
---|
1024 | |
---|
1025 | c write (6,'(a)') '100.*f:' |
---|
1026 | c write (6,'(8f7.2)') |
---|
1027 | c & (100.*fluxtop(j,ibox),ibox=1,ncolprint) |
---|
1028 | |
---|
1029 | c write (6,'(a)') 'total_trans:' |
---|
1030 | c write (6,'(8f7.2)') |
---|
1031 | c & (trans_layers_above(j,ibox),ibox=1,ncolprint) |
---|
1032 | c enddo |
---|
1033 | c endif |
---|
1034 | |
---|
1035 | enddo ! ilev |
---|
1036 | |
---|
1037 | |
---|
1038 | do j=1,npoints |
---|
1039 | !add in surface emission |
---|
1040 | bb(j)=1/( exp(1307.27/skt(j)) - 1. ) |
---|
1041 | !bb(j)=5.67e-8*skt(j)**4 |
---|
1042 | end do |
---|
1043 | |
---|
1044 | do ibox=1,ncol |
---|
1045 | do j=1,npoints |
---|
1046 | |
---|
1047 | !add in surface emission |
---|
1048 | |
---|
1049 | fluxtop(j,ibox) = fluxtop(j,ibox) |
---|
1050 | & + emsfc_lw * bb(j) |
---|
1051 | & * trans_layers_above(j,ibox) |
---|
1052 | |
---|
1053 | end do |
---|
1054 | end do |
---|
1055 | |
---|
1056 | cIM |
---|
1057 | c if (ncolprint.ne.0) then |
---|
1058 | |
---|
1059 | c do j=1,npoints ,1000 |
---|
1060 | c write(6,'(a10)') 'j=' |
---|
1061 | c write(6,'(8I10)') j |
---|
1062 | c write (6,'(a)') 'id:' |
---|
1063 | c write (6,'(a)') 'surface' |
---|
1064 | |
---|
1065 | c write (6,'(a)') 'emiss_layer:' |
---|
1066 | c write (6,'(8f7.2)') (dem(1,ibox),ibox=1,ncolprint) |
---|
1067 | |
---|
1068 | c write (6,'(a)') '100.*bb(j):' |
---|
1069 | c write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint) |
---|
1070 | |
---|
1071 | c write (6,'(a)') '100.*f:' |
---|
1072 | c write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint) |
---|
1073 | c end do |
---|
1074 | c endif |
---|
1075 | |
---|
1076 | !now that you have the top of atmosphere radiance account |
---|
1077 | !for ISCCP procedures to determine cloud top temperature |
---|
1078 | |
---|
1079 | !account for partially transmitting cloud recompute flux |
---|
1080 | !ISCCP would see assuming a single layer cloud |
---|
1081 | !note choice here of 2.13, as it is primarily ice |
---|
1082 | !clouds which have partial emissivity and need the |
---|
1083 | !adjustment performed in this section |
---|
1084 | ! |
---|
1085 | !If it turns out that the cloud brightness temperature |
---|
1086 | !is greater than 260K, then the liquid cloud conversion |
---|
1087 | !factor of 2.56 is used. |
---|
1088 | ! |
---|
1089 | !Note that this is discussed on pages 85-87 of |
---|
1090 | !the ISCCP D level documentation (Rossow et al. 1996) |
---|
1091 | |
---|
1092 | do j=1,npoints |
---|
1093 | !compute minimum brightness temperature and optical depth |
---|
1094 | btcmin(j) = 1. / ( exp(1307.27/(attrop(j)-5.)) - 1. ) |
---|
1095 | enddo |
---|
1096 | do ibox=1,ncol |
---|
1097 | do j=1,npoints |
---|
1098 | transmax(j) = (fluxtop(j,ibox)-btcmin(j)) |
---|
1099 | & /(fluxtop_clrsky(j)-btcmin(j)) |
---|
1100 | !note that the initial setting of tauir(j) is needed so that |
---|
1101 | !tauir(j) has a realistic value should the next if block be |
---|
1102 | !bypassed |
---|
1103 | tauir(j) = tau(j,ibox) * rec2p13 |
---|
1104 | taumin(j) = -1. * log(max(min(transmax(j),0.9999999),0.001)) |
---|
1105 | |
---|
1106 | enddo |
---|
1107 | |
---|
1108 | if (top_height .eq. 1) then |
---|
1109 | do j=1,npoints |
---|
1110 | if (transmax(j) .gt. 0.001 .and. |
---|
1111 | & transmax(j) .le. 0.9999999) then |
---|
1112 | fluxtopinit(j) = fluxtop(j,ibox) |
---|
1113 | tauir(j) = tau(j,ibox) *rec2p13 |
---|
1114 | endif |
---|
1115 | enddo |
---|
1116 | do icycle=1,2 |
---|
1117 | do j=1,npoints |
---|
1118 | if (tau(j,ibox) .gt. (tauchk )) then |
---|
1119 | if (transmax(j) .gt. 0.001 .and. |
---|
1120 | & transmax(j) .le. 0.9999999) then |
---|
1121 | emcld(j,ibox) = 1. - exp(-1. * tauir(j) ) |
---|
1122 | fluxtop(j,ibox) = fluxtopinit(j) - |
---|
1123 | & ((1.-emcld(j,ibox))*fluxtop_clrsky(j)) |
---|
1124 | fluxtop(j,ibox)=max(1.E-06, |
---|
1125 | & (fluxtop(j,ibox)/emcld(j,ibox))) |
---|
1126 | tb(j,ibox)= 1307.27 |
---|
1127 | & / (log(1. + (1./fluxtop(j,ibox)))) |
---|
1128 | if (tb(j,ibox) .gt. 260.) then |
---|
1129 | tauir(j) = tau(j,ibox) / 2.56 |
---|
1130 | end if |
---|
1131 | end if |
---|
1132 | end if |
---|
1133 | enddo |
---|
1134 | enddo |
---|
1135 | |
---|
1136 | endif |
---|
1137 | |
---|
1138 | do j=1,npoints |
---|
1139 | if (tau(j,ibox) .gt. (tauchk )) then |
---|
1140 | !cloudy box |
---|
1141 | tb(j,ibox)= 1307.27/ (log(1. + (1./fluxtop(j,ibox)))) |
---|
1142 | if (top_height.eq.1.and.tauir(j).lt.taumin(j)) then |
---|
1143 | tb(j,ibox) = attrop(j) - 5. |
---|
1144 | tau(j,ibox) = 2.13*taumin(j) |
---|
1145 | end if |
---|
1146 | else |
---|
1147 | !clear sky brightness temperature |
---|
1148 | tb(j,ibox) = 1307.27/(log(1.+(1./fluxtop_clrsky(j)))) |
---|
1149 | end if |
---|
1150 | enddo ! j |
---|
1151 | enddo ! ibox |
---|
1152 | |
---|
1153 | cIM |
---|
1154 | c if (ncolprint.ne.0) then |
---|
1155 | |
---|
1156 | c do j=1,npoints,1000 |
---|
1157 | c write(6,'(a10)') 'j=' |
---|
1158 | c write(6,'(8I10)') j |
---|
1159 | |
---|
1160 | c write (6,'(a)') 'attrop:' |
---|
1161 | c write (6,'(8f7.2)') (attrop(j)) |
---|
1162 | |
---|
1163 | c write (6,'(a)') 'btcmin:' |
---|
1164 | c write (6,'(8f7.2)') (btcmin(j)) |
---|
1165 | |
---|
1166 | c write (6,'(a)') 'fluxtop_clrsky*100:' |
---|
1167 | c write (6,'(8f7.2)') |
---|
1168 | c & (100.*fluxtop_clrsky(j)) |
---|
1169 | |
---|
1170 | c write (6,'(a)') '100.*f_adj:' |
---|
1171 | c write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint) |
---|
1172 | |
---|
1173 | c write (6,'(a)') 'transmax:' |
---|
1174 | c write (6,'(8f7.2)') (transmax(ibox),ibox=1,ncolprint) |
---|
1175 | |
---|
1176 | c write (6,'(a)') 'tau:' |
---|
1177 | c write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint) |
---|
1178 | |
---|
1179 | c write (6,'(a)') 'emcld:' |
---|
1180 | c write (6,'(8f7.2)') (emcld(j,ibox),ibox=1,ncolprint) |
---|
1181 | |
---|
1182 | c write (6,'(a)') 'total_trans:' |
---|
1183 | c write (6,'(8f7.2)') |
---|
1184 | c & (trans_layers_above(j,ibox),ibox=1,ncolprint) |
---|
1185 | |
---|
1186 | c write (6,'(a)') 'total_emiss:' |
---|
1187 | c write (6,'(8f7.2)') |
---|
1188 | c & (1.0-trans_layers_above(j,ibox),ibox=1,ncolprint) |
---|
1189 | |
---|
1190 | c write (6,'(a)') 'total_trans:' |
---|
1191 | c write (6,'(8f7.2)') |
---|
1192 | c & (trans_layers_above(j,ibox),ibox=1,ncolprint) |
---|
1193 | |
---|
1194 | c write (6,'(a)') 'ppout:' |
---|
1195 | c write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint) |
---|
1196 | c enddo ! j |
---|
1197 | c endif |
---|
1198 | |
---|
1199 | end if |
---|
1200 | |
---|
1201 | ! ---------------------------------------------------! |
---|
1202 | |
---|
1203 | ! |
---|
1204 | ! ---------------------------------------------------! |
---|
1205 | ! DETERMINE CLOUD TOP PRESSURE |
---|
1206 | ! |
---|
1207 | ! again the 2 methods differ according to whether |
---|
1208 | ! or not you use the physical cloud top pressure (top_height = 2) |
---|
1209 | ! or the radiatively determined cloud top pressure (top_height = 1 or 3) |
---|
1210 | ! |
---|
1211 | |
---|
1212 | !compute cloud top pressure |
---|
1213 | do 30 ibox=1,ncol |
---|
1214 | !segregate according to optical thickness |
---|
1215 | if (top_height .eq. 1 .or. top_height .eq. 3) then |
---|
1216 | !find level whose temperature |
---|
1217 | !most closely matches brightness temperature |
---|
1218 | do j=1,npoints |
---|
1219 | nmatch(j)=0 |
---|
1220 | enddo |
---|
1221 | do 29 ilev=1,nlev-1 |
---|
1222 | !cdir nodep |
---|
1223 | do j=1,npoints |
---|
1224 | if ((at(j,ilev) .ge. tb(j,ibox) .and. |
---|
1225 | & at(j,ilev+1) .lt. tb(j,ibox)) .or. |
---|
1226 | & (at(j,ilev) .le. tb(j,ibox) .and. |
---|
1227 | & at(j,ilev+1) .gt. tb(j,ibox))) then |
---|
1228 | |
---|
1229 | nmatch(j)=nmatch(j)+1 |
---|
1230 | if(abs(at(j,ilev)-tb(j,ibox)) .lt. |
---|
1231 | & abs(at(j,ilev+1)-tb(j,ibox))) then |
---|
1232 | match(j,nmatch(j))=ilev |
---|
1233 | else |
---|
1234 | match(j,nmatch(j))=ilev+1 |
---|
1235 | end if |
---|
1236 | end if |
---|
1237 | enddo |
---|
1238 | 29 continue |
---|
1239 | |
---|
1240 | do j=1,npoints |
---|
1241 | if (nmatch(j) .ge. 1) then |
---|
1242 | ptop(j,ibox)=pfull(j,match(j,nmatch(j))) |
---|
1243 | levmatch(j,ibox)=match(j,nmatch(j)) |
---|
1244 | else |
---|
1245 | if (tb(j,ibox) .lt. atmin(j)) then |
---|
1246 | ptop(j,ibox)=ptrop(j) |
---|
1247 | levmatch(j,ibox)=itrop(j) |
---|
1248 | end if |
---|
1249 | if (tb(j,ibox) .gt. atmax(j)) then |
---|
1250 | ptop(j,ibox)=pfull(j,nlev) |
---|
1251 | levmatch(j,ibox)=nlev |
---|
1252 | end if |
---|
1253 | end if |
---|
1254 | enddo ! j |
---|
1255 | |
---|
1256 | else ! if (top_height .eq. 1 .or. top_height .eq. 3) |
---|
1257 | |
---|
1258 | do j=1,npoints |
---|
1259 | ptop(j,ibox)=0. |
---|
1260 | enddo |
---|
1261 | do ilev=1,nlev |
---|
1262 | do j=1,npoints |
---|
1263 | if ((ptop(j,ibox) .eq. 0. ) |
---|
1264 | cIM & .and.(frac_out(j,ibox,ilev) .ne. 0)) then |
---|
1265 | & .and.(frac_out(j,ibox,ilev) .ne. 0.0)) then |
---|
1266 | ptop(j,ibox)=pfull(j,ilev) |
---|
1267 | levmatch(j,ibox)=ilev |
---|
1268 | end if |
---|
1269 | end do |
---|
1270 | end do |
---|
1271 | end if |
---|
1272 | |
---|
1273 | do j=1,npoints |
---|
1274 | if (tau(j,ibox) .le. (tauchk )) then |
---|
1275 | ptop(j,ibox)=0. |
---|
1276 | levmatch(j,ibox)=0 |
---|
1277 | endif |
---|
1278 | enddo |
---|
1279 | |
---|
1280 | 30 continue |
---|
1281 | |
---|
1282 | ! |
---|
1283 | ! |
---|
1284 | ! ---------------------------------------------------! |
---|
1285 | |
---|
1286 | |
---|
1287 | ! |
---|
1288 | ! ---------------------------------------------------! |
---|
1289 | ! DETERMINE ISCCP CLOUD TYPE FREQUENCIES |
---|
1290 | ! |
---|
1291 | ! Now that ptop and tau have been determined, |
---|
1292 | ! determine amount of each of the 49 ISCCP cloud |
---|
1293 | ! types |
---|
1294 | ! |
---|
1295 | ! Also compute grid box mean cloud top pressure and |
---|
1296 | ! optical thickness. The mean cloud top pressure and |
---|
1297 | ! optical thickness are averages over the cloudy |
---|
1298 | ! area only. The mean cloud top pressure is a linear |
---|
1299 | ! average of the cloud top pressures. The mean cloud |
---|
1300 | ! optical thickness is computed by converting optical |
---|
1301 | ! thickness to an albedo, averaging in albedo units, |
---|
1302 | ! then converting the average albedo back to a mean |
---|
1303 | ! optical thickness. |
---|
1304 | ! |
---|
1305 | |
---|
1306 | !compute isccp frequencies |
---|
1307 | |
---|
1308 | !reset frequencies |
---|
1309 | do 38 ilev=1,7 |
---|
1310 | do 38 ilev2=1,7 |
---|
1311 | do j=1,npoints ! |
---|
1312 | fq_isccp(j,ilev,ilev2)=0. |
---|
1313 | enddo |
---|
1314 | 38 continue |
---|
1315 | |
---|
1316 | !reset variables need for averaging cloud properties |
---|
1317 | do j=1,npoints |
---|
1318 | totalcldarea(j) = 0. |
---|
1319 | meanalbedocld(j) = 0. |
---|
1320 | meanptop(j) = 0. |
---|
1321 | meantaucld(j) = 0. |
---|
1322 | enddo ! j |
---|
1323 | |
---|
1324 | boxarea = 1./real(ncol) |
---|
1325 | |
---|
1326 | !determine optical depth category |
---|
1327 | cIM do 39 j=1,npoints |
---|
1328 | cIM do ibox=1,ncol |
---|
1329 | do 39 ibox=1,ncol |
---|
1330 | do j=1,npoints |
---|
1331 | |
---|
1332 | cIM |
---|
1333 | c CALL CPU_time(t1) |
---|
1334 | cIM |
---|
1335 | |
---|
1336 | if (tau(j,ibox) .gt. (tauchk ) |
---|
1337 | & .and. ptop(j,ibox) .gt. 0.) then |
---|
1338 | box_cloudy(j,ibox)=.true. |
---|
1339 | endif |
---|
1340 | |
---|
1341 | cIM |
---|
1342 | c CALL CPU_time(t2) |
---|
1343 | c print*,'IF tau t2 - t1',t2 - t1 |
---|
1344 | |
---|
1345 | c CALL CPU_time(t1) |
---|
1346 | cIM |
---|
1347 | |
---|
1348 | if (box_cloudy(j,ibox)) then |
---|
1349 | |
---|
1350 | ! totalcldarea always diagnosed day or night |
---|
1351 | totalcldarea(j) = totalcldarea(j) + boxarea |
---|
1352 | |
---|
1353 | if (sunlit(j).eq.1) then |
---|
1354 | |
---|
1355 | ! tau diagnostics only with sunlight |
---|
1356 | |
---|
1357 | boxtau(j,ibox) = tau(j,ibox) |
---|
1358 | |
---|
1359 | !convert optical thickness to albedo |
---|
1360 | albedocld(j,ibox) |
---|
1361 | & =real(invtau(min(nint(100.*tau(j,ibox)),45000))) |
---|
1362 | |
---|
1363 | !contribute to averaging |
---|
1364 | meanalbedocld(j) = meanalbedocld(j) |
---|
1365 | & +albedocld(j,ibox)*boxarea |
---|
1366 | |
---|
1367 | endif |
---|
1368 | |
---|
1369 | endif |
---|
1370 | |
---|
1371 | cIM |
---|
1372 | c CALL CPU_time(t2) |
---|
1373 | c print*,'IF box_cloudy t2 - t1',t2 - t1 |
---|
1374 | |
---|
1375 | c CALL CPU_time(t1) |
---|
1376 | cIM BEG |
---|
1377 | cIM !convert ptop to millibars |
---|
1378 | ptop(j,ibox)=ptop(j,ibox) / 100. |
---|
1379 | |
---|
1380 | cIM !save for output cloud top pressure and optical thickness |
---|
1381 | boxptop(j,ibox) = ptop(j,ibox) |
---|
1382 | cIM END |
---|
1383 | |
---|
1384 | cIM BEG |
---|
1385 | !reset itau(j), ipres(j) |
---|
1386 | itau(j) = 0 |
---|
1387 | ipres(j) = 0 |
---|
1388 | |
---|
1389 | if (tau(j,ibox) .lt. isccp_taumin) then |
---|
1390 | itau(j)=1 |
---|
1391 | else if (tau(j,ibox) .ge. isccp_taumin |
---|
1392 | & |
---|
1393 | & .and. tau(j,ibox) .lt. 1.3) then |
---|
1394 | itau(j)=2 |
---|
1395 | else if (tau(j,ibox) .ge. 1.3 |
---|
1396 | & .and. tau(j,ibox) .lt. 3.6) then |
---|
1397 | itau(j)=3 |
---|
1398 | else if (tau(j,ibox) .ge. 3.6 |
---|
1399 | & .and. tau(j,ibox) .lt. 9.4) then |
---|
1400 | itau(j)=4 |
---|
1401 | else if (tau(j,ibox) .ge. 9.4 |
---|
1402 | & .and. tau(j,ibox) .lt. 23.) then |
---|
1403 | itau(j)=5 |
---|
1404 | else if (tau(j,ibox) .ge. 23. |
---|
1405 | & .and. tau(j,ibox) .lt. 60.) then |
---|
1406 | itau(j)=6 |
---|
1407 | else if (tau(j,ibox) .ge. 60.) then |
---|
1408 | itau(j)=7 |
---|
1409 | end if |
---|
1410 | |
---|
1411 | !determine cloud top pressure category |
---|
1412 | if ( ptop(j,ibox) .gt. 0. |
---|
1413 | & .and.ptop(j,ibox) .lt. 180.) then |
---|
1414 | ipres(j)=1 |
---|
1415 | else if(ptop(j,ibox) .ge. 180. |
---|
1416 | & .and.ptop(j,ibox) .lt. 310.) then |
---|
1417 | ipres(j)=2 |
---|
1418 | else if(ptop(j,ibox) .ge. 310. |
---|
1419 | & .and.ptop(j,ibox) .lt. 440.) then |
---|
1420 | ipres(j)=3 |
---|
1421 | else if(ptop(j,ibox) .ge. 440. |
---|
1422 | & .and.ptop(j,ibox) .lt. 560.) then |
---|
1423 | ipres(j)=4 |
---|
1424 | else if(ptop(j,ibox) .ge. 560. |
---|
1425 | & .and.ptop(j,ibox) .lt. 680.) then |
---|
1426 | ipres(j)=5 |
---|
1427 | else if(ptop(j,ibox) .ge. 680. |
---|
1428 | & .and.ptop(j,ibox) .lt. 800.) then |
---|
1429 | ipres(j)=6 |
---|
1430 | else if(ptop(j,ibox) .ge. 800.) then |
---|
1431 | ipres(j)=7 |
---|
1432 | end if |
---|
1433 | cIM END |
---|
1434 | |
---|
1435 | if (sunlit(j).eq.1 .or. top_height .eq. 3) then |
---|
1436 | |
---|
1437 | cIM !convert ptop to millibars |
---|
1438 | cIM ptop(j,ibox)=ptop(j,ibox) / 100. |
---|
1439 | |
---|
1440 | cIM !save for output cloud top pressure and optical thickness |
---|
1441 | cIM boxptop(j,ibox) = ptop(j,ibox) |
---|
1442 | |
---|
1443 | if (box_cloudy(j,ibox)) then |
---|
1444 | |
---|
1445 | meanptop(j) = meanptop(j) + ptop(j,ibox)*boxarea |
---|
1446 | |
---|
1447 | cIM !reset itau(j), ipres(j) |
---|
1448 | cIM itau(j) = 0 |
---|
1449 | cIM ipres(j) = 0 |
---|
1450 | |
---|
1451 | c if (tau(j,ibox) .lt. isccp_taumin) then |
---|
1452 | c itau(j)=1 |
---|
1453 | c else if (tau(j,ibox) .ge. isccp_taumin |
---|
1454 | c & |
---|
1455 | c & .and. tau(j,ibox) .lt. 1.3) then |
---|
1456 | c itau(j)=2 |
---|
1457 | c else if (tau(j,ibox) .ge. 1.3 |
---|
1458 | c & .and. tau(j,ibox) .lt. 3.6) then |
---|
1459 | c itau(j)=3 |
---|
1460 | c else if (tau(j,ibox) .ge. 3.6 |
---|
1461 | c & .and. tau(j,ibox) .lt. 9.4) then |
---|
1462 | c itau(j)=4 |
---|
1463 | c else if (tau(j,ibox) .ge. 9.4 |
---|
1464 | c & .and. tau(j,ibox) .lt. 23.) then |
---|
1465 | c itau(j)=5 |
---|
1466 | c else if (tau(j,ibox) .ge. 23. |
---|
1467 | c & .and. tau(j,ibox) .lt. 60.) then |
---|
1468 | c itau(j)=6 |
---|
1469 | c else if (tau(j,ibox) .ge. 60.) then |
---|
1470 | c itau(j)=7 |
---|
1471 | c end if |
---|
1472 | |
---|
1473 | c !determine cloud top pressure category |
---|
1474 | c if ( ptop(j,ibox) .gt. 0. |
---|
1475 | c & .and.ptop(j,ibox) .lt. 180.) then |
---|
1476 | c ipres(j)=1 |
---|
1477 | c else if(ptop(j,ibox) .ge. 180. |
---|
1478 | c & .and.ptop(j,ibox) .lt. 310.) then |
---|
1479 | c ipres(j)=2 |
---|
1480 | c else if(ptop(j,ibox) .ge. 310. |
---|
1481 | c & .and.ptop(j,ibox) .lt. 440.) then |
---|
1482 | c ipres(j)=3 |
---|
1483 | c else if(ptop(j,ibox) .ge. 440. |
---|
1484 | c & .and.ptop(j,ibox) .lt. 560.) then |
---|
1485 | c ipres(j)=4 |
---|
1486 | c else if(ptop(j,ibox) .ge. 560. |
---|
1487 | c & .and.ptop(j,ibox) .lt. 680.) then |
---|
1488 | c ipres(j)=5 |
---|
1489 | c else if(ptop(j,ibox) .ge. 680. |
---|
1490 | c & .and.ptop(j,ibox) .lt. 800.) then |
---|
1491 | c ipres(j)=6 |
---|
1492 | c else if(ptop(j,ibox) .ge. 800.) then |
---|
1493 | c ipres(j)=7 |
---|
1494 | c end if |
---|
1495 | |
---|
1496 | !update frequencies |
---|
1497 | if(ipres(j) .gt. 0.and.itau(j) .gt. 0) then |
---|
1498 | fq_isccp(j,itau(j),ipres(j))= |
---|
1499 | & fq_isccp(j,itau(j),ipres(j))+ boxarea |
---|
1500 | end if |
---|
1501 | |
---|
1502 | cIM calcul stats regime dynamique BEG |
---|
1503 | ! iw(j) = int((w(j)-wmin)/pas_w) +1 |
---|
1504 | ! pctj(itau(j),ipres(j),iw(j))=.FALSE. |
---|
1505 | ! !update frequencies W500 |
---|
1506 | ! if (pct_ocean(j)) then |
---|
1507 | ! if (ipres(j) .gt. 0.and.itau(j) .gt. 0) then |
---|
1508 | ! if (iw(j) .gt. int(wmin).and.iw(j) .le. iwmx) then |
---|
1509 | c print*,' ISCCP iw=',iw(j),j |
---|
1510 | ! fq_dynreg(itau(j),ipres(j),iw(j))= |
---|
1511 | ! & fq_dynreg(itau(j),ipres(j),iw(j))+ |
---|
1512 | ! & boxarea |
---|
1513 | c & fq_isccp(j,itau(j),ipres(j)) |
---|
1514 | ! pctj(itau(j),ipres(j),iw(j))=.TRUE. |
---|
1515 | c nfq_dynreg(itau(j),ipres(j),iw(j))= |
---|
1516 | c & nfq_dynreg(itau(j),ipres(j),iw(j))+1. |
---|
1517 | ! end if |
---|
1518 | ! end if |
---|
1519 | ! end if |
---|
1520 | cIM calcul stats regime dynamique END |
---|
1521 | end if !IM boxcloudy |
---|
1522 | |
---|
1523 | end if !IM sunlit |
---|
1524 | |
---|
1525 | cIM |
---|
1526 | c CALL CPU_time(t2) |
---|
1527 | c print*,'IF sunlit boxcloudy t2 - t1',t2 - t1 |
---|
1528 | cIM |
---|
1529 | enddo !IM ibox/j |
---|
1530 | |
---|
1531 | |
---|
1532 | cIM ajout stats s/ W500 BEG |
---|
1533 | cIM ajout stats s/ W500 END |
---|
1534 | |
---|
1535 | c if(j.EQ.6722) then |
---|
1536 | c print*,' ISCCP',w(j),iw(j),ipres(j),itau(j) |
---|
1537 | c endif |
---|
1538 | |
---|
1539 | ! if (pct_ocean(j)) then |
---|
1540 | ! if (ipres(j) .gt. 0.and.itau(j) .gt. 0) then |
---|
1541 | ! if (iw(j) .gt. int(wmin).and.iw(j) .le. iwmx) then |
---|
1542 | ! if(pctj(itau(j),ipres(j),iw(j))) THEN |
---|
1543 | ! nfq_dynreg(itau(j),ipres(j),iw(j))= |
---|
1544 | ! & nfq_dynreg(itau(j),ipres(j),iw(j))+1. |
---|
1545 | c if(itau(j).EQ.4.AND.ipres(j).EQ.2.AND. |
---|
1546 | c & iw(j).EQ.10) then |
---|
1547 | c PRINT*,' isccp AVANT', |
---|
1548 | c & nfq_dynreg(itau(j),ipres(j),iw(j)), |
---|
1549 | c & fq_dynreg(itau(j),ipres(j),iw(j)) |
---|
1550 | c endif |
---|
1551 | ! endif |
---|
1552 | ! endif |
---|
1553 | ! endif |
---|
1554 | ! endif |
---|
1555 | 39 continue !IM j/ibox |
---|
1556 | |
---|
1557 | !compute mean cloud properties |
---|
1558 | do j=1,npoints |
---|
1559 | if (totalcldarea(j) .gt. 0.) then |
---|
1560 | meanptop(j) = meanptop(j) / totalcldarea(j) |
---|
1561 | if (sunlit(j).eq.1) then |
---|
1562 | meanalbedocld(j) = meanalbedocld(j) / totalcldarea(j) |
---|
1563 | meantaucld(j) = tautab(min(255,max(1,nint(meanalbedocld(j))))) |
---|
1564 | end if |
---|
1565 | end if |
---|
1566 | enddo ! j |
---|
1567 | ! |
---|
1568 | cIM ajout stats s/ W500 BEG |
---|
1569 | ! do nw = 1, iwmx |
---|
1570 | ! do l = 1, 7 |
---|
1571 | ! do k = 1, 7 |
---|
1572 | ! if (nfq_dynreg(k,l,nw).GT.0.) then |
---|
1573 | ! fq_dynreg(k,l,nw) = fq_dynreg(k,l,nw)/nfq_dynreg(k,l,nw) |
---|
1574 | c if(k.EQ.4.AND.l.EQ.2.AND.nw.EQ.10) then |
---|
1575 | c print*,' isccp APRES',nfq_dynreg(k,l,nw), |
---|
1576 | c & fq_dynreg(k,l,nw) |
---|
1577 | c endif |
---|
1578 | ! else |
---|
1579 | ! if(fq_dynreg(k,l,nw).NE.0.) then |
---|
1580 | ! print*,'nfq_dynreg = 0 tau,pc,nw',k,l,nw,fq_dynreg(k,l,nw) |
---|
1581 | ! endif |
---|
1582 | c fq_dynreg(k,l,nw) = -1.E+20 |
---|
1583 | c nfq_dynreg(k,l,nw) = 1.E+20 |
---|
1584 | ! end if |
---|
1585 | ! enddo !k |
---|
1586 | ! enddo !l |
---|
1587 | ! enddo !nw |
---|
1588 | cIM ajout stats s/ W500 END |
---|
1589 | ! ---------------------------------------------------! |
---|
1590 | |
---|
1591 | ! ---------------------------------------------------! |
---|
1592 | ! OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM |
---|
1593 | ! |
---|
1594 | !cIM |
---|
1595 | c if (debugcol.ne.0) then |
---|
1596 | ! |
---|
1597 | c do j=1,npoints,debugcol |
---|
1598 | |
---|
1599 | c !produce character output |
---|
1600 | c do ilev=1,nlev |
---|
1601 | c do ibox=1,ncol |
---|
1602 | c acc(ilev,ibox)=0 |
---|
1603 | c enddo |
---|
1604 | c enddo |
---|
1605 | |
---|
1606 | c do ilev=1,nlev |
---|
1607 | c do ibox=1,ncol |
---|
1608 | c acc(ilev,ibox)=frac_out(j,ibox,ilev)*2 |
---|
1609 | c if (levmatch(j,ibox) .eq. ilev) |
---|
1610 | c & acc(ilev,ibox)=acc(ilev,ibox)+1 |
---|
1611 | c enddo |
---|
1612 | c enddo |
---|
1613 | |
---|
1614 | !print test |
---|
1615 | |
---|
1616 | c write(ftn09,11) j |
---|
1617 | c11 format('ftn09.',i4.4) |
---|
1618 | c open(9, FILE=ftn09, FORM='FORMATTED') |
---|
1619 | |
---|
1620 | c write(9,'(a1)') ' ' |
---|
1621 | c write(9,'(10i5)') |
---|
1622 | c & (ilev,ilev=5,nlev,5) |
---|
1623 | c write(9,'(a1)') ' ' |
---|
1624 | |
---|
1625 | c do ibox=1,ncol |
---|
1626 | c write(9,'(40(a1),1x,40(a1))') |
---|
1627 | c & (cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev) |
---|
1628 | c & ,(cchar(acc(ilev,ibox)+1),ilev=1,nlev) |
---|
1629 | c end do |
---|
1630 | c close(9) |
---|
1631 | c |
---|
1632 | cIM |
---|
1633 | c if (ncolprint.ne.0) then |
---|
1634 | c write(6,'(a1)') ' ' |
---|
1635 | c write(6,'(a2,1X,5(a7,1X),a50)') |
---|
1636 | c & 'ilev', |
---|
1637 | c & 'pfull','at', |
---|
1638 | c & 'cc*100','dem_s','dtau_s', |
---|
1639 | c & 'cchar' |
---|
1640 | |
---|
1641 | ! do 4012 ilev=1,nlev |
---|
1642 | ! write(6,'(60i2)') (box(i,ilev),i=1,ncolprint) |
---|
1643 | ! write(6,'(i2,1X,5(f7.2,1X),50(a1))') |
---|
1644 | ! & ilev, |
---|
1645 | ! & pfull(j,ilev)/100.,at(j,ilev), |
---|
1646 | ! & cc(j,ilev)*100.0,dem_s(j,ilev),dtau_s(j,ilev) |
---|
1647 | ! & ,(cchar(acc(ilev,ibox)+1),ibox=1,ncolprint) |
---|
1648 | !4012 continue |
---|
1649 | c write (6,'(a)') 'skt(j):' |
---|
1650 | c write (6,'(8f7.2)') skt(j) |
---|
1651 | |
---|
1652 | c write (6,'(8I7)') (ibox,ibox=1,ncolprint) |
---|
1653 | |
---|
1654 | c write (6,'(a)') 'tau:' |
---|
1655 | c write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint) |
---|
1656 | |
---|
1657 | c write (6,'(a)') 'tb:' |
---|
1658 | c write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint) |
---|
1659 | |
---|
1660 | c write (6,'(a)') 'ptop:' |
---|
1661 | c write (6,'(8f7.2)') (ptop(j,ibox),ibox=1,ncolprint) |
---|
1662 | c endif |
---|
1663 | |
---|
1664 | c enddo |
---|
1665 | |
---|
1666 | c end if |
---|
1667 | |
---|
1668 | return |
---|
1669 | end |
---|