Context Navigation

aeropacity_mod.F @ 4063

Last change on this file since 4063 was 4030, checked in by emillour, 4 weeks ago

Mars PCM:

add possibility to read MY37 dust scenario and the possibility to read in a scenario provided in IR extinction opacity (at 9.3um) identified as such by being "cdodext" in the input file (whereas "cdod" is for IR absorption).
add possibility to load EUV daily E10.7 input for MY37.

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.37))
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