Context Navigation

mo_simple_plumes.F90 @ 5171

Last change on this file since 5171 was 5084, checked in by Laurent Fairhead, 12 months ago
Reverting to r4065. Updating fortran standard broke too much stuff. Will do it by smaller chunks AB, LF
File size: 22.4 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 :: &
37	sp_initialized = .FALSE. !< parameter determining whether input needs to be read
38	!$OMP THREADPRIVATE(sp_initialized)
39
40	REAL, SAVE :: &
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
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)
60
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
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	!
75	USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
76	USE mod_phys_lmdz_omp_data, ONLY: is_omp_root
77	USE mod_phys_lmdz_transfert_para, ONLY: bcast
78	!
79	! ----------
80	!
81	INTEGER :: iret, ncid, DimID, VarID, xdmy
82	CHARACTER (len = 50) :: modname = 'mo_simple_plumes.sp_setup'
83	CHARACTER (len = 80) :: abort_message
84	!
85	! ----------
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)
90	IF (iret /= NF90_NOERR) THEN
91	abort_message='NetCDF File not opened'
92	CALL abort_physic(modname,abort_message,1)
93	ENDIF
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)
99	IF (xdmy /= nplumes) THEN
100	abort_message='NetCDF improperly dimensioned -- plume_number'
101	CALL abort_physic(modname,abort_message,1)
102	ENDIF
103	!
104	iret = nf90_inq_dimid(ncid, "plume_feature", DimId)
105	iret = nf90_inquire_dimension(ncid, DimId, len = xdmy)
106	IF (xdmy /= nfeatures) THEN
107	abort_message='NetCDF improperly dimensioned -- plume_feature'
108	CALL abort_physic(modname,abort_message,1)
109	ENDIF
110	!
111	iret = nf90_inq_dimid(ncid, "year_fr" , DimId)
112	iret = nf90_inquire_dimension(ncid, DimID, len = xdmy)
113	IF (xdmy /= ntimes) THEN
114	abort_message='NetCDF improperly dimensioned -- year_fr'
115	CALL abort_physic(modname,abort_message,1)
116	ENDIF
117	!
118	iret = nf90_inq_dimid(ncid, "years" , DimId)
119	iret = nf90_inquire_dimension(ncid, DimID, len = xdmy)
120	IF (xdmy /= nyears) THEN
121	abort_message='NetCDF improperly dimensioned -- years'
122	CALL abort_physic(modname,abort_message,1)
123	ENDIF
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/))
129	IF (iret /= NF90_NOERR) THEN
130	abort_message='NetCDF Error reading plume_lat'
131	CALL abort_physic(modname,abort_message,1)
132	ENDIF
133	!
134	iret = nf90_inq_varid(ncid, "plume_lon", VarId)
135	iret = nf90_get_var(ncid, VarID, plume_lon(:), start=(/1/),count=(/nplumes/))
136	IF (iret /= NF90_NOERR) THEN
137	abort_message='NetCDF Error reading plume_lon'
138	CALL abort_physic(modname,abort_message,1)
139	ENDIF
140	!
141	iret = nf90_inq_varid(ncid, "beta_a" , VarId)
142	iret = nf90_get_var(ncid, VarID, beta_a(:) , start=(/1/),count=(/nplumes/))
143	IF (iret /= NF90_NOERR) THEN
144	abort_message='NetCDF Error reading beta_a'
145	CALL abort_physic(modname,abort_message,1)
146	ENDIF
147	!
148	iret = nf90_inq_varid(ncid, "beta_b" , VarId)
149	iret = nf90_get_var(ncid, VarID, beta_b(:) , start=(/1/),count=(/nplumes/))
150	IF (iret /= NF90_NOERR) THEN
151	abort_message='NetCDF Error reading beta_b'
152	CALL abort_physic(modname,abort_message,1)
153	ENDIF
154	!
155	iret = nf90_inq_varid(ncid, "aod_spmx" , VarId)
156	iret = nf90_get_var(ncid, VarID, aod_spmx(:) , start=(/1/),count=(/nplumes/))
157	IF (iret /= NF90_NOERR) THEN
158	abort_message='NetCDF Error reading aod_spmx'
159	CALL abort_physic(modname,abort_message,1)
160	ENDIF
161	!
162	iret = nf90_inq_varid(ncid, "aod_fmbg" , VarId)
163	iret = nf90_get_var(ncid, VarID, aod_fmbg(:) , start=(/1/),count=(/nplumes/))
164	IF (iret /= NF90_NOERR) THEN
165	abort_message='NetCDF Error reading aod_fmbg'
166	CALL abort_physic(modname,abort_message,1)
167	ENDIF
168	!
169	iret = nf90_inq_varid(ncid, "ssa550" , VarId)
170	iret = nf90_get_var(ncid, VarID, ssa550(:) , start=(/1/),count=(/nplumes/))
171	IF (iret /= NF90_NOERR) THEN
172	abort_message='NetCDF Error reading ssa550'
173	CALL abort_physic(modname,abort_message,1)
174	ENDIF
175	!
176	iret = nf90_inq_varid(ncid, "asy550" , VarId)
177	iret = nf90_get_var(ncid, VarID, asy550(:) , start=(/1/),count=(/nplumes/))
178	IF (iret /= NF90_NOERR) THEN
179	abort_message='NetCDF Error reading asy550'
180	CALL abort_physic(modname,abort_message,1)
181	ENDIF
182	!
183	iret = nf90_inq_varid(ncid, "angstrom" , VarId)
184	iret = nf90_get_var(ncid, VarID, angstrom(:), start=(/1/),count=(/nplumes/))
185	IF (iret /= NF90_NOERR) THEN
186	abort_message='NetCDF Error reading angstrom'
187	CALL abort_physic(modname,abort_message,1)
188	ENDIF
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/))
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
196	!
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/))
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
203	!
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/))
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
210	!
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/))
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
217	!
218	iret = nf90_inq_varid(ncid, "theta" , VarId)
219	iret = nf90_get_var(ncid, VarID, theta(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
220	IF (iret /= NF90_NOERR) THEN
221	abort_message='NetCDF Error reading theta'
222	CALL abort_physic(modname,abort_message,1)
223	ENDIF
224	!
225	iret = nf90_inq_varid(ncid, "ftr_weight" , VarId)
226	iret = nf90_get_var(ncid, VarID, ftr_weight(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
227	IF (iret /= NF90_NOERR) THEN
228	abort_message='NetCDF Error reading plume_lat'
229	CALL abort_physic(modname,abort_message,1)
230	ENDIF
231	!
232	iret = nf90_inq_varid(ncid, "year_weight" , VarId)
233	iret = nf90_get_var(ncid, VarID, year_weight(:,:) , start=(/1,1/),count=(/nyears,nplumes /))
234	IF (iret /= NF90_NOERR) THEN
235	abort_message='NetCDF Error reading year_weight'
236	CALL abort_physic(modname,abort_message,1)
237	ENDIF
238	!
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/))
241	IF (iret /= NF90_NOERR) THEN
242	abort_message='NetCDF Error reading ann_cycle'
243	CALL abort_physic(modname,abort_message,1)
244	ENDIF
245	!
246	iret = nf90_close(ncid)
247	!
248	ENDIF !--root processor
249	!
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)
267	!
268	sp_initialized = .TRUE.
269	!
270	RETURN
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
298	CALL abort_physic('set_time_weight','Time out of bounds',1)
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)
306	ENDDO
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.))
393	ENDDO
394	ENDDO
395	DO icol=1,ncol
396	dNovrN(icol) = 1.0
397	caod_sp(icol) = 0.0
398	caod_bg(icol) = 0.02
399	ENDDO
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.
409	ENDDO
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)
414	ENDDO
415	ENDDO
416	DO k=1,nlevels
417	DO icol=1,ncol
418	prof(icol,k) = ( prof(icol,k) / beta_sum(icol) ) * z_beta(icol,k)
419	ENDDO
420	ENDDO
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)
461	ENDDO
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
476	ENDDO
477	ENDDO
478	ENDDO
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.))
486	ENDDO
487	ENDDO
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)
493	ENDDO
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