cosp_stats.F90

Go to the documentation of this file.

! (c) British Crown Copyright 2008, 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.

!
! History:
! Jul 2007 - A. Bodas-Salcedo - Initial version
! Jul 2008 - A. Bodas-Salcedo - Added capability of producing outputs in standard grid
! Oct 2008 - J.-L. Dufresne   - Bug fixed. Assignment of Npoints,Nlevels,Nhydro,Ncolumns in COSP_STATS
! Oct 2008 - H. Chepfer       - Added PARASOL reflectance arguments
! Jun 2010 - T. Yokohata, T. Nishimura and K. Ogochi - Added NEC SXs optimisations
!
!
MODULE mod_cosp_stats
  USE mod_cosp_constants
  USE mod_cosp_types
  USE mod_llnl_stats
  USE mod_lmd_ipsl_stats
  IMPLICIT NONE

CONTAINS

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!------------------- SUBROUTINE COSP_STATS ------------------------
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SUBROUTINE cosp_stats(gbx,sgx,cfg,sgradar,sglidar,vgrid,stradar,stlidar)

   ! Input arguments
   type(cosp_gridbox),intent(in) :: gbx
   type(cosp_subgrid),intent(in) :: sgx
   type(cosp_config),intent(in)  :: cfg
   type(cosp_sgradar),intent(in) :: sgradar
   type(cosp_sglidar),intent(in) :: sglidar
   type(cosp_vgrid),intent(in)   :: vgrid
   ! Output arguments
   type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics for radar
   type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics for lidar

   ! Local variables
   integer :: Npoints  !# of grid points
   integer :: Nlevels  !# of levels
   integer :: Nhydro   !# of hydrometeors
   integer :: Ncolumns !# of columns
   integer :: Nlr
   logical :: ok_lidar_cfad = .false.
   real,dimension(:,:,:),allocatable :: Ze_out,betatot_out,betamol_in,betamol_out,ph_in,ph_out
   real,dimension(:,:),allocatable :: ph_c,betamol_c

   npoints  = gbx%Npoints
   nlevels  = gbx%Nlevels
   nhydro   = gbx%Nhydro
   ncolumns = gbx%Ncolumns
   nlr      = vgrid%Nlvgrid

   if (cfg%Lcfad_Lidarsr532) ok_lidar_cfad=.true.

   if (vgrid%use_vgrid) then ! Statistics in a different vertical grid
        allocate(ze_out(npoints,ncolumns,nlr),betatot_out(npoints,ncolumns,nlr), &
                 betamol_in(npoints,1,nlevels),betamol_out(npoints,1,nlr),betamol_c(npoints,nlr), &
                 ph_in(npoints,1,nlevels),ph_out(npoints,1,nlr),ph_c(npoints,nlr))
        ze_out = 0.0
        betatot_out  = 0.0
        betamol_out= 0.0
        betamol_c  = 0.0
        ph_in(:,1,:)  = gbx%ph(:,:)
        ph_out  = 0.0
        ph_c    = 0.0
        !++++++++++++ Radar CFAD ++++++++++++++++
        if (cfg%Lradar_sim) then
            call cosp_change_vertical_grid(npoints,ncolumns,nlevels,gbx%zlev,gbx%zlev_half,sgradar%Ze_tot, &
                                           nlr,vgrid%zl,vgrid%zu,ze_out,log_units=.true.)
            stradar%cfad_ze = cosp_cfad(npoints,ncolumns,nlr,dbze_bins,ze_out, &
                                        dbze_min,dbze_max,cfad_ze_min,cfad_ze_width)
        endif
        !++++++++++++ Lidar CFAD ++++++++++++++++
        if (cfg%Llidar_sim) then
            betamol_in(:,1,:) = sglidar%beta_mol(:,:)
            call cosp_change_vertical_grid(npoints,1,nlevels,gbx%zlev,gbx%zlev_half,betamol_in, &
                                           nlr,vgrid%zl,vgrid%zu,betamol_out)
            call cosp_change_vertical_grid(npoints,ncolumns,nlevels,gbx%zlev,gbx%zlev_half,sglidar%beta_tot, &
                                           nlr,vgrid%zl,vgrid%zu,betatot_out)
            call cosp_change_vertical_grid(npoints,1,nlevels,gbx%zlev,gbx%zlev_half,ph_in, &
                                           nlr,vgrid%zl,vgrid%zu,ph_out)
            ph_c(:,:) = ph_out(:,1,:)
            betamol_c(:,:) = betamol_out(:,1,:)
            ! Stats from lidar_stat_summary
            call diag_lidar(npoints,ncolumns,nlr,sr_bins,parasol_nrefl &
                            ,betatot_out,betamol_c,sglidar%refl,gbx%land,ph_c &
                            ,lidar_undef,ok_lidar_cfad &
                            ,stlidar%cfad_sr,stlidar%srbval &
                            ,lidar_ncat,stlidar%lidarcld,stlidar%cldlayer,stlidar%parasolrefl)
        endif
        !++++++++++++ Lidar-only cloud amount and lidar&radar total cloud mount ++++++++++++++++
        if (cfg%Lradar_sim.and.cfg%Llidar_sim) call cosp_lidar_only_cloud(npoints,ncolumns,nlr, &
                                    betatot_out,betamol_c,ze_out, &
                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)
        ! Deallocate arrays at coarse resolution
        deallocate(ze_out,betatot_out,betamol_in,betamol_out,betamol_c,ph_in,ph_out,ph_c)
   else ! Statistics in model levels
        !++++++++++++ Radar CFAD ++++++++++++++++
        if (cfg%Lradar_sim) stradar%cfad_ze = cosp_cfad(npoints,ncolumns,nlr,dbze_bins,sgradar%Ze_tot, &
                                        dbze_min,dbze_max,cfad_ze_min,cfad_ze_width)
        !++++++++++++ Lidar CFAD ++++++++++++++++
        ! Stats from lidar_stat_summary
        if (cfg%Llidar_sim) call diag_lidar(npoints,ncolumns,nlr,sr_bins,parasol_nrefl &
                        ,sglidar%beta_tot,sglidar%beta_mol,sglidar%refl,gbx%land,gbx%ph &
                        ,lidar_undef,ok_lidar_cfad &
                        ,stlidar%cfad_sr,stlidar%srbval &
                        ,lidar_ncat,stlidar%lidarcld,stlidar%cldlayer,stlidar%parasolrefl)
        !++++++++++++ Lidar-only cloud amount and lidar&radar total cloud mount ++++++++++++++++
        if (cfg%Lradar_sim.and.cfg%Llidar_sim) call cosp_lidar_only_cloud(npoints,ncolumns,nlr, &
                                    sglidar%beta_tot,sglidar%beta_mol,sgradar%Ze_tot, &
                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)
   endif
   ! Replace undef
   where (stlidar%cfad_sr   == lidar_undef) stlidar%cfad_sr   = r_undef
   where (stlidar%lidarcld  == lidar_undef) stlidar%lidarcld  = r_undef
   where (stlidar%cldlayer  == lidar_undef) stlidar%cldlayer  = r_undef
   where (stlidar%parasolrefl == lidar_undef) stlidar%parasolrefl = r_undef

END SUBROUTINE cosp_stats


!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!---------- SUBROUTINE COSP_CHANGE_VERTICAL_GRID ----------------
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SUBROUTINE cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,zfull,zhalf,y,M,zl,zu,r,log_units)
   implicit none
   ! Input arguments
   integer,intent(in) :: Npoints  !# of grid points
   integer,intent(in) :: Nlevels  !# of levels
   integer,intent(in) :: Ncolumns !# of columns
   real,dimension(Npoints,Nlevels),intent(in) :: zfull ! Height at model levels [m] (Bottom of model layer)
   real,dimension(Npoints,Nlevels),intent(in) :: zhalf ! Height at half model levels [m] (Bottom of model layer)
   real,dimension(Npoints,Ncolumns,Nlevels),intent(in) :: y     ! Variable to be changed to a different grid
   integer,intent(in) :: M  !# levels in the new grid
   real,dimension(M),intent(in) :: zl ! Lower boundary of new levels  [m]
   real,dimension(M),intent(in) :: zu ! Upper boundary of new levels  [m]
   logical,optional,intent(in) :: log_units ! log units, need to convert to linear units
   ! Output
   real,dimension(Npoints,Ncolumns,M),intent(out) :: r ! Variable on new grid

   ! Local variables
   integer :: i,j,k
   logical :: lunits

   integer :: l
   real,dimension(Npoints) :: ws,sumwyp
   real,dimension(Npoints,Nlevels) :: xl,xu
   real,dimension(Npoints,Nlevels) :: w
   real,dimension(Npoints,Ncolumns,Nlevels) :: yp

   lunits=.false.
   if (present(log_units)) lunits=log_units

   r(:,:,:) = r_ground
   ! Vertical grid at that point
   xl(:,:) = zhalf(:,:)
   xu(:,1:nlevels-1) = xl(:,2:nlevels)
   xu(:,nlevels) = zfull(:,nlevels) +  zfull(:,nlevels) - zhalf(:,nlevels) ! Top level symmetric
   yp(:,:,:) = y(:,:,:) ! Temporary variable to regrid
   ! Check for dBZ and change if necessary
   if (lunits) then
     where (y /= r_undef)
       yp = 10.0**(y/10.0)
     elsewhere
       yp = 0.0
     end where
   endif
   do k=1,m
     ! Find weights
     w(:,:) = 0.0
     do j=1,nlevels
       do i=1,npoints
         if ((xl(i,j) < zl(k)).and.(xu(i,j) > zl(k)).and.(xu(i,j) <= zu(k))) then
           !xl(j)-----------------xu(j)
           !      zl(k)------------------------------zu(k)
           w(i,j) = xu(i,j) - zl(k)
         else if ((xl(i,j) >= zl(k)).and.(xu(i,j) <= zu(k))) then
           !           xl(j)-----------------xu(j)
           !      zl(k)------------------------------zu(k)
           w(i,j) = xu(i,j) - xl(i,j)
         else if ((xl(i,j) >= zl(k)).and.(xl(i,j) < zu(k)).and.(xu(i,j) >= zu(k))) then
           !                           xl(j)-----------------xu(j)
           !      zl(k)------------------------------zu(k)
           w(i,j) = zu(k) - xl(i,j)
         else if ((xl(i,j) <= zl(k)).and.(xu(i,j) >= zu(k))) then
           !  xl(j)---------------------------xu(j)
           !        zl(k)--------------zu(k)
           w(i,j) = zu(k) - zl(k)
         endif
       enddo
     enddo
     ! Do the weighted mean
     do j=1,ncolumns
       ws(:) = 0.0
       sumwyp(:) = 0.0
       do l=1,nlevels
         do i=1,npoints
           if (zu(k) > zhalf(i,1)) then ! Level above model bottom level
             ws(i) = ws(i) + w(i,l)
             sumwyp(i) = sumwyp(i) + w(i,l)*yp(i,j,l)
           endif
         enddo
       enddo
       do i=1,npoints
         if (zu(k) > zhalf(i,1)) then ! Level above model bottom level
           if (ws(i) > 0.0) r(i,j,k) = sumwyp(i)/ws(i)
         endif
       enddo
     enddo
   enddo
   ! Check for dBZ and change if necessary
   if (lunits) then
     do k=1,m
       do j=1,ncolumns
         do i=1,npoints
           if (zu(k) > zhalf(i,1)) then ! Level above model bottom level
             if (r(i,j,k) <= 0.0) then
                 r(i,j,k) = r_undef
             else
                 r(i,j,k) = 10.0*log10(r(i,j,k))
             endif
           endif
         enddo
       enddo
     enddo
   endif



END SUBROUTINE cosp_change_vertical_grid

END MODULE mod_cosp_stats