Context Navigation

mo_simple_plumes.F90 @ 3654

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

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

Download in other formats:

Original Format