Context Navigation

mo_simple_plumes.f90 @ 5408

Last change on this file since 5408 was 5395, checked in by evignon, 9 days ago
suite precedent commit
File size: 22.5 KB

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

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

Download in other formats:

Original Format