source: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/mo_simple_plumes.F90

Last change on this file was 3318, checked in by jghattas, 7 years ago

Integration des developpements fait sur la trunk par O. Boucher concernant le MACspV2 aerosol plume climatology(nouveau option flag_aersol=7). Les commits suivants fait sur le trunk sont ici merge : [3274], [3279], [3287], [3288], [3290], [3295], [3296], [3297].
Tout les modifications dans newmicro.f90 ne sont pas retenu mais les changemnets lie au flag_aerosol=7 sont prise.
File size: 20.2 KB
RevLine 
[3318]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    !
81    ! ----------
82    !--only one processor reads the input data
83    IF (is_mpi_root.AND.is_omp_root) THEN
84       !   
85       iret = nf90_open("MACv2.0-SP_v1.nc", NF90_NOWRITE, ncid)
86       IF (iret /= NF90_NOERR) STOP 'NetCDF File not opened'
87       !
88       ! read dimensions and make sure file conforms to expected size
89       !
90       iret = nf90_inq_dimid(ncid, "plume_number"  , DimId)
91       iret = nf90_inquire_dimension(ncid, DimId, len = xdmy)
92       IF (xdmy /= nplumes) STOP 'NetCDF improperly dimensioned -- plume_number'
93       !
94       iret = nf90_inq_dimid(ncid, "plume_feature", DimId)
95       iret = nf90_inquire_dimension(ncid, DimId, len = xdmy)
96       IF (xdmy /= nfeatures) STOP 'NetCDF improperly dimensioned -- plume_feature'
97       !
98       iret = nf90_inq_dimid(ncid, "year_fr"   , DimId)
99       iret = nf90_inquire_dimension(ncid, DimID, len = xdmy)
100       IF (xdmy /= ntimes) STOP 'NetCDF improperly dimensioned -- year_fr'
101       !
102       iret = nf90_inq_dimid(ncid, "years"   , DimId)
103       iret = nf90_inquire_dimension(ncid, DimID, len = xdmy)
104       IF (xdmy /= nyears) STOP 'NetCDF improperly dimensioned -- years'
105       !
106       ! read variables that define the simple plume climatology
107       !
108       iret = nf90_inq_varid(ncid, "plume_lat", VarId)
109       iret = nf90_get_var(ncid, VarID, plume_lat(:), start=(/1/),count=(/nplumes/))
110       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading plume_lat'
111       iret = nf90_inq_varid(ncid, "plume_lon", VarId)
112       iret = nf90_get_var(ncid, VarID, plume_lon(:), start=(/1/),count=(/nplumes/))
113       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading plume_lon'
114       iret = nf90_inq_varid(ncid, "beta_a"   , VarId)
115       iret = nf90_get_var(ncid, VarID, beta_a(:)   , start=(/1/),count=(/nplumes/))
116       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading beta_a'
117       iret = nf90_inq_varid(ncid, "beta_b"   , VarId)
118       iret = nf90_get_var(ncid, VarID, beta_b(:)   , start=(/1/),count=(/nplumes/))
119       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading beta_b'
120       iret = nf90_inq_varid(ncid, "aod_spmx" , VarId)
121       iret = nf90_get_var(ncid, VarID, aod_spmx(:)  , start=(/1/),count=(/nplumes/))
122       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading aod_spmx'
123       iret = nf90_inq_varid(ncid, "aod_fmbg" , VarId)
124       iret = nf90_get_var(ncid, VarID, aod_fmbg(:)  , start=(/1/),count=(/nplumes/))
125       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading aod_fmbg'
126       iret = nf90_inq_varid(ncid, "ssa550"   , VarId)
127       iret = nf90_get_var(ncid, VarID, ssa550(:)  , start=(/1/),count=(/nplumes/))
128       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading ssa550'
129       iret = nf90_inq_varid(ncid, "asy550"   , VarId)
130       iret = nf90_get_var(ncid, VarID, asy550(:)  , start=(/1/),count=(/nplumes/))
131       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading asy550'
132       iret = nf90_inq_varid(ncid, "angstrom" , VarId)
133       iret = nf90_get_var(ncid, VarID, angstrom(:), start=(/1/),count=(/nplumes/))
134       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading angstrom'
135       !
136       iret = nf90_inq_varid(ncid, "sig_lat_W"     , VarId)
137       iret = nf90_get_var(ncid, VarID, sig_lat_W(:,:)    , start=(/1,1/),count=(/nfeatures,nplumes/))
138       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lat_W'
139       iret = nf90_inq_varid(ncid, "sig_lat_E"     , VarId)
140       iret = nf90_get_var(ncid, VarID, sig_lat_E(:,:)    , start=(/1,1/),count=(/nfeatures,nplumes/))
141       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lat_E'
142       iret = nf90_inq_varid(ncid, "sig_lon_E"     , VarId)
143       iret = nf90_get_var(ncid, VarID, sig_lon_E(:,:)    , start=(/1,1/),count=(/nfeatures,nplumes/))
144       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lon_E'
145       iret = nf90_inq_varid(ncid, "sig_lon_W"     , VarId)
146       iret = nf90_get_var(ncid, VarID, sig_lon_W(:,:)    , start=(/1,1/),count=(/nfeatures,nplumes/))
147       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lon_W'
148       iret = nf90_inq_varid(ncid, "theta"         , VarId)
149       iret = nf90_get_var(ncid, VarID, theta(:,:)        , start=(/1,1/),count=(/nfeatures,nplumes/))
150       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading theta'
151       iret = nf90_inq_varid(ncid, "ftr_weight"    , VarId)
152       iret = nf90_get_var(ncid, VarID, ftr_weight(:,:)   , start=(/1,1/),count=(/nfeatures,nplumes/))
153       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading plume_lat'
154       iret = nf90_inq_varid(ncid, "year_weight"   , VarId)
155       iret = nf90_get_var(ncid, VarID, year_weight(:,:)  , start=(/1,1/),count=(/nyears,nplumes   /))
156       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading year_weight'
157       iret = nf90_inq_varid(ncid, "ann_cycle"     , VarId)
158       iret = nf90_get_var(ncid, VarID, ann_cycle(:,:,:)  , start=(/1,1,1/),count=(/nfeatures,ntimes,nplumes/))
159       IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading ann_cycle'
160       !
161       iret = nf90_close(ncid)
162       !
163    ENDIF !--root processor
164    !
165    CALL bcast(plume_lat)
166    CALL bcast(plume_lon)
167    CALL bcast(beta_a)
168    CALL bcast(beta_b)
169    CALL bcast(aod_spmx)
170    CALL bcast(aod_fmbg)
171    CALL bcast(asy550)
172    CALL bcast(ssa550)
173    CALL bcast(angstrom)
174    CALL bcast(sig_lon_E)
175    CALL bcast(sig_lon_W)
176    CALL bcast(sig_lat_E)
177    CALL bcast(sig_lat_W)
178    CALL bcast(theta)
179    CALL bcast(ftr_weight)
180    CALL bcast(year_weight)
181    CALL bcast(ann_cycle)
182    !
183    sp_initialized = .TRUE.
184    !
185    RETURN
186    !
187  END SUBROUTINE sp_setup
188  !
189  ! ------------------------------------------------------------------------------------------------------------------------
190  ! SET_TIME_WEIGHT:  The simple plume model assumes that meteorology constrains plume shape and that only source strength
191  ! influences the amplitude of a plume associated with a given source region.   This routine retrieves the temporal weights
192  ! for the plumes.  Each plume feature has its own temporal weights which varies yearly.  The annual cycle is indexed by
193  ! week in the year and superimposed on the yearly mean value of the weight.
194  !
195  SUBROUTINE set_time_weight(year_fr)
196    !
197    ! ----------
198    !
199    REAL, INTENT(IN) ::  &
200         year_fr           !< Fractional Year (1850.0 - 2100.99)
201
202    INTEGER          ::  &
203         iyear          ,& !< Integer year values between 1 and 156 (1850-2100)
204         iweek          ,& !< Integer index (between 1 and ntimes); for ntimes=52 this corresponds to weeks (roughly)
205         iplume            ! plume number
206    !
207    ! ----------
208    !
209    iyear = FLOOR(year_fr) - 1849
210    iweek = FLOOR((year_fr - FLOOR(year_fr)) * ntimes) + 1
211
212    IF ((iweek > ntimes) .OR. (iweek < 1) .OR. (iyear > nyears) .OR. (iyear < 1)) THEN
213      CALL abort_physic('set_time_weight','Time out of bounds')
214    ENDIF
215
216    DO iplume=1,nplumes
217      time_weight(1,iplume) = year_weight(iyear,iplume) * ann_cycle(1,iweek,iplume)
218      time_weight(2,iplume) = year_weight(iyear,iplume) * ann_cycle(2,iweek,iplume)
219      time_weight_bg(1,iplume) = ann_cycle(1,iweek,iplume)
220      time_weight_bg(2,iplume) = ann_cycle(2,iweek,iplume)
221    END DO
222   
223    RETURN
224  END SUBROUTINE set_time_weight
225  !
226  ! ------------------------------------------------------------------------------------------------------------------------
227  ! SP_AOP_PROFILE:  This subroutine calculates the simple plume aerosol and cloud active optical properties based on the
228  ! the simple plume fit to the MPI Aerosol Climatology (Version 2).  It sums over nplumes to provide a profile of aerosol
229  ! optical properties on a host models vertical grid.
230  !
231  SUBROUTINE sp_aop_profile                                                                           ( &
232       nlevels        ,ncol           ,lambda         ,oro            ,lon            ,lat            , &
233       year_fr        ,z              ,dz             ,dNovrN         ,aod_prof       ,ssa_prof       , &
234       asy_prof       )
235    !
236    ! ----------
237    !
238    INTEGER, INTENT(IN)        :: &
239         nlevels,                 & !< number of levels
240         ncol                       !< number of columns
241
242    REAL, INTENT(IN)           :: &
243         lambda,                  & !< wavelength
244         year_fr,                 & !< Fractional Year (1903.0 is the 0Z on the first of January 1903, Gregorian)
245         oro(ncol),               & !< orographic height (m)
246         lon(ncol),               & !< longitude
247         lat(ncol),               & !< latitude
248         z (ncol,nlevels),        & !< height above sea-level (m)
249         dz(ncol,nlevels)           !< level thickness (difference between half levels) (m)
250
251    REAL, INTENT(OUT)          :: &
252         dNovrN(ncol)           , & !< anthropogenic increase in cloud drop number concentration (factor)
253         aod_prof(ncol,nlevels) , & !< profile of aerosol optical depth
254         ssa_prof(ncol,nlevels) , & !< profile of single scattering albedo
255         asy_prof(ncol,nlevels)     !< profile of asymmetry parameter
256
257    INTEGER                    :: iplume, icol, k
258
259    REAL                       ::  &
260         eta(ncol,nlevels),        & !< normalized height (by 15 km)
261         z_beta(ncol,nlevels),     & !< profile for scaling column optical depth
262         prof(ncol,nlevels),       & !< scaled profile (by beta function)
263         beta_sum(ncol),           & !< vertical sum of beta function
264         ssa(ncol),                & !< single scattering albedo
265         asy(ncol),                & !< asymmetry parameter
266         cw_an(ncol),              & !< column weight for simple plume (anthropogenic) AOD at 550 nm
267         cw_bg(ncol),              & !< column weight for fine-mode natural background AOD at 550 nm
268         caod_sp(ncol),            & !< column simple plume anthropogenic AOD at 550 nm
269         caod_bg(ncol),            & !< column fine-mode natural background AOD at 550 nm
270         a_plume1,                 & !< gaussian longitude factor for feature 1
271         a_plume2,                 & !< gaussian longitude factor for feature 2
272         b_plume1,                 & !< gaussian latitude factor for feature 1
273         b_plume2,                 & !< gaussian latitude factor for feature 2
274         delta_lat,                & !< latitude offset
275         delta_lon,                & !< longitude offset
276         delta_lon_t,              & !< threshold for maximum longitudinal plume extent used in transition from 360 to 0 degrees
277         lon1,                     & !< rotated longitude for feature 1
278         lat1,                     & !< rotated latitude for feature 2
279         lon2,                     & !< rotated longitude for feature 1
280         lat2,                     & !< rotated latitude for feature 2
281         f1,                       & !< contribution from feature 1
282         f2,                       & !< contribution from feature 2
283         f3,                       & !< contribution from feature 1 in natural background of Twomey effect
284         f4,                       & !< contribution from feature 2 in natural background of Twomey effect
285         aod_550,                  & !< aerosol optical depth at 550nm
286         aod_lmd,                  & !< aerosol optical depth at input wavelength
287         lfactor                     !< factor to compute wavelength dependence of optical properties
288    !
289    ! ----------
290    !
291    ! initialize input data (by calling setup at first instance)
292    !
293    IF (.NOT.sp_initialized) CALL sp_setup
294    !
295    ! get time weights
296    !
297    CALL set_time_weight(year_fr)
298    !
299    ! initialize variables, including output
300    !
301    DO k=1,nlevels
302      DO icol=1,ncol
303        aod_prof(icol,k) = 0.0
304        ssa_prof(icol,k) = 0.0
305        asy_prof(icol,k) = 0.0
306        z_beta(icol,k)   = MERGE(1.0, 0.0, z(icol,k) >= oro(icol))
307        eta(icol,k)      = MAX(0.0,MIN(1.0,z(icol,k)/15000.))
308      END DO
309    END DO
310    DO icol=1,ncol
311      dNovrN(icol)   = 1.0
312      caod_sp(icol)  = 0.0
313      caod_bg(icol)  = 0.02
314    END DO
315    !
316    ! sum contribution from plumes to construct composite profiles of aerosol optical properties
317    !
318    DO iplume=1,nplumes
319      !
320      ! calculate vertical distribution function from parameters of beta distribution
321      !
322      DO icol=1,ncol
323        beta_sum(icol) = 0.
324      END DO
325      DO k=1,nlevels
326        DO icol=1,ncol
327          prof(icol,k)   = (eta(icol,k)**(beta_a(iplume)-1.) * (1.-eta(icol,k))**(beta_b(iplume)-1.)) * dz(icol,k)
328          beta_sum(icol) = beta_sum(icol) + prof(icol,k)
329        END DO
330      END DO
331      DO k=1,nlevels
332        DO icol=1,ncol
333          prof(icol,k)   = ( prof(icol,k) / beta_sum(icol) ) * z_beta(icol,k)
334        END DO
335      END DO
336      !
337      ! calculate plume weights
338      !
339      DO icol=1,ncol
340        !
341        ! get plume-center relative spatial parameters for specifying amplitude of plume at given lat and lon
342        !
343        delta_lat   = lat(icol) - plume_lat(iplume)
344        delta_lon   = lon(icol) - plume_lon(iplume)
345        delta_lon_t = MERGE (260., 180., iplume == 1)
346        delta_lon   = MERGE ( delta_lon-SIGN(360.,delta_lon) , delta_lon , ABS(delta_lon) > delta_lon_t)
347
348        a_plume1  = 0.5 / (MERGE(sig_lon_E(1,iplume), sig_lon_W(1,iplume), delta_lon > 0)**2)
349        b_plume1  = 0.5 / (MERGE(sig_lat_E(1,iplume), sig_lat_W(1,iplume), delta_lon > 0)**2)
350        a_plume2  = 0.5 / (MERGE(sig_lon_E(2,iplume), sig_lon_W(2,iplume), delta_lon > 0)**2)
351        b_plume2  = 0.5 / (MERGE(sig_lat_E(2,iplume), sig_lat_W(2,iplume), delta_lon > 0)**2)
352        !
353        ! adjust for a plume specific rotation which helps match plume state to climatology.
354        !
355        lon1 =   COS(theta(1,iplume))*(delta_lon) + SIN(theta(1,iplume))*(delta_lat)
356        lat1 = - SIN(theta(1,iplume))*(delta_lon) + COS(theta(1,iplume))*(delta_lat)
357        lon2 =   COS(theta(2,iplume))*(delta_lon) + SIN(theta(2,iplume))*(delta_lat)
358        lat2 = - SIN(theta(2,iplume))*(delta_lon) + COS(theta(2,iplume))*(delta_lat)
359        !
360        ! calculate contribution to plume from its different features, to get a column weight for the anthropogenic
361        ! (cw_an) and the fine-mode natural background aerosol (cw_bg)
362        !
363        f1 = time_weight(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)**2) + (b_plume1 * ((lat1)**2))))
364        f2 = time_weight(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)**2) + (b_plume2 * ((lat2)**2))))
365        f3 = time_weight_bg(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)**2) + (b_plume1 * ((lat1)**2))))
366        f4 = time_weight_bg(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)**2) + (b_plume2 * ((lat2)**2))))
367
368        cw_an(icol) = f1 * aod_spmx(iplume) + f2 * aod_spmx(iplume) 
369        cw_bg(icol) = f3 * aod_fmbg(iplume) + f4 * aod_fmbg(iplume)
370        !
371        ! calculate wavelength-dependent scattering properties
372        !
373        lfactor   = MIN(1.0,700.0/lambda)
374        ssa(icol) = (ssa550(iplume) * lfactor**4) / ((ssa550(iplume) * lfactor**4) + ((1-ssa550(iplume)) * lfactor))
375        asy(icol) =  asy550(iplume) * SQRT(lfactor)
376      END DO
377      !
378      ! distribute plume optical properties across its vertical profile weighting by optical depth and scaling for
379      ! wavelength using the angstrom parameter.
380      !     
381      lfactor = EXP(-angstrom(iplume) * LOG(lambda/550.0))
382      DO k=1,nlevels
383        DO icol = 1,ncol
384          aod_550          = prof(icol,k)     * cw_an(icol)
385          aod_lmd          = aod_550          * lfactor
386          caod_sp(icol)    = caod_sp(icol)    + aod_550
387          caod_bg(icol)    = caod_bg(icol)    + prof(icol,k) * cw_bg(icol)
388          asy_prof(icol,k) = asy_prof(icol,k) + aod_lmd * ssa(icol) * asy(icol)
389          ssa_prof(icol,k) = ssa_prof(icol,k) + aod_lmd * ssa(icol)
390          aod_prof(icol,k) = aod_prof(icol,k) + aod_lmd
391        END DO
392      END DO
393    END DO
394    !
395    ! complete optical depth weighting
396    !
397    DO k=1,nlevels
398      DO icol = 1,ncol
399        asy_prof(icol,k) = MERGE(asy_prof(icol,k)/ssa_prof(icol,k), 0.0, ssa_prof(icol,k) > TINY(1.))
400        ssa_prof(icol,k) = MERGE(ssa_prof(icol,k)/aod_prof(icol,k), 1.0, aod_prof(icol,k) > TINY(1.))
401      END DO
402    END DO
403    !
404    ! calculate effective radius normalization (divisor) factor
405    !
406    DO icol=1,ncol
407      dNovrN(icol) = LOG((1000.0 * (caod_sp(icol) + caod_bg(icol))) + 1.0)/LOG((1000.0 * caod_bg(icol)) + 1.0)
408    END DO
409
410    RETURN
411  END SUBROUTINE sp_aop_profile
412 
413END MODULE MO_SIMPLE_PLUMES
Note: See TracBrowser for help on using the repository browser.