Context Navigation

icarus.F90 @ 5134

Last change on this file since 5134 was 3358, checked in by idelkadi, 6 years ago
Implementation de la nouvelle version COSPv2 dans LMDZ. Pour compiler avec makelmdz_fcma utiliser l'option "-cosp2 true"
File size: 30.1 KB

Rev	Line
[3358]	1	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	2	! Copyright (c) 2009, Lawrence Livemore National Security Limited Liability
	3	! All rights reserved.
	4	!
	5	! Redistribution and use in source and binary forms, with or without modification, are
	6	! permitted provided that the following conditions are met:
	7	!
	8	! 1. Redistributions of source code must retain the above copyright notice, this list of
	9	! conditions and the following disclaimer.
	10	!
	11	! 2. Redistributions in binary form must reproduce the above copyright notice, this list
	12	! of conditions and the following disclaimer in the documentation and/or other
	13	! materials provided with the distribution.
	14	!
	15	! 3. Neither the name of the copyright holder nor the names of its contributors may be
	16	! used to endorse or promote products derived from this software without specific prior
	17	! written permission.
	18	!
	19	! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
	20	! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
	21	! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
	22	! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
	23	! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
	24	! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
	25	! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
	26	! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
	27	! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
	28	!
	29	! History
	30	! May 2015 - D. Swales - Modified for COSPv2.0
	31	! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	32	MODULE MOD_ICARUS
	33	USE COSP_KINDS, ONLY: wp
	34	USE COSP_PHYS_CONSTANTS, ONLY: amd,amw,avo,grav
	35	use MOD_COSP_STATS, ONLY: hist2D
	36	USE MOD_COSP_CONFIG, ONLY: R_UNDEF,numISCCPTauBins,numISCCPPresBins,isccp_histTau, &
	37	isccp_histPres
	38	implicit none
	39
	40	! Shared Parameters
	41	integer,parameter :: &
	42	ncolprint = 0 ! Flag for debug printing (set as parameter to increase performance)
	43
	44	! Cloud-top height determination
	45	integer :: &
	46	isccp_top_height, & ! Top height adjustment method
	47	isccp_top_height_direction ! Direction for finding atmosphere pressure level
	48
	49	! Parameters used by icarus
	50	real(wp),parameter :: &
	51	tauchk = -1._wp*log(0.9999999_wp), & ! Lower limit on optical depth
	52	isccp_taumin = 0.3_wp, & ! Minimum optical depth for joint-hostogram
	53	pstd = 1013250._wp, & ! Mean sea-level pressure (Pa)
	54	isccp_t0 = 296._wp, & ! Mean surface temperature (K)
	55	output_missing_value = -1.E+30 ! Missing values
	56
	57	contains
	58	! ##########################################################################
	59	! ##########################################################################
	60	SUBROUTINE ICARUS(debug,debugcol,npoints,sunlit,nlev,ncol,pfull, &
	61	phalf,qv,cc,conv,dtau_s,dtau_c,th,thd,frac_out,skt,emsfc_lw,at,&
	62	dem_s,dem_c,fq_isccp,totalcldarea, meanptop,meantaucld, &
	63	meanalbedocld, meantb,meantbclr,boxtau,boxptop,levmatch)
	64
	65	! INPUTS
	66	INTEGER,intent(in) :: & !
	67	npoints, & ! Number of model points in the horizontal
	68	nlev, & ! Number of model levels in column
	69	ncol, & ! Number of subcolumns
	70	debug, & ! Debug flag
	71	debugcol ! Debug column flag
	72	INTEGER,intent(in),dimension(npoints) :: & !
	73	sunlit ! 1 for day points, 0 for night time
	74	REAL(WP),intent(in) :: & !
	75	emsfc_lw ! 10.5 micron emissivity of surface (fraction)
	76	REAL(WP),intent(in),dimension(npoints) :: & !
	77	skt ! Skin Temperature (K)
	78	REAL(WP),intent(in),dimension(npoints,ncol,nlev) :: & !
	79	frac_out ! Boxes gridbox divided up into subcolumns
	80	REAL(WP),intent(in),dimension(npoints,nlev) :: & !
	81	pfull, & ! Pressure of full model levels (Pascals)
	82	qv, & ! Water vapor specific humidity (kg vapor/ kg air)
	83	cc, & ! Cloud cover in each model level (fraction)
	84	conv, & ! Convective cloud cover in each model
	85	at, & ! Temperature in each model level (K)
	86	dem_c, & ! Emissivity for convective clouds
	87	dem_s, & ! Emissivity for stratiform clouds
	88	dtau_c, & ! Optical depth for convective clouds
	89	dtau_s ! Optical depth for stratiform clouds
	90	REAL(WP),intent(in),dimension(npoints,nlev+1) :: & !
	91	phalf ! Pressure of half model levels (Pascals)!
	92	integer,intent(in) :: th,thd
	93
	94	! OUTPUTS
	95	REAL(WP),intent(out),dimension(npoints,7,7) :: &
	96	fq_isccp ! The fraction of the model grid box covered by clouds
	97	REAL(WP),intent(out),dimension(npoints) :: &
	98	totalcldarea, & ! The fraction of model grid box columns with cloud present
	99	meanptop, & ! Mean cloud top pressure (mb) - linear averaging
	100	meantaucld, & ! Mean optical thickness
	101	meanalbedocld, & ! Mean cloud albedo
	102	meantb, & ! Mean all-sky 10.5 micron brightness temperature
	103	meantbclr ! Mean clear-sky 10.5 micron brightness temperature
	104	REAL(WP),intent(out),dimension(npoints,ncol) :: &
	105	boxtau, & ! Optical thickness in each column
	106	boxptop ! Cloud top pressure (mb) in each column
	107	INTEGER,intent(out),dimension(npoints,ncol) :: &
	108	levmatch ! Used for icarus unit testing only
	109
	110
	111	! INTERNAL VARIABLES
	112	CHARACTER(len=10) :: ftn09
	113	REAL(WP),dimension(npoints,ncol) :: boxttop
	114	REAL(WP),dimension(npoints,ncol,nlev) :: dtau,demIN
	115	INTEGER :: j,ilev,ibox
	116	INTEGER,dimension(nlev,ncol ) :: acc
	117
	118	! PARAMETERS
	119	character ,parameter, dimension(6) :: cchar=(/' ','-','1','+','I','+'/)
	120	character(len=1),parameter,dimension(6) :: cchar_realtops=(/ ' ',' ','1','1','I','I'/)
	121	! ##########################################################################
	122
	123	call cosp_simulator_optics(npoints,ncol,nlev,frac_out,dem_c,dem_s,demIN)
	124	call cosp_simulator_optics(npoints,ncol,nlev,frac_out,dtau_c,dtau_s,dtau)
	125
	126	call ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demIN,skt,emsfc_lw,qv,at, &
	127	pfull,phalf,frac_out,levmatch,boxtau,boxptop,boxttop,meantbclr)
	128
	129	call ICARUS_COLUMN(npoints,ncol,boxtau,boxptop/100._wp,sunlit,boxttop,&
	130	fq_isccp,meanalbedocld,meanptop,meantaucld,totalcldarea,meantb)
	131
	132	! ##########################################################################
	133	! OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM
	134	! ##########################################################################
	135
	136	if (debugcol.ne.0) then
	137	do j=1,npoints,debugcol
	138
	139	! Produce character output
	140	do ilev=1,nlev
	141	acc(ilev,1:ncol)=frac_out(j,1:ncol,ilev)*2
	142	where(levmatch(j,1:ncol) .eq. ilev) acc(ilev,1:ncol)=acc(ilev,1:ncol)+1
	143	enddo
	144
	145	write(ftn09,11) j
	146	11 format('ftn09.',i4.4)
	147	open(9, FILE=ftn09, FORM='FORMATTED')
	148
	149	write(9,'(a1)') ' '
	150	write(9,'(10i5)') (ilev,ilev=5,nlev,5)
	151	write(9,'(a1)') ' '
	152
	153	do ibox=1,ncol
	154	write(9,'(40(a1),1x,40(a1))') &
	155	(cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev),&
	156	(cchar(acc(ilev,ibox)+1),ilev=1,nlev)
	157	end do
	158	close(9)
	159
	160	enddo
	161	end if
	162
	163	return
	164	end SUBROUTINE ICARUS
	165
	166	! ############################################################################
	167	! ############################################################################
	168	! ############################################################################
	169	SUBROUTINE ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demiN,skt,emsfc_lw,qv,at, &
	170	pfull,phalf,frac_out,levmatch,boxtau,boxptop,boxttop,meantbclr)
	171	! Inputs
	172	INTEGER, intent(in) :: &
	173	ncol, & ! Number of subcolumns
	174	npoints, & ! Number of horizontal gridpoints
	175	nlev ! Number of vertical levels
	176	INTEGER, intent(in), dimension(npoints) :: &
	177	sunlit ! 1=day 0=night
	178	REAL(WP),intent(in) :: &
	179	emsfc_lw ! 10.5 micron emissivity of surface (fraction)
	180	REAL(WP),intent(in), dimension(npoints) :: &
	181	skt ! Skin temperature
	182	REAL(WP),intent(in), dimension(npoints,nlev) :: &
	183	at, & ! Temperature
	184	pfull, & ! Presure
	185	qv ! Specific humidity
	186	REAL(WP),intent(in), dimension(npoints,ncol,nlev) :: &
	187	frac_out, & ! Subcolumn cloud cover
	188	dtau, & ! Subcolumn optical thickness
	189	demIN ! Subcolumn emissivity
	190	REAL(WP),intent(in), dimension(npoints,nlev+1) :: &
	191	phalf ! Pressure at model half levels
	192
	193	! Outputs
	194	REAL(WP),intent(inout),dimension(npoints) :: &
	195	meantbclr ! Mean clear-sky 10.5 micron brightness temperature
	196	REAL(WP),intent(inout),dimension(npoints,ncol) :: &
	197	boxtau, & ! Optical thickness in each column
	198	boxptop, & ! Cloud top pressure (mb) in each column
	199	boxttop ! Cloud top temperature in each column
	200	INTEGER, intent(inout),dimension(npoints,ncol) :: levmatch
	201
	202	! Local Variables
	203	INTEGER :: &
	204	j,ibox,ilev,k1,k2,icycle
	205	INTEGER,dimension(npoints) :: &
	206	nmatch,itrop
	207	INTEGER,dimension(npoints,nlev-1) :: &
	208	match
	209	REAL(WP) :: &
	210	logp,logp1,logp2,atd
	211	REAL(WP),dimension(npoints) :: &
	212	bb,attropmin,attrop,ptrop,atmax,btcmin,transmax,tauir,taumin,fluxtopinit,press, &
	213	dpress,atmden,rvh20,rhoave,rh20s,rfrgn,tmpexp,tauwv,wk,trans_layers_above_clrsky, &
	214	fluxtop_clrsky
	215	REAL(WP),dimension(npoints,nlev) :: &
	216	dem_wv
	217	REAL(WP),dimension(npoints,ncol) :: &
	218	trans_layers_above,dem,tb,emcld,fluxtop,tau,ptop
	219
	220	! ####################################################################################
	221	! Compute cloud optical depth for each column by summing up subcolumns
	222	tau(1:npoints,1:ncol) = 0._wp
	223	tau(1:npoints,1:ncol) = sum(dtau,dim=3)
	224
	225	! Set tropopause values
	226	if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then
	227	ptrop(1:npoints) = 5000._wp
	228	attropmin(1:npoints) = 400._wp
	229	atmax(1:npoints) = 0._wp
	230	attrop(1:npoints) = 120._wp
	231	itrop(1:npoints) = 1
	232
	233	do ilev=1,nlev
	234	where(pfull(1:npoints,ilev) .lt. 40000. .and. &
	235	pfull(1:npoints,ilev) .gt. 5000. .and. &
	236	at(1:npoints,ilev) .lt. attropmin(1:npoints))
	237	ptrop(1:npoints) = pfull(1:npoints,ilev)
	238	attropmin(1:npoints) = at(1:npoints,ilev)
	239	attrop(1:npoints) = attropmin(1:npoints)
	240	itrop = ilev
	241	endwhere
	242	enddo
	243
	244	do ilev=1,nlev
	245	atmax(1:npoints) = merge(at(1:npoints,ilev),atmax(1:npoints),&
	246	at(1:npoints,ilev) .gt. atmax(1:npoints) .and. ilev .ge. itrop(1:npoints))
	247	enddo
	248	end if
	249
	250	if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then
	251	! ############################################################################
	252	! Clear-sky radiance calculation
	253	!
	254	! Compute water vapor continuum emissivity this treatment follows Schwarkzopf
	255	! and Ramasamy JGR 1999,vol 104, pages 9467-9499. The emissivity is calculated
	256	! at a wavenumber of 955 cm-1, or 10.47 microns
	257	! ############################################################################
	258	do ilev=1,nlev
	259	press(1:npoints) = pfull(1:npoints,ilev)*10._wp
	260	dpress(1:npoints) = (phalf(1:npoints,ilev+1)-phalf(1:npoints,ilev))*10
	261	atmden(1:npoints) = dpress(1:npoints)/(grav*100._wp)
	262	rvh20(1:npoints) = qv(1:npoints,ilev)*amd/amw
	263	wk(1:npoints) = rvh20(1:npoints)avoatmden/amd
	264	rhoave(1:npoints) = (press(1:npoints)/pstd)*(isccp_t0/at(1:npoints,ilev))
	265	rh20s(1:npoints) = rvh20(1:npoints)*rhoave(1:npoints)
	266	rfrgn(1:npoints) = rhoave(1:npoints)-rh20s(1:npoints)
	267	tmpexp(1:npoints) = exp(-0.02_wp*(at(1:npoints,ilev)-isccp_t0))
	268	tauwv(1:npoints) = wk(1:npoints)1.e-20((0.0224697_wprh20s(1:npoints) &
	269	tmpexp(1:npoints))+(3.41817e-7rfrgn(1:npoints)))0.98_wp
	270	dem_wv(1:npoints,ilev) = 1._wp - exp( -1._wp * tauwv(1:npoints))
	271	enddo
	272
	273	fluxtop_clrsky(1:npoints) = 0._wp
	274	trans_layers_above_clrsky(1:npoints) = 1._wp
	275	do ilev=1,nlev
	276	! Black body emission at temperature of the layer
	277	bb(1:npoints) = 1._wp / ( exp(1307.27_wp/at(1:npoints,ilev)) - 1._wp )
	278
	279	! Increase TOA flux by flux emitted from layer times total transmittance in layers above
	280	fluxtop_clrsky(1:npoints) = fluxtop_clrsky(1:npoints) + &
	281	dem_wv(1:npoints,ilev)bb(1:npoints)trans_layers_above_clrsky(1:npoints)
	282
	283	! Update trans_layers_above with transmissivity from this layer for next time around loop
	284	trans_layers_above_clrsky(1:npoints) = trans_layers_above_clrsky(1:npoints)*&
	285	(1.-dem_wv(1:npoints,ilev))
	286	enddo
	287
	288	! Add in surface emission
	289	bb(1:npoints) = 1._wp/( exp(1307.27_wp/skt(1:npoints)) - 1._wp )
	290	fluxtop_clrsky(1:npoints) = fluxtop_clrsky(1:npoints) + &
	291	emsfc_lw * bb(1:npoints)*trans_layers_above_clrsky(1:npoints)
	292
	293	! Clear Sky brightness temperature
	294	meantbclr(1:npoints) = 1307.27_wp/(log(1._wp+(1._wp/fluxtop_clrsky(1:npoints))))
	295
	296	! #################################################################################
	297	! All-sky radiance calculation
	298	! #################################################################################
	299
	300	fluxtop(1:npoints,1:ncol) = 0._wp
	301	trans_layers_above(1:npoints,1:ncol) = 1._wp
	302	do ilev=1,nlev
	303	! Black body emission at temperature of the layer
	304	bb=1._wp/(exp(1307.27_wp/at(1:npoints,ilev)) - 1._wp)
	305
	306	do ibox=1,ncol
	307	! Emissivity
	308	dem(1:npoints,ibox) = merge(dem_wv(1:npoints,ilev), &
	309	1._wp-(1._wp-demIN(1:npoints,ibox,ilev))*(1._wp-dem_wv(1:npoints,ilev)), &
	310	demIN(1:npoints,ibox,ilev) .eq. 0)
	311
	312	! Increase TOA flux emitted from layer
	313	fluxtop(1:npoints,ibox) = fluxtop(1:npoints,ibox) + dem(1:npoints,ibox)bbtrans_layers_above(1:npoints,ibox)
	314
	315	! Update trans_layer by emitted layer from above
	316	trans_layers_above(1:npoints,ibox) = trans_layers_above(1:npoints,ibox)*(1._wp-dem(1:npoints,ibox))
	317	enddo
	318	enddo
	319
	320	! Add in surface emission
	321	bb(1:npoints)=1._wp/( exp(1307.27_wp/skt(1:npoints)) - 1._wp )
	322	do ibox=1,ncol
	323	fluxtop(1:npoints,ibox) = fluxtop(1:npoints,ibox) + emsfc_lwbb(1:npoints)trans_layers_above(1:npoints,ibox)
	324	end do
	325
	326	! All Sky brightness temperature
	327	boxttop(1:npoints,1:ncol) = 1307.27_wp/(log(1._wp+(1._wp/fluxtop(1:npoints,1:ncol))))
	328
	329	! #################################################################################
	330	! Cloud-Top Temperature
	331	!
	332	! Now that you have the top of atmosphere radiance, account for ISCCP
	333	! procedures to determine cloud top temperature account for partially
	334	! transmitting cloud recompute flux ISCCP would see assuming a single layer
	335	! cloud. NOTE choice here of 2.13, as it is primarily ice clouds which have
	336	! partial emissivity and need the adjustment performed in this section. If it
	337	! turns out that the cloud brightness temperature is greater than 260K, then
	338	! the liquid cloud conversion factor of 2.56 is used. NOTE that this is
	339	! discussed on pages 85-87 of the ISCCP D level documentation
	340	! (Rossow et al. 1996)
	341	! #################################################################################
	342
	343	! Compute minimum brightness temperature and optical depth
	344	btcmin(1:npoints) = 1._wp / ( exp(1307.27_wp/(attrop(1:npoints)-5._wp)) - 1._wp )
	345
	346	do ibox=1,ncol
	347	transmax(1:npoints) = (fluxtop(1:npoints,ibox)-btcmin) /(fluxtop_clrsky(1:npoints)-btcmin(1:npoints))
	348	tauir(1:npoints) = tau(1:npoints,ibox)/2.13_wp
	349	taumin(1:npoints) = -log(max(min(transmax(1:npoints),0.9999999_wp),0.001_wp))
	350	if (isccp_top_height .eq. 1) then
	351	do j=1,npoints
	352	if (transmax(j) .gt. 0.001 .and. transmax(j) .le. 0.9999999) then
	353	fluxtopinit(j) = fluxtop(j,ibox)
	354	tauir(j) = tau(j,ibox)/2.13_wp
	355	endif
	356	enddo
	357	do icycle=1,2
	358	do j=1,npoints
	359	if (tau(j,ibox) .gt. (tauchk)) then
	360	if (transmax(j) .gt. 0.001 .and. transmax(j) .le. 0.9999999) then
	361	emcld(j,ibox) = 1._wp - exp(-1._wp * tauir(j) )
	362	fluxtop(j,ibox) = fluxtopinit(j) - ((1.-emcld(j,ibox))*fluxtop_clrsky(j))
	363	fluxtop(j,ibox)=max(1.E-06_wp,(fluxtop(j,ibox)/emcld(j,ibox)))
	364	tb(j,ibox)= 1307.27_wp / (log(1._wp + (1._wp/fluxtop(j,ibox))))
	365	if (tb(j,ibox) .gt. 260.) then
	366	tauir(j) = tau(j,ibox) / 2.56_wp
	367	end if
	368	end if
	369	end if
	370	enddo
	371	enddo
	372	endif
	373
	374	! Cloud-top temperature
	375	where(tau(1:npoints,ibox) .gt. tauchk)
	376	tb(1:npoints,ibox)= 1307.27_wp/ (log(1. + (1._wp/fluxtop(1:npoints,ibox))))
	377	where (isccp_top_height .eq. 1 .and. tauir(1:npoints) .lt. taumin(1:npoints))
	378	tb(1:npoints,ibox) = attrop(1:npoints) - 5._wp
	379	tau(1:npoints,ibox) = 2.13_wp*taumin(1:npoints)
	380	endwhere
	381	endwhere
	382
	383	! Clear-sky brightness temperature
	384	where(tau(1:npoints,ibox) .le. tauchk)
	385	tb(1:npoints,ibox) = meantbclr(1:npoints)
	386	endwhere
	387	enddo
	388	else
	389	meantbclr(1:npoints) = output_missing_value
	390	end if
	391
	392	! ####################################################################################
	393	! Cloud-Top Pressure
	394	!
	395	! The 2 methods differ according to whether or not you use the physical cloud
	396	! top pressure (isccp_top_height = 2) or the radiatively determined cloud top
	397	! pressure (isccp_top_height = 1 or 3)
	398	! ####################################################################################
	399	do ibox=1,ncol
	400	!segregate according to optical thickness
	401	if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then
	402
	403	! Find level whose temperature most closely matches brightness temperature
	404	nmatch(1:npoints)=0
	405	do k1=1,nlev-1
	406	ilev = merge(nlev-k1,k1,isccp_top_height_direction .eq. 2)
	407	do j=1,npoints
	408	if (ilev .ge. itrop(j) .and. &
	409	((at(j,ilev) .ge. tb(j,ibox) .and. &
	410	at(j,ilev+1) .le. tb(j,ibox)) .or. &
	411	(at(j,ilev) .le. tb(j,ibox) .and. &
	412	at(j,ilev+1) .ge. tb(j,ibox)))) then
	413	nmatch(j)=nmatch(j)+1
	414	match(j,nmatch(j))=ilev
	415	endif
	416	enddo
	417	enddo
	418
	419	do j=1,npoints
	420	if (nmatch(j) .ge. 1) then
	421	k1 = match(j,nmatch(j))
	422	k2 = k1 + 1
	423	logp1 = log(pfull(j,k1))
	424	logp2 = log(pfull(j,k2))
	425	atd = max(tauchk,abs(at(j,k2) - at(j,k1)))
	426	logp=logp1+(logp2-logp1)*abs(tb(j,ibox)-at(j,k1))/atd
	427	ptop(j,ibox) = exp(logp)
	428	levmatch(j,ibox) = merge(k1,k2,abs(pfull(j,k1)-ptop(j,ibox)) .lt. abs(pfull(j,k2)-ptop(j,ibox)))
	429	else
	430	if (tb(j,ibox) .le. attrop(j)) then
	431	ptop(j,ibox)=ptrop(j)
	432	levmatch(j,ibox)=itrop(j)
	433	end if
	434	if (tb(j,ibox) .ge. atmax(j)) then
	435	ptop(j,ibox)=pfull(j,nlev)
	436	levmatch(j,ibox)=nlev
	437	end if
	438	end if
	439	enddo
	440	else
	441	ptop(1:npoints,ibox)=0.
	442	do ilev=1,nlev
	443	where((ptop(1:npoints,ibox) .eq. 0. ) .and.(frac_out(1:npoints,ibox,ilev) .ne. 0))
	444	ptop(1:npoints,ibox)=phalf(1:npoints,ilev)
	445	levmatch(1:npoints,ibox)=ilev
	446	endwhere
	447	end do
	448	end if
	449	where(tau(1:npoints,ibox) .le. tauchk)
	450	ptop(1:npoints,ibox)=0._wp
	451	levmatch(1:npoints,ibox)=0._wp
	452	endwhere
	453	enddo
	454
	455	! ####################################################################################
	456	! Compute subcolumn pressure and optical depth
	457	! ####################################################################################
	458	boxtau(1:npoints,1:ncol) = output_missing_value
	459	boxptop(1:npoints,1:ncol) = output_missing_value
	460	do ibox=1,ncol
	461	do j=1,npoints
	462	if (tau(j,ibox) .gt. (tauchk) .and. ptop(j,ibox) .gt. 0.) then
	463	if (sunlit(j).eq.1 .or. isccp_top_height .eq. 3) then
	464	boxtau(j,ibox) = tau(j,ibox)
	465	boxptop(j,ibox) = ptop(j,ibox)!/100._wp
	466	endif
	467	endif
	468	enddo
	469	enddo
	470
	471	end SUBROUTINE ICARUS_SUBCOLUMN
	472
	473	! ######################################################################################
	474	! SUBROUTINE icarus_column
	475	! ######################################################################################
	476	SUBROUTINE ICARUS_column(npoints,ncol,boxtau,boxptop,sunlit,boxttop,fq_isccp, &
	477	meanalbedocld,meanptop,meantaucld,totalcldarea,meantb)
	478	! Inputs
	479	INTEGER, intent(in) :: &
	480	ncol, & ! Number of subcolumns
	481	npoints ! Number of horizontal gridpoints
	482	INTEGER, intent(in),dimension(npoints) :: &
	483	sunlit ! day=1 night=0
	484	REAL(WP),intent(in),dimension(npoints,ncol) :: &
	485	boxttop, & ! Subcolumn top temperature
	486	boxptop, & ! Subcolumn cloud top pressure
	487	boxtau ! Subcolumn optical depth
	488
	489	! Outputs
	490	REAL(WP),intent(inout),dimension(npoints) :: &
	491	meanalbedocld, & ! Gridmean cloud albedo
	492	meanptop, & ! Gridmean cloud top pressure (mb) - linear averaging
	493	meantaucld, & ! Gridmean optical thickness
	494	totalcldarea, & ! The fraction of model grid box columns with cloud present
	495	meantb ! Gridmean all-sky 10.5 micron brightness temperature
	496	REAL(WP),intent(inout),dimension(npoints,7,7) :: &
	497	fq_isccp ! The fraction of the model grid box covered by clouds
	498
	499	! Local Variables
	500	INTEGER :: j,ilev,ilev2
	501	REAL(WP),dimension(npoints,ncol) :: albedocld
	502	LOGICAL, dimension(npoints,ncol) :: box_cloudy
	503
	504	! Variables for new joint-histogram implementation
	505	logical,dimension(ncol) :: box_cloudy2
	506
	507	! ####################################################################################
	508	! Brightness Temperature
	509	! ####################################################################################
	510	if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then
	511	meantb(1:npoints)=sum(boxttop,2)/ncol
	512	else
	513	meantb(1:npoints) = output_missing_value
	514	endif
	515
	516	! ####################################################################################
	517	! Determines ISCCP cloud type frequencies
	518	!
	519	! Now that boxptop and boxtau have been determined, determine amount of each of the
	520	! 49 ISCCP cloud types. Also compute grid box mean cloud top pressure and
	521	! optical thickness. The mean cloud top pressure and optical thickness are
	522	! averages over the cloudy area only. The mean cloud top pressure is a linear
	523	! average of the cloud top pressures. The mean cloud optical thickness is
	524	! computed by converting optical thickness to an albedo, averaging in albedo
	525	! units, then converting the average albedo back to a mean optical thickness.
	526	! ####################################################################################
	527
	528	! Initialize
	529	albedocld(1:npoints,1:ncol) = 0._wp
	530	box_cloudy(1:npoints,1:ncol) = .false.
	531
	532	! Reset frequencies
	533	!fq_isccp = spread(spread(merge(0._wp,output_missing_value,sunlit .eq. 1 .or. isccp_top_height .eq. 3),2,7),2,7)
	534	do ilev=1,7
	535	do ilev2=1,7
	536	do j=1,npoints !
	537	if (sunlit(j).eq.1 .or. isccp_top_height .eq. 3) then
	538	fq_isccp(j,ilev,ilev2)= 0.
	539	else
	540	fq_isccp(j,ilev,ilev2)= output_missing_value
	541	end if
	542	enddo
	543	enddo
	544	enddo
	545
	546
	547	! Reset variables need for averaging cloud properties
	548	where(sunlit .eq. 1 .or. isccp_top_height .eq. 3)
	549	totalcldarea(1:npoints) = 0._wp
	550	meanalbedocld(1:npoints) = 0._wp
	551	meanptop(1:npoints) = 0._wp
	552	meantaucld(1:npoints) = 0._wp
	553	elsewhere
	554	totalcldarea(1:npoints) = output_missing_value
	555	meanalbedocld(1:npoints) = output_missing_value
	556	meanptop(1:npoints) = output_missing_value
	557	meantaucld(1:npoints) = output_missing_value
	558	endwhere
	559
	560	! Compute column quantities and joint-histogram
	561	do j=1,npoints
	562	! Subcolumns that are cloudy(true) and not(false)
	563	box_cloudy2(1:ncol) = merge(.true.,.false.,boxtau(j,1:ncol) .gt. tauchk .and. boxptop(j,1:ncol) .gt. 0.)
	564
	565	! Compute joint histogram and column quantities for points that are sunlit and cloudy
	566	if (sunlit(j) .eq.1 .or. isccp_top_height .eq. 3) then
	567	! Joint-histogram
	568	call hist2D(boxtau(j,1:ncol),boxptop(j,1:ncol),ncol,isccp_histTau,numISCCPTauBins, &
	569	isccp_histPres,numISCCPPresBins,fq_isccp(j,1:numISCCPTauBins,1:numISCCPPresBins))
	570	fq_isccp(j,1:numISCCPTauBins,1:numISCCPPresBins) = &
	571	fq_isccp(j,1:numISCCPTauBins,1:numISCCPPresBins)/ncol
	572
	573	! Column cloud area
	574	totalcldarea(j) = real(count(box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin))/ncol
	575
	576	! Subcolumn cloud albedo
	577	!albedocld(j,1:ncol) = merge((boxtau(j,1:ncol)0.895_wp)/((boxtau(j,1:ncol)0.895_wp)+6.82_wp),&
	578	! 0._wp,box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin)
	579	where(box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin)
	580	albedocld(j,1:ncol) = (boxtau(j,1:ncol)0.895_wp)/((boxtau(j,1:ncol)0.895_wp)+6.82_wp)
	581	elsewhere
	582	albedocld(j,1:ncol) = 0._wp
	583	endwhere
	584
	585	! Column cloud albedo
	586	meanalbedocld(j) = sum(albedocld(j,1:ncol))/ncol
	587
	588	! Column cloud top pressure
	589	meanptop(j) = sum(boxptop(j,1:ncol),box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin)/ncol
	590	endif
	591	enddo
	592
	593	! Compute mean cloud properties. Set to mssing value in the event that totalcldarea=0
	594	where(totalcldarea(1:npoints) .gt. 0)
	595	meanptop(1:npoints) = 100._wp*meanptop(1:npoints)/totalcldarea(1:npoints)
	596	meanalbedocld(1:npoints) = meanalbedocld(1:npoints)/totalcldarea(1:npoints)
	597	meantaucld(1:npoints) = (6.82_wp/((1._wp/meanalbedocld(1:npoints))-1.))**(1._wp/0.895_wp)
	598	elsewhere
	599	meanptop(1:nPoints) = output_missing_value
	600	meanalbedocld(1:nPoints) = output_missing_value
	601	meantaucld(1:nPoints) = output_missing_value
	602	endwhere
	603	!meanptop(1:npoints) = merge(100._wp*meanptop(1:npoints)/totalcldarea(1:npoints),&
	604	! output_missing_value,totalcldarea(1:npoints) .gt. 0)
	605	!meanalbedocld(1:npoints) = merge(meanalbedocld(1:npoints)/totalcldarea(1:npoints), &
	606	! output_missing_value,totalcldarea(1:npoints) .gt. 0)
	607	!meantaucld(1:npoints) = merge((6.82_wp/((1._wp/meanalbedocld(1:npoints))-1.))**(1._wp/0.895_wp), &
	608	! output_missing_value,totalcldarea(1:npoints) .gt. 0)
	609
	610	! Represent in percent
	611	where(totalcldarea .ne. output_missing_value) totalcldarea = totalcldarea*100._wp
	612	where(fq_isccp .ne. output_missing_value) fq_isccp = fq_isccp*100._wp
	613
	614
	615	end SUBROUTINE ICARUS_column
	616
	617	subroutine cosp_simulator_optics(dim1,dim2,dim3,flag,varIN1,varIN2,varOUT)
	618	! INPUTS
	619	integer,intent(in) :: &
	620	dim1, & ! Dimension 1 extent (Horizontal)
	621	dim2, & ! Dimension 2 extent (Subcolumn)
	622	dim3 ! Dimension 3 extent (Vertical)
	623	real(wp),intent(in),dimension(dim1,dim2,dim3) :: &
	624	flag ! Logical to determine the of merge var1IN and var2IN
	625	real(wp),intent(in),dimension(dim1, dim3) :: &
	626	varIN1, & ! Input field 1
	627	varIN2 ! Input field 2
	628	! OUTPUTS
	629	real(wp),intent(out),dimension(dim1,dim2,dim3) :: &
	630	varOUT ! Merged output field
	631	! LOCAL VARIABLES
	632	integer :: j
	633
	634	varOUT(1:dim1,1:dim2,1:dim3) = 0._wp
	635	do j=1,dim2
	636	where(flag(:,j,:) .eq. 1)
	637	varOUT(:,j,:) = varIN2
	638	endwhere
	639	where(flag(:,j,:) .eq. 2)
	640	varOUT(:,j,:) = varIN1
	641	endwhere
	642	enddo
	643	end subroutine cosp_simulator_optics
	644	end module MOD_ICARUS
	645

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

Download in other formats:

Original Format