Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

isccp_cloud_types.F90 @ 3627

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

Converted to free source form files in libf/phylmd which were still in
fixed source form. The conversion was done using the polish mode of
the NAG Fortran Compiler.

In addition to converting to free source form, the processing of the
files also:

-- indented the code (including comments);

-- set Fortran keywords to uppercase, and set all other identifiers
to lower case;

-- added qualifiers to end statements (for example "end subroutine
conflx", instead of "end");

-- changed the terminating statements of all DO loops so that each
loop ends with an ENDDO statement (instead of a labeled continue).

-- replaced #include by include.

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: 45.5 KB

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

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

Download in other formats:

Original Format