Context Navigation

aeropacity.F90 @ 1493

Last change on this file since 1493 was 1397, checked in by milmd, 10 years ago
In LMDZ.GENERIC replacement of all phystd .h files by module files.
File size: 16.2 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
[787]	6	USE comgeomfi_h
[858]	7	USE tracer_h, only: noms,rho_co2,rho_ice
[1384]	8	use comcstfi_mod, only: g
[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, &
	12	pres_bottom_strato,pres_top_strato,obs_tau_col_strato
[135]	13
	14	implicit none
	15
	16	!==================================================================
	17	!
	18	! Purpose
	19	! -------
	20	! Compute aerosol optical depth in each gridbox.
	21	!
	22	! Authors
	23	! -------
	24	! F. Forget
	25	! F. Montmessin (water ice scheme)
	26	! update J.-B. Madeleine (2008)
	27	! dust removal, simplification by Robin Wordsworth (2009)
	28	!
	29	! Input
	30	! -----
	31	! ngrid Number of horizontal gridpoints
	32	! nlayer Number of layers
	33	! nq Number of tracers
	34	! pplev Pressure (Pa) at each layer boundary
	35	! pq Aerosol mixing ratio
	36	! reffrad(ngrid,nlayer,naerkind) Aerosol effective radius
[1308]	37	! QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
	38	! QREFir3d(ngrid,nlayer,naerkind) / at reference wavelengths
[135]	39	!
	40	! Output
	41	! ------
	42	! aerosol Aerosol optical depth in layer l, grid point ig
	43	! tau_col Total column optical depth at grid point ig
	44	!
	45	!=======================================================================
	46
[858]	47	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
	48	INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
	49	INTEGER,INTENT(IN) :: nq ! number of tracers
	50	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
	51	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
	52	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
	53	REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth
	54	REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) ! aerosol effective radius
[1308]	55	REAL,INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind) ! extinction coefficient in the visible
	56	REAL,INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
[858]	57	REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth
	58	! BENJAMIN MODIFS
[1308]	59	real,intent(in) :: cloudfrac(ngrid,nlayer) ! cloud fraction
[858]	60	real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
	61	logical,intent(in) :: clearsky
[135]	62
[253]	63	real aerosol0
	64
[135]	65	INTEGER l,ig,iq,iaer
	66
[858]	67	LOGICAL,SAVE :: firstcall=.true.
[1315]	68	!$OMP THREADPRIVATE(firstcall)
[135]	69	REAL CBRT
	70	EXTERNAL CBRT
	71
	72	INTEGER,SAVE :: i_co2ice=0 ! co2 ice
	73	INTEGER,SAVE :: i_h2oice=0 ! water ice
[1315]	74	!$OMP THREADPRIVATE(i_co2ice,i_h2oice)
[135]	75	CHARACTER(LEN=20) :: tracername ! to temporarily store text
	76
[253]	77	! for fixed dust profiles
	78	real topdust, expfactor, zp
[787]	79	REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
	80	REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling
[858]	81
[728]	82	real CLFtot
[253]	83
	84	! identify tracers
[135]	85	IF (firstcall) THEN
	86
[726]	87	write(,) "Tracers found in aeropacity:"
[787]	88	do iq=1,nq
[135]	89	tracername=noms(iq)
	90	if (tracername.eq."co2_ice") then
	91	i_co2ice=iq
[726]	92	write(,) "i_co2ice=",i_co2ice
	93
[135]	94	endif
	95	if (tracername.eq."h2o_ice") then
	96	i_h2oice=iq
[726]	97	write(,) "i_h2oice=",i_h2oice
[135]	98	endif
	99	enddo
	100
[741]	101	if (noaero) then
	102	print*, "No active aerosols found in aeropacity"
	103	else
[726]	104	print*, "If you would like to use aerosols, make sure any old"
	105	print*, "start files are updated in newstart using the option"
	106	print*, "q=0"
[741]	107	write(,) "Active aerosols found in aeropacity:"
[726]	108	endif
[253]	109
[741]	110	if ((iaero_co2.ne.0).and.(.not.noaero)) then
[726]	111	print*, 'iaero_co2= ',iaero_co2
	112	endif
	113	if (iaero_h2o.ne.0) then
	114	print*,'iaero_h2o= ',iaero_h2o
	115	endif
	116	if (iaero_dust.ne.0) then
	117	print*,'iaero_dust= ',iaero_dust
	118	endif
	119	if (iaero_h2so4.ne.0) then
	120	print*,'iaero_h2so4= ',iaero_h2so4
	121	endif
[1026]	122	if (iaero_back2lay.ne.0) then
	123	print*,'iaero_back2lay= ',iaero_back2lay
	124	endif
[726]	125
[135]	126	firstcall=.false.
	127	ENDIF ! of IF (firstcall)
	128
	129
	130	! ---------------------------------------------------------
	131	!==================================================================
[726]	132	! CO2 ice aerosols
[135]	133	!==================================================================
	134
[728]	135	if (iaero_co2.ne.0) then
[726]	136	iaer=iaero_co2
[135]	137	! 1. Initialization
[1308]	138	aerosol(1:ngrid,1:nlayer,iaer)=0.0
[135]	139	! 2. Opacity calculation
[804]	140	if (noaero) then ! aerosol set to zero
[1308]	141	aerosol(1:ngrid,1:nlayer,iaer)=0.0
[741]	142	elseif (aerofixco2.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed
[1308]	143	aerosol(1:ngrid,1:nlayer,iaer)=1.e-9
[787]	144	!aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option
[253]	145	else
	146	DO ig=1, ngrid
	147	DO l=1,nlayer-1 ! to stop the rad tran bug
[135]	148
[253]	149	aerosol0 = &
	150	( 0.75 * QREFvis3d(ig,l,iaer) / &
	151	( rho_co2 * reffrad(ig,l,iaer) ) ) * &
	152	( pq(ig,l,i_co2ice) + 1.E-9 ) * &
	153	( pplev(ig,l) - pplev(ig,l+1) ) / g
	154	aerosol0 = max(aerosol0,1.e-9)
	155	aerosol0 = min(aerosol0,L_TAUMAX)
	156	aerosol(ig,l,iaer) = aerosol0
	157	! aerosol(ig,l,iaer) = 0.0
[726]	158	! print*, aerosol(ig,l,iaer)
[253]	159	! using cloud fraction
	160	! aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF))
	161	! aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX)
	162
	163
	164	ENDDO
	165	ENDDO
[726]	166	end if ! if fixed or varying
[728]	167	end if ! if CO2 aerosols
[135]	168	!==================================================================
[726]	169	! Water ice / liquid
[135]	170	!==================================================================
	171
[728]	172	if (iaero_h2o.ne.0) then
[726]	173	iaer=iaero_h2o
[135]	174	! 1. Initialization
[1308]	175	aerosol(1:ngrid,1:nlayer,iaer)=0.0
[135]	176	! 2. Opacity calculation
[726]	177	if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
[1308]	178	aerosol(1:ngrid,1:nlayer,iaer) =1.e-9
[305]	179
	180	! put cloud at cloudlvl
	181	if(kastprof.and.(cloudlvl.ne.0.0))then
	182	ig=1
	183	do l=1,nlayer
	184	if(int(cloudlvl).eq.l)then
	185	!if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
	186	print*,'Inserting cloud at level ',l
	187	!aerosol(ig,l,iaer)=10.0
	188
	189	rho_ice=920.0
	190
	191	! the Kasting approximation
	192	aerosol(ig,l,iaer) = &
	193	( 0.75 * QREFvis3d(ig,l,iaer) / &
	194	( rho_ice * reffrad(ig,l,iaer) ) ) * &
	195	!( pq(ig,l,i_h2oice) + 1.E-9 ) * &
	196	( 4.0e-4 + 1.E-9 ) * &
	197	( pplev(ig,l) - pplev(ig,l+1) ) / g
	198
	199
	200	open(115,file='clouds.out',form='formatted')
	201	write(115,*) l,aerosol(ig,l,iaer)
	202	close(115)
	203
	204	return
	205	endif
	206	end do
	207
	208	call abort
	209	endif
	210
[253]	211	else
[728]	212
[253]	213	do ig=1, ngrid
	214	do l=1,nlayer-1 ! to stop the rad tran bug
[135]	215
[728]	216	aerosol(ig,l,iaer) = & !modification by BC
[253]	217	( 0.75 * QREFvis3d(ig,l,iaer) / &
	218	( rho_ice * reffrad(ig,l,iaer) ) ) * &
[716]	219	! pq(ig,l,i_h2oice) * & !JL I dropped the +1e-9 here to have the same
[728]	220	!( pplev(ig,l) - pplev(ig,l+1) ) / g ! opacity in the clearsky=true and the
	221	! clear=false/pq=0 case
[716]	222	( pq(ig,l,i_h2oice) + 1.E-9 ) * & ! Doing this makes the code unstable, so I have restored it (RW)
[728]	223	( pplev(ig,l) - pplev(ig,l+1) ) / g
[135]	224
[253]	225	enddo
	226	enddo
[135]	227
[728]	228	if(CLFvarying)then
[1308]	229	call totalcloudfrac(ngrid,nlayer,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1,1,iaer))
[728]	230	do ig=1, ngrid
	231	do l=1,nlayer-1 ! to stop the rad tran bug
	232	CLFtot = max(totcloudfrac(ig),0.01)
	233	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
	234	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
	235	enddo
	236	enddo
	237	else
	238	do ig=1, ngrid
	239	do l=1,nlayer-1 ! to stop the rad tran bug
	240	CLFtot = CLFfixval
	241	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
	242	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
	243	enddo
	244	enddo
	245	end if!(CLFvarying)
	246	endif !(aerofixed.or.(i_h2oice.eq.0).or.clearsky)
	247
	248	end if ! End if h2o aerosol
	249
[726]	250	!==================================================================
	251	! Dust
	252	!==================================================================
[728]	253	if (iaero_dust.ne.0) then
[726]	254	iaer=iaero_dust
	255	! 1. Initialization
[1308]	256	aerosol(1:ngrid,1:nlayer,iaer)=0.0
[726]	257
[728]	258	topdust=30.0 ! km (used to be 10.0 km) LK
[726]	259
	260	! 2. Opacity calculation
	261
	262	! expfactor=0.
	263	DO l=1,nlayer-1
	264	DO ig=1,ngrid
	265	! Typical mixing ratio profile
	266
[1308]	267	zp=(pplev(ig,1)/pplay(ig,l))**(70./topdust)
[726]	268	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
	269
	270	! Vertical scaling function
	271	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
	272	*expfactor
	273
	274
	275	ENDDO
	276	ENDDO
	277
	278	! Rescaling each layer to reproduce the choosen (or assimilated)
	279	! dust extinction opacity at visible reference wavelength, which
[1308]	280	! is scaled to the surface pressure pplev(ig,1)
[726]	281
	282	taudusttmp(1:ngrid)=0.
	283	DO l=1,nlayer
	284	DO ig=1,ngrid
	285	taudusttmp(ig) = taudusttmp(ig) &
	286	+ aerosol(ig,l,iaer)
	287	ENDDO
	288	ENDDO
	289	DO l=1,nlayer-1
	290	DO ig=1,ngrid
	291	aerosol(ig,l,iaer) = max(1E-20, &
	292	dusttau &
[1308]	293	* pplev(ig,1) / pplev(ig,1) &
[726]	294	* aerosol(ig,l,iaer) &
	295	/ taudusttmp(ig))
	296
	297	ENDDO
	298	ENDDO
[728]	299	end if ! If dust aerosol
[726]	300
[135]	301	!==================================================================
[726]	302	! H2SO4
[253]	303	!==================================================================
[726]	304	! added by LK
	305	if (iaero_h2so4.ne.0) then
[728]	306	iaer=iaero_h2so4
[135]	307
[253]	308	! 1. Initialization
[1308]	309	aerosol(1:ngrid,1:nlayer,iaer)=0.0
[253]	310
	311
	312	! 2. Opacity calculation
	313
[726]	314	! expfactor=0.
[728]	315	DO l=1,nlayer-1
	316	DO ig=1,ngrid
	317	! Typical mixing ratio profile
[253]	318
[1308]	319	zp=(pplev(ig,1)/pplay(ig,l))**(70./30) !emulating topdust
[728]	320	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
[135]	321
[726]	322	! Vertical scaling function
[728]	323	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor
[135]	324
[728]	325	ENDDO
	326	ENDDO
	327	tauh2so4tmp(1:ngrid)=0.
	328	DO l=1,nlayer
	329	DO ig=1,ngrid
	330	tauh2so4tmp(ig) = tauh2so4tmp(ig) + aerosol(ig,l,iaer)
	331	ENDDO
	332	ENDDO
	333	DO l=1,nlayer-1
	334	DO ig=1,ngrid
	335	aerosol(ig,l,iaer) = max(1E-20, &
[726]	336	1 &
[1308]	337	* pplev(ig,1) / pplev(ig,1) &
[726]	338	* aerosol(ig,l,iaer) &
	339	/ tauh2so4tmp(ig))
	340
	341	ENDDO
[728]	342	ENDDO
[726]	343
	344	! 1/700. is assuming a "sulfurtau" of 1
	345	! Sulfur aerosol routine to be improved.
	346	! aerosol0 = &
	347	! ( 0.75 * QREFvis3d(ig,l,iaer) / &
	348	! ( rho_h2so4 * reffrad(ig,l,iaer) ) ) * &
	349	! ( pq(ig,l,i_h2so4) + 1.E-9 ) * &
	350	! ( pplev(ig,l) - pplev(ig,l+1) ) / g
	351	! aerosol0 = max(aerosol0,1.e-9)
	352	! aerosol0 = min(aerosol0,L_TAUMAX)
	353	! aerosol(ig,l,iaer) = aerosol0
	354
	355	! ENDDO
	356	! ENDDO
	357	end if
	358
	359
[1026]	360	! ---------------------------------------------------------
	361	!==================================================================
	362	! Two-layer aerosols (unknown composition)
	363	! S. Guerlet (2013)
	364	!==================================================================
[726]	365
[1026]	366	if (iaero_back2lay .ne.0) then
	367	iaer=iaero_back2lay
	368	! 1. Initialization
[1308]	369	aerosol(1:ngrid,1:nlayer,iaer)=0.0
[1026]	370	! 2. Opacity calculation
	371	DO ig=1,ngrid
	372	DO l=1,nlayer-1
	373	aerosol(ig,l,iaer) = ( pplev(ig,l) - pplev(ig,l+1) )
	374	!! 1. below tropospheric layer: no aerosols
	375	IF (pplev(ig,l) .gt. pres_bottom_tropo) THEN
	376	aerosol(ig,l,iaer) = 0.*aerosol(ig,l,iaer)
	377	!! 2. tropo layer
	378	ELSEIF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN
	379	aerosol(ig,l,iaer) = obs_tau_col_tropo*aerosol(ig,l,iaer)
	380	!! 3. linear transition
	381	ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN
	382	expfactor=log(obs_tau_col_strato/obs_tau_col_tropo)/log(pres_bottom_strato/pres_top_tropo)
[1132]	383	aerosol(ig,l,iaer)= obs_tau_col_tropo((pplev(ig,l)/pres_top_tropo)expfactor)aerosol(ig,l,iaer)/1.5
[1026]	384	!! 4. strato layer
	385	ELSEIF (pplev(ig,l) .le. pres_bottom_strato .and. pplev(ig,l) .gt. pres_top_strato) THEN
	386	aerosol(ig,l,iaer)= obs_tau_col_strato*aerosol(ig,l,iaer)
	387	!! 5. above strato layer: no aerosols
	388	ELSEIF (pplev(ig,l) .lt. pres_top_strato) THEN
	389	aerosol(ig,l,iaer) = 0.*aerosol(ig,l,iaer)
	390	ENDIF
	391	ENDDO
	392	ENDDO
	393
	394	! 3. Re-normalize to observed total column
	395	tau_col(:)=0.0
	396	DO l=1,nlayer
	397	DO ig=1,ngrid
	398	tau_col(ig) = tau_col(ig) &
	399	+ aerosol(ig,l,iaer)/(obs_tau_col_tropo+obs_tau_col_strato)
	400	ENDDO
	401	ENDDO
	402
	403	DO ig=1,ngrid
	404	DO l=1,nlayer-1
	405	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig)
	406	ENDDO
	407	ENDDO
	408
	409
	410	end if ! if Two-layer aerosols
	411
	412
[135]	413	! --------------------------------------------------------------------------
	414	! Column integrated visible optical depth in each point (used for diagnostic)
	415
[253]	416	tau_col(:)=0.0
[135]	417	do iaer = 1, naerkind
	418	do l=1,nlayer
	419	do ig=1,ngrid
	420	tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
	421	end do
	422	end do
	423	end do
	424
[787]	425	do ig=1,ngrid
[253]	426	do l=1,nlayer
	427	do iaer = 1, naerkind
	428	if(aerosol(ig,l,iaer).gt.1.e3)then
	429	print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
	430	print*,'at ig=',ig,', l=',l,', iaer=',iaer
	431	print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
	432	print*,'reffrad=',reffrad(ig,l,iaer)
	433	endif
	434	end do
	435	end do
	436	end do
	437
[787]	438	do ig=1,ngrid
[253]	439	if(tau_col(ig).gt.1.e3)then
	440	print*,'WARNING: tau_col=',tau_col(ig)
	441	print*,'at ig=',ig
	442	print*,'aerosol=',aerosol(ig,:,:)
	443	print*,'QREFvis3d=',QREFvis3d(ig,:,:)
	444	print*,'reffrad=',reffrad(ig,:,:)
	445	endif
	446	end do
[135]	447	return
	448	end subroutine aeropacity
	449

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

Download in other formats:

Original Format