| 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: 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 | |
|---|
| 47 | contains |
|---|
| 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 |
|---|
| 442 | end module cosp_optics |
|---|