source: LMDZ6/branches/Amaury_dev/libf/phylmd/cosp/mod_cosp_radar.F90 @ 5136

! (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.

MODULE MOD_COSP_RADAR
  USE MOD_COSP_CONSTANTS
  USE MOD_COSP_TYPES
  USE MOD_COSP_UTILS
  use radar_simulator_types
  use array_lib
  use atmos_lib
  use format_input
  IMPLICIT NONE

  INTERFACE
    subroutine radar_simulator(hp,nprof,ngate,undef, &
        hgt_matrix,hm_matrix,re_matrix,Np_matrix, &
        p_matrix,t_matrix,rh_matrix, &
        Ze_non,Ze_ray,g_to_vol,a_to_vol,dBZe, &
        g_to_vol_in,g_to_vol_out)

        use m_mrgrnk
        use array_lib
        use math_lib
        use optics_lib
        use radar_simulator_types
        implicit none

        ! ----- INPUTS -----  
        type(class_param) :: hp

        integer, intent(in) :: nprof,ngate

        real undef
        real*8, dimension(nprof,ngate), intent(in) :: hgt_matrix, p_matrix, &
            t_matrix,rh_matrix
        real*8, dimension(hp%nhclass,nprof,ngate), intent(in) :: hm_matrix
        real*8, dimension(hp%nhclass,nprof,ngate), intent(inout) :: re_matrix
        real*8, dimension(hp%nhclass,nprof,ngate), intent(inout) :: Np_matrix

        ! ----- OUTPUTS -----
        real*8, dimension(nprof,ngate), intent(out) :: Ze_non,Ze_ray, &
            g_to_vol,dBZe,a_to_vol
        ! ----- OPTIONAL -----
        real*8, optional, dimension(nprof,ngate) :: &
            g_to_vol_in,g_to_vol_out
     end subroutine radar_simulator
  END INTERFACE

CONTAINS

!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!------------------- SUBROUTINE COSP_RADAR ------------------------
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
SUBROUTINE COSP_RADAR(gbx,sgx,sghydro,z)
  IMPLICIT NONE

  ! Arguments
  type(cosp_gridbox),intent(inout) :: gbx  ! Gridbox info
  type(cosp_subgrid),intent(in) :: sgx  ! Subgrid info
  type(cosp_sghydro),intent(in) :: sghydro  ! Subgrid info for hydrometeors
  type(cosp_sgradar),intent(inout) :: z ! Output from simulator, subgrid

  ! Local variables 
  integer :: & 
  nsizes            ! num of discrete drop sizes

  real*8, dimension(:,:), allocatable :: &
  g_to_vol ! integrated atten due to gases, r>v (dB)

  real*8, dimension(:,:), allocatable :: &
  Ze_non, &         ! radar reflectivity withOUT attenuation (dBZ)
  Ze_ray, &         ! Rayleigh reflectivity (dBZ)
  h_atten_to_vol, &     ! attenuation by hydromets, radar to vol (dB)
  g_atten_to_vol, &     ! gaseous atteunation, radar to vol (dB)
  dBZe, &           ! effective radar reflectivity factor (dBZ)
  hgt_matrix, &         ! height of hydrometeors (km)
  t_matrix, &                   !temperature (k)
  p_matrix, &                   !pressure (hPa)
  rh_matrix                     !relative humidity (%)

  real*8, dimension(:,:,:), allocatable :: &
  hm_matrix, &          ! hydrometeor mixing ratio (g kg^-1)
  re_matrix, &          ! effective radius (microns).   Optional. 0 ==> use Np_matrix or defaults
  Np_matrix         ! total number concentration (kg^-1).   Optional 0==> use defaults 

  integer, parameter :: one = 1
  ! logical :: hgt_reversed
  logical :: hgt_descending
  integer :: pr,i,j,k,unt,ngate

! ----- main program settings ------

  ! Inputs to Quickbeam
  allocate(hgt_matrix(gbx%Npoints,gbx%Nlevels),p_matrix(gbx%Npoints,gbx%Nlevels), &
           t_matrix(gbx%Npoints,gbx%Nlevels),rh_matrix(gbx%Npoints,gbx%Nlevels))
  allocate(hm_matrix(gbx%Nhydro,gbx%Npoints,gbx%Nlevels))
  allocate(re_matrix(gbx%Nhydro,gbx%Npoints,gbx%Nlevels))
  allocate(Np_matrix(gbx%Nhydro,gbx%Npoints,gbx%Nlevels))

  ! Outputs from Quickbeam
  allocate(Ze_non(gbx%Npoints,gbx%Nlevels))
  allocate(Ze_ray(gbx%Npoints,gbx%Nlevels))
  allocate(h_atten_to_vol(gbx%Npoints,gbx%Nlevels))
  allocate(g_atten_to_vol(gbx%Npoints,gbx%Nlevels))
  allocate(dBZe(gbx%Npoints,gbx%Nlevels))

  ! Optional argument. It is computed and returned in the first call to
  ! radar_simulator, and passed as input in the rest
  allocate(g_to_vol(gbx%Npoints,gbx%Nlevels))

  ! Even if there is no unit conversion, they are needed for type conversion
  p_matrix   = gbx%p/100.0     ! From Pa to hPa
  hgt_matrix = gbx%zlev/1000.0 ! From m to km
  t_matrix   = gbx%T
  rh_matrix  = gbx%q
  re_matrix  = 0.0


  ! set flag denoting position of radar relative to hgt_matrix orientation
          ngate = size(hgt_matrix,2)

          hgt_descending = hgt_matrix(1,1) > hgt_matrix(1,ngate)

          if ( &
             (gbx%surface_radar == 1 .and. hgt_descending) .or.  &
             (gbx%surface_radar == 0 .and. (.not. hgt_descending)) &
             ) &
          then
            gbx%hp%radar_at_layer_one = .false.
          else
            gbx%hp%radar_at_layer_one = .true.
          endif

  ! ----- loop over subcolumns -----
  do pr=1,sgx%Ncolumns

      !  NOTE:
      !  atmospheric profiles are the same within the same gridbox
      !  only hydrometeor profiles will be different for each subgridbox

         do i=1,gbx%Nhydro
            hm_matrix(i,:,:) = sghydro%mr_hydro(:,pr,:,i)*1000.0 ! Units from kg/kg to g/kg
            if (gbx%use_reff) then
              re_matrix(i,:,:) = sghydro%Reff(:,pr,:,i)*1.e6       ! Units from m to micron
              Np_matrix(i,:,:) = sghydro%Np(:,pr,:,i)              ! Units [#/kg]
            endif
         enddo

      !   ----- call radar simulator -----
      if (pr == 1) then ! Compute gaseous attenuation for all profiles
         call radar_simulator(gbx%hp,gbx%Npoints,gbx%Nlevels,R_UNDEF, &
           hgt_matrix,hm_matrix,re_matrix,Np_matrix, &
           p_matrix,t_matrix,rh_matrix, &
           Ze_non,Ze_ray,h_atten_to_vol,g_atten_to_vol,dBZe,g_to_vol_out=g_to_vol)
      else ! Use gaseous atteunuation for pr = 1
         call radar_simulator(gbx%hp,gbx%Npoints,gbx%Nlevels,R_UNDEF, &
           hgt_matrix,hm_matrix,re_matrix,Np_matrix, &
           p_matrix,t_matrix,rh_matrix, &
           Ze_non,Ze_ray,h_atten_to_vol,g_atten_to_vol,dBZe,g_to_vol_in=g_to_vol)
      endif

      ! store caluculated dBZe values for later output/processing
      z%Ze_tot(:,pr,:)=dBZe(:,:)
  enddo !pr

  deallocate(hgt_matrix,p_matrix,t_matrix,rh_matrix)
  deallocate(hm_matrix,re_matrix, &
      Ze_non,Ze_ray,h_atten_to_vol,g_atten_to_vol,dBZe)
  deallocate(g_to_vol)
END SUBROUTINE COSP_RADAR

END MODULE MOD_COSP_RADAR