Index: LMDZ5/trunk/libf/phylmd/cosp/calc_Re.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/cosp/calc_Re.F90 (revision 2432) +++ LMDZ5/trunk/libf/phylmd/cosp/calc_Re.F90 (revision 2432) @@ -0,0 +1,262 @@ + subroutine calc_Re(Q,Np,rho_a, & + dtype,dmin,dmax,apm,bpm,rho_c,p1,p2,p3, & + Re) + use math_lib + implicit none + +! Purpose: +! Calculates Effective Radius (1/2 distribution 3rd moment / 2nd moment). +! +! For some distribution types, the total number concentration (per kg), Np +! may be optionally specified. Should be set to zero, otherwise. +! +! Roj Marchand July 2010 + + +! Inputs: +! +! [Q] hydrometeor mixing ratio (g/kg) ! not needed for some distribution types +! [Np] Optional Total number concentration (per kg). 0 = use defaults (p1, p2, p3) +! +! [rho_a] ambient air density (kg m^-3) +! +! Distribution parameters as per quickbeam documentation. +! [dtype] distribution type +! [dmin] minimum size cutoff (um) +! [dmax] maximum size cutoff (um) +! [apm] a parameter for mass (kg m^[-bpm]) +! [bmp] b params for mass +! [p1],[p2],[p3] distribution parameters +! +! +! Outputs: +! [Re] Effective radius, 1/2 the 3rd moment/2nd moment (um) +! +! Created: +! July 2010 Roj Marchand +! + + +! ----- INPUTS ----- + + real*8, intent(in) :: Q,Np,rho_a + + integer, intent(in):: dtype + real*8, intent(in) :: dmin,dmax,rho_c,p1,p2,p3 + + real*8, intent(inout) :: apm,bpm + +! ----- OUTPUTS ----- + + real*8, intent(out) :: Re + +! ----- INTERNAL ----- + + integer :: local_dtype + real*8 :: local_p3,local_Np + + real*8 :: pi, & + N0,D0,vu,dm,ld, & ! gamma, exponential variables + rg,log_sigma_g + + real*8 :: tmp1,tmp2 + + pi = acos(-1.0) + + ! // if density is constant, set equivalent values for apm and bpm + if ((rho_c > 0) .and. (apm < 0)) then + apm = (pi/6)*rho_c + bpm = 3. + endif + + ! Exponential is same as modified gamma with vu =1 + ! if Np is specified then we will just treat as modified gamma + if(dtype.eq.2 .and. Np>0) then + local_dtype=1; + local_p3=1; + else + local_dtype=dtype; + local_p3=p3; + endif + + select case(local_dtype) + +! ---------------------------------------------------------! +! // modified gamma ! +! ---------------------------------------------------------! +! :: Np = total number concentration (1/kg) = Nt / rho_a +! :: D0 = characteristic diameter (um) +! :: dm = mean diameter (um) - first moment over zeroth moment +! :: vu = distribution width parameter + + case(1) + + if( abs(local_p3+2) < 1E-8) then + + if(Np>1E-30) then + ! Morrison scheme with Martin 1994 shape parameter (NOTE: vu = pc +1) + ! fixed Roj. Dec. 2010 -- after comment by S. Mcfarlane + vu = (1/(0.2714 + 0.00057145*Np*rho_a*1E-6))**2 ! units of Nt = Np*rhoa = #/cm^3 + else + print *, 'Error: Must specify a value for Np in each volume', & + ' with Morrison/Martin Scheme.' + stop + endif + + elseif (abs(local_p3+1) > 1E-8) then + + ! vu is fixed in hp structure + vu = local_p3 + + else + + ! vu isn't specified + + print *, 'Error: Must specify a value for vu for Modified Gamma distribution' + stop + + endif + + + if( Np.eq.0 .and. p2+1 > 1E-8) then ! use default value for MEAN diameter as first default + + dm = p2 ! by definition, should have units of microns + D0 = gamma(vu)/gamma(vu+1)*dm + + else ! use value of Np + + if(Np.eq.0) then + + if( abs(p1+1) > 1E-8 ) then ! use default number concentration + + local_Np = p1 ! total number concentration / pa --- units kg^-1 + else + print *, 'Error: Must specify Np or default value ', & + '(p1=Dm [um] or p2=Np [1/kg]) for ', & + 'Modified Gamma distribution' + stop + endif + else + local_Np=Np; + endif + + D0 = 1E6 * ( Q*1E-3*gamma(vu)/(apm*local_Np*gamma(vu+bpm)) )**(1/bpm) ! units = microns + + endif + + Re = 0.5*D0*gamma(vu+3)/gamma(vu+2) + + +! ---------------------------------------------------------! +! // exponential ! +! ---------------------------------------------------------! +! :: N0 = intercept parameter (m^-4) +! :: ld = slope parameter (um) + + case(2) + + ! Np not specified (see if statement above) + + if((abs(p1+1) > 1E-8) ) then ! N0 has been specified, determine ld + + N0 = p1 + tmp1 = 1./(1.+bpm) + ld = ((apm*gamma(1.+bpm)*N0)/(rho_a*Q*1E-3))**tmp1 + ld = ld/1E6 ! set units to microns^-1 + + elseif (abs(p2+1) > 1E-8) then ! lambda=ld has been specified as default + + ld = p2 ! should have units of microns^-1 + + else + + print *, 'Error: Must specify Np or default value ', & + '(p1=No or p2=lambda) for Exponential distribution' + stop + + endif + + Re = 1.5/ld ; + +! ---------------------------------------------------------! +! // power law ! +! ---------------------------------------------------------! +! :: ahp = Ar parameter (m^-4 mm^-bhp) +! :: bhp = br parameter +! :: dmin_mm = lower bound (mm) +! :: dmax_mm = upper bound (mm) + + case(3) + + Re=0; ! Not supporting LUT approach for power-law ... + + if(Np>0) then + + print *, 'Variable Np not supported for ', & + 'Power Law distribution' + stop + endif +! ---------------------------------------------------------! +! // monodisperse ! +! ---------------------------------------------------------! +! :: D0 = particle diameter (um) == Re + + case(4) + + Re = p1 + + if(Np>0) then + print *, 'Variable Np not supported for ', & + 'Monodispersed distribution' + stop + endif + +! ---------------------------------------------------------! +! // lognormal ! +! ---------------------------------------------------------! +! :: N0 = total number concentration (m^-3) +! :: np = fixed number concentration (kg^-1) +! :: rg = mean radius (um) +! :: log_sigma_g = ln(geometric standard deviation) + + case(5) + + if( abs(local_p3+1) > 1E-8 ) then + + !set natural log width + log_sigma_g = local_p3 + else + print *, 'Error: Must specify a value for sigma_g ', & + 'when using a Log-Normal distribution' + stop + endif + + ! get rg ... + if( Np.eq.0 .and. (abs(p2+1) > 1E-8) ) then ! use default value of rg + + rg = p2 + + else + if(Np>0) then + + local_Np=Np; + + elseif(abs(p2+1) < 1E-8) then + + local_Np=p1 + else + print *, 'Error: Must specify Np or default value ', & + '(p2=Rg or p1=Np) for Log-Normal distribution' + endif + + log_sigma_g = p3 + tmp1 = (Q*1E-3)/(2.**bpm*apm*local_Np) + tmp2 = exp(0.5*bpm**2.*(log_sigma_g))**2. + rg = ((tmp1/tmp2)**(1/bpm))*1E6 + endif + + Re = rg*exp(+2.5*(log_sigma_g**2)) + + end select + + end subroutine calc_Re Index: LMDZ5/trunk/libf/phylmd/cosp/cosp_modis_simulator.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/cosp/cosp_modis_simulator.F90 (revision 2432) +++ LMDZ5/trunk/libf/phylmd/cosp/cosp_modis_simulator.F90 (revision 2432) @@ -0,0 +1,453 @@ +! (c) 2009, Regents of the Unversity of Colorado +! Author: Robert Pincus, Cooperative Institute for Research in the Environmental Sciences +! All rights reserved. +! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $ +! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/cosp_modis_simulator.F90 $ +! +! Redistribution and use in source and binary forms, with or without modification, are permitted +! provided that the following conditions are met: +! +! * Redistributions of source code must retain the above copyright notice, this list +! of conditions and the following disclaimer. +! * Redistributions in binary form must reproduce the above copyright notice, this list +! of conditions and the following disclaimer in the documentation and/or other materials +! provided with the distribution. +! * Neither the name of the Met Office nor the names of its contributors may be used +! to endorse or promote products derived from this software without specific prior written +! permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR +! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR +! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER +! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! + +! +! History: +! May 2009 - Robert Pincus - Initial version +! Dec 2009 - Robert Pincus - Tiny revisions +! +MODULE MOD_COSP_Modis_Simulator + USE MOD_COSP_CONSTANTS + USE MOD_COSP_TYPES + use mod_modis_sim, numModisTauBins => numTauHistogramBins, & + numModisPressureBins => numPressureHistogramBins, & + MODIS_TAU => nominalTauHistogramCenters, & + MODIS_TAU_BNDS => nominalTauHistogramBoundaries, & + MODIS_PC => nominalPressureHistogramCenters, & + MODIS_PC_BNDS => nominalPressureHistogramBoundaries + implicit none + !------------------------------------------------------------------------------------------------ + ! Public type + ! + ! Summary statistics from MODIS retrievals + type COSP_MODIS + ! Dimensions + integer :: Npoints ! Number of gridpoints + + ! + ! Grid means; dimension nPoints + ! + real, dimension(:), pointer :: & + Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, & + Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, & + Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, & + Optical_Thickness_Total_LogMean, Optical_Thickness_Water_LogMean, Optical_Thickness_Ice_LogMean, & + Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, & + Cloud_Top_Pressure_Total_Mean, & + Liquid_Water_Path_Mean, Ice_Water_Path_Mean + ! + ! Also need the ISCCP-type optical thickness/cloud top pressure histogram + ! + real, dimension(:, :, :), pointer :: Optical_Thickness_vs_Cloud_Top_Pressure + end type COSP_MODIS + +contains + !------------------------------------------------------------------------------------------------ + subroutine COSP_Modis_Simulator(gridBox, subCols, subcolHydro, isccpSim, modisSim) + ! Arguments + type(cosp_gridbox), intent(in ) :: gridBox ! Gridbox info + type(cosp_subgrid), intent(in ) :: subCols ! subCol indicators of convective/stratiform + type(cosp_sghydro), intent(in ) :: subcolHydro ! subcol hydrometeor contens + type(cosp_isccp), intent(in ) :: isccpSim ! ISCCP simulator output + type(cosp_modis), intent( out) :: modisSim ! MODIS simulator subcol output + + ! ------------------------------------------------------------ + ! Local variables + ! Leave space only for sunlit points + + integer :: nPoints, nSubCols, nLevels, nSunlit, i, j, k + + ! Grid-mean quanties; dimensions nPoints, nLevels + real, & + dimension(count(gridBox%sunlit(:) > 0), gridBox%nLevels) :: & + temperature, pressureLayers + real, & + dimension(count(gridBox%sunlit(:) > 0), gridBox%nLevels + 1) :: & + pressureLevels + + ! Subcol quantities, dimension nPoints, nSubCols, nLevels + real, & + dimension(count(gridBox%sunlit(:) > 0), subCols%nColumns, gridBox%nLevels) :: & + opticalThickness, cloudWater, cloudIce, waterSize, iceSize + + ! Vertically-integrated subcol quantities; dimensions nPoints, nSubcols + integer, & + dimension(count(gridBox%sunlit(:) > 0), subCols%nColumns) :: & + retrievedPhase + real, & + dimension(count(gridBox%sunlit(:) > 0), subCols%nColumns) :: & + isccpTau, isccpCloudTopPressure, retrievedCloudTopPressure, retrievedTau, retrievedSize + + ! Vertically-integrated results + real, dimension(count(gridBox%sunlit(:) > 0)) :: & + cfTotal, cfLiquid, cfIce, & + cfHigh, cfMid, cfLow, & + meanTauTotal, meanTauLiquid, meanTauIce, & + meanLogTauTotal, meanLogTauLiquid, meanLogTauIce , & + meanSizeLiquid, meanSizeIce, & + meanCloudTopPressure, & + meanLiquidWaterPath, meanIceWaterPath + + real, dimension(count(gridBox%sunlit(:) > 0), numModisTauBins, numModisPressureBins) :: & + jointHistogram + + integer, dimension(count(gridBox%sunlit(:) > 0)) :: sunlit + integer, dimension(count(gridBox%sunlit(:) <= 0)) :: notSunlit + ! ------------------------------------------------------------ + + ! + ! Are there any sunlit points? + ! + nSunlit = count(gridBox%sunlit(:) > 0) + if(nSunlit > 0) then + nLevels = gridBox%Nlevels + nPoints = gridBox%Npoints + nSubCols = subCols%Ncolumns + ! + ! This is a vector index indicating which points are sunlit + ! + sunlit(:) = pack((/ (i, i = 1, nPoints ) /), mask = gridBox%sunlit(:) > 0) + notSunlit(:) = pack((/ (i, i = 1, nPoints ) /), mask = .not. gridBox%sunlit(:) > 0) + + ! + ! Copy needed quantities, reversing vertical order and removing points with no sunlight + ! + pressureLevels(:, 1) = 0.0 ! Top of model, following ISCCP sim + temperature(:, :) = gridBox%T (sunlit(:), nLevels:1:-1) + pressureLayers(:, :) = gridBox%p (sunlit(:), nLevels:1:-1) + pressureLevels(:, 2:) = gridBox%ph(sunlit(:), nLevels:1:-1) + + ! + ! Subcolumn properties - first stratiform cloud... + ! + where(subCols%frac_out(sunlit(:), :, :) == I_LSC) + !opticalThickness(:, :, :) = & + ! spread(gridBox%dtau_s (sunlit(:), :), dim = 2, nCopies = nSubCols) + cloudWater(:, :, :) = subcolHydro%mr_hydro(sunlit(:), :, :, I_LSCLIQ) + waterSize (:, :, :) = subcolHydro%reff (sunlit(:), :, :, I_LSCLIQ) + cloudIce (:, :, :) = subcolHydro%mr_hydro(sunlit(:), :, :, I_LSCICE) + iceSize (:, :, :) = subcolHydro%reff (sunlit(:), :, :, I_LSCICE) + elsewhere + opticalThickness(:, :, :) = 0. + cloudWater (:, :, :) = 0. + cloudIce (:, :, :) = 0. + waterSize (:, :, :) = 0. + iceSize (:, :, :) = 0. + end where + + ! Loop version of spread above - intrinsic doesn't work on certain platforms. + do k = 1, nLevels + do j = 1, nSubCols + do i = 1, nSunlit + if(subCols%frac_out(sunlit(i), j, k) == I_LSC) then + opticalThickness(i, j, k) = gridBox%dtau_s(sunlit(i), k) + else + opticalThickness(i, j, k) = 0. + end if + end do + end do + end do + + ! + ! .. then add convective cloud... + ! + where(subCols%frac_out(sunlit(:), :, :) == I_CVC) + !opticalThickness(:, :, :) = & + ! spread(gridBox%dtau_c( sunlit(:), :), dim = 2, nCopies = nSubCols) + cloudWater(:, :, :) = subcolHydro%mr_hydro(sunlit(:), :, :, I_CVCLIQ) + waterSize (:, :, :) = subcolHydro%reff (sunlit(:), :, :, I_CVCLIQ) + cloudIce (:, :, :) = subcolHydro%mr_hydro(sunlit(:), :, :, I_CVCICE) + iceSize (:, :, :) = subcolHydro%reff (sunlit(:), :, :, I_CVCICE) + end where + + ! Loop version of spread above - intrinsic doesn't work on certain platforms. + do k = 1, nLevels + do j = 1, nSubCols + do i = 1, nSunlit + if(subCols%frac_out(sunlit(i), j, k) == I_CVC) opticalThickness(i, j, k) = gridBox%dtau_c(sunlit(i), k) + end do + end do + end do + + ! + ! Reverse vertical order + ! + opticalThickness(:, :, :) = opticalThickness(:, :, nLevels:1:-1) + cloudWater (:, :, :) = cloudWater (:, :, nLevels:1:-1) + waterSize (:, :, :) = waterSize (:, :, nLevels:1:-1) + cloudIce (:, :, :) = cloudIce (:, :, nLevels:1:-1) + iceSize (:, :, :) = iceSize (:, :, nLevels:1:-1) + + isccpTau(:, :) = isccpSim%boxtau (sunlit(:), :) + isccpCloudTopPressure(:, :) = isccpSim%boxptop(sunlit(:), :) + + do i = 1, nSunlit + call modis_L2_simulator(temperature(i, :), pressureLayers(i, :), pressureLevels(i, :), & + opticalThickness(i, :, :), cloudWater(i, :, :), cloudIce(i, :, :), & + waterSize(i, :, :), iceSize(i, :, :), & + isccpTau(i, :), isccpCloudTopPressure(i, :), & + retrievedPhase(i, :), retrievedCloudTopPressure(i, :), & + retrievedTau(i, :), retrievedSize(i, :)) + end do + call modis_L3_simulator(retrievedPhase, & + retrievedCloudTopPressure, & + retrievedTau, retrievedSize, & + cfTotal, cfLiquid, cfIce, & + cfHigh, cfMid, cfLow, & + meanTauTotal, meanTauLiquid, meanTauIce, & + meanLogTauTotal, meanLogTauLiquid, meanLogTauIce, & + meanSizeLiquid, meanSizeIce, & + meanCloudTopPressure, & + meanLiquidWaterPath, meanIceWaterPath, & + jointHistogram) + ! + ! Copy results into COSP structure + ! + modisSim%Cloud_Fraction_Total_Mean(sunlit(:)) = cfTotal(:) + modisSim%Cloud_Fraction_Water_Mean(sunlit(:)) = cfLiquid + modisSim%Cloud_Fraction_Ice_Mean (sunlit(:)) = cfIce + + modisSim%Cloud_Fraction_High_Mean(sunlit(:)) = cfHigh + modisSim%Cloud_Fraction_Mid_Mean (sunlit(:)) = cfMid + modisSim%Cloud_Fraction_Low_Mean (sunlit(:)) = cfLow + + modisSim%Optical_Thickness_Total_Mean(sunlit(:)) = meanTauTotal + modisSim%Optical_Thickness_Water_Mean(sunlit(:)) = meanTauLiquid + modisSim%Optical_Thickness_Ice_Mean (sunlit(:)) = meanTauIce + + modisSim%Optical_Thickness_Total_LogMean(sunlit(:)) = meanLogTauTotal + modisSim%Optical_Thickness_Water_LogMean(sunlit(:)) = meanLogTauLiquid + modisSim%Optical_Thickness_Ice_LogMean (sunlit(:)) = meanLogTauIce + + modisSim%Cloud_Particle_Size_Water_Mean(sunlit(:)) = meanSizeLiquid + modisSim%Cloud_Particle_Size_Ice_Mean (sunlit(:)) = meanSizeIce + + modisSim%Cloud_Top_Pressure_Total_Mean(sunlit(:)) = meanCloudTopPressure + + modisSim%Liquid_Water_Path_Mean(sunlit(:)) = meanLiquidWaterPath + modisSim%Ice_Water_Path_Mean (sunlit(:)) = meanIceWaterPath + + modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(sunlit(:), 2:numModisTauBins+1, :) = jointHistogram(:, :, :) + ! + ! Reorder pressure bins in joint histogram to go from surface to TOA + ! + modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(:,:,:) = & + modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(:, :, numModisPressureBins:1:-1) + if(nSunlit < nPoints) then + ! + ! Where it's night and we haven't done the retrievals the values are undefined + ! + modisSim%Cloud_Fraction_Total_Mean(notSunlit(:)) = R_UNDEF + modisSim%Cloud_Fraction_Water_Mean(notSunlit(:)) = R_UNDEF + modisSim%Cloud_Fraction_Ice_Mean (notSunlit(:)) = R_UNDEF + + modisSim%Cloud_Fraction_High_Mean(notSunlit(:)) = R_UNDEF + modisSim%Cloud_Fraction_Mid_Mean (notSunlit(:)) = R_UNDEF + modisSim%Cloud_Fraction_Low_Mean (notSunlit(:)) = R_UNDEF + + modisSim%Optical_Thickness_Total_Mean(notSunlit(:)) = R_UNDEF + modisSim%Optical_Thickness_Water_Mean(notSunlit(:)) = R_UNDEF + modisSim%Optical_Thickness_Ice_Mean (notSunlit(:)) = R_UNDEF + + modisSim%Optical_Thickness_Total_LogMean(notSunlit(:)) = R_UNDEF + modisSim%Optical_Thickness_Water_LogMean(notSunlit(:)) = R_UNDEF + modisSim%Optical_Thickness_Ice_LogMean (notSunlit(:)) = R_UNDEF + + modisSim%Cloud_Particle_Size_Water_Mean(notSunlit(:)) = R_UNDEF + modisSim%Cloud_Particle_Size_Ice_Mean (notSunlit(:)) = R_UNDEF + + modisSim%Cloud_Top_Pressure_Total_Mean(notSunlit(:)) = R_UNDEF + + modisSim%Liquid_Water_Path_Mean(notSunlit(:)) = R_UNDEF + modisSim%Ice_Water_Path_Mean (notSunlit(:)) = R_UNDEF + + modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(notSunlit(:), :, :) = R_UNDEF + end if + else + ! + ! It's nightime everywhere - everything is undefined + ! + modisSim%Cloud_Fraction_Total_Mean(:) = R_UNDEF + modisSim%Cloud_Fraction_Water_Mean(:) = R_UNDEF + modisSim%Cloud_Fraction_Ice_Mean (:) = R_UNDEF + + modisSim%Cloud_Fraction_High_Mean(:) = R_UNDEF + modisSim%Cloud_Fraction_Mid_Mean (:) = R_UNDEF + modisSim%Cloud_Fraction_Low_Mean (:) = R_UNDEF + + modisSim%Optical_Thickness_Total_Mean(:) = R_UNDEF + modisSim%Optical_Thickness_Water_Mean(:) = R_UNDEF + modisSim%Optical_Thickness_Ice_Mean (:) = R_UNDEF + + modisSim%Optical_Thickness_Total_LogMean(:) = R_UNDEF + modisSim%Optical_Thickness_Water_LogMean(:) = R_UNDEF + modisSim%Optical_Thickness_Ice_LogMean (:) = R_UNDEF + + modisSim%Cloud_Particle_Size_Water_Mean(:) = R_UNDEF + modisSim%Cloud_Particle_Size_Ice_Mean (:) = R_UNDEF + + modisSim%Cloud_Top_Pressure_Total_Mean(:) = R_UNDEF + + modisSim%Liquid_Water_Path_Mean(:) = R_UNDEF + modisSim%Ice_Water_Path_Mean (:) = R_UNDEF + + modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(:, :, :) = R_UNDEF + end if + + end subroutine COSP_Modis_Simulator + !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + !------------- SUBROUTINE CONSTRUCT_COSP_MODIS ------------------ + !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + SUBROUTINE CONSTRUCT_COSP_MODIS(cfg, nPoints, x) + type(cosp_config), intent(in) :: cfg ! Configuration options + integer, intent(in) :: Npoints ! Number of sampled points + type(cosp_MODIS), intent(out) :: x + ! + ! Allocate minumum storage if simulator not used + ! + if (cfg%LMODIS_sim) then + x%nPoints = nPoints + else + x%Npoints = 1 + endif + + ! --- Allocate arrays --- + allocate(x%Cloud_Fraction_Total_Mean(x%nPoints)) + allocate(x%Cloud_Fraction_Water_Mean(x%nPoints)) + allocate(x%Cloud_Fraction_Ice_Mean(x%nPoints)) + + allocate(x%Cloud_Fraction_High_Mean(x%nPoints)) + allocate(x%Cloud_Fraction_Mid_Mean(x%nPoints)) + allocate(x%Cloud_Fraction_Low_Mean(x%nPoints)) + + allocate(x%Optical_Thickness_Total_Mean(x%nPoints)) + allocate(x%Optical_Thickness_Water_Mean(x%nPoints)) + allocate(x%Optical_Thickness_Ice_Mean(x%nPoints)) + + allocate(x%Optical_Thickness_Total_LogMean(x%nPoints)) + allocate(x%Optical_Thickness_Water_LogMean(x%nPoints)) + allocate(x%Optical_Thickness_Ice_LogMean(x%nPoints)) + + allocate(x%Cloud_Particle_Size_Water_Mean(x%nPoints)) + allocate(x%Cloud_Particle_Size_Ice_Mean(x%nPoints)) + + allocate(x%Cloud_Top_Pressure_Total_Mean(x%nPoints)) + + allocate(x%Liquid_Water_Path_Mean(x%nPoints)) + allocate(x%Ice_Water_Path_Mean(x%nPoints)) + + allocate(x%Optical_Thickness_vs_Cloud_Top_Pressure(nPoints, numModisTauBins+1, numModisPressureBins)) + x%Optical_Thickness_vs_Cloud_Top_Pressure(:, :, :) = R_UNDEF + END SUBROUTINE CONSTRUCT_COSP_MODIS + + !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + !------------- SUBROUTINE FREE_COSP_MODIS ----------------------- + !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + SUBROUTINE FREE_COSP_MODIS(x) + type(cosp_MODIS),intent(inout) :: x + ! + ! Free space used by cosp_modis variable. + ! + + if(associated(x%Cloud_Fraction_Total_Mean)) deallocate(x%Cloud_Fraction_Total_Mean) + if(associated(x%Cloud_Fraction_Water_Mean)) deallocate(x%Cloud_Fraction_Water_Mean) + if(associated(x%Cloud_Fraction_Ice_Mean )) deallocate(x%Cloud_Fraction_Ice_Mean) + + if(associated(x%Cloud_Fraction_High_Mean)) deallocate(x%Cloud_Fraction_High_Mean) + if(associated(x%Cloud_Fraction_Mid_Mean )) deallocate(x%Cloud_Fraction_Mid_Mean) + if(associated(x%Cloud_Fraction_Low_Mean )) deallocate(x%Cloud_Fraction_Low_Mean) + + if(associated(x%Optical_Thickness_Total_Mean)) deallocate(x%Optical_Thickness_Total_Mean) + if(associated(x%Optical_Thickness_Water_Mean)) deallocate(x%Optical_Thickness_Water_Mean) + if(associated(x%Optical_Thickness_Ice_Mean )) deallocate(x%Optical_Thickness_Ice_Mean) + + if(associated(x%Optical_Thickness_Total_LogMean)) deallocate(x%Optical_Thickness_Total_LogMean) + if(associated(x%Optical_Thickness_Water_LogMean)) deallocate(x%Optical_Thickness_Water_LogMean) + if(associated(x%Optical_Thickness_Ice_LogMean )) deallocate(x%Optical_Thickness_Ice_LogMean) + + if(associated(x%Cloud_Particle_Size_Water_Mean)) deallocate(x%Cloud_Particle_Size_Water_Mean) + if(associated(x%Cloud_Particle_Size_Ice_Mean )) deallocate(x%Cloud_Particle_Size_Ice_Mean) + + if(associated(x%Cloud_Top_Pressure_Total_Mean )) deallocate(x%Cloud_Top_Pressure_Total_Mean ) + + if(associated(x%Liquid_Water_Path_Mean)) deallocate(x%Liquid_Water_Path_Mean ) + if(associated(x%Ice_Water_Path_Mean )) deallocate(x%Ice_Water_Path_Mean ) + + if(associated(x%Optical_Thickness_vs_Cloud_Top_Pressure)) deallocate(x%Optical_Thickness_vs_Cloud_Top_Pressure ) + END SUBROUTINE FREE_COSP_MODIS + ! ----------------------------------------------------- + + !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + !------------- SUBROUTINE COSP_MODIS_CPSECTION ----------------- + !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + SUBROUTINE COSP_MODIS_CPSECTION(ix, iy, orig, copy) + integer, dimension(2), intent(in) :: ix, iy + type(cosp_modis), intent(in ) :: orig + type(cosp_modis), intent( out) :: copy + ! + ! Copy a set of grid points from one cosp_modis variable to another. + ! Should test to be sure ix and iy refer to the same number of grid points + ! + integer :: orig_start, orig_end, copy_start, copy_end + + orig_start = ix(1); orig_end = ix(2) + copy_start = iy(1); copy_end = iy(2) + + copy%Cloud_Fraction_Total_Mean(copy_start:copy_end) = orig%Cloud_Fraction_Total_Mean(orig_start:orig_end) + copy%Cloud_Fraction_Water_Mean(copy_start:copy_end) = orig%Cloud_Fraction_Water_Mean(orig_start:orig_end) + copy%Cloud_Fraction_Ice_Mean (copy_start:copy_end) = orig%Cloud_Fraction_Ice_Mean (orig_start:orig_end) + + copy%Cloud_Fraction_High_Mean(copy_start:copy_end) = orig%Cloud_Fraction_High_Mean(orig_start:orig_end) + copy%Cloud_Fraction_Mid_Mean (copy_start:copy_end) = orig%Cloud_Fraction_Mid_Mean (orig_start:orig_end) + copy%Cloud_Fraction_Low_Mean (copy_start:copy_end) = orig%Cloud_Fraction_Low_Mean (orig_start:orig_end) + + copy%Optical_Thickness_Total_Mean(copy_start:copy_end) = orig%Optical_Thickness_Total_Mean(orig_start:orig_end) + copy%Optical_Thickness_Water_Mean(copy_start:copy_end) = orig%Optical_Thickness_Water_Mean(orig_start:orig_end) + copy%Optical_Thickness_Ice_Mean (copy_start:copy_end) = orig%Optical_Thickness_Ice_Mean (orig_start:orig_end) + + copy%Optical_Thickness_Total_LogMean(copy_start:copy_end) = & + orig%Optical_Thickness_Total_LogMean(orig_start:orig_end) + copy%Optical_Thickness_Water_LogMean(copy_start:copy_end) = & + orig%Optical_Thickness_Water_LogMean(orig_start:orig_end) + copy%Optical_Thickness_Ice_LogMean (copy_start:copy_end) = & + orig%Optical_Thickness_Ice_LogMean (orig_start:orig_end) + + copy%Cloud_Particle_Size_Water_Mean(copy_start:copy_end) = orig%Cloud_Particle_Size_Water_Mean(orig_start:orig_end) + copy%Cloud_Particle_Size_Ice_Mean (copy_start:copy_end) = orig%Cloud_Particle_Size_Ice_Mean (orig_start:orig_end) + + copy%Cloud_Top_Pressure_Total_Mean(copy_start:copy_end) = orig%Cloud_Top_Pressure_Total_Mean(orig_start:orig_end) + + copy%Liquid_Water_Path_Mean(copy_start:copy_end) = orig%Liquid_Water_Path_Mean(orig_start:orig_end) + copy%Ice_Water_Path_Mean (copy_start:copy_end) = orig%Ice_Water_Path_Mean (orig_start:orig_end) + + copy%Optical_Thickness_vs_Cloud_Top_Pressure(copy_start:copy_end, :, :) = & + orig%Optical_Thickness_vs_Cloud_Top_Pressure(orig_start:orig_end, :, :) + END SUBROUTINE COSP_MODIS_CPSECTION + ! ----------------------------------------------------- + +END MODULE MOD_COSP_Modis_Simulator Index: LMDZ5/trunk/libf/phylmd/cosp/isccp_cloud_types.F =================================================================== --- LMDZ5/trunk/libf/phylmd/cosp/isccp_cloud_types.F (revision 2432) +++ LMDZ5/trunk/libf/phylmd/cosp/isccp_cloud_types.F (revision 2432) @@ -0,0 +1,347 @@ +! $Revision: 23 $, $Date: 2011-03-31 15:41:37 +0200 (jeu. 31 mars 2011) $ +! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/icarus-scops-4.1-bsd/isccp_cloud_types.f $ + SUBROUTINE ISCCP_CLOUD_TYPES( + & debug, + & debugcol, + & npoints, + & sunlit, + & nlev, + & ncol, + & seed, + & pfull, + & phalf, + & qv, + & cc, + & conv, + & dtau_s, + & dtau_c, + & top_height, + & top_height_direction, + & overlap, + & frac_out, + & skt, + & emsfc_lw, + & at, + & dem_s, + & dem_c, + & fq_isccp, + & totalcldarea, + & meanptop, + & meantaucld, + & meanalbedocld, + & meantb, + & meantbclr, + & boxtau, + & boxptop + &) + +!$Id: isccp_cloud_types.f,v 4.0 2009/03/06 11:05:11 hadmw Exp $ + +! *****************************COPYRIGHT**************************** +! (c) British Crown Copyright 2009, the Met Office. +! All rights reserved. +! +! Redistribution and use in source and binary forms, with or without +! modification, are permitted provided that the +! following conditions are met: +! +! * Redistributions of source code must retain the above +! copyright notice, this list of conditions and the following +! disclaimer. +! * Redistributions in binary form must reproduce the above +! copyright notice, this list of conditions and the following +! disclaimer in the documentation and/or other materials +! provided with the distribution. +! * Neither the name of the Met Office nor the names of its +! contributors may be used to endorse or promote products +! derived from this software without specific prior written +! permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! *****************************COPYRIGHT******************************* +! *****************************COPYRIGHT******************************* +! *****************************COPYRIGHT******************************* + + implicit none + +! NOTE: the maximum number of levels and columns is set by +! the following parameter statement + + INTEGER ncolprint + +! ----- +! Input +! ----- + + INTEGER npoints ! number of model points in the horizontal + INTEGER nlev ! number of model levels in column + INTEGER ncol ! number of subcolumns + + INTEGER sunlit(npoints) ! 1 for day points, 0 for night time + + INTEGER seed(npoints) + ! seed values for marsaglia random number generator + ! It is recommended that the seed is set + ! to a different value for each model + ! gridbox it is called on, as it is + ! possible that the choice of the same + ! seed value every time may introduce some + ! statistical bias in the results, particularly + ! for low values of NCOL. + + REAL pfull(npoints,nlev) + ! pressure of full model levels (Pascals) + ! pfull(npoints,1) is top level of model + ! pfull(npoints,nlev) is bot of model + + REAL phalf(npoints,nlev+1) + ! pressure of half model levels (Pascals) + ! phalf(npoints,1) is top of model + ! phalf(npoints,nlev+1) is the surface pressure + + REAL qv(npoints,nlev) + ! water vapor specific humidity (kg vapor/ kg air) + ! on full model levels + + REAL cc(npoints,nlev) + ! input cloud cover in each model level (fraction) + ! NOTE: This is the HORIZONTAL area of each + ! grid box covered by clouds + + REAL conv(npoints,nlev) + ! input convective cloud cover in each model + ! level (fraction) + ! NOTE: This is the HORIZONTAL area of each + ! grid box covered by convective clouds + + REAL dtau_s(npoints,nlev) + ! mean 0.67 micron optical depth of stratiform + ! clouds in each model level + ! NOTE: this the cloud optical depth of only the + ! cloudy part of the grid box, it is not weighted + ! with the 0 cloud optical depth of the clear + ! part of the grid box + + REAL dtau_c(npoints,nlev) + ! mean 0.67 micron optical depth of convective + ! clouds in each + ! model level. Same note applies as in dtau_s. + + INTEGER overlap ! overlap type + ! 1=max + ! 2=rand + ! 3=max/rand + + INTEGER top_height ! 1 = adjust top height using both a computed + ! infrared brightness temperature and the visible + ! optical depth to adjust cloud top pressure. Note + ! that this calculation is most appropriate to compare + ! to ISCCP data during sunlit hours. + ! 2 = do not adjust top height, that is cloud top + ! pressure is the actual cloud top pressure + ! in the model + ! 3 = adjust top height using only the computed + ! infrared brightness temperature. Note that this + ! calculation is most appropriate to compare to ISCCP + ! IR only algortihm (i.e. you can compare to nighttime + ! ISCCP data with this option) + + INTEGER top_height_direction ! direction for finding atmosphere pressure level + ! with interpolated temperature equal to the radiance + ! determined cloud-top temperature + ! + ! 1 = find the *lowest* altitude (highest pressure) level + ! with interpolated temperature equal to the radiance + ! determined cloud-top temperature + ! + ! 2 = find the *highest* altitude (lowest pressure) level + ! with interpolated temperature equal to the radiance + ! determined cloud-top temperature + ! + ! ONLY APPLICABLE IF top_height EQUALS 1 or 3 + ! + ! 1 = old setting: matches all versions of + ! ISCCP simulator with versions numbers 3.5.1 and lower + ! + ! 2 = default setting: for version numbers 4.0 and higher +! +! The following input variables are used only if top_height = 1 or top_height = 3 +! + REAL skt(npoints) ! skin Temperature (K) + REAL emsfc_lw ! 10.5 micron emissivity of surface (fraction) + REAL at(npoints,nlev) ! temperature in each model level (K) + REAL dem_s(npoints,nlev) ! 10.5 micron longwave emissivity of stratiform + ! clouds in each + ! model level. Same note applies as in dtau_s. + REAL dem_c(npoints,nlev) ! 10.5 micron longwave emissivity of convective + ! clouds in each + ! model level. Same note applies as in dtau_s. + + REAL frac_out(npoints,ncol,nlev) ! boxes gridbox divided up into + ! Equivalent of BOX in original version, but + ! indexed by column then row, rather than + ! by row then column + + + +! ------ +! Output +! ------ + + REAL fq_isccp(npoints,7,7) ! the fraction of the model grid box covered by + ! each of the 49 ISCCP D level cloud types + + REAL totalcldarea(npoints) ! the fraction of model grid box columns + ! with cloud somewhere in them. NOTE: This diagnostic + ! does not count model clouds with tau < isccp_taumin + ! Thus this diagnostic does not equal the sum over all entries of fq_isccp. + ! However, this diagnostic does equal the sum over entries of fq_isccp with + ! itau = 2:7 (omitting itau = 1) + + + ! The following three means are averages only over the cloudy areas with tau > isccp_taumin. + ! If no clouds with tau > isccp_taumin are in grid box all three quantities should equal zero. + + REAL meanptop(npoints) ! mean cloud top pressure (mb) - linear averaging + ! in cloud top pressure. + + REAL meantaucld(npoints) ! mean optical thickness + ! linear averaging in albedo performed. + + real meanalbedocld(npoints) ! mean cloud albedo + ! linear averaging in albedo performed + + real meantb(npoints) ! mean all-sky 10.5 micron brightness temperature + + real meantbclr(npoints) ! mean clear-sky 10.5 micron brightness temperature + + REAL boxtau(npoints,ncol) ! optical thickness in each column + + REAL boxptop(npoints,ncol) ! cloud top pressure (mb) in each column + + +! +! ------ +! Working variables added when program updated to mimic Mark Webb's PV-Wave code +! ------ + + REAL dem(npoints,ncol),bb(npoints) ! working variables for 10.5 micron longwave + ! emissivity in part of + ! gridbox under consideration + + REAL ptrop(npoints) + REAL attrop(npoints) + REAL attropmin (npoints) + REAL atmax(npoints) + REAL atmin(npoints) + REAL btcmin(npoints) + REAL transmax(npoints) + + INTEGER i,j,ilev,ibox,itrop(npoints) + INTEGER ipres(npoints) + INTEGER itau(npoints),ilev2 + INTEGER acc(nlev,ncol) + INTEGER match(npoints,nlev-1) + INTEGER nmatch(npoints) + INTEGER levmatch(npoints,ncol) + + !variables needed for water vapor continuum absorption + real fluxtop_clrsky(npoints),trans_layers_above_clrsky(npoints) + real taumin(npoints) + real dem_wv(npoints,nlev), wtmair, wtmh20, Navo, grav, pstd, t0 + real press(npoints), dpress(npoints), atmden(npoints) + real rvh20(npoints), wk(npoints), rhoave(npoints) + real rh20s(npoints), rfrgn(npoints) + real tmpexp(npoints),tauwv(npoints) + + character*1 cchar(6),cchar_realtops(6) + integer icycle + REAL tau(npoints,ncol) + LOGICAL box_cloudy(npoints,ncol) + REAL tb(npoints,ncol) + REAL ptop(npoints,ncol) + REAL emcld(npoints,ncol) + REAL fluxtop(npoints,ncol) + REAL trans_layers_above(npoints,ncol) + real isccp_taumin,fluxtopinit(npoints),tauir(npoints) + REAL albedocld(npoints,ncol) + real boxarea + integer debug ! set to non-zero value to print out inputs + ! with step debug + integer debugcol ! set to non-zero value to print out column + ! decomposition with step debugcol + integer rangevec(npoints),rangeerror + + integer index1(npoints),num1,jj,k1,k2 + real rec2p13,tauchk,logp,logp1,logp2,atd + + character*10 ftn09 + + DATA isccp_taumin / 0.3 / + DATA cchar / ' ','-','1','+','I','+'/ + DATA cchar_realtops / ' ',' ','1','1','I','I'/ + +! ------ End duplicate definitions common to wrapper routine + + ncolprint=0 + + CALL SCOPS( + & npoints, + & nlev, + & ncol, + & seed, + & cc, + & conv, + & overlap, + & frac_out, + & ncolprint + &) + + CALL ICARUS( + & debug, + & debugcol, + & npoints, + & sunlit, + & nlev, + & ncol, + & pfull, + & phalf, + & qv, + & cc, + & conv, + & dtau_s, + & dtau_c, + & top_height, + & top_height_direction, + & overlap, + & frac_out, + & skt, + & emsfc_lw, + & at, + & dem_s, + & dem_c, + & fq_isccp, + & totalcldarea, + & meanptop, + & meantaucld, + & meanalbedocld, + & meantb, + & meantbclr, + & boxtau, + & boxptop + &) + + return + end + Index: LMDZ5/trunk/libf/phylmd/cosp/modis_simulator.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/cosp/modis_simulator.F90 (revision 2432) +++ LMDZ5/trunk/libf/phylmd/cosp/modis_simulator.F90 (revision 2432) @@ -0,0 +1,1020 @@ +! (c) 2009-2010, Regents of the Unversity of Colorado +! Author: Robert Pincus, Cooperative Institute for Research in the Environmental Sciences +! All rights reserved. +! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $ +! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/MODIS_simulator/modis_simulator.F90 $ +! +! Redistribution and use in source and binary forms, with or without modification, are permitted +! provided that the following conditions are met: +! +! * Redistributions of source code must retain the above copyright notice, this list +! of conditions and the following disclaimer. +! * Redistributions in binary form must reproduce the above copyright notice, this list +! of conditions and the following disclaimer in the documentation and/or other materials +! provided with the distribution. +! * Neither the name of the Met Office nor the names of its contributors may be used +! to endorse or promote products derived from this software without specific prior written +! permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR +! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND +! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR +! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER +! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! + +! +! History: +! May 2009 - Robert Pincus - Initial version +! June 2009 - Steve Platnick and Robert Pincus - Simple radiative transfer for size retrievals +! August 2009 - Robert Pincus - Consistency and bug fixes suggested by Rick Hemler (GFDL) +! November 2009 - Robert Pincus - Bux fixes and speed-ups after experience with Rick Hemler using AM2 (GFDL) +! January 2010 - Robert Pincus - Added high, middle, low cloud fractions +! + +! +! Notes on using the MODIS simulator: +! *) You may provide either layer-by-layer values of optical thickness at 0.67 and 2.1 microns, or +! optical thickness at 0.67 microns and ice- and liquid-water contents (in consistent units of +! your choosing) +! *) Required input also includes the optical thickness and cloud top pressure +! derived from the ISCCP simulator run with parameter top_height = 1. +! *) Cloud particle sizes are specified as radii, measured in meters, though within the module we +! use units of microns. Where particle sizes are outside the bounds used in the MODIS retrieval +! libraries (parameters re_water_min, re_ice_min, etc.) the simulator returns missing values (re_fill) + +! +! When error conditions are encountered this code calls the function complain_and_die, supplied at the +! bottom of this module. Users probably want to replace this with something more graceful. +! +module mod_modis_sim + USE MOD_COSP_TYPES, only: R_UNDEF + implicit none + ! ------------------------------ + ! Algorithmic parameters + ! + + real, parameter :: ice_density = 0.93 ! liquid density is 1. + ! + ! Retrieval parameters + ! + real, parameter :: min_OpticalThickness = 0.3, & ! Minimum detectable optical thickness + CO2Slicing_PressureLimit = 700. * 100., & ! Cloud with higher pressures use thermal methods, units Pa + CO2Slicing_TauLimit = 1., & ! How deep into the cloud does CO2 slicing see? + phase_TauLimit = 1., & ! How deep into the cloud does the phase detection see? + size_TauLimit = 2., & ! Depth of the re retreivals + phaseDiscrimination_Threshold = 0.7 ! What fraction of total extincton needs to be + ! in a single category to make phase discrim. work? + real, parameter :: re_fill= -999. + integer, parameter :: phaseIsNone = 0, phaseIsLiquid = 1, phaseIsIce = 2, phaseIsUndetermined = 3 + + logical, parameter :: useSimpleReScheme = .false. + ! + ! These are the limits of the libraries for the MODIS collection 5 algorithms + ! They are also the limits used in the fits for g and w0 + ! + real, parameter :: re_water_min= 4., re_water_max= 30., re_ice_min= 5., re_ice_max= 90. + integer, parameter :: num_trial_res = 15 ! increase to make the linear pseudo-retrieval of size more accurate + logical, parameter :: use_two_re_iterations = .false. ! do two retrieval iterations? + + ! + ! Precompute near-IR optical params vs size for retrieval scheme + ! + integer, private :: i + real, dimension(num_trial_res), parameter :: & + trial_re_w = re_water_min + (re_water_max - re_water_min)/(num_trial_res-1) * (/ (i - 1, i = 1, num_trial_res) /), & + trial_re_i = re_ice_min + (re_ice_max - re_ice_min) /(num_trial_res-1) * (/ (i - 1, i = 1, num_trial_res) /) + + ! Can't initialze these during compilation, but do in before looping columns in retrievals + real, dimension(num_trial_res) :: g_w, g_i, w0_w, w0_i + ! ------------------------------ + ! Bin boundaries for the joint optical thickness/cloud top pressure histogram + ! + integer, parameter :: numTauHistogramBins = 6, numPressureHistogramBins = 7 + + real, private :: dummy_real + real, dimension(numTauHistogramBins + 1), parameter :: & + tauHistogramBoundaries = (/ min_OpticalThickness, 1.3, 3.6, 9.4, 23., 60., 10000. /) + real, dimension(numPressureHistogramBins + 1), parameter :: & ! Units Pa + pressureHistogramBoundaries = (/ 0., 18000., 31000., 44000., 56000., 68000., 80000., 1000000. /) + real, parameter :: highCloudPressureLimit = 440. * 100., lowCloudPressureLimit = 680. * 100. + ! + ! For output - nominal bin centers and bin boundaries. On output pressure bins are highest to lowest. + ! + integer, private :: k, l + real, parameter, dimension(2, numTauHistogramBins) :: & + nominalTauHistogramBoundaries = & + reshape(source = (/ tauHistogramBoundaries(1), & + ((tauHistogramBoundaries(k), l = 1, 2), k = 2, numTauHistogramBins), & + 100000. /), & + shape = (/2, numTauHistogramBins /) ) + real, parameter, dimension(numTauHistogramBins) :: & + nominalTauHistogramCenters = (nominalTauHistogramBoundaries(1, :) + & + nominalTauHistogramBoundaries(2, :) ) / 2. + + real, parameter, dimension(2, numPressureHistogramBins) :: & + nominalPressureHistogramBoundaries = & + reshape(source = (/ 100000., & + ((pressureHistogramBoundaries(k), l = 1, 2), k = numPressureHistogramBins, 2, -1), & + 0. /), & + shape = (/2, numPressureHistogramBins /) ) + real, parameter, dimension(numPressureHistogramBins) :: & + nominalPressureHistogramCenters = (nominalPressureHistogramBoundaries(1, :) + & + nominalPressureHistogramBoundaries(2, :) ) / 2. + ! ------------------------------ + ! There are two ways to call the MODIS simulator: + ! 1) Provide total optical thickness and liquid/ice water content and we'll partition tau in + ! subroutine modis_L2_simulator_oneTau, or + ! 2) Provide ice and liquid optical depths in each layer + ! + interface modis_L2_simulator + module procedure modis_L2_simulator_oneTau, modis_L2_simulator_twoTaus + end interface +contains + !------------------------------------------------------------------------------------------------ + ! MODIS simulator using specified liquid and ice optical thickness in each layer + ! + ! Note: this simulator operates on all points; to match MODIS itself night-time + ! points should be excluded + ! + ! Note: the simulator requires as input the optical thickness and cloud top pressure + ! derived from the ISCCP simulator run with parameter top_height = 1. + ! If cloud top pressure is higher than about 700 mb, MODIS can't use CO2 slicing + ! and reverts to a thermal algorithm much like ISCCP's. Rather than replicate that + ! alogrithm in this simulator we simply report the values from the ISCCP simulator. + ! + subroutine modis_L2_simulator_twoTaus( & + temp, pressureLayers, pressureLevels, & + liquid_opticalThickness, ice_opticalThickness, & + waterSize, iceSize, & + isccpTau, isccpCloudTopPressure, & + retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize) + + ! Grid-mean quantities at layer centers, starting at the model top + ! dimension nLayers + real, dimension(:), intent(in ) :: temp, & ! Temperature, K + pressureLayers, & ! Pressure, Pa + pressureLevels ! Pressure at layer edges, Pa (dimension nLayers + 1) + ! Sub-column quantities + ! dimension nSubcols, nLayers + real, dimension(:, :), intent(in ) :: liquid_opticalThickness, & ! Layer optical thickness @ 0.67 microns due to liquid + ice_opticalThickness ! ditto, due to ice + real, dimension(:, :), intent(in ) :: waterSize, & ! Cloud drop effective radius, microns + iceSize ! Cloud ice effective radius, microns + + ! Cloud properties retrieved from ISCCP using top_height = 1 + ! dimension nSubcols + real, dimension(:), intent(in ) :: isccpTau, & ! Column-integrated optical thickness + isccpCloudTopPressure ! ISCCP-retrieved cloud top pressure (Pa) + + ! Properties retrieved by MODIS + ! dimension nSubcols + integer, dimension(:), intent(out) :: retrievedPhase ! liquid/ice/other - integer, defined in module header + real, dimension(:), intent(out) :: retrievedCloudTopPressure, & ! units of pressureLayers + retrievedTau, & ! unitless + retrievedSize ! microns + ! --------------------------------------------------- + ! Local variables + logical, dimension(size(retrievedTau)) :: cloudMask + real, dimension(size(waterSize, 1), size(waterSize, 2)) :: tauLiquidFraction, tauTotal + real :: integratedLiquidFraction + integer :: i, nSubcols, nLevels + + ! --------------------------------------------------- + nSubcols = size(liquid_opticalThickness, 1) + nLevels = size(liquid_opticalThickness, 2) + + ! + ! Initial error checks + ! + if(any((/ size(ice_opticalThickness, 1), size(waterSize, 1), size(iceSize, 1), & + size(isccpTau), size(isccpCloudTopPressure), & + size(retrievedPhase), size(retrievedCloudTopPressure), & + size(retrievedTau), size(retrievedSize) /) /= nSubcols )) & + call complain_and_die("Differing number of subcolumns in one or more arrays") + + if(any((/ size(ice_opticalThickness, 2), size(waterSize, 2), size(iceSize, 2), & + size(temp), size(pressureLayers), size(pressureLevels)-1 /) /= nLevels )) & + call complain_and_die("Differing number of levels in one or more arrays") + + + if(any( (/ any(temp <= 0.), any(pressureLayers <= 0.), & + any(liquid_opticalThickness < 0.), & + any(ice_opticalThickness < 0.), & + any(waterSize < 0.), any(iceSize < 0.) /) )) & + call complain_and_die("Input values out of bounds") + + ! --------------------------------------------------- + ! + ! Compute the total optical thickness and the proportion due to liquid in each cell + ! + where(liquid_opticalThickness(:, :) + ice_opticalThickness(:, :) > 0.) + tauLiquidFraction(:, :) = liquid_opticalThickness(:, :)/(liquid_opticalThickness(:, :) + ice_opticalThickness(:, :)) + elsewhere + tauLiquidFraction(:, :) = 0. + end where + tauTotal(:, :) = liquid_opticalThickness(:, :) + ice_opticalThickness(:, :) + + ! + ! Optical depth retrieval + ! This is simply a sum over the optical thickness in each layer + ! It should agree with the ISCCP values after min values have been excluded + ! + retrievedTau(:) = sum(tauTotal(:, :), dim = 2) + + ! + ! Cloud detection - does optical thickness exceed detection threshold? + ! + cloudMask = retrievedTau(:) >= min_OpticalThickness + + ! + ! Initialize initial estimates for size retrievals + ! + if(any(cloudMask) .and. .not. useSimpleReScheme) then + g_w(:) = get_g_nir( phaseIsLiquid, trial_re_w(:)) + w0_w(:) = get_ssa_nir(phaseIsLiquid, trial_re_w(:)) + g_i(:) = get_g_nir( phaseIsIce, trial_re_i(:)) + w0_i(:) = get_ssa_nir(phaseIsIce, trial_re_i(:)) + end if + + do i = 1, nSubCols + if(cloudMask(i)) then + ! + ! Cloud top pressure determination + ! MODIS uses CO2 slicing for clouds with tops above about 700 mb and thermal methods for clouds + ! lower than that. + ! For CO2 slicing we report the optical-depth weighted pressure, integrating to a specified + ! optical depth + ! This assumes linear variation in p between levels. Linear in ln(p) is probably better, + ! though we'd need to deal with the lowest pressure gracefully. + ! + retrievedCloudTopPressure(i) = cloud_top_pressure((/ 0., tauTotal(i, :) /), & + pressureLevels, & + CO2Slicing_TauLimit) + + + ! + ! Phase determination - determine fraction of total tau that's liquid + ! When ice and water contribute about equally to the extinction we can't tell + ! what the phase is + ! + integratedLiquidFraction = weight_by_extinction(tauTotal(i, :), & + tauLiquidFraction(i, :), & + phase_TauLimit) + if(integratedLiquidFraction >= phaseDiscrimination_Threshold) then + retrievedPhase(i) = phaseIsLiquid + else if (integratedLiquidFraction <= 1.- phaseDiscrimination_Threshold) then + retrievedPhase(i) = phaseIsIce + else + retrievedPhase(i) = phaseIsUndetermined + end if + + ! + ! Size determination + ! + if(useSimpleReScheme) then + ! This is the extinction-weighted size considering only the phase we've chosen + ! + if(retrievedPhase(i) == phaseIsIce) then + retrievedSize(i) = weight_by_extinction(ice_opticalThickness(i, :), & + iceSize(i, :), & + (1. - integratedLiquidFraction) * size_TauLimit) + + else if(retrievedPhase(i) == phaseIsLiquid) then + retrievedSize(i) = weight_by_extinction(liquid_opticalThickness(i, :), & + waterSize(i, :), & + integratedLiquidFraction * size_TauLimit) + + else + retrievedSize(i) = 0. + end if + else + retrievedSize(i) = 1.0e-06*retrieve_re(retrievedPhase(i), retrievedTau(i), & + obs_Refl_nir = compute_nir_reflectance(liquid_opticalThickness(i, :), waterSize(i, :)*1.0e6, & + ice_opticalThickness(i, :), iceSize(i, :)*1.0e6)) + end if + else + ! + ! Values when we don't think there's a cloud. + ! + retrievedCloudTopPressure(i) = R_UNDEF + retrievedPhase(i) = phaseIsNone + retrievedSize(i) = R_UNDEF + retrievedTau(i) = R_UNDEF + end if + end do + where((retrievedSize(:) < 0.).and.(retrievedSize(:) /= R_UNDEF)) retrievedSize(:) = 1.0e-06*re_fill + + ! We use the ISCCP-derived CTP for low clouds, since the ISCCP simulator ICARUS + ! mimics what MODIS does to first order. + ! Of course, ISCCP cloud top pressures are in mb. + ! + where(cloudMask(:) .and. retrievedCloudTopPressure(:) > CO2Slicing_PressureLimit) & + retrievedCloudTopPressure(:) = isccpCloudTopPressure * 100. + + end subroutine modis_L2_simulator_twoTaus + !------------------------------------------------------------------------------------------------ + ! + ! MODIS simulator: provide a single optical thickness and the cloud ice and liquid contents; + ! we'll partition this into ice and liquid optical thickness and call the full MODIS simulator + ! + subroutine modis_L2_simulator_oneTau( & + temp, pressureLayers, pressureLevels, & + opticalThickness, cloudWater, cloudIce, & + waterSize, iceSize, & + isccpTau, isccpCloudTopPressure, & + retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize) + ! Grid-mean quantities at layer centers, + ! dimension nLayers + real, dimension(:), intent(in ) :: temp, & ! Temperature, K + pressureLayers, & ! Pressure, Pa + pressureLevels ! Pressure at layer edges, Pa (dimension nLayers + 1) + ! Sub-column quantities + ! dimension nLayers, nSubcols + real, dimension(:, :), intent(in ) :: opticalThickness, & ! Layer optical thickness @ 0.67 microns + cloudWater, & ! Cloud water content, arbitrary units + cloudIce ! Cloud water content, same units as cloudWater + real, dimension(:, :), intent(in ) :: waterSize, & ! Cloud drop effective radius, microns + iceSize ! Cloud ice effective radius, microns + + ! Cloud properties retrieved from ISCCP using top_height = 1 + ! dimension nSubcols + + real, dimension(:), intent(in ) :: isccpTau, & ! Column-integrated optical thickness + isccpCloudTopPressure ! ISCCP-retrieved cloud top pressure (Pa) + + ! Properties retrieved by MODIS + ! dimension nSubcols + integer, dimension(:), intent(out) :: retrievedPhase ! liquid/ice/other - integer + real, dimension(:), intent(out) :: retrievedCloudTopPressure, & ! units of pressureLayers + retrievedTau, & ! unitless + retrievedSize ! microns (or whatever units + ! waterSize and iceSize are supplied in) + ! --------------------------------------------------- + ! Local variables + real, dimension(size(opticalThickness, 1), size(opticalThickness, 2)) :: & + liquid_opticalThickness, ice_opticalThickness, tauLiquidFraction + + ! --------------------------------------------------- + + where(cloudIce(:, :) <= 0.) + tauLiquidFraction(:, :) = 1. + elsewhere + where (cloudWater(:, :) <= 0.) + tauLiquidFraction(:, :) = 0. + elsewhere + ! + ! Geometic optics limit - tau as LWP/re (proportional to LWC/re) + ! + tauLiquidFraction(:, :) = (cloudWater(:, :)/waterSize(:, :)) / & + (cloudWater(:, :)/waterSize(:, :) + cloudIce(:, :)/(ice_density * iceSize(:, :)) ) + end where + end where + liquid_opticalThickness(:, :) = tauLiquidFraction(:, :) * opticalThickness(:, :) + ice_opticalThickness (:, :) = opticalThickness(:, :) - liquid_opticalThickness(:, :) + + call modis_L2_simulator_twoTaus(temp, pressureLayers, pressureLevels, & + liquid_opticalThickness, ice_opticalThickness, & + waterSize, iceSize, & + isccpTau, isccpCloudTopPressure, & + retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize) + + end subroutine modis_L2_simulator_oneTau + !------------------------------------------------------------------------------------------------ + subroutine modis_L3_simulator(phase, cloud_top_pressure, optical_thickness, particle_size, & + Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, & + Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, & + Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, & + Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10, & + Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, & + Cloud_Top_Pressure_Total_Mean, & + Liquid_Water_Path_Mean, Ice_Water_Path_Mean, & + Optical_Thickness_vs_Cloud_Top_Pressure) + ! + ! Inputs; dimension nPoints, nSubcols + ! + integer, dimension(:, :), intent(in) :: phase + real, dimension(:, :), intent(in) :: cloud_top_pressure, optical_thickness, particle_size + ! + ! Outputs; dimension nPoints + ! + real, dimension(:), intent(out) :: & + Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, & + Cloud_Fraction_High_Mean, Cloud_Fraction_Mid_Mean, Cloud_Fraction_Low_Mean, & + Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, & + Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10, & + Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, & + Cloud_Top_Pressure_Total_Mean, & + Liquid_Water_Path_Mean, Ice_Water_Path_Mean + ! tau/ctp histogram; dimensions nPoints, numTauHistogramBins , numPressureHistogramBins + real, dimension(:, :, :), intent(out) :: Optical_Thickness_vs_Cloud_Top_Pressure + ! --------------------------- + ! Local variables + ! + real, parameter :: LWP_conversion = 2./3. * 1000. ! MKS units + integer :: i, j + integer :: nPoints, nSubcols + logical, dimension(size(phase, 1), size(phase, 2)) :: & + cloudMask, waterCloudMask, iceCloudMask, validRetrievalMask + logical, dimension(size(phase, 1), size(phase, 2), numTauHistogramBins ) :: tauMask + logical, dimension(size(phase, 1), size(phase, 2), numPressureHistogramBins) :: pressureMask + ! --------------------------- + + nPoints = size(phase, 1) + nSubcols = size(phase, 2) + ! + ! Array conformance checks + ! + if(any( (/ size(cloud_top_pressure, 1), size(optical_thickness, 1), size(particle_size, 1), & + size(Cloud_Fraction_Total_Mean), size(Cloud_Fraction_Water_Mean), size(Cloud_Fraction_Ice_Mean), & + size(Cloud_Fraction_High_Mean), size(Cloud_Fraction_Mid_Mean), size(Cloud_Fraction_Low_Mean), & + size(Optical_Thickness_Total_Mean), size(Optical_Thickness_Water_Mean), size(Optical_Thickness_Ice_Mean), & + size(Optical_Thickness_Total_MeanLog10), size(Optical_Thickness_Water_MeanLog10), & + size(Optical_Thickness_Ice_MeanLog10), size(Cloud_Particle_Size_Water_Mean), & + size(Cloud_Particle_Size_Ice_Mean), size(Cloud_Top_Pressure_Total_Mean), & + size(Liquid_Water_Path_Mean), size(Ice_Water_Path_Mean) /) /= nPoints)) & + call complain_and_die("Some L3 arrays have wrong number of grid points") + if(any( (/ size(cloud_top_pressure, 2), size(optical_thickness, 2), size(particle_size, 2) /) /= nSubcols)) & + call complain_and_die("Some L3 arrays have wrong number of subcolumns") + + + ! + ! Include only those pixels with successful retrievals in the statistics + ! + validRetrievalMask(:, :) = particle_size(:, :) > 0. + cloudMask = phase(:, :) /= phaseIsNone .and. validRetrievalMask(:, :) + waterCloudMask = phase(:, :) == phaseIsLiquid .and. validRetrievalMask(:, :) + iceCloudMask = phase(:, :) == phaseIsIce .and. validRetrievalMask(:, :) + ! + ! Use these as pixel counts at first + ! + Cloud_Fraction_Total_Mean(:) = real(count(cloudMask, dim = 2)) + Cloud_Fraction_Water_Mean(:) = real(count(waterCloudMask, dim = 2)) + Cloud_Fraction_Ice_Mean(:) = real(count(iceCloudMask, dim = 2)) + + Cloud_Fraction_High_Mean(:) = real(count(cloudMask .and. cloud_top_pressure <= highCloudPressureLimit, dim = 2)) + Cloud_Fraction_Low_Mean(:) = real(count(cloudMask .and. cloud_top_pressure > lowCloudPressureLimit, dim = 2)) + Cloud_Fraction_Mid_Mean(:) = Cloud_Fraction_Total_Mean(:) - Cloud_Fraction_High_Mean(:) - Cloud_Fraction_Low_Mean(:) + + ! + ! Don't want to divide by 0, even though the sums will be 0 where the pixel counts are 0. + ! + where (Cloud_Fraction_Total_Mean == 0) Cloud_Fraction_Total_Mean = -1. + where (Cloud_Fraction_Water_Mean == 0) Cloud_Fraction_Water_Mean = -1. + where (Cloud_Fraction_Ice_Mean == 0) Cloud_Fraction_Ice_Mean = -1. + + Optical_Thickness_Total_Mean = sum(optical_thickness, mask = cloudMask, dim = 2) / Cloud_Fraction_Total_Mean(:) + Optical_Thickness_Water_Mean = sum(optical_thickness, mask = waterCloudMask, dim = 2) / Cloud_Fraction_Water_Mean(:) + Optical_Thickness_Ice_Mean = sum(optical_thickness, mask = iceCloudMask, dim = 2) / Cloud_Fraction_Ice_Mean(:) + + ! We take the absolute value of optical thickness here to satisfy compilers that complains when we + ! evaluate the logarithm of a negative number, even though it's not included in the sum. + Optical_Thickness_Total_MeanLog10 = sum(log10(abs(optical_thickness)), mask = cloudMask, dim = 2) / & + Cloud_Fraction_Total_Mean(:) + Optical_Thickness_Water_MeanLog10 = sum(log10(abs(optical_thickness)), mask = waterCloudMask, dim = 2) / & + Cloud_Fraction_Water_Mean(:) + Optical_Thickness_Ice_MeanLog10 = sum(log10(abs(optical_thickness)), mask = iceCloudMask, dim = 2) / & + Cloud_Fraction_Ice_Mean(:) + + Cloud_Particle_Size_Water_Mean = sum(particle_size, mask = waterCloudMask, dim = 2) / Cloud_Fraction_Water_Mean(:) + Cloud_Particle_Size_Ice_Mean = sum(particle_size, mask = iceCloudMask, dim = 2) / Cloud_Fraction_Ice_Mean(:) + + Cloud_Top_Pressure_Total_Mean = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / max(1, count(cloudMask, dim = 2)) + + Liquid_Water_Path_Mean = LWP_conversion & + * sum(particle_size * optical_thickness, mask = waterCloudMask, dim = 2) & + / Cloud_Fraction_Water_Mean(:) + Ice_Water_Path_Mean = LWP_conversion * ice_density & + * sum(particle_size * optical_thickness, mask = iceCloudMask, dim = 2) & + / Cloud_Fraction_Ice_Mean(:) + + ! + ! Normalize pixel counts to fraction + ! The first three cloud fractions have been set to -1 in cloud-free areas, so set those places to 0. + ! + Cloud_Fraction_Total_Mean(:) = max(0., Cloud_Fraction_Total_Mean(:)/nSubcols) + Cloud_Fraction_Water_Mean(:) = max(0., Cloud_Fraction_Water_Mean(:)/nSubcols) + Cloud_Fraction_Ice_Mean(:) = max(0., Cloud_Fraction_Ice_Mean(:) /nSubcols) + + Cloud_Fraction_High_Mean(:) = Cloud_Fraction_High_Mean(:) /nSubcols + Cloud_Fraction_Mid_Mean(:) = Cloud_Fraction_Mid_Mean(:) /nSubcols + Cloud_Fraction_Low_Mean(:) = Cloud_Fraction_Low_Mean(:) /nSubcols + + ! ---- + ! Joint histogram + ! + do i = 1, numTauHistogramBins + where(cloudMask(:, :)) + tauMask(:, :, i) = optical_thickness(:, :) >= tauHistogramBoundaries(i) .and. & + optical_thickness(:, :) < tauHistogramBoundaries(i+1) + elsewhere + tauMask(:, :, i) = .false. + end where + end do + + do i = 1, numPressureHistogramBins + where(cloudMask(:, :)) + pressureMask(:, :, i) = cloud_top_pressure(:, :) >= pressureHistogramBoundaries(i) .and. & + cloud_top_pressure(:, :) < pressureHistogramBoundaries(i+1) + elsewhere + pressureMask(:, :, i) = .false. + end where + end do + + do i = 1, numPressureHistogramBins + do j = 1, numTauHistogramBins + Optical_Thickness_vs_Cloud_Top_Pressure(:, j, i) = & + real(count(tauMask(:, :, j) .and. pressureMask(:, :, i), dim = 2)) / real(nSubcols) + end do + end do + + end subroutine modis_L3_simulator + !------------------------------------------------------------------------------------------------ + function cloud_top_pressure(tauIncrement, pressure, tauLimit) + real, dimension(:), intent(in) :: tauIncrement, pressure + real, intent(in) :: tauLimit + real :: cloud_top_pressure + ! + ! Find the extinction-weighted pressure. Assume that pressure varies linearly between + ! layers and use the trapezoidal rule. + ! + + real :: deltaX, totalTau, totalProduct + integer :: i + + totalTau = 0.; totalProduct = 0. + do i = 2, size(tauIncrement) + if(totalTau + tauIncrement(i) > tauLimit) then + deltaX = tauLimit - totalTau + totalTau = totalTau + deltaX + ! + ! Result for trapezoidal rule when you take less than a full step + ! tauIncrement is a layer-integrated value + ! + totalProduct = totalProduct & + + pressure(i-1) * deltaX & + + (pressure(i) - pressure(i-1)) * deltaX**2/(2. * tauIncrement(i)) + else + totalTau = totalTau + tauIncrement(i) + totalProduct = totalProduct + tauIncrement(i) * (pressure(i) + pressure(i-1)) / 2. + end if + if(totalTau >= tauLimit) exit + end do + cloud_top_pressure = totalProduct/totalTau + end function cloud_top_pressure + !------------------------------------------------------------------------------------------------ + function weight_by_extinction(tauIncrement, f, tauLimit) + real, dimension(:), intent(in) :: tauIncrement, f + real, intent(in) :: tauLimit + real :: weight_by_extinction + ! + ! Find the extinction-weighted value of f(tau), assuming constant f within each layer + ! + + real :: deltaX, totalTau, totalProduct + integer :: i + + totalTau = 0.; totalProduct = 0. + do i = 1, size(tauIncrement) + if(totalTau + tauIncrement(i) > tauLimit) then + deltaX = tauLimit - totalTau + totalTau = totalTau + deltaX + totalProduct = totalProduct + deltaX * f(i) + else + totalTau = totalTau + tauIncrement(i) + totalProduct = totalProduct + tauIncrement(i) * f(i) + end if + if(totalTau >= tauLimit) exit + end do + weight_by_extinction = totalProduct/totalTau + end function weight_by_extinction + !------------------------------------------------------------------------------------------------ + pure function compute_nir_reflectance(water_tau, water_size, ice_tau, ice_size) + real, dimension(:), intent(in) :: water_tau, water_size, ice_tau, ice_size + real :: compute_nir_reflectance + + real, dimension(size(water_tau)) :: water_g, water_w0, ice_g, ice_w0, & + tau, g, w0 + !---------------------------------------- + water_g(:) = get_g_nir( phaseIsLiquid, water_size) + water_w0(:) = get_ssa_nir(phaseIsLiquid, water_size) + ice_g(:) = get_g_nir( phaseIsIce, ice_size) + ice_w0(:) = get_ssa_nir(phaseIsIce, ice_size) + ! + ! Combine ice and water optical properties + ! + g(:) = 0; w0(:) = 0. + tau(:) = ice_tau(:) + water_tau(:) + where (tau(:) > 0) + g(:) = (water_tau(:) * water_g(:) + ice_tau(:) * ice_g(:) ) / & + tau(:) + w0(:) = (water_tau(:) * water_g(:) * water_w0(:) + ice_tau(:) * ice_g(:) * ice_w0(:)) / & + (g(:) * tau(:)) + end where + + compute_nir_reflectance = compute_toa_reflectace(tau, g, w0) + end function compute_nir_reflectance + !------------------------------------------------------------------------------------------------ + ! Retreivals + !------------------------------------------------------------------------------------------------ + elemental function retrieve_re (phase, tau, obs_Refl_nir) + integer, intent(in) :: phase + real, intent(in) :: tau, obs_Refl_nir + real :: retrieve_re + ! + ! Finds the re that produces the minimum mis-match between predicted and observed reflectance in + ! MODIS band 7 (near IR) + ! Uses + ! fits for asymmetry parameter g(re) and single scattering albedo w0(re) based on MODIS tables + ! two-stream for layer reflectance and transmittance as a function of optical thickness tau, g, and w0 + ! adding-doubling for total reflectance + ! + ! + ! + ! Local variables + ! + real, parameter :: min_distance_to_boundary = 0.01 + real :: re_min, re_max, delta_re + integer :: i + + real, dimension(num_trial_res) :: trial_re, g, w0, predicted_Refl_nir + ! -------------------------- + + if(any(phase == (/ phaseIsLiquid, phaseIsUndetermined, phaseIsIce /))) then + if (phase == phaseIsLiquid .OR. phase == phaseIsUndetermined) then + re_min = re_water_min + re_max = re_water_max + trial_re(:) = trial_re_w + g(:) = g_w(:) + w0(:) = w0_w(:) + else + re_min = re_ice_min + re_max = re_ice_max + trial_re(:) = trial_re_i + g(:) = g_i(:) + w0(:) = w0_i(:) + end if + ! + ! 1st attempt at index: w/coarse re resolution + ! + predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:)) + retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir) + ! + ! If first retrieval works, can try 2nd iteration using greater re resolution + ! + if(use_two_re_iterations .and. retrieve_re > 0.) then + re_min = retrieve_re - delta_re + re_max = retrieve_re + delta_re + delta_re = (re_max - re_min)/real(num_trial_res-1) + + trial_re(:) = re_min + delta_re * (/ (i - 1, i = 1, num_trial_res) /) + g(:) = get_g_nir( phase, trial_re(:)) + w0(:) = get_ssa_nir(phase, trial_re(:)) + predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:)) + retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir) + end if + else + retrieve_re = re_fill + end if + + end function retrieve_re + ! -------------------------------------------- + pure function interpolate_to_min(x, y, yobs) + real, dimension(:), intent(in) :: x, y + real, intent(in) :: yobs + real :: interpolate_to_min + ! + ! Given a set of values of y as y(x), find the value of x that minimizes abs(y - yobs) + ! y must be monotonic in x + ! + real, dimension(size(x)) :: diff + integer :: nPoints, minDiffLoc, lowerBound, upperBound + ! --------------------------------- + nPoints = size(y) + diff(:) = y(:) - yobs + minDiffLoc = minloc(abs(diff), dim = 1) + + if(minDiffLoc == 1) then + lowerBound = minDiffLoc + upperBound = minDiffLoc + 1 + else if(minDiffLoc == nPoints) then + lowerBound = minDiffLoc - 1 + upperBound = minDiffLoc + else + if(diff(minDiffLoc-1) * diff(minDiffLoc) < 0) then + lowerBound = minDiffLoc-1 + upperBound = minDiffLoc + else + lowerBound = minDiffLoc + upperBound = minDiffLoc + 1 + end if + end if + + if(diff(lowerBound) * diff(upperBound) < 0) then + ! + ! Interpolate the root position linearly if we bracket the root + ! + interpolate_to_min = x(upperBound) - & + diff(upperBound) * (x(upperBound) - x(lowerBound)) / (diff(upperBound) - diff(lowerBound)) + else + interpolate_to_min = re_fill + end if + + + end function interpolate_to_min + ! -------------------------------------------- + ! Optical properties + ! -------------------------------------------- + elemental function get_g_nir (phase, re) + ! + ! Polynomial fit for asummetry parameter g in MODIS band 7 (near IR) as a function + ! of size for ice and water + ! Fits from Steve Platnick + ! + + integer, intent(in) :: phase + real, intent(in) :: re + real :: get_g_nir + + real, dimension(3), parameter :: ice_coefficients = (/ 0.7432, 4.5563e-3, -2.8697e-5 /), & + small_water_coefficients = (/ 0.8027, -1.0496e-2, 1.7071e-3 /), & + big_water_coefficients = (/ 0.7931, 5.3087e-3, -7.4995e-5 /) + + ! approx. fits from MODIS Collection 5 LUT scattering calculations + if(phase == phaseIsLiquid) then + if(re < 8.) then + get_g_nir = fit_to_quadratic(re, small_water_coefficients) + if(re < re_water_min) get_g_nir = fit_to_quadratic(re_water_min, small_water_coefficients) + else + get_g_nir = fit_to_quadratic(re, big_water_coefficients) + if(re > re_water_max) get_g_nir = fit_to_quadratic(re_water_max, big_water_coefficients) + end if + else + get_g_nir = fit_to_quadratic(re, ice_coefficients) + if(re < re_ice_min) get_g_nir = fit_to_quadratic(re_ice_min, ice_coefficients) + if(re > re_ice_max) get_g_nir = fit_to_quadratic(re_ice_max, ice_coefficients) + end if + + end function get_g_nir + + ! -------------------------------------------- + elemental function get_ssa_nir (phase, re) + integer, intent(in) :: phase + real, intent(in) :: re + real :: get_ssa_nir + ! + ! Polynomial fit for single scattering albedo in MODIS band 7 (near IR) as a function + ! of size for ice and water + ! Fits from Steve Platnick + ! + + real, dimension(4), parameter :: ice_coefficients = (/ 0.9994, -4.5199e-3, 3.9370e-5, -1.5235e-7 /) + real, dimension(3), parameter :: water_coefficients = (/ 1.0008, -2.5626e-3, 1.6024e-5 /) + + ! approx. fits from MODIS Collection 5 LUT scattering calculations + if(phase == phaseIsLiquid) then + get_ssa_nir = fit_to_quadratic(re, water_coefficients) + if(re < re_water_min) get_ssa_nir = fit_to_quadratic(re_water_min, water_coefficients) + if(re > re_water_max) get_ssa_nir = fit_to_quadratic(re_water_max, water_coefficients) + else + get_ssa_nir = fit_to_cubic(re, ice_coefficients) + if(re < re_ice_min) get_ssa_nir = fit_to_cubic(re_ice_min, ice_coefficients) + if(re > re_ice_max) get_ssa_nir = fit_to_cubic(re_ice_max, ice_coefficients) + end if + + end function get_ssa_nir + ! -------------------------------------------- + pure function fit_to_cubic(x, coefficients) + real, intent(in) :: x + real, dimension(:), intent(in) :: coefficients + real :: fit_to_cubic + + + fit_to_cubic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3) + x * coefficients(4))) + end function fit_to_cubic + ! -------------------------------------------- + pure function fit_to_quadratic(x, coefficients) + real, intent(in) :: x + real, dimension(:), intent(in) :: coefficients + real :: fit_to_quadratic + + + fit_to_quadratic = coefficients(1) + x * (coefficients(2) + x * (coefficients(3))) + end function fit_to_quadratic + ! -------------------------------------------- + ! Radiative transfer + ! -------------------------------------------- + pure function compute_toa_reflectace(tau, g, w0) + real, dimension(:), intent(in) :: tau, g, w0 + real :: compute_toa_reflectace + + logical, dimension(size(tau)) :: cloudMask + integer, dimension(count(tau(:) > 0)) :: cloudIndicies + real, dimension(count(tau(:) > 0)) :: Refl, Trans + real :: Refl_tot, Trans_tot + integer :: i + ! --------------------------------------- + ! + ! This wrapper reports reflectance only and strips out non-cloudy elements from the calculation + ! + cloudMask = tau(:) > 0. + cloudIndicies = pack((/ (i, i = 1, size(tau)) /), mask = cloudMask) + do i = 1, size(cloudIndicies) + call two_stream(tau(cloudIndicies(i)), g(cloudIndicies(i)), w0(cloudIndicies(i)), Refl(i), Trans(i)) + end do + + call adding_doubling(Refl(:), Trans(:), Refl_tot, Trans_tot) + + compute_toa_reflectace = Refl_tot + + end function compute_toa_reflectace + ! -------------------------------------------- + pure subroutine two_stream(tauint, gint, w0int, ref, tra) + real, intent(in) :: tauint, gint, w0int + real, intent(out) :: ref, tra + ! + ! Compute reflectance in a single layer using the two stream approximation + ! The code itself is from Lazaros Oreopoulos via Steve Platnick + ! + ! ------------------------ + ! Local variables + ! for delta Eddington code + ! xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1) + integer, parameter :: beam = 2 + real, parameter :: xmu = 0.866, minConservativeW0 = 0.9999999 + real :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, & + rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den, th + ! + ! Compute reflectance and transmittance in a single layer using the two stream approximation + ! The code itself is from Lazaros Oreopoulos via Steve Platnick + ! + f = gint**2 + tau = (1 - w0int * f) * tauint + w0 = (1 - f) * w0int / (1 - w0int * f) + g = (gint - f) / (1 - f) + + ! delta-Eddington (Joseph et al. 1976) + gamma1 = (7 - w0* (4 + 3 * g)) / 4.0 + gamma2 = -(1 - w0* (4 - 3 * g)) / 4.0 + gamma3 = (2 - 3*g*xmu) / 4.0 + gamma4 = 1 - gamma3 + + if (w0int > minConservativeW0) then + ! Conservative scattering + if (beam == 1) then + rh = (gamma1*tau+(gamma3-gamma1*xmu)*(1-exp(-tau/xmu))) + ref = rh / (1 + gamma1 * tau) + tra = 1 - ref + else if(beam == 2) then + ref = gamma1*tau/(1 + gamma1*tau) + tra = 1 - ref + endif + else + ! Non-conservative scattering + a1 = gamma1 * gamma4 + gamma2 * gamma3 + a2 = gamma1 * gamma3 + gamma2 * gamma4 + + rk = sqrt(gamma1**2 - gamma2**2) + + r1 = (1 - rk * xmu) * (a2 + rk * gamma3) + r2 = (1 + rk * xmu) * (a2 - rk * gamma3) + r3 = 2 * rk *(gamma3 - a2 * xmu) + r4 = (1 - (rk * xmu)**2) * (rk + gamma1) + r5 = (1 - (rk * xmu)**2) * (rk - gamma1) + + t1 = (1 + rk * xmu) * (a1 + rk * gamma4) + t2 = (1 - rk * xmu) * (a1 - rk * gamma4) + t3 = 2 * rk * (gamma4 + a1 * xmu) + t4 = r4 + t5 = r5 + + beta = -r5 / r4 + + e1 = min(rk * tau, 500.) + e2 = min(tau / xmu, 500.) + + if (beam == 1) then + den = r4 * exp(e1) + r5 * exp(-e1) + ref = w0*(r1*exp(e1)-r2*exp(-e1)-r3*exp(-e2))/den + den = t4 * exp(e1) + t5 * exp(-e1) + th = exp(-e2) + tra = th-th*w0*(t1*exp(e1)-t2*exp(-e1)-t3*exp(e2))/den + elseif (beam == 2) then + ef1 = exp(-e1) + ef2 = exp(-2*e1) + ref = (gamma2*(1-ef2))/((rk+gamma1)*(1-beta*ef2)) + tra = (2*rk*ef1)/((rk+gamma1)*(1-beta*ef2)) + endif + end if + end subroutine two_stream + ! -------------------------------------------------- + elemental function two_stream_reflectance(tauint, gint, w0int) + real, intent(in) :: tauint, gint, w0int + real :: two_stream_reflectance + ! + ! Compute reflectance in a single layer using the two stream approximation + ! The code itself is from Lazaros Oreopoulos via Steve Platnick + ! + ! ------------------------ + ! Local variables + ! for delta Eddington code + ! xmu, gamma3, and gamma4 only used for collimated beam approximation (i.e., beam=1) + integer, parameter :: beam = 2 + real, parameter :: xmu = 0.866, minConservativeW0 = 0.9999999 + real :: tau, w0, g, f, gamma1, gamma2, gamma3, gamma4, & + rh, a1, a2, rk, r1, r2, r3, r4, r5, t1, t2, t3, t4, t5, beta, e1, e2, ef1, ef2, den + ! ------------------------ + + + f = gint**2 + tau = (1 - w0int * f) * tauint + w0 = (1 - f) * w0int / (1 - w0int * f) + g = (gint - f) / (1 - f) + + ! delta-Eddington (Joseph et al. 1976) + gamma1 = (7 - w0* (4 + 3 * g)) / 4.0 + gamma2 = -(1 - w0* (4 - 3 * g)) / 4.0 + gamma3 = (2 - 3*g*xmu) / 4.0 + gamma4 = 1 - gamma3 + + if (w0int > minConservativeW0) then + ! Conservative scattering + if (beam == 1) then + rh = (gamma1*tau+(gamma3-gamma1*xmu)*(1-exp(-tau/xmu))) + two_stream_reflectance = rh / (1 + gamma1 * tau) + elseif (beam == 2) then + two_stream_reflectance = gamma1*tau/(1 + gamma1*tau) + endif + + else ! + + ! Non-conservative scattering + a1 = gamma1 * gamma4 + gamma2 * gamma3 + a2 = gamma1 * gamma3 + gamma2 * gamma4 + + rk = sqrt(gamma1**2 - gamma2**2) + + r1 = (1 - rk * xmu) * (a2 + rk * gamma3) + r2 = (1 + rk * xmu) * (a2 - rk * gamma3) + r3 = 2 * rk *(gamma3 - a2 * xmu) + r4 = (1 - (rk * xmu)**2) * (rk + gamma1) + r5 = (1 - (rk * xmu)**2) * (rk - gamma1) + + t1 = (1 + rk * xmu) * (a1 + rk * gamma4) + t2 = (1 - rk * xmu) * (a1 - rk * gamma4) + t3 = 2 * rk * (gamma4 + a1 * xmu) + t4 = r4 + t5 = r5 + + beta = -r5 / r4 + + e1 = min(rk * tau, 500.) + e2 = min(tau / xmu, 500.) + + if (beam == 1) then + den = r4 * exp(e1) + r5 * exp(-e1) + two_stream_reflectance = w0*(r1*exp(e1)-r2*exp(-e1)-r3*exp(-e2))/den + elseif (beam == 2) then + ef1 = exp(-e1) + ef2 = exp(-2*e1) + two_stream_reflectance = (gamma2*(1-ef2))/((rk+gamma1)*(1-beta*ef2)) + endif + + end if + end function two_stream_reflectance + ! -------------------------------------------- + pure subroutine adding_doubling (Refl, Tran, Refl_tot, Tran_tot) + real, dimension(:), intent(in) :: Refl, Tran + real, intent(out) :: Refl_tot, Tran_tot + ! + ! Use adding/doubling formulas to compute total reflectance and transmittance from layer values + ! + + integer :: i + real, dimension(size(Refl)) :: Refl_cumulative, Tran_cumulative + + Refl_cumulative(1) = Refl(1); Tran_cumulative(1) = Tran(1) + + do i=2, size(Refl) + ! place (add) previous combined layer(s) reflectance on top of layer i, w/black surface (or ignoring surface): + Refl_cumulative(i) = Refl_cumulative(i-1) + Refl(i)*(Tran_cumulative(i-1)**2)/(1 - Refl_cumulative(i-1) * Refl(i)) + Tran_cumulative(i) = (Tran_cumulative(i-1)*Tran(i)) / (1 - Refl_cumulative(i-1) * Refl(i)) + end do + + Refl_tot = Refl_cumulative(size(Refl)) + Tran_tot = Tran_cumulative(size(Refl)) + + end subroutine adding_doubling + ! -------------------------------------------------- + subroutine complain_and_die(message) + character(len = *), intent(in) :: message + + write(6, *) "Failure in MODIS simulator" + write(6, *) trim(message) + stop + end subroutine complain_and_die + !------------------------------------------------------------------------------------------------ +end module mod_modis_sim Index: LMDZ5/trunk/libf/phylmd/cosp/mrgrnk.F90.new =================================================================== --- LMDZ5/trunk/libf/phylmd/cosp/mrgrnk.F90.new (revision 2432) +++ LMDZ5/trunk/libf/phylmd/cosp/mrgrnk.F90.new (revision 2432) @@ -0,0 +1,611 @@ +! $Revision: 23 $, $Date: 2011-03-31 15:41:37 +0200 (jeu. 31 mars 2011) $ +! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/quickbeam/mrgrnk.f90 $ +Module m_mrgrnk +Integer, Parameter :: kdp = selected_real_kind(15) +public :: mrgrnk +private :: kdp +private :: R_mrgrnk, I_mrgrnk, D_mrgrnk +interface mrgrnk + module procedure D_mrgrnk, R_mrgrnk, I_mrgrnk +end interface mrgrnk +contains + +Subroutine D_mrgrnk (XDONT, IRNGT) +! __________________________________________________________ +! MRGRNK = Merge-sort ranking of an array +! For performance reasons, the first 2 passes are taken +! out of the standard loop, and use dedicated coding. +! __________________________________________________________ +! __________________________________________________________ + Real (kind=kdp), Dimension (:), Intent (In) :: XDONT + Integer*4, Dimension (:), Intent (Out) :: IRNGT +! __________________________________________________________ + Real (kind=kdp) :: XVALA, XVALB +! + Integer, Dimension (SIZE(IRNGT)) :: JWRKT + Integer :: LMTNA, LMTNC, IRNG1, IRNG2 + Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB +! + NVAL = Min (SIZE(XDONT), SIZE(IRNGT)) + Select Case (NVAL) + Case (:0) + Return + Case (1) + IRNGT (1) = 1 + Return + Case Default + Continue + End Select +! +! Fill-in the index array, creating ordered couples +! + Do IIND = 2, NVAL, 2 + If (XDONT(IIND-1) <= XDONT(IIND)) Then + IRNGT (IIND-1) = IIND - 1 + IRNGT (IIND) = IIND + Else + IRNGT (IIND-1) = IIND + IRNGT (IIND) = IIND - 1 + End If + End Do + If (Modulo(NVAL, 2) /= 0) Then + IRNGT (NVAL) = NVAL + End If +! +! We will now have ordered subsets A - B - A - B - ... +! and merge A and B couples into C - C - ... +! + LMTNA = 2 + LMTNC = 4 +! +! First iteration. The length of the ordered subsets goes from 2 to 4 +! + Do + If (NVAL <= 2) Exit +! +! Loop on merges of A and B into C +! + Do IWRKD = 0, NVAL - 1, 4 + If ((IWRKD+4) > NVAL) Then + If ((IWRKD+2) >= NVAL) Exit +! +! 1 2 3 +! + If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit +! +! 1 3 2 +! + If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then + IRNG2 = IRNGT (IWRKD+2) + IRNGT (IWRKD+2) = IRNGT (IWRKD+3) + IRNGT (IWRKD+3) = IRNG2 +! +! 3 1 2 +! + Else + IRNG1 = IRNGT (IWRKD+1) + IRNGT (IWRKD+1) = IRNGT (IWRKD+3) + IRNGT (IWRKD+3) = IRNGT (IWRKD+2) + IRNGT (IWRKD+2) = IRNG1 + End If + Exit + End If +! +! 1 2 3 4 +! + If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle +! +! 1 3 x x +! + If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then + IRNG2 = IRNGT (IWRKD+2) + IRNGT (IWRKD+2) = IRNGT (IWRKD+3) + If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then +! 1 3 2 4 + IRNGT (IWRKD+3) = IRNG2 + Else +! 1 3 4 2 + IRNGT (IWRKD+3) = IRNGT (IWRKD+4) + IRNGT (IWRKD+4) = IRNG2 + End If +! +! 3 x x x +! + Else + IRNG1 = IRNGT (IWRKD+1) + IRNG2 = IRNGT (IWRKD+2) + IRNGT (IWRKD+1) = IRNGT (IWRKD+3) + If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then + IRNGT (IWRKD+2) = IRNG1 + If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then +! 3 1 2 4 + IRNGT (IWRKD+3) = IRNG2 + Else +! 3 1 4 2 + IRNGT (IWRKD+3) = IRNGT (IWRKD+4) + IRNGT (IWRKD+4) = IRNG2 + End If + Else +! 3 4 1 2 + IRNGT (IWRKD+2) = IRNGT (IWRKD+4) + IRNGT (IWRKD+3) = IRNG1 + IRNGT (IWRKD+4) = IRNG2 + End If + End If + End Do +! +! The Cs become As and Bs +! + LMTNA = 4 + Exit + End Do +! +! Iteration loop. Each time, the length of the ordered subsets +! is doubled. +! + Do + If (LMTNA >= NVAL) Exit + IWRKF = 0 + LMTNC = 2 * LMTNC +! +! Loop on merges of A and B into C +! + Do + IWRK = IWRKF + IWRKD = IWRKF + 1 + JINDA = IWRKF + LMTNA + IWRKF = IWRKF + LMTNC + If (IWRKF >= NVAL) Then + If (JINDA >= NVAL) Exit + IWRKF = NVAL + End If + IINDA = 1 + IINDB = JINDA + 1 +! +! Shortcut for the case when the max of A is smaller +! than the min of B. This line may be activated when the +! initial set is already close to sorted. +! +! IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE +! +! One steps in the C subset, that we build in the final rank array +! +! Make a copy of the rank array for the merge iteration +! + JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA) +! + XVALA = XDONT (JWRKT(IINDA)) + XVALB = XDONT (IRNGT(IINDB)) +! + Do + IWRK = IWRK + 1 +! +! We still have unprocessed values in both A and B +! + If (XVALA > XVALB) Then + IRNGT (IWRK) = IRNGT (IINDB) + IINDB = IINDB + 1 + If (IINDB > IWRKF) Then +! Only A still with unprocessed values + IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA) + Exit + End If + XVALB = XDONT (IRNGT(IINDB)) + Else + IRNGT (IWRK) = JWRKT (IINDA) + IINDA = IINDA + 1 + If (IINDA > LMTNA) Exit! Only B still with unprocessed values + XVALA = XDONT (JWRKT(IINDA)) + End If +! + End Do + End Do +! +! The Cs become As and Bs +! + LMTNA = 2 * LMTNA + End Do +! + Return +! +End Subroutine D_mrgrnk + +Subroutine R_mrgrnk (XDONT, IRNGT) +! __________________________________________________________ +! MRGRNK = Merge-sort ranking of an array +! For performance reasons, the first 2 passes are taken +! out of the standard loop, and use dedicated coding. +! __________________________________________________________ +! _________________________________________________________ + Real, Dimension (:), Intent (In) :: XDONT + Integer*4, Dimension (:), Intent (Out) :: IRNGT +! __________________________________________________________ + Real :: XVALA, XVALB +! + Integer, Dimension (SIZE(IRNGT)) :: JWRKT + Integer :: LMTNA, LMTNC, IRNG1, IRNG2 + Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB +! + NVAL = Min (SIZE(XDONT), SIZE(IRNGT)) + Select Case (NVAL) + Case (:0) + Return + Case (1) + IRNGT (1) = 1 + Return + Case Default + Continue + End Select +! +! Fill-in the index array, creating ordered couples +! + Do IIND = 2, NVAL, 2 + If (XDONT(IIND-1) <= XDONT(IIND)) Then + IRNGT (IIND-1) = IIND - 1 + IRNGT (IIND) = IIND + Else + IRNGT (IIND-1) = IIND + IRNGT (IIND) = IIND - 1 + End If + End Do + If (Modulo(NVAL, 2) /= 0) Then + IRNGT (NVAL) = NVAL + End If +! +! We will now have ordered subsets A - B - A - B - ... +! and merge A and B couples into C - C - ... +! + LMTNA = 2 + LMTNC = 4 +! +! First iteration. The length of the ordered subsets goes from 2 to 4 +! + Do + If (NVAL <= 2) Exit +! +! Loop on merges of A and B into C +! + Do IWRKD = 0, NVAL - 1, 4 + If ((IWRKD+4) > NVAL) Then + If ((IWRKD+2) >= NVAL) Exit +! +! 1 2 3 +! + If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit +! +! 1 3 2 +! + If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then + IRNG2 = IRNGT (IWRKD+2) + IRNGT (IWRKD+2) = IRNGT (IWRKD+3) + IRNGT (IWRKD+3) = IRNG2 +! +! 3 1 2 +! + Else + IRNG1 = IRNGT (IWRKD+1) + IRNGT (IWRKD+1) = IRNGT (IWRKD+3) + IRNGT (IWRKD+3) = IRNGT (IWRKD+2) + IRNGT (IWRKD+2) = IRNG1 + End If + Exit + End If +! +! 1 2 3 4 +! + If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle +! +! 1 3 x x +! + If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then + IRNG2 = IRNGT (IWRKD+2) + IRNGT (IWRKD+2) = IRNGT (IWRKD+3) + If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then +! 1 3 2 4 + IRNGT (IWRKD+3) = IRNG2 + Else +! 1 3 4 2 + IRNGT (IWRKD+3) = IRNGT (IWRKD+4) + IRNGT (IWRKD+4) = IRNG2 + End If +! +! 3 x x x +! + Else + IRNG1 = IRNGT (IWRKD+1) + IRNG2 = IRNGT (IWRKD+2) + IRNGT (IWRKD+1) = IRNGT (IWRKD+3) + If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then + IRNGT (IWRKD+2) = IRNG1 + If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then +! 3 1 2 4 + IRNGT (IWRKD+3) = IRNG2 + Else +! 3 1 4 2 + IRNGT (IWRKD+3) = IRNGT (IWRKD+4) + IRNGT (IWRKD+4) = IRNG2 + End If + Else +! 3 4 1 2 + IRNGT (IWRKD+2) = IRNGT (IWRKD+4) + IRNGT (IWRKD+3) = IRNG1 + IRNGT (IWRKD+4) = IRNG2 + End If + End If + End Do +! +! The Cs become As and Bs +! + LMTNA = 4 + Exit + End Do +! +! Iteration loop. Each time, the length of the ordered subsets +! is doubled. +! + Do + If (LMTNA >= NVAL) Exit + IWRKF = 0 + LMTNC = 2 * LMTNC +! +! Loop on merges of A and B into C +! + Do + IWRK = IWRKF + IWRKD = IWRKF + 1 + JINDA = IWRKF + LMTNA + IWRKF = IWRKF + LMTNC + If (IWRKF >= NVAL) Then + If (JINDA >= NVAL) Exit + IWRKF = NVAL + End If + IINDA = 1 + IINDB = JINDA + 1 +! +! Shortcut for the case when the max of A is smaller +! than the min of B. This line may be activated when the +! initial set is already close to sorted. +! +! IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE +! +! One steps in the C subset, that we build in the final rank array +! +! Make a copy of the rank array for the merge iteration +! + JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA) +! + XVALA = XDONT (JWRKT(IINDA)) + XVALB = XDONT (IRNGT(IINDB)) +! + Do + IWRK = IWRK + 1 +! +! We still have unprocessed values in both A and B +! + If (XVALA > XVALB) Then + IRNGT (IWRK) = IRNGT (IINDB) + IINDB = IINDB + 1 + If (IINDB > IWRKF) Then +! Only A still with unprocessed values + IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA) + Exit + End If + XVALB = XDONT (IRNGT(IINDB)) + Else + IRNGT (IWRK) = JWRKT (IINDA) + IINDA = IINDA + 1 + If (IINDA > LMTNA) Exit! Only B still with unprocessed values + XVALA = XDONT (JWRKT(IINDA)) + End If +! + End Do + End Do +! +! The Cs become As and Bs +! + LMTNA = 2 * LMTNA + End Do +! + Return +! +End Subroutine R_mrgrnk +Subroutine I_mrgrnk (XDONT, IRNGT) +! __________________________________________________________ +! MRGRNK = Merge-sort ranking of an array +! For performance reasons, the first 2 passes are taken +! out of the standard loop, and use dedicated coding. +! __________________________________________________________ +! __________________________________________________________ + Integer, Dimension (:), Intent (In) :: XDONT + Integer*4, Dimension (:), Intent (Out) :: IRNGT +! __________________________________________________________ + Integer :: XVALA, XVALB +! + Integer, Dimension (SIZE(IRNGT)) :: JWRKT + Integer :: LMTNA, LMTNC, IRNG1, IRNG2 + Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB +! + NVAL = Min (SIZE(XDONT), SIZE(IRNGT)) + Select Case (NVAL) + Case (:0) + Return + Case (1) + IRNGT (1) = 1 + Return + Case Default + Continue + End Select +! +! Fill-in the index array, creating ordered couples +! + Do IIND = 2, NVAL, 2 + If (XDONT(IIND-1) <= XDONT(IIND)) Then + IRNGT (IIND-1) = IIND - 1 + IRNGT (IIND) = IIND + Else + IRNGT (IIND-1) = IIND + IRNGT (IIND) = IIND - 1 + End If + End Do + If (Modulo(NVAL, 2) /= 0) Then + IRNGT (NVAL) = NVAL + End If +! +! We will now have ordered subsets A - B - A - B - ... +! and merge A and B couples into C - C - ... +! + LMTNA = 2 + LMTNC = 4 +! +! First iteration. The length of the ordered subsets goes from 2 to 4 +! + Do + If (NVAL <= 2) Exit +! +! Loop on merges of A and B into C +! + Do IWRKD = 0, NVAL - 1, 4 + If ((IWRKD+4) > NVAL) Then + If ((IWRKD+2) >= NVAL) Exit +! +! 1 2 3 +! + If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit +! +! 1 3 2 +! + If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then + IRNG2 = IRNGT (IWRKD+2) + IRNGT (IWRKD+2) = IRNGT (IWRKD+3) + IRNGT (IWRKD+3) = IRNG2 +! +! 3 1 2 +! + Else + IRNG1 = IRNGT (IWRKD+1) + IRNGT (IWRKD+1) = IRNGT (IWRKD+3) + IRNGT (IWRKD+3) = IRNGT (IWRKD+2) + IRNGT (IWRKD+2) = IRNG1 + End If + Exit + End If +! +! 1 2 3 4 +! + If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle +! +! 1 3 x x +! + If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then + IRNG2 = IRNGT (IWRKD+2) + IRNGT (IWRKD+2) = IRNGT (IWRKD+3) + If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then +! 1 3 2 4 + IRNGT (IWRKD+3) = IRNG2 + Else +! 1 3 4 2 + IRNGT (IWRKD+3) = IRNGT (IWRKD+4) + IRNGT (IWRKD+4) = IRNG2 + End If +! +! 3 x x x +! + Else + IRNG1 = IRNGT (IWRKD+1) + IRNG2 = IRNGT (IWRKD+2) + IRNGT (IWRKD+1) = IRNGT (IWRKD+3) + If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then + IRNGT (IWRKD+2) = IRNG1 + If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then +! 3 1 2 4 + IRNGT (IWRKD+3) = IRNG2 + Else +! 3 1 4 2 + IRNGT (IWRKD+3) = IRNGT (IWRKD+4) + IRNGT (IWRKD+4) = IRNG2 + End If + Else +! 3 4 1 2 + IRNGT (IWRKD+2) = IRNGT (IWRKD+4) + IRNGT (IWRKD+3) = IRNG1 + IRNGT (IWRKD+4) = IRNG2 + End If + End If + End Do +! +! The Cs become As and Bs +! + LMTNA = 4 + Exit + End Do +! +! Iteration loop. Each time, the length of the ordered subsets +! is doubled. +! + Do + If (LMTNA >= NVAL) Exit + IWRKF = 0 + LMTNC = 2 * LMTNC +! +! Loop on merges of A and B into C +! + Do + IWRK = IWRKF + IWRKD = IWRKF + 1 + JINDA = IWRKF + LMTNA + IWRKF = IWRKF + LMTNC + If (IWRKF >= NVAL) Then + If (JINDA >= NVAL) Exit + IWRKF = NVAL + End If + IINDA = 1 + IINDB = JINDA + 1 +! +! Shortcut for the case when the max of A is smaller +! than the min of B. This line may be activated when the +! initial set is already close to sorted. +! +! IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE +! +! One steps in the C subset, that we build in the final rank array +! +! Make a copy of the rank array for the merge iteration +! + JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA) +! + XVALA = XDONT (JWRKT(IINDA)) + XVALB = XDONT (IRNGT(IINDB)) +! + Do + IWRK = IWRK + 1 +! +! We still have unprocessed values in both A and B +! + If (XVALA > XVALB) Then + IRNGT (IWRK) = IRNGT (IINDB) + IINDB = IINDB + 1 + If (IINDB > IWRKF) Then +! Only A still with unprocessed values + IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA) + Exit + End If + XVALB = XDONT (IRNGT(IINDB)) + Else + IRNGT (IWRK) = JWRKT (IINDA) + IINDA = IINDA + 1 + If (IINDA > LMTNA) Exit! Only B still with unprocessed values + XVALA = XDONT (JWRKT(IINDA)) + End If +! + End Do + End Do +! +! The Cs become As and Bs +! + LMTNA = 2 * LMTNA + End Do +! + Return +! +End Subroutine I_mrgrnk +end module m_mrgrnk Index: LMDZ5/trunk/libf/phylmd/cosp/predict_mom07.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/cosp/predict_mom07.F90 (revision 2432) +++ LMDZ5/trunk/libf/phylmd/cosp/predict_mom07.F90 (revision 2432) @@ -0,0 +1,41 @@ +! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $ +! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/predict_mom07.F90 $ + subroutine predict_mom07(m2,tc,n,m) + + ! subroutine to predict nth moment m given m2, tc, n + ! if m2=-9999 then the routine will predict m2 given m,n,tc + implicit none + + real*8 :: a1,a2,a3,b1,b2,b3,c1,c2,c3 + real*8 :: m2,tc,n,m,a_,b_,c_,A,B,C,n2 + + a1= 13.6078 + a2= -7.76062 + a3= 0.478694 + b1= -0.0360722 + b2= 0.0150830 + b3= 0.00149453 + c1= 0.806856 + c2= 0.00581022 + c3= 0.0456723 + + n2=n*n + a_=a1+a2*n+a3*n2 + b_=b1+b2*n+b3*n2 + c_=c1+c2*n+c3*n2 + + A=exp(a_) + B=b_ + C=c_ + + ! predict m from m2 and tc + if(m2.ne.-9999) then + m=A*exp(B*tc)*m2**C + endif + ! get m2 if mass-dimension relationship not proportional to D**2 + if(m2.eq.-9999) then + m2=(m/(A*exp(B*tc)))**(1.0/C) + endif + + return + end subroutine predict_mom07 Index: LMDZ5/trunk/libf/phylmd/cosp/radar_simulator_init.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/cosp/radar_simulator_init.F90 (revision 2432) +++ LMDZ5/trunk/libf/phylmd/cosp/radar_simulator_init.F90 (revision 2432) @@ -0,0 +1,156 @@ + subroutine radar_simulator_init(freq,k2,use_gas_abs,do_ray,undef, & + nhclass, & + hclass_type,hclass_phase, & + hclass_dmin,hclass_dmax, & + hclass_apm,hclass_bpm,hclass_rho, & + hclass_p1,hclass_p2,hclass_p3, & + load_scale_LUTs_flag,update_scale_LUTs_flag,LUT_file_name, & + hp & ! output + ) + use radar_simulator_types + implicit none + +! Purpose: +! +! Initialize variables used by the radar simulator. +! Part of QuickBeam v3.0 by John Haynes and Roj Marchand +! +! +! Inputs: +! [] from data in hydrometeor class input +! +! [freq] radar frequency (GHz) +! +! [k2] |K|^2, the dielectric constant, set to -1 to use the +! frequency dependent default +! +! [use_gas_abs] 1=do gaseous abs calcs, 0=no gasesous absorbtion calculated, +! 2=calculate (and use) absorption for first profile on all profiles +! +! [undef] mising data value +! [nhclass] number of hydrometeor types +! +! For each hydrometero type: +! hclass_type Type of distribution (see quickbeam documentation) +! hclass_phase 1==ice, 0=liquid +! +! hclass_dmin minimum diameter allowed is drop size distribution N(DDmax)=0 +! +! hclass_apm,hclass_bpm Density of partical apm*D^bpm or constant = rho +! hclass_rho, +! +! hclass_p1,hclass_p2, Default values of DSD parameters (see quickbeam documentation) +! hclass_p3, +! +! load_scale_LUTs_flag Flag, load scale factors Look Up Table from file at start of run +! update_scale_LUTs_flag Flag, save new scale factors calculated during this run to LUT +! LUT_file_name Name of file containing LUT +! +! Outputs: +! [hp] structure that define hydrometeor types +! +! Modified: +! 08/23/2006 placed into subroutine form (Roger Marchand) +! June 2010 New interface to support "radar_simulator_params" structure + +! ----- INPUT ----- + + real, intent(in) :: freq,k2 + integer, intent(in) :: nhclass ! number of hydrometeor classes in use + integer, intent(in) :: use_gas_abs,do_ray + real :: undef + real,dimension(nhclass),intent(in) :: hclass_dmin,hclass_dmax, & + hclass_apm,hclass_bpm,hclass_rho, & + hclass_p1,hclass_p2,hclass_p3 + integer,dimension(nhclass),intent(in) :: hclass_type,hclass_phase + + logical, intent(in) :: load_scale_LUTs_flag,update_scale_LUTs_flag + character*240, intent(in) :: LUT_file_name + +! ----- OUTPUTS ----- + type(class_param), intent(out) :: hp + +! ----- INTERNAL ----- + integer :: i,j + real*8 :: delt, deltp + + ! + ! set radar simulation properites + ! + hp%freq=freq + hp%k2=k2 + hp%use_gas_abs=use_gas_abs + hp%do_ray=do_ray + hp%nhclass=nhclass + + hp%load_scale_LUTs=load_scale_LUTs_flag + hp%update_scale_LUTs=update_scale_LUTs_flag + hp%scale_LUT_file_name=LUT_file_name + + ! + ! Store settings for hydrometeor types in hp (class_parameter) structure. + ! + do i = 1,nhclass + hp%dtype(i) = hclass_type(i) + hp%phase(i) = hclass_phase(i) + hp%dmin(i) = hclass_dmin(i) + hp%dmax(i) = hclass_dmax(i) + hp%apm(i) = hclass_apm(i) + hp%bpm(i) = hclass_bpm(i) + hp%rho(i) = hclass_rho(i) + hp%p1(i) = hclass_p1(i) + hp%p2(i) = hclass_p2(i) + hp%p3(i) = hclass_p3(i) + enddo + + ! + ! initialize scaling array + ! + hp%N_scale_flag = .false. + hp%fc = undef + hp%rho_eff = undef + + hp%Z_scale_flag = .false. + hp%Ze_scaled = 0.0 + hp%Zr_scaled = 0.0 + hp%kr_scaled = 0.0 + + ! + ! set up Re bin "structure" for z_scaling + ! + hp%base_list(1)=0; + do j=1,Re_MAX_BIN + hp%step_list(j)=0.1+0.1*((j-1)**1.5); + if(hp%step_list(j)>Re_BIN_LENGTH) then + hp%step_list(j)=Re_BIN_LENGTH; + endif + if(j>1) then + hp%base_list(j)=hp%base_list(j-1)+floor(Re_BIN_LENGTH/hp%step_list(j-1)); + endif + enddo + + ! + ! set up Temperature bin structure used for z scaling + ! + do i=1,cnt_ice + hp%mt_tti(i)=(i-1)*5-90 + 273.15 + enddo + + do i=1,cnt_liq + hp%mt_ttl(i)=(i-1)*5-60 + 273.15 + enddo + + ! + ! define array discrete diameters used in mie calculations + ! + delt = (log(dmax)-log(dmin))/(nd-1) + deltp = exp(delt) + + hp%D(1) = dmin + do i=2,nd + hp%D(i) = hp%D(i-1)*deltp + enddo + + + end subroutine radar_simulator_init Index: LMDZ5/trunk/libf/phylmd/cosp/scale_LUTs_io.F90 =================================================================== --- LMDZ5/trunk/libf/phylmd/cosp/scale_LUTs_io.F90 (revision 2432) +++ LMDZ5/trunk/libf/phylmd/cosp/scale_LUTs_io.F90 (revision 2432) @@ -0,0 +1,137 @@ + ! scale_LUT_io: Contains subroutines to load and save scaling Look Up Tables (LUTs) to a file + ! + ! June 2010 Written by Roj Marchand + + module scale_LUTs_io + implicit none + + contains + + subroutine load_scale_LUTs(hp) + + use radar_simulator_types + + type(class_param), intent(inout) :: hp + + logical :: LUT_file_exists + integer :: i,j,k,ind + + ! + ! load scale LUT from file + ! + inquire(file=trim(hp%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat', & + exist=LUT_file_exists) + + if(.not.LUT_file_exists) then + + write(*,*) '*************************************************' + write(*,*) 'Warning: Could NOT FIND radar LUT file: ', & + trim(hp%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat' + write(*,*) 'Will calculated LUT values as needed' + write(*,*) '*************************************************' + + return + else + + OPEN(unit=12,file=trim(hp%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat',& + form='unformatted', & + err= 89, & + access='DIRECT',& + recl=28) + + write(*,*) 'Loading radar LUT file: ', & + trim(hp%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat' + + do i=1,maxhclass + do j=1,mt_ntt + do k=1,nRe_types + + ind = i+(j-1)*maxhclass+(k-1)*(nRe_types*mt_ntt) + + read(12,rec=ind) hp%Z_scale_flag(i,j,k), & + hp%Ze_scaled(i,j,k), & + hp%Zr_scaled(i,j,k), & + hp%kr_scaled(i,j,k) + + ! if(ind==1482329) then + ! write (*,*) ind, hp%Z_scale_flag(i,j,k), & + ! hp%Ze_scaled(i,j,k), & + ! hp%Zr_scaled(i,j,k), & + ! hp%kr_scaled(i,j,k) + !endif + + enddo + enddo + enddo + + close(unit=12) + return + endif + + 89 write(*,*) 'Error: Found but could NOT READ radar LUT file: ', & + trim(hp%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat' + stop + + end subroutine load_scale_LUTs + + subroutine save_scale_LUTs(hp) + + use radar_simulator_types + + type(class_param), intent(inout) :: hp + + logical :: LUT_file_exists + integer :: i,j,k,ind + + inquire(file=trim(hp%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat', & + exist=LUT_file_exists) + + OPEN(unit=12,file=trim(hp%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat',& + form='unformatted', & + err= 99, & + access='DIRECT',& + recl=28) + + write(*,*) 'Creating or Updating radar LUT file: ', & + trim(hp%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat' + + do i=1,maxhclass + do j=1,mt_ntt + do k=1,nRe_types + + ind = i+(j-1)*maxhclass+(k-1)*(nRe_types*mt_ntt) + + if(.not.LUT_file_exists .or. hp%Z_scale_added_flag(i,j,k)) then + + hp%Z_scale_added_flag(i,j,k)=.false. + + write(12,rec=ind) hp%Z_scale_flag(i,j,k), & + hp%Ze_scaled(i,j,k), & + hp%Zr_scaled(i,j,k), & + hp%kr_scaled(i,j,k) + + ! 1482329 T 0.170626345026495 0.00000000000000 1.827402935860823E-003 + + !if(ind==1482329) then + ! write (*,*) ind, hp%Z_scale_flag(i,j,k), & + ! hp%Ze_scaled(i,j,k), & + ! hp%Zr_scaled(i,j,k), & + ! hp%kr_scaled(i,j,k) + !endif + endif + enddo + enddo + enddo + + close(unit=12) + return + + 99 write(*,*) 'Error: Unable to create/update radar LUT file: ', & + trim(hp%scale_LUT_file_name) // '_radar_Z_scale_LUT.dat' + return + + end subroutine save_scale_LUTs + + + end module scale_LUTs_io + Index: LMDZ5/trunk/libf/phylmd/cosp/test_modis_simulator.F90.s =================================================================== --- LMDZ5/trunk/libf/phylmd/cosp/test_modis_simulator.F90.s (revision 2432) +++ LMDZ5/trunk/libf/phylmd/cosp/test_modis_simulator.F90.s (revision 2432) @@ -0,0 +1,154 @@ +! $Revision: 23 $, $Date: 2011-03-31 15:41:37 +0200 (jeu. 31 mars 2011) $ +! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/MODIS_simulator/test_modis_simulator.F90 $ +program test_modis_simulator + use mod_modis_sim + implicit none + ! + ! Tests cases for MODIS L2 (pixel) and L3 (grid-scale) simulators + ! + + ! + ! Inputs + ! + integer, parameter :: nLayers = 6, nSubCols = 10 + real, dimension(1, nLayers) :: temp, pressureLayers + real, dimension(1, nLayers+1) :: pressureLevels + real, dimension(1, nSubCols, nLayers) :: & + opticalThickness, cloudWater, cloudIce, waterSize, iceSize, & + liquid_opticalThickness, ice_opticalThickness + ! + ! Pixel-scale retreivals + ! + integer, dimension(1, nSubCols) :: retrievedPhase + real, dimension(1, nSubCols) :: isccpTau, isccpCloudTopPressure, & + retrievedCloudTopPressure, retrievedTau, retrievedSize + integer :: i, k + character(len = 6) :: phase + + ! + ! Grid mean properties + ! + real, dimension(1) :: & + Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, & + Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, & + Optical_Thickness_Total_LogMean, Optical_Thickness_Water_LogMean, Optical_Thickness_Ice_LogMean, & + Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, & + Cloud_Top_Pressure_Total_Mean, & + Liquid_Water_Path_Mean, Ice_Water_Path_Mean + real, dimension(1, numTauHistogramBins, numPressureHistogramBins) :: Optical_Thickness_vs_Cloud_Top_Pressure + + ! --------------------------------------------------------------------------------------------- + pressureLevels(1, :) = (/ 0., 200., 400., 600., 800., 1000., 1200. /) + pressureLayers(1, :) = (pressureLevels(1, 2:) + pressureLevels(1, 1:nLayers)) / 2 + temp(1, :) = (/ 200., 220., 240., 260., 280., 300.0 /) + isccpCloudTopPressure(:, :) = 999. + + + ! + ! Tests: these 9 from Steve Platnick + ! + ice_opticalThickness = 0.; liquid_opticalThickness = 0. + iceSize = 0.; waterSize = 0. + + ice_opticalThickness (1, 1, 3) = 0.2 + iceSize (1, 1, 3) = 30 + + ice_opticalThickness (1, 2, 3) = 1. + iceSize (1, 2, 3) = 30 + + liquid_opticalThickness(1, 3, 5) = 0.2 + waterSize (1, 3, 5) = 10 + + liquid_opticalThickness(1, 4, 5) = 1. + waterSize (1, 4, 5) = 10 + + ice_opticalThickness (1, 5, 3) = 0.2 + iceSize (1, 5, 3) = 30 + liquid_opticalThickness(1, 5, 5) = 5. + waterSize (1, 5, 5) = 10 + + ice_opticalThickness (1, 6, 3) = 1. + iceSize (1, 6, 3) = 30 + liquid_opticalThickness(1, 6, 5) = 5. + waterSize (1, 6, 5) = 10 + + ice_opticalThickness (1, 7, 3) = 10. + iceSize (1, 7, 3) = 30 + liquid_opticalThickness(1, 7, 5) = 5. + waterSize (1, 7, 5) = 10 + + ice_opticalThickness (1, 8, 2) = 1. + iceSize (1, 8, 2) = 15. + ice_opticalThickness (1, 8, 3) = 1. + iceSize (1, 8, 3) = 30. + ice_opticalThickness (1, 8, 4) = 1. + iceSize (1, 8, 4) = 60. + + liquid_opticalThickness(1, 9, 4) = 1. + waterSize (1, 9, 4) = 15 + liquid_opticalThickness(1, 9, 5) = 3. + waterSize (1, 9, 5) = 10. + liquid_opticalThickness(1, 9, 6) = 1. + waterSize (1, 9, 6) = 50. + + call modis_l2_simulator(temp(1,:), pressureLayers(1,:), pressureLevels(1,:), & + liquid_opticalThickness(1,:, :), ice_opticalThickness(1,:, :), & + waterSize(1,:, :), iceSize(1,:, :), & + isccpTau(1, :), isccpCloudTopPressure(1, :), & + retrievedPhase(1, :), retrievedCloudTopPressure(1, :), retrievedTau(1, :), retrievedSize(1, :)) + do i = 1, nSubcols + print * + print *, "Profile ", i + if(any(ice_opticalThickness(1, i, :) > 0.) .or. any(liquid_opticalThickness(1, i, :) > 0.)) & + print *, "center p ice: tau size water: tau size" + do k = 1, nLayers + if(ice_opticalThickness(1, i, k) > 0. .or. liquid_opticalThickness(1, i, k) > 0.) & + write(*, '(3x, f5.0, 7x, 2(f4.1, 3x, f4.1, 10x))'), & + pressureLayers(1, k), ice_opticalThickness(1, i, k), iceSize(1, i, k), liquid_opticalThickness(1, i, k), waterSize(1, i, k) + end do + select case(retrievedPhase(1, i)) + case(phaseIsNone) + phase = "None" + case(phaseIsIce) + phase = "Ice" + case(phaseIsLiquid) + phase = "Liquid" + case(phaseIsUndetermined) + phase = "Undet." + case default + phase = "???" + end select + print *, "Retrievals" + write(*, '(" Phase: ", a6, 2x, "CTP: ", f5.0, ", tau = ", f4.1, ", size = ", f4.1)') & + phase, retrievedCloudTopPressure(1, i), retrievedTau(1, i), retrievedSize(1, i) + end do + + call modis_L3_simulator(retrievedPhase, retrievedCloudTopPressure, retrievedTau, retrievedSize, & + Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean, & + Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean, & + Optical_Thickness_Total_LogMean, Optical_Thickness_Water_LogMean, Optical_Thickness_Ice_LogMean, & + Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean, & + Cloud_Top_Pressure_Total_Mean, & + Liquid_Water_Path_Mean, Ice_Water_Path_Mean, & + Optical_Thickness_vs_Cloud_Top_Pressure) + + print *; print *; print *, "Grid means" + do i = 1, nSubcols + write(*, '(" Phase: ", i2, 2x, "CTP: ", f5.0, ", tau = ", f4.1, ", size = ", f4.1)') & + retrievedPhase(1, i), retrievedCloudTopPressure(1, i), retrievedTau(1, i), retrievedSize(1, i) + end do + print *, " Total Liquid Ice" + write(*, '("Fraction: ", 3(f4.3, 6x))') Cloud_Fraction_Total_Mean, Cloud_Fraction_Water_Mean, Cloud_Fraction_Ice_Mean + write(*, '("Tau: ", 3(f4.1, 6x))') Optical_Thickness_Total_Mean, Optical_Thickness_Water_Mean, Optical_Thickness_Ice_Mean + write(*, '("TauLog: ", 3(f4.1, 6x))') Optical_Thickness_Total_LogMean, Optical_Thickness_Water_LogMean, Optical_Thickness_Ice_LogMean + write(*, '("Size: ", 10x, 2(f4.1, 6x))') Cloud_Particle_Size_Water_Mean, Cloud_Particle_Size_Ice_Mean + write(*, '("WaterPath:", 9x, 2(f5.1, 5x))') Liquid_Water_Path_Mean, Ice_Water_Path_Mean + write(*, '("CTP: ", f5.1)') Cloud_Top_Pressure_Total_Mean + print *; print *, "Histogram" + write(*, '(8x, 7(f4.1, 3x))') tauHistogramBoundaries(:) + do k = 1, numPressureHistogramBins + write(*, '(f5.0, 3x, 7(f4.2, 3x))') pressureHistogramBoundaries(k), Optical_Thickness_vs_Cloud_Top_Pressure(1, :, k) + end do + print *, "Total fraction from histogram:", sum(Optical_Thickness_vs_Cloud_Top_Pressure(1, :, :)) +end program test_modis_simulator +