Context Navigation

aeropacity_mod.F @ 2642

Last change on this file since 2642 was 2634, checked in by abierjon, 4 years ago

Mars GCM:
Changes on dust_rad_adjust for GCM6:

put an upper limit on dust_rad_adjust to avoid skyrocketing values when there are strong diurnal variations of opacity
move the condition to compute dust_rad_adjust only once per time step before the whole call to compute_dust_rad_adjust, instead of repetitive local checks
some cleaning

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

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

Download in other formats:

Original Format