Context Navigation

aeropacity.F90 @ 2613

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

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

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

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

+Extra info :

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

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

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

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

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

File size: 26.6 KB

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

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

Download in other formats:

Original Format