Context Navigation

aeropacity.F90 @ 735

Last change on this file since 735 was 728, checked in by jleconte, 13 years ago

18/07/2012 == JL

New water cycle scheme:
- largescale now in F90. Robustness increased by i) including evap inside largescale ii) computing the

condensed water amount iteratively

same improvements in moistadj.
Water thermodynamical data and saturation curves centralized in module watercommn_h
- The saturation curves used are now Tetens formula as they are analyticaly inversible (Ts(P)-> Ps(T)).

New saturation curve yields very good agreement with the former one.

Saturation curves are now generalized for arbitrary water amount (not just q<<1)
The old watersat should be removed soon.

The effect of water vapor on total (surface) pressure can be taken into account by setting

mass_redistrib=.true. in callphys.def (routine mass_redistribution inspired from co2_condense in martian
model but with a different scheme as many routines evaporate/condense water vapor).

New cloud and precipitation scheme (JL + BC):
- The default recovery assumption for computing the total cloud fraction has been changed (total random gave too

large cloud fractions). See totalcloudfrac.F90 for details and to change this.

Totalcloudfraction now set the total cloud fraction to the fraction of the

optically thickest cloud and totalcloudfrac is thus called in aeropacity.

Only the total cloud fraction is used to compute optical depth in aeropacity (no more effective

optical depth with exponential formula).

4 precipitation schemes are now available (see rain.F90 for details). The choice can be made using precip_scheme

in callphys.def. Usage of the more physically based model of Boucher et al 95 (precip_scheme=4) is recommended.
default behavior is set to the former "simple scheme" (precip_scheme=1).

See rain.f90 to determine the parameter to be defined in callphys.def as a function of the precipitation scheme used.

Physiq.F90 now written in a matricial (more F90) way.
Radii (H2O and CO2 cloud particles, aerosols, duts, ...) calculations now centralized in module radii_mod.F90

and work with the new aerosol scheme implemented by Laura K. Some inconsistency may remain in callsedim.

Implementation compiled with ifort and pgf90.
gcm.e runs in Earth and Early Mars case with CO2 and H2O cycle + dust.

File size: 13.6 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
[726]	118	if (naerkind.ne.0) then
	119	print*, "If you would like to use aerosols, make sure any old"
	120	print*, "start files are updated in newstart using the option"
	121	print*, "q=0"
	122	write(,) "Aerosols found in aeropacity:"
	123	endif
[253]	124
[726]	125	if (iaero_co2.ne.0) then
	126	print*, 'iaero_co2= ',iaero_co2
	127	endif
	128	if (iaero_h2o.ne.0) then
	129	print*,'iaero_h2o= ',iaero_h2o
	130	endif
	131	if (iaero_dust.ne.0) then
	132	print*,'iaero_dust= ',iaero_dust
	133	endif
	134	if (iaero_h2so4.ne.0) then
	135	print*,'iaero_h2so4= ',iaero_h2so4
	136	endif
	137
[135]	138	firstcall=.false.
	139	ENDIF ! of IF (firstcall)
	140
	141
	142	! ---------------------------------------------------------
	143	!==================================================================
[726]	144	! CO2 ice aerosols
[135]	145	!==================================================================
	146
[728]	147	if (iaero_co2.ne.0) then
[726]	148	iaer=iaero_co2
[135]	149	! 1. Initialization
[728]	150	aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
[135]	151	! 2. Opacity calculation
[726]	152	if (aerofixco2.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed
[728]	153	aerosol(1:ngridmx,1:nlayermx,iaer)=1.e-9
	154	!aerosol(1:ngridmx,12,iaer)=4.0 ! single cloud layer option
[253]	155	else
	156	DO ig=1, ngrid
	157	DO l=1,nlayer-1 ! to stop the rad tran bug
[135]	158
[253]	159	aerosol0 = &
	160	( 0.75 * QREFvis3d(ig,l,iaer) / &
	161	( rho_co2 * reffrad(ig,l,iaer) ) ) * &
	162	( pq(ig,l,i_co2ice) + 1.E-9 ) * &
	163	( pplev(ig,l) - pplev(ig,l+1) ) / g
	164	aerosol0 = max(aerosol0,1.e-9)
	165	aerosol0 = min(aerosol0,L_TAUMAX)
	166	aerosol(ig,l,iaer) = aerosol0
	167	! aerosol(ig,l,iaer) = 0.0
[726]	168	! print*, aerosol(ig,l,iaer)
[253]	169	! using cloud fraction
	170	! aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF))
	171	! aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX)
	172
	173
	174	ENDDO
	175	ENDDO
[726]	176	end if ! if fixed or varying
[728]	177	end if ! if CO2 aerosols
[135]	178	!==================================================================
[726]	179	! Water ice / liquid
[135]	180	!==================================================================
	181
[728]	182	if (iaero_h2o.ne.0) then
[726]	183	iaer=iaero_h2o
[135]	184	! 1. Initialization
[728]	185	aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
[135]	186	! 2. Opacity calculation
[726]	187	if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
[728]	188	aerosol(1:ngridmx,1:nlayermx,iaer) =1.e-9
[305]	189
	190	! put cloud at cloudlvl
	191	if(kastprof.and.(cloudlvl.ne.0.0))then
	192	ig=1
	193	do l=1,nlayer
	194	if(int(cloudlvl).eq.l)then
	195	!if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
	196	print*,'Inserting cloud at level ',l
	197	!aerosol(ig,l,iaer)=10.0
	198
	199	rho_ice=920.0
	200
	201	! the Kasting approximation
	202	aerosol(ig,l,iaer) = &
	203	( 0.75 * QREFvis3d(ig,l,iaer) / &
	204	( rho_ice * reffrad(ig,l,iaer) ) ) * &
	205	!( pq(ig,l,i_h2oice) + 1.E-9 ) * &
	206	( 4.0e-4 + 1.E-9 ) * &
	207	( pplev(ig,l) - pplev(ig,l+1) ) / g
	208
	209
	210	open(115,file='clouds.out',form='formatted')
	211	write(115,*) l,aerosol(ig,l,iaer)
	212	close(115)
	213
	214	return
	215	endif
	216	end do
	217
	218	call abort
	219	endif
	220
[253]	221	else
[728]	222
[253]	223	do ig=1, ngrid
	224	do l=1,nlayer-1 ! to stop the rad tran bug
[135]	225
[728]	226	aerosol(ig,l,iaer) = & !modification by BC
[253]	227	( 0.75 * QREFvis3d(ig,l,iaer) / &
	228	( rho_ice * reffrad(ig,l,iaer) ) ) * &
[716]	229	! pq(ig,l,i_h2oice) * & !JL I dropped the +1e-9 here to have the same
[728]	230	!( pplev(ig,l) - pplev(ig,l+1) ) / g ! opacity in the clearsky=true and the
	231	! clear=false/pq=0 case
[716]	232	( pq(ig,l,i_h2oice) + 1.E-9 ) * & ! Doing this makes the code unstable, so I have restored it (RW)
[728]	233	( pplev(ig,l) - pplev(ig,l+1) ) / g
[135]	234
[253]	235	enddo
	236	enddo
[135]	237
[728]	238	if(CLFvarying)then
	239	call totalcloudfrac(cloudfrac,totcloudfrac,pplev,pq,aerosol(1:ngridmx,1:nlayermx,iaer))
	240	do ig=1, ngrid
	241	do l=1,nlayer-1 ! to stop the rad tran bug
	242	CLFtot = max(totcloudfrac(ig),0.01)
	243	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
	244	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
	245	enddo
	246	enddo
	247	else
	248	do ig=1, ngrid
	249	do l=1,nlayer-1 ! to stop the rad tran bug
	250	CLFtot = CLFfixval
	251	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
	252	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
	253	enddo
	254	enddo
	255	end if!(CLFvarying)
	256	endif !(aerofixed.or.(i_h2oice.eq.0).or.clearsky)
	257
	258	end if ! End if h2o aerosol
	259
[726]	260	!==================================================================
	261	! Dust
	262	!==================================================================
[728]	263	if (iaero_dust.ne.0) then
[726]	264	iaer=iaero_dust
	265	! 1. Initialization
[728]	266	aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
[726]	267
[728]	268	topdust=30.0 ! km (used to be 10.0 km) LK
[726]	269
	270	! 2. Opacity calculation
	271
	272	! expfactor=0.
	273	DO l=1,nlayer-1
	274	DO ig=1,ngrid
	275	! Typical mixing ratio profile
	276
	277	zp=(preff/pplay(ig,l))**(70./topdust)
	278	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
	279
	280	! Vertical scaling function
	281	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
	282	*expfactor
	283
	284
	285	ENDDO
	286	ENDDO
	287
	288	! Rescaling each layer to reproduce the choosen (or assimilated)
	289	! dust extinction opacity at visible reference wavelength, which
	290	! is scaled to the "preff" reference surface pressure available in comvert.h
	291	! and stored in startfi.nc
	292
	293	taudusttmp(1:ngrid)=0.
	294	DO l=1,nlayer
	295	DO ig=1,ngrid
	296	taudusttmp(ig) = taudusttmp(ig) &
	297	+ aerosol(ig,l,iaer)
	298	ENDDO
	299	ENDDO
	300	DO l=1,nlayer-1
	301	DO ig=1,ngrid
	302	aerosol(ig,l,iaer) = max(1E-20, &
	303	dusttau &
	304	* pplev(ig,1) / preff &
	305	* aerosol(ig,l,iaer) &
	306	/ taudusttmp(ig))
	307
	308	ENDDO
	309	ENDDO
[728]	310	end if ! If dust aerosol
[726]	311
[135]	312	!==================================================================
[726]	313	! H2SO4
[253]	314	!==================================================================
[726]	315	! added by LK
	316	if (iaero_h2so4.ne.0) then
[728]	317	iaer=iaero_h2so4
[135]	318
[253]	319	! 1. Initialization
[728]	320	aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
[253]	321
	322
	323	! 2. Opacity calculation
	324
[726]	325	! expfactor=0.
[728]	326	DO l=1,nlayer-1
	327	DO ig=1,ngrid
	328	! Typical mixing ratio profile
[253]	329
[728]	330	zp=(preff/pplay(ig,l))**(70./30) !emulating topdust
	331	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
[135]	332
[726]	333	! Vertical scaling function
[728]	334	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor
[135]	335
[728]	336	ENDDO
	337	ENDDO
	338	tauh2so4tmp(1:ngrid)=0.
	339	DO l=1,nlayer
	340	DO ig=1,ngrid
	341	tauh2so4tmp(ig) = tauh2so4tmp(ig) + aerosol(ig,l,iaer)
	342	ENDDO
	343	ENDDO
	344	DO l=1,nlayer-1
	345	DO ig=1,ngrid
	346	aerosol(ig,l,iaer) = max(1E-20, &
[726]	347	1 &
	348	* pplev(ig,1) / preff &
	349	* aerosol(ig,l,iaer) &
	350	/ tauh2so4tmp(ig))
	351
	352	ENDDO
[728]	353	ENDDO
[726]	354
	355	! 1/700. is assuming a "sulfurtau" of 1
	356	! Sulfur aerosol routine to be improved.
	357	! aerosol0 = &
	358	! ( 0.75 * QREFvis3d(ig,l,iaer) / &
	359	! ( rho_h2so4 * reffrad(ig,l,iaer) ) ) * &
	360	! ( pq(ig,l,i_h2so4) + 1.E-9 ) * &
	361	! ( pplev(ig,l) - pplev(ig,l+1) ) / g
	362	! aerosol0 = max(aerosol0,1.e-9)
	363	! aerosol0 = min(aerosol0,L_TAUMAX)
	364	! aerosol(ig,l,iaer) = aerosol0
	365
	366	! ENDDO
	367	! ENDDO
	368	end if
	369
	370
	371
[135]	372	! --------------------------------------------------------------------------
	373	! Column integrated visible optical depth in each point (used for diagnostic)
	374
[253]	375	tau_col(:)=0.0
[135]	376	do iaer = 1, naerkind
	377	do l=1,nlayer
	378	do ig=1,ngrid
	379	tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
	380	end do
	381	end do
	382	end do
	383
[253]	384	do ig=1, ngrid
	385	do l=1,nlayer
	386	do iaer = 1, naerkind
	387	if(aerosol(ig,l,iaer).gt.1.e3)then
	388	print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
	389	print*,'at ig=',ig,', l=',l,', iaer=',iaer
	390	print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
	391	print*,'reffrad=',reffrad(ig,l,iaer)
	392	endif
	393	end do
	394	end do
	395	end do
	396
	397	do ig=1, ngrid
	398	if(tau_col(ig).gt.1.e3)then
	399	print*,'WARNING: tau_col=',tau_col(ig)
	400	print*,'at ig=',ig
	401	print*,'aerosol=',aerosol(ig,:,:)
	402	print*,'QREFvis3d=',QREFvis3d(ig,:,:)
	403	print*,'reffrad=',reffrad(ig,:,:)
	404	endif
	405	end do
[135]	406	return
	407	end subroutine aeropacity
	408

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

Download in other formats:

Original Format