Context Navigation

icarus.F @ 5423

Last change on this file since 5423 was 1414, checked in by idelkadi, 14 years ago
Passage a la version cosp.v1.3 pour le Lidar et ISCCP Corrections de bugs pour ISCCP et optimisation pour le Lidar
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	call flush(6)
453	STOP
454	endif
455	enddo
456
457	!
458	! ---------------------------------------------------!
459
460
461	!
462	! ---------------------------------------------------!
463	! COMPUTE CLOUD OPTICAL DEPTH FOR EACH COLUMN and
464	! put into vector tau
465
466	!initialize tau and albedocld to zero
467	do 15 ibox=1,ncol
468	do j=1,npoints
469	tau(j,ibox)=0.
470	albedocld(j,ibox)=0.
471	boxtau(j,ibox)=output_missing_value
472	boxptop(j,ibox)=output_missing_value
473	box_cloudy(j,ibox)=.false.
474	enddo
475	15 continue
476
477	!compute total cloud optical depth for each column
478	do ilev=1,nlev
479	!increment tau for each of the boxes
480	do ibox=1,ncol
481	do j=1,npoints
482	if (frac_out(j,ibox,ilev).eq.1) then
483	tau(j,ibox)=tau(j,ibox)
484	& + dtau_s(j,ilev)
485	endif
486	if (frac_out(j,ibox,ilev).eq.2) then
487	tau(j,ibox)=tau(j,ibox)
488	& + dtau_c(j,ilev)
489	end if
490	enddo
491	enddo ! ibox
492	enddo ! ilev
493	if (ncolprint.ne.0) then
494
495	do j=1,npoints ,1000
496	write(6,'(a10)') 'j='
497	write(6,'(8I10)') j
498	write(6,'(i2,1X,8(f7.2,1X))')
499	& ilev,
500	& (tau(j,ibox),ibox=1,ncolprint)
501	enddo
502	endif
503	!
504	! ---------------------------------------------------!
505
506
507
508	!
509	! ---------------------------------------------------!
510	! COMPUTE INFRARED BRIGHTNESS TEMPERUATRES
511	! AND CLOUD TOP TEMPERATURE SATELLITE SHOULD SEE
512	!
513	! again this is only done if top_height = 1 or 3
514	!
515	! fluxtop is the 10.5 micron radiance at the top of the
516	! atmosphere
517	! trans_layers_above is the total transmissivity in the layers
518	! above the current layer
519	! fluxtop_clrsky(j) and trans_layers_above_clrsky(j) are the clear
520	! sky versions of these quantities.
521
522	if (top_height .eq. 1 .or. top_height .eq. 3) then
523
524
525	!----------------------------------------------------------------------
526	!
527	! DO CLEAR SKY RADIANCE CALCULATION FIRST
528	!
529	!compute water vapor continuum emissivity
530	!this treatment follows Schwarkzopf and Ramasamy
531	!JGR 1999,vol 104, pages 9467-9499.
532	!the emissivity is calculated at a wavenumber of 955 cm-1,
533	!or 10.47 microns
534	wtmair = 28.9644
535	wtmh20 = 18.01534
536	Navo = 6.023E+23
537	grav = 9.806650E+02
538	pstd = 1.013250E+06
539	t0 = 296.
540	if (ncolprint .ne. 0)
541	& write(6,*) 'ilev pw (kg/m2) tauwv(j) dem_wv'
542	do 125 ilev=1,nlev
543	do j=1,npoints
544	!press and dpress are dyne/cm2 = Pascals *10
545	press(j) = pfull(j,ilev)*10.
546	dpress(j) = (phalf(j,ilev+1)-phalf(j,ilev))*10
547	!atmden = g/cm2 = kg/m2 / 10
548	atmden(j) = dpress(j)/grav
549	rvh20(j) = qv(j,ilev)*wtmair/wtmh20
550	wk(j) = rvh20(j)Navoatmden(j)/wtmair
551	rhoave(j) = (press(j)/pstd)*(t0/at(j,ilev))
552	rh20s(j) = rvh20(j)*rhoave(j)
553	rfrgn(j) = rhoave(j)-rh20s(j)
554	tmpexp(j) = exp(-0.02*(at(j,ilev)-t0))
555	tauwv(j) = wk(j)1.e-20(
556	& (0.0224697rh20s(j)tmpexp(j)) +
557	& (3.41817e-7rfrgn(j)) )0.98
558	dem_wv(j,ilev) = 1. - exp( -1. * tauwv(j))
559	enddo
560	if (ncolprint .ne. 0) then
561	do j=1,npoints ,1000
562	write(6,'(a10)') 'j='
563	write(6,'(8I10)') j
564	write(6,'(i2,1X,3(f8.3,3X))') ilev,
565	& qv(j,ilev)*(phalf(j,ilev+1)-phalf(j,ilev))/(grav/100.),
566	& tauwv(j),dem_wv(j,ilev)
567	enddo
568	endif
569	125 continue
570
571	!initialize variables
572	do j=1,npoints
573	fluxtop_clrsky(j) = 0.
574	trans_layers_above_clrsky(j)=1.
575	enddo
576
577	do ilev=1,nlev
578	do j=1,npoints
579
580	! Black body emission at temperature of the layer
581
582	bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. )
583	!bb(j)= 5.67e-8at(j,ilev)*4
584
585	! increase TOA flux by flux emitted from layer
586	! times total transmittance in layers above
587
588	fluxtop_clrsky(j) = fluxtop_clrsky(j)
589	& + dem_wv(j,ilev)bb(j)trans_layers_above_clrsky(j)
590
591	! update trans_layers_above with transmissivity
592	! from this layer for next time around loop
593
594	trans_layers_above_clrsky(j)=
595	& trans_layers_above_clrsky(j)*(1.-dem_wv(j,ilev))
596
597
598	enddo
599	if (ncolprint.ne.0) then
600	do j=1,npoints ,1000
601	write(6,'(a10)') 'j='
602	write(6,'(8I10)') j
603	write (6,'(a)') 'ilev:'
604	write (6,'(I2)') ilev
605
606	write (6,'(a)')
607	& 'emiss_layer,100.bb(j),100.f,total_trans:'
608	write (6,'(4(f7.2,1X))') dem_wv(j,ilev),100.*bb(j),
609	& 100.*fluxtop_clrsky(j),trans_layers_above_clrsky(j)
610	enddo
611	endif
612
613	enddo !loop over level
614
615	do j=1,npoints
616	!add in surface emission
617	bb(j)=1/( exp(1307.27/skt(j)) - 1. )
618	!bb(j)=5.67e-8skt(j)*4
619
620	fluxtop_clrsky(j) = fluxtop_clrsky(j) + emsfc_lw * bb(j)
621	& * trans_layers_above_clrsky(j)
622
623	!clear sky brightness temperature
624	meantbclr(j) = 1307.27/(log(1.+(1./fluxtop_clrsky(j))))
625
626	enddo
627
628	if (ncolprint.ne.0) then
629	do j=1,npoints ,1000
630	write(6,'(a10)') 'j='
631	write(6,'(8I10)') j
632	write (6,'(a)') 'id:'
633	write (6,'(a)') 'surface'
634
635	write (6,'(a)') 'emsfc,100.bb(j),100.f,total_trans:'
636	write (6,'(5(f7.2,1X))') emsfc_lw,100.*bb(j),
637	& 100.*fluxtop_clrsky(j),
638	& trans_layers_above_clrsky(j), meantbclr(j)
639	enddo
640	endif
641
642
643	!
644	! END OF CLEAR SKY CALCULATION
645	!
646	!----------------------------------------------------------------
647
648
649
650	if (ncolprint.ne.0) then
651
652	do j=1,npoints ,1000
653	write(6,'(a10)') 'j='
654	write(6,'(8I10)') j
655	write (6,'(a)') 'ts:'
656	write (6,'(8f7.2)') (skt(j),ibox=1,ncolprint)
657
658	write (6,'(a)') 'ta_rev:'
659	write (6,'(8f7.2)')
660	& ((at(j,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
661
662	enddo
663	endif
664	!loop over columns
665	do ibox=1,ncol
666	do j=1,npoints
667	fluxtop(j,ibox)=0.
668	trans_layers_above(j,ibox)=1.
669	enddo
670	enddo
671
672	do ilev=1,nlev
673	do j=1,npoints
674	! Black body emission at temperature of the layer
675
676	bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. )
677	!bb(j)= 5.67e-8at(j,ilev)*4
678	enddo
679
680	do ibox=1,ncol
681	do j=1,npoints
682
683	! emissivity for point in this layer
684	if (frac_out(j,ibox,ilev).eq.1) then
685	dem(j,ibox)= 1. -
686	& ( (1. - dem_wv(j,ilev)) * (1. - dem_s(j,ilev)) )
687	else if (frac_out(j,ibox,ilev).eq.2) then
688	dem(j,ibox)= 1. -
689	& ( (1. - dem_wv(j,ilev)) * (1. - dem_c(j,ilev)) )
690	else
691	dem(j,ibox)= dem_wv(j,ilev)
692	end if
693
694
695	! increase TOA flux by flux emitted from layer
696	! times total transmittance in layers above
697
698	fluxtop(j,ibox) = fluxtop(j,ibox)
699	& + dem(j,ibox) * bb(j)
700	& * trans_layers_above(j,ibox)
701
702	! update trans_layers_above with transmissivity
703	! from this layer for next time around loop
704
705	trans_layers_above(j,ibox)=
706	& trans_layers_above(j,ibox)*(1.-dem(j,ibox))
707
708	enddo ! j
709	enddo ! ibox
710
711	if (ncolprint.ne.0) then
712	do j=1,npoints,1000
713	write (6,'(a)') 'ilev:'
714	write (6,'(I2)') ilev
715
716	write(6,'(a10)') 'j='
717	write(6,'(8I10)') j
718	write (6,'(a)') 'emiss_layer:'
719	write (6,'(8f7.2)') (dem(j,ibox),ibox=1,ncolprint)
720
721	write (6,'(a)') '100.*bb(j):'
722	write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
723
724	write (6,'(a)') '100.*f:'
725	write (6,'(8f7.2)')
726	& (100.*fluxtop(j,ibox),ibox=1,ncolprint)
727
728	write (6,'(a)') 'total_trans:'
729	write (6,'(8f7.2)')
730	& (trans_layers_above(j,ibox),ibox=1,ncolprint)
731	enddo
732	endif
733
734	enddo ! ilev
735
736
737	do j=1,npoints
738	!add in surface emission
739	bb(j)=1/( exp(1307.27/skt(j)) - 1. )
740	!bb(j)=5.67e-8skt(j)*4
741	end do
742
743	do ibox=1,ncol
744	do j=1,npoints
745
746	!add in surface emission
747
748	fluxtop(j,ibox) = fluxtop(j,ibox)
749	& + emsfc_lw * bb(j)
750	& * trans_layers_above(j,ibox)
751
752	end do
753	end do
754
755	!calculate mean infrared brightness temperature
756	do ibox=1,ncol
757	do j=1,npoints
758	meantb(j) = meantb(j)+1307.27/(log(1.+(1./fluxtop(j,ibox))))
759	end do
760	end do
761	do j=1, npoints
762	meantb(j) = meantb(j) / real(ncol)
763	end do
764
765	if (ncolprint.ne.0) then
766
767	do j=1,npoints ,1000
768	write(6,'(a10)') 'j='
769	write(6,'(8I10)') j
770	write (6,'(a)') 'id:'
771	write (6,'(a)') 'surface'
772
773	write (6,'(a)') 'emiss_layer:'
774	write (6,'(8f7.2)') (dem(1,ibox),ibox=1,ncolprint)
775
776	write (6,'(a)') '100.*bb(j):'
777	write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
778
779	write (6,'(a)') '100.*f:'
780	write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint)
781
782	write (6,'(a)') 'meantb(j):'
783	write (6,'(8f7.2)') (meantb(j),ibox=1,ncolprint)
784
785	end do
786	endif
787
788	!now that you have the top of atmosphere radiance account
789	!for ISCCP procedures to determine cloud top temperature
790
791	!account for partially transmitting cloud recompute flux
792	!ISCCP would see assuming a single layer cloud
793	!note choice here of 2.13, as it is primarily ice
794	!clouds which have partial emissivity and need the
795	!adjustment performed in this section
796	!
797	!If it turns out that the cloud brightness temperature
798	!is greater than 260K, then the liquid cloud conversion
799	!factor of 2.56 is used.
800	!
801	!Note that this is discussed on pages 85-87 of
802	!the ISCCP D level documentation (Rossow et al. 1996)
803
804	do j=1,npoints
805	!compute minimum brightness temperature and optical depth
806	btcmin(j) = 1. / ( exp(1307.27/(attrop(j)-5.)) - 1. )
807	enddo
808	do ibox=1,ncol
809	do j=1,npoints
810	transmax(j) = (fluxtop(j,ibox)-btcmin(j))
811	& /(fluxtop_clrsky(j)-btcmin(j))
812	!note that the initial setting of tauir(j) is needed so that
813	!tauir(j) has a realistic value should the next if block be
814	!bypassed
815	tauir(j) = tau(j,ibox) * rec2p13
816	taumin(j) = -1. * log(max(min(transmax(j),0.9999999),0.001))
817
818	enddo
819
820	if (top_height .eq. 1) then
821	do j=1,npoints
822	if (transmax(j) .gt. 0.001 .and.
823	& transmax(j) .le. 0.9999999) then
824	fluxtopinit(j) = fluxtop(j,ibox)
825	tauir(j) = tau(j,ibox) *rec2p13
826	endif
827	enddo
828	do icycle=1,2
829	do j=1,npoints
830	if (tau(j,ibox) .gt. (tauchk )) then
831	if (transmax(j) .gt. 0.001 .and.
832	& transmax(j) .le. 0.9999999) then
833	emcld(j,ibox) = 1. - exp(-1. * tauir(j) )
834	fluxtop(j,ibox) = fluxtopinit(j) -
835	& ((1.-emcld(j,ibox))*fluxtop_clrsky(j))
836	fluxtop(j,ibox)=max(1.E-06,
837	& (fluxtop(j,ibox)/emcld(j,ibox)))
838	tb(j,ibox)= 1307.27
839	& / (log(1. + (1./fluxtop(j,ibox))))
840	if (tb(j,ibox) .gt. 260.) then
841	tauir(j) = tau(j,ibox) / 2.56
842	end if
843	end if
844	end if
845	enddo
846	enddo
847
848	endif
849
850	do j=1,npoints
851	if (tau(j,ibox) .gt. (tauchk )) then
852	!cloudy box
853	!NOTE: tb is the cloud-top temperature not infrared brightness temperature
854	!at this point in the code
855	tb(j,ibox)= 1307.27/ (log(1. + (1./fluxtop(j,ibox))))
856	if (top_height.eq.1.and.tauir(j).lt.taumin(j)) then
857	tb(j,ibox) = attrop(j) - 5.
858	tau(j,ibox) = 2.13*taumin(j)
859	end if
860	else
861	!clear sky brightness temperature
862	tb(j,ibox) = meantbclr(j)
863	end if
864	enddo ! j
865	enddo ! ibox
866
867	if (ncolprint.ne.0) then
868
869	do j=1,npoints,1000
870	write(6,'(a10)') 'j='
871	write(6,'(8I10)') j
872
873	write (6,'(a)') 'attrop:'
874	write (6,'(8f7.2)') (attrop(j))
875
876	write (6,'(a)') 'btcmin:'
877	write (6,'(8f7.2)') (btcmin(j))
878
879	write (6,'(a)') 'fluxtop_clrsky*100:'
880	write (6,'(8f7.2)')
881	& (100.*fluxtop_clrsky(j))
882
883	write (6,'(a)') '100.*f_adj:'
884	write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint)
885
886	write (6,'(a)') 'transmax:'
887	write (6,'(8f7.2)') (transmax(ibox),ibox=1,ncolprint)
888
889	write (6,'(a)') 'tau:'
890	write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
891
892	write (6,'(a)') 'emcld:'
893	write (6,'(8f7.2)') (emcld(j,ibox),ibox=1,ncolprint)
894
895	write (6,'(a)') 'total_trans:'
896	write (6,'(8f7.2)')
897	& (trans_layers_above(j,ibox),ibox=1,ncolprint)
898
899	write (6,'(a)') 'total_emiss:'
900	write (6,'(8f7.2)')
901	& (1.0-trans_layers_above(j,ibox),ibox=1,ncolprint)
902
903	write (6,'(a)') 'total_trans:'
904	write (6,'(8f7.2)')
905	& (trans_layers_above(j,ibox),ibox=1,ncolprint)
906
907	write (6,'(a)') 'ppout:'
908	write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
909	enddo ! j
910	endif
911
912	end if
913
914	! ---------------------------------------------------!
915
916	!
917	! ---------------------------------------------------!
918	! DETERMINE CLOUD TOP PRESSURE
919	!
920	! again the 2 methods differ according to whether
921	! or not you use the physical cloud top pressure (top_height = 2)
922	! or the radiatively determined cloud top pressure (top_height = 1 or 3)
923	!
924
925	!compute cloud top pressure
926	do 30 ibox=1,ncol
927	!segregate according to optical thickness
928	if (top_height .eq. 1 .or. top_height .eq. 3) then
929	!find level whose temperature
930	!most closely matches brightness temperature
931	do j=1,npoints
932	nmatch(j)=0
933	enddo
934	do 29 k1=1,nlev-1
935	if (top_height_direction .eq. 2) then
936	ilev = nlev - k1
937	else
938	ilev = k1
939	end if
940	!cdir nodep
941	do j=1,npoints
942	if (ilev .ge. itrop(j)) then
943	if ((at(j,ilev) .ge. tb(j,ibox) .and.
944	& at(j,ilev+1) .le. tb(j,ibox)) .or.
945	& (at(j,ilev) .le. tb(j,ibox) .and.
946	& at(j,ilev+1) .ge. tb(j,ibox))) then
947	nmatch(j)=nmatch(j)+1
948	match(j,nmatch(j))=ilev
949	end if
950	end if
951	enddo
952	29 continue
953
954	do j=1,npoints
955	if (nmatch(j) .ge. 1) then
956	k1 = match(j,nmatch(j))
957	k2 = k1 + 1
958	logp1 = log(pfull(j,k1))
959	logp2 = log(pfull(j,k2))
960	atd = max(tauchk,abs(at(j,k2) - at(j,k1)))
961	logp=logp1+(logp2-logp1)*abs(tb(j,ibox)-at(j,k1))/atd
962	ptop(j,ibox) = exp(logp)
963	if(abs(pfull(j,k1)-ptop(j,ibox)) .lt.
964	& abs(pfull(j,k2)-ptop(j,ibox))) then
965	levmatch(j,ibox)=k1
966	else
967	levmatch(j,ibox)=k2
968	end if
969	else
970	if (tb(j,ibox) .le. attrop(j)) then
971	ptop(j,ibox)=ptrop(j)
972	levmatch(j,ibox)=itrop(j)
973	end if
974	if (tb(j,ibox) .ge. atmax(j)) then
975	ptop(j,ibox)=pfull(j,nlev)
976	levmatch(j,ibox)=nlev
977	end if
978	end if
979	enddo ! j
980
981	else ! if (top_height .eq. 1 .or. top_height .eq. 3)
982
983	do j=1,npoints
984	ptop(j,ibox)=0.
985	enddo
986	do ilev=1,nlev
987	do j=1,npoints
988	if ((ptop(j,ibox) .eq. 0. )
989	& .and.(frac_out(j,ibox,ilev) .ne. 0)) then
990	ptop(j,ibox)=phalf(j,ilev)
991	levmatch(j,ibox)=ilev
992	end if
993	end do
994	end do
995	end if
996
997	do j=1,npoints
998	if (tau(j,ibox) .le. (tauchk )) then
999	ptop(j,ibox)=0.
1000	levmatch(j,ibox)=0
1001	endif
1002	enddo
1003
1004	30 continue
1005
1006	!
1007	!
1008	! ---------------------------------------------------!
1009
1010
1011	!
1012	! ---------------------------------------------------!
1013	! DETERMINE ISCCP CLOUD TYPE FREQUENCIES
1014	!
1015	! Now that ptop and tau have been determined,
1016	! determine amount of each of the 49 ISCCP cloud
1017	! types
1018	!
1019	! Also compute grid box mean cloud top pressure and
1020	! optical thickness. The mean cloud top pressure and
1021	! optical thickness are averages over the cloudy
1022	! area only. The mean cloud top pressure is a linear
1023	! average of the cloud top pressures. The mean cloud
1024	! optical thickness is computed by converting optical
1025	! thickness to an albedo, averaging in albedo units,
1026	! then converting the average albedo back to a mean
1027	! optical thickness.
1028	!
1029
1030	!compute isccp frequencies
1031
1032	!reset frequencies
1033	do 38 ilev=1,7
1034	do 38 ilev2=1,7
1035	do j=1,npoints !
1036	if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1037	fq_isccp(j,ilev,ilev2)= 0.
1038	else
1039	fq_isccp(j,ilev,ilev2)= output_missing_value
1040	end if
1041	enddo
1042	38 continue
1043
1044	!reset variables need for averaging cloud properties
1045	do j=1,npoints
1046	if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1047	totalcldarea(j) = 0.
1048	meanalbedocld(j) = 0.
1049	meanptop(j) = 0.
1050	meantaucld(j) = 0.
1051	else
1052	totalcldarea(j) = output_missing_value
1053	meanalbedocld(j) = output_missing_value
1054	meanptop(j) = output_missing_value
1055	meantaucld(j) = output_missing_value
1056	end if
1057	enddo ! j
1058
1059	boxarea = 1./real(ncol)
1060
1061	do 39 ibox=1,ncol
1062	do j=1,npoints
1063
1064	if (tau(j,ibox) .gt. (tauchk )
1065	& .and. ptop(j,ibox) .gt. 0.) then
1066	box_cloudy(j,ibox)=.true.
1067	endif
1068
1069	if (box_cloudy(j,ibox)) then
1070
1071	if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1072
1073	boxtau(j,ibox) = tau(j,ibox)
1074
1075	if (tau(j,ibox) .ge. isccp_taumin) then
1076	totalcldarea(j) = totalcldarea(j) + boxarea
1077
1078	!convert optical thickness to albedo
1079	albedocld(j,ibox)
1080	& = (tau(j,ibox)0.895)/((tau(j,ibox)0.895)+6.82)
1081
1082	!contribute to averaging
1083	meanalbedocld(j) = meanalbedocld(j)
1084	& +albedocld(j,ibox)*boxarea
1085
1086	end if
1087
1088	endif
1089
1090	endif
1091
1092	if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1093
1094	if (box_cloudy(j,ibox)) then
1095
1096	!convert ptop to millibars
1097	ptop(j,ibox)=ptop(j,ibox) / 100.
1098
1099	!save for output cloud top pressure and optical thickness
1100	boxptop(j,ibox) = ptop(j,ibox)
1101
1102	if (tau(j,ibox) .ge. isccp_taumin) then
1103	meanptop(j) = meanptop(j) + ptop(j,ibox)*boxarea
1104	end if
1105
1106	!reset itau(j), ipres(j)
1107	itau(j) = 0
1108	ipres(j) = 0
1109
1110	!determine optical depth category
1111	if (tau(j,ibox) .lt. isccp_taumin) then
1112	itau(j)=1
1113	else if (tau(j,ibox) .ge. isccp_taumin
1114	&
1115	& .and. tau(j,ibox) .lt. 1.3) then
1116	itau(j)=2
1117	else if (tau(j,ibox) .ge. 1.3
1118	& .and. tau(j,ibox) .lt. 3.6) then
1119	itau(j)=3
1120	else if (tau(j,ibox) .ge. 3.6
1121	& .and. tau(j,ibox) .lt. 9.4) then
1122	itau(j)=4
1123	else if (tau(j,ibox) .ge. 9.4
1124	& .and. tau(j,ibox) .lt. 23.) then
1125	itau(j)=5
1126	else if (tau(j,ibox) .ge. 23.
1127	& .and. tau(j,ibox) .lt. 60.) then
1128	itau(j)=6
1129	else if (tau(j,ibox) .ge. 60.) then
1130	itau(j)=7
1131	end if
1132
1133	!determine cloud top pressure category
1134	if ( ptop(j,ibox) .gt. 0.
1135	& .and.ptop(j,ibox) .lt. 180.) then
1136	ipres(j)=1
1137	else if(ptop(j,ibox) .ge. 180.
1138	& .and.ptop(j,ibox) .lt. 310.) then
1139	ipres(j)=2
1140	else if(ptop(j,ibox) .ge. 310.
1141	& .and.ptop(j,ibox) .lt. 440.) then
1142	ipres(j)=3
1143	else if(ptop(j,ibox) .ge. 440.
1144	& .and.ptop(j,ibox) .lt. 560.) then
1145	ipres(j)=4
1146	else if(ptop(j,ibox) .ge. 560.
1147	& .and.ptop(j,ibox) .lt. 680.) then
1148	ipres(j)=5
1149	else if(ptop(j,ibox) .ge. 680.
1150	& .and.ptop(j,ibox) .lt. 800.) then
1151	ipres(j)=6
1152	else if(ptop(j,ibox) .ge. 800.) then
1153	ipres(j)=7
1154	end if
1155
1156	!update frequencies
1157	if(ipres(j) .gt. 0.and.itau(j) .gt. 0) then
1158	fq_isccp(j,itau(j),ipres(j))=
1159	& fq_isccp(j,itau(j),ipres(j))+ boxarea
1160	end if
1161
1162	end if
1163
1164	end if
1165
1166	enddo ! j
1167	39 continue
1168
1169	!compute mean cloud properties
1170	do j=1,npoints
1171	if (totalcldarea(j) .gt. 0.) then
1172	! code above guarantees that totalcldarea > 0
1173	! only if sunlit .eq. 1 .or. top_height = 3
1174	! and applies only to clouds with tau > isccp_taumin
1175	meanptop(j) = meanptop(j) / totalcldarea(j)
1176	meanalbedocld(j) = meanalbedocld(j) / totalcldarea(j)
1177	meantaucld(j) = (6.82/((1./meanalbedocld(j))-1.))**(1./0.895)
1178	else
1179	! this code is necessary so that in the case that totalcldarea = 0.,
1180	! that these variables, which are in-cloud averages, are set to missing
1181	! note that totalcldarea will be 0. if all the clouds in the grid box have
1182	! tau < isccp_taumin
1183	meanptop(j) = output_missing_value
1184	meanalbedocld(j) = output_missing_value
1185	meantaucld(j) = output_missing_value
1186	end if
1187	enddo ! j
1188	!
1189	! ---------------------------------------------------!
1190
1191	! ---------------------------------------------------!
1192	! OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM
1193	!
1194	if (debugcol.ne.0) then
1195	!
1196	do j=1,npoints,debugcol
1197
1198	!produce character output
1199	do ilev=1,nlev
1200	do ibox=1,ncol
1201	acc(ilev,ibox)=0
1202	enddo
1203	enddo
1204
1205	do ilev=1,nlev
1206	do ibox=1,ncol
1207	acc(ilev,ibox)=frac_out(j,ibox,ilev)*2
1208	if (levmatch(j,ibox) .eq. ilev)
1209	& acc(ilev,ibox)=acc(ilev,ibox)+1
1210	enddo
1211	enddo
1212
1213	!print test
1214
1215	write(ftn09,11) j
1216	11 format('ftn09.',i4.4)
1217	open(9, FILE=ftn09, FORM='FORMATTED')
1218
1219	write(9,'(a1)') ' '
1220	write(9,'(10i5)')
1221	& (ilev,ilev=5,nlev,5)
1222	write(9,'(a1)') ' '
1223
1224	do ibox=1,ncol
1225	write(9,'(40(a1),1x,40(a1))')
1226	& (cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev)
1227	& ,(cchar(acc(ilev,ibox)+1),ilev=1,nlev)
1228	end do
1229	close(9)
1230
1231	if (ncolprint.ne.0) then
1232	write(6,'(a1)') ' '
1233	write(6,'(a2,1X,5(a7,1X),a50)')
1234	& 'ilev',
1235	& 'pfull','at',
1236	& 'cc*100','dem_s','dtau_s',
1237	& 'cchar'
1238
1239	! do 4012 ilev=1,nlev
1240	! write(6,'(60i2)') (box(i,ilev),i=1,ncolprint)
1241	! write(6,'(i2,1X,5(f7.2,1X),50(a1))')
1242	! & ilev,
1243	! & pfull(j,ilev)/100.,at(j,ilev),
1244	! & cc(j,ilev)*100.0,dem_s(j,ilev),dtau_s(j,ilev)
1245	! & ,(cchar(acc(ilev,ibox)+1),ibox=1,ncolprint)
1246	!4012 continue
1247	write (6,'(a)') 'skt(j):'
1248	write (6,'(8f7.2)') skt(j)
1249
1250	write (6,'(8I7)') (ibox,ibox=1,ncolprint)
1251
1252	write (6,'(a)') 'tau:'
1253	write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
1254
1255	write (6,'(a)') 'tb:'
1256	write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
1257
1258	write (6,'(a)') 'ptop:'
1259	write (6,'(8f7.2)') (ptop(j,ibox),ibox=1,ncolprint)
1260	endif
1261
1262	enddo
1263
1264	end if
1265
1266	return
1267	end
1268
1269

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

Download in other formats:

Original Format