| [3491] | 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 | ! May 2009: Robert Pincus - Initial version |
|---|
| 31 | ! June 2009: Steve Platnick and Robert Pincus - Simple radiative transfer for size |
|---|
| 32 | ! retrievals |
|---|
| 33 | ! August 2009: Robert Pincus - Consistency and bug fixes suggested by Rick Hemler (GFDL) |
|---|
| 34 | ! November 2009: Robert Pincus - Bux fixes and speed-ups after experience with Rick Hemler |
|---|
| 35 | ! using AM2 (GFDL) |
|---|
| 36 | ! January 2010: Robert Pincus - Added high, middle, low cloud fractions |
|---|
| 37 | ! May 2015: Dustin Swales - Modified for COSPv2.0 |
|---|
| 38 | ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|---|
| 39 | ! |
|---|
| 40 | ! Notes on using the MODIS simulator: |
|---|
| 41 | ! *) You may provide either layer-by-layer values of optical thickness at 0.67 and 2.1 |
|---|
| 42 | ! microns, or optical thickness at 0.67 microns and ice- and liquid-water contents |
|---|
| 43 | ! (in consistent units of your choosing) |
|---|
| 44 | ! *) Required input also includes the optical thickness and cloud top pressure |
|---|
| 45 | ! derived from the ISCCP simulator run with parameter top_height = 1. |
|---|
| 46 | ! *) Cloud particle sizes are specified as radii, measured in meters, though within the |
|---|
| 47 | ! module we use units of microns. Where particle sizes are outside the bounds used in |
|---|
| 48 | ! the MODIS retrieval libraries (parameters re_water_min, re_ice_min, etc.) the |
|---|
| 49 | ! simulator returns missing values (re_fill) |
|---|
| 50 | ! |
|---|
| 51 | ! When error conditions are encountered this code calls the function complain_and_die, |
|---|
| 52 | ! supplied at the bottom of this module. Users probably want to replace this with |
|---|
| 53 | ! something more graceful. |
|---|
| 54 | ! |
|---|
| 55 | module mod_modis_sim |
|---|
| 56 | USE MOD_COSP_CONFIG, only: R_UNDEF,modis_histTau,modis_histPres,numMODISTauBins, & |
|---|
| 57 | numMODISPresBins,numMODISReffIceBins,numMODISReffLiqBins, & |
|---|
| 58 | modis_histReffIce,modis_histReffLiq |
|---|
| 59 | USE COSP_KINDS, ONLY: wp |
|---|
| 60 | use MOD_COSP_STATS, ONLY: hist2D |
|---|
| 61 | |
|---|
| 62 | implicit none |
|---|
| 63 | ! ########################################################################## |
|---|
| 64 | ! Retrieval parameters |
|---|
| 65 | integer, parameter :: & |
|---|
| 66 | num_trial_res = 15 ! Increase to make the linear pseudo-retrieval of size more accurate |
|---|
| 67 | |
|---|
| 68 | real(wp) :: & |
|---|
| 69 | min_OpticalThickness, & ! Minimum detectable optical thickness |
|---|
| 70 | CO2Slicing_PressureLimit, & ! Cloud with higher pressures use thermal methods, units Pa |
|---|
| 71 | CO2Slicing_TauLimit, & ! How deep into the cloud does CO2 slicing see? |
|---|
| 72 | phase_TauLimit, & ! How deep into the cloud does the phase detection see? |
|---|
| 73 | size_TauLimit, & ! Depth of the re retreivals |
|---|
| 74 | phaseDiscrimination_Threshold, & ! What fraction of total extincton needs to be in a single |
|---|
| 75 | ! category to make phase discrim. work? |
|---|
| 76 | re_fill, & ! |
|---|
| 77 | re_water_min, & ! Minimum effective radius (liquid) |
|---|
| 78 | re_water_max, & ! Maximum effective radius (liquid) |
|---|
| 79 | re_ice_min, & ! Minimum effective radius (ice) |
|---|
| 80 | re_ice_max, & ! Minimum effective radius (ice) |
|---|
| 81 | highCloudPressureLimit, & ! High cloud pressure limit (Pa) |
|---|
| 82 | lowCloudPressureLimit ! Low cloud pressure limit (Pa) |
|---|
| 83 | integer :: & |
|---|
| 84 | phaseIsNone, & ! |
|---|
| 85 | phaseIsLiquid, & ! |
|---|
| 86 | phaseIsIce, & ! |
|---|
| 87 | phaseIsUndetermined ! |
|---|
| 88 | |
|---|
| 89 | real(wp),dimension(num_trial_res) :: & |
|---|
| 90 | trial_re_w, & ! Near-IR optical params vs size for retrieval scheme (liquid) |
|---|
| 91 | trial_re_i ! Near-IR optical params vs size for retrieval scheme (ice) |
|---|
| 92 | real(wp),dimension(num_trial_res) :: & |
|---|
| 93 | g_w, & ! Assymettry parameter for size retrieval (liquid) |
|---|
| 94 | g_i, & ! Assymettry parameter for size retrieval (ice) |
|---|
| 95 | w0_w, & ! Single-scattering albedo for size retrieval (liquid) |
|---|
| 96 | w0_i ! Single-scattering albedo for size retrieval (ice) |
|---|
| 97 | ! Algorithmic parameters |
|---|
| 98 | real(wp),parameter :: & |
|---|
| 99 | ice_density = 0.93_wp ! Liquid density is 1. |
|---|
| 100 | |
|---|
| 101 | contains |
|---|
| 102 | ! ######################################################################################## |
|---|
| 103 | ! MODIS simulator using specified liquid and ice optical thickness in each layer |
|---|
| 104 | ! |
|---|
| 105 | ! Note: this simulator operates on all points; to match MODIS itself night-time |
|---|
| 106 | ! points should be excluded |
|---|
| 107 | ! |
|---|
| 108 | ! Note: the simulator requires as input the optical thickness and cloud top pressure |
|---|
| 109 | ! derived from the ISCCP simulator run with parameter top_height = 1. |
|---|
| 110 | ! If cloud top pressure is higher than about 700 mb, MODIS can't use CO2 slicing |
|---|
| 111 | ! and reverts to a thermal algorithm much like ISCCP's. Rather than replicate that |
|---|
| 112 | ! alogrithm in this simulator we simply report the values from the ISCCP simulator. |
|---|
| 113 | ! ######################################################################################## |
|---|
| 114 | subroutine modis_subcolumn(nSubCols, nLevels, pressureLevels, optical_thickness, & |
|---|
| 115 | tauLiquidFraction, g, w0,isccpCloudTopPressure, & |
|---|
| 116 | retrievedPhase, retrievedCloudTopPressure, & |
|---|
| 117 | retrievedTau, retrievedSize) |
|---|
| 118 | |
|---|
| 119 | ! INPUTS |
|---|
| 120 | integer,intent(in) :: & |
|---|
| 121 | nSubCols, & ! Number of subcolumns |
|---|
| 122 | nLevels ! Number of levels |
|---|
| 123 | real(wp),dimension(nLevels+1),intent(in) :: & |
|---|
| 124 | pressureLevels ! Gridmean pressure at layer edges (Pa) |
|---|
| 125 | real(wp),dimension(nSubCols,nLevels),intent(in) :: & |
|---|
| 126 | optical_thickness, & ! Subcolumn optical thickness @ 0.67 microns. |
|---|
| 127 | tauLiquidFraction, & ! Liquid water fraction |
|---|
| 128 | g, & ! Subcolumn assymetry parameter |
|---|
| 129 | w0 ! Subcolumn single-scattering albedo |
|---|
| 130 | real(wp),dimension(nSubCols),intent(in) :: & |
|---|
| 131 | isccpCloudTopPressure ! ISCCP retrieved cloud top pressure (Pa) |
|---|
| 132 | |
|---|
| 133 | ! OUTPUTS |
|---|
| 134 | integer, dimension(nSubCols), intent(inout) :: & |
|---|
| 135 | retrievedPhase ! MODIS retrieved phase (liquid/ice/other) |
|---|
| 136 | real(wp),dimension(nSubCols), intent(inout) :: & |
|---|
| 137 | retrievedCloudTopPressure, & ! MODIS retrieved CTP (Pa) |
|---|
| 138 | retrievedTau, & ! MODIS retrieved optical depth (unitless) |
|---|
| 139 | retrievedSize ! MODIS retrieved particle size (microns) |
|---|
| 140 | |
|---|
| 141 | ! LOCAL VARIABLES |
|---|
| 142 | logical, dimension(nSubCols) :: & |
|---|
| 143 | cloudMask |
|---|
| 144 | real(wp) :: & |
|---|
| 145 | integratedLiquidFraction, & |
|---|
| 146 | obs_Refl_nir |
|---|
| 147 | real(wp),dimension(num_trial_res) :: & |
|---|
| 148 | predicted_Refl_nir |
|---|
| 149 | integer :: & |
|---|
| 150 | i |
|---|
| 151 | |
|---|
| 152 | ! ######################################################################################## |
|---|
| 153 | ! Optical depth retrieval |
|---|
| 154 | ! This is simply a sum over the optical thickness in each layer. |
|---|
| 155 | ! It should agree with the ISCCP values after min values have been excluded. |
|---|
| 156 | ! ######################################################################################## |
|---|
| 157 | retrievedTau(1:nSubCols) = sum(optical_thickness(1:nSubCols,1:nLevels), dim = 2) |
|---|
| 158 | |
|---|
| 159 | ! ######################################################################################## |
|---|
| 160 | ! Cloud detection |
|---|
| 161 | ! does optical thickness exceed detection threshold? |
|---|
| 162 | ! ######################################################################################## |
|---|
| 163 | cloudMask = retrievedTau(1:nSubCols) >= min_OpticalThickness |
|---|
| 164 | |
|---|
| 165 | do i = 1, nSubCols |
|---|
| 166 | if(cloudMask(i)) then |
|---|
| 167 | ! ################################################################################## |
|---|
| 168 | ! Cloud top pressure determination |
|---|
| 169 | ! MODIS uses CO2 slicing for clouds with tops above about 700 mb and thermal |
|---|
| 170 | ! methods for clouds lower than that. For CO2 slicing we report the optical-depth |
|---|
| 171 | ! weighted pressure, integrating to a specified optical depth. |
|---|
| 172 | ! This assumes linear variation in p between levels. Linear in ln(p) is probably |
|---|
| 173 | ! better, though we'd need to deal with the lowest pressure gracefully. |
|---|
| 174 | ! ################################################################################## |
|---|
| 175 | retrievedCloudTopPressure(i) = cloud_top_pressure(nLevels,(/ 0._wp, optical_thickness(i,1:nLevels) /), & |
|---|
| 176 | pressureLevels(1:nLevels),CO2Slicing_TauLimit) |
|---|
| 177 | |
|---|
| 178 | ! ################################################################################## |
|---|
| 179 | ! Phase determination |
|---|
| 180 | ! Determine fraction of total tau that's liquid when ice and water contribute about |
|---|
| 181 | ! equally to the extinction we can't tell what the phase is. |
|---|
| 182 | ! ################################################################################## |
|---|
| 183 | integratedLiquidFraction = weight_by_extinction(nLevels,optical_thickness(i,1:nLevels), & |
|---|
| 184 | tauLiquidFraction(i, 1:nLevels), & |
|---|
| 185 | phase_TauLimit) |
|---|
| 186 | if(integratedLiquidFraction >= phaseDiscrimination_Threshold) then |
|---|
| 187 | retrievedPhase(i) = phaseIsLiquid |
|---|
| 188 | else if (integratedLiquidFraction <= 1._wp- phaseDiscrimination_Threshold) then |
|---|
| 189 | retrievedPhase(i) = phaseIsIce |
|---|
| 190 | else |
|---|
| 191 | retrievedPhase(i) = phaseIsUndetermined |
|---|
| 192 | end if |
|---|
| 193 | |
|---|
| 194 | ! ################################################################################## |
|---|
| 195 | ! Size determination |
|---|
| 196 | ! ################################################################################## |
|---|
| 197 | |
|---|
| 198 | ! Compute observed reflectance |
|---|
| 199 | obs_Refl_nir = compute_toa_reflectace(nLevels,optical_thickness(i,1:nLevels), g(i,1:nLevels), w0(i,1:nLevels)) |
|---|
| 200 | |
|---|
| 201 | ! Compute predicted reflectance |
|---|
| 202 | if(any(retrievedPhase(i) == (/ phaseIsLiquid, phaseIsUndetermined, phaseIsIce /))) then |
|---|
| 203 | if (retrievedPhase(i) == phaseIsLiquid .OR. retrievedPhase(i) == phaseIsUndetermined) then |
|---|
| 204 | predicted_Refl_nir(1:num_trial_res) = two_stream_reflectance(retrievedTau(i), & |
|---|
| 205 | g_w(1:num_trial_res), w0_w(1:num_trial_res)) |
|---|
| 206 | retrievedSize(i) = 1.0e-06_wp*interpolate_to_min(trial_re_w(1:num_trial_res), & |
|---|
| 207 | predicted_Refl_nir(1:num_trial_res), obs_Refl_nir) |
|---|
| 208 | else |
|---|
| 209 | predicted_Refl_nir(1:num_trial_res) = two_stream_reflectance(retrievedTau(i), & |
|---|
| 210 | g_i(1:num_trial_res), w0_i(1:num_trial_res)) |
|---|
| 211 | retrievedSize(i) = 1.0e-06_wp*interpolate_to_min(trial_re_i(1:num_trial_res), & |
|---|
| 212 | predicted_Refl_nir(1:num_trial_res), obs_Refl_nir) |
|---|
| 213 | endif |
|---|
| 214 | else |
|---|
| 215 | retrievedSize(i) = re_fill |
|---|
| 216 | endif |
|---|
| 217 | else |
|---|
| 218 | ! Values when we don't think there's a cloud. |
|---|
| 219 | retrievedCloudTopPressure(i) = R_UNDEF |
|---|
| 220 | retrievedPhase(i) = phaseIsNone |
|---|
| 221 | retrievedSize(i) = R_UNDEF |
|---|
| 222 | retrievedTau(i) = R_UNDEF |
|---|
| 223 | end if |
|---|
| 224 | end do |
|---|
| 225 | where((retrievedSize(1:nSubCols) < 0.).and.(retrievedSize(1:nSubCols) /= R_UNDEF)) & |
|---|
| 226 | retrievedSize(1:nSubCols) = 1.0e-06_wp*re_fill |
|---|
| 227 | |
|---|
| 228 | ! We use the ISCCP-derived CTP for low clouds, since the ISCCP simulator ICARUS |
|---|
| 229 | ! mimics what MODIS does to first order. |
|---|
| 230 | ! Of course, ISCCP cloud top pressures are in mb. |
|---|
| 231 | where(cloudMask(1:nSubCols) .and. retrievedCloudTopPressure(1:nSubCols) > CO2Slicing_PressureLimit) & |
|---|
| 232 | retrievedCloudTopPressure(1:nSubCols) = isccpCloudTopPressure! * 100._wp |
|---|
| 233 | |
|---|
| 234 | end subroutine modis_subcolumn |
|---|
| 235 | |
|---|
| 236 | ! ######################################################################################## |
|---|
| 237 | subroutine modis_column(nPoints,nSubCols,phase, cloud_top_pressure, optical_thickness, particle_size, & |
|---|
| 238 | Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, & |
|---|
| 239 | Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, & |
|---|
| 240 | Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, & |
|---|
| 241 | Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10,& |
|---|
| 242 | Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, Cloud_Top_Pressure_Total_Mean, & |
|---|
| 243 | Liquid_Water_Path_Mean, Ice_Water_Path_Mean, & |
|---|
| 244 | Optical_Thickness_vs_Cloud_Top_Pressure,Optical_Thickness_vs_ReffIce,Optical_Thickness_vs_ReffLiq) |
|---|
| 245 | |
|---|
| 246 | ! INPUTS |
|---|
| 247 | integer,intent(in) :: & |
|---|
| 248 | nPoints, & ! Number of horizontal gridpoints |
|---|
| 249 | nSubCols ! Number of subcolumns |
|---|
| 250 | integer,intent(in), dimension(nPoints, nSubCols) :: & |
|---|
| 251 | phase |
|---|
| 252 | real(wp),intent(in),dimension(nPoints, nSubCols) :: & |
|---|
| 253 | cloud_top_pressure, & |
|---|
| 254 | optical_thickness, & |
|---|
| 255 | particle_size |
|---|
| 256 | |
|---|
| 257 | ! OUTPUTS |
|---|
| 258 | real(wp),intent(inout),dimension(nPoints) :: & ! |
|---|
| 259 | Cloud_Fraction_Total_Mean, & ! |
|---|
| 260 | Cloud_Fraction_Water_Mean, & ! |
|---|
| 261 | Cloud_Fraction_Ice_Mean, & ! |
|---|
| 262 | Cloud_Fraction_High_Mean, & ! |
|---|
| 263 | Cloud_Fraction_Mid_Mean, & ! |
|---|
| 264 | Cloud_Fraction_Low_Mean, & ! |
|---|
| 265 | Optical_Thickness_Total_Mean, & ! |
|---|
| 266 | Optical_Thickness_Water_Mean, & ! |
|---|
| 267 | Optical_Thickness_Ice_Mean, & ! |
|---|
| 268 | Optical_Thickness_Total_MeanLog10, & ! |
|---|
| 269 | Optical_Thickness_Water_MeanLog10, & ! |
|---|
| 270 | Optical_Thickness_Ice_MeanLog10, & ! |
|---|
| 271 | Cloud_Particle_Size_Water_Mean, & ! |
|---|
| 272 | Cloud_Particle_Size_Ice_Mean, & ! |
|---|
| 273 | Cloud_Top_Pressure_Total_Mean, & ! |
|---|
| 274 | Liquid_Water_Path_Mean, & ! |
|---|
| 275 | Ice_Water_Path_Mean ! |
|---|
| 276 | real(wp),intent(inout),dimension(nPoints,numMODISTauBins,numMODISPresBins) :: & |
|---|
| 277 | Optical_Thickness_vs_Cloud_Top_Pressure |
|---|
| 278 | real(wp),intent(inout),dimension(nPoints,numMODISTauBins,numMODISReffIceBins) :: & |
|---|
| 279 | Optical_Thickness_vs_ReffIce |
|---|
| 280 | real(wp),intent(inout),dimension(nPoints,numMODISTauBins,numMODISReffLiqBins) :: & |
|---|
| 281 | Optical_Thickness_vs_ReffLiq |
|---|
| 282 | |
|---|
| 283 | ! LOCAL VARIABLES |
|---|
| 284 | real(wp), parameter :: & |
|---|
| 285 | LWP_conversion = 2._wp/3._wp * 1000._wp ! MKS units |
|---|
| 286 | integer :: j |
|---|
| 287 | logical, dimension(nPoints,nSubCols) :: & |
|---|
| 288 | cloudMask, & |
|---|
| 289 | waterCloudMask, & |
|---|
| 290 | iceCloudMask, & |
|---|
| 291 | validRetrievalMask |
|---|
| 292 | real(wp),dimension(nPoints,nSubCols) :: & |
|---|
| 293 | tauWRK,ctpWRK,reffIceWRK,reffLiqWRK |
|---|
| 294 | |
|---|
| 295 | ! ######################################################################################## |
|---|
| 296 | ! Include only those pixels with successful retrievals in the statistics |
|---|
| 297 | ! ######################################################################################## |
|---|
| 298 | validRetrievalMask(1:nPoints,1:nSubCols) = particle_size(1:nPoints,1:nSubCols) > 0. |
|---|
| 299 | cloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) /= phaseIsNone .and. & |
|---|
| 300 | validRetrievalMask(1:nPoints,1:nSubCols) |
|---|
| 301 | waterCloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) == phaseIsLiquid .and. & |
|---|
| 302 | validRetrievalMask(1:nPoints,1:nSubCols) |
|---|
| 303 | iceCloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) == phaseIsIce .and. & |
|---|
| 304 | validRetrievalMask(1:nPoints,1:nSubCols) |
|---|
| 305 | |
|---|
| 306 | ! ######################################################################################## |
|---|
| 307 | ! Use these as pixel counts at first |
|---|
| 308 | ! ######################################################################################## |
|---|
| 309 | Cloud_Fraction_Total_Mean(1:nPoints) = real(count(cloudMask, dim = 2)) |
|---|
| 310 | Cloud_Fraction_Water_Mean(1:nPoints) = real(count(waterCloudMask, dim = 2)) |
|---|
| 311 | Cloud_Fraction_Ice_Mean(1:nPoints) = real(count(iceCloudMask, dim = 2)) |
|---|
| 312 | Cloud_Fraction_High_Mean(1:nPoints) = real(count(cloudMask .and. cloud_top_pressure <= & |
|---|
| 313 | highCloudPressureLimit, dim = 2)) |
|---|
| 314 | Cloud_Fraction_Low_Mean(1:nPoints) = real(count(cloudMask .and. cloud_top_pressure > & |
|---|
| 315 | lowCloudPressureLimit, dim = 2)) |
|---|
| 316 | Cloud_Fraction_Mid_Mean(1:nPoints) = Cloud_Fraction_Total_Mean(1:nPoints) - Cloud_Fraction_High_Mean(1:nPoints)& |
|---|
| 317 | - Cloud_Fraction_Low_Mean(1:nPoints) |
|---|
| 318 | |
|---|
| 319 | ! ######################################################################################## |
|---|
| 320 | ! Compute column amounts. |
|---|
| 321 | ! ######################################################################################## |
|---|
| 322 | where(Cloud_Fraction_Total_Mean(1:nPoints) > 0) |
|---|
| 323 | Optical_Thickness_Total_Mean(1:nPoints) = sum(optical_thickness, mask = cloudMask, dim = 2) / & |
|---|
| 324 | Cloud_Fraction_Total_Mean(1:nPoints) |
|---|
| 325 | Optical_Thickness_Total_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = cloudMask, & |
|---|
| 326 | dim = 2) / Cloud_Fraction_Total_Mean(1:nPoints) |
|---|
| 327 | elsewhere |
|---|
| 328 | Optical_Thickness_Total_Mean = R_UNDEF |
|---|
| 329 | Optical_Thickness_Total_MeanLog10 = R_UNDEF |
|---|
| 330 | endwhere |
|---|
| 331 | where(Cloud_Fraction_Water_Mean(1:nPoints) > 0) |
|---|
| 332 | Optical_Thickness_Water_Mean(1:nPoints) = sum(optical_thickness, mask = waterCloudMask, dim = 2) / & |
|---|
| 333 | Cloud_Fraction_Water_Mean(1:nPoints) |
|---|
| 334 | Liquid_Water_Path_Mean(1:nPoints) = LWP_conversion*sum(particle_size*optical_thickness, & |
|---|
| 335 | mask=waterCloudMask,dim=2)/Cloud_Fraction_Water_Mean(1:nPoints) |
|---|
| 336 | Optical_Thickness_Water_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = waterCloudMask,& |
|---|
| 337 | dim = 2) / Cloud_Fraction_Water_Mean(1:nPoints) |
|---|
| 338 | Cloud_Particle_Size_Water_Mean(1:nPoints) = sum(particle_size, mask = waterCloudMask, dim = 2) / & |
|---|
| 339 | Cloud_Fraction_Water_Mean(1:nPoints) |
|---|
| 340 | elsewhere |
|---|
| 341 | Optical_Thickness_Water_Mean = R_UNDEF |
|---|
| 342 | Optical_Thickness_Water_MeanLog10 = R_UNDEF |
|---|
| 343 | Cloud_Particle_Size_Water_Mean = R_UNDEF |
|---|
| 344 | Liquid_Water_Path_Mean = R_UNDEF |
|---|
| 345 | endwhere |
|---|
| 346 | where(Cloud_Fraction_Ice_Mean(1:nPoints) > 0) |
|---|
| 347 | Optical_Thickness_Ice_Mean(1:nPoints) = sum(optical_thickness, mask = iceCloudMask, dim = 2) / & |
|---|
| 348 | Cloud_Fraction_Ice_Mean(1:nPoints) |
|---|
| 349 | Ice_Water_Path_Mean(1:nPoints) = LWP_conversion * ice_density*sum(particle_size*optical_thickness,& |
|---|
| 350 | mask=iceCloudMask,dim = 2) /Cloud_Fraction_Ice_Mean(1:nPoints) |
|---|
| 351 | Optical_Thickness_Ice_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = iceCloudMask,& |
|---|
| 352 | dim = 2) / Cloud_Fraction_Ice_Mean(1:nPoints) |
|---|
| 353 | Cloud_Particle_Size_Ice_Mean(1:nPoints) = sum(particle_size, mask = iceCloudMask, dim = 2) / & |
|---|
| 354 | Cloud_Fraction_Ice_Mean(1:nPoints) |
|---|
| 355 | elsewhere |
|---|
| 356 | Optical_Thickness_Ice_Mean = R_UNDEF |
|---|
| 357 | Optical_Thickness_Ice_MeanLog10 = R_UNDEF |
|---|
| 358 | Cloud_Particle_Size_Ice_Mean = R_UNDEF |
|---|
| 359 | Ice_Water_Path_Mean = R_UNDEF |
|---|
| 360 | endwhere |
|---|
| 361 | Cloud_Top_Pressure_Total_Mean = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / & |
|---|
| 362 | max(1, count(cloudMask, dim = 2)) |
|---|
| 363 | |
|---|
| 364 | ! ######################################################################################## |
|---|
| 365 | ! Normalize pixel counts to fraction. |
|---|
| 366 | ! ######################################################################################## |
|---|
| 367 | Cloud_Fraction_High_Mean(1:nPoints) = Cloud_Fraction_High_Mean(1:nPoints) /nSubcols |
|---|
| 368 | Cloud_Fraction_Mid_Mean(1:nPoints) = Cloud_Fraction_Mid_Mean(1:nPoints) /nSubcols |
|---|
| 369 | Cloud_Fraction_Low_Mean(1:nPoints) = Cloud_Fraction_Low_Mean(1:nPoints) /nSubcols |
|---|
| 370 | Cloud_Fraction_Total_Mean(1:nPoints) = Cloud_Fraction_Total_Mean(1:nPoints) /nSubcols |
|---|
| 371 | Cloud_Fraction_Ice_Mean(1:nPoints) = Cloud_Fraction_Ice_Mean(1:nPoints) /nSubcols |
|---|
| 372 | Cloud_Fraction_Water_Mean(1:nPoints) = Cloud_Fraction_Water_Mean(1:nPoints) /nSubcols |
|---|
| 373 | |
|---|
| 374 | ! ######################################################################################## |
|---|
| 375 | ! Joint histograms |
|---|
| 376 | ! ######################################################################################## |
|---|
| 377 | ! Loop over all points |
|---|
| 378 | tauWRK(1:nPoints,1:nSubCols) = optical_thickness(1:nPoints,1:nSubCols) |
|---|
| 379 | ctpWRK(1:nPoints,1:nSubCols) = cloud_top_pressure(1:nPoints,1:nSubCols) |
|---|
| 380 | reffIceWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,iceCloudMask) |
|---|
| 381 | reffLiqWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,waterCloudMask) |
|---|
| 382 | do j=1,nPoints |
|---|
| 383 | |
|---|
| 384 | ! Fill clear and optically thin subcolumns with fill |
|---|
| 385 | where(.not. cloudMask(j,1:nSubCols)) |
|---|
| 386 | tauWRK(j,1:nSubCols) = -999._wp |
|---|
| 387 | ctpWRK(j,1:nSubCols) = -999._wp |
|---|
| 388 | endwhere |
|---|
| 389 | ! Joint histogram of tau/CTP |
|---|
| 390 | call hist2D(tauWRK(j,1:nSubCols),ctpWRK(j,1:nSubCols),nSubCols,& |
|---|
| 391 | modis_histTau,numMODISTauBins,& |
|---|
| 392 | modis_histPres,numMODISPresBins,& |
|---|
| 393 | Optical_Thickness_vs_Cloud_Top_Pressure(j,1:numMODISTauBins,1:numMODISPresBins)) |
|---|
| 394 | ! Joint histogram of tau/ReffICE |
|---|
| 395 | call hist2D(tauWRK(j,1:nSubCols),reffIceWrk(j,1:nSubCols),nSubCols, & |
|---|
| 396 | modis_histTau,numMODISTauBins,modis_histReffIce, & |
|---|
| 397 | numMODISReffIceBins, Optical_Thickness_vs_ReffIce(j,1:numMODISTauBins,1:numMODISReffIceBins)) |
|---|
| 398 | ! Joint histogram of tau/ReffLIQ |
|---|
| 399 | call hist2D(tauWRK(j,1:nSubCols),reffLiqWrk(j,1:nSubCols),nSubCols, & |
|---|
| 400 | modis_histTau,numMODISTauBins,modis_histReffLiq, & |
|---|
| 401 | numMODISReffLiqBins, Optical_Thickness_vs_ReffLiq(j,1:numMODISTauBins,1:numMODISReffLiqBins)) |
|---|
| 402 | |
|---|
| 403 | enddo |
|---|
| 404 | Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numMODISTauBins,1:numMODISPresBins) = & |
|---|
| 405 | Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numMODISTauBins,1:numMODISPresBins)/nSubCols |
|---|
| 406 | Optical_Thickness_vs_ReffIce(1:nPoints,1:numMODISTauBins,1:numMODISReffIceBins) = & |
|---|
| 407 | Optical_Thickness_vs_ReffIce(1:nPoints,1:numMODISTauBins,1:numMODISReffIceBins)/nSubCols |
|---|
| 408 | Optical_Thickness_vs_ReffLiq(1:nPoints,1:numMODISTauBins,1:numMODISReffLiqBins) = & |
|---|
| 409 | Optical_Thickness_vs_ReffLiq(1:nPoints,1:numMODISTauBins,1:numMODISReffLiqBins)/nSubCols |
|---|
| 410 | |
|---|
| 411 | |
|---|
| 412 | ! Unit conversion |
|---|
| 413 | where(Optical_Thickness_vs_Cloud_Top_Pressure /= R_UNDEF) & |
|---|
| 414 | Optical_Thickness_vs_Cloud_Top_Pressure = Optical_Thickness_vs_Cloud_Top_Pressure*100._wp |
|---|
| 415 | where(Optical_Thickness_vs_ReffIce /= R_UNDEF) Optical_Thickness_vs_ReffIce = Optical_Thickness_vs_ReffIce*100._wp |
|---|
| 416 | where(Optical_Thickness_vs_ReffLiq /= R_UNDEF) Optical_Thickness_vs_ReffLiq = Optical_Thickness_vs_ReffLiq*100._wp |
|---|
| 417 | where(Cloud_Fraction_Total_Mean /= R_UNDEF) Cloud_Fraction_Total_Mean = Cloud_Fraction_Total_Mean*100._wp |
|---|
| 418 | where(Cloud_Fraction_Water_Mean /= R_UNDEF) Cloud_Fraction_Water_Mean = Cloud_Fraction_Water_Mean*100._wp |
|---|
| 419 | where(Cloud_Fraction_Ice_Mean /= R_UNDEF) Cloud_Fraction_Ice_Mean = Cloud_Fraction_Ice_Mean*100._wp |
|---|
| 420 | where(Cloud_Fraction_High_Mean /= R_UNDEF) Cloud_Fraction_High_Mean = Cloud_Fraction_High_Mean*100._wp |
|---|
| 421 | where(Cloud_Fraction_Mid_Mean /= R_UNDEF) Cloud_Fraction_Mid_Mean = Cloud_Fraction_Mid_Mean*100._wp |
|---|
| 422 | where(Cloud_Fraction_Low_Mean /= R_UNDEF) Cloud_Fraction_Low_Mean = Cloud_Fraction_Low_Mean*100._wp |
|---|
| 423 | |
|---|
| 424 | end subroutine modis_column |
|---|
| 425 | |
|---|
| 426 | ! ######################################################################################## |
|---|
| 427 | function cloud_top_pressure(nLevels,tauIncrement, pressure, tauLimit) |
|---|
| 428 | ! INPUTS |
|---|
| 429 | integer, intent(in) :: nLevels |
|---|
| 430 | real(wp),intent(in),dimension(nLevels) :: tauIncrement, pressure |
|---|
| 431 | real(wp),intent(in) :: tauLimit |
|---|
| 432 | ! OUTPUTS |
|---|
| 433 | real(wp) :: cloud_top_pressure |
|---|
| 434 | ! LOCAL VARIABLES |
|---|
| 435 | real(wp) :: deltaX, totalTau, totalProduct |
|---|
| 436 | integer :: i |
|---|
| 437 | |
|---|
| 438 | ! Find the extinction-weighted pressure. Assume that pressure varies linearly between |
|---|
| 439 | ! layers and use the trapezoidal rule. |
|---|
| 440 | totalTau = 0._wp; totalProduct = 0._wp |
|---|
| 441 | do i = 2, size(tauIncrement) |
|---|
| 442 | if(totalTau + tauIncrement(i) > tauLimit) then |
|---|
| 443 | deltaX = tauLimit - totalTau |
|---|
| 444 | totalTau = totalTau + deltaX |
|---|
| 445 | ! |
|---|
| 446 | ! Result for trapezoidal rule when you take less than a full step |
|---|
| 447 | ! tauIncrement is a layer-integrated value |
|---|
| 448 | ! |
|---|
| 449 | totalProduct = totalProduct & |
|---|
| 450 | + pressure(i-1) * deltaX & |
|---|
| 451 | + (pressure(i) - pressure(i-1)) * deltaX**2/(2._wp * tauIncrement(i)) |
|---|
| 452 | else |
|---|
| 453 | totalTau = totalTau + tauIncrement(i) |
|---|
| 454 | totalProduct = totalProduct + tauIncrement(i) * (pressure(i) + pressure(i-1)) / 2._wp |
|---|
| 455 | end if |
|---|
| 456 | if(totalTau >= tauLimit) exit |
|---|
| 457 | end do |
|---|
| 458 | |
|---|
| 459 | if (totalTau > 0._wp) then |
|---|
| 460 | cloud_top_pressure = totalProduct/totalTau |
|---|
| 461 | else |
|---|
| 462 | cloud_top_pressure = 0._wp |
|---|
| 463 | endif |
|---|
| 464 | |
|---|
| 465 | end function cloud_top_pressure |
|---|
| 466 | |
|---|
| 467 | ! ######################################################################################## |
|---|
| 468 | function weight_by_extinction(nLevels,tauIncrement, f, tauLimit) |
|---|
| 469 | ! INPUTS |
|---|
| 470 | integer, intent(in) :: nLevels |
|---|
| 471 | real(wp),intent(in),dimension(nLevels) :: tauIncrement, f |
|---|
| 472 | real(wp),intent(in) :: tauLimit |
|---|
| 473 | ! OUTPUTS |
|---|
| 474 | real(wp) :: weight_by_extinction |
|---|
| 475 | ! LOCAL VARIABLES |
|---|
| 476 | real(wp) :: deltaX, totalTau, totalProduct |
|---|
| 477 | integer :: i |
|---|
| 478 | |
|---|
| 479 | ! Find the extinction-weighted value of f(tau), assuming constant f within each layer |
|---|
| 480 | totalTau = 0._wp; totalProduct = 0._wp |
|---|
| 481 | do i = 1, size(tauIncrement) |
|---|
| 482 | if(totalTau + tauIncrement(i) > tauLimit) then |
|---|
| 483 | deltaX = tauLimit - totalTau |
|---|
| 484 | totalTau = totalTau + deltaX |
|---|
| 485 | totalProduct = totalProduct + deltaX * f(i) |
|---|
| 486 | else |
|---|
| 487 | totalTau = totalTau + tauIncrement(i) |
|---|
| 488 | totalProduct = totalProduct + tauIncrement(i) * f(i) |
|---|
| 489 | end if |
|---|
| 490 | if(totalTau >= tauLimit) exit |
|---|
| 491 | end do |
|---|
| 492 | |
|---|
| 493 | if (totalTau > 0._wp) then |
|---|
| 494 | weight_by_extinction = totalProduct/totalTau |
|---|
| 495 | else |
|---|
| 496 | weight_by_extinction = 0._wp |
|---|
| 497 | endif |
|---|
| 498 | |
|---|
| 499 | end function weight_by_extinction |
|---|
| 500 | |
|---|
| 501 | ! ######################################################################################## |
|---|
| 502 | pure function interpolate_to_min(x, y, yobs) |
|---|
| 503 | ! INPUTS |
|---|
| 504 | real(wp),intent(in),dimension(num_trial_res) :: x, y |
|---|
| 505 | real(wp),intent(in) :: yobs |
|---|
| 506 | ! OUTPUTS |
|---|
| 507 | real(wp) :: interpolate_to_min |
|---|
| 508 | ! LOCAL VARIABLES |
|---|
| 509 | real(wp), dimension(num_trial_res) :: diff |
|---|
| 510 | integer :: nPoints, minDiffLoc, lowerBound, upperBound |
|---|
| 511 | |
|---|
| 512 | ! Given a set of values of y as y(x), find the value of x that minimizes abs(y - yobs) |
|---|
| 513 | ! y must be monotonic in x |
|---|
| 514 | |
|---|
| 515 | nPoints = size(y) |
|---|
| 516 | diff(1:num_trial_res) = y(1:num_trial_res) - yobs |
|---|
| 517 | minDiffLoc = minloc(abs(diff), dim = 1) |
|---|
| 518 | |
|---|
| 519 | if(minDiffLoc == 1) then |
|---|
| 520 | lowerBound = minDiffLoc |
|---|
| 521 | upperBound = minDiffLoc + 1 |
|---|
| 522 | else if(minDiffLoc == nPoints) then |
|---|
| 523 | lowerBound = minDiffLoc - 1 |
|---|
| 524 | upperBound = minDiffLoc |
|---|
| 525 | else |
|---|
| 526 | if(diff(minDiffLoc-1) * diff(minDiffLoc) < 0) then |
|---|
| 527 | lowerBound = minDiffLoc-1 |
|---|
| 528 | upperBound = minDiffLoc |
|---|
| 529 | else |
|---|
| 530 | lowerBound = minDiffLoc |
|---|
| 531 | upperBound = minDiffLoc + 1 |
|---|
| 532 | end if |
|---|
| 533 | end if |
|---|
| 534 | |
|---|
| 535 | if(diff(lowerBound) * diff(upperBound) < 0) then |
|---|
| 536 | ! |
|---|
| 537 | ! Interpolate the root position linearly if we bracket the root |
|---|
| 538 | ! |
|---|
| 539 | interpolate_to_min = x(upperBound) - & |
|---|
| 540 | diff(upperBound) * (x(upperBound) - x(lowerBound)) / (diff(upperBound) - diff(lowerBound)) |
|---|
| 541 | else |
|---|
| 542 | interpolate_to_min = re_fill |
|---|
| 543 | end if |
|---|
| 544 | |
|---|
| 545 | |
|---|
| 546 | end function interpolate_to_min |
|---|
| 547 | |
|---|
| 548 | ! ######################################################################################## |
|---|
| 549 | ! Optical properties |
|---|
| 550 | ! ######################################################################################## |
|---|
| 551 | elemental function get_g_nir_old (phase, re) |
|---|
| 552 | ! Polynomial fit for asummetry parameter g in MODIS band 7 (near IR) as a function |
|---|
| 553 | ! of size for ice and water |
|---|
| 554 | ! Fits from Steve Platnick |
|---|
| 555 | |
|---|
| 556 | ! INPUTS |
|---|
| 557 | integer, intent(in) :: phase |
|---|
| 558 | real(wp),intent(in) :: re |
|---|
| 559 | ! OUTPUTS |
|---|
| 560 | real(wp) :: get_g_nir_old |
|---|
| 561 | ! LOCAL VARIABLES(parameters) |
|---|
| 562 | real(wp), dimension(3), parameter :: & |
|---|
| 563 | ice_coefficients = (/ 0.7432, 4.5563e-3, -2.8697e-5 /), & |
|---|
| 564 | small_water_coefficients = (/ 0.8027, -1.0496e-2, 1.7071e-3 /), & |
|---|
| 565 | big_water_coefficients = (/ 0.7931, 5.3087e-3, -7.4995e-5 /) |
|---|
| 566 | |
|---|
| 567 | ! approx. fits from MODIS Collection 5 LUT scattering calculations |
|---|
| 568 | if(phase == phaseIsLiquid) then |
|---|
| 569 | if(re < 8.) then |
|---|
| 570 | get_g_nir_old = fit_to_quadratic(re, small_water_coefficients) |
|---|
| 571 | if(re < re_water_min) get_g_nir_old = fit_to_quadratic(re_water_min, small_water_coefficients) |
|---|
| 572 | else |
|---|
| 573 | get_g_nir_old = fit_to_quadratic(re, big_water_coefficients) |
|---|
| 574 | if(re > re_water_max) get_g_nir_old = fit_to_quadratic(re_water_max, big_water_coefficients) |
|---|
| 575 | end if |
|---|
| 576 | else |
|---|
| 577 | get_g_nir_old = fit_to_quadratic(re, ice_coefficients) |
|---|
| 578 | if(re < re_ice_min) get_g_nir_old = fit_to_quadratic(re_ice_min, ice_coefficients) |
|---|
| 579 | if(re > re_ice_max) get_g_nir_old = fit_to_quadratic(re_ice_max, ice_coefficients) |
|---|
| 580 | end if |
|---|
| 581 | |
|---|
| 582 | end function get_g_nir_old |
|---|
| 583 | |
|---|
| 584 | ! ######################################################################################## |
|---|
| 585 | elemental function get_ssa_nir_old (phase, re) |
|---|
| 586 | ! Polynomial fit for single scattering albedo in MODIS band 7 (near IR) as a function |
|---|
| 587 | ! of size for ice and water |
|---|
| 588 | ! Fits from Steve Platnick |
|---|
| 589 | |
|---|
| 590 | ! INPUTS |
|---|
| 591 | integer, intent(in) :: phase |
|---|
| 592 | real(wp),intent(in) :: re |
|---|
| 593 | ! OUTPUTS |
|---|
| 594 | real(wp) :: get_ssa_nir_old |
|---|
| 595 | ! LOCAL VARIABLES (parameters) |
|---|
| 596 | real(wp), dimension(4), parameter :: ice_coefficients = (/ 0.9994, -4.5199e-3, 3.9370e-5, -1.5235e-7 /) |
|---|
| 597 | real(wp), dimension(3), parameter :: water_coefficients = (/ 1.0008, -2.5626e-3, 1.6024e-5 /) |
|---|
| 598 | |
|---|
| 599 | ! approx. fits from MODIS Collection 5 LUT scattering calculations |
|---|
| 600 | if(phase == phaseIsLiquid) then |
|---|
| 601 | get_ssa_nir_old = fit_to_quadratic(re, water_coefficients) |
|---|
| 602 | if(re < re_water_min) get_ssa_nir_old = fit_to_quadratic(re_water_min, water_coefficients) |
|---|
| 603 | if(re > re_water_max) get_ssa_nir_old = fit_to_quadratic(re_water_max, water_coefficients) |
|---|
| 604 | else |
|---|
| 605 | get_ssa_nir_old = fit_to_cubic(re, ice_coefficients) |
|---|
| 606 | if(re < re_ice_min) get_ssa_nir_old = fit_to_cubic(re_ice_min, ice_coefficients) |
|---|
| 607 | if(re > re_ice_max) get_ssa_nir_old = fit_to_cubic(re_ice_max, ice_coefficients) |
|---|
| 608 | end if |
|---|
| 609 | |
|---|
| 610 | end function get_ssa_nir_old |
|---|
| 611 | |
|---|
| 612 | elemental function get_g_nir (phase, re) |
|---|
| 613 | ! |
|---|
| 614 | ! Polynomial fit for asummetry parameter g in MODIS band 7 (near IR) as a function |
|---|
| 615 | ! of size for ice and water |
|---|
| 616 | ! Fits from Steve Platnick |
|---|
| 617 | ! |
|---|
| 618 | |
|---|
| 619 | integer, intent(in) :: phase |
|---|
| 620 | real(wp), intent(in) :: re |
|---|
| 621 | real(wp) :: get_g_nir |
|---|
| 622 | |
|---|
| 623 | real(wp), dimension(3), parameter :: ice_coefficients = (/ 0.7490, 6.5153e-3, -5.4136e-5 /), & |
|---|
| 624 | small_water_coefficients = (/ 1.0364, -8.8800e-2, 7.0000e-3 /) |
|---|
| 625 | real(wp), dimension(4), parameter :: big_water_coefficients = (/ 0.6035, 2.8993e-2, -1.1051e-3, 1.5134e-5 /) |
|---|
| 626 | |
|---|
| 627 | ! approx. fits from MODIS Collection 6 LUT scattering calculations for 3.7 µm channel size retrievals |
|---|
| 628 | if(phase == phaseIsLiquid) then |
|---|
| 629 | if(re < 7.) then |
|---|
| 630 | get_g_nir = fit_to_quadratic(re, small_water_coefficients) |
|---|
| 631 | if(re < re_water_min) get_g_nir = fit_to_quadratic(re_water_min, small_water_coefficients) |
|---|
| 632 | else |
|---|
| 633 | get_g_nir = fit_to_cubic(re, big_water_coefficients) |
|---|
| 634 | if(re > re_water_max) get_g_nir = fit_to_cubic(re_water_max, big_water_coefficients) |
|---|
| 635 | end if |
|---|
| 636 | else |
|---|
| 637 | get_g_nir = fit_to_quadratic(re, ice_coefficients) |
|---|
| 638 | if(re < re_ice_min) get_g_nir = fit_to_quadratic(re_ice_min, ice_coefficients) |
|---|
| 639 | if(re > re_ice_max) get_g_nir = fit_to_quadratic(re_ice_max, ice_coefficients) |
|---|
| 640 | end if |
|---|
| 641 | |
|---|
| 642 | end function get_g_nir |
|---|
| 643 | |
|---|
| 644 | ! -------------------------------------------- |
|---|
| 645 | elemental function get_ssa_nir (phase, re) |
|---|
| 646 | integer, intent(in) :: phase |
|---|
| 647 | real(wp), intent(in) :: re |
|---|
| 648 | real(wp) :: get_ssa_nir |
|---|
| 649 | ! |
|---|
| 650 | ! Polynomial fit for single scattering albedo in MODIS band 7 (near IR) as a function |
|---|
| 651 | ! of size for ice and water |
|---|
| 652 | ! Fits from Steve Platnick |
|---|
| 653 | ! |
|---|
| 654 | real(wp), dimension(4), parameter :: ice_coefficients = (/ 0.9625, -1.8069e-2, 3.3281e-4,-2.2865e-6/) |
|---|
| 655 | real(wp), dimension(3), parameter :: water_coefficients = (/ 1.0044, -1.1397e-2, 1.3300e-4 /) |
|---|
| 656 | |
|---|
| 657 | ! approx. fits from MODIS Collection 6 LUT scattering calculations |
|---|
| 658 | if(phase == phaseIsLiquid) then |
|---|
| 659 | get_ssa_nir = fit_to_quadratic(re, water_coefficients) |
|---|
| 660 | if(re < re_water_min) get_ssa_nir = fit_to_quadratic(re_water_min, water_coefficients) |
|---|
| 661 | if(re > re_water_max) get_ssa_nir = fit_to_quadratic(re_water_max, water_coefficients) |
|---|
| 662 | else |
|---|
| 663 | get_ssa_nir = fit_to_cubic(re, ice_coefficients) |
|---|
| 664 | if(re < re_ice_min) get_ssa_nir = fit_to_cubic(re_ice_min, ice_coefficients) |
|---|
| 665 | if(re > re_ice_max) get_ssa_nir = fit_to_cubic(re_ice_max, ice_coefficients) |
|---|
| 666 | end if |
|---|
| 667 | |
|---|
| 668 | end function get_ssa_nir |
|---|
| 669 | |
|---|
| 670 | |
|---|
| 671 | |
|---|
| 672 | ! ######################################################################################## |
|---|
| 673 | pure function fit_to_cubic(x, coefficients) |
|---|
| 674 | ! INPUTS |
|---|
| 675 | real(wp), intent(in) :: x |
|---|
| 676 | real(wp), dimension(4), intent(in) :: coefficients |
|---|
| 677 | ! OUTPUTS |
|---|
| 678 | real(wp) :: fit_to_cubic |
|---|
| 679 | |
|---|
| 680 | fit_to_cubic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3) + x * coefficients(4))) |
|---|
| 681 | end function fit_to_cubic |
|---|
| 682 | |
|---|
| 683 | ! ######################################################################################## |
|---|
| 684 | pure function fit_to_quadratic(x, coefficients) |
|---|
| 685 | ! INPUTS |
|---|
| 686 | real(wp), intent(in) :: x |
|---|
| 687 | real(wp), dimension(3), intent(in) :: coefficients |
|---|
| 688 | ! OUTPUTS |
|---|
| 689 | real(wp) :: fit_to_quadratic |
|---|
| 690 | |
|---|
| 691 | fit_to_quadratic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3))) |
|---|
| 692 | end function fit_to_quadratic |
|---|
| 693 | |
|---|
| 694 | ! ######################################################################################## |
|---|
| 695 | ! Radiative transfer |
|---|
| 696 | ! ######################################################################################## |
|---|
| 697 | pure function compute_toa_reflectace(nLevels,tau, g, w0) |
|---|
| 698 | ! This wrapper reports reflectance only and strips out non-cloudy elements from the |
|---|
| 699 | ! calculation |
|---|
| 700 | |
|---|
| 701 | ! INPUTS |
|---|
| 702 | integer,intent(in) :: nLevels |
|---|
| 703 | real(wp),intent(in),dimension(nLevels) :: tau, g, w0 |
|---|
| 704 | ! OUTPUTS |
|---|
| 705 | real(wp) :: compute_toa_reflectace |
|---|
| 706 | ! LOCAL VARIABLES |
|---|
| 707 | logical, dimension(nLevels) :: cloudMask |
|---|
| 708 | integer, dimension(count(tau(1:nLevels) > 0)) :: cloudIndicies |
|---|
| 709 | real(wp),dimension(count(tau(1:nLevels) > 0)) :: Refl,Trans |
|---|
| 710 | real(wp) :: Refl_tot, Trans_tot |
|---|
| 711 | integer :: i |
|---|
| 712 | |
|---|
| 713 | cloudMask(1:nLevels) = tau(1:nLevels) > 0. |
|---|
| 714 | cloudIndicies = pack((/ (i, i = 1, nLevels) /), mask = cloudMask) |
|---|
| 715 | do i = 1, size(cloudIndicies) |
|---|
| 716 | call two_stream(tau(cloudIndicies(i)), g(cloudIndicies(i)), w0(cloudIndicies(i)), Refl(i), Trans(i)) |
|---|
| 717 | end do |
|---|
| 718 | |
|---|
| 719 | call adding_doubling(count(tau(1:nLevels) > 0),Refl(:), Trans(:), Refl_tot, Trans_tot) |
|---|
| 720 | |
|---|
| 721 | compute_toa_reflectace = Refl_tot |
|---|
| 722 | |
|---|
| 723 | end function compute_toa_reflectace |
|---|
| 724 | |
|---|
| 725 | ! ######################################################################################## |
|---|
| 726 | pure subroutine two_stream(tauint, gint, w0int, ref, tra) |
|---|
| 727 | ! Compute reflectance in a single layer using the two stream approximation |
|---|
| 728 | ! The code itself is from Lazaros Oreopoulos via Steve Platnick |
|---|
| 729 | ! INPUTS |
|---|
| 730 | real(wp), intent(in) :: tauint, gint, w0int |
|---|
| 731 | ! OUTPUTS |
|---|
| 732 | real(wp), intent(out) :: ref, tra |
|---|
| 733 | ! LOCAL VARIABLES |
|---|
| 734 | ! for delta Eddington code |
|---|
| 735 | ! xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1) |
|---|
| 736 | integer, parameter :: beam = 2 |
|---|
| 737 | real(wp),parameter :: xmu = 0.866, minConservativeW0 = 0.9999999 |
|---|
| 738 | real(wp) :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, & |
|---|
| 739 | rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den, th |
|---|
| 740 | |
|---|
| 741 | ! Compute reflectance and transmittance in a single layer using the two stream approximation |
|---|
| 742 | ! The code itself is from Lazaros Oreopoulos via Steve Platnick |
|---|
| 743 | f = gint**2 |
|---|
| 744 | tau = (1._wp - w0int * f) * tauint |
|---|
| 745 | w0 = (1._wp - f) * w0int / (1._wp - w0int * f) |
|---|
| 746 | g = (gint - f) / (1._wp - f) |
|---|
| 747 | |
|---|
| 748 | ! delta-Eddington (Joseph et al. 1976) |
|---|
| 749 | gamma1 = (7._wp - w0* (4._wp + 3._wp * g)) / 4._wp |
|---|
| 750 | gamma2 = -(1._wp - w0* (4._wp - 3._wp * g)) / 4._wp |
|---|
| 751 | gamma3 = (2._wp - 3._wp*g*xmu) / 4._wp |
|---|
| 752 | gamma4 = 1._wp - gamma3 |
|---|
| 753 | |
|---|
| 754 | if (w0int > minConservativeW0) then |
|---|
| 755 | ! Conservative scattering |
|---|
| 756 | if (beam == 1) then |
|---|
| 757 | rh = (gamma1*tau+(gamma3-gamma1*xmu)*(1-exp(-tau/xmu))) |
|---|
| 758 | |
|---|
| 759 | ref = rh / (1._wp + gamma1 * tau) |
|---|
| 760 | tra = 1._wp - ref |
|---|
| 761 | else if(beam == 2) then |
|---|
| 762 | ref = gamma1*tau/(1._wp + gamma1*tau) |
|---|
| 763 | tra = 1._wp - ref |
|---|
| 764 | endif |
|---|
| 765 | else |
|---|
| 766 | ! Non-conservative scattering |
|---|
| 767 | a1 = gamma1 * gamma4 + gamma2 * gamma3 |
|---|
| 768 | a2 = gamma1 * gamma3 + gamma2 * gamma4 |
|---|
| 769 | |
|---|
| 770 | rk = sqrt(gamma1**2 - gamma2**2) |
|---|
| 771 | |
|---|
| 772 | r1 = (1._wp - rk * xmu) * (a2 + rk * gamma3) |
|---|
| 773 | r2 = (1._wp + rk * xmu) * (a2 - rk * gamma3) |
|---|
| 774 | r3 = 2._wp * rk *(gamma3 - a2 * xmu) |
|---|
| 775 | r4 = (1._wp - (rk * xmu)**2) * (rk + gamma1) |
|---|
| 776 | r5 = (1._wp - (rk * xmu)**2) * (rk - gamma1) |
|---|
| 777 | |
|---|
| 778 | t1 = (1._wp + rk * xmu) * (a1 + rk * gamma4) |
|---|
| 779 | t2 = (1._wp - rk * xmu) * (a1 - rk * gamma4) |
|---|
| 780 | t3 = 2._wp * rk * (gamma4 + a1 * xmu) |
|---|
| 781 | t4 = r4 |
|---|
| 782 | t5 = r5 |
|---|
| 783 | |
|---|
| 784 | beta = -r5 / r4 |
|---|
| 785 | |
|---|
| 786 | e1 = min(rk * tau, 500._wp) |
|---|
| 787 | e2 = min(tau / xmu, 500._wp) |
|---|
| 788 | |
|---|
| 789 | if (beam == 1) then |
|---|
| 790 | den = r4 * exp(e1) + r5 * exp(-e1) |
|---|
| 791 | ref = w0*(r1*exp(e1)-r2*exp(-e1)-r3*exp(-e2))/den |
|---|
| 792 | den = t4 * exp(e1) + t5 * exp(-e1) |
|---|
| 793 | th = exp(-e2) |
|---|
| 794 | tra = th-th*w0*(t1*exp(e1)-t2*exp(-e1)-t3*exp(e2))/den |
|---|
| 795 | elseif (beam == 2) then |
|---|
| 796 | ef1 = exp(-e1) |
|---|
| 797 | ef2 = exp(-2*e1) |
|---|
| 798 | ref = (gamma2*(1._wp-ef2))/((rk+gamma1)*(1._wp-beta*ef2)) |
|---|
| 799 | tra = (2._wp*rk*ef1)/((rk+gamma1)*(1._wp-beta*ef2)) |
|---|
| 800 | endif |
|---|
| 801 | end if |
|---|
| 802 | end subroutine two_stream |
|---|
| 803 | |
|---|
| 804 | ! ######################################################################################## |
|---|
| 805 | elemental function two_stream_reflectance(tauint, gint, w0int) |
|---|
| 806 | ! Compute reflectance in a single layer using the two stream approximation |
|---|
| 807 | ! The code itself is from Lazaros Oreopoulos via Steve Platnick |
|---|
| 808 | |
|---|
| 809 | ! INPUTS |
|---|
| 810 | real(wp), intent(in) :: tauint, gint, w0int |
|---|
| 811 | ! OUTPUTS |
|---|
| 812 | real(wp) :: two_stream_reflectance |
|---|
| 813 | ! LOCAL VARIABLES |
|---|
| 814 | ! for delta Eddington code |
|---|
| 815 | ! xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1) |
|---|
| 816 | integer, parameter :: beam = 2 |
|---|
| 817 | real(wp),parameter :: xmu = 0.866, minConservativeW0 = 0.9999999 |
|---|
| 818 | real(wp) :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, & |
|---|
| 819 | rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den |
|---|
| 820 | |
|---|
| 821 | f = gint**2 |
|---|
| 822 | tau = (1._wp - w0int * f) * tauint |
|---|
| 823 | w0 = (1._wp - f) * w0int / (1._wp - w0int * f) |
|---|
| 824 | g = (gint - f) / (1._wp - f) |
|---|
| 825 | |
|---|
| 826 | ! delta-Eddington (Joseph et al. 1976) |
|---|
| 827 | gamma1 = (7._wp - w0* (4._wp + 3._wp * g)) / 4._wp |
|---|
| 828 | gamma2 = -(1._wp - w0* (4._wp - 3._wp * g)) / 4._wp |
|---|
| 829 | gamma3 = (2._wp - 3._wp*g*xmu) / 4._wp |
|---|
| 830 | gamma4 = 1._wp - gamma3 |
|---|
| 831 | |
|---|
| 832 | if (w0int > minConservativeW0) then |
|---|
| 833 | ! Conservative scattering |
|---|
| 834 | if (beam == 1) then |
|---|
| 835 | rh = (gamma1*tau+(gamma3-gamma1*xmu)*(1-exp(-tau/xmu))) |
|---|
| 836 | two_stream_reflectance = rh / (1._wp + gamma1 * tau) |
|---|
| 837 | elseif (beam == 2) then |
|---|
| 838 | two_stream_reflectance = gamma1*tau/(1._wp + gamma1*tau) |
|---|
| 839 | endif |
|---|
| 840 | |
|---|
| 841 | else ! |
|---|
| 842 | |
|---|
| 843 | ! Non-conservative scattering |
|---|
| 844 | a1 = gamma1 * gamma4 + gamma2 * gamma3 |
|---|
| 845 | a2 = gamma1 * gamma3 + gamma2 * gamma4 |
|---|
| 846 | |
|---|
| 847 | rk = sqrt(gamma1**2 - gamma2**2) |
|---|
| 848 | |
|---|
| 849 | r1 = (1._wp - rk * xmu) * (a2 + rk * gamma3) |
|---|
| 850 | r2 = (1._wp + rk * xmu) * (a2 - rk * gamma3) |
|---|
| 851 | r3 = 2._wp * rk *(gamma3 - a2 * xmu) |
|---|
| 852 | r4 = (1._wp - (rk * xmu)**2) * (rk + gamma1) |
|---|
| 853 | r5 = (1._wp - (rk * xmu)**2) * (rk - gamma1) |
|---|
| 854 | |
|---|
| 855 | t1 = (1._wp + rk * xmu) * (a1 + rk * gamma4) |
|---|
| 856 | t2 = (1._wp - rk * xmu) * (a1 - rk * gamma4) |
|---|
| 857 | t3 = 2._wp * rk * (gamma4 + a1 * xmu) |
|---|
| 858 | t4 = r4 |
|---|
| 859 | t5 = r5 |
|---|
| 860 | |
|---|
| 861 | beta = -r5 / r4 |
|---|
| 862 | |
|---|
| 863 | e1 = min(rk * tau, 500._wp) |
|---|
| 864 | e2 = min(tau / xmu, 500._wp) |
|---|
| 865 | |
|---|
| 866 | if (beam == 1) then |
|---|
| 867 | den = r4 * exp(e1) + r5 * exp(-e1) |
|---|
| 868 | two_stream_reflectance = w0*(r1*exp(e1)-r2*exp(-e1)-r3*exp(-e2))/den |
|---|
| 869 | elseif (beam == 2) then |
|---|
| 870 | ef1 = exp(-e1) |
|---|
| 871 | ef2 = exp(-2*e1) |
|---|
| 872 | two_stream_reflectance = (gamma2*(1._wp-ef2))/((rk+gamma1)*(1._wp-beta*ef2)) |
|---|
| 873 | endif |
|---|
| 874 | |
|---|
| 875 | end if |
|---|
| 876 | end function two_stream_reflectance |
|---|
| 877 | |
|---|
| 878 | ! ######################################################################################## |
|---|
| 879 | pure subroutine adding_doubling (npts,Refl, Tran, Refl_tot, Tran_tot) |
|---|
| 880 | ! Use adding/doubling formulas to compute total reflectance and transmittance from |
|---|
| 881 | ! layer values |
|---|
| 882 | |
|---|
| 883 | ! INPUTS |
|---|
| 884 | integer,intent(in) :: npts |
|---|
| 885 | real(wp),intent(in),dimension(npts) :: Refl,Tran |
|---|
| 886 | ! OUTPUTS |
|---|
| 887 | real(wp),intent(out) :: Refl_tot, Tran_tot |
|---|
| 888 | ! LOCAL VARIABLES |
|---|
| 889 | integer :: i |
|---|
| 890 | real(wp), dimension(npts) :: Refl_cumulative, Tran_cumulative |
|---|
| 891 | |
|---|
| 892 | Refl_cumulative(1) = Refl(1) |
|---|
| 893 | Tran_cumulative(1) = Tran(1) |
|---|
| 894 | |
|---|
| 895 | do i=2, npts |
|---|
| 896 | ! place (add) previous combined layer(s) reflectance on top of layer i, w/black surface (or ignoring surface): |
|---|
| 897 | Refl_cumulative(i) = Refl_cumulative(i-1) + Refl(i)*(Tran_cumulative(i-1)**2)/(1._wp - Refl_cumulative(i-1) * Refl(i)) |
|---|
| 898 | Tran_cumulative(i) = (Tran_cumulative(i-1)*Tran(i)) / (1._wp - Refl_cumulative(i-1) * Refl(i)) |
|---|
| 899 | end do |
|---|
| 900 | |
|---|
| 901 | Refl_tot = Refl_cumulative(size(Refl)) |
|---|
| 902 | Tran_tot = Tran_cumulative(size(Refl)) |
|---|
| 903 | |
|---|
| 904 | end subroutine adding_doubling |
|---|
| 905 | |
|---|
| 906 | end module mod_modis_sim |
|---|
