GCC Code Coverage Report


Directory: ./
File: phys/isccp_cloud_types.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 303 0.0%
Branches: 0 300 0.0%

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