Context Navigation

aeropacity_mod.F @ 2156

Last change on this file since 2156 was 2137, checked in by emillour, 6 years ago
Mars GCM: Updates in code to be able to read in the MY34 dust scenario and the MY34 solar EUV input. EM
File size: 30.5 KB

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

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

Download in other formats:

Original Format