Context Navigation

aeropacity.F90 @ 2176

Last change on this file since 2176 was 1988, checked in by jleconte, 6 years ago

28/08/2018 == JL

We now shift the radiative model top from p=0 to the middle of the last physical layer. This is done by changing pmid and plevrad in callcorrk and some corrections need to be done in gfluxv.
This seems to get rid of the aratic temperature behavior in the last two layers of the model (especially on the night side on synchronous planets).
Additional speedup corrections have been made in gfluxi that change nothing to the result.
Finally, 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 and changes in last commit.)
This has been done for water ice in aeropacity, but same correction should probably be done for other aerosol types.

File size: 20.8 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, &
	13	tau_nh3_cloud, pres_nh3_cloud
[135]	14
	15	implicit none
	16
	17	!==================================================================
	18	!
	19	! Purpose
	20	! -------
	21	! Compute aerosol optical depth in each gridbox.
	22	!
	23	! Authors
	24	! -------
	25	! F. Forget
	26	! F. Montmessin (water ice scheme)
	27	! update J.-B. Madeleine (2008)
	28	! dust removal, simplification by Robin Wordsworth (2009)
	29	!
	30	! Input
	31	! -----
	32	! ngrid Number of horizontal gridpoints
	33	! nlayer Number of layers
	34	! nq Number of tracers
	35	! pplev Pressure (Pa) at each layer boundary
	36	! pq Aerosol mixing ratio
	37	! reffrad(ngrid,nlayer,naerkind) Aerosol effective radius
[1308]	38	! QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
	39	! QREFir3d(ngrid,nlayer,naerkind) / at reference wavelengths
[135]	40	!
	41	! Output
	42	! ------
	43	! aerosol Aerosol optical depth in layer l, grid point ig
	44	! tau_col Total column optical depth at grid point ig
	45	!
	46	!=======================================================================
	47
[858]	48	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
	49	INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
	50	INTEGER,INTENT(IN) :: nq ! number of tracers
	51	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
	52	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
	53	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
	54	REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth
	55	REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) ! aerosol effective radius
[1308]	56	REAL,INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind) ! extinction coefficient in the visible
	57	REAL,INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
[858]	58	REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth
	59	! BENJAMIN MODIFS
[1308]	60	real,intent(in) :: cloudfrac(ngrid,nlayer) ! cloud fraction
[858]	61	real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
	62	logical,intent(in) :: clearsky
[135]	63
[1677]	64	real aerosol0, obs_tau_col_aurora, pm, pente_cloud
[253]	65
[135]	66	INTEGER l,ig,iq,iaer
	67
[858]	68	LOGICAL,SAVE :: firstcall=.true.
[1315]	69	!$OMP THREADPRIVATE(firstcall)
[135]	70	REAL CBRT
	71	EXTERNAL CBRT
	72
	73	INTEGER,SAVE :: i_co2ice=0 ! co2 ice
	74	INTEGER,SAVE :: i_h2oice=0 ! water ice
[1315]	75	!$OMP THREADPRIVATE(i_co2ice,i_h2oice)
[135]	76	CHARACTER(LEN=20) :: tracername ! to temporarily store text
	77
[253]	78	! for fixed dust profiles
	79	real topdust, expfactor, zp
[787]	80	REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
	81	REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling
[858]	82
[728]	83	real CLFtot
[253]	84
	85	! identify tracers
[135]	86	IF (firstcall) THEN
	87
[726]	88	write(,) "Tracers found in aeropacity:"
[787]	89	do iq=1,nq
[135]	90	tracername=noms(iq)
	91	if (tracername.eq."co2_ice") then
	92	i_co2ice=iq
[726]	93	write(,) "i_co2ice=",i_co2ice
	94
[135]	95	endif
	96	if (tracername.eq."h2o_ice") then
	97	i_h2oice=iq
[726]	98	write(,) "i_h2oice=",i_h2oice
[135]	99	endif
	100	enddo
	101
[741]	102	if (noaero) then
	103	print*, "No active aerosols found in aeropacity"
	104	else
[726]	105	print*, "If you would like to use aerosols, make sure any old"
	106	print*, "start files are updated in newstart using the option"
	107	print*, "q=0"
[741]	108	write(,) "Active aerosols found in aeropacity:"
[726]	109	endif
[253]	110
[741]	111	if ((iaero_co2.ne.0).and.(.not.noaero)) then
[726]	112	print*, 'iaero_co2= ',iaero_co2
	113	endif
	114	if (iaero_h2o.ne.0) then
	115	print*,'iaero_h2o= ',iaero_h2o
	116	endif
	117	if (iaero_dust.ne.0) then
	118	print*,'iaero_dust= ',iaero_dust
	119	endif
	120	if (iaero_h2so4.ne.0) then
	121	print*,'iaero_h2so4= ',iaero_h2so4
	122	endif
[1026]	123	if (iaero_back2lay.ne.0) then
	124	print*,'iaero_back2lay= ',iaero_back2lay
	125	endif
[1677]	126	if (iaero_nh3.ne.0) then
	127	print*,'iaero_nh3= ',iaero_nh3
	128	endif
	129	if (iaero_aurora.ne.0) then
	130	print*,'iaero_aurora= ',iaero_aurora
	131	endif
[726]	132
[135]	133	firstcall=.false.
	134	ENDIF ! of IF (firstcall)
	135
	136
	137	! ---------------------------------------------------------
	138	!==================================================================
[726]	139	! CO2 ice aerosols
[135]	140	!==================================================================
	141
[728]	142	if (iaero_co2.ne.0) then
[726]	143	iaer=iaero_co2
[135]	144	! 1. Initialization
[1308]	145	aerosol(1:ngrid,1:nlayer,iaer)=0.0
[135]	146	! 2. Opacity calculation
[804]	147	if (noaero) then ! aerosol set to zero
[1308]	148	aerosol(1:ngrid,1:nlayer,iaer)=0.0
[741]	149	elseif (aerofixco2.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed
[1308]	150	aerosol(1:ngrid,1:nlayer,iaer)=1.e-9
[787]	151	!aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option
[253]	152	else
	153	DO ig=1, ngrid
	154	DO l=1,nlayer-1 ! to stop the rad tran bug
[135]	155
[253]	156	aerosol0 = &
	157	( 0.75 * QREFvis3d(ig,l,iaer) / &
	158	( rho_co2 * reffrad(ig,l,iaer) ) ) * &
	159	( pq(ig,l,i_co2ice) + 1.E-9 ) * &
	160	( pplev(ig,l) - pplev(ig,l+1) ) / g
	161	aerosol0 = max(aerosol0,1.e-9)
	162	aerosol0 = min(aerosol0,L_TAUMAX)
	163	aerosol(ig,l,iaer) = aerosol0
	164	! aerosol(ig,l,iaer) = 0.0
[726]	165	! print*, aerosol(ig,l,iaer)
[253]	166	! using cloud fraction
	167	! aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF))
	168	! aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX)
	169
	170
	171	ENDDO
	172	ENDDO
[726]	173	end if ! if fixed or varying
[728]	174	end if ! if CO2 aerosols
[135]	175	!==================================================================
[726]	176	! Water ice / liquid
[135]	177	!==================================================================
	178
[728]	179	if (iaero_h2o.ne.0) then
[726]	180	iaer=iaero_h2o
[135]	181	! 1. Initialization
[1308]	182	aerosol(1:ngrid,1:nlayer,iaer)=0.0
[135]	183	! 2. Opacity calculation
[726]	184	if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
[1308]	185	aerosol(1:ngrid,1:nlayer,iaer) =1.e-9
[305]	186
	187	! put cloud at cloudlvl
	188	if(kastprof.and.(cloudlvl.ne.0.0))then
	189	ig=1
	190	do l=1,nlayer
	191	if(int(cloudlvl).eq.l)then
	192	!if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
	193	print*,'Inserting cloud at level ',l
	194	!aerosol(ig,l,iaer)=10.0
	195
	196	rho_ice=920.0
	197
	198	! the Kasting approximation
	199	aerosol(ig,l,iaer) = &
	200	( 0.75 * QREFvis3d(ig,l,iaer) / &
	201	( rho_ice * reffrad(ig,l,iaer) ) ) * &
	202	!( pq(ig,l,i_h2oice) + 1.E-9 ) * &
	203	( 4.0e-4 + 1.E-9 ) * &
	204	( pplev(ig,l) - pplev(ig,l+1) ) / g
	205
	206
	207	open(115,file='clouds.out',form='formatted')
	208	write(115,*) l,aerosol(ig,l,iaer)
	209	close(115)
	210
	211	return
	212	endif
	213	end do
	214
	215	call abort
	216	endif
	217
[253]	218	else
[728]	219
[253]	220	do ig=1, ngrid
[1988]	221	!do l=1,nlayer-1 ! to stop the rad tran bug
	222	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)
	223	! same correction should b-probably be done for other aerosol types.
[728]	224	aerosol(ig,l,iaer) = & !modification by BC
[253]	225	( 0.75 * QREFvis3d(ig,l,iaer) / &
	226	( rho_ice * reffrad(ig,l,iaer) ) ) * &
[716]	227	! pq(ig,l,i_h2oice) * & !JL I dropped the +1e-9 here to have the same
[728]	228	!( pplev(ig,l) - pplev(ig,l+1) ) / g ! opacity in the clearsky=true and the
	229	! clear=false/pq=0 case
[716]	230	( pq(ig,l,i_h2oice) + 1.E-9 ) * & ! Doing this makes the code unstable, so I have restored it (RW)
[728]	231	( pplev(ig,l) - pplev(ig,l+1) ) / g
[135]	232
[253]	233	enddo
	234	enddo
[135]	235
[728]	236	if(CLFvarying)then
[1308]	237	call totalcloudfrac(ngrid,nlayer,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1,1,iaer))
[728]	238	do ig=1, ngrid
[1988]	239	!do l=1,nlayer-1 ! to stop the rad tran bug
	240	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]	241	CLFtot = max(totcloudfrac(ig),0.01)
	242	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
	243	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
	244	enddo
	245	enddo
	246	else
	247	do ig=1, ngrid
[1988]	248	!do l=1,nlayer-1 ! to stop the rad tran bug
	249	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]	250	CLFtot = CLFfixval
	251	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
	252	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
	253	enddo
	254	enddo
	255	end if!(CLFvarying)
	256	endif !(aerofixed.or.(i_h2oice.eq.0).or.clearsky)
	257
	258	end if ! End if h2o aerosol
	259
[726]	260	!==================================================================
	261	! Dust
	262	!==================================================================
[728]	263	if (iaero_dust.ne.0) then
[726]	264	iaer=iaero_dust
	265	! 1. Initialization
[1308]	266	aerosol(1:ngrid,1:nlayer,iaer)=0.0
[726]	267
[728]	268	topdust=30.0 ! km (used to be 10.0 km) LK
[726]	269
	270	! 2. Opacity calculation
	271
	272	! expfactor=0.
	273	DO l=1,nlayer-1
	274	DO ig=1,ngrid
	275	! Typical mixing ratio profile
	276
[1308]	277	zp=(pplev(ig,1)/pplay(ig,l))**(70./topdust)
[726]	278	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
	279
	280	! Vertical scaling function
	281	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
	282	*expfactor
	283
	284
	285	ENDDO
	286	ENDDO
	287
	288	! Rescaling each layer to reproduce the choosen (or assimilated)
	289	! dust extinction opacity at visible reference wavelength, which
[1308]	290	! is scaled to the surface pressure pplev(ig,1)
[726]	291
	292	taudusttmp(1:ngrid)=0.
	293	DO l=1,nlayer
	294	DO ig=1,ngrid
	295	taudusttmp(ig) = taudusttmp(ig) &
	296	+ aerosol(ig,l,iaer)
	297	ENDDO
	298	ENDDO
	299	DO l=1,nlayer-1
	300	DO ig=1,ngrid
	301	aerosol(ig,l,iaer) = max(1E-20, &
	302	dusttau &
[1308]	303	* pplev(ig,1) / pplev(ig,1) &
[726]	304	* aerosol(ig,l,iaer) &
	305	/ taudusttmp(ig))
	306
	307	ENDDO
	308	ENDDO
[728]	309	end if ! If dust aerosol
[726]	310
[135]	311	!==================================================================
[726]	312	! H2SO4
[253]	313	!==================================================================
[726]	314	! added by LK
	315	if (iaero_h2so4.ne.0) then
[728]	316	iaer=iaero_h2so4
[135]	317
[253]	318	! 1. Initialization
[1308]	319	aerosol(1:ngrid,1:nlayer,iaer)=0.0
[253]	320
	321
	322	! 2. Opacity calculation
	323
[726]	324	! expfactor=0.
[728]	325	DO l=1,nlayer-1
	326	DO ig=1,ngrid
	327	! Typical mixing ratio profile
[253]	328
[1308]	329	zp=(pplev(ig,1)/pplay(ig,l))**(70./30) !emulating topdust
[728]	330	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
[135]	331
[726]	332	! Vertical scaling function
[728]	333	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor
[135]	334
[728]	335	ENDDO
	336	ENDDO
	337	tauh2so4tmp(1:ngrid)=0.
	338	DO l=1,nlayer
	339	DO ig=1,ngrid
	340	tauh2so4tmp(ig) = tauh2so4tmp(ig) + aerosol(ig,l,iaer)
	341	ENDDO
	342	ENDDO
	343	DO l=1,nlayer-1
	344	DO ig=1,ngrid
	345	aerosol(ig,l,iaer) = max(1E-20, &
[726]	346	1 &
[1308]	347	* pplev(ig,1) / pplev(ig,1) &
[726]	348	* aerosol(ig,l,iaer) &
	349	/ tauh2so4tmp(ig))
	350
	351	ENDDO
[728]	352	ENDDO
[726]	353
	354	! 1/700. is assuming a "sulfurtau" of 1
	355	! Sulfur aerosol routine to be improved.
	356	! aerosol0 = &
	357	! ( 0.75 * QREFvis3d(ig,l,iaer) / &
	358	! ( rho_h2so4 * reffrad(ig,l,iaer) ) ) * &
	359	! ( pq(ig,l,i_h2so4) + 1.E-9 ) * &
	360	! ( pplev(ig,l) - pplev(ig,l+1) ) / g
	361	! aerosol0 = max(aerosol0,1.e-9)
	362	! aerosol0 = min(aerosol0,L_TAUMAX)
	363	! aerosol(ig,l,iaer) = aerosol0
	364
	365	! ENDDO
	366	! ENDDO
	367	end if
	368
	369
[1026]	370	! ---------------------------------------------------------
	371	!==================================================================
	372	! Two-layer aerosols (unknown composition)
	373	! S. Guerlet (2013)
	374	!==================================================================
[726]	375
[1026]	376	if (iaero_back2lay .ne.0) then
	377	iaer=iaero_back2lay
	378	! 1. Initialization
[1308]	379	aerosol(1:ngrid,1:nlayer,iaer)=0.0
[1026]	380	! 2. Opacity calculation
	381	DO ig=1,ngrid
	382	DO l=1,nlayer-1
	383	aerosol(ig,l,iaer) = ( pplev(ig,l) - pplev(ig,l+1) )
	384	!! 1. below tropospheric layer: no aerosols
	385	IF (pplev(ig,l) .gt. pres_bottom_tropo) THEN
	386	aerosol(ig,l,iaer) = 0.*aerosol(ig,l,iaer)
	387	!! 2. tropo layer
	388	ELSEIF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN
	389	aerosol(ig,l,iaer) = obs_tau_col_tropo*aerosol(ig,l,iaer)
	390	!! 3. linear transition
	391	ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN
	392	expfactor=log(obs_tau_col_strato/obs_tau_col_tropo)/log(pres_bottom_strato/pres_top_tropo)
[1132]	393	aerosol(ig,l,iaer)= obs_tau_col_tropo((pplev(ig,l)/pres_top_tropo)expfactor)aerosol(ig,l,iaer)/1.5
[1026]	394	!! 4. strato layer
	395	ELSEIF (pplev(ig,l) .le. pres_bottom_strato .and. pplev(ig,l) .gt. pres_top_strato) THEN
	396	aerosol(ig,l,iaer)= obs_tau_col_strato*aerosol(ig,l,iaer)
	397	!! 5. above strato layer: no aerosols
	398	ELSEIF (pplev(ig,l) .lt. pres_top_strato) THEN
	399	aerosol(ig,l,iaer) = 0.*aerosol(ig,l,iaer)
	400	ENDIF
	401	ENDDO
	402	ENDDO
	403
	404	! 3. Re-normalize to observed total column
	405	tau_col(:)=0.0
	406	DO l=1,nlayer
	407	DO ig=1,ngrid
	408	tau_col(ig) = tau_col(ig) &
	409	+ aerosol(ig,l,iaer)/(obs_tau_col_tropo+obs_tau_col_strato)
	410	ENDDO
	411	ENDDO
	412
	413	DO ig=1,ngrid
	414	DO l=1,nlayer-1
	415	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig)
	416	ENDDO
	417	ENDDO
	418
	419
	420	end if ! if Two-layer aerosols
	421
[1677]	422	!==================================================================
	423	! Saturn/Jupiter ammonia cloud = thin cloud (scale height 0.2 hard coded...)
	424	! S. Guerlet (2013)
	425	!==================================================================
[1026]	426
[1677]	427	if (iaero_nh3 .ne.0) then
	428	iaer=iaero_nh3
	429	! 1. Initialization
	430	aerosol(1:ngrid,1:nlayer,iaer)=0.D0
	431	! 2. Opacity calculation
	432	DO ig=1,ngrid
	433
	434	DO l=1,nlayer-1
	435	!! 1. below cloud layer: no opacity
	436
	437	IF (pplev(ig,l) .gt. pres_nh3_cloud ) THEN
	438	aerosol(ig,l,iaer) = 0.D0
	439
	440	ELSEIF (pplev(ig,l) .le. pres_nh3_cloud ) THEN
	441	pente_cloud=5. !!(hard-coded, correspond to scale height 0.2)
	442	aerosol(ig,l,iaer) = ((pplev(ig,l)/pres_nh3_cloud)*(pente_cloud))tau_nh3_cloud
	443
	444	ENDIF
	445	ENDDO
	446
	447	END DO
	448
	449	! 3. Re-normalize to observed total column
	450	tau_col(:)=0.0
	451	DO l=1,nlayer
	452	DO ig=1,ngrid
	453	tau_col(ig) = tau_col(ig) &
	454	+ aerosol(ig,l,iaer)/tau_nh3_cloud
	455	ENDDO
	456	ENDDO
	457
	458	DO ig=1,ngrid
	459	DO l=1,nlayer-1
	460	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig)
	461	ENDDO
	462	ENDDO
	463
	464	end if ! if NH3 cloud
	465	!==================================================================
	466	! Jovian auroral aerosols (unknown composition) NON-GENERIC: vertical and meridional profile tuned to observations
	467	! S. Guerlet (2015)
	468	!==================================================================
	469
	470
	471	if (iaero_aurora .ne.0) then
	472	iaer=iaero_aurora
	473	! 1. Initialization
	474	aerosol(1:ngrid,1:nlayer,iaer)=0.D0
	475	pm = 2000. !!case study: maxi aerosols at 20 hPa
	476	! 2. Opacity calculation
	477	DO ig=1,ngrid
	478
	479	!! Test Jupiter (based on Zhang et al 2013 observations, but a bit different), decembre 2015
	480	DO l=1,nlayer
	481	aerosol(ig,l,iaer) = (pplev(ig,l)/pm)*2 exp(-(pplev(ig,l)/pm)**2)
	482	ENDDO
	483	ENDDO
	484
	485	! 3. Meridional distribution, and re-normalize to observed total column
	486	tau_col(:)=0.D0
	487	DO ig=1,ngrid
	488	!!Jupiter
	489	!!Hem sud:
	490	IF (latitude(ig)180.D0/pi .lt. -45.D0 .and. latitude(ig)180.D0/pi .gt. -70.) THEN
	491	obs_tau_col_aurora= 10.D0*(-0.06D0latitude(ig)*180.D0/pi-3.4D0)
	492	ELSEIF (latitude(ig)180.D0/pi .lt. -37.D0 .and. latitude(ig)180.D0/pi .ge. -45.) THEN
	493	obs_tau_col_aurora= 10.D0*(-0.3D0latitude(ig)*180.D0/pi-14.3D0)
	494	ELSEIF (latitude(ig)*180./pi .le. -70. ) THEN
	495	obs_tau_col_aurora= 10*(0.0670.-3.4D0)
	496	!!Hem Nord:
	497	ELSEIF (latitude(ig)180.D0/pi .gt. 30.D0 .and. latitude(ig)180.D0/pi .lt. 70.) THEN
	498	obs_tau_col_aurora= 10.D0*(0.03D0latitude(ig)*180.D0/pi-1.17D0)
	499	ELSEIF (latitude(ig)180.D0/pi .gt. 22.D0 .and. latitude(ig)180.D0/pi .le. 30.) THEN
	500	obs_tau_col_aurora= 10.D0*(0.3D0latitude(ig)*180.D0/pi-9.4D0)
	501	ELSEIF (latitude(ig)*180.D0/pi .ge. 70.) THEN
	502	obs_tau_col_aurora= 10*(0.0370.-1.17D0)
	503	ELSEIF (latitude(ig)180.D0/pi .ge. -37. .and. latitude(ig)180.D0/pi .le. 22.) THEN
	504	obs_tau_col_aurora = 0.001D0 !!Jupiter: mini pas a zero
	505	ENDIF
	506
	507	DO l=1,nlayer
	508	tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)/obs_tau_col_aurora
	509	ENDDO
	510	ENDDO
	511
	512	DO ig=1,ngrid
	513	DO l=1,nlayer-1
	514	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig)
	515	ENDDO
	516	ENDDO
	517
	518
	519	end if ! if Auroral aerosols
	520
	521
[135]	522	! --------------------------------------------------------------------------
	523	! Column integrated visible optical depth in each point (used for diagnostic)
	524
[253]	525	tau_col(:)=0.0
[135]	526	do iaer = 1, naerkind
	527	do l=1,nlayer
	528	do ig=1,ngrid
	529	tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
	530	end do
	531	end do
	532	end do
	533
[787]	534	do ig=1,ngrid
[253]	535	do l=1,nlayer
	536	do iaer = 1, naerkind
	537	if(aerosol(ig,l,iaer).gt.1.e3)then
	538	print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
	539	print*,'at ig=',ig,', l=',l,', iaer=',iaer
	540	print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
	541	print*,'reffrad=',reffrad(ig,l,iaer)
	542	endif
	543	end do
	544	end do
	545	end do
	546
[787]	547	do ig=1,ngrid
[253]	548	if(tau_col(ig).gt.1.e3)then
	549	print*,'WARNING: tau_col=',tau_col(ig)
	550	print*,'at ig=',ig
	551	print*,'aerosol=',aerosol(ig,:,:)
	552	print*,'QREFvis3d=',QREFvis3d(ig,:,:)
	553	print*,'reffrad=',reffrad(ig,:,:)
	554	endif
	555	end do
[135]	556	return
	557	end subroutine aeropacity
	558

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

Download in other formats:

Original Format