Context Navigation

aeropacity_mod.F @ 2616

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

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

Download in other formats:

Original Format