Context Navigation

aeropacity.F90 @ 537

Last change on this file since 537 was 526, checked in by jleconte, 13 years ago

13/02/2012 == JL + AS

All outputs are now in netCDF format. Even in 1D (No more G1D)
Clean up of the call to callcorrk when CLFvarying=true
Corrects a bug in writediagspecIR/VI. Output are now in W/m2/cm-1 as a function of the wavenumber in cm-1
Enable writediagspecIR/V to work in the CLFvarying=true case (output now done in Physiq after writediagfi)
Add a simple treatment for the supersaturation of CO2 (see forget et al 2012)
corrects a small bug when no clouds are present in aeropacity

File size: 10.5 KB

Rev	Line
[253]	1	Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pq, &
	2	aerosol,reffrad,QREFvis3d,QREFir3d,tau_col,cloudfrac,totcloudfrac,clearsky)
[135]	3
[253]	4	use radinc_h, only : naerkind, L_TAUMAX
[135]	5
	6	implicit none
	7
	8	!==================================================================
	9	!
	10	! Purpose
	11	! -------
	12	! Compute aerosol optical depth in each gridbox.
	13	!
	14	! Authors
	15	! -------
	16	! F. Forget
	17	! F. Montmessin (water ice scheme)
	18	! update J.-B. Madeleine (2008)
	19	! dust removal, simplification by Robin Wordsworth (2009)
	20	!
	21	! Input
	22	! -----
	23	! ngrid Number of horizontal gridpoints
	24	! nlayer Number of layers
	25	! nq Number of tracers
	26	! pplev Pressure (Pa) at each layer boundary
	27	! pq Aerosol mixing ratio
	28	! reffrad(ngrid,nlayer,naerkind) Aerosol effective radius
	29	! QREFvis3d(ngridmx,nlayermx,naerkind) \ 3d extinction coefficients
	30	! QREFir3d(ngridmx,nlayermx,naerkind) / at reference wavelengths
	31	!
	32	! Output
	33	! ------
	34	! aerosol Aerosol optical depth in layer l, grid point ig
	35	! tau_col Total column optical depth at grid point ig
	36	!
	37	!=======================================================================
	38
	39	#include "dimensions.h"
	40	#include "dimphys.h"
	41	#include "callkeys.h"
	42	#include "comcstfi.h"
	43	#include "comgeomfi.h"
	44	#include "tracer.h"
	45
	46	INTEGER ngrid,nlayer,nq
	47
[253]	48	REAL pplay(ngrid,nlayer)
[135]	49	REAL pplev(ngrid,nlayer+1)
	50	REAL pq(ngrid,nlayer,nq)
	51	REAL aerosol(ngrid,nlayer,naerkind)
	52	REAL reffrad(ngrid,nlayer,naerkind)
	53	REAL QREFvis3d(ngridmx,nlayermx,naerkind)
	54	REAL QREFir3d(ngridmx,nlayermx,naerkind)
	55
[253]	56	REAL tau_col(ngrid)
	57	! REAL tauref(ngrid), tau_col(ngrid)
[135]	58
[253]	59	real cloudfrac(ngridmx,nlayermx)
	60	real aerosol0
	61
[135]	62	INTEGER l,ig,iq,iaer
	63
	64	LOGICAL firstcall
	65	DATA firstcall/.true./
	66	SAVE firstcall
	67
	68	REAL CBRT
	69	EXTERNAL CBRT
	70
	71	INTEGER,SAVE :: i_co2ice=0 ! co2 ice
	72	INTEGER,SAVE :: i_h2oice=0 ! water ice
	73	CHARACTER(LEN=20) :: tracername ! to temporarily store text
	74
[253]	75	! for fixed dust profiles
	76	real topdust, expfactor, zp
	77	REAL taudusttmp(ngridmx) ! Temporary dust opacity used before scaling
[135]	78
[253]	79	! BENJAMIN MODIFS
	80	real CLFtot,CLF1,CLF2
	81	real totcloudfrac(ngridmx)
	82	logical clearsky
	83
	84	! identify tracers
[135]	85	IF (firstcall) THEN
	86
[253]	87	! are these tests of any real use ?
[135]	88	IF(ngrid.NE.ngridmx) THEN
	89	PRINT*,'STOP in aeropacity!'
	90	PRINT*,'problem of dimensions:'
	91	PRINT*,'ngrid =',ngrid
	92	PRINT*,'ngridmx =',ngridmx
	93	STOP
	94	ENDIF
	95
	96	if (nq.gt.nqmx) then
	97	write(,) 'STOP in aeropacity: (nq .gt. nqmx)!'
	98	write(,) 'nq=',nq,' nqmx=',nqmx
	99	stop
	100	endif
	101
	102	do iq=1,nqmx
	103	tracername=noms(iq)
	104	if (tracername.eq."co2_ice") then
	105	i_co2ice=iq
	106	endif
	107	if (tracername.eq."h2o_ice") then
	108	i_h2oice=iq
	109	endif
	110	enddo
	111
	112	write(,) "aeropacity: i_co2ice=",i_co2ice
	113	write(,) " i_h2oice=",i_h2oice
	114
[253]	115	if(watercond.and.(.not.aerofixed))then
	116	if(naerkind.lt.2)then
	117	print*,'Cannot have active H2O clouds with naerkind less than 2!'
	118	call abort
	119	endif
	120	endif
[135]	121
[253]	122	if(dusttau.gt.0.0)then
	123	if(naerkind.lt.3)then
	124	print*,'Cannot have active dust with naerkind less than 3!'
	125	call abort
	126	endif
	127	endif
	128
[135]	129	firstcall=.false.
	130	ENDIF ! of IF (firstcall)
	131
	132
	133	DO iaer = 1, naerkind ! Loop on aerosol kind (NOT tracer)
	134	! ---------------------------------------------------------
[253]	135	aerkind: SELECT CASE (iaer)
[135]	136	!==================================================================
[253]	137	CASE(1) aerkind ! CO2 ice (iaer=1)
[135]	138	!==================================================================
	139
	140	! 1. Initialization
[253]	141	aerosol(:,:,iaer)=0.0
[135]	142
	143	! 2. Opacity calculation
[253]	144	if (aerofixed.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed
	145	do ig=1, ngrid
	146	do l=1,nlayer
	147	aerosol(ig,l,iaer)=1.e-9
	148	end do
	149	!aerosol(ig,12,iaer)=4.0 ! single cloud layer option
	150	end do
	151	else
	152	DO ig=1, ngrid
	153	DO l=1,nlayer-1 ! to stop the rad tran bug
[135]	154
[253]	155	aerosol0 = &
	156	( 0.75 * QREFvis3d(ig,l,iaer) / &
	157	( rho_co2 * reffrad(ig,l,iaer) ) ) * &
	158	( pq(ig,l,i_co2ice) + 1.E-9 ) * &
	159	( pplev(ig,l) - pplev(ig,l+1) ) / g
	160	aerosol0 = max(aerosol0,1.e-9)
	161	aerosol0 = min(aerosol0,L_TAUMAX)
	162	aerosol(ig,l,iaer) = aerosol0
	163	! aerosol(ig,l,iaer) = 0.0
	164
	165	! using cloud fraction
	166	! aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF))
	167	! aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX)
	168
	169
	170	ENDDO
	171	ENDDO
	172	end if
	173
[135]	174	!==================================================================
[253]	175	CASE(2) aerkind ! Water ice / liquid (iaer=2)
[135]	176	!==================================================================
	177
	178	! 1. Initialization
[253]	179	aerosol(:,:,iaer)=0.0
[135]	180
	181	! 2. Opacity calculation
[253]	182	if (aerofixed.or.(i_h2oice.eq.0).or.clearsky) then
[305]	183	do ig=1,ngrid ! to stop the rad tran bug
[253]	184	do l=1,nlayer
	185	aerosol(ig,l,iaer) =1.e-9
	186	end do
	187	end do
[305]	188
	189	! put cloud at cloudlvl
	190	if(kastprof.and.(cloudlvl.ne.0.0))then
	191	ig=1
	192	do l=1,nlayer
	193	if(int(cloudlvl).eq.l)then
	194	!if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
	195	print*,'Inserting cloud at level ',l
	196	!aerosol(ig,l,iaer)=10.0
	197
	198	rho_ice=920.0
	199
	200	! the Kasting approximation
	201	aerosol(ig,l,iaer) = &
	202	( 0.75 * QREFvis3d(ig,l,iaer) / &
	203	( rho_ice * reffrad(ig,l,iaer) ) ) * &
	204	!( pq(ig,l,i_h2oice) + 1.E-9 ) * &
	205	( 4.0e-4 + 1.E-9 ) * &
	206	( pplev(ig,l) - pplev(ig,l+1) ) / g
	207
	208
	209	open(115,file='clouds.out',form='formatted')
	210	write(115,*) l,aerosol(ig,l,iaer)
	211	close(115)
	212
	213	return
	214	endif
	215	end do
	216
	217	call abort
	218	endif
	219
[253]	220	else
	221	do ig=1, ngrid
	222	do l=1,nlayer-1 ! to stop the rad tran bug
[135]	223
[253]	224	if(CLFvarying)then
	225	CLF1 = max(cloudfrac(ig,l),0.01)
	226	CLFtot = max(totcloudfrac(ig),0.01)
	227	CLF2 = CLF1/CLFtot
	228	CLF2 = min(1.,CLF2)
	229	else
	230	CLF1=1.0
	231	CLF2=CLFfixval
	232	endif
	233
	234	aerosol0 = &
	235	( 0.75 * QREFvis3d(ig,l,iaer) / &
	236	( rho_ice * reffrad(ig,l,iaer) ) ) * &
[526]	237	pq(ig,l,i_h2oice) * & !JL I dropped the +1e-9 here to have the same
	238	( pplev(ig,l) - pplev(ig,l+1) ) / g / & ! opacity in the clearsky=true and the
	239	CLF1 ! clear=false/pq=0 case
[135]	240
[253]	241	aerosol(ig,l,iaer) = -log(1 - CLF2 + CLF2*exp(-aerosol0))
	242	! why no division in exponential?
	243	! because its already done in aerosol0
	244
	245	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
	246	aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),aerosol0)
	247
	248	enddo
	249	enddo
	250	end if
[135]	251
[253]	252
[135]	253	!==================================================================
[253]	254	CASE(3) aerkind ! Dust (iaer=3)
	255	!==================================================================
[135]	256
[253]	257	! 1. Initialization
	258	aerosol(:,:,iaer)=0.0
	259
	260	topdust=10.0 ! km
	261
	262	print*,'WARNING, dust is experimental for Early Mars'
	263
	264	! 2. Opacity calculation
	265
[305]	266	do ig=1,ngrid
	267	do l=1,nlayer
[253]	268	zp=700./pplay(ig,l)
[486]	269	aerosol(ig,l,iaer)=(dusttau/700.)*(pplev(ig,l)-pplev(ig,l+1)) &
[253]	270	max( exp(.03(1.-max(zp,1.))) , 1.E-3 ) &
	271	*QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
	272	! takes into account particle variation
	273	enddo
	274	enddo
	275
	276	!==================================================================
	277	END SELECT aerkind
[135]	278	ENDDO ! iaer (loop on aerosol kind)
	279
	280
	281	! --------------------------------------------------------------------------
	282	! Column integrated visible optical depth in each point (used for diagnostic)
	283
[253]	284	tau_col(:)=0.0
[135]	285	do iaer = 1, naerkind
	286	do l=1,nlayer
	287	do ig=1,ngrid
	288	tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
	289	end do
	290	end do
	291	end do
	292
[253]	293	do ig=1, ngrid
	294	do l=1,nlayer
	295	do iaer = 1, naerkind
	296	if(aerosol(ig,l,iaer).gt.1.e3)then
	297	print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
	298	print*,'at ig=',ig,', l=',l,', iaer=',iaer
	299	print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
	300	print*,'reffrad=',reffrad(ig,l,iaer)
	301	endif
	302	end do
	303	end do
	304	end do
	305
	306	do ig=1, ngrid
	307	if(tau_col(ig).gt.1.e3)then
	308	print*,'WARNING: tau_col=',tau_col(ig)
	309	print*,'at ig=',ig
	310	print*,'aerosol=',aerosol(ig,:,:)
	311	print*,'QREFvis3d=',QREFvis3d(ig,:,:)
	312	print*,'reffrad=',reffrad(ig,:,:)
	313	endif
	314	end do
	315
[135]	316	return
	317	end subroutine aeropacity
	318

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

Download in other formats:

Original Format