source: LMDZ6/branches/Amaury_dev/libf/phylmd/cospv2/cosp_optics.F90 @ 5218

Last change on this file since 5218 was 5185, checked in by abarral, 3 months ago

Replace REPROBUS CPP KEY by logical using handmade wonky wrapper
File size: 23.2 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! 04/04/18  Rodrigo Guzman- Added CALIOP-like Ground LIDar routines (GLID)
32! 10/04/18  Rodrigo Guzman- Added ATLID-like (EarthCare) lidar routines (ATLID)
33
34! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35module cosp_optics
36  USE COSP_KINDS, ONLY: wp,dp
37  USE COSP_MATH_CONSTANTS,  ONLY: pi
38  USE COSP_PHYS_CONSTANTS,  ONLY: rholiq,km,rd,grav
39  USE MOD_MODIS_SIM,        ONLY: get_g_nir,get_ssa_nir,phaseIsLiquid,phaseIsIce
40  implicit none
41 
42  real(wp),parameter ::        & !
43       ice_density   = 0.93_wp   ! Ice density used in MODIS phase partitioning
44
45  interface cosp_simulator_optics
46     module procedure cosp_simulator_optics2D, cosp_simulator_optics3D
47  end interface cosp_simulator_optics
48 
49contains
50  ! ##########################################################################
51  !                          COSP_SIMULATOR_OPTICS
52
53  ! Used by: ISCCP, MISR and MODIS simulators
54  ! ##########################################################################
55  subroutine cosp_simulator_optics2D(dim1,dim2,dim3,flag,varIN1,varIN2,varOUT)
56    ! INPUTS
57    integer,intent(in) :: &
58         dim1,   & ! Dimension 1 extent (Horizontal)
59         dim2,   & ! Dimension 2 extent (Subcolumn)
60         dim3      ! Dimension 3 extent (Vertical)
61    real(wp),intent(in),dimension(dim1,dim2,dim3) :: &
62         flag      ! Logical to determine the of merge var1IN and var2IN
63    real(wp),intent(in),dimension(dim1,     dim3) :: &
64         varIN1, & ! Input field 1
65         varIN2    ! Input field 2
66    ! OUTPUTS
67    real(wp),intent(out),dimension(dim1,dim2,dim3) :: &
68         varOUT    ! Merged output field
69    ! LOCAL VARIABLES
70    integer :: j
71   
72    varOUT(1:dim1,1:dim2,1:dim3) = 0._wp
73    DO j=1,dim2
74       where(flag(:,j,:) .eq. 1)
75          varOUT(:,j,:) = varIN2
76       endwhere
77       where(flag(:,j,:) .eq. 2)
78          varOUT(:,j,:) = varIN1
79       endwhere
80    enddo
81  end subroutine cosp_simulator_optics2D
82  subroutine cosp_simulator_optics3D(dim1,dim2,dim3,flag,varIN1,varIN2,varOUT)
83    ! INPUTS
84    integer,intent(in) :: &
85         dim1,   & ! Dimension 1 extent (Horizontal)
86         dim2,   & ! Dimension 2 extent (Subcolumn)
87         dim3      ! Dimension 3 extent (Vertical)
88    real(wp),intent(in),dimension(dim1,dim2,dim3) :: &
89         flag      ! Logical to determine the of merge var1IN and var2IN
90    real(wp),intent(in),dimension(dim1,dim2,dim3) :: &
91         varIN1, & ! Input field 1
92         varIN2    ! Input field 2
93    ! OUTPUTS
94    real(wp),intent(out),dimension(dim1,dim2,dim3) :: &
95         varOUT    ! Merged output field
96   
97    varOUT(1:dim1,1:dim2,1:dim3) = 0._wp
98   where(flag(:,:,:) .eq. 1)
99       varOUT(:,:,:) = varIN2
100    endwhere
101    where(flag(:,:,:) .eq. 2)
102       varOUT(:,:,:) = varIN1
103    endwhere
104   
105  end subroutine cosp_simulator_optics3D
106 
107  ! ##############################################################################
108  !                           MODIS_OPTICS_PARTITION
109
110  ! For the MODIS simulator, there are times when only a sinlge optical depth
111  ! profile, cloud-ice and cloud-water are provided. In this case, the optical
112  ! depth is partitioned by phase.
113  ! ##############################################################################
114  subroutine MODIS_OPTICS_PARTITION(npoints,nlev,ncolumns,cloudWater,cloudIce,waterSize, &
115                                    iceSize,tau,tauL,tauI)
116    ! INPUTS
117    INTEGER,intent(in) :: &
118         npoints,   & ! Number of horizontal gridpoints
119         nlev,      & ! Number of levels
120         ncolumns     ! Number of subcolumns
121    REAL(wp),intent(in),dimension(npoints,nlev,ncolumns) :: &
122         cloudWater, & ! Subcolumn cloud water content
123         cloudIce,   & ! Subcolumn cloud ice content
124         waterSize,  & ! Subcolumn cloud water effective radius
125         iceSize,    & ! Subcolumn cloud ice effective radius
126         tau           ! Optical thickness
127   
128    ! OUTPUTS
129    real(wp),intent(out),dimension(npoints,nlev,ncolumns) :: &
130         tauL,       & ! Partitioned liquid optical thickness.
131         tauI          ! Partitioned ice optical thickness.
132    ! LOCAL VARIABLES
133    real(wp),dimension(nlev,ncolumns) :: fracL
134    integer                           :: i
135   
136   
137    DO i=1,npoints
138       where(cloudIce(i,:, :) <= 0.)
139          fracL(:, :) = 1._wp
140       elsewhere
141          where (cloudWater(i,:, :) <= 0.)
142             fracL(:, :) = 0._wp
143          elsewhere
144             ! Geometic optics limit - tau as LWP/re  (proportional to LWC/re)
145             fracL(:, :) = (cloudWater(i,:, :)/waterSize(i,:, :)) / &
146                  (cloudWater(i,:, :)/waterSize(i,:, :) + cloudIce(i,:, :)/(ice_density * iceSize(i,:, :)) )
147          end where
148       end where
149       tauL(i,:, :) = fracL(:, :) * tau(i,:, :)
150       tauI(i,:, :) = tau(i,:, :) - tauL(i,:, :)
151    enddo
152   
153  end subroutine MODIS_OPTICS_PARTITION
154  ! ########################################################################################
155  !                                   MODIS_OPTICS
156
157  ! ########################################################################################
158  subroutine modis_optics(nPoints,nLevels,nSubCols,tauLIQ,sizeLIQ,tauICE,sizeICE,fracLIQ, g, w0)
159    ! INPUTS
160    integer, intent(in)                                      :: nPoints,nLevels,nSubCols
161    real(wp),intent(in),dimension(nPoints,nSubCols,nLevels)  :: tauLIQ, sizeLIQ, tauICE, sizeICE
162    ! OUTPUTS
163    real(wp),intent(out),dimension(nPoints,nSubCols,nLevels) :: g,w0,fracLIQ
164    ! LOCAL VARIABLES
165    real(wp), dimension(nLevels)            :: water_g, water_w0, ice_g, ice_w0,tau
166    integer :: i,j
167   
168    ! Initialize
169    g(1:nPoints,1:nSubCols,1:nLevels)  = 0._wp
170    w0(1:nPoints,1:nSubCols,1:nLevels) = 0._wp
171   
172    DO j =1,nPoints
173       DO i=1,nSubCols
174          water_g(1:nLevels)  = get_g_nir(  phaseIsLiquid, sizeLIQ(j,i,1:nLevels))
175          water_w0(1:nLevels) = get_ssa_nir(phaseIsLiquid, sizeLIQ(j,i,1:nLevels))
176          ice_g(1:nLevels)    = get_g_nir(  phaseIsIce,    sizeICE(j,i,1:nLevels))
177          ice_w0(1:nLevels)   = get_ssa_nir(phaseIsIce,    sizeICE(j,i,1:nLevels))
178         
179          ! Combine ice and water optical properties
180          tau(1:nLevels) = tauICE(j,i,1:nLevels) + tauLIQ(j,i,1:nLevels)
181          where (tau(1:nLevels) > 0)
182             w0(j,i,1:nLevels) = (tauLIQ(j,i,1:nLevels)*water_w0(1:nLevels) + tauICE(j,i,1:nLevels) *ice_w0(1:nLevels)) / &
183                                  (tau(1:nLevels))
184             g(j,i,1:nLevels) = (tauLIQ(j,i,1:nLevels)*water_g(1:nLevels)*water_w0(1:nLevels) + tauICE(j,i,1:nLevels) * &
185                  ice_g(1:nLevels) * ice_w0(1:nLevels)) / (w0(j,i,1:nLevels) * tau(1:nLevels))
186          end where
187       enddo
188    enddo
189   
190    ! Compute the total optical thickness and the proportion due to liquid in each cell
191    DO i=1,npoints
192       where(tauLIQ(i,1:nSubCols,1:nLevels) + tauICE(i,1:nSubCols,1:nLevels) > 0.)
193          fracLIQ(i,1:nSubCols,1:nLevels) = tauLIQ(i,1:nSubCols,1:nLevels)/ &
194               (tauLIQ(i,1:nSubCols,1:nLevels) + tauICE(i,1:nSubCols,1:nLevels))
195       elsewhere
196          fracLIQ(i,1:nSubCols,1:nLevels) = 0._wp
197       end  where
198    enddo
199   
200  end subroutine modis_optics
201 
202  ! ######################################################################################
203  ! SUBROUTINE lidar_optics
204  ! ######################################################################################
205  subroutine lidar_optics(npoints, ncolumns, nlev, npart, ice_type, lidar_freq, lground, &
206       q_lsliq, q_lsice, q_cvliq, q_cvice, ls_radliq, ls_radice, cv_radliq, cv_radice,   &
207       pres, presf, temp, beta_mol, betatot, tau_mol, tautot, tautot_S_liq, tautot_S_ice,&
208       betatot_ice, betatot_liq, tautot_ice, tautot_liq)
209
210    ! ####################################################################################
211    ! NOTE: Using "grav" from cosp_constants.f90, instead of grav=9.81, introduces
212    ! changes of up to 2% in atb532 adn 0.003% in parasolRefl and lidarBetaMol532.
213    ! This also results in  small changes in the joint-histogram, cfadLidarsr532.
214    ! ####################################################################################
215   
216    ! INPUTS
217    INTEGER,intent(in) :: &
218         npoints,      & ! Number of gridpoints
219         ncolumns,     & ! Number of subcolumns
220         nlev,         & ! Number of levels
221         npart,        & ! Number of cloud meteors (stratiform_liq, stratiform_ice, conv_liq, conv_ice).
222         ice_type,     & ! Ice particle shape hypothesis (0 for spheres, 1 for non-spherical)
223         lidar_freq      ! Lidar frequency (nm). Use to change between lidar platforms
224    logical,intent(in) :: &
225         lground         ! True for ground-based lidar
226    REAL(WP),intent(in),dimension(npoints,nlev) :: &
227         temp,         & ! Temperature of layer k
228         pres,         & ! Pressure at full levels
229         ls_radliq,    & ! Effective radius of LS liquid particles (meters)
230         ls_radice,    & ! Effective radius of LS ice particles (meters)
231         cv_radliq,    & ! Effective radius of CONV liquid particles (meters)
232         cv_radice       ! Effective radius of CONV ice particles (meters)
233    REAL(WP),intent(in),dimension(npoints,ncolumns,nlev) :: &
234         q_lsliq,      & ! LS sub-column liquid water mixing ratio (kg/kg)
235         q_lsice,      & ! LS sub-column ice water mixing ratio (kg/kg)
236         q_cvliq,      & ! CONV sub-column liquid water mixing ratio (kg/kg)
237         q_cvice         ! CONV sub-column ice water mixing ratio (kg/kg)
238    REAL(WP),intent(in),dimension(npoints,nlev+1) :: &
239         presf           ! Pressure at half levels
240   
241    ! OUTPUTS
242    REAL(WP),intent(out),dimension(npoints,ncolumns,nlev)       :: &
243         betatot,        & !
244         tautot            ! Optical thickess integrated from top
245    REAL(WP),intent(out),dimension(npoints,nlev) :: &
246         beta_mol,       & ! Molecular backscatter coefficient
247         tau_mol           ! Molecular optical depth
248    ! OUTPUTS (optional)
249    REAL(WP),optional,intent(out),dimension(npoints,ncolumns) :: &
250         tautot_S_liq,   & ! TOA optical depth for liquid
251         tautot_S_ice      ! TOA optical depth for ice
252    REAL(WP),optional,intent(out),dimension(npoints,ncolumns,nlev)       :: &
253         betatot_ice,    & ! Backscatter coefficient for ice particles
254         betatot_liq,    & ! Backscatter coefficient for liquid particles
255         tautot_ice,     & ! Total optical thickness of ice
256         tautot_liq        ! Total optical thickness of liq
257   
258    ! LOCAL VARIABLES
259    REAL(WP),dimension(npart)              :: rhopart
260    REAL(WP),dimension(npart,5)            :: polpart
261    REAL(WP),dimension(npoints,nlev)       :: rhoair,alpha_mol
262    REAL(WP),dimension(npoints,nlev+1)     :: zheight         
263    REAL(WP),dimension(npoints,nlev,npart) :: rad_part,kp_part,qpart,alpha_part,tau_part
264    real(wp)                               :: Cmol,rdiffm
265    logical                                :: lparasol,lphaseoptics
266    INTEGER                                :: i,k,icol,zi,zf,zinc,zoffset
267   
268    ! Local data
269    REAL(WP),PARAMETER :: rhoice     = 0.5e+03    ! Density of ice (kg/m3)
270    REAL(WP),PARAMETER :: Cmol_532nm = 6.2446e-32 ! Wavelength dependent
271    REAL(WP),PARAMETER :: Cmol_355nm = 3.2662e-31! Wavelength dependent
272    REAL(WP),PARAMETER :: rdiffm_532nm = 0.7_wp     ! Multiple scattering correction parameter
273    REAL(WP),PARAMETER :: rdiffm_355nm = 0.6_wp     ! Multiple scattering correction parameter
274    REAL(WP),PARAMETER :: Qscat      = 2.0_wp     ! Particle scattering efficiency at 532 nm
275    ! Local indicies for large-scale and convective ice and liquid
276    INTEGER,PARAMETER  :: INDX_LSLIQ  = 1
277    INTEGER,PARAMETER  :: INDX_LSICE  = 2
278    INTEGER,PARAMETER  :: INDX_CVLIQ  = 3
279    INTEGER,PARAMETER  :: INDX_CVICE  = 4
280   
281    ! Polarized optics parameterization
282    ! Polynomial coefficients for spherical liq/ice particles derived from Mie theory.
283    ! Polynomial coefficients for non spherical particles derived from a composite of
284    ! Ray-tracing theory for large particles (e.g. Noel et al., Appl. Opt., 2001)
285    ! and FDTD theory for very small particles (Yang et al., JQSRT, 2003).
286    ! We repeat the same coefficients for LS and CONV cloud to make code more readable
287    REAL(WP),PARAMETER,dimension(5) :: &
288         polpartCVLIQ  = (/ 2.6980e-8_wp,  -3.7701e-6_wp,  1.6594e-4_wp,    -0.0024_wp,    0.0626_wp/), &
289         polpartLSLIQ  = (/ 2.6980e-8_wp,  -3.7701e-6_wp,  1.6594e-4_wp,    -0.0024_wp,    0.0626_wp/), &
290         polpartCVICE0 = (/-1.0176e-8_wp,   1.7615e-6_wp, -1.0480e-4_wp,     0.0019_wp,    0.0460_wp/), &
291         polpartLSICE0 = (/-1.0176e-8_wp,   1.7615e-6_wp, -1.0480e-4_wp,     0.0019_wp,    0.0460_wp/), &
292         polpartCVICE1 = (/ 1.3615e-8_wp, -2.04206e-6_wp, 7.51799e-5_wp, 0.00078213_wp, 0.0182131_wp/), &
293         polpartLSICE1 = (/ 1.3615e-8_wp, -2.04206e-6_wp, 7.51799e-5_wp, 0.00078213_wp, 0.0182131_wp/)
294    ! ##############################################################################
295
296    ! Which LIDAR frequency are we using?
297    if (lidar_freq .eq. 355) then
298       Cmol   = Cmol_355nm
299       rdiffm = rdiffm_355nm
300    endif
301    if (lidar_freq .eq. 532) then
302       Cmol   = Cmol_532nm
303       rdiffm = rdiffm_532nm
304    endif
305   
306    ! Do we need to generate optical inputs for Parasol simulator?
307    lparasol = .false.
308    if (present(tautot_S_liq) .AND. present(tautot_S_ice)) lparasol = .true.
309   
310    ! Are optical-depths and backscatter coefficients for ice and liquid requested?
311    lphaseoptics=.false.
312    if (present(betatot_ice) .AND. present(betatot_liq) .AND. present(tautot_liq) .AND. &
313         present(tautot_ice)) lphaseoptics=.true.
314
315    ! Is this lidar spaceborne (default) or ground-based (lground=.true.)?
316    zi   = 2
317    zf   = nlev
318    zinc = 1
319    zoffset = -1
320    if (lground) then
321       zi   = nlev-1
322       zf   = 1
323       zinc = -1
324       zoffset = 1
325    endif
326   
327    ! Liquid/ice particles
328    rhopart(INDX_LSLIQ)  = rholiq
329    rhopart(INDX_LSICE)  = rhoice
330    rhopart(INDX_CVLIQ)  = rholiq
331    rhopart(INDX_CVICE)  = rhoice
332   
333    ! LS and CONV Liquid water coefficients
334    polpart(INDX_LSLIQ,1:5)  = polpartLSLIQ
335    polpart(INDX_CVLIQ,1:5)  = polpartCVLIQ
336   
337    ! LS and CONV Ice water coefficients
338    if (ice_type .eq. 0) then
339       polpart(INDX_LSICE,1:5) = polpartLSICE0
340       polpart(INDX_CVICE,1:5) = polpartCVICE0
341    endif
342    if (ice_type .eq. 1) then
343       polpart(INDX_LSICE,1:5) = polpartLSICE1
344       polpart(INDX_CVICE,1:5) = polpartCVICE1
345    endif
346   
347    ! Effective radius particles:
348    rad_part(1:npoints,1:nlev,INDX_LSLIQ)  = ls_radliq(1:npoints,1:nlev)
349    rad_part(1:npoints,1:nlev,INDX_LSICE)  = ls_radice(1:npoints,1:nlev)
350    rad_part(1:npoints,1:nlev,INDX_CVLIQ)  = cv_radliq(1:npoints,1:nlev)
351    rad_part(1:npoints,1:nlev,INDX_CVICE)  = cv_radice(1:npoints,1:nlev)   
352    rad_part(1:npoints,1:nlev,1:npart)     = MAX(rad_part(1:npoints,1:nlev,1:npart),0._wp)
353    rad_part(1:npoints,1:nlev,1:npart)     = MIN(rad_part(1:npoints,1:nlev,1:npart),70.0e-6_wp)
354   
355    ! Density (clear-sky air)
356    rhoair(1:npoints,1:nlev) = pres(1:npoints,1:nlev)/(rd*temp(1:npoints,1:nlev))
357   
358    ! Altitude at half pressure levels:
359    zheight(1:npoints,nlev+1) = 0._wp
360    DO k=nlev,1,-1
361       zheight(1:npoints,k) = zheight(1:npoints,k+1) &
362            -(presf(1:npoints,k)-presf(1:npoints,k+1))/(rhoair(1:npoints,k)*grav)
363    enddo
364   
365    ! ##############################################################################
366    ! *) Molecular alpha, beta and optical thickness
367    ! ##############################################################################
368   
369    beta_mol(1:npoints,1:nlev)  = pres(1:npoints,1:nlev)/km/temp(1:npoints,1:nlev)*Cmol
370    alpha_mol(1:npoints,1:nlev) = 8._wp*pi/3._wp * beta_mol(1:npoints,1:nlev)
371   
372    ! Optical thickness of each layer (molecular) 
373    tau_mol(1:npoints,1:nlev) = alpha_mol(1:npoints,1:nlev)*(zheight(1:npoints,1:nlev)-&
374         zheight(1:npoints,2:nlev+1))
375             
376    ! Optical thickness from TOA to layer k (molecular)
377    DO k = zi,zf,zinc
378       tau_mol(1:npoints,k) = tau_mol(1:npoints,k) + tau_mol(1:npoints,k+zoffset)
379    ENDDO   
380
381    betatot    (1:npoints,1:ncolumns,1:nlev) = spread(beta_mol(1:npoints,1:nlev), dim=2, NCOPIES=ncolumns)
382    tautot     (1:npoints,1:ncolumns,1:nlev) = spread(tau_mol (1:npoints,1:nlev), dim=2, NCOPIES=ncolumns)
383    if (lphaseoptics) then
384       betatot_liq(1:npoints,1:ncolumns,1:nlev) = betatot(1:npoints,1:ncolumns,1:nlev)
385       betatot_ice(1:npoints,1:ncolumns,1:nlev) = betatot(1:npoints,1:ncolumns,1:nlev)
386       tautot_liq (1:npoints,1:ncolumns,1:nlev) = tautot(1:npoints,1:ncolumns,1:nlev)
387       tautot_ice (1:npoints,1:ncolumns,1:nlev) = tautot(1:npoints,1:ncolumns,1:nlev)     
388    endif
389
390    ! ##############################################################################
391    ! *) Particles alpha, beta and optical thickness
392    ! ##############################################################################
393    ! Polynomials kp_lidar derived from Mie theory
394    DO i = 1, npart
395       where (rad_part(1:npoints,1:nlev,i) .gt. 0.0)
396          kp_part(1:npoints,1:nlev,i) = &
397               polpart(i,1)*(rad_part(1:npoints,1:nlev,i)*1e6)**4 &
398               + polpart(i,2)*(rad_part(1:npoints,1:nlev,i)*1e6)**3 &
399               + polpart(i,3)*(rad_part(1:npoints,1:nlev,i)*1e6)**2 &
400               + polpart(i,4)*(rad_part(1:npoints,1:nlev,i)*1e6) &
401               + polpart(i,5)
402       elsewhere
403          kp_part(1:npoints,1:nlev,i) = 0._wp
404       endwhere
405    enddo   
406
407    ! Initialize (if necessary)
408    if (lparasol) then
409       tautot_S_liq(1:npoints,1:ncolumns) = 0._wp
410       tautot_S_ice(1:npoints,1:ncolumns) = 0._wp
411    endif
412
413    ! Loop over all subcolumns
414    DO icol=1,ncolumns
415       ! ##############################################################################
416       ! Mixing ratio particles in each subcolum
417       ! ##############################################################################
418       qpart(1:npoints,1:nlev,INDX_LSLIQ) =  q_lsliq(1:npoints,icol,1:nlev)
419       qpart(1:npoints,1:nlev,INDX_LSICE)  = q_lsice(1:npoints,icol,1:nlev)
420       qpart(1:npoints,1:nlev,INDX_CVLIQ)  = q_cvliq(1:npoints,icol,1:nlev)
421       qpart(1:npoints,1:nlev,INDX_CVICE)  = q_cvice(1:npoints,icol,1:nlev)
422
423       ! ##############################################################################
424       ! Alpha and optical thickness (particles)
425       ! ##############################################################################
426       ! Alpha of particles in each subcolumn:
427       DO i = 1, npart
428          where (rad_part(1:npoints,1:nlev,i) .gt. 0.0)
429             alpha_part(1:npoints,1:nlev,i) = 3._wp/4._wp * Qscat &
430                  * rhoair(1:npoints,1:nlev) * qpart(1:npoints,1:nlev,i) &
431                  / (rhopart(i) * rad_part(1:npoints,1:nlev,i) )
432          elsewhere
433             alpha_part(1:npoints,1:nlev,i) = 0._wp
434          endwhere
435       enddo
436       
437       ! Optical thicknes
438       tau_part(1:npoints,1:nlev,1:npart) = rdiffm * alpha_part(1:npoints,1:nlev,1:npart)
439       DO i = 1, npart
440          ! Optical thickness of each layer (particles)
441          tau_part(1:npoints,1:nlev,i) = tau_part(1:npoints,1:nlev,i) &
442               & * (zheight(1:npoints,1:nlev)-zheight(1:npoints,2:nlev+1) )
443          ! Optical thickness from TOA to layer k (particles)
444          DO k=zi,zf,zinc
445             tau_part(1:npoints,k,i) = tau_part(1:npoints,k,i) + tau_part(1:npoints,k+zoffset,i)
446          enddo
447       enddo
448 
449       ! ##############################################################################
450       ! Beta and optical thickness (total=molecular + particules)
451       ! ##############################################################################
452       
453       DO i = 1, npart
454          betatot(1:npoints,icol,1:nlev) = betatot(1:npoints,icol,1:nlev) + &
455               kp_part(1:npoints,1:nlev,i)*alpha_part(1:npoints,1:nlev,i)
456          tautot(1:npoints,icol,1:nlev) = tautot(1:npoints,icol,1:nlev)  + &
457              tau_part(1:npoints,1:nlev,i)
458       ENDDO
459       
460       ! ##############################################################################
461       ! Beta and optical thickness (liquid/ice)
462       ! ##############################################################################
463       if (lphaseoptics) then
464          ! Ice
465          betatot_ice(1:npoints,icol,1:nlev) = betatot_ice(1:npoints,icol,1:nlev)+ &
466               kp_part(1:npoints,1:nlev,INDX_LSICE)*alpha_part(1:npoints,1:nlev,INDX_LSICE)+ &
467               kp_part(1:npoints,1:nlev,INDX_CVICE)*alpha_part(1:npoints,1:nlev,INDX_CVICE)
468          tautot_ice(1:npoints,icol,1:nlev) = tautot_ice(1:npoints,icol,1:nlev)  + &
469               tau_part(1:npoints,1:nlev,INDX_LSICE) + &
470               tau_part(1:npoints,1:nlev,INDX_CVICE)
471         
472          ! Liquid
473          betatot_liq(1:npoints,icol,1:nlev) = betatot_liq(1:npoints,icol,1:nlev)+ &
474               kp_part(1:npoints,1:nlev,INDX_LSLIQ)*alpha_part(1:npoints,1:nlev,INDX_LSLIQ)+ &
475               kp_part(1:npoints,1:nlev,INDX_CVLIQ)*alpha_part(1:npoints,1:nlev,INDX_CVLIQ)
476          tautot_liq(1:npoints,icol,1:nlev) = tautot_liq(1:npoints,icol,1:nlev)  + &
477               tau_part(1:npoints,1:nlev,INDX_LSLIQ) + &
478               tau_part(1:npoints,1:nlev,INDX_CVLIQ)
479       endif
480
481       ! ##############################################################################   
482       ! Optical depths used by the PARASOL simulator
483       ! ##############################################################################             
484       if (lparasol) then
485          tautot_S_liq(:,icol) = tau_part(:,nlev,1)+tau_part(:,nlev,3)
486          tautot_S_ice(:,icol) = tau_part(:,nlev,2)+tau_part(:,nlev,4)             
487       endif
488    enddo
489   
490  end subroutine lidar_optics
491
492end module cosp_optics
Note: See TracBrowser for help on using the repository browser.