Changeset 2713 for LMDZ5/trunk/libf


Ignore:
Timestamp:
Nov 24, 2016, 4:02:04 PM (8 years ago)
Author:
idelkadi
Message:

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

Location:
LMDZ5/trunk/libf/phylmd/cosp
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/cosp/cosp_constants.F90

    r2571 r2713  
    3535#include "cosp_defs.h"
    3636MODULE MOD_COSP_CONSTANTS
    37 
    38     use netcdf, only: nf90_fill_real
    3937    IMPLICIT NONE
    4038
    … …  
    5452    ! Missing value
    5553    real,parameter :: R_UNDEF = -1.0E30
    56 !    real,parameter :: R_UNDEF = nf90_fill_real
    5754
    5855    ! Number of possible output variables
    59     integer,parameter :: N_OUT_LIST = 63
    60     integer,parameter :: N3D = 8
     56    integer,parameter :: N_OUT_LIST = 65
     57    integer,parameter :: N3D = 10
    6158    integer,parameter :: N2D = 14
    6259    integer,parameter :: N1D = 40
    … …  
    108105                   -31.5,-28.5,-25.5,-22.5,-19.5,-16.5,-13.5,-10.5, -7.5, -4.5, &
    109106                    -1.5,  1.5,  4.5,  7.5, 10.5, 13.5, 16.5, 19.5, 22.5, 25.5/)
    110     real,parameter,dimension(2,LIDAR_NTEMP) :: LIDAR_PHASE_TEMP_BNDS=reshape(source=&
    111                     (/-273.15,-90.,-90.,-87.,-87.,-84.,-84.,-81.,-81.,-78., &
     107    real,parameter,dimension(2,LIDAR_NTEMP) :: LIDAR_PHASE_TEMP_BNDS=reshape(source=(/-273.15,-90.,-90.,-87.,-87.,-84.,-84.,-81.,-81.,-78., &
    112108                   -78.,-75.,-75.,-72.,-72.,-69.,-69.,-66.,-66.,-63., &
    113109                   -63.,-60.,-60.,-57.,-57.,-54.,-54.,-51.,-51.,-48., &

  • LMDZ5/trunk/libf/phylmd/cosp/cosp_modis_simulator.F90

    r2432 r2713  
    22!   Author: Robert Pincus, Cooperative Institute for Research in the Environmental Sciences
    33! All rights reserved.
    4 ! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $
     4! $Revision: 88 $, $Date: 2013-11-13 07:08:38 -0700 (Wed, 13 Nov 2013) $
    55! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/cosp_modis_simulator.F90 $
    66!
    … …  
    6565     !
    6666     real, dimension(:, :, :), pointer :: Optical_Thickness_vs_Cloud_Top_Pressure
     67     real, dimension(:, :, :), pointer :: Optical_Thickness_vs_ReffICE
     68     real, dimension(:, :, :), pointer :: Optical_Thickness_vs_ReffLIQ
    6769  end type COSP_MODIS
    6870 
    … …  
    115117       
    116118    real, dimension(count(gridBox%sunlit(:) > 0), numModisTauBins, numModisPressureBins) :: &
    117        jointHistogram
     119         jointHistogram
     120    real, dimension(count(gridBox%sunlit(:) > 0), numModisTauBins, numMODISReffIceBins) :: &
     121         jointHistogram2
     122    real, dimension(count(gridBox%sunlit(:) > 0), numModisTauBins, numMODISReffLiqBins) :: &
     123         jointHistogram3
    118124   
    119125    integer, dimension(count(gridBox%sunlit(:) >  0)) :: sunlit
    … …  
    214220                                retrievedPhase(i, :), retrievedCloudTopPressure(i, :),      &
    215221                                retrievedTau(i, :), retrievedSize(i, :))
    216       end do
    217       call modis_L3_simulator(retrievedPhase,              &
    218                               retrievedCloudTopPressure,   &
    219                               retrievedTau, retrievedSize, &
    220                               cfTotal,         cfLiquid,         cfIce,         &
    221                               cfHigh,          cfMid,            cfLow,         &
    222                               meanTauTotal,    meanTauLiquid,    meanTauIce,    &
    223                               meanLogTauTotal, meanLogTauLiquid, meanLogTauIce, &
    224                                                meanSizeLiquid,   meanSizeIce,   &
    225                               meanCloudTopPressure,                             &
    226                                                meanLiquidWaterPath, meanIceWaterPath, &
    227                               jointHistogram)
     222     end do
     223     
     224      ! DJS2015: Call L3 modis simulator used by cospv2.0
     225     ! call modis_L3_simulator(retrievedPhase,              &
     226     !                         retrievedCloudTopPressure,   &
     227     !                         retrievedTau, retrievedSize, &
     228     !                         cfTotal,         cfLiquid,         cfIce,         &
     229     !                         cfHigh,          cfMid,            cfLow,         &
     230     !                         meanTauTotal,    meanTauLiquid,    meanTauIce,    &
     231     !                         meanLogTauTotal, meanLogTauLiquid, meanLogTauIce, &
     232     !                         meanSizeLiquid,   meanSizeIce,   &
     233     !                         meanCloudTopPressure,                             &
     234     !                         meanLiquidWaterPath, meanIceWaterPath, &
     235     !                         jointHistogram)
     236     call modis_column(nSunlit,nSubcols,retrievedPhase,retrievedCloudTopPressure,   &
     237                        retrievedTau,retrievedSize,cfTotal,cfLiquid,cfIce,cfHigh,    &
     238                        cfMid,cfLow,meanTauTotal,meanTauLiquid,meanTauIce,           &
     239                        meanLogTauTotal,meanLogTauLiquid,meanLogTauIce,              &
     240                        meanSizeLiquid,meanSizeIce,meanCloudTopPressure,             &
     241                        meanLiquidWaterPath, meanIceWaterPath,                       &
     242                        jointHistogram,jointHistogram2,jointHistogram3)
     243      ! DJS2015: END
     244     
    228245      !
    229246      ! Copy results into COSP structure
    … …  
    254271     
    255272      modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(sunlit(:), 2:numModisTauBins+1, :) = jointHistogram(:, :, :)
     273      modisSim%Optical_Thickness_vs_ReffICE(sunlit(:),2:numModisTauBins+1,:)              = jointHistogram2(:, :, :)
     274      modisSim%Optical_Thickness_vs_ReffLIQ(sunlit(:),2:numModisTauBins+1,:)              = jointHistogram3(:, :, :)
    256275      !
    257276      ! Reorder pressure bins in joint histogram to go from surface to TOA
    258277      !
    259       modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(:,:,:) = &
    260         modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(:, :, numModisPressureBins:1:-1)
     278      modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(:,:,:) = modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(:, :, numModisPressureBins:1:-1)
    261279      if(nSunlit < nPoints) then
    262280        !
    … …  
    288306 
    289307        modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(notSunlit(:), :, :) = R_UNDEF
     308        modisSim%Optical_Thickness_vs_ReffICE(notSunlit(:), :, :) = R_UNDEF
     309        modisSim%Optical_Thickness_vs_ReffLIQ(notSunlit(:), :, :) = R_UNDEF
    290310      end if
    291311    else
    … …  
    318338 
    319339      modisSim%Optical_Thickness_vs_Cloud_Top_Pressure(:, :, :) = R_UNDEF
     340      modisSim%Optical_Thickness_vs_ReffICE(:, :, :) = R_UNDEF
     341      modisSim%Optical_Thickness_vs_ReffLIQ(:, :, :) = R_UNDEF
    320342    end if
    321343
    … …  
    363385     
    364386    allocate(x%Optical_Thickness_vs_Cloud_Top_Pressure(nPoints, numModisTauBins+1, numModisPressureBins))
     387    allocate(x%Optical_Thickness_vs_ReffICE(nPoints, numModisTauBins+1, numModisReffIceBins))
     388    allocate(x%Optical_Thickness_vs_ReffLIQ(nPoints, numModisTauBins+1, numModisReffLiqBins))
    365389    x%Optical_Thickness_vs_Cloud_Top_Pressure(:, :, :) = R_UNDEF
    366390  END SUBROUTINE CONSTRUCT_COSP_MODIS
    … …  
    400424   
    401425    if(associated(x%Optical_Thickness_vs_Cloud_Top_Pressure)) deallocate(x%Optical_Thickness_vs_Cloud_Top_Pressure   )
     426    if(associated(x%Optical_Thickness_vs_ReffIce)) deallocate(x%Optical_Thickness_vs_ReffIce)
     427    if(associated(x%Optical_Thickness_vs_ReffLiq)) deallocate(x%Optical_Thickness_vs_ReffLiq)
    402428  END SUBROUTINE FREE_COSP_MODIS
    403429  ! -----------------------------------------------------
    … …  
    447473   
    448474    copy%Optical_Thickness_vs_Cloud_Top_Pressure(copy_start:copy_end, :, :) = &
    449                           orig%Optical_Thickness_vs_Cloud_Top_Pressure(orig_start:orig_end, :, :)
     475         orig%Optical_Thickness_vs_Cloud_Top_Pressure(orig_start:orig_end, :, :)
     476    copy%Optical_Thickness_vs_ReffIce(copy_start:copy_end, :, :) = &
     477         orig%Optical_Thickness_vs_ReffIce(orig_start:orig_end, :, :)
     478    copy%Optical_Thickness_vs_ReffLiq(copy_start:copy_end, :, :) = &
     479         orig%Optical_Thickness_vs_ReffLiq(orig_start:orig_end, :, :)
     480
    450481  END SUBROUTINE COSP_MODIS_CPSECTION
    451482  ! -----------------------------------------------------

  • LMDZ5/trunk/libf/phylmd/cosp/cosp_output_mod.F90

    r2652 r2713  
    99  USE MOD_COSP_TYPES
    1010  use MOD_COSP_Modis_Simulator, only : cosp_modis
     11  use mod_modis_sim, only : numMODISReffIceBins, reffICE_binCenters, &
     12                            numMODISReffLiqBins, reffLIQ_binCenters
    1113
    1214! cosp_output_mod
    … …  
    1719!$OMP THREADPRIVATE(cosp_outfilekeys, cosp_nidfiles)
    1820      INTEGER, DIMENSION(3), SAVE  :: nhoricosp,nvert,nvertmcosp,nvertcol,nvertbze, &
    19                                       nvertsratio,nvertisccp,nvertp,nverttemp,nvertmisr
     21                                      nvertsratio,nvertisccp,nvertp,nverttemp,nvertmisr, &
     22                                      nvertReffIce,nvertReffLiq
    2023      REAL, DIMENSION(3), SAVE                :: zoutm_cosp
    2124!$OMP THREADPRIVATE(nhoricosp, nvert,nvertmcosp,nvertcol,nvertsratio,nvertbze,nvertisccp,nvertp,zoutm_cosp,nverttemp,nvertmisr)
     25!$OMP THREADPRIVATE(nvertReffIce,nvertReffLiq)
    2226      REAL, SAVE                   :: zdtimemoy_cosp
    2327!$OMP THREADPRIVATE(zdtimemoy_cosp)
    … …  
    176180  TYPE(ctrl_outcosp), SAVE :: o_clmodis = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
    177181         "clmodis", "MODIS Cloud Area Fraction", "%", (/ ('', i=1, 3) /))
     182  TYPE(ctrl_outcosp), SAVE :: o_crimodis = ctrl_outcosp((/ .TRUE., .TRUE., .TRUE. /), &
     183         "crimodis", "Optical_Thickness_vs_ReffIce from Modis", "%", (/ ('',i=1, 3) /))         
     184  TYPE(ctrl_outcosp), SAVE :: o_crlmodis = ctrl_outcosp((/ .TRUE., .TRUE.,.TRUE. /), &
     185         "crlmodis", "Optical_Thickness_vs_ReffLiq from Modis", "%", (/ ('',i=1, 3) /))         
    178186
    179187! Rttovs simulator
    … …  
    325333!   CALL wxios_add_vaxis("dbze", DBZE_BINS, dbze_ax)
    326334!   CALL wxios_add_vaxis("scatratio", SR_BINS, sratio_ax)
     335   CALL wxios_add_vaxis("ReffIce", numMODISReffIceBins, reffICE_binCenters)
     336   CALL wxios_add_vaxis("ReffLiq", numMODISReffLiqBins, reffLIQ_binCenters)
     337
    327338#endif
    328339   
    … …  
    366377
    367378      CALL histvert(cosp_nidfiles(iff),"cth16","altitude","m",MISR_N_CTH,MISR_CTH,nvertmisr(iff))
     379 
     380      CALL histvert(cosp_nidfiles(iff),"ReffIce","Effective_particle_size_Ice","microns",numMODISReffIceBins, reffICE_binCenters, &
     381                    nvertReffIce(iff))                                         
     382     
     383      CALL histvert(cosp_nidfiles(iff),"ReffLiq","Effective_particle_size_Liq","microns",numMODISReffLiqBins, reffLIQ_binCenters, &                                 
     384                    nvertReffLiq(iff))
    368385
    369386!      CALL histvert(cosp_nidfiles(iff),"dbze","equivalent_reflectivity_factor","dBZ",DBZE_BINS,dbze_ax,nvertbze(iff))

  • LMDZ5/trunk/libf/phylmd/cosp/cosp_output_write_mod.F90

    r2560 r2713  
    372372 endif
    373373
     374    where(modis%Optical_Thickness_vs_ReffIce == R_UNDEF) &
     375          modis%Optical_Thickness_vs_ReffIce = Cosp_fill_value
     376
     377    where(modis%Optical_Thickness_vs_ReffLiq == R_UNDEF) &
     378          modis%Optical_Thickness_vs_ReffLiq = Cosp_fill_value
     379
     380    do icl=1,7
     381   CALL histwrite3d_cosp(o_crimodis, &
     382     modis%Optical_Thickness_vs_ReffIce(:,icl,:),nvertReffIce,icl)
     383   CALL histwrite3d_cosp(o_crlmodis, &
     384     modis%Optical_Thickness_vs_ReffLiq(:,icl,:),nvertReffLiq,icl)
     385    enddo
     386
    374387 IF(.NOT.cosp_varsdefined) THEN
    375388!$OMP MASTER
    … …  
    521534      ELSE IF (nvertsave.eq.nvertmisr(iff)) THEN
    522535          klevs=MISR_N_CTH
    523            nam_axvert="cth16"
     536          nam_axvert="cth16"
     537      ELSE IF (nvertsave.eq.nvertReffIce(iff)) THEN
     538          klevs= numMODISReffIceBins
     539          nam_axvert="ReffIce"
     540      ELSE IF (nvertsave.eq.nvertReffLiq(iff)) THEN
     541          klevs= numMODISReffLiqBins
     542          nam_axvert="ReffLiq"
    524543      ELSE
    525544           klevs=Nlevout

  • LMDZ5/trunk/libf/phylmd/cosp/cosp_types.F90

    r2428 r2713  
    5151                Lcltmodis,Lclwmodis,Lclimodis,Lclhmodis,Lclmmodis,Lcllmodis,Ltautmodis,Ltauwmodis,Ltauimodis,Ltautlogmodis, &
    5252                Ltauwlogmodis,Ltauilogmodis,Lreffclwmodis,Lreffclimodis,Lpctmodis,Llwpmodis, &
    53                 Liwpmodis,Lclmodis
     53                Liwpmodis,Lclmodis,Lcrimodis,Lcrlmodis
    5454
    5555     character(len=32) :: out_list(N_OUT_LIST)

  • LMDZ5/trunk/libf/phylmd/cosp/modis_simulator.F90

    r2432 r2713  
    22!   Author: Robert Pincus, Cooperative Institute for Research in the Environmental Sciences
    33! All rights reserved.
    4 ! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $
     4! $Revision: 88 $, $Date: 2013-11-13 07:08:38 -0700 (Wed, 13 Nov 2013) $
    55! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/MODIS_simulator/modis_simulator.F90 $
    66!
    … …  
    7979  real,    parameter :: re_water_min= 4., re_water_max= 30., re_ice_min= 5., re_ice_max= 90.
    8080  integer, parameter :: num_trial_res = 15             ! increase to make the linear pseudo-retrieval of size more accurate
    81   logical, parameter :: use_two_re_iterations = .false. ! do two retrieval iterations?
    82  
     81! DJS2015: Remove unused parameter
     82!  logical, parameter :: use_two_re_iterations = .false. ! do two retrieval iterations?
     83! DJS2015 END 
    8384  !
    8485  ! Precompute near-IR optical params vs size for retrieval scheme
    … …  
    125126    nominalPressureHistogramCenters = (nominalPressureHistogramBoundaries(1, :) + &
    126127                                       nominalPressureHistogramBoundaries(2, :) ) / 2.
     128  ! DJS2015 START: Add bin descriptions for joint-histograms of partice-sizes and optical depth. This is
     129  !                identical to what is done in COSPv.2.0.0 for histogram bin initialization.
     130  integer :: j
     131  integer,parameter :: &
     132       numMODISReffLiqBins = 6, & ! Number of bins for tau/ReffLiq joint-histogram
     133       numMODISReffIceBins = 6    ! Number of bins for tau/ReffICE joint-histogram
     134  real,parameter,dimension(numMODISReffLiqBins+1) :: &
     135       reffLIQ_binBounds = (/0., 8e-6, 1.0e-5, 1.3e-5, 1.5e-5, 2.0e-5, 3.0e-5/)
     136  real,parameter,dimension(numMODISReffIceBins+1) :: &
     137       reffICE_binBounds = (/0., 1.0e-5, 2.0e-5, 3.0e-5, 4.0e-5, 6.0e-5, 9.0e-5/)
     138  real,parameter,dimension(2,numMODISReffIceBins) :: &
     139       reffICE_binEdges = reshape(source=(/reffICE_binBounds(1),((reffICE_binBounds(k),  &
     140                                  l=1,2),k=2,numMODISReffIceBins),reffICE_binBounds(numMODISReffIceBins+1)/),  &
     141                                  shape = (/2,numMODISReffIceBins/))
     142  real,parameter,dimension(2,numMODISReffLiqBins) :: &
     143       reffLIQ_binEdges = reshape(source=(/reffLIQ_binBounds(1),((reffLIQ_binBounds(k),  &
     144                                  l=1,2),k=2,numMODISReffLiqBins),reffLIQ_binBounds(numMODISReffIceBins+1)/),  &
     145                                  shape = (/2,numMODISReffLiqBins/))             
     146  real,parameter,dimension(numMODISReffIceBins) :: &
     147       reffICE_binCenters = (reffICE_binEdges(1,:)+reffICE_binEdges(2,:))/2.
     148  real,parameter,dimension(numMODISReffLiqBins) :: &
     149       reffLIQ_binCenters = (reffLIQ_binEdges(1,:)+reffLIQ_binEdges(2,:))/2.
     150  ! DJS2015 END
     151
    127152  ! ------------------------------
    128153  ! There are two ways to call the MODIS simulator:
    … …  
    384409                               
    385410  end subroutine modis_L2_simulator_oneTau
     411
     412  ! ########################################################################################
     413  subroutine modis_column(nPoints,nSubCols,phase, cloud_top_pressure, optical_thickness, particle_size,      &
     414       Cloud_Fraction_Total_Mean,         Cloud_Fraction_Water_Mean,         Cloud_Fraction_Ice_Mean,        &
     415       Cloud_Fraction_High_Mean,          Cloud_Fraction_Mid_Mean,           Cloud_Fraction_Low_Mean,        &
     416       Optical_Thickness_Total_Mean,      Optical_Thickness_Water_Mean,      Optical_Thickness_Ice_Mean,     &
     417       Optical_Thickness_Total_MeanLog10, Optical_Thickness_Water_MeanLog10, Optical_Thickness_Ice_MeanLog10,&
     418       Cloud_Particle_Size_Water_Mean,    Cloud_Particle_Size_Ice_Mean,      Cloud_Top_Pressure_Total_Mean,  &
     419       Liquid_Water_Path_Mean,            Ice_Water_Path_Mean,                                               &   
     420       Optical_Thickness_vs_Cloud_Top_Pressure,Optical_Thickness_vs_ReffIce,Optical_Thickness_vs_ReffLiq)
     421   
     422    ! INPUTS
     423    integer,intent(in) :: &
     424         nPoints,                           & ! Number of horizontal gridpoints
     425         nSubCols                             ! Number of subcolumns
     426    integer,intent(in), dimension(:,:) ::  &
     427!ds    integer,intent(in), dimension(nPoints, nSubCols) ::  &
     428         phase                             
     429    real,intent(in),dimension(:,:) ::  &
     430!ds    real,intent(in),dimension(nPoints, nSubCols) ::  &
     431         cloud_top_pressure,                &
     432         optical_thickness,                 &
     433         particle_size
     434 
     435    ! OUTPUTS
     436    real,intent(inout),dimension(:)  ::   & !
     437!ds    real,intent(inout),dimension(nPoints)  ::   & !
     438         Cloud_Fraction_Total_Mean,         & !
     439         Cloud_Fraction_Water_Mean,         & !
     440         Cloud_Fraction_Ice_Mean,           & !
     441         Cloud_Fraction_High_Mean,          & !
     442         Cloud_Fraction_Mid_Mean,           & !
     443         Cloud_Fraction_Low_Mean,           & !
     444         Optical_Thickness_Total_Mean,      & !
     445         Optical_Thickness_Water_Mean,      & !
     446         Optical_Thickness_Ice_Mean,        & !
     447         Optical_Thickness_Total_MeanLog10, & !
     448         Optical_Thickness_Water_MeanLog10, & !
     449         Optical_Thickness_Ice_MeanLog10,   & !
     450         Cloud_Particle_Size_Water_Mean,    & !
     451         Cloud_Particle_Size_Ice_Mean,      & !
     452         Cloud_Top_Pressure_Total_Mean,     & !
     453         Liquid_Water_Path_Mean,            & !
     454         Ice_Water_Path_Mean                  !
     455    real,intent(inout),dimension(:,:,:) :: &
     456!ds    real,intent(inout),dimension(nPoints,numTauHistogramBins,numPressureHistogramBins) :: &
     457         Optical_Thickness_vs_Cloud_Top_Pressure
     458    real,intent(inout),dimension(:,:,:) :: &   
     459!ds    real,intent(inout),dimension(nPoints,numTauHistogramBins,numMODISReffIceBins) :: &   
     460         Optical_Thickness_vs_ReffIce
     461    real,intent(inout),dimension(:,:,:) :: &   
     462!ds    real,intent(inout),dimension(nPoints,numTauHistogramBins,numMODISReffLiqBins) :: &   
     463         Optical_Thickness_vs_ReffLiq         
     464
     465    ! LOCAL VARIABLES
     466    real, parameter :: &
     467         LWP_conversion = 2./3. * 1000. ! MKS units 
     468    integer :: i, j
     469    logical, dimension(nPoints,nSubCols) :: &
     470         cloudMask,      &
     471         waterCloudMask, &
     472         iceCloudMask,   &
     473         validRetrievalMask
     474    real,dimension(nPoints,nSubCols) :: &
     475         tauWRK,ctpWRK,reffIceWRK,reffLiqWRK
     476   
     477    ! ########################################################################################
     478    ! Include only those pixels with successful retrievals in the statistics
     479    ! ########################################################################################
     480    validRetrievalMask(1:nPoints,1:nSubCols) = particle_size(1:nPoints,1:nSubCols) > 0.
     481    cloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) /= phaseIsNone .and.       &
     482         validRetrievalMask(1:nPoints,1:nSubCols)
     483    waterCloudMask(1:nPoints,1:nSubCols) = phase(1:nPoints,1:nSubCols) == phaseIsLiquid .and. &
     484         validRetrievalMask(1:nPoints,1:nSubCols)
     485    iceCloudMask(1:nPoints,1:nSubCols)   = phase(1:nPoints,1:nSubCols) == phaseIsIce .and.    &
     486         validRetrievalMask(1:nPoints,1:nSubCols)
     487
     488    ! ########################################################################################
     489    ! Use these as pixel counts at first
     490    ! ########################################################################################
     491    Cloud_Fraction_Total_Mean(1:nPoints) = real(count(cloudMask,      dim = 2))
     492    Cloud_Fraction_Water_Mean(1:nPoints) = real(count(waterCloudMask, dim = 2))
     493    Cloud_Fraction_Ice_Mean(1:nPoints)   = real(count(iceCloudMask,   dim = 2))
     494    Cloud_Fraction_High_Mean(1:nPoints)  = real(count(cloudMask .and. cloud_top_pressure <=          &
     495                                           highCloudPressureLimit, dim = 2))
     496    Cloud_Fraction_Low_Mean(1:nPoints)   = real(count(cloudMask .and. cloud_top_pressure >           &
     497                                           lowCloudPressureLimit,  dim = 2))
     498    Cloud_Fraction_Mid_Mean(1:nPoints)   = Cloud_Fraction_Total_Mean(1:nPoints) - Cloud_Fraction_High_Mean(1:nPoints)&
     499                                           - Cloud_Fraction_Low_Mean(1:nPoints)
     500
     501    ! ########################################################################################
     502    ! Compute mean optical thickness.
     503    ! ########################################################################################
     504    Optical_Thickness_Total_Mean(1:nPoints) = sum(optical_thickness, mask = cloudMask,      dim = 2) / &
     505                                              Cloud_Fraction_Total_Mean(1:nPoints)
     506    Optical_Thickness_Water_Mean(1:nPoints) = sum(optical_thickness, mask = waterCloudMask, dim = 2) / &
     507                                              Cloud_Fraction_Water_Mean(1:nPoints)
     508    Optical_Thickness_Ice_Mean(1:nPoints)   = sum(optical_thickness, mask = iceCloudMask,   dim = 2) / &
     509                                              Cloud_Fraction_Ice_Mean(1:nPoints)
     510       
     511    ! ########################################################################################
     512    ! We take the absolute value of optical thickness here to satisfy compilers that complains
     513    ! when we evaluate the logarithm of a negative number, even though it's not included in
     514    ! the sum.
     515    ! ########################################################################################
     516    Optical_Thickness_Total_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = cloudMask, &
     517         dim = 2) / Cloud_Fraction_Total_Mean(1:nPoints)
     518    Optical_Thickness_Water_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = waterCloudMask,&
     519         dim = 2) / Cloud_Fraction_Water_Mean(1:nPoints)
     520    Optical_Thickness_Ice_MeanLog10(1:nPoints) = sum(log10(abs(optical_thickness)), mask = iceCloudMask,&
     521         dim = 2) / Cloud_Fraction_Ice_Mean(1:nPoints)
     522    Cloud_Particle_Size_Water_Mean(1:nPoints) = sum(particle_size, mask = waterCloudMask, dim = 2) / &
     523         Cloud_Fraction_Water_Mean(1:nPoints)
     524    Cloud_Particle_Size_Ice_Mean(1:nPoints) = sum(particle_size, mask = iceCloudMask,   dim = 2) / &
     525         Cloud_Fraction_Ice_Mean(1:nPoints)
     526    Cloud_Top_Pressure_Total_Mean(1:nPoints) = sum(cloud_top_pressure, mask = cloudMask, dim = 2) / &
     527         max(1, count(cloudMask, dim = 2))
     528    Liquid_Water_Path_Mean(1:nPoints) = LWP_conversion*sum(particle_size*optical_thickness, &
     529         mask=waterCloudMask,dim=2)/Cloud_Fraction_Water_Mean(1:nPoints)
     530    Ice_Water_Path_Mean(1:nPoints) = LWP_conversion * ice_density*sum(particle_size*optical_thickness,&
     531         mask=iceCloudMask,dim = 2) /Cloud_Fraction_Ice_Mean(1:nPoints)
     532
     533    ! ########################################################################################
     534    ! Normalize pixel counts to fraction. The first three cloud fractions have been set to -1
     535    ! in cloud-free areas, so set those places to 0.
     536    ! ########################################################################################
     537    Cloud_Fraction_High_Mean(1:nPoints)  = Cloud_Fraction_High_Mean(1:nPoints) /nSubcols
     538    Cloud_Fraction_Mid_Mean(1:nPoints)   = Cloud_Fraction_Mid_Mean(1:nPoints)  /nSubcols
     539    Cloud_Fraction_Low_Mean(1:nPoints)   = Cloud_Fraction_Low_Mean(1:nPoints)  /nSubcols
     540
     541    ! ########################################################################################
     542    ! Set clear-scenes to undefined
     543    ! ########################################################################################
     544    where (Cloud_Fraction_Total_Mean == 0)
     545       Optical_Thickness_Total_Mean      = R_UNDEF
     546       Optical_Thickness_Total_MeanLog10 = R_UNDEF
     547       Cloud_Top_Pressure_Total_Mean     = R_UNDEF
     548    endwhere
     549    where (Cloud_Fraction_Water_Mean == 0)
     550       Optical_Thickness_Water_Mean      = R_UNDEF
     551       Optical_Thickness_Water_MeanLog10 = R_UNDEF
     552       Cloud_Particle_Size_Water_Mean    = R_UNDEF
     553       Liquid_Water_Path_Mean            = R_UNDEF
     554    endwhere
     555    where (Cloud_Fraction_Ice_Mean == 0)
     556       Optical_Thickness_Ice_Mean        = R_UNDEF
     557       Optical_Thickness_Ice_MeanLog10   = R_UNDEF
     558       Cloud_Particle_Size_Ice_Mean      = R_UNDEF
     559       Ice_Water_Path_Mean               = R_UNDEF
     560    endwhere
     561    where (Cloud_Fraction_High_Mean == 0)  Cloud_Fraction_High_Mean = R_UNDEF
     562    where (Cloud_Fraction_Mid_Mean == 0)   Cloud_Fraction_Mid_Mean = R_UNDEF
     563    where (Cloud_Fraction_Low_Mean == 0)   Cloud_Fraction_Low_Mean = R_UNDEF
     564
     565    ! ########################################################################################
     566    ! Joint histogram 
     567    ! ########################################################################################
     568
     569    ! Loop over all points
     570    tauWRK(1:nPoints,1:nSubCols)     = optical_thickness(1:nPoints,1:nSubCols)
     571    ctpWRK(1:nPoints,1:nSubCols)     = cloud_top_pressure(1:nPoints,1:nSubCols)
     572    reffIceWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,iceCloudMask)
     573    reffLiqWRK(1:nPoints,1:nSubCols) = merge(particle_size,R_UNDEF,waterCloudMask)
     574    do j=1,nPoints
     575
     576       ! Fill clear and optically thin subcolumns with fill
     577       where(.not. cloudMask(j,1:nSubCols))
     578          tauWRK(j,1:nSubCols) = -999.
     579          ctpWRK(j,1:nSubCols) = -999.
     580       endwhere
     581       ! Joint histogram of tau/CTP
     582       call hist2D(tauWRK(j,1:nSubCols),ctpWRK(j,1:nSubCols),nSubCols,&
     583                   tauHistogramBoundaries,numTauHistogramBins,&
     584                   pressureHistogramBoundaries,numPressureHistogramBins,&
     585                   Optical_Thickness_vs_Cloud_Top_Pressure(j,1:numTauHistogramBins,1:numPressureHistogramBins))
     586       ! Joint histogram of tau/ReffICE
     587       call hist2D(tauWRK(j,1:nSubCols),reffIceWrk(j,1:nSubCols),nSubCols,               &
     588                   tauHistogramBoundaries,numTauHistogramBins,reffICE_binBounds,         &
     589                   numMODISReffIceBins, Optical_Thickness_vs_ReffIce(j,1:numTauHistogramBins,1:numMODISReffIceBins))
     590       ! Joint histogram of tau/ReffLIQ
     591       call hist2D(tauWRK(j,1:nSubCols),reffLiqWrk(j,1:nSubCols),nSubCols,               &
     592                   tauHistogramBoundaries,numTauHistogramBins,reffLIQ_binBounds,         &
     593                   numMODISReffLiqBins, Optical_Thickness_vs_ReffLiq(j,1:numTauHistogramBins,1:numMODISReffLiqBins))                   
     594
     595    enddo   
     596    Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numTauHistogramBins,1:numPressureHistogramBins) = &
     597         Optical_Thickness_vs_Cloud_Top_Pressure(1:nPoints,1:numTauHistogramBins,1:numPressureHistogramBins)/nSubCols
     598    Optical_Thickness_vs_ReffIce(1:nPoints,1:numTauHistogramBins,1:numMODISReffIceBins) = &
     599         Optical_Thickness_vs_ReffIce(1:nPoints,1:numTauHistogramBins,1:numMODISReffIceBins)/nSubCols
     600    Optical_Thickness_vs_ReffLiq(1:nPoints,1:numTauHistogramBins,1:numMODISReffLiqBins) = &
     601         Optical_Thickness_vs_ReffLiq(1:nPoints,1:numTauHistogramBins,1:numMODISReffLiqBins)/nSubCols
     602
     603  end subroutine modis_column
     604  ! ######################################################################################
     605  ! SUBROUTINE hist2D
     606  ! ######################################################################################
     607  subroutine hist2D(var1,var2,npts,bin1,nbin1,bin2,nbin2,jointHist)
     608    implicit none
     609   
     610    ! INPUTS
     611    integer, intent(in) :: &
     612         npts,  & ! Number of data points to be sorted
     613         nbin1, & ! Number of bins in histogram direction 1
     614         nbin2    ! Number of bins in histogram direction 2
     615    real,intent(in),dimension(npts) :: &
     616         var1,  & ! Variable 1 to be sorted into bins
     617         var2     ! variable 2 to be sorted into bins
     618    real,intent(in),dimension(nbin1+1) :: &
     619         bin1     ! Histogram bin 1 boundaries
     620    real,intent(in),dimension(nbin2+1) :: &
     621         bin2     ! Histogram bin 2 boundaries
     622    ! OUTPUTS
     623    real,intent(out),dimension(nbin1,nbin2) :: &
     624         jointHist
     625   
     626    ! LOCAL VARIABLES
     627    integer :: ij,ik
     628   
     629    do ij=2,nbin1+1
     630       do ik=2,nbin2+1
     631          jointHist(ij-1,ik-1)=count(var1 .ge. bin1(ij-1) .and. var1 .lt. bin1(ij) .and. &
     632               var2 .ge. bin2(ik-1) .and. var2 .lt. bin2(ik))       
     633       enddo
     634    enddo
     635  end subroutine hist2D
     636 
    386637  !------------------------------------------------------------------------------------------------
    387638  subroutine modis_L3_simulator(phase, cloud_top_pressure, optical_thickness, particle_size,            &
    … …  
    666917      ! If first retrieval works, can try 2nd iteration using greater re resolution
    667918      !
    668       if(use_two_re_iterations .and. retrieve_re > 0.) then
    669         re_min = retrieve_re - delta_re
    670         re_max = retrieve_re + delta_re
    671         delta_re = (re_max - re_min)/real(num_trial_res-1)
    672  
    673         trial_re(:) = re_min + delta_re * (/ (i - 1, i = 1, num_trial_res) /)
    674         g(:)  = get_g_nir(  phase, trial_re(:))
    675         w0(:) = get_ssa_nir(phase, trial_re(:))
    676         predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:))
    677         retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir)
    678       end if
     919! DJS2015: Remove unused piece of code     
     920!      if(use_two_re_iterations .and. retrieve_re > 0.) then
     921!        re_min = retrieve_re - delta_re
     922!        re_max = retrieve_re + delta_re
     923!        delta_re = (re_max - re_min)/real(num_trial_res-1)
     924! 
     925!        trial_re(:) = re_min + delta_re * (/ (i - 1, i = 1, num_trial_res) /)
     926!        g(:)  = get_g_nir(  phase, trial_re(:))
     927!        w0(:) = get_ssa_nir(phase, trial_re(:))
     928!        predicted_Refl_nir(:) = two_stream_reflectance(tau, g(:), w0(:))
     929!        retrieve_re = interpolate_to_min(trial_re(:), predicted_Refl_nir(:), obs_Refl_nir)
     930!      end if
     931! DJS2015 END
    679932    else
    680933      retrieve_re = re_fill
    … …  
    739992    real,    intent(in) :: re
    740993    real :: get_g_nir
    741    
    742     real, dimension(3), parameter :: ice_coefficients   = (/ 0.7432,  4.5563e-3, -2.8697e-5 /), &
    743                                small_water_coefficients = (/ 0.8027, -1.0496e-2,  1.7071e-3 /), &
    744                                  big_water_coefficients = (/ 0.7931,  5.3087e-3, -7.4995e-5 /)
    745    
    746     ! approx. fits from MODIS Collection 5 LUT scattering calculations
    747     if(phase == phaseIsLiquid) then
    748       if(re < 8.) then
    749         get_g_nir = fit_to_quadratic(re, small_water_coefficients)
    750         if(re < re_water_min) get_g_nir = fit_to_quadratic(re_water_min, small_water_coefficients)
    751       else
    752         get_g_nir = fit_to_quadratic(re,  big_water_coefficients)
    753         if(re > re_water_max) get_g_nir = fit_to_quadratic(re_water_max, big_water_coefficients)
    754       end if
     994
     995    real, dimension(3), parameter :: ice_coefficients         = (/ 0.7490, 6.5153e-3, -5.4136e-5 /), &
     996                                     small_water_coefficients = (/ 1.0364, -8.8800e-2, 7.0000e-3 /)
     997    real, dimension(4), parameter :: big_water_coefficients   = (/ 0.6035, 2.8993e-2, -1.1051e-3, 1.5134e-5 /)
     998
     999    ! approx. fits from MODIS Collection 6 LUT scattering calculations for 3.7 µm channel size retrievals
     1000    if(phase == phaseIsLiquid) then 
     1001       if(re < 7.) then
     1002          get_g_nir = fit_to_quadratic(re, small_water_coefficients)
     1003          if(re < re_water_min) get_g_nir = fit_to_quadratic(re_water_min, small_water_coefficients)
     1004       else
     1005          get_g_nir = fit_to_cubic(re, big_water_coefficients)
     1006          if(re > re_water_max) get_g_nir = fit_to_cubic(re_water_max, big_water_coefficients)
     1007       end if
    7551008    else
    756       get_g_nir = fit_to_quadratic(re, ice_coefficients)
     1009       get_g_nir = fit_to_quadratic(re, ice_coefficients)
    7571010      if(re < re_ice_min) get_g_nir = fit_to_quadratic(re_ice_min, ice_coefficients)
    7581011      if(re > re_ice_max) get_g_nir = fit_to_quadratic(re_ice_max, ice_coefficients)
    … …  
    7711024        ! Fits from Steve Platnick
    7721025        !
     1026        real, dimension(4), parameter :: ice_coefficients   = (/ 0.9625, -1.8069e-2, 3.3281e-4,-2.2865e-6/)
     1027        real, dimension(3), parameter :: water_coefficients = (/ 1.0044, -1.1397e-2, 1.3300e-4 /)
    7731028       
    774         real, dimension(4), parameter :: ice_coefficients   = (/ 0.9994, -4.5199e-3, 3.9370e-5, -1.5235e-7 /)
    775         real, dimension(3), parameter :: water_coefficients = (/ 1.0008, -2.5626e-3, 1.6024e-5 /)
    776        
    777         ! approx. fits from MODIS Collection 5 LUT scattering calculations
     1029        ! approx. fits from MODIS Collection 6 LUT scattering calculations
    7781030        if(phase == phaseIsLiquid) then
    7791031          get_ssa_nir = fit_to_quadratic(re, water_coefficients)

  • LMDZ5/trunk/libf/phylmd/cosp/read_cosp_output_nl.F90

    r2428 r2713  
    4141             Lcltmodis,Lclwmodis,Lclimodis,Lclhmodis,Lclmmodis,Lcllmodis,Ltautmodis,Ltauwmodis,Ltauimodis,Ltautlogmodis, &
    4242             Ltauwlogmodis,Ltauilogmodis,Lreffclwmodis,Lreffclimodis,Lpctmodis,Llwpmodis, &
    43              Liwpmodis,Lclmodis
     43             Liwpmodis,Lclmodis,Lcrimodis,Lcrlmodis
    4444
    4545  namelist/COSP_OUTPUT/Lradar_sim,Llidar_sim,Lisccp_sim,Lmodis_sim,Lmisr_sim,Lrttov_sim, &
    … …  
    5757             Lcltmodis,Lclwmodis,Lclimodis,Lclhmodis,Lclmmodis,Lcllmodis,Ltautmodis,Ltauwmodis,Ltauimodis,Ltautlogmodis, &
    5858             Ltauwlogmodis,Ltauilogmodis,Lreffclwmodis,Lreffclimodis,Lpctmodis,Llwpmodis, &
    59              Liwpmodis,Lclmodis
     59             Liwpmodis,Lclmodis,Lcrimodis,Lcrlmodis
    6060   
    6161  do i=1,N_OUT_LIST
    … …  
    137137  CALL bcast(Lclmodis)
    138138  CALL bcast(Ltbrttov)
     139  CALL bcast(Lcrimodis)
     140  CALL bcast(Lcrlmodis)
     141
    139142!$OMP BARRIER
    140143
    … …  
    223226    Liwpmodis=.false.
    224227    Lclmodis=.false.
     228    Lcrimodis=.false.
     229    Lcrlmodis=.false.
    225230  endif
    226231  if (Lmodis_sim) Lisccp_sim = .true.
    … …  
    381386  i = i+1
    382387  if (Lclmodis)         cfg%out_list(i) = 'clmodis'
     388  i = i+1
     389  if (Lcrimodis)         cfg%out_list(i) = 'crimodis'
     390  i = i+1
     391  if (Lcrlmodis)         cfg%out_list(i) = 'crlmodis'
    383392   
    384393  if (i /= N_OUT_LIST) then
    … …  
    459468  cfg%Liwpmodis=Liwpmodis
    460469  cfg%Lclmodis=Lclmodis
    461  
     470  cfg%Lcrimodis=Lcrimodis
     471  cfg%Lcrlmodis=Lcrlmodis
     472
    462473 END SUBROUTINE READ_COSP_OUTPUT_NL
    463474
Note: See TracChangeset for help on using the changeset viewer.