source: LMDZ6/trunk/libf/phylmd/cosp/mod_modis_sim.F90 @ 4514

Last change on this file since 4514 was 3286, checked in by idelkadi, 7 years ago
  • La cles ok_debug_cosp pour tourner en mode debug est mise a true par defaut
  • Correction dans le simulateur Modis
File size: 64.2 KB
Line 
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.
4! $Revision: 88 $, $Date: 2013-11-13 07:08:38 -0700 (Wed, 13 Nov 2013) $
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!
53module mod_modis_sim
54  USE MOD_COSP_TYPES, only: R_UNDEF, ok_debug_cosp
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
81! DJS2015: Remove unused parameter
82!  logical, parameter :: use_two_re_iterations = .false. ! do two retrieval iterations?
83! DJS2015 END 
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.
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
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
161contains
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    real :: seuil
387    ! ---------------------------------------------------
388   
389    if (ok_debug_cosp) then
390       seuil=1.e-9
391    else
392       seuil=0.0
393    endif
394
395    where(cloudIce(:, :) <= 0.)
396      tauLiquidFraction(:, :) = 1.
397    elsewhere
398      where (cloudWater(:, :) <= 0.)
399        tauLiquidFraction(:, :) = 0.
400      elsewhere
401        !
402        ! Geometic optics limit - tau as LWP/re  (proportional to LWC/re)
403        !
404! Modif AI 02 2018
405        tauLiquidFraction(:, :) = (cloudWater(:, :)/max(waterSize(:, :), seuil) ) / &
406                                  (cloudWater(:, :)/max(waterSize(:, :), seuil) + &
407                                   cloudIce(:, :)/(ice_density * max(iceSize(:, :), seuil)) )
408      end where
409    end where
410    liquid_opticalThickness(:, :) = tauLiquidFraction(:, :) * opticalThickness(:, :)
411    ice_opticalThickness   (:, :) = opticalThickness(:, :) - liquid_opticalThickness(:, :)
412   
413    call modis_L2_simulator_twoTaus(temp, pressureLayers, pressureLevels,          &
414                                    liquid_opticalThickness, ice_opticalThickness, &
415                                    waterSize, iceSize,                            &
416                                    isccpTau, isccpCloudTopPressure,               &
417                                    retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize)
418                               
419  end subroutine modis_L2_simulator_oneTau
420
421  ! ########################################################################################
422  subroutine modis_column(nPoints,nSubCols,phase, cloud_top_pressure, optical_thickness, particle_size,      &
423       Cloud_Fraction_Total_Mean,         Cloud_Fraction_Water_Mean,         Cloud_Fraction_Ice_Mean,        &
424       Cloud_Fraction_High_Mean,          Cloud_Fraction_Mid_Mean,           Cloud_Fraction_Low_Mean,        &
425       Optical_Thickness_Total_Mean,      Optical_Thickness_Water_Mean,      Optical_Thickness_Ice_Mean,     &
426       Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10,&
427       Cloud_Particle_Size_Water_Mean,    Cloud_Particle_Size_Ice_Mean,      Cloud_Top_Pressure_Total_Mean,  &
428       Liquid_Water_Path_Mean,            Ice_Water_Path_Mean,                                               &   
429       Optical_Thickness_vs_Cloud_Top_Pressure,Optical_Thickness_vs_ReffIce,Optical_Thickness_vs_ReffLiq)
430   
431    ! INPUTS
432    integer,intent(in) :: &
433         nPoints,                           & ! Number of horizontal gridpoints
434         nSubCols                             ! Number of subcolumns
435    integer,intent(in), dimension(:,:) ::  &
436!ds    integer,intent(in), dimension(nPoints, nSubCols) ::  &
437         phase                             
438    real,intent(in),dimension(:,:) ::  &
439!ds    real,intent(in),dimension(nPoints, nSubCols) ::  &
440         cloud_top_pressure,                &
441         optical_thickness,                 &
442         particle_size
443 
444    ! OUTPUTS
445    real,intent(inout),dimension(:)  ::   & !
446!ds    real,intent(inout),dimension(nPoints)  ::   & !
447         Cloud_Fraction_Total_Mean,         & !
448         Cloud_Fraction_Water_Mean,         & !
449         Cloud_Fraction_Ice_Mean,           & !
450         Cloud_Fraction_High_Mean,          & !
451         Cloud_Fraction_Mid_Mean,           & !
452         Cloud_Fraction_Low_Mean,           & !
453         Optical_Thickness_Total_Mean,      & !
454         Optical_Thickness_Water_Mean,      & !
455         Optical_Thickness_Ice_Mean,        & !
456         Optical_Thickness_Total_MeanLog10, & !
457         Optical_Thickness_Water_MeanLog10, & !
458         Optical_Thickness_Ice_MeanLog10,   & !
459         Cloud_Particle_Size_Water_Mean,    & !
460         Cloud_Particle_Size_Ice_Mean,      & !
461         Cloud_Top_Pressure_Total_Mean,     & !
462         Liquid_Water_Path_Mean,            & !
463         Ice_Water_Path_Mean                  !
464    real,intent(inout),dimension(:,:,:) :: &
465!ds    real,intent(inout),dimension(nPoints,numTauHistogramBins,numPressureHistogramBins) :: &
466         Optical_Thickness_vs_Cloud_Top_Pressure
467    real,intent(inout),dimension(:,:,:) :: &   
468!ds    real,intent(inout),dimension(nPoints,numTauHistogramBins,numMODISReffIceBins) :: &   
469         Optical_Thickness_vs_ReffIce
470    real,intent(inout),dimension(:,:,:) :: &   
471!ds    real,intent(inout),dimension(nPoints,numTauHistogramBins,numMODISReffLiqBins) :: &   
472         Optical_Thickness_vs_ReffLiq         
473
474    ! LOCAL VARIABLES
475    real, parameter :: &
476         LWP_conversion = 2./3. * 1000. ! MKS units 
477    integer :: i, j
478    logical, dimension(nPoints,nSubCols) :: &
479         cloudMask,      &
480         waterCloudMask, &
481         iceCloudMask,   &
482         validRetrievalMask
483    real,dimension(nPoints,nSubCols) :: &
484         tauWRK,ctpWRK,reffIceWRK,reffLiqWRK
485   
486    ! ########################################################################################
487    ! Include only those pixels with successful retrievals in the statistics
488    ! ########################################################################################
489    validRetrievalMask(1:nPoints,1:nSubCols) = particle_size(1:nPoints,1:nSubCols) > 0.
490    cloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) /= phaseIsNone .and.       &
491         validRetrievalMask(1:nPoints,1:nSubCols)
492    waterCloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) == phaseIsLiquid .and. &
493         validRetrievalMask(1:nPoints,1:nSubCols)
494    iceCloudMask(1:nPoints,1:nSubCols)   = phase(1:nPoints,1:nSubCols) == phaseIsIce .and.    &
495         validRetrievalMask(1:nPoints,1:nSubCols)
496
497    ! ########################################################################################
498    ! Use these as pixel counts at first
499    ! ########################################################################################
500    Cloud_Fraction_Total_Mean(1:nPoints) = real(count(cloudMask,      dim = 2))
501    Cloud_Fraction_Water_Mean(1:nPoints) = real(count(waterCloudMask, dim = 2))
502    Cloud_Fraction_Ice_Mean(1:nPoints)   = real(count(iceCloudMask,   dim = 2))
503    Cloud_Fraction_High_Mean(1:nPoints)  = real(count(cloudMask .and. cloud_top_pressure <=          &
504                                           highCloudPressureLimit, dim = 2))
505    Cloud_Fraction_Low_Mean(1:nPoints)   = real(count(cloudMask .and. cloud_top_pressure >           &
506                                           lowCloudPressureLimit,  dim = 2))
507    Cloud_Fraction_Mid_Mean(1:nPoints)   = Cloud_Fraction_Total_Mean(1:nPoints) - Cloud_Fraction_High_Mean(1:nPoints)&
508                                           - Cloud_Fraction_Low_Mean(1:nPoints)
509
510    ! ########################################################################################
511    ! Compute mean optical thickness.
512    ! ########################################################################################
513    where (Cloud_Fraction_Total_Mean == 0)
514       Optical_Thickness_Total_Mean      = R_UNDEF
515    elsewhere
516       Optical_Thickness_Total_Mean(1:nPoints) = sum(optical_thickness, mask = cloudMask,      dim = 2) / &
517                                                 Cloud_Fraction_Total_Mean(1:nPoints)
518    endwhere
519    where (Cloud_Fraction_Water_Mean == 0)
520       Optical_Thickness_Water_Mean      = R_UNDEF
521    elsewhere
522       Optical_Thickness_Water_Mean(1:nPoints) = sum(optical_thickness, mask = waterCloudMask, dim = 2) / &
523                                                 Cloud_Fraction_Water_Mean(1:nPoints)
524    endwhere
525    where (Cloud_Fraction_Ice_Mean == 0)
526       Optical_Thickness_Ice_Mean        = R_UNDEF
527    elsewhere
528       Optical_Thickness_Ice_Mean(1:nPoints)   = sum(optical_thickness, mask = iceCloudMask,   dim = 2) / &
529                                                 Cloud_Fraction_Ice_Mean(1:nPoints)
530    endwhere
531   
532    ! ########################################################################################
533    ! We take the absolute value of optical thickness here to satisfy compilers that complains
534    ! when we evaluate the logarithm of a negative number, even though it's not included in
535    ! the sum.
536    ! ########################################################################################
537    where (Cloud_Fraction_Total_Mean == 0)
538       Optical_Thickness_Total_MeanLog10 = R_UNDEF
539    elsewhere
540       Optical_Thickness_Total_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = cloudMask, &
541           dim = 2) / Cloud_Fraction_Total_Mean(1:nPoints)
542    endwhere
543    where (Cloud_Fraction_Water_Mean == 0)
544       Optical_Thickness_Water_MeanLog10 = R_UNDEF
545    elsewhere
546       Optical_Thickness_Water_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = waterCloudMask,&
547            dim = 2) / Cloud_Fraction_Water_Mean(1:nPoints)
548    endwhere
549    where (Cloud_Fraction_Ice_Mean == 0)
550       Optical_Thickness_Ice_MeanLog10   = R_UNDEF
551    elsewhere
552       Optical_Thickness_Ice_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = iceCloudMask,&
553            dim = 2) / Cloud_Fraction_Ice_Mean(1:nPoints)
554    endwhere
555    where (Cloud_Fraction_Water_Mean == 0)
556       Cloud_Particle_Size_Water_Mean    = R_UNDEF
557    elsewhere
558       Cloud_Particle_Size_Water_Mean(1:nPoints) = sum(particle_size, mask = waterCloudMask, dim = 2) / &
559            Cloud_Fraction_Water_Mean(1:nPoints)
560    endwhere
561    where (Cloud_Fraction_Ice_Mean == 0)
562       Cloud_Particle_Size_Ice_Mean      = R_UNDEF 
563    elsewhere
564       Cloud_Particle_Size_Ice_Mean(1:nPoints) = sum(particle_size, mask = iceCloudMask,   dim = 2) / &
565            Cloud_Fraction_Ice_Mean(1:nPoints)
566    endwhere
567    where (Cloud_Fraction_Total_Mean == 0)
568       Cloud_Top_Pressure_Total_Mean     = R_UNDEF
569    elsewhere
570       Cloud_Top_Pressure_Total_Mean(1:nPoints) = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / &
571            max(1, count(cloudMask, dim = 2))
572    endwhere
573    where (Cloud_Fraction_Water_Mean == 0)
574       Liquid_Water_Path_Mean            = R_UNDEF
575    elsewhere
576       Liquid_Water_Path_Mean(1:nPoints) = LWP_conversion*sum(particle_size*optical_thickness, &
577            mask=waterCloudMask,dim=2)/Cloud_Fraction_Water_Mean(1:nPoints)
578    endwhere
579    where (Cloud_Fraction_Ice_Mean == 0)
580       Ice_Water_Path_Mean               = R_UNDEF
581    elsewhere
582       Ice_Water_Path_Mean(1:nPoints) = LWP_conversion * ice_density*sum(particle_size*optical_thickness,&
583             mask=iceCloudMask,dim = 2) /Cloud_Fraction_Ice_Mean(1:nPoints)
584    endwhere
585    ! ########################################################################################
586    ! Normalize pixel counts to fraction.
587    ! ########################################################################################
588    Cloud_Fraction_High_Mean(1:nPoints)  = Cloud_Fraction_High_Mean(1:nPoints)  /nSubcols
589    Cloud_Fraction_Mid_Mean(1:nPoints)   = Cloud_Fraction_Mid_Mean(1:nPoints)   /nSubcols
590    Cloud_Fraction_Low_Mean(1:nPoints)   = Cloud_Fraction_Low_Mean(1:nPoints)   /nSubcols
591    Cloud_Fraction_Total_Mean(1:nPoints) = Cloud_Fraction_Total_Mean(1:nPoints) /nSubcols
592    Cloud_Fraction_Ice_Mean(1:nPoints)   = Cloud_Fraction_Ice_Mean(1:nPoints)   /nSubcols
593    Cloud_Fraction_Water_Mean(1:nPoints) = Cloud_Fraction_Water_Mean(1:nPoints) /nSubcols
594   
595    ! ########################################################################################
596    ! Set clear-scenes to undefined
597    ! ########################################################################################
598!    where (Cloud_Fraction_High_Mean == 0)  Cloud_Fraction_High_Mean = R_UNDEF
599!    where (Cloud_Fraction_Mid_Mean == 0)   Cloud_Fraction_Mid_Mean = R_UNDEF
600!    where (Cloud_Fraction_Low_Mean == 0)   Cloud_Fraction_Low_Mean = R_UNDEF
601
602    ! ########################################################################################
603    ! Joint histogram 
604    ! ########################################################################################
605
606    ! Loop over all points
607    tauWRK(1:nPoints,1:nSubCols)     = optical_thickness(1:nPoints,1:nSubCols)
608    ctpWRK(1:nPoints,1:nSubCols)     = cloud_top_pressure(1:nPoints,1:nSubCols)
609    reffIceWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,iceCloudMask)
610    reffLiqWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,waterCloudMask)
611    do j=1,nPoints
612
613       ! Fill clear and optically thin subcolumns with fill
614       where(.not. cloudMask(j,1:nSubCols))
615          tauWRK(j,1:nSubCols) = -999.
616          ctpWRK(j,1:nSubCols) = -999.
617       endwhere
618       ! Joint histogram of tau/CTP
619       call hist2D(tauWRK(j,1:nSubCols),ctpWRK(j,1:nSubCols),nSubCols,&
620                   tauHistogramBoundaries,numTauHistogramBins,&
621                   pressureHistogramBoundaries,numPressureHistogramBins,&
622                   Optical_Thickness_vs_Cloud_Top_Pressure(j,1:numTauHistogramBins,1:numPressureHistogramBins))
623       ! Joint histogram of tau/ReffICE
624       call hist2D(tauWRK(j,1:nSubCols),reffIceWrk(j,1:nSubCols),nSubCols,               &
625                   tauHistogramBoundaries,numTauHistogramBins,reffICE_binBounds,         &
626                   numMODISReffIceBins, Optical_Thickness_vs_ReffIce(j,1:numTauHistogramBins,1:numMODISReffIceBins))
627       ! Joint histogram of tau/ReffLIQ
628       call hist2D(tauWRK(j,1:nSubCols),reffLiqWrk(j,1:nSubCols),nSubCols,               &
629                   tauHistogramBoundaries,numTauHistogramBins,reffLIQ_binBounds,         &
630                   numMODISReffLiqBins, Optical_Thickness_vs_ReffLiq(j,1:numTauHistogramBins,1:numMODISReffLiqBins))                   
631
632    enddo   
633    Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numTauHistogramBins,1:numPressureHistogramBins) = &
634         Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numTauHistogramBins,1:numPressureHistogramBins)/nSubCols
635    Optical_Thickness_vs_ReffIce(1:nPoints,1:numTauHistogramBins,1:numMODISReffIceBins) = &
636         Optical_Thickness_vs_ReffIce(1:nPoints,1:numTauHistogramBins,1:numMODISReffIceBins)/nSubCols
637    Optical_Thickness_vs_ReffLiq(1:nPoints,1:numTauHistogramBins,1:numMODISReffLiqBins) = &
638         Optical_Thickness_vs_ReffLiq(1:nPoints,1:numTauHistogramBins,1:numMODISReffLiqBins)/nSubCols
639
640  end subroutine modis_column
641  ! ######################################################################################
642  ! SUBROUTINE hist2D
643  ! ######################################################################################
644  subroutine hist2D(var1,var2,npts,bin1,nbin1,bin2,nbin2,jointHist)
645    implicit none
646   
647    ! INPUTS
648    integer, intent(in) :: &
649         npts,  & ! Number of data points to be sorted
650         nbin1, & ! Number of bins in histogram direction 1
651         nbin2    ! Number of bins in histogram direction 2
652    real,intent(in),dimension(npts) :: &
653         var1,  & ! Variable 1 to be sorted into bins
654         var2     ! variable 2 to be sorted into bins
655    real,intent(in),dimension(nbin1+1) :: &
656         bin1     ! Histogram bin 1 boundaries
657    real,intent(in),dimension(nbin2+1) :: &
658         bin2     ! Histogram bin 2 boundaries
659    ! OUTPUTS
660    real,intent(out),dimension(nbin1,nbin2) :: &
661         jointHist
662   
663    ! LOCAL VARIABLES
664    integer :: ij,ik
665   
666    do ij=2,nbin1+1
667       do ik=2,nbin2+1
668          jointHist(ij-1,ik-1)=count(var1 .ge. bin1(ij-1) .and. var1 .lt. bin1(ij) .and. &
669               var2 .ge. bin2(ik-1) .and. var2 .lt. bin2(ik))       
670       enddo
671    enddo
672  end subroutine hist2D
673 
674  !------------------------------------------------------------------------------------------------
675  subroutine modis_L3_simulator(phase, cloud_top_pressure, optical_thickness, particle_size,            &
676       Cloud_Fraction_Total_Mean,       Cloud_Fraction_Water_Mean,       Cloud_Fraction_Ice_Mean,       &
677       Cloud_Fraction_High_Mean,        Cloud_Fraction_Mid_Mean,         Cloud_Fraction_Low_Mean,       &
678       Optical_Thickness_Total_Mean,    Optical_Thickness_Water_Mean,    Optical_Thickness_Ice_Mean,    &
679       Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10, &
680                                        Cloud_Particle_Size_Water_Mean,  Cloud_Particle_Size_Ice_Mean,  &
681       Cloud_Top_Pressure_Total_Mean,                                                                   &
682                                        Liquid_Water_Path_Mean,          Ice_Water_Path_Mean,           &   
683       Optical_Thickness_vs_Cloud_Top_Pressure)
684    !
685    ! Inputs; dimension nPoints, nSubcols
686    !
687    integer, dimension(:, :),   intent(in)  :: phase
688    real,    dimension(:, :),   intent(in)  :: cloud_top_pressure, optical_thickness, particle_size
689    !
690    ! Outputs; dimension nPoints
691    !
692    real,    dimension(:),      intent(out) :: &
693       Cloud_Fraction_Total_Mean,       Cloud_Fraction_Water_Mean,       Cloud_Fraction_Ice_Mean,       &
694       Cloud_Fraction_High_Mean,        Cloud_Fraction_Mid_Mean,         Cloud_Fraction_Low_Mean,       &
695       Optical_Thickness_Total_Mean,    Optical_Thickness_Water_Mean,    Optical_Thickness_Ice_Mean,    &
696       Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10, &
697                                        Cloud_Particle_Size_Water_Mean,  Cloud_Particle_Size_Ice_Mean,  &
698       Cloud_Top_Pressure_Total_Mean,                                                                   &
699                                        Liquid_Water_Path_Mean,          Ice_Water_Path_Mean
700    ! tau/ctp histogram; dimensions nPoints, numTauHistogramBins , numPressureHistogramBins
701    real,    dimension(:, :, :), intent(out) :: Optical_Thickness_vs_Cloud_Top_Pressure
702    ! ---------------------------
703    ! Local variables
704    !
705    real, parameter :: LWP_conversion = 2./3. * 1000. ! MKS units 
706    integer :: i, j
707    integer :: nPoints, nSubcols
708    logical, dimension(size(phase, 1), size(phase, 2)) :: &
709      cloudMask, waterCloudMask, iceCloudMask, validRetrievalMask
710    logical, dimension(size(phase, 1), size(phase, 2), numTauHistogramBins     ) :: tauMask
711    logical, dimension(size(phase, 1), size(phase, 2), numPressureHistogramBins) :: pressureMask
712    ! ---------------------------
713   
714    nPoints  = size(phase, 1)
715    nSubcols = size(phase, 2)
716    !
717    ! Array conformance checks
718    !
719    if(any( (/ size(cloud_top_pressure, 1), size(optical_thickness, 1), size(particle_size, 1),                                &
720               size(Cloud_Fraction_Total_Mean),       size(Cloud_Fraction_Water_Mean),       size(Cloud_Fraction_Ice_Mean),    &
721               size(Cloud_Fraction_High_Mean),        size(Cloud_Fraction_Mid_Mean),         size(Cloud_Fraction_Low_Mean),    &
722               size(Optical_Thickness_Total_Mean),    size(Optical_Thickness_Water_Mean),    size(Optical_Thickness_Ice_Mean), &
723               size(Optical_Thickness_Total_MeanLog10), size(Optical_Thickness_Water_MeanLog10), &
724               size(Optical_Thickness_Ice_MeanLog10),   size(Cloud_Particle_Size_Water_Mean),    &
725               size(Cloud_Particle_Size_Ice_Mean),      size(Cloud_Top_Pressure_Total_Mean),     &
726               size(Liquid_Water_Path_Mean),          size(Ice_Water_Path_Mean) /) /= nPoints))  &
727      call complain_and_die("Some L3 arrays have wrong number of grid points")
728    if(any( (/ size(cloud_top_pressure, 2), size(optical_thickness, 2), size(particle_size, 2) /)  /= nSubcols)) &
729      call complain_and_die("Some L3 arrays have wrong number of subcolumns")
730   
731   
732    !
733    ! Include only those pixels with successful retrievals in the statistics
734    !
735    validRetrievalMask(:, :) = particle_size(:, :) > 0.
736    cloudMask      = phase(:, :) /= phaseIsNone   .and. validRetrievalMask(:, :)
737    waterCloudMask = phase(:, :) == phaseIsLiquid .and. validRetrievalMask(:, :)
738    iceCloudMask   = phase(:, :) == phaseIsIce    .and. validRetrievalMask(:, :)
739    !
740    ! Use these as pixel counts at first
741    !
742    Cloud_Fraction_Total_Mean(:) = real(count(cloudMask,      dim = 2))
743    Cloud_Fraction_Water_Mean(:) = real(count(waterCloudMask, dim = 2))
744    Cloud_Fraction_Ice_Mean(:)   = real(count(iceCloudMask,   dim = 2))
745   
746    Cloud_Fraction_High_Mean(:) = real(count(cloudMask .and. cloud_top_pressure <= highCloudPressureLimit, dim = 2))
747    Cloud_Fraction_Low_Mean(:)  = real(count(cloudMask .and. cloud_top_pressure >  lowCloudPressureLimit,  dim = 2))
748    Cloud_Fraction_Mid_Mean(:)  = Cloud_Fraction_Total_Mean(:) - Cloud_Fraction_High_Mean(:) - Cloud_Fraction_Low_Mean(:)
749   
750    !
751    ! Don't want to divide by 0, even though the sums will be 0 where the pixel counts are 0.
752    !
753    where (Cloud_Fraction_Total_Mean == 0) Cloud_Fraction_Total_Mean = -1.
754    where (Cloud_Fraction_Water_Mean == 0) Cloud_Fraction_Water_Mean = -1.
755    where (Cloud_Fraction_Ice_Mean   == 0) Cloud_Fraction_Ice_Mean   = -1.
756   
757    Optical_Thickness_Total_Mean = sum(optical_thickness, mask = cloudMask,      dim = 2) / Cloud_Fraction_Total_Mean(:)
758    Optical_Thickness_Water_Mean = sum(optical_thickness, mask = waterCloudMask, dim = 2) / Cloud_Fraction_Water_Mean(:)
759    Optical_Thickness_Ice_Mean   = sum(optical_thickness, mask = iceCloudMask,   dim = 2) / Cloud_Fraction_Ice_Mean(:)
760   
761    ! We take the absolute value of optical thickness here to satisfy compilers that complains when we
762    !   evaluate the logarithm of a negative number, even though it's not included in the sum.
763    Optical_Thickness_Total_MeanLog10 = sum(log10(abs(optical_thickness)), mask = cloudMask,      dim = 2) / &
764                                        Cloud_Fraction_Total_Mean(:)
765    Optical_Thickness_Water_MeanLog10 = sum(log10(abs(optical_thickness)), mask = waterCloudMask, dim = 2) / &
766                                        Cloud_Fraction_Water_Mean(:)
767    Optical_Thickness_Ice_MeanLog10   = sum(log10(abs(optical_thickness)), mask = iceCloudMask,   dim = 2) / &
768                                        Cloud_Fraction_Ice_Mean(:)
769   
770    Cloud_Particle_Size_Water_Mean = sum(particle_size, mask = waterCloudMask, dim = 2) / Cloud_Fraction_Water_Mean(:)
771    Cloud_Particle_Size_Ice_Mean   = sum(particle_size, mask = iceCloudMask,   dim = 2) / Cloud_Fraction_Ice_Mean(:)
772   
773    Cloud_Top_Pressure_Total_Mean = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / max(1, count(cloudMask, dim = 2))
774   
775    Liquid_Water_Path_Mean = LWP_conversion &
776                             * sum(particle_size * optical_thickness, mask = waterCloudMask, dim = 2) &
777                             / Cloud_Fraction_Water_Mean(:)
778    Ice_Water_Path_Mean    = LWP_conversion * ice_density &
779                             * sum(particle_size * optical_thickness, mask = iceCloudMask,   dim = 2) &
780                             / Cloud_Fraction_Ice_Mean(:)
781
782    !
783    ! Normalize pixel counts to fraction
784    !   The first three cloud fractions have been set to -1 in cloud-free areas, so set those places to 0.
785    !
786    Cloud_Fraction_Total_Mean(:) = max(0., Cloud_Fraction_Total_Mean(:)/nSubcols)
787    Cloud_Fraction_Water_Mean(:) = max(0., Cloud_Fraction_Water_Mean(:)/nSubcols)
788    Cloud_Fraction_Ice_Mean(:)   = max(0., Cloud_Fraction_Ice_Mean(:)  /nSubcols)
789   
790    Cloud_Fraction_High_Mean(:)  = Cloud_Fraction_High_Mean(:) /nSubcols
791    Cloud_Fraction_Mid_Mean(:)   = Cloud_Fraction_Mid_Mean(:)  /nSubcols
792    Cloud_Fraction_Low_Mean(:)   = Cloud_Fraction_Low_Mean(:)  /nSubcols
793   
794    ! ----
795    ! Joint histogram
796    !
797    do i = 1, numTauHistogramBins
798      where(cloudMask(:, :))
799        tauMask(:, :, i) = optical_thickness(:, :) >= tauHistogramBoundaries(i) .and. &
800                           optical_thickness(:, :) <  tauHistogramBoundaries(i+1)
801      elsewhere
802        tauMask(:, :, i) = .false.
803      end where
804    end do
805
806    do i = 1, numPressureHistogramBins
807      where(cloudMask(:, :))
808        pressureMask(:, :, i) = cloud_top_pressure(:, :) >= pressureHistogramBoundaries(i) .and. &
809                                cloud_top_pressure(:, :) <  pressureHistogramBoundaries(i+1)
810      elsewhere
811        pressureMask(:, :, i) = .false.
812      end where
813    end do
814   
815    do i = 1, numPressureHistogramBins
816      do j = 1, numTauHistogramBins
817        Optical_Thickness_vs_Cloud_Top_Pressure(:, j, i) = &
818          real(count(tauMask(:, :, j) .and. pressureMask(:, :, i), dim = 2)) / real(nSubcols)
819      end do
820    end do
821   
822  end subroutine modis_L3_simulator
823  !------------------------------------------------------------------------------------------------
824  function cloud_top_pressure(tauIncrement, pressure, tauLimit)
825    real, dimension(:), intent(in) :: tauIncrement, pressure
826    real,               intent(in) :: tauLimit
827    real                           :: cloud_top_pressure
828    !
829    ! Find the extinction-weighted pressure. Assume that pressure varies linearly between
830    !   layers and use the trapezoidal rule.
831    !
832   
833    real :: deltaX, totalTau, totalProduct
834    integer :: i
835   
836    totalTau = 0.; totalProduct = 0.
837    do i = 2, size(tauIncrement)
838      if(totalTau + tauIncrement(i) > tauLimit) then
839        deltaX = tauLimit - totalTau
840        totalTau = totalTau + deltaX
841        !
842        ! Result for trapezoidal rule when you take less than a full step
843        !   tauIncrement is a layer-integrated value
844        !
845        totalProduct = totalProduct           &
846                     + pressure(i-1) * deltaX &
847                     + (pressure(i) - pressure(i-1)) * deltaX**2/(2. * tauIncrement(i))
848      else
849        totalTau =     totalTau     + tauIncrement(i)
850        totalProduct = totalProduct + tauIncrement(i) * (pressure(i) + pressure(i-1)) / 2.
851      end if
852      if(totalTau >= tauLimit) exit
853    end do
854    cloud_top_pressure = totalProduct/totalTau
855  end function cloud_top_pressure
856  !------------------------------------------------------------------------------------------------
857  function weight_by_extinction(tauIncrement, f, tauLimit)
858    real, dimension(:), intent(in) :: tauIncrement, f
859    real,               intent(in) :: tauLimit
860    real                           :: weight_by_extinction
861    !
862    ! Find the extinction-weighted value of f(tau), assuming constant f within each layer
863    !
864   
865    real    :: deltaX, totalTau, totalProduct
866    integer :: i
867   
868    totalTau = 0.; totalProduct = 0.
869    do i = 1, size(tauIncrement)
870      if(totalTau + tauIncrement(i) > tauLimit) then
871        deltaX       = tauLimit - totalTau
872        totalTau     = totalTau     + deltaX
873        totalProduct = totalProduct + deltaX * f(i)
874      else
875        totalTau     = totalTau     + tauIncrement(i)
876        totalProduct = totalProduct + tauIncrement(i) * f(i)
877      end if
878      if(totalTau >= tauLimit) exit
879    end do
880    weight_by_extinction = totalProduct/totalTau
881  end function weight_by_extinction
882  !------------------------------------------------------------------------------------------------
883  pure function compute_nir_reflectance(water_tau, water_size, ice_tau, ice_size)
884    real, dimension(:), intent(in) :: water_tau, water_size, ice_tau, ice_size
885    real                           :: compute_nir_reflectance
886   
887    real, dimension(size(water_tau)) :: water_g, water_w0, ice_g, ice_w0, &
888                                        tau, g, w0
889    !----------------------------------------
890    water_g(:)  = get_g_nir(  phaseIsLiquid, water_size)
891    water_w0(:) = get_ssa_nir(phaseIsLiquid, water_size)
892    ice_g(:)    = get_g_nir(  phaseIsIce,    ice_size)
893    ice_w0(:)   = get_ssa_nir(phaseIsIce,    ice_size)
894    !
895    ! Combine ice and water optical properties
896    !
897    g(:) = 0; w0(:) = 0.
898    tau(:) = ice_tau(:) + water_tau(:)
899    where (tau(:) > 0)
900      w0(:) = (water_tau(:) * water_w0(:)  + ice_tau(:) * ice_w0(:)) / &
901              tau(:)
902      g(:) = (water_tau(:) * water_g(:) * water_w0(:)  + ice_tau(:) * ice_g(:) * ice_w0(:)) / &
903             (w0(:) * tau(:))
904    end where
905   
906    compute_nir_reflectance = compute_toa_reflectace(tau, g, w0)
907  end function compute_nir_reflectance
908  !------------------------------------------------------------------------------------------------
909  ! Retreivals
910  !------------------------------------------------------------------------------------------------
911  elemental function retrieve_re (phase, tau, obs_Refl_nir)
912      integer, intent(in) :: phase
913      real,    intent(in) :: tau, obs_Refl_nir
914      real                :: retrieve_re
915      !
916      ! Finds the re that produces the minimum mis-match between predicted and observed reflectance in
917      !   MODIS band 7 (near IR)
918      ! Uses
919      !  fits for asymmetry parameter g(re) and single scattering albedo w0(re) based on MODIS tables
920      !  two-stream for layer reflectance and transmittance as a function of optical thickness tau, g, and w0
921      !  adding-doubling for total reflectance
922      ! 
923      !
924      !
925      ! Local variables
926      !
927      real, parameter :: min_distance_to_boundary = 0.01
928      real    :: re_min, re_max, delta_re
929      integer :: i
930     
931      real, dimension(num_trial_res) :: trial_re, g, w0, predicted_Refl_nir
932      ! --------------------------
933   
934    if(any(phase == (/ phaseIsLiquid, phaseIsUndetermined, phaseIsIce /))) then
935      if (phase == phaseIsLiquid .OR. phase == phaseIsUndetermined) then
936        re_min = re_water_min
937        re_max = re_water_max
938        trial_re(:) = trial_re_w
939        g(:)   = g_w(:)
940        w0(:)  = w0_w(:)
941      else
942        re_min = re_ice_min
943        re_max = re_ice_max
944        trial_re(:) = trial_re_i
945        g(:)   = g_i(:)
946        w0(:)  = w0_i(:)
947      end if
948      !
949      ! 1st attempt at index: w/coarse re resolution
950      !
951      predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:))
952      retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir)
953      !
954      ! If first retrieval works, can try 2nd iteration using greater re resolution
955      !
956! DJS2015: Remove unused piece of code     
957!      if(use_two_re_iterations .and. retrieve_re > 0.) then
958!        re_min = retrieve_re - delta_re
959!        re_max = retrieve_re + delta_re
960!        delta_re = (re_max - re_min)/real(num_trial_res-1)
961
962!        trial_re(:) = re_min + delta_re * (/ (i - 1, i = 1, num_trial_res) /)
963!        g(:)  = get_g_nir(  phase, trial_re(:))
964!        w0(:) = get_ssa_nir(phase, trial_re(:))
965!        predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:))
966!        retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir)
967!      end if
968! DJS2015 END
969    else
970      retrieve_re = re_fill
971    end if
972   
973  end function retrieve_re
974  ! --------------------------------------------
975  pure function interpolate_to_min(x, y, yobs)
976    real, dimension(:), intent(in) :: x, y
977    real,               intent(in) :: yobs
978    real                           :: interpolate_to_min
979    !
980    ! Given a set of values of y as y(x), find the value of x that minimizes abs(y - yobs)
981    !   y must be monotonic in x
982    !
983    real, dimension(size(x)) :: diff
984    integer                  :: nPoints, minDiffLoc, lowerBound, upperBound
985    ! ---------------------------------
986    nPoints = size(y)
987    diff(:) = y(:) - yobs
988    minDiffLoc = minloc(abs(diff), dim = 1)
989   
990    if(minDiffLoc == 1) then
991      lowerBound = minDiffLoc
992      upperBound = minDiffLoc + 1
993    else if(minDiffLoc == nPoints) then
994      lowerBound = minDiffLoc - 1
995      upperBound = minDiffLoc
996    else
997      if(diff(minDiffLoc-1) * diff(minDiffLoc) < 0) then
998        lowerBound = minDiffLoc-1
999        upperBound = minDiffLoc
1000      else
1001        lowerBound = minDiffLoc
1002        upperBound = minDiffLoc + 1
1003      end if
1004    end if
1005   
1006    if(diff(lowerBound) * diff(upperBound) < 0) then     
1007      !
1008      ! Interpolate the root position linearly if we bracket the root
1009      !
1010      interpolate_to_min = x(upperBound) - &
1011                           diff(upperBound) * (x(upperBound) - x(lowerBound)) / (diff(upperBound) - diff(lowerBound))
1012    else
1013      interpolate_to_min = re_fill
1014    end if
1015   
1016
1017  end function interpolate_to_min
1018  ! --------------------------------------------
1019  ! Optical properties
1020  ! --------------------------------------------
1021  elemental function get_g_nir (phase, re)
1022    !
1023    ! Polynomial fit for asummetry parameter g in MODIS band 7 (near IR) as a function
1024    !   of size for ice and water
1025    ! Fits from Steve Platnick
1026    !
1027
1028    integer, intent(in) :: phase
1029    real,    intent(in) :: re
1030    real :: get_g_nir
1031
1032    real, dimension(3), parameter :: ice_coefficients         = (/ 0.7490, 6.5153e-3, -5.4136e-5 /), &
1033                                     small_water_coefficients = (/ 1.0364, -8.8800e-2, 7.0000e-3 /)
1034    real, dimension(4), parameter :: big_water_coefficients   = (/ 0.6035, 2.8993e-2, -1.1051e-3, 1.5134e-5 /)
1035
1036    ! approx. fits from MODIS Collection 6 LUT scattering calculations for 3.7 µm channel size retrievals
1037    if(phase == phaseIsLiquid) then
1038       if(re < 7.) then
1039          get_g_nir = fit_to_quadratic(re, small_water_coefficients)
1040          if(re < re_water_min) get_g_nir = fit_to_quadratic(re_water_min, small_water_coefficients)
1041       else
1042          get_g_nir = fit_to_cubic(re, big_water_coefficients)
1043          if(re > re_water_max) get_g_nir = fit_to_cubic(re_water_max, big_water_coefficients)
1044       end if
1045    else
1046       get_g_nir = fit_to_quadratic(re, ice_coefficients)
1047      if(re < re_ice_min) get_g_nir = fit_to_quadratic(re_ice_min, ice_coefficients)
1048      if(re > re_ice_max) get_g_nir = fit_to_quadratic(re_ice_max, ice_coefficients)
1049    end if
1050   
1051  end function get_g_nir
1052
1053  ! --------------------------------------------
1054    elemental function get_ssa_nir (phase, re)
1055        integer, intent(in) :: phase
1056        real,    intent(in) :: re
1057        real                :: get_ssa_nir
1058        !
1059        ! Polynomial fit for single scattering albedo in MODIS band 7 (near IR) as a function
1060        !   of size for ice and water
1061        ! Fits from Steve Platnick
1062        !
1063        real, dimension(4), parameter :: ice_coefficients   = (/ 0.9625, -1.8069e-2, 3.3281e-4,-2.2865e-6/)
1064        real, dimension(3), parameter :: water_coefficients = (/ 1.0044, -1.1397e-2, 1.3300e-4 /)
1065       
1066        ! approx. fits from MODIS Collection 6 LUT scattering calculations
1067        if(phase == phaseIsLiquid) then
1068          get_ssa_nir = fit_to_quadratic(re, water_coefficients)
1069          if(re < re_water_min) get_ssa_nir = fit_to_quadratic(re_water_min, water_coefficients)
1070          if(re > re_water_max) get_ssa_nir = fit_to_quadratic(re_water_max, water_coefficients)
1071        else
1072          get_ssa_nir = fit_to_cubic(re, ice_coefficients)
1073          if(re < re_ice_min) get_ssa_nir = fit_to_cubic(re_ice_min, ice_coefficients)
1074          if(re > re_ice_max) get_ssa_nir = fit_to_cubic(re_ice_max, ice_coefficients)
1075        end if
1076
1077    end function get_ssa_nir
1078   ! --------------------------------------------
1079  pure function fit_to_cubic(x, coefficients)
1080    real,               intent(in) :: x
1081    real, dimension(:), intent(in) :: coefficients
1082    real                           :: fit_to_cubic
1083   
1084   
1085    fit_to_cubic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3) + x * coefficients(4)))
1086 end function fit_to_cubic
1087   ! --------------------------------------------
1088  pure function fit_to_quadratic(x, coefficients)
1089    real,               intent(in) :: x
1090    real, dimension(:), intent(in) :: coefficients
1091    real                           :: fit_to_quadratic
1092   
1093   
1094    fit_to_quadratic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3)))
1095 end function fit_to_quadratic
1096  ! --------------------------------------------
1097  ! Radiative transfer
1098  ! --------------------------------------------
1099  pure function compute_toa_reflectace(tau, g, w0)
1100    real, dimension(:), intent(in) :: tau, g, w0
1101    real                           :: compute_toa_reflectace
1102   
1103    logical, dimension(size(tau))         :: cloudMask
1104    integer, dimension(count(tau(:) > 0)) :: cloudIndicies
1105    real,    dimension(count(tau(:) > 0)) :: Refl,     Trans
1106    real                                  :: Refl_tot, Trans_tot
1107    integer                               :: i
1108    ! ---------------------------------------
1109    !
1110    ! This wrapper reports reflectance only and strips out non-cloudy elements from the calculation
1111    !
1112    cloudMask = tau(:) > 0.
1113    cloudIndicies = pack((/ (i, i = 1, size(tau)) /), mask = cloudMask)
1114    do i = 1, size(cloudIndicies)
1115      call two_stream(tau(cloudIndicies(i)), g(cloudIndicies(i)), w0(cloudIndicies(i)), Refl(i), Trans(i))
1116    end do
1117                   
1118    call adding_doubling(Refl(:), Trans(:), Refl_tot, Trans_tot) 
1119   
1120    compute_toa_reflectace = Refl_tot
1121   
1122  end function compute_toa_reflectace
1123  ! --------------------------------------------
1124  pure subroutine two_stream(tauint, gint, w0int, ref, tra)
1125    real, intent(in)  :: tauint, gint, w0int
1126    real, intent(out) :: ref, tra
1127    !
1128    ! Compute reflectance in a single layer using the two stream approximation
1129    !   The code itself is from Lazaros Oreopoulos via Steve Platnick
1130    !
1131    ! ------------------------
1132    ! Local variables
1133    !   for delta Eddington code
1134    !   xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1)
1135    integer, parameter :: beam = 2
1136    real,    parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
1137    real :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
1138            rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den, th
1139    !
1140    ! Compute reflectance and transmittance in a single layer using the two stream approximation
1141    !   The code itself is from Lazaros Oreopoulos via Steve Platnick
1142    !
1143    f   = gint**2
1144    tau = (1 - w0int * f) * tauint
1145    w0  = (1 - f) * w0int / (1 - w0int * f)
1146    g   = (gint - f) / (1 - f)
1147
1148    ! delta-Eddington (Joseph et al. 1976)
1149    gamma1 =  (7 - w0* (4 + 3 * g)) / 4.0
1150    gamma2 = -(1 - w0* (4 - 3 * g)) / 4.0
1151    gamma3 =  (2 - 3*g*xmu) / 4.0
1152    gamma4 =   1 - gamma3
1153
1154    if (w0int > minConservativeW0) then
1155      ! Conservative scattering
1156      if (beam == 1) then
1157          rh = (gamma1*tau+(gamma3-gamma1*xmu)*(1-exp(-tau/xmu)))
1158          ref = rh / (1 + gamma1 * tau)
1159          tra = 1 - ref       
1160      else if(beam == 2) then
1161          ref = gamma1*tau/(1 + gamma1*tau)
1162          tra = 1 - ref
1163      endif
1164    else
1165      ! Non-conservative scattering
1166      a1 = gamma1 * gamma4 + gamma2 * gamma3
1167      a2 = gamma1 * gamma3 + gamma2 * gamma4
1168
1169      rk = sqrt(gamma1**2 - gamma2**2)
1170     
1171      r1 = (1 - rk * xmu) * (a2 + rk * gamma3)
1172      r2 = (1 + rk * xmu) * (a2 - rk * gamma3)
1173      r3 = 2 * rk *(gamma3 - a2 * xmu)
1174      r4 = (1 - (rk * xmu)**2) * (rk + gamma1)
1175      r5 = (1 - (rk * xmu)**2) * (rk - gamma1)
1176     
1177      t1 = (1 + rk * xmu) * (a1 + rk * gamma4)
1178      t2 = (1 - rk * xmu) * (a1 - rk * gamma4)
1179      t3 = 2 * rk * (gamma4 + a1 * xmu)
1180      t4 = r4
1181      t5 = r5
1182
1183      beta = -r5 / r4         
1184     
1185      e1 = min(rk * tau, 500.)
1186      e2 = min(tau / xmu, 500.)
1187     
1188      if (beam == 1) then
1189         den = r4 * exp(e1) + r5 * exp(-e1)
1190         ref  = w0*(r1*exp(e1)-r2*exp(-e1)-r3*exp(-e2))/den
1191         den = t4 * exp(e1) + t5 * exp(-e1)
1192         th  = exp(-e2)
1193         tra = th-th*w0*(t1*exp(e1)-t2*exp(-e1)-t3*exp(e2))/den
1194      elseif (beam == 2) then
1195         ef1 = exp(-e1)
1196         ef2 = exp(-2*e1)
1197         ref = (gamma2*(1-ef2))/((rk+gamma1)*(1-beta*ef2))
1198         tra = (2*rk*ef1)/((rk+gamma1)*(1-beta*ef2))
1199      endif
1200    end if
1201  end subroutine two_stream
1202  ! --------------------------------------------------
1203  elemental function two_stream_reflectance(tauint, gint, w0int)
1204    real, intent(in) :: tauint, gint, w0int
1205    real             :: two_stream_reflectance
1206    !
1207    ! Compute reflectance in a single layer using the two stream approximation
1208    !   The code itself is from Lazaros Oreopoulos via Steve Platnick
1209    !
1210    ! ------------------------
1211    ! Local variables
1212    !   for delta Eddington code
1213    !   xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1)
1214    integer, parameter :: beam = 2
1215    real,    parameter :: xmu = 0.866, minConservativeW0 = 0.9999999
1216    real :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, &
1217            rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den
1218    ! ------------------------
1219
1220
1221    f   = gint**2
1222    tau = (1 - w0int * f) * tauint
1223    w0  = (1 - f) * w0int / (1 - w0int * f)
1224    g   = (gint - f) / (1 - f)
1225
1226    ! delta-Eddington (Joseph et al. 1976)
1227    gamma1 =  (7 - w0* (4 + 3 * g)) / 4.0
1228    gamma2 = -(1 - w0* (4 - 3 * g)) / 4.0
1229    gamma3 =  (2 - 3*g*xmu) / 4.0
1230    gamma4 =   1 - gamma3
1231
1232    if (w0int > minConservativeW0) then
1233      ! Conservative scattering
1234      if (beam == 1) then
1235          rh = (gamma1*tau+(gamma3-gamma1*xmu)*(1-exp(-tau/xmu)))
1236          two_stream_reflectance = rh / (1 + gamma1 * tau)
1237      elseif (beam == 2) then
1238          two_stream_reflectance = gamma1*tau/(1 + gamma1*tau)
1239      endif
1240       
1241    else    !
1242
1243        ! Non-conservative scattering
1244         a1 = gamma1 * gamma4 + gamma2 * gamma3
1245         a2 = gamma1 * gamma3 + gamma2 * gamma4
1246
1247         rk = sqrt(gamma1**2 - gamma2**2)
1248         
1249         r1 = (1 - rk * xmu) * (a2 + rk * gamma3)
1250         r2 = (1 + rk * xmu) * (a2 - rk * gamma3)
1251         r3 = 2 * rk *(gamma3 - a2 * xmu)
1252         r4 = (1 - (rk * xmu)**2) * (rk + gamma1)
1253         r5 = (1 - (rk * xmu)**2) * (rk - gamma1)
1254         
1255         t1 = (1 + rk * xmu) * (a1 + rk * gamma4)
1256         t2 = (1 - rk * xmu) * (a1 - rk * gamma4)
1257         t3 = 2 * rk * (gamma4 + a1 * xmu)
1258         t4 = r4
1259         t5 = r5
1260
1261         beta = -r5 / r4         
1262         
1263         e1 = min(rk * tau, 500.)
1264         e2 = min(tau / xmu, 500.)
1265         
1266         if (beam == 1) then
1267           den = r4 * exp(e1) + r5 * exp(-e1)
1268           two_stream_reflectance  = w0*(r1*exp(e1)-r2*exp(-e1)-r3*exp(-e2))/den
1269         elseif (beam == 2) then
1270           ef1 = exp(-e1)
1271           ef2 = exp(-2*e1)
1272           two_stream_reflectance = (gamma2*(1-ef2))/((rk+gamma1)*(1-beta*ef2))
1273         endif
1274           
1275      end if
1276  end function two_stream_reflectance
1277  ! --------------------------------------------
1278    pure subroutine adding_doubling (Refl, Tran, Refl_tot, Tran_tot)     
1279      real,    dimension(:), intent(in)  :: Refl,     Tran
1280      real,                  intent(out) :: Refl_tot, Tran_tot
1281      !
1282      ! Use adding/doubling formulas to compute total reflectance and transmittance from layer values
1283      !
1284     
1285      integer :: i
1286      real, dimension(size(Refl)) :: Refl_cumulative, Tran_cumulative
1287     
1288      Refl_cumulative(1) = Refl(1); Tran_cumulative(1) = Tran(1)   
1289     
1290      do i=2, size(Refl)
1291          ! place (add) previous combined layer(s) reflectance on top of layer i, w/black surface (or ignoring surface):
1292          Refl_cumulative(i) = Refl_cumulative(i-1) + Refl(i)*(Tran_cumulative(i-1)**2)/(1 - Refl_cumulative(i-1) * Refl(i))
1293          Tran_cumulative(i) = (Tran_cumulative(i-1)*Tran(i)) / (1 - Refl_cumulative(i-1) * Refl(i))
1294      end do
1295     
1296      Refl_tot = Refl_cumulative(size(Refl))
1297      Tran_tot = Tran_cumulative(size(Refl))
1298
1299    end subroutine adding_doubling
1300  ! --------------------------------------------------
1301  subroutine complain_and_die(message)
1302    character(len = *), intent(in) :: message
1303   
1304    write(6, *) "Failure in MODIS simulator"
1305    write(6, *)  trim(message)
1306    stop
1307  end subroutine complain_and_die
1308  !------------------------------------------------------------------------------------------------
1309end module mod_modis_sim
Note: See TracBrowser for help on using the repository browser.