Context Navigation

aeropacity.F90 @ 778

Last change on this file since 778 was 741, checked in by jleconte, 13 years ago
Model now works without aerosols at all. You still need to compile it with -s 1 option. LK12
File size: 13.8 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
[726]	4	use radinc_h, only : L_TAUMAX,naerkind
	5	use aerosol_mod
[135]	6
	7	implicit none
	8
	9	!==================================================================
	10	!
	11	! Purpose
	12	! -------
	13	! Compute aerosol optical depth in each gridbox.
	14	!
	15	! Authors
	16	! -------
	17	! F. Forget
	18	! F. Montmessin (water ice scheme)
	19	! update J.-B. Madeleine (2008)
	20	! dust removal, simplification by Robin Wordsworth (2009)
	21	!
	22	! Input
	23	! -----
	24	! ngrid Number of horizontal gridpoints
	25	! nlayer Number of layers
	26	! nq Number of tracers
	27	! pplev Pressure (Pa) at each layer boundary
	28	! pq Aerosol mixing ratio
	29	! reffrad(ngrid,nlayer,naerkind) Aerosol effective radius
	30	! QREFvis3d(ngridmx,nlayermx,naerkind) \ 3d extinction coefficients
	31	! QREFir3d(ngridmx,nlayermx,naerkind) / at reference wavelengths
	32	!
	33	! Output
	34	! ------
	35	! aerosol Aerosol optical depth in layer l, grid point ig
	36	! tau_col Total column optical depth at grid point ig
	37	!
	38	!=======================================================================
	39
	40	#include "dimensions.h"
	41	#include "dimphys.h"
	42	#include "callkeys.h"
	43	#include "comcstfi.h"
	44	#include "comgeomfi.h"
	45	#include "tracer.h"
[726]	46	#include "comvert.h"
[135]	47
	48	INTEGER ngrid,nlayer,nq
	49
[253]	50	REAL pplay(ngrid,nlayer)
[135]	51	REAL pplev(ngrid,nlayer+1)
	52	REAL pq(ngrid,nlayer,nq)
	53	REAL aerosol(ngrid,nlayer,naerkind)
	54	REAL reffrad(ngrid,nlayer,naerkind)
	55	REAL QREFvis3d(ngridmx,nlayermx,naerkind)
	56	REAL QREFir3d(ngridmx,nlayermx,naerkind)
	57
[253]	58	REAL tau_col(ngrid)
	59	! REAL tauref(ngrid), tau_col(ngrid)
[135]	60
[253]	61	real cloudfrac(ngridmx,nlayermx)
	62	real aerosol0
	63
[135]	64	INTEGER l,ig,iq,iaer
	65
	66	LOGICAL firstcall
	67	DATA firstcall/.true./
	68	SAVE firstcall
	69
	70	REAL CBRT
	71	EXTERNAL CBRT
	72
	73	INTEGER,SAVE :: i_co2ice=0 ! co2 ice
	74	INTEGER,SAVE :: i_h2oice=0 ! water ice
	75	CHARACTER(LEN=20) :: tracername ! to temporarily store text
	76
[253]	77	! for fixed dust profiles
	78	real topdust, expfactor, zp
	79	REAL taudusttmp(ngridmx) ! Temporary dust opacity used before scaling
[726]	80	REAL tauh2so4tmp(ngridmx) ! Temporary h2so4 opacity used before scaling
[253]	81	! BENJAMIN MODIFS
[728]	82	real CLFtot
[253]	83	real totcloudfrac(ngridmx)
	84	logical clearsky
	85
	86	! identify tracers
[135]	87	IF (firstcall) THEN
	88
[253]	89	! are these tests of any real use ?
[135]	90	IF(ngrid.NE.ngridmx) THEN
	91	PRINT*,'STOP in aeropacity!'
	92	PRINT*,'problem of dimensions:'
	93	PRINT*,'ngrid =',ngrid
	94	PRINT*,'ngridmx =',ngridmx
	95	STOP
	96	ENDIF
	97
	98	if (nq.gt.nqmx) then
	99	write(,) 'STOP in aeropacity: (nq .gt. nqmx)!'
	100	write(,) 'nq=',nq,' nqmx=',nqmx
	101	stop
	102	endif
	103
[726]	104	write(,) "Tracers found in aeropacity:"
[135]	105	do iq=1,nqmx
	106	tracername=noms(iq)
	107	if (tracername.eq."co2_ice") then
	108	i_co2ice=iq
[726]	109	write(,) "i_co2ice=",i_co2ice
	110
[135]	111	endif
	112	if (tracername.eq."h2o_ice") then
	113	i_h2oice=iq
[726]	114	write(,) "i_h2oice=",i_h2oice
[135]	115	endif
	116	enddo
	117
[741]	118	if (noaero) then
	119	print*, "No active aerosols found in aeropacity"
	120	else
[726]	121	print*, "If you would like to use aerosols, make sure any old"
	122	print*, "start files are updated in newstart using the option"
	123	print*, "q=0"
[741]	124	write(,) "Active aerosols found in aeropacity:"
[726]	125	endif
[253]	126
[741]	127	if ((iaero_co2.ne.0).and.(.not.noaero)) then
[726]	128	print*, 'iaero_co2= ',iaero_co2
	129	endif
	130	if (iaero_h2o.ne.0) then
	131	print*,'iaero_h2o= ',iaero_h2o
	132	endif
	133	if (iaero_dust.ne.0) then
	134	print*,'iaero_dust= ',iaero_dust
	135	endif
	136	if (iaero_h2so4.ne.0) then
	137	print*,'iaero_h2so4= ',iaero_h2so4
	138	endif
	139
[135]	140	firstcall=.false.
	141	ENDIF ! of IF (firstcall)
	142
	143
	144	! ---------------------------------------------------------
	145	!==================================================================
[726]	146	! CO2 ice aerosols
[135]	147	!==================================================================
	148
[728]	149	if (iaero_co2.ne.0) then
[726]	150	iaer=iaero_co2
[135]	151	! 1. Initialization
[728]	152	aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
[135]	153	! 2. Opacity calculation
[741]	154	if (noaero) then ! aerosol set to a very low value
	155	aerosol(1:ngridmx,1:nlayermx,iaer)=1.e-9
	156	elseif (aerofixco2.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed
[728]	157	aerosol(1:ngridmx,1:nlayermx,iaer)=1.e-9
	158	!aerosol(1:ngridmx,12,iaer)=4.0 ! single cloud layer option
[253]	159	else
	160	DO ig=1, ngrid
	161	DO l=1,nlayer-1 ! to stop the rad tran bug
[135]	162
[253]	163	aerosol0 = &
	164	( 0.75 * QREFvis3d(ig,l,iaer) / &
	165	( rho_co2 * reffrad(ig,l,iaer) ) ) * &
	166	( pq(ig,l,i_co2ice) + 1.E-9 ) * &
	167	( pplev(ig,l) - pplev(ig,l+1) ) / g
	168	aerosol0 = max(aerosol0,1.e-9)
	169	aerosol0 = min(aerosol0,L_TAUMAX)
	170	aerosol(ig,l,iaer) = aerosol0
	171	! aerosol(ig,l,iaer) = 0.0
[726]	172	! print*, aerosol(ig,l,iaer)
[253]	173	! using cloud fraction
	174	! aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF))
	175	! aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX)
	176
	177
	178	ENDDO
	179	ENDDO
[726]	180	end if ! if fixed or varying
[728]	181	end if ! if CO2 aerosols
[135]	182	!==================================================================
[726]	183	! Water ice / liquid
[135]	184	!==================================================================
	185
[728]	186	if (iaero_h2o.ne.0) then
[726]	187	iaer=iaero_h2o
[135]	188	! 1. Initialization
[728]	189	aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
[135]	190	! 2. Opacity calculation
[726]	191	if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
[728]	192	aerosol(1:ngridmx,1:nlayermx,iaer) =1.e-9
[305]	193
	194	! put cloud at cloudlvl
	195	if(kastprof.and.(cloudlvl.ne.0.0))then
	196	ig=1
	197	do l=1,nlayer
	198	if(int(cloudlvl).eq.l)then
	199	!if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
	200	print*,'Inserting cloud at level ',l
	201	!aerosol(ig,l,iaer)=10.0
	202
	203	rho_ice=920.0
	204
	205	! the Kasting approximation
	206	aerosol(ig,l,iaer) = &
	207	( 0.75 * QREFvis3d(ig,l,iaer) / &
	208	( rho_ice * reffrad(ig,l,iaer) ) ) * &
	209	!( pq(ig,l,i_h2oice) + 1.E-9 ) * &
	210	( 4.0e-4 + 1.E-9 ) * &
	211	( pplev(ig,l) - pplev(ig,l+1) ) / g
	212
	213
	214	open(115,file='clouds.out',form='formatted')
	215	write(115,*) l,aerosol(ig,l,iaer)
	216	close(115)
	217
	218	return
	219	endif
	220	end do
	221
	222	call abort
	223	endif
	224
[253]	225	else
[728]	226
[253]	227	do ig=1, ngrid
	228	do l=1,nlayer-1 ! to stop the rad tran bug
[135]	229
[728]	230	aerosol(ig,l,iaer) = & !modification by BC
[253]	231	( 0.75 * QREFvis3d(ig,l,iaer) / &
	232	( rho_ice * reffrad(ig,l,iaer) ) ) * &
[716]	233	! pq(ig,l,i_h2oice) * & !JL I dropped the +1e-9 here to have the same
[728]	234	!( pplev(ig,l) - pplev(ig,l+1) ) / g ! opacity in the clearsky=true and the
	235	! clear=false/pq=0 case
[716]	236	( pq(ig,l,i_h2oice) + 1.E-9 ) * & ! Doing this makes the code unstable, so I have restored it (RW)
[728]	237	( pplev(ig,l) - pplev(ig,l+1) ) / g
[135]	238
[253]	239	enddo
	240	enddo
[135]	241
[728]	242	if(CLFvarying)then
	243	call totalcloudfrac(cloudfrac,totcloudfrac,pplev,pq,aerosol(1:ngridmx,1:nlayermx,iaer))
	244	do ig=1, ngrid
	245	do l=1,nlayer-1 ! to stop the rad tran bug
	246	CLFtot = max(totcloudfrac(ig),0.01)
	247	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
	248	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
	249	enddo
	250	enddo
	251	else
	252	do ig=1, ngrid
	253	do l=1,nlayer-1 ! to stop the rad tran bug
	254	CLFtot = CLFfixval
	255	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
	256	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
	257	enddo
	258	enddo
	259	end if!(CLFvarying)
	260	endif !(aerofixed.or.(i_h2oice.eq.0).or.clearsky)
	261
	262	end if ! End if h2o aerosol
	263
[726]	264	!==================================================================
	265	! Dust
	266	!==================================================================
[728]	267	if (iaero_dust.ne.0) then
[726]	268	iaer=iaero_dust
	269	! 1. Initialization
[728]	270	aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
[726]	271
[728]	272	topdust=30.0 ! km (used to be 10.0 km) LK
[726]	273
	274	! 2. Opacity calculation
	275
	276	! expfactor=0.
	277	DO l=1,nlayer-1
	278	DO ig=1,ngrid
	279	! Typical mixing ratio profile
	280
	281	zp=(preff/pplay(ig,l))**(70./topdust)
	282	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
	283
	284	! Vertical scaling function
	285	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
	286	*expfactor
	287
	288
	289	ENDDO
	290	ENDDO
	291
	292	! Rescaling each layer to reproduce the choosen (or assimilated)
	293	! dust extinction opacity at visible reference wavelength, which
	294	! is scaled to the "preff" reference surface pressure available in comvert.h
	295	! and stored in startfi.nc
	296
	297	taudusttmp(1:ngrid)=0.
	298	DO l=1,nlayer
	299	DO ig=1,ngrid
	300	taudusttmp(ig) = taudusttmp(ig) &
	301	+ aerosol(ig,l,iaer)
	302	ENDDO
	303	ENDDO
	304	DO l=1,nlayer-1
	305	DO ig=1,ngrid
	306	aerosol(ig,l,iaer) = max(1E-20, &
	307	dusttau &
	308	* pplev(ig,1) / preff &
	309	* aerosol(ig,l,iaer) &
	310	/ taudusttmp(ig))
	311
	312	ENDDO
	313	ENDDO
[728]	314	end if ! If dust aerosol
[726]	315
[135]	316	!==================================================================
[726]	317	! H2SO4
[253]	318	!==================================================================
[726]	319	! added by LK
	320	if (iaero_h2so4.ne.0) then
[728]	321	iaer=iaero_h2so4
[135]	322
[253]	323	! 1. Initialization
[728]	324	aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
[253]	325
	326
	327	! 2. Opacity calculation
	328
[726]	329	! expfactor=0.
[728]	330	DO l=1,nlayer-1
	331	DO ig=1,ngrid
	332	! Typical mixing ratio profile
[253]	333
[728]	334	zp=(preff/pplay(ig,l))**(70./30) !emulating topdust
	335	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
[135]	336
[726]	337	! Vertical scaling function
[728]	338	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor
[135]	339
[728]	340	ENDDO
	341	ENDDO
	342	tauh2so4tmp(1:ngrid)=0.
	343	DO l=1,nlayer
	344	DO ig=1,ngrid
	345	tauh2so4tmp(ig) = tauh2so4tmp(ig) + aerosol(ig,l,iaer)
	346	ENDDO
	347	ENDDO
	348	DO l=1,nlayer-1
	349	DO ig=1,ngrid
	350	aerosol(ig,l,iaer) = max(1E-20, &
[726]	351	1 &
	352	* pplev(ig,1) / preff &
	353	* aerosol(ig,l,iaer) &
	354	/ tauh2so4tmp(ig))
	355
	356	ENDDO
[728]	357	ENDDO
[726]	358
	359	! 1/700. is assuming a "sulfurtau" of 1
	360	! Sulfur aerosol routine to be improved.
	361	! aerosol0 = &
	362	! ( 0.75 * QREFvis3d(ig,l,iaer) / &
	363	! ( rho_h2so4 * reffrad(ig,l,iaer) ) ) * &
	364	! ( pq(ig,l,i_h2so4) + 1.E-9 ) * &
	365	! ( pplev(ig,l) - pplev(ig,l+1) ) / g
	366	! aerosol0 = max(aerosol0,1.e-9)
	367	! aerosol0 = min(aerosol0,L_TAUMAX)
	368	! aerosol(ig,l,iaer) = aerosol0
	369
	370	! ENDDO
	371	! ENDDO
	372	end if
	373
	374
	375
[135]	376	! --------------------------------------------------------------------------
	377	! Column integrated visible optical depth in each point (used for diagnostic)
	378
[253]	379	tau_col(:)=0.0
[135]	380	do iaer = 1, naerkind
	381	do l=1,nlayer
	382	do ig=1,ngrid
	383	tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
	384	end do
	385	end do
	386	end do
	387
[253]	388	do ig=1, ngrid
	389	do l=1,nlayer
	390	do iaer = 1, naerkind
	391	if(aerosol(ig,l,iaer).gt.1.e3)then
	392	print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
	393	print*,'at ig=',ig,', l=',l,', iaer=',iaer
	394	print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
	395	print*,'reffrad=',reffrad(ig,l,iaer)
	396	endif
	397	end do
	398	end do
	399	end do
	400
	401	do ig=1, ngrid
	402	if(tau_col(ig).gt.1.e3)then
	403	print*,'WARNING: tau_col=',tau_col(ig)
	404	print*,'at ig=',ig
	405	print*,'aerosol=',aerosol(ig,:,:)
	406	print*,'QREFvis3d=',QREFvis3d(ig,:,:)
	407	print*,'reffrad=',reffrad(ig,:,:)
	408	endif
	409	end do
[135]	410	return
	411	end subroutine aeropacity
	412

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

Download in other formats:

Original Format