scale_grid_real.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
module scale_grid_real
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  use scale_prof
  use scale_grid_index
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: real_setup
  public :: real_update_z
  public :: real_calc_areavol

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  real(RP), public, allocatable :: real_lon(:,:)
  real(RP), public, allocatable :: real_lat(:,:)
  real(RP), public, allocatable :: real_cz (:,:,:)
  real(RP), public, allocatable :: real_fz (:,:,:)

  real(RP), public              :: real_basepoint_lon
  real(RP), public              :: real_basepoint_lat

  real(RP), public, allocatable :: real_lonx (:,:)
  real(RP), public, allocatable :: real_lony (:,:)
  real(RP), public, allocatable :: real_lonxy(:,:)
  real(RP), public, allocatable :: real_latx (:,:)
  real(RP), public, allocatable :: real_laty (:,:)
  real(RP), public, allocatable :: real_latxy(:,:)
  real(RP), public, allocatable :: real_dlon (:,:)
  real(RP), public, allocatable :: real_dlat (:,:)

  real(RP), public, allocatable :: real_z1  (:,:)
  real(RP), public              :: real_aspect_max
  real(RP), public              :: real_aspect_min

  real(RP), public, allocatable :: real_phi (:,:,:)

  real(RP), public, allocatable :: real_area(:,:)
  real(RP), public, allocatable :: real_vol (:,:,:)

  real(RP), public, allocatable :: real_domain_catalogue(:,:,:)

  real(RP), public :: real_totarea
  real(RP), public :: real_totvol

  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: real_calc_latlon
  private :: real_calc_z

  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine real_setup
    use scale_process, only: &
       prc_nprocs,  &
       prc_mpistop
    use scale_grid, only: &
       grid_domain_center_x, &
       grid_domain_center_y
    use scale_mapproj, only: &
       mprj_setup
    use scale_fileio, only: &
       fileio_set_coordinates
    implicit none

    !---------------------------------------------------------------------------

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '+++ Module[REAL]/Categ[GRID]'

    allocate( real_lon(ia,ja) )
    allocate( real_lat(ia,ja) )
    allocate( real_lonx(ia,ja) )
    allocate( real_lony(ia,ja) )
    allocate( real_lonxy(ia,ja) )
    allocate( real_latx(ia,ja) )
    allocate( real_laty(ia,ja) )
    allocate( real_latxy(ia,ja) )
    allocate( real_dlon(ia,ja) )
    allocate( real_dlat(ia,ja) )

    allocate( real_cz(  ka,ia,ja) )
    allocate( real_fz(0:ka,ia,ja) )
    allocate( real_z1(     ia,ja) )
    allocate( real_phi(  ka,ia,ja) )

    allocate( real_area(   ia,ja) )
    allocate( real_vol(ka,ia,ja) )

    allocate( real_domain_catalogue(prc_nprocs,4,2) )

    ! setup map projection
    call mprj_setup( grid_domain_center_x, grid_domain_center_y )

    ! calc longitude & latitude
    call real_calc_latlon

    ! calc real height
    call real_calc_z

    ! calc control area & volume
    ! call REAL_calc_areavol ! must be called after GTRANS_setup

    ! set latlon and z to fileio module
    call fileio_set_coordinates( real_lon, real_lonx, real_lony, real_lonxy, &
                                 real_lat, real_latx, real_laty, real_latxy, &
                                 real_cz,  real_fz    )

    return
  end subroutine real_setup

  !-----------------------------------------------------------------------------
  subroutine real_update_z
    use scale_process, only: &
       prc_mpistop
    use scale_fileio, only: &
       fileio_set_coordinates
    implicit none
    !---------------------------------------------------------------------------

    ! calc real height
    call real_calc_z

    ! set latlon and z to fileio module
    call fileio_set_coordinates( real_lon, real_lonx, real_lony, real_lonxy, &
                                 real_lat, real_latx, real_laty, real_latxy, &
                                 real_cz,  real_fz    )

    return
  end subroutine real_update_z

  !-----------------------------------------------------------------------------
  subroutine real_calc_latlon
    use scale_const, only: &
       d2r => const_d2r
    use scale_grid, only: &
       grid_domain_center_x, &
       grid_domain_center_y, &
       grid_cx, &
       grid_cy, &
       grid_fx, &
       grid_fy
    use scale_mapproj, only: &
       mprj_basepoint_lon, &
       mprj_basepoint_lat, &
       mprj_xy2lonlat
    use scale_process, only: &
       prc_mpistop, &
       prc_nprocs,  &
       prc_ismaster
    use scale_comm, only: &
       comm_gather,  &
       comm_bcast
    implicit none

    real(RP), allocatable :: mine(:,:)
    real(RP), allocatable :: whole(:,:)

    real(RP) :: CX, CY, FX, FY

    integer, parameter :: I_LON  = 1
    integer, parameter :: I_LAT  = 2
    integer, parameter :: I_NW   = 1
    integer, parameter :: I_NE   = 2
    integer, parameter :: I_SW   = 3
    integer, parameter :: I_SE   = 4

    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 :: i, j
    integer :: fid, ierr
    !---------------------------------------------------------------------------

    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_domain_catalogue,iostat=ierr)
    if( ierr < 0 ) then !--- missing
       if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       write(*,*) 'xxx Not appropriate names in namelist PARAM_DOMAIN_CATALOGUE. Check!'
       call prc_mpistop
    endif
    if( io_lnml ) write(io_fid_log,nml=param_domain_catalogue)

    real_basepoint_lon = mprj_basepoint_lon * d2r
    real_basepoint_lat = mprj_basepoint_lat * d2r

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '*** Base position in the global domain'
    if( io_l ) write(io_fid_log,*) '->(',real_basepoint_lon/d2r,',',real_basepoint_lat/d2r,')'

    do j = 1, ja
    do i = 1, ia
       ! apply offset
       cx = grid_cx(i) - grid_domain_center_x
       cy = grid_cy(j) - grid_domain_center_y
       fx = grid_fx(i) - grid_domain_center_x
       fy = grid_fy(j) - grid_domain_center_y

       call mprj_xy2lonlat( grid_cx(i), grid_cy(j), real_lon(i,j), real_lat(i,j) )
       call mprj_xy2lonlat( grid_fx(i), grid_cy(j), real_lonx(i,j), real_latx(i,j) )
       call mprj_xy2lonlat( grid_cx(i), grid_fy(j), real_lony(i,j), real_laty(i,j) )
       call mprj_xy2lonlat( grid_fx(i), grid_fy(j), real_lonxy(i,j), real_latxy(i,j) )
    enddo
    enddo

    real_dlon(:,:) = 0.0_rp
    real_dlat(:,:) = 0.0_rp
    do j = js, je
    do i = is, ie
       real_dlon(i,j) = real_lonx(i,j) - real_lonx(i-1,j)
       real_dlat(i,j) = real_laty(i,j) - real_laty(i,j-1)
    enddo
    enddo

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '*** Position on the earth (Local)'
    if( io_l ) write(io_fid_log,'(1x,A,F10.5,A,F9.5,A,A,F10.5,A,F9.5,A)') &
                                'NW(',real_lon(is,je)/d2r,',',real_lat(is,je)/d2r,')-', &
                                'NE(',real_lon(ie,je)/d2r,',',real_lat(ie,je)/d2r,')'
    if( io_l ) write(io_fid_log,'(1x,A)') &
                                '            |                       |'
    if( io_l ) write(io_fid_log,'(1x,A,F10.5,A,F9.5,A,A,F10.5,A,F9.5,A)') &
                                'SW(',real_lon(is,js)/d2r,',',real_lat(is,js)/d2r,')-', &
                                'SE(',real_lon(ie,js)/d2r,',',real_lat(ie,js)/d2r,')'


    allocate( mine(4, 2           ) )
    allocate( whole(4, 2*prc_nprocs) )

    mine(i_nw,i_lon) = real_lonxy(is-1,je  )/d2r
    mine(i_ne,i_lon) = real_lonxy(ie  ,je  )/d2r
    mine(i_sw,i_lon) = real_lonxy(is-1,js-1)/d2r
    mine(i_se,i_lon) = real_lonxy(ie  ,js-1)/d2r
    mine(i_nw,i_lat) = real_latxy(is-1,je  )/d2r
    mine(i_ne,i_lat) = real_latxy(ie  ,je  )/d2r
    mine(i_sw,i_lat) = real_latxy(is-1,js-1)/d2r
    mine(i_se,i_lat) = real_latxy(ie  ,js-1)/d2r

    call comm_gather( whole, mine, 4, 2 ) ! everytime do for online nesting

    if( prc_ismaster )then
       if( domain_catalogue_output ) then
          fid = io_get_available_fid()
          open( fid,                                    &
                file   = trim(domain_catalogue_fname), &
                form   = 'formatted',                  &
                status = 'replace',                    &
                iostat = ierr                          )

          if ( ierr /= 0 ) then
             if( io_l ) write(*,*) 'xxx cannot create latlon-catalogue file!'
             call prc_mpistop
          endif

          do i = 1, prc_nprocs ! for offline nesting
             write(fid,'(i8,8f32.24)',iostat=ierr) i, whole(i_nw,i_lon+2*(i-1)), whole(i_ne,i_lon+2*(i-1)), & ! LON: NW, NE
                                                      whole(i_sw,i_lon+2*(i-1)), whole(i_se,i_lon+2*(i-1)), & ! LON: SW, SE
                                                      whole(i_nw,i_lat+2*(i-1)), whole(i_ne,i_lat+2*(i-1)), & ! LAT: NW, NE
                                                      whole(i_sw,i_lat+2*(i-1)), whole(i_se,i_lat+2*(i-1))    ! LAT: SW, SE
             if ( ierr /= 0 ) exit
          enddo
          close(fid)
       endif

       do i = 1, prc_nprocs ! for online nesting
          real_domain_catalogue(i,i_nw,i_lon) = whole(i_nw,i_lon+2*(i-1)) ! LON: NW
          real_domain_catalogue(i,i_ne,i_lon) = whole(i_ne,i_lon+2*(i-1)) ! LON: NE
          real_domain_catalogue(i,i_sw,i_lon) = whole(i_sw,i_lon+2*(i-1)) ! LON: SW
          real_domain_catalogue(i,i_se,i_lon) = whole(i_se,i_lon+2*(i-1)) ! LON: SE
          real_domain_catalogue(i,i_nw,i_lat) = whole(i_nw,i_lat+2*(i-1)) ! LAT: NW
          real_domain_catalogue(i,i_ne,i_lat) = whole(i_ne,i_lat+2*(i-1)) ! LAT: NE
          real_domain_catalogue(i,i_sw,i_lat) = whole(i_sw,i_lat+2*(i-1)) ! LAT: SW
          real_domain_catalogue(i,i_se,i_lat) = whole(i_se,i_lat+2*(i-1)) ! LAT: SE
       enddo
    endif

    call comm_bcast( real_domain_catalogue, prc_nprocs, 4, 2 )

    return
  end subroutine real_calc_latlon

  !-----------------------------------------------------------------------------
  subroutine real_calc_z
    use scale_const, only: &
       grav => const_grav
    use scale_grid, only: &
       grid_cz,  &
       grid_fz,  &
       grid_cdx, &
       grid_cdy
    use scale_topography, only: &
       zsfc => topo_zsfc
    implicit none

    real(RP) :: Htop
    real(RP) :: DFZ

    integer  :: k, i, j
    !---------------------------------------------------------------------------

    htop = grid_fz(ke) - grid_fz(ks-1)

    do j = 1, ja
    do i = 1, ia
    do k = 1, ka
       real_cz(k,i,j) = ( htop - zsfc(i,j) ) / htop * grid_cz(k) + zsfc(i,j)
    enddo
    enddo
    enddo

    do j = 1, ja
    do i = 1, ia
    do k = 0, ka
       real_fz(k,i,j) = ( htop - zsfc(i,j) ) / htop * grid_fz(k) + zsfc(i,j)
    enddo
    enddo
    enddo

    real_z1(:,:) = real_cz(ks,:,:) - zsfc(:,:)

    real_phi(:,:,:) = grav * real_cz(:,:,:)

    real_aspect_max = -1.e+30_rp
    real_aspect_min =  1.e+30_rp
    do j = js, je
    do i = is, ie
    do k = ks, ke
       dfz = real_fz(k,i,j) - real_fz(k-1,i,j)
       real_aspect_max = max( real_aspect_max, grid_cdx(i) / dfz, grid_cdy(j) / dfz )
       real_aspect_min = min( real_aspect_min, grid_cdx(i) / dfz, grid_cdy(j) / dfz )
    enddo
    enddo
    enddo

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '*** Minimum & maximum aspect ratio'
    if( io_l ) write(io_fid_log,*) '->(',real_aspect_min,',',real_aspect_max,')'

    return
  end subroutine real_calc_z

  !-----------------------------------------------------------------------------
  subroutine real_calc_areavol( &
       MAPF )
    use scale_const, only: &
       radius => const_radius
    use scale_grid, only: &
       dz, &
       dx, &
       dy
    implicit none
    real(RP), intent(in) :: MAPF(ia,ja,2)

    integer :: k, i, j
    !---------------------------------------------------------------------------

    real_totarea   = 0.0_rp
    real_area(:,:) = 0.0_rp
    do j = js, je
    do i = is, ie
       real_area(i,j) = dx * dy / ( mapf(i,j,1) * mapf(i,j,2) )
       real_totarea = real_totarea + real_area(i,j)
    enddo
    enddo

    real_totvol     = 0.0_rp
    real_vol(:,:,:) = 0.0_rp
    do j = js, je
    do i = is, ie
    do k = ks, ke
       real_vol(k,i,j) = ( real_fz(k,i,j) - real_fz(k-1,i,j) ) * real_area(i,j)
       real_totvol = real_totvol + real_vol(k,i,j)
    enddo
    enddo
    enddo

    return
  end subroutine real_calc_areavol

end module scale_grid_real