Context Navigation

aeropacity_mod.F @ 2628

Last change on this file since 2628 was 2628, checked in by abierjon, 3 years ago

Mars GCM:
Big changes on mountain top dust flows for GCM6:

the scheme now activates only in grid meshes that contain a mountain among a hard-written list, instead of every meshes. This is done to prevent strong artificial reinjections of dust in places that don't present huge converging slopes enabling the concentration of dust (ex: Valles Marineris, Hellas). Topdust is now also detrained as soon as it leaves the column it originated from.
the list of the 19 allowed mountains is used by the subroutine topmons_setup in module topmons_mod, to compute a logical array contains_mons(ngrid). alpha_hmons and hsummit are also set up once and for all by this subroutine, which is called in physiq_mod's firstcall.
contains_mons, alpha_hmons and hsummit are now saved variables of the module surfdat_h, and are called as such and not as arguments in the subroutines using them.
the logical flag "slpwind" and the comments in the code have also been updated to the new terminology "mountain top dust flows", accordingly to ticket #71. The new flag read in callphys.def is "topflows".

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

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

Download in other formats:

Original Format