Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

icarus.F @ 3663

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

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

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

Download in other formats:

Original Format