Context Navigation

aeropacity_mod.F @ 3026

Last change on this file since 3026 was 2921, checked in by emillour, 22 months ago
Mars PCM: Add the possibility to use MY36 dust scenario and MY36 EUV inputs EM
File size: 39.0 KB

Rev	Line
[1711]	1	MODULE aeropacity_mod
	2
	3	IMPLICIT NONE
	4
[2409]	5	INTEGER :: iddist ! flag for vertical dust ditribution type (when imposed)
	6	! 0: Pollack90, 1: top set by "topdustref"
	7	! 2: Viking scenario; =3 MGS scenario
	8	REAL :: topdustref ! Dust top altitude (km); only matters only if iddist=1)
[1711]	9	CONTAINS
	10
[38]	11	SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
[2643]	12	& pq,pt,tauscaling,dust_rad_adjust,IRtoVIScoef,tau_pref_scenario,
	13	& tau_pref_gcm,tau,taucloudtes,aerosol,dsodust,reffrad,
	14	& QREFvis3d,QREFir3d,omegaREFir3d,
	15	& totstormfract,clearatm,dsords,dsotop,
	16	& nohmons,
	17	& clearsky,totcloudfrac)
[1974]	18
[2304]	19	use ioipsl_getin_p_mod, only: getin_p
[1036]	20	use tracer_mod, only: noms, igcm_h2o_ice, igcm_dust_mass,
[1224]	21	& igcm_dust_submicron, rho_dust, rho_ice,
[2199]	22	& nqdust, igcm_stormdust_mass,
[2447]	23	& igcm_topdust_mass, igcm_co2_ice
[1543]	24	use geometry_mod, only: latitude ! grid point latitudes (rad)
[1541]	25	use comgeomfi_h, only: sinlat ! sines of grid point latitudes
[1375]	26	#ifdef DUSTSTORM
[1543]	27	use geometry_mod, only: longitude
[1375]	28	use tracer_mod, only: r3n_q, ref_r0, igcm_dust_number
	29	#endif
[2304]	30	use comcstfi_h, only: g, pi
[1246]	31	use dimradmars_mod, only: naerkind, name_iaer,
	32	& iaerdust,tauvis,
	33	& iaer_dust_conrath,iaer_dust_doubleq,
[1974]	34	& iaer_dust_submicron,iaer_h2o_ice,
[2199]	35	& iaer_stormdust_doubleq,
	36	& iaer_topdust_doubleq
[2643]	37	use dust_param_mod, only: odpref, freedust,
	38	& reff_driven_IRtoVIS_scenario
[2417]	39	use dust_scaling_mod, only: compute_dustscaling
[2494]	40	use density_co2_ice_mod, only: density_co2_ice
[2628]	41	use surfdat_h,only: alpha_hmons,contains_mons
[2643]	42	use read_dust_scenario_mod, only: read_dust_scenario
[2628]	43
[38]	44	IMPLICIT NONE
	45	c=======================================================================
	46	c subject:
	47	c --------
	48	c Computing aerosol optical depth in each gridbox.
	49	c
	50	c author: F.Forget
	51	c ------
	52	c update F. Montmessin (water ice scheme)
	53	c and S. Lebonnois (12/06/2003) compatibility dust/ice/chemistry
	54	c update J.-B. Madeleine 2008-2009:
	55	c - added 3D scattering by aerosols;
	56	c - dustopacity transferred from physiq.F to callradite.F,
	57	c and renamed into aeropacity.F;
[607]	58	c update E. Millour, march 2012:
	59	c - reference pressure is now set to 610Pa (not 700Pa)
[38]	60	c
	61	c=======================================================================
[1974]	62	include "callkeys.h"
[38]	63
	64	c-----------------------------------------------------------------------
	65	c
	66	c Declarations :
	67	c --------------
	68	c
	69	c Input/Output
	70	c ------------
[2415]	71	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
	72	INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
	73	INTEGER,INTENT(IN) :: nq ! number of tracers
	74	REAL,INTENT(IN) :: ls ! Solar Longitude (rad)
	75	REAL,INTENT(IN) :: zday ! date (in martian sols) since Ls=0
	76	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! pressure (Pa) in the middle of
	77	! each atmospheric layer
	78	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! pressure (Pa) at the boundaries
	79	! of the atmospheric layers
	80	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers
[2494]	81	REAL,INTENT(IN) :: pt(ngrid,nlayer) !temperature
[2415]	82	REAL,INTENT(OUT) :: tau_pref_scenario(ngrid) ! prescribed dust column
	83	! visible opacity at odpref from scenario
	84	REAL,INTENT(OUT) :: tau_pref_gcm(ngrid) ! computed dust column
	85	! visible opacity at odpref in the GCM
	86	REAL,INTENT(OUT) :: tau(ngrid,naerkind) ! column total visible
	87	! optical depth of each aerosol
	88	REAL,INTENT(OUT) :: taucloudtes(ngrid)! Water ice cloud opacity at
	89	! infrared reference wavelength using
	90	! Qabs instead of Qext
	91	! (for direct comparison with TES)
	92	REAL, INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! optical
[2643]	93	! depth of each aerosol in each layer
[2415]	94	REAL, INTENT(OUT) :: dsodust(ngrid,nlayer) ! density scaled opacity
	95	! of (background) dust
[2413]	96	REAL, INTENT(OUT) :: dsords(ngrid,nlayer) !dso of stormdust
	97	REAL, INTENT(OUT) :: dsotop(ngrid,nlayer) !dso of topdust
[2415]	98	REAL, INTENT(INOUT) :: reffrad(ngrid,nlayer,naerkind) ! effective radius
	99	! of the aerosols in the grid boxes
	100	REAL, INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind) ! 3D extinction
	101	! coefficients (in the visible) of aerosols
	102	REAL, INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind) ! 3D extinction
	103	! coefficients (in the infra-red) of aerosols
	104	REAL, INTENT(IN) :: omegaREFir3d(ngrid,nlayer,naerkind) ! at the
	105	! reference wavelengths
	106	LOGICAL, INTENT(IN) :: clearatm ! true to compute RT without stormdust
	107	! and false to compute RT in rocket dust storms
	108	REAL, INTENT(IN) :: totstormfract(ngrid) ! mesh fraction with a rocket
	109	! dust storm
[2628]	110	LOGICAL, INTENT(IN) :: nohmons ! true to compute RT without topdust,
	111	! false to compute RT in the topdust
[2417]	112	REAL,INTENT(OUT) :: tauscaling(ngrid) ! Scaling factor for qdust and Ndust
[2634]	113	REAL,INTENT(INOUT) :: dust_rad_adjust(ngrid) ! Radiative adjustment
[2643]	114	! factor for dust
	115	REAL,INTENT(INOUT) :: IRtoVIScoef(ngrid) ! conversion coefficient to apply on
	116	! scenario absorption IR (9.3um) CDOD
	117	! = tau_pref_gcm_VIS / tau_pref_gcm_IR
[2415]	118	REAL,INTENT(IN) :: totcloudfrac(ngrid) ! total water ice cloud fraction
	119	LOGICAL,INTENT(IN) :: clearsky ! true to compute RT without water ice clouds
	120	! false to compute RT with clouds (total or sub-grid clouds)
[38]	121	c
	122	c Local variables :
	123	c -----------------
[1974]	124	REAL CLFtot ! total cloud fraction
	125	real expfactor
[38]	126	INTEGER l,ig,iq,i,j
	127	INTEGER iaer ! Aerosol index
[1047]	128	real topdust(ngrid)
[38]	129	real zlsconst, zp
	130	real taueq,tauS,tauN
	131	c Mean Qext(vis)/Qext(ir) profile
[1047]	132	real msolsir(nlayer,naerkind)
[38]	133	c Mean Qext(ir)/Qabs(ir) profile
[1047]	134	real mqextsqabs(nlayer,naerkind)
[38]	135	c Variables used when multiple particle sizes are used
	136	c for dust or water ice particles in the radiative transfer
	137	c (see callradite.F for more information).
[1047]	138	REAL taucloudvis(ngrid)! Cloud opacity at visible
[38]	139	! reference wavelength
[1224]	140	REAL topdust0(ngrid)
[2643]	141
	142	REAL aerosol_IRabs(ngrid,nlayer)
	143	REAL taudust_IRabs(ngrid)
	144	REAL taudust_VISext(ngrid)
[83]	145
[2494]	146	! -- CO2 clouds
	147	real CLFtotco2
	148	real taucloudco2vis(ngrid)
	149	real taucloudco2tes(ngrid)
	150	real totcloudco2frac(ngrid) ! a mettre en (in) [CM]
	151	double precision :: rho_ice_co2
[1375]	152	#ifdef DUSTSTORM
	153	!! Local dust storms
	154	logical localstorm ! =true to create a local dust storm
	155	real taulocref,ztoploc,radloc,lonloc,latloc ! local dust storm parameters
	156	real reffstorm, yeah
	157	REAL ray(ngrid) ! distance from dust storm center
	158	REAL tauuser(ngrid) ! opacity perturbation due to dust storm
	159	REAL more_dust(ngrid,nlayer,2) ! Mass mixing ratio perturbation due to the dust storm
	160	REAL int_factor(ngrid) ! useful factor to compute mmr perturbation
	161	real l_top ! layer of the storm's top
	162	REAL zalt(ngrid, nlayer) ! useful factor to compute l_top
	163	#endif
	164
[38]	165	c local saved variables
	166	c ---------------------
	167
	168	c Level under which the dust mixing ratio is held constant
	169	c when computing the dust opacity in each layer
	170	c (this applies when doubleq and active are true)
[1208]	171	INTEGER, PARAMETER :: cstdustlevel0 = 7
	172	INTEGER, SAVE :: cstdustlevel
[38]	173
[607]	174	LOGICAL,SAVE :: firstcall=.true.
[38]	175
	176	! indexes of water ice and dust tracers:
	177	INTEGER,SAVE :: i_ice=0 ! water ice
	178	CHARACTER(LEN=20) :: txt ! to temporarly store text
	179	CHARACTER(LEN=1) :: txt2 ! to temporarly store text
[2447]	180	! indexes of co2 ice :
	181	INTEGER,SAVE :: i_co2ice=0 ! co2 ice
[38]	182	! indexes of dust scatterers:
	183	INTEGER,SAVE :: naerdust ! number of dust scatterers
	184
[2584]	185	!$OMP THREADPRIVATE(cstdustlevel,firstcall,i_ice,
	186	!$OMP& i_co2ice,naerdust)
	187
[2252]	188	! initializations
[38]	189	tau(1:ngrid,1:naerkind)=0
	190
	191	! identify tracers
	192
[1775]	193	!! AS: firstcall OK absolute
[38]	194	IF (firstcall) THEN
	195	! identify scatterers that are dust
	196	naerdust=0
[2494]	197	iaerdust(1:naerkind) = 0
	198	nqdust(1:nq) = 0
[38]	199	DO iaer=1,naerkind
	200	txt=name_iaer(iaer)
[1974]	201	! CW17: choice tauscaling for stormdust or not
[2199]	202	IF ((txt(1:4).eq."dust").OR.(txt(1:5).eq."storm")
	203	& .OR.(txt(1:3).eq."top")) THEN !MV19: topdust tracer
[38]	204	naerdust=naerdust+1
	205	iaerdust(naerdust)=iaer
	206	ENDIF
	207	ENDDO
	208	! identify tracers which are dust
	209	i=0
	210	DO iq=1,nq
	211	txt=noms(iq)
	212	IF (txt(1:4).eq."dust") THEN
	213	i=i+1
	214	nqdust(i)=iq
	215	ENDIF
	216	ENDDO
	217	IF (water.AND.activice) THEN
	218	i_ice=igcm_h2o_ice
	219	write(,) "aeropacity: i_ice=",i_ice
	220	ENDIF
	221
[2447]	222	IF (co2clouds.AND.activeco2ice) THEN
	223	i_co2ice=igcm_co2_ice
	224	write(,) "aeropacity: i_co2ice =",i_co2ice
	225	ENDIF
	226
[38]	227	c typical profile of solsir and (1-w)^(-1):
[1775]	228	c --- purely for diagnostics and printing
[38]	229	msolsir(1:nlayer,1:naerkind)=0
	230	mqextsqabs(1:nlayer,1:naerkind)=0
[222]	231	WRITE(,) "Typical profiles of Qext(vis)/Qext(IR)"
	232	WRITE(,) " and Qext(IR)/Qabs(IR):"
[38]	233	DO iaer = 1, naerkind ! Loop on aerosol kind
	234	WRITE(,) "Aerosol # ",iaer
	235	DO l=1,nlayer
[1047]	236	DO ig=1,ngrid
[38]	237	msolsir(l,iaer)=msolsir(l,iaer)+
	238	& QREFvis3d(ig,l,iaer)/
	239	& QREFir3d(ig,l,iaer)
	240	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)+
	241	& (1.E0-omegaREFir3d(ig,l,iaer))**(-1)
	242	ENDDO
[1047]	243	msolsir(l,iaer)=msolsir(l,iaer)/REAL(ngrid)
	244	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)/REAL(ngrid)
[38]	245	ENDDO
	246	WRITE(,) "solsir: ",msolsir(:,iaer)
	247	WRITE(,) "Qext/Qabs(IR): ",mqextsqabs(:,iaer)
	248	ENDDO
	249
	250	! load value of tauvis from callphys.def (if given there,
	251	! otherwise default value read from starfi.nc file will be used)
[2304]	252	call getin_p("tauvis",tauvis)
[38]	253
[1974]	254	IF (freedust.or.rdstorm) THEN ! if rdstorm no need to held opacity constant at the first levels
[1208]	255	cstdustlevel = 1
	256	ELSE
[1974]	257	cstdustlevel = cstdustlevel0 !Opacity in the first levels is held constant to
	258	!avoid unrealistic values due to constant lifting
[1208]	259	ENDIF
	260
[1375]	261	#ifndef DUSTSTORM
[38]	262	firstcall=.false.
[1375]	263	#endif
[38]	264
[2252]	265	END IF ! end of if firstcall
[38]	266
[2415]	267	! 1. Get prescribed tau_pref_scenario, Dust column optical depth at "odpref" Pa
	268	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
[1088]	269
[2415]	270	IF(iaervar.eq.1) THEN
[1047]	271	do ig=1, ngrid
[2415]	272	tau_pref_scenario(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
[38]	273	! or read in starfi
	274	end do
	275	ELSE IF (iaervar.eq.2) THEN ! << "Viking" Scenario>>
	276
[2415]	277	tau_pref_scenario(1) = 0.7+.3cos(ls+80.pi/180.) ! like seen by VL1
[38]	278	do ig=2,ngrid
[2415]	279	tau_pref_scenario(ig) = tau_pref_scenario(1)
[38]	280	end do
	281
	282	ELSE IF (iaervar.eq.3) THEN ! << "MGS" scenario >>
	283
	284	taueq= 0.2 +(0.5-0.2) (cos(0.5(ls-4.363)))**14
	285	tauS= 0.1 +(0.5-0.1) (cos(0.5(ls-4.363)))**14
	286	tauN = 0.1
[1047]	287	do ig=1,ngrid
[1541]	288	if (latitude(ig).ge.0) then
[1047]	289	! Northern hemisphere
[2415]	290	tau_pref_scenario(ig)= tauN +
[1541]	291	& (taueq-tauN)0.5(1+tanh((45-latitude(ig)180./pi)6/60))
[1047]	292	else
	293	! Southern hemisphere
[2415]	294	tau_pref_scenario(ig)= tauS +
[1541]	295	& (taueq-tauS)0.5(1+tanh((45+latitude(ig)180./pi)6/60))
[1047]	296	endif
	297	enddo ! of do ig=1,ngrid
[38]	298	ELSE IF (iaervar.eq.5) THEN ! << Escalier Scenario>>
[2415]	299	tau_pref_scenario(1) = 2.5
[38]	300	if ((ls.ge.30.pi/180.).and.(ls.le.150.pi/180.))
[2415]	301	& tau_pref_scenario(1) = .2
[38]	302
	303	do ig=2,ngrid
[2415]	304	tau_pref_scenario(ig) = tau_pref_scenario(1)
[38]	305	end do
[2643]	306	!!!!!!!!!!!!!!!!!!!!!!!!!!!
	307	! NB: here, IRtoVIScoef=2.6
	308	! ( useful to be here only if iddist=0 (Pollack90 vertical distribution) )
[1278]	309	ELSE IF ((iaervar.ge.6).and.(iaervar.le.8)) THEN
	310	! clim, cold or warm synthetic scenarios
[2415]	311	call read_dust_scenario(ngrid,nlayer,zday,pplev,
[2643]	312	& IRtoVIScoef,tau_pref_scenario)
[2921]	313	ELSE IF ((iaervar.ge.24).and.(iaervar.le.36))
[607]	314	& THEN ! << MY... dust scenarios >>
[2415]	315	call read_dust_scenario(ngrid,nlayer,zday,pplev,
[2643]	316	& IRtoVIScoef,tau_pref_scenario)
[607]	317	ELSE IF ((iaervar.eq.4).or.
	318	& ((iaervar.ge.124).and.(iaervar.le.126))) THEN
	319	! "old" TES assimation dust scenario (values at 700Pa in files!)
[2415]	320	call read_dust_scenario(ngrid,nlayer,zday,pplev,
[2643]	321	& IRtoVIScoef,tau_pref_scenario)
	322	!!!!!!!!!!!!!!!!!!!!!!!!!!!
[38]	323	ELSE
[2304]	324	call abort_physic("aeropacity","wrong value for iaervar",1)
[38]	325	ENDIF
	326
[2413]	327	! -----------------------------------------------------------------
	328	! 2. Compute/set the opacity of each aerosol in each layer
	329	! -----------------------------------------------------------------
[38]	330
[2413]	331	DO iaer = 1, naerkind ! Loop on all aerosols
[38]	332	c --------------------------------------------
	333	aerkind: SELECT CASE (name_iaer(iaer))
	334	c==================================================================
	335	CASE("dust_conrath") aerkind ! Typical dust profile
	336	c==================================================================
	337
	338	c Altitude of the top of the dust layer
	339	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	340	zlsconst=SIN(ls-2.76)
	341	if (iddist.eq.1) then
	342	do ig=1,ngrid
	343	topdust(ig)=topdustref ! constant dust layer top
	344	end do
	345
	346	else if (iddist.eq.2) then ! "Viking" scenario
	347	do ig=1,ngrid
[1224]	348	! altitude of the top of the aerosol layer (km) at Ls=2.76rad:
	349	! in the Viking year scenario
	350	topdust0(ig)=60. -22.sinlat(ig)*2
[38]	351	topdust(ig)=topdust0(ig)+18.*zlsconst
	352	end do
	353
	354	else if(iddist.eq.3) then !"MGS" scenario
	355	do ig=1,ngrid
	356	topdust(ig)=60.+18.*zlsconst
[1541]	357	& -(32+18zlsconst)sin(latitude(ig))**4
	358	& - 8zlsconst(sin(latitude(ig)))**5
[38]	359	end do
	360	endif
	361
	362	c Optical depth in each layer :
	363	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	364	if(iddist.ge.1) then
	365
	366	expfactor=0.
	367	DO l=1,nlayer
	368	DO ig=1,ngrid
	369	c Typical mixing ratio profile
[607]	370	if(pplay(ig,l).gt.odpref
[38]	371	$ /(988.**(topdust(ig)/70.))) then
[607]	372	zp=(odpref/pplay(ig,l))**(70./topdust(ig))
[38]	373	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
	374	else
	375	expfactor=1.e-3
	376	endif
	377	c Vertical scaling function
	378	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) *
	379	& expfactor *
	380	& QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
	381	ENDDO
	382	ENDDO
	383
	384	else if(iddist.eq.0) then
	385	c old dust vertical distribution function (pollack90)
	386	DO l=1,nlayer
	387	DO ig=1,ngrid
[607]	388	zp=odpref/pplay(ig,l)
[2415]	389	aerosol(ig,l,1)= tau_pref_scenario(ig)/odpref *
[38]	390	s (pplev(ig,l)-pplev(ig,l+1))
	391	s max( exp(.03(1.-max(zp,1.))) , 1.E-3 )
	392	ENDDO
	393	ENDDO
	394	end if
	395
	396	c==================================================================
[2246]	397	CASE("dust_doubleq") aerkind! Two-moment scheme for background dust
[38]	398	c (transport of mass and number mixing ratio)
	399	c==================================================================
[2643]	400	! Some initialisations for the IRtoVIScoef
	401	aerosol_IRabs(:,:)=0.
	402	taudust_IRabs(:)=0.
	403	taudust_VISext(:)=0.
	404
[38]	405	DO l=1,nlayer
	406	IF (l.LE.cstdustlevel) THEN
	407	c Opacity in the first levels is held constant to
	408	c avoid unrealistic values due to constant lifting:
	409	DO ig=1,ngrid
[2414]	410	! OPTICAL DEPTH used in the radiative transfer
[2246]	411	! => visible wavelength
[38]	412	aerosol(ig,l,iaer) =
	413	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
	414	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
	415	& pq(ig,cstdustlevel,igcm_dust_mass) *
	416	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
[2246]	417	! DENSITY SCALED OPACITY :
[2414]	418	! Diagnostic output to be compared with observations
[2246]	419	! => infrared wavelength
[2161]	420	dsodust(ig,l) =
	421	& ( 0.75 * QREFir3d(ig,cstdustlevel,iaer) /
	422	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
	423	& pq(ig,cstdustlevel,igcm_dust_mass)
[2643]	424
	425	if (reff_driven_IRtoVIS_scenario) then
	426	if ((clearatm).and.(nohmons)) then ! the IRtoVIScoef is computed only during the first call to the RT
	427	! OPTICAL DEPTH in IR absorption to compute the IRtoVIScoef
	428	aerosol_IRabs(ig,l) =
	429	& ( 0.75 * QREFir3d(ig,cstdustlevel,iaer) /
	430	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
	431	& ( 1. - omegaREFir3d(ig,cstdustlevel,iaer) ) *
	432	& pq(ig,cstdustlevel,igcm_dust_mass) *
	433	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	434	endif
	435	endif
[38]	436	ENDDO
	437	ELSE
	438	DO ig=1,ngrid
[2414]	439	! OPTICAL DEPTH used in the radiative transfer
[2246]	440	! => visible wavelength
	441	aerosol(ig,l,iaer) =
[38]	442	& ( 0.75 * QREFvis3d(ig,l,iaer) /
	443	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
	444	& pq(ig,l,igcm_dust_mass) *
	445	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
[2246]	446	! DENSITY SCALED OPACITY :
[2414]	447	! Diagnostic output to be compared with observations
[2246]	448	! => infrared wavelength
[2161]	449	dsodust(ig,l) =
	450	& ( 0.75 * QREFir3d(ig,l,iaer) /
	451	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
	452	& pq(ig,l,igcm_dust_mass)
[2643]	453
	454	if (reff_driven_IRtoVIS_scenario) then
	455	if ((clearatm).and.(nohmons)) then ! the IRtoVIScoef is computed only during the first call to the RT
	456	! OPTICAL DEPTH in IR absorption to compute the IRtoVIScoef
	457	aerosol_IRabs(ig,l) =
	458	& ( 0.75 * QREFir3d(ig,l,iaer) /
	459	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
	460	& ( 1. - omegaREFir3d(ig,l,iaer) ) *
	461	& pq(ig,l,igcm_dust_mass) *
	462	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	463	endif
	464	endif
[38]	465	ENDDO
	466	ENDIF
[2643]	467	if (reff_driven_IRtoVIS_scenario) then
	468	if ((clearatm).and.(nohmons)) then ! the IRtoVIScoef is computed only during the first call to the RT
	469	taudust_VISext(:) = taudust_VISext(:) + aerosol(:,l,iaer)
	470	taudust_IRabs(:) = taudust_IRabs(:) + aerosol_IRabs(:,l)
	471	endif
	472	endif
[38]	473	ENDDO
[2643]	474
	475	if (reff_driven_IRtoVIS_scenario) then
	476	if ((clearatm).and.(nohmons)) then ! the IRtoVIScoef is computed only during the first call to the RT
	477	IRtoVIScoef(:) = taudust_VISext(:) / taudust_IRabs(:)
	478	endif
	479	endif
[38]	480
	481	c==================================================================
	482	CASE("dust_submicron") aerkind ! Small dust population
	483	c==================================================================
	484
	485	DO l=1,nlayer
	486	IF (l.LE.cstdustlevel) THEN
	487	c Opacity in the first levels is held constant to
	488	c avoid unrealistic values due to constant lifting:
	489	DO ig=1,ngrid
	490	aerosol(ig,l,iaer) =
	491	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
	492	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
	493	& pq(ig,cstdustlevel,igcm_dust_submicron) *
	494	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	495	ENDDO
	496	ELSE
	497	DO ig=1,ngrid
	498	aerosol(ig,l,iaer) =
	499	& ( 0.75 * QREFvis3d(ig,l,iaer) /
	500	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
	501	& pq(ig,l,igcm_dust_submicron) *
	502	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	503	ENDDO
	504	ENDIF
	505	ENDDO
	506
	507	c==================================================================
	508	CASE("h2o_ice") aerkind ! Water ice crystals
	509	c==================================================================
	510
	511	c 1. Initialization
	512	aerosol(1:ngrid,1:nlayer,iaer) = 0.
	513	taucloudvis(1:ngrid) = 0.
	514	taucloudtes(1:ngrid) = 0.
	515	c 2. Opacity calculation
[1711]	516	! NO CLOUDS
	517	IF (clearsky) THEN
	518	aerosol(1:ngrid,1:nlayer,iaer) =1.e-9
	519	! CLOUDSs
	520	ELSE ! else (clearsky)
	521	DO ig=1, ngrid
	522	DO l=1,nlayer
	523	aerosol(ig,l,iaer) = max(1E-20,
	524	& ( 0.75 * QREFvis3d(ig,l,iaer) /
	525	& ( rho_ice * reffrad(ig,l,iaer) ) ) *
	526	& pq(ig,l,i_ice) *
	527	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
[38]	528	& )
[1711]	529	taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
	530	taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
	531	& QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
	532	& ( 1.E0 - omegaREFir3d(ig,l,iaer) )
	533	ENDDO
[38]	534	ENDDO
[1711]	535	! SUB-GRID SCALE CLOUDS
	536	IF (CLFvarying) THEN
	537	DO ig=1, ngrid
	538	DO l=1,nlayer-1
	539	CLFtot = max(totcloudfrac(ig),0.01)
	540	aerosol(ig,l,iaer)=
	541	& aerosol(ig,l,iaer)/CLFtot
	542	aerosol(ig,l,iaer) =
	543	& max(aerosol(ig,l,iaer),1.e-9)
	544	ENDDO
	545	ENDDO
	546	ENDIF ! end (CLFvarying)
	547	ENDIF ! end (clearsky)
	548
[38]	549	c==================================================================
[2447]	550	CASE("co2_ice") aerkind ! CO2 ice crystals
	551	c==================================================================
	552
	553	c 1. Initialization
	554	aerosol(1:ngrid,1:nlayer,iaer) = 0.
[2494]	555	taucloudco2vis(1:ngrid) = 0.
	556	taucloudco2tes(1:ngrid) = 0.
[2447]	557	c 2. Opacity calculation
	558	! NO CLOUDS
	559	IF (clearsky) THEN
[2494]	560	aerosol(1:ngrid,1:nlayer,iaer) = 1.e-9
[2447]	561	! CLOUDSs
	562	ELSE ! else (clearsky)
[2494]	563	DO ig = 1, ngrid
	564	DO l = 1, nlayer
	565	call density_co2_ice(dble(pt(ig,l)), rho_ice_co2)
	566
[2447]	567	aerosol(ig,l,iaer) = max(1E-20,
	568	& ( 0.75 * QREFvis3d(ig,l,iaer) /
[2494]	569	& ( rho_ice_co2 * reffrad(ig,l,iaer) ) ) *
[2447]	570	& pq(ig,l,i_co2ice) *
	571	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	572	& )
[2494]	573	taucloudco2vis(ig) = taucloudco2vis(ig)
	574	& + aerosol(ig,l,iaer)
	575	taucloudco2tes(ig) = taucloudco2tes(ig)
	576	& + aerosol(ig,l,iaer) *
[2447]	577	& QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
	578	& ( 1.E0 - omegaREFir3d(ig,l,iaer) )
	579	ENDDO
	580	ENDDO
	581	! SUB-GRID SCALE CLOUDS
[2494]	582	IF (CLFvaryingCO2) THEN
[2447]	583	DO ig=1, ngrid
[2494]	584	DO l= 1, nlayer-1
	585	CLFtotco2 = max(totcloudco2frac(ig),0.01)
[2447]	586	aerosol(ig,l,iaer)=
[2494]	587	& aerosol(ig,l,iaer)/CLFtotco2
[2447]	588	aerosol(ig,l,iaer) =
	589	& max(aerosol(ig,l,iaer),1.e-9)
	590	ENDDO
	591	ENDDO
[2494]	592	ENDIF ! end (CLFvaryingCO2)
[2447]	593	ENDIF ! end (clearsky)
	594
	595	c==================================================================
[1974]	596	CASE("stormdust_doubleq") aerkind ! CW17 : Two-moment scheme for
	597	c stormdust (transport of mass and number mixing ratio)
	598	c==================================================================
	599	c aerosol is calculated twice : once within the dust storm (clearatm=false)
	600	c and once in the part of the mesh without dust storm (clearatm=true)
	601	aerosol(1:ngrid,1:nlayer,iaer) = 0.
	602	IF (clearatm) THEN ! considering part of the mesh without storm
	603	aerosol(1:ngrid,1:nlayer,iaer)=1.e-25
	604	ELSE ! part of the mesh with concentred dust storm
	605	DO l=1,nlayer
	606	IF (l.LE.cstdustlevel) THEN
	607	c Opacity in the first levels is held constant to
	608	c avoid unrealistic values due to constant lifting:
	609	DO ig=1,ngrid
[2414]	610	! OPTICAL DEPTH used in the radiative transfer
[2246]	611	! => visible wavelength
	612	aerosol(ig,l,iaer) =
[1974]	613	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
	614	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
	615	& pq(ig,cstdustlevel,igcm_stormdust_mass) *
	616	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
[2246]	617	! DENSITY SCALED OPACITY :
[2414]	618	! Diagnostic output to be compared with observations
[2246]	619	! => infrared wavelength
	620	dsords(ig,l) =
	621	& ( 0.75 * QREFir3d(ig,cstdustlevel,iaer) /
	622	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
	623	& pq(ig,cstdustlevel,igcm_stormdust_mass)
[1974]	624	ENDDO
	625	ELSE
[2246]	626	DO ig=1,ngrid
[2414]	627	! OPTICAL DEPTH used in the radiative transfer
[2246]	628	! => visible wavelength
[1974]	629	aerosol(ig,l,iaer) =
	630	& ( 0.75 * QREFvis3d(ig,l,iaer) /
	631	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
	632	& pq(ig,l,igcm_stormdust_mass) *
	633	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
[2246]	634	! DENSITY SCALED OPACITY :
[2414]	635	! Diagnostic output to be compared with observations
[2246]	636	! => infrared wavelength
	637	dsords(ig,l) =
	638	& ( 0.75 * QREFir3d(ig,l,iaer) /
	639	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
	640	& pq(ig,l,igcm_stormdust_mass)
	641	ENDDO
[1974]	642	ENDIF
	643	ENDDO
	644	ENDIF
	645	c==================================================================
[2199]	646	CASE("topdust_doubleq") aerkind ! MV18 : Two-moment scheme for
	647	c topdust (transport of mass and number mixing ratio)
	648	c==================================================================
	649	c aerosol is calculated twice : once "above" the sub-grid mountain (nohmons=false)
	650	c and once in the part of the mesh without the sub-grid mountain (nohmons=true)
	651	aerosol(1:ngrid,1:nlayer,iaer) = 0.
[2634]	652	IF (nohmons) THEN ! considering part of the mesh without top dust flows
[2199]	653	aerosol(1:ngrid,1:nlayer,iaer)=1.e-25
[2628]	654	ELSE ! part of the mesh with concentrated dust flows
[2199]	655	DO l=1,nlayer
[2246]	656	IF (l.LE.cstdustlevel) THEN
	657	c Opacity in the first levels is held constant to
	658	c avoid unrealistic values due to constant lifting:
	659	DO ig=1,ngrid
[2414]	660	! OPTICAL DEPTH used in the radiative transfer
	661	! => visible wavelength
[2246]	662	aerosol(ig,l,iaer) =
	663	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
	664	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
	665	& pq(ig,cstdustlevel,igcm_topdust_mass) *
	666	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	667	! DENSITY SCALED OPACITY :
[2414]	668	! Diagnostic output to be compared with observations
[2246]	669	! => infrared wavelength
	670	dsotop(ig,l) =
	671	& ( 0.75 * QREFir3d(ig,cstdustlevel,iaer) /
	672	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
	673	& pq(ig,cstdustlevel,igcm_topdust_mass)
	674	ENDDO
	675	ELSE
	676	DO ig=1,ngrid
[2414]	677	! OPTICAL DEPTH used in the radiative transfer
[2246]	678	! => visible wavelength
	679	aerosol(ig,l,iaer) =
[2199]	680	& ( 0.75 * QREFvis3d(ig,l,iaer) /
	681	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
	682	& pq(ig,l,igcm_topdust_mass) *
	683	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
[2246]	684	! DENSITY SCALED OPACITY :
[2414]	685	! Diagnostic output to be compared with observations
[2246]	686	! => infrared wavelength
	687	dsotop(ig,l) =
	688	& ( 0.75 * QREFir3d(ig,l,iaer) /
	689	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
	690	& pq(ig,l,igcm_topdust_mass)
	691	ENDDO
	692	ENDIF
[2199]	693	ENDDO
	694	ENDIF
	695	c==================================================================
[38]	696	END SELECT aerkind
	697	c -----------------------------------
	698	ENDDO ! iaer (loop on aerosol kind)
	699
[2413]	700	! 3. Specific treatments for the dust aerosols
[38]	701
[2643]	702	! here IRtoVIScoef has been updated, we can call again read_dust_scenario
	703	if (reff_driven_IRtoVIS_scenario) then
	704	IF ((iaervar.ge.6).and.(iaervar.le.8)) THEN
	705	! clim, cold or warm synthetic scenarios
	706	call read_dust_scenario(ngrid,nlayer,zday,pplev,
	707	& IRtoVIScoef,tau_pref_scenario)
[2921]	708	ELSE IF ((iaervar.ge.24).and.(iaervar.le.36))
[2643]	709	& THEN ! << MY... dust scenarios >>
	710	call read_dust_scenario(ngrid,nlayer,zday,pplev,
	711	& IRtoVIScoef,tau_pref_scenario)
	712	ELSE IF ((iaervar.eq.4).or.
	713	& ((iaervar.ge.124).and.(iaervar.le.126))) THEN
	714	! "old" TES assimation dust scenario (values at 700Pa in files!)
	715	call read_dust_scenario(ngrid,nlayer,zday,pplev,
	716	& IRtoVIScoef,tau_pref_scenario)
	717	ENDIF
	718	endif
	719
[1375]	720	#ifdef DUSTSTORM
	721	c -----------------------------------------------------------------
[1410]	722	! Calculate reference opacity without perturbation
[1375]	723	c -----------------------------------------------------------------
	724	IF (firstcall) THEN
	725	DO iaer=1,naerdust
	726	DO l=1,nlayer
	727	DO ig=1,ngrid
[2415]	728	tau_pref_gcm(ig) = tau_pref_gcm(ig) +
[1375]	729	& aerosol(ig,l,iaerdust(iaer))
	730	ENDDO
	731	ENDDO
	732	ENDDO
[2415]	733	tau_pref_gcm(:) = tau_pref_gcm(:) * odpref / pplev(:,1)
[1410]	734
[1375]	735	c--------------------------------------------------
[1410]	736	c Get parameters of the opacity perturbation
[1375]	737	c--------------------------------------------------
[1410]	738	iaer=1 ! just change dust
[1375]	739
	740	write(,) "Add a local storm ?"
	741	localstorm=.true. ! default value
[2304]	742	call getin_p("localstorm",localstorm)
[1375]	743	write(,) " localstorm = ",localstorm
	744
	745	IF (localstorm) THEN
	746	WRITE(,) "********************"
	747	WRITE(,) "ADDING A LOCAL STORM"
	748	WRITE(,) "********************"
	749
	750	write(,) "ref opacity of local dust storm"
	751	taulocref = 4.25 ! default value
[2304]	752	call getin_p("taulocref",taulocref)
[1375]	753	write(,) " taulocref = ",taulocref
	754
	755	write(,) "target altitude of local storm (km)"
	756	ztoploc = 10.0 ! default value
[2304]	757	call getin_p("ztoploc",ztoploc)
[1375]	758	write(,) " ztoploc = ",ztoploc
	759
	760	write(,) "radius of dust storm (degree)"
	761	radloc = 0.5 ! default value
[2304]	762	call getin_p("radloc",radloc)
[1375]	763	write(,) " radloc = ",radloc
	764
	765	write(,) "center longitude of storm (deg)"
	766	lonloc = 25.0 ! default value
[2304]	767	call getin_p("lonloc",lonloc)
[1375]	768	write(,) " lonloc = ",lonloc
	769
	770	write(,) "center latitude of storm (deg)"
	771	latloc = -2.5 ! default value
[2304]	772	call getin_p("latloc",latloc)
[1375]	773	write(,) " latloc = ",latloc
	774
	775	write(,) "reff storm (mic) 0. for background"
	776	reffstorm = 0.0 ! default value
[2304]	777	call getin_p("reffstorm",reffstorm)
[1375]	778	write(,) " reffstorm = ",reffstorm
	779
[1410]	780	!! LOOP: modify opacity
[1375]	781	DO ig=1,ngrid
	782
[1410]	783	!! distance to the center:
[1541]	784	ray(ig)=SQRT((latitude(ig)180./pi-latloc)*2 +
	785	& (longitude(ig)180./pi -lonloc)*2)
[1375]	786
	787	!! transition factor for storm
[1410]	788	!! factor is hardcoded -- increase it to steepen
[1375]	789	yeah = (TANH(2.+(radloc-ray(ig))*10.)+1.)/2.
	790
[1410]	791	!! new opacity field
	792	!! -- add an opacity set to taulocref
	793	!! -- the additional reference opacity will
	794	!! thus be taulocref*odpref/pplev
[2415]	795	tauuser(ig)=max( tau_pref_gcm(ig) * pplev(ig,1) /odpref ,
[1410]	796	& taulocref * yeah )
[1375]	797
[1410]	798	!! compute l_top
[1375]	799	DO l=1,nlayer
	800	zalt(ig,l) = LOG( pplev(ig,1)/pplev(ig,l) )
	801	& / g / 44.01
	802	& * 8.31 * 210.
	803	IF ( (ztoploc .lt. zalt(ig,l) )
	804	& .and. (ztoploc .gt. zalt(ig,l-1)) ) l_top=l-1
	805	ENDDO
	806
[1410]	807	!! change reffrad if ever needed
[1375]	808	IF (reffstorm .gt. 0.) THEN
	809	DO l=1,nlayer
	810	IF (l .lt. l_top+1) THEN
	811	reffrad(ig,l,iaer) = max( reffrad(ig,l,iaer), reffstorm
	812	& * 1.e-6 * yeah )
	813	ENDIF
	814	ENDDO
	815	ENDIF
	816
[1410]	817	ENDDO
	818	!! END LOOP
[1375]	819
[1410]	820	!! compute perturbation in each layer (equation 8 in Spiga et al. JGR 2013)
[1375]	821	DO ig=1,ngrid
	822	int_factor(ig)=0.
	823	DO l=1,nlayer
	824	IF (l .lt. l_top+1) THEN
	825	int_factor(ig) =
	826	& int_factor(ig) +
	827	& ( 0.75 * QREFvis3d(ig,l,iaer) /
	828	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
	829	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	830	ENDIF
	831	ENDDO
	832	DO l=1, nlayer
[1410]	833	!! Mass mixing ratio perturbation due to local dust storm in each layer
[1375]	834	more_dust(ig,l,1)=
[2415]	835	& (tauuser(ig)-(tau_pref_gcm(ig)
[1375]	836	& * pplev(ig,1) /odpref)) /
	837	& int_factor(ig)
	838	more_dust(ig,l,2)=
[2415]	839	& (tauuser(ig)-(tau_pref_gcm(ig) *
[1375]	840	& pplev(ig,1) /odpref))
	841	& / int_factor(ig) *
	842	& ((ref_r0/reffrad(ig,l,iaer))**3)
	843	& * r3n_q
	844	ENDDO
	845	ENDDO
	846
[1410]	847	!! quantity of dust for each layer with the addition of the perturbation
	848	DO l=1, l_top
[1376]	849	pq(:,l,igcm_dust_mass)= pq(:,l,igcm_dust_mass)
[1410]	850	. + more_dust(:,l,1)
[1376]	851	pq(:,l,igcm_dust_number)= pq(:,l,igcm_dust_number)
[1410]	852	. + more_dust(:,l,2)
	853	ENDDO
	854	ENDIF !! IF (localstorm)
[2415]	855	tau_pref_gcm(:)=0.
[1410]	856	ENDIF !! IF (firstcall)
[1375]	857	#endif
	858
[2413]	859	!
[2417]	860	! 3.1. Compute "tauscaling" and "dust_rad_adjust", the dust rescaling
	861	! coefficients and adjust aerosol() dust opacities accordingly
	862	call compute_dustscaling(ngrid,nlayer,naerkind,naerdust,zday,pplev
[2643]	863	& ,tau_pref_scenario,IRtoVIScoef,
	864	& tauscaling,dust_rad_adjust,aerosol)
[358]	865
[2415]	866	! 3.2. Recompute tau_pref_gcm, the reference dust opacity, based on dust tracer
[2413]	867	! mixing ratios and their optical properties
[358]	868
[1088]	869	IF (freedust) THEN
[2415]	870	! Initialisation :
	871	tau_pref_gcm(:)=0
[1208]	872	DO iaer=1,naerdust
	873	DO l=1,nlayer
	874	DO ig=1,ngrid
[1410]	875	#ifdef DUSTSTORM
	876	!! recalculate opacity because storm perturbation has been added
	877	IF (firstcall) THEN
	878	aerosol(ig,l,iaer) =
	879	& ( 0.75 * QREFvis3d(ig,l,iaer) /
	880	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
	881	& pq(ig,l,igcm_dust_mass) *
	882	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	883	ENDIF
	884	#endif
[2415]	885	c MV19: tau_pref_gcm must ALWAYS contain the opacity of all dust tracers
	886	! GCM DUST OPTICAL DEPTH tau_pref_gcm is to be compared
	887	! with the observation CDOD tau_pref_scenario
[2414]	888	! => visible wavelength
[2199]	889	IF (name_iaer(iaerdust(iaer)).eq."dust_doubleq") THEN
[2415]	890	tau_pref_gcm(ig) = tau_pref_gcm(ig) +
[2199]	891	& ( 0.75 * QREFvis3d(ig,l,iaerdust(iaer)) /
	892	& ( rho_dust * reffrad(ig,l,iaerdust(iaer)) ) ) *
	893	& pq(ig,l,igcm_dust_mass) *
	894	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	895	ELSE IF (name_iaer(iaerdust(iaer)).eq."stormdust_doubleq") THEN
[2415]	896	tau_pref_gcm(ig) = tau_pref_gcm(ig) +
[2199]	897	& ( 0.75 * QREFvis3d(ig,l,iaerdust(iaer)) /
	898	& ( rho_dust * reffrad(ig,l,iaerdust(iaer)) ) ) *
	899	& pq(ig,l,igcm_stormdust_mass) *
	900	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	901	ELSE IF (name_iaer(iaerdust(iaer)).eq."topdust_doubleq") THEN
[2415]	902	tau_pref_gcm(ig) = tau_pref_gcm(ig) +
[2199]	903	& ( 0.75 * QREFvis3d(ig,l,iaerdust(iaer)) /
	904	& ( rho_dust * reffrad(ig,l,iaerdust(iaer)) ) ) *
	905	& pq(ig,l,igcm_topdust_mass) *
	906	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
	907	ENDIF
	908
[1208]	909	ENDDO
[1088]	910	ENDDO
	911	ENDDO
[2415]	912	tau_pref_gcm(:) = tau_pref_gcm(:) * odpref / pplev(:,1)
	913	ELSE
	914	! dust opacity strictly follows what is imposed by the dust scenario
	915	tau_pref_gcm(:)=tau_pref_scenario(:)
[2413]	916	ENDIF ! of IF (freedust)
[1974]	917
[2413]	918	! -----------------------------------------------------------------
	919	! 4. Total integrated visible optical depth of aerosols in each column
	920	! -----------------------------------------------------------------
[38]	921	DO iaer=1,naerkind
	922	do l=1,nlayer
	923	do ig=1,ngrid
	924	tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer)
	925	end do
	926	end do
	927	ENDDO
[1375]	928
[1974]	929
[1375]	930	#ifdef DUSTSTORM
	931	IF (firstcall) THEN
	932	firstcall=.false.
	933	ENDIF
	934	#endif
	935
[2413]	936	!
	937	! 5. Adapt aerosol() for the radiative transfer (i.e. account for
	938	! cases when it refers to a fraction of the global mesh)
	939	!
[38]	940
	941	c -----------------------------------------------------------------
[1974]	942	c aerosol/X for stormdust to prepare calculation of radiative transfer
	943	c -----------------------------------------------------------------
[2199]	944	IF (rdstorm) THEN
[1974]	945	DO l=1,nlayer
	946	DO ig=1,ngrid
[2199]	947	! stormdust: opacity relative to the storm fraction (stormdust/x)
[1974]	948	aerosol(ig,l,iaer_stormdust_doubleq) =
	949	& aerosol(ig,l,iaer_stormdust_doubleq)/totstormfract(ig)
	950	ENDDO
	951	ENDDO
[2199]	952	ENDIF
[1711]	953
[2199]	954	c -----------------------------------------------------------------
	955	c aerosol/X for topdust to prepare calculation of radiative transfer
	956	c -----------------------------------------------------------------
[2628]	957	IF (topflows) THEN
[2199]	958	DO ig=1,ngrid
[2628]	959	IF (contains_mons(ig)) THEN ! contains_mons=True ensures that alpha_hmons>0
[2199]	960	DO l=1,nlayer
[2628]	961	! topdust: opacity relative to the mons fraction (topdust/x)
[2199]	962	aerosol(ig,l,iaer_topdust_doubleq) =
	963	& aerosol(ig,l,iaer_topdust_doubleq)/alpha_hmons(ig)
	964	ENDDO
	965	ENDIF
	966	ENDDO
	967	ENDIF
[1974]	968
[1711]	969	END SUBROUTINE aeropacity
	970
	971	END MODULE aeropacity_mod

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

Download in other formats:

Original Format