Context Navigation

aeropacity_mod.F @ 2472

Last change on this file since 2472 was 2449, checked in by emillour, 5 years ago
Mars GCM: Minor cleanup and update in aeropacity_mod.F and read_dust_scenario.F90 to be able to use MY35 dust scenario. EM
File size: 34.7 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, igcm_co2_ice
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 co2 ice :
163	INTEGER,SAVE :: i_co2ice=0 ! co2 ice
164	! indexes of dust scatterers:
165	INTEGER,SAVE :: naerdust ! number of dust scatterers
166
167	! initializations
168	tau(1:ngrid,1:naerkind)=0
169
170	! identify tracers
171
172	!! AS: firstcall OK absolute
173	IF (firstcall) THEN
174	! identify scatterers that are dust
175	naerdust=0
176	DO iaer=1,naerkind
177	txt=name_iaer(iaer)
178	! CW17: choice tauscaling for stormdust or not
179	IF ((txt(1:4).eq."dust").OR.(txt(1:5).eq."storm")
180	& .OR.(txt(1:3).eq."top")) THEN !MV19: topdust tracer
181	naerdust=naerdust+1
182	iaerdust(naerdust)=iaer
183	ENDIF
184	ENDDO
185	! identify tracers which are dust
186	i=0
187	DO iq=1,nq
188	txt=noms(iq)
189	IF (txt(1:4).eq."dust") THEN
190	i=i+1
191	nqdust(i)=iq
192	ENDIF
193	ENDDO
194
195	IF (water.AND.activice) THEN
196	i_ice=igcm_h2o_ice
197	write(,) "aeropacity: i_ice=",i_ice
198	ENDIF
199
200	IF (co2clouds.AND.activeco2ice) THEN
201	i_co2ice=igcm_co2_ice
202	write(,) "aeropacity: i_co2ice =",i_co2ice
203	ENDIF
204
205	c typical profile of solsir and (1-w)^(-1):
206	c --- purely for diagnostics and printing
207	msolsir(1:nlayer,1:naerkind)=0
208	mqextsqabs(1:nlayer,1:naerkind)=0
209	WRITE(,) "Typical profiles of Qext(vis)/Qext(IR)"
210	WRITE(,) " and Qext(IR)/Qabs(IR):"
211	DO iaer = 1, naerkind ! Loop on aerosol kind
212	WRITE(,) "Aerosol # ",iaer
213	DO l=1,nlayer
214	DO ig=1,ngrid
215	msolsir(l,iaer)=msolsir(l,iaer)+
216	& QREFvis3d(ig,l,iaer)/
217	& QREFir3d(ig,l,iaer)
218	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)+
219	& (1.E0-omegaREFir3d(ig,l,iaer))**(-1)
220	ENDDO
221	msolsir(l,iaer)=msolsir(l,iaer)/REAL(ngrid)
222	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)/REAL(ngrid)
223	ENDDO
224	WRITE(,) "solsir: ",msolsir(:,iaer)
225	WRITE(,) "Qext/Qabs(IR): ",mqextsqabs(:,iaer)
226	ENDDO
227
228	! load value of tauvis from callphys.def (if given there,
229	! otherwise default value read from starfi.nc file will be used)
230	call getin_p("tauvis",tauvis)
231
232	IF (freedust.or.rdstorm) THEN ! if rdstorm no need to held opacity constant at the first levels
233	cstdustlevel = 1
234	ELSE
235	cstdustlevel = cstdustlevel0 !Opacity in the first levels is held constant to
236	!avoid unrealistic values due to constant lifting
237	ENDIF
238
239
240	#ifndef DUSTSTORM
241	firstcall=.false.
242	#endif
243
244	END IF ! end of if firstcall
245
246	! 1. Get prescribed tau_pref_scenario, Dust column optical depth at "odpref" Pa
247	!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
248
249	IF(iaervar.eq.1) THEN
250	do ig=1, ngrid
251	tau_pref_scenario(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
252	! or read in starfi
253	end do
254	ELSE IF (iaervar.eq.2) THEN ! << "Viking" Scenario>>
255
256	tau_pref_scenario(1) = 0.7+.3cos(ls+80.pi/180.) ! like seen by VL1
257	do ig=2,ngrid
258	tau_pref_scenario(ig) = tau_pref_scenario(1)
259	end do
260
261	ELSE IF (iaervar.eq.3) THEN ! << "MGS" scenario >>
262
263	taueq= 0.2 +(0.5-0.2) (cos(0.5(ls-4.363)))**14
264	tauS= 0.1 +(0.5-0.1) (cos(0.5(ls-4.363)))**14
265	tauN = 0.1
266	do ig=1,ngrid
267	if (latitude(ig).ge.0) then
268	! Northern hemisphere
269	tau_pref_scenario(ig)= tauN +
270	& (taueq-tauN)0.5(1+tanh((45-latitude(ig)180./pi)6/60))
271	else
272	! Southern hemisphere
273	tau_pref_scenario(ig)= tauS +
274	& (taueq-tauS)0.5(1+tanh((45+latitude(ig)180./pi)6/60))
275	endif
276	enddo ! of do ig=1,ngrid
277	ELSE IF (iaervar.eq.5) THEN ! << Escalier Scenario>>
278	tau_pref_scenario(1) = 2.5
279	if ((ls.ge.30.pi/180.).and.(ls.le.150.pi/180.))
280	& tau_pref_scenario(1) = .2
281
282	do ig=2,ngrid
283	tau_pref_scenario(ig) = tau_pref_scenario(1)
284	end do
285	ELSE IF ((iaervar.ge.6).and.(iaervar.le.8)) THEN
286	! clim, cold or warm synthetic scenarios
287	call read_dust_scenario(ngrid,nlayer,zday,pplev,
288	& tau_pref_scenario)
289	ELSE IF ((iaervar.ge.24).and.(iaervar.le.35))
290	& THEN ! << MY... dust scenarios >>
291	call read_dust_scenario(ngrid,nlayer,zday,pplev,
292	& tau_pref_scenario)
293	ELSE IF ((iaervar.eq.4).or.
294	& ((iaervar.ge.124).and.(iaervar.le.126))) THEN
295	! "old" TES assimation dust scenario (values at 700Pa in files!)
296	call read_dust_scenario(ngrid,nlayer,zday,pplev,
297	& tau_pref_scenario)
298	ELSE
299	call abort_physic("aeropacity","wrong value for iaervar",1)
300	ENDIF
301
302	! -----------------------------------------------------------------
303	! 2. Compute/set the opacity of each aerosol in each layer
304	! -----------------------------------------------------------------
305
306	DO iaer = 1, naerkind ! Loop on all aerosols
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
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
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
332	& -(32+18zlsconst)sin(latitude(ig))**4
333	& - 8zlsconst(sin(latitude(ig)))**5
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
345	if(pplay(ig,l).gt.odpref
346	$ /(988.**(topdust(ig)/70.))) then
347	zp=(odpref/pplay(ig,l))**(70./topdust(ig))
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
363	zp=odpref/pplay(ig,l)
364	aerosol(ig,l,1)= tau_pref_scenario(ig)/odpref *
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 background 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	! OPTICAL DEPTH used in the radiative transfer
382	! => visible wavelength
383	aerosol(ig,l,iaer) =
384	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
385	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
386	& pq(ig,cstdustlevel,igcm_dust_mass) *
387	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
388	! DENSITY SCALED OPACITY :
389	! Diagnostic output to be compared with observations
390	! => infrared wavelength
391	dsodust(ig,l) =
392	& ( 0.75 * QREFir3d(ig,cstdustlevel,iaer) /
393	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
394	& pq(ig,cstdustlevel,igcm_dust_mass)
395	ENDDO
396	ELSE
397	DO ig=1,ngrid
398	! OPTICAL DEPTH used in the radiative transfer
399	! => visible wavelength
400	aerosol(ig,l,iaer) =
401	& ( 0.75 * QREFvis3d(ig,l,iaer) /
402	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
403	& pq(ig,l,igcm_dust_mass) *
404	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
405	! DENSITY SCALED OPACITY :
406	! Diagnostic output to be compared with observations
407	! => infrared wavelength
408	dsodust(ig,l) =
409	& ( 0.75 * QREFir3d(ig,l,iaer) /
410	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
411	& pq(ig,l,igcm_dust_mass)
412	ENDDO
413	ENDIF
414	ENDDO
415
416	c==================================================================
417	CASE("dust_submicron") aerkind ! Small dust population
418	c==================================================================
419
420	DO l=1,nlayer
421	IF (l.LE.cstdustlevel) THEN
422	c Opacity in the first levels is held constant to
423	c avoid unrealistic values due to constant lifting:
424	DO ig=1,ngrid
425	aerosol(ig,l,iaer) =
426	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
427	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
428	& pq(ig,cstdustlevel,igcm_dust_submicron) *
429	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
430	ENDDO
431	ELSE
432	DO ig=1,ngrid
433	aerosol(ig,l,iaer) =
434	& ( 0.75 * QREFvis3d(ig,l,iaer) /
435	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
436	& pq(ig,l,igcm_dust_submicron) *
437	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
438	ENDDO
439	ENDIF
440	ENDDO
441
442	c==================================================================
443	CASE("h2o_ice") aerkind ! Water ice crystals
444	c==================================================================
445
446	c 1. Initialization
447	aerosol(1:ngrid,1:nlayer,iaer) = 0.
448	taucloudvis(1:ngrid) = 0.
449	taucloudtes(1:ngrid) = 0.
450	c 2. Opacity calculation
451	! NO CLOUDS
452	IF (clearsky) THEN
453	aerosol(1:ngrid,1:nlayer,iaer) =1.e-9
454	! CLOUDSs
455	ELSE ! else (clearsky)
456	DO ig=1, ngrid
457	DO l=1,nlayer
458	aerosol(ig,l,iaer) = max(1E-20,
459	& ( 0.75 * QREFvis3d(ig,l,iaer) /
460	& ( rho_ice * reffrad(ig,l,iaer) ) ) *
461	& pq(ig,l,i_ice) *
462	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
463	& )
464	taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
465	taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
466	& QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
467	& ( 1.E0 - omegaREFir3d(ig,l,iaer) )
468	ENDDO
469	ENDDO
470	! SUB-GRID SCALE CLOUDS
471	IF (CLFvarying) THEN
472	DO ig=1, ngrid
473	DO l=1,nlayer-1
474	CLFtot = max(totcloudfrac(ig),0.01)
475	aerosol(ig,l,iaer)=
476	& aerosol(ig,l,iaer)/CLFtot
477	aerosol(ig,l,iaer) =
478	& max(aerosol(ig,l,iaer),1.e-9)
479	ENDDO
480	ENDDO
481	ENDIF ! end (CLFvarying)
482	ENDIF ! end (clearsky)
483
484	c==================================================================
485	CASE("co2_ice") aerkind ! CO2 ice crystals
486	c==================================================================
487
488	c 1. Initialization
489	aerosol(1:ngrid,1:nlayer,iaer) = 0.
490	taucloudvis(1:ngrid) = 0.
491	taucloudtes(1:ngrid) = 0.
492	c 2. Opacity calculation
493	! NO CLOUDS
494	IF (clearsky) THEN
495	aerosol(1:ngrid,1:nlayer,iaer) =1.e-9
496	! CLOUDSs
497	ELSE ! else (clearsky)
498	DO ig=1, ngrid
499	DO l=1,nlayer
500	aerosol(ig,l,iaer) = max(1E-20,
501	& ( 0.75 * QREFvis3d(ig,l,iaer) /
502	& ( rho_ice * reffrad(ig,l,iaer) ) ) *
503	& pq(ig,l,i_co2ice) *
504	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
505	& )
506	taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
507	taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
508	& QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
509	& ( 1.E0 - omegaREFir3d(ig,l,iaer) )
510	ENDDO
511	ENDDO
512	! SUB-GRID SCALE CLOUDS
513	IF (CLFvarying) THEN
514	DO ig=1, ngrid
515	DO l=1,nlayer-1
516	CLFtot = max(totcloudfrac(ig),0.01)
517	aerosol(ig,l,iaer)=
518	& aerosol(ig,l,iaer)/CLFtot
519	aerosol(ig,l,iaer) =
520	& max(aerosol(ig,l,iaer),1.e-9)
521	ENDDO
522	ENDDO
523	ENDIF ! end (CLFvarying)
524	ENDIF ! end (clearsky)
525
526	c==================================================================
527	CASE("stormdust_doubleq") aerkind ! CW17 : Two-moment scheme for
528	c stormdust (transport of mass and number mixing ratio)
529	c==================================================================
530	c aerosol is calculated twice : once within the dust storm (clearatm=false)
531	c and once in the part of the mesh without dust storm (clearatm=true)
532	aerosol(1:ngrid,1:nlayer,iaer) = 0.
533	IF (clearatm) THEN ! considering part of the mesh without storm
534	aerosol(1:ngrid,1:nlayer,iaer)=1.e-25
535	ELSE ! part of the mesh with concentred dust storm
536	DO l=1,nlayer
537	IF (l.LE.cstdustlevel) THEN
538	c Opacity in the first levels is held constant to
539	c avoid unrealistic values due to constant lifting:
540	DO ig=1,ngrid
541	! OPTICAL DEPTH used in the radiative transfer
542	! => visible wavelength
543	aerosol(ig,l,iaer) =
544	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
545	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
546	& pq(ig,cstdustlevel,igcm_stormdust_mass) *
547	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
548	! DENSITY SCALED OPACITY :
549	! Diagnostic output to be compared with observations
550	! => infrared wavelength
551	dsords(ig,l) =
552	& ( 0.75 * QREFir3d(ig,cstdustlevel,iaer) /
553	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
554	& pq(ig,cstdustlevel,igcm_stormdust_mass)
555	ENDDO
556	ELSE
557	DO ig=1,ngrid
558	! OPTICAL DEPTH used in the radiative transfer
559	! => visible wavelength
560	aerosol(ig,l,iaer) =
561	& ( 0.75 * QREFvis3d(ig,l,iaer) /
562	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
563	& pq(ig,l,igcm_stormdust_mass) *
564	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
565	! DENSITY SCALED OPACITY :
566	! Diagnostic output to be compared with observations
567	! => infrared wavelength
568	dsords(ig,l) =
569	& ( 0.75 * QREFir3d(ig,l,iaer) /
570	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
571	& pq(ig,l,igcm_stormdust_mass)
572	ENDDO
573	ENDIF
574	ENDDO
575	ENDIF
576	c==================================================================
577	CASE("topdust_doubleq") aerkind ! MV18 : Two-moment scheme for
578	c topdust (transport of mass and number mixing ratio)
579	c==================================================================
580	c aerosol is calculated twice : once "above" the sub-grid mountain (nohmons=false)
581	c and once in the part of the mesh without the sub-grid mountain (nohmons=true)
582	aerosol(1:ngrid,1:nlayer,iaer) = 0.
583	IF (nohmons) THEN ! considering part of the mesh without storm
584	aerosol(1:ngrid,1:nlayer,iaer)=1.e-25
585	ELSE ! part of the mesh with concentred dust storm
586	DO l=1,nlayer
587	IF (l.LE.cstdustlevel) THEN
588	c Opacity in the first levels is held constant to
589	c avoid unrealistic values due to constant lifting:
590	DO ig=1,ngrid
591	! OPTICAL DEPTH used in the radiative transfer
592	! => visible wavelength
593	aerosol(ig,l,iaer) =
594	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
595	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
596	& pq(ig,cstdustlevel,igcm_topdust_mass) *
597	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
598	! DENSITY SCALED OPACITY :
599	! Diagnostic output to be compared with observations
600	! => infrared wavelength
601	dsotop(ig,l) =
602	& ( 0.75 * QREFir3d(ig,cstdustlevel,iaer) /
603	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
604	& pq(ig,cstdustlevel,igcm_topdust_mass)
605	ENDDO
606	ELSE
607	DO ig=1,ngrid
608	! OPTICAL DEPTH used in the radiative transfer
609	! => visible wavelength
610	aerosol(ig,l,iaer) =
611	& ( 0.75 * QREFvis3d(ig,l,iaer) /
612	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
613	& pq(ig,l,igcm_topdust_mass) *
614	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
615	! DENSITY SCALED OPACITY :
616	! Diagnostic output to be compared with observations
617	! => infrared wavelength
618	dsotop(ig,l) =
619	& ( 0.75 * QREFir3d(ig,l,iaer) /
620	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
621	& pq(ig,l,igcm_topdust_mass)
622	ENDDO
623	ENDIF
624	ENDDO
625	ENDIF
626	c==================================================================
627	END SELECT aerkind
628	c -----------------------------------
629	ENDDO ! iaer (loop on aerosol kind)
630
631	! 3. Specific treatments for the dust aerosols
632
633	#ifdef DUSTSTORM
634	c -----------------------------------------------------------------
635	! Calculate reference opacity without perturbation
636	c -----------------------------------------------------------------
637	IF (firstcall) THEN
638	DO iaer=1,naerdust
639	DO l=1,nlayer
640	DO ig=1,ngrid
641	tau_pref_gcm(ig) = tau_pref_gcm(ig) +
642	& aerosol(ig,l,iaerdust(iaer))
643	ENDDO
644	ENDDO
645	ENDDO
646	tau_pref_gcm(:) = tau_pref_gcm(:) * odpref / pplev(:,1)
647
648	c--------------------------------------------------
649	c Get parameters of the opacity perturbation
650	c--------------------------------------------------
651	iaer=1 ! just change dust
652
653	write(,) "Add a local storm ?"
654	localstorm=.true. ! default value
655	call getin_p("localstorm",localstorm)
656	write(,) " localstorm = ",localstorm
657
658	IF (localstorm) THEN
659	WRITE(,) "********************"
660	WRITE(,) "ADDING A LOCAL STORM"
661	WRITE(,) "********************"
662
663	write(,) "ref opacity of local dust storm"
664	taulocref = 4.25 ! default value
665	call getin_p("taulocref",taulocref)
666	write(,) " taulocref = ",taulocref
667
668	write(,) "target altitude of local storm (km)"
669	ztoploc = 10.0 ! default value
670	call getin_p("ztoploc",ztoploc)
671	write(,) " ztoploc = ",ztoploc
672
673	write(,) "radius of dust storm (degree)"
674	radloc = 0.5 ! default value
675	call getin_p("radloc",radloc)
676	write(,) " radloc = ",radloc
677
678	write(,) "center longitude of storm (deg)"
679	lonloc = 25.0 ! default value
680	call getin_p("lonloc",lonloc)
681	write(,) " lonloc = ",lonloc
682
683	write(,) "center latitude of storm (deg)"
684	latloc = -2.5 ! default value
685	call getin_p("latloc",latloc)
686	write(,) " latloc = ",latloc
687
688	write(,) "reff storm (mic) 0. for background"
689	reffstorm = 0.0 ! default value
690	call getin_p("reffstorm",reffstorm)
691	write(,) " reffstorm = ",reffstorm
692
693	!! LOOP: modify opacity
694	DO ig=1,ngrid
695
696	!! distance to the center:
697	ray(ig)=SQRT((latitude(ig)180./pi-latloc)*2 +
698	& (longitude(ig)180./pi -lonloc)*2)
699
700	!! transition factor for storm
701	!! factor is hardcoded -- increase it to steepen
702	yeah = (TANH(2.+(radloc-ray(ig))*10.)+1.)/2.
703
704	!! new opacity field
705	!! -- add an opacity set to taulocref
706	!! -- the additional reference opacity will
707	!! thus be taulocref*odpref/pplev
708	tauuser(ig)=max( tau_pref_gcm(ig) * pplev(ig,1) /odpref ,
709	& taulocref * yeah )
710
711	!! compute l_top
712	DO l=1,nlayer
713	zalt(ig,l) = LOG( pplev(ig,1)/pplev(ig,l) )
714	& / g / 44.01
715	& * 8.31 * 210.
716	IF ( (ztoploc .lt. zalt(ig,l) )
717	& .and. (ztoploc .gt. zalt(ig,l-1)) ) l_top=l-1
718	ENDDO
719
720	!! change reffrad if ever needed
721	IF (reffstorm .gt. 0.) THEN
722	DO l=1,nlayer
723	IF (l .lt. l_top+1) THEN
724	reffrad(ig,l,iaer) = max( reffrad(ig,l,iaer), reffstorm
725	& * 1.e-6 * yeah )
726	ENDIF
727	ENDDO
728	ENDIF
729
730	ENDDO
731	!! END LOOP
732
733	!! compute perturbation in each layer (equation 8 in Spiga et al. JGR 2013)
734	DO ig=1,ngrid
735	int_factor(ig)=0.
736	DO l=1,nlayer
737	IF (l .lt. l_top+1) THEN
738	int_factor(ig) =
739	& int_factor(ig) +
740	& ( 0.75 * QREFvis3d(ig,l,iaer) /
741	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
742	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
743	ENDIF
744	ENDDO
745	DO l=1, nlayer
746	!! Mass mixing ratio perturbation due to local dust storm in each layer
747	more_dust(ig,l,1)=
748	& (tauuser(ig)-(tau_pref_gcm(ig)
749	& * pplev(ig,1) /odpref)) /
750	& int_factor(ig)
751	more_dust(ig,l,2)=
752	& (tauuser(ig)-(tau_pref_gcm(ig) *
753	& pplev(ig,1) /odpref))
754	& / int_factor(ig) *
755	& ((ref_r0/reffrad(ig,l,iaer))**3)
756	& * r3n_q
757	ENDDO
758	ENDDO
759
760	!! quantity of dust for each layer with the addition of the perturbation
761	DO l=1, l_top
762	pq(:,l,igcm_dust_mass)= pq(:,l,igcm_dust_mass)
763	. + more_dust(:,l,1)
764	pq(:,l,igcm_dust_number)= pq(:,l,igcm_dust_number)
765	. + more_dust(:,l,2)
766	ENDDO
767	ENDIF !! IF (localstorm)
768	tau_pref_gcm(:)=0.
769	ENDIF !! IF (firstcall)
770	#endif
771
772	!
773	! 3.1. Compute "tauscaling" and "dust_rad_adjust", the dust rescaling
774	! coefficients and adjust aerosol() dust opacities accordingly
775	call compute_dustscaling(ngrid,nlayer,naerkind,naerdust,zday,pplev
776	& ,tau_pref_scenario,tauscaling,
777	& dust_rad_adjust,aerosol)
778
779	! 3.2. Recompute tau_pref_gcm, the reference dust opacity, based on dust tracer
780	! mixing ratios and their optical properties
781
782	IF (freedust) THEN
783	! Initialisation :
784	tau_pref_gcm(:)=0
785	DO iaer=1,naerdust
786	DO l=1,nlayer
787	DO ig=1,ngrid
788	#ifdef DUSTSTORM
789	!! recalculate opacity because storm perturbation has been added
790	IF (firstcall) THEN
791	aerosol(ig,l,iaer) =
792	& ( 0.75 * QREFvis3d(ig,l,iaer) /
793	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
794	& pq(ig,l,igcm_dust_mass) *
795	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
796	ENDIF
797	#endif
798	c MV19: tau_pref_gcm must ALWAYS contain the opacity of all dust tracers
799	! GCM DUST OPTICAL DEPTH tau_pref_gcm is to be compared
800	! with the observation CDOD tau_pref_scenario
801	! => visible wavelength
802	IF (name_iaer(iaerdust(iaer)).eq."dust_doubleq") THEN
803	tau_pref_gcm(ig) = tau_pref_gcm(ig) +
804	& ( 0.75 * QREFvis3d(ig,l,iaerdust(iaer)) /
805	& ( rho_dust * reffrad(ig,l,iaerdust(iaer)) ) ) *
806	& pq(ig,l,igcm_dust_mass) *
807	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
808	ELSE IF (name_iaer(iaerdust(iaer)).eq."stormdust_doubleq") THEN
809	tau_pref_gcm(ig) = tau_pref_gcm(ig) +
810	& ( 0.75 * QREFvis3d(ig,l,iaerdust(iaer)) /
811	& ( rho_dust * reffrad(ig,l,iaerdust(iaer)) ) ) *
812	& pq(ig,l,igcm_stormdust_mass) *
813	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
814	ELSE IF (name_iaer(iaerdust(iaer)).eq."topdust_doubleq") THEN
815	tau_pref_gcm(ig) = tau_pref_gcm(ig) +
816	& ( 0.75 * QREFvis3d(ig,l,iaerdust(iaer)) /
817	& ( rho_dust * reffrad(ig,l,iaerdust(iaer)) ) ) *
818	& pq(ig,l,igcm_topdust_mass) *
819	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
820	ENDIF
821
822	ENDDO
823	ENDDO
824	ENDDO
825	tau_pref_gcm(:) = tau_pref_gcm(:) * odpref / pplev(:,1)
826	ELSE
827	! dust opacity strictly follows what is imposed by the dust scenario
828	tau_pref_gcm(:)=tau_pref_scenario(:)
829	ENDIF ! of IF (freedust)
830
831	! -----------------------------------------------------------------
832	! 4. Total integrated visible optical depth of aerosols in each column
833	! -----------------------------------------------------------------
834	DO iaer=1,naerkind
835	do l=1,nlayer
836	do ig=1,ngrid
837	tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer)
838	end do
839	end do
840	ENDDO
841
842
843	#ifdef DUSTSTORM
844	IF (firstcall) THEN
845	firstcall=.false.
846	ENDIF
847	#endif
848
849	!
850	! 5. Adapt aerosol() for the radiative transfer (i.e. account for
851	! cases when it refers to a fraction of the global mesh)
852	!
853
854	c -----------------------------------------------------------------
855	c aerosol/X for stormdust to prepare calculation of radiative transfer
856	c -----------------------------------------------------------------
857	IF (rdstorm) THEN
858	DO l=1,nlayer
859	DO ig=1,ngrid
860	! stormdust: opacity relative to the storm fraction (stormdust/x)
861	aerosol(ig,l,iaer_stormdust_doubleq) =
862	& aerosol(ig,l,iaer_stormdust_doubleq)/totstormfract(ig)
863	ENDDO
864	ENDDO
865	ENDIF
866
867	c -----------------------------------------------------------------
868	c aerosol/X for topdust to prepare calculation of radiative transfer
869	c -----------------------------------------------------------------
870	IF (slpwind) THEN
871	DO ig=1,ngrid
872	IF (alpha_hmons(ig) .gt. 0.) THEN
873	DO l=1,nlayer
874	! topdust: opacity relative to the storm fraction (topdust/x)
875	aerosol(ig,l,iaer_topdust_doubleq) =
876	& aerosol(ig,l,iaer_topdust_doubleq)/alpha_hmons(ig)
877	ENDDO
878	ENDIF
879	ENDDO
880	ENDIF
881
882	END SUBROUTINE aeropacity
883
884	END MODULE aeropacity_mod

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

Download in other formats:

Original Format