Context Navigation

icarus.F @ 4100

Last change on this file since 4100 was 1454, checked in by Laurent Fairhead, 14 years ago

Merge of LMDZ5V1.0-dev branch r1453 into LMDZ5 trunk r1434

Fusion entre la version r1453 de la branche de développement LMDZ5V1.0-dev
et le tronc LMDZ5 (r1434)

File size: 42.8 KB

Line
1	SUBROUTINE ICARUS(
2	& debug,
3	& debugcol,
4	& npoints,
5	& sunlit,
6	& nlev,
7	& ncol,
8	& pfull,
9	& phalf,
10	& qv,
11	& cc,
12	& conv,
13	& dtau_s,
14	& dtau_c,
15	& top_height,
16	& top_height_direction,
17	& overlap,
18	& frac_out,
19	& skt,
20	& emsfc_lw,
21	& at,
22	& dem_s,
23	& dem_c,
24	& fq_isccp,
25	& totalcldarea,
26	& meanptop,
27	& meantaucld,
28	& meanalbedocld,
29	& meantb,
30	& meantbclr,
31	& boxtau,
32	& boxptop
33	&)
34
35	!$Id: icarus.f,v 4.1 2010/05/27 16:30:18 hadmw Exp $
36
37	! ***************************COPYRIGHT**************************
38	! (c) 2009, Lawrence Livermore National Security Limited Liability
39	! Corporation.
40	! All rights reserved.
41	!
42	! Redistribution and use in source and binary forms, with or without
43	! modification, are permitted provided that the
44	! following conditions are met:
45	!
46	! * Redistributions of source code must retain the above
47	! copyright notice, this list of conditions and the following
48	! disclaimer.
49	! * Redistributions in binary form must reproduce the above
50	! copyright notice, this list of conditions and the following
51	! disclaimer in the documentation and/or other materials
52	! provided with the distribution.
53	! * Neither the name of the Lawrence Livermore National Security
54	! Limited Liability Corporation nor the names of its
55	! contributors may be used to endorse or promote products
56	! derived from this software without specific prior written
57	! permission.
58	!
59	! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
60	! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
61	! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
62	! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
63	! OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
64	! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
65	! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
66	! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
67	! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
68	! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
69	! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
70	!
71	! ***************************COPYRIGHT*****************************
72	! ***************************COPYRIGHT*****************************
73	! ***************************COPYRIGHT*****************************
74	! ***************************COPYRIGHT*****************************
75
76	implicit none
77
78	! NOTE: the maximum number of levels and columns is set by
79	! the following parameter statement
80
81	INTEGER ncolprint
82
83	! -----
84	! Input
85	! -----
86
87	INTEGER npoints ! number of model points in the horizontal
88	INTEGER nlev ! number of model levels in column
89	INTEGER ncol ! number of subcolumns
90
91	INTEGER sunlit(npoints) ! 1 for day points, 0 for night time
92
93	REAL pfull(npoints,nlev)
94	! pressure of full model levels (Pascals)
95	! pfull(npoints,1) is top level of model
96	! pfull(npoints,nlev) is bot of model
97
98	REAL phalf(npoints,nlev+1)
99	! pressure of half model levels (Pascals)
100	! phalf(npoints,1) is top of model
101	! phalf(npoints,nlev+1) is the surface pressure
102
103	REAL qv(npoints,nlev)
104	! water vapor specific humidity (kg vapor/ kg air)
105	! on full model levels
106
107	REAL cc(npoints,nlev)
108	! input cloud cover in each model level (fraction)
109	! NOTE: This is the HORIZONTAL area of each
110	! grid box covered by clouds
111
112	REAL conv(npoints,nlev)
113	! input convective cloud cover in each model
114	! level (fraction)
115	! NOTE: This is the HORIZONTAL area of each
116	! grid box covered by convective clouds
117
118	REAL dtau_s(npoints,nlev)
119	! mean 0.67 micron optical depth of stratiform
120	! clouds in each model level
121	! NOTE: this the cloud optical depth of only the
122	! cloudy part of the grid box, it is not weighted
123	! with the 0 cloud optical depth of the clear
124	! part of the grid box
125
126	REAL dtau_c(npoints,nlev)
127	! mean 0.67 micron optical depth of convective
128	! clouds in each
129	! model level. Same note applies as in dtau_s.
130
131	INTEGER overlap ! overlap type
132	! 1=max
133	! 2=rand
134	! 3=max/rand
135
136	INTEGER top_height ! 1 = adjust top height using both a computed
137	! infrared brightness temperature and the visible
138	! optical depth to adjust cloud top pressure. Note
139	! that this calculation is most appropriate to compare
140	! to ISCCP data during sunlit hours.
141	! 2 = do not adjust top height, that is cloud top
142	! pressure is the actual cloud top pressure
143	! in the model
144	! 3 = adjust top height using only the computed
145	! infrared brightness temperature. Note that this
146	! calculation is most appropriate to compare to ISCCP
147	! IR only algortihm (i.e. you can compare to nighttime
148	! ISCCP data with this option)
149
150	INTEGER top_height_direction ! direction for finding atmosphere pressure level
151	! with interpolated temperature equal to the radiance
152	! determined cloud-top temperature
153	!
154	! 1 = find the lowest altitude (highest pressure) level
155	! with interpolated temperature equal to the radiance
156	! determined cloud-top temperature
157	!
158	! 2 = find the highest altitude (lowest pressure) level
159	! with interpolated temperature equal to the radiance
160	! determined cloud-top temperature
161	!
162	! ONLY APPLICABLE IF top_height EQUALS 1 or 3
163	! !
164	! 1 = old setting: matches all versions of
165	! ISCCP simulator with versions numbers 3.5.1 and lower
166	!
167	! 2 = default setting: for version numbers 4.0 and higher
168	!
169	! The following input variables are used only if top_height = 1 or top_height = 3
170	!
171	REAL skt(npoints) ! skin Temperature (K)
172	REAL emsfc_lw ! 10.5 micron emissivity of surface (fraction)
173	REAL at(npoints,nlev) ! temperature in each model level (K)
174	REAL dem_s(npoints,nlev) ! 10.5 micron longwave emissivity of stratiform
175	! clouds in each
176	! model level. Same note applies as in dtau_s.
177	REAL dem_c(npoints,nlev) ! 10.5 micron longwave emissivity of convective
178	! clouds in each
179	! model level. Same note applies as in dtau_s.
180
181	REAL frac_out(npoints,ncol,nlev) ! boxes gridbox divided up into
182	! Equivalent of BOX in original version, but
183	! indexed by column then row, rather than
184	! by row then column
185
186
187
188	! ------
189	! Output
190	! ------
191
192	REAL fq_isccp(npoints,7,7) ! the fraction of the model grid box covered by
193	! each of the 49 ISCCP D level cloud types
194
195	REAL totalcldarea(npoints) ! the fraction of model grid box columns
196	! with cloud somewhere in them. NOTE: This diagnostic
197	! does not count model clouds with tau < isccp_taumin
198	! Thus this diagnostic does not equal the sum over all entries of fq_isccp.
199	! However, this diagnostic does equal the sum over entries of fq_isccp with
200	! itau = 2:7 (omitting itau = 1)
201
202
203	! The following three means are averages only over the cloudy areas with tau > isccp_taumin.
204	! If no clouds with tau > isccp_taumin are in grid box all three quantities should equal zero.
205
206	REAL meanptop(npoints) ! mean cloud top pressure (mb) - linear averaging
207	! in cloud top pressure.
208
209	REAL meantaucld(npoints) ! mean optical thickness
210	! linear averaging in albedo performed.
211
212	real meanalbedocld(npoints) ! mean cloud albedo
213	! linear averaging in albedo performed
214
215	real meantb(npoints) ! mean all-sky 10.5 micron brightness temperature
216
217	real meantbclr(npoints) ! mean clear-sky 10.5 micron brightness temperature
218
219	REAL boxtau(npoints,ncol) ! optical thickness in each column
220
221	REAL boxptop(npoints,ncol) ! cloud top pressure (mb) in each column
222
223
224	!
225	! ------
226	! Working variables added when program updated to mimic Mark Webb's PV-Wave code
227	! ------
228
229	REAL dem(npoints,ncol),bb(npoints) ! working variables for 10.5 micron longwave
230	! emissivity in part of
231	! gridbox under consideration
232
233	REAL ptrop(npoints)
234	REAL attrop(npoints)
235	REAL attropmin (npoints)
236	REAL atmax(npoints)
237	REAL btcmin(npoints)
238	REAL transmax(npoints)
239
240	INTEGER i,j,ilev,ibox,itrop(npoints)
241	INTEGER ipres(npoints)
242	INTEGER itau(npoints),ilev2
243	INTEGER acc(nlev,ncol)
244	INTEGER match(npoints,nlev-1)
245	INTEGER nmatch(npoints)
246	INTEGER levmatch(npoints,ncol)
247
248	!variables needed for water vapor continuum absorption
249	real fluxtop_clrsky(npoints),trans_layers_above_clrsky(npoints)
250	real taumin(npoints)
251	real dem_wv(npoints,nlev), wtmair, wtmh20, Navo, grav, pstd, t0
252	real press(npoints), dpress(npoints), atmden(npoints)
253	real rvh20(npoints), wk(npoints), rhoave(npoints)
254	real rh20s(npoints), rfrgn(npoints)
255	real tmpexp(npoints),tauwv(npoints)
256
257	character*1 cchar(6),cchar_realtops(6)
258	integer icycle
259	REAL tau(npoints,ncol)
260	LOGICAL box_cloudy(npoints,ncol)
261	REAL tb(npoints,ncol)
262	REAL ptop(npoints,ncol)
263	REAL emcld(npoints,ncol)
264	REAL fluxtop(npoints,ncol)
265	REAL trans_layers_above(npoints,ncol)
266	real isccp_taumin,fluxtopinit(npoints),tauir(npoints)
267	REAL albedocld(npoints,ncol)
268	real boxarea
269	integer debug ! set to non-zero value to print out inputs
270	! with step debug
271	integer debugcol ! set to non-zero value to print out column
272	! decomposition with step debugcol
273	integer rangevec(npoints),rangeerror
274
275	integer index1(npoints),num1,jj,k1,k2
276	real rec2p13,tauchk,logp,logp1,logp2,atd
277	real output_missing_value
278
279	character*10 ftn09
280
281	DATA isccp_taumin / 0.3 /
282	DATA output_missing_value / -1.E+30 /
283	DATA cchar / ' ','-','1','+','I','+'/
284	DATA cchar_realtops / ' ',' ','1','1','I','I'/
285
286	! ------ End duplicate definitions common to wrapper routine
287
288	tauchk = -1.*log(0.9999999)
289	rec2p13=1./2.13
290
291	ncolprint=0
292
293	if ( debug.ne.0 ) then
294	j=1
295	write(6,'(a10)') 'j='
296	write(6,'(8I10)') j
297	write(6,'(a10)') 'debug='
298	write(6,'(8I10)') debug
299	write(6,'(a10)') 'debugcol='
300	write(6,'(8I10)') debugcol
301	write(6,'(a10)') 'npoints='
302	write(6,'(8I10)') npoints
303	write(6,'(a10)') 'nlev='
304	write(6,'(8I10)') nlev
305	write(6,'(a10)') 'ncol='
306	write(6,'(8I10)') ncol
307	write(6,'(a11)') 'top_height='
308	write(6,'(8I10)') top_height
309	write(6,'(a21)') 'top_height_direction='
310	write(6,'(8I10)') top_height_direction
311	write(6,'(a10)') 'overlap='
312	write(6,'(8I10)') overlap
313	write(6,'(a10)') 'emsfc_lw='
314	write(6,'(8f10.2)') emsfc_lw
315	do j=1,npoints,debug
316	write(6,'(a10)') 'j='
317	write(6,'(8I10)') j
318	write(6,'(a10)') 'sunlit='
319	write(6,'(8I10)') sunlit(j)
320	write(6,'(a10)') 'pfull='
321	write(6,'(8f10.2)') (pfull(j,i),i=1,nlev)
322	write(6,'(a10)') 'phalf='
323	write(6,'(8f10.2)') (phalf(j,i),i=1,nlev+1)
324	write(6,'(a10)') 'qv='
325	write(6,'(8f10.3)') (qv(j,i),i=1,nlev)
326	write(6,'(a10)') 'cc='
327	write(6,'(8f10.3)') (cc(j,i),i=1,nlev)
328	write(6,'(a10)') 'conv='
329	write(6,'(8f10.2)') (conv(j,i),i=1,nlev)
330	write(6,'(a10)') 'dtau_s='
331	write(6,'(8g12.5)') (dtau_s(j,i),i=1,nlev)
332	write(6,'(a10)') 'dtau_c='
333	write(6,'(8f10.2)') (dtau_c(j,i),i=1,nlev)
334	write(6,'(a10)') 'skt='
335	write(6,'(8f10.2)') skt(j)
336	write(6,'(a10)') 'at='
337	write(6,'(8f10.2)') (at(j,i),i=1,nlev)
338	write(6,'(a10)') 'dem_s='
339	write(6,'(8f10.3)') (dem_s(j,i),i=1,nlev)
340	write(6,'(a10)') 'dem_c='
341	write(6,'(8f10.3)') (dem_c(j,i),i=1,nlev)
342	enddo
343	endif
344
345	! ---------------------------------------------------!
346
347	if (ncolprint.ne.0) then
348	do j=1,npoints,1000
349	write(6,'(a10)') 'j='
350	write(6,'(8I10)') j
351	enddo
352	endif
353
354	if (top_height .eq. 1 .or. top_height .eq. 3) then
355
356	do j=1,npoints
357	ptrop(j)=5000.
358	attropmin(j) = 400.
359	atmax(j) = 0.
360	attrop(j) = 120.
361	itrop(j) = 1
362	enddo
363
364	do 12 ilev=1,nlev
365	do j=1,npoints
366	if (pfull(j,ilev) .lt. 40000. .and.
367	& pfull(j,ilev) .gt. 5000. .and.
368	& at(j,ilev) .lt. attropmin(j)) then
369	ptrop(j) = pfull(j,ilev)
370	attropmin(j) = at(j,ilev)
371	attrop(j) = attropmin(j)
372	itrop(j)=ilev
373	end if
374	enddo
375	12 continue
376
377	do 13 ilev=1,nlev
378	do j=1,npoints
379	if (at(j,ilev) .gt. atmax(j) .and.
380	& ilev .ge. itrop(j)) atmax(j)=at(j,ilev)
381	enddo
382	13 continue
383
384	end if
385
386
387	if (top_height .eq. 1 .or. top_height .eq. 3) then
388	do j=1,npoints
389	meantb(j) = 0.
390	meantbclr(j) = 0.
391	end do
392	else
393	do j=1,npoints
394	meantb(j) = output_missing_value
395	meantbclr(j) = output_missing_value
396	end do
397	end if
398
399	! -----------------------------------------------------!
400
401	! ---------------------------------------------------!
402
403	do ilev=1,nlev
404	do j=1,npoints
405
406	rangevec(j)=0
407
408	if (cc(j,ilev) .lt. 0. .or. cc(j,ilev) .gt. 1.) then
409	! error = cloud fraction less than zero
410	! error = cloud fraction greater than 1
411	rangevec(j)=rangevec(j)+1
412	endif
413
414	if (conv(j,ilev) .lt. 0. .or. conv(j,ilev) .gt. 1.) then
415	! ' error = convective cloud fraction less than zero'
416	! ' error = convective cloud fraction greater than 1'
417	rangevec(j)=rangevec(j)+2
418	endif
419
420	if (dtau_s(j,ilev) .lt. 0.) then
421	! ' error = stratiform cloud opt. depth less than zero'
422	rangevec(j)=rangevec(j)+4
423	endif
424
425	if (dtau_c(j,ilev) .lt. 0.) then
426	! ' error = convective cloud opt. depth less than zero'
427	rangevec(j)=rangevec(j)+8
428	endif
429
430	if (dem_s(j,ilev) .lt. 0. .or. dem_s(j,ilev) .gt. 1.) then
431	! ' error = stratiform cloud emissivity less than zero'
432	! ' error = stratiform cloud emissivity greater than 1'
433	rangevec(j)=rangevec(j)+16
434	endif
435
436	if (dem_c(j,ilev) .lt. 0. .or. dem_c(j,ilev) .gt. 1.) then
437	! ' error = convective cloud emissivity less than zero'
438	! ' error = convective cloud emissivity greater than 1'
439	rangevec(j)=rangevec(j)+32
440	endif
441	enddo
442
443	rangeerror=0
444	do j=1,npoints
445	rangeerror=rangeerror+rangevec(j)
446	enddo
447
448	if (rangeerror.ne.0) then
449	write (6,*) 'Input variable out of range'
450	write (6,*) 'rangevec:'
451	write (6,*) rangevec
452	STOP
453	endif
454	enddo
455
456	!
457	! ---------------------------------------------------!
458
459
460	!
461	! ---------------------------------------------------!
462	! COMPUTE CLOUD OPTICAL DEPTH FOR EACH COLUMN and
463	! put into vector tau
464
465	!initialize tau and albedocld to zero
466	do 15 ibox=1,ncol
467	do j=1,npoints
468	tau(j,ibox)=0.
469	albedocld(j,ibox)=0.
470	boxtau(j,ibox)=output_missing_value
471	boxptop(j,ibox)=output_missing_value
472	box_cloudy(j,ibox)=.false.
473	enddo
474	15 continue
475
476	!compute total cloud optical depth for each column
477	do ilev=1,nlev
478	!increment tau for each of the boxes
479	do ibox=1,ncol
480	do j=1,npoints
481	if (frac_out(j,ibox,ilev).eq.1) then
482	tau(j,ibox)=tau(j,ibox)
483	& + dtau_s(j,ilev)
484	endif
485	if (frac_out(j,ibox,ilev).eq.2) then
486	tau(j,ibox)=tau(j,ibox)
487	& + dtau_c(j,ilev)
488	end if
489	enddo
490	enddo ! ibox
491	enddo ! ilev
492	if (ncolprint.ne.0) then
493
494	do j=1,npoints ,1000
495	write(6,'(a10)') 'j='
496	write(6,'(8I10)') j
497	write(6,'(i2,1X,8(f7.2,1X))')
498	& ilev,
499	& (tau(j,ibox),ibox=1,ncolprint)
500	enddo
501	endif
502	!
503	! ---------------------------------------------------!
504
505
506
507	!
508	! ---------------------------------------------------!
509	! COMPUTE INFRARED BRIGHTNESS TEMPERUATRES
510	! AND CLOUD TOP TEMPERATURE SATELLITE SHOULD SEE
511	!
512	! again this is only done if top_height = 1 or 3
513	!
514	! fluxtop is the 10.5 micron radiance at the top of the
515	! atmosphere
516	! trans_layers_above is the total transmissivity in the layers
517	! above the current layer
518	! fluxtop_clrsky(j) and trans_layers_above_clrsky(j) are the clear
519	! sky versions of these quantities.
520
521	if (top_height .eq. 1 .or. top_height .eq. 3) then
522
523
524	!----------------------------------------------------------------------
525	!
526	! DO CLEAR SKY RADIANCE CALCULATION FIRST
527	!
528	!compute water vapor continuum emissivity
529	!this treatment follows Schwarkzopf and Ramasamy
530	!JGR 1999,vol 104, pages 9467-9499.
531	!the emissivity is calculated at a wavenumber of 955 cm-1,
532	!or 10.47 microns
533	wtmair = 28.9644
534	wtmh20 = 18.01534
535	Navo = 6.023E+23
536	grav = 9.806650E+02
537	pstd = 1.013250E+06
538	t0 = 296.
539	if (ncolprint .ne. 0)
540	& write(6,*) 'ilev pw (kg/m2) tauwv(j) dem_wv'
541	do 125 ilev=1,nlev
542	do j=1,npoints
543	!press and dpress are dyne/cm2 = Pascals *10
544	press(j) = pfull(j,ilev)*10.
545	dpress(j) = (phalf(j,ilev+1)-phalf(j,ilev))*10
546	!atmden = g/cm2 = kg/m2 / 10
547	atmden(j) = dpress(j)/grav
548	rvh20(j) = qv(j,ilev)*wtmair/wtmh20
549	wk(j) = rvh20(j)Navoatmden(j)/wtmair
550	rhoave(j) = (press(j)/pstd)*(t0/at(j,ilev))
551	rh20s(j) = rvh20(j)*rhoave(j)
552	rfrgn(j) = rhoave(j)-rh20s(j)
553	tmpexp(j) = exp(-0.02*(at(j,ilev)-t0))
554	tauwv(j) = wk(j)1.e-20(
555	& (0.0224697rh20s(j)tmpexp(j)) +
556	& (3.41817e-7rfrgn(j)) )0.98
557	dem_wv(j,ilev) = 1. - exp( -1. * tauwv(j))
558	enddo
559	if (ncolprint .ne. 0) then
560	do j=1,npoints ,1000
561	write(6,'(a10)') 'j='
562	write(6,'(8I10)') j
563	write(6,'(i2,1X,3(f8.3,3X))') ilev,
564	& qv(j,ilev)*(phalf(j,ilev+1)-phalf(j,ilev))/(grav/100.),
565	& tauwv(j),dem_wv(j,ilev)
566	enddo
567	endif
568	125 continue
569
570	!initialize variables
571	do j=1,npoints
572	fluxtop_clrsky(j) = 0.
573	trans_layers_above_clrsky(j)=1.
574	enddo
575
576	do ilev=1,nlev
577	do j=1,npoints
578
579	! Black body emission at temperature of the layer
580
581	bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. )
582	!bb(j)= 5.67e-8at(j,ilev)*4
583
584	! increase TOA flux by flux emitted from layer
585	! times total transmittance in layers above
586
587	fluxtop_clrsky(j) = fluxtop_clrsky(j)
588	& + dem_wv(j,ilev)bb(j)trans_layers_above_clrsky(j)
589
590	! update trans_layers_above with transmissivity
591	! from this layer for next time around loop
592
593	trans_layers_above_clrsky(j)=
594	& trans_layers_above_clrsky(j)*(1.-dem_wv(j,ilev))
595
596
597	enddo
598	if (ncolprint.ne.0) then
599	do j=1,npoints ,1000
600	write(6,'(a10)') 'j='
601	write(6,'(8I10)') j
602	write (6,'(a)') 'ilev:'
603	write (6,'(I2)') ilev
604
605	write (6,'(a)')
606	& 'emiss_layer,100.bb(j),100.f,total_trans:'
607	write (6,'(4(f7.2,1X))') dem_wv(j,ilev),100.*bb(j),
608	& 100.*fluxtop_clrsky(j),trans_layers_above_clrsky(j)
609	enddo
610	endif
611
612	enddo !loop over level
613
614	do j=1,npoints
615	!add in surface emission
616	bb(j)=1/( exp(1307.27/skt(j)) - 1. )
617	!bb(j)=5.67e-8skt(j)*4
618
619	fluxtop_clrsky(j) = fluxtop_clrsky(j) + emsfc_lw * bb(j)
620	& * trans_layers_above_clrsky(j)
621
622	!clear sky brightness temperature
623	meantbclr(j) = 1307.27/(log(1.+(1./fluxtop_clrsky(j))))
624
625	enddo
626
627	if (ncolprint.ne.0) then
628	do j=1,npoints ,1000
629	write(6,'(a10)') 'j='
630	write(6,'(8I10)') j
631	write (6,'(a)') 'id:'
632	write (6,'(a)') 'surface'
633
634	write (6,'(a)') 'emsfc,100.bb(j),100.f,total_trans:'
635	write (6,'(5(f7.2,1X))') emsfc_lw,100.*bb(j),
636	& 100.*fluxtop_clrsky(j),
637	& trans_layers_above_clrsky(j), meantbclr(j)
638	enddo
639	endif
640
641
642	!
643	! END OF CLEAR SKY CALCULATION
644	!
645	!----------------------------------------------------------------
646
647
648
649	if (ncolprint.ne.0) then
650
651	do j=1,npoints ,1000
652	write(6,'(a10)') 'j='
653	write(6,'(8I10)') j
654	write (6,'(a)') 'ts:'
655	write (6,'(8f7.2)') (skt(j),ibox=1,ncolprint)
656
657	write (6,'(a)') 'ta_rev:'
658	write (6,'(8f7.2)')
659	& ((at(j,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
660
661	enddo
662	endif
663	!loop over columns
664	do ibox=1,ncol
665	do j=1,npoints
666	fluxtop(j,ibox)=0.
667	trans_layers_above(j,ibox)=1.
668	enddo
669	enddo
670
671	do ilev=1,nlev
672	do j=1,npoints
673	! Black body emission at temperature of the layer
674
675	bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. )
676	!bb(j)= 5.67e-8at(j,ilev)*4
677	enddo
678
679	do ibox=1,ncol
680	do j=1,npoints
681
682	! emissivity for point in this layer
683	if (frac_out(j,ibox,ilev).eq.1) then
684	dem(j,ibox)= 1. -
685	& ( (1. - dem_wv(j,ilev)) * (1. - dem_s(j,ilev)) )
686	else if (frac_out(j,ibox,ilev).eq.2) then
687	dem(j,ibox)= 1. -
688	& ( (1. - dem_wv(j,ilev)) * (1. - dem_c(j,ilev)) )
689	else
690	dem(j,ibox)= dem_wv(j,ilev)
691	end if
692
693
694	! increase TOA flux by flux emitted from layer
695	! times total transmittance in layers above
696
697	fluxtop(j,ibox) = fluxtop(j,ibox)
698	& + dem(j,ibox) * bb(j)
699	& * trans_layers_above(j,ibox)
700
701	! update trans_layers_above with transmissivity
702	! from this layer for next time around loop
703
704	trans_layers_above(j,ibox)=
705	& trans_layers_above(j,ibox)*(1.-dem(j,ibox))
706
707	enddo ! j
708	enddo ! ibox
709
710	if (ncolprint.ne.0) then
711	do j=1,npoints,1000
712	write (6,'(a)') 'ilev:'
713	write (6,'(I2)') ilev
714
715	write(6,'(a10)') 'j='
716	write(6,'(8I10)') j
717	write (6,'(a)') 'emiss_layer:'
718	write (6,'(8f7.2)') (dem(j,ibox),ibox=1,ncolprint)
719
720	write (6,'(a)') '100.*bb(j):'
721	write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
722
723	write (6,'(a)') '100.*f:'
724	write (6,'(8f7.2)')
725	& (100.*fluxtop(j,ibox),ibox=1,ncolprint)
726
727	write (6,'(a)') 'total_trans:'
728	write (6,'(8f7.2)')
729	& (trans_layers_above(j,ibox),ibox=1,ncolprint)
730	enddo
731	endif
732
733	enddo ! ilev
734
735
736	do j=1,npoints
737	!add in surface emission
738	bb(j)=1/( exp(1307.27/skt(j)) - 1. )
739	!bb(j)=5.67e-8skt(j)*4
740	end do
741
742	do ibox=1,ncol
743	do j=1,npoints
744
745	!add in surface emission
746
747	fluxtop(j,ibox) = fluxtop(j,ibox)
748	& + emsfc_lw * bb(j)
749	& * trans_layers_above(j,ibox)
750
751	end do
752	end do
753
754	!calculate mean infrared brightness temperature
755	do ibox=1,ncol
756	do j=1,npoints
757	meantb(j) = meantb(j)+1307.27/(log(1.+(1./fluxtop(j,ibox))))
758	end do
759	end do
760	do j=1, npoints
761	meantb(j) = meantb(j) / real(ncol)
762	end do
763
764	if (ncolprint.ne.0) then
765
766	do j=1,npoints ,1000
767	write(6,'(a10)') 'j='
768	write(6,'(8I10)') j
769	write (6,'(a)') 'id:'
770	write (6,'(a)') 'surface'
771
772	write (6,'(a)') 'emiss_layer:'
773	write (6,'(8f7.2)') (dem(1,ibox),ibox=1,ncolprint)
774
775	write (6,'(a)') '100.*bb(j):'
776	write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
777
778	write (6,'(a)') '100.*f:'
779	write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint)
780
781	write (6,'(a)') 'meantb(j):'
782	write (6,'(8f7.2)') (meantb(j),ibox=1,ncolprint)
783
784	end do
785	endif
786
787	!now that you have the top of atmosphere radiance account
788	!for ISCCP procedures to determine cloud top temperature
789
790	!account for partially transmitting cloud recompute flux
791	!ISCCP would see assuming a single layer cloud
792	!note choice here of 2.13, as it is primarily ice
793	!clouds which have partial emissivity and need the
794	!adjustment performed in this section
795	!
796	!If it turns out that the cloud brightness temperature
797	!is greater than 260K, then the liquid cloud conversion
798	!factor of 2.56 is used.
799	!
800	!Note that this is discussed on pages 85-87 of
801	!the ISCCP D level documentation (Rossow et al. 1996)
802
803	do j=1,npoints
804	!compute minimum brightness temperature and optical depth
805	btcmin(j) = 1. / ( exp(1307.27/(attrop(j)-5.)) - 1. )
806	enddo
807	do ibox=1,ncol
808	do j=1,npoints
809	transmax(j) = (fluxtop(j,ibox)-btcmin(j))
810	& /(fluxtop_clrsky(j)-btcmin(j))
811	!note that the initial setting of tauir(j) is needed so that
812	!tauir(j) has a realistic value should the next if block be
813	!bypassed
814	tauir(j) = tau(j,ibox) * rec2p13
815	taumin(j) = -1. * log(max(min(transmax(j),0.9999999),0.001))
816
817	enddo
818
819	if (top_height .eq. 1) then
820	do j=1,npoints
821	if (transmax(j) .gt. 0.001 .and.
822	& transmax(j) .le. 0.9999999) then
823	fluxtopinit(j) = fluxtop(j,ibox)
824	tauir(j) = tau(j,ibox) *rec2p13
825	endif
826	enddo
827	do icycle=1,2
828	do j=1,npoints
829	if (tau(j,ibox) .gt. (tauchk )) then
830	if (transmax(j) .gt. 0.001 .and.
831	& transmax(j) .le. 0.9999999) then
832	emcld(j,ibox) = 1. - exp(-1. * tauir(j) )
833	fluxtop(j,ibox) = fluxtopinit(j) -
834	& ((1.-emcld(j,ibox))*fluxtop_clrsky(j))
835	fluxtop(j,ibox)=max(1.E-06,
836	& (fluxtop(j,ibox)/emcld(j,ibox)))
837	tb(j,ibox)= 1307.27
838	& / (log(1. + (1./fluxtop(j,ibox))))
839	if (tb(j,ibox) .gt. 260.) then
840	tauir(j) = tau(j,ibox) / 2.56
841	end if
842	end if
843	end if
844	enddo
845	enddo
846
847	endif
848
849	do j=1,npoints
850	if (tau(j,ibox) .gt. (tauchk )) then
851	!cloudy box
852	!NOTE: tb is the cloud-top temperature not infrared brightness temperature
853	!at this point in the code
854	tb(j,ibox)= 1307.27/ (log(1. + (1./fluxtop(j,ibox))))
855	if (top_height.eq.1.and.tauir(j).lt.taumin(j)) then
856	tb(j,ibox) = attrop(j) - 5.
857	tau(j,ibox) = 2.13*taumin(j)
858	end if
859	else
860	!clear sky brightness temperature
861	tb(j,ibox) = meantbclr(j)
862	end if
863	enddo ! j
864	enddo ! ibox
865
866	if (ncolprint.ne.0) then
867
868	do j=1,npoints,1000
869	write(6,'(a10)') 'j='
870	write(6,'(8I10)') j
871
872	write (6,'(a)') 'attrop:'
873	write (6,'(8f7.2)') (attrop(j))
874
875	write (6,'(a)') 'btcmin:'
876	write (6,'(8f7.2)') (btcmin(j))
877
878	write (6,'(a)') 'fluxtop_clrsky*100:'
879	write (6,'(8f7.2)')
880	& (100.*fluxtop_clrsky(j))
881
882	write (6,'(a)') '100.*f_adj:'
883	write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint)
884
885	write (6,'(a)') 'transmax:'
886	write (6,'(8f7.2)') (transmax(ibox),ibox=1,ncolprint)
887
888	write (6,'(a)') 'tau:'
889	write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
890
891	write (6,'(a)') 'emcld:'
892	write (6,'(8f7.2)') (emcld(j,ibox),ibox=1,ncolprint)
893
894	write (6,'(a)') 'total_trans:'
895	write (6,'(8f7.2)')
896	& (trans_layers_above(j,ibox),ibox=1,ncolprint)
897
898	write (6,'(a)') 'total_emiss:'
899	write (6,'(8f7.2)')
900	& (1.0-trans_layers_above(j,ibox),ibox=1,ncolprint)
901
902	write (6,'(a)') 'total_trans:'
903	write (6,'(8f7.2)')
904	& (trans_layers_above(j,ibox),ibox=1,ncolprint)
905
906	write (6,'(a)') 'ppout:'
907	write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
908	enddo ! j
909	endif
910
911	end if
912
913	! ---------------------------------------------------!
914
915	!
916	! ---------------------------------------------------!
917	! DETERMINE CLOUD TOP PRESSURE
918	!
919	! again the 2 methods differ according to whether
920	! or not you use the physical cloud top pressure (top_height = 2)
921	! or the radiatively determined cloud top pressure (top_height = 1 or 3)
922	!
923
924	!compute cloud top pressure
925	do 30 ibox=1,ncol
926	!segregate according to optical thickness
927	if (top_height .eq. 1 .or. top_height .eq. 3) then
928	!find level whose temperature
929	!most closely matches brightness temperature
930	do j=1,npoints
931	nmatch(j)=0
932	enddo
933	do 29 k1=1,nlev-1
934	if (top_height_direction .eq. 2) then
935	ilev = nlev - k1
936	else
937	ilev = k1
938	end if
939	!cdir nodep
940	do j=1,npoints
941	if (ilev .ge. itrop(j)) then
942	if ((at(j,ilev) .ge. tb(j,ibox) .and.
943	& at(j,ilev+1) .le. tb(j,ibox)) .or.
944	& (at(j,ilev) .le. tb(j,ibox) .and.
945	& at(j,ilev+1) .ge. tb(j,ibox))) then
946	nmatch(j)=nmatch(j)+1
947	match(j,nmatch(j))=ilev
948	end if
949	end if
950	enddo
951	29 continue
952
953	do j=1,npoints
954	if (nmatch(j) .ge. 1) then
955	k1 = match(j,nmatch(j))
956	k2 = k1 + 1
957	logp1 = log(pfull(j,k1))
958	logp2 = log(pfull(j,k2))
959	atd = max(tauchk,abs(at(j,k2) - at(j,k1)))
960	logp=logp1+(logp2-logp1)*abs(tb(j,ibox)-at(j,k1))/atd
961	ptop(j,ibox) = exp(logp)
962	if(abs(pfull(j,k1)-ptop(j,ibox)) .lt.
963	& abs(pfull(j,k2)-ptop(j,ibox))) then
964	levmatch(j,ibox)=k1
965	else
966	levmatch(j,ibox)=k2
967	end if
968	else
969	if (tb(j,ibox) .le. attrop(j)) then
970	ptop(j,ibox)=ptrop(j)
971	levmatch(j,ibox)=itrop(j)
972	end if
973	if (tb(j,ibox) .ge. atmax(j)) then
974	ptop(j,ibox)=pfull(j,nlev)
975	levmatch(j,ibox)=nlev
976	end if
977	end if
978	enddo ! j
979
980	else ! if (top_height .eq. 1 .or. top_height .eq. 3)
981
982	do j=1,npoints
983	ptop(j,ibox)=0.
984	enddo
985	do ilev=1,nlev
986	do j=1,npoints
987	if ((ptop(j,ibox) .eq. 0. )
988	& .and.(frac_out(j,ibox,ilev) .ne. 0)) then
989	ptop(j,ibox)=phalf(j,ilev)
990	levmatch(j,ibox)=ilev
991	end if
992	end do
993	end do
994	end if
995
996	do j=1,npoints
997	if (tau(j,ibox) .le. (tauchk )) then
998	ptop(j,ibox)=0.
999	levmatch(j,ibox)=0
1000	endif
1001	enddo
1002
1003	30 continue
1004
1005	!
1006	!
1007	! ---------------------------------------------------!
1008
1009
1010	!
1011	! ---------------------------------------------------!
1012	! DETERMINE ISCCP CLOUD TYPE FREQUENCIES
1013	!
1014	! Now that ptop and tau have been determined,
1015	! determine amount of each of the 49 ISCCP cloud
1016	! types
1017	!
1018	! Also compute grid box mean cloud top pressure and
1019	! optical thickness. The mean cloud top pressure and
1020	! optical thickness are averages over the cloudy
1021	! area only. The mean cloud top pressure is a linear
1022	! average of the cloud top pressures. The mean cloud
1023	! optical thickness is computed by converting optical
1024	! thickness to an albedo, averaging in albedo units,
1025	! then converting the average albedo back to a mean
1026	! optical thickness.
1027	!
1028
1029	!compute isccp frequencies
1030
1031	!reset frequencies
1032	do 38 ilev=1,7
1033	do 38 ilev2=1,7
1034	do j=1,npoints !
1035	if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1036	fq_isccp(j,ilev,ilev2)= 0.
1037	else
1038	fq_isccp(j,ilev,ilev2)= output_missing_value
1039	end if
1040	enddo
1041	38 continue
1042
1043	!reset variables need for averaging cloud properties
1044	do j=1,npoints
1045	if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1046	totalcldarea(j) = 0.
1047	meanalbedocld(j) = 0.
1048	meanptop(j) = 0.
1049	meantaucld(j) = 0.
1050	else
1051	totalcldarea(j) = output_missing_value
1052	meanalbedocld(j) = output_missing_value
1053	meanptop(j) = output_missing_value
1054	meantaucld(j) = output_missing_value
1055	end if
1056	enddo ! j
1057
1058	boxarea = 1./real(ncol)
1059
1060	do 39 ibox=1,ncol
1061	do j=1,npoints
1062
1063	if (tau(j,ibox) .gt. (tauchk )
1064	& .and. ptop(j,ibox) .gt. 0.) then
1065	box_cloudy(j,ibox)=.true.
1066	endif
1067
1068	if (box_cloudy(j,ibox)) then
1069
1070	if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1071
1072	boxtau(j,ibox) = tau(j,ibox)
1073
1074	if (tau(j,ibox) .ge. isccp_taumin) then
1075	totalcldarea(j) = totalcldarea(j) + boxarea
1076
1077	!convert optical thickness to albedo
1078	albedocld(j,ibox)
1079	& = (tau(j,ibox)0.895)/((tau(j,ibox)0.895)+6.82)
1080
1081	!contribute to averaging
1082	meanalbedocld(j) = meanalbedocld(j)
1083	& +albedocld(j,ibox)*boxarea
1084
1085	end if
1086
1087	endif
1088
1089	endif
1090
1091	if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1092
1093	if (box_cloudy(j,ibox)) then
1094
1095	!convert ptop to millibars
1096	ptop(j,ibox)=ptop(j,ibox) / 100.
1097
1098	!save for output cloud top pressure and optical thickness
1099	boxptop(j,ibox) = ptop(j,ibox)
1100
1101	if (tau(j,ibox) .ge. isccp_taumin) then
1102	meanptop(j) = meanptop(j) + ptop(j,ibox)*boxarea
1103	end if
1104
1105	!reset itau(j), ipres(j)
1106	itau(j) = 0
1107	ipres(j) = 0
1108
1109	!determine optical depth category
1110	if (tau(j,ibox) .lt. isccp_taumin) then
1111	itau(j)=1
1112	else if (tau(j,ibox) .ge. isccp_taumin
1113	&
1114	& .and. tau(j,ibox) .lt. 1.3) then
1115	itau(j)=2
1116	else if (tau(j,ibox) .ge. 1.3
1117	& .and. tau(j,ibox) .lt. 3.6) then
1118	itau(j)=3
1119	else if (tau(j,ibox) .ge. 3.6
1120	& .and. tau(j,ibox) .lt. 9.4) then
1121	itau(j)=4
1122	else if (tau(j,ibox) .ge. 9.4
1123	& .and. tau(j,ibox) .lt. 23.) then
1124	itau(j)=5
1125	else if (tau(j,ibox) .ge. 23.
1126	& .and. tau(j,ibox) .lt. 60.) then
1127	itau(j)=6
1128	else if (tau(j,ibox) .ge. 60.) then
1129	itau(j)=7
1130	end if
1131
1132	!determine cloud top pressure category
1133	if ( ptop(j,ibox) .gt. 0.
1134	& .and.ptop(j,ibox) .lt. 180.) then
1135	ipres(j)=1
1136	else if(ptop(j,ibox) .ge. 180.
1137	& .and.ptop(j,ibox) .lt. 310.) then
1138	ipres(j)=2
1139	else if(ptop(j,ibox) .ge. 310.
1140	& .and.ptop(j,ibox) .lt. 440.) then
1141	ipres(j)=3
1142	else if(ptop(j,ibox) .ge. 440.
1143	& .and.ptop(j,ibox) .lt. 560.) then
1144	ipres(j)=4
1145	else if(ptop(j,ibox) .ge. 560.
1146	& .and.ptop(j,ibox) .lt. 680.) then
1147	ipres(j)=5
1148	else if(ptop(j,ibox) .ge. 680.
1149	& .and.ptop(j,ibox) .lt. 800.) then
1150	ipres(j)=6
1151	else if(ptop(j,ibox) .ge. 800.) then
1152	ipres(j)=7
1153	end if
1154
1155	!update frequencies
1156	if(ipres(j) .gt. 0.and.itau(j) .gt. 0) then
1157	fq_isccp(j,itau(j),ipres(j))=
1158	& fq_isccp(j,itau(j),ipres(j))+ boxarea
1159	end if
1160
1161	end if
1162
1163	end if
1164
1165	enddo ! j
1166	39 continue
1167
1168	!compute mean cloud properties
1169	do j=1,npoints
1170	if (totalcldarea(j) .gt. 0.) then
1171	! code above guarantees that totalcldarea > 0
1172	! only if sunlit .eq. 1 .or. top_height = 3
1173	! and applies only to clouds with tau > isccp_taumin
1174	meanptop(j) = meanptop(j) / totalcldarea(j)
1175	meanalbedocld(j) = meanalbedocld(j) / totalcldarea(j)
1176	meantaucld(j) = (6.82/((1./meanalbedocld(j))-1.))**(1./0.895)
1177	else
1178	! this code is necessary so that in the case that totalcldarea = 0.,
1179	! that these variables, which are in-cloud averages, are set to missing
1180	! note that totalcldarea will be 0. if all the clouds in the grid box have
1181	! tau < isccp_taumin
1182	meanptop(j) = output_missing_value
1183	meanalbedocld(j) = output_missing_value
1184	meantaucld(j) = output_missing_value
1185	end if
1186	enddo ! j
1187	!
1188	! ---------------------------------------------------!
1189
1190	! ---------------------------------------------------!
1191	! OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM
1192	!
1193	if (debugcol.ne.0) then
1194	!
1195	do j=1,npoints,debugcol
1196
1197	!produce character output
1198	do ilev=1,nlev
1199	do ibox=1,ncol
1200	acc(ilev,ibox)=0
1201	enddo
1202	enddo
1203
1204	do ilev=1,nlev
1205	do ibox=1,ncol
1206	acc(ilev,ibox)=frac_out(j,ibox,ilev)*2
1207	if (levmatch(j,ibox) .eq. ilev)
1208	& acc(ilev,ibox)=acc(ilev,ibox)+1
1209	enddo
1210	enddo
1211
1212	!print test
1213
1214	write(ftn09,11) j
1215	11 format('ftn09.',i4.4)
1216	open(9, FILE=ftn09, FORM='FORMATTED')
1217
1218	write(9,'(a1)') ' '
1219	write(9,'(10i5)')
1220	& (ilev,ilev=5,nlev,5)
1221	write(9,'(a1)') ' '
1222
1223	do ibox=1,ncol
1224	write(9,'(40(a1),1x,40(a1))')
1225	& (cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev)
1226	& ,(cchar(acc(ilev,ibox)+1),ilev=1,nlev)
1227	end do
1228	close(9)
1229
1230	if (ncolprint.ne.0) then
1231	write(6,'(a1)') ' '
1232	write(6,'(a2,1X,5(a7,1X),a50)')
1233	& 'ilev',
1234	& 'pfull','at',
1235	& 'cc*100','dem_s','dtau_s',
1236	& 'cchar'
1237
1238	! do 4012 ilev=1,nlev
1239	! write(6,'(60i2)') (box(i,ilev),i=1,ncolprint)
1240	! write(6,'(i2,1X,5(f7.2,1X),50(a1))')
1241	! & ilev,
1242	! & pfull(j,ilev)/100.,at(j,ilev),
1243	! & cc(j,ilev)*100.0,dem_s(j,ilev),dtau_s(j,ilev)
1244	! & ,(cchar(acc(ilev,ibox)+1),ibox=1,ncolprint)
1245	!4012 continue
1246	write (6,'(a)') 'skt(j):'
1247	write (6,'(8f7.2)') skt(j)
1248
1249	write (6,'(8I7)') (ibox,ibox=1,ncolprint)
1250
1251	write (6,'(a)') 'tau:'
1252	write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
1253
1254	write (6,'(a)') 'tb:'
1255	write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
1256
1257	write (6,'(a)') 'ptop:'
1258	write (6,'(8f7.2)') (ptop(j,ibox),ibox=1,ncolprint)
1259	endif
1260
1261	enddo
1262
1263	end if
1264
1265	return
1266	end
1267
1268

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

Download in other formats:

Original Format