Context Navigation

aeropacity_mod.F @ 2325

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

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

Download in other formats:

Original Format