Context Navigation

aeropacity_mod.F @ 2246

Last change on this file since 2246 was 2246, checked in by abierjon, 5 years ago
Mars GCM: Resolved ticket #32 : 1) dsodust is now calculated only once in the InfraRed? by aeropacity_mod (used to be wrongly calculated twice, such as dsodust=IR_part+Visible_part) ; 2) dsords is now calculated in the IR by aeropacity_mod (used to be calculated in the Visible) ; 3) dsotop is added and calculated in the IR in aeropacity_mod + some cleaning and commenting of the code AB
File size: 34.9 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,dsotop,
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(OUT) :: dsotop(ngrid,nlayer) !dso of topdust
94	REAL, INTENT(INOUT) :: reffrad(ngrid,nlayer,naerkind)
95	REAL, INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind)
96	REAL, INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
97	REAL, INTENT(IN) :: omegaREFir3d(ngrid,nlayer,naerkind)
98	LOGICAL, INTENT(IN) :: clearatm
99	REAL, INTENT(IN) :: totstormfract(ngrid)
100	LOGICAL, INTENT(IN) :: nohmons
101	REAL, INTENT(IN) :: alpha_hmons(ngrid)
102	REAL, INTENT(OUT) :: tauscaling(ngrid) ! Scaling factor for qdust and Ndust
103	REAL,INTENT(IN) :: totcloudfrac(ngrid) ! total cloud fraction
104	LOGICAL,INTENT(IN) :: clearsky ! true for part without clouds,false for part with clouds (total or sub-grid clouds)
105	c
106	c Local variables :
107	c -----------------
108	REAL CLFtot ! total cloud fraction
109	real expfactor
110	INTEGER l,ig,iq,i,j
111	INTEGER iaer ! Aerosol index
112	real topdust(ngrid)
113	real zlsconst, zp
114	real taueq,tauS,tauN
115	c Mean Qext(vis)/Qext(ir) profile
116	real msolsir(nlayer,naerkind)
117	c Mean Qext(ir)/Qabs(ir) profile
118	real mqextsqabs(nlayer,naerkind)
119	c Variables used when multiple particle sizes are used
120	c for dust or water ice particles in the radiative transfer
121	c (see callradite.F for more information).
122	REAL taudusttmp(ngrid)! Temporary dust opacity used before scaling
123	REAL taubackdusttmp(ngrid)! Temporary background dust opacity used before scaling
124	REAL taualldust(ngrid)! dust opacity all dust
125	REAL taudust(ngrid)! dust opacity dust doubleq
126	REAL taustormdust(ngrid)! dust opacity stormdust doubleq
127	REAL taustormdusttmp(ngrid)! dust opacity stormdust doubleq before tauscaling
128	REAL taudustvis(ngrid) ! Dust opacity after scaling
129	REAL taudusttes(ngrid) ! Dust opacity at IR ref. wav. as
130	! "seen" by the GCM.
131	REAL taucloudvis(ngrid)! Cloud opacity at visible
132	! reference wavelength
133	REAL taucloudtes(ngrid)! Cloud opacity at infrared
134	! reference wavelength using
135	! Qabs instead of Qext
136	! (direct comparison with TES)
137	REAL topdust0(ngrid)
138
139	#ifdef DUSTSTORM
140	!! Local dust storms
141	logical localstorm ! =true to create a local dust storm
142	real taulocref,ztoploc,radloc,lonloc,latloc ! local dust storm parameters
143	real reffstorm, yeah
144	REAL ray(ngrid) ! distance from dust storm center
145	REAL tauuser(ngrid) ! opacity perturbation due to dust storm
146	REAL more_dust(ngrid,nlayer,2) ! Mass mixing ratio perturbation due to the dust storm
147	REAL int_factor(ngrid) ! useful factor to compute mmr perturbation
148	real l_top ! layer of the storm's top
149	REAL zalt(ngrid, nlayer) ! useful factor to compute l_top
150	#endif
151
152	c local saved variables
153	c ---------------------
154
155	c Level under which the dust mixing ratio is held constant
156	c when computing the dust opacity in each layer
157	c (this applies when doubleq and active are true)
158	INTEGER, PARAMETER :: cstdustlevel0 = 7
159	INTEGER, SAVE :: cstdustlevel
160
161	LOGICAL,SAVE :: firstcall=.true.
162
163	! indexes of water ice and dust tracers:
164	INTEGER,SAVE :: i_ice=0 ! water ice
165	real,parameter :: odpref=610. ! DOD reference pressure (Pa)
166	CHARACTER(LEN=20) :: txt ! to temporarly store text
167	CHARACTER(LEN=1) :: txt2 ! to temporarly store text
168	! indexes of dust scatterers:
169	INTEGER,SAVE :: naerdust ! number of dust scatterers
170
171	tau(1:ngrid,1:naerkind)=0
172	dsords(:,:)=0. !CW17: initialize dsords
173	dsodust(:,:)=0.
174	dsotop(:,:)=0.
175
176	! identify tracers
177
178	!! AS: firstcall OK absolute
179	IF (firstcall) THEN
180	! identify scatterers that are dust
181	naerdust=0
182	DO iaer=1,naerkind
183	txt=name_iaer(iaer)
184	! CW17: choice tauscaling for stormdust or not
185	IF ((txt(1:4).eq."dust").OR.(txt(1:5).eq."storm")
186	& .OR.(txt(1:3).eq."top")) THEN !MV19: topdust tracer
187	naerdust=naerdust+1
188	iaerdust(naerdust)=iaer
189	ENDIF
190	ENDDO
191	! identify tracers which are dust
192	i=0
193	DO iq=1,nq
194	txt=noms(iq)
195	IF (txt(1:4).eq."dust") THEN
196	i=i+1
197	nqdust(i)=iq
198	ENDIF
199	ENDDO
200
201	IF (water.AND.activice) THEN
202	i_ice=igcm_h2o_ice
203	write(,) "aeropacity: i_ice=",i_ice
204	ENDIF
205
206	c typical profile of solsir and (1-w)^(-1):
207	c --- purely for diagnostics and printing
208	msolsir(1:nlayer,1:naerkind)=0
209	mqextsqabs(1:nlayer,1:naerkind)=0
210	WRITE(,) "Typical profiles of Qext(vis)/Qext(IR)"
211	WRITE(,) " and Qext(IR)/Qabs(IR):"
212	DO iaer = 1, naerkind ! Loop on aerosol kind
213	WRITE(,) "Aerosol # ",iaer
214	DO l=1,nlayer
215	DO ig=1,ngrid
216	msolsir(l,iaer)=msolsir(l,iaer)+
217	& QREFvis3d(ig,l,iaer)/
218	& QREFir3d(ig,l,iaer)
219	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)+
220	& (1.E0-omegaREFir3d(ig,l,iaer))**(-1)
221	ENDDO
222	msolsir(l,iaer)=msolsir(l,iaer)/REAL(ngrid)
223	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)/REAL(ngrid)
224	ENDDO
225	WRITE(,) "solsir: ",msolsir(:,iaer)
226	WRITE(,) "Qext/Qabs(IR): ",mqextsqabs(:,iaer)
227	ENDDO
228
229	! load value of tauvis from callphys.def (if given there,
230	! otherwise default value read from starfi.nc file will be used)
231	call getin("tauvis",tauvis)
232
233	IF (freedust.or.rdstorm) THEN ! if rdstorm no need to held opacity constant at the first levels
234	cstdustlevel = 1
235	ELSE
236	cstdustlevel = cstdustlevel0 !Opacity in the first levels is held constant to
237	!avoid unrealistic values due to constant lifting
238	ENDIF
239
240
241	#ifndef DUSTSTORM
242	firstcall=.false.
243	#endif
244
245	END IF
246
247	c Vertical column optical depth at "odpref" Pa
248	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
249	IF(freedust) THEN
250	tauref(:) = 0. ! tauref is computed after, instead of being forced
251
252	ELSE IF(iaervar.eq.1) THEN
253	do ig=1, ngrid
254	tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
255	! or read in starfi
256	end do
257	ELSE IF (iaervar.eq.2) THEN ! << "Viking" Scenario>>
258
259	tauref(1) = 0.7+.3cos(ls+80.pi/180.) ! like seen by VL1
260	do ig=2,ngrid
261	tauref(ig) = tauref(1)
262	end do
263
264	ELSE IF (iaervar.eq.3) THEN ! << "MGS" scenario >>
265
266	taueq= 0.2 +(0.5-0.2) (cos(0.5(ls-4.363)))**14
267	tauS= 0.1 +(0.5-0.1) (cos(0.5(ls-4.363)))**14
268	tauN = 0.1
269	c if (peri_day.eq.150) then
270	c tauS=0.1
271	c tauN=0.1 +(0.5-0.1) (cos(0.5(ls+pi-4.363)))**14
272	c taueq= 0.2 +(0.5-0.2) (cos(0.5(ls+pi-4.363)))**14
273	c endif
274	do ig=1,ngrid
275	if (latitude(ig).ge.0) then
276	! Northern hemisphere
277	tauref(ig)= tauN +
278	& (taueq-tauN)0.5(1+tanh((45-latitude(ig)180./pi)6/60))
279	else
280	! Southern hemisphere
281	tauref(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	c tauref(1) = 0.2
287	c if ((ls.ge.210.pi/180.).and.(ls.le.330.pi/180.))
288	c & tauref(1) = 2.5
289	tauref(1) = 2.5
290	if ((ls.ge.30.pi/180.).and.(ls.le.150.pi/180.))
291	& tauref(1) = .2
292
293	do ig=2,ngrid
294	tauref(ig) = tauref(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,tauref)
299	ELSE IF ((iaervar.ge.24).and.(iaervar.le.34))
300	& THEN ! << MY... dust scenarios >>
301	call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
302	ELSE IF ((iaervar.eq.4).or.
303	& ((iaervar.ge.124).and.(iaervar.le.126))) THEN
304	! "old" TES assimation dust scenario (values at 700Pa in files!)
305	call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
306	ELSE
307	stop 'problem with iaervar in aeropacity.F'
308	ENDIF
309
310	c -----------------------------------------------------------------
311	c Computing the opacity in each layer
312	c -----------------------------------------------------------------
313
314	DO iaer = 1, naerkind ! Loop on aerosol kind
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)= tauref(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 for the computation of tauref,
390	! which is to be compared with tauref_scenario
391	! => visible wavelength
392	aerosol(ig,l,iaer) =
393	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
394	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
395	& pq(ig,cstdustlevel,igcm_dust_mass) *
396	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
397	! DENSITY SCALED OPACITY :
398	! GCM output to be compared with observations
399	! => infrared wavelength
400	dsodust(ig,l) =
401	& ( 0.75 * QREFir3d(ig,cstdustlevel,iaer) /
402	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
403	& pq(ig,cstdustlevel,igcm_dust_mass)
404	ENDDO
405	ELSE
406	DO ig=1,ngrid
407	! OPTICAL DEPTH for the computation of tauref,
408	! which is to be compared with tauref_scenario
409	! => visible wavelength
410	aerosol(ig,l,iaer) =
411	& ( 0.75 * QREFvis3d(ig,l,iaer) /
412	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
413	& pq(ig,l,igcm_dust_mass) *
414	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
415	! DENSITY SCALED OPACITY :
416	! GCM output to be compared with observations
417	! => infrared wavelength
418	dsodust(ig,l) =
419	& ( 0.75 * QREFir3d(ig,l,iaer) /
420	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
421	& pq(ig,l,igcm_dust_mass)
422	ENDDO
423	ENDIF
424	ENDDO
425
426	c==================================================================
427	CASE("dust_submicron") aerkind ! Small dust population
428	c==================================================================
429
430	DO l=1,nlayer
431	IF (l.LE.cstdustlevel) THEN
432	c Opacity in the first levels is held constant to
433	c avoid unrealistic values due to constant lifting:
434	DO ig=1,ngrid
435	aerosol(ig,l,iaer) =
436	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
437	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
438	& pq(ig,cstdustlevel,igcm_dust_submicron) *
439	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
440	ENDDO
441	ELSE
442	DO ig=1,ngrid
443	aerosol(ig,l,iaer) =
444	& ( 0.75 * QREFvis3d(ig,l,iaer) /
445	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
446	& pq(ig,l,igcm_dust_submicron) *
447	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
448	ENDDO
449	ENDIF
450	ENDDO
451
452	c==================================================================
453	CASE("h2o_ice") aerkind ! Water ice crystals
454	c==================================================================
455
456	c 1. Initialization
457	aerosol(1:ngrid,1:nlayer,iaer) = 0.
458	taucloudvis(1:ngrid) = 0.
459	taucloudtes(1:ngrid) = 0.
460	c 2. Opacity calculation
461	! NO CLOUDS
462	IF (clearsky) THEN
463	aerosol(1:ngrid,1:nlayer,iaer) =1.e-9
464	! CLOUDSs
465	ELSE ! else (clearsky)
466	DO ig=1, ngrid
467	DO l=1,nlayer
468	aerosol(ig,l,iaer) = max(1E-20,
469	& ( 0.75 * QREFvis3d(ig,l,iaer) /
470	& ( rho_ice * reffrad(ig,l,iaer) ) ) *
471	& pq(ig,l,i_ice) *
472	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
473	& )
474	taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
475	taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
476	& QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
477	& ( 1.E0 - omegaREFir3d(ig,l,iaer) )
478	ENDDO
479	ENDDO
480	! SUB-GRID SCALE CLOUDS
481	IF (CLFvarying) THEN
482	DO ig=1, ngrid
483	DO l=1,nlayer-1
484	CLFtot = max(totcloudfrac(ig),0.01)
485	aerosol(ig,l,iaer)=
486	& aerosol(ig,l,iaer)/CLFtot
487	aerosol(ig,l,iaer) =
488	& max(aerosol(ig,l,iaer),1.e-9)
489	ENDDO
490	ENDDO
491	! ELSE ! else (CLFvarying)
492	! DO ig=1, ngrid
493	! DO l=1,nlayer-1 ! to stop the rad tran bug
494	! CLFtot = CLFfixval
495	! aerosol(ig,l,iaer)=
496	! & aerosol(ig,l,iaer)/CLFtot
497	! aerosol(ig,l,iaer) =
498	! & max(aerosol(ig,l,iaer),1.e-9)
499	! ENDDO
500	! ENDDO
501	ENDIF ! end (CLFvarying)
502	ENDIF ! end (clearsky)
503
504	c==================================================================
505	CASE("stormdust_doubleq") aerkind ! CW17 : Two-moment scheme for
506	c stormdust (transport of mass and number mixing ratio)
507	c==================================================================
508	c aerosol is calculated twice : once within the dust storm (clearatm=false)
509	c and once in the part of the mesh without dust storm (clearatm=true)
510	aerosol(1:ngrid,1:nlayer,iaer) = 0.
511	IF (clearatm) THEN ! considering part of the mesh without storm
512	aerosol(1:ngrid,1:nlayer,iaer)=1.e-25
513	ELSE ! part of the mesh with concentred dust storm
514	DO l=1,nlayer
515	IF (l.LE.cstdustlevel) THEN
516	c Opacity in the first levels is held constant to
517	c avoid unrealistic values due to constant lifting:
518	DO ig=1,ngrid
519	! OPTICAL DEPTH for the computation of tauref,
520	! which is to be compared with tauref_scenario
521	! => visible wavelength
522	aerosol(ig,l,iaer) =
523	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
524	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
525	& pq(ig,cstdustlevel,igcm_stormdust_mass) *
526	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
527	! DENSITY SCALED OPACITY :
528	! GCM output to be compared with observations
529	! => infrared wavelength
530	dsords(ig,l) =
531	& ( 0.75 * QREFir3d(ig,cstdustlevel,iaer) /
532	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
533	& pq(ig,cstdustlevel,igcm_stormdust_mass)
534	ENDDO
535	ELSE
536	DO ig=1,ngrid
537	! OPTICAL DEPTH for the computation of tauref,
538	! which is to be compared with tauref_scenario
539	! => visible wavelength
540	aerosol(ig,l,iaer) =
541	& ( 0.75 * QREFvis3d(ig,l,iaer) /
542	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
543	& pq(ig,l,igcm_stormdust_mass) *
544	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
545	! DENSITY SCALED OPACITY :
546	! GCM output to be compared with observations
547	! => infrared wavelength
548	dsords(ig,l) =
549	& ( 0.75 * QREFir3d(ig,l,iaer) /
550	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
551	& pq(ig,l,igcm_stormdust_mass)
552	ENDDO
553	ENDIF
554	ENDDO
555	ENDIF
556	c==================================================================
557	CASE("topdust_doubleq") aerkind ! MV18 : Two-moment scheme for
558	c topdust (transport of mass and number mixing ratio)
559	c==================================================================
560	c aerosol is calculated twice : once "above" the sub-grid mountain (nohmons=false)
561	c and once in the part of the mesh without the sub-grid mountain (nohmons=true)
562	aerosol(1:ngrid,1:nlayer,iaer) = 0.
563	IF (nohmons) THEN ! considering part of the mesh without storm
564	aerosol(1:ngrid,1:nlayer,iaer)=1.e-25
565	ELSE ! part of the mesh with concentred dust storm
566	DO l=1,nlayer
567	IF (l.LE.cstdustlevel) THEN
568	c Opacity in the first levels is held constant to
569	c avoid unrealistic values due to constant lifting:
570	DO ig=1,ngrid
571	aerosol(ig,l,iaer) =
572	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
573	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
574	& pq(ig,cstdustlevel,igcm_topdust_mass) *
575	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
576	! DENSITY SCALED OPACITY :
577	! GCM output to be compared with observations
578	! => infrared wavelength
579	dsotop(ig,l) =
580	& ( 0.75 * QREFir3d(ig,cstdustlevel,iaer) /
581	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
582	& pq(ig,cstdustlevel,igcm_topdust_mass)
583	ENDDO
584	ELSE
585	DO ig=1,ngrid
586	! OPTICAL DEPTH for the computation of tauref,
587	! which is to be compared with tauref_scenario
588	! => visible wavelength
589	aerosol(ig,l,iaer) =
590	& ( 0.75 * QREFvis3d(ig,l,iaer) /
591	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
592	& pq(ig,l,igcm_topdust_mass) *
593	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
594	! DENSITY SCALED OPACITY :
595	! GCM output to be compared with observations
596	! => infrared wavelength
597	dsotop(ig,l) =
598	& ( 0.75 * QREFir3d(ig,l,iaer) /
599	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
600	& pq(ig,l,igcm_topdust_mass)
601	ENDDO
602	ENDIF
603	ENDDO
604	ENDIF
605	c==================================================================
606	END SELECT aerkind
607	c -----------------------------------
608	ENDDO ! iaer (loop on aerosol kind)
609
610	c -----------------------------------------------------------------
611	c Rescaling each layer to reproduce the choosen (or assimilated)
612	c dust extinction opacity at visible reference wavelength, which
613	c is originally scaled to an equivalent odpref Pa pressure surface.
614	c -----------------------------------------------------------------
615
616
617	#ifdef DUSTSTORM
618	c -----------------------------------------------------------------
619	! Calculate reference opacity without perturbation
620	c -----------------------------------------------------------------
621	IF (firstcall) THEN
622	DO iaer=1,naerdust
623	DO l=1,nlayer
624	DO ig=1,ngrid
625	tauref(ig) = tauref(ig) +
626	& aerosol(ig,l,iaerdust(iaer))
627	ENDDO
628	ENDDO
629	ENDDO
630	tauref(:) = tauref(:) * odpref / pplev(:,1)
631
632	c--------------------------------------------------
633	c Get parameters of the opacity perturbation
634	c--------------------------------------------------
635	iaer=1 ! just change dust
636
637	write(,) "Add a local storm ?"
638	localstorm=.true. ! default value
639	call getin("localstorm",localstorm)
640	write(,) " localstorm = ",localstorm
641
642	IF (localstorm) THEN
643	WRITE(,) "********************"
644	WRITE(,) "ADDING A LOCAL STORM"
645	WRITE(,) "********************"
646
647	write(,) "ref opacity of local dust storm"
648	taulocref = 4.25 ! default value
649	call getin("taulocref",taulocref)
650	write(,) " taulocref = ",taulocref
651
652	write(,) "target altitude of local storm (km)"
653	ztoploc = 10.0 ! default value
654	call getin("ztoploc",ztoploc)
655	write(,) " ztoploc = ",ztoploc
656
657	write(,) "radius of dust storm (degree)"
658	radloc = 0.5 ! default value
659	call getin("radloc",radloc)
660	write(,) " radloc = ",radloc
661
662	write(,) "center longitude of storm (deg)"
663	lonloc = 25.0 ! default value
664	call getin("lonloc",lonloc)
665	write(,) " lonloc = ",lonloc
666
667	write(,) "center latitude of storm (deg)"
668	latloc = -2.5 ! default value
669	call getin("latloc",latloc)
670	write(,) " latloc = ",latloc
671
672	write(,) "reff storm (mic) 0. for background"
673	reffstorm = 0.0 ! default value
674	call getin("reffstorm",reffstorm)
675	write(,) " reffstorm = ",reffstorm
676
677	!! LOOP: modify opacity
678	DO ig=1,ngrid
679
680	!! distance to the center:
681	ray(ig)=SQRT((latitude(ig)180./pi-latloc)*2 +
682	& (longitude(ig)180./pi -lonloc)*2)
683
684	!! transition factor for storm
685	!! factor is hardcoded -- increase it to steepen
686	yeah = (TANH(2.+(radloc-ray(ig))*10.)+1.)/2.
687
688	!! new opacity field
689	!! -- add an opacity set to taulocref
690	!! -- the additional reference opacity will
691	!! thus be taulocref*odpref/pplev
692	tauuser(ig)=max( tauref(ig) * pplev(ig,1) /odpref ,
693	& taulocref * yeah )
694
695	!! compute l_top
696	DO l=1,nlayer
697	zalt(ig,l) = LOG( pplev(ig,1)/pplev(ig,l) )
698	& / g / 44.01
699	& * 8.31 * 210.
700	IF ( (ztoploc .lt. zalt(ig,l) )
701	& .and. (ztoploc .gt. zalt(ig,l-1)) ) l_top=l-1
702	ENDDO
703
704	!! change reffrad if ever needed
705	IF (reffstorm .gt. 0.) THEN
706	DO l=1,nlayer
707	IF (l .lt. l_top+1) THEN
708	reffrad(ig,l,iaer) = max( reffrad(ig,l,iaer), reffstorm
709	& * 1.e-6 * yeah )
710	ENDIF
711	ENDDO
712	ENDIF
713
714	ENDDO
715	!! END LOOP
716
717	!! compute perturbation in each layer (equation 8 in Spiga et al. JGR 2013)
718	DO ig=1,ngrid
719	int_factor(ig)=0.
720	DO l=1,nlayer
721	IF (l .lt. l_top+1) THEN
722	int_factor(ig) =
723	& int_factor(ig) +
724	& ( 0.75 * QREFvis3d(ig,l,iaer) /
725	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
726	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
727	ENDIF
728	ENDDO
729	DO l=1, nlayer
730	!! Mass mixing ratio perturbation due to local dust storm in each layer
731	more_dust(ig,l,1)=
732	& (tauuser(ig)-(tauref(ig)
733	& * pplev(ig,1) /odpref)) /
734	& int_factor(ig)
735	more_dust(ig,l,2)=
736	& (tauuser(ig)-(tauref(ig) *
737	& pplev(ig,1) /odpref))
738	& / int_factor(ig) *
739	& ((ref_r0/reffrad(ig,l,iaer))**3)
740	& * r3n_q
741	ENDDO
742	ENDDO
743
744	!! quantity of dust for each layer with the addition of the perturbation
745	DO l=1, l_top
746	pq(:,l,igcm_dust_mass)= pq(:,l,igcm_dust_mass)
747	. + more_dust(:,l,1)
748	pq(:,l,igcm_dust_number)= pq(:,l,igcm_dust_number)
749	. + more_dust(:,l,2)
750	ENDDO
751	ENDIF !! IF (localstorm)
752	tauref(:)=0.
753	ENDIF !! IF (firstcall)
754	#endif
755
756	IF (freedust) THEN
757	tauscaling(:) = 1.
758	c opacity obtained with stormdust
759	IF (rdstorm) THEN
760	taustormdusttmp(1:ngrid)=0.
761	DO l=1,nlayer
762	DO ig=1,ngrid
763	taustormdusttmp(ig) = taustormdusttmp(ig)+
764	& aerosol(ig,l,iaerdust(2))
765	ENDDO
766	ENDDO
767	!opacity obtained with background dust only
768	taubackdusttmp(1:ngrid)=0.
769	DO l=1,nlayer
770	DO ig=1,ngrid
771	taubackdusttmp(ig) = taubackdusttmp(ig)+
772	& aerosol(ig,l,iaerdust(1))
773	ENDDO
774	ENDDO
775	ENDIF !rdsstorm
776	ELSE
777	c Temporary scaling factor
778	taudusttmp(1:ngrid)=0.
779	DO iaer=1,naerdust
780	DO l=1,nlayer
781	DO ig=1,ngrid
782	c Scaling factor
783	taudusttmp(ig) = taudusttmp(ig) +
784	& aerosol(ig,l,iaerdust(iaer))
785	ENDDO
786	ENDDO
787	ENDDO
788
789	c Saved scaling factor
790	DO ig=1,ngrid
791	tauscaling(ig) = tauref(ig) *
792	& pplev(ig,1) / odpref / taudusttmp(ig)
793	ENDDO
794
795	ENDIF ! IF (freedust)
796
797	c Opacity computation
798	DO iaer=1,naerdust
799	DO l=1,nlayer
800	DO ig=1,ngrid
801	aerosol(ig,l,iaerdust(iaer)) = max(1E-20,
802	& aerosol(ig,l,iaerdust(iaer))* tauscaling(ig))
803	ENDDO
804	ENDDO
805	ENDDO
806
807	IF (freedust) THEN
808	! tauref has been initialized to 0 before.
809	DO iaer=1,naerdust
810	DO l=1,nlayer
811	DO ig=1,ngrid
812	#ifdef DUSTSTORM
813	!! recalculate opacity because storm perturbation has been added
814	IF (firstcall) THEN
815	aerosol(ig,l,iaer) =
816	& ( 0.75 * QREFvis3d(ig,l,iaer) /
817	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
818	& pq(ig,l,igcm_dust_mass) *
819	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
820	ENDIF
821	#endif
822	! tauref(ig) = tauref(ig) +
823	! & aerosol(ig,l,iaerdust(iaer))
824	c MV19: tauref must ALWAYS contain the opacity of all dust tracers
825	IF (name_iaer(iaerdust(iaer)).eq."dust_doubleq") THEN
826	tauref(ig) = tauref(ig) +
827	& ( 0.75 * QREFvis3d(ig,l,iaerdust(iaer)) /
828	& ( rho_dust * reffrad(ig,l,iaerdust(iaer)) ) ) *
829	& pq(ig,l,igcm_dust_mass) *
830	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
831	ELSE IF (name_iaer(iaerdust(iaer)).eq."stormdust_doubleq") THEN
832	tauref(ig) = tauref(ig) +
833	& ( 0.75 * QREFvis3d(ig,l,iaerdust(iaer)) /
834	& ( rho_dust * reffrad(ig,l,iaerdust(iaer)) ) ) *
835	& pq(ig,l,igcm_stormdust_mass) *
836	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
837	ELSE IF (name_iaer(iaerdust(iaer)).eq."topdust_doubleq") THEN
838	tauref(ig) = tauref(ig) +
839	& ( 0.75 * QREFvis3d(ig,l,iaerdust(iaer)) /
840	& ( rho_dust * reffrad(ig,l,iaerdust(iaer)) ) ) *
841	& pq(ig,l,igcm_topdust_mass) *
842	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
843	ENDIF
844
845	ENDDO
846	ENDDO
847	ENDDO
848	tauref(:) = tauref(:) * odpref / pplev(:,1)
849	ENDIF
850
851	c -----------------------------------------------------------------
852	c Column integrated visible optical depth in each point
853	c -----------------------------------------------------------------
854	DO iaer=1,naerkind
855	do l=1,nlayer
856	do ig=1,ngrid
857	tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer)
858	end do
859	end do
860	ENDDO
861
862	c for diagnostics: opacity for all dust scatterers stormdust included
863	taualldust(1:ngrid)=0.
864	DO iaer=1,naerdust
865	DO l=1,nlayer
866	DO ig=1,ngrid
867	taualldust(ig) = taualldust(ig) +
868	& aerosol(ig,l,iaerdust(iaer))
869	ENDDO
870	ENDDO
871	ENDDO
872
873	IF (rdstorm) THEN
874
875	c for diagnostics: opacity for dust in background only
876	taudust(1:ngrid)=0.
877	DO l=1,nlayer
878	DO ig=1,ngrid
879	taudust(ig) = taudust(ig) +
880	& aerosol(ig,l,iaer_dust_doubleq)
881	ENDDO
882	ENDDO
883
884	c for diagnostics: opacity for dust in storm only
885	taustormdust(1:ngrid)=0.
886	DO l=1,nlayer
887	DO ig=1,ngrid
888	taustormdust(ig) = taustormdust(ig) +
889	& aerosol(ig,l,iaer_stormdust_doubleq)
890	ENDDO
891	ENDDO
892
893	ENDIF
894
895
896	#ifdef DUSTSTORM
897	IF (firstcall) THEN
898	firstcall=.false.
899	ENDIF
900	#endif
901
902
903	c -----------------------------------------------------------------
904	c aerosol/X for stormdust to prepare calculation of radiative transfer
905	c -----------------------------------------------------------------
906	IF (rdstorm) THEN
907	DO l=1,nlayer
908	DO ig=1,ngrid
909	! stormdust: opacity relative to the storm fraction (stormdust/x)
910	aerosol(ig,l,iaer_stormdust_doubleq) =
911	& aerosol(ig,l,iaer_stormdust_doubleq)/totstormfract(ig)
912	ENDDO
913	ENDDO
914	ENDIF
915
916	c -----------------------------------------------------------------
917	c aerosol/X for topdust to prepare calculation of radiative transfer
918	c -----------------------------------------------------------------
919	IF (slpwind) THEN
920	DO ig=1,ngrid
921	IF (alpha_hmons(ig) .gt. 0.) THEN
922	DO l=1,nlayer
923	! topdust: opacity relative to the storm fraction (topdust/x)
924	aerosol(ig,l,iaer_topdust_doubleq) =
925	& aerosol(ig,l,iaer_topdust_doubleq)/alpha_hmons(ig)
926	ENDDO
927	ENDIF
928	ENDDO
929	ENDIF
930
931	END SUBROUTINE aeropacity
932
933	END MODULE aeropacity_mod

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

Download in other formats:

Original Format