Context Navigation

aeropacity_mod.F @ 2437

Last change on this file since 2437 was 2417, checked in by emillour, 5 years ago

Mars GCM:
Add a new scheme to handle the dust and its radiative impact. Triggered using
a new flag dustscaling_mode=2 (dustscaling_mod=0: no rescaling at all, and
dustscaling_mode=1: full rescaling using tauscaling, GCMv5.3 style). Rescaling
is then only done on the radiative impact (see dust_scaling_mod.F90) of dust.
Moreover the scaling factor "dust_rad_adjust" is evaluated using the target dust
scenario opacity for the next sol and left to evolve linearly until then to not
impose the diurnal evolution of dust.
In practice, main changes or additions in the code are:

renamed flag "tauscaling_mode" as "dustscaling_mode"
moved parameter "t_scenario_sol" to "dust_param_mod"
adapted "compute_dustscaling" routine in "dust_scaling_mod"
added module "dust_rad_adjust_mod"
2D fields "dust_rad_adjust_prev" and "dust_rad_adjust_next" required to compute coefficient "dust_rad_adjust" need to be stored in (re)startfi files

File size: 32.9 KB

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

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

Download in other formats:

Original Format