Context Navigation

icarus.f90 @ 5310

Last change on this file since 5310 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

Rev	Line
[2428]	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 $
[5248]	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	)
[1262]	36
[5248]	37	!$Id: icarus.f,v 4.1 2010/05/27 16:30:18 hadmw Exp $
[1262]	38
[5248]	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*****************************
[1262]	77
[5248]	78	implicit none
[1262]	79
[5248]	80	! NOTE: the maximum number of levels and columns is set by
	81	! the following parameter statement
[1262]	82
[5248]	83	INTEGER :: ncolprint
[1262]	84
[5248]	85	! -----
	86	! Input
	87	! -----
[1262]	88
[5248]	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
[1262]	92
[5248]	93	INTEGER :: sunlit(npoints) ! 1 for day points, 0 for night time
[1262]	94
[5248]	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
[1262]	99
[5248]	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
[1262]	104
[5248]	105	REAL :: qv(npoints,nlev)
	106	! ! water vapor specific humidity (kg vapor/ kg air)
	107	! ! on full model levels
[1262]	108
[5248]	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
[1262]	113
[5248]	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
[1262]	119
[5248]	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
[1262]	127
[5248]	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.
[1262]	132
[5248]	133	INTEGER :: overlap ! overlap type
	134	! ! 1=max
	135	! ! 2=rand
	136	! ! 3=max/rand
[1262]	137
[5248]	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)
[1262]	151
[5248]	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.
[1262]	182
[5248]	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
[1262]	187
	188
	189
[5248]	190	! ------
	191	! Output
	192	! ------
[1262]	193
[5248]	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
[1262]	196
[5248]	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)
[1262]	203
	204
[5248]	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.
[1262]	207
[5248]	208	REAL :: meanptop(npoints) ! mean cloud top pressure (mb) - linear averaging
	209	! ! in cloud top pressure.
[1262]	210
[5248]	211	REAL :: meantaucld(npoints) ! mean optical thickness
	212	! ! linear averaging in albedo performed.
[1262]	213
[5248]	214	real :: meanalbedocld(npoints) ! mean cloud albedo
	215	! ! linear averaging in albedo performed
[1262]	216
[5248]	217	real :: meantb(npoints) ! mean all-sky 10.5 micron brightness temperature
[1262]	218
[5248]	219	real :: meantbclr(npoints) ! mean clear-sky 10.5 micron brightness temperature
[1262]	220
[5248]	221	REAL :: boxtau(npoints,ncol) ! optical thickness in each column
[1262]	222
[5248]	223	REAL :: boxptop(npoints,ncol) ! cloud top pressure (mb) in each column
[1262]	224
	225
[5248]	226	!
	227	! ------
	228	! Working variables added when program updated to mimic Mark Webb's PV-Wave code
	229	! ------
[1262]	230
[5248]	231	REAL :: dem(npoints,ncol),bb(npoints) ! working variables for 10.5 micron longwave
	232	! ! emissivity in part of
	233	! ! gridbox under consideration
[1262]	234
[5248]	235	REAL :: ptrop(npoints)
	236	REAL :: attrop(npoints)
	237	REAL :: attropmin (npoints)
	238	REAL :: atmax(npoints)
	239	REAL :: btcmin(npoints)
	240	REAL :: transmax(npoints)
[1262]	241
[5248]	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)
[1414]	249
[5248]	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)
[1262]	258
[5248]	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
[1262]	276
[5248]	277	integer :: index1(npoints),num1,jj,k1,k2
	278	real :: rec2p13,tauchk,logp,logp1,logp2,atd
	279	real :: output_missing_value
[1262]	280
[5248]	281	character(len=10) :: ftn09
[1262]	282
[5248]	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'/
[1262]	287
[5248]	288	! ------ End duplicate definitions common to wrapper routine
[1262]	289
[5248]	290	tauchk = -1.*log(0.9999999)
	291	rec2p13=1./2.13
[1262]	292
[5248]	293	ncolprint=0
[1262]	294
[5248]	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
[1262]	346
[5248]	347	! ---------------------------------------------------!
[1262]	348
[5248]	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
[1262]	355
[5248]	356	if (top_height .eq. 1 .or. top_height .eq. 3) then
[1262]	357
[5248]	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
[1262]	365
[5248]	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
[1262]	378
[5248]	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
[1262]	385
[5248]	386	end if
[1262]	387
	388
[5248]	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
[1262]	400
[5248]	401	! -----------------------------------------------------!
[1262]	402
[5248]	403	! ---------------------------------------------------!
[1262]	404
[5248]	405	do ilev=1,nlev
	406	do j=1,npoints
[1262]	407
[5248]	408	rangevec(j)=0
[1262]	409
[5248]	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
[1262]	415
[5248]	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
[1262]	421
[5248]	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
[1262]	426
[5248]	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
[1262]	431
[5248]	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
[1262]	437
[5248]	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
[1262]	444
[5248]	445	rangeerror=0
	446	do j=1,npoints
	447	rangeerror=rangeerror+rangevec(j)
	448	enddo
[1262]	449
[5248]	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
[1262]	457
[5248]	458	!
	459	! ---------------------------------------------------!
[1262]	460
	461
[5248]	462	!
	463	! ---------------------------------------------------!
	464	! COMPUTE CLOUD OPTICAL DEPTH FOR EACH COLUMN and
	465	! put into vector tau
[1262]	466
[5248]	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
[1262]	477
[5248]	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
[1262]	503	endif
[5248]	504	!
	505	! ---------------------------------------------------!
[1262]	506
	507
	508
[5248]	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.
[1262]	522
[5248]	523	if (top_height .eq. 1 .or. top_height .eq. 3) then
[1262]	524
	525
[5248]	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
[1262]	571
[5248]	572	! !initialize variables
	573	do j=1,npoints
	574	fluxtop_clrsky(j) = 0.
	575	trans_layers_above_clrsky(j)=1.
	576	enddo
[1262]	577
[5248]	578	do ilev=1,nlev
	579	do j=1,npoints
[1262]	580
[5248]	581	! ! Black body emission at temperature of the layer
[1262]	582
[5248]	583	bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. )
	584	! !bb(j)= 5.67e-8at(j,ilev)*4
[1262]	585
[5248]	586	! ! increase TOA flux by flux emitted from layer
	587	! ! times total transmittance in layers above
[1262]	588
[5248]	589	fluxtop_clrsky(j) = fluxtop_clrsky(j) &
	590	+ dem_wv(j,ilev)bb(j)trans_layers_above_clrsky(j)
[1262]	591
[5248]	592	! ! update trans_layers_above with transmissivity
	593	! ! from this layer for next time around loop
[1262]	594
[5248]	595	trans_layers_above_clrsky(j)= &
	596	trans_layers_above_clrsky(j)*(1.-dem_wv(j,ilev))
[1262]	597
	598
[5248]	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
[1262]	606
[5248]	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
[1262]	613
[5248]	614	enddo !loop over level
[1262]	615
[5248]	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
[1262]	620
[5248]	621	fluxtop_clrsky(j) = fluxtop_clrsky(j) + emsfc_lw * bb(j) &
	622	* trans_layers_above_clrsky(j)
[1262]	623
[5248]	624	! !clear sky brightness temperature
	625	meantbclr(j) = 1307.27/(log(1.+(1./fluxtop_clrsky(j))))
[1262]	626
[5248]	627	enddo
[1262]	628
[5248]	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'
[1262]	635
[5248]	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
[1262]	642
	643
[5248]	644	! !
	645	! ! END OF CLEAR SKY CALCULATION
	646	! !
	647	! !----------------------------------------------------------------
[1262]	648
	649
[5248]	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
[1262]	674	do j=1,npoints
[5248]	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)) )
[1262]	691	else
[5248]	692	dem(j,ibox)= dem_wv(j,ilev)
[1262]	693	end if
[5248]	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
[1262]	709	enddo ! j
	710	enddo ! ibox
	711
	712	if (ncolprint.ne.0) then
[5248]	713	do j=1,npoints,1000
	714	write (6,'(a)') 'ilev:'
	715	write (6,'(I2)') ilev
[1262]	716
	717	write(6,'(a10)') 'j='
	718	write(6,'(8I10)') j
[5248]	719	write (6,'(a)') 'emiss_layer:'
	720	write (6,'(8f7.2)') (dem(j,ibox),ibox=1,ncolprint)
[1262]	721
[5248]	722	write (6,'(a)') '100.*bb(j):'
	723	write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)
[1262]	724
[5248]	725	write (6,'(a)') '100.*f:'
	726	write (6,'(8f7.2)') &
	727	(100.*fluxtop(j,ibox),ibox=1,ncolprint)
	728
[1262]	729	write (6,'(a)') 'total_trans:'
[5248]	730	write (6,'(8f7.2)') &
	731	(trans_layers_above(j,ibox),ibox=1,ncolprint)
	732	enddo
[1262]	733	endif
	734
[5248]	735	enddo ! ilev
[1262]	736
	737
[5248]	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
[1262]	743
[5248]	744	do ibox=1,ncol
	745	do j=1,npoints
[1262]	746
[5248]	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
[1262]	843	end if
	844	end if
[5248]	845	end if
[1262]	846	enddo
	847	enddo
	848
[5248]	849	endif
[1262]	850
[5248]	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
[1262]	867
[5248]	868	if (ncolprint.ne.0) then
[1262]	869
[5248]	870	do j=1,npoints,1000
	871	write(6,'(a10)') 'j='
	872	write(6,'(8I10)') j
[1262]	873
[5248]	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
[1262]	952	enddo
[5248]	953	end do
[1262]	954
[5248]	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
[1262]	980	enddo ! j
	981
[5248]	982	else ! if (top_height .eq. 1 .or. top_height .eq. 3)
[1262]	983
[5248]	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
[1262]	997
[5248]	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
[1262]	1004
[5248]	1005	end do
[1262]	1006
[5248]	1007	!
	1008	!
	1009	! ---------------------------------------------------!
[1262]	1010
	1011
[5248]	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	!
[1262]	1030
[5248]	1031	! !compute isccp frequencies
[1262]	1032
[5248]	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
[1262]	1045
[5248]	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
[1262]	1060
[5248]	1061	boxarea = 1./real(ncol)
[1262]	1062
[5248]	1063	do ibox=1,ncol
	1064	do j=1,npoints
[1262]	1065
[5248]	1066	if (tau(j,ibox) .gt. (tauchk ) &
	1067	.and. ptop(j,ibox) .gt. 0.) then
	1068	box_cloudy(j,ibox)=.true.
	1069	endif
[1262]	1070
[5248]	1071	if (box_cloudy(j,ibox)) then
[1262]	1072
[5248]	1073	if (sunlit(j).eq.1 .or. top_height .eq. 3) then
[1262]	1074
[5248]	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
[1262]	1088	end if
	1089
[5248]	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
[1262]	1106	end if
	1107
[5248]	1108	! !reset itau(j), ipres(j)
	1109	itau(j) = 0
	1110	ipres(j) = 0
[1262]	1111
[5248]	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
[1262]	1133
[5248]	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
[1262]	1156
[5248]	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
[1262]	1162
[5248]	1163	end if
[1262]	1164
[5248]	1165	end if
[1262]	1166
[5248]	1167	enddo ! j
	1168	end do
[1262]	1169
[5248]	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
[1262]	1204	enddo
	1205
[5248]	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
[1262]	1213
[5248]	1214	! !print test
[1262]	1215
[5248]	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