MOM_transcribe_grid.F90

!> Module with routines for copying information from a shared dynamic horizontal
!! grid to an ocean-specific horizontal grid and the reverse.
module mom_transcribe_grid

! This file is part of MOM6. See LICENSE.md for the license.

use mom_array_transform, only: rotate_array, rotate_array_pair
use mom_domains, only : pass_var, pass_vector
use mom_domains, only : to_all, scalar_pair, cgrid_ne, agrid, bgrid_ne, corner
use mom_dyn_horgrid, only : dyn_horgrid_type, set_derived_dyn_horgrid
use mom_error_handler, only : mom_error, mom_mesg, fatal, warning
use mom_grid, only : ocean_grid_type, set_derived_metrics
use mom_unit_scaling, only : unit_scale_type


implicit none ; private

public copy_dyngrid_to_mom_grid, copy_mom_grid_to_dyngrid, rotate_dyngrid

contains

!> Copies information from a dynamic (shared) horizontal grid type into an
!! ocean_grid_type.
subroutine copy_dyngrid_to_mom_grid(dG, oG, US)
  type(dyn_horgrid_type), intent(in)    :: dg  !< Common horizontal grid type
  type(ocean_grid_type),  intent(inout) :: og  !< Ocean grid type
  type(unit_scale_type),  intent(in)    :: us  !< A dimensional unit scaling type

  integer :: isd, ied, jsd, jed      ! Common data domains.
  integer :: isdb, iedb, jsdb, jedb  ! Common data domains.
  integer :: ido, jdo, ido2, jdo2    ! Indexing offsets between the grids.
  integer :: igst, jgst              ! Global starting indices.
  integer :: i, j

  ! MOM_grid_init and create_dyn_horgrid are called outside of this routine.
  ! This routine copies over the fields that were set by MOM_initialized_fixed.

  ! Determine the indexing offsets between the grids.
  ido = dg%idg_offset - og%idg_offset
  jdo = dg%jdg_offset - og%jdg_offset

  isd = max(og%isd, dg%isd+ido) ; jsd = max(og%jsd, dg%jsd+jdo)
  ied = min(og%ied, dg%ied+ido) ; jed = min(og%jed, dg%jed+jdo)
  isdb = max(og%IsdB, dg%IsdB+ido) ; jsdb = max(og%JsdB, dg%JsdB+jdo)
  iedb = min(og%IedB, dg%IedB+ido) ; jedb = min(og%JedB, dg%JedB+jdo)

  ! Check that the grids conform.
  if ((isd > og%isc) .or. (ied < og%ied) .or. (jsd > og%jsc) .or. (jed > og%jed)) &
    call mom_error(fatal, "copy_dyngrid_to_MOM_grid called with incompatible grids.")

  do i=isd,ied ; do j=jsd,jed
    og%geoLonT(i,j) = dg%geoLonT(i+ido,j+jdo)
    og%geoLatT(i,j) = dg%geoLatT(i+ido,j+jdo)
    og%dxT(i,j) = dg%dxT(i+ido,j+jdo)
    og%dyT(i,j) = dg%dyT(i+ido,j+jdo)
    og%areaT(i,j) = dg%areaT(i+ido,j+jdo)
    og%bathyT(i,j) = dg%bathyT(i+ido,j+jdo)

    og%dF_dx(i,j) = dg%dF_dx(i+ido,j+jdo)
    og%dF_dy(i,j) = dg%dF_dy(i+ido,j+jdo)
    og%sin_rot(i,j) = dg%sin_rot(i+ido,j+jdo)
    og%cos_rot(i,j) = dg%cos_rot(i+ido,j+jdo)
    og%mask2dT(i,j) = dg%mask2dT(i+ido,j+jdo)
  enddo ; enddo

  do i=isdb,iedb ; do j=jsd,jed
    og%geoLonCu(i,j) = dg%geoLonCu(i+ido,j+jdo)
    og%geoLatCu(i,j) = dg%geoLatCu(i+ido,j+jdo)
    og%dxCu(i,j) = dg%dxCu(i+ido,j+jdo)
    og%dyCu(i,j) = dg%dyCu(i+ido,j+jdo)
    og%dy_Cu(i,j) = dg%dy_Cu(i+ido,j+jdo)

    og%mask2dCu(i,j) = dg%mask2dCu(i+ido,j+jdo)
    og%areaCu(i,j) = dg%areaCu(i+ido,j+jdo)
    og%IareaCu(i,j) = dg%IareaCu(i+ido,j+jdo)
  enddo ; enddo

  do i=isd,ied ; do j=jsdb,jedb
    og%geoLonCv(i,j) = dg%geoLonCv(i+ido,j+jdo)
    og%geoLatCv(i,j) = dg%geoLatCv(i+ido,j+jdo)
    og%dxCv(i,j) = dg%dxCv(i+ido,j+jdo)
    og%dyCv(i,j) = dg%dyCv(i+ido,j+jdo)
    og%dx_Cv(i,j) = dg%dx_Cv(i+ido,j+jdo)

    og%mask2dCv(i,j) = dg%mask2dCv(i+ido,j+jdo)
    og%areaCv(i,j) = dg%areaCv(i+ido,j+jdo)
    og%IareaCv(i,j) = dg%IareaCv(i+ido,j+jdo)
  enddo ; enddo

  do i=isdb,iedb ; do j=jsdb,jedb
    og%geoLonBu(i,j) = dg%geoLonBu(i+ido,j+jdo)
    og%geoLatBu(i,j) = dg%geoLatBu(i+ido,j+jdo)
    og%dxBu(i,j) = dg%dxBu(i+ido,j+jdo)
    og%dyBu(i,j) = dg%dyBu(i+ido,j+jdo)
    og%areaBu(i,j) = dg%areaBu(i+ido,j+jdo)
    og%CoriolisBu(i,j) = dg%CoriolisBu(i+ido,j+jdo)
    og%mask2dBu(i,j) = dg%mask2dBu(i+ido,j+jdo)
  enddo ; enddo

  og%bathymetry_at_vel = dg%bathymetry_at_vel
  if (og%bathymetry_at_vel) then
    do i=isdb,iedb ; do j=jsd,jed
      og%Dblock_u(i,j) = dg%Dblock_u(i+ido,j+jdo)
      og%Dopen_u(i,j) = dg%Dopen_u(i+ido,j+jdo)
    enddo ; enddo
    do i=isd,ied ; do j=jsdb,jedb
      og%Dblock_v(i,j) = dg%Dblock_v(i+ido,j+jdo)
      og%Dopen_v(i,j) = dg%Dopen_v(i+ido,j+jdo)
    enddo ; enddo
  endif

  og%gridLonT(og%isg:og%ieg) = dg%gridLonT(dg%isg:dg%ieg)
  og%gridLatT(og%jsg:og%jeg) = dg%gridLatT(dg%jsg:dg%jeg)
  ! The more complicated logic here avoids segmentation faults if one grid uses
  ! global symmetric memory while the other does not.  Because a northeast grid
  ! convention is being used, the upper bounds for each array correspond.
  !   Note that the dynamic grid always uses symmetric memory.
  ido2 = dg%IegB-og%IegB ; igst = max(og%IsgB, (dg%isg-1)-ido2)
  jdo2 = dg%JegB-og%JegB ; jgst = max(og%JsgB, (dg%jsg-1)-jdo2)
  do i=igst,og%IegB ; og%gridLonB(i) = dg%gridLonB(i+ido2) ; enddo
  do j=jgst,og%JegB ; og%gridLatB(j) = dg%gridLatB(j+jdo2) ; enddo

  ! Copy various scalar variables and strings.
  og%x_axis_units = dg%x_axis_units ; og%y_axis_units = dg%y_axis_units
  og%areaT_global = dg%areaT_global ; og%IareaT_global = dg%IareaT_global
  og%south_lat = dg%south_lat ; og%west_lon  = dg%west_lon
  og%len_lat = dg%len_lat ; og%len_lon = dg%len_lon
  og%Rad_Earth = dg%Rad_Earth ; og%max_depth = dg%max_depth

! Update the halos in case the dynamic grid has smaller halos than the ocean grid.
  call pass_var(og%areaT, og%Domain)
  call pass_var(og%bathyT, og%Domain)
  call pass_var(og%geoLonT, og%Domain)
  call pass_var(og%geoLatT, og%Domain)
  call pass_vector(og%dxT, og%dyT, og%Domain, to_all+scalar_pair, agrid)
  call pass_vector(og%dF_dx, og%dF_dy, og%Domain, to_all, agrid)
  call pass_vector(og%cos_rot, og%sin_rot, og%Domain, to_all, agrid)
  call pass_var(og%mask2dT, og%Domain)

  call pass_vector(og%areaCu, og%areaCv, og%Domain, to_all+scalar_pair, cgrid_ne)
  call pass_vector(og%dyCu, og%dxCv, og%Domain, to_all+scalar_pair, cgrid_ne)
  call pass_vector(og%dxCu, og%dyCv, og%Domain, to_all+scalar_pair, cgrid_ne)
  call pass_vector(og%dy_Cu, og%dx_Cv, og%Domain, to_all+scalar_pair, cgrid_ne)
  call pass_vector(og%mask2dCu, og%mask2dCv, og%Domain, to_all+scalar_pair, cgrid_ne)
  call pass_vector(og%IareaCu, og%IareaCv, og%Domain, to_all+scalar_pair, cgrid_ne)
  call pass_vector(og%IareaCu, og%IareaCv, og%Domain, to_all+scalar_pair, cgrid_ne)
  call pass_vector(og%geoLatCu, og%geoLatCv, og%Domain, to_all+scalar_pair, cgrid_ne)

  call pass_var(og%areaBu, og%Domain, position=corner)
  call pass_var(og%geoLonBu, og%Domain, position=corner, inner_halo=og%isc-isd)
  call pass_var(og%geoLatBu, og%Domain, position=corner)
  call pass_vector(og%dxBu, og%dyBu, og%Domain, to_all+scalar_pair, bgrid_ne)
  call pass_var(og%CoriolisBu, og%Domain, position=corner)
  call pass_var(og%mask2dBu, og%Domain, position=corner)

  if (og%bathymetry_at_vel) then
    call pass_vector(og%Dblock_u, og%Dblock_v, og%Domain, to_all+scalar_pair, cgrid_ne)
    call pass_vector(og%Dopen_u, og%Dopen_v, og%Domain, to_all+scalar_pair, cgrid_ne)
  endif

  call set_derived_metrics(og, us)

end subroutine copy_dyngrid_to_mom_grid


!> Copies information from an ocean_grid_type into a dynamic (shared)
!! horizontal grid type.
subroutine copy_mom_grid_to_dyngrid(oG, dG, US)
  type(ocean_grid_type),  intent(in)    :: og  !< Ocean grid type
  type(dyn_horgrid_type), intent(inout) :: dg  !< Common horizontal grid type
  type(unit_scale_type), optional, intent(in) :: us !< A dimensional unit scaling type

  integer :: isd, ied, jsd, jed      ! Common data domains.
  integer :: isdb, iedb, jsdb, jedb  ! Common data domains.
  integer :: ido, jdo, ido2, jdo2    ! Indexing offsets between the grids.
  integer :: igst, jgst              ! Global starting indices.
  integer :: i, j

  ! MOM_grid_init and create_dyn_horgrid are called outside of this routine.
  ! This routine copies over the fields that were set by MOM_initialized_fixed.

  ! Determine the indexing offsets between the grids.
  ido = og%idG_offset - dg%idG_offset
  jdo = og%jdG_offset - dg%jdG_offset

  isd = max(dg%isd, og%isd+ido) ; jsd = max(dg%jsd, og%jsd+jdo)
  ied = min(dg%ied, og%ied+ido) ; jed = min(dg%jed, og%jed+jdo)
  isdb = max(dg%IsdB, og%IsdB+ido) ; jsdb = max(dg%JsdB, og%JsdB+jdo)
  iedb = min(dg%IedB, og%IedB+ido) ; jedb = min(dg%JedB, og%JedB+jdo)

  ! Check that the grids conform.
  if ((isd > dg%isc) .or. (ied < dg%ied) .or. (jsd > dg%jsc) .or. (jed > dg%jed)) &
    call mom_error(fatal, "copy_dyngrid_to_MOM_grid called with incompatible grids.")

  do i=isd,ied ; do j=jsd,jed
    dg%geoLonT(i,j) = og%geoLonT(i+ido,j+jdo)
    dg%geoLatT(i,j) = og%geoLatT(i+ido,j+jdo)
    dg%dxT(i,j) = og%dxT(i+ido,j+jdo)
    dg%dyT(i,j) = og%dyT(i+ido,j+jdo)
    dg%areaT(i,j) = og%areaT(i+ido,j+jdo)
    dg%bathyT(i,j) = og%bathyT(i+ido,j+jdo)

    dg%dF_dx(i,j) = og%dF_dx(i+ido,j+jdo)
    dg%dF_dy(i,j) = og%dF_dy(i+ido,j+jdo)
    dg%sin_rot(i,j) = og%sin_rot(i+ido,j+jdo)
    dg%cos_rot(i,j) = og%cos_rot(i+ido,j+jdo)
    dg%mask2dT(i,j) = og%mask2dT(i+ido,j+jdo)
  enddo ; enddo

  do i=isdb,iedb ; do j=jsd,jed
    dg%geoLonCu(i,j) = og%geoLonCu(i+ido,j+jdo)
    dg%geoLatCu(i,j) = og%geoLatCu(i+ido,j+jdo)
    dg%dxCu(i,j) = og%dxCu(i+ido,j+jdo)
    dg%dyCu(i,j) = og%dyCu(i+ido,j+jdo)
    dg%dy_Cu(i,j) = og%dy_Cu(i+ido,j+jdo)

    dg%mask2dCu(i,j) = og%mask2dCu(i+ido,j+jdo)
    dg%areaCu(i,j) = og%areaCu(i+ido,j+jdo)
    dg%IareaCu(i,j) = og%IareaCu(i+ido,j+jdo)
  enddo ; enddo

  do i=isd,ied ; do j=jsdb,jedb
    dg%geoLonCv(i,j) = og%geoLonCv(i+ido,j+jdo)
    dg%geoLatCv(i,j) = og%geoLatCv(i+ido,j+jdo)
    dg%dxCv(i,j) = og%dxCv(i+ido,j+jdo)
    dg%dyCv(i,j) = og%dyCv(i+ido,j+jdo)
    dg%dx_Cv(i,j) = og%dx_Cv(i+ido,j+jdo)

    dg%mask2dCv(i,j) = og%mask2dCv(i+ido,j+jdo)
    dg%areaCv(i,j) = og%areaCv(i+ido,j+jdo)
    dg%IareaCv(i,j) = og%IareaCv(i+ido,j+jdo)
  enddo ; enddo

  do i=isdb,iedb ; do j=jsdb,jedb
    dg%geoLonBu(i,j) = og%geoLonBu(i+ido,j+jdo)
    dg%geoLatBu(i,j) = og%geoLatBu(i+ido,j+jdo)
    dg%dxBu(i,j) = og%dxBu(i+ido,j+jdo)
    dg%dyBu(i,j) = og%dyBu(i+ido,j+jdo)
    dg%areaBu(i,j) = og%areaBu(i+ido,j+jdo)
    dg%CoriolisBu(i,j) = og%CoriolisBu(i+ido,j+jdo)
    dg%mask2dBu(i,j) = og%mask2dBu(i+ido,j+jdo)
  enddo ; enddo

  dg%bathymetry_at_vel = og%bathymetry_at_vel
  if (dg%bathymetry_at_vel) then
    do i=isdb,iedb ; do j=jsd,jed
      dg%Dblock_u(i,j) = og%Dblock_u(i+ido,j+jdo)
      dg%Dopen_u(i,j) = og%Dopen_u(i+ido,j+jdo)
    enddo ; enddo
    do i=isd,ied ; do j=jsdb,jedb
      dg%Dblock_v(i,j) = og%Dblock_v(i+ido,j+jdo)
      dg%Dopen_v(i,j) = og%Dopen_v(i+ido,j+jdo)
    enddo ; enddo
  endif

  dg%gridLonT(dg%isg:dg%ieg) = og%gridLonT(og%isg:og%ieg)
  dg%gridLatT(dg%jsg:dg%jeg) = og%gridLatT(og%jsg:og%jeg)

  ! The more complicated logic here avoids segmentation faults if one grid uses
  ! global symmetric memory while the other does not.  Because a northeast grid
  ! convention is being used, the upper bounds for each array correspond.
  !   Note that the dynamic grid always uses symmetric memory.
  ido2 = og%IegB-dg%IegB ; igst = max(dg%isg-1, og%IsgB-ido2)
  jdo2 = og%JegB-dg%JegB ; jgst = max(dg%jsg-1, og%JsgB-jdo2)
  do i=igst,dg%IegB ; dg%gridLonB(i) = og%gridLonB(i+ido2) ; enddo
  do j=jgst,dg%JegB ; dg%gridLatB(j) = og%gridLatB(j+jdo2) ; enddo

  ! Copy various scalar variables and strings.
  dg%x_axis_units = og%x_axis_units ; dg%y_axis_units = og%y_axis_units
  dg%areaT_global = og%areaT_global ; dg%IareaT_global = og%IareaT_global
  dg%south_lat = og%south_lat ; dg%west_lon  = og%west_lon
  dg%len_lat = og%len_lat ; dg%len_lon = og%len_lon
  dg%Rad_Earth = og%Rad_Earth ; dg%max_depth = og%max_depth

! Update the halos in case the dynamic grid has smaller halos than the ocean grid.
  call pass_var(dg%areaT, dg%Domain)
  call pass_var(dg%bathyT, dg%Domain)
  call pass_var(dg%geoLonT, dg%Domain)
  call pass_var(dg%geoLatT, dg%Domain)
  call pass_vector(dg%dxT, dg%dyT, dg%Domain, to_all+scalar_pair, agrid)
  call pass_vector(dg%dF_dx, dg%dF_dy, dg%Domain, to_all, agrid)
  call pass_vector(dg%cos_rot, dg%sin_rot, dg%Domain, to_all, agrid)
  call pass_var(dg%mask2dT, dg%Domain)

  call pass_vector(dg%areaCu, dg%areaCv, dg%Domain, to_all+scalar_pair, cgrid_ne)
  call pass_vector(dg%dyCu, dg%dxCv, dg%Domain, to_all+scalar_pair, cgrid_ne)
  call pass_vector(dg%dxCu, dg%dyCv, dg%Domain, to_all+scalar_pair, cgrid_ne)
  call pass_vector(dg%dy_Cu, dg%dx_Cv, dg%Domain, to_all+scalar_pair, cgrid_ne)
  call pass_vector(dg%mask2dCu, dg%mask2dCv, dg%Domain, to_all+scalar_pair, cgrid_ne)
  call pass_vector(dg%IareaCu, dg%IareaCv, dg%Domain, to_all+scalar_pair, cgrid_ne)
  call pass_vector(dg%IareaCu, dg%IareaCv, dg%Domain, to_all+scalar_pair, cgrid_ne)
  call pass_vector(dg%geoLatCu, dg%geoLatCv, dg%Domain, to_all+scalar_pair, cgrid_ne)

  call pass_var(dg%areaBu, dg%Domain, position=corner)
  call pass_var(dg%geoLonBu, dg%Domain, position=corner, inner_halo=dg%isc-isd)
  call pass_var(dg%geoLatBu, dg%Domain, position=corner)
  call pass_vector(dg%dxBu, dg%dyBu, dg%Domain, to_all+scalar_pair, bgrid_ne)
  call pass_var(dg%CoriolisBu, dg%Domain, position=corner)
  call pass_var(dg%mask2dBu, dg%Domain, position=corner)

  if (dg%bathymetry_at_vel) then
    call pass_vector(dg%Dblock_u, dg%Dblock_v, dg%Domain, to_all+scalar_pair, cgrid_ne)
    call pass_vector(dg%Dopen_u, dg%Dopen_v, dg%Domain, to_all+scalar_pair, cgrid_ne)
  endif

  call  set_derived_dyn_horgrid(dg, us)

end subroutine copy_mom_grid_to_dyngrid

subroutine rotate_dyngrid(G_in, G, US, turns)
  type(dyn_horgrid_type), intent(in)    :: g_in   !< Common horizontal grid type
  type(dyn_horgrid_type), intent(inout) :: g      !< Ocean grid type
  type(unit_scale_type),  intent(in)    :: us     !< A dimensional unit scaling type
  integer, intent(in) :: turns                    !< Number of quarter turns

  integer :: jsc, jec, jscb, jecb
  integer :: qturn

  ! Center point
  call rotate_array(g_in%geoLonT, turns, g%geoLonT)
  call rotate_array(g_in%geoLatT, turns, g%geoLatT)
  call rotate_array_pair(g_in%dxT, g_in%dyT, turns, g%dxT, g%dyT)
  call rotate_array(g_in%areaT, turns, g%areaT)
  call rotate_array(g_in%bathyT, turns, g%bathyT)

  call rotate_array_pair(g_in%df_dx, g_in%df_dy, turns, g%df_dx, g%df_dy)
  call rotate_array(g_in%sin_rot, turns, g%sin_rot)
  call rotate_array(g_in%cos_rot, turns, g%cos_rot)
  call rotate_array(g_in%mask2dT, turns, g%mask2dT)

  ! Face point
  call rotate_array_pair(g_in%geoLonCu, g_in%geoLonCv, turns, &
      g%geoLonCu, g%geoLonCv)
  call rotate_array_pair(g_in%geoLatCu, g_in%geoLatCv, turns, &
      g%geoLatCu, g%geoLatCv)
  call rotate_array_pair(g_in%dxCu, g_in%dyCv, turns, g%dxCu, g%dyCv)
  call rotate_array_pair(g_in%dxCv, g_in%dyCu, turns, g%dxCv, g%dyCu)
  call rotate_array_pair(g_in%dx_Cv, g_in%dy_Cu, turns, g%dx_Cv, g%dy_Cu)

  call rotate_array_pair(g_in%mask2dCu, g_in%mask2dCv, turns, &
      g%mask2dCu, g%mask2dCv)
  call rotate_array_pair(g_in%areaCu, g_in%areaCv, turns, &
      g%areaCu, g%areaCv)
  call rotate_array_pair(g_in%IareaCu, g_in%IareaCv, turns, &
      g%IareaCu, g%IareaCv)

  ! Vertex point
  call rotate_array(g_in%geoLonBu, turns, g%geoLonBu)
  call rotate_array(g_in%geoLatBu, turns, g%geoLatBu)
  call rotate_array_pair(g_in%dxBu, g_in%dyBu, turns, g%dxBu, g%dyBu)
  call rotate_array(g_in%areaBu, turns, g%areaBu)
  call rotate_array(g_in%CoriolisBu, turns, g%CoriolisBu)
  call rotate_array(g_in%mask2dBu, turns, g%mask2dBu)

  ! Topographic
  g%bathymetry_at_vel = g_in%bathymetry_at_vel
  if (g%bathymetry_at_vel) then
    call rotate_array_pair(g_in%Dblock_u, g_in%Dblock_v, turns, &
        g%Dblock_u, g%Dblock_v)
    call rotate_array_pair(g_in%Dopen_u, g_in%Dopen_v, turns, &
        g%Dopen_u, g%Dopen_v)
  endif

  ! Nominal grid axes
  ! TODO: We should not assign lat values to the lon axis, and vice versa.
  !   We temporarily copy lat <-> lon since several components still expect
  !   lat and lon sizes to match the first and second dimension sizes.
  !   But we ought to instead leave them unchanged and adjust the references to
  !   these axes.
  if (modulo(turns, 2) /= 0) then
    g%gridLonT(:) = g_in%gridLatT(g_in%jeg:g_in%jsg:-1)
    g%gridLatT(:) = g_in%gridLonT(:)
    g%gridLonB(:) = g_in%gridLatB(g_in%jeg:(g_in%jsg-1):-1)
    g%gridLatB(:) = g_in%gridLonB(:)
  else
    g%gridLonT(:) = g_in%gridLonT(:)
    g%gridLatT(:) = g_in%gridLatT(:)
    g%gridLonB(:) = g_in%gridLonB(:)
    g%gridLatB(:) = g_in%gridLatB(:)
  endif

  g%x_axis_units = g_in%y_axis_units
  g%y_axis_units = g_in%x_axis_units
  g%south_lat = g_in%south_lat
  g%west_lon = g_in%west_lon
  g%len_lat = g_in%len_lat
  g%len_lon = g_in%len_lon

  ! Rotation-invariant fields
  g%areaT_global = g_in%areaT_global
  g%IareaT_global = g_in%IareaT_global
  g%Rad_Earth = g_in%Rad_Earth
  g%max_depth = g_in%max_depth

  call set_derived_dyn_horgrid(g, us)
end subroutine rotate_dyngrid

end module mom_transcribe_grid