source: LMDZ6/trunk/libf/phylmd/cosp2/modis_simulator.F90 @ 5069

Last change on this file since 5069 was 3358, checked in by idelkadi, 6 years ago

Implementation de la nouvelle version COSPv2 dans LMDZ.
Pour compiler avec makelmdz_fcma utiliser l'option "-cosp2 true"

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