Context Navigation

aeropacity.F90 @ 3436

Last change on this file since 3436 was 2899, checked in by emillour, 22 months ago
Generic PCM: More code tidying: turn aeropacity, aeroptproperties, gfluxi, gfluxv, sfluxi and sfluxv into modules. EM
File size: 37.9 KB

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

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

Download in other formats:

Original Format