Context Navigation

aeropacity.F90 @ 2889

Last change on this file since 2889 was 2831, checked in by emillour, 3 years ago
Generic PCM: Add the possibility to include Venus-like aerosols (triggered by option aerovenus=.true. in callphys.def); baseline is to use 5 distinct scatterers but each may be turned on/off (via aerovenus1, aerovenus2, aerovenus2p, aerovenus3, aerovenusUV flags which may be specified in callphys.def). GG
File size: 37.8 KB

Line
1	Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pt, pq, &
2	aerosol,reffrad,nueffrad, QREFvis3d,QREFir3d,tau_col, &
3	cloudfrac,totcloudfrac,clearsky)
4
5	use radinc_h, only : L_TAUMAX,naerkind
6	use aerosol_mod, only: iaero_nlay, iaero_generic, &
7	iaero_aurora, iaero_back2lay, iaero_co2, &
8	iaero_dust, iaero_h2o, iaero_h2so4, &
9	iaero_nh3, i_rgcs_ice, noaero, &
10	iaero_venus1, iaero_venus2, iaero_venus2p, &
11	iaero_venus3, iaero_venusUV
12	USE tracer_h, only: noms,rho_co2,rho_ice,rho_q,mmol
13	use comcstfi_mod, only: g, pi, mugaz, avocado
14	use geometry_mod, only: latitude
15	use callkeys_mod, only: aerofixco2,aerofixh2o,kastprof,cloudlvl, &
16	CLFvarying,CLFfixval,dusttau, &
17	pres_bottom_tropo,pres_top_tropo,obs_tau_col_tropo, &
18	pres_bottom_strato,pres_top_strato,obs_tau_col_strato, &
19	tau_nh3_cloud, pres_nh3_cloud, &
20	nlayaero, aeronlay_tauref, aeronlay_choice, &
21	aeronlay_pbot, aeronlay_ptop, aeronlay_sclhght, &
22	aerogeneric
23	use generic_tracer_index_mod, only: generic_tracer_index
24	implicit none
25
26	!==================================================================
27	!
28	! Purpose
29	! -------
30	! Compute aerosol optical depth in each gridbox.
31	!
32	! Authors
33	! -------
34	! F. Forget
35	! F. Montmessin (water ice scheme)
36	! update J.-B. Madeleine (2008)
37	! dust removal, simplification by Robin Wordsworth (2009)
38	! Generic n-layer aerosol - J. Vatant d'Ollone (2020)
39	! Radiative Generic Condensable Species - Lucas Teinturier (2022)
40	!
41	! Input
42	! -----
43	! ngrid Number of horizontal gridpoints
44	! nlayer Number of layers
45	! nq Number of tracers
46	! pplev Pressure (Pa) at each layer boundary
47	! pq Aerosol mixing ratio
48	! reffrad(ngrid,nlayer,naerkind) Aerosol effective radius
49	! QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
50	! QREFir3d(ngrid,nlayer,naerkind) / at reference wavelengths
51	!
52	! Output
53	! ------
54	! aerosol Aerosol optical depth in layer l, grid point ig
55	! tau_col Total column optical depth at grid point ig
56	!
57	!=======================================================================
58
59	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
60	INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
61	INTEGER,INTENT(IN) :: nq ! number of tracers
62	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
63	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
64	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
65	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! mid-layer temperature (K)
66	REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth
67	REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) ! aerosol effective radius
68	REAL,INTENT(IN) :: nueffrad(ngrid,nlayer,naerkind) ! aerosol effective variance
69	REAL,INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind) ! extinction coefficient in the visible
70	REAL,INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
71	REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth
72	! BENJAMIN MODIFS
73	real,intent(in) :: cloudfrac(ngrid,nlayer) ! cloud fraction
74	real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
75	logical,intent(in) :: clearsky
76
77	real aerosol0, obs_tau_col_aurora, pm
78	real pcloud_deck, cloud_slope
79
80	real dp_strato(ngrid)
81	real dp_tropo(ngrid)
82	real dp_layer(ngrid)
83
84	INTEGER l,ig,iq,iaer,ia
85
86	LOGICAL,SAVE :: firstcall=.true.
87	!$OMP THREADPRIVATE(firstcall)
88	REAL CBRT
89	EXTERNAL CBRT
90
91	INTEGER,SAVE :: i_co2ice=0 ! co2 ice
92	INTEGER,SAVE :: i_h2oice=0 ! water ice
93	!$OMP THREADPRIVATE(i_co2ice,i_h2oice)
94	CHARACTER(LEN=20) :: tracername ! to temporarily store text
95
96	! for fixed dust profiles
97	real topdust, expfactor, zp
98	REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
99	REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling
100
101	real CLFtot
102	integer igen_ice,igen_vap ! to store the index of generic tracer
103	logical dummy_bool ! dummy boolean just in case we need one
104	! integer i_rgcs_ice(aerogeneric)
105	! for venus clouds
106	real :: p_bot,p_top,h_bot,h_top,mode_dens,h_lay
107
108	! identify tracers
109	IF (firstcall) THEN
110	ia =0
111	write(,) "Tracers found in aeropacity:"
112	do iq=1,nq
113	tracername=noms(iq)
114	if (tracername.eq."co2_ice") then
115	i_co2ice=iq
116	write(,) "i_co2ice=",i_co2ice
117
118	endif
119	if (tracername.eq."h2o_ice") then
120	i_h2oice=iq
121	write(,) "i_h2oice=",i_h2oice
122	endif
123	enddo
124
125	if (noaero) then
126	print*, "No active aerosols found in aeropacity"
127	else
128	print*, "If you would like to use aerosols, make sure any old"
129	print*, "start files are updated in newstart using the option"
130	print*, "q=0"
131	write(,) "Active aerosols found in aeropacity:"
132	endif
133
134	if ((iaero_co2.ne.0).and.(.not.noaero)) then
135	print*, 'iaero_co2= ',iaero_co2
136	endif
137	if (iaero_h2o.ne.0) then
138	print*,'iaero_h2o= ',iaero_h2o
139	endif
140	if (iaero_dust.ne.0) then
141	print*,'iaero_dust= ',iaero_dust
142	endif
143	if (iaero_h2so4.ne.0) then
144	print*,'iaero_h2so4= ',iaero_h2so4
145	endif
146	if (iaero_back2lay.ne.0) then
147	print*,'iaero_back2lay= ',iaero_back2lay
148	endif
149	if (iaero_nh3.ne.0) then
150	print*,'iaero_nh3= ',iaero_nh3
151	endif
152	if (iaero_nlay(1).ne.0) then
153	print*,'iaero_nlay= ',iaero_nlay(:)
154	endif
155	if (iaero_aurora.ne.0) then
156	print*,'iaero_aurora= ',iaero_aurora
157	endif
158
159	if (iaero_venus1.ne.0) then
160	print*,'iaero_venus1= ',iaero_venus1
161	endif
162	if (iaero_venus2.ne.0) then
163	print*,'iaero_venus2= ',iaero_venus2
164	endif
165	if (iaero_venus2p.ne.0) then
166	print*,'iaero_venus2p= ',iaero_venus2p
167	endif
168	if (iaero_venus3.ne.0) then
169	print*,'iaero_venus3= ',iaero_venus3
170	endif
171	if (iaero_venusUV.ne.0) then
172	print*,'iaero_venusUV= ',iaero_venusUV
173	endif
174
175	if (aerogeneric .ne. 0) then
176	print*,"iaero_generic= ",iaero_generic(:)
177	endif
178	firstcall=.false.
179	ENDIF ! of IF (firstcall)
180
181
182	! ---------------------------------------------------------
183	!==================================================================
184	! CO2 ice aerosols
185	!==================================================================
186
187	if (iaero_co2.ne.0) then
188	iaer=iaero_co2
189	! 1. Initialization
190	aerosol(1:ngrid,1:nlayer,iaer)=0.0
191	! 2. Opacity calculation
192	if (noaero) then ! aerosol set to zero
193	aerosol(1:ngrid,1:nlayer,iaer)=0.0
194	elseif (aerofixco2.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed
195	aerosol(1:ngrid,1:nlayer,iaer)=1.e-9
196	!aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option
197	else
198	DO ig=1, ngrid
199	DO l=1,nlayer-1 ! to stop the rad tran bug
200
201	aerosol0 = &
202	( 0.75 * QREFvis3d(ig,l,iaer) / &
203	( rho_co2 * reffrad(ig,l,iaer) ) ) * &
204	( pq(ig,l,i_co2ice) + 1.E-9 ) * &
205	( pplev(ig,l) - pplev(ig,l+1) ) / g
206	aerosol0 = max(aerosol0,1.e-9)
207	aerosol0 = min(aerosol0,L_TAUMAX)
208	aerosol(ig,l,iaer) = aerosol0
209	! aerosol(ig,l,iaer) = 0.0
210	! print*, aerosol(ig,l,iaer)
211	! using cloud fraction
212	! aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF))
213	! aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX)
214
215
216	ENDDO
217	ENDDO
218	end if ! if fixed or varying
219	end if ! if CO2 aerosols
220	!==================================================================
221	! Water ice / liquid
222	!==================================================================
223
224	if (iaero_h2o.ne.0) then
225	iaer=iaero_h2o
226	! 1. Initialization
227	aerosol(1:ngrid,1:nlayer,iaer)=0.0
228	! 2. Opacity calculation
229	if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
230	aerosol(1:ngrid,1:nlayer,iaer) =1.e-9
231
232	! put cloud at cloudlvl
233	if(kastprof.and.(cloudlvl.ne.0.0))then
234	ig=1
235	do l=1,nlayer
236	if(int(cloudlvl).eq.l)then
237	!if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
238	print*,'Inserting cloud at level ',l
239	!aerosol(ig,l,iaer)=10.0
240
241	rho_ice=920.0
242
243	! the Kasting approximation
244	aerosol(ig,l,iaer) = &
245	( 0.75 * QREFvis3d(ig,l,iaer) / &
246	( rho_ice * reffrad(ig,l,iaer) ) ) * &
247	!( pq(ig,l,i_h2oice) + 1.E-9 ) * &
248	( 4.0e-4 + 1.E-9 ) * &
249	( pplev(ig,l) - pplev(ig,l+1) ) / g
250
251
252	open(115,file='clouds.out',form='formatted')
253	write(115,*) l,aerosol(ig,l,iaer)
254	close(115)
255
256	return
257	endif
258	end do
259
260	call abort
261	endif
262
263	else
264
265	do ig=1, ngrid
266	!do l=1,nlayer-1 ! to stop the rad tran bug
267	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)
268	! same correction should b-probably be done for other aerosol types.
269	aerosol(ig,l,iaer) = & !modification by BC
270	( 0.75 * QREFvis3d(ig,l,iaer) / &
271	( rho_ice * reffrad(ig,l,iaer) ) ) * &
272	! pq(ig,l,i_h2oice) * & !JL I dropped the +1e-9 here to have the same
273	!( pplev(ig,l) - pplev(ig,l+1) ) / g ! opacity in the clearsky=true and the
274	! clear=false/pq=0 case
275	( pq(ig,l,i_h2oice) + 1.E-9 ) * & ! Doing this makes the code unstable, so I have restored it (RW)
276	( pplev(ig,l) - pplev(ig,l+1) ) / g
277
278	enddo
279	enddo
280
281	if(CLFvarying)then
282	call totalcloudfrac(ngrid,nlayer,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1,1,iaer))
283	do ig=1, ngrid
284	!do l=1,nlayer-1 ! to stop the rad tran bug
285	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)
286	CLFtot = max(totcloudfrac(ig),0.01)
287	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
288	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
289	enddo
290	enddo
291	else
292	do ig=1, ngrid
293	!do l=1,nlayer-1 ! to stop the rad tran bug
294	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)
295	CLFtot = CLFfixval
296	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
297	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
298	enddo
299	enddo
300	end if!(CLFvarying)
301	endif !(aerofixed.or.(i_h2oice.eq.0).or.clearsky)
302
303	end if ! End if h2o aerosol
304
305	!==================================================================
306	! Dust
307	!==================================================================
308	if (iaero_dust.ne.0) then
309	iaer=iaero_dust
310	! 1. Initialization
311	aerosol(1:ngrid,1:nlayer,iaer)=0.0
312
313	topdust=30.0 ! km (used to be 10.0 km) LK
314
315	! 2. Opacity calculation
316
317	! expfactor=0.
318	DO l=1,nlayer-1
319	DO ig=1,ngrid
320	! Typical mixing ratio profile
321
322	zp=(pplev(ig,1)/pplay(ig,l))**(70./topdust)
323	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
324
325	! Vertical scaling function
326	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
327	*expfactor
328
329
330	ENDDO
331	ENDDO
332
333	! Rescaling each layer to reproduce the choosen (or assimilated)
334	! dust extinction opacity at visible reference wavelength, which
335	! is scaled to the surface pressure pplev(ig,1)
336
337	taudusttmp(1:ngrid)=0.
338	DO l=1,nlayer
339	DO ig=1,ngrid
340	taudusttmp(ig) = taudusttmp(ig) &
341	+ aerosol(ig,l,iaer)
342	ENDDO
343	ENDDO
344	DO l=1,nlayer-1
345	DO ig=1,ngrid
346	aerosol(ig,l,iaer) = max(1E-20, &
347	dusttau &
348	* pplev(ig,1) / pplev(ig,1) &
349	* aerosol(ig,l,iaer) &
350	/ taudusttmp(ig))
351
352	ENDDO
353	ENDDO
354	end if ! If dust aerosol
355
356	!==================================================================
357	! H2SO4
358	!==================================================================
359	! added by LK
360	if (iaero_h2so4.ne.0) then
361	iaer=iaero_h2so4
362
363	! 1. Initialization
364	aerosol(1:ngrid,1:nlayer,iaer)=0.0
365
366
367	! 2. Opacity calculation
368
369	! expfactor=0.
370	DO l=1,nlayer-1
371	DO ig=1,ngrid
372	! Typical mixing ratio profile
373
374	zp=(pplev(ig,1)/pplay(ig,l))**(70./30) !emulating topdust
375	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
376
377	! Vertical scaling function
378	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor
379
380	ENDDO
381	ENDDO
382	tauh2so4tmp(1:ngrid)=0.
383	DO l=1,nlayer
384	DO ig=1,ngrid
385	tauh2so4tmp(ig) = tauh2so4tmp(ig) + aerosol(ig,l,iaer)
386	ENDDO
387	ENDDO
388	DO l=1,nlayer-1
389	DO ig=1,ngrid
390	aerosol(ig,l,iaer) = max(1E-20, &
391	1 &
392	* pplev(ig,1) / pplev(ig,1) &
393	* aerosol(ig,l,iaer) &
394	/ tauh2so4tmp(ig))
395
396	ENDDO
397	ENDDO
398
399	! 1/700. is assuming a "sulfurtau" of 1
400	! Sulfur aerosol routine to be improved.
401	! aerosol0 = &
402	! ( 0.75 * QREFvis3d(ig,l,iaer) / &
403	! ( rho_h2so4 * reffrad(ig,l,iaer) ) ) * &
404	! ( pq(ig,l,i_h2so4) + 1.E-9 ) * &
405	! ( pplev(ig,l) - pplev(ig,l+1) ) / g
406	! aerosol0 = max(aerosol0,1.e-9)
407	! aerosol0 = min(aerosol0,L_TAUMAX)
408	! aerosol(ig,l,iaer) = aerosol0
409
410	! ENDDO
411	! ENDDO
412	end if
413
414
415	! ---------------------------------------------------------
416	!==================================================================
417	! Two-layer aerosols (unknown composition)
418	! S. Guerlet (2013) - Modif by J. Vatant d'Ollone (2020)
419	!
420	! This scheme is deprecated and left for retrocompatibility
421	! You should use the n-layer scheme below !
422	!
423	!==================================================================
424
425	if (iaero_back2lay .ne.0) then
426	iaer=iaero_back2lay
427	! 1. Initialization
428	aerosol(1:ngrid,1:nlayer,iaer)=0.0
429	! 2. Opacity calculation
430
431
432	! JVO 20 : Modif to have each of the layers (strato and tropo) correctly normalized
433	! Otherwise we previously had the total optical depth correct but for each
434	! separately, so it didn't match the input values + what's more normalizing
435	! to the sum was making them non-independent : eg changing tau_tropo was
436	! affecting stratopsheric values of optical depth ...
437	!
438	! Note that the main consequence of the former version bug was (in most cases)
439	! to strongly underestimate the stratospheric optical depths compared to the
440	! required values, eg, with tau_tropo=10 and tau_strato=0.1, you actually ended
441	! with an actual tau_strato of 1E-4 ... !
442	!
443	! NB : Because of the extra transition opacity if the layers are non contiguous,
444	! be aware that at the the bottom we have tau > tau_strato + tau_tropo
445
446	DO ig=1,ngrid
447	dp_tropo(ig) = 0.D0
448	dp_strato(ig) = 0.D0
449	DO l=1,nlayer-1
450	aerosol(ig,l,iaer) = ( pplev(ig,l) - pplev(ig,l+1) )
451	!! 1. below tropospheric layer: no aerosols
452	IF (pplev(ig,l) .gt. pres_bottom_tropo) THEN
453	aerosol(ig,l,iaer) = 0.D0
454	!! 2. tropo layer
455	ELSEIF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN
456	dp_tropo(ig) = dp_tropo(ig) + aerosol(ig,l,iaer)
457	!! 3. linear transition
458	! JVO 20 : This interpolation needs to be done AFTER we set strato and tropo (see below)
459	!! 4. strato layer
460	ELSEIF (pplev(ig,l) .le. pres_bottom_strato .and. pplev(ig,l) .ge. pres_top_strato) THEN
461	dp_strato(ig) = dp_strato(ig) + aerosol(ig,l,iaer)
462	!! 5. above strato layer: no aerosols
463	ELSEIF (pplev(ig,l) .lt. pres_top_strato) THEN
464	aerosol(ig,l,iaer) = 0.D0
465	ENDIF
466	ENDDO
467	ENDDO
468
469	! 3. Re-normalize to the (input) observed (total) column (for each of the layers)
470
471	DO ig=1,ngrid
472	DO l=1,nlayer-1
473	IF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN
474	aerosol(ig,l,iaer) = obs_tau_col_tropo*aerosol(ig,l,iaer)/dp_tropo(ig)
475	ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN
476	expfactor=log(pplev(ig,l)/pres_top_tropo)/log(pres_bottom_strato/pres_top_tropo)
477	aerosol(ig,l,iaer) = (obs_tau_col_strato/dp_strato(ig))**expfactor &
478	* (obs_tau_col_tropo/dp_tropo(ig))**(1.0-expfactor) &
479	* aerosol(ig,l,iaer)
480	ELSEIF (pplev(ig,l) .le. pres_bottom_strato .and. pplev(ig,l) .ge. pres_top_strato) THEN
481	aerosol(ig,l,iaer) = obs_tau_col_strato*aerosol(ig,l,iaer)/dp_strato(ig)
482	ENDIF
483	ENDDO
484	ENDDO
485
486
487	end if ! if Two-layer aerosols
488
489	!==================================================================
490	! Saturn/Jupiter ammonia cloud = thin cloud (scale height 0.2 hard coded...)
491	! S. Guerlet (2013)
492	! JVO 20 : You should now use the generic n-layer scheme below
493	!==================================================================
494
495	if (iaero_nh3 .ne.0) then
496	iaer=iaero_nh3
497	! 1. Initialization
498	aerosol(1:ngrid,1:nlayer,iaer)=0.D0
499	! 2. Opacity calculation
500	DO ig=1,ngrid
501
502	DO l=1,nlayer-1
503	!! 1. below cloud layer: no opacity
504
505	IF (pplev(ig,l) .gt. pres_nh3_cloud ) THEN
506	aerosol(ig,l,iaer) = 0.D0
507
508	ELSEIF (pplev(ig,l) .le. pres_nh3_cloud ) THEN
509	cloud_slope=5. !!(hard-coded, correspond to scale height 0.2)
510	aerosol(ig,l,iaer) = ((pplev(ig,l)/pres_nh3_cloud)*(cloud_slope))tau_nh3_cloud
511
512	ENDIF
513	ENDDO
514
515	END DO
516
517	! 3. Re-normalize to observed total column
518	dp_layer(:)=0.0
519	DO l=1,nlayer
520	DO ig=1,ngrid
521	dp_layer(ig) = dp_layer(ig) &
522	+ aerosol(ig,l,iaer)/tau_nh3_cloud
523	ENDDO
524	ENDDO
525
526	DO ig=1,ngrid
527	DO l=1,nlayer-1
528	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/dp_layer(ig)
529	ENDDO
530	ENDDO
531
532	end if ! if NH3 cloud
533
534	!=========================================================================================================
535	! Generic N-layers aerosols/clouds
536	! Author : J. Vatant d'Ollone (2020)
537	!
538	! Purpose: Replaces and extents the former buggy 2-layer scheme as well as hard-coded NH3 cloud
539	!
540	! + Each layer can have different optical properties, size of particle ...
541	! + Enables up to n=4 layers as we apparently cannot run with more scatterers (could be worth checking...)
542	! + You have different choices for vertical profile of the aerosol layers :
543	! * aeronlay_choice = 1 : Layer tau is spread between ptop and pbot following atm scale height.
544	! * aeronlay_choice = 2 : Layer tau follows its own scale height above cloud deck (pbot).
545	! In this case ptop is dummy and sclhght gives the ratio H_cl/H_atm.
546	! * aeronlay_choice = ... feel free to add more cases !
547	! + Layers can overlap if needed (if you want a 'transition layer' as in the 2-scheme, just add it)
548	!
549	!=========================================================================================================
550
551	if (iaero_nlay(1) .ne.0) then
552
553	DO ia=1,nlayaero
554	iaer=iaero_nlay(ia)
555
556	! a. Initialization
557	aerosol(1:ngrid,1:nlayer,iaer)=0.D0
558
559	! b. Opacity calculation
560
561	! Case 1 : Follows atmospheric scale height between boundaries pressures
562	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
563	IF (aeronlay_choice(ia).eq.1) THEN
564
565	dp_layer(:)=0.D0
566	DO ig=1,ngrid
567	DO l=1,nlayer-1
568	!! i. Opacity follows scale height
569	IF ( pplev(ig,l).le.aeronlay_pbot(ia) .AND. &
570	pplev(ig,l).ge.aeronlay_ptop(ia) ) THEN
571	aerosol(ig,l,iaer) = ( pplev(ig,l) - pplev(ig,l+1) )
572	dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)
573	!! ii. Outside aerosol layer boundaries: no aerosols
574	ELSE
575	aerosol(ig,l,iaer) = 0.D0
576	ENDIF
577	ENDDO
578	ENDDO
579	! iii. Re-normalize to required total opacity
580	DO ig=1,ngrid
581	DO l=1,nlayer-1
582	IF ( pplev(ig,l).le.aeronlay_pbot(ia) .AND. &
583	pplev(ig,l).ge.aeronlay_ptop(ia) ) THEN
584	aerosol(ig,l,iaer) = aerosol(ig,l,iaer) / dp_layer(ig) &
585	* aeronlay_tauref(ia)
586	ENDIF
587	ENDDO
588	ENDDO
589
590	! Case 2 : Follows input scale height
591	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
592	ELSE IF (aeronlay_choice(ia).eq.2) THEN
593
594	cloud_slope = 1.D0/aeronlay_sclhght(ia)
595	pcloud_deck = aeronlay_pbot(ia)
596	dp_layer(:) = 0.D0
597
598	DO ig=1,ngrid
599	DO l=1,nlayer-1
600	!! i. Below cloud layer: no opacity
601	IF (pplev(ig,l) .gt. pcloud_deck) THEN
602	aerosol(ig,l,iaer) = 0.D0
603	!! ii. Follows scale height above cloud deck
604	ELSEIF (pplev(ig,l) .le. pcloud_deck) THEN
605	aerosol(ig,l,iaer) = ((pplev(ig,l)/pcloud_deck)**(cloud_slope))
606	dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)
607	ENDIF
608	ENDDO
609	ENDDO
610	! iii. Re-normalize to required total opacity
611	DO ig=1,ngrid
612	DO l=1,nlayer-1
613	IF (pplev(ig,l) .le. pcloud_deck) THEN
614	aerosol(ig,l,iaer) = aerosol(ig,l,iaer) / dp_layer(ig) &
615	* aeronlay_tauref(ia)
616	ENDIF
617	ENDDO
618	ENDDO
619
620	ENDIF ! aeronlay_choice
621
622	ENDDO ! loop on n aerosol layers
623
624	end if ! if N-layer aerosols
625
626	!==================================================================
627	! Jovian auroral aerosols (unknown composition) NON-GENERIC: vertical and meridional profile tuned to observations
628	! S. Guerlet (2015)
629	!==================================================================
630
631
632	if (iaero_aurora .ne.0) then
633	iaer=iaero_aurora
634	! 1. Initialization
635	aerosol(1:ngrid,1:nlayer,iaer)=0.D0
636	pm = 2000. !!case study: maxi aerosols at 20 hPa
637	! 2. Opacity calculation
638	DO ig=1,ngrid
639
640	!! Test Jupiter (based on Zhang et al 2013 observations, but a bit different), decembre 2015
641	DO l=1,nlayer
642	aerosol(ig,l,iaer) = (pplev(ig,l)/pm)*2 exp(-(pplev(ig,l)/pm)**2)
643	ENDDO
644	ENDDO
645
646	! 3. Meridional distribution, and re-normalize to observed total column
647	dp_layer(:)=0.D0
648	DO ig=1,ngrid
649	!!Jupiter
650	!!Hem sud:
651	IF (latitude(ig)180.D0/pi .lt. -45.D0 .and. latitude(ig)180.D0/pi .gt. -70.) THEN
652	obs_tau_col_aurora= 10.D0*(-0.06D0latitude(ig)*180.D0/pi-3.4D0)
653	ELSEIF (latitude(ig)180.D0/pi .lt. -37.D0 .and. latitude(ig)180.D0/pi .ge. -45.) THEN
654	obs_tau_col_aurora= 10.D0*(-0.3D0latitude(ig)*180.D0/pi-14.3D0)
655	ELSEIF (latitude(ig)*180./pi .le. -70. ) THEN
656	obs_tau_col_aurora= 10*(0.0670.-3.4D0)
657	!!Hem Nord:
658	ELSEIF (latitude(ig)180.D0/pi .gt. 30.D0 .and. latitude(ig)180.D0/pi .lt. 70.) THEN
659	obs_tau_col_aurora= 10.D0*(0.03D0latitude(ig)*180.D0/pi-1.17D0)
660	ELSEIF (latitude(ig)180.D0/pi .gt. 22.D0 .and. latitude(ig)180.D0/pi .le. 30.) THEN
661	obs_tau_col_aurora= 10.D0*(0.3D0latitude(ig)*180.D0/pi-9.4D0)
662	ELSEIF (latitude(ig)*180.D0/pi .ge. 70.) THEN
663	obs_tau_col_aurora= 10*(0.0370.-1.17D0)
664	ELSEIF (latitude(ig)180.D0/pi .ge. -37. .and. latitude(ig)180.D0/pi .le. 22.) THEN
665	obs_tau_col_aurora = 0.001D0 !!Jupiter: mini pas a zero
666	ENDIF
667
668	DO l=1,nlayer
669	dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)/obs_tau_col_aurora
670	ENDDO
671	ENDDO
672
673	DO ig=1,ngrid
674	DO l=1,nlayer-1
675	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/dp_layer(ig)
676	ENDDO
677	ENDDO
678
679
680	end if ! if Auroral aerosols
681	!===========================================================================
682	! Radiative Generic Condensable aerosols scheme
683	! Only used when we give aerogeneric != 0 in callphys.def
684	! Computes the generic aerosols' opacity in the same fashion as water of
685	! dust, using the QREFvis3d of the concerned specie
686	! Lucas Teinturier (2022)
687	!===========================================================================
688	if (aerogeneric .ne. 0) then ! we enter the scheme
689	do ia=1,aerogeneric
690	iaer = iaero_generic(ia)
691	! Initialization
692	aerosol(1:ngrid,1:nlayer,iaer) = 0.D0
693	igen_ice = i_rgcs_ice(ia)
694	! Let's loop on the horizontal and vertical grid
695	do ig=1,ngrid
696	do l=1,nlayer
697	aerosol(ig,l,iaer) = ( 0.75*QREFvis3d(ig,l,iaer) / &
698	(rho_q(igen_ice) * reffrad(ig,l,iaer)) ) * &
699	(pq(ig,l,igen_ice)+1E-9 ) * &
700	(pplev(ig,l) - pplev(ig,l+1)) /g
701	enddo !l=1,nlayer
702	enddo !ig=1,ngrid
703	enddo !ia=1,aerogeneric
704	endif !aerogeneric .ne. 0
705
706	!==================================================================
707	! Venus clouds (4 modes)
708	! S. Lebonnois (jan 2016)
709	!==================================================================
710	! distributions from Haus et al, 2013
711	! mode 1 2 2p 3
712	! r (microns) 0.30 1.05 1.40 3.65
713	! sigma 1.56 1.29 1.23 1.28
714	! reff (microns) 0.49 1.23 1.56 4.25
715	! nueff 0.21 0.067 0.044 0.063
716	! (nueff=exp(ln^2 sigma)-1)
717	!
718	! p_bot <=> zb ; p_top <=> zb+zc ; h_bot <=> Hlo ; h_top <=> Hup
719	! p<p_top: N=No(p/p_top)*(h_lay/h_top) h_lay=RT/g (in m)
720	! p>p_bot: N=No(p_bot/p)*(h_lay/h_bot) R=8.31/mu (mu in kg/mol)
721	! N is in m-3
722	!
723	! dTau = Qext[pireff*2exp(-3ln(1+nueff))]Nh_lay(-dp)/p
724
725	! Mode 1
726	if (iaero_venus1 .ne.0) then
727	iaer=iaero_venus1
728
729	! 1. Initialization
730	aerosol(1:ngrid,1:nlayer,iaer)=0.0
731	p_bot = 1.e5
732	p_top = 1.e4
733	h_bot = 1.0e3 ! m
734	h_top = 5.0e3
735
736	! 2. Opacity calculation
737
738	DO ig=1,ngrid
739	DO l=1,nlayer-1
740
741	h_lay=8.31pt(ig,l)/(g0.044)
742
743	!! 1. below 2e5 Pa: no aerosols
744	IF (pplay(ig,l) .gt. 2.e5) THEN
745	mode_dens = 0.
746
747	!! 2. cloud layer
748	ELSEIF (pplay(ig,l) .le. 2.e5 .and. pplay(ig,l) .gt. p_bot) THEN
749	mode_dens = 1.81e8(p_bot/pplay(ig,l))*(h_lay/h_bot)
750
751	ELSEIF (pplay(ig,l) .le. p_bot .and. pplay(ig,l) .gt. p_top) THEN
752	mode_dens = 1.81e8 ! m-3
753
754	ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. 1.e2) THEN
755	mode_dens = 1.81e8(pplay(ig,l)/p_top)*(h_lay/h_top)
756
757	!! 3. above 1.e2 Pa: no aerosols
758	ELSEIF (pplay(ig,l) .le. 1.e2) THEN
759	mode_dens = 0.
760	ENDIF
761
762	aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)* &
763	pi(reffrad(ig,l,iaer))2exp(-3.log(1+nueffrad(ig,l,iaer))) &
764	mode_densh_lay(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
765
766	ENDDO
767	ENDDO
768
769	end if ! mode 1
770
771	! Mode 2
772	if (iaero_venus2 .ne.0) then
773	iaer=iaero_venus2
774
775	! 1. Initialization
776	aerosol(1:ngrid,1:nlayer,iaer)=0.0
777	p_bot = 1.1e4
778	p_top = 1.e4
779	h_bot = 3.0e3
780	h_top = 3.5e3
781
782	! 2. Opacity calculation
783
784	DO ig=1,ngrid
785	DO l=1,nlayer-1
786
787	h_lay=8.31pt(ig,l)/(g0.044)
788
789	!! 1. below 2e5 Pa: no aerosols
790	IF (pplay(ig,l) .gt. 2.e5) THEN
791	mode_dens = 0.
792
793	!! 2. cloud layer
794	ELSEIF (pplay(ig,l) .le. 2.e5 .and. pplay(ig,l) .gt. p_bot) THEN
795	mode_dens = 1.00e8(p_bot/pplay(ig,l))*(h_lay/h_bot)
796
797	ELSEIF (pplay(ig,l) .le. p_bot .and. pplay(ig,l) .gt. p_top) THEN
798	mode_dens = 1.00e8
799
800	ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. 1.e2) THEN
801	mode_dens = 1.00e8(pplay(ig,l)/p_top)*(h_lay/h_top)
802
803	!! 3. above 1.e2 Pa: no aerosols
804	ELSEIF (pplay(ig,l) .le. 1.e2) THEN
805	mode_dens = 0.
806	ENDIF
807
808	aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)* &
809	pi(reffrad(ig,l,iaer))2exp(-3.log(1+nueffrad(ig,l,iaer))) &
810	mode_densh_lay(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
811
812	ENDDO
813	ENDDO
814
815	end if ! mode 2
816
817	! Mode 2p
818	if (iaero_venus2p .ne.0) then
819	iaer=iaero_venus2p
820
821	! 1. Initialization
822	aerosol(1:ngrid,1:nlayer,iaer)=0.0
823	p_bot = 1.e5
824	p_top = 2.3e4
825	h_bot = 0.1e3
826	h_top = 1.0e3
827
828	! 2. Opacity calculation
829
830	DO ig=1,ngrid
831	DO l=1,nlayer-1
832
833	h_lay=8.31pt(ig,l)/(g0.044)
834
835	!! 1. below 2e5 Pa: no aerosols
836	IF (pplay(ig,l) .gt. 2.e5) THEN
837	mode_dens = 0.
838
839	!! 2. cloud layer
840	ELSEIF (pplay(ig,l) .le. 2.e5 .and. pplay(ig,l) .gt. p_bot) THEN
841	mode_dens = 5.00e7(p_bot/pplay(ig,l))*(h_lay/h_bot)
842
843	ELSEIF (pplay(ig,l) .le. p_bot .and. pplay(ig,l) .gt. p_top) THEN
844	mode_dens = 5.00e7
845
846	ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. 1.e2) THEN
847	mode_dens = 5.00e7(pplay(ig,l)/p_top)*(h_lay/h_top)
848
849	!! 3. above 1.e2 Pa: no aerosols
850	ELSEIF (pplay(ig,l) .le. 1.e2) THEN
851	mode_dens = 0.
852	ENDIF
853
854	aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)* &
855	pi(reffrad(ig,l,iaer))2exp(-3.log(1+nueffrad(ig,l,iaer))) &
856	mode_densh_lay(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
857
858	ENDDO
859	ENDDO
860
861	end if ! mode 2p
862
863	! Mode 3
864	if (iaero_venus3 .ne.0) then
865	iaer=iaero_venus3
866
867	! 1. Initialization
868	aerosol(1:ngrid,1:nlayer,iaer)=0.0
869	p_bot = 1.e5
870	p_top = 4.e4
871	h_bot = 0.5e3
872	h_top = 1.0e3
873
874	! 2. Opacity calculation
875
876	DO ig=1,ngrid
877	DO l=1,nlayer-1
878
879	h_lay=8.31pt(ig,l)/(g0.044)
880
881	!! 1. below 2e5 Pa: no aerosols
882	IF (pplay(ig,l) .gt. 2.e5) THEN
883	mode_dens = 0.
884
885	!! 2. cloud layer
886	ELSEIF (pplay(ig,l) .le. 2.e5 .and. pplay(ig,l) .gt. p_bot) THEN
887	mode_dens = 1.40e7(p_bot/pplay(ig,l))*(h_lay/h_bot)
888
889	ELSEIF (pplay(ig,l) .le. p_bot .and. pplay(ig,l) .gt. p_top) THEN
890	mode_dens = 1.40e7
891
892	ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. 1.e2) THEN
893	mode_dens = 1.40e7(pplay(ig,l)/p_top)*(h_lay/h_top)
894
895	!! 3. above 1.e2 Pa: no aerosols
896	ELSEIF (pplay(ig,l) .le. 1.e2) THEN
897	mode_dens = 0.
898	ENDIF
899
900	aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)* &
901	pi(reffrad(ig,l,iaer))2exp(-3.log(1+nueffrad(ig,l,iaer))) &
902	mode_densh_lay(pplev(ig,l)-pplev(ig,l+1))/pplay(ig,l)
903
904	ENDDO
905	ENDDO
906
907	end if ! mode 3
908
909	! UV absorber
910	if (iaero_venusUV .ne.0) then
911	iaer=iaero_venusUV
912
913	! 1. Initialization
914	aerosol(1:ngrid,1:nlayer,iaer)=0.0
915	p_bot = 3.3e4 ! 58 km
916	p_top = 3.7e3 ! 70 km
917	h_bot = 1.0e3
918	h_top = 1.0e3
919
920	! 2. Opacity calculation
921
922	DO ig=1,ngrid
923	DO l=1,nlayer-1
924
925	h_lay=8.31pt(ig,l)/(g0.044)
926
927	!! 1. below 7e4 Pa: no aerosols
928	IF (pplay(ig,l) .gt. 7.e4) THEN
929	mode_dens = 0.
930
931	!! 2. cloud layer
932	ELSEIF (pplay(ig,l) .le. 7.e4 .and. pplay(ig,l) .gt. p_bot) THEN
933	mode_dens = 1.00e7(p_bot/pplay(ig,l))*(h_lay/h_bot)
934
935	ELSEIF (pplay(ig,l) .le. p_bot .and. pplay(ig,l) .gt. p_top) THEN
936	mode_dens = 1.00e7
937
938	ELSEIF (pplay(ig,l) .le. p_top .and. pplay(ig,l) .gt. 1.e3) THEN
939	mode_dens = 1.00e7(pplay(ig,l)/p_top)*(h_lay/h_top)
940
941	!! 3. above 1.e3 Pa: no aerosols
942	ELSEIF (pplay(ig,l) .le. 1.e3) THEN
943	mode_dens = 0.
944	ENDIF
945
946	! normalized to 0.35 microns (peak of absorption)
947	aerosol(ig,l,iaer) = QREFvis3d(ig,l,iaer)*mode_dens
948
949	ENDDO
950	ENDDO
951
952	! 3. Re-normalize to Haus et al 2015 total column optical depth
953	tau_col(:)=0.0
954	DO l=1,nlayer
955	DO ig=1,ngrid
956	tau_col(ig) = tau_col(ig) &
957	+ aerosol(ig,l,iaer)
958	ENDDO
959	ENDDO
960	DO ig=1,ngrid
961	DO l=1,nlayer-1
962	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)*0.205/tau_col(ig)
963	ENDDO
964	ENDDO
965
966	end if ! UV absorber
967
968	!==================================================================
969	! ig=10
970	! do l=1,nlayer
971	! print,8.31pt(ig,l)/(g*0.044),pplay(ig,l),aerosol(ig,l,1),aerosol(ig,l,2),aerosol(ig,l,3),aerosol(ig,l,4)
972	! print*,l,pplay(ig,l),aerosol(ig,l,5)
973	! enddo
974	! stop
975	!==================================================================
976
977
978	! --------------------------------------------------------------------------
979	! Column integrated visible optical depth in each point (used for diagnostic)
980
981	tau_col(:)=0.0
982	do iaer = 1, naerkind
983	do l=1,nlayer
984	do ig=1,ngrid
985	tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
986	end do
987	end do
988	end do
989
990	! do ig=1,ngrid
991	! do l=1,nlayer
992	! do iaer = 1, naerkind
993	! if(aerosol(ig,l,iaer).gt.1.e3)then
994	! print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
995	! print*,'at ig=',ig,', l=',l,', iaer=',iaer
996	! print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
997	! print*,'reffrad=',reffrad(ig,l,iaer)
998	! endif
999	! end do
1000	! end do
1001	! end do
1002
1003	! do ig=1,ngrid
1004	! if(tau_col(ig).gt.1.e3)then
1005	! print*,'WARNING: tau_col=',tau_col(ig)
1006	! print*,'at ig=',ig
1007	! print*,'aerosol=',aerosol(ig,:,:)
1008	! print*,'QREFvis3d=',QREFvis3d(ig,:,:)
1009	! print*,'reffrad=',reffrad(ig,:,:)
1010	! endif
1011	! end do
1012
1013	end subroutine aeropacity
1014

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

Download in other formats:

Original Format