Context Navigation

mo_simple_plumes.f90 @ 5394

Last change on this file since 5394 was 5268, checked in by abarral, 7 weeks ago
.f90 <-> .F90 depending on cpp key use
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