source: LMDZ5/trunk/libf/phylmd/cosp/modis_simulator.F90 @ 2713

Last change on this file since 2713 was 2713, checked in by idelkadi, 8 years ago

Mise a jour du simulateur COSP.
Passage de la version 1.4.0 a la version 1.4.1
(Version suggeree pour CFMIP3/CMIP6)

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