Context Navigation

isccp_cloud_types.F @ 486

Last change on this file since 486 was 486, checked in by lmdzadmin, 20 years ago
Phasage avec la version de Ionela IM/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.5 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 486 2003-12-15 17:50:41Z lmdzadmin $
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	c write(6,'(a10)') 'top_height='
284	c write(6,'(8I10)') top_height
285	c write(6,'(a10)') 'overlap='
286	c 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	call ran0_vec(npoints,seed,ran)
659
660	do j=1,npoints
661	threshold(j,ibox)=
662	!if max overlapped conv cloud
663	& maxocc(j,ibox) * (
664	& boxpos(j,ibox)
665	& ) +
666	!else
667	& (1-maxocc(j,ibox)) * (
668	!if max overlapped strat cloud
669	& (maxosc(j,ibox)) * (
670	!threshold=boxpos
671	& threshold(j,ibox)
672	& ) +
673	!else
674	& (1-maxosc(j,ibox)) * (
675	!threshold_min=random[thrmin,1]
676	& threshold_min(j,ibox)+
677	& (1-threshold_min(j,ibox))*ran(j)
678	& )
679	& )
680	enddo
681
682	ENDDO ! ibox
683
684	! Fill frac_out with 1's where tca is greater than the threshold
685
686	DO ibox=1,ncol
687	do j=1,npoints
688	if (tca(j,ilev).gt.threshold(j,ibox)) then
689	cIM REAL frac_out(j,ibox,ilev)=1
690	frac_out(j,ibox,ilev)=1.0
691	else
692	cIM REAL frac_out(j,ibox,ilev)=0
693	frac_out(j,ibox,ilev)=0.0
694	end if
695	enddo
696	ENDDO
697
698	! Code to partition boxes into startiform and convective parts
699	! goes here
700
701	DO ibox=1,ncol
702	do j=1,npoints
703	if (threshold(j,ibox).le.cca(j,ilev)) then
704	! = 2 IF threshold le cca(j)
705	cIM REAL frac_out(j,ibox,ilev) = 2
706	frac_out(j,ibox,ilev) = 2.0
707	else
708	! = the same IF NOT threshold le cca(j)
709	frac_out(j,ibox,ilev) = frac_out(j,ibox,ilev)
710	end if
711	enddo
712	ENDDO
713
714	! Set last_frac to tca at this level, so as to be tca
715	! from last level next time round
716
717	cIM
718	c if (ncolprint.ne.0) then
719
720	c do j=1,npoints ,1000
721	c write(6,'(a10)') 'j='
722	c write(6,'(8I10)') j
723	c write (6,'(a)') 'last_frac:'
724	c write (6,'(8f5.2)') (tca(j,ilev-1))
725
726	c write (6,'(a)') 'cca:'
727	c write (6,'(8f5.2)') (cca(j,ilev),ibox=1,ncolprint)
728
729	c write (6,'(a)') 'max_overlap_cc:'
730	c write (6,'(8f5.2)') (maxocc(j,ibox),ibox=1,ncolprint)
731
732	c write (6,'(a)') 'max_overlap_sc:'
733	c write (6,'(8f5.2)') (maxosc(j,ibox),ibox=1,ncolprint)
734
735	c write (6,'(a)') 'threshold_min_nsf2:'
736	c write (6,'(8f5.2)') (threshold_min(j,ibox),ibox=1,ncolprint)
737
738	c write (6,'(a)') 'threshold_nsf2:'
739	c write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)
740
741	c write (6,'(a)') 'frac_out_pp_rev:'
742	c write (6,'(8f5.2)')
743	c & ((frac_out(j,ibox,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
744	c enddo
745	c endif
746
747	200 CONTINUE !loop over nlev
748
749	!
750	! ---------------------------------------------------!
751
752
753	!
754	! ---------------------------------------------------!
755	! COMPUTE CLOUD OPTICAL DEPTH FOR EACH COLUMN and
756	! put into vector tau
757
758	!initialize tau and albedocld to zero
759	do 15 ibox=1,ncol
760	do j=1,npoints
761	tau(j,ibox)=0.
762	albedocld(j,ibox)=0.
763	boxtau(j,ibox)=0.
764	boxptop(j,ibox)=0.
765	box_cloudy(j,ibox)=.false.
766	enddo
767	15 continue
768
769	!compute total cloud optical depth for each column
770	do ilev=1,nlev
771	!increment tau for each of the boxes
772	do ibox=1,ncol
773	do j=1,npoints
774	cIM REAL if (frac_out(j,ibox,ilev).eq.1) then
775	if (frac_out(j,ibox,ilev).eq.1.0) then
776	tau(j,ibox)=tau(j,ibox)
777	& + dtau_s(j,ilev)
778	endif
779	cIM REAL if (frac_out(j,ibox,ilev).eq.2) then
780	if (frac_out(j,ibox,ilev).eq.2.0) then
781	tau(j,ibox)=tau(j,ibox)
782	& + dtau_c(j,ilev)
783	end if
784	enddo
785	enddo ! ibox
786	enddo ! ilev
787	cIM
788	c if (ncolprint.ne.0) then
789
790	c do j=1,npoints ,1000
791	c write(6,'(a10)') 'j='
792	c write(6,'(8I10)') j
793	c write(6,'(i2,1X,8(f7.2,1X))')
794	c & ilev,
795	c & (tau(j,ibox),ibox=1,ncolprint)
796	c enddo
797	c endif
798	!
799	! ---------------------------------------------------!
800
801
802
803	!
804	! ---------------------------------------------------!
805	! COMPUTE INFRARED BRIGHTNESS TEMPERUATRES
806	! AND CLOUD TOP TEMPERATURE SATELLITE SHOULD SEE
807	!
808	! again this is only done if top_height = 1 or 3
809	!
810	! fluxtop is the 10.5 micron radiance at the top of the
811	! atmosphere
812	! trans_layers_above is the total transmissivity in the layers
813	! above the current layer
814	! fluxtop_clrsky(j) and trans_layers_above_clrsky(j) are the clear
815	! sky versions of these quantities.
816
817	if (top_height .eq. 1 .or. top_height .eq. 3) then
818
819
820	!----------------------------------------------------------------------
821	!
822	! DO CLEAR SKY RADIANCE CALCULATION FIRST
823	!
824	!compute water vapor continuum emissivity
825	!this treatment follows Schwarkzopf and Ramasamy
826	!JGR 1999,vol 104, pages 9467-9499.
827	!the emissivity is calculated at a wavenumber of 955 cm-1,
828	!or 10.47 microns
829	wtmair = 28.9644
830	wtmh20 = 18.01534
831	Navo = 6.023E+23
832	grav = 9.806650E+02
833	pstd = 1.013250E+06
834	t0 = 296.
835	cIM
836	c if (ncolprint .ne. 0)
837	c & write(6,*) 'ilev pw (kg/m2) tauwv(j) dem_wv'
838	do 125 ilev=1,nlev
839	do j=1,npoints
840	!press and dpress are dyne/cm2 = Pascals *10
841	press(j) = pfull(j,ilev)*10.
842	dpress(j) = (phalf(j,ilev+1)-phalf(j,ilev))*10
843	!atmden = g/cm2 = kg/m2 / 10
844	atmden(j) = dpress(j)/grav
845	rvh20(j) = qv(j,ilev)*wtmair/wtmh20
846	wk(j) = rvh20(j)Navoatmden(j)/wtmair
847	rhoave(j) = (press(j)/pstd)*(t0/at(j,ilev))
848	rh20s(j) = rvh20(j)*rhoave(j)
849	rfrgn(j) = rhoave(j)-rh20s(j)
850	tmpexp(j) = exp(-0.02*(at(j,ilev)-t0))
851	tauwv(j) = wk(j)1.e-20(
852	& (0.0224697rh20s(j)tmpexp(j)) +
853	& (3.41817e-7rfrgn(j)) )0.98
854	dem_wv(j,ilev) = 1. - exp( -1. * tauwv(j))
855	enddo
856	cIM
857	c if (ncolprint .ne. 0) then
858	c do j=1,npoints ,1000
859	c write(6,'(a10)') 'j='
860	c write(6,'(8I10)') j
861	c write(6,'(i2,1X,3(f8.3,3X))') ilev,
862	c & qv(j,ilev)*(phalf(j,ilev+1)-phalf(j,ilev))/(grav/100.),
863	c & tauwv(j),dem_wv(j,ilev)
864	c enddo
865	c endif
866	125 continue
867
868	!initialize variables
869	do j=1,npoints
870	fluxtop_clrsky(j) = 0.
871	trans_layers_above_clrsky(j)=1.
872	enddo
873
874	do ilev=1,nlev
875	do j=1,npoints
876
877	! Black body emission at temperature of the layer
878
879	bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. )
880	!bb(j)= 5.67e-8at(j,ilev)*4
881
882	! increase TOA flux by flux emitted from layer
883	! times total transmittance in layers above
884
885	fluxtop_clrsky(j) = fluxtop_clrsky(j)
886	& + dem_wv(j,ilev)bb(j)trans_layers_above_clrsky(j)
887
888	! update trans_layers_above with transmissivity
889	! from this layer for next time around loop
890
891	trans_layers_above_clrsky(j)=
892	& trans_layers_above_clrsky(j)*(1.-dem_wv(j,ilev))
893
894
895	enddo
896	cIM
897	c if (ncolprint.ne.0) then
898	c do j=1,npoints ,1000
899	c write(6,'(a10)') 'j='
900	c write(6,'(8I10)') j
901	c write (6,'(a)') 'ilev:'
902	c write (6,'(I2)') ilev
903
904	c write (6,'(a)')
905	c & 'emiss_layer,100.bb(j),100.f,total_trans:'
906	c write (6,'(4(f7.2,1X))') dem_wv(j,ilev),100.*bb(j),
907	c & 100.*fluxtop_clrsky(j),trans_layers_above_clrsky(j)
908	c enddo
909	c endif
910
911	enddo !loop over level
912
913	do j=1,npoints
914	!add in surface emission
915	bb(j)=1/( exp(1307.27/skt(j)) - 1. )
916	!bb(j)=5.67e-8skt(j)*4
917
918	fluxtop_clrsky(j) = fluxtop_clrsky(j) + emsfc_lw * bb(j)
919	& * trans_layers_above_clrsky(j)
920	enddo
921
922	cIM
923	c if (ncolprint.ne.0) then
924	c do j=1,npoints ,1000
925	c write(6,'(a10)') 'j='
926	c write(6,'(8I10)') j
927	c write (6,'(a)') 'id:'
928	c write (6,'(a)') 'surface'
929
930	c write (6,'(a)') 'emsfc,100.bb(j),100.f,total_trans:'
931	c write (6,'(4(f7.2,1X))') emsfc_lw,100.*bb(j),
932	c & 100.*fluxtop_clrsky(j),
933	c & trans_layers_above_clrsky(j)
934	c enddo
935	c endif
936
937
938	!
939	! END OF CLEAR SKY CALCULATION
940	!
941	!----------------------------------------------------------------
942
943
944	cIM
945	c if (ncolprint.ne.0) then
946
947	c do j=1,npoints ,1000
948	c write(6,'(a10)') 'j='
949	c write(6,'(8I10)') j
950	c write (6,'(a)') 'ts:'
951	c write (6,'(8f7.2)') (skt(j),ibox=1,ncolprint)
952
953	c write (6,'(a)') 'ta_rev:'
954	c write (6,'(8f7.2)')
955	c & ((at(j,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
956
957	c enddo
958	c endif
959	!loop over columns
960	do ibox=1,ncol
961	do j=1,npoints
962	fluxtop(j,ibox)=0.
963	trans_layers_above(j,ibox)=1.
964	enddo
965	enddo
966
967	do ilev=1,nlev
968	do j=1,npoints
969	! Black body emission at temperature of the layer
970
971	bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. )
972	!bb(j)= 5.67e-8at(j,ilev)*4
973	enddo
974
975	do ibox=1,ncol
976	do j=1,npoints
977
978	! emissivity for point in this layer
979	cIM REAL if (frac_out(j,ibox,ilev).eq.1) then
980	if (frac_out(j,ibox,ilev).eq.1.0) then
981	dem(j,ibox)= 1. -
982	& ( (1. - dem_wv(j,ilev)) * (1. - dem_s(j,ilev)) )
983	cIM REAL else if (frac_out(j,ibox,ilev).eq.2) then
984	else if (frac_out(j,ibox,ilev).eq.2.0) then
985	dem(j,ibox)= 1. -
986	& ( (1. - dem_wv(j,ilev)) * (1. - dem_c(j,ilev)) )
987	else
988	dem(j,ibox)= dem_wv(j,ilev)
989	end if
990
991
992	! increase TOA flux by flux emitted from layer
993	! times total transmittance in layers above
994
995	fluxtop(j,ibox) = fluxtop(j,ibox)
996	& + dem(j,ibox) * bb(j)
997	& * trans_layers_above(j,ibox)
998
999	! update trans_layers_above with transmissivity
1000	! from this layer for next time around loop
1001
1002	trans_layers_above(j,ibox)=
1003	& trans_layers_above(j,ibox)*(1.-dem(j,ibox))
1004
1005	enddo ! j
1006	enddo ! ibox
1007
1008	cIM
1009	c if (ncolprint.ne.0) then
1010	c do j=1,npoints,1000
1011	c write (6,'(a)') 'ilev:'
1012	c write (6,'(I2)') ilev
1013
1014	c write(6,'(a10)') 'j='
1015	c write(6,'(8I10)') j
1016	c write (6,'(a)') 'emiss_layer:'
1017	c write (6,'(8f7.2)') (dem(j,ibox),ibox=1,ncolprint)
1018
1019	c write (6,'(a)') '100.*bb(j):'
1020	c write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
1021
1022	c write (6,'(a)') '100.*f:'
1023	c write (6,'(8f7.2)')
1024	c & (100.*fluxtop(j,ibox),ibox=1,ncolprint)
1025
1026	c write (6,'(a)') 'total_trans:'
1027	c write (6,'(8f7.2)')
1028	c & (trans_layers_above(j,ibox),ibox=1,ncolprint)
1029	c enddo
1030	c endif
1031
1032	enddo ! ilev
1033
1034
1035	do j=1,npoints
1036	!add in surface emission
1037	bb(j)=1/( exp(1307.27/skt(j)) - 1. )
1038	!bb(j)=5.67e-8skt(j)*4
1039	end do
1040
1041	do ibox=1,ncol
1042	do j=1,npoints
1043
1044	!add in surface emission
1045
1046	fluxtop(j,ibox) = fluxtop(j,ibox)
1047	& + emsfc_lw * bb(j)
1048	& * trans_layers_above(j,ibox)
1049
1050	end do
1051	end do
1052
1053	cIM
1054	c if (ncolprint.ne.0) then
1055
1056	c do j=1,npoints ,1000
1057	c write(6,'(a10)') 'j='
1058	c write(6,'(8I10)') j
1059	c write (6,'(a)') 'id:'
1060	c write (6,'(a)') 'surface'
1061
1062	c write (6,'(a)') 'emiss_layer:'
1063	c write (6,'(8f7.2)') (dem(1,ibox),ibox=1,ncolprint)
1064
1065	c write (6,'(a)') '100.*bb(j):'
1066	c write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
1067
1068	c write (6,'(a)') '100.*f:'
1069	c write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint)
1070	c end do
1071	c endif
1072
1073	!now that you have the top of atmosphere radiance account
1074	!for ISCCP procedures to determine cloud top temperature
1075
1076	!account for partially transmitting cloud recompute flux
1077	!ISCCP would see assuming a single layer cloud
1078	!note choice here of 2.13, as it is primarily ice
1079	!clouds which have partial emissivity and need the
1080	!adjustment performed in this section
1081	!
1082	!If it turns out that the cloud brightness temperature
1083	!is greater than 260K, then the liquid cloud conversion
1084	!factor of 2.56 is used.
1085	!
1086	!Note that this is discussed on pages 85-87 of
1087	!the ISCCP D level documentation (Rossow et al. 1996)
1088
1089	do j=1,npoints
1090	!compute minimum brightness temperature and optical depth
1091	btcmin(j) = 1. / ( exp(1307.27/(attrop(j)-5.)) - 1. )
1092	enddo
1093	do ibox=1,ncol
1094	do j=1,npoints
1095	transmax(j) = (fluxtop(j,ibox)-btcmin(j))
1096	& /(fluxtop_clrsky(j)-btcmin(j))
1097	!note that the initial setting of tauir(j) is needed so that
1098	!tauir(j) has a realistic value should the next if block be
1099	!bypassed
1100	tauir(j) = tau(j,ibox) * rec2p13
1101	taumin(j) = -1. * log(max(min(transmax(j),0.9999999),0.001))
1102
1103	enddo
1104
1105	if (top_height .eq. 1) then
1106	do j=1,npoints
1107	if (transmax(j) .gt. 0.001 .and.
1108	& transmax(j) .le. 0.9999999) then
1109	fluxtopinit(j) = fluxtop(j,ibox)
1110	tauir(j) = tau(j,ibox) *rec2p13
1111	endif
1112	enddo
1113	do icycle=1,2
1114	do j=1,npoints
1115	if (tau(j,ibox) .gt. (tauchk )) then
1116	if (transmax(j) .gt. 0.001 .and.
1117	& transmax(j) .le. 0.9999999) then
1118	emcld(j,ibox) = 1. - exp(-1. * tauir(j) )
1119	fluxtop(j,ibox) = fluxtopinit(j) -
1120	& ((1.-emcld(j,ibox))*fluxtop_clrsky(j))
1121	fluxtop(j,ibox)=max(1.E-06,
1122	& (fluxtop(j,ibox)/emcld(j,ibox)))
1123	tb(j,ibox)= 1307.27
1124	& / (log(1. + (1./fluxtop(j,ibox))))
1125	if (tb(j,ibox) .gt. 260.) then
1126	tauir(j) = tau(j,ibox) / 2.56
1127	end if
1128	end if
1129	end if
1130	enddo
1131	enddo
1132
1133	endif
1134
1135	do j=1,npoints
1136	if (tau(j,ibox) .gt. (tauchk )) then
1137	!cloudy box
1138	tb(j,ibox)= 1307.27/ (log(1. + (1./fluxtop(j,ibox))))
1139	if (top_height.eq.1.and.tauir(j).lt.taumin(j)) then
1140	tb(j,ibox) = attrop(j) - 5.
1141	tau(j,ibox) = 2.13*taumin(j)
1142	end if
1143	else
1144	!clear sky brightness temperature
1145	tb(j,ibox) = 1307.27/(log(1.+(1./fluxtop_clrsky(j))))
1146	end if
1147	enddo ! j
1148	enddo ! ibox
1149
1150	cIM
1151	c if (ncolprint.ne.0) then
1152
1153	c do j=1,npoints,1000
1154	c write(6,'(a10)') 'j='
1155	c write(6,'(8I10)') j
1156
1157	c write (6,'(a)') 'attrop:'
1158	c write (6,'(8f7.2)') (attrop(j))
1159
1160	c write (6,'(a)') 'btcmin:'
1161	c write (6,'(8f7.2)') (btcmin(j))
1162
1163	c write (6,'(a)') 'fluxtop_clrsky*100:'
1164	c write (6,'(8f7.2)')
1165	c & (100.*fluxtop_clrsky(j))
1166
1167	c write (6,'(a)') '100.*f_adj:'
1168	c write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint)
1169
1170	c write (6,'(a)') 'transmax:'
1171	c write (6,'(8f7.2)') (transmax(ibox),ibox=1,ncolprint)
1172
1173	c write (6,'(a)') 'tau:'
1174	c write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
1175
1176	c write (6,'(a)') 'emcld:'
1177	c write (6,'(8f7.2)') (emcld(j,ibox),ibox=1,ncolprint)
1178
1179	c write (6,'(a)') 'total_trans:'
1180	c write (6,'(8f7.2)')
1181	c & (trans_layers_above(j,ibox),ibox=1,ncolprint)
1182
1183	c write (6,'(a)') 'total_emiss:'
1184	c write (6,'(8f7.2)')
1185	c & (1.0-trans_layers_above(j,ibox),ibox=1,ncolprint)
1186
1187	c write (6,'(a)') 'total_trans:'
1188	c write (6,'(8f7.2)')
1189	c & (trans_layers_above(j,ibox),ibox=1,ncolprint)
1190
1191	c write (6,'(a)') 'ppout:'
1192	c write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
1193	c enddo ! j
1194	c endif
1195
1196	end if
1197
1198	! ---------------------------------------------------!
1199
1200	!
1201	! ---------------------------------------------------!
1202	! DETERMINE CLOUD TOP PRESSURE
1203	!
1204	! again the 2 methods differ according to whether
1205	! or not you use the physical cloud top pressure (top_height = 2)
1206	! or the radiatively determined cloud top pressure (top_height = 1 or 3)
1207	!
1208
1209	!compute cloud top pressure
1210	do 30 ibox=1,ncol
1211	!segregate according to optical thickness
1212	if (top_height .eq. 1 .or. top_height .eq. 3) then
1213	!find level whose temperature
1214	!most closely matches brightness temperature
1215	do j=1,npoints
1216	nmatch(j)=0
1217	enddo
1218	do 29 ilev=1,nlev-1
1219	!cdir nodep
1220	do j=1,npoints
1221	if ((at(j,ilev) .ge. tb(j,ibox) .and.
1222	& at(j,ilev+1) .lt. tb(j,ibox)) .or.
1223	& (at(j,ilev) .le. tb(j,ibox) .and.
1224	& at(j,ilev+1) .gt. tb(j,ibox))) then
1225
1226	nmatch(j)=nmatch(j)+1
1227	if(abs(at(j,ilev)-tb(j,ibox)) .lt.
1228	& abs(at(j,ilev+1)-tb(j,ibox))) then
1229	match(j,nmatch(j))=ilev
1230	else
1231	match(j,nmatch(j))=ilev+1
1232	end if
1233	end if
1234	enddo
1235	29 continue
1236
1237	do j=1,npoints
1238	if (nmatch(j) .ge. 1) then
1239	ptop(j,ibox)=pfull(j,match(j,nmatch(j)))
1240	levmatch(j,ibox)=match(j,nmatch(j))
1241	else
1242	if (tb(j,ibox) .lt. atmin(j)) then
1243	ptop(j,ibox)=ptrop(j)
1244	levmatch(j,ibox)=itrop(j)
1245	end if
1246	if (tb(j,ibox) .gt. atmax(j)) then
1247	ptop(j,ibox)=pfull(j,nlev)
1248	levmatch(j,ibox)=nlev
1249	end if
1250	end if
1251	enddo ! j
1252
1253	else ! if (top_height .eq. 1 .or. top_height .eq. 3)
1254
1255	do j=1,npoints
1256	ptop(j,ibox)=0.
1257	enddo
1258	do ilev=1,nlev
1259	do j=1,npoints
1260	if ((ptop(j,ibox) .eq. 0. )
1261	cIM & .and.(frac_out(j,ibox,ilev) .ne. 0)) then
1262	& .and.(frac_out(j,ibox,ilev) .ne. 0.0)) then
1263	ptop(j,ibox)=pfull(j,ilev)
1264	levmatch(j,ibox)=ilev
1265	end if
1266	end do
1267	end do
1268	end if
1269
1270	do j=1,npoints
1271	if (tau(j,ibox) .le. (tauchk )) then
1272	ptop(j,ibox)=0.
1273	levmatch(j,ibox)=0
1274	endif
1275	enddo
1276
1277	30 continue
1278
1279	!
1280	!
1281	! ---------------------------------------------------!
1282
1283
1284	!
1285	! ---------------------------------------------------!
1286	! DETERMINE ISCCP CLOUD TYPE FREQUENCIES
1287	!
1288	! Now that ptop and tau have been determined,
1289	! determine amount of each of the 49 ISCCP cloud
1290	! types
1291	!
1292	! Also compute grid box mean cloud top pressure and
1293	! optical thickness. The mean cloud top pressure and
1294	! optical thickness are averages over the cloudy
1295	! area only. The mean cloud top pressure is a linear
1296	! average of the cloud top pressures. The mean cloud
1297	! optical thickness is computed by converting optical
1298	! thickness to an albedo, averaging in albedo units,
1299	! then converting the average albedo back to a mean
1300	! optical thickness.
1301	!
1302
1303	!compute isccp frequencies
1304
1305	!reset frequencies
1306	do 38 ilev=1,7
1307	do 38 ilev2=1,7
1308	do j=1,npoints !
1309	fq_isccp(j,ilev,ilev2)=0.
1310	enddo
1311	38 continue
1312
1313	!reset variables need for averaging cloud properties
1314	do j=1,npoints
1315	totalcldarea(j) = 0.
1316	meanalbedocld(j) = 0.
1317	meanptop(j) = 0.
1318	meantaucld(j) = 0.
1319	enddo ! j
1320
1321	boxarea = 1./real(ncol)
1322
1323	!determine optical depth category
1324	cIM do 39 j=1,npoints
1325	cIM do ibox=1,ncol
1326	do 39 ibox=1,ncol
1327	do j=1,npoints
1328
1329	cIM
1330	c CALL CPU_time(t1)
1331	cIM
1332
1333	if (tau(j,ibox) .gt. (tauchk )
1334	& .and. ptop(j,ibox) .gt. 0.) then
1335	box_cloudy(j,ibox)=.true.
1336	endif
1337
1338	cIM
1339	c CALL CPU_time(t2)
1340	c print*,'IF tau t2 - t1',t2 - t1
1341
1342	c CALL CPU_time(t1)
1343	cIM
1344
1345	if (box_cloudy(j,ibox)) then
1346
1347	! totalcldarea always diagnosed day or night
1348	totalcldarea(j) = totalcldarea(j) + boxarea
1349
1350	if (sunlit(j).eq.1) then
1351
1352	! tau diagnostics only with sunlight
1353
1354	boxtau(j,ibox) = tau(j,ibox)
1355
1356	!convert optical thickness to albedo
1357	albedocld(j,ibox)
1358	& =real(invtau(min(nint(100.*tau(j,ibox)),45000)))
1359
1360	!contribute to averaging
1361	meanalbedocld(j) = meanalbedocld(j)
1362	& +albedocld(j,ibox)*boxarea
1363
1364	endif
1365
1366	endif
1367
1368	cIM
1369	c CALL CPU_time(t2)
1370	c print*,'IF box_cloudy t2 - t1',t2 - t1
1371
1372	c CALL CPU_time(t1)
1373	cIM BEG
1374	cIM !convert ptop to millibars
1375	ptop(j,ibox)=ptop(j,ibox) / 100.
1376
1377	cIM !save for output cloud top pressure and optical thickness
1378	boxptop(j,ibox) = ptop(j,ibox)
1379	cIM END
1380
1381	cIM BEG
1382	!reset itau(j), ipres(j)
1383	itau(j) = 0
1384	ipres(j) = 0
1385
1386	if (tau(j,ibox) .lt. isccp_taumin) then
1387	itau(j)=1
1388	else if (tau(j,ibox) .ge. isccp_taumin
1389	&
1390	& .and. tau(j,ibox) .lt. 1.3) then
1391	itau(j)=2
1392	else if (tau(j,ibox) .ge. 1.3
1393	& .and. tau(j,ibox) .lt. 3.6) then
1394	itau(j)=3
1395	else if (tau(j,ibox) .ge. 3.6
1396	& .and. tau(j,ibox) .lt. 9.4) then
1397	itau(j)=4
1398	else if (tau(j,ibox) .ge. 9.4
1399	& .and. tau(j,ibox) .lt. 23.) then
1400	itau(j)=5
1401	else if (tau(j,ibox) .ge. 23.
1402	& .and. tau(j,ibox) .lt. 60.) then
1403	itau(j)=6
1404	else if (tau(j,ibox) .ge. 60.) then
1405	itau(j)=7
1406	end if
1407
1408	!determine cloud top pressure category
1409	if ( ptop(j,ibox) .gt. 0.
1410	& .and.ptop(j,ibox) .lt. 180.) then
1411	ipres(j)=1
1412	else if(ptop(j,ibox) .ge. 180.
1413	& .and.ptop(j,ibox) .lt. 310.) then
1414	ipres(j)=2
1415	else if(ptop(j,ibox) .ge. 310.
1416	& .and.ptop(j,ibox) .lt. 440.) then
1417	ipres(j)=3
1418	else if(ptop(j,ibox) .ge. 440.
1419	& .and.ptop(j,ibox) .lt. 560.) then
1420	ipres(j)=4
1421	else if(ptop(j,ibox) .ge. 560.
1422	& .and.ptop(j,ibox) .lt. 680.) then
1423	ipres(j)=5
1424	else if(ptop(j,ibox) .ge. 680.
1425	& .and.ptop(j,ibox) .lt. 800.) then
1426	ipres(j)=6
1427	else if(ptop(j,ibox) .ge. 800.) then
1428	ipres(j)=7
1429	end if
1430	cIM END
1431
1432	if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1433
1434	cIM !convert ptop to millibars
1435	cIM ptop(j,ibox)=ptop(j,ibox) / 100.
1436
1437	cIM !save for output cloud top pressure and optical thickness
1438	cIM boxptop(j,ibox) = ptop(j,ibox)
1439
1440	if (box_cloudy(j,ibox)) then
1441
1442	meanptop(j) = meanptop(j) + ptop(j,ibox)*boxarea
1443
1444	cIM !reset itau(j), ipres(j)
1445	cIM itau(j) = 0
1446	cIM ipres(j) = 0
1447
1448	c if (tau(j,ibox) .lt. isccp_taumin) then
1449	c itau(j)=1
1450	c else if (tau(j,ibox) .ge. isccp_taumin
1451	c &
1452	c & .and. tau(j,ibox) .lt. 1.3) then
1453	c itau(j)=2
1454	c else if (tau(j,ibox) .ge. 1.3
1455	c & .and. tau(j,ibox) .lt. 3.6) then
1456	c itau(j)=3
1457	c else if (tau(j,ibox) .ge. 3.6
1458	c & .and. tau(j,ibox) .lt. 9.4) then
1459	c itau(j)=4
1460	c else if (tau(j,ibox) .ge. 9.4
1461	c & .and. tau(j,ibox) .lt. 23.) then
1462	c itau(j)=5
1463	c else if (tau(j,ibox) .ge. 23.
1464	c & .and. tau(j,ibox) .lt. 60.) then
1465	c itau(j)=6
1466	c else if (tau(j,ibox) .ge. 60.) then
1467	c itau(j)=7
1468	c end if
1469
1470	c !determine cloud top pressure category
1471	c if ( ptop(j,ibox) .gt. 0.
1472	c & .and.ptop(j,ibox) .lt. 180.) then
1473	c ipres(j)=1
1474	c else if(ptop(j,ibox) .ge. 180.
1475	c & .and.ptop(j,ibox) .lt. 310.) then
1476	c ipres(j)=2
1477	c else if(ptop(j,ibox) .ge. 310.
1478	c & .and.ptop(j,ibox) .lt. 440.) then
1479	c ipres(j)=3
1480	c else if(ptop(j,ibox) .ge. 440.
1481	c & .and.ptop(j,ibox) .lt. 560.) then
1482	c ipres(j)=4
1483	c else if(ptop(j,ibox) .ge. 560.
1484	c & .and.ptop(j,ibox) .lt. 680.) then
1485	c ipres(j)=5
1486	c else if(ptop(j,ibox) .ge. 680.
1487	c & .and.ptop(j,ibox) .lt. 800.) then
1488	c ipres(j)=6
1489	c else if(ptop(j,ibox) .ge. 800.) then
1490	c ipres(j)=7
1491	c end if
1492
1493	!update frequencies
1494	if(ipres(j) .gt. 0.and.itau(j) .gt. 0) then
1495	fq_isccp(j,itau(j),ipres(j))=
1496	& fq_isccp(j,itau(j),ipres(j))+ boxarea
1497	end if
1498
1499	cIM calcul stats regime dynamique BEG
1500	! iw(j) = int((w(j)-wmin)/pas_w) +1
1501	! pctj(itau(j),ipres(j),iw(j))=.FALSE.
1502	! !update frequencies W500
1503	! if (pct_ocean(j)) then
1504	! if (ipres(j) .gt. 0.and.itau(j) .gt. 0) then
1505	! if (iw(j) .gt. int(wmin).and.iw(j) .le. iwmx) then
1506	c print*,' ISCCP iw=',iw(j),j
1507	! fq_dynreg(itau(j),ipres(j),iw(j))=
1508	! & fq_dynreg(itau(j),ipres(j),iw(j))+
1509	! & boxarea
1510	c & fq_isccp(j,itau(j),ipres(j))
1511	! pctj(itau(j),ipres(j),iw(j))=.TRUE.
1512	c nfq_dynreg(itau(j),ipres(j),iw(j))=
1513	c & nfq_dynreg(itau(j),ipres(j),iw(j))+1.
1514	! end if
1515	! end if
1516	! end if
1517	cIM calcul stats regime dynamique END
1518	end if !IM boxcloudy
1519
1520	end if !IM sunlit
1521
1522	cIM
1523	c CALL CPU_time(t2)
1524	c print*,'IF sunlit boxcloudy t2 - t1',t2 - t1
1525	cIM
1526	enddo !IM ibox/j
1527
1528
1529	cIM ajout stats s/ W500 BEG
1530	cIM ajout stats s/ W500 END
1531
1532	c if(j.EQ.6722) then
1533	c print*,' ISCCP',w(j),iw(j),ipres(j),itau(j)
1534	c endif
1535
1536	! if (pct_ocean(j)) then
1537	! if (ipres(j) .gt. 0.and.itau(j) .gt. 0) then
1538	! if (iw(j) .gt. int(wmin).and.iw(j) .le. iwmx) then
1539	! if(pctj(itau(j),ipres(j),iw(j))) THEN
1540	! nfq_dynreg(itau(j),ipres(j),iw(j))=
1541	! & nfq_dynreg(itau(j),ipres(j),iw(j))+1.
1542	c if(itau(j).EQ.4.AND.ipres(j).EQ.2.AND.
1543	c & iw(j).EQ.10) then
1544	c PRINT*,' isccp AVANT',
1545	c & nfq_dynreg(itau(j),ipres(j),iw(j)),
1546	c & fq_dynreg(itau(j),ipres(j),iw(j))
1547	c endif
1548	! endif
1549	! endif
1550	! endif
1551	! endif
1552	39 continue !IM j/ibox
1553
1554	!compute mean cloud properties
1555	do j=1,npoints
1556	if (totalcldarea(j) .gt. 0.) then
1557	meanptop(j) = meanptop(j) / totalcldarea(j)
1558	if (sunlit(j).eq.1) then
1559	meanalbedocld(j) = meanalbedocld(j) / totalcldarea(j)
1560	meantaucld(j) = tautab(min(255,max(1,nint(meanalbedocld(j)))))
1561	end if
1562	end if
1563	enddo ! j
1564	!
1565	cIM ajout stats s/ W500 BEG
1566	! do nw = 1, iwmx
1567	! do l = 1, 7
1568	! do k = 1, 7
1569	! if (nfq_dynreg(k,l,nw).GT.0.) then
1570	! fq_dynreg(k,l,nw) = fq_dynreg(k,l,nw)/nfq_dynreg(k,l,nw)
1571	c if(k.EQ.4.AND.l.EQ.2.AND.nw.EQ.10) then
1572	c print*,' isccp APRES',nfq_dynreg(k,l,nw),
1573	c & fq_dynreg(k,l,nw)
1574	c endif
1575	! else
1576	! if(fq_dynreg(k,l,nw).NE.0.) then
1577	! print*,'nfq_dynreg = 0 tau,pc,nw',k,l,nw,fq_dynreg(k,l,nw)
1578	! endif
1579	c fq_dynreg(k,l,nw) = -1.E+20
1580	c nfq_dynreg(k,l,nw) = 1.E+20
1581	! end if
1582	! enddo !k
1583	! enddo !l
1584	! enddo !nw
1585	cIM ajout stats s/ W500 END
1586	! ---------------------------------------------------!
1587
1588	! ---------------------------------------------------!
1589	! OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM
1590	!
1591	!cIM
1592	c if (debugcol.ne.0) then
1593	!
1594	c do j=1,npoints,debugcol
1595
1596	c !produce character output
1597	c do ilev=1,nlev
1598	c do ibox=1,ncol
1599	c acc(ilev,ibox)=0
1600	c enddo
1601	c enddo
1602
1603	c do ilev=1,nlev
1604	c do ibox=1,ncol
1605	c acc(ilev,ibox)=frac_out(j,ibox,ilev)*2
1606	c if (levmatch(j,ibox) .eq. ilev)
1607	c & acc(ilev,ibox)=acc(ilev,ibox)+1
1608	c enddo
1609	c enddo
1610
1611	!print test
1612
1613	c write(ftn09,11) j
1614	c11 format('ftn09.',i4.4)
1615	c open(9, FILE=ftn09, FORM='FORMATTED')
1616
1617	c write(9,'(a1)') ' '
1618	c write(9,'(10i5)')
1619	c & (ilev,ilev=5,nlev,5)
1620	c write(9,'(a1)') ' '
1621
1622	c do ibox=1,ncol
1623	c write(9,'(40(a1),1x,40(a1))')
1624	c & (cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev)
1625	c & ,(cchar(acc(ilev,ibox)+1),ilev=1,nlev)
1626	c end do
1627	c close(9)
1628	c
1629	cIM
1630	c if (ncolprint.ne.0) then
1631	c write(6,'(a1)') ' '
1632	c write(6,'(a2,1X,5(a7,1X),a50)')
1633	c & 'ilev',
1634	c & 'pfull','at',
1635	c & 'cc*100','dem_s','dtau_s',
1636	c & 'cchar'
1637
1638	! do 4012 ilev=1,nlev
1639	! write(6,'(60i2)') (box(i,ilev),i=1,ncolprint)
1640	! write(6,'(i2,1X,5(f7.2,1X),50(a1))')
1641	! & ilev,
1642	! & pfull(j,ilev)/100.,at(j,ilev),
1643	! & cc(j,ilev)*100.0,dem_s(j,ilev),dtau_s(j,ilev)
1644	! & ,(cchar(acc(ilev,ibox)+1),ibox=1,ncolprint)
1645	!4012 continue
1646	c write (6,'(a)') 'skt(j):'
1647	c write (6,'(8f7.2)') skt(j)
1648
1649	c write (6,'(8I7)') (ibox,ibox=1,ncolprint)
1650
1651	c write (6,'(a)') 'tau:'
1652	c write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
1653
1654	c write (6,'(a)') 'tb:'
1655	c write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
1656
1657	c write (6,'(a)') 'ptop:'
1658	c write (6,'(8f7.2)') (ptop(j,ibox),ibox=1,ncolprint)
1659	c endif
1660
1661	c enddo
1662
1663	c end if
1664
1665	return
1666	end

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

Download in other formats:

Original Format