Context Navigation

aeropacity.F90 @ 2613

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

Add a generic n-layer aerosol scheme to replace the former buggy 2-layer scheme as well as the hard-coded NH3 cloud.

It can be called using 'aeronlay=.true.' in callphys.def, and set the number of layers (up to 4) with 'nlayaero'.
Then, the following parameters are read as arrays of size nlayaero in callphys.def (separated by blank space)

*aeronlay_tauref (Optical depth of aerosol layer at ref wavelenght)
*aeronlay_lamref (Ref wavelenght (m))
*aeronlay_choice (Choice for vertical profile - 1:tau follows atm scale height btwn top and bottom - 2:tau follows it own scale height)
*aeronlay_pbot (Bottom pressure (Pa))
*aeronlay_ptop (Top pressure (Pa) - useful only if choice=1)
*aeronlay_sclhght (Ratio of aerosol layer scale height / atmospheric scale height - useful only if choice=2 )
*aeronlay_size (Particle size (m))
*optprop_aeronlay_vis File for VIS opt properties.
*optprop_aeronlay_ir File for IR opt properties.

+Extra info :

+ In addition of solving the bug from 2-layer it enables different optical properties.
+ The former scheme are left for retrocompatibility (for now) but you should use the new one.
+ See aeropacity.F90 for the calculations

+ Each layer can have different optical properties, size of particle ...
+ You have different choices for vertical profile of the aerosol layers :

aeronlay_choice = 1 : Layer tau is spread between ptop and pbot following atm scale height.
aeronlay_choice = 2 : Layer tau follows its own scale height above cloud deck (pbot).

In this case ptop is dummy and sclhght gives the ratio H_cl/H_atm.

+ The reference wavelenght for input optical depth is now read as input (aeronlay_lamref)
+ Following the last point some comment is added in suaer_corrk about the 'not-really-dummy'ness of IR lamref..

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

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

Download in other formats:

Original Format