Context Navigation

isccp_cloud_types.F @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 55.3 KB

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

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

Download in other formats:

Original Format