Context Navigation

isccp_cloud_types.F @ 466

Last change on this file since 466 was 466, checked in by (none), 21 years ago
This commit was manufactured by cvs2svn to create branch 'rel-LF'.
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 53.6 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format