source: LMDZ6/branches/blowing_snow/libf/phylmd/cosp2/cosp_optics.F90 @ 5322

Last change on this file since 5322 was 3358, checked in by idelkadi, 6 years ago

Implementation de la nouvelle version COSPv2 dans LMDZ.
Pour compiler avec makelmdz_fcma utiliser l'option "-cosp2 true"
File size: 21.8 KB
Line 
1! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2! Copyright (c) 2015, Regents of the University of Colorado
3! All rights reserved.
4!
5! Redistribution and use in source and binary forms, with or without modification, are
6! permitted provided that the following conditions are met:
7!
8! 1. Redistributions of source code must retain the above copyright notice, this list of
9!    conditions and the following disclaimer.
10!
11! 2. Redistributions in binary form must reproduce the above copyright notice, this list
12!    of conditions and the following disclaimer in the documentation and/or other
13!    materials provided with the distribution.
14!
15! 3. Neither the name of the copyright holder nor the names of its contributors may be
16!    used to endorse or promote products derived from this software without specific prior
17!    written permission.
18!
19! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
20! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
22! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
24! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28!
29! History:
30! 05/01/15  Dustin Swales - Original version
31!
32! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33module cosp_optics
34  USE COSP_KINDS, ONLY: wp,dp
35  USE COSP_MATH_CONSTANTS,  ONLY: pi
36  USE COSP_PHYS_CONSTANTS,  ONLY: rholiq,km,rd,grav
37  USE MOD_MODIS_SIM,        ONLY: get_g_nir,get_ssa_nir,phaseIsLiquid,phaseIsIce
38  implicit none
39 
40  real(wp),parameter ::        & !
41       ice_density   = 0.93_wp   ! Ice density used in MODIS phase partitioning
42
43  interface cosp_simulator_optics
44     module procedure cosp_simulator_optics2D, cosp_simulator_optics3D
45  end interface cosp_simulator_optics
46 
47contains
48  ! ##########################################################################
49  !                          COSP_SIMULATOR_OPTICS
50  !
51  ! Used by: ISCCP, MISR and MODIS simulators
52  ! ##########################################################################
53  subroutine cosp_simulator_optics2D(dim1,dim2,dim3,flag,varIN1,varIN2,varOUT)
54    ! INPUTS
55    integer,intent(in) :: &
56         dim1,   & ! Dimension 1 extent (Horizontal)
57         dim2,   & ! Dimension 2 extent (Subcolumn)
58         dim3      ! Dimension 3 extent (Vertical)
59    real(wp),intent(in),dimension(dim1,dim2,dim3) :: &
60         flag      ! Logical to determine the of merge var1IN and var2IN
61    real(wp),intent(in),dimension(dim1,     dim3) :: &
62         varIN1, & ! Input field 1
63         varIN2    ! Input field 2
64    ! OUTPUTS
65    real(wp),intent(out),dimension(dim1,dim2,dim3) :: &
66         varOUT    ! Merged output field
67    ! LOCAL VARIABLES
68    integer :: j
69   
70    varOUT(1:dim1,1:dim2,1:dim3) = 0._wp
71    do j=1,dim2
72       where(flag(:,j,:) .eq. 1)
73          varOUT(:,j,:) = varIN2
74       endwhere
75       where(flag(:,j,:) .eq. 2)
76          varOUT(:,j,:) = varIN1
77       endwhere
78    enddo
79  end subroutine cosp_simulator_optics2D
80  subroutine cosp_simulator_optics3D(dim1,dim2,dim3,flag,varIN1,varIN2,varOUT)
81    ! INPUTS
82    integer,intent(in) :: &
83         dim1,   & ! Dimension 1 extent (Horizontal)
84         dim2,   & ! Dimension 2 extent (Subcolumn)
85         dim3      ! Dimension 3 extent (Vertical)
86    real(wp),intent(in),dimension(dim1,dim2,dim3) :: &
87         flag      ! Logical to determine the of merge var1IN and var2IN
88    real(wp),intent(in),dimension(dim1,dim2,dim3) :: &
89         varIN1, & ! Input field 1
90         varIN2    ! Input field 2
91    ! OUTPUTS
92    real(wp),intent(out),dimension(dim1,dim2,dim3) :: &
93         varOUT    ! Merged output field
94   
95    varOUT(1:dim1,1:dim2,1:dim3) = 0._wp
96   where(flag(:,:,:) .eq. 1)
97       varOUT(:,:,:) = varIN2
98    endwhere
99    where(flag(:,:,:) .eq. 2)
100       varOUT(:,:,:) = varIN1
101    endwhere
102   
103  end subroutine cosp_simulator_optics3D
104 
105  ! ##############################################################################
106  !                           MODIS_OPTICS_PARTITION
107  !
108  ! For the MODIS simulator, there are times when only a sinlge optical depth
109  ! profile, cloud-ice and cloud-water are provided. In this case, the optical
110  ! depth is partitioned by phase.
111  ! ##############################################################################
112  subroutine MODIS_OPTICS_PARTITION(npoints,nlev,ncolumns,cloudWater,cloudIce,waterSize, &
113                                    iceSize,tau,tauL,tauI)
114    ! INPUTS
115    INTEGER,intent(in) :: &
116         npoints,   & ! Number of horizontal gridpoints
117         nlev,      & ! Number of levels
118         ncolumns     ! Number of subcolumns
119    REAL(wp),intent(in),dimension(npoints,nlev,ncolumns) :: &
120         cloudWater, & ! Subcolumn cloud water content
121         cloudIce,   & ! Subcolumn cloud ice content
122         waterSize,  & ! Subcolumn cloud water effective radius
123         iceSize,    & ! Subcolumn cloud ice effective radius
124         tau           ! Optical thickness
125   
126    ! OUTPUTS
127    real(wp),intent(out),dimension(npoints,nlev,ncolumns) :: &
128         tauL,       & ! Partitioned liquid optical thickness.
129         tauI          ! Partitioned ice optical thickness.
130    ! LOCAL VARIABLES
131    real(wp),dimension(nlev,ncolumns) :: fracL
132    integer                           :: i
133   
134   
135    do i=1,npoints
136       where(cloudIce(i,:, :) <= 0.)
137          fracL(:, :) = 1._wp
138       elsewhere
139          where (cloudWater(i,:, :) <= 0.)
140             fracL(:, :) = 0._wp
141          elsewhere
142             ! Geometic optics limit - tau as LWP/re  (proportional to LWC/re)
143             fracL(:, :) = (cloudWater(i,:, :)/waterSize(i,:, :)) / &
144                  (cloudWater(i,:, :)/waterSize(i,:, :) + cloudIce(i,:, :)/(ice_density * iceSize(i,:, :)) )
145          end where
146       end where
147       tauL(i,:, :) = fracL(:, :) * tau(i,:, :)
148       tauI(i,:, :) = tau(i,:, :) - tauL(i,:, :)
149    enddo
150   
151  end subroutine MODIS_OPTICS_PARTITION
152  ! ########################################################################################
153  !                                   MODIS_OPTICS
154  !
155  ! ########################################################################################
156  subroutine modis_optics(nPoints,nLevels,nSubCols,tauLIQ,sizeLIQ,tauICE,sizeICE,fracLIQ, g, w0)
157    ! INPUTS
158    integer, intent(in)                                      :: nPoints,nLevels,nSubCols
159    real(wp),intent(in),dimension(nPoints,nSubCols,nLevels)  :: tauLIQ, sizeLIQ, tauICE, sizeICE
160    ! OUTPUTS
161    real(wp),intent(out),dimension(nPoints,nSubCols,nLevels) :: g,w0,fracLIQ
162    ! LOCAL VARIABLES
163    real(wp), dimension(nLevels)            :: water_g, water_w0, ice_g, ice_w0,tau
164    integer :: i,j
165   
166    ! Initialize
167    g(1:nPoints,1:nSubCols,1:nLevels)  = 0._wp
168    w0(1:nPoints,1:nSubCols,1:nLevels) = 0._wp
169   
170    do j =1,nPoints
171       do i=1,nSubCols
172          water_g(1:nLevels)  = get_g_nir(  phaseIsLiquid, sizeLIQ(j,i,1:nLevels))
173          water_w0(1:nLevels) = get_ssa_nir(phaseIsLiquid, sizeLIQ(j,i,1:nLevels))
174          ice_g(1:nLevels)    = get_g_nir(  phaseIsIce,    sizeICE(j,i,1:nLevels))
175          ice_w0(1:nLevels)   = get_ssa_nir(phaseIsIce,    sizeICE(j,i,1:nLevels))
176         
177          ! Combine ice and water optical properties
178          tau(1:nLevels) = tauICE(j,i,1:nLevels) + tauLIQ(j,i,1:nLevels)
179          where (tau(1:nLevels) > 0)
180             g(j,i,1:nLevels)  = (tauLIQ(j,i,1:nLevels)*water_g(1:nLevels) + tauICE(j,i,1:nLevels)*ice_g(1:nLevels)) / &
181                  tau(1:nLevels)
182             w0(j,i,1:nLevels) = (tauLIQ(j,i,1:nLevels)*water_g(1:nLevels)*water_w0(1:nLevels) + tauICE(j,i,1:nLevels) * &
183                  ice_g(1:nLevels) * ice_w0(1:nLevels)) / (g(j,i,1:nLevels) * tau(1:nLevels))
184          end where
185       enddo
186    enddo
187   
188    ! Compute the total optical thickness and the proportion due to liquid in each cell
189    do i=1,npoints
190       where(tauLIQ(i,1:nSubCols,1:nLevels) + tauICE(i,1:nSubCols,1:nLevels) > 0.)
191          fracLIQ(i,1:nSubCols,1:nLevels) = tauLIQ(i,1:nSubCols,1:nLevels)/ &
192               (tauLIQ(i,1:nSubCols,1:nLevels) + tauICE(i,1:nSubCols,1:nLevels))
193       elsewhere
194          fracLIQ(i,1:nSubCols,1:nLevels) = 0._wp
195       end  where
196    enddo
197   
198  end subroutine modis_optics
199 
200  ! ######################################################################################
201  ! SUBROUTINE lidar_optics
202  ! ######################################################################################
203  subroutine lidar_optics(npoints,ncolumns,nlev,npart,ice_type,q_lsliq, q_lsice,     &
204                              q_cvliq, q_cvice,ls_radliq,ls_radice,cv_radliq,cv_radice,  &
205                              pres,presf,temp,beta_mol,betatot,tau_mol,tautot,  &
206                              tautot_S_liq,tautot_S_ice,betatot_ice,betatot_liq,         &
207                              tautot_ice,tautot_liq)
208    ! ####################################################################################
209    ! NOTE: Using "grav" from cosp_constants.f90, instead of grav=9.81, introduces
210    ! changes of up to 2% in atb532 adn 0.003% in parasolRefl and lidarBetaMol532.
211    ! This also results in  small changes in the joint-histogram, cfadLidarsr532.
212    ! ####################################################################################
213   
214    ! INPUTS
215    INTEGER,intent(in) :: &
216         npoints,      & ! Number of gridpoints
217         ncolumns,     & ! Number of subcolumns
218         nlev,         & ! Number of levels
219         npart,        & ! Number of cloud meteors (stratiform_liq, stratiform_ice, conv_liq, conv_ice).
220         ice_type        ! Ice particle shape hypothesis (0 for spheres, 1 for non-spherical)
221    REAL(WP),intent(in),dimension(npoints,nlev) :: &
222         temp,         & ! Temperature of layer k
223         pres,         & ! Pressure at full levels
224         ls_radliq,    & ! Effective radius of LS liquid particles (meters)
225         ls_radice,    & ! Effective radius of LS ice particles (meters)
226         cv_radliq,    & ! Effective radius of CONV liquid particles (meters)
227         cv_radice       ! Effective radius of CONV ice particles (meters)
228    REAL(WP),intent(in),dimension(npoints,ncolumns,nlev) :: &
229         q_lsliq,      & ! LS sub-column liquid water mixing ratio (kg/kg)
230         q_lsice,      & ! LS sub-column ice water mixing ratio (kg/kg)
231         q_cvliq,      & ! CONV sub-column liquid water mixing ratio (kg/kg)
232         q_cvice         ! CONV sub-column ice water mixing ratio (kg/kg)
233    REAL(WP),intent(in),dimension(npoints,nlev+1) :: &
234         presf           ! Pressure at half levels
235   
236    ! OUTPUTS
237    REAL(WP),intent(out),dimension(npoints,ncolumns,nlev)       :: &
238         betatot,        & !
239         tautot            ! Optical thickess integrated from top
240    REAL(WP),intent(out),dimension(npoints,ncolumns,nlev)       :: &
241         betatot_ice,    & ! Backscatter coefficient for ice particles
242         betatot_liq,    & ! Backscatter coefficient for liquid particles
243         tautot_ice,     & ! Total optical thickness of ice
244         tautot_liq        ! Total optical thickness of liq
245    REAL(WP),intent(out),dimension(npoints,nlev) :: &
246         beta_mol,       & ! Molecular backscatter coefficient
247         tau_mol           ! Molecular optical depth
248    REAL(WP),intent(out),dimension(npoints,ncolumns) :: &
249         tautot_S_liq,   & ! TOA optical depth for liquid
250         tautot_S_ice      ! TOA optical depth for ice
251   
252    ! LOCAL VARIABLES
253    REAL(WP),dimension(npart)                       :: rhopart
254    REAL(WP),dimension(npart,5)                     :: polpart
255    REAL(WP),dimension(npoints,nlev)                :: rhoair,alpha_mol
256    REAL(WP),dimension(npoints,nlev+1)              :: zheight         
257    REAL(WP),dimension(npoints,nlev,npart)          :: rad_part,kp_part,qpart
258    REAL(WP),dimension(npoints,ncolumns,nlev,npart) :: alpha_part,tau_part
259    INTEGER                                         :: i,k,icol
260   
261    ! Local data
262    REAL(WP),PARAMETER :: rhoice     = 0.5e+03    ! Density of ice (kg/m3)
263    REAL(WP),PARAMETER :: Cmol       = 6.2446e-32 ! Wavelength dependent
264    REAL(WP),PARAMETER :: rdiffm     = 0.7_wp     ! Multiple scattering correction parameter
265    REAL(WP),PARAMETER :: Qscat      = 2.0_wp     ! Particle scattering efficiency at 532 nm
266    ! Local indicies for large-scale and convective ice and liquid
267    INTEGER,PARAMETER  :: INDX_LSLIQ = 1
268    INTEGER,PARAMETER  :: INDX_LSICE = 2
269    INTEGER,PARAMETER  :: INDX_CVLIQ = 3
270    INTEGER,PARAMETER  :: INDX_CVICE = 4
271   
272    ! Polarized optics parameterization
273    ! Polynomial coefficients for spherical liq/ice particles derived from Mie theory.
274    ! Polynomial coefficients for non spherical particles derived from a composite of
275    ! Ray-tracing theory for large particles (e.g. Noel et al., Appl. Opt., 2001)
276    ! and FDTD theory for very small particles (Yang et al., JQSRT, 2003).
277    ! We repeat the same coefficients for LS and CONV cloud to make code more readable
278    REAL(WP),PARAMETER,dimension(5) :: &
279         polpartCVLIQ  = (/ 2.6980e-8_wp,  -3.7701e-6_wp,  1.6594e-4_wp,    -0.0024_wp,    0.0626_wp/), &
280         polpartLSLIQ  = (/ 2.6980e-8_wp,  -3.7701e-6_wp,  1.6594e-4_wp,    -0.0024_wp,    0.0626_wp/), &
281         polpartCVICE0 = (/-1.0176e-8_wp,   1.7615e-6_wp, -1.0480e-4_wp,     0.0019_wp,    0.0460_wp/), &
282         polpartLSICE0 = (/-1.0176e-8_wp,   1.7615e-6_wp, -1.0480e-4_wp,     0.0019_wp,    0.0460_wp/), &
283         polpartCVICE1 = (/ 1.3615e-8_wp, -2.04206e-6_wp, 7.51799e-5_wp, 0.00078213_wp, 0.0182131_wp/), &
284         polpartLSICE1 = (/ 1.3615e-8_wp, -2.04206e-6_wp, 7.51799e-5_wp, 0.00078213_wp, 0.0182131_wp/)
285    ! ##############################################################################
286   
287    ! Liquid/ice particles
288    rhopart(INDX_LSLIQ) = rholiq
289    rhopart(INDX_LSICE) = rhoice
290    rhopart(INDX_CVLIQ) = rholiq
291    rhopart(INDX_CVICE) = rhoice
292   
293    ! LS and CONV Liquid water coefficients
294    polpart(INDX_LSLIQ,1:5) = polpartLSLIQ
295    polpart(INDX_CVLIQ,1:5) = polpartCVLIQ
296    ! LS and CONV Ice water coefficients
297    if (ice_type .eq. 0) then
298       polpart(INDX_LSICE,1:5) = polpartLSICE0
299       polpart(INDX_CVICE,1:5) = polpartCVICE0
300    endif
301    if (ice_type .eq. 1) then
302       polpart(INDX_LSICE,1:5) = polpartLSICE1
303       polpart(INDX_CVICE,1:5) = polpartCVICE1
304    endif
305   
306    ! Effective radius particles:
307    rad_part(1:npoints,1:nlev,INDX_LSLIQ) = ls_radliq(1:npoints,1:nlev)
308    rad_part(1:npoints,1:nlev,INDX_LSICE) = ls_radice(1:npoints,1:nlev)
309    rad_part(1:npoints,1:nlev,INDX_CVLIQ) = cv_radliq(1:npoints,1:nlev)
310    rad_part(1:npoints,1:nlev,INDX_CVICE) = cv_radice(1:npoints,1:nlev)   
311    rad_part(1:npoints,1:nlev,1:npart)    = MAX(rad_part(1:npoints,1:nlev,1:npart),0._wp)
312    rad_part(1:npoints,1:nlev,1:npart)    = MIN(rad_part(1:npoints,1:nlev,1:npart),70.0e-6_wp)
313   
314    ! Density (clear-sky air)
315    rhoair(1:npoints,1:nlev) = pres(1:npoints,1:nlev)/(rd*temp(1:npoints,1:nlev))
316   
317    ! Altitude at half pressure levels:
318    zheight(1:npoints,nlev+1) = 0._wp
319    do k=nlev,1,-1
320       zheight(1:npoints,k) = zheight(1:npoints,k+1) &
321            -(presf(1:npoints,k)-presf(1:npoints,k+1))/(rhoair(1:npoints,k)*grav)
322    enddo
323   
324    ! ##############################################################################
325    ! *) Molecular alpha, beta and optical thickness
326    ! ##############################################################################
327   
328    beta_mol(1:npoints,1:nlev)  = pres(1:npoints,1:nlev)/km/temp(1:npoints,1:nlev)*Cmol
329    alpha_mol(1:npoints,1:nlev) = 8._wp*pi/3._wp * beta_mol(1:npoints,1:nlev)
330   
331    ! Optical thickness of each layer (molecular) 
332    tau_mol(1:npoints,1:nlev) = alpha_mol(1:npoints,1:nlev)*(zheight(1:npoints,1:nlev)-&
333         zheight(1:npoints,2:nlev+1))
334   
335    ! Optical thickness from TOA to layer k (molecular)
336    DO k = 2,nlev
337       tau_mol(1:npoints,k) = tau_mol(1:npoints,k) + tau_mol(1:npoints,k-1)
338    ENDDO
339   
340    betatot    (1:npoints,1:ncolumns,1:nlev) = spread(beta_mol(1:npoints,1:nlev), dim=2, NCOPIES=ncolumns)
341    tautot     (1:npoints,1:ncolumns,1:nlev) = spread(tau_mol (1:npoints,1:nlev), dim=2, NCOPIES=ncolumns)
342    betatot_liq(1:npoints,1:ncolumns,1:nlev) = betatot(1:npoints,1:ncolumns,1:nlev)
343    betatot_ice(1:npoints,1:ncolumns,1:nlev) = betatot(1:npoints,1:ncolumns,1:nlev)
344    tautot_liq (1:npoints,1:ncolumns,1:nlev) = tautot(1:npoints,1:ncolumns,1:nlev)
345    tautot_ice (1:npoints,1:ncolumns,1:nlev) = tautot(1:npoints,1:ncolumns,1:nlev)
346   
347    ! ##############################################################################
348    ! *) Particles alpha, beta and optical thickness
349    ! ##############################################################################
350    ! Polynomials kp_lidar derived from Mie theory
351    do i = 1, npart
352       where (rad_part(1:npoints,1:nlev,i) .gt. 0.0)
353          kp_part(1:npoints,1:nlev,i) = &
354               polpart(i,1)*(rad_part(1:npoints,1:nlev,i)*1e6)**4 &
355               + polpart(i,2)*(rad_part(1:npoints,1:nlev,i)*1e6)**3 &
356               + polpart(i,3)*(rad_part(1:npoints,1:nlev,i)*1e6)**2 &
357               + polpart(i,4)*(rad_part(1:npoints,1:nlev,i)*1e6) &
358               + polpart(i,5)
359       elsewhere
360          kp_part(1:npoints,1:nlev,i) = 0._wp
361       endwhere
362    enddo
363   
364    ! Loop over all subcolumns
365    do icol=1,ncolumns
366       ! ##############################################################################
367       ! Mixing ratio particles in each subcolum
368       ! ##############################################################################
369       qpart(1:npoints,1:nlev,INDX_LSLIQ) = q_lsliq(1:npoints,icol,1:nlev)
370       qpart(1:npoints,1:nlev,INDX_LSICE) = q_lsice(1:npoints,icol,1:nlev)
371       qpart(1:npoints,1:nlev,INDX_CVLIQ) = q_cvliq(1:npoints,icol,1:nlev)
372       qpart(1:npoints,1:nlev,INDX_CVICE) = q_cvice(1:npoints,icol,1:nlev)
373       
374       ! ##############################################################################
375       ! Alpha and optical thickness (particles)
376       ! ##############################################################################
377       ! Alpha of particles in each subcolumn:
378       do i = 1, npart
379          where (rad_part(1:npoints,1:nlev,i) .gt. 0.0)
380             alpha_part(1:npoints,icol,1:nlev,i) = 3._wp/4._wp * Qscat &
381                  * rhoair(1:npoints,1:nlev) * qpart(1:npoints,1:nlev,i) &
382                  / (rhopart(i) * rad_part(1:npoints,1:nlev,i) )
383          elsewhere
384             alpha_part(1:npoints,icol,1:nlev,i) = 0._wp
385          endwhere
386       enddo
387       
388       ! Optical thicknes
389       tau_part(1:npoints,icol,1:nlev,1:npart) = rdiffm * alpha_part(1:npoints,icol,1:nlev,1:npart)
390       do i = 1, npart
391          ! Optical thickness of each layer (particles)
392          tau_part(1:npoints,icol,1:nlev,i) = tau_part(1:npoints,icol,1:nlev,i) &
393               & * (zheight(1:npoints,1:nlev)-zheight(1:npoints,2:nlev+1) )
394          ! Optical thickness from TOA to layer k (particles)
395          do k=2,nlev
396             tau_part(1:npoints,icol,k,i) = tau_part(1:npoints,icol,k,i) + tau_part(1:npoints,icol,k-1,i)
397          enddo
398       enddo
399       
400       ! ##############################################################################
401       ! Beta and optical thickness (total=molecular + particules)
402       ! ##############################################################################
403       
404       DO i = 1, npart
405          betatot(1:npoints,icol,1:nlev) = betatot(1:npoints,icol,1:nlev) + &
406               kp_part(1:npoints,1:nlev,i)*alpha_part(1:npoints,icol,1:nlev,i)
407          tautot(1:npoints,icol,1:nlev) = tautot(1:npoints,icol,1:nlev)  + &
408               tau_part(1:npoints,icol,1:nlev,i)
409       ENDDO
410       
411       ! ##############################################################################
412       ! Beta and optical thickness (liquid/ice)
413       ! ##############################################################################
414       ! Ice
415       betatot_ice(1:npoints,icol,1:nlev) = betatot_ice(1:npoints,icol,1:nlev)+ &
416            kp_part(1:npoints,1:nlev,INDX_LSICE)*alpha_part(1:npoints,icol,1:nlev,INDX_LSICE)+ &
417            kp_part(1:npoints,1:nlev,INDX_CVICE)*alpha_part(1:npoints,icol,1:nlev,INDX_CVICE)
418       tautot_ice(1:npoints,icol,1:nlev) = tautot_ice(1:npoints,icol,1:nlev)  + &
419            tau_part(1:npoints,icol,1:nlev,INDX_LSICE) + &
420            tau_part(1:npoints,icol,1:nlev,INDX_CVICE)
421       
422       ! Liquid
423       betatot_liq(1:npoints,icol,1:nlev) = betatot_liq(1:npoints,icol,1:nlev)+ &
424            kp_part(1:npoints,1:nlev,INDX_LSLIQ)*alpha_part(1:npoints,icol,1:nlev,INDX_LSLIQ)+ &
425            kp_part(1:npoints,1:nlev,INDX_CVLIQ)*alpha_part(1:npoints,icol,1:nlev,INDX_CVLIQ)
426       tautot_liq(1:npoints,icol,1:nlev) = tautot_liq(1:npoints,icol,1:nlev)  + &
427            tau_part(1:npoints,icol,1:nlev,INDX_LSLIQ) + &
428            tau_part(1:npoints,icol,1:nlev,INDX_CVLIQ)
429    enddo
430   
431    ! ##############################################################################   
432    ! Optical depths used by the PARASOL simulator
433    ! ##############################################################################   
434    tautot_S_liq(1:npoints,1:ncolumns) = 0._wp
435    tautot_S_ice(1:npoints,1:ncolumns) = 0._wp
436    do icol=1,ncolumns   
437       tautot_S_liq(1:npoints,icol) = tautot_S_liq(1:npoints,icol)+tau_part(1:npoints,icol,nlev,1)+tau_part(1:npoints,icol,nlev,3)
438       tautot_S_ice(1:npoints,icol) = tautot_S_ice(1:npoints,icol)+tau_part(1:npoints,icol,nlev,2)+tau_part(1:npoints,icol,nlev,4)
439    enddo
440   
441  end subroutine lidar_optics
442end module cosp_optics
Note: See TracBrowser for help on using the repository browser.