Context Navigation

aeropacity.F @ 1635

Last change on this file since 1635 was 1543, checked in by emillour, 9 years ago

All models: Further adaptations to keep up with changes in LMDZ5 concerning
physics/dynamics separation:

dyn3d:
adapted gcm.F so that all physics initializations are now done in iniphysiq.

dyn3dpar:
adapted gcm.F so that all physics initializations are now done in iniphysiq.
updated calfis_p.F to follow up with changes.
copied over updated "bands.F90" from LMDZ5.

dynphy_lonlat:
calfis_p.F90, mod_interface_dyn_phys.F90, follow up of changes in phy_common/mod_* routines

phy_common:
added "geometry_mod.F90" to store information about the grid (replaces phy*/comgeomphy.F90) and give variables friendlier names: rlond => longitude , rlatd => latitude, airephy => cell_area, cuphy => dx , cvphy => dy
added "physics_distribution_mod.F90"
updated "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_mpi_data.F90", "mod_phys_lmdz_para.F90", "mod_phys_lmdz_mpi_transfert.F90", "mod_grid_phy_lmdz.F90", "mod_phys_lmdz_omp_data.F90", "mod_phys_lmdz_omp_transfert.F90", "write_field_phy.F90" and "ioipsl_getin_p_mod.F90" to LMDZ5 versions.

phy[venus/titan/mars/std]:
removed "init_phys_lmdz.F90", "comgeomphy.F90"; adapted routines to use geometry_mod (longitude, latitude, cell_area, etc.)

File size: 25.9 KB

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

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

Download in other formats:

Original Format