Context Navigation

aeropacity_mod.F @ 2001

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

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

Download in other formats:

Original Format