Context Navigation

aeropacity.F90 @ 3981

Last change on this file since 3981 was 3922, checked in by mmaurice, 2 months ago

Generic PCM:

Add time-dependent dust scenario for (present-day) Mars from Montmessin
et al., 2004.

File size: 40.6 KB

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

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

Download in other formats:

Original Format