Context Navigation

aeropacity_mod.F @ 2246

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

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

Download in other formats:

Original Format