[2432] | 1 | ! (c) 2009-2010, Regents of the Unversity of Colorado |
---|
| 2 | ! Author: Robert Pincus, Cooperative Institute for Research in the Environmental Sciences |
---|
| 3 | ! All rights reserved. |
---|
[2720] | 4 | ! $Revision: 88 $, $Date: 2013-11-13 07:08:38 -0700 (Wed, 13 Nov 2013) $ |
---|
[2432] | 5 | ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/MODIS_simulator/modis_simulator.F90 $ |
---|
| 6 | ! |
---|
| 7 | ! Redistribution and use in source and binary forms, with or without modification, are permitted |
---|
| 8 | ! provided that the following conditions are met: |
---|
| 9 | ! |
---|
| 10 | ! * Redistributions of source code must retain the above copyright notice, this list |
---|
| 11 | ! of conditions and the following disclaimer. |
---|
| 12 | ! * Redistributions in binary form must reproduce the above copyright notice, this list |
---|
| 13 | ! of conditions and the following disclaimer in the documentation and/or other materials |
---|
| 14 | ! provided with the distribution. |
---|
| 15 | ! * Neither the name of the Met Office nor the names of its contributors may be used |
---|
| 16 | ! to endorse or promote products derived from this software without specific prior written |
---|
| 17 | ! permission. |
---|
| 18 | ! |
---|
| 19 | ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR |
---|
| 20 | ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND |
---|
| 21 | ! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR |
---|
| 22 | ! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL |
---|
| 23 | ! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, |
---|
| 24 | ! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER |
---|
| 25 | ! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT |
---|
| 26 | ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
---|
| 27 | ! |
---|
| 28 | |
---|
| 29 | ! |
---|
| 30 | ! History: |
---|
| 31 | ! May 2009 - Robert Pincus - Initial version |
---|
| 32 | ! June 2009 - Steve Platnick and Robert Pincus - Simple radiative transfer for size 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 using AM2 (GFDL) |
---|
| 35 | ! January 2010 - Robert Pincus - Added high, middle, low cloud fractions |
---|
| 36 | ! |
---|
| 37 | |
---|
| 38 | ! |
---|
| 39 | ! Notes on using the MODIS simulator: |
---|
| 40 | ! *) You may provide either layer-by-layer values of optical thickness at 0.67 and 2.1 microns, or |
---|
| 41 | ! optical thickness at 0.67 microns and ice- and liquid-water contents (in consistent units of |
---|
| 42 | ! your choosing) |
---|
| 43 | ! *) Required input also includes the optical thickness and cloud top pressure |
---|
| 44 | ! derived from the ISCCP simulator run with parameter top_height = 1. |
---|
| 45 | ! *) Cloud particle sizes are specified as radii, measured in meters, though within the module we |
---|
| 46 | ! use units of microns. Where particle sizes are outside the bounds used in the MODIS retrieval |
---|
| 47 | ! libraries (parameters re_water_min, re_ice_min, etc.) the simulator returns missing values (re_fill) |
---|
| 48 | |
---|
| 49 | ! |
---|
| 50 | ! When error conditions are encountered this code calls the function complain_and_die, supplied at the |
---|
| 51 | ! bottom of this module. Users probably want to replace this with something more graceful. |
---|
| 52 | ! |
---|
| 53 | module mod_modis_sim |
---|
| 54 | USE MOD_COSP_TYPES, only: R_UNDEF |
---|
| 55 | implicit none |
---|
| 56 | ! ------------------------------ |
---|
| 57 | ! Algorithmic parameters |
---|
| 58 | ! |
---|
| 59 | |
---|
| 60 | real, parameter :: ice_density = 0.93 ! liquid density is 1. |
---|
| 61 | ! |
---|
| 62 | ! Retrieval parameters |
---|
| 63 | ! |
---|
| 64 | real, parameter :: min_OpticalThickness = 0.3, & ! Minimum detectable optical thickness |
---|
| 65 | CO2Slicing_PressureLimit = 700. * 100., & ! Cloud with higher pressures use thermal methods, units Pa |
---|
| 66 | CO2Slicing_TauLimit = 1., & ! How deep into the cloud does CO2 slicing see? |
---|
| 67 | phase_TauLimit = 1., & ! How deep into the cloud does the phase detection see? |
---|
| 68 | size_TauLimit = 2., & ! Depth of the re retreivals |
---|
| 69 | phaseDiscrimination_Threshold = 0.7 ! What fraction of total extincton needs to be |
---|
| 70 | ! in a single category to make phase discrim. work? |
---|
| 71 | real, parameter :: re_fill= -999. |
---|
| 72 | integer, parameter :: phaseIsNone = 0, phaseIsLiquid = 1, phaseIsIce = 2, phaseIsUndetermined = 3 |
---|
| 73 | |
---|
| 74 | logical, parameter :: useSimpleReScheme = .false. |
---|
| 75 | ! |
---|
| 76 | ! These are the limits of the libraries for the MODIS collection 5 algorithms |
---|
| 77 | ! They are also the limits used in the fits for g and w0 |
---|
| 78 | ! |
---|
| 79 | real, parameter :: re_water_min= 4., re_water_max= 30., re_ice_min= 5., re_ice_max= 90. |
---|
| 80 | integer, parameter :: num_trial_res = 15 ! increase to make the linear pseudo-retrieval of size more accurate |
---|
[2720] | 81 | ! DJS2015: Remove unused parameter |
---|
| 82 | ! logical, parameter :: use_two_re_iterations = .false. ! do two retrieval iterations? |
---|
| 83 | ! DJS2015 END |
---|
[2432] | 84 | ! |
---|
| 85 | ! Precompute near-IR optical params vs size for retrieval scheme |
---|
| 86 | ! |
---|
| 87 | integer, private :: i |
---|
| 88 | real, dimension(num_trial_res), parameter :: & |
---|
| 89 | trial_re_w = re_water_min + (re_water_max - re_water_min)/(num_trial_res-1) * (/ (i - 1, i = 1, num_trial_res) /), & |
---|
| 90 | trial_re_i = re_ice_min + (re_ice_max - re_ice_min) /(num_trial_res-1) * (/ (i - 1, i = 1, num_trial_res) /) |
---|
| 91 | |
---|
| 92 | ! Can't initialze these during compilation, but do in before looping columns in retrievals |
---|
| 93 | real, dimension(num_trial_res) :: g_w, g_i, w0_w, w0_i |
---|
| 94 | ! ------------------------------ |
---|
| 95 | ! Bin boundaries for the joint optical thickness/cloud top pressure histogram |
---|
| 96 | ! |
---|
| 97 | integer, parameter :: numTauHistogramBins = 6, numPressureHistogramBins = 7 |
---|
| 98 | |
---|
| 99 | real, private :: dummy_real |
---|
| 100 | real, dimension(numTauHistogramBins + 1), parameter :: & |
---|
| 101 | tauHistogramBoundaries = (/ min_OpticalThickness, 1.3, 3.6, 9.4, 23., 60., 10000. /) |
---|
| 102 | real, dimension(numPressureHistogramBins + 1), parameter :: & ! Units Pa |
---|
| 103 | pressureHistogramBoundaries = (/ 0., 18000., 31000., 44000., 56000., 68000., 80000., 1000000. /) |
---|
| 104 | real, parameter :: highCloudPressureLimit = 440. * 100., lowCloudPressureLimit = 680. * 100. |
---|
| 105 | ! |
---|
| 106 | ! For output - nominal bin centers and bin boundaries. On output pressure bins are highest to lowest. |
---|
| 107 | ! |
---|
| 108 | integer, private :: k, l |
---|
| 109 | real, parameter, dimension(2, numTauHistogramBins) :: & |
---|
| 110 | nominalTauHistogramBoundaries = & |
---|
| 111 | reshape(source = (/ tauHistogramBoundaries(1), & |
---|
| 112 | ((tauHistogramBoundaries(k), l = 1, 2), k = 2, numTauHistogramBins), & |
---|
| 113 | 100000. /), & |
---|
| 114 | shape = (/2, numTauHistogramBins /) ) |
---|
| 115 | real, parameter, dimension(numTauHistogramBins) :: & |
---|
| 116 | nominalTauHistogramCenters = (nominalTauHistogramBoundaries(1, :) + & |
---|
| 117 | nominalTauHistogramBoundaries(2, :) ) / 2. |
---|
| 118 | |
---|
| 119 | real, parameter, dimension(2, numPressureHistogramBins) :: & |
---|
| 120 | nominalPressureHistogramBoundaries = & |
---|
| 121 | reshape(source = (/ 100000., & |
---|
| 122 | ((pressureHistogramBoundaries(k), l = 1, 2), k = numPressureHistogramBins, 2, -1), & |
---|
| 123 | 0. /), & |
---|
| 124 | shape = (/2, numPressureHistogramBins /) ) |
---|
| 125 | real, parameter, dimension(numPressureHistogramBins) :: & |
---|
| 126 | nominalPressureHistogramCenters = (nominalPressureHistogramBoundaries(1, :) + & |
---|
| 127 | nominalPressureHistogramBoundaries(2, :) ) / 2. |
---|
[2720] | 128 | ! DJS2015 START: Add bin descriptions for joint-histograms of partice-sizes and optical depth. This is |
---|
| 129 | ! identical to what is done in COSPv.2.0.0 for histogram bin initialization. |
---|
| 130 | integer :: j |
---|
| 131 | integer,parameter :: & |
---|
| 132 | numMODISReffLiqBins = 6, & ! Number of bins for tau/ReffLiq joint-histogram |
---|
| 133 | numMODISReffIceBins = 6 ! Number of bins for tau/ReffICE joint-histogram |
---|
| 134 | real,parameter,dimension(numMODISReffLiqBins+1) :: & |
---|
| 135 | reffLIQ_binBounds = (/0., 8e-6, 1.0e-5, 1.3e-5, 1.5e-5, 2.0e-5, 3.0e-5/) |
---|
| 136 | real,parameter,dimension(numMODISReffIceBins+1) :: & |
---|
| 137 | reffICE_binBounds = (/0., 1.0e-5, 2.0e-5, 3.0e-5, 4.0e-5, 6.0e-5, 9.0e-5/) |
---|
| 138 | real,parameter,dimension(2,numMODISReffIceBins) :: & |
---|
| 139 | reffICE_binEdges = reshape(source=(/reffICE_binBounds(1),((reffICE_binBounds(k), & |
---|
| 140 | l=1,2),k=2,numMODISReffIceBins),reffICE_binBounds(numMODISReffIceBins+1)/), & |
---|
| 141 | shape = (/2,numMODISReffIceBins/)) |
---|
| 142 | real,parameter,dimension(2,numMODISReffLiqBins) :: & |
---|
| 143 | reffLIQ_binEdges = reshape(source=(/reffLIQ_binBounds(1),((reffLIQ_binBounds(k), & |
---|
| 144 | l=1,2),k=2,numMODISReffLiqBins),reffLIQ_binBounds(numMODISReffIceBins+1)/), & |
---|
| 145 | shape = (/2,numMODISReffLiqBins/)) |
---|
| 146 | real,parameter,dimension(numMODISReffIceBins) :: & |
---|
| 147 | reffICE_binCenters = (reffICE_binEdges(1,:)+reffICE_binEdges(2,:))/2. |
---|
| 148 | real,parameter,dimension(numMODISReffLiqBins) :: & |
---|
| 149 | reffLIQ_binCenters = (reffLIQ_binEdges(1,:)+reffLIQ_binEdges(2,:))/2. |
---|
| 150 | ! DJS2015 END |
---|
| 151 | |
---|
[2432] | 152 | ! ------------------------------ |
---|
| 153 | ! There are two ways to call the MODIS simulator: |
---|
| 154 | ! 1) Provide total optical thickness and liquid/ice water content and we'll partition tau in |
---|
| 155 | ! subroutine modis_L2_simulator_oneTau, or |
---|
| 156 | ! 2) Provide ice and liquid optical depths in each layer |
---|
| 157 | ! |
---|
| 158 | interface modis_L2_simulator |
---|
| 159 | module procedure modis_L2_simulator_oneTau, modis_L2_simulator_twoTaus |
---|
| 160 | end interface |
---|
| 161 | contains |
---|
| 162 | !------------------------------------------------------------------------------------------------ |
---|
| 163 | ! MODIS simulator using specified liquid and ice optical thickness in each layer |
---|
| 164 | ! |
---|
| 165 | ! Note: this simulator operates on all points; to match MODIS itself night-time |
---|
| 166 | ! points should be excluded |
---|
| 167 | ! |
---|
| 168 | ! Note: the simulator requires as input the optical thickness and cloud top pressure |
---|
| 169 | ! derived from the ISCCP simulator run with parameter top_height = 1. |
---|
| 170 | ! If cloud top pressure is higher than about 700 mb, MODIS can't use CO2 slicing |
---|
| 171 | ! and reverts to a thermal algorithm much like ISCCP's. Rather than replicate that |
---|
| 172 | ! alogrithm in this simulator we simply report the values from the ISCCP simulator. |
---|
| 173 | ! |
---|
| 174 | subroutine modis_L2_simulator_twoTaus( & |
---|
| 175 | temp, pressureLayers, pressureLevels, & |
---|
| 176 | liquid_opticalThickness, ice_opticalThickness, & |
---|
| 177 | waterSize, iceSize, & |
---|
| 178 | isccpTau, isccpCloudTopPressure, & |
---|
| 179 | retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize) |
---|
| 180 | |
---|
| 181 | ! Grid-mean quantities at layer centers, starting at the model top |
---|
| 182 | ! dimension nLayers |
---|
| 183 | real, dimension(:), intent(in ) :: temp, & ! Temperature, K |
---|
| 184 | pressureLayers, & ! Pressure, Pa |
---|
| 185 | pressureLevels ! Pressure at layer edges, Pa (dimension nLayers + 1) |
---|
| 186 | ! Sub-column quantities |
---|
| 187 | ! dimension nSubcols, nLayers |
---|
| 188 | real, dimension(:, :), intent(in ) :: liquid_opticalThickness, & ! Layer optical thickness @ 0.67 microns due to liquid |
---|
| 189 | ice_opticalThickness ! ditto, due to ice |
---|
| 190 | real, dimension(:, :), intent(in ) :: waterSize, & ! Cloud drop effective radius, microns |
---|
| 191 | iceSize ! Cloud ice effective radius, microns |
---|
| 192 | |
---|
| 193 | ! Cloud properties retrieved from ISCCP using top_height = 1 |
---|
| 194 | ! dimension nSubcols |
---|
| 195 | real, dimension(:), intent(in ) :: isccpTau, & ! Column-integrated optical thickness |
---|
| 196 | isccpCloudTopPressure ! ISCCP-retrieved cloud top pressure (Pa) |
---|
| 197 | |
---|
| 198 | ! Properties retrieved by MODIS |
---|
| 199 | ! dimension nSubcols |
---|
| 200 | integer, dimension(:), intent(out) :: retrievedPhase ! liquid/ice/other - integer, defined in module header |
---|
| 201 | real, dimension(:), intent(out) :: retrievedCloudTopPressure, & ! units of pressureLayers |
---|
| 202 | retrievedTau, & ! unitless |
---|
| 203 | retrievedSize ! microns |
---|
| 204 | ! --------------------------------------------------- |
---|
| 205 | ! Local variables |
---|
| 206 | logical, dimension(size(retrievedTau)) :: cloudMask |
---|
| 207 | real, dimension(size(waterSize, 1), size(waterSize, 2)) :: tauLiquidFraction, tauTotal |
---|
| 208 | real :: integratedLiquidFraction |
---|
| 209 | integer :: i, nSubcols, nLevels |
---|
| 210 | |
---|
| 211 | ! --------------------------------------------------- |
---|
| 212 | nSubcols = size(liquid_opticalThickness, 1) |
---|
| 213 | nLevels = size(liquid_opticalThickness, 2) |
---|
| 214 | |
---|
| 215 | ! |
---|
| 216 | ! Initial error checks |
---|
| 217 | ! |
---|
| 218 | if(any((/ size(ice_opticalThickness, 1), size(waterSize, 1), size(iceSize, 1), & |
---|
| 219 | size(isccpTau), size(isccpCloudTopPressure), & |
---|
| 220 | size(retrievedPhase), size(retrievedCloudTopPressure), & |
---|
| 221 | size(retrievedTau), size(retrievedSize) /) /= nSubcols )) & |
---|
| 222 | call complain_and_die("Differing number of subcolumns in one or more arrays") |
---|
| 223 | |
---|
| 224 | if(any((/ size(ice_opticalThickness, 2), size(waterSize, 2), size(iceSize, 2), & |
---|
| 225 | size(temp), size(pressureLayers), size(pressureLevels)-1 /) /= nLevels )) & |
---|
| 226 | call complain_and_die("Differing number of levels in one or more arrays") |
---|
| 227 | |
---|
| 228 | |
---|
| 229 | if(any( (/ any(temp <= 0.), any(pressureLayers <= 0.), & |
---|
| 230 | any(liquid_opticalThickness < 0.), & |
---|
| 231 | any(ice_opticalThickness < 0.), & |
---|
| 232 | any(waterSize < 0.), any(iceSize < 0.) /) )) & |
---|
| 233 | call complain_and_die("Input values out of bounds") |
---|
| 234 | |
---|
| 235 | ! --------------------------------------------------- |
---|
| 236 | ! |
---|
| 237 | ! Compute the total optical thickness and the proportion due to liquid in each cell |
---|
| 238 | ! |
---|
| 239 | where(liquid_opticalThickness(:, :) + ice_opticalThickness(:, :) > 0.) |
---|
| 240 | tauLiquidFraction(:, :) = liquid_opticalThickness(:, :)/(liquid_opticalThickness(:, :) + ice_opticalThickness(:, :)) |
---|
| 241 | elsewhere |
---|
| 242 | tauLiquidFraction(:, :) = 0. |
---|
| 243 | end where |
---|
| 244 | tauTotal(:, :) = liquid_opticalThickness(:, :) + ice_opticalThickness(:, :) |
---|
| 245 | |
---|
| 246 | ! |
---|
| 247 | ! Optical depth retrieval |
---|
| 248 | ! This is simply a sum over the optical thickness in each layer |
---|
| 249 | ! It should agree with the ISCCP values after min values have been excluded |
---|
| 250 | ! |
---|
| 251 | retrievedTau(:) = sum(tauTotal(:, :), dim = 2) |
---|
| 252 | |
---|
| 253 | ! |
---|
| 254 | ! Cloud detection - does optical thickness exceed detection threshold? |
---|
| 255 | ! |
---|
| 256 | cloudMask = retrievedTau(:) >= min_OpticalThickness |
---|
| 257 | |
---|
| 258 | ! |
---|
| 259 | ! Initialize initial estimates for size retrievals |
---|
| 260 | ! |
---|
| 261 | if(any(cloudMask) .and. .not. useSimpleReScheme) then |
---|
| 262 | g_w(:) = get_g_nir( phaseIsLiquid, trial_re_w(:)) |
---|
| 263 | w0_w(:) = get_ssa_nir(phaseIsLiquid, trial_re_w(:)) |
---|
| 264 | g_i(:) = get_g_nir( phaseIsIce, trial_re_i(:)) |
---|
| 265 | w0_i(:) = get_ssa_nir(phaseIsIce, trial_re_i(:)) |
---|
| 266 | end if |
---|
| 267 | |
---|
| 268 | do i = 1, nSubCols |
---|
| 269 | if(cloudMask(i)) then |
---|
| 270 | ! |
---|
| 271 | ! Cloud top pressure determination |
---|
| 272 | ! MODIS uses CO2 slicing for clouds with tops above about 700 mb and thermal methods for clouds |
---|
| 273 | ! lower than that. |
---|
| 274 | ! For CO2 slicing we report the optical-depth weighted pressure, integrating to a specified |
---|
| 275 | ! optical depth |
---|
| 276 | ! This assumes linear variation in p between levels. Linear in ln(p) is probably better, |
---|
| 277 | ! though we'd need to deal with the lowest pressure gracefully. |
---|
| 278 | ! |
---|
| 279 | retrievedCloudTopPressure(i) = cloud_top_pressure((/ 0., tauTotal(i, :) /), & |
---|
| 280 | pressureLevels, & |
---|
| 281 | CO2Slicing_TauLimit) |
---|
| 282 | |
---|
| 283 | |
---|
| 284 | ! |
---|
| 285 | ! Phase determination - determine fraction of total tau that's liquid |
---|
| 286 | ! When ice and water contribute about equally to the extinction we can't tell |
---|
| 287 | ! what the phase is |
---|
| 288 | ! |
---|
| 289 | integratedLiquidFraction = weight_by_extinction(tauTotal(i, :), & |
---|
| 290 | tauLiquidFraction(i, :), & |
---|
| 291 | phase_TauLimit) |
---|
| 292 | if(integratedLiquidFraction >= phaseDiscrimination_Threshold) then |
---|
| 293 | retrievedPhase(i) = phaseIsLiquid |
---|
| 294 | else if (integratedLiquidFraction <= 1.- phaseDiscrimination_Threshold) then |
---|
| 295 | retrievedPhase(i) = phaseIsIce |
---|
| 296 | else |
---|
| 297 | retrievedPhase(i) = phaseIsUndetermined |
---|
| 298 | end if |
---|
| 299 | |
---|
| 300 | ! |
---|
| 301 | ! Size determination |
---|
| 302 | ! |
---|
| 303 | if(useSimpleReScheme) then |
---|
| 304 | ! This is the extinction-weighted size considering only the phase we've chosen |
---|
| 305 | ! |
---|
| 306 | if(retrievedPhase(i) == phaseIsIce) then |
---|
| 307 | retrievedSize(i) = weight_by_extinction(ice_opticalThickness(i, :), & |
---|
| 308 | iceSize(i, :), & |
---|
| 309 | (1. - integratedLiquidFraction) * size_TauLimit) |
---|
| 310 | |
---|
| 311 | else if(retrievedPhase(i) == phaseIsLiquid) then |
---|
| 312 | retrievedSize(i) = weight_by_extinction(liquid_opticalThickness(i, :), & |
---|
| 313 | waterSize(i, :), & |
---|
| 314 | integratedLiquidFraction * size_TauLimit) |
---|
| 315 | |
---|
| 316 | else |
---|
| 317 | retrievedSize(i) = 0. |
---|
| 318 | end if |
---|
| 319 | else |
---|
| 320 | retrievedSize(i) = 1.0e-06*retrieve_re(retrievedPhase(i), retrievedTau(i), & |
---|
| 321 | obs_Refl_nir = compute_nir_reflectance(liquid_opticalThickness(i, :), waterSize(i, :)*1.0e6, & |
---|
| 322 | ice_opticalThickness(i, :), iceSize(i, :)*1.0e6)) |
---|
| 323 | end if |
---|
| 324 | else |
---|
| 325 | ! |
---|
| 326 | ! Values when we don't think there's a cloud. |
---|
| 327 | ! |
---|
| 328 | retrievedCloudTopPressure(i) = R_UNDEF |
---|
| 329 | retrievedPhase(i) = phaseIsNone |
---|
| 330 | retrievedSize(i) = R_UNDEF |
---|
| 331 | retrievedTau(i) = R_UNDEF |
---|
| 332 | end if |
---|
| 333 | end do |
---|
| 334 | where((retrievedSize(:) < 0.).and.(retrievedSize(:) /= R_UNDEF)) retrievedSize(:) = 1.0e-06*re_fill |
---|
| 335 | |
---|
| 336 | ! We use the ISCCP-derived CTP for low clouds, since the ISCCP simulator ICARUS |
---|
| 337 | ! mimics what MODIS does to first order. |
---|
| 338 | ! Of course, ISCCP cloud top pressures are in mb. |
---|
| 339 | ! |
---|
| 340 | where(cloudMask(:) .and. retrievedCloudTopPressure(:) > CO2Slicing_PressureLimit) & |
---|
| 341 | retrievedCloudTopPressure(:) = isccpCloudTopPressure * 100. |
---|
| 342 | |
---|
| 343 | end subroutine modis_L2_simulator_twoTaus |
---|
| 344 | !------------------------------------------------------------------------------------------------ |
---|
| 345 | ! |
---|
| 346 | ! MODIS simulator: provide a single optical thickness and the cloud ice and liquid contents; |
---|
| 347 | ! we'll partition this into ice and liquid optical thickness and call the full MODIS simulator |
---|
| 348 | ! |
---|
| 349 | subroutine modis_L2_simulator_oneTau( & |
---|
| 350 | temp, pressureLayers, pressureLevels, & |
---|
| 351 | opticalThickness, cloudWater, cloudIce, & |
---|
| 352 | waterSize, iceSize, & |
---|
| 353 | isccpTau, isccpCloudTopPressure, & |
---|
| 354 | retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize) |
---|
| 355 | ! Grid-mean quantities at layer centers, |
---|
| 356 | ! dimension nLayers |
---|
| 357 | real, dimension(:), intent(in ) :: temp, & ! Temperature, K |
---|
| 358 | pressureLayers, & ! Pressure, Pa |
---|
| 359 | pressureLevels ! Pressure at layer edges, Pa (dimension nLayers + 1) |
---|
| 360 | ! Sub-column quantities |
---|
| 361 | ! dimension nLayers, nSubcols |
---|
| 362 | real, dimension(:, :), intent(in ) :: opticalThickness, & ! Layer optical thickness @ 0.67 microns |
---|
| 363 | cloudWater, & ! Cloud water content, arbitrary units |
---|
| 364 | cloudIce ! Cloud water content, same units as cloudWater |
---|
| 365 | real, dimension(:, :), intent(in ) :: waterSize, & ! Cloud drop effective radius, microns |
---|
| 366 | iceSize ! Cloud ice effective radius, microns |
---|
| 367 | |
---|
| 368 | ! Cloud properties retrieved from ISCCP using top_height = 1 |
---|
| 369 | ! dimension nSubcols |
---|
| 370 | |
---|
| 371 | real, dimension(:), intent(in ) :: isccpTau, & ! Column-integrated optical thickness |
---|
| 372 | isccpCloudTopPressure ! ISCCP-retrieved cloud top pressure (Pa) |
---|
| 373 | |
---|
| 374 | ! Properties retrieved by MODIS |
---|
| 375 | ! dimension nSubcols |
---|
| 376 | integer, dimension(:), intent(out) :: retrievedPhase ! liquid/ice/other - integer |
---|
| 377 | real, dimension(:), intent(out) :: retrievedCloudTopPressure, & ! units of pressureLayers |
---|
| 378 | retrievedTau, & ! unitless |
---|
| 379 | retrievedSize ! microns (or whatever units |
---|
| 380 | ! waterSize and iceSize are supplied in) |
---|
| 381 | ! --------------------------------------------------- |
---|
| 382 | ! Local variables |
---|
| 383 | real, dimension(size(opticalThickness, 1), size(opticalThickness, 2)) :: & |
---|
| 384 | liquid_opticalThickness, ice_opticalThickness, tauLiquidFraction |
---|
| 385 | |
---|
| 386 | ! --------------------------------------------------- |
---|
| 387 | |
---|
| 388 | where(cloudIce(:, :) <= 0.) |
---|
| 389 | tauLiquidFraction(:, :) = 1. |
---|
| 390 | elsewhere |
---|
| 391 | where (cloudWater(:, :) <= 0.) |
---|
| 392 | tauLiquidFraction(:, :) = 0. |
---|
| 393 | elsewhere |
---|
| 394 | ! |
---|
| 395 | ! Geometic optics limit - tau as LWP/re (proportional to LWC/re) |
---|
| 396 | ! |
---|
| 397 | tauLiquidFraction(:, :) = (cloudWater(:, :)/waterSize(:, :)) / & |
---|
| 398 | (cloudWater(:, :)/waterSize(:, :) + cloudIce(:, :)/(ice_density * iceSize(:, :)) ) |
---|
| 399 | end where |
---|
| 400 | end where |
---|
| 401 | liquid_opticalThickness(:, :) = tauLiquidFraction(:, :) * opticalThickness(:, :) |
---|
| 402 | ice_opticalThickness (:, :) = opticalThickness(:, :) - liquid_opticalThickness(:, :) |
---|
| 403 | |
---|
| 404 | call modis_L2_simulator_twoTaus(temp, pressureLayers, pressureLevels, & |
---|
| 405 | liquid_opticalThickness, ice_opticalThickness, & |
---|
| 406 | waterSize, iceSize, & |
---|
| 407 | isccpTau, isccpCloudTopPressure, & |
---|
| 408 | retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize) |
---|
| 409 | |
---|
| 410 | end subroutine modis_L2_simulator_oneTau |
---|
[2720] | 411 | |
---|
| 412 | ! ######################################################################################## |
---|
| 413 | subroutine modis_column(nPoints,nSubCols,phase, cloud_top_pressure, optical_thickness, particle_size, & |
---|
| 414 | Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, & |
---|
| 415 | Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, & |
---|
| 416 | Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, & |
---|
| 417 | Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10,& |
---|
| 418 | Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, Cloud_Top_Pressure_Total_Mean, & |
---|
| 419 | Liquid_Water_Path_Mean, Ice_Water_Path_Mean, & |
---|
| 420 | Optical_Thickness_vs_Cloud_Top_Pressure,Optical_Thickness_vs_ReffIce,Optical_Thickness_vs_ReffLiq) |
---|
| 421 | |
---|
| 422 | ! INPUTS |
---|
| 423 | integer,intent(in) :: & |
---|
| 424 | nPoints, & ! Number of horizontal gridpoints |
---|
| 425 | nSubCols ! Number of subcolumns |
---|
| 426 | integer,intent(in), dimension(:,:) :: & |
---|
| 427 | !ds integer,intent(in), dimension(nPoints, nSubCols) :: & |
---|
| 428 | phase |
---|
| 429 | real,intent(in),dimension(:,:) :: & |
---|
| 430 | !ds real,intent(in),dimension(nPoints, nSubCols) :: & |
---|
| 431 | cloud_top_pressure, & |
---|
| 432 | optical_thickness, & |
---|
| 433 | particle_size |
---|
| 434 | |
---|
| 435 | ! OUTPUTS |
---|
| 436 | real,intent(inout),dimension(:) :: & ! |
---|
| 437 | !ds real,intent(inout),dimension(nPoints) :: & ! |
---|
| 438 | Cloud_Fraction_Total_Mean, & ! |
---|
| 439 | Cloud_Fraction_Water_Mean, & ! |
---|
| 440 | Cloud_Fraction_Ice_Mean, & ! |
---|
| 441 | Cloud_Fraction_High_Mean, & ! |
---|
| 442 | Cloud_Fraction_Mid_Mean, & ! |
---|
| 443 | Cloud_Fraction_Low_Mean, & ! |
---|
| 444 | Optical_Thickness_Total_Mean, & ! |
---|
| 445 | Optical_Thickness_Water_Mean, & ! |
---|
| 446 | Optical_Thickness_Ice_Mean, & ! |
---|
| 447 | Optical_Thickness_Total_MeanLog10, & ! |
---|
| 448 | Optical_Thickness_Water_MeanLog10, & ! |
---|
| 449 | Optical_Thickness_Ice_MeanLog10, & ! |
---|
| 450 | Cloud_Particle_Size_Water_Mean, & ! |
---|
| 451 | Cloud_Particle_Size_Ice_Mean, & ! |
---|
| 452 | Cloud_Top_Pressure_Total_Mean, & ! |
---|
| 453 | Liquid_Water_Path_Mean, & ! |
---|
| 454 | Ice_Water_Path_Mean ! |
---|
| 455 | real,intent(inout),dimension(:,:,:) :: & |
---|
| 456 | !ds real,intent(inout),dimension(nPoints,numTauHistogramBins,numPressureHistogramBins) :: & |
---|
| 457 | Optical_Thickness_vs_Cloud_Top_Pressure |
---|
| 458 | real,intent(inout),dimension(:,:,:) :: & |
---|
| 459 | !ds real,intent(inout),dimension(nPoints,numTauHistogramBins,numMODISReffIceBins) :: & |
---|
| 460 | Optical_Thickness_vs_ReffIce |
---|
| 461 | real,intent(inout),dimension(:,:,:) :: & |
---|
| 462 | !ds real,intent(inout),dimension(nPoints,numTauHistogramBins,numMODISReffLiqBins) :: & |
---|
| 463 | Optical_Thickness_vs_ReffLiq |
---|
| 464 | |
---|
| 465 | ! LOCAL VARIABLES |
---|
| 466 | real, parameter :: & |
---|
| 467 | LWP_conversion = 2./3. * 1000. ! MKS units |
---|
| 468 | integer :: i, j |
---|
| 469 | logical, dimension(nPoints,nSubCols) :: & |
---|
| 470 | cloudMask, & |
---|
| 471 | waterCloudMask, & |
---|
| 472 | iceCloudMask, & |
---|
| 473 | validRetrievalMask |
---|
| 474 | real,dimension(nPoints,nSubCols) :: & |
---|
| 475 | tauWRK,ctpWRK,reffIceWRK,reffLiqWRK |
---|
| 476 | |
---|
| 477 | ! ######################################################################################## |
---|
| 478 | ! Include only those pixels with successful retrievals in the statistics |
---|
| 479 | ! ######################################################################################## |
---|
| 480 | validRetrievalMask(1:nPoints,1:nSubCols) = particle_size(1:nPoints,1:nSubCols) > 0. |
---|
| 481 | cloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) /= phaseIsNone .and. & |
---|
| 482 | validRetrievalMask(1:nPoints,1:nSubCols) |
---|
| 483 | waterCloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) == phaseIsLiquid .and. & |
---|
| 484 | validRetrievalMask(1:nPoints,1:nSubCols) |
---|
| 485 | iceCloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) == phaseIsIce .and. & |
---|
| 486 | validRetrievalMask(1:nPoints,1:nSubCols) |
---|
| 487 | |
---|
| 488 | ! ######################################################################################## |
---|
| 489 | ! Use these as pixel counts at first |
---|
| 490 | ! ######################################################################################## |
---|
| 491 | Cloud_Fraction_Total_Mean(1:nPoints) = real(count(cloudMask, dim = 2)) |
---|
| 492 | Cloud_Fraction_Water_Mean(1:nPoints) = real(count(waterCloudMask, dim = 2)) |
---|
| 493 | Cloud_Fraction_Ice_Mean(1:nPoints) = real(count(iceCloudMask, dim = 2)) |
---|
| 494 | Cloud_Fraction_High_Mean(1:nPoints) = real(count(cloudMask .and. cloud_top_pressure <= & |
---|
| 495 | highCloudPressureLimit, dim = 2)) |
---|
| 496 | Cloud_Fraction_Low_Mean(1:nPoints) = real(count(cloudMask .and. cloud_top_pressure > & |
---|
| 497 | lowCloudPressureLimit, dim = 2)) |
---|
| 498 | Cloud_Fraction_Mid_Mean(1:nPoints) = Cloud_Fraction_Total_Mean(1:nPoints) - Cloud_Fraction_High_Mean(1:nPoints)& |
---|
| 499 | - Cloud_Fraction_Low_Mean(1:nPoints) |
---|
| 500 | |
---|
| 501 | ! ######################################################################################## |
---|
| 502 | ! Compute mean optical thickness. |
---|
| 503 | ! ######################################################################################## |
---|
| 504 | Optical_Thickness_Total_Mean(1:nPoints) = sum(optical_thickness, mask = cloudMask, dim = 2) / & |
---|
| 505 | Cloud_Fraction_Total_Mean(1:nPoints) |
---|
| 506 | Optical_Thickness_Water_Mean(1:nPoints) = sum(optical_thickness, mask = waterCloudMask, dim = 2) / & |
---|
| 507 | Cloud_Fraction_Water_Mean(1:nPoints) |
---|
| 508 | Optical_Thickness_Ice_Mean(1:nPoints) = sum(optical_thickness, mask = iceCloudMask, dim = 2) / & |
---|
| 509 | Cloud_Fraction_Ice_Mean(1:nPoints) |
---|
| 510 | |
---|
| 511 | ! ######################################################################################## |
---|
| 512 | ! We take the absolute value of optical thickness here to satisfy compilers that complains |
---|
| 513 | ! when we evaluate the logarithm of a negative number, even though it's not included in |
---|
| 514 | ! the sum. |
---|
| 515 | ! ######################################################################################## |
---|
| 516 | Optical_Thickness_Total_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = cloudMask, & |
---|
| 517 | dim = 2) / Cloud_Fraction_Total_Mean(1:nPoints) |
---|
| 518 | Optical_Thickness_Water_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = waterCloudMask,& |
---|
| 519 | dim = 2) / Cloud_Fraction_Water_Mean(1:nPoints) |
---|
| 520 | Optical_Thickness_Ice_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = iceCloudMask,& |
---|
| 521 | dim = 2) / Cloud_Fraction_Ice_Mean(1:nPoints) |
---|
| 522 | Cloud_Particle_Size_Water_Mean(1:nPoints) = sum(particle_size, mask = waterCloudMask, dim = 2) / & |
---|
| 523 | Cloud_Fraction_Water_Mean(1:nPoints) |
---|
| 524 | Cloud_Particle_Size_Ice_Mean(1:nPoints) = sum(particle_size, mask = iceCloudMask, dim = 2) / & |
---|
| 525 | Cloud_Fraction_Ice_Mean(1:nPoints) |
---|
| 526 | Cloud_Top_Pressure_Total_Mean(1:nPoints) = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / & |
---|
| 527 | max(1, count(cloudMask, dim = 2)) |
---|
| 528 | Liquid_Water_Path_Mean(1:nPoints) = LWP_conversion*sum(particle_size*optical_thickness, & |
---|
| 529 | mask=waterCloudMask,dim=2)/Cloud_Fraction_Water_Mean(1:nPoints) |
---|
| 530 | Ice_Water_Path_Mean(1:nPoints) = LWP_conversion * ice_density*sum(particle_size*optical_thickness,& |
---|
| 531 | mask=iceCloudMask,dim = 2) /Cloud_Fraction_Ice_Mean(1:nPoints) |
---|
| 532 | |
---|
| 533 | ! ######################################################################################## |
---|
[2839] | 534 | ! Normalize pixel counts to fraction. |
---|
[2720] | 535 | ! ######################################################################################## |
---|
[2839] | 536 | Cloud_Fraction_High_Mean(1:nPoints) = Cloud_Fraction_High_Mean(1:nPoints) /nSubcols |
---|
| 537 | Cloud_Fraction_Mid_Mean(1:nPoints) = Cloud_Fraction_Mid_Mean(1:nPoints) /nSubcols |
---|
| 538 | Cloud_Fraction_Low_Mean(1:nPoints) = Cloud_Fraction_Low_Mean(1:nPoints) /nSubcols |
---|
| 539 | Cloud_Fraction_Total_Mean(1:nPoints) = Cloud_Fraction_Total_Mean(1:nPoints) /nSubcols |
---|
| 540 | Cloud_Fraction_Ice_Mean(1:nPoints) = Cloud_Fraction_Ice_Mean(1:nPoints) /nSubcols |
---|
| 541 | Cloud_Fraction_Water_Mean(1:nPoints) = Cloud_Fraction_Water_Mean(1:nPoints) /nSubcols |
---|
| 542 | |
---|
[2720] | 543 | ! ######################################################################################## |
---|
| 544 | ! Set clear-scenes to undefined |
---|
| 545 | ! ######################################################################################## |
---|
| 546 | where (Cloud_Fraction_Total_Mean == 0) |
---|
| 547 | Optical_Thickness_Total_Mean = R_UNDEF |
---|
| 548 | Optical_Thickness_Total_MeanLog10 = R_UNDEF |
---|
| 549 | Cloud_Top_Pressure_Total_Mean = R_UNDEF |
---|
| 550 | endwhere |
---|
| 551 | where (Cloud_Fraction_Water_Mean == 0) |
---|
| 552 | Optical_Thickness_Water_Mean = R_UNDEF |
---|
| 553 | Optical_Thickness_Water_MeanLog10 = R_UNDEF |
---|
| 554 | Cloud_Particle_Size_Water_Mean = R_UNDEF |
---|
| 555 | Liquid_Water_Path_Mean = R_UNDEF |
---|
| 556 | endwhere |
---|
| 557 | where (Cloud_Fraction_Ice_Mean == 0) |
---|
| 558 | Optical_Thickness_Ice_Mean = R_UNDEF |
---|
| 559 | Optical_Thickness_Ice_MeanLog10 = R_UNDEF |
---|
| 560 | Cloud_Particle_Size_Ice_Mean = R_UNDEF |
---|
| 561 | Ice_Water_Path_Mean = R_UNDEF |
---|
| 562 | endwhere |
---|
| 563 | where (Cloud_Fraction_High_Mean == 0) Cloud_Fraction_High_Mean = R_UNDEF |
---|
| 564 | where (Cloud_Fraction_Mid_Mean == 0) Cloud_Fraction_Mid_Mean = R_UNDEF |
---|
| 565 | where (Cloud_Fraction_Low_Mean == 0) Cloud_Fraction_Low_Mean = R_UNDEF |
---|
| 566 | |
---|
| 567 | ! ######################################################################################## |
---|
| 568 | ! Joint histogram |
---|
| 569 | ! ######################################################################################## |
---|
| 570 | |
---|
| 571 | ! Loop over all points |
---|
| 572 | tauWRK(1:nPoints,1:nSubCols) = optical_thickness(1:nPoints,1:nSubCols) |
---|
| 573 | ctpWRK(1:nPoints,1:nSubCols) = cloud_top_pressure(1:nPoints,1:nSubCols) |
---|
| 574 | reffIceWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,iceCloudMask) |
---|
| 575 | reffLiqWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,waterCloudMask) |
---|
| 576 | do j=1,nPoints |
---|
| 577 | |
---|
| 578 | ! Fill clear and optically thin subcolumns with fill |
---|
| 579 | where(.not. cloudMask(j,1:nSubCols)) |
---|
| 580 | tauWRK(j,1:nSubCols) = -999. |
---|
| 581 | ctpWRK(j,1:nSubCols) = -999. |
---|
| 582 | endwhere |
---|
| 583 | ! Joint histogram of tau/CTP |
---|
| 584 | call hist2D(tauWRK(j,1:nSubCols),ctpWRK(j,1:nSubCols),nSubCols,& |
---|
| 585 | tauHistogramBoundaries,numTauHistogramBins,& |
---|
| 586 | pressureHistogramBoundaries,numPressureHistogramBins,& |
---|
| 587 | Optical_Thickness_vs_Cloud_Top_Pressure(j,1:numTauHistogramBins,1:numPressureHistogramBins)) |
---|
| 588 | ! Joint histogram of tau/ReffICE |
---|
| 589 | call hist2D(tauWRK(j,1:nSubCols),reffIceWrk(j,1:nSubCols),nSubCols, & |
---|
| 590 | tauHistogramBoundaries,numTauHistogramBins,reffICE_binBounds, & |
---|
| 591 | numMODISReffIceBins, Optical_Thickness_vs_ReffIce(j,1:numTauHistogramBins,1:numMODISReffIceBins)) |
---|
| 592 | ! Joint histogram of tau/ReffLIQ |
---|
| 593 | call hist2D(tauWRK(j,1:nSubCols),reffLiqWrk(j,1:nSubCols),nSubCols, & |
---|
| 594 | tauHistogramBoundaries,numTauHistogramBins,reffLIQ_binBounds, & |
---|
| 595 | numMODISReffLiqBins, Optical_Thickness_vs_ReffLiq(j,1:numTauHistogramBins,1:numMODISReffLiqBins)) |
---|
| 596 | |
---|
| 597 | enddo |
---|
| 598 | Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numTauHistogramBins,1:numPressureHistogramBins) = & |
---|
| 599 | Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numTauHistogramBins,1:numPressureHistogramBins)/nSubCols |
---|
| 600 | Optical_Thickness_vs_ReffIce(1:nPoints,1:numTauHistogramBins,1:numMODISReffIceBins) = & |
---|
| 601 | Optical_Thickness_vs_ReffIce(1:nPoints,1:numTauHistogramBins,1:numMODISReffIceBins)/nSubCols |
---|
| 602 | Optical_Thickness_vs_ReffLiq(1:nPoints,1:numTauHistogramBins,1:numMODISReffLiqBins) = & |
---|
| 603 | Optical_Thickness_vs_ReffLiq(1:nPoints,1:numTauHistogramBins,1:numMODISReffLiqBins)/nSubCols |
---|
| 604 | |
---|
| 605 | end subroutine modis_column |
---|
| 606 | ! ###################################################################################### |
---|
| 607 | ! SUBROUTINE hist2D |
---|
| 608 | ! ###################################################################################### |
---|
| 609 | subroutine hist2D(var1,var2,npts,bin1,nbin1,bin2,nbin2,jointHist) |
---|
| 610 | implicit none |
---|
| 611 | |
---|
| 612 | ! INPUTS |
---|
| 613 | integer, intent(in) :: & |
---|
| 614 | npts, & ! Number of data points to be sorted |
---|
| 615 | nbin1, & ! Number of bins in histogram direction 1 |
---|
| 616 | nbin2 ! Number of bins in histogram direction 2 |
---|
| 617 | real,intent(in),dimension(npts) :: & |
---|
| 618 | var1, & ! Variable 1 to be sorted into bins |
---|
| 619 | var2 ! variable 2 to be sorted into bins |
---|
| 620 | real,intent(in),dimension(nbin1+1) :: & |
---|
| 621 | bin1 ! Histogram bin 1 boundaries |
---|
| 622 | real,intent(in),dimension(nbin2+1) :: & |
---|
| 623 | bin2 ! Histogram bin 2 boundaries |
---|
| 624 | ! OUTPUTS |
---|
| 625 | real,intent(out),dimension(nbin1,nbin2) :: & |
---|
| 626 | jointHist |
---|
| 627 | |
---|
| 628 | ! LOCAL VARIABLES |
---|
| 629 | integer :: ij,ik |
---|
| 630 | |
---|
| 631 | do ij=2,nbin1+1 |
---|
| 632 | do ik=2,nbin2+1 |
---|
| 633 | jointHist(ij-1,ik-1)=count(var1 .ge. bin1(ij-1) .and. var1 .lt. bin1(ij) .and. & |
---|
| 634 | var2 .ge. bin2(ik-1) .and. var2 .lt. bin2(ik)) |
---|
| 635 | enddo |
---|
| 636 | enddo |
---|
| 637 | end subroutine hist2D |
---|
| 638 | |
---|
[2432] | 639 | !------------------------------------------------------------------------------------------------ |
---|
| 640 | subroutine modis_L3_simulator(phase, cloud_top_pressure, optical_thickness, particle_size, & |
---|
| 641 | Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, & |
---|
| 642 | Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, & |
---|
| 643 | Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, & |
---|
| 644 | Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10, & |
---|
| 645 | Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, & |
---|
| 646 | Cloud_Top_Pressure_Total_Mean, & |
---|
| 647 | Liquid_Water_Path_Mean, Ice_Water_Path_Mean, & |
---|
| 648 | Optical_Thickness_vs_Cloud_Top_Pressure) |
---|
| 649 | ! |
---|
| 650 | ! Inputs; dimension nPoints, nSubcols |
---|
| 651 | ! |
---|
| 652 | integer, dimension(:, :), intent(in) :: phase |
---|
| 653 | real, dimension(:, :), intent(in) :: cloud_top_pressure, optical_thickness, particle_size |
---|
| 654 | ! |
---|
| 655 | ! Outputs; dimension nPoints |
---|
| 656 | ! |
---|
| 657 | real, dimension(:), intent(out) :: & |
---|
| 658 | Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, & |
---|
| 659 | Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, & |
---|
| 660 | Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, & |
---|
| 661 | Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10, & |
---|
| 662 | Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, & |
---|
| 663 | Cloud_Top_Pressure_Total_Mean, & |
---|
| 664 | Liquid_Water_Path_Mean, Ice_Water_Path_Mean |
---|
| 665 | ! tau/ctp histogram; dimensions nPoints, numTauHistogramBins , numPressureHistogramBins |
---|
| 666 | real, dimension(:, :, :), intent(out) :: Optical_Thickness_vs_Cloud_Top_Pressure |
---|
| 667 | ! --------------------------- |
---|
| 668 | ! Local variables |
---|
| 669 | ! |
---|
| 670 | real, parameter :: LWP_conversion = 2./3. * 1000. ! MKS units |
---|
| 671 | integer :: i, j |
---|
| 672 | integer :: nPoints, nSubcols |
---|
| 673 | logical, dimension(size(phase, 1), size(phase, 2)) :: & |
---|
| 674 | cloudMask, waterCloudMask, iceCloudMask, validRetrievalMask |
---|
| 675 | logical, dimension(size(phase, 1), size(phase, 2), numTauHistogramBins ) :: tauMask |
---|
| 676 | logical, dimension(size(phase, 1), size(phase, 2), numPressureHistogramBins) :: pressureMask |
---|
| 677 | ! --------------------------- |
---|
| 678 | |
---|
| 679 | nPoints = size(phase, 1) |
---|
| 680 | nSubcols = size(phase, 2) |
---|
| 681 | ! |
---|
| 682 | ! Array conformance checks |
---|
| 683 | ! |
---|
| 684 | if(any( (/ size(cloud_top_pressure, 1), size(optical_thickness, 1), size(particle_size, 1), & |
---|
| 685 | size(Cloud_Fraction_Total_Mean), size(Cloud_Fraction_Water_Mean), size(Cloud_Fraction_Ice_Mean), & |
---|
| 686 | size(Cloud_Fraction_High_Mean), size(Cloud_Fraction_Mid_Mean), size(Cloud_Fraction_Low_Mean), & |
---|
| 687 | size(Optical_Thickness_Total_Mean), size(Optical_Thickness_Water_Mean), size(Optical_Thickness_Ice_Mean), & |
---|
| 688 | size(Optical_Thickness_Total_MeanLog10), size(Optical_Thickness_Water_MeanLog10), & |
---|
| 689 | size(Optical_Thickness_Ice_MeanLog10), size(Cloud_Particle_Size_Water_Mean), & |
---|
| 690 | size(Cloud_Particle_Size_Ice_Mean), size(Cloud_Top_Pressure_Total_Mean), & |
---|
| 691 | size(Liquid_Water_Path_Mean), size(Ice_Water_Path_Mean) /) /= nPoints)) & |
---|
| 692 | call complain_and_die("Some L3 arrays have wrong number of grid points") |
---|
| 693 | if(any( (/ size(cloud_top_pressure, 2), size(optical_thickness, 2), size(particle_size, 2) /) /= nSubcols)) & |
---|
| 694 | call complain_and_die("Some L3 arrays have wrong number of subcolumns") |
---|
| 695 | |
---|
| 696 | |
---|
| 697 | ! |
---|
| 698 | ! Include only those pixels with successful retrievals in the statistics |
---|
| 699 | ! |
---|
| 700 | validRetrievalMask(:, :) = particle_size(:, :) > 0. |
---|
| 701 | cloudMask = phase(:, :) /= phaseIsNone .and. validRetrievalMask(:, :) |
---|
| 702 | waterCloudMask = phase(:, :) == phaseIsLiquid .and. validRetrievalMask(:, :) |
---|
| 703 | iceCloudMask = phase(:, :) == phaseIsIce .and. validRetrievalMask(:, :) |
---|
| 704 | ! |
---|
| 705 | ! Use these as pixel counts at first |
---|
| 706 | ! |
---|
| 707 | Cloud_Fraction_Total_Mean(:) = real(count(cloudMask, dim = 2)) |
---|
| 708 | Cloud_Fraction_Water_Mean(:) = real(count(waterCloudMask, dim = 2)) |
---|
| 709 | Cloud_Fraction_Ice_Mean(:) = real(count(iceCloudMask, dim = 2)) |
---|
| 710 | |
---|
| 711 | Cloud_Fraction_High_Mean(:) = real(count(cloudMask .and. cloud_top_pressure <= highCloudPressureLimit, dim = 2)) |
---|
| 712 | Cloud_Fraction_Low_Mean(:) = real(count(cloudMask .and. cloud_top_pressure > lowCloudPressureLimit, dim = 2)) |
---|
| 713 | Cloud_Fraction_Mid_Mean(:) = Cloud_Fraction_Total_Mean(:) - Cloud_Fraction_High_Mean(:) - Cloud_Fraction_Low_Mean(:) |
---|
| 714 | |
---|
| 715 | ! |
---|
| 716 | ! Don't want to divide by 0, even though the sums will be 0 where the pixel counts are 0. |
---|
| 717 | ! |
---|
| 718 | where (Cloud_Fraction_Total_Mean == 0) Cloud_Fraction_Total_Mean = -1. |
---|
| 719 | where (Cloud_Fraction_Water_Mean == 0) Cloud_Fraction_Water_Mean = -1. |
---|
| 720 | where (Cloud_Fraction_Ice_Mean == 0) Cloud_Fraction_Ice_Mean = -1. |
---|
| 721 | |
---|
| 722 | Optical_Thickness_Total_Mean = sum(optical_thickness, mask = cloudMask, dim = 2) / Cloud_Fraction_Total_Mean(:) |
---|
| 723 | Optical_Thickness_Water_Mean = sum(optical_thickness, mask = waterCloudMask, dim = 2) / Cloud_Fraction_Water_Mean(:) |
---|
| 724 | Optical_Thickness_Ice_Mean = sum(optical_thickness, mask = iceCloudMask, dim = 2) / Cloud_Fraction_Ice_Mean(:) |
---|
| 725 | |
---|
| 726 | ! We take the absolute value of optical thickness here to satisfy compilers that complains when we |
---|
| 727 | ! evaluate the logarithm of a negative number, even though it's not included in the sum. |
---|
| 728 | Optical_Thickness_Total_MeanLog10 = sum(log10(abs(optical_thickness)), mask = cloudMask, dim = 2) / & |
---|
| 729 | Cloud_Fraction_Total_Mean(:) |
---|
| 730 | Optical_Thickness_Water_MeanLog10 = sum(log10(abs(optical_thickness)), mask = waterCloudMask, dim = 2) / & |
---|
| 731 | Cloud_Fraction_Water_Mean(:) |
---|
| 732 | Optical_Thickness_Ice_MeanLog10 = sum(log10(abs(optical_thickness)), mask = iceCloudMask, dim = 2) / & |
---|
| 733 | Cloud_Fraction_Ice_Mean(:) |
---|
| 734 | |
---|
| 735 | Cloud_Particle_Size_Water_Mean = sum(particle_size, mask = waterCloudMask, dim = 2) / Cloud_Fraction_Water_Mean(:) |
---|
| 736 | Cloud_Particle_Size_Ice_Mean = sum(particle_size, mask = iceCloudMask, dim = 2) / Cloud_Fraction_Ice_Mean(:) |
---|
| 737 | |
---|
| 738 | Cloud_Top_Pressure_Total_Mean = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / max(1, count(cloudMask, dim = 2)) |
---|
| 739 | |
---|
| 740 | Liquid_Water_Path_Mean = LWP_conversion & |
---|
| 741 | * sum(particle_size * optical_thickness, mask = waterCloudMask, dim = 2) & |
---|
| 742 | / Cloud_Fraction_Water_Mean(:) |
---|
| 743 | Ice_Water_Path_Mean = LWP_conversion * ice_density & |
---|
| 744 | * sum(particle_size * optical_thickness, mask = iceCloudMask, dim = 2) & |
---|
| 745 | / Cloud_Fraction_Ice_Mean(:) |
---|
| 746 | |
---|
| 747 | ! |
---|
| 748 | ! Normalize pixel counts to fraction |
---|
| 749 | ! The first three cloud fractions have been set to -1 in cloud-free areas, so set those places to 0. |
---|
| 750 | ! |
---|
| 751 | Cloud_Fraction_Total_Mean(:) = max(0., Cloud_Fraction_Total_Mean(:)/nSubcols) |
---|
| 752 | Cloud_Fraction_Water_Mean(:) = max(0., Cloud_Fraction_Water_Mean(:)/nSubcols) |
---|
| 753 | Cloud_Fraction_Ice_Mean(:) = max(0., Cloud_Fraction_Ice_Mean(:) /nSubcols) |
---|
| 754 | |
---|
| 755 | Cloud_Fraction_High_Mean(:) = Cloud_Fraction_High_Mean(:) /nSubcols |
---|
| 756 | Cloud_Fraction_Mid_Mean(:) = Cloud_Fraction_Mid_Mean(:) /nSubcols |
---|
| 757 | Cloud_Fraction_Low_Mean(:) = Cloud_Fraction_Low_Mean(:) /nSubcols |
---|
| 758 | |
---|
| 759 | ! ---- |
---|
| 760 | ! Joint histogram |
---|
| 761 | ! |
---|
| 762 | do i = 1, numTauHistogramBins |
---|
| 763 | where(cloudMask(:, :)) |
---|
| 764 | tauMask(:, :, i) = optical_thickness(:, :) >= tauHistogramBoundaries(i) .and. & |
---|
| 765 | optical_thickness(:, :) < tauHistogramBoundaries(i+1) |
---|
| 766 | elsewhere |
---|
| 767 | tauMask(:, :, i) = .false. |
---|
| 768 | end where |
---|
| 769 | end do |
---|
| 770 | |
---|
| 771 | do i = 1, numPressureHistogramBins |
---|
| 772 | where(cloudMask(:, :)) |
---|
| 773 | pressureMask(:, :, i) = cloud_top_pressure(:, :) >= pressureHistogramBoundaries(i) .and. & |
---|
| 774 | cloud_top_pressure(:, :) < pressureHistogramBoundaries(i+1) |
---|
| 775 | elsewhere |
---|
| 776 | pressureMask(:, :, i) = .false. |
---|
| 777 | end where |
---|
| 778 | end do |
---|
| 779 | |
---|
| 780 | do i = 1, numPressureHistogramBins |
---|
| 781 | do j = 1, numTauHistogramBins |
---|
| 782 | Optical_Thickness_vs_Cloud_Top_Pressure(:, j, i) = & |
---|
| 783 | real(count(tauMask(:, :, j) .and. pressureMask(:, :, i), dim = 2)) / real(nSubcols) |
---|
| 784 | end do |
---|
| 785 | end do |
---|
| 786 | |
---|
| 787 | end subroutine modis_L3_simulator |
---|
| 788 | !------------------------------------------------------------------------------------------------ |
---|
| 789 | function cloud_top_pressure(tauIncrement, pressure, tauLimit) |
---|
| 790 | real, dimension(:), intent(in) :: tauIncrement, pressure |
---|
| 791 | real, intent(in) :: tauLimit |
---|
| 792 | real :: cloud_top_pressure |
---|
| 793 | ! |
---|
| 794 | ! Find the extinction-weighted pressure. Assume that pressure varies linearly between |
---|
| 795 | ! layers and use the trapezoidal rule. |
---|
| 796 | ! |
---|
| 797 | |
---|
| 798 | real :: deltaX, totalTau, totalProduct |
---|
| 799 | integer :: i |
---|
| 800 | |
---|
| 801 | totalTau = 0.; totalProduct = 0. |
---|
| 802 | do i = 2, size(tauIncrement) |
---|
| 803 | if(totalTau + tauIncrement(i) > tauLimit) then |
---|
| 804 | deltaX = tauLimit - totalTau |
---|
| 805 | totalTau = totalTau + deltaX |
---|
| 806 | ! |
---|
| 807 | ! Result for trapezoidal rule when you take less than a full step |
---|
| 808 | ! tauIncrement is a layer-integrated value |
---|
| 809 | ! |
---|
| 810 | totalProduct = totalProduct & |
---|
| 811 | + pressure(i-1) * deltaX & |
---|
| 812 | + (pressure(i) - pressure(i-1)) * deltaX**2/(2. * tauIncrement(i)) |
---|
| 813 | else |
---|
| 814 | totalTau = totalTau + tauIncrement(i) |
---|
| 815 | totalProduct = totalProduct + tauIncrement(i) * (pressure(i) + pressure(i-1)) / 2. |
---|
| 816 | end if |
---|
| 817 | if(totalTau >= tauLimit) exit |
---|
| 818 | end do |
---|
| 819 | cloud_top_pressure = totalProduct/totalTau |
---|
| 820 | end function cloud_top_pressure |
---|
| 821 | !------------------------------------------------------------------------------------------------ |
---|
| 822 | function weight_by_extinction(tauIncrement, f, tauLimit) |
---|
| 823 | real, dimension(:), intent(in) :: tauIncrement, f |
---|
| 824 | real, intent(in) :: tauLimit |
---|
| 825 | real :: weight_by_extinction |
---|
| 826 | ! |
---|
| 827 | ! Find the extinction-weighted value of f(tau), assuming constant f within each layer |
---|
| 828 | ! |
---|
| 829 | |
---|
| 830 | real :: deltaX, totalTau, totalProduct |
---|
| 831 | integer :: i |
---|
| 832 | |
---|
| 833 | totalTau = 0.; totalProduct = 0. |
---|
| 834 | do i = 1, size(tauIncrement) |
---|
| 835 | if(totalTau + tauIncrement(i) > tauLimit) then |
---|
| 836 | deltaX = tauLimit - totalTau |
---|
| 837 | totalTau = totalTau + deltaX |
---|
| 838 | totalProduct = totalProduct + deltaX * f(i) |
---|
| 839 | else |
---|
| 840 | totalTau = totalTau + tauIncrement(i) |
---|
| 841 | totalProduct = totalProduct + tauIncrement(i) * f(i) |
---|
| 842 | end if |
---|
| 843 | if(totalTau >= tauLimit) exit |
---|
| 844 | end do |
---|
| 845 | weight_by_extinction = totalProduct/totalTau |
---|
| 846 | end function weight_by_extinction |
---|
| 847 | !------------------------------------------------------------------------------------------------ |
---|
| 848 | pure function compute_nir_reflectance(water_tau, water_size, ice_tau, ice_size) |
---|
| 849 | real, dimension(:), intent(in) :: water_tau, water_size, ice_tau, ice_size |
---|
| 850 | real :: compute_nir_reflectance |
---|
| 851 | |
---|
| 852 | real, dimension(size(water_tau)) :: water_g, water_w0, ice_g, ice_w0, & |
---|
| 853 | tau, g, w0 |
---|
| 854 | !---------------------------------------- |
---|
| 855 | water_g(:) = get_g_nir( phaseIsLiquid, water_size) |
---|
| 856 | water_w0(:) = get_ssa_nir(phaseIsLiquid, water_size) |
---|
| 857 | ice_g(:) = get_g_nir( phaseIsIce, ice_size) |
---|
| 858 | ice_w0(:) = get_ssa_nir(phaseIsIce, ice_size) |
---|
| 859 | ! |
---|
| 860 | ! Combine ice and water optical properties |
---|
| 861 | ! |
---|
| 862 | g(:) = 0; w0(:) = 0. |
---|
| 863 | tau(:) = ice_tau(:) + water_tau(:) |
---|
| 864 | where (tau(:) > 0) |
---|
| 865 | g(:) = (water_tau(:) * water_g(:) + ice_tau(:) * ice_g(:) ) / & |
---|
| 866 | tau(:) |
---|
| 867 | w0(:) = (water_tau(:) * water_g(:) * water_w0(:) + ice_tau(:) * ice_g(:) * ice_w0(:)) / & |
---|
| 868 | (g(:) * tau(:)) |
---|
| 869 | end where |
---|
| 870 | |
---|
| 871 | compute_nir_reflectance = compute_toa_reflectace(tau, g, w0) |
---|
| 872 | end function compute_nir_reflectance |
---|
| 873 | !------------------------------------------------------------------------------------------------ |
---|
| 874 | ! Retreivals |
---|
| 875 | !------------------------------------------------------------------------------------------------ |
---|
| 876 | elemental function retrieve_re (phase, tau, obs_Refl_nir) |
---|
| 877 | integer, intent(in) :: phase |
---|
| 878 | real, intent(in) :: tau, obs_Refl_nir |
---|
| 879 | real :: retrieve_re |
---|
| 880 | ! |
---|
| 881 | ! Finds the re that produces the minimum mis-match between predicted and observed reflectance in |
---|
| 882 | ! MODIS band 7 (near IR) |
---|
| 883 | ! Uses |
---|
| 884 | ! fits for asymmetry parameter g(re) and single scattering albedo w0(re) based on MODIS tables |
---|
| 885 | ! two-stream for layer reflectance and transmittance as a function of optical thickness tau, g, and w0 |
---|
| 886 | ! adding-doubling for total reflectance |
---|
| 887 | ! |
---|
| 888 | ! |
---|
| 889 | ! |
---|
| 890 | ! Local variables |
---|
| 891 | ! |
---|
| 892 | real, parameter :: min_distance_to_boundary = 0.01 |
---|
| 893 | real :: re_min, re_max, delta_re |
---|
| 894 | integer :: i |
---|
| 895 | |
---|
| 896 | real, dimension(num_trial_res) :: trial_re, g, w0, predicted_Refl_nir |
---|
| 897 | ! -------------------------- |
---|
| 898 | |
---|
| 899 | if(any(phase == (/ phaseIsLiquid, phaseIsUndetermined, phaseIsIce /))) then |
---|
| 900 | if (phase == phaseIsLiquid .OR. phase == phaseIsUndetermined) then |
---|
| 901 | re_min = re_water_min |
---|
| 902 | re_max = re_water_max |
---|
| 903 | trial_re(:) = trial_re_w |
---|
| 904 | g(:) = g_w(:) |
---|
| 905 | w0(:) = w0_w(:) |
---|
| 906 | else |
---|
| 907 | re_min = re_ice_min |
---|
| 908 | re_max = re_ice_max |
---|
| 909 | trial_re(:) = trial_re_i |
---|
| 910 | g(:) = g_i(:) |
---|
| 911 | w0(:) = w0_i(:) |
---|
| 912 | end if |
---|
| 913 | ! |
---|
| 914 | ! 1st attempt at index: w/coarse re resolution |
---|
| 915 | ! |
---|
| 916 | predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:)) |
---|
| 917 | retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir) |
---|
| 918 | ! |
---|
| 919 | ! If first retrieval works, can try 2nd iteration using greater re resolution |
---|
| 920 | ! |
---|
[2720] | 921 | ! DJS2015: Remove unused piece of code |
---|
| 922 | ! if(use_two_re_iterations .and. retrieve_re > 0.) then |
---|
| 923 | ! re_min = retrieve_re - delta_re |
---|
| 924 | ! re_max = retrieve_re + delta_re |
---|
| 925 | ! delta_re = (re_max - re_min)/real(num_trial_res-1) |
---|
| 926 | ! |
---|
| 927 | ! trial_re(:) = re_min + delta_re * (/ (i - 1, i = 1, num_trial_res) /) |
---|
| 928 | ! g(:) = get_g_nir( phase, trial_re(:)) |
---|
| 929 | ! w0(:) = get_ssa_nir(phase, trial_re(:)) |
---|
| 930 | ! predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:)) |
---|
| 931 | ! retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir) |
---|
| 932 | ! end if |
---|
| 933 | ! DJS2015 END |
---|
[2432] | 934 | else |
---|
| 935 | retrieve_re = re_fill |
---|
| 936 | end if |
---|
| 937 | |
---|
| 938 | end function retrieve_re |
---|
| 939 | ! -------------------------------------------- |
---|
| 940 | pure function interpolate_to_min(x, y, yobs) |
---|
| 941 | real, dimension(:), intent(in) :: x, y |
---|
| 942 | real, intent(in) :: yobs |
---|
| 943 | real :: interpolate_to_min |
---|
| 944 | ! |
---|
| 945 | ! Given a set of values of y as y(x), find the value of x that minimizes abs(y - yobs) |
---|
| 946 | ! y must be monotonic in x |
---|
| 947 | ! |
---|
| 948 | real, dimension(size(x)) :: diff |
---|
| 949 | integer :: nPoints, minDiffLoc, lowerBound, upperBound |
---|
| 950 | ! --------------------------------- |
---|
| 951 | nPoints = size(y) |
---|
| 952 | diff(:) = y(:) - yobs |
---|
| 953 | minDiffLoc = minloc(abs(diff), dim = 1) |
---|
| 954 | |
---|
| 955 | if(minDiffLoc == 1) then |
---|
| 956 | lowerBound = minDiffLoc |
---|
| 957 | upperBound = minDiffLoc + 1 |
---|
| 958 | else if(minDiffLoc == nPoints) then |
---|
| 959 | lowerBound = minDiffLoc - 1 |
---|
| 960 | upperBound = minDiffLoc |
---|
| 961 | else |
---|
| 962 | if(diff(minDiffLoc-1) * diff(minDiffLoc) < 0) then |
---|
| 963 | lowerBound = minDiffLoc-1 |
---|
| 964 | upperBound = minDiffLoc |
---|
| 965 | else |
---|
| 966 | lowerBound = minDiffLoc |
---|
| 967 | upperBound = minDiffLoc + 1 |
---|
| 968 | end if |
---|
| 969 | end if |
---|
| 970 | |
---|
| 971 | if(diff(lowerBound) * diff(upperBound) < 0) then |
---|
| 972 | ! |
---|
| 973 | ! Interpolate the root position linearly if we bracket the root |
---|
| 974 | ! |
---|
| 975 | interpolate_to_min = x(upperBound) - & |
---|
| 976 | diff(upperBound) * (x(upperBound) - x(lowerBound)) / (diff(upperBound) - diff(lowerBound)) |
---|
| 977 | else |
---|
| 978 | interpolate_to_min = re_fill |
---|
| 979 | end if |
---|
| 980 | |
---|
| 981 | |
---|
| 982 | end function interpolate_to_min |
---|
| 983 | ! -------------------------------------------- |
---|
| 984 | ! Optical properties |
---|
| 985 | ! -------------------------------------------- |
---|
| 986 | elemental function get_g_nir (phase, re) |
---|
| 987 | ! |
---|
| 988 | ! Polynomial fit for asummetry parameter g in MODIS band 7 (near IR) as a function |
---|
| 989 | ! of size for ice and water |
---|
| 990 | ! Fits from Steve Platnick |
---|
| 991 | ! |
---|
| 992 | |
---|
| 993 | integer, intent(in) :: phase |
---|
| 994 | real, intent(in) :: re |
---|
| 995 | real :: get_g_nir |
---|
[2720] | 996 | |
---|
| 997 | real, dimension(3), parameter :: ice_coefficients = (/ 0.7490, 6.5153e-3, -5.4136e-5 /), & |
---|
| 998 | small_water_coefficients = (/ 1.0364, -8.8800e-2, 7.0000e-3 /) |
---|
| 999 | real, dimension(4), parameter :: big_water_coefficients = (/ 0.6035, 2.8993e-2, -1.1051e-3, 1.5134e-5 /) |
---|
| 1000 | |
---|
| 1001 | ! approx. fits from MODIS Collection 6 LUT scattering calculations for 3.7 µm channel size retrievals |
---|
| 1002 | if(phase == phaseIsLiquid) then |
---|
| 1003 | if(re < 7.) then |
---|
| 1004 | get_g_nir = fit_to_quadratic(re, small_water_coefficients) |
---|
| 1005 | if(re < re_water_min) get_g_nir = fit_to_quadratic(re_water_min, small_water_coefficients) |
---|
| 1006 | else |
---|
| 1007 | get_g_nir = fit_to_cubic(re, big_water_coefficients) |
---|
| 1008 | if(re > re_water_max) get_g_nir = fit_to_cubic(re_water_max, big_water_coefficients) |
---|
| 1009 | end if |
---|
[2432] | 1010 | else |
---|
[2720] | 1011 | get_g_nir = fit_to_quadratic(re, ice_coefficients) |
---|
[2432] | 1012 | if(re < re_ice_min) get_g_nir = fit_to_quadratic(re_ice_min, ice_coefficients) |
---|
| 1013 | if(re > re_ice_max) get_g_nir = fit_to_quadratic(re_ice_max, ice_coefficients) |
---|
| 1014 | end if |
---|
| 1015 | |
---|
| 1016 | end function get_g_nir |
---|
| 1017 | |
---|
| 1018 | ! -------------------------------------------- |
---|
| 1019 | elemental function get_ssa_nir (phase, re) |
---|
| 1020 | integer, intent(in) :: phase |
---|
| 1021 | real, intent(in) :: re |
---|
| 1022 | real :: get_ssa_nir |
---|
| 1023 | ! |
---|
| 1024 | ! Polynomial fit for single scattering albedo in MODIS band 7 (near IR) as a function |
---|
| 1025 | ! of size for ice and water |
---|
| 1026 | ! Fits from Steve Platnick |
---|
| 1027 | ! |
---|
[2720] | 1028 | real, dimension(4), parameter :: ice_coefficients = (/ 0.9625, -1.8069e-2, 3.3281e-4,-2.2865e-6/) |
---|
| 1029 | real, dimension(3), parameter :: water_coefficients = (/ 1.0044, -1.1397e-2, 1.3300e-4 /) |
---|
[2432] | 1030 | |
---|
[2720] | 1031 | ! approx. fits from MODIS Collection 6 LUT scattering calculations |
---|
[2432] | 1032 | if(phase == phaseIsLiquid) then |
---|
| 1033 | get_ssa_nir = fit_to_quadratic(re, water_coefficients) |
---|
| 1034 | if(re < re_water_min) get_ssa_nir = fit_to_quadratic(re_water_min, water_coefficients) |
---|
| 1035 | if(re > re_water_max) get_ssa_nir = fit_to_quadratic(re_water_max, water_coefficients) |
---|
| 1036 | else |
---|
| 1037 | get_ssa_nir = fit_to_cubic(re, ice_coefficients) |
---|
| 1038 | if(re < re_ice_min) get_ssa_nir = fit_to_cubic(re_ice_min, ice_coefficients) |
---|
| 1039 | if(re > re_ice_max) get_ssa_nir = fit_to_cubic(re_ice_max, ice_coefficients) |
---|
| 1040 | end if |
---|
| 1041 | |
---|
| 1042 | end function get_ssa_nir |
---|
| 1043 | ! -------------------------------------------- |
---|
| 1044 | pure function fit_to_cubic(x, coefficients) |
---|
| 1045 | real, intent(in) :: x |
---|
| 1046 | real, dimension(:), intent(in) :: coefficients |
---|
| 1047 | real :: fit_to_cubic |
---|
| 1048 | |
---|
| 1049 | |
---|
| 1050 | fit_to_cubic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3) + x * coefficients(4))) |
---|
| 1051 | end function fit_to_cubic |
---|
| 1052 | ! -------------------------------------------- |
---|
| 1053 | pure function fit_to_quadratic(x, coefficients) |
---|
| 1054 | real, intent(in) :: x |
---|
| 1055 | real, dimension(:), intent(in) :: coefficients |
---|
| 1056 | real :: fit_to_quadratic |
---|
| 1057 | |
---|
| 1058 | |
---|
| 1059 | fit_to_quadratic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3))) |
---|
| 1060 | end function fit_to_quadratic |
---|
| 1061 | ! -------------------------------------------- |
---|
| 1062 | ! Radiative transfer |
---|
| 1063 | ! -------------------------------------------- |
---|
| 1064 | pure function compute_toa_reflectace(tau, g, w0) |
---|
| 1065 | real, dimension(:), intent(in) :: tau, g, w0 |
---|
| 1066 | real :: compute_toa_reflectace |
---|
| 1067 | |
---|
| 1068 | logical, dimension(size(tau)) :: cloudMask |
---|
| 1069 | integer, dimension(count(tau(:) > 0)) :: cloudIndicies |
---|
| 1070 | real, dimension(count(tau(:) > 0)) :: Refl, Trans |
---|
| 1071 | real :: Refl_tot, Trans_tot |
---|
| 1072 | integer :: i |
---|
| 1073 | ! --------------------------------------- |
---|
| 1074 | ! |
---|
| 1075 | ! This wrapper reports reflectance only and strips out non-cloudy elements from the calculation |
---|
| 1076 | ! |
---|
| 1077 | cloudMask = tau(:) > 0. |
---|
| 1078 | cloudIndicies = pack((/ (i, i = 1, size(tau)) /), mask = cloudMask) |
---|
| 1079 | do i = 1, size(cloudIndicies) |
---|
| 1080 | call two_stream(tau(cloudIndicies(i)), g(cloudIndicies(i)), w0(cloudIndicies(i)), Refl(i), Trans(i)) |
---|
| 1081 | end do |
---|
| 1082 | |
---|
| 1083 | call adding_doubling(Refl(:), Trans(:), Refl_tot, Trans_tot) |
---|
| 1084 | |
---|
| 1085 | compute_toa_reflectace = Refl_tot |
---|
| 1086 | |
---|
| 1087 | end function compute_toa_reflectace |
---|
| 1088 | ! -------------------------------------------- |
---|
| 1089 | pure subroutine two_stream(tauint, gint, w0int, ref, tra) |
---|
| 1090 | real, intent(in) :: tauint, gint, w0int |
---|
| 1091 | real, intent(out) :: ref, tra |
---|
| 1092 | ! |
---|
| 1093 | ! Compute reflectance in a single layer using the two stream approximation |
---|
| 1094 | ! The code itself is from Lazaros Oreopoulos via Steve Platnick |
---|
| 1095 | ! |
---|
| 1096 | ! ------------------------ |
---|
| 1097 | ! Local variables |
---|
| 1098 | ! for delta Eddington code |
---|
| 1099 | ! xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1) |
---|
| 1100 | integer, parameter :: beam = 2 |
---|
| 1101 | real, parameter :: xmu = 0.866, minConservativeW0 = 0.9999999 |
---|
| 1102 | real :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, & |
---|
| 1103 | rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den, th |
---|
| 1104 | ! |
---|
| 1105 | ! Compute reflectance and transmittance in a single layer using the two stream approximation |
---|
| 1106 | ! The code itself is from Lazaros Oreopoulos via Steve Platnick |
---|
| 1107 | ! |
---|
| 1108 | f = gint**2 |
---|
| 1109 | tau = (1 - w0int * f) * tauint |
---|
| 1110 | w0 = (1 - f) * w0int / (1 - w0int * f) |
---|
| 1111 | g = (gint - f) / (1 - f) |
---|
| 1112 | |
---|
| 1113 | ! delta-Eddington (Joseph et al. 1976) |
---|
| 1114 | gamma1 = (7 - w0* (4 + 3 * g)) / 4.0 |
---|
| 1115 | gamma2 = -(1 - w0* (4 - 3 * g)) / 4.0 |
---|
| 1116 | gamma3 = (2 - 3*g*xmu) / 4.0 |
---|
| 1117 | gamma4 = 1 - gamma3 |
---|
| 1118 | |
---|
| 1119 | if (w0int > minConservativeW0) then |
---|
| 1120 | ! Conservative scattering |
---|
| 1121 | if (beam == 1) then |
---|
| 1122 | rh = (gamma1*tau+(gamma3-gamma1*xmu)*(1-exp(-tau/xmu))) |
---|
| 1123 | ref = rh / (1 + gamma1 * tau) |
---|
| 1124 | tra = 1 - ref |
---|
| 1125 | else if(beam == 2) then |
---|
| 1126 | ref = gamma1*tau/(1 + gamma1*tau) |
---|
| 1127 | tra = 1 - ref |
---|
| 1128 | endif |
---|
| 1129 | else |
---|
| 1130 | ! Non-conservative scattering |
---|
| 1131 | a1 = gamma1 * gamma4 + gamma2 * gamma3 |
---|
| 1132 | a2 = gamma1 * gamma3 + gamma2 * gamma4 |
---|
| 1133 | |
---|
| 1134 | rk = sqrt(gamma1**2 - gamma2**2) |
---|
| 1135 | |
---|
| 1136 | r1 = (1 - rk * xmu) * (a2 + rk * gamma3) |
---|
| 1137 | r2 = (1 + rk * xmu) * (a2 - rk * gamma3) |
---|
| 1138 | r3 = 2 * rk *(gamma3 - a2 * xmu) |
---|
| 1139 | r4 = (1 - (rk * xmu)**2) * (rk + gamma1) |
---|
| 1140 | r5 = (1 - (rk * xmu)**2) * (rk - gamma1) |
---|
| 1141 | |
---|
| 1142 | t1 = (1 + rk * xmu) * (a1 + rk * gamma4) |
---|
| 1143 | t2 = (1 - rk * xmu) * (a1 - rk * gamma4) |
---|
| 1144 | t3 = 2 * rk * (gamma4 + a1 * xmu) |
---|
| 1145 | t4 = r4 |
---|
| 1146 | t5 = r5 |
---|
| 1147 | |
---|
| 1148 | beta = -r5 / r4 |
---|
| 1149 | |
---|
| 1150 | e1 = min(rk * tau, 500.) |
---|
| 1151 | e2 = min(tau / xmu, 500.) |
---|
| 1152 | |
---|
| 1153 | if (beam == 1) then |
---|
| 1154 | den = r4 * exp(e1) + r5 * exp(-e1) |
---|
| 1155 | ref = w0*(r1*exp(e1)-r2*exp(-e1)-r3*exp(-e2))/den |
---|
| 1156 | den = t4 * exp(e1) + t5 * exp(-e1) |
---|
| 1157 | th = exp(-e2) |
---|
| 1158 | tra = th-th*w0*(t1*exp(e1)-t2*exp(-e1)-t3*exp(e2))/den |
---|
| 1159 | elseif (beam == 2) then |
---|
| 1160 | ef1 = exp(-e1) |
---|
| 1161 | ef2 = exp(-2*e1) |
---|
| 1162 | ref = (gamma2*(1-ef2))/((rk+gamma1)*(1-beta*ef2)) |
---|
| 1163 | tra = (2*rk*ef1)/((rk+gamma1)*(1-beta*ef2)) |
---|
| 1164 | endif |
---|
| 1165 | end if |
---|
| 1166 | end subroutine two_stream |
---|
| 1167 | ! -------------------------------------------------- |
---|
| 1168 | elemental function two_stream_reflectance(tauint, gint, w0int) |
---|
| 1169 | real, intent(in) :: tauint, gint, w0int |
---|
| 1170 | real :: two_stream_reflectance |
---|
| 1171 | ! |
---|
| 1172 | ! Compute reflectance in a single layer using the two stream approximation |
---|
| 1173 | ! The code itself is from Lazaros Oreopoulos via Steve Platnick |
---|
| 1174 | ! |
---|
| 1175 | ! ------------------------ |
---|
| 1176 | ! Local variables |
---|
| 1177 | ! for delta Eddington code |
---|
| 1178 | ! xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1) |
---|
| 1179 | integer, parameter :: beam = 2 |
---|
| 1180 | real, parameter :: xmu = 0.866, minConservativeW0 = 0.9999999 |
---|
| 1181 | real :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, & |
---|
| 1182 | rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den |
---|
| 1183 | ! ------------------------ |
---|
| 1184 | |
---|
| 1185 | |
---|
| 1186 | f = gint**2 |
---|
| 1187 | tau = (1 - w0int * f) * tauint |
---|
| 1188 | w0 = (1 - f) * w0int / (1 - w0int * f) |
---|
| 1189 | g = (gint - f) / (1 - f) |
---|
| 1190 | |
---|
| 1191 | ! delta-Eddington (Joseph et al. 1976) |
---|
| 1192 | gamma1 = (7 - w0* (4 + 3 * g)) / 4.0 |
---|
| 1193 | gamma2 = -(1 - w0* (4 - 3 * g)) / 4.0 |
---|
| 1194 | gamma3 = (2 - 3*g*xmu) / 4.0 |
---|
| 1195 | gamma4 = 1 - gamma3 |
---|
| 1196 | |
---|
| 1197 | if (w0int > minConservativeW0) then |
---|
| 1198 | ! Conservative scattering |
---|
| 1199 | if (beam == 1) then |
---|
| 1200 | rh = (gamma1*tau+(gamma3-gamma1*xmu)*(1-exp(-tau/xmu))) |
---|
| 1201 | two_stream_reflectance = rh / (1 + gamma1 * tau) |
---|
| 1202 | elseif (beam == 2) then |
---|
| 1203 | two_stream_reflectance = gamma1*tau/(1 + gamma1*tau) |
---|
| 1204 | endif |
---|
| 1205 | |
---|
| 1206 | else ! |
---|
| 1207 | |
---|
| 1208 | ! Non-conservative scattering |
---|
| 1209 | a1 = gamma1 * gamma4 + gamma2 * gamma3 |
---|
| 1210 | a2 = gamma1 * gamma3 + gamma2 * gamma4 |
---|
| 1211 | |
---|
| 1212 | rk = sqrt(gamma1**2 - gamma2**2) |
---|
| 1213 | |
---|
| 1214 | r1 = (1 - rk * xmu) * (a2 + rk * gamma3) |
---|
| 1215 | r2 = (1 + rk * xmu) * (a2 - rk * gamma3) |
---|
| 1216 | r3 = 2 * rk *(gamma3 - a2 * xmu) |
---|
| 1217 | r4 = (1 - (rk * xmu)**2) * (rk + gamma1) |
---|
| 1218 | r5 = (1 - (rk * xmu)**2) * (rk - gamma1) |
---|
| 1219 | |
---|
| 1220 | t1 = (1 + rk * xmu) * (a1 + rk * gamma4) |
---|
| 1221 | t2 = (1 - rk * xmu) * (a1 - rk * gamma4) |
---|
| 1222 | t3 = 2 * rk * (gamma4 + a1 * xmu) |
---|
| 1223 | t4 = r4 |
---|
| 1224 | t5 = r5 |
---|
| 1225 | |
---|
| 1226 | beta = -r5 / r4 |
---|
| 1227 | |
---|
| 1228 | e1 = min(rk * tau, 500.) |
---|
| 1229 | e2 = min(tau / xmu, 500.) |
---|
| 1230 | |
---|
| 1231 | if (beam == 1) then |
---|
| 1232 | den = r4 * exp(e1) + r5 * exp(-e1) |
---|
| 1233 | two_stream_reflectance = w0*(r1*exp(e1)-r2*exp(-e1)-r3*exp(-e2))/den |
---|
| 1234 | elseif (beam == 2) then |
---|
| 1235 | ef1 = exp(-e1) |
---|
| 1236 | ef2 = exp(-2*e1) |
---|
| 1237 | two_stream_reflectance = (gamma2*(1-ef2))/((rk+gamma1)*(1-beta*ef2)) |
---|
| 1238 | endif |
---|
| 1239 | |
---|
| 1240 | end if |
---|
| 1241 | end function two_stream_reflectance |
---|
| 1242 | ! -------------------------------------------- |
---|
| 1243 | pure subroutine adding_doubling (Refl, Tran, Refl_tot, Tran_tot) |
---|
| 1244 | real, dimension(:), intent(in) :: Refl, Tran |
---|
| 1245 | real, intent(out) :: Refl_tot, Tran_tot |
---|
| 1246 | ! |
---|
| 1247 | ! Use adding/doubling formulas to compute total reflectance and transmittance from layer values |
---|
| 1248 | ! |
---|
| 1249 | |
---|
| 1250 | integer :: i |
---|
| 1251 | real, dimension(size(Refl)) :: Refl_cumulative, Tran_cumulative |
---|
| 1252 | |
---|
| 1253 | Refl_cumulative(1) = Refl(1); Tran_cumulative(1) = Tran(1) |
---|
| 1254 | |
---|
| 1255 | do i=2, size(Refl) |
---|
| 1256 | ! place (add) previous combined layer(s) reflectance on top of layer i, w/black surface (or ignoring surface): |
---|
| 1257 | Refl_cumulative(i) = Refl_cumulative(i-1) + Refl(i)*(Tran_cumulative(i-1)**2)/(1 - Refl_cumulative(i-1) * Refl(i)) |
---|
| 1258 | Tran_cumulative(i) = (Tran_cumulative(i-1)*Tran(i)) / (1 - Refl_cumulative(i-1) * Refl(i)) |
---|
| 1259 | end do |
---|
| 1260 | |
---|
| 1261 | Refl_tot = Refl_cumulative(size(Refl)) |
---|
| 1262 | Tran_tot = Tran_cumulative(size(Refl)) |
---|
| 1263 | |
---|
| 1264 | end subroutine adding_doubling |
---|
| 1265 | ! -------------------------------------------------- |
---|
| 1266 | subroutine complain_and_die(message) |
---|
| 1267 | character(len = *), intent(in) :: message |
---|
| 1268 | |
---|
| 1269 | write(6, *) "Failure in MODIS simulator" |
---|
| 1270 | write(6, *) trim(message) |
---|
| 1271 | stop |
---|
| 1272 | end subroutine complain_and_die |
---|
| 1273 | !------------------------------------------------------------------------------------------------ |
---|
[2720] | 1274 | end module mod_modis_sim |
---|