Context Navigation

aeropacity_mod.F @ 2266

Last change on this file since 2266 was 2252, checked in by abierjon, 5 years ago
Mars GCM: Bug fix following r2248 in aeropacity_mod and topmons_mod : since dsodust, dsords and dsotop are diagnostic physiq_mod variables, we don't want them to be reinitialized at each call of aeropacity_mod and topmons_mod, but we initialize them once and for all at the beginning of physiq_mod instead. AB
File size: 34.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,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(INOUT) :: dsodust(ngrid,nlayer)
92	REAL, INTENT(INOUT) :: dsords(ngrid,nlayer) !dso of stormdust
93	REAL, INTENT(INOUT) :: 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	! initializations
172	tau(1:ngrid,1:naerkind)=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 ! end of if firstcall
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 background 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	! OPTICAL DEPTH for the computation of tauref,
388	! which is to be compared with tauref_scenario
389	! => visible wavelength
390	aerosol(ig,l,iaer) =
391	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
392	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
393	& pq(ig,cstdustlevel,igcm_dust_mass) *
394	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
395	! DENSITY SCALED OPACITY :
396	! GCM output to be compared with observations
397	! => infrared wavelength
398	dsodust(ig,l) =
399	& ( 0.75 * QREFir3d(ig,cstdustlevel,iaer) /
400	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
401	& pq(ig,cstdustlevel,igcm_dust_mass)
402	ENDDO
403	ELSE
404	DO ig=1,ngrid
405	! OPTICAL DEPTH for the computation of tauref,
406	! which is to be compared with tauref_scenario
407	! => visible wavelength
408	aerosol(ig,l,iaer) =
409	& ( 0.75 * QREFvis3d(ig,l,iaer) /
410	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
411	& pq(ig,l,igcm_dust_mass) *
412	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
413	! DENSITY SCALED OPACITY :
414	! GCM output to be compared with observations
415	! => infrared wavelength
416	dsodust(ig,l) =
417	& ( 0.75 * QREFir3d(ig,l,iaer) /
418	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
419	& pq(ig,l,igcm_dust_mass)
420	ENDDO
421	ENDIF
422	ENDDO
423
424	c==================================================================
425	CASE("dust_submicron") aerkind ! Small dust population
426	c==================================================================
427
428	DO l=1,nlayer
429	IF (l.LE.cstdustlevel) THEN
430	c Opacity in the first levels is held constant to
431	c avoid unrealistic values due to constant lifting:
432	DO ig=1,ngrid
433	aerosol(ig,l,iaer) =
434	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
435	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
436	& pq(ig,cstdustlevel,igcm_dust_submicron) *
437	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
438	ENDDO
439	ELSE
440	DO ig=1,ngrid
441	aerosol(ig,l,iaer) =
442	& ( 0.75 * QREFvis3d(ig,l,iaer) /
443	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
444	& pq(ig,l,igcm_dust_submicron) *
445	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
446	ENDDO
447	ENDIF
448	ENDDO
449
450	c==================================================================
451	CASE("h2o_ice") aerkind ! Water ice crystals
452	c==================================================================
453
454	c 1. Initialization
455	aerosol(1:ngrid,1:nlayer,iaer) = 0.
456	taucloudvis(1:ngrid) = 0.
457	taucloudtes(1:ngrid) = 0.
458	c 2. Opacity calculation
459	! NO CLOUDS
460	IF (clearsky) THEN
461	aerosol(1:ngrid,1:nlayer,iaer) =1.e-9
462	! CLOUDSs
463	ELSE ! else (clearsky)
464	DO ig=1, ngrid
465	DO l=1,nlayer
466	aerosol(ig,l,iaer) = max(1E-20,
467	& ( 0.75 * QREFvis3d(ig,l,iaer) /
468	& ( rho_ice * reffrad(ig,l,iaer) ) ) *
469	& pq(ig,l,i_ice) *
470	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
471	& )
472	taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
473	taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
474	& QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
475	& ( 1.E0 - omegaREFir3d(ig,l,iaer) )
476	ENDDO
477	ENDDO
478	! SUB-GRID SCALE CLOUDS
479	IF (CLFvarying) THEN
480	DO ig=1, ngrid
481	DO l=1,nlayer-1
482	CLFtot = max(totcloudfrac(ig),0.01)
483	aerosol(ig,l,iaer)=
484	& aerosol(ig,l,iaer)/CLFtot
485	aerosol(ig,l,iaer) =
486	& max(aerosol(ig,l,iaer),1.e-9)
487	ENDDO
488	ENDDO
489	! ELSE ! else (CLFvarying)
490	! DO ig=1, ngrid
491	! DO l=1,nlayer-1 ! to stop the rad tran bug
492	! CLFtot = CLFfixval
493	! aerosol(ig,l,iaer)=
494	! & aerosol(ig,l,iaer)/CLFtot
495	! aerosol(ig,l,iaer) =
496	! & max(aerosol(ig,l,iaer),1.e-9)
497	! ENDDO
498	! ENDDO
499	ENDIF ! end (CLFvarying)
500	ENDIF ! end (clearsky)
501
502	c==================================================================
503	CASE("stormdust_doubleq") aerkind ! CW17 : Two-moment scheme for
504	c stormdust (transport of mass and number mixing ratio)
505	c==================================================================
506	c aerosol is calculated twice : once within the dust storm (clearatm=false)
507	c and once in the part of the mesh without dust storm (clearatm=true)
508	aerosol(1:ngrid,1:nlayer,iaer) = 0.
509	IF (clearatm) THEN ! considering part of the mesh without storm
510	aerosol(1:ngrid,1:nlayer,iaer)=1.e-25
511	ELSE ! part of the mesh with concentred dust storm
512	DO l=1,nlayer
513	IF (l.LE.cstdustlevel) THEN
514	c Opacity in the first levels is held constant to
515	c avoid unrealistic values due to constant lifting:
516	DO ig=1,ngrid
517	! OPTICAL DEPTH for the computation of tauref,
518	! which is to be compared with tauref_scenario
519	! => visible wavelength
520	aerosol(ig,l,iaer) =
521	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
522	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
523	& pq(ig,cstdustlevel,igcm_stormdust_mass) *
524	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
525	! DENSITY SCALED OPACITY :
526	! GCM output to be compared with observations
527	! => infrared wavelength
528	dsords(ig,l) =
529	& ( 0.75 * QREFir3d(ig,cstdustlevel,iaer) /
530	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
531	& pq(ig,cstdustlevel,igcm_stormdust_mass)
532	ENDDO
533	ELSE
534	DO ig=1,ngrid
535	! OPTICAL DEPTH for the computation of tauref,
536	! which is to be compared with tauref_scenario
537	! => visible wavelength
538	aerosol(ig,l,iaer) =
539	& ( 0.75 * QREFvis3d(ig,l,iaer) /
540	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
541	& pq(ig,l,igcm_stormdust_mass) *
542	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
543	! DENSITY SCALED OPACITY :
544	! GCM output to be compared with observations
545	! => infrared wavelength
546	dsords(ig,l) =
547	& ( 0.75 * QREFir3d(ig,l,iaer) /
548	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
549	& pq(ig,l,igcm_stormdust_mass)
550	ENDDO
551	ENDIF
552	ENDDO
553	ENDIF
554	c==================================================================
555	CASE("topdust_doubleq") aerkind ! MV18 : Two-moment scheme for
556	c topdust (transport of mass and number mixing ratio)
557	c==================================================================
558	c aerosol is calculated twice : once "above" the sub-grid mountain (nohmons=false)
559	c and once in the part of the mesh without the sub-grid mountain (nohmons=true)
560	aerosol(1:ngrid,1:nlayer,iaer) = 0.
561	IF (nohmons) THEN ! considering part of the mesh without storm
562	aerosol(1:ngrid,1:nlayer,iaer)=1.e-25
563	ELSE ! part of the mesh with concentred dust storm
564	DO l=1,nlayer
565	IF (l.LE.cstdustlevel) THEN
566	c Opacity in the first levels is held constant to
567	c avoid unrealistic values due to constant lifting:
568	DO ig=1,ngrid
569	aerosol(ig,l,iaer) =
570	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
571	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
572	& pq(ig,cstdustlevel,igcm_topdust_mass) *
573	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
574	! DENSITY SCALED OPACITY :
575	! GCM output to be compared with observations
576	! => infrared wavelength
577	dsotop(ig,l) =
578	& ( 0.75 * QREFir3d(ig,cstdustlevel,iaer) /
579	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
580	& pq(ig,cstdustlevel,igcm_topdust_mass)
581	ENDDO
582	ELSE
583	DO ig=1,ngrid
584	! OPTICAL DEPTH for the computation of tauref,
585	! which is to be compared with tauref_scenario
586	! => visible wavelength
587	aerosol(ig,l,iaer) =
588	& ( 0.75 * QREFvis3d(ig,l,iaer) /
589	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
590	& pq(ig,l,igcm_topdust_mass) *
591	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
592	! DENSITY SCALED OPACITY :
593	! GCM output to be compared with observations
594	! => infrared wavelength
595	dsotop(ig,l) =
596	& ( 0.75 * QREFir3d(ig,l,iaer) /
597	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
598	& pq(ig,l,igcm_topdust_mass)
599	ENDDO
600	ENDIF
601	ENDDO
602	ENDIF
603	c==================================================================
604	END SELECT aerkind
605	c -----------------------------------
606	ENDDO ! iaer (loop on aerosol kind)
607
608	c -----------------------------------------------------------------
609	c Rescaling each layer to reproduce the choosen (or assimilated)
610	c dust extinction opacity at visible reference wavelength, which
611	c is originally scaled to an equivalent odpref Pa pressure surface.
612	c -----------------------------------------------------------------
613
614
615	#ifdef DUSTSTORM
616	c -----------------------------------------------------------------
617	! Calculate reference opacity without perturbation
618	c -----------------------------------------------------------------
619	IF (firstcall) THEN
620	DO iaer=1,naerdust
621	DO l=1,nlayer
622	DO ig=1,ngrid
623	tauref(ig) = tauref(ig) +
624	& aerosol(ig,l,iaerdust(iaer))
625	ENDDO
626	ENDDO
627	ENDDO
628	tauref(:) = tauref(:) * odpref / pplev(:,1)
629
630	c--------------------------------------------------
631	c Get parameters of the opacity perturbation
632	c--------------------------------------------------
633	iaer=1 ! just change dust
634
635	write(,) "Add a local storm ?"
636	localstorm=.true. ! default value
637	call getin("localstorm",localstorm)
638	write(,) " localstorm = ",localstorm
639
640	IF (localstorm) THEN
641	WRITE(,) "********************"
642	WRITE(,) "ADDING A LOCAL STORM"
643	WRITE(,) "********************"
644
645	write(,) "ref opacity of local dust storm"
646	taulocref = 4.25 ! default value
647	call getin("taulocref",taulocref)
648	write(,) " taulocref = ",taulocref
649
650	write(,) "target altitude of local storm (km)"
651	ztoploc = 10.0 ! default value
652	call getin("ztoploc",ztoploc)
653	write(,) " ztoploc = ",ztoploc
654
655	write(,) "radius of dust storm (degree)"
656	radloc = 0.5 ! default value
657	call getin("radloc",radloc)
658	write(,) " radloc = ",radloc
659
660	write(,) "center longitude of storm (deg)"
661	lonloc = 25.0 ! default value
662	call getin("lonloc",lonloc)
663	write(,) " lonloc = ",lonloc
664
665	write(,) "center latitude of storm (deg)"
666	latloc = -2.5 ! default value
667	call getin("latloc",latloc)
668	write(,) " latloc = ",latloc
669
670	write(,) "reff storm (mic) 0. for background"
671	reffstorm = 0.0 ! default value
672	call getin("reffstorm",reffstorm)
673	write(,) " reffstorm = ",reffstorm
674
675	!! LOOP: modify opacity
676	DO ig=1,ngrid
677
678	!! distance to the center:
679	ray(ig)=SQRT((latitude(ig)180./pi-latloc)*2 +
680	& (longitude(ig)180./pi -lonloc)*2)
681
682	!! transition factor for storm
683	!! factor is hardcoded -- increase it to steepen
684	yeah = (TANH(2.+(radloc-ray(ig))*10.)+1.)/2.
685
686	!! new opacity field
687	!! -- add an opacity set to taulocref
688	!! -- the additional reference opacity will
689	!! thus be taulocref*odpref/pplev
690	tauuser(ig)=max( tauref(ig) * pplev(ig,1) /odpref ,
691	& taulocref * yeah )
692
693	!! compute l_top
694	DO l=1,nlayer
695	zalt(ig,l) = LOG( pplev(ig,1)/pplev(ig,l) )
696	& / g / 44.01
697	& * 8.31 * 210.
698	IF ( (ztoploc .lt. zalt(ig,l) )
699	& .and. (ztoploc .gt. zalt(ig,l-1)) ) l_top=l-1
700	ENDDO
701
702	!! change reffrad if ever needed
703	IF (reffstorm .gt. 0.) THEN
704	DO l=1,nlayer
705	IF (l .lt. l_top+1) THEN
706	reffrad(ig,l,iaer) = max( reffrad(ig,l,iaer), reffstorm
707	& * 1.e-6 * yeah )
708	ENDIF
709	ENDDO
710	ENDIF
711
712	ENDDO
713	!! END LOOP
714
715	!! compute perturbation in each layer (equation 8 in Spiga et al. JGR 2013)
716	DO ig=1,ngrid
717	int_factor(ig)=0.
718	DO l=1,nlayer
719	IF (l .lt. l_top+1) THEN
720	int_factor(ig) =
721	& int_factor(ig) +
722	& ( 0.75 * QREFvis3d(ig,l,iaer) /
723	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
724	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
725	ENDIF
726	ENDDO
727	DO l=1, nlayer
728	!! Mass mixing ratio perturbation due to local dust storm in each layer
729	more_dust(ig,l,1)=
730	& (tauuser(ig)-(tauref(ig)
731	& * pplev(ig,1) /odpref)) /
732	& int_factor(ig)
733	more_dust(ig,l,2)=
734	& (tauuser(ig)-(tauref(ig) *
735	& pplev(ig,1) /odpref))
736	& / int_factor(ig) *
737	& ((ref_r0/reffrad(ig,l,iaer))**3)
738	& * r3n_q
739	ENDDO
740	ENDDO
741
742	!! quantity of dust for each layer with the addition of the perturbation
743	DO l=1, l_top
744	pq(:,l,igcm_dust_mass)= pq(:,l,igcm_dust_mass)
745	. + more_dust(:,l,1)
746	pq(:,l,igcm_dust_number)= pq(:,l,igcm_dust_number)
747	. + more_dust(:,l,2)
748	ENDDO
749	ENDIF !! IF (localstorm)
750	tauref(:)=0.
751	ENDIF !! IF (firstcall)
752	#endif
753
754	IF (freedust) THEN
755	tauscaling(:) = 1.
756	c opacity obtained with stormdust
757	IF (rdstorm) THEN
758	taustormdusttmp(1:ngrid)=0.
759	DO l=1,nlayer
760	DO ig=1,ngrid
761	taustormdusttmp(ig) = taustormdusttmp(ig)+
762	& aerosol(ig,l,iaerdust(2))
763	ENDDO
764	ENDDO
765	!opacity obtained with background dust only
766	taubackdusttmp(1:ngrid)=0.
767	DO l=1,nlayer
768	DO ig=1,ngrid
769	taubackdusttmp(ig) = taubackdusttmp(ig)+
770	& aerosol(ig,l,iaerdust(1))
771	ENDDO
772	ENDDO
773	ENDIF !rdsstorm
774	ELSE
775	c Temporary scaling factor
776	taudusttmp(1:ngrid)=0.
777	DO iaer=1,naerdust
778	DO l=1,nlayer
779	DO ig=1,ngrid
780	c Scaling factor
781	taudusttmp(ig) = taudusttmp(ig) +
782	& aerosol(ig,l,iaerdust(iaer))
783	ENDDO
784	ENDDO
785	ENDDO
786
787	c Saved scaling factor
788	DO ig=1,ngrid
789	tauscaling(ig) = tauref(ig) *
790	& pplev(ig,1) / odpref / taudusttmp(ig)
791	ENDDO
792
793	ENDIF ! IF (freedust)
794
795	c Opacity computation
796	DO iaer=1,naerdust
797	DO l=1,nlayer
798	DO ig=1,ngrid
799	aerosol(ig,l,iaerdust(iaer)) = max(1E-20,
800	& aerosol(ig,l,iaerdust(iaer))* tauscaling(ig))
801	ENDDO
802	ENDDO
803	ENDDO
804
805	IF (freedust) THEN
806	! tauref has been initialized to 0 before.
807	DO iaer=1,naerdust
808	DO l=1,nlayer
809	DO ig=1,ngrid
810	#ifdef DUSTSTORM
811	!! recalculate opacity because storm perturbation has been added
812	IF (firstcall) THEN
813	aerosol(ig,l,iaer) =
814	& ( 0.75 * QREFvis3d(ig,l,iaer) /
815	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
816	& pq(ig,l,igcm_dust_mass) *
817	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
818	ENDIF
819	#endif
820	! tauref(ig) = tauref(ig) +
821	! & aerosol(ig,l,iaerdust(iaer))
822	c MV19: tauref must ALWAYS contain the opacity of all dust tracers
823	IF (name_iaer(iaerdust(iaer)).eq."dust_doubleq") THEN
824	tauref(ig) = tauref(ig) +
825	& ( 0.75 * QREFvis3d(ig,l,iaerdust(iaer)) /
826	& ( rho_dust * reffrad(ig,l,iaerdust(iaer)) ) ) *
827	& pq(ig,l,igcm_dust_mass) *
828	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
829	ELSE IF (name_iaer(iaerdust(iaer)).eq."stormdust_doubleq") THEN
830	tauref(ig) = tauref(ig) +
831	& ( 0.75 * QREFvis3d(ig,l,iaerdust(iaer)) /
832	& ( rho_dust * reffrad(ig,l,iaerdust(iaer)) ) ) *
833	& pq(ig,l,igcm_stormdust_mass) *
834	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
835	ELSE IF (name_iaer(iaerdust(iaer)).eq."topdust_doubleq") THEN
836	tauref(ig) = tauref(ig) +
837	& ( 0.75 * QREFvis3d(ig,l,iaerdust(iaer)) /
838	& ( rho_dust * reffrad(ig,l,iaerdust(iaer)) ) ) *
839	& pq(ig,l,igcm_topdust_mass) *
840	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
841	ENDIF
842
843	ENDDO
844	ENDDO
845	ENDDO
846	tauref(:) = tauref(:) * odpref / pplev(:,1)
847	ENDIF
848
849	c -----------------------------------------------------------------
850	c Column integrated visible optical depth in each point
851	c -----------------------------------------------------------------
852	DO iaer=1,naerkind
853	do l=1,nlayer
854	do ig=1,ngrid
855	tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer)
856	end do
857	end do
858	ENDDO
859
860	c for diagnostics: opacity for all dust scatterers stormdust included
861	taualldust(1:ngrid)=0.
862	DO iaer=1,naerdust
863	DO l=1,nlayer
864	DO ig=1,ngrid
865	taualldust(ig) = taualldust(ig) +
866	& aerosol(ig,l,iaerdust(iaer))
867	ENDDO
868	ENDDO
869	ENDDO
870
871	IF (rdstorm) THEN
872
873	c for diagnostics: opacity for dust in background only
874	taudust(1:ngrid)=0.
875	DO l=1,nlayer
876	DO ig=1,ngrid
877	taudust(ig) = taudust(ig) +
878	& aerosol(ig,l,iaer_dust_doubleq)
879	ENDDO
880	ENDDO
881
882	c for diagnostics: opacity for dust in storm only
883	taustormdust(1:ngrid)=0.
884	DO l=1,nlayer
885	DO ig=1,ngrid
886	taustormdust(ig) = taustormdust(ig) +
887	& aerosol(ig,l,iaer_stormdust_doubleq)
888	ENDDO
889	ENDDO
890
891	ENDIF
892
893
894	#ifdef DUSTSTORM
895	IF (firstcall) THEN
896	firstcall=.false.
897	ENDIF
898	#endif
899
900
901	c -----------------------------------------------------------------
902	c aerosol/X for stormdust to prepare calculation of radiative transfer
903	c -----------------------------------------------------------------
904	IF (rdstorm) THEN
905	DO l=1,nlayer
906	DO ig=1,ngrid
907	! stormdust: opacity relative to the storm fraction (stormdust/x)
908	aerosol(ig,l,iaer_stormdust_doubleq) =
909	& aerosol(ig,l,iaer_stormdust_doubleq)/totstormfract(ig)
910	ENDDO
911	ENDDO
912	ENDIF
913
914	c -----------------------------------------------------------------
915	c aerosol/X for topdust to prepare calculation of radiative transfer
916	c -----------------------------------------------------------------
917	IF (slpwind) THEN
918	DO ig=1,ngrid
919	IF (alpha_hmons(ig) .gt. 0.) THEN
920	DO l=1,nlayer
921	! topdust: opacity relative to the storm fraction (topdust/x)
922	aerosol(ig,l,iaer_topdust_doubleq) =
923	& aerosol(ig,l,iaer_topdust_doubleq)/alpha_hmons(ig)
924	ENDDO
925	ENDIF
926	ENDDO
927	ENDIF
928
929	END SUBROUTINE aeropacity
930
931	END MODULE aeropacity_mod

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

Download in other formats:

Original Format