Context Navigation

aeropacity_mod.F @ 3026

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

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

Download in other formats:

Original Format