Context Navigation

aeropacity_mod.F @ 2240

Last change on this file since 2240 was 2199, checked in by mvals, 6 years ago
Mars GCM: Implementation of a new parametrization of the dust entrainment by slope winds above the sub-grid scale topography. The parametrization is activated with the flag slpwind=.true. (set to "false" by default) in callphys.def. The new parametrization involves the new tracers topdust_mass and topdust_number. MV
File size: 33.8 KB

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

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

Download in other formats:

Original Format