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

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

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '++++++ Module[GRID_REAL] / Categ[ATMOS-RM GRID] / Origin[SCALElib]'

    !--- 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_nml ) write(io_fid_nml,nml=param_domain_catalogue)

    allocate( real_lon(  ia,  ja) )
    allocate( real_lat(  ia,  ja) )
    allocate( real_lonx(0:ia,  ja) )
    allocate( real_lony(  ia,0:ja) )
    allocate( real_lonxy(0:ia,0:ja) )
    allocate( real_latx(0:ia,  ja) )
    allocate( real_laty(  ia,0:ja) )
    allocate( real_latxy(0:ia,0: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( domain_catalogue_fname, domain_catalogue_output )

    ! 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( &
       catalogue_fname, &
       catalogue_output )
    use scale_process, only: &
       prc_mpistop, &
       prc_nprocs,  &
       prc_ismaster
    use scale_const, only: &
       pi  => const_pi, &
       d2r => const_d2r
    use scale_grid, only: &
       grid_cx, &
       grid_cy, &
       grid_fx, &
       grid_fy
    use scale_comm, only: &
       comm_gather,  &
       comm_bcast
    use scale_mapproj, only: &
       mprj_basepoint_lon, &
       mprj_basepoint_lat, &
       mprj_xy2lonlat
    implicit none

    character(len=*), intent(in) :: catalogue_fname
    logical,          intent(in) :: catalogue_output

    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

    real(RP) :: mine (4,2)
    real(RP) :: whole(4,2*prc_nprocs)

    integer  :: i, j
    integer  :: fid, ierr
    !---------------------------------------------------------------------------

    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 (lat,lon)'
    if( io_l ) write(io_fid_log,*) '*** ->(',real_basepoint_lon/d2r,',',real_basepoint_lat/d2r,')'

    do j = 1, ja
    do i = 1, ia
       call mprj_xy2lonlat( grid_cx(i), grid_cy(j), real_lon(i,j), real_lat(i,j) )
    enddo
    enddo

    do j = 1, ja
    do i = 0, ia
       call mprj_xy2lonlat( grid_fx(i), grid_cy(j), real_lonx(i,j), real_latx(i,j) )
    enddo
    enddo

    do j = 0, ja
    do i = 1, ia
       call mprj_xy2lonlat( grid_cx(i), grid_fy(j), real_lony(i,j), real_laty(i,j) )
    enddo
    enddo

    do j = 0, ja
    do i = 0, ia
       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)
       if( real_dlon(i,j) < 0.0_rp ) real_dlon(i,j) = real_dlon(i,j) + 2.0_rp*pi

       if (      real_dlon(i,j) == 0.0_rp &
            .OR. real_dlat(i,j) == 0.0_rp ) then
          write(*,*) 'xxx Invalid grid distance in lat-lon! i,j=', i,j
          write(*,*) 'xxx Lon(i-1),Lon(i),dlon=', real_lonx(i-1,j)/d2r,real_lonx(i,j)/d2r,real_dlon(i,j)/d2r
          write(*,*) 'xxx Lat(j-1),Lat(j),dlat=', real_laty(i,j-1)/d2r,real_laty(i,j)/d2r,real_dlat(i,j)/d2r
          call prc_mpistop
       endif
    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,')'

    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 ( catalogue_output ) then

          fid = io_get_available_fid()
          open( fid,                            &
                file   = trim(catalogue_fname), &
                form   = 'formatted',           &
                status = 'replace',             &
                iostat = ierr                   )

          if ( ierr /= 0 ) then
             write(*,*) 'xxx [REAL_calc_latlon] 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