Context Navigation

isccp_cloud_types.F @ 594

Last change on this file since 594 was 524, checked in by lmdzadmin, 20 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 53.5 KB

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

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

Download in other formats:

Original Format