Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

icarus.F @ 3331

Last change on this file since 3331 was 3331, checked in by acozic, 6 years ago
Add modification for isotopes
Property svn:executable set to ``*
File size: 42.9 KB

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

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

Download in other formats:

Original Format