Context Navigation

aeropacity.F90 @ 255

Last change on this file since 255 was 253, checked in by emillour, 14 years ago

Generic GCM

Massive update to version 0.7

EM+RW

File size: 9.2 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
	183	do ig=1, ngrid ! to stop the rad tran bug
	184	do l=1,nlayer
	185	aerosol(ig,l,iaer) =1.e-9
	186	end do
	187	end do
	188	else
	189	do ig=1, ngrid
	190	do l=1,nlayer-1 ! to stop the rad tran bug
[135]	191
[253]	192	if(CLFvarying)then
	193	CLF1 = max(cloudfrac(ig,l),0.01)
	194	CLFtot = max(totcloudfrac(ig),0.01)
	195	CLF2 = CLF1/CLFtot
	196	CLF2 = min(1.,CLF2)
	197	else
	198	CLF1=1.0
	199	CLF2=CLFfixval
	200	endif
	201
	202	aerosol0 = &
	203	( 0.75 * QREFvis3d(ig,l,iaer) / &
	204	( rho_ice * reffrad(ig,l,iaer) ) ) * &
	205	( pq(ig,l,i_h2oice) + 1.E-9 ) * &
	206	( pplev(ig,l) - pplev(ig,l+1) ) / g / &
	207	CLF1
[135]	208
[253]	209	aerosol(ig,l,iaer) = -log(1 - CLF2 + CLF2*exp(-aerosol0))
	210	! why no division in exponential?
	211	! because its already done in aerosol0
	212
	213	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
	214	aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),aerosol0)
	215
	216	enddo
	217	enddo
	218	end if
[135]	219
[253]	220
[135]	221	!==================================================================
[253]	222	CASE(3) aerkind ! Dust (iaer=3)
	223	!==================================================================
[135]	224
[253]	225	! 1. Initialization
	226	aerosol(:,:,iaer)=0.0
	227
	228	topdust=10.0 ! km
	229
	230	print*,'WARNING, dust is experimental for Early Mars'
	231
	232	! 2. Opacity calculation
	233
	234	do l=1,nlayer
	235	do ig=1,ngrid-1 ! to stop the rad tran bug
	236	zp=700./pplay(ig,l)
	237	aerosol(ig,l,1)=(dusttau/700.)*(pplev(ig,l)-pplev(ig,l+1)) &
	238	max( exp(.03(1.-max(zp,1.))) , 1.E-3 ) &
	239	*QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
	240	! takes into account particle variation
	241	enddo
	242	enddo
	243
	244	!==================================================================
	245	END SELECT aerkind
[135]	246	ENDDO ! iaer (loop on aerosol kind)
	247
	248
	249	! --------------------------------------------------------------------------
	250	! Column integrated visible optical depth in each point (used for diagnostic)
	251
[253]	252	tau_col(:)=0.0
[135]	253	do iaer = 1, naerkind
	254	do l=1,nlayer
	255	do ig=1,ngrid
	256	tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
	257	end do
	258	end do
	259	end do
	260
[253]	261	do ig=1, ngrid
	262	do l=1,nlayer
	263	do iaer = 1, naerkind
	264	if(aerosol(ig,l,iaer).gt.1.e3)then
	265	print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
	266	print*,'at ig=',ig,', l=',l,', iaer=',iaer
	267	print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
	268	print*,'reffrad=',reffrad(ig,l,iaer)
	269	endif
	270	end do
	271	end do
	272	end do
	273
	274	do ig=1, ngrid
	275	if(tau_col(ig).gt.1.e3)then
	276	print*,'WARNING: tau_col=',tau_col(ig)
	277	print*,'at ig=',ig
	278	print*,'aerosol=',aerosol(ig,:,:)
	279	print*,'QREFvis3d=',QREFvis3d(ig,:,:)
	280	print*,'reffrad=',reffrad(ig,:,:)
	281	endif
	282	end do
	283
[135]	284	return
	285	end subroutine aeropacity
	286

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

Download in other formats:

Original Format