Context Navigation

aeropacity.F90 @ 3300

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

Rev	Line
[2899]	1	module aeropacity_mod
	2
	3	implicit none
	4
	5	contains
	6
[2831]	7	Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pt, pq, &
	8	aerosol,reffrad,nueffrad, QREFvis3d,QREFir3d,tau_col, &
	9	cloudfrac,totcloudfrac,clearsky)
[135]	10
[726]	11	use radinc_h, only : L_TAUMAX,naerkind
[2810]	12	use aerosol_mod, only: iaero_nlay, iaero_generic, &
	13	iaero_aurora, iaero_back2lay, iaero_co2, &
	14	iaero_dust, iaero_h2o, iaero_h2so4, &
[2831]	15	iaero_nh3, i_rgcs_ice, noaero, &
	16	iaero_venus1, iaero_venus2, iaero_venus2p, &
	17	iaero_venus3, iaero_venusUV
[2803]	18	USE tracer_h, only: noms,rho_co2,rho_ice,rho_q,mmol
	19	use comcstfi_mod, only: g, pi, mugaz, avocado
[1677]	20	use geometry_mod, only: latitude
[1397]	21	use callkeys_mod, only: aerofixco2,aerofixh2o,kastprof,cloudlvl, &
	22	CLFvarying,CLFfixval,dusttau, &
	23	pres_bottom_tropo,pres_top_tropo,obs_tau_col_tropo, &
[1677]	24	pres_bottom_strato,pres_top_strato,obs_tau_col_strato, &
[2297]	25	tau_nh3_cloud, pres_nh3_cloud, &
	26	nlayaero, aeronlay_tauref, aeronlay_choice, &
[2803]	27	aeronlay_pbot, aeronlay_ptop, aeronlay_sclhght, &
	28	aerogeneric
	29	use generic_tracer_index_mod, only: generic_tracer_index
[135]	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)
[2297]	44	! Generic n-layer aerosol - J. Vatant d'Ollone (2020)
[2804]	45	! Radiative Generic Condensable Species - Lucas Teinturier (2022)
[135]	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
[1308]	55	! QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
	56	! QREFir3d(ngrid,nlayer,naerkind) / at reference wavelengths
[135]	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
[858]	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)
[2831]	71	REAL,INTENT(IN) :: pt(ngrid,nlayer) ! mid-layer temperature (K)
[858]	72	REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth
	73	REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) ! aerosol effective radius
[2831]	74	REAL,INTENT(IN) :: nueffrad(ngrid,nlayer,naerkind) ! aerosol effective variance
[1308]	75	REAL,INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind) ! extinction coefficient in the visible
	76	REAL,INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
[858]	77	REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth
	78	! BENJAMIN MODIFS
[1308]	79	real,intent(in) :: cloudfrac(ngrid,nlayer) ! cloud fraction
[858]	80	real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
	81	logical,intent(in) :: clearsky
[135]	82
[2297]	83	real aerosol0, obs_tau_col_aurora, pm
	84	real pcloud_deck, cloud_slope
	85
[2254]	86	real dp_strato(ngrid)
	87	real dp_tropo(ngrid)
[2297]	88	real dp_layer(ngrid)
[253]	89
[2297]	90	INTEGER l,ig,iq,iaer,ia
[135]	91
[858]	92	LOGICAL,SAVE :: firstcall=.true.
[1315]	93	!$OMP THREADPRIVATE(firstcall)
[135]	94	REAL CBRT
	95	EXTERNAL CBRT
	96
	97	INTEGER,SAVE :: i_co2ice=0 ! co2 ice
	98	INTEGER,SAVE :: i_h2oice=0 ! water ice
[1315]	99	!$OMP THREADPRIVATE(i_co2ice,i_h2oice)
[135]	100	CHARACTER(LEN=20) :: tracername ! to temporarily store text
	101
[253]	102	! for fixed dust profiles
	103	real topdust, expfactor, zp
[787]	104	REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
	105	REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling
[858]	106
[728]	107	real CLFtot
[2803]	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
[2804]	110	! integer i_rgcs_ice(aerogeneric)
[2831]	111	! for venus clouds
	112	real :: p_bot,p_top,h_bot,h_top,mode_dens,h_lay
	113
[253]	114	! identify tracers
[135]	115	IF (firstcall) THEN
[2804]	116	ia =0
[726]	117	write(,) "Tracers found in aeropacity:"
[787]	118	do iq=1,nq
[135]	119	tracername=noms(iq)
	120	if (tracername.eq."co2_ice") then
	121	i_co2ice=iq
[726]	122	write(,) "i_co2ice=",i_co2ice
	123
[135]	124	endif
	125	if (tracername.eq."h2o_ice") then
	126	i_h2oice=iq
[726]	127	write(,) "i_h2oice=",i_h2oice
[135]	128	endif
	129	enddo
	130
[741]	131	if (noaero) then
	132	print*, "No active aerosols found in aeropacity"
	133	else
[726]	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"
[741]	137	write(,) "Active aerosols found in aeropacity:"
[726]	138	endif
[253]	139
[741]	140	if ((iaero_co2.ne.0).and.(.not.noaero)) then
[726]	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
[1026]	152	if (iaero_back2lay.ne.0) then
	153	print*,'iaero_back2lay= ',iaero_back2lay
	154	endif
[1677]	155	if (iaero_nh3.ne.0) then
	156	print*,'iaero_nh3= ',iaero_nh3
	157	endif
[2297]	158	if (iaero_nlay(1).ne.0) then
	159	print*,'iaero_nlay= ',iaero_nlay(:)
	160	endif
[1677]	161	if (iaero_aurora.ne.0) then
	162	print*,'iaero_aurora= ',iaero_aurora
	163	endif
[2831]	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
[2810]	181	if (aerogeneric .ne. 0) then
[2803]	182	print*,"iaero_generic= ",iaero_generic(:)
	183	endif
[135]	184	firstcall=.false.
	185	ENDIF ! of IF (firstcall)
	186
	187
	188	! ---------------------------------------------------------
	189	!==================================================================
[726]	190	! CO2 ice aerosols
[135]	191	!==================================================================
	192
[728]	193	if (iaero_co2.ne.0) then
[726]	194	iaer=iaero_co2
[135]	195	! 1. Initialization
[1308]	196	aerosol(1:ngrid,1:nlayer,iaer)=0.0
[135]	197	! 2. Opacity calculation
[804]	198	if (noaero) then ! aerosol set to zero
[1308]	199	aerosol(1:ngrid,1:nlayer,iaer)=0.0
[741]	200	elseif (aerofixco2.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed
[1308]	201	aerosol(1:ngrid,1:nlayer,iaer)=1.e-9
[787]	202	!aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option
[253]	203	else
	204	DO ig=1, ngrid
	205	DO l=1,nlayer-1 ! to stop the rad tran bug
[135]	206
[253]	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
[726]	216	! print*, aerosol(ig,l,iaer)
[253]	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
[726]	224	end if ! if fixed or varying
[728]	225	end if ! if CO2 aerosols
[135]	226	!==================================================================
[726]	227	! Water ice / liquid
[135]	228	!==================================================================
	229
[728]	230	if (iaero_h2o.ne.0) then
[726]	231	iaer=iaero_h2o
[135]	232	! 1. Initialization
[1308]	233	aerosol(1:ngrid,1:nlayer,iaer)=0.0
[135]	234	! 2. Opacity calculation
[726]	235	if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
[1308]	236	aerosol(1:ngrid,1:nlayer,iaer) =1.e-9
[305]	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
[253]	269	else
[728]	270
[253]	271	do ig=1, ngrid
[1988]	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.
[728]	275	aerosol(ig,l,iaer) = & !modification by BC
[253]	276	( 0.75 * QREFvis3d(ig,l,iaer) / &
	277	( rho_ice * reffrad(ig,l,iaer) ) ) * &
[716]	278	! pq(ig,l,i_h2oice) * & !JL I dropped the +1e-9 here to have the same
[728]	279	!( pplev(ig,l) - pplev(ig,l+1) ) / g ! opacity in the clearsky=true and the
	280	! clear=false/pq=0 case
[716]	281	( pq(ig,l,i_h2oice) + 1.E-9 ) * & ! Doing this makes the code unstable, so I have restored it (RW)
[728]	282	( pplev(ig,l) - pplev(ig,l+1) ) / g
[135]	283
[253]	284	enddo
	285	enddo
[135]	286
[728]	287	if(CLFvarying)then
[1308]	288	call totalcloudfrac(ngrid,nlayer,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1,1,iaer))
[728]	289	do ig=1, ngrid
[1988]	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)
[728]	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
[1988]	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)
[728]	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
[726]	311	!==================================================================
	312	! Dust
	313	!==================================================================
[728]	314	if (iaero_dust.ne.0) then
[726]	315	iaer=iaero_dust
	316	! 1. Initialization
[1308]	317	aerosol(1:ngrid,1:nlayer,iaer)=0.0
[726]	318
[728]	319	topdust=30.0 ! km (used to be 10.0 km) LK
[726]	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
[1308]	328	zp=(pplev(ig,1)/pplay(ig,l))**(70./topdust)
[726]	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
[1308]	341	! is scaled to the surface pressure pplev(ig,1)
[726]	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 &
[1308]	354	* pplev(ig,1) / pplev(ig,1) &
[726]	355	* aerosol(ig,l,iaer) &
	356	/ taudusttmp(ig))
	357
	358	ENDDO
	359	ENDDO
[728]	360	end if ! If dust aerosol
[726]	361
[135]	362	!==================================================================
[726]	363	! H2SO4
[253]	364	!==================================================================
[726]	365	! added by LK
	366	if (iaero_h2so4.ne.0) then
[728]	367	iaer=iaero_h2so4
[135]	368
[253]	369	! 1. Initialization
[1308]	370	aerosol(1:ngrid,1:nlayer,iaer)=0.0
[253]	371
	372
	373	! 2. Opacity calculation
	374
[726]	375	! expfactor=0.
[728]	376	DO l=1,nlayer-1
	377	DO ig=1,ngrid
	378	! Typical mixing ratio profile
[253]	379
[1308]	380	zp=(pplev(ig,1)/pplay(ig,l))**(70./30) !emulating topdust
[728]	381	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
[135]	382
[726]	383	! Vertical scaling function
[728]	384	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor
[135]	385
[728]	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, &
[726]	397	1 &
[1308]	398	* pplev(ig,1) / pplev(ig,1) &
[726]	399	* aerosol(ig,l,iaer) &
	400	/ tauh2so4tmp(ig))
	401
	402	ENDDO
[728]	403	ENDDO
[2803]	404
[726]	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
[1026]	421	! ---------------------------------------------------------
	422	!==================================================================
	423	! Two-layer aerosols (unknown composition)
[2254]	424	! S. Guerlet (2013) - Modif by J. Vatant d'Ollone (2020)
[2297]	425	!
	426	! This scheme is deprecated and left for retrocompatibility
	427	! You should use the n-layer scheme below !
	428	!
[1026]	429	!==================================================================
[726]	430
[1026]	431	if (iaero_back2lay .ne.0) then
	432	iaer=iaero_back2lay
	433	! 1. Initialization
[1308]	434	aerosol(1:ngrid,1:nlayer,iaer)=0.0
[1026]	435	! 2. Opacity calculation
[2254]	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
[1026]	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
[2254]	459	aerosol(ig,l,iaer) = 0.D0
[1026]	460	!! 2. tropo layer
	461	ELSEIF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN
[2254]	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)
[1026]	468	!! 5. above strato layer: no aerosols
	469	ELSEIF (pplev(ig,l) .lt. pres_top_strato) THEN
[2254]	470	aerosol(ig,l,iaer) = 0.D0
[1026]	471	ENDIF
[2254]	472	ENDDO
[1026]	473	ENDDO
	474
[2254]	475	! 3. Re-normalize to the (input) observed (total) column (for each of the layers)
	476
[1026]	477	DO ig=1,ngrid
[2254]	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
[1026]	490	ENDDO
	491
	492
	493	end if ! if Two-layer aerosols
	494
[1677]	495	!==================================================================
	496	! Saturn/Jupiter ammonia cloud = thin cloud (scale height 0.2 hard coded...)
	497	! S. Guerlet (2013)
[2297]	498	! JVO 20 : You should now use the generic n-layer scheme below
[1677]	499	!==================================================================
[1026]	500
[1677]	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
[2297]	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
[1677]	517
	518	ENDIF
	519	ENDDO
	520
	521	END DO
	522
	523	! 3. Re-normalize to observed total column
[2297]	524	dp_layer(:)=0.0
[1677]	525	DO l=1,nlayer
	526	DO ig=1,ngrid
[2297]	527	dp_layer(ig) = dp_layer(ig) &
[1677]	528	+ aerosol(ig,l,iaer)/tau_nh3_cloud
	529	ENDDO
	530	ENDDO
	531
	532	DO ig=1,ngrid
	533	DO l=1,nlayer-1
[2297]	534	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/dp_layer(ig)
[1677]	535	ENDDO
	536	ENDDO
	537
	538	end if ! if NH3 cloud
[2297]	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
[1677]	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
[2297]	653	dp_layer(:)=0.D0
[1677]	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
[2297]	675	dp_layer(ig) = dp_layer(ig) + aerosol(ig,l,iaer)/obs_tau_col_aurora
[1677]	676	ENDDO
	677	ENDDO
	678
	679	DO ig=1,ngrid
	680	DO l=1,nlayer-1
[2297]	681	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/dp_layer(ig)
[1677]	682	ENDDO
	683	ENDDO
	684
	685
	686	end if ! if Auroral aerosols
[2803]	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	!===========================================================================
[2810]	694	if (aerogeneric .ne. 0) then ! we enter the scheme
[2803]	695	do ia=1,aerogeneric
	696	iaer = iaero_generic(ia)
	697	! Initialization
	698	aerosol(1:ngrid,1:nlayer,iaer) = 0.D0
[2804]	699	igen_ice = i_rgcs_ice(ia)
[2803]	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
[2810]	710	endif !aerogeneric .ne. 0
[1677]	711
[2831]	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
[135]	984	! --------------------------------------------------------------------------
	985	! Column integrated visible optical depth in each point (used for diagnostic)
	986
[253]	987	tau_col(:)=0.0
[135]	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
[2804]	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
[253]	1008
[2804]	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
[2831]	1018
[135]	1019	end subroutine aeropacity
	1020
[2899]	1021	end module aeropacity_mod

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

Download in other formats:

Original Format