Context Navigation

icarus.F @ 4024

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

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

Download in other formats:

Original Format