Context Navigation

aeropacity_mod.F @ 2526

Last change on this file since 2526 was 2494, checked in by cmathe, 4 years ago

Mars GCM:
co2_ice as scatterer in radiative transfert. Need co2clouds and

activeco2ice .eqv. True. Files involved:

aeropacity_mod.F
callradite_mod.F
physiq_mod.F
updatereffrad_mod.F
suaer.F90

determine co2_ice density from temperature. Used in riceco2 computation.

Files involved:

co2cloud.F90
improvedco2clouds_mod.F90
updaterad.F90
updatereffrad_mod.F

co2condens_mod4micro.F: variable initialization
initracer.F: add nuiceco2_ref = 0.2
phyredem.F: remove co2_ice from qsurf since co2_ice => co2ice
watercloud_mod.F: tiny typo

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

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

Download in other formats:

Original Format