Context Navigation

mo_simple_plumes.F90 @ 3991

Last change on this file since 3991 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

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 :: &
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
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)
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	! ----------
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)
89	IF (iret /= NF90_NOERR) THEN
90	abort_message='NetCDF File not opened'
91	CALL abort_physic(modname,abort_message,1)
92	ENDIF
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)
98	IF (xdmy /= nplumes) THEN
99	abort_message='NetCDF improperly dimensioned -- plume_number'
100	CALL abort_physic(modname,abort_message,1)
101	ENDIF
102	!
103	iret = nf90_inq_dimid(ncid, "plume_feature", DimId)
104	iret = nf90_inquire_dimension(ncid, DimId, len = xdmy)
105	IF (xdmy /= nfeatures) THEN
106	abort_message='NetCDF improperly dimensioned -- plume_feature'
107	CALL abort_physic(modname,abort_message,1)
108	ENDIF
109	!
110	iret = nf90_inq_dimid(ncid, "year_fr" , DimId)
111	iret = nf90_inquire_dimension(ncid, DimID, len = xdmy)
112	IF (xdmy /= ntimes) THEN
113	abort_message='NetCDF improperly dimensioned -- year_fr'
114	CALL abort_physic(modname,abort_message,1)
115	ENDIF
116	!
117	iret = nf90_inq_dimid(ncid, "years" , DimId)
118	iret = nf90_inquire_dimension(ncid, DimID, len = xdmy)
119	IF (xdmy /= nyears) THEN
120	abort_message='NetCDF improperly dimensioned -- years'
121	CALL abort_physic(modname,abort_message,1)
122	ENDIF
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/))
128	IF (iret /= NF90_NOERR) THEN
129	abort_message='NetCDF Error reading plume_lat'
130	CALL abort_physic(modname,abort_message,1)
131	ENDIF
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	iret = nf90_inq_varid(ncid, "beta_a" , VarId)
139	iret = nf90_get_var(ncid, VarID, beta_a(:) , start=(/1/),count=(/nplumes/))
140	IF (iret /= NF90_NOERR) THEN
141	abort_message='NetCDF Error reading beta_a'
142	CALL abort_physic(modname,abort_message,1)
143	ENDIF
144	iret = nf90_inq_varid(ncid, "beta_b" , VarId)
145	iret = nf90_get_var(ncid, VarID, beta_b(:) , start=(/1/),count=(/nplumes/))
146	IF (iret /= NF90_NOERR) THEN
147	abort_message='NetCDF Error reading beta_b'
148	CALL abort_physic(modname,abort_message,1)
149	ENDIF
150	iret = nf90_inq_varid(ncid, "aod_spmx" , VarId)
151	iret = nf90_get_var(ncid, VarID, aod_spmx(:) , start=(/1/),count=(/nplumes/))
152	IF (iret /= NF90_NOERR) THEN
153	abort_message='NetCDF Error reading aod_spmx'
154	CALL abort_physic(modname,abort_message,1)
155	ENDIF
156	iret = nf90_inq_varid(ncid, "aod_fmbg" , VarId)
157	iret = nf90_get_var(ncid, VarID, aod_fmbg(:) , start=(/1/),count=(/nplumes/))
158	IF (iret /= NF90_NOERR) THEN
159	abort_message='NetCDF Error reading aod_fmbg'
160	CALL abort_physic(modname,abort_message,1)
161	ENDIF
162	iret = nf90_inq_varid(ncid, "ssa550" , VarId)
163	iret = nf90_get_var(ncid, VarID, ssa550(:) , start=(/1/),count=(/nplumes/))
164	IF (iret /= NF90_NOERR) THEN
165	abort_message='NetCDF Error reading ssa550'
166	CALL abort_physic(modname,abort_message,1)
167	ENDIF
168	iret = nf90_inq_varid(ncid, "asy550" , VarId)
169	iret = nf90_get_var(ncid, VarID, asy550(:) , start=(/1/),count=(/nplumes/))
170	IF (iret /= NF90_NOERR) THEN
171	abort_message='NetCDF Error reading asy550'
172	CALL abort_physic(modname,abort_message,1)
173	ENDIF
174	iret = nf90_inq_varid(ncid, "angstrom" , VarId)
175	iret = nf90_get_var(ncid, VarID, angstrom(:), start=(/1/),count=(/nplumes/))
176	IF (iret /= NF90_NOERR) THEN
177	abort_message='NetCDF Error reading angstrom'
178	CALL abort_physic(modname,abort_message,1)
179	ENDIF
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/))
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
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/))
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
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/))
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
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/))
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
205	iret = nf90_inq_varid(ncid, "theta" , VarId)
206	iret = nf90_get_var(ncid, VarID, theta(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
207	IF (iret /= NF90_NOERR) THEN
208	abort_message='NetCDF Error reading theta'
209	CALL abort_physic(modname,abort_message,1)
210	ENDIF
211	iret = nf90_inq_varid(ncid, "ftr_weight" , VarId)
212	iret = nf90_get_var(ncid, VarID, ftr_weight(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/))
213	IF (iret /= NF90_NOERR) THEN
214	abort_message='NetCDF Error reading plume_lat'
215	CALL abort_physic(modname,abort_message,1)
216	ENDIF
217	iret = nf90_inq_varid(ncid, "year_weight" , VarId)
218	iret = nf90_get_var(ncid, VarID, year_weight(:,:) , start=(/1,1/),count=(/nyears,nplumes /))
219	IF (iret /= NF90_NOERR) THEN
220	abort_message='NetCDF Error reading year_weight'
221	CALL abort_physic(modname,abort_message,1)
222	ENDIF
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/))
225	IF (iret /= NF90_NOERR) THEN
226	abort_message='NetCDF Error reading ann_cycle'
227	CALL abort_physic(modname,abort_message,1)
228	ENDIF
229	!
230	iret = nf90_close(ncid)
231	!
232	ENDIF !--root processor
233	!
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)
251	!
252	sp_initialized = .TRUE.
253	!
254	RETURN
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