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

Last change on this file since 3709 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!
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 ::                                      &
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
65CONTAINS
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 
482END MODULE MO_SIMPLE_PLUMES
Note: See TracBrowser for help on using the repository browser.