Context Navigation

isccp_cloud_types.F90 @ 5103

Last change on this file since 5103 was 5103, checked in by abarral, 8 weeks ago
Handle CPP_INLANDSIS in lmdz_cppkeys_wrapper.F90 Remove obsolete key wrgrads_thermcell, _ADV_HALO, _ADV_HALLO, isminmax Remove redundant uses of CPPKEY_INCA (thanks acozic) Remove obsolete misc/write_field.F90 Remove unused ioipsl_* wrappers Remove calls to WriteField_u with wrong signature Convert .F -> .[fF]90 (lint) uppercase fortran operators [note: 1d and iso still broken - working on it]
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	! END OF CLEAR SKY CALCULATION
901
902	!----------------------------------------------------------------
903
904
905	! IM
906	! if (ncolprint.ne.0) then
907
908	! do j=1,npoints ,1000
909	! write(6,'(a10)') 'j='
910	! write(6,'(8I10)') j
911	! write (6,'(a)') 'ts:'
912	! write (6,'(8f7.2)') (skt(j),ibox=1,ncolprint)
913
914	! write (6,'(a)') 'ta_rev:'
915	! write (6,'(8f7.2)')
916	! & ((at(j,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
917
918	! enddo
919	! endif
920	!loop over columns
921	DO ibox = 1, ncol
922	DO j = 1, npoints
923	fluxtop(j, ibox) = 0.
924	trans_layers_above(j, ibox) = 1.
925	END DO
926	END DO
927
928	DO ilev = 1, nlev
929	DO j = 1, npoints
930	! Black body emission at temperature of the layer
931
932	bb(j) = 1/(exp(1307.27/at(j,ilev))-1.)
933	!bb(j)= 5.67e-8at(j,ilev)*4
934	END DO
935
936	DO ibox = 1, ncol
937	DO j = 1, npoints
938
939	! emissivity for point in this layer
940	! IM REAL if (frac_out(j,ibox,ilev).eq.1) then
941	IF (frac_out(j,ibox,ilev)==1.0) THEN
942	dem(j, ibox) = 1. - ((1.-dem_wv(j,ilev))*(1.-dem_s(j,ilev)))
943	! IM REAL else if (frac_out(j,ibox,ilev).eq.2) then
944	ELSE IF (frac_out(j,ibox,ilev)==2.0) THEN
945	dem(j, ibox) = 1. - ((1.-dem_wv(j,ilev))*(1.-dem_c(j,ilev)))
946	ELSE
947	dem(j, ibox) = dem_wv(j, ilev)
948	END IF
949
950
951	! increase TOA flux by flux emitted from layer
952	! times total transmittance in layers above
953
954	fluxtop(j, ibox) = fluxtop(j, ibox) + dem(j, ibox)bb(j) &
955	trans_layers_above(j, ibox)
956
957	! update trans_layers_above with transmissivity
958	! from this layer for next time around loop
959
960	trans_layers_above(j, ibox) = trans_layers_above(j, ibox)* &
961	(1.-dem(j,ibox))
962
963	END DO ! j
964	END DO ! ibox
965
966	! IM
967	! if (ncolprint.ne.0) then
968	! do j=1,npoints,1000
969	! write (6,'(a)') 'ilev:'
970	! write (6,'(I2)') ilev
971
972	! write(6,'(a10)') 'j='
973	! write(6,'(8I10)') j
974	! write (6,'(a)') 'emiss_layer:'
975	! write (6,'(8f7.2)') (dem(j,ibox),ibox=1,ncolprint)
976
977	! write (6,'(a)') '100.*bb(j):'
978	! write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
979
980	! write (6,'(a)') '100.*f:'
981	! write (6,'(8f7.2)')
982	! & (100.*fluxtop(j,ibox),ibox=1,ncolprint)
983
984	! write (6,'(a)') 'total_trans:'
985	! write (6,'(8f7.2)')
986	! & (trans_layers_above(j,ibox),ibox=1,ncolprint)
987	! enddo
988	! endif
989
990	END DO ! ilev
991
992
993	DO j = 1, npoints
994	!add in surface emission
995	bb(j) = 1/(exp(1307.27/skt(j))-1.)
996	!bb(j)=5.67e-8skt(j)*4
997	END DO
998
999	DO ibox = 1, ncol
1000	DO j = 1, npoints
1001
1002	!add in surface emission
1003
1004	fluxtop(j, ibox) = fluxtop(j, ibox) + emsfc_lwbb(j) &
1005	trans_layers_above(j, ibox)
1006
1007	END DO
1008	END DO
1009
1010	! IM
1011	! if (ncolprint.ne.0) then
1012
1013	! do j=1,npoints ,1000
1014	! write(6,'(a10)') 'j='
1015	! write(6,'(8I10)') j
1016	! write (6,'(a)') 'id:'
1017	! write (6,'(a)') 'surface'
1018
1019	! write (6,'(a)') 'emiss_layer:'
1020	! write (6,'(8f7.2)') (dem(1,ibox),ibox=1,ncolprint)
1021
1022	! write (6,'(a)') '100.*bb(j):'
1023	! write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
1024
1025	! write (6,'(a)') '100.*f:'
1026	! write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint)
1027	! END DO
1028	! endif
1029
1030	!now that you have the top of atmosphere radiance account
1031	!for ISCCP procedures to determine cloud top temperature
1032
1033	!account for partially transmitting cloud recompute flux
1034	!ISCCP would see assuming a single layer cloud
1035	!note choice here of 2.13, as it is primarily ice
1036	!clouds which have partial emissivity and need the
1037	!adjustment performed in this section
1038
1039	!If it turns out that the cloud brightness temperature
1040	!is greater than 260K, then the liquid cloud conversion
1041	!factor of 2.56 is used.
1042
1043	!Note that this is discussed on pages 85-87 of
1044	!the ISCCP D level documentation (Rossow et al. 1996)
1045
1046	DO j = 1, npoints
1047	!compute minimum brightness temperature and optical depth
1048	btcmin(j) = 1./(exp(1307.27/(attrop(j)-5.))-1.)
1049	END DO
1050	DO ibox = 1, ncol
1051	DO j = 1, npoints
1052	transmax(j) = (fluxtop(j,ibox)-btcmin(j))/(fluxtop_clrsky(j)-btcmin(j &
1053	))
1054	!note that the initial setting of tauir(j) is needed so that
1055	!tauir(j) has a realistic value should the next if block be
1056	!bypassed
1057	tauir(j) = tau(j, ibox)*rec2p13
1058	taumin(j) = -1.*log(max(min(transmax(j),0.9999999),0.001))
1059
1060	END DO
1061
1062	IF (top_height==1) THEN
1063	DO j = 1, npoints
1064	IF (transmax(j)>0.001 .AND. transmax(j)<=0.9999999) THEN
1065	fluxtopinit(j) = fluxtop(j, ibox)
1066	tauir(j) = tau(j, ibox)*rec2p13
1067	END IF
1068	END DO
1069	DO icycle = 1, 2
1070	DO j = 1, npoints
1071	IF (tau(j,ibox)>(tauchk)) THEN
1072	IF (transmax(j)>0.001 .AND. transmax(j)<=0.9999999) THEN
1073	emcld(j, ibox) = 1. - exp(-1.*tauir(j))
1074	fluxtop(j, ibox) = fluxtopinit(j) - ((1.-emcld(j, &
1075	ibox))*fluxtop_clrsky(j))
1076	fluxtop(j, ibox) = max(1.E-06, (fluxtop(j,ibox)/emcld(j, &
1077	ibox)))
1078	tb(j, ibox) = 1307.27/(log(1.+(1./fluxtop(j,ibox))))
1079	IF (tb(j,ibox)>260.) THEN
1080	tauir(j) = tau(j, ibox)/2.56
1081	END IF
1082	END IF
1083	END IF
1084	END DO
1085	END DO
1086
1087	END IF
1088
1089	DO j = 1, npoints
1090	IF (tau(j,ibox)>(tauchk)) THEN
1091	!cloudy box
1092	tb(j, ibox) = 1307.27/(log(1.+(1./fluxtop(j,ibox))))
1093	IF (top_height==1 .AND. tauir(j)<taumin(j)) THEN
1094	tb(j, ibox) = attrop(j) - 5.
1095	tau(j, ibox) = 2.13*taumin(j)
1096	END IF
1097	ELSE
1098	!clear sky brightness temperature
1099	tb(j, ibox) = 1307.27/(log(1.+(1./fluxtop_clrsky(j))))
1100	END IF
1101	END DO ! j
1102	END DO ! ibox
1103
1104	! IM
1105	! if (ncolprint.ne.0) then
1106
1107	! do j=1,npoints,1000
1108	! write(6,'(a10)') 'j='
1109	! write(6,'(8I10)') j
1110
1111	! write (6,'(a)') 'attrop:'
1112	! write (6,'(8f7.2)') (attrop(j))
1113
1114	! write (6,'(a)') 'btcmin:'
1115	! write (6,'(8f7.2)') (btcmin(j))
1116
1117	! write (6,'(a)') 'fluxtop_clrsky*100:'
1118	! write (6,'(8f7.2)')
1119	! & (100.*fluxtop_clrsky(j))
1120
1121	! write (6,'(a)') '100.*f_adj:'
1122	! write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint)
1123
1124	! write (6,'(a)') 'transmax:'
1125	! write (6,'(8f7.2)') (transmax(ibox),ibox=1,ncolprint)
1126
1127	! write (6,'(a)') 'tau:'
1128	! write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
1129
1130	! write (6,'(a)') 'emcld:'
1131	! write (6,'(8f7.2)') (emcld(j,ibox),ibox=1,ncolprint)
1132
1133	! write (6,'(a)') 'total_trans:'
1134	! write (6,'(8f7.2)')
1135	! & (trans_layers_above(j,ibox),ibox=1,ncolprint)
1136
1137	! write (6,'(a)') 'total_emiss:'
1138	! write (6,'(8f7.2)')
1139	! & (1.0-trans_layers_above(j,ibox),ibox=1,ncolprint)
1140
1141	! write (6,'(a)') 'total_trans:'
1142	! write (6,'(8f7.2)')
1143	! & (trans_layers_above(j,ibox),ibox=1,ncolprint)
1144
1145	! write (6,'(a)') 'ppout:'
1146	! write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
1147	! enddo ! j
1148	! endif
1149
1150	END IF
1151
1152	! ---------------------------------------------------!
1153
1154
1155	! ---------------------------------------------------!
1156	! DETERMINE CLOUD TOP PRESSURE
1157
1158	! again the 2 methods differ according to whether
1159	! or not you use the physical cloud top pressure (top_height = 2)
1160	! or the radiatively determined cloud top pressure (top_height = 1 or 3)
1161
1162
1163	!compute cloud top pressure
1164	DO ibox = 1, ncol
1165	!segregate according to optical thickness
1166	IF (top_height==1 .OR. top_height==3) THEN
1167	!find level whose temperature
1168	!most closely matches brightness temperature
1169	DO j = 1, npoints
1170	nmatch(j) = 0
1171	END DO
1172	DO ilev = 1, nlev - 1
1173	! cdir nodep
1174	DO j = 1, npoints
1175	IF ((at(j,ilev)>=tb(j,ibox) .AND. at(j,ilev+1)<tb(j, &
1176	ibox)) .OR. (at(j,ilev)<=tb(j,ibox) .AND. at(j,ilev+1)>tb(j, &
1177	ibox))) THEN
1178
1179	nmatch(j) = nmatch(j) + 1
1180	IF (abs(at(j,ilev)-tb(j,ibox))<abs(at(j,ilev+1)-tb(j,ibox))) THEN
1181	match(j, nmatch(j)) = ilev
1182	ELSE
1183	match(j, nmatch(j)) = ilev + 1
1184	END IF
1185	END IF
1186	END DO
1187	END DO
1188
1189	DO j = 1, npoints
1190	IF (nmatch(j)>=1) THEN
1191	ptop(j, ibox) = pfull(j, match(j,nmatch(j)))
1192	levmatch(j, ibox) = match(j, nmatch(j))
1193	ELSE
1194	IF (tb(j,ibox)<atmin(j)) THEN
1195	ptop(j, ibox) = ptrop(j)
1196	levmatch(j, ibox) = itrop(j)
1197	END IF
1198	IF (tb(j,ibox)>atmax(j)) THEN
1199	ptop(j, ibox) = pfull(j, nlev)
1200	levmatch(j, ibox) = nlev
1201	END IF
1202	END IF
1203	END DO ! j
1204
1205	ELSE ! if (top_height .eq. 1 .or. top_height .eq. 3)
1206
1207	DO j = 1, npoints
1208	ptop(j, ibox) = 0.
1209	END DO
1210	DO ilev = 1, nlev
1211	DO j = 1, npoints
1212	IF ((ptop(j,ibox)==0.) & ! IM &
1213	! .and.(frac_out(j,ibox,ilev) .ne. 0))
1214	! then
1215	.AND. (frac_out(j,ibox,ilev)/=0.0)) THEN
1216	ptop(j, ibox) = pfull(j, ilev)
1217	levmatch(j, ibox) = ilev
1218	END IF
1219	END DO
1220	END DO
1221	END IF
1222
1223	DO j = 1, npoints
1224	IF (tau(j,ibox)<=(tauchk)) THEN
1225	ptop(j, ibox) = 0.
1226	levmatch(j, ibox) = 0
1227	END IF
1228	END DO
1229
1230	END DO
1231
1232
1233
1234	! ---------------------------------------------------!
1235
1236
1237
1238	! ---------------------------------------------------!
1239	! DETERMINE ISCCP CLOUD TYPE FREQUENCIES
1240
1241	! Now that ptop and tau have been determined,
1242	! determine amount of each of the 49 ISCCP cloud
1243	! types
1244
1245	! Also compute grid box mean cloud top pressure and
1246	! optical thickness. The mean cloud top pressure and
1247	! optical thickness are averages over the cloudy
1248	! area only. The mean cloud top pressure is a linear
1249	! average of the cloud top pressures. The mean cloud
1250	! optical thickness is computed by converting optical
1251	! thickness to an albedo, averaging in albedo units,
1252	! then converting the average albedo back to a mean
1253	! optical thickness.
1254
1255
1256	!compute isccp frequencies
1257
1258	!reset frequencies
1259	DO ilev = 1, 7
1260	DO ilev2 = 1, 7
1261	DO j = 1, npoints !
1262	fq_isccp(j, ilev, ilev2) = 0.
1263	END DO
1264	END DO
1265	END DO
1266
1267	!reset variables need for averaging cloud properties
1268	DO j = 1, npoints
1269	totalcldarea(j) = 0.
1270	meanalbedocld(j) = 0.
1271	meanptop(j) = 0.
1272	meantaucld(j) = 0.
1273	END DO ! j
1274
1275	boxarea = 1./real(ncol)
1276
1277	!determine optical depth category
1278	! IM do 39 j=1,npoints
1279	! IM do ibox=1,ncol
1280	DO ibox = 1, ncol
1281	DO j = 1, npoints
1282
1283	! IM
1284	! CALL CPU_time(t1)
1285	! IM
1286
1287	IF (tau(j,ibox)>(tauchk) .AND. ptop(j,ibox)>0.) THEN
1288	box_cloudy(j, ibox) = .TRUE.
1289	END IF
1290
1291	! IM
1292	! CALL CPU_time(t2)
1293	! PRINT*,'IF tau t2 - t1',t2 - t1
1294
1295	! CALL CPU_time(t1)
1296	! IM
1297
1298	IF (box_cloudy(j,ibox)) THEN
1299
1300	! totalcldarea always diagnosed day or night
1301	totalcldarea(j) = totalcldarea(j) + boxarea
1302
1303	IF (sunlit(j)==1) THEN
1304
1305	! tau diagnostics only with sunlight
1306
1307	boxtau(j, ibox) = tau(j, ibox)
1308
1309	!convert optical thickness to albedo
1310	albedocld(j, ibox) = real(invtau(min(nint(100.*tau(j,ibox)), &
1311	45000)))
1312
1313	!contribute to averaging
1314	meanalbedocld(j) = meanalbedocld(j) + albedocld(j, ibox)*boxarea
1315
1316	END IF
1317
1318	END IF
1319
1320	! IM
1321	! CALL CPU_time(t2)
1322	! PRINT*,'IF box_cloudy t2 - t1',t2 - t1
1323
1324	! CALL CPU_time(t1)
1325	! IM BEG
1326	! IM !convert ptop to millibars
1327	ptop(j, ibox) = ptop(j, ibox)/100.
1328
1329	! IM !save for output cloud top pressure and optical
1330	! thickness
1331	boxptop(j, ibox) = ptop(j, ibox)
1332	! IM END
1333
1334	! IM BEG
1335	!reset itau(j), ipres(j)
1336	itau(j) = 0
1337	ipres(j) = 0
1338
1339	IF (tau(j,ibox)<isccp_taumin) THEN
1340	itau(j) = 1
1341	ELSE IF (tau(j,ibox)>=isccp_taumin .AND. tau(j,ibox)<1.3) THEN
1342	itau(j) = 2
1343	ELSE IF (tau(j,ibox)>=1.3 .AND. tau(j,ibox)<3.6) THEN
1344	itau(j) = 3
1345	ELSE IF (tau(j,ibox)>=3.6 .AND. tau(j,ibox)<9.4) THEN
1346	itau(j) = 4
1347	ELSE IF (tau(j,ibox)>=9.4 .AND. tau(j,ibox)<23.) THEN
1348	itau(j) = 5
1349	ELSE IF (tau(j,ibox)>=23. .AND. tau(j,ibox)<60.) THEN
1350	itau(j) = 6
1351	ELSE IF (tau(j,ibox)>=60.) THEN
1352	itau(j) = 7
1353	END IF
1354
1355	!determine cloud top pressure category
1356	IF (ptop(j,ibox)>0. .AND. ptop(j,ibox)<180.) THEN
1357	ipres(j) = 1
1358	ELSE IF (ptop(j,ibox)>=180. .AND. ptop(j,ibox)<310.) THEN
1359	ipres(j) = 2
1360	ELSE IF (ptop(j,ibox)>=310. .AND. ptop(j,ibox)<440.) THEN
1361	ipres(j) = 3
1362	ELSE IF (ptop(j,ibox)>=440. .AND. ptop(j,ibox)<560.) THEN
1363	ipres(j) = 4
1364	ELSE IF (ptop(j,ibox)>=560. .AND. ptop(j,ibox)<680.) THEN
1365	ipres(j) = 5
1366	ELSE IF (ptop(j,ibox)>=680. .AND. ptop(j,ibox)<800.) THEN
1367	ipres(j) = 6
1368	ELSE IF (ptop(j,ibox)>=800.) THEN
1369	ipres(j) = 7
1370	END IF
1371	! IM END
1372
1373	IF (sunlit(j)==1 .OR. top_height==3) THEN
1374
1375	! IM !convert ptop to millibars
1376	! IM ptop(j,ibox)=ptop(j,ibox) / 100.
1377
1378	! IM !save for output cloud top pressure and optical
1379	! thickness
1380	! IM boxptop(j,ibox) = ptop(j,ibox)
1381
1382	IF (box_cloudy(j,ibox)) THEN
1383
1384	meanptop(j) = meanptop(j) + ptop(j, ibox)*boxarea
1385
1386	! IM !reset itau(j), ipres(j)
1387	! IM itau(j) = 0
1388	! IM ipres(j) = 0
1389
1390	! if (tau(j,ibox) .lt. isccp_taumin) then
1391	! itau(j)=1
1392	! else if (tau(j,ibox) .ge. isccp_taumin
1393	! &
1394	! & .and. tau(j,ibox) .lt. 1.3) then
1395	! itau(j)=2
1396	! else if (tau(j,ibox) .ge. 1.3
1397	! & .and. tau(j,ibox) .lt. 3.6) then
1398	! itau(j)=3
1399	! else if (tau(j,ibox) .ge. 3.6
1400	! & .and. tau(j,ibox) .lt. 9.4) then
1401	! itau(j)=4
1402	! else if (tau(j,ibox) .ge. 9.4
1403	! & .and. tau(j,ibox) .lt. 23.) then
1404	! itau(j)=5
1405	! else if (tau(j,ibox) .ge. 23.
1406	! & .and. tau(j,ibox) .lt. 60.) then
1407	! itau(j)=6
1408	! else if (tau(j,ibox) .ge. 60.) then
1409	! itau(j)=7
1410	! end if
1411
1412	! !determine cloud top pressure category
1413	! if ( ptop(j,ibox) .gt. 0.
1414	! & .and.ptop(j,ibox) .lt. 180.) then
1415	! ipres(j)=1
1416	! else if(ptop(j,ibox) .ge. 180.
1417	! & .and.ptop(j,ibox) .lt. 310.) then
1418	! ipres(j)=2
1419	! else if(ptop(j,ibox) .ge. 310.
1420	! & .and.ptop(j,ibox) .lt. 440.) then
1421	! ipres(j)=3
1422	! else if(ptop(j,ibox) .ge. 440.
1423	! & .and.ptop(j,ibox) .lt. 560.) then
1424	! ipres(j)=4
1425	! else if(ptop(j,ibox) .ge. 560.
1426	! & .and.ptop(j,ibox) .lt. 680.) then
1427	! ipres(j)=5
1428	! else if(ptop(j,ibox) .ge. 680.
1429	! & .and.ptop(j,ibox) .lt. 800.) then
1430	! ipres(j)=6
1431	! else if(ptop(j,ibox) .ge. 800.) then
1432	! ipres(j)=7
1433	! end if
1434
1435	!update frequencies
1436	IF (ipres(j)>0 .AND. itau(j)>0) THEN
1437	fq_isccp(j, itau(j), ipres(j)) = fq_isccp(j, itau(j), ipres(j)) + &
1438	boxarea
1439	END IF
1440
1441	! IM calcul stats regime dynamique BEG
1442	! iw(j) = int((w(j)-wmin)/pas_w) +1
1443	! pctj(itau(j),ipres(j),iw(j))=.FALSE.
1444	! !update frequencies W500
1445	! if (pct_ocean(j)) then
1446	! if (ipres(j) .gt. 0.and.itau(j) .gt. 0) then
1447	! if (iw(j) .gt. int(wmin).and.iw(j) .le. iwmx) then
1448	! PRINT*,' ISCCP iw=',iw(j),j
1449	! fq_dynreg(itau(j),ipres(j),iw(j))=
1450	! & fq_dynreg(itau(j),ipres(j),iw(j))+
1451	! & boxarea
1452	! & fq_isccp(j,itau(j),ipres(j))
1453	! pctj(itau(j),ipres(j),iw(j))=.TRUE.
1454	! nfq_dynreg(itau(j),ipres(j),iw(j))=
1455	! & nfq_dynreg(itau(j),ipres(j),iw(j))+1.
1456	! end if
1457	! end if
1458	! end if
1459	! IM calcul stats regime dynamique END
1460	END IF !IM boxcloudy
1461
1462	END IF !IM sunlit
1463
1464	! IM
1465	! CALL CPU_time(t2)
1466	! PRINT*,'IF sunlit boxcloudy t2 - t1',t2 - t1
1467	! IM
1468	END DO !IM ibox/j
1469
1470
1471	! IM ajout stats s/ W500 BEG
1472	! IM ajout stats s/ W500 END
1473
1474	! if(j.EQ.6722) then
1475	! PRINT*,' ISCCP',w(j),iw(j),ipres(j),itau(j)
1476	! endif
1477
1478	! if (pct_ocean(j)) then
1479	! if (ipres(j) .gt. 0.and.itau(j) .gt. 0) then
1480	! if (iw(j) .gt. int(wmin).and.iw(j) .le. iwmx) then
1481	! if(pctj(itau(j),ipres(j),iw(j))) THEN
1482	! nfq_dynreg(itau(j),ipres(j),iw(j))=
1483	! & nfq_dynreg(itau(j),ipres(j),iw(j))+1.
1484	! if(itau(j).EQ.4.AND.ipres(j).EQ.2.AND.
1485	! & iw(j).EQ.10) then
1486	! PRINT*,' isccp AVANT',
1487	! & nfq_dynreg(itau(j),ipres(j),iw(j)),
1488	! & fq_dynreg(itau(j),ipres(j),iw(j))
1489	! endif
1490	! endif
1491	! endif
1492	! endif
1493	! endif
1494
1495	END DO
1496	!compute mean cloud properties
1497	! IM j/ibox
1498	DO j = 1, npoints
1499	IF (totalcldarea(j)>0.) THEN
1500	meanptop(j) = meanptop(j)/totalcldarea(j)
1501	IF (sunlit(j)==1) THEN
1502	meanalbedocld(j) = meanalbedocld(j)/totalcldarea(j)
1503	meantaucld(j) = tautab(min(255,max(1,nint(meanalbedocld(j)))))
1504	END IF
1505	END IF
1506	END DO ! j
1507
1508	! IM ajout stats s/ W500 BEG
1509	! do nw = 1, iwmx
1510	! do l = 1, 7
1511	! do k = 1, 7
1512	! if (nfq_dynreg(k,l,nw).GT.0.) then
1513	! fq_dynreg(k,l,nw) = fq_dynreg(k,l,nw)/nfq_dynreg(k,l,nw)
1514	! if(k.EQ.4.AND.l.EQ.2.AND.nw.EQ.10) then
1515	! PRINT*,' isccp APRES',nfq_dynreg(k,l,nw),
1516	! & fq_dynreg(k,l,nw)
1517	! endif
1518	! else
1519	! if(fq_dynreg(k,l,nw).NE.0.) then
1520	! PRINT*,'nfq_dynreg = 0 tau,pc,nw',k,l,nw,fq_dynreg(k,l,nw)
1521	! endif
1522	! fq_dynreg(k,l,nw) = -1.E+20
1523	! nfq_dynreg(k,l,nw) = 1.E+20
1524	! end if
1525	! enddo !k
1526	! enddo !l
1527	! enddo !nw
1528	! IM ajout stats s/ W500 END
1529	! ---------------------------------------------------!
1530
1531	! ---------------------------------------------------!
1532	! OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM
1533
1534	! cIM
1535	! if (debugcol.ne.0) then
1536
1537	! do j=1,npoints,debugcol
1538
1539	! !produce character output
1540	! do ilev=1,nlev
1541	! do ibox=1,ncol
1542	! acc(ilev,ibox)=0
1543	! enddo
1544	! enddo
1545
1546	! do ilev=1,nlev
1547	! do ibox=1,ncol
1548	! acc(ilev,ibox)=frac_out(j,ibox,ilev)*2
1549	! if (levmatch(j,ibox) .eq. ilev)
1550	! & acc(ilev,ibox)=acc(ilev,ibox)+1
1551	! enddo
1552	! enddo
1553
1554	!print test
1555
1556	! write(ftn09,11) j
1557	! 11 format('ftn09.',i4.4)
1558	! open(9, FILE=ftn09, FORM='FORMATTED')
1559
1560	! write(9,'(a1)') ' '
1561	! write(9,'(10i5)')
1562	! & (ilev,ilev=5,nlev,5)
1563	! write(9,'(a1)') ' '
1564
1565	! do ibox=1,ncol
1566	! write(9,'(40(a1),1x,40(a1))')
1567	! & (cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev)
1568	! & ,(cchar(acc(ilev,ibox)+1),ilev=1,nlev)
1569	! END DO
1570	! close(9)
1571
1572	! IM
1573	! if (ncolprint.ne.0) then
1574	! write(6,'(a1)') ' '
1575	! write(6,'(a2,1X,5(a7,1X),a50)')
1576	! & 'ilev',
1577	! & 'pfull','at',
1578	! & 'cc*100','dem_s','dtau_s',
1579	! & 'cchar'
1580
1581	! do 4012 ilev=1,nlev
1582	! write(6,'(60i2)') (box(i,ilev),i=1,ncolprint)
1583	! write(6,'(i2,1X,5(f7.2,1X),50(a1))')
1584	! & ilev,
1585	! & pfull(j,ilev)/100.,at(j,ilev),
1586	! & cc(j,ilev)*100.0,dem_s(j,ilev),dtau_s(j,ilev)
1587	! & ,(cchar(acc(ilev,ibox)+1),ibox=1,ncolprint)
1588	! 4012 continue
1589	! write (6,'(a)') 'skt(j):'
1590	! write (6,'(8f7.2)') skt(j)
1591
1592	! write (6,'(8I7)') (ibox,ibox=1,ncolprint)
1593
1594	! write (6,'(a)') 'tau:'
1595	! write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
1596
1597	! write (6,'(a)') 'tb:'
1598	! write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
1599
1600	! write (6,'(a)') 'ptop:'
1601	! write (6,'(8f7.2)') (ptop(j,ibox),ibox=1,ncolprint)
1602	! endif
1603
1604	! enddo
1605
1606	! end if
1607
1608	RETURN
1609	END SUBROUTINE isccp_cloud_types

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

Download in other formats:

Original Format