Context Navigation

aeropacity_mod.F @ 2212

Last change on this file since 2212 was 2199, checked in by mvals, 6 years ago
Mars GCM: Implementation of a new parametrization of the dust entrainment by slope winds above the sub-grid scale topography. The parametrization is activated with the flag slpwind=.true. (set to "false" by default) in callphys.def. The new parametrization involves the new tracers topdust_mass and topdust_number. MV
File size: 33.8 KB

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

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

Download in other formats:

Original Format