Context Navigation

mo_simple_plumes.F90 @ 5049

Last change on this file since 5049 was 4164, checked in by oboucher, 2 years ago
add a SAVE argument to variables read in simple aerosol plume model
File size: 22.4 KB

Rev	Line
[3274]	1	!>
	2	!!
	3	!! @brief Module MO_SIMPLE_PLUMES: provides anthropogenic aerosol optical properties as a function of lat, lon
	4	!! height, time, and wavelength
	5	!!
	6	!! @remarks
	7	!!
	8	!! @author Bjorn Stevens, Stephanie Fiedler and Karsten Peters MPI-Met, Hamburg (v1 release 2016-11-10)
	9	!!
	10	!! @change-log:
	11	!! - 2016-12-05: beta release (BS, SF and KP, MPI-Met)
	12	!! - 2016-09-28: revised representation of Twomey effect (SF, MPI-Met)
	13	!! - 2015-09-28: bug fixes (SF, MPI-Met)
	14	!! - 2016-10-12: revised maximum longitudinal extent of European plume (KP, SF, MPI-Met)
	15	!! $ID: n/a$
	16	!!
	17	!! @par Origin
	18	!! Based on code originally developed at the MPI-Met by Karsten Peters, Bjorn Stevens, Stephanie Fiedler
	19	!! and Stefan Kinne with input from Thorsten Mauritsen and Robert Pincus
	20	!!
	21	!! @par Copyright
	22	!!
	23	!
	24	MODULE MO_SIMPLE_PLUMES
	25
	26	USE netcdf
	27
	28	IMPLICIT NONE
	29
	30	INTEGER, PARAMETER :: &
	31	nplumes = 9 ,& !< Number of plumes
	32	nfeatures = 2 ,& !< Number of features per plume
	33	ntimes = 52 ,& !< Number of times resolved per year (52 => weekly resolution)
	34	nyears = 251 !< Number of years of available forcing
	35
	36	LOGICAL, SAVE :: &
	37	sp_initialized = .FALSE. !< parameter determining whether input needs to be read
[3295]	38	!$OMP THREADPRIVATE(sp_initialized)
[3274]	39
[4164]	40	REAL, SAVE :: &
[3274]	41	plume_lat (nplumes) ,& !< latitude of plume center (AOD maximum)
	42	plume_lon (nplumes) ,& !< longitude of plume center (AOD maximum)
	43	beta_a (nplumes) ,& !< parameter a for beta function vertical profile
	44	beta_b (nplumes) ,& !< parameter b for beta function vertical profile
	45	aod_spmx (nplumes) ,& !< anthropogenic AOD maximum at 550 for plumes
	46	aod_fmbg (nplumes) ,& !< anthropogenic AOD at 550 for fine-mode natural background (idealized to mimic Twomey effect)
	47	asy550 (nplumes) ,& !< asymmetry parameter at 550nm for plume
	48	ssa550 (nplumes) ,& !< single scattering albedo at 550nm for plume
	49	angstrom (nplumes) ,& !< Angstrom parameter for plume
	50	sig_lon_E (nfeatures,nplumes) ,& !< Eastward extent of plume feature
	51	sig_lon_W (nfeatures,nplumes) ,& !< Westward extent of plume feature
	52	sig_lat_E (nfeatures,nplumes) ,& !< Southward extent of plume feature
	53	sig_lat_W (nfeatures,nplumes) ,& !< Northward extent of plume feature
	54	theta (nfeatures,nplumes) ,& !< Rotation angle of plume feature
	55	ftr_weight (nfeatures,nplumes) ,& !< Feature weights
	56	year_weight (nyears,nplumes) ,& !< Yearly weight for plume
	57	ann_cycle (nfeatures,ntimes,nplumes) !< annual cycle for plume feature
[3295]	58	!$OMP THREADPRIVATE(plume_lat,plume_lon,beta_a,beta_b,aod_spmx,aod_fmbg,asy550,ssa550,angstrom)
	59	!$OMP THREADPRIVATE(sig_lon_E,sig_lon_W,sig_lat_E,sig_lat_W,theta,ftr_weight,year_weight,ann_cycle)
[3274]	60
[4164]	61	REAL :: &
	62	time_weight (nfeatures,nplumes) ,& !< Time weights
	63	time_weight_bg (nfeatures,nplumes) !< as time_weight but for natural background in Twomey effect
	64
[3274]	65	PUBLIC sp_aop_profile
	66
	67	CONTAINS
	68	!
	69	! ------------------------------------------------------------------------------------------------------------------------
	70	! SP_SETUP: This subroutine should be called at initialization to read the netcdf data that describes the simple plume
	71	! climatology. The information needs to be either read by each processor or distributed to processors.
	72	!
	73	SUBROUTINE sp_setup
	74	!
[3295]	75	USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
	76	USE mod_phys_lmdz_omp_data, ONLY: is_omp_root
[3297]	77	USE mod_phys_lmdz_transfert_para, ONLY: bcast
[3295]	78	!
[3274]	79	! ----------
	80	!
	81	INTEGER :: iret, ncid, DimID, VarID, xdmy
[3531]	82	CHARACTER (len = 50) :: modname = 'mo_simple_plumes.sp_setup'
	83	CHARACTER (len = 80) :: abort_message
[3274]	84	!
	85	! ----------
[3295]	86	!--only one processor reads the input data
	87	IF (is_mpi_root.AND.is_omp_root) THEN
	88	!
	89	iret = nf90_open("MACv2.0-SP_v1.nc", NF90_NOWRITE, ncid)
[3531]	90	IF (iret /= NF90_NOERR) THEN
	91	abort_message='NetCDF File not opened'
	92	CALL abort_physic(modname,abort_message,1)
	93	ENDIF
[3295]	94	!
	95	! read dimensions and make sure file conforms to expected size
	96	!
	97	iret = nf90_inq_dimid(ncid, "plume_number" , DimId)
	98	iret = nf90_inquire_dimension(ncid, DimId, len = xdmy)
[3531]	99	IF (xdmy /= nplumes) THEN
	100	abort_message='NetCDF improperly dimensioned -- plume_number'
	101	CALL abort_physic(modname,abort_message,1)
	102	ENDIF
[3295]	103	!
	104	iret = nf90_inq_dimid(ncid, "plume_feature", DimId)
	105	iret = nf90_inquire_dimension(ncid, DimId, len = xdmy)
[3531]	106	IF (xdmy /= nfeatures) THEN
	107	abort_message='NetCDF improperly dimensioned -- plume_feature'
	108	CALL abort_physic(modname,abort_message,1)
	109	ENDIF
[3295]	110	!
	111	iret = nf90_inq_dimid(ncid, "year_fr" , DimId)
	112	iret = nf90_inquire_dimension(ncid, DimID, len = xdmy)
[3531]	113	IF (xdmy /= ntimes) THEN
	114	abort_message='NetCDF improperly dimensioned -- year_fr'
	115	CALL abort_physic(modname,abort_message,1)
	116	ENDIF
[3295]	117	!
	118	iret = nf90_inq_dimid(ncid, "years" , DimId)
	119	iret = nf90_inquire_dimension(ncid, DimID, len = xdmy)
[3531]	120	IF (xdmy /= nyears) THEN
	121	abort_message='NetCDF improperly dimensioned -- years'
	122	CALL abort_physic(modname,abort_message,1)
	123	ENDIF
[3295]	124	!
	125	! read variables that define the simple plume climatology
	126	!
	127	iret = nf90_inq_varid(ncid, "plume_lat", VarId)
	128	iret = nf90_get_var(ncid, VarID, plume_lat(:), start=(/1/),count=(/nplumes/))
[3531]	129	IF (iret /= NF90_NOERR) THEN
	130	abort_message='NetCDF Error reading plume_lat'
	131	CALL abort_physic(modname,abort_message,1)
	132	ENDIF
[4164]	133	!
[3295]	134	iret = nf90_inq_varid(ncid, "plume_lon", VarId)
	135	iret = nf90_get_var(ncid, VarID, plume_lon(:), start=(/1/),count=(/nplumes/))
[3531]	136	IF (iret /= NF90_NOERR) THEN
	137	abort_message='NetCDF Error reading plume_lon'
	138	CALL abort_physic(modname,abort_message,1)
	139	ENDIF
[4164]	140	!
[3295]	141	iret = nf90_inq_varid(ncid, "beta_a" , VarId)
	142	iret = nf90_get_var(ncid, VarID, beta_a(:) , start=(/1/),count=(/nplumes/))
[3531]	143	IF (iret /= NF90_NOERR) THEN
	144	abort_message='NetCDF Error reading beta_a'
	145	CALL abort_physic(modname,abort_message,1)
	146	ENDIF
[4164]	147	!
[3295]	148	iret = nf90_inq_varid(ncid, "beta_b" , VarId)
	149	iret = nf90_get_var(ncid, VarID, beta_b(:) , start=(/1/),count=(/nplumes/))
[3531]	150	IF (iret /= NF90_NOERR) THEN
	151	abort_message='NetCDF Error reading beta_b'
	152	CALL abort_physic(modname,abort_message,1)
	153	ENDIF
[4164]	154	!
[3295]	155	iret = nf90_inq_varid(ncid, "aod_spmx" , VarId)
	156	iret = nf90_get_var(ncid, VarID, aod_spmx(:) , start=(/1/),count=(/nplumes/))
[3531]	157	IF (iret /= NF90_NOERR) THEN
	158	abort_message='NetCDF Error reading aod_spmx'
	159	CALL abort_physic(modname,abort_message,1)
	160	ENDIF
[4164]	161	!
[3295]	162	iret = nf90_inq_varid(ncid, "aod_fmbg" , VarId)
	163	iret = nf90_get_var(ncid, VarID, aod_fmbg(:) , start=(/1/),count=(/nplumes/))
[3531]	164	IF (iret /= NF90_NOERR) THEN
	165	abort_message='NetCDF Error reading aod_fmbg'
	166	CALL abort_physic(modname,abort_message,1)
	167	ENDIF
[4164]	168	!
[3295]	169	iret = nf90_inq_varid(ncid, "ssa550" , VarId)
	170	iret = nf90_get_var(ncid, VarID, ssa550(:) , start=(/1/),count=(/nplumes/))
[3531]	171	IF (iret /= NF90_NOERR) THEN
	172	abort_message='NetCDF Error reading ssa550'
	173	CALL abort_physic(modname,abort_message,1)
	174	ENDIF
[4164]	175	!
[3295]	176	iret = nf90_inq_varid(ncid, "asy550" , VarId)
	177	iret = nf90_get_var(ncid, VarID, asy550(:) , start=(/1/),count=(/nplumes/))
[3531]	178	IF (iret /= NF90_NOERR) THEN
	179	abort_message='NetCDF Error reading asy550'
	180	CALL abort_physic(modname,abort_message,1)
	181	ENDIF
[4164]	182	!
[3295]	183	iret = nf90_inq_varid(ncid, "angstrom" , VarId)
	184	iret = nf90_get_var(ncid, VarID, angstrom(:), start=(/1/),count=(/nplumes/))
[3531]	185	IF (iret /= NF90_NOERR) THEN
	186	abort_message='NetCDF Error reading angstrom'
	187	CALL abort_physic(modname,abort_message,1)
	188	ENDIF
[3295]	189	!
	190	iret = nf90_inq_varid(ncid, "sig_lat_W" , VarId)
	191	iret = nf90_get_var(ncid, VarID, sig_lat_W(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
[3531]	192	IF (iret /= NF90_NOERR) THEN
	193	abort_message='NetCDF Error reading sig_lat_W'
	194	CALL abort_physic(modname,abort_message,1)
	195	ENDIF
[4164]	196	!
[3295]	197	iret = nf90_inq_varid(ncid, "sig_lat_E" , VarId)
	198	iret = nf90_get_var(ncid, VarID, sig_lat_E(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
[3531]	199	IF (iret /= NF90_NOERR) THEN
	200	abort_message='NetCDF Error reading sig_lat_E'
	201	CALL abort_physic(modname,abort_message,1)
	202	ENDIF
[4164]	203	!
[3295]	204	iret = nf90_inq_varid(ncid, "sig_lon_E" , VarId)
	205	iret = nf90_get_var(ncid, VarID, sig_lon_E(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
[3531]	206	IF (iret /= NF90_NOERR) THEN
	207	abort_message='NetCDF Error reading sig_lon_E'
	208	CALL abort_physic(modname,abort_message,1)
	209	ENDIF
[4164]	210	!
[3295]	211	iret = nf90_inq_varid(ncid, "sig_lon_W" , VarId)
	212	iret = nf90_get_var(ncid, VarID, sig_lon_W(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
[3531]	213	IF (iret /= NF90_NOERR) THEN
	214	abort_message='NetCDF Error reading sig_lon_W'
	215	CALL abort_physic(modname,abort_message,1)
	216	ENDIF
[4164]	217	!
[3295]	218	iret = nf90_inq_varid(ncid, "theta" , VarId)
	219	iret = nf90_get_var(ncid, VarID, theta(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
[3531]	220	IF (iret /= NF90_NOERR) THEN
	221	abort_message='NetCDF Error reading theta'
	222	CALL abort_physic(modname,abort_message,1)
	223	ENDIF
[4164]	224	!
[3295]	225	iret = nf90_inq_varid(ncid, "ftr_weight" , VarId)
	226	iret = nf90_get_var(ncid, VarID, ftr_weight(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
[3531]	227	IF (iret /= NF90_NOERR) THEN
	228	abort_message='NetCDF Error reading plume_lat'
	229	CALL abort_physic(modname,abort_message,1)
	230	ENDIF
[4164]	231	!
[3295]	232	iret = nf90_inq_varid(ncid, "year_weight" , VarId)
	233	iret = nf90_get_var(ncid, VarID, year_weight(:,:) , start=(/1,1/),count=(/nyears,nplumes /))
[3531]	234	IF (iret /= NF90_NOERR) THEN
	235	abort_message='NetCDF Error reading year_weight'
	236	CALL abort_physic(modname,abort_message,1)
	237	ENDIF
[4164]	238	!
[3295]	239	iret = nf90_inq_varid(ncid, "ann_cycle" , VarId)
	240	iret = nf90_get_var(ncid, VarID, ann_cycle(:,:,:) , start=(/1,1,1/),count=(/nfeatures,ntimes,nplumes/))
[3531]	241	IF (iret /= NF90_NOERR) THEN
	242	abort_message='NetCDF Error reading ann_cycle'
	243	CALL abort_physic(modname,abort_message,1)
	244	ENDIF
[3295]	245	!
	246	iret = nf90_close(ncid)
	247	!
	248	ENDIF !--root processor
[3274]	249	!
[3295]	250	CALL bcast(plume_lat)
	251	CALL bcast(plume_lon)
	252	CALL bcast(beta_a)
	253	CALL bcast(beta_b)
	254	CALL bcast(aod_spmx)
	255	CALL bcast(aod_fmbg)
	256	CALL bcast(asy550)
	257	CALL bcast(ssa550)
	258	CALL bcast(angstrom)
	259	CALL bcast(sig_lon_E)
	260	CALL bcast(sig_lon_W)
	261	CALL bcast(sig_lat_E)
	262	CALL bcast(sig_lat_W)
	263	CALL bcast(theta)
	264	CALL bcast(ftr_weight)
	265	CALL bcast(year_weight)
	266	CALL bcast(ann_cycle)
[3274]	267	!
[3295]	268	sp_initialized = .TRUE.
[3274]	269	!
[3295]	270	RETURN
[3274]	271	!
	272	END SUBROUTINE sp_setup
	273	!
	274	! ------------------------------------------------------------------------------------------------------------------------
	275	! SET_TIME_WEIGHT: The simple plume model assumes that meteorology constrains plume shape and that only source strength
	276	! influences the amplitude of a plume associated with a given source region. This routine retrieves the temporal weights
	277	! for the plumes. Each plume feature has its own temporal weights which varies yearly. The annual cycle is indexed by
	278	! week in the year and superimposed on the yearly mean value of the weight.
	279	!
	280	SUBROUTINE set_time_weight(year_fr)
	281	!
	282	! ----------
	283	!
	284	REAL, INTENT(IN) :: &
	285	year_fr !< Fractional Year (1850.0 - 2100.99)
	286
	287	INTEGER :: &
	288	iyear ,& !< Integer year values between 1 and 156 (1850-2100)
	289	iweek ,& !< Integer index (between 1 and ntimes); for ntimes=52 this corresponds to weeks (roughly)
	290	iplume ! plume number
	291	!
	292	! ----------
	293	!
	294	iyear = FLOOR(year_fr) - 1849
	295	iweek = FLOOR((year_fr - FLOOR(year_fr)) * ntimes) + 1
	296
	297	IF ((iweek > ntimes) .OR. (iweek < 1) .OR. (iyear > nyears) .OR. (iyear < 1)) THEN
[4103]	298	CALL abort_physic('set_time_weight','Time out of bounds',1)
[3274]	299	ENDIF
	300
	301	DO iplume=1,nplumes
	302	time_weight(1,iplume) = year_weight(iyear,iplume) * ann_cycle(1,iweek,iplume)
	303	time_weight(2,iplume) = year_weight(iyear,iplume) * ann_cycle(2,iweek,iplume)
	304	time_weight_bg(1,iplume) = ann_cycle(1,iweek,iplume)
	305	time_weight_bg(2,iplume) = ann_cycle(2,iweek,iplume)
[4164]	306	ENDDO
[3274]	307
	308	RETURN
	309	END SUBROUTINE set_time_weight
	310	!
	311	! ------------------------------------------------------------------------------------------------------------------------
	312	! SP_AOP_PROFILE: This subroutine calculates the simple plume aerosol and cloud active optical properties based on the
	313	! the simple plume fit to the MPI Aerosol Climatology (Version 2). It sums over nplumes to provide a profile of aerosol
	314	! optical properties on a host models vertical grid.
	315	!
	316	SUBROUTINE sp_aop_profile ( &
	317	nlevels ,ncol ,lambda ,oro ,lon ,lat , &
	318	year_fr ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , &
	319	asy_prof )
	320	!
	321	! ----------
	322	!
	323	INTEGER, INTENT(IN) :: &
	324	nlevels, & !< number of levels
	325	ncol !< number of columns
	326
	327	REAL, INTENT(IN) :: &
	328	lambda, & !< wavelength
	329	year_fr, & !< Fractional Year (1903.0 is the 0Z on the first of January 1903, Gregorian)
	330	oro(ncol), & !< orographic height (m)
	331	lon(ncol), & !< longitude
	332	lat(ncol), & !< latitude
	333	z (ncol,nlevels), & !< height above sea-level (m)
	334	dz(ncol,nlevels) !< level thickness (difference between half levels) (m)
	335
	336	REAL, INTENT(OUT) :: &
	337	dNovrN(ncol) , & !< anthropogenic increase in cloud drop number concentration (factor)
	338	aod_prof(ncol,nlevels) , & !< profile of aerosol optical depth
	339	ssa_prof(ncol,nlevels) , & !< profile of single scattering albedo
	340	asy_prof(ncol,nlevels) !< profile of asymmetry parameter
	341
	342	INTEGER :: iplume, icol, k
	343
	344	REAL :: &
	345	eta(ncol,nlevels), & !< normalized height (by 15 km)
	346	z_beta(ncol,nlevels), & !< profile for scaling column optical depth
	347	prof(ncol,nlevels), & !< scaled profile (by beta function)
	348	beta_sum(ncol), & !< vertical sum of beta function
	349	ssa(ncol), & !< single scattering albedo
	350	asy(ncol), & !< asymmetry parameter
	351	cw_an(ncol), & !< column weight for simple plume (anthropogenic) AOD at 550 nm
	352	cw_bg(ncol), & !< column weight for fine-mode natural background AOD at 550 nm
	353	caod_sp(ncol), & !< column simple plume anthropogenic AOD at 550 nm
	354	caod_bg(ncol), & !< column fine-mode natural background AOD at 550 nm
	355	a_plume1, & !< gaussian longitude factor for feature 1
	356	a_plume2, & !< gaussian longitude factor for feature 2
	357	b_plume1, & !< gaussian latitude factor for feature 1
	358	b_plume2, & !< gaussian latitude factor for feature 2
	359	delta_lat, & !< latitude offset
	360	delta_lon, & !< longitude offset
	361	delta_lon_t, & !< threshold for maximum longitudinal plume extent used in transition from 360 to 0 degrees
	362	lon1, & !< rotated longitude for feature 1
	363	lat1, & !< rotated latitude for feature 2
	364	lon2, & !< rotated longitude for feature 1
	365	lat2, & !< rotated latitude for feature 2
	366	f1, & !< contribution from feature 1
	367	f2, & !< contribution from feature 2
	368	f3, & !< contribution from feature 1 in natural background of Twomey effect
	369	f4, & !< contribution from feature 2 in natural background of Twomey effect
	370	aod_550, & !< aerosol optical depth at 550nm
	371	aod_lmd, & !< aerosol optical depth at input wavelength
	372	lfactor !< factor to compute wavelength dependence of optical properties
	373	!
	374	! ----------
	375	!
	376	! initialize input data (by calling setup at first instance)
	377	!
	378	IF (.NOT.sp_initialized) CALL sp_setup
	379	!
	380	! get time weights
	381	!
	382	CALL set_time_weight(year_fr)
	383	!
	384	! initialize variables, including output
	385	!
	386	DO k=1,nlevels
	387	DO icol=1,ncol
	388	aod_prof(icol,k) = 0.0
	389	ssa_prof(icol,k) = 0.0
	390	asy_prof(icol,k) = 0.0
	391	z_beta(icol,k) = MERGE(1.0, 0.0, z(icol,k) >= oro(icol))
	392	eta(icol,k) = MAX(0.0,MIN(1.0,z(icol,k)/15000.))
[4164]	393	ENDDO
	394	ENDDO
[3274]	395	DO icol=1,ncol
	396	dNovrN(icol) = 1.0
	397	caod_sp(icol) = 0.0
	398	caod_bg(icol) = 0.02
[4164]	399	ENDDO
[3274]	400	!
	401	! sum contribution from plumes to construct composite profiles of aerosol optical properties
	402	!
	403	DO iplume=1,nplumes
	404	!
	405	! calculate vertical distribution function from parameters of beta distribution
	406	!
	407	DO icol=1,ncol
	408	beta_sum(icol) = 0.
[4164]	409	ENDDO
[3274]	410	DO k=1,nlevels
	411	DO icol=1,ncol
	412	prof(icol,k) = (eta(icol,k)*(beta_a(iplume)-1.) (1.-eta(icol,k))*(beta_b(iplume)-1.)) dz(icol,k)
	413	beta_sum(icol) = beta_sum(icol) + prof(icol,k)
[4164]	414	ENDDO
	415	ENDDO
[3274]	416	DO k=1,nlevels
	417	DO icol=1,ncol
	418	prof(icol,k) = ( prof(icol,k) / beta_sum(icol) ) * z_beta(icol,k)
[4164]	419	ENDDO
	420	ENDDO
[3274]	421	!
	422	! calculate plume weights
	423	!
	424	DO icol=1,ncol
	425	!
	426	! get plume-center relative spatial parameters for specifying amplitude of plume at given lat and lon
	427	!
	428	delta_lat = lat(icol) - plume_lat(iplume)
	429	delta_lon = lon(icol) - plume_lon(iplume)
	430	delta_lon_t = MERGE (260., 180., iplume == 1)
	431	delta_lon = MERGE ( delta_lon-SIGN(360.,delta_lon) , delta_lon , ABS(delta_lon) > delta_lon_t)
	432
	433	a_plume1 = 0.5 / (MERGE(sig_lon_E(1,iplume), sig_lon_W(1,iplume), delta_lon > 0)**2)
	434	b_plume1 = 0.5 / (MERGE(sig_lat_E(1,iplume), sig_lat_W(1,iplume), delta_lon > 0)**2)
	435	a_plume2 = 0.5 / (MERGE(sig_lon_E(2,iplume), sig_lon_W(2,iplume), delta_lon > 0)**2)
	436	b_plume2 = 0.5 / (MERGE(sig_lat_E(2,iplume), sig_lat_W(2,iplume), delta_lon > 0)**2)
	437	!
	438	! adjust for a plume specific rotation which helps match plume state to climatology.
	439	!
	440	lon1 = COS(theta(1,iplume))(delta_lon) + SIN(theta(1,iplume))(delta_lat)
	441	lat1 = - SIN(theta(1,iplume))(delta_lon) + COS(theta(1,iplume))(delta_lat)
	442	lon2 = COS(theta(2,iplume))(delta_lon) + SIN(theta(2,iplume))(delta_lat)
	443	lat2 = - SIN(theta(2,iplume))(delta_lon) + COS(theta(2,iplume))(delta_lat)
	444	!
	445	! calculate contribution to plume from its different features, to get a column weight for the anthropogenic
	446	! (cw_an) and the fine-mode natural background aerosol (cw_bg)
	447	!
	448	f1 = time_weight(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)*2) + (b_plume1 ((lat1)**2))))
	449	f2 = time_weight(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)*2) + (b_plume2 ((lat2)**2))))
	450	f3 = time_weight_bg(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)*2) + (b_plume1 ((lat1)**2))))
	451	f4 = time_weight_bg(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)*2) + (b_plume2 ((lat2)**2))))
	452
	453	cw_an(icol) = f1 * aod_spmx(iplume) + f2 * aod_spmx(iplume)
	454	cw_bg(icol) = f3 * aod_fmbg(iplume) + f4 * aod_fmbg(iplume)
	455	!
	456	! calculate wavelength-dependent scattering properties
	457	!
	458	lfactor = MIN(1.0,700.0/lambda)
	459	ssa(icol) = (ssa550(iplume) * lfactor*4) / ((ssa550(iplume) lfactor*4) + ((1-ssa550(iplume)) lfactor))
	460	asy(icol) = asy550(iplume) * SQRT(lfactor)
[4164]	461	ENDDO
[3274]	462	!
	463	! distribute plume optical properties across its vertical profile weighting by optical depth and scaling for
	464	! wavelength using the angstrom parameter.
	465	!
	466	lfactor = EXP(-angstrom(iplume) * LOG(lambda/550.0))
	467	DO k=1,nlevels
	468	DO icol = 1,ncol
	469	aod_550 = prof(icol,k) * cw_an(icol)
	470	aod_lmd = aod_550 * lfactor
	471	caod_sp(icol) = caod_sp(icol) + aod_550
	472	caod_bg(icol) = caod_bg(icol) + prof(icol,k) * cw_bg(icol)
	473	asy_prof(icol,k) = asy_prof(icol,k) + aod_lmd * ssa(icol) * asy(icol)
	474	ssa_prof(icol,k) = ssa_prof(icol,k) + aod_lmd * ssa(icol)
	475	aod_prof(icol,k) = aod_prof(icol,k) + aod_lmd
[4164]	476	ENDDO
	477	ENDDO
	478	ENDDO
[3274]	479	!
	480	! complete optical depth weighting
	481	!
	482	DO k=1,nlevels
	483	DO icol = 1,ncol
	484	asy_prof(icol,k) = MERGE(asy_prof(icol,k)/ssa_prof(icol,k), 0.0, ssa_prof(icol,k) > TINY(1.))
	485	ssa_prof(icol,k) = MERGE(ssa_prof(icol,k)/aod_prof(icol,k), 1.0, aod_prof(icol,k) > TINY(1.))
[4164]	486	ENDDO
	487	ENDDO
[3274]	488	!
	489	! calculate effective radius normalization (divisor) factor
	490	!
	491	DO icol=1,ncol
	492	dNovrN(icol) = LOG((1000.0 * (caod_sp(icol) + caod_bg(icol))) + 1.0)/LOG((1000.0 * caod_bg(icol)) + 1.0)
[4164]	493	ENDDO
[3274]	494
	495	RETURN
	496	END SUBROUTINE sp_aop_profile
	497
	498	END MODULE MO_SIMPLE_PLUMES

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

Download in other formats:

Original Format