Context Navigation

icarus.f90 @ 5312

Last change on this file since 5312 was 5263, checked in by abarral, 7 weeks ago
fix cospv1 compilation
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: 39.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(len=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(len=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 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	end do
378
379	do 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	end do
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 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	end do
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 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	end do
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 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 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	end do
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	end do
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 ilev=1,7
1035	do 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	end do
1044	end do
1045
1046	! !reset variables need for averaging cloud properties
1047	do j=1,npoints
1048	if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1049	totalcldarea(j) = 0.
1050	meanalbedocld(j) = 0.
1051	meanptop(j) = 0.
1052	meantaucld(j) = 0.
1053	else
1054	totalcldarea(j) = output_missing_value
1055	meanalbedocld(j) = output_missing_value
1056	meanptop(j) = output_missing_value
1057	meantaucld(j) = output_missing_value
1058	end if
1059	enddo ! j
1060
1061	boxarea = 1./real(ncol)
1062
1063	do ibox=1,ncol
1064	do j=1,npoints
1065
1066	if (tau(j,ibox) .gt. (tauchk ) &
1067	.and. ptop(j,ibox) .gt. 0.) then
1068	box_cloudy(j,ibox)=.true.
1069	endif
1070
1071	if (box_cloudy(j,ibox)) then
1072
1073	if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1074
1075	boxtau(j,ibox) = tau(j,ibox)
1076
1077	if (tau(j,ibox) .ge. isccp_taumin) then
1078	totalcldarea(j) = totalcldarea(j) + boxarea
1079
1080	! !convert optical thickness to albedo
1081	albedocld(j,ibox) &
1082	= (tau(j,ibox)0.895)/((tau(j,ibox)0.895)+6.82)
1083
1084	! !contribute to averaging
1085	meanalbedocld(j) = meanalbedocld(j) &
1086	+albedocld(j,ibox)*boxarea
1087
1088	end if
1089
1090	endif
1091
1092	endif
1093
1094	if (sunlit(j).eq.1 .or. top_height .eq. 3) then
1095
1096	if (box_cloudy(j,ibox)) then
1097
1098	! !convert ptop to millibars
1099	ptop(j,ibox)=ptop(j,ibox) / 100.
1100
1101	! !save for output cloud top pressure and optical thickness
1102	boxptop(j,ibox) = ptop(j,ibox)
1103
1104	if (tau(j,ibox) .ge. isccp_taumin) then
1105	meanptop(j) = meanptop(j) + ptop(j,ibox)*boxarea
1106	end if
1107
1108	! !reset itau(j), ipres(j)
1109	itau(j) = 0
1110	ipres(j) = 0
1111
1112	! !determine optical depth category
1113	if (tau(j,ibox) .lt. isccp_taumin) then
1114	itau(j)=1
1115	else if (tau(j,ibox) .ge. isccp_taumin &
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	end do
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 subroutine icarus
1269
1270

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

Download in other formats:

Original Format