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