Context Navigation

aeropacity_mod.F @ 2472

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

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

Download in other formats:

Original Format