Context Navigation

aeropacity.F90 @ 2254

Last change on this file since 2254 was 2254, checked in by jvatant, 5 years ago

Fix a bug in the 2-layers aerosol scheme in aeropacity.F90
-> It appeared that the normalization for total column optical depth was a

wrong way to go, the actual strato and tropo values weren't matcihing the
input values, and one of the main consequence was to strongly underestimate
the actual strato optical depth compared to input value in most cases.

-> Now we ensure the normalization to the input values for both layers, while

keeping a dependence in dP for each GCM layer.

-> Also add a sanity check ensuring pres_top_tropo > pres_bottom_strato.
--JVO

File size: 22.2 KB

Line
1	Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pq, &
2	aerosol,reffrad,QREFvis3d,QREFir3d,tau_col,cloudfrac,totcloudfrac,clearsky)
3
4	use radinc_h, only : L_TAUMAX,naerkind
5	use aerosol_mod
6	USE tracer_h, only: noms,rho_co2,rho_ice
7	use comcstfi_mod, only: g, pi
8	use geometry_mod, only: latitude
9	use callkeys_mod, only: aerofixco2,aerofixh2o,kastprof,cloudlvl, &
10	CLFvarying,CLFfixval,dusttau, &
11	pres_bottom_tropo,pres_top_tropo,obs_tau_col_tropo, &
12	pres_bottom_strato,pres_top_strato,obs_tau_col_strato, &
13	tau_nh3_cloud, pres_nh3_cloud
14
15	implicit none
16
17	!==================================================================
18	!
19	! Purpose
20	! -------
21	! Compute aerosol optical depth in each gridbox.
22	!
23	! Authors
24	! -------
25	! F. Forget
26	! F. Montmessin (water ice scheme)
27	! update J.-B. Madeleine (2008)
28	! dust removal, simplification by Robin Wordsworth (2009)
29	!
30	! Input
31	! -----
32	! ngrid Number of horizontal gridpoints
33	! nlayer Number of layers
34	! nq Number of tracers
35	! pplev Pressure (Pa) at each layer boundary
36	! pq Aerosol mixing ratio
37	! reffrad(ngrid,nlayer,naerkind) Aerosol effective radius
38	! QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
39	! QREFir3d(ngrid,nlayer,naerkind) / at reference wavelengths
40	!
41	! Output
42	! ------
43	! aerosol Aerosol optical depth in layer l, grid point ig
44	! tau_col Total column optical depth at grid point ig
45	!
46	!=======================================================================
47
48	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
49	INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
50	INTEGER,INTENT(IN) :: nq ! number of tracers
51	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
52	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
53	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
54	REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth
55	REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) ! aerosol effective radius
56	REAL,INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind) ! extinction coefficient in the visible
57	REAL,INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
58	REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth
59	! BENJAMIN MODIFS
60	real,intent(in) :: cloudfrac(ngrid,nlayer) ! cloud fraction
61	real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
62	logical,intent(in) :: clearsky
63
64	real aerosol0, obs_tau_col_aurora, pm, pente_cloud
65
66	real dp_strato(ngrid)
67	real dp_tropo(ngrid)
68
69	INTEGER l,ig,iq,iaer
70
71	LOGICAL,SAVE :: firstcall=.true.
72	!$OMP THREADPRIVATE(firstcall)
73	REAL CBRT
74	EXTERNAL CBRT
75
76	INTEGER,SAVE :: i_co2ice=0 ! co2 ice
77	INTEGER,SAVE :: i_h2oice=0 ! water ice
78	!$OMP THREADPRIVATE(i_co2ice,i_h2oice)
79	CHARACTER(LEN=20) :: tracername ! to temporarily store text
80
81	! for fixed dust profiles
82	real topdust, expfactor, zp
83	REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
84	REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling
85
86	real CLFtot
87
88	! identify tracers
89	IF (firstcall) THEN
90
91	write(,) "Tracers found in aeropacity:"
92	do iq=1,nq
93	tracername=noms(iq)
94	if (tracername.eq."co2_ice") then
95	i_co2ice=iq
96	write(,) "i_co2ice=",i_co2ice
97
98	endif
99	if (tracername.eq."h2o_ice") then
100	i_h2oice=iq
101	write(,) "i_h2oice=",i_h2oice
102	endif
103	enddo
104
105	if (noaero) then
106	print*, "No active aerosols found in aeropacity"
107	else
108	print*, "If you would like to use aerosols, make sure any old"
109	print*, "start files are updated in newstart using the option"
110	print*, "q=0"
111	write(,) "Active aerosols found in aeropacity:"
112	endif
113
114	if ((iaero_co2.ne.0).and.(.not.noaero)) then
115	print*, 'iaero_co2= ',iaero_co2
116	endif
117	if (iaero_h2o.ne.0) then
118	print*,'iaero_h2o= ',iaero_h2o
119	endif
120	if (iaero_dust.ne.0) then
121	print*,'iaero_dust= ',iaero_dust
122	endif
123	if (iaero_h2so4.ne.0) then
124	print*,'iaero_h2so4= ',iaero_h2so4
125	endif
126	if (iaero_back2lay.ne.0) then
127	print*,'iaero_back2lay= ',iaero_back2lay
128	endif
129	if (iaero_nh3.ne.0) then
130	print*,'iaero_nh3= ',iaero_nh3
131	endif
132	if (iaero_aurora.ne.0) then
133	print*,'iaero_aurora= ',iaero_aurora
134	endif
135
136	firstcall=.false.
137	ENDIF ! of IF (firstcall)
138
139
140	! ---------------------------------------------------------
141	!==================================================================
142	! CO2 ice aerosols
143	!==================================================================
144
145	if (iaero_co2.ne.0) then
146	iaer=iaero_co2
147	! 1. Initialization
148	aerosol(1:ngrid,1:nlayer,iaer)=0.0
149	! 2. Opacity calculation
150	if (noaero) then ! aerosol set to zero
151	aerosol(1:ngrid,1:nlayer,iaer)=0.0
152	elseif (aerofixco2.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed
153	aerosol(1:ngrid,1:nlayer,iaer)=1.e-9
154	!aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option
155	else
156	DO ig=1, ngrid
157	DO l=1,nlayer-1 ! to stop the rad tran bug
158
159	aerosol0 = &
160	( 0.75 * QREFvis3d(ig,l,iaer) / &
161	( rho_co2 * reffrad(ig,l,iaer) ) ) * &
162	( pq(ig,l,i_co2ice) + 1.E-9 ) * &
163	( pplev(ig,l) - pplev(ig,l+1) ) / g
164	aerosol0 = max(aerosol0,1.e-9)
165	aerosol0 = min(aerosol0,L_TAUMAX)
166	aerosol(ig,l,iaer) = aerosol0
167	! aerosol(ig,l,iaer) = 0.0
168	! print*, aerosol(ig,l,iaer)
169	! using cloud fraction
170	! aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF))
171	! aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX)
172
173
174	ENDDO
175	ENDDO
176	end if ! if fixed or varying
177	end if ! if CO2 aerosols
178	!==================================================================
179	! Water ice / liquid
180	!==================================================================
181
182	if (iaero_h2o.ne.0) then
183	iaer=iaero_h2o
184	! 1. Initialization
185	aerosol(1:ngrid,1:nlayer,iaer)=0.0
186	! 2. Opacity calculation
187	if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
188	aerosol(1:ngrid,1:nlayer,iaer) =1.e-9
189
190	! put cloud at cloudlvl
191	if(kastprof.and.(cloudlvl.ne.0.0))then
192	ig=1
193	do l=1,nlayer
194	if(int(cloudlvl).eq.l)then
195	!if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
196	print*,'Inserting cloud at level ',l
197	!aerosol(ig,l,iaer)=10.0
198
199	rho_ice=920.0
200
201	! the Kasting approximation
202	aerosol(ig,l,iaer) = &
203	( 0.75 * QREFvis3d(ig,l,iaer) / &
204	( rho_ice * reffrad(ig,l,iaer) ) ) * &
205	!( pq(ig,l,i_h2oice) + 1.E-9 ) * &
206	( 4.0e-4 + 1.E-9 ) * &
207	( pplev(ig,l) - pplev(ig,l+1) ) / g
208
209
210	open(115,file='clouds.out',form='formatted')
211	write(115,*) l,aerosol(ig,l,iaer)
212	close(115)
213
214	return
215	endif
216	end do
217
218	call abort
219	endif
220
221	else
222
223	do ig=1, ngrid
224	!do l=1,nlayer-1 ! to stop the rad tran bug
225	do l=1,nlayer !JL18 if aerosols are present in the last layer we must account for them. Provides better upper boundary condition in the IR. They must however be put to zero in the sw (see optcv)
226	! same correction should b-probably be done for other aerosol types.
227	aerosol(ig,l,iaer) = & !modification by BC
228	( 0.75 * QREFvis3d(ig,l,iaer) / &
229	( rho_ice * reffrad(ig,l,iaer) ) ) * &
230	! pq(ig,l,i_h2oice) * & !JL I dropped the +1e-9 here to have the same
231	!( pplev(ig,l) - pplev(ig,l+1) ) / g ! opacity in the clearsky=true and the
232	! clear=false/pq=0 case
233	( pq(ig,l,i_h2oice) + 1.E-9 ) * & ! Doing this makes the code unstable, so I have restored it (RW)
234	( pplev(ig,l) - pplev(ig,l+1) ) / g
235
236	enddo
237	enddo
238
239	if(CLFvarying)then
240	call totalcloudfrac(ngrid,nlayer,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1,1,iaer))
241	do ig=1, ngrid
242	!do l=1,nlayer-1 ! to stop the rad tran bug
243	do l=1,nlayer !JL18 if aerosols are present in the last layer we must account for them. Provides better upper boundary condition in the IR. They must however be put to zero in the sw (see optcv)
244	CLFtot = max(totcloudfrac(ig),0.01)
245	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
246	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
247	enddo
248	enddo
249	else
250	do ig=1, ngrid
251	!do l=1,nlayer-1 ! to stop the rad tran bug
252	do l=1,nlayer !JL18 if aerosols are present in the last layer we must account for them. Provides better upper boundary condition in the IR. They must however be put to zero in the sw (see optcv)
253	CLFtot = CLFfixval
254	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
255	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
256	enddo
257	enddo
258	end if!(CLFvarying)
259	endif !(aerofixed.or.(i_h2oice.eq.0).or.clearsky)
260
261	end if ! End if h2o aerosol
262
263	!==================================================================
264	! Dust
265	!==================================================================
266	if (iaero_dust.ne.0) then
267	iaer=iaero_dust
268	! 1. Initialization
269	aerosol(1:ngrid,1:nlayer,iaer)=0.0
270
271	topdust=30.0 ! km (used to be 10.0 km) LK
272
273	! 2. Opacity calculation
274
275	! expfactor=0.
276	DO l=1,nlayer-1
277	DO ig=1,ngrid
278	! Typical mixing ratio profile
279
280	zp=(pplev(ig,1)/pplay(ig,l))**(70./topdust)
281	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
282
283	! Vertical scaling function
284	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
285	*expfactor
286
287
288	ENDDO
289	ENDDO
290
291	! Rescaling each layer to reproduce the choosen (or assimilated)
292	! dust extinction opacity at visible reference wavelength, which
293	! is scaled to the surface pressure pplev(ig,1)
294
295	taudusttmp(1:ngrid)=0.
296	DO l=1,nlayer
297	DO ig=1,ngrid
298	taudusttmp(ig) = taudusttmp(ig) &
299	+ aerosol(ig,l,iaer)
300	ENDDO
301	ENDDO
302	DO l=1,nlayer-1
303	DO ig=1,ngrid
304	aerosol(ig,l,iaer) = max(1E-20, &
305	dusttau &
306	* pplev(ig,1) / pplev(ig,1) &
307	* aerosol(ig,l,iaer) &
308	/ taudusttmp(ig))
309
310	ENDDO
311	ENDDO
312	end if ! If dust aerosol
313
314	!==================================================================
315	! H2SO4
316	!==================================================================
317	! added by LK
318	if (iaero_h2so4.ne.0) then
319	iaer=iaero_h2so4
320
321	! 1. Initialization
322	aerosol(1:ngrid,1:nlayer,iaer)=0.0
323
324
325	! 2. Opacity calculation
326
327	! expfactor=0.
328	DO l=1,nlayer-1
329	DO ig=1,ngrid
330	! Typical mixing ratio profile
331
332	zp=(pplev(ig,1)/pplay(ig,l))**(70./30) !emulating topdust
333	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
334
335	! Vertical scaling function
336	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor
337
338	ENDDO
339	ENDDO
340	tauh2so4tmp(1:ngrid)=0.
341	DO l=1,nlayer
342	DO ig=1,ngrid
343	tauh2so4tmp(ig) = tauh2so4tmp(ig) + aerosol(ig,l,iaer)
344	ENDDO
345	ENDDO
346	DO l=1,nlayer-1
347	DO ig=1,ngrid
348	aerosol(ig,l,iaer) = max(1E-20, &
349	1 &
350	* pplev(ig,1) / pplev(ig,1) &
351	* aerosol(ig,l,iaer) &
352	/ tauh2so4tmp(ig))
353
354	ENDDO
355	ENDDO
356
357	! 1/700. is assuming a "sulfurtau" of 1
358	! Sulfur aerosol routine to be improved.
359	! aerosol0 = &
360	! ( 0.75 * QREFvis3d(ig,l,iaer) / &
361	! ( rho_h2so4 * reffrad(ig,l,iaer) ) ) * &
362	! ( pq(ig,l,i_h2so4) + 1.E-9 ) * &
363	! ( pplev(ig,l) - pplev(ig,l+1) ) / g
364	! aerosol0 = max(aerosol0,1.e-9)
365	! aerosol0 = min(aerosol0,L_TAUMAX)
366	! aerosol(ig,l,iaer) = aerosol0
367
368	! ENDDO
369	! ENDDO
370	end if
371
372
373	! ---------------------------------------------------------
374	!==================================================================
375	! Two-layer aerosols (unknown composition)
376	! S. Guerlet (2013) - Modif by J. Vatant d'Ollone (2020)
377	!==================================================================
378
379	if (iaero_back2lay .ne.0) then
380	iaer=iaero_back2lay
381	! 1. Initialization
382	aerosol(1:ngrid,1:nlayer,iaer)=0.0
383	! 2. Opacity calculation
384
385
386	! JVO 20 : Modif to have each of the layers (strato and tropo) correctly normalized
387	! Otherwise we previously had the total optical depth correct but for each
388	! separately, so it didn't match the input values + what's more normalizing
389	! to the sum was making them non-independent : eg changing tau_tropo was
390	! affecting stratopsheric values of optical depth ...
391	!
392	! Note that the main consequence of the former version bug was (in most cases)
393	! to strongly underestimate the stratospheric optical depths compared to the
394	! required values, eg, with tau_tropo=10 and tau_strato=0.1, you actually ended
395	! with an actual tau_strato of 1E-4 ... !
396	!
397	! NB : Because of the extra transition opacity if the layers are non contiguous,
398	! be aware that at the the bottom we have tau > tau_strato + tau_tropo
399
400	DO ig=1,ngrid
401	dp_tropo(ig) = 0.D0
402	dp_strato(ig) = 0.D0
403	DO l=1,nlayer-1
404	aerosol(ig,l,iaer) = ( pplev(ig,l) - pplev(ig,l+1) )
405	!! 1. below tropospheric layer: no aerosols
406	IF (pplev(ig,l) .gt. pres_bottom_tropo) THEN
407	aerosol(ig,l,iaer) = 0.D0
408	!! 2. tropo layer
409	ELSEIF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN
410	dp_tropo(ig) = dp_tropo(ig) + aerosol(ig,l,iaer)
411	!! 3. linear transition
412	! JVO 20 : This interpolation needs to be done AFTER we set strato and tropo (see below)
413	!! 4. strato layer
414	ELSEIF (pplev(ig,l) .le. pres_bottom_strato .and. pplev(ig,l) .ge. pres_top_strato) THEN
415	dp_strato(ig) = dp_strato(ig) + aerosol(ig,l,iaer)
416	!! 5. above strato layer: no aerosols
417	ELSEIF (pplev(ig,l) .lt. pres_top_strato) THEN
418	aerosol(ig,l,iaer) = 0.D0
419	ENDIF
420	ENDDO
421	ENDDO
422
423	! 3. Re-normalize to the (input) observed (total) column (for each of the layers)
424
425	DO ig=1,ngrid
426	DO l=1,nlayer-1
427	IF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN
428	aerosol(ig,l,iaer) = obs_tau_col_tropo*aerosol(ig,l,iaer)/dp_tropo(ig)
429	ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN
430	expfactor=log(pplev(ig,l)/pres_top_tropo)/log(pres_bottom_strato/pres_top_tropo)
431	aerosol(ig,l,iaer) = (obs_tau_col_strato/dp_strato(ig))**expfactor &
432	* (obs_tau_col_tropo/dp_tropo(ig))**(1.0-expfactor) &
433	* aerosol(ig,l,iaer)
434	ELSEIF (pplev(ig,l) .le. pres_bottom_strato .and. pplev(ig,l) .ge. pres_top_strato) THEN
435	aerosol(ig,l,iaer) = obs_tau_col_strato*aerosol(ig,l,iaer)/dp_strato(ig)
436	ENDIF
437	write(,), aerosol(ig,l,iaer)
438	ENDDO
439	ENDDO
440
441
442	end if ! if Two-layer aerosols
443
444	!==================================================================
445	! Saturn/Jupiter ammonia cloud = thin cloud (scale height 0.2 hard coded...)
446	! S. Guerlet (2013)
447	!==================================================================
448
449	if (iaero_nh3 .ne.0) then
450	iaer=iaero_nh3
451	! 1. Initialization
452	aerosol(1:ngrid,1:nlayer,iaer)=0.D0
453	! 2. Opacity calculation
454	DO ig=1,ngrid
455
456	DO l=1,nlayer-1
457	!! 1. below cloud layer: no opacity
458
459	IF (pplev(ig,l) .gt. pres_nh3_cloud ) THEN
460	aerosol(ig,l,iaer) = 0.D0
461
462	ELSEIF (pplev(ig,l) .le. pres_nh3_cloud ) THEN
463	pente_cloud=5. !!(hard-coded, correspond to scale height 0.2)
464	aerosol(ig,l,iaer) = ((pplev(ig,l)/pres_nh3_cloud)*(pente_cloud))tau_nh3_cloud
465
466	ENDIF
467	ENDDO
468
469	END DO
470
471	! 3. Re-normalize to observed total column
472	tau_col(:)=0.0
473	DO l=1,nlayer
474	DO ig=1,ngrid
475	tau_col(ig) = tau_col(ig) &
476	+ aerosol(ig,l,iaer)/tau_nh3_cloud
477	ENDDO
478	ENDDO
479
480	DO ig=1,ngrid
481	DO l=1,nlayer-1
482	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig)
483	ENDDO
484	ENDDO
485
486	end if ! if NH3 cloud
487	!==================================================================
488	! Jovian auroral aerosols (unknown composition) NON-GENERIC: vertical and meridional profile tuned to observations
489	! S. Guerlet (2015)
490	!==================================================================
491
492
493	if (iaero_aurora .ne.0) then
494	iaer=iaero_aurora
495	! 1. Initialization
496	aerosol(1:ngrid,1:nlayer,iaer)=0.D0
497	pm = 2000. !!case study: maxi aerosols at 20 hPa
498	! 2. Opacity calculation
499	DO ig=1,ngrid
500
501	!! Test Jupiter (based on Zhang et al 2013 observations, but a bit different), decembre 2015
502	DO l=1,nlayer
503	aerosol(ig,l,iaer) = (pplev(ig,l)/pm)*2 exp(-(pplev(ig,l)/pm)**2)
504	ENDDO
505	ENDDO
506
507	! 3. Meridional distribution, and re-normalize to observed total column
508	tau_col(:)=0.D0
509	DO ig=1,ngrid
510	!!Jupiter
511	!!Hem sud:
512	IF (latitude(ig)180.D0/pi .lt. -45.D0 .and. latitude(ig)180.D0/pi .gt. -70.) THEN
513	obs_tau_col_aurora= 10.D0*(-0.06D0latitude(ig)*180.D0/pi-3.4D0)
514	ELSEIF (latitude(ig)180.D0/pi .lt. -37.D0 .and. latitude(ig)180.D0/pi .ge. -45.) THEN
515	obs_tau_col_aurora= 10.D0*(-0.3D0latitude(ig)*180.D0/pi-14.3D0)
516	ELSEIF (latitude(ig)*180./pi .le. -70. ) THEN
517	obs_tau_col_aurora= 10*(0.0670.-3.4D0)
518	!!Hem Nord:
519	ELSEIF (latitude(ig)180.D0/pi .gt. 30.D0 .and. latitude(ig)180.D0/pi .lt. 70.) THEN
520	obs_tau_col_aurora= 10.D0*(0.03D0latitude(ig)*180.D0/pi-1.17D0)
521	ELSEIF (latitude(ig)180.D0/pi .gt. 22.D0 .and. latitude(ig)180.D0/pi .le. 30.) THEN
522	obs_tau_col_aurora= 10.D0*(0.3D0latitude(ig)*180.D0/pi-9.4D0)
523	ELSEIF (latitude(ig)*180.D0/pi .ge. 70.) THEN
524	obs_tau_col_aurora= 10*(0.0370.-1.17D0)
525	ELSEIF (latitude(ig)180.D0/pi .ge. -37. .and. latitude(ig)180.D0/pi .le. 22.) THEN
526	obs_tau_col_aurora = 0.001D0 !!Jupiter: mini pas a zero
527	ENDIF
528
529	DO l=1,nlayer
530	tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)/obs_tau_col_aurora
531	ENDDO
532	ENDDO
533
534	DO ig=1,ngrid
535	DO l=1,nlayer-1
536	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig)
537	ENDDO
538	ENDDO
539
540
541	end if ! if Auroral aerosols
542
543
544	! --------------------------------------------------------------------------
545	! Column integrated visible optical depth in each point (used for diagnostic)
546
547	tau_col(:)=0.0
548	do iaer = 1, naerkind
549	do l=1,nlayer
550	do ig=1,ngrid
551	tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
552	end do
553	end do
554	end do
555
556	do ig=1,ngrid
557	do l=1,nlayer
558	do iaer = 1, naerkind
559	if(aerosol(ig,l,iaer).gt.1.e3)then
560	print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
561	print*,'at ig=',ig,', l=',l,', iaer=',iaer
562	print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
563	print*,'reffrad=',reffrad(ig,l,iaer)
564	endif
565	end do
566	end do
567	end do
568
569	do ig=1,ngrid
570	if(tau_col(ig).gt.1.e3)then
571	print*,'WARNING: tau_col=',tau_col(ig)
572	print*,'at ig=',ig
573	print*,'aerosol=',aerosol(ig,:,:)
574	print*,'QREFvis3d=',QREFvis3d(ig,:,:)
575	print*,'reffrad=',reffrad(ig,:,:)
576	endif
577	end do
578	return
579	end subroutine aeropacity
580

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

Download in other formats:

Original Format