source: LMDZ6/trunk/libf/phylmd/mo_simple_plumes.F90 @ 4639

Last change on this file since 4639 was 4164, checked in by oboucher, 2 years ago

add a SAVE argument to variables read in simple aerosol plume model

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!
24MODULE 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
67CONTAINS
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 
498END MODULE MO_SIMPLE_PLUMES
Note: See TracBrowser for help on using the repository browser.