Context Navigation

icarus.F @ 3241

Last change on this file since 3241 was 1603, checked in by idelkadi, 13 years ago
Appel de la fonction flush desactive
File size: 42.8 KB

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

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

Download in other formats:

Original Format