scale_atmos_grid_cartesC_real.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_atmos_grid_cartesc_real
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  use scale_atmos_grid_cartesc_index
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: atmos_grid_cartesc_real_setup
  public :: atmos_grid_cartesc_real_calc_z
  public :: atmos_grid_cartesc_real_calc_areavol
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  real(rp), public              :: atmos_grid_cartesc_real_basepoint_lon
  real(rp), public              :: atmos_grid_cartesc_real_basepoint_lat
 
  real(rp), public, allocatable :: atmos_grid_cartesc_real_cz  (:,:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_czuy(:,:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_czxv(:,:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_czuv(:,:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_fz  (:,:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_fzuy(:,:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_fzxv(:,:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_fzuv(:,:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_f2h (:,:,:,:)
 
  real(rp), public, allocatable :: atmos_grid_cartesc_real_lon  (:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_lonuy(:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_lonxv(:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_lonuv(:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_lat  (:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_latuy(:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_latxv(:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_latuv(:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_dlon (:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_dlat (:,:)
 
  real(rp), public, allocatable :: atmos_grid_cartesc_real_z1  (:,:)
  real(rp), public              :: atmos_grid_cartesc_real_aspect_max
  real(rp), public              :: atmos_grid_cartesc_real_aspect_min
 
  real(rp), public, allocatable :: atmos_grid_cartesc_real_phi (:,:,:)
 
  real(rp), public, allocatable :: atmos_grid_cartesc_real_area     (:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_areazuy_x(:,:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_areazxv_y(:,:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_areawuy_x(:,:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_areawxv_y(:,:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_areauy   (:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_areazxy_x(:,:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_areazuv_y(:,:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_areaxv   (:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_areazuv_x(:,:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_areazxy_y(:,:,:)
 
  real(rp), public              :: atmos_grid_cartesc_real_totarea
  real(rp), public              :: atmos_grid_cartesc_real_totareauy
  real(rp), public              :: atmos_grid_cartesc_real_totareaxv
  real(rp), public, allocatable :: atmos_grid_cartesc_real_totareazuy_x(:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_totareazxv_y(:)
 
  real(rp), public, allocatable :: atmos_grid_cartesc_real_vol   (:,:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_volwxy(:,:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_volzuy(:,:,:)
  real(rp), public, allocatable :: atmos_grid_cartesc_real_volzxv(:,:,:)
  real(rp), public              :: atmos_grid_cartesc_real_totvol
  real(rp), public              :: atmos_grid_cartesc_real_totvolwxy
  real(rp), public              :: atmos_grid_cartesc_real_totvolzuy
  real(rp), public              :: atmos_grid_cartesc_real_totvolzxv
 
  real(rp), public, allocatable :: atmos_grid_cartesc_real_domain_catalogue(:,:,:)
 
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: atmos_grid_cartesc_real_calc_latlon
 
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine atmos_grid_cartesc_real_setup
    use scale_prc, only: &
       prc_nprocs,  &
       prc_abort
    use scale_atmos_grid_cartesc, only: &
       atmos_grid_cartesc_domain_center_x, &
       atmos_grid_cartesc_domain_center_y, &
       atmos_grid_cartesc_cz,              &
       atmos_grid_cartesc_fz
    use scale_topography, only: &
       topography_exist
    use scale_mapprojection, only: &
       mapprojection_setup
    use scale_interp_vert, only: &
       interp_vert_setcoef
    implicit none
 
    character(len=H_LONG) :: domain_catalogue_fname  = 'latlon_domain_catalogue.txt'
    logical               :: domain_catalogue_output = .false.
 
    namelist / param_domain_catalogue / &
       domain_catalogue_fname,  &
       domain_catalogue_output
 
    integer :: ierr
    !---------------------------------------------------------------------------
 
    log_newline
    log_info("ATMOS_GRID_CARTESC_REAL_setup",*) 'Setup'
 
    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_domain_catalogue,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       log_info("ATMOS_GRID_CARTESC_REAL_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("ATMOS_GRID_CARTESC_REAL_setup",*) 'Not appropriate names in namelist PARAM_DOMAIN_CATALOGUE. Check!'
       call prc_abort
    endif
    log_nml(param_domain_catalogue)
 
    allocate( atmos_grid_cartesc_real_lon(  ia,  ja) )
    allocate( atmos_grid_cartesc_real_lat(  ia,  ja) )
    allocate( atmos_grid_cartesc_real_lonuy(0:ia,  ja) )
    allocate( atmos_grid_cartesc_real_lonxv(  ia,0:ja) )
    allocate( atmos_grid_cartesc_real_lonuv(0:ia,0:ja) )
    allocate( atmos_grid_cartesc_real_latuy(0:ia,  ja) )
    allocate( atmos_grid_cartesc_real_latxv(  ia,0:ja) )
    allocate( atmos_grid_cartesc_real_latuv(0:ia,0:ja) )
    allocate( atmos_grid_cartesc_real_dlon(  ia,  ja) )
    allocate( atmos_grid_cartesc_real_dlat(  ia,  ja) )
 
    allocate( atmos_grid_cartesc_real_cz(  ka,ia,ja) )
    allocate( atmos_grid_cartesc_real_czuy(  ka,ia,ja) )
    allocate( atmos_grid_cartesc_real_czxv(  ka,ia,ja) )
    allocate( atmos_grid_cartesc_real_czuv(  ka,ia,ja) )
    allocate( atmos_grid_cartesc_real_fz(0:ka,ia,ja) )
    allocate( atmos_grid_cartesc_real_fzuy(0:ka,ia,ja) )
    allocate( atmos_grid_cartesc_real_fzxv(0:ka,ia,ja) )
    allocate( atmos_grid_cartesc_real_fzuv(0:ka,ia,ja) )
    allocate( atmos_grid_cartesc_real_f2h(ka,2,ia,ja) )
    allocate( atmos_grid_cartesc_real_z1(     ia,ja) )
    allocate( atmos_grid_cartesc_real_phi(  ka,ia,ja) )
 
    allocate( atmos_grid_cartesc_real_area(     ia,ja) )
    allocate( atmos_grid_cartesc_real_areazuy_x(ka  ,ia,ja) )
    allocate( atmos_grid_cartesc_real_areazxv_y(ka  ,ia,ja) )
    allocate( atmos_grid_cartesc_real_areawuy_x(ka+1,ia,ja) )
    allocate( atmos_grid_cartesc_real_areawxv_y(ka+1,ia,ja) )
    allocate( atmos_grid_cartesc_real_areauy(     ia,ja) )
    allocate( atmos_grid_cartesc_real_areazxy_x(ka,  ia,ja) )
    allocate( atmos_grid_cartesc_real_areazuv_y(ka,  ia,ja) )
    allocate( atmos_grid_cartesc_real_areaxv(     ia,ja) )
    allocate( atmos_grid_cartesc_real_areazuv_x(ka,  ia,ja) )
    allocate( atmos_grid_cartesc_real_areazxy_y(ka,  ia,ja) )
 
    allocate( atmos_grid_cartesc_real_totareazuy_x(ia) )
    allocate( atmos_grid_cartesc_real_totareazxv_y(ja) )
 
    allocate( atmos_grid_cartesc_real_vol(  ka,ia,ja) )
    allocate( atmos_grid_cartesc_real_volwxy(0:ka,ia,ja) )
    allocate( atmos_grid_cartesc_real_volzuy(  ka,ia,ja) )
    allocate( atmos_grid_cartesc_real_volzxv(  ka,ia,ja) )
 
    allocate( atmos_grid_cartesc_real_domain_catalogue(prc_nprocs,2,2) )
 
    ! setup map projection
    call mapprojection_setup( atmos_grid_cartesc_domain_center_x, atmos_grid_cartesc_domain_center_y )
 
    ! calc longitude & latitude
    call atmos_grid_cartesc_real_calc_latlon( domain_catalogue_fname, domain_catalogue_output )
 
    ! calc real height
    call atmos_grid_cartesc_real_calc_z
 
    call interp_vert_setcoef( ka, ks, ke, ia, 1,  ia, ja, 1,  ja, & ! [IN]
                              topography_exist,                  & ! [IN]
                              atmos_grid_cartesc_cz(:),          & ! [IN]
                              atmos_grid_cartesc_fz(:),          & ! [IN]
                              atmos_grid_cartesc_real_cz(:,:,:), & ! [IN]
                              atmos_grid_cartesc_real_fz(:,:,:)  ) ! [IN]
 
    return
  end subroutine atmos_grid_cartesc_real_setup
 
  !-----------------------------------------------------------------------------
  subroutine atmos_grid_cartesc_real_calc_latlon( &
       catalogue_fname, &
       catalogue_output )
    use scale_prc, only: &
       prc_abort, &
       prc_nprocs,  &
       prc_ismaster
    use scale_const, only: &
       pi  => const_pi, &
       d2r => const_d2r
    use scale_atmos_grid_cartesc, only: &
       atmos_grid_cartesc_cx, &
       atmos_grid_cartesc_cy, &
       atmos_grid_cartesc_fx, &
       atmos_grid_cartesc_fy
    use scale_comm_cartesc, only: &
       comm_gather,  &
       comm_bcast
    use scale_mapprojection, only: &
       mapprojection_basepoint_lon, &
       mapprojection_basepoint_lat, &
       mapprojection_xy2lonlat
    implicit none
 
    character(len=*), intent(in) :: catalogue_fname
    logical,          intent(in) :: catalogue_output
 
    integer, parameter :: i_min  = 1
    integer, parameter :: i_max  = 2
    integer, parameter :: i_lon  = 1
    integer, parameter :: i_lat  = 2
 
    real(rp) :: mine (2,2)
    real(rp) :: whole(2,2,prc_nprocs)
 
    integer  :: i, j
    integer  :: fid, ierr
    !---------------------------------------------------------------------------
 
    atmos_grid_cartesc_real_basepoint_lon = mapprojection_basepoint_lon * d2r
    atmos_grid_cartesc_real_basepoint_lat = mapprojection_basepoint_lat * d2r
 
    log_newline
    log_info("ATMOS_GRID_CARTESC_REAL_calc_latlon",*) 'Base position in the global domain (lat,lon)'
    log_info_cont(*) '-> (',atmos_grid_cartesc_real_basepoint_lon/d2r,',',atmos_grid_cartesc_real_basepoint_lat/d2r,')'
 
    !$omp parallel do collapse(2)
    do j = 1, ja
    do i = 1, ia
       call mapprojection_xy2lonlat( atmos_grid_cartesc_cx(i), atmos_grid_cartesc_cy(j), atmos_grid_cartesc_real_lon(i,j), atmos_grid_cartesc_real_lat(i,j) )
    enddo
    enddo
 
    !$omp parallel do collapse(2)
    do j = 1, ja
    do i = 0, ia
       call mapprojection_xy2lonlat( atmos_grid_cartesc_fx(i), atmos_grid_cartesc_cy(j), atmos_grid_cartesc_real_lonuy(i,j), atmos_grid_cartesc_real_latuy(i,j) )
    enddo
    enddo
 
    !$omp parallel do collapse(2)
    do j = 0, ja
    do i = 1, ia
       call mapprojection_xy2lonlat( atmos_grid_cartesc_cx(i), atmos_grid_cartesc_fy(j), atmos_grid_cartesc_real_lonxv(i,j), atmos_grid_cartesc_real_latxv(i,j) )
    enddo
    enddo
 
    !$omp parallel do collapse(2)
    do j = 0, ja
    do i = 0, ia
       call mapprojection_xy2lonlat( atmos_grid_cartesc_fx(i), atmos_grid_cartesc_fy(j), atmos_grid_cartesc_real_lonuv(i,j), atmos_grid_cartesc_real_latuv(i,j) )
    enddo
    enddo
 
    !$omp workshare
!OCL ZFILL
    atmos_grid_cartesc_real_dlon(:,:) = 0.0_rp
!OCL ZFILL
    atmos_grid_cartesc_real_dlat(:,:) = 0.0_rp
    !$omp end workshare
 
    !$omp parallel do
    do j = js, je
    do i = is, ie
       atmos_grid_cartesc_real_dlon(i,j) = abs( atmos_grid_cartesc_real_lonuy(i,j) - atmos_grid_cartesc_real_lonuy(i-1,j) )
       atmos_grid_cartesc_real_dlat(i,j) = abs( atmos_grid_cartesc_real_latxv(i,j) - atmos_grid_cartesc_real_latxv(i,j-1) )
       if( atmos_grid_cartesc_real_dlon(i,j) > 2.0_rp*pi - atmos_grid_cartesc_real_dlon(i,j)  ) atmos_grid_cartesc_real_dlon(i,j) = 2.0_rp*pi - atmos_grid_cartesc_real_dlon(i,j)
 
       if (      atmos_grid_cartesc_real_dlon(i,j) == 0.0_rp &
            .OR. atmos_grid_cartesc_real_dlat(i,j) == 0.0_rp ) then
          log_error("ATMOS_GRID_CARTESC_REAL_calc_latlon",*) 'Invalid grid distance in lat-lon! i,j=', i,j
          log_error_cont(*) 'Lon(i-1),Lon(i),dlon=', atmos_grid_cartesc_real_lonuy(i-1,j)/d2r,atmos_grid_cartesc_real_lonuy(i,j)/d2r,atmos_grid_cartesc_real_dlon(i,j)/d2r
          log_error_cont(*) 'Lat(j-1),Lat(j),dlat=', atmos_grid_cartesc_real_latxv(i,j-1)/d2r,atmos_grid_cartesc_real_latxv(i,j)/d2r,atmos_grid_cartesc_real_dlat(i,j)/d2r
          call prc_abort
       endif
    enddo
    enddo
 
    log_newline
    log_info("ATMOS_GRID_CARTESC_REAL_calc_latlon",*) 'Position on the earth (Local)'
    log_info_cont('(1x,A,F10.5,A,F9.5,A,A,F10.5,A,F9.5,A)') &
                               'NW(',atmos_grid_cartesc_real_lon(is,je)/d2r,',',atmos_grid_cartesc_real_lat(is,je)/d2r,')', &
                            ' - NE(',atmos_grid_cartesc_real_lon(ie,je)/d2r,',',atmos_grid_cartesc_real_lat(ie,je)/d2r,')'
 
    log_info_cont('(1x,A)') '             |                          |'
    log_info_cont('(1x,A,F10.5,A,F9.5,A,A,F10.5,A,F9.5,A)') &
                               'SW(',atmos_grid_cartesc_real_lon(is,js)/d2r,',',atmos_grid_cartesc_real_lat(is,js)/d2r,')', &
                            ' - SE(',atmos_grid_cartesc_real_lon(ie,js)/d2r,',',atmos_grid_cartesc_real_lat(ie,js)/d2r,')'
 
    mine(i_min,i_lon) = minval(atmos_grid_cartesc_real_lonuv(:,:)) / d2r
    mine(i_max,i_lon) = maxval(atmos_grid_cartesc_real_lonuv(:,:)) / d2r
    mine(i_min,i_lat) = minval(atmos_grid_cartesc_real_latuv(:,:)) / d2r
    mine(i_max,i_lat) = maxval(atmos_grid_cartesc_real_latuv(:,:)) / d2r
 
    call comm_gather( whole(:,:,:), mine(:,:), 2, 2 ) ! everytime do for online nesting
 
    if ( prc_ismaster ) then
       if ( catalogue_output ) then
 
          fid = io_get_available_fid()
          open( fid,                            &
                file   = trim(catalogue_fname), &
                form   = 'formatted',           &
                status = 'replace',             &
                iostat = ierr                   )
 
          if ( ierr /= 0 ) then
             log_error("ATMOS_GRID_CARTESC_REAL_calc_latlon",*) 'cannot create latlon-catalogue file!'
             call prc_abort
          endif
 
          do i = 1, prc_nprocs ! for offline nesting
             write(fid,'(I8,8F32.24)',iostat=ierr) i, whole(i_min,i_lon,i), whole(i_max,i_lon,i), & ! LON: MIN, MAX
                                                      whole(i_min,i_lat,i), whole(i_max,i_lat,i)    ! LAT: MIN, MAX
             if ( ierr /= 0 ) exit
          enddo
 
          close(fid)
 
       endif
 
       do i = 1, prc_nprocs ! for online nesting
          atmos_grid_cartesc_real_domain_catalogue(i,i_min,i_lon) = whole(i_min,i_lon,i)
          atmos_grid_cartesc_real_domain_catalogue(i,i_max,i_lon) = whole(i_max,i_lon,i)
          atmos_grid_cartesc_real_domain_catalogue(i,i_min,i_lat) = whole(i_min,i_lat,i)
          atmos_grid_cartesc_real_domain_catalogue(i,i_max,i_lat) = whole(i_max,i_lat,i)
       enddo
    endif
 
    call comm_bcast( atmos_grid_cartesc_real_domain_catalogue(:,:,:), prc_nprocs, 2, 2 )
 
    return
  end subroutine atmos_grid_cartesc_real_calc_latlon
 
  !-----------------------------------------------------------------------------
  subroutine atmos_grid_cartesc_real_calc_z
    use scale_const, only: &
       grav => const_grav
    use scale_atmos_grid_cartesc, only: &
       atmos_grid_cartesc_cz,  &
       atmos_grid_cartesc_fz,  &
       atmos_grid_cartesc_cdx, &
       atmos_grid_cartesc_cdy
    use scale_file_cartesc, only: &
       file_cartesc_set_coordinates_atmos
    use scale_topography, only: &
       zsfc => topography_zsfc
    use scale_landuse, only: &
       landuse_frac_land
    implicit none
 
    real(dp) :: htop
    real(rp) :: zs
    real(rp) :: dfz
 
    real(rp) :: dz1, dz2
 
    integer  :: k, i, j
    !---------------------------------------------------------------------------
 
    htop = atmos_grid_cartesc_fz(ke) - atmos_grid_cartesc_fz(ks-1)
 
    !$omp parallel do private(zs) collapse(2)
    do j = 1, ja
    do i = 1, ia
       zs = zsfc(i,j)
       do k = 1, ka
          atmos_grid_cartesc_real_cz(k,i,j) = ( htop - zs ) / htop * atmos_grid_cartesc_cz(k) + zs
       enddo
    enddo
    enddo
 
    !$omp parallel do private(zs) collapse(2)
    do j = 1, ja
    do i = 1, ia-1
       zs = ( zsfc(i,j) + zsfc(i+1,j) ) * 0.5_rp
       do k = 1, ka
          atmos_grid_cartesc_real_czuy(k,i,j) = ( htop - zs ) / htop * atmos_grid_cartesc_cz(k) + zs
       enddo
    enddo
    enddo
    !$omp parallel do private(zs)
    do j = 1, ja
       zs = zsfc(ia,j)
       do k = 1, ka
          atmos_grid_cartesc_real_czuy(k,ia,j) = ( htop - zs ) / htop * atmos_grid_cartesc_cz(k) + zs
       enddo
    enddo
 
    !$omp parallel do private(zs) collapse(2)
    do j = 1, ja-1
    do i = 1, ia
       zs = ( zsfc(i,j) + zsfc(i,j+1) ) * 0.5_rp
       do k = 1, ka
          atmos_grid_cartesc_real_czxv(k,i,j) = ( htop - zs ) / htop * atmos_grid_cartesc_cz(k) + zs
       enddo
    enddo
    enddo
    do i = 1, ia
       zs = zsfc(i,ja)
       do k = 1, ka
          atmos_grid_cartesc_real_czxv(k,i,ja) = ( htop - zs ) / htop * atmos_grid_cartesc_cz(k) + zs
       enddo
    enddo
 
    !$omp parallel do private(zs) collapse(2)
    do j = 1, ja-1
    do i = 1, ia-1
       zs = ( zsfc(i,j) + zsfc(i+1,j) + zsfc(i,j+1) + zsfc(i+1,j+1) ) * 0.25_rp
       do k = 1, ka
          atmos_grid_cartesc_real_czuv(k,i,j) = ( htop - zs ) / htop * atmos_grid_cartesc_cz(k) + zs
       enddo
    enddo
    enddo
    !$omp parallel do private(zs)
    do j = 1, ja-1
       zs = ( zsfc(ia,j) + zsfc(ia,j+1) ) * 0.5_rp
       do k = 1, ka
          atmos_grid_cartesc_real_czuv(k,ia,j) = ( htop - zs ) / htop * atmos_grid_cartesc_cz(k) + zs
       enddo
    enddo
    do i = 1, ia-1
       zs = ( zsfc(i,ja) + zsfc(i+1,ja) ) * 0.5_rp
       do k = 1, ka
          atmos_grid_cartesc_real_czuv(k,i,ja) = ( htop - zs ) / htop * atmos_grid_cartesc_cz(k) + zs
       enddo
    enddo
    zs = zsfc(ia,ja)
    do k = 1, ka
       atmos_grid_cartesc_real_czuv(k,ia,ja) = ( htop - zs ) / htop * atmos_grid_cartesc_cz(k) + zs
    enddo
 
 
    !$omp parallel do private(zs) collapse(2)
    do j = 1, ja
    do i = 1, ia
       zs = zsfc(i,j)
       do k = 0, ka
          atmos_grid_cartesc_real_fz(k,i,j) = ( htop - zs ) / htop * atmos_grid_cartesc_fz(k) + zs
       end do
    end do
    end do
 
    !$omp parallel do private(zs) collapse(2)
    do j = 1, ja
    do i = 1, ia-1
       zs = ( zsfc(i,j) + zsfc(i+1,j) ) * 0.5_rp
       do k = 0, ka
          atmos_grid_cartesc_real_fzuy(k,i,j) = ( htop - zs ) / htop * atmos_grid_cartesc_fz(k) + zs
       end do
    end do
    end do
    !$omp parallel do private(zs)
    do j = 1, ja
       zs = zsfc(ia,j)
       do k = 0, ka
          atmos_grid_cartesc_real_fzuy(k,ia,j) = ( htop - zs ) / htop * atmos_grid_cartesc_fz(k) + zs
       end do
    end do
 
    !$omp parallel do private(zs) collapse(2)
    do j = 1, ja-1
    do i = 1, ia
       zs = ( zsfc(i,j) + zsfc(i,j+1) ) * 0.5_rp
       do k = 0, ka
          atmos_grid_cartesc_real_fzxv(k,i,j) = ( htop - zs ) / htop * atmos_grid_cartesc_fz(k) + zs
       enddo
    enddo
    enddo
    !$omp parallel do private(zs)
    do i = 1, ia
       zs = zsfc(i,ja)
       do k = 0, ka
          atmos_grid_cartesc_real_fzxv(k,i,ja) = ( htop - zs ) / htop * atmos_grid_cartesc_fz(k) + zs
       enddo
    enddo
 
    !$omp parallel do private(zs) collapse(2)
    do j = 1, ja-1
    do i = 1, ia-1
       zs = ( zsfc(i,j) + zsfc(i+1,j) + zsfc(i,j+1) + zsfc(i+1,j+1) ) * 0.25_rp
       do k = 0, ka
          atmos_grid_cartesc_real_fzuv(k,i,j) = ( htop - zs ) / htop * atmos_grid_cartesc_fz(k) + zs
       enddo
    enddo
    enddo
    !$omp parallel do private(zs)
    do j = 1, ja-1
       zs = ( zsfc(ia,j) + zsfc(ia,j+1) ) * 0.5_rp
       do k = 0, ka
          atmos_grid_cartesc_real_fzuv(k,ia,j) = ( htop - zs ) / htop * atmos_grid_cartesc_fz(k) + zs
       enddo
    enddo
    !$omp parallel do private(zs)
    do i = 1, ia-1
       zs = ( zsfc(i,ja) + zsfc(i+1,ja) ) * 0.5_rp
       do k = 0, ka
          atmos_grid_cartesc_real_fzuv(k,i,ja) = ( htop - zs ) / htop * atmos_grid_cartesc_fz(k) + zs
       enddo
    enddo
    zs = zsfc(ia,ja)
    do k = 0, ka
       atmos_grid_cartesc_real_fzuv(k,ia,ja) = ( htop - zs ) / htop * atmos_grid_cartesc_fz(k) + zs
    enddo
 
    !$omp parallel do private(dz1,dz2) collapse(2)
    do j = 1, ja
    do i = 1, ia
       do k = ks, ke-1
          dz1 = atmos_grid_cartesc_real_fz(k+1,i,j) - atmos_grid_cartesc_real_fz(k  ,i,j)
          dz2 = atmos_grid_cartesc_real_fz(k  ,i,j) - atmos_grid_cartesc_real_fz(k-1,i,j)
          atmos_grid_cartesc_real_f2h(k,1,i,j) = dz2 / ( dz1 + dz2 )
          atmos_grid_cartesc_real_f2h(k,2,i,j) = dz1 / ( dz1 + dz2 )
       end do
       atmos_grid_cartesc_real_f2h(1:ks-1,1,i,j) = 0.5_rp
       atmos_grid_cartesc_real_f2h(1:ks-1,2,i,j) = 0.5_rp
       atmos_grid_cartesc_real_f2h(ke:ka ,1,i,j) = 0.5_rp
       atmos_grid_cartesc_real_f2h(ke:ka ,2,i,j) = 0.5_rp
    end do
    end do
 
 
    !$omp workshare
!OCL ZFILL
    atmos_grid_cartesc_real_z1(:,:) = atmos_grid_cartesc_real_cz(ks,:,:) - zsfc(:,:)
 
!OCL ZFILL
    atmos_grid_cartesc_real_phi(:,:,:) = grav * atmos_grid_cartesc_real_cz(:,:,:)
    !$omp end workshare
 
    atmos_grid_cartesc_real_aspect_max = -1.e+30_rp
    atmos_grid_cartesc_real_aspect_min =  1.e+30_rp
 
    !$omp parallel do private(dfz) collapse(2) &
    !$omp reduction(max:ATMOS_GRID_CARTESC_REAL_ASPECT_MAX) &
    !$omp reduction(min:ATMOS_GRID_CARTESC_REAL_ASPECT_MIN)
    do j = js, je
    do i = is, ie
    do k = ks, ke
       dfz = atmos_grid_cartesc_real_fz(k,i,j) - atmos_grid_cartesc_real_fz(k-1,i,j)
       atmos_grid_cartesc_real_aspect_max = max( atmos_grid_cartesc_real_aspect_max, atmos_grid_cartesc_cdx(i) / dfz, atmos_grid_cartesc_cdy(j) / dfz )
       atmos_grid_cartesc_real_aspect_min = min( atmos_grid_cartesc_real_aspect_min, atmos_grid_cartesc_cdx(i) / dfz, atmos_grid_cartesc_cdy(j) / dfz )
    enddo
    enddo
    enddo
 
    log_newline
    log_info("ATMOS_GRID_CARTESC_REAL_calc_Z",*) 'Minimum & maximum lowermost CZ'
    log_info_cont(*) '-> (',minval( atmos_grid_cartesc_real_cz(ks,:,:) ),',',maxval( atmos_grid_cartesc_real_cz(ks,:,:) ),')'
    log_info("ATMOS_GRID_CARTESC_REAL_calc_Z",*) 'Minimum & maximum aspect ratio'
    log_info_cont(*) '-> (',atmos_grid_cartesc_real_aspect_min,',',atmos_grid_cartesc_real_aspect_max,')'
 
    ! set latlon and z to fileio module
    call file_cartesc_set_coordinates_atmos( atmos_grid_cartesc_real_cz,  atmos_grid_cartesc_real_fz,                                                                      & ! [IN]
                                             atmos_grid_cartesc_real_lon, atmos_grid_cartesc_real_lonuy, atmos_grid_cartesc_real_lonxv, atmos_grid_cartesc_real_lonuv,     & ! [IN]
                                             atmos_grid_cartesc_real_lat, atmos_grid_cartesc_real_latuy, atmos_grid_cartesc_real_latxv, atmos_grid_cartesc_real_latuv,     & ! [IN]
                                             zsfc, landuse_frac_land                                                                                                      ) ! [IN]
 
    return
  end subroutine atmos_grid_cartesc_real_calc_z
 
  !-----------------------------------------------------------------------------
  subroutine atmos_grid_cartesc_real_calc_areavol( &
       MAPF )
    use scale_prc_cartesc, only: &
       prc_twod
    use scale_atmos_grid_cartesc, only: &
       atmos_grid_cartesc_cdx, &
       atmos_grid_cartesc_fdx, &
       atmos_grid_cartesc_cdy, &
       atmos_grid_cartesc_fdy
    use scale_comm_cartesc, only: &
       comm_vars8, &
       comm_wait
    use scale_file_cartesc, only: &
       file_cartesc_set_areavol_atmos
    use scale_topography, only: &
       topography_zsfc
    use scale_landuse, only: &
       landuse_frac_land
    implicit none
 
    real(rp), intent(in) :: mapf(ia,ja,2,4)
 
    real(rp) :: areauv(ia,ja)
 
    integer :: k, i, j
    !---------------------------------------------------------------------------
 
    atmos_grid_cartesc_real_area(:,:)   = 0.0_rp
    atmos_grid_cartesc_real_areazuy_x(:,:,:) = 0.0_rp
    atmos_grid_cartesc_real_areazxv_y(:,:,:) = 0.0_rp
    atmos_grid_cartesc_real_areawuy_x(:,:,:) = 0.0_rp
    atmos_grid_cartesc_real_areawxv_y(:,:,:) = 0.0_rp
    atmos_grid_cartesc_real_areauy(:,:)   = 0.0_rp
    atmos_grid_cartesc_real_areazxy_x(:,:,:) = 0.0_rp
    atmos_grid_cartesc_real_areazuv_y(:,:,:) = 0.0_rp
    atmos_grid_cartesc_real_areaxv(:,:)   = 0.0_rp
    atmos_grid_cartesc_real_areazuv_y(:,:,:) = 0.0_rp
    atmos_grid_cartesc_real_areazxy_y(:,:,:) = 0.0_rp
 
    atmos_grid_cartesc_real_totarea         = 0.0_rp
    atmos_grid_cartesc_real_totareauy       = 0.0_rp
    atmos_grid_cartesc_real_totareaxv       = 0.0_rp
    atmos_grid_cartesc_real_totareazuy_x(:) = 0.0_rp
    atmos_grid_cartesc_real_totareazxv_y(:) = 0.0_rp
 
    atmos_grid_cartesc_real_vol(:,:,:) = 0.0_rp
    atmos_grid_cartesc_real_volwxy(:,:,:) = 0.0_rp
    atmos_grid_cartesc_real_volzuy(:,:,:) = 0.0_rp
    atmos_grid_cartesc_real_volzxv(:,:,:) = 0.0_rp
    atmos_grid_cartesc_real_totvol    = 0.0_rp
    atmos_grid_cartesc_real_totvolwxy = 0.0_rp
    atmos_grid_cartesc_real_totvolzuy = 0.0_rp
    atmos_grid_cartesc_real_totvolzxv = 0.0_rp
 
    !$omp parallel do
    do j = js, je
    do i = is, ie
       atmos_grid_cartesc_real_area(i,j) = atmos_grid_cartesc_cdx(i) * atmos_grid_cartesc_cdy(j) / ( mapf(i,j,1,i_xy) * mapf(i,j,2,i_xy) )
       atmos_grid_cartesc_real_areaxv(i,j) = atmos_grid_cartesc_cdx(i) * atmos_grid_cartesc_fdy(j) / ( mapf(i,j,1,i_xv) * mapf(i,j,2,i_xv) )
    end do
    end do
    if ( prc_twod ) then
       !$omp parallel do
       do j = js, je
       do i = is, ie
          atmos_grid_cartesc_real_areauy(i,j) = atmos_grid_cartesc_real_area(i,j)
                                  areauv(i,j) = atmos_grid_cartesc_real_areaxv(i,j)
 
       end do
       end do
    else
       !$omp parallel do
       do j = js, je
       do i = is, ie
          atmos_grid_cartesc_real_areauy(i,j) = atmos_grid_cartesc_fdx(i) * atmos_grid_cartesc_cdy(j) / ( mapf(i,j,1,i_uy) * mapf(i,j,2,i_uy) )
                                  areauv(i,j) = atmos_grid_cartesc_fdx(i) * atmos_grid_cartesc_fdy(j) / ( mapf(i,j,1,i_uv) * mapf(i,j,2,i_uv) )
 
       end do
       end do
    end if
 
    call comm_vars8( atmos_grid_cartesc_real_area(:,:), 1 )
    call comm_vars8( atmos_grid_cartesc_real_areaxv(:,:), 2 )
    call comm_vars8( atmos_grid_cartesc_real_areauy(:,:), 3 )
    call comm_vars8(                         areauv(:,:), 4 )
 
    !$omp parallel do &
    !$omp reduction(+:ATMOS_GRID_CARTESC_REAL_TOTAREA,ATMOS_GRID_CARTESC_REAL_TOTAREAXV,ATMOS_GRID_CARTESC_REAL_TOTAREAUY)
    do j = js, je
    do i = is, ie
       atmos_grid_cartesc_real_totarea   = atmos_grid_cartesc_real_totarea   + atmos_grid_cartesc_real_area(i,j)
       atmos_grid_cartesc_real_totareaxv = atmos_grid_cartesc_real_totareaxv + atmos_grid_cartesc_real_areaxv(i,j)
       atmos_grid_cartesc_real_totareauy = atmos_grid_cartesc_real_totareauy + atmos_grid_cartesc_real_areauy(i,j)
    end do
    end do
 
    !$omp parallel do collapse(2)
    do j = 1, ja
    do i = 1, ia
       do k = ks, ke
          atmos_grid_cartesc_real_areazxv_y(k,i,j) = atmos_grid_cartesc_cdx(i) / mapf(i,j,1,i_xv) * ( atmos_grid_cartesc_real_fzxv(k,i,j) - atmos_grid_cartesc_real_fzxv(k-1,i,j) )
       end do
       do k = ks-1, ke
          atmos_grid_cartesc_real_areawxv_y(k,i,j) = atmos_grid_cartesc_cdx(i) / mapf(i,j,1,i_xv) * ( atmos_grid_cartesc_real_czxv(k+1,i,j) - atmos_grid_cartesc_real_czxv(k,i,j) )
       end do
       do k = ks, ke
          atmos_grid_cartesc_real_areazxy_x(k,i,j) = atmos_grid_cartesc_cdy(j) / mapf(i,j,2,i_xy) * ( atmos_grid_cartesc_real_fz(k,i,j) - atmos_grid_cartesc_real_fz(k-1,i,j) )
          atmos_grid_cartesc_real_areazxy_y(k,i,j) = atmos_grid_cartesc_cdx(i) / mapf(i,j,1,i_xy) * ( atmos_grid_cartesc_real_fz(k,i,j) - atmos_grid_cartesc_real_fz(k-1,i,j) )
       end do
    end do
    end do
 
    if ( .not. prc_twod ) then
       !$omp parallel do collapse(2)
       do j = 1, ja
       do i = 1, ia
          do k = ks, ke
             atmos_grid_cartesc_real_areazuy_x(k,i,j) = atmos_grid_cartesc_cdy(j) / mapf(i,j,2,i_uy) * ( atmos_grid_cartesc_real_fzuy(k,i,j) - atmos_grid_cartesc_real_fzuy(k-1,i,j) )
          end do
          do k = ks-1, ke
             atmos_grid_cartesc_real_areawuy_x(k,i,j) = atmos_grid_cartesc_cdy(j) / mapf(i,j,2,i_uy) * ( atmos_grid_cartesc_real_czuy(k+1,i,j) - atmos_grid_cartesc_real_czuy(k,i,j) )
          end do
          do k = ks, ke
             atmos_grid_cartesc_real_areazuv_y(k,i,j) = atmos_grid_cartesc_cdx(i) / mapf(i,j,1,i_uv) * ( atmos_grid_cartesc_real_fzuv(k,i,j) - atmos_grid_cartesc_real_fzuv(k-1,i,j) )
             atmos_grid_cartesc_real_areazuv_x(k,i,j) = atmos_grid_cartesc_cdy(j) / mapf(i,j,2,i_uv) * ( atmos_grid_cartesc_real_fzuv(k,i,j) - atmos_grid_cartesc_real_fzuv(k-1,i,j) )
          end do
       end do
       end do
    end if
 
    call comm_wait( atmos_grid_cartesc_real_area(:,:), 1 )
    call comm_wait( atmos_grid_cartesc_real_areaxv(:,:), 2 )
    call comm_wait( atmos_grid_cartesc_real_areauy(:,:), 3 )
    call comm_wait(                         areauv(:,:), 4 )
 
 
    !$omp parallel do collapse(2)
    do j = 1,  ja
    do i = is, ie
    do k = ks, ke
       atmos_grid_cartesc_real_totareazxv_y(j) = atmos_grid_cartesc_real_totareazxv_y(j) + atmos_grid_cartesc_real_areazxv_y(k,i,j)
    end do
    end do
    end do
    !$omp parallel do collapse(2)
    do j = js, je
    do i = 1,  ia
    do k = ks, ke
       atmos_grid_cartesc_real_totareazuy_x(i) = atmos_grid_cartesc_real_totareazuy_x(i) + atmos_grid_cartesc_real_areazuy_x(k,i,j)
    end do
    end do
    end do
 
 
    !$omp parallel do collapse(2)
    do j = 1, ja
    do i = 1, ia
       do k = ks, ke
          atmos_grid_cartesc_real_vol(k,i,j) = ( atmos_grid_cartesc_real_fz(k,i,j) - atmos_grid_cartesc_real_fz(k-1,i,j) ) * atmos_grid_cartesc_real_area(i,j)
          atmos_grid_cartesc_real_volzxv(k,i,j) = ( atmos_grid_cartesc_real_fzxv(k,i,j) - atmos_grid_cartesc_real_fzxv(k-1,i,j) ) * atmos_grid_cartesc_real_areaxv(i,j)
       end do
       do k = ks-1, ke
          atmos_grid_cartesc_real_volwxy(k,i,j) = ( atmos_grid_cartesc_real_cz(k+1,i,j) - atmos_grid_cartesc_real_cz(k,i,j) ) * atmos_grid_cartesc_real_area(i,j)
       end do
    end do
    end do
    if ( prc_twod ) then
       !$omp parallel do collapse(2)
       do j = 1, ja
       do i = 1, ia
       do k = ks, ke
          atmos_grid_cartesc_real_volzuy(k,i,j) = atmos_grid_cartesc_real_vol(k,i,j)
       end do
       end do
       end do
    else
       !$omp parallel do collapse(2)
       do j = 1, ja
       do i = 1, ia
       do k = ks, ke
          atmos_grid_cartesc_real_volzuy(k,i,j) = ( atmos_grid_cartesc_real_fzuy(k,i,j) - atmos_grid_cartesc_real_fzuy(k-1,i,j) ) * atmos_grid_cartesc_real_areauy(i,j)
       end do
       end do
       end do
    end if
 
    !$omp parallel do &
    !$omp reduction(+:ATMOS_GRID_CARTESC_REAL_TOTVOL,ATMOS_GRID_CARTESC_REAL_TOTVOLZXV,ATMOS_GRID_CARTESC_REAL_TOTVOLWXY,ATMOS_GRID_CARTESC_REAL_TOTVOLZUY)
    do j = js, je
    do i = is, ie
    do k = ks, ke
       atmos_grid_cartesc_real_totvol    = atmos_grid_cartesc_real_totvol    + atmos_grid_cartesc_real_vol(k,i,j)
       atmos_grid_cartesc_real_totvolzxv = atmos_grid_cartesc_real_totvolzxv + atmos_grid_cartesc_real_volzxv(k,i,j)
       atmos_grid_cartesc_real_totvolwxy = atmos_grid_cartesc_real_totvolwxy + atmos_grid_cartesc_real_volwxy(k,i,j)
       atmos_grid_cartesc_real_totvolzuy = atmos_grid_cartesc_real_totvolzuy + atmos_grid_cartesc_real_volzuy(k,i,j)
    enddo
    enddo
    enddo
 
    ! set latlon and z to fileio module
    call file_cartesc_set_areavol_atmos( atmos_grid_cartesc_real_area,   atmos_grid_cartesc_real_areazuy_x, atmos_grid_cartesc_real_areazxv_y,                         & ! [IN]
                                                                         atmos_grid_cartesc_real_areawuy_x, atmos_grid_cartesc_real_areawxv_y,                         & ! [IN]
                                         atmos_grid_cartesc_real_areauy, atmos_grid_cartesc_real_areazxy_x, atmos_grid_cartesc_real_areazuv_y,                         & ! [IN]
                                         atmos_grid_cartesc_real_areaxv, atmos_grid_cartesc_real_areazuv_x, atmos_grid_cartesc_real_areazxy_y,                         & ! [IN]
                                         atmos_grid_cartesc_real_vol, atmos_grid_cartesc_real_volwxy, atmos_grid_cartesc_real_volzuy, atmos_grid_cartesc_real_volzxv   ) ! [IN]
 
    return
  end subroutine atmos_grid_cartesc_real_calc_areavol
 
end module scale_atmos_grid_cartesc_real