GCC Code Coverage Report


Directory: ./
File: phys/mo_simple_plumes.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 217 0.0%
Branches: 0 110 0.0%

Line Branch Exec Source
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
483