Context Navigation

aeropacity_mod.F @ 3757

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

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

Download in other formats:

Original Format