Context Navigation

aeropacity_mod.F @ 1921

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

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

Download in other formats:

Original Format