mod_cnv2d.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module mod_cnv2d
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_atmos_grid_cartesc_index
  use scale_tracer
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: cnv2d_setup
  public :: cnv2d_tile_init
  public :: cnv2d_grads_init
  public :: cnv2d_exec
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: cnv2d_init
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  integer, parameter :: I_TILE  = 1
  integer, parameter :: I_GrADS = 2
  integer :: CNV2D_ftype
 
  integer  :: GLOBAL_IA
 
  integer               :: nLON, nLAT
  real(RP), allocatable :: data_org(:,:)
  real(RP), allocatable :: LAT_org (:,:)
  real(RP), allocatable :: LON_org (:,:)
  real(RP), allocatable :: LAT_1d  (:)
  real(RP), allocatable :: LON_1d  (:)
 
  ! interpolation
  integer               :: itp_lev
  integer,  allocatable :: idx_i(:,:,:)
  integer,  allocatable :: idx_j(:,:,:)
  real(RP), allocatable :: hfact(:,:,:)
  logical :: zonal, pole
 
  logical :: first = .true.
 
  ! TILE data
  character(len=H_SHORT) :: CNV2D_tile_dtype
  character(len=H_LONG)  :: CNV2D_tile_dir
  real(RP)               :: TILE_DLAT, TILE_DLON
  integer                :: TILE_nlim
  integer                :: TILE_nmax
  character(len=H_LONG), allocatable :: TILE_fname(:)
  logical,               allocatable :: TILE_hit  (:)
  integer,               allocatable :: TILE_JS   (:)
  integer,               allocatable :: TILE_JE   (:)
  integer,               allocatable :: TILE_IS   (:)
  integer,               allocatable :: TILE_IE   (:)
 
  integer :: dom_is, dom_ie, dom_js, dom_je
 
 
  ! GrADS data
  integer :: CNV2D_grads_fid
  integer :: CNV2D_grads_vid
 
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine cnv2d_setup
    use scale_prc, only: &
       prc_abort
    implicit none
 
    !---------------------------------------------------------------------------
 
    log_newline
    log_info("CNV2D_setup",*) 'Setup'
 
    !--- read namelist
!!$    rewind(IO_FID_CONF)
!!$    read(IO_FID_CONF,nml=PARAM_CNV2D,iostat=ierr)
!!$    if( ierr < 0 ) then !--- missing
!!$       LOG_INFO("CNV2D_setup",*) 'Not found namelist. Default used.'
!!$    elseif( ierr > 0 ) then !--- fatal error
!!$       LOG_ERROR("CNV2D_setup",*) 'Not appropriate names in namelist PARAM_CNV2D. Check!'
!!$       call PRC_abort
!!$    endif
!!$    LOG_NML(PARAM_CNV2D)
 
 
    return
  end subroutine cnv2d_setup
 
  subroutine cnv2d_tile_init( &
       dtype,        &
       dlat, dlon,   &
       dir,          &
       catalogue,    &
       interp_type,  &
       interp_level, &
       nmax          )
    use scale_prc, only: &
       prc_abort
    use scale_const, only: &
       undef  => const_undef, &
       d2r    => const_d2r
    use scale_file_tiledata, only: &
       file_tiledata_get_info,   &
       file_tiledata_get_latlon, &
       file_tiledata_get_data
    use scale_atmos_grid_cartesc_real, only: &
       latxv => atmos_grid_cartesc_real_latxv, &
       lonuy => atmos_grid_cartesc_real_lonuy
    implicit none
    character(len=*), intent(in) :: dtype
    real(rp),         intent(in) :: dlat, dlon
    character(len=*), intent(in) :: dir
    character(len=*), intent(in) :: catalogue
    character(len=*), intent(in) :: interp_type
 
    integer, intent(in), optional :: interp_level
    integer, intent(in), optional :: nmax
 
    real(rp) :: domain_lats, domain_late
    real(rp) :: domain_lons, domain_lone
 
    character(len=H_LONG) :: fname
 
    integer :: ierr
    integer :: i, j
    !---------------------------------------------------------------------------
 
    if ( present(nmax) ) then
       tile_nlim = nmax
    else
       tile_nlim = 1000
    end if
 
    domain_lats = minval( latxv(:,:) )
    domain_late = maxval( latxv(:,:) )
    domain_lons = minval( lonuy(:,:) )
    domain_lone = maxval( lonuy(:,:) )
 
    log_info("CNV2D_setup",*) 'Domain Information'
    log_info_cont(*) 'Domain (LAT)    :', domain_lats/d2r, domain_late/d2r
    log_info_cont(*) '       (LON)    :', domain_lons/d2r, domain_lone/d2r
 
    tile_dlat = dlat * d2r
    tile_dlon = dlon * d2r
 
    if ( .not. first ) then
       deallocate( tile_fname, tile_hit )
       deallocate( tile_js, tile_je, tile_is, tile_ie )
       deallocate( lat_1d, lon_1d )
    end if
 
    allocate( tile_fname(tile_nlim), tile_hit(tile_nlim) )
    allocate( tile_js(tile_nlim), tile_je(tile_nlim), tile_is(tile_nlim), tile_ie(tile_nlim) )
 
    ! catalogue file
    fname = trim(dir)//'/'//trim(catalogue)
 
    call file_tiledata_get_info( tile_nlim,                                          & ! [IN]
                                 tile_dlat, tile_dlon,                               & ! [IN]
                                 domain_lats, domain_late, domain_lons, domain_lone, & ! [IN]
                                 fname,                                              & ! [IN]
                                 global_ia,                                          & ! [OUT]
                                 tile_nmax,                                          & ! [OUT]
                                 tile_fname(:), tile_hit(:),                         & ! [OUT]
                                 tile_js(:), tile_je(:), tile_is(:), tile_ie(:),     & ! [OUT]
                                 nlat, nlon, dom_js, dom_je, dom_is, dom_ie,         & ! [OUT]
                                 zonal, pole                                         ) ! [OUT]
 
    allocate( lat_1d(nlat) )
    allocate( lon_1d(nlon) )
 
    call file_tiledata_get_latlon( nlat, nlon,           & ! [IN]
                                   dom_js, dom_is,       & ! [IN]
                                   tile_dlat, tile_dlon, & ! [IN]
                                   lat_1d(:), lon_1d(:)  ) ! [OUT]
 
    cnv2d_ftype = i_tile
    cnv2d_tile_dtype = dtype
    cnv2d_tile_dir   = dir
 
    call cnv2d_init( interp_type, &
                     interp_level = interp_level, &
                     ll_struct = .true.           )
 
    return
  end subroutine cnv2d_tile_init
 
 
  subroutine cnv2d_grads_init( &
       FILE_NAME,    &
       VAR_NAME,     &
       LAT_NAME,     &
       LON_NAME,     &
       interp_type,  &
       interp_level, &
       search_limit, &
       POSTFIX       )
    use scale_file_grads, only: &
       file_grads_open,      &
       file_grads_get_shape, &
       file_grads_varid,     &
       file_grads_isoned,    &
       file_grads_read
    use scale_const, only: &
       d2r => const_d2r
    implicit none
    character(len=*), intent(in) :: file_name
    character(len=*), intent(in) :: var_name
    character(len=*), intent(in) :: lat_name
    character(len=*), intent(in) :: lon_name
    character(len=*), intent(in) :: interp_type
    integer,          intent(in), optional :: interp_level
    real(rp),         intent(in), optional :: search_limit
    character(len=*), intent(in), optional :: postfix
 
    integer :: file_id, var_id
    integer :: shape(2)
 
    integer :: i, j
 
    call file_grads_open( file_name, & ! [IN]
                          file_id    ) ! [OUT]
 
    call file_grads_get_shape( file_id, var_name, &
                               shape(:)           )
    nlon = shape(1)
    nlat = shape(2)
 
    if ( .not. first ) deallocate( lat_org, lon_org, lat_1d, lon_1d )
    allocate( lat_org(nlon,nlat), lon_org(nlon,nlat) )
    allocate( lat_1d(nlat), lon_1d(nlon) )
 
    ! lat
    call file_grads_varid( file_id, lat_name, & ! (in)
                           var_id             ) ! (out)
    if ( file_grads_isoned( file_id, var_id ) ) then
       call file_grads_read( file_id, var_id, & ! (in)
                             lat_1d(:)        ) ! (out)
       !$omp parallel do collapse(2)
       do j = 1, nlat
       do i = 1, nlon
          lat_org(i,j) = lat_1d(j) * d2r
       end do
       end do
    else
       call file_grads_read( file_id, var_id,  & ! (in)
                             lat_org(:,:),     & ! (out)
                             postfix = postfix ) ! (in)
       !$omp parallel do collapse(2)
       do j = 1, nlat
       do i = 1, nlon
          lat_org(i,j) = lat_org(i,j) * d2r
       end do
       end do
    end if
 
    ! lon
    call file_grads_varid( file_id, lon_name, & ! (in)
                           var_id             ) ! (out)
    if ( file_grads_isoned( file_id, var_id ) ) then
       call file_grads_read( file_id, var_id, & ! (in)
                             lon_1d(:)         ) ! (out)
       !$omp parallel do collapse(2)
       do j = 1, nlat
       do i = 1, nlon
          lon_org(i,j) = lon_1d(i) * d2r
       end do
       end do
    else
       call file_grads_read( file_id, var_id,  & ! (in)
                             lon_org(:,:),     & ! (out)
                             postfix = postfix ) ! (in)
       !$omp parallel do collapse(2)
       do j = 1, nlat
       do i = 1, nlon
          lon_org(i,j) = lon_org(i,j) * d2r
       end do
       end do
    end if
 
    call file_grads_varid( file_id, var_name, & ! [IN]
                           var_id             ) ! [OUT]
 
    cnv2d_ftype = i_grads
    cnv2d_grads_fid = file_id
    cnv2d_grads_vid = var_id
 
    call cnv2d_init( interp_type, &
                     interp_level = interp_level, &
                     search_limit = search_limit, &
                     ll_struct = .false.          )
 
    return
  end subroutine cnv2d_grads_init
 
  subroutine cnv2d_exec( &
       var,       &
       step,      &
       min_value, &
       yrevers    )
    use scale_file_tiledata, only: &
       file_tiledata_get_data
    use scale_file_grads, only: &
       file_grads_read
    use scale_interp, only: &
       interp_interp2d
    implicit none
    real(rp), intent(out) :: var(ia,ja)
    integer,  intent(in), optional :: step
    real(rp), intent(in), optional :: min_value
    logical,  intent(in), optional :: yrevers
 
    select case ( cnv2d_ftype )
    case ( i_tile )
       call file_tiledata_get_data( nlat, nlon,                                     & ! [IN]
                                    cnv2d_tile_dir,                                 & ! [IN]
                                    global_ia,                                      & ! [IN]
                                    tile_nmax,                                      & ! [IN]
                                    tile_dlat, tile_dlon,                           & ! [IN]
                                    tile_fname(:), tile_hit(:),                     & ! [IN]
                                    tile_js(:), tile_je(:), tile_is(:), tile_ie(:), & ! [IN]
                                    dom_js, dom_je, dom_is, dom_ie,                 & ! [IN]
                                    cnv2d_tile_dtype,                               & ! [IN]
                                    data_org(:,:),                                  & ! [OUT]
                                    step = step,                                    & ! [IN]
                                    min_value = min_value, yrevers = yrevers        ) ! [IN]
    case ( i_grads )
       call file_grads_read( cnv2d_grads_fid, cnv2d_grads_vid, & ! [IN]
                             data_org(:,:),                    & ! [OUT]
                             step = step                       ) ! [IN]
    end select
 
 
    call interp_interp2d( itp_lev,                    & ! [IN]
                          nlon,nlat,                  & ! [IN]
                          ia, ja,                     & ! [IN]
                          idx_i(:,:,:), idx_j(:,:,:), & ! [IN]
                          hfact(:,:,:),               & ! [IN]
                          data_org(:,:),              & ! [IN]
                          var(:,:)                    ) ! [OUT]
 
    return
  end subroutine cnv2d_exec
 
 
  ! private
 
  subroutine cnv2d_init( &
       interp_type, &
       interp_level, &
       search_limit, &
       ll_struct     )
    use scale_prc, only: &
       prc_abort
    use scale_const, only: &
       pi => const_pi
    use scale_interp, only: &
       interp_factor2d_linear_latlon, &
       interp_factor2d_linear_xy,     &
       interp_factor2d_weight
    use scale_mapprojection, only: &
       mapprojection_lonlat2xy
    use scale_atmos_grid_cartesc, only: &
       cx => atmos_grid_cartesc_cx, &
       cy => atmos_grid_cartesc_cy
    use scale_atmos_grid_cartesc_real, only: &
       lat => atmos_grid_cartesc_real_lat, &
       lon => atmos_grid_cartesc_real_lon
    implicit none
 
    character(len=*), intent(in) :: interp_type
    integer,          intent(in), optional :: interp_level
    real(rp),         intent(in), optional :: search_limit
    logical,          intent(in), optional :: ll_struct
 
    real(rp), allocatable :: x_org(:,:), y_org(:,:)
    !---------------------------------------------------------------------------
 
    if ( .not. first ) deallocate( idx_i, idx_j, hfact, data_org )
 
    select case ( interp_type )
    case ( 'LINEAR' )
       itp_lev = 4
    case ( 'DIST-WEIGHT' )
       if ( present(interp_level) ) then
          itp_lev = interp_level
       else
          log_error('CNV2D_init',*) 'INTERP_LEVEL is required for the DIST-WEIGHT interpolation'
          call prc_abort
       end if
    case default
       log_error('CNV2D_init',*) 'INTERP_TYPE is invalid: ', trim(interp_type)
       call prc_abort
    end select
 
    ! interporation
    allocate( idx_i(ia,ja,itp_lev) )
    allocate( idx_j(ia,ja,itp_lev) )
    allocate( hfact(ia,ja,itp_lev) )
 
    select case ( interp_type )
    case ( 'LINEAR' )
       if ( ll_struct ) then
          call interp_factor2d_linear_latlon( nlon, nlat,                 & ! [IN]
                                              ia, ja,                     & ! [IN]
                                              lon_1d(:), lat_1d(:),       & ! [IN]
                                              lon(:,:), lat(:,:),         & ! [IN]
                                              idx_i(:,:,:), idx_j(:,:,:), & ! [OUT]
                                              hfact(:,:,:)                ) ! [OUT]
       else
          allocate( x_org(nlon,nlat), y_org(nlon,nlat) )
          call mapprojection_lonlat2xy( nlon, 1, nlon, nlat, 1, nlat, & ! [IN]
                                        lon_org(:,:), lat_org(:,:),   & ! [IN]
                                        x_org(:,:), y_org(:,:)        ) ! [OUT]
 
          zonal = ( maxval(lon_org) - minval(lat_org) ) > 2.0_rp * pi * 0.9_rp
          pole = ( maxval(lat_org) > pi * 0.5_rp * 0.9_rp ) .or. ( minval(lat_org) < - pi * 0.5_rp * 0.9_rp )
          call interp_factor2d_linear_xy( nlon, nlat,                 & ! [IN]
                                          ia, ja,                     & ! [IN]
                                          x_org(:,:), y_org(:,:),     & ! [IN]
                                          cx(:), cy(:),               & ! [IN]
                                          idx_i(:,:,:), idx_j(:,:,:), & ! [OUT]
                                          hfact(:,:,:),               & ! [OUT]
                                          zonal = zonal, pole = pole  ) ! [IN]
          deallocate( x_org, y_org )
       end if
    case ( 'DIST-WEIGHT' )
       call interp_factor2d_weight( itp_lev,       & ! [IN]
                                    nlon, nlat,                   & ! [IN]
                                    ia, ja,                       & ! [IN]
                                    lon_org(:,:), lat_org(:,:),   & ! [IN]
                                    lon(:,:), lat(:,:),           & ! [IN]
                                    idx_i(:,:,:), idx_j(:,:,:),   & ! [OUT]
                                    hfact(:,:,:),                 & ! [OUT]
                                    latlon_structure = ll_struct, & ! [IN]
                                    lon_1d = lon_1d(:),           & ! [IN]
                                    lat_1d = lat_1d(:),           & ! [IN]
                                    search_limit = search_limit   ) ! [IN]
    end select
 
    allocate( data_org(nlon,nlat) )
 
    first = .false.
 
    return
  end subroutine cnv2d_init
 
  ! private
 
end module mod_cnv2d