source: LMDZ6/branches/DYNAMICO-conv-GC/libf/phylmd/cosp/modis_simulator.F90 @ 4543

Last change on this file since 4543 was 2822, checked in by idelkadi, 8 years ago
  1. Ecriture des champs de sorties Cosp 4D dans les fichiers Netcdf en utilisant les fonctionalites Xios.
  • ajout d'une routine histwrite_cosp4D appelee si les sorties sont gerees par Xios (dans cosp_output_write.F90)
  • definitions des differents axes et grilles 4D pour Cosp (dans context_lmdz.xml)
  1. Traitement des valeurs indefinies par Xios et mise a zero si l'on utilise Ioipsl.
  1. Correction d'erreur sur les sorties Modis
File size: 63.3 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
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    ! ---------------------------------------------------
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
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    ! ########################################################################################
534    ! Normalize pixel counts to fraction.
535    ! ########################################################################################
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   
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 
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      !
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
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
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
1010    else
1011       get_g_nir = fit_to_quadratic(re, ice_coefficients)
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        !
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 /)
1030       
1031        ! approx. fits from MODIS Collection 6 LUT scattering calculations
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  !------------------------------------------------------------------------------------------------
1274end module mod_modis_sim
Note: See TracBrowser for help on using the repository browser.