Context Navigation

mo_simple_plumes.F90 @ 5077

Last change on this file since 5077 was 5075, checked in by abarral, 5 months ago

[continued & end] replace netcdf by lmdz_netcdf.F90 wrapper
"use netcdf" is now only used in lmdz_netcdf.F90 (except ecrad and obsolete/)
<include "netcdf.inc"> is now likewise only used in lmdz_netcdf.F90.

systematically specify explicitely <USE lmdz_netcdf, ONLY:> (probably left some missing, to correct later on)

Further replacement of nf_put_* by nf90_put_* (same for _get_)

[minor] replace deprecated boolean operators along the way

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

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

Download in other formats:

Original Format