Changeset 1414 for LMDZ4/trunk


Ignore:
Timestamp:
Jul 15, 2010, 5:21:22 PM (14 years ago)
Author:
idelkadi
Message:

Passage a la version cosp.v1.3 pour le Lidar et ISCCP
Corrections de bugs pour ISCCP et optimisation pour le Lidar

Location:
LMDZ4/trunk/libf/cosp
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • LMDZ4/trunk/libf/cosp/cosp.F90

    r1327 r1414  
    2323! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
    2424
     25!!#include "cosp_defs.h"
    2526MODULE MOD_COSP
    2627  USE MOD_COSP_TYPES
    … …  
    128129      !     and reff_zero    == .false.  Reff use in lidar and set to 0 for radar
    129130  endif
    130   if ((gbx%use_reff) .and. (reff_zero)) then ! Inconsistent choice. Want to use Reff but not inputs passed
    131         print *, '---------- COSP ERROR ------------'
    132         print *, ''
    133         print *, 'use_reff==.true. but Reff is always zero'
    134         print *, ''
    135         print *, '----------------------------------'
    136         stop
    137   endif
     131!  if ((gbx%use_reff) .and. (reff_zero)) then ! Inconsistent choice. Want to use Reff but not inputs passed
     132!        print *, '---------- COSP ERROR ------------'
     133!        print *, ''
     134!        print *, 'use_reff==.true. but Reff is always zero'
     135!        print *, ''
     136!        print *, '----------------------------------'
     137!        stop
     138!  endif
    138139  if ((.not. gbx%use_reff) .and. (reff_zero)) then ! No Reff in radar. Default in lidar
    139140        gbx%Reff = DEFAULT_LIDAR_REFF
    … …  
    322323  integer :: Niter     ! Number of calls to cosp_simulator
    323324  integer :: i,j,k
     325  integer :: I_HYDRO
    324326  real,dimension(:,:),pointer :: column_frac_out ! Array with one column of frac_out
    325327  integer,parameter :: scops_debug=0    !  set to non-zero value to print out inputs for debugging in SCOPS
    … …  
    330332  real,dimension(:,:),allocatable :: frac_ls,prec_ls,frac_cv,prec_cv ! Cloud/Precipitation fraction in each model level
    331333                                                                     ! Levels are from SURFACE to TOA
     334  real,dimension(:,:),allocatable :: rho ! (Npoints, Nlevels). Atmospheric dens
    332335  type(cosp_sghydro) :: sghydro   ! Subgrid info for hydrometeors en each iteration
    333336
    … …  
    378381        do k=1,Nlevels,1
    379382            do i=1,Ncolumns,1
    380                 if (sgx%frac_out (j,i,Nlevels+1-k) .eq. 1) frac_ls(j,k)=frac_ls(j,k)+1.
    381                 if (sgx%frac_out (j,i,Nlevels+1-k) .eq. 2) frac_cv(j,k)=frac_cv(j,k)+1.
     383                if (sgx%frac_out (j,i,Nlevels+1-k) == I_LSC) frac_ls(j,k)=frac_ls(j,k)+1.
     384                if (sgx%frac_out (j,i,Nlevels+1-k) == I_CVC) frac_cv(j,k)=frac_cv(j,k)+1.
    382385                if (sgx%prec_frac(j,i,Nlevels+1-k) .eq. 1) prec_ls(j,k)=prec_ls(j,k)+1.
    383386                if (sgx%prec_frac(j,i,Nlevels+1-k) .eq. 2) prec_cv(j,k)=prec_cv(j,k)+1.
    … …  
    414417            !--------- Mixing ratios for clouds and Reff for Clouds and precip -------
    415418            column_frac_out => sgx%frac_out(:,k,:)
    416             where (column_frac_out == 1)     !+++++++++++ LS clouds ++++++++
     419            where (column_frac_out == I_LSC)     !+++++++++++ LS clouds ++++++++
    417420                sghydro%mr_hydro(:,k,:,I_LSCLIQ) = gbx%mr_hydro(:,:,I_LSCLIQ)
    418421                sghydro%mr_hydro(:,k,:,I_LSCICE) = gbx%mr_hydro(:,:,I_LSCICE)
    … …  
    423426                sghydro%Reff(:,k,:,I_LSSNOW)     = gbx%Reff(:,:,I_LSSNOW)
    424427                sghydro%Reff(:,k,:,I_LSGRPL)     = gbx%Reff(:,:,I_LSGRPL)
    425             elsewhere (column_frac_out == 2) !+++++++++++ CONV clouds ++++++++
     428            elsewhere (column_frac_out == I_CVC) !+++++++++++ CONV clouds ++++++++
    426429                sghydro%mr_hydro(:,k,:,I_CVCLIQ) = gbx%mr_hydro(:,:,I_CVCLIQ)
    427430                sghydro%mr_hydro(:,k,:,I_CVCICE) = gbx%mr_hydro(:,:,I_CVCICE)
    … …  
    434437            !--------- Precip -------
    435438            if (.not. gbx%use_precipitation_fluxes) then
    436                 where (column_frac_out == 1)  !+++++++++++ LS Precipitation ++++++++
     439                where (column_frac_out == I_LSC)  !+++++++++++ LS Precipitation ++++++++
    437440                    sghydro%mr_hydro(:,k,:,I_LSRAIN) = gbx%mr_hydro(:,:,I_LSRAIN)
    438441                    sghydro%mr_hydro(:,k,:,I_LSSNOW) = gbx%mr_hydro(:,:,I_LSSNOW)
    439442                    sghydro%mr_hydro(:,k,:,I_LSGRPL) = gbx%mr_hydro(:,:,I_LSGRPL)
    440                 elsewhere (column_frac_out == 2) !+++++++++++ CONV Precipitation ++++++++
     443                elsewhere (column_frac_out == I_CVC) !+++++++++++ CONV Precipitation ++++++++
    441444                    sghydro%mr_hydro(:,k,:,I_CVRAIN) = gbx%mr_hydro(:,:,I_CVRAIN)
    442445                    sghydro%mr_hydro(:,k,:,I_CVSNOW) = gbx%mr_hydro(:,:,I_CVSNOW)
    … …  
    488491                        sghydro%mr_hydro(:,:,:,I_LSRAIN),sghydro%mr_hydro(:,:,:,I_LSSNOW),sghydro%mr_hydro(:,:,:,I_LSGRPL), &
    489492                        sghydro%mr_hydro(:,:,:,I_CVRAIN),sghydro%mr_hydro(:,:,:,I_CVSNOW))
    490         endif
     493       endif
    491494   !++++++++++ CRM mode ++++++++++
    492495   else

  • LMDZ4/trunk/libf/cosp/cosp_constants.F90

    r1279 r1414  
    5151    ! Number of possible output variables
    5252    integer,parameter :: N_OUT_LIST = 27
     53    ! Value for forward model result from a level that is under the ground
     54    real,parameter :: R_GROUND = -1.0E20
     55
     56    ! Stratiform and convective clouds in frac_out
     57    integer, parameter :: I_LSC = 1, & ! Large-scale clouds
     58                          I_CVC = 2    ! Convective clouds
    5359   
    5460    !--- Radar constants
    … …  
    5662    integer,parameter :: DBZE_BINS     =   15   ! Number of dBZe bins in histogram (cfad)
    5763    real,parameter    :: DBZE_MIN      = -100.0 ! Minimum value for radar reflectivity
    58     real,parameter    :: DBZE_MAX      =   30.0 ! Maximum value for radar reflectivity
     64    real,parameter    :: DBZE_MAX      =   80.0 ! Maximum value for radar reflectivity
    5965    real,parameter    :: CFAD_ZE_MIN   =  -50.0 ! Lower value of the first CFAD Ze bin
    6066    real,parameter    :: CFAD_ZE_WIDTH =    5.0 ! Bin width (dBZe)
    … …  
    6975    integer,parameter :: LIDAR_NCAT    =   4
    7076    integer,parameter :: PARASOL_NREFL =   5 ! parasol
    71 !     real,parameter,dimension(PARASOL_NREFL) :: PARASOL_SZA = (/0.0, 15.0, 30.0, 45.0, 60.0/)
    72     real,parameter,dimension(PARASOL_NREFL) :: PARASOL_SZA = (/1.0, 2.0, 3.0, 4.0, 5.0/)
     77    real,parameter,dimension(PARASOL_NREFL) :: PARASOL_SZA = (/0.0, 20.0, 40.0, 6.0, 80.0/)
     78!    real,parameter,dimension(PARASOL_NREFL) :: PARASOL_SZA = (/1.0, 2.0, 3.0, 4.0, 5.0/)
    7379    real,parameter    :: DEFAULT_LIDAR_REFF = 30.0e-6 ! Default lidar effective radius
    7480   

  • LMDZ4/trunk/libf/cosp/cosp_isccp_simulator.F90

    r1279 r1414  
    8585     
    8686  ! Change boxptop from hPa to Pa. This avoids using UDUNITS in CMOR
    87   y%boxptop = y%boxptop*100.0
     87!  y%boxptop = y%boxptop*100.0
    8888 
    8989  ! Check if there is any value slightly greater than 1

  • LMDZ4/trunk/libf/cosp/cosp_simulator.F90

    r1279 r1414  
    6565  !+++++++++ Radar model ++++++++++ 
    6666  if (cfg%Lradar_sim) then
    67     call system_clock(t0,count_rate,count_max)
    6867    call cosp_radar(gbx,sgx,sghydro,sgradar)
    69     call system_clock(t1,count_rate,count_max)
    70     print *, '%%%%%%  Radar:', (t1-t0)*1.0/count_rate, ' s'
    71   else
    72     print *, '%%%%%%  Radar not used'
    7368  endif
    7469 
    7570  !+++++++++ Lidar model ++++++++++
    7671  if (cfg%Llidar_sim) then
    77     call system_clock(t0,count_rate,count_max)
    7872    call cosp_lidar(gbx,sgx,sghydro,sglidar)
    79     call system_clock(t1,count_rate,count_max)
    80     print *, '%%%%%%  Lidar:', (t1-t0)*1.0/count_rate, ' s'
    81   else
    82     print *, '%%%%%%  Lidar not used'
    8373  endif
    8474
    … …  
    8676  !+++++++++ ISCCP simulator ++++++++++
    8777  if (cfg%Lisccp_sim) then
    88     call system_clock(t0,count_rate,count_max)
    8978    call cosp_isccp_simulator(gbx,sgx,isccp)
    90     call system_clock(t1,count_rate,count_max)
    91     print *, '%%%%%%  ISCCP:', (t1-t0)*1.0/count_rate, ' s'
    92   else
    93     print *, '%%%%%%  ISCCP not used'
    9479  endif
    9580 
    9681  !+++++++++ MISR simulator ++++++++++
    9782  if (cfg%Lmisr_sim) then
    98     call system_clock(t0,count_rate,count_max)
    9983    call cosp_misr_simulator(gbx,sgx,misr)
    100     call system_clock(t1,count_rate,count_max)
    101     print *, '%%%%%%  MISR:', (t1-t0)*1.0/count_rate, ' s'
    102   else
    103     print *, '%%%%%%  MISR not used'
    10484  endif
    10585 
    10686
    10787  !+++++++++++ Summary statistics +++++++++++
    108 !   write(*,*) 'Stats:'
    109 !   read(*,*) c
    11088  if (cfg%Lstats) then
    111     call system_clock(t0,count_rate,count_max)
    11289    call cosp_stats(gbx,sgx,cfg,sgradar,sglidar,vgrid,stradar,stlidar)
    113     call system_clock(t1,count_rate,count_max)
    114     print *, '%%%%%%  Stats:', (t1-t0)*1.0/count_rate, ' s'
    115   endif
    116   !+++++++++++ change of units after computation of statistics +++++++++++
    117   if (cfg%Llidar_sim) then
    118     where((sglidar%beta_tot > 0.0) .and. (sglidar%beta_tot /= R_UNDEF))
    119         sglidar%beta_tot = log10(sglidar%beta_tot)
    120     elsewhere
    121         sglidar%beta_tot = R_UNDEF
    122     end where
     90!    print *, '%%%%%%  Stats:', (t1-t0)*1.0/count_rate, ' s'
    12391  endif
    12492
     93 
    12594END SUBROUTINE COSP_SIMULATOR
    12695

  • LMDZ4/trunk/libf/cosp/cosp_stats.F90

    r1396 r1414  
    11! (c) British Crown Copyright 2008, the Met Office.
    22! All rights reserved.
    3 ! 
    4 ! Redistribution and use in source and binary forms, with or without modification, are permitted 
     3!
     4! Redistribution and use in source and binary forms, with or without modification, are permitted
    55! provided that the following conditions are met:
    6 ! 
    7 !     * Redistributions of source code must retain the above copyright notice, this list 
     6!
     7!     * Redistributions of source code must retain the above copyright notice, this list
    88!       of conditions and the following disclaimer.
    99!     * Redistributions in binary form must reproduce the above copyright notice, this list
    10 !       of conditions and the following disclaimer in the documentation and/or other materials 
     10!       of conditions and the following disclaimer in the documentation and/or other materials
    1111!       provided with the distribution.
    12 !     * Neither the name of the Met Office nor the names of its contributors may be used 
    13 !       to endorse or promote products derived from this software without specific prior written 
     12!     * Neither the name of the Met Office nor the names of its contributors may be used
     13!       to endorse or promote products derived from this software without specific prior written
    1414!       permission.
    15 ! 
    16 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 
    17 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND 
    18 ! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 
    19 ! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
    20 ! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
    21 ! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER 
    22 ! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT 
     15!
     16! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
     17! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
     18! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
     19! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
     20! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
     21! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
     22! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
    2323! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
    2424
    … …  
    2929! Oct 2008 - J.-L. Dufresne   - Bug fixed. Assignment of Npoints,Nlevels,Nhydro,Ncolumns in COSP_STATS
    3030! Oct 2008 - H. Chepfer       - Added PARASOL reflectance arguments
    31 ! May 2010 - L. Fairhead      - Optimisation of COSP_CHANGE_VERTICAL_GRID routine for NEC SX8 computer
    32 !
    33 ! 
     31! Jun 2010 - T. Yokohata, T. Nishimura and K. Ogochi - Added NEC SXs optimisations
     32!
     33!
    3434MODULE MOD_COSP_STATS
    3535  USE MOD_COSP_CONSTANTS
    … …  
    4545!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    4646SUBROUTINE COSP_STATS(gbx,sgx,cfg,sgradar,sglidar,vgrid,stradar,stlidar)
    47  
     47
    4848   ! Input arguments
    4949   type(cosp_gridbox),intent(in) :: gbx
    … …  
    5555   ! Output arguments
    5656   type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics for radar
    57    type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics for lidar 
    58    
    59    ! Local variables 
     57   type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics for lidar
     58
     59   ! Local variables
    6060   integer :: Npoints  !# of grid points
    6161   integer :: Nlevels  !# of levels
    … …  
    6666   real,dimension(:,:,:),allocatable :: Ze_out,betatot_out,betamol_in,betamol_out,ph_in,ph_out
    6767   real,dimension(:,:),allocatable :: ph_c,betamol_c
    68  
     68
    6969   Npoints  = gbx%Npoints
    7070   Nlevels  = gbx%Nlevels
    … …  
    7373   Nlr      = vgrid%Nlvgrid
    7474
    75    if (cfg%Lcfad_lidarsr532) ok_lidar_cfad=.true.
     75   if (cfg%Lcfad_Lidarsr532) ok_lidar_cfad=.true.
    7676
    7777   if (vgrid%use_vgrid) then ! Statistics in a different vertical grid
    … …  
    8181        Ze_out = 0.0
    8282        betatot_out  = 0.0
    83         betamol_in(:,1,:) = sglidar%beta_mol(:,:)
    8483        betamol_out= 0.0
    8584        betamol_c  = 0.0
    … …  
    8988        !++++++++++++ Radar CFAD ++++++++++++++++
    9089        if (cfg%Lradar_sim) then
    91            call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sgradar%Ze_tot, &
     90            call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sgradar%Ze_tot, &
    9291                                           Nlr,vgrid%zl,vgrid%zu,Ze_out,log_units=.true.)
    9392            stradar%cfad_ze = cosp_cfad(Npoints,Ncolumns,Nlr,DBZE_BINS,Ze_out, &
    … …  
    9695        !++++++++++++ Lidar CFAD ++++++++++++++++
    9796        if (cfg%Llidar_sim) then
     97            betamol_in(:,1,:) = sglidar%beta_mol(:,:)
    9898            call cosp_change_vertical_grid(Npoints,1,Nlevels,gbx%zlev,gbx%zlev_half,betamol_in, &
    9999                                           Nlr,vgrid%zl,vgrid%zu,betamol_out)
    … …  
    114114        if (cfg%Lradar_sim.and.cfg%Llidar_sim) call cosp_lidar_only_cloud(Npoints,Ncolumns,Nlr, &
    115115                                    betatot_out,betamol_c,Ze_out, &
    116                                     stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)   
     116                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)
    117117        ! Deallocate arrays at coarse resolution
    118118        deallocate(Ze_out,betatot_out,betamol_in,betamol_out,betamol_c,ph_in,ph_out,ph_c)
    … …  
    131131        if (cfg%Lradar_sim.and.cfg%Llidar_sim) call cosp_lidar_only_cloud(Npoints,Ncolumns,Nlr, &
    132132                                    sglidar%beta_tot,sglidar%beta_mol,sgradar%Ze_tot, &
    133                                     stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)   
     133                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)
    134134   endif
    135135   ! Replace undef
    136    where (stlidar%cfad_sr   == LIDAR_UNDEF) stlidar%cfad_sr   = R_UNDEF 
    137    where (stlidar%lidarcld  == LIDAR_UNDEF) stlidar%lidarcld  = R_UNDEF 
    138    where (stlidar%cldlayer  == LIDAR_UNDEF) stlidar%cldlayer  = R_UNDEF 
    139    where (stlidar%parasolrefl == LIDAR_UNDEF) stlidar%parasolrefl = R_UNDEF 
     136   where (stlidar%cfad_sr   == LIDAR_UNDEF) stlidar%cfad_sr   = R_UNDEF
     137   where (stlidar%lidarcld  == LIDAR_UNDEF) stlidar%lidarcld  = R_UNDEF
     138   where (stlidar%cldlayer  == LIDAR_UNDEF) stlidar%cldlayer  = R_UNDEF
     139   where (stlidar%parasolrefl == LIDAR_UNDEF) stlidar%parasolrefl = R_UNDEF
    140140
    141141END SUBROUTINE COSP_STATS
    … …  
    162162
    163163   ! Local variables
    164    integer :: i,j,k,c
     164   integer :: i,j,k
    165165   logical :: lunits
    166    real :: ws(Npoints,M), ws_temp(Npoints,M)
    167    real,dimension(Npoints,Nlevels) :: xl, xu ! Lower and upper boundaries of model grid
    168    real,dimension(M) :: dz          ! Layer depth
    169    real,dimension(Npoints,Nlevels,M) :: w   ! Weights to do the mean at each point
    170    real,dimension(Npoints,Ncolumns,Nlevels) :: yp  ! Variable to be changed to a different grid.
    171                                            ! Local copy at a particular point.
    172                                            ! This allows for change of units.
    173  
     166
     167   integer :: l
     168   real,dimension(Npoints) :: ws,sumwyp
     169   real,dimension(Npoints,Nlevels) :: xl,xu
     170   real,dimension(Npoints,Nlevels) :: w
     171   real,dimension(Npoints,Ncolumns,Nlevels) :: yp
    174172
    175173   lunits=.false.
    176174   if (present(log_units)) lunits=log_units
    177175
    178    r(:,:,:) = 0.0
    179    yp(:,:,:) = y(:,:,:)
    180    w(:,:,:) = 0.0
    181    ws(:,:) = 0.0
    182    ws_temp(:,:) = 0.0
    183    dz = zu - zl
    184 
     176   r(:,:,:) = R_GROUND
     177   ! Vertical grid at that point
     178   xl(:,:) = zhalf(:,:)
     179   xu(:,1:Nlevels-1) = xl(:,2:Nlevels)
     180   xu(:,Nlevels) = zfull(:,Nlevels) +  zfull(:,Nlevels) - zhalf(:,Nlevels) ! Top level symmetric
     181   yp(:,:,:) = y(:,:,:) ! Temporary variable to regrid
    185182   ! Check for dBZ and change if necessary
    186183   if (lunits) then
    187184     where (y /= R_UNDEF)
    188         yp = 10.0**(y/10.0)
     185       yp = 10.0**(y/10.0)
     186     elsewhere
     187       yp = 0.0
    189188     end where
    190189   endif
    191 
    192    ! Vertical grids
    193    xl(:,:) = zhalf(:,:)
    194    xu(:,1:Nlevels-1) = zhalf(:,2:Nlevels)
    195    xu(:,Nlevels) = zfull(:,Nlevels) +  zfull(:,Nlevels) - zhalf(:,Nlevels) ! Top level symmetric
    196   ! Find weights
    197    do k=1, M
    198      do j=1, Nlevels
    199        do i=1, Npoints
     190   do k=1,M
     191     ! Find weights
     192     w(:,:) = 0.0
     193     do j=1,Nlevels
     194       do i=1,Npoints
    200195         if ((xl(i,j) < zl(k)).and.(xu(i,j) > zl(k)).and.(xu(i,j) <= zu(k))) then
    201196           !xl(j)-----------------xu(j)
    202197           !      zl(k)------------------------------zu(k)
    203            w(i,j,k) = xu(i,j) - zl(k)
     198           w(i,j) = xu(i,j) - zl(k)
    204199         else if ((xl(i,j) >= zl(k)).and.(xu(i,j) <= zu(k))) then
    205200           !           xl(j)-----------------xu(j)
    206201           !      zl(k)------------------------------zu(k)
    207            w(i,j,k) =  xl(i,j+1) - xl(i,j)
     202           w(i,j) = xu(i,j) - xl(i,j)
    208203         else if ((xl(i,j) >= zl(k)).and.(xl(i,j) < zu(k)).and.(xu(i,j) >= zu(k))) then
    209204           !                           xl(j)-----------------xu(j)
    210205           !      zl(k)------------------------------zu(k)
    211            w(i,j,k) = zu(k) - xl(i,j)
     206           w(i,j) = zu(k) - xl(i,j)
    212207         else if ((xl(i,j) <= zl(k)).and.(xu(i,j) >= zu(k))) then
    213208           !  xl(j)---------------------------xu(j)
    214209           !        zl(k)--------------zu(k)
    215            w(i,j,k) = dz(j)
     210           w(i,j) = zu(k) - zl(k)
     211         endif
     212       enddo
     213     enddo
     214     ! Do the weighted mean
     215     do j=1,Ncolumns
     216       ws    (:) = 0.0
     217       sumwyp(:) = 0.0
     218       do l=1,Nlevels
     219         do i=1,Npoints
     220           if (zu(k) > zhalf(i,1)) then ! Level above model bottom level
     221             ws    (i) = ws    (i) + w(i,l)
     222             sumwyp(i) = sumwyp(i) + w(i,l)*yp(i,j,l)
     223           endif
     224         enddo
     225       enddo
     226       do i=1,Npoints
     227         if (zu(k) > zhalf(i,1)) then ! Level above model bottom level
     228           if (ws(i) > 0.0) r(i,j,k) = sumwyp(i)/ws(i)
    216229         endif
    217230       enddo
    218231     enddo
    219232   enddo
    220 
    221    ! Do the weighted mean
    222    do k=1,M
    223      do i = 1, Npoints
    224         ws(i,k) = sum(w(i,:,k))
    225      enddo
    226    enddo
    227 
    228    ws_temp = 1.
    229    where (ws(:,:) > 0.0) ws_temp(:,:)=ws(:,:)
    230 
    231    do c=1,Ncolumns
    232      do k=1,M
    233        do i = 1, Npoints
    234           r(i,c,k) = sum(w(i,:,k)*yp(i,c,:))/ws_temp(i,k)
    235        enddo
    236      enddo
    237    enddo
    238 
    239    do k=1,M
    240      do i = 1, Npoints
    241         if (ws(i,k) <= 0.0) r(i,:,k)=0.0
    242      enddo
    243    enddo
    244  
    245233   ! Check for dBZ and change if necessary
    246234   if (lunits) then
    247       where (r(:,:,:) <= 0.0)
    248         r(:,:,:) = R_UNDEF
    249       elsewhere
    250         r(:,:,:) = 10.0*log10(r(:,:,:))
    251       end where
     235     do k=1,M
     236       do j=1,Ncolumns
     237         do i=1,Npoints
     238           if (zu(k) > zhalf(i,1)) then ! Level above model bottom level
     239             if (r(i,j,k) <= 0.0) then
     240                 r(i,j,k) = R_UNDEF
     241             else
     242                 r(i,j,k) = 10.0*log10(r(i,j,k))
     243             endif
     244           endif
     245         enddo
     246       enddo
     247     enddo
    252248   endif
    253    
    254 END SUBROUTINE COSP_CHANGE_VERTICAL_GRID
     249
     250
     251
     252END SUBROUTINE COSP_CHANGE_VERTICAL_GRID
    255253
    256254END MODULE MOD_COSP_STATS

  • LMDZ4/trunk/libf/cosp/cosp_utils.F90

    r1279 r1414  
    4040  END INTERFACE
    4141CONTAINS
     42
     43!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     44!------------------- SUBROUTINE COSP_PRECIP_MXRATIO --------------
     45!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     46SUBROUTINE COSP_PRECIP_MXRATIO(Npoints,Nlevels,Ncolumns,p,T,prec_frac,prec_type, &
     47                          n_ax,n_bx,alpha_x,c_x,d_x,g_x,a_x,b_x,gamma1,gamma2, &
     48                          flux,mxratio)
     49
     50    ! Input arguments, (IN)
     51    integer,intent(in) :: Npoints,Nlevels,Ncolumns
     52    real,intent(in),dimension(Npoints,Nlevels) :: p,T,flux
     53    real,intent(in),dimension(Npoints,Ncolumns,Nlevels) :: prec_frac
     54    real,intent(in) :: n_ax,n_bx,alpha_x,c_x,d_x,g_x,a_x,b_x,gamma1,gamma2,prec_type
     55    ! Input arguments, (OUT)
     56    real,intent(out),dimension(Npoints,Ncolumns,Nlevels) :: mxratio
     57    ! Local variables
     58    integer :: i,j,k
     59    real :: sigma,one_over_xip1,xi,rho0,rho
     60   
     61    mxratio = 0.0
     62
     63    if (n_ax >= 0.0) then ! N_ax is used to control which hydrometeors need to be computed
     64        !gamma1  = gamma(alpha_x + b_x + d_x + 1.0)
     65        !gamma2  = gamma(alpha_x + b_x + 1.0)
     66        xi      = d_x/(alpha_x + b_x - n_bx + 1.0)
     67        rho0    = 1.29
     68        sigma   = (gamma2/(gamma1*c_x))*(n_ax*a_x*gamma2)**xi
     69        one_over_xip1 = 1.0/(xi + 1.0)
     70       
     71        do k=1,Nlevels
     72            do j=1,Ncolumns
     73                do i=1,Npoints
     74                    if ((prec_frac(i,j,k)==prec_type).or.(prec_frac(i,j,k)==3.)) then
     75                        rho = p(i,k)/(287.05*T(i,k))
     76                        mxratio(i,j,k)=(flux(i,k)*((rho/rho0)**g_x)*sigma)**one_over_xip1
     77                        mxratio(i,j,k)=mxratio(i,j,k)/rho
     78                    endif
     79                enddo
     80            enddo
     81        enddo
     82    endif
     83END SUBROUTINE COSP_PRECIP_MXRATIO
    4284
    4385

  • LMDZ4/trunk/libf/cosp/icarus.F

    r1279 r1414  
    3333     &)
    3434
    35 !Id: icarus.f,v 4.0 2009/02/12 13:59:20 hadmw Exp $
     35!$Id: icarus.f,v 4.1 2010/05/27 16:30:18 hadmw Exp $
    3636
    3737! *****************************COPYRIGHT****************************
    … …  
    7272! *****************************COPYRIGHT*******************************
    7373! *****************************COPYRIGHT*******************************
     74! *****************************COPYRIGHT*******************************
    7475
    7576      implicit none
    … …  
    234235      REAL attropmin (npoints)
    235236      REAL atmax(npoints)
    236       REAL atmin(npoints)
    237237      REAL btcmin(npoints)
    238238      REAL transmax(npoints)
    … …  
    356356      do j=1,npoints
    357357          ptrop(j)=5000.
    358           atmin(j) = 400.
    359358          attropmin(j) = 400.
    360359          atmax(j) = 0.
    … …  
    373372                itrop(j)=ilev
    374373           end if
    375            if (at(j,ilev) .gt. atmax(j)) atmax(j)=at(j,ilev)
    376            if (at(j,ilev) .lt. atmin(j)) atmin(j)=at(j,ilev)
    377374        enddo
    37837512    continue
     376
     377      do 13 ilev=1,nlev
     378        do j=1,npoints
     379           if (at(j,ilev) .gt. atmax(j) .and.
     380     &              ilev  .ge. itrop(j)) atmax(j)=at(j,ilev)
     381        enddo
     38213    continue
    379383
    380384      end if

  • LMDZ4/trunk/libf/cosp/llnl_stats.F90

    r1279 r1414  
    2424
    2525MODULE MOD_LLNL_STATS
     26  USE MOD_COSP_CONSTANTS
    2627  IMPLICIT NONE
    2728
    … …  
    6263   do j = 1, Nlevels, 1
    6364      do k = 1, Ncolumns, 1
    64          do i = 1, Npoints, 1
    65             if ((x(i,k,j) >= xmin) .and. (x(i,k,j) <= xmax)) then
     65         do i = 1, Npoints, 1
     66            if (x(i,k,j) == R_GROUND) then
     67               cosp_cfad(i,:,j) = R_UNDEF
     68            elseif ((x(i,k,j) >= xmin) .and. (x(i,k,j) <= xmax)) then
    6669               ibin = ceiling((x(i,k,j) - bmin)/bwidth)
    6770               if (ibin > Nbins) ibin = Nbins
    … …  
    7275      enddo  !k
    7376   enddo  !j
    74    cosp_cfad = cosp_cfad / Ncolumns
     77   where ((cosp_cfad /= R_UNDEF).and.(cosp_cfad /= 0.0)) cosp_cfad = cosp_cfad / Ncolumns
    7578END FUNCTION COSP_CFAD
    7679
    … …  
    98101   integer :: pr,i,j
    99102   
    100 !    lidar_only_freq_cloud = 0.0
    101 !    tcc = 0.0
     103   lidar_only_freq_cloud = 0.0
     104   tcc = 0.0
    102105   do pr=1,Npoints
    103106     do i=1,Ncolumns

  • LMDZ4/trunk/libf/cosp/lmd_ipsl_stats.F90

    r1279 r1414  
    11! Copyright (c) 2009, Centre National de la Recherche Scientifique
    22! All rights reserved.
    3 ! 
    4 ! Redistribution and use in source and binary forms, with or without modification, are permitted 
     3!
     4! Redistribution and use in source and binary forms, with or without modification, are permitted
    55! provided that the following conditions are met:
    6 ! 
    7 !     * Redistributions of source code must retain the above copyright notice, this list 
     6!
     7!     * Redistributions of source code must retain the above copyright notice, this list
    88!       of conditions and the following disclaimer.
    99!     * Redistributions in binary form must reproduce the above copyright notice, this list
    10 !       of conditions and the following disclaimer in the documentation and/or other materials 
     10!       of conditions and the following disclaimer in the documentation and/or other materials
    1111!       provided with the distribution.
    1212!     * Neither the name of the LMD/IPSL/CNRS/UPMC nor the names of its
    13 !       contributors may be used to endorse or promote products derived from this software without 
     13!       contributors may be used to endorse or promote products derived from this software without
    1414!       specific prior written permission.
    15 ! 
    16 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 
    17 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND 
    18 ! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 
    19 ! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
    20 ! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
    21 ! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER 
    22 ! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT 
     15!
     16! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
     17! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
     18! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
     19! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
     20! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
     21! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
     22! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
    2323! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
    2424
    … …  
    3939! -----------------------------------------------------------------------------------
    4040! Lidar outputs :
    41 ! 
     41!
    4242! Diagnose cloud fraction (3D cloud fraction + low/middle/high/total cloud fraction
    4343! from the lidar signals (ATB and molecular ATB) computed from model outputs
    4444!      +
    4545! Compute CFADs of lidar scattering ratio SR and of depolarization index
    46 ! 
     46!
    4747! Authors: Sandrine Bony and Helene Chepfer (LMD/IPSL, CNRS, UPMC, France).
    4848!
    49 ! December 2008, S. Bony,  H. Chepfer and J-L. Dufresne : 
     49! December 2008, S. Bony,  H. Chepfer and J-L. Dufresne :
    5050! - change of the cloud detection threshold S_cld from 3 to 5, for better
    5151! with both day and night observations. The optical thinest clouds are missed.
    … …  
    5353! December 2008, A. Bodas-Salcedo:
    5454! - Dimensions of pmol reduced to (npoints,llm)
     55! August 2009, A. Bodas-Salcedo:
     56! - Warning message regarding PARASOL being valid only over ocean deleted.
     57! February 2010, A. Bodas-Salcedo:
     58! - Undef passed into cosp_cfad_sr
     59! June 2010, T. Yokohata, T. Nishimura and K. Ogochi
     60! Optimisation of COSP_CFAD_SR
    5561!
    5662! Version 1.0 (June 2007)
    … …  
    7076
    7177      real undef                    ! undefined value
    72       real pnorm(npoints,ncol,llm)  ! lidar ATB 
     78      real pnorm(npoints,ncol,llm)  ! lidar ATB
    7379      real pmol(npoints,llm)        ! molecular ATB
    74       real land(npoints)            ! Landmask [0 - Ocean, 1 - Land]   
     80      real land(npoints)            ! Landmask [0 - Ocean, 1 - Land]
    7581      real pplay(npoints,llm)       ! pressure on model levels (Pa)
    7682      logical ok_lidar_cfad         ! true if lidar CFAD diagnostics need to be computed
    … …  
    7884
    7985! c outputs :
    80       real lidarcld(npoints,llm)     ! 3D "lidar" cloud fraction 
     86      real lidarcld(npoints,llm)     ! 3D "lidar" cloud fraction
    8187      real cldlayer(npoints,ncat)    ! "lidar" cloud fraction (low, mid, high, total)
    82       real cfad2(npoints,max_bin,llm) ! CFADs of SR 
    83       real srbval(max_bin)           ! SR bins in CFADs 
     88      real cfad2(npoints,max_bin,llm) ! CFADs of SR
     89      real srbval(max_bin)           ! SR bins in CFADs
    8490      real parasolrefl(npoints,nrefl)! grid-averaged parasol reflectance
    8591
    8692! c threshold for cloud detection :
    87       real S_clr 
    88       parameter (S_clr = 1.2) 
     93      real S_clr
     94      parameter (S_clr = 1.2)
    8995      real S_cld
    9096!      parameter (S_cld = 3.0)  ! Previous thresold for cloud detection
    … …  
    102108! c 0- Initializations
    103109! c -------------------------------------------------------
    104 ! Parasol reflectance algorithm is not valid over land. Write
    105 ! a warning if there is no land. Landmask [0 - Ocean, 1 - Land]
    106       IF ( MAXVAL(land(:)) .EQ. 0.0) THEN
    107           WRITE (*,*) 'WARNING. PARASOL reflectance is not valid over land' &
    108              & ,' and there is only land'
    109       END IF
     110!
    110111
    111112!  Should be modified in future version
    … …  
    116117! c -------------------------------------------------------
    117118!
    118 !       where ((pnorm.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 )) 
     119!       where ((pnorm.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
    119120!          x3d = pnorm/pmol
    120121!       elsewhere
    … …  
    124125      do ic = 1, ncol
    125126        pnorm_c = pnorm(:,ic,:)
    126         where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 )) 
     127        where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
    127128            x3d_c = pnorm_c/pmol
    128129        elsewhere
    … …  
    142143
    143144! c -------------------------------------------------------
    144 ! c 3- CFADs 
     145! c 3- CFADs
    145146! c -------------------------------------------------------
    146147      if (ok_lidar_cfad) then
    … …  
    148149! c CFADs of subgrid-scale lidar scattering ratios :
    149150! c -------------------------------------------------------
    150       CALL COSP_CFAD_SR(npoints,ncol,llm,max_bin, &
     151      CALL COSP_CFAD_SR(npoints,ncol,llm,max_bin,undef, &
    151152                 x3d, &
    152153                 S_att,S_clr,xmax,cfad2,srbval)
    … …  
    172173! if land=0 -> parasolrefl=parasolrefl
    173174        parasolrefl(:,k) = parasolrefl(:,k) * MAX(1.0-land(:),0.0) &
    174                            + (1.0 - MAX(1.0-land(:),0.0))*undef 
     175                           + (1.0 - MAX(1.0-land(:),0.0))*undef
    175176      enddo
    176177
    177178      RETURN
    178179      END SUBROUTINE diag_lidar
    179          
    180          
     180
     181
    181182!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    182183!-------------------- FUNCTION COSP_CFAD_SR ------------------------
    183184! Author: Sandrine Bony (LMD/IPSL, CNRS, Paris)
    184185!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    185       SUBROUTINE COSP_CFAD_SR(Npoints,Ncolumns,Nlevels,Nbins, &
     186      SUBROUTINE COSP_CFAD_SR(Npoints,Ncolumns,Nlevels,Nbins,undef, &
    186187                      x,S_att,S_clr,xmax,cfad,srbval)
    187188      IMPLICIT NONE
    … …  
    205206! Input arguments
    206207      integer Npoints,Ncolumns,Nlevels,Nbins
    207       real xmax,S_att,S_clr
     208      real xmax,S_att,S_clr,undef
    208209! Input-output arguments
    209210      real x(Npoints,Ncolumns,Nlevels)
    … …  
    213214! Local variables
    214215      integer i, j, k, ib
     216      real srbval_ext(0:Nbins)
    215217
    216218! c -------------------------------------------------------
    … …  
    235237      cfad(:,:,:) = 0.0
    236238
    237 
     239      srbval_ext(1:Nbins) = srbval
     240      srbval_ext(0) = -1.0
    238241! c -------------------------------------------------------
    239242! c c- Compute CFAD
    240243! c -------------------------------------------------------
    241244
    242         do j = Nlevels, 1, -1
    243           do k = 1, Ncolumns
    244               where ( x(:,k,j).le.srbval(1) ) &
    245                         cfad(:,1,j) = cfad(:,1,j) + 1.0
    246           enddo  !k
    247         enddo  !j
    248 
    249       do ib = 2, Nbins
    250         do j = Nlevels, 1, -1
    251           do k = 1, Ncolumns
    252               where ( ( x(:,k,j).gt.srbval(ib-1) ) .and. ( x(:,k,j).le.srbval(ib) ) ) &
    253                         cfad(:,ib,j) = cfad(:,ib,j) + 1.0
    254           enddo  !k
    255         enddo  !j
    256       enddo  !ib
    257 
    258       cfad(:,:,:) = cfad(:,:,:) / float(Ncolumns)
     245      do j = 1, Nlevels
     246         do ib = 1, Nbins
     247            do k = 1, Ncolumns
     248               do i = 1, Npoints
     249                  if (x(i,k,j) /= undef) then
     250                     if ((x(i,k,j).gt.srbval_ext(ib-1)).and.(x(i,k,j).le.srbval_ext(ib))) &
     251                          cfad(i,ib,j) = cfad(i,ib,j) + 1.0
     252                  else
     253                     cfad(i,ib,j) = undef
     254                  endif
     255               enddo
     256            enddo
     257         enddo
     258      enddo
     259
     260      where (cfad .ne. undef)  cfad = cfad / float(Ncolumns)
    259261
    260262! c -------------------------------------------------------
    … …  
    264266!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    265267!-------------------- SUBROUTINE COSP_CLDFRAC -------------------
    266 ! c Purpose: Cloud fraction diagnosed from lidar measurements 
     268! c Purpose: Cloud fraction diagnosed from lidar measurements
    267269!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    268270      SUBROUTINE COSP_CLDFRAC(Npoints,Ncolumns,Nlevels,Ncat, &
    … …  
    288290      real nsub(Npoints,Nlevels)
    289291
    290 
    291 ! ---------------------------------------------------------------
    292 ! 1- initialization
     292      real cldlay1(Npoints,Ncolumns)
     293      real cldlay2(Npoints,Ncolumns)
     294      real cldlay3(Npoints,Ncolumns)
     295      real nsublay1(Npoints,Ncolumns)
     296      real nsublay2(Npoints,Ncolumns)
     297      real nsublay3(Npoints,Ncolumns)
     298
     299
     300! ---------------------------------------------------------------
     301! 1- initialization
    293302! ---------------------------------------------------------------
    294303
    … …  
    317326
    318327! number of usefull sub-columns:
    319          where ( (x(:,:,k).gt.S_att) .and. (x(:,:,k).ne. undef)  ) 
     328         where ( (x(:,:,k).gt.S_att) .and. (x(:,:,k).ne. undef)  )
    320329           srok(:,:,k)=1.0
    321330         elsewhere
    … …  
    329338! categories) :
    330339! ---------------------------------------------------------------
    331 ! Test abderr
     340      lidarcld = 0.0
     341      nsub = 0.0
     342
     343!! XXX: Use cldlay[1-3] and nsublay[1-3] to avoid bank-conflicts.
     344      cldlay1 = 0.0
     345      cldlay2 = 0.0
     346      cldlay3 = 0.0
     347      cldlay(:,:,4) = 0.0 !! XXX: Ncat == 4
     348      nsublay1 = 0.0
     349      nsublay2 = 0.0
     350      nsublay3 = 0.0
     351      nsublay(:,:,4) = 0.0
    332352      do k = Nlevels, 1, -1
    333353       do ic = 1, Ncolumns
    334354        do ip = 1, Npoints
    335           iz=1
    336           p1 = pplay(ip,k)
    337           if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds
    338             iz=3
    339           else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then  ! mid clouds
    340             iz=2
     355         p1 = pplay(ip,k)
     356
     357         if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds
     358            cldlay3(ip,ic) = MAX(cldlay3(ip,ic), cldy(ip,ic,k))
     359            nsublay3(ip,ic) = MAX(nsublay3(ip,ic), srok(ip,ic,k))
     360         else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then  ! mid clouds
     361            cldlay2(ip,ic) = MAX(cldlay2(ip,ic), cldy(ip,ic,k))
     362            nsublay2(ip,ic) = MAX(nsublay2(ip,ic), srok(ip,ic,k))
     363         else
     364            cldlay1(ip,ic) = MAX(cldlay1(ip,ic), cldy(ip,ic,k))
     365            nsublay1(ip,ic) = MAX(nsublay1(ip,ic), srok(ip,ic,k))
    341366         endif
    342367
    343          cldlay(ip,ic,iz) = MAX(cldlay(ip,ic,iz),cldy(ip,ic,k))
    344          cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4),cldy(ip,ic,k))
     368         cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4), cldy(ip,ic,k))
    345369         lidarcld(ip,k)=lidarcld(ip,k) + cldy(ip,ic,k)
    346 
    347          nsublay(ip,ic,iz) = MAX(nsublay(ip,ic,iz),srok(ip,ic,k))
    348370         nsublay(ip,ic,4) = MAX(nsublay(ip,ic,4),srok(ip,ic,k))
    349371         nsub(ip,k)=nsub(ip,k) + srok(ip,ic,k)
    350 
    351372        enddo
    352373       enddo
    353374      enddo
     375      cldlay(:,:,1) = cldlay1
     376      cldlay(:,:,2) = cldlay2
     377      cldlay(:,:,3) = cldlay3
     378      nsublay(:,:,1) = nsublay1
     379      nsublay(:,:,2) = nsublay2
     380      nsublay(:,:,3) = nsublay3
    354381
    355382! -- grid-box 3D cloud fraction
    … …  
    369396       do ic = 1, Ncolumns
    370397
    371           cldlayer(:,iz)=cldlayer(:,iz) + cldlay(:,ic,iz)   
    372           nsublayer(:,iz)=nsublayer(:,iz) + nsublay(:,ic,iz) 
     398          cldlayer(:,iz)=cldlayer(:,iz) + cldlay(:,ic,iz)
     399          nsublayer(:,iz)=nsublayer(:,iz) + nsublay(:,ic,iz)
    373400
    374401       enddo
    … …  
    383410      END SUBROUTINE COSP_CLDFRAC
    384411! ---------------------------------------------------------------
    385          
     412
    386413END MODULE MOD_LMD_IPSL_STATS
Note: See TracChangeset for help on using the changeset viewer.