Line data Source code
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 : ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33 : module 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: phaseIsLiquid,phaseIsIce,get_g_nir,get_ssa_nir
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 :
47 : contains
48 : ! ##########################################################################
49 : ! COSP_SIMULATOR_OPTICS
50 : !
51 : ! Used by: ISCCP, MISR and MODIS simulators
52 : ! ##########################################################################
53 18576 : 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 262126800 : varOUT(1:dim1,1:dim2,1:dim3) = 0._wp
71 204336 : do j=1,dim2
72 260733600 : where(flag(:,j,:) .eq. 1)
73 185760 : varOUT(:,j,:) = varIN2
74 : endwhere
75 260752176 : where(flag(:,j,:) .eq. 2)
76 : varOUT(:,j,:) = varIN1
77 : endwhere
78 : enddo
79 18576 : end subroutine cosp_simulator_optics2D
80 37152 : 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 524253600 : varOUT(1:dim1,1:dim2,1:dim3) = 0._wp
96 524253600 : where(flag(:,:,:) .eq. 1)
97 : varOUT(:,:,:) = varIN2
98 : endwhere
99 524253600 : where(flag(:,:,:) .eq. 2)
100 : varOUT(:,:,:) = varIN1
101 : endwhere
102 :
103 37152 : 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 9288 : subroutine MODIS_OPTICS_PARTITION(npoints,nlev,ncolumns,cloudWater,cloudIce,cloudSnow, &
113 9288 : waterSize,iceSize,snowSize,tau,tauL,tauI,tauS)
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 : cloudSnow, & ! Subcolumn cloud snow content
123 : waterSize, & ! Subcolumn cloud water effective radius
124 : iceSize, & ! Subcolumn cloud ice effective radius
125 : snowSize, & ! Subcolumn cloud snow 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 : tauS ! Partitioned snow optical thickness.
133 :
134 : ! LOCAL VARIABLES
135 18576 : real(wp),dimension(nlev,ncolumns) :: totalExtinction
136 : integer :: i,j,k
137 :
138 155088 : do i=1,npoints
139 124075800 : totalExtinction(:,:) = 0._wp
140 124075800 : where(waterSize(i,:,:) > 0.)
141 145800 : totalExtinction(:,:) = cloudWater(i,:,:)/waterSize(i,:,:)
142 : elsewhere
143 : totalExtinction(:,:) = 0.
144 : end where
145 :
146 : where(iceSize(i,:,:) > 0.) totalExtinction(:,:) = totalExtinction(:,:) + &
147 124075800 : cloudIce(i,:,:)/(ice_density * iceSize(i,:,:))
148 : where(snowSize(i,:,:) > 0.) totalExtinction(:,:) = totalExtinction(:,:) + &
149 124075800 : cloudSnow(i,:,:)/(ice_density * snowSize(i,:,:))
150 :
151 124075800 : where((waterSize(i,:, :) > 0.) .and. (totalExtinction(:,:) > 0.))
152 : tauL(i,:,:) = tau(i,:,:) * cloudWater(i,:,:) / &
153 : ( waterSize(i,:,:) * totalExtinction(:,:))
154 : elsewhere
155 : tauL(i,:,:) = 0.
156 : end where
157 :
158 124075800 : where( (iceSize(i,:,:) > 0.) .and. (totalExtinction(:,:) > 0.))
159 : tauI(i,:,:) = tau(i,:,:) * cloudIce(i,:,:) / &
160 : (ice_density * iceSize(i,:,:) * totalExtinction(:,:))
161 : elsewhere
162 : tauI(i,:,:) = 0.
163 : end where
164 :
165 124085088 : where( (snowSize(i,:,:) > 0.) .and. (totalExtinction(:,:) > 0.))
166 : tauS(i,:,:) = tau(i,:,:) * cloudSnow(i,:,:) / &
167 : (ice_density * snowSize(i,:,:) * totalExtinction(:,:))
168 : elsewhere
169 : tauS(i,:,:) = 0.
170 : end where
171 : enddo
172 :
173 9288 : end subroutine MODIS_OPTICS_PARTITION
174 : ! ########################################################################################
175 : ! SUBROUTINE MODIS_OPTICS
176 : ! ########################################################################################
177 9288 : subroutine modis_optics(nPoints, nLevels, nSubCols, tauLIQ, sizeLIQ, tauICE, sizeICE, &
178 9288 : tauSNOW, sizeSNOW, fracLIQ, g, w0)
179 : ! INPUTS
180 : integer, intent(in) :: nPoints,nLevels,nSubCols
181 : real(wp),intent(in),dimension(nPoints,nSubCols,nLevels) :: tauLIQ, sizeLIQ, tauICE, sizeICE,tauSNOW,sizeSNOW
182 :
183 : ! OUTPUTS
184 : real(wp),intent(out),dimension(nPoints,nSubCols,nLevels) :: g,w0,fracLIQ
185 : ! LOCAL VARIABLES
186 9288 : real(wp), dimension(nLevels) :: water_g, water_w0, ice_g, ice_w0, tau, snow_g, snow_w0
187 : integer :: i,j
188 :
189 : ! Initialize
190 131063400 : g(1:nPoints,1:nSubCols,1:nLevels) = 0._wp
191 131063400 : w0(1:nPoints,1:nSubCols,1:nLevels) = 0._wp
192 :
193 155088 : do j =1,nPoints
194 1613088 : do i=1,nSubCols
195 123930000 : water_g(1:nLevels) = get_g_nir( phaseIsLiquid, sizeLIQ(j,i,1:nLevels))
196 123930000 : water_w0(1:nLevels) = get_ssa_nir(phaseIsLiquid, sizeLIQ(j,i,1:nLevels))
197 123930000 : ice_g(1:nLevels) = get_g_nir( phaseIsIce, sizeICE(j,i,1:nLevels))
198 123930000 : ice_w0(1:nLevels) = get_ssa_nir(phaseIsIce, sizeICE(j,i,1:nLevels))
199 123930000 : snow_g(1:nLevels) = get_g_nir( phaseIsIce, sizeSNOW(j,i,1:nLevels))
200 123930000 : snow_w0(1:nLevels) = get_ssa_nir(phaseIsIce, sizeSNOW(j,i,1:nLevels))
201 :
202 : ! Combine ice, snow and water optical properties
203 123930000 : tau(1:nLevels) = tauICE(j,i,1:nLevels) + tauLIQ(j,i,1:nLevels) + tauSNOW(j,i,1:nLevels)
204 369019800 : where (tau(1:nLevels) > 0)
205 : w0(j,i,1:nLevels) = (tauLIQ(j,i,1:nLevels)*water_w0(1:nLevels) + &
206 : tauICE(j,i,1:nLevels)*ice_w0(1:nLevels) + &
207 : tauSNOW(j,i,1:nLevels)*snow_w0(1:nLevels)) / &
208 : tau(1:nLevels)
209 : g(j,i,1:nLevels) = (tauLIQ(j,i,1:nLevels)*water_g(1:nLevels)*water_w0(1:nLevels) + &
210 : tauICE(j,i,1:nLevels)*ice_g(1:nLevels)*ice_w0(1:nLevels) + &
211 : tauSNOW(j,i,1:nLevels)*snow_g(1:nLevels)*snow_w0(1:nLevels)) / &
212 : (w0(j,i,1:nLevels) * tau(1:nLevels))
213 : end where
214 : enddo
215 : enddo
216 :
217 : ! Compute the total optical thickness and the proportion due to liquid in each cell
218 155088 : do i=1,npoints
219 134874288 : where(tauLIQ(i,1:nSubCols,1:nLevels) + tauICE(i,1:nSubCols,1:nLevels) + tauSNOW(i,1:nSubCols,1:nLevels) > 0.)
220 : fracLIQ(i,1:nSubCols,1:nLevels) = tauLIQ(i,1:nSubCols,1:nLevels)/ &
221 : (tauLIQ(i,1:nSubCols,1:nLevels) + tauICE(i,1:nSubCols,1:nLevels) + tauSNOW(i,1:nSubCols,1:nLevels) )
222 : elsewhere
223 145800 : fracLIQ(i,1:nSubCols,1:nLevels) = 0._wp
224 : end where
225 : enddo
226 :
227 9288 : end subroutine modis_optics
228 :
229 : ! ######################################################################################
230 : ! SUBROUTINE lidar_optics
231 : ! ######################################################################################
232 18576 : subroutine lidar_optics(npoints,ncolumns,nlev,npart,ice_type,q_lsliq,q_lsice,q_cvliq, &
233 18576 : q_cvice,q_lssnow,ls_radliq,ls_radice,cv_radliq,cv_radice,ls_radsnow, &
234 9288 : pres,presf,temp,beta_mol,betatot,tau_mol,tautot, &
235 18576 : tautot_S_liq,tautot_S_ice,betatot_ice,betatot_liq, &
236 18576 : tautot_ice,tautot_liq)
237 :
238 : ! ####################################################################################
239 : ! NOTE: Using "grav" from cosp_constants.f90, instead of grav=9.81, introduces
240 : ! changes of up to 2% in atb532 adn 0.003% in parasolRefl and lidarBetaMol532.
241 : ! This also results in small changes in the joint-histogram, cfadLidarsr532.
242 : ! ####################################################################################
243 :
244 : ! INPUTS
245 : INTEGER,intent(in) :: &
246 : npoints, & ! Number of gridpoints
247 : ncolumns, & ! Number of subcolumns
248 : nlev, & ! Number of levels
249 : npart, & ! Number of cloud meteors (stratiform_liq, stratiform_ice, conv_liq, conv_ice).
250 : ice_type ! Ice particle shape hypothesis (0 for spheres, 1 for non-spherical)
251 : REAL(WP),intent(in),dimension(npoints,nlev) :: &
252 : temp, & ! Temperature of layer k
253 : pres, & ! Pressure at full levels
254 : ls_radliq, & ! Effective radius of LS liquid particles (meters)
255 : ls_radice, & ! Effective radius of LS ice particles (meters)
256 : cv_radliq, & ! Effective radius of CONV liquid particles (meters)
257 : cv_radice ! Effective radius of CONV ice particles (meters)
258 : REAL(WP),intent(inout),dimension(npoints,nlev) :: &
259 : ls_radsnow ! Effective radius of LS snow particles (meters)
260 : REAL(WP),intent(in),dimension(npoints,ncolumns,nlev) :: &
261 : q_lsliq, & ! LS sub-column liquid water mixing ratio (kg/kg)
262 : q_lsice, & ! LS sub-column ice water mixing ratio (kg/kg)
263 : q_cvliq, & ! CONV sub-column liquid water mixing ratio (kg/kg)
264 : q_cvice, & ! CONV sub-column ice water mixing ratio (kg/kg)
265 : q_lssnow ! LS sub-column snow mixing ratio (kg/kg)
266 : REAL(WP),intent(in),dimension(npoints,nlev+1) :: &
267 : presf ! Pressure at half levels
268 :
269 : ! OUTPUTS
270 : REAL(WP),intent(out),dimension(npoints,ncolumns,nlev) :: &
271 : betatot, & !
272 : tautot ! Optical thickess integrated from top
273 : REAL(WP),optional,intent(out),dimension(npoints,ncolumns,nlev) :: &
274 : betatot_ice, & ! Backscatter coefficient for ice particles
275 : betatot_liq, & ! Backscatter coefficient for liquid particles
276 : tautot_ice, & ! Total optical thickness of ice
277 : tautot_liq ! Total optical thickness of liq
278 : REAL(WP),intent(out),dimension(npoints,nlev) :: &
279 : beta_mol, & ! Molecular backscatter coefficient
280 : tau_mol ! Molecular optical depth
281 : REAL(WP),intent(out),dimension(npoints,ncolumns) :: &
282 : tautot_S_liq, & ! TOA optical depth for liquid
283 : tautot_S_ice ! TOA optical depth for ice
284 :
285 : ! LOCAL VARIABLES
286 18576 : REAL(WP),dimension(npart) :: rhopart
287 18576 : REAL(WP),dimension(npart,5) :: polpart
288 18576 : REAL(WP),dimension(npoints,nlev) :: rhoair,alpha_mol
289 18576 : REAL(WP),dimension(npoints,nlev+1) :: zheight
290 18576 : REAL(WP),dimension(npoints,nlev,npart) :: rad_part,kp_part,qpart,alpha_part,tau_part
291 :
292 : INTEGER :: i,k,icol
293 :
294 : ! Local data
295 : REAL(WP),PARAMETER :: rhoice = 0.5e+03 ! Density of ice (kg/m3)
296 : REAL(WP),PARAMETER :: Cmol = 6.2446e-32 ! Wavelength dependent
297 : REAL(WP),PARAMETER :: rdiffm = 0.7_wp ! Multiple scattering correction parameter
298 : REAL(WP),PARAMETER :: Qscat = 2.0_wp ! Particle scattering efficiency at 532 nm
299 : ! Local indicies for large-scale and convective ice and liquid
300 : INTEGER,PARAMETER :: INDX_LSLIQ = 1
301 : INTEGER,PARAMETER :: INDX_LSICE = 2
302 : INTEGER,PARAMETER :: INDX_CVLIQ = 3
303 : INTEGER,PARAMETER :: INDX_CVICE = 4
304 : INTEGER,PARAMETER :: INDX_LSSNOW = 5
305 :
306 : ! Polarized optics parameterization
307 : ! Polynomial coefficients for spherical liq/ice particles derived from Mie theory.
308 : ! Polynomial coefficients for non spherical particles derived from a composite of
309 : ! Ray-tracing theory for large particles (e.g. Noel et al., Appl. Opt., 2001)
310 : ! and FDTD theory for very small particles (Yang et al., JQSRT, 2003).
311 : ! We repeat the same coefficients for LS and CONV cloud to make code more readable
312 : REAL(WP),PARAMETER,dimension(5) :: &
313 : polpartCVLIQ = (/ 2.6980e-8_wp, -3.7701e-6_wp, 1.6594e-4_wp, -0.0024_wp, 0.0626_wp/), &
314 : polpartLSLIQ = (/ 2.6980e-8_wp, -3.7701e-6_wp, 1.6594e-4_wp, -0.0024_wp, 0.0626_wp/), &
315 : polpartCVICE0 = (/-1.0176e-8_wp, 1.7615e-6_wp, -1.0480e-4_wp, 0.0019_wp, 0.0460_wp/), &
316 : polpartLSICE0 = (/-1.0176e-8_wp, 1.7615e-6_wp, -1.0480e-4_wp, 0.0019_wp, 0.0460_wp/), &
317 : polpartCVICE1 = (/ 1.3615e-8_wp, -2.04206e-6_wp, 7.51799e-5_wp, 0.00078213_wp, 0.0182131_wp/), &
318 : polpartLSICE1 = (/ 1.3615e-8_wp, -2.04206e-6_wp, 7.51799e-5_wp, 0.00078213_wp, 0.0182131_wp/), &
319 : polpartLSSNOW = (/ 1.3615e-8_wp, -2.04206e-6_wp, 7.51799e-5_wp, 0.00078213_wp, 0.0182131_wp/)
320 : ! ##############################################################################
321 :
322 : ! Liquid/ice particles
323 9288 : rhopart(INDX_LSLIQ) = rholiq
324 9288 : rhopart(INDX_LSICE) = rhoice
325 9288 : rhopart(INDX_CVLIQ) = rholiq
326 9288 : rhopart(INDX_CVICE) = rhoice
327 9288 : rhopart(INDX_LSSNOW) = rhoice/2._wp
328 :
329 : ! LS and CONV Liquid water coefficients
330 55728 : polpart(INDX_LSLIQ,1:5) = polpartLSLIQ
331 55728 : polpart(INDX_CVLIQ,1:5) = polpartCVLIQ
332 55728 : polpart(INDX_LSSNOW,1:5) = polpartLSSNOW
333 : ! LS and CONV Ice water coefficients
334 9288 : if (ice_type .eq. 0) then
335 55728 : polpart(INDX_LSICE,1:5) = polpartLSICE0
336 55728 : polpart(INDX_CVICE,1:5) = polpartCVICE0
337 : endif
338 9288 : if (ice_type .eq. 1) then
339 0 : polpart(INDX_LSICE,1:5) = polpartLSICE1
340 0 : polpart(INDX_CVICE,1:5) = polpartCVICE1
341 : endif
342 :
343 : ! Effective radius particles:
344 13036680 : rad_part(1:npoints,1:nlev,INDX_LSLIQ) = ls_radliq(1:npoints,1:nlev)
345 13036680 : rad_part(1:npoints,1:nlev,INDX_LSICE) = ls_radice(1:npoints,1:nlev)
346 13036680 : rad_part(1:npoints,1:nlev,INDX_CVLIQ) = cv_radliq(1:npoints,1:nlev)
347 13036680 : rad_part(1:npoints,1:nlev,INDX_CVICE) = cv_radice(1:npoints,1:nlev)
348 13036680 : rad_part(1:npoints,1:nlev,INDX_LSSNOW) = ls_radsnow(1:npoints,1:nlev)
349 65192688 : rad_part(:,:,:) = MAX(rad_part(:,:,:),0._wp)
350 65192688 : rad_part(:,:,:) = MIN(rad_part(:,:,:),70.0e-6_wp)
351 13036680 : ls_radsnow(:,:) = MAX(ls_radsnow(:,:),0._wp)
352 13036680 : ls_radsnow(:,:) = MIN(ls_radsnow(:,:),1000.e-6_wp)
353 :
354 : ! Density (clear-sky air)
355 13036680 : rhoair(1:npoints,1:nlev) = pres(1:npoints,1:nlev)/(rd*temp(1:npoints,1:nlev))
356 :
357 : ! Altitude at half pressure levels:
358 155088 : zheight(1:npoints,nlev+1) = 0._wp
359 789480 : do k=nlev,1,-1
360 780192 : zheight(1:npoints,k) = zheight(1:npoints,k+1) &
361 13816872 : -(presf(1:npoints,k)-presf(1:npoints,k+1))/(rhoair(1:npoints,k)*grav)
362 : enddo
363 :
364 : ! ##############################################################################
365 : ! *) Molecular alpha, beta and optical thickness
366 : ! ##############################################################################
367 :
368 13036680 : beta_mol(1:npoints,1:nlev) = pres(1:npoints,1:nlev)/km/temp(1:npoints,1:nlev)*Cmol
369 13036680 : alpha_mol(1:npoints,1:nlev) = 8._wp*pi/3._wp * beta_mol(1:npoints,1:nlev)
370 :
371 : ! Optical thickness of each layer (molecular)
372 9288 : tau_mol(1:npoints,1:nlev) = alpha_mol(1:npoints,1:nlev)*(zheight(1:npoints,1:nlev)-&
373 13045968 : zheight(1:npoints,2:nlev+1))
374 :
375 : ! Optical thickness from TOA to layer k (molecular)
376 780192 : DO k = 2,nlev
377 12881592 : tau_mol(1:npoints,k) = tau_mol(1:npoints,k) + tau_mol(1:npoints,k-1)
378 : ENDDO
379 :
380 9288 : betatot (1:npoints,1:ncolumns,1:nlev) = spread(beta_mol(1:npoints,1:nlev), dim=2, NCOPIES=ncolumns)
381 9288 : tautot (1:npoints,1:ncolumns,1:nlev) = spread(tau_mol (1:npoints,1:nlev), dim=2, NCOPIES=ncolumns)
382 131072688 : betatot_liq(1:npoints,1:ncolumns,1:nlev) = betatot(1:npoints,1:ncolumns,1:nlev)
383 131072688 : betatot_ice(1:npoints,1:ncolumns,1:nlev) = betatot(1:npoints,1:ncolumns,1:nlev)
384 131072688 : tautot_liq (1:npoints,1:ncolumns,1:nlev) = tautot(1:npoints,1:ncolumns,1:nlev)
385 131072688 : tautot_ice (1:npoints,1:ncolumns,1:nlev) = tautot(1:npoints,1:ncolumns,1:nlev)
386 :
387 : ! ##############################################################################
388 : ! *) Particles alpha, beta and optical thickness
389 : ! ##############################################################################
390 : ! Polynomials kp_lidar derived from Mie theory
391 55728 : do i = 1, npart
392 65192688 : where (rad_part(1:npoints,1:nlev,i) .gt. 0.0)
393 : kp_part(1:npoints,1:nlev,i) = &
394 46440 : polpart(i,1)*(rad_part(1:npoints,1:nlev,i)*1e6)**4 &
395 : + polpart(i,2)*(rad_part(1:npoints,1:nlev,i)*1e6)**3 &
396 : + polpart(i,3)*(rad_part(1:npoints,1:nlev,i)*1e6)**2 &
397 : + polpart(i,4)*(rad_part(1:npoints,1:nlev,i)*1e6) &
398 : + polpart(i,5)
399 : elsewhere
400 : kp_part(1:npoints,1:nlev,i) = 0._wp
401 : endwhere
402 : enddo
403 :
404 : ! Loop over all subcolumns
405 1560168 : tautot_S_liq(1:npoints,1:ncolumns) = 0._wp
406 1560168 : tautot_S_ice(1:npoints,1:ncolumns) = 0._wp
407 102168 : do icol=1,ncolumns
408 : ! ##############################################################################
409 : ! Mixing ratio particles in each subcolum
410 : ! ##############################################################################
411 130366800 : qpart(1:npoints,1:nlev,INDX_LSLIQ) = q_lsliq(1:npoints,icol,1:nlev)
412 130366800 : qpart(1:npoints,1:nlev,INDX_LSICE) = q_lsice(1:npoints,icol,1:nlev)
413 130366800 : qpart(1:npoints,1:nlev,INDX_CVLIQ) = q_cvliq(1:npoints,icol,1:nlev)
414 130366800 : qpart(1:npoints,1:nlev,INDX_CVICE) = q_cvice(1:npoints,icol,1:nlev)
415 130366800 : qpart(1:npoints,1:nlev,INDX_LSSNOW) = q_lssnow(1:npoints,icol,1:nlev)
416 :
417 : ! ##############################################################################
418 : ! Alpha and optical thickness (particles)
419 : ! ##############################################################################
420 : ! Alpha of particles in each subcolumn:
421 464400 : do i = 1, npart-1
422 521560080 : where (rad_part(1:npoints,1:nlev,i) .gt. 0.0)
423 : alpha_part(1:npoints,1:nlev,i) = 3._wp/4._wp * Qscat &
424 : * rhoair(1:npoints,1:nlev) * qpart(1:npoints,1:nlev,i) &
425 371520 : / (rhopart(i) * rad_part(1:npoints,1:nlev,i) )
426 : elsewhere
427 : alpha_part(1:npoints,1:nlev,i) = 0._wp
428 : endwhere
429 : enddo
430 130366800 : where ( ls_radsnow(:,:).gt.0.0)
431 : alpha_part(:,:,5) = 3.0/4.0 * Qscat * rhoair(:,:) * qpart(:,:,5)/(rhopart(5) * ls_radsnow(:,:) )
432 : elsewhere
433 : alpha_part(:,:,5) = 0.
434 : endwhere
435 :
436 : ! Optical thicknes
437 651926880 : tau_part(1:npoints,1:nlev,1:npart) = rdiffm * alpha_part(1:npoints,1:nlev,1:npart)
438 557280 : do i = 1, npart
439 : ! Optical thickness of each layer (particles)
440 464400 : tau_part(1:npoints,1:nlev,i) = tau_part(1:npoints,1:nlev,i) &
441 652298400 : & * (zheight(1:npoints,1:nlev)-zheight(1:npoints,2:nlev+1) )
442 : ! Optical thickness from TOA to layer k (particles)
443 39102480 : do k=2,nlev
444 644079600 : tau_part(1:npoints,k,i) = tau_part(1:npoints,k,i) + tau_part(1:npoints,k-1,i)
445 : enddo
446 : enddo
447 :
448 : ! ##############################################################################
449 : ! Beta and optical thickness (total=molecular + particules)
450 : ! ##############################################################################
451 :
452 557280 : DO i = 1, npart
453 : betatot(1:npoints,icol,1:nlev) = betatot(1:npoints,icol,1:nlev) + &
454 651834000 : kp_part(1:npoints,1:nlev,i)*alpha_part(1:npoints,1:nlev,i)
455 : tautot(1:npoints,icol,1:nlev) = tautot(1:npoints,icol,1:nlev) + &
456 651926880 : tau_part(1:npoints,1:nlev,i)
457 : ENDDO
458 :
459 : ! ##############################################################################
460 : ! Beta and optical thickness (liquid/ice)
461 : ! ##############################################################################
462 : ! Ice
463 371520 : betatot_ice(1:npoints,icol,1:nlev) = betatot_ice(1:npoints,icol,1:nlev)+ &
464 92880 : kp_part(1:npoints,1:nlev,INDX_LSICE)*alpha_part(1:npoints,1:nlev,INDX_LSICE)+ &
465 : kp_part(1:npoints,1:nlev,INDX_CVICE)*alpha_part(1:npoints,1:nlev,INDX_CVICE)+ &
466 130738320 : kp_part(1:npoints,1:nlev,INDX_LSSNOW)*alpha_part(1:npoints,1:nlev,INDX_LSSNOW)
467 371520 : tautot_ice(1:npoints,icol,1:nlev) = tautot_ice(1:npoints,icol,1:nlev) + &
468 92880 : tau_part(1:npoints,1:nlev,INDX_LSICE) + &
469 : tau_part(1:npoints,1:nlev,INDX_CVICE) + &
470 130738320 : tau_part(1:npoints,1:nlev,INDX_LSSNOW)
471 :
472 : ! Liquid
473 371520 : betatot_liq(1:npoints,icol,1:nlev) = betatot_liq(1:npoints,icol,1:nlev)+ &
474 92880 : kp_part(1:npoints,1:nlev,INDX_LSLIQ)*alpha_part(1:npoints,1:nlev,INDX_LSLIQ)+ &
475 130738320 : kp_part(1:npoints,1:nlev,INDX_CVLIQ)*alpha_part(1:npoints,1:nlev,INDX_CVLIQ)
476 371520 : tautot_liq(1:npoints,icol,1:nlev) = tautot_liq(1:npoints,icol,1:nlev) + &
477 92880 : tau_part(1:npoints,1:nlev,INDX_LSLIQ) + &
478 130738320 : tau_part(1:npoints,1:nlev,INDX_CVLIQ)
479 :
480 : ! ##############################################################################
481 : ! Optical depths used by the PARASOL simulator
482 : ! ##############################################################################
483 1550880 : tautot_S_liq(:,icol) = tau_part(:,nlev,1)+tau_part(:,nlev,3)
484 1560168 : tautot_S_ice(:,icol) = tau_part(:,nlev,2)+tau_part(:,nlev,4)+tau_part(:,nlev,5)
485 : enddo
486 :
487 37152 : end subroutine lidar_optics
488 :
489 :
490 : end module cosp_optics
|
