Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

icarus.F

Last change on this file was 5185, checked in by abarral, 2 months ago
Replace REPROBUS CPP KEY by logical using handmade wonky wrapper
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 42.7 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	! Working variables added when program updated to mimic Mark Webb's PV-Wave code
227	! ------
228
229	REAL dem(npoints,ncol),bb(npoints) ! working variables for 10.5 micron longwave
230	! emissivity in part of
231	! gridbox under consideration
232
233	REAL ptrop(npoints)
234	REAL attrop(npoints)
235	REAL attropmin (npoints)
236	REAL atmax(npoints)
237	REAL btcmin(npoints)
238	REAL transmax(npoints)
239
240	INTEGER i,j,ilev,ibox,itrop(npoints)
241	INTEGER ipres(npoints)
242	INTEGER itau(npoints),ilev2
243	INTEGER acc(nlev,ncol)
244	INTEGER match(npoints,nlev-1)
245	INTEGER nmatch(npoints)
246	INTEGER levmatch(npoints,ncol)
247
248	!variables needed for water vapor continuum absorption
249	real fluxtop_clrsky(npoints),trans_layers_above_clrsky(npoints)
250	real taumin(npoints)
251	real dem_wv(npoints,nlev), wtmair, wtmh20, Navo, grav, pstd, t0
252	real press(npoints), dpress(npoints), atmden(npoints)
253	real rvh20(npoints), wk(npoints), rhoave(npoints)
254	real rh20s(npoints), rfrgn(npoints)
255	real tmpexp(npoints),tauwv(npoints)
256
257	character*1 cchar(6),cchar_realtops(6)
258	integer icycle
259	REAL tau(npoints,ncol)
260	LOGICAL box_cloudy(npoints,ncol)
261	REAL tb(npoints,ncol)
262	REAL ptop(npoints,ncol)
263	REAL emcld(npoints,ncol)
264	REAL fluxtop(npoints,ncol)
265	REAL trans_layers_above(npoints,ncol)
266	real isccp_taumin,fluxtopinit(npoints),tauir(npoints)
267	REAL albedocld(npoints,ncol)
268	real boxarea
269	integer debug ! set to non-zero value to print out inputs
270	! with step debug
271	integer debugcol ! set to non-zero value to print out column
272	! decomposition with step debugcol
273	integer rangevec(npoints),rangeerror
274
275	integer index1(npoints),num1,jj,k1,k2
276	real rec2p13,tauchk,logp,logp1,logp2,atd
277	real output_missing_value
278
279	character*10 ftn09
280
281	DATA isccp_taumin / 0.3 /
282	DATA output_missing_value / -1.E+30 /
283	DATA cchar / ' ','-','1','+','I','+'/
284	DATA cchar_realtops / ' ',' ','1','1','I','I'/
285
286	! ------ End duplicate definitions common to wrapper routine
287
288	tauchk = -1.*log(0.9999999)
289	rec2p13=1./2.13
290
291	ncolprint=0
292
293	if ( debug.ne.0 ) then
294	j=1
295	write(6,'(a10)') 'j='
296	write(6,'(8I10)') j
297	write(6,'(a10)') 'debug='
298	write(6,'(8I10)') debug
299	write(6,'(a10)') 'debugcol='
300	write(6,'(8I10)') debugcol
301	write(6,'(a10)') 'npoints='
302	write(6,'(8I10)') npoints
303	write(6,'(a10)') 'nlev='
304	write(6,'(8I10)') nlev
305	write(6,'(a10)') 'ncol='
306	write(6,'(8I10)') ncol
307	write(6,'(a11)') 'top_height='
308	write(6,'(8I10)') top_height
309	write(6,'(a21)') 'top_height_direction='
310	write(6,'(8I10)') top_height_direction
311	write(6,'(a10)') 'overlap='
312	write(6,'(8I10)') overlap
313	write(6,'(a10)') 'emsfc_lw='
314	write(6,'(8f10.2)') emsfc_lw
315	do j=1,npoints,debug
316	write(6,'(a10)') 'j='
317	write(6,'(8I10)') j
318	write(6,'(a10)') 'sunlit='
319	write(6,'(8I10)') sunlit(j)
320	write(6,'(a10)') 'pfull='
321	write(6,'(8f10.2)') (pfull(j,i),i=1,nlev)
322	write(6,'(a10)') 'phalf='
323	write(6,'(8f10.2)') (phalf(j,i),i=1,nlev+1)
324	write(6,'(a10)') 'qv='
325	write(6,'(8f10.3)') (qv(j,i),i=1,nlev)
326	write(6,'(a10)') 'cc='
327	write(6,'(8f10.3)') (cc(j,i),i=1,nlev)
328	write(6,'(a10)') 'conv='
329	write(6,'(8f10.2)') (conv(j,i),i=1,nlev)
330	write(6,'(a10)') 'dtau_s='
331	write(6,'(8g12.5)') (dtau_s(j,i),i=1,nlev)
332	write(6,'(a10)') 'dtau_c='
333	write(6,'(8f10.2)') (dtau_c(j,i),i=1,nlev)
334	write(6,'(a10)') 'skt='
335	write(6,'(8f10.2)') skt(j)
336	write(6,'(a10)') 'at='
337	write(6,'(8f10.2)') (at(j,i),i=1,nlev)
338	write(6,'(a10)') 'dem_s='
339	write(6,'(8f10.3)') (dem_s(j,i),i=1,nlev)
340	write(6,'(a10)') 'dem_c='
341	write(6,'(8f10.3)') (dem_c(j,i),i=1,nlev)
342	enddo
343	endif
344
345	! ---------------------------------------------------!
346
347	if (ncolprint.ne.0) then
348	do j=1,npoints,1000
349	write(6,'(a10)') 'j='
350	write(6,'(8I10)') j
351	enddo
352	endif
353
354	if (top_height .eq. 1 .or. top_height .eq. 3) then
355
356	do j=1,npoints
357	ptrop(j)=5000.
358	attropmin(j) = 400.
359	atmax(j) = 0.
360	attrop(j) = 120.
361	itrop(j) = 1
362	enddo
363
364	do 12 ilev=1,nlev
365	do j=1,npoints
366	if (pfull(j,ilev) .lt. 40000. .AND.
367	& pfull(j,ilev) .gt. 5000. .AND.
368	& at(j,ilev) .lt. attropmin(j)) then
369	ptrop(j) = pfull(j,ilev)
370	attropmin(j) = at(j,ilev)
371	attrop(j) = attropmin(j)
372	itrop(j)=ilev
373	end if
374	enddo
375	12 continue
376
377	do 13 ilev=1,nlev
378	do j=1,npoints
379	if (at(j,ilev) .gt. atmax(j) .AND.
380	& ilev .ge. itrop(j)) atmax(j)=at(j,ilev)
381	enddo
382	13 continue
383
384	end if
385
386
387	if (top_height .eq. 1 .or. top_height .eq. 3) then
388	do j=1,npoints
389	meantb(j) = 0.
390	meantbclr(j) = 0.
391	end do
392	else
393	do j=1,npoints
394	meantb(j) = output_missing_value
395	meantbclr(j) = output_missing_value
396	end do
397	end if
398
399	! -----------------------------------------------------!
400
401	! ---------------------------------------------------!
402
403	do ilev=1,nlev
404	do j=1,npoints
405
406	rangevec(j)=0
407
408	if (cc(j,ilev) .lt. 0. .or. cc(j,ilev) .gt. 1.) then
409	! error = cloud fraction less than zero
410	! error = cloud fraction greater than 1
411	rangevec(j)=rangevec(j)+1
412	endif
413
414	if (conv(j,ilev) .lt. 0. .or. conv(j,ilev) .gt. 1.) then
415	! ' error = convective cloud fraction less than zero'
416	! ' error = convective cloud fraction greater than 1'
417	rangevec(j)=rangevec(j)+2
418	endif
419
420	if (dtau_s(j,ilev) .lt. 0.) then
421	! ' error = stratiform cloud opt. depth less than zero'
422	rangevec(j)=rangevec(j)+4
423	endif
424
425	if (dtau_c(j,ilev) .lt. 0.) then
426	! ' error = convective cloud opt. depth less than zero'
427	rangevec(j)=rangevec(j)+8
428	endif
429
430	if (dem_s(j,ilev) .lt. 0. .or. dem_s(j,ilev) .gt. 1.) then
431	! ' error = stratiform cloud emissivity less than zero'
432	! ' error = stratiform cloud emissivity greater than 1'
433	rangevec(j)=rangevec(j)+16
434	endif
435
436	if (dem_c(j,ilev) .lt. 0. .or. dem_c(j,ilev) .gt. 1.) then
437	! ' error = convective cloud emissivity less than zero'
438	! ' error = convective cloud emissivity greater than 1'
439	rangevec(j)=rangevec(j)+32
440	endif
441	enddo
442
443	rangeerror=0
444	do j=1,npoints
445	rangeerror=rangeerror+rangevec(j)
446	enddo
447
448	if (rangeerror.ne.0) then
449	write (6,*) 'Input variable out of range'
450	write (6,*) 'rangevec:'
451	write (6,*) rangevec
452	STOP
453	endif
454	enddo
455
456	! ---------------------------------------------------!
457
458	! ---------------------------------------------------!
459	! 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	! COMPUTE INFRARED BRIGHTNESS TEMPERUATRES
504	! AND CLOUD TOP TEMPERATURE SATELLITE SHOULD SEE
505
506	! again this is only done if top_height = 1 or 3
507
508	! fluxtop is the 10.5 micron radiance at the top of the
509	! atmosphere
510	! trans_layers_above is the total transmissivity in the layers
511	! above the current layer
512	! fluxtop_clrsky(j) and trans_layers_above_clrsky(j) are the clear
513	! sky versions of these quantities.
514
515	if (top_height .eq. 1 .or. top_height .eq. 3) then
516
517
518	!----------------------------------------------------------------------
519
520	! DO CLEAR SKY RADIANCE CALCULATION FIRST
521
522	!compute water vapor continuum emissivity
523	!this treatment follows Schwarkzopf and Ramasamy
524	!JGR 1999,vol 104, pages 9467-9499.
525	!the emissivity is calculated at a wavenumber of 955 cm-1,
526	!or 10.47 microns
527	wtmair = 28.9644
528	wtmh20 = 18.01534
529	Navo = 6.023E+23
530	grav = 9.806650E+02
531	pstd = 1.013250E+06
532	t0 = 296.
533	if (ncolprint .ne. 0)
534	& write(6,*) 'ilev pw (kg/m2) tauwv(j) dem_wv'
535	do 125 ilev=1,nlev
536	do j=1,npoints
537	!press and dpress are dyne/cm2 = Pascals *10
538	press(j) = pfull(j,ilev)*10.
539	dpress(j) = (phalf(j,ilev+1)-phalf(j,ilev))*10
540	!atmden = g/cm2 = kg/m2 / 10
541	atmden(j) = dpress(j)/grav
542	rvh20(j) = qv(j,ilev)*wtmair/wtmh20
543	wk(j) = rvh20(j)Navoatmden(j)/wtmair
544	rhoave(j) = (press(j)/pstd)*(t0/at(j,ilev))
545	rh20s(j) = rvh20(j)*rhoave(j)
546	rfrgn(j) = rhoave(j)-rh20s(j)
547	tmpexp(j) = exp(-0.02*(at(j,ilev)-t0))
548	tauwv(j) = wk(j)1.e-20(
549	& (0.0224697rh20s(j)tmpexp(j)) +
550	& (3.41817e-7rfrgn(j)) )0.98
551	dem_wv(j,ilev) = 1. - exp( -1. * tauwv(j))
552	enddo
553	if (ncolprint .ne. 0) then
554	do j=1,npoints ,1000
555	write(6,'(a10)') 'j='
556	write(6,'(8I10)') j
557	write(6,'(i2,1X,3(f8.3,3X))') ilev,
558	& qv(j,ilev)*(phalf(j,ilev+1)-phalf(j,ilev))/(grav/100.),
559	& tauwv(j),dem_wv(j,ilev)
560	enddo
561	endif
562	125 continue
563
564	!initialize variables
565	do j=1,npoints
566	fluxtop_clrsky(j) = 0.
567	trans_layers_above_clrsky(j)=1.
568	enddo
569
570	do ilev=1,nlev
571	do j=1,npoints
572
573	! Black body emission at temperature of the layer
574
575	bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. )
576	!bb(j)= 5.67e-8at(j,ilev)*4
577
578	! increase TOA flux by flux emitted from layer
579	! times total transmittance in layers above
580
581	fluxtop_clrsky(j) = fluxtop_clrsky(j)
582	& + dem_wv(j,ilev)bb(j)trans_layers_above_clrsky(j)
583
584	! update trans_layers_above with transmissivity
585	! from this layer for next time around loop
586
587	trans_layers_above_clrsky(j)=
588	& trans_layers_above_clrsky(j)*(1.-dem_wv(j,ilev))
589
590
591	enddo
592	if (ncolprint.ne.0) then
593	do j=1,npoints ,1000
594	write(6,'(a10)') 'j='
595	write(6,'(8I10)') j
596	write (6,'(a)') 'ilev:'
597	write (6,'(I2)') ilev
598
599	write (6,'(a)')
600	& 'emiss_layer,100.bb(j),100.f,total_trans:'
601	write (6,'(4(f7.2,1X))') dem_wv(j,ilev),100.*bb(j),
602	& 100.*fluxtop_clrsky(j),trans_layers_above_clrsky(j)
603	enddo
604	endif
605
606	enddo !loop over level
607
608	do j=1,npoints
609	!add in surface emission
610	bb(j)=1/( exp(1307.27/skt(j)) - 1. )
611	!bb(j)=5.67e-8skt(j)*4
612
613	fluxtop_clrsky(j) = fluxtop_clrsky(j) + emsfc_lw * bb(j)
614	& * trans_layers_above_clrsky(j)
615
616	!clear sky brightness temperature
617	meantbclr(j) = 1307.27/(log(1.+(1./fluxtop_clrsky(j))))
618
619	enddo
620
621	if (ncolprint.ne.0) then
622	do j=1,npoints ,1000
623	write(6,'(a10)') 'j='
624	write(6,'(8I10)') j
625	write (6,'(a)') 'id:'
626	write (6,'(a)') 'surface'
627
628	write (6,'(a)') 'emsfc,100.bb(j),100.f,total_trans:'
629	write (6,'(5(f7.2,1X))') emsfc_lw,100.*bb(j),
630	& 100.*fluxtop_clrsky(j),
631	& trans_layers_above_clrsky(j), meantbclr(j)
632	enddo
633	endif
634
635	! END OF CLEAR SKY CALCULATION
636
637	!----------------------------------------------------------------
638
639
640
641	if (ncolprint.ne.0) then
642
643	do j=1,npoints ,1000
644	write(6,'(a10)') 'j='
645	write(6,'(8I10)') j
646	write (6,'(a)') 'ts:'
647	write (6,'(8f7.2)') (skt(j),ibox=1,ncolprint)
648
649	write (6,'(a)') 'ta_rev:'
650	write (6,'(8f7.2)')
651	& ((at(j,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
652
653	enddo
654	endif
655	!loop over columns
656	do ibox=1,ncol
657	do j=1,npoints
658	fluxtop(j,ibox)=0.
659	trans_layers_above(j,ibox)=1.
660	enddo
661	enddo
662
663	do ilev=1,nlev
664	do j=1,npoints
665	! Black body emission at temperature of the layer
666
667	bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. )
668	!bb(j)= 5.67e-8at(j,ilev)*4
669	enddo
670
671	do ibox=1,ncol
672	do j=1,npoints
673
674	! emissivity for point in this layer
675	if (frac_out(j,ibox,ilev).eq.1) then
676	dem(j,ibox)= 1. -
677	& ( (1. - dem_wv(j,ilev)) * (1. - dem_s(j,ilev)) )
678	else if (frac_out(j,ibox,ilev).eq.2) then
679	dem(j,ibox)= 1. -
680	& ( (1. - dem_wv(j,ilev)) * (1. - dem_c(j,ilev)) )
681	else
682	dem(j,ibox)= dem_wv(j,ilev)
683	end if
684
685
686	! increase TOA flux by flux emitted from layer
687	! times total transmittance in layers above
688
689	fluxtop(j,ibox) = fluxtop(j,ibox)
690	& + dem(j,ibox) * bb(j)
691	& * trans_layers_above(j,ibox)
692
693	! update trans_layers_above with transmissivity
694	! from this layer for next time around loop
695
696	trans_layers_above(j,ibox)=
697	& trans_layers_above(j,ibox)*(1.-dem(j,ibox))
698
699	enddo ! j
700	enddo ! ibox
701
702	if (ncolprint.ne.0) then
703	do j=1,npoints,1000
704	write (6,'(a)') 'ilev:'
705	write (6,'(I2)') ilev
706
707	write(6,'(a10)') 'j='
708	write(6,'(8I10)') j
709	write (6,'(a)') 'emiss_layer:'
710	write (6,'(8f7.2)') (dem(j,ibox),ibox=1,ncolprint)
711
712	write (6,'(a)') '100.*bb(j):'
713	write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
714
715	write (6,'(a)') '100.*f:'
716	write (6,'(8f7.2)')
717	& (100.*fluxtop(j,ibox),ibox=1,ncolprint)
718
719	write (6,'(a)') 'total_trans:'
720	write (6,'(8f7.2)')
721	& (trans_layers_above(j,ibox),ibox=1,ncolprint)
722	enddo
723	endif
724
725	enddo ! ilev
726
727
728	do j=1,npoints
729	!add in surface emission
730	bb(j)=1/( exp(1307.27/skt(j)) - 1. )
731	!bb(j)=5.67e-8skt(j)*4
732	end do
733
734	do ibox=1,ncol
735	do j=1,npoints
736
737	!add in surface emission
738
739	fluxtop(j,ibox) = fluxtop(j,ibox)
740	& + emsfc_lw * bb(j)
741	& * trans_layers_above(j,ibox)
742
743	end do
744	end do
745
746	!calculate mean infrared brightness temperature
747	do ibox=1,ncol
748	do j=1,npoints
749	meantb(j) = meantb(j)+1307.27/(log(1.+(1./fluxtop(j,ibox))))
750	end do
751	end do
752	do j=1, npoints
753	meantb(j) = meantb(j) / real(ncol)
754	end do
755
756	if (ncolprint.ne.0) then
757
758	do j=1,npoints ,1000
759	write(6,'(a10)') 'j='
760	write(6,'(8I10)') j
761	write (6,'(a)') 'id:'
762	write (6,'(a)') 'surface'
763
764	write (6,'(a)') 'emiss_layer:'
765	write (6,'(8f7.2)') (dem(1,ibox),ibox=1,ncolprint)
766
767	write (6,'(a)') '100.*bb(j):'
768	write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
769
770	write (6,'(a)') '100.*f:'
771	write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint)
772
773	write (6,'(a)') 'meantb(j):'
774	write (6,'(8f7.2)') (meantb(j),ibox=1,ncolprint)
775
776	end do
777	endif
778
779	!now that you have the top of atmosphere radiance account
780	!for ISCCP procedures to determine cloud top temperature
781
782	!account for partially transmitting cloud recompute flux
783	!ISCCP would see assuming a single layer cloud
784	!note choice here of 2.13, as it is primarily ice
785	!clouds which have partial emissivity and need the
786	!adjustment performed in this section
787
788	!If it turns out that the cloud brightness temperature
789	!is greater than 260K, then the liquid cloud conversion
790	!factor of 2.56 is used.
791
792	!Note that this is discussed on pages 85-87 of
793	!the ISCCP D level documentation (Rossow et al. 1996)
794
795	do j=1,npoints
796	!compute minimum brightness temperature and optical depth
797	btcmin(j) = 1. / ( exp(1307.27/(attrop(j)-5.)) - 1. )
798	enddo
799	do ibox=1,ncol
800	do j=1,npoints
801	transmax(j) = (fluxtop(j,ibox)-btcmin(j))
802	& /(fluxtop_clrsky(j)-btcmin(j))
803	!note that the initial setting of tauir(j) is needed so that
804	!tauir(j) has a realistic value should the next if block be
805	!bypassed
806	tauir(j) = tau(j,ibox) * rec2p13
807	taumin(j) = -1. * log(max(min(transmax(j),0.9999999),0.001))
808
809	enddo
810
811	if (top_height .eq. 1) then
812	do j=1,npoints
813	if (transmax(j) .gt. 0.001 .AND.
814	& transmax(j) .le. 0.9999999) then
815	fluxtopinit(j) = fluxtop(j,ibox)
816	tauir(j) = tau(j,ibox) *rec2p13
817	endif
818	enddo
819	do icycle=1,2
820	do j=1,npoints
821	if (tau(j,ibox) .gt. (tauchk )) then
822	if (transmax(j) .gt. 0.001 .AND.
823	& transmax(j) .le. 0.9999999) then
824	emcld(j,ibox) = 1. - exp(-1. * tauir(j) )
825	fluxtop(j,ibox) = fluxtopinit(j) -
826	& ((1.-emcld(j,ibox))*fluxtop_clrsky(j))
827	fluxtop(j,ibox)=max(1.E-06,
828	& (fluxtop(j,ibox)/emcld(j,ibox)))
829	tb(j,ibox)= 1307.27
830	& / (log(1. + (1./fluxtop(j,ibox))))
831	if (tb(j,ibox) .gt. 260.) then
832	tauir(j) = tau(j,ibox) / 2.56
833	end if
834	end if
835	end if
836	enddo
837	enddo
838
839	endif
840
841	do j=1,npoints
842	if (tau(j,ibox) .gt. (tauchk )) then
843	!cloudy box
844	!NOTE: tb is the cloud-top temperature not infrared brightness temperature
845	!at this point in the code
846	tb(j,ibox)= 1307.27/ (log(1. + (1./fluxtop(j,ibox))))
847	if (top_height.eq.1.AND.tauir(j).lt.taumin(j)) then
848	tb(j,ibox) = attrop(j) - 5.
849	tau(j,ibox) = 2.13*taumin(j)
850	end if
851	else
852	!clear sky brightness temperature
853	tb(j,ibox) = meantbclr(j)
854	end if
855	enddo ! j
856	enddo ! ibox
857
858	if (ncolprint.ne.0) then
859
860	do j=1,npoints,1000
861	write(6,'(a10)') 'j='
862	write(6,'(8I10)') j
863
864	write (6,'(a)') 'attrop:'
865	write (6,'(8f7.2)') (attrop(j))
866
867	write (6,'(a)') 'btcmin:'
868	write (6,'(8f7.2)') (btcmin(j))
869
870	write (6,'(a)') 'fluxtop_clrsky*100:'
871	write (6,'(8f7.2)')
872	& (100.*fluxtop_clrsky(j))
873
874	write (6,'(a)') '100.*f_adj:'
875	write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint)
876
877	write (6,'(a)') 'transmax:'
878	write (6,'(8f7.2)') (transmax(ibox),ibox=1,ncolprint)
879
880	write (6,'(a)') 'tau:'
881	write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
882
883	write (6,'(a)') 'emcld:'
884	write (6,'(8f7.2)') (emcld(j,ibox),ibox=1,ncolprint)
885
886	write (6,'(a)') 'total_trans:'
887	write (6,'(8f7.2)')
888	& (trans_layers_above(j,ibox),ibox=1,ncolprint)
889
890	write (6,'(a)') 'total_emiss:'
891	write (6,'(8f7.2)')
892	& (1.0-trans_layers_above(j,ibox),ibox=1,ncolprint)
893
894	write (6,'(a)') 'total_trans:'
895	write (6,'(8f7.2)')
896	& (trans_layers_above(j,ibox),ibox=1,ncolprint)
897
898	write (6,'(a)') 'ppout:'
899	write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
900	enddo ! j
901	endif
902
903	end if
904
905	! ---------------------------------------------------!
906
907	! ---------------------------------------------------!
908	! DETERMINE CLOUD TOP PRESSURE
909
910	! again the 2 methods differ according to whether
911	! or not you use the physical cloud top pressure (top_height = 2)
912	! or the radiatively determined cloud top pressure (top_height = 1 or 3)
913
914	!compute cloud top pressure
915	do 30 ibox=1,ncol
916	!segregate according to optical thickness
917	if (top_height .eq. 1 .or. top_height .eq. 3) then
918	!find level whose temperature
919	!most closely matches brightness temperature
920	do j=1,npoints
921	nmatch(j)=0
922	enddo
923	do 29 k1=1,nlev-1
924	if (top_height_direction .eq. 2) then
925	ilev = nlev - k1
926	else
927	ilev = k1
928	end if
929	!cdir nodep
930	do j=1,npoints
931	if (ilev .ge. itrop(j)) then
932	if ((at(j,ilev) .ge. tb(j,ibox) .AND.
933	& at(j,ilev+1) .le. tb(j,ibox)) .or.
934	& (at(j,ilev) .le. tb(j,ibox) .AND.
935	& at(j,ilev+1) .ge. tb(j,ibox))) then
936	nmatch(j)=nmatch(j)+1
937	match(j,nmatch(j))=ilev
938	end if
939	end if
940	enddo
941	29 continue
942
943	do j=1,npoints
944	if (nmatch(j) .ge. 1) then
945	k1 = match(j,nmatch(j))
946	k2 = k1 + 1
947	logp1 = log(pfull(j,k1))
948	logp2 = log(pfull(j,k2))
949	atd = max(tauchk,abs(at(j,k2) - at(j,k1)))
950	logp=logp1+(logp2-logp1)*abs(tb(j,ibox)-at(j,k1))/atd
951	ptop(j,ibox) = exp(logp)
952	if(abs(pfull(j,k1)-ptop(j,ibox)) .lt.
953	& abs(pfull(j,k2)-ptop(j,ibox))) then
954	levmatch(j,ibox)=k1
955	else
956	levmatch(j,ibox)=k2
957	end if
958	else
959	if (tb(j,ibox) .le. attrop(j)) then
960	ptop(j,ibox)=ptrop(j)
961	levmatch(j,ibox)=itrop(j)
962	end if
963	if (tb(j,ibox) .ge. atmax(j)) then
964	ptop(j,ibox)=pfull(j,nlev)
965	levmatch(j,ibox)=nlev
966	end if
967	end if
968	enddo ! j
969
970	else ! if (top_height .eq. 1 .or. top_height .eq. 3)
971
972	do j=1,npoints
973	ptop(j,ibox)=0.
974	enddo
975	do ilev=1,nlev
976	do j=1,npoints
977	if ((ptop(j,ibox) .eq. 0. )
978	& .AND.(frac_out(j,ibox,ilev) .ne. 0)) then
979	ptop(j,ibox)=phalf(j,ilev)
980	levmatch(j,ibox)=ilev
981	end if
982	end do
983	end do
984	end if
985
986	do j=1,npoints
987	if (tau(j,ibox) .le. (tauchk )) then
988	ptop(j,ibox)=0.
989	levmatch(j,ibox)=0
990	endif
991	enddo
992
993	30 continue
994
995
996	! ---------------------------------------------------!
997
998	! ---------------------------------------------------!
999	! DETERMINE ISCCP CLOUD TYPE FREQUENCIES
1000
1001	! Now that ptop and tau have been determined,
1002	! determine amount of each of the 49 ISCCP cloud
1003	! types
1004
1005	! Also compute grid box mean cloud top pressure and
1006	! optical thickness. The mean cloud top pressure and
1007	! optical thickness are averages over the cloudy
1008	! area only. The mean cloud top pressure is a linear
1009	! average of the cloud top pressures. The mean cloud
1010	! optical thickness is computed by converting optical
1011	! thickness to an albedo, averaging in albedo units,
1012	! then converting the average albedo back to a mean
1013	! optical thickness.
1014
1015	!compute isccp frequencies
1016
1017	!reset frequencies
1018	do 38 ilev=1,7
1019	do 38 ilev2=1,7
1020	do j=1,npoints !
1021	if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1022	fq_isccp(j,ilev,ilev2)= 0.
1023	else
1024	fq_isccp(j,ilev,ilev2)= output_missing_value
1025	end if
1026	enddo
1027	38 continue
1028
1029	!reset variables need for averaging cloud properties
1030	do j=1,npoints
1031	if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1032	totalcldarea(j) = 0.
1033	meanalbedocld(j) = 0.
1034	meanptop(j) = 0.
1035	meantaucld(j) = 0.
1036	else
1037	totalcldarea(j) = output_missing_value
1038	meanalbedocld(j) = output_missing_value
1039	meanptop(j) = output_missing_value
1040	meantaucld(j) = output_missing_value
1041	end if
1042	enddo ! j
1043
1044	boxarea = 1./real(ncol)
1045
1046	do 39 ibox=1,ncol
1047	do j=1,npoints
1048
1049	if (tau(j,ibox) .gt. (tauchk )
1050	& .AND. ptop(j,ibox) .gt. 0.) then
1051	box_cloudy(j,ibox)=.true.
1052	endif
1053
1054	if (box_cloudy(j,ibox)) then
1055
1056	if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1057
1058	boxtau(j,ibox) = tau(j,ibox)
1059
1060	if (tau(j,ibox) .ge. isccp_taumin) then
1061	totalcldarea(j) = totalcldarea(j) + boxarea
1062
1063	!convert optical thickness to albedo
1064	albedocld(j,ibox)
1065	& = (tau(j,ibox)0.895)/((tau(j,ibox)0.895)+6.82)
1066
1067	!contribute to averaging
1068	meanalbedocld(j) = meanalbedocld(j)
1069	& +albedocld(j,ibox)*boxarea
1070
1071	end if
1072
1073	endif
1074
1075	endif
1076
1077	if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1078
1079	if (box_cloudy(j,ibox)) then
1080
1081	!convert ptop to millibars
1082	ptop(j,ibox)=ptop(j,ibox) / 100.
1083
1084	!save for output cloud top pressure and optical thickness
1085	boxptop(j,ibox) = ptop(j,ibox)
1086
1087	if (tau(j,ibox) .ge. isccp_taumin) then
1088	meanptop(j) = meanptop(j) + ptop(j,ibox)*boxarea
1089	end if
1090
1091	!reset itau(j), ipres(j)
1092	itau(j) = 0
1093	ipres(j) = 0
1094
1095	!determine optical depth category
1096	if (tau(j,ibox) .lt. isccp_taumin) then
1097	itau(j)=1
1098	else if (tau(j,ibox) .ge. isccp_taumin
1099	&
1100	& .AND. tau(j,ibox) .lt. 1.3) then
1101	itau(j)=2
1102	else if (tau(j,ibox) .ge. 1.3
1103	& .AND. tau(j,ibox) .lt. 3.6) then
1104	itau(j)=3
1105	else if (tau(j,ibox) .ge. 3.6
1106	& .AND. tau(j,ibox) .lt. 9.4) then
1107	itau(j)=4
1108	else if (tau(j,ibox) .ge. 9.4
1109	& .AND. tau(j,ibox) .lt. 23.) then
1110	itau(j)=5
1111	else if (tau(j,ibox) .ge. 23.
1112	& .AND. tau(j,ibox) .lt. 60.) then
1113	itau(j)=6
1114	else if (tau(j,ibox) .ge. 60.) then
1115	itau(j)=7
1116	end if
1117
1118	!determine cloud top pressure category
1119	if ( ptop(j,ibox) .gt. 0.
1120	& .AND.ptop(j,ibox) .lt. 180.) then
1121	ipres(j)=1
1122	else if(ptop(j,ibox) .ge. 180.
1123	& .AND.ptop(j,ibox) .lt. 310.) then
1124	ipres(j)=2
1125	else if(ptop(j,ibox) .ge. 310.
1126	& .AND.ptop(j,ibox) .lt. 440.) then
1127	ipres(j)=3
1128	else if(ptop(j,ibox) .ge. 440.
1129	& .AND.ptop(j,ibox) .lt. 560.) then
1130	ipres(j)=4
1131	else if(ptop(j,ibox) .ge. 560.
1132	& .AND.ptop(j,ibox) .lt. 680.) then
1133	ipres(j)=5
1134	else if(ptop(j,ibox) .ge. 680.
1135	& .AND.ptop(j,ibox) .lt. 800.) then
1136	ipres(j)=6
1137	else if(ptop(j,ibox) .ge. 800.) then
1138	ipres(j)=7
1139	end if
1140
1141	!update frequencies
1142	if(ipres(j) .gt. 0.AND.itau(j) .gt. 0) then
1143	fq_isccp(j,itau(j),ipres(j))=
1144	& fq_isccp(j,itau(j),ipres(j))+ boxarea
1145	end if
1146
1147	end if
1148
1149	end if
1150
1151	enddo ! j
1152	39 continue
1153
1154	!compute mean cloud properties
1155	do j=1,npoints
1156	if (totalcldarea(j) .gt. 0.) then
1157	! code above guarantees that totalcldarea > 0
1158	! only if sunlit .eq. 1 .or. top_height = 3
1159	! and applies only to clouds with tau > isccp_taumin
1160	meanptop(j) = meanptop(j) / totalcldarea(j)
1161	meanalbedocld(j) = meanalbedocld(j) / totalcldarea(j)
1162	meantaucld(j) = (6.82/((1./meanalbedocld(j))-1.))**(1./0.895)
1163	else
1164	! this code is necessary so that in the case that totalcldarea = 0.,
1165	! that these variables, which are in-cloud averages, are set to missing
1166	! note that totalcldarea will be 0. if all the clouds in the grid box have
1167	! tau < isccp_taumin
1168	meanptop(j) = output_missing_value
1169	meanalbedocld(j) = output_missing_value
1170	meantaucld(j) = output_missing_value
1171	end if
1172	enddo ! j
1173
1174	! ---------------------------------------------------!
1175
1176	! ---------------------------------------------------!
1177	! OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM
1178
1179	if (debugcol.ne.0) then
1180
1181	do j=1,npoints,debugcol
1182
1183	!produce character output
1184	do ilev=1,nlev
1185	do ibox=1,ncol
1186	acc(ilev,ibox)=0
1187	enddo
1188	enddo
1189
1190	do ilev=1,nlev
1191	do ibox=1,ncol
1192	acc(ilev,ibox)=frac_out(j,ibox,ilev)*2
1193	if (levmatch(j,ibox) .eq. ilev)
1194	& acc(ilev,ibox)=acc(ilev,ibox)+1
1195	enddo
1196	enddo
1197
1198	!print test
1199
1200	write(ftn09,11) j
1201	11 format('ftn09.',i4.4)
1202	open(9, FILE=ftn09, FORM='FORMATTED')
1203
1204	write(9,'(a1)') ' '
1205	write(9,'(10i5)')
1206	& (ilev,ilev=5,nlev,5)
1207	write(9,'(a1)') ' '
1208
1209	do ibox=1,ncol
1210	write(9,'(40(a1),1x,40(a1))')
1211	& (cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev)
1212	& ,(cchar(acc(ilev,ibox)+1),ilev=1,nlev)
1213	end do
1214	close(9)
1215
1216	if (ncolprint.ne.0) then
1217	write(6,'(a1)') ' '
1218	write(6,'(a2,1X,5(a7,1X),a50)')
1219	& 'ilev',
1220	& 'pfull','at',
1221	& 'cc*100','dem_s','dtau_s',
1222	& 'cchar'
1223
1224	! do 4012 ilev=1,nlev
1225	! write(6,'(60i2)') (box(i,ilev),i=1,ncolprint)
1226	! write(6,'(i2,1X,5(f7.2,1X),50(a1))')
1227	! & ilev,
1228	! & pfull(j,ilev)/100.,at(j,ilev),
1229	! & cc(j,ilev)*100.0,dem_s(j,ilev),dtau_s(j,ilev)
1230	! & ,(cchar(acc(ilev,ibox)+1),ibox=1,ncolprint)
1231	!4012 continue
1232	write (6,'(a)') 'skt(j):'
1233	write (6,'(8f7.2)') skt(j)
1234
1235	write (6,'(8I7)') (ibox,ibox=1,ncolprint)
1236
1237	write (6,'(a)') 'tau:'
1238	write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
1239
1240	write (6,'(a)') 'tb:'
1241	write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
1242
1243	write (6,'(a)') 'ptop:'
1244	write (6,'(8f7.2)') (ptop(j,ibox),ibox=1,ncolprint)
1245	endif
1246
1247	enddo
1248
1249	end if
1250
1251	return
1252	end
1253
1254

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

Download in other formats:

Original Format