Context Navigation

aeropacity_mod.F @ 2514

Last change on this file since 2514 was 2494, checked in by cmathe, 4 years ago

Mars GCM:
co2_ice as scatterer in radiative transfert. Need co2clouds and

activeco2ice .eqv. True. Files involved:

aeropacity_mod.F
callradite_mod.F
physiq_mod.F
updatereffrad_mod.F
suaer.F90

determine co2_ice density from temperature. Used in riceco2 computation.

Files involved:

co2cloud.F90
improvedco2clouds_mod.F90
updaterad.F90
updatereffrad_mod.F

co2condens_mod4micro.F: variable initialization
initracer.F: add nuiceco2_ref = 0.2
phyredem.F: remove co2_ice from qsurf since co2_ice => co2ice
watercloud_mod.F: tiny typo

File size: 35.2 KB

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

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

Download in other formats:

Original Format