Context Navigation

aeropacity_mod.F @ 1769

Last change on this file since 1769 was 1711, checked in by mvals, 8 years ago

Mars GCM

implementation of a sub-grid water cloud fraction scheme (by A. Pottier): If CLFvarying is set to true in callphys.def, then the sub-grid cloud fraction routine is applied in watercloud.F and aeropacity.F.
accordingly modified files are: watercloud.F, aeropacity.F, callradite.F, conf_phys.F, phyetat0.F, phyredem.F90, physiq_mod.F, testphys1d.F, callkeys.h, newtsart.F, start2archive.F, lect_start_archive.F
added file: tcondwater.F90, used by watercloud.F to calculate the condensation temperature of water
watercloud.F, aeropacity.F, callradite.F are converted to module files as watercloud_mod.F, aeropacity_mod.F, callradite_mod.F

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

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

Download in other formats:

Original Format