diff --git a/src/core/MOM_dyn_horgrid.F90 b/src/core/MOM_dyn_horgrid.F90 new file mode 100644 index 0000000000..17bb0bbcb7 --- /dev/null +++ b/src/core/MOM_dyn_horgrid.F90 @@ -0,0 +1,349 @@ +module MOM_dyn_horgrid + +!*********************************************************************** +!* GNU General Public License * +!* This file is a part of MOM. * +!* * +!* MOM is free software; you can redistribute it and/or modify it and * +!* are expected to follow the terms of the GNU General Public License * +!* as published by the Free Software Foundation; either version 2 of * +!* the License, or (at your option) any later version. * +!* * +!* MOM is distributed in the hope that it will be useful, but WITHOUT * +!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * +!* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public * +!* License for more details. * +!* * +!* For the full text of the GNU General Public License, * +!* write to: Free Software Foundation, Inc., * +!* 675 Mass Ave, Cambridge, MA 02139, USA. * +!* or see: http://www.gnu.org/licenses/gpl.html * +!*********************************************************************** + +use MOM_hor_index, only : hor_index_type +use MOM_domains, only : MOM_domain_type +use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL + +implicit none ; private + +public create_dyn_horgrid, destroy_dyn_horgrid, set_derived_dyn_horgrid + +type, public :: dyn_horgrid_type + type(MOM_domain_type), pointer :: Domain => NULL() + type(MOM_domain_type), pointer :: Domain_aux => NULL() ! A non-symmetric auxiliary domain type. + + ! These elements can be copied from a provided hor_index_type. + type(hor_index_type), pointer :: HI + integer :: isc, iec, jsc, jec ! The range of the computational domain indices + integer :: isd, ied, jsd, jed ! and data domain indices at tracer cell centers. + integer :: isg, ieg, jsg, jeg ! The range of the global domain tracer cell indices. + integer :: IscB, IecB, JscB, JecB ! The range of the computational domain indices + integer :: IsdB, IedB, JsdB, JedB ! and data domain indices at tracer cell vertices. + integer :: IsgB, IegB, JsgB, JegB ! The range of the global domain vertex indices. + integer :: isd_global ! The values of isd and jsd in the global + integer :: jsd_global ! (decomposition invariant) index space. + integer :: idg_offset ! The offset between the corresponding global + integer :: jdg_offset ! and local array indices. + logical :: symmetric ! True if symmetric memory is used. + + logical :: nonblocking_updates ! If true, non-blocking halo updates are + ! allowed. The default is .false. (for now). + integer :: first_direction ! An integer that indicates which direction is + ! to be updated first in directionally split + ! parts of the calculation. This can be altered + ! during the course of the run via calls to + ! set_first_direction. + + real, allocatable, dimension(:,:) :: & + mask2dT, & ! 0 for land points and 1 for ocean points on the h-grid. Nd. + geoLatT, & ! The geographic latitude at q points in degrees of latitude or m. + geoLonT, & ! The geographic longitude at q points in degrees of longitude or m. + dxT, IdxT, & ! dxT is delta x at h points, in m, and IdxT is 1/dxT in m-1. + dyT, IdyT, & ! dyT is delta y at h points, in m, and IdyT is 1/dyT in m-1. + areaT, & ! areaT is the area of an h-cell, in m2. + IareaT, & ! IareaT = 1/areaT, in m-2. + sin_rot, & ! The sine and cosine of the angular rotation between the local + cos_rot ! model grid's northward and the true northward directions. + + real, allocatable, dimension(:,:) :: & + mask2dCu, & ! 0 for boundary points and 1 for ocean points on the u grid. Nondim. + geoLatCu, & ! The geographic latitude at u points in degrees of latitude or m. + geoLonCu, & ! The geographic longitude at u points in degrees of longitude or m. + dxCu, IdxCu, & ! dxCu is delta x at u points, in m, and IdxCu is 1/dxCu in m-1. + dyCu, IdyCu, & ! dyCu is delta y at u points, in m, and IdyCu is 1/dyCu in m-1. + dy_Cu, & ! The unblocked lengths of the u-faces of the h-cell in m. + dy_Cu_obc, & ! The unblocked lengths of the u-faces of the h-cell in m for OBC. + IareaCu, & ! The masked inverse areas of u-grid cells in m2. + areaCu ! The areas of the u-grid cells in m2. + + real, allocatable, dimension(:,:) :: & + mask2dCv, & ! 0 for boundary points and 1 for ocean points on the v grid. Nondim. + geoLatCv, & ! The geographic latitude at v points in degrees of latitude or m. + geoLonCv, & ! The geographic longitude at v points in degrees of longitude or m. + dxCv, IdxCv, & ! dxCv is delta x at v points, in m, and IdxCv is 1/dxCv in m-1. + dyCv, IdyCv, & ! dyCv is delta y at v points, in m, and IdyCv is 1/dyCv in m-1. + dx_Cv, & ! The unblocked lengths of the v-faces of the h-cell in m. + dx_Cv_obc, & ! The unblocked lengths of the v-faces of the h-cell in m for OBC. + IareaCv, & ! The masked inverse areas of v-grid cells in m2. + areaCv ! The areas of the v-grid cells in m2. + + real, allocatable, dimension(:,:) :: & + mask2dBu, & ! 0 for boundary points and 1 for ocean points on the q grid. Nondim. + geoLatBu, & ! The geographic latitude at q points in degrees of latitude or m. + geoLonBu, & ! The geographic longitude at q points in degrees of longitude or m. + dxBu, IdxBu, & ! dxBu is delta x at q points, in m, and IdxBu is 1/dxBu in m-1. + dyBu, IdyBu, & ! dyBu is delta y at q points, in m, and IdyBu is 1/dyBu in m-1. + areaBu, & ! areaBu is the area of a q-cell, in m2 + IareaBu ! IareaBu = 1/areaBu in m-2. + + real, pointer, dimension(:) :: & + gridLatT => NULL(), gridLatB => NULL() ! The latitude of T or B points for + ! the purpose of labeling the output axes. + ! On many grids these are the same as geoLatT & geoLatBu. + real, pointer, dimension(:) :: & + gridLonT => NULL(), gridLonB => NULL() ! The longitude of T or B points for + ! the purpose of labeling the output axes. + ! On many grids these are the same as geoLonT & geoLonBu. + character(len=40) :: & + x_axis_units, & ! The units that are used in labeling the coordinate + y_axis_units ! axes. Except on a Cartesian grid, these are usually + ! some variant of "degrees". + + real, allocatable, dimension(:,:) :: & + bathyT ! Ocean bottom depth at tracer points, in m. + + logical :: bathymetry_at_vel ! If true, there are separate values for the + ! basin depths at velocity points. Otherwise the effects of + ! of topography are entirely determined from thickness points. + real, allocatable, dimension(:,:) :: & + Dblock_u, & ! Topographic depths at u-points at which the flow is blocked + Dopen_u ! (Dblock_u) and open at width dy_Cu (Dopen_u), both in m. + real, allocatable, dimension(:,:) :: & + Dblock_v, & ! Topographic depths at v-points at which the flow is blocked + Dopen_v ! (Dblock_v) and open at width dx_Cv (Dopen_v), both in m. + real, allocatable, dimension(:,:) :: & + CoriolisBu ! The Coriolis parameter at corner points, in s-1. + real, allocatable, dimension(:,:) :: & + dF_dx, dF_dy ! Derivatives of f (Coriolis parameter) at h-points, in s-1 m-1. + real :: g_Earth ! The gravitational acceleration in m s-2. + + ! These variables are global sums that are useful for 1-d diagnostics + real :: areaT_global ! Global sum of h-cell area in m2 + real :: IareaT_global ! Global sum of inverse h-cell area (1/areaT_global) + ! in m2 + + ! These parameters are run-time parameters that are used during some + ! initialization routines (but not all) + real :: south_lat ! The latitude (or y-coordinate) of the first v-line + real :: west_lon ! The longitude (or x-coordinate) of the first u-line + real :: len_lat = 0. ! The latitudinal (or y-coord) extent of physical domain + real :: len_lon = 0. ! The longitudinal (or x-coord) extent of physical domain + real :: Rad_Earth = 6.378e6 ! The radius of the planet in meters. + real :: max_depth ! The maximum depth of the ocean in meters. +end type dyn_horgrid_type + +contains + +!--------------------------------------------------------------------- +!> Allocate memory used by the dyn_horgrid_type and related structures. +subroutine create_dyn_horgrid(G, HI) + type(dyn_horgrid_type), pointer :: G !< A pointer to the dynamic horizontal grid type + type(hor_index_type), intent(in) :: HI !< A hor_index_type for array extents + integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isg, ieg, jsg, jeg + + ! This subroutine allocates the lateral elements of the dyn_horgrid_type that + ! are always used and zeros them out. + + if (associated(G)) then + call MOM_error(FATAL, "destroy_dyn_horgrid called with an unassociated horgrid_type.") + endif + allocate(G) + + G%HI = HI + + G%isc = HI%isc ; G%iec = HI%iec ; G%jsc = HI%jsc ; G%jec = HI%jec + G%isd = HI%isd ; G%ied = HI%ied ; G%jsd = HI%jsd ; G%jed = HI%jed + G%isg = HI%isg ; G%ieg = HI%ieg ; G%jsg = HI%jsg ; G%jeg = HI%jeg + + G%IscB = HI%IscB ; G%IecB = HI%IecB ; G%JscB = HI%JscB ; G%JecB = HI%JecB + G%IsdB = HI%IsdB ; G%IedB = HI%IedB ; G%JsdB = HI%JsdB ; G%JedB = HI%JedB + G%IsgB = HI%IsgB ; G%IegB = HI%IegB ; G%JsgB = HI%JsgB ; G%JegB = HI%JegB + + G%idg_offset = HI%idg_offset ; G%jdg_offset = HI%jdg_offset + G%isd_global = G%isd + HI%idg_offset ; G%jsd_global = G%jsd + HI%jdg_offset + G%symmetric = HI%symmetric + + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + isg = G%isg ; ieg = G%ieg ; jsg = G%jsg ; jeg = G%jeg + + allocate(G%dxT(isd:ied,jsd:jed)) ; G%dxT(:,:) = 0.0 + allocate(G%dxCu(IsdB:IedB,jsd:jed)) ; G%dxCu(:,:) = 0.0 + allocate(G%dxCv(isd:ied,JsdB:JedB)) ; G%dxCv(:,:) = 0.0 + allocate(G%dxBu(IsdB:IedB,JsdB:JedB)) ; G%dxBu(:,:) = 0.0 + allocate(G%IdxT(isd:ied,jsd:jed)) ; G%IdxT(:,:) = 0.0 + allocate(G%IdxCu(IsdB:IedB,jsd:jed)) ; G%IdxCu(:,:) = 0.0 + allocate(G%IdxCv(isd:ied,JsdB:JedB)) ; G%IdxCv(:,:) = 0.0 + allocate(G%IdxBu(IsdB:IedB,JsdB:JedB)) ; G%IdxBu(:,:) = 0.0 + + allocate(G%dyT(isd:ied,jsd:jed)) ; G%dyT(:,:) = 0.0 + allocate(G%dyCu(IsdB:IedB,jsd:jed)) ; G%dyCu(:,:) = 0.0 + allocate(G%dyCv(isd:ied,JsdB:JedB)) ; G%dyCv(:,:) = 0.0 + allocate(G%dyBu(IsdB:IedB,JsdB:JedB)) ; G%dyBu(:,:) = 0.0 + allocate(G%IdyT(isd:ied,jsd:jed)) ; G%IdyT(:,:) = 0.0 + allocate(G%IdyCu(IsdB:IedB,jsd:jed)) ; G%IdyCu(:,:) = 0.0 + allocate(G%IdyCv(isd:ied,JsdB:JedB)) ; G%IdyCv(:,:) = 0.0 + allocate(G%IdyBu(IsdB:IedB,JsdB:JedB)) ; G%IdyBu(:,:) = 0.0 + + allocate(G%areaT(isd:ied,jsd:jed)) ; G%areaT(:,:) = 0.0 + allocate(G%IareaT(isd:ied,jsd:jed)) ; G%IareaT(:,:) = 0.0 + allocate(G%areaBu(IsdB:IedB,JsdB:JedB)) ; G%areaBu(:,:) = 0.0 + allocate(G%IareaBu(IsdB:IedB,JsdB:JedB)) ; G%IareaBu(:,:) = 0.0 + + allocate(G%mask2dT(isd:ied,jsd:jed)) ; G%mask2dT(:,:) = 0.0 + allocate(G%mask2dCu(IsdB:IedB,jsd:jed)) ; G%mask2dCu(:,:) = 0.0 + allocate(G%mask2dCv(isd:ied,JsdB:JedB)) ; G%mask2dCv(:,:) = 0.0 + allocate(G%mask2dBu(IsdB:IedB,JsdB:JedB)) ; G%mask2dBu(:,:) = 0.0 + allocate(G%geoLatT(isd:ied,jsd:jed)) ; G%geoLatT(:,:) = 0.0 + allocate(G%geoLatCu(IsdB:IedB,jsd:jed)) ; G%geoLatCu(:,:) = 0.0 + allocate(G%geoLatCv(isd:ied,JsdB:JedB)) ; G%geoLatCv(:,:) = 0.0 + allocate(G%geoLatBu(IsdB:IedB,JsdB:JedB)) ; G%geoLatBu(:,:) = 0.0 + allocate(G%geoLonT(isd:ied,jsd:jed)) ; G%geoLonT(:,:) = 0.0 + allocate(G%geoLonCu(IsdB:IedB,jsd:jed)) ; G%geoLonCu(:,:) = 0.0 + allocate(G%geoLonCv(isd:ied,JsdB:JedB)) ; G%geoLonCv(:,:) = 0.0 + allocate(G%geoLonBu(IsdB:IedB,JsdB:JedB)) ; G%geoLonBu(:,:) = 0.0 + + allocate(G%dx_Cv(isd:ied,JsdB:JedB)) ; G%dx_Cv(:,:) = 0.0 + allocate(G%dy_Cu(IsdB:IedB,jsd:jed)) ; G%dy_Cu(:,:) = 0.0 + allocate(G%dx_Cv_obc(isd:ied,JsdB:JedB)) ; G%dx_Cv_obc(:,:) = 0.0 + allocate(G%dy_Cu_obc(IsdB:IedB,jsd:jed)) ; G%dy_Cu_obc(:,:) = 0.0 + + allocate(G%areaCu(IsdB:IedB,jsd:jed)) ; G%areaCu(:,:) = 0.0 + allocate(G%areaCv(isd:ied,JsdB:JedB)) ; G%areaCv(:,:) = 0.0 + allocate(G%IareaCu(IsdB:IedB,jsd:jed)) ; G%IareaCu(:,:) = 0.0 + allocate(G%IareaCv(isd:ied,JsdB:JedB)) ; G%IareaCv(:,:) = 0.0 + + allocate(G%bathyT(isd:ied, jsd:jed)) ; G%bathyT(:,:) = 0.0 + allocate(G%CoriolisBu(IsdB:IedB, JsdB:JedB)) ; G%CoriolisBu(:,:) = 0.0 + allocate(G%dF_dx(isd:ied, jsd:jed)) ; G%dF_dx(:,:) = 0.0 + allocate(G%dF_dy(isd:ied, jsd:jed)) ; G%dF_dy(:,:) = 0.0 + + allocate(G%sin_rot(isd:ied,jsd:jed)) ; G%sin_rot(:,:) = 0.0 + allocate(G%cos_rot(isd:ied,jsd:jed)) ; G%cos_rot(:,:) = 1.0 + + allocate(G%Dblock_u(IsdB:IedB, jsd:jed)) ; G%Dblock_u(:,:) = 0.0 + allocate(G%Dopen_u(IsdB:IedB, jsd:jed)) ; G%Dopen_u(:,:) = 0.0 + allocate(G%Dblock_v(isd:ied, JsdB:JedB)) ; G%Dblock_v(:,:) = 0.0 + allocate(G%Dopen_v(isd:ied, JsdB:JedB)) ; G%Dopen_v(:,:) = 0.0 + + allocate(G%gridLonT(isg:ieg)) ; G%gridLonT(:) = 0.0 + allocate(G%gridLonB(G%IsgB:G%IegB)) ; G%gridLonB(:) = 0.0 + allocate(G%gridLatT(jsg:jeg)) ; G%gridLatT(:) = 0.0 + allocate(G%gridLatB(G%JsgB:G%JegB)) ; G%gridLatB(:) = 0.0 + +end subroutine create_dyn_horgrid + +!> set_derived_dyn_horgrid calculates metric terms that are derived from other metrics. +subroutine set_derived_dyn_horgrid(G) + type(dyn_horgrid_type), pointer :: G !< The dynamic horizontal grid type +! Various inverse grid spacings and derived areas are calculated within this +! subroutine. + integer :: i, j, isd, ied, jsd, jed + integer :: IsdB, IedB, JsdB, JedB + + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + + do j=jsd,jed ; do i=isd,ied + if (G%dxT(i,j) < 0.0) G%dxT(i,j) = 0.0 + if (G%dyT(i,j) < 0.0) G%dyT(i,j) = 0.0 + G%IdxT(i,j) = Adcroft_reciprocal(G%dxT(i,j)) + G%IdyT(i,j) = Adcroft_reciprocal(G%dyT(i,j)) + G%IareaT(i,j) = Adcroft_reciprocal(G%areaT(i,j)) + enddo ; enddo + + do j=jsd,jed ; do I=IsdB,IedB + if (G%dxCu(I,j) < 0.0) G%dxCu(I,j) = 0.0 + if (G%dyCu(I,j) < 0.0) G%dyCu(I,j) = 0.0 + G%IdxCu(I,j) = Adcroft_reciprocal(G%dxCu(I,j)) + G%IdyCu(I,j) = Adcroft_reciprocal(G%dyCu(I,j)) + enddo ; enddo + + do J=JsdB,JedB ; do i=isd,ied + if (G%dxCv(i,J) < 0.0) G%dxCv(i,J) = 0.0 + if (G%dyCv(i,J) < 0.0) G%dyCv(i,J) = 0.0 + G%IdxCv(i,J) = Adcroft_reciprocal(G%dxCv(i,J)) + G%IdyCv(i,J) = Adcroft_reciprocal(G%dyCv(i,J)) + enddo ; enddo + + do J=JsdB,JedB ; do I=IsdB,IedB + if (G%dxBu(I,J) < 0.0) G%dxBu(I,J) = 0.0 + if (G%dyBu(I,J) < 0.0) G%dyBu(I,J) = 0.0 + + G%IdxBu(I,J) = Adcroft_reciprocal(G%dxBu(I,J)) + G%IdyBu(I,J) = Adcroft_reciprocal(G%dyBu(I,J)) + ! areaBu has usually been set to a positive area elsewhere. + if (G%areaBu(I,J) <= 0.0) G%areaBu(I,J) = G%dxBu(I,J) * G%dyBu(I,J) + G%IareaBu(I,J) = Adcroft_reciprocal(G%areaBu(I,J)) + enddo ; enddo + +end subroutine set_derived_dyn_horgrid + +!> Adcroft_reciprocal(x) = 1/x for |x|>0 or 0 for x=0. +function Adcroft_reciprocal(val) result(I_val) + real, intent(in) :: val !< The value being inverted. + real :: I_val !< The Adcroft reciprocal of val. + + I_val = 0.0 ; if (val /= 0.0) I_val = 1.0/val +end function Adcroft_reciprocal + +!--------------------------------------------------------------------- +!> Release memory used by the dyn_horgrid_type and related structures. +subroutine destroy_dyn_horgrid(G) + type(dyn_horgrid_type), pointer :: G !< The dynamic horizontal grid type + + if (.not.associated(G)) then + call MOM_error(FATAL, "destroy_dyn_horgrid called with an unassociated horgrid_type.") + endif + + deallocate(G%dxT) ; deallocate(G%dxCu) ; deallocate(G%dxCv) ; deallocate(G%dxBu) + deallocate(G%IdxT) ; deallocate(G%IdxCu) ; deallocate(G%IdxCv) ; deallocate(G%IdxBu) + + deallocate(G%dyT) ; deallocate(G%dyCu) ; deallocate(G%dyCv) ; deallocate(G%dyBu) + deallocate(G%IdyT) ; deallocate(G%IdyCu) ; deallocate(G%IdyCv) ; deallocate(G%IdyBu) + + deallocate(G%areaT) ; deallocate(G%IareaT) + deallocate(G%areaBu) ; deallocate(G%IareaBu) + deallocate(G%areaCu) ; deallocate(G%IareaCu) + deallocate(G%areaCv) ; deallocate(G%IareaCv) + + deallocate(G%mask2dT) ; deallocate(G%mask2dCu) + deallocate(G%mask2dCv) ; deallocate(G%mask2dBu) + + deallocate(G%geoLatT) ; deallocate(G%geoLatCu) + deallocate(G%geoLatCv) ; deallocate(G%geoLatBu) + deallocate(G%geoLonT) ; deallocate(G%geoLonCu) + deallocate(G%geoLonCv) ; deallocate(G%geoLonBu) + + deallocate(G%dx_Cv) ; deallocate(G%dy_Cu) + deallocate(G%dx_Cv_obc) ; deallocate(G%dy_Cu_obc) + + deallocate(G%bathyT) ; deallocate(G%CoriolisBu) + deallocate(G%dF_dx) ; deallocate(G%dF_dy) + deallocate(G%sin_rot) ; deallocate(G%cos_rot) + + deallocate(G%Dblock_u) ; deallocate(G%Dopen_u) + deallocate(G%Dblock_v) ; deallocate(G%Dopen_v) + + deallocate(G%gridLonT) ; deallocate(G%gridLatT) + deallocate(G%gridLonB) ; deallocate(G%gridLatB) + + ! Do not deallocate G%Domain%mpp_domain or deallocate(G%Domain) because + ! these are pointers to types used elsewhere. + G%Domain => NULL() + + deallocate(G) + +end subroutine destroy_dyn_horgrid + +end module MOM_dyn_horgrid diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index bf978cac76..fd8f252eb4 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -29,7 +29,7 @@ module MOM_grid #include -public MOM_grid_init, MOM_grid_end, set_first_direction +public MOM_grid_init, MOM_grid_end, set_derived_metrics, set_first_direction public isPointInCell, hor_index_type type, public :: ocean_grid_type @@ -320,6 +320,60 @@ subroutine MOM_grid_init(G, param_file) end subroutine MOM_grid_init + +!> set_derived_metrics calculates metric terms that are derived from other metrics. +subroutine set_derived_metrics(G) + type(ocean_grid_type), intent(inout) :: G !< The horizontal grid structure +! Various inverse grid spacings and derived areas are calculated within this +! subroutine. + integer :: i, j, isd, ied, jsd, jed + integer :: IsdB, IedB, JsdB, JedB + + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + + do j=jsd,jed ; do i=isd,ied + if (G%dxT(i,j) < 0.0) G%dxT(i,j) = 0.0 + if (G%dyT(i,j) < 0.0) G%dyT(i,j) = 0.0 + G%IdxT(i,j) = Adcroft_reciprocal(G%dxT(i,j)) + G%IdyT(i,j) = Adcroft_reciprocal(G%dyT(i,j)) + G%IareaT(i,j) = Adcroft_reciprocal(G%areaT(i,j)) + enddo ; enddo + + do j=jsd,jed ; do I=IsdB,IedB + if (G%dxCu(I,j) < 0.0) G%dxCu(I,j) = 0.0 + if (G%dyCu(I,j) < 0.0) G%dyCu(I,j) = 0.0 + G%IdxCu(I,j) = Adcroft_reciprocal(G%dxCu(I,j)) + G%IdyCu(I,j) = Adcroft_reciprocal(G%dyCu(I,j)) + enddo ; enddo + + do J=JsdB,JedB ; do i=isd,ied + if (G%dxCv(i,J) < 0.0) G%dxCv(i,J) = 0.0 + if (G%dyCv(i,J) < 0.0) G%dyCv(i,J) = 0.0 + G%IdxCv(i,J) = Adcroft_reciprocal(G%dxCv(i,J)) + G%IdyCv(i,J) = Adcroft_reciprocal(G%dyCv(i,J)) + enddo ; enddo + + do J=JsdB,JedB ; do I=IsdB,IedB + if (G%dxBu(I,J) < 0.0) G%dxBu(I,J) = 0.0 + if (G%dyBu(I,J) < 0.0) G%dyBu(I,J) = 0.0 + + G%IdxBu(I,J) = Adcroft_reciprocal(G%dxBu(I,J)) + G%IdyBu(I,J) = Adcroft_reciprocal(G%dyBu(I,J)) + ! areaBu has usually been set to a positive area elsewhere. + if (G%areaBu(I,J) <= 0.0) G%areaBu(I,J) = G%dxBu(I,J) * G%dyBu(I,J) + G%IareaBu(I,J) = Adcroft_reciprocal(G%areaBu(I,J)) + enddo ; enddo +end subroutine set_derived_metrics + +!> Adcroft_reciprocal(x) = 1/x for |x|>0 or 0 for x=0. +function Adcroft_reciprocal(val) result(I_val) + real, intent(in) :: val !< The value being inverted. + real :: I_val !< The Adcroft reciprocal of val. + + I_val = 0.0 ; if (val /= 0.0) I_val = 1.0/val +end function Adcroft_reciprocal + !> Returns true if the coordinates (x,y) are within the h-cell (i,j) logical function isPointInCell(G, i, j, x, y) type(ocean_grid_type), intent(in) :: G !< Grid type diff --git a/src/core/MOM_transcribe_grid.F90 b/src/core/MOM_transcribe_grid.F90 new file mode 100644 index 0000000000..db8aa57c50 --- /dev/null +++ b/src/core/MOM_transcribe_grid.F90 @@ -0,0 +1,190 @@ +module MOM_transcribe_grid + +!*********************************************************************** +!* GNU General Public License * +!* This file is a part of MOM. * +!* * +!* MOM is free software; you can redistribute it and/or modify it and * +!* are expected to follow the terms of the GNU General Public License * +!* as published by the Free Software Foundation; either version 2 of * +!* the License, or (at your option) any later version. * +!* * +!* MOM is distributed in the hope that it will be useful, but WITHOUT * +!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY * +!* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public * +!* License for more details. * +!* * +!* For the full text of the GNU General Public License, * +!* write to: Free Software Foundation, Inc., * +!* 675 Mass Ave, Cambridge, MA 02139, USA. * +!* or see: http://www.gnu.org/licenses/gpl.html * +!*********************************************************************** + +use MOM_hor_index, only : hor_index_type, hor_index_init +use MOM_domains, only : MOM_domain_type, get_domain_extent, compute_block_extent +use MOM_domains, only : pass_var, pass_vector, sum_across_PEs, broadcast +use MOM_domains, only : root_PE, To_All, SCALAR_PAIR, CGRID_NE, AGRID, BGRID_NE, CORNER +use MOM_dyn_horgrid, only : dyn_horgrid_type +use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_grid, only : ocean_grid_type, set_derived_metrics +! use MOM_grid, only : ! , set_first_direction + +implicit none ; private + +#include + +contains + +!> Copies information from a dynamic (shared) horizontal grid type into an +!! ocean_grid_type. +subroutine copy_dyngrid_to_MOM_grid(dG, G) + type(dyn_horgrid_type), intent(in) :: dG !< Common horizontal grid type + type(ocean_grid_type), intent(inout) :: G !< Ocean grid type + + integer :: isd, ied, jsd, jed, nk ! Common data domains. + integer :: IsdB, IedB, JsdB, JedB ! Common data domains. + integer :: ido, jdo ! Indexing offsets between the grids. + integer :: i, j, k + + ! MOM_grid_init and create_dyn_horgrid are called outside of this routine. + ! This routine copies over the fields that were set by MOM_initilized_fixed. + + ! Determine the indexing offsets between the grids. + ido = dG%idg_offset - G%idg_offset + jdo = dG%jdg_offset - G%jdg_offset + + isd = max(G%isd, dG%isd+ido) ; jsd = max(G%jsd, dG%jsd+jdo) + ied = min(G%ied, dG%ied+ido) ; jed = min(G%jed, dG%jed+jdo) + IsdB = max(G%IsdB, dG%IsdB+ido) ; JsdB = max(G%JsdB, dG%JsdB+jdo) + IedB = min(G%IedB, dG%IedB+ido) ; JedB = min(G%JedB, dG%JedB+jdo) + + ! Check that the grids conform. + if ((isd > G%isc) .or. (ied < G%ied) .or. (jsd > G%jsc) .or. (jed > G%jed)) & + call MOM_error(FATAL, "copy_dyngrid_to_MOM_grid called with incompatible grids.") + + do i=isd,ied ; do j=jsd,jed + G%geoLonT(i,j) = dG%geoLonT(i+ido,j+jdo) + G%geoLatT(i,j) = dG%geoLatT(i+ido,j+jdo) + G%dxT(i,j) = dG%dxT(i+ido,j+jdo) + G%dyT(i,j) = dG%dyT(i+ido,j+jdo) + G%areaT(i,j) = dG%areaT(i+ido,j+jdo) + G%bathyT(i,j) = dG%bathyT(i+ido,j+jdo) + enddo ; enddo + + do I=IsdB,IedB ; do j=jsd,jed + G%geoLonCu(I,j) = dG%geoLonCu(I+ido,j+jdo) + G%geoLatCu(I,j) = dG%geoLatCu(I+ido,j+jdo) + G%dxCu(I,j) = dG%dxCu(I+ido,j+jdo) + G%dyCu(I,j) = dG%dyCu(I+ido,j+jdo) + enddo ; enddo + + do i=isd,ied ; do J=JsdB,JedB + G%geoLonCv(i,J) = dG%geoLonCv(i+ido,J+jdo) + G%geoLatCv(i,J) = dG%geoLatCv(i+ido,J+jdo) + G%dxCv(i,J) = dG%dxCv(i+ido,J+jdo) + G%dyCv(i,J) = dG%dyCv(i+ido,J+jdo) + enddo ; enddo + + do I=IsdB,IedB ; do J=JsdB,JedB + G%geoLonBu(I,J) = dG%geoLonBu(I+ido,J+jdo) + G%geoLatBu(I,J) = dG%geoLatBu(I+ido,J+jdo) + G%dxBu(I,J) = dG%dxBu(I+ido,J+jdo) + G%dyBu(I,J) = dG%dyBu(I+ido,J+jdo) + G%areaBu(i,j) = dG%areaBu(I+ido,J+jdo) + enddo ; enddo + +! Update the halos in case the dynamic grid has smaller halos than the ocean grid. + call pass_vector(G%dyCu, G%dxCv, G%Domain, To_All+Scalar_Pair, CGRID_NE) + call pass_vector(G%dxCu, G%dyCv, G%Domain, To_All+Scalar_Pair, CGRID_NE) + call pass_vector(G%dxBu, G%dyBu, G%Domain, To_All+Scalar_Pair, BGRID_NE) + call pass_var(G%areaT, G%Domain) + call pass_var(G%areaBu, G%Domain, position=CORNER) + + call set_derived_metrics(G) + +end subroutine copy_dyngrid_to_MOM_grid + +! subroutine set_first_direction(G, y_first) +! type(ocean_grid_type), intent(inout) :: G +! integer, intent(in) :: y_first + +! G%first_direction = y_first +! end subroutine set_first_direction + +!--------------------------------------------------------------------- +!> Allocate memory used by the ocean_grid_type and related structures. +subroutine allocate_metrics(G) + type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type + integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isg, ieg, jsg, jeg + + ! This subroutine allocates the lateral elements of the ocean_grid_type that + ! are always used and zeros them out. + + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + isg = G%isg ; ieg = G%ieg ; jsg = G%jsg ; jeg = G%jeg + + ALLOC_(G%dxT(isd:ied,jsd:jed)) ; G%dxT(:,:) = 0.0 + ALLOC_(G%dxCu(IsdB:IedB,jsd:jed)) ; G%dxCu(:,:) = 0.0 + ALLOC_(G%dxCv(isd:ied,JsdB:JedB)) ; G%dxCv(:,:) = 0.0 + ALLOC_(G%dxBu(IsdB:IedB,JsdB:JedB)) ; G%dxBu(:,:) = 0.0 + ALLOC_(G%IdxT(isd:ied,jsd:jed)) ; G%IdxT(:,:) = 0.0 + ALLOC_(G%IdxCu(IsdB:IedB,jsd:jed)) ; G%IdxCu(:,:) = 0.0 + ALLOC_(G%IdxCv(isd:ied,JsdB:JedB)) ; G%IdxCv(:,:) = 0.0 + ALLOC_(G%IdxBu(IsdB:IedB,JsdB:JedB)) ; G%IdxBu(:,:) = 0.0 + + ALLOC_(G%dyT(isd:ied,jsd:jed)) ; G%dyT(:,:) = 0.0 + ALLOC_(G%dyCu(IsdB:IedB,jsd:jed)) ; G%dyCu(:,:) = 0.0 + ALLOC_(G%dyCv(isd:ied,JsdB:JedB)) ; G%dyCv(:,:) = 0.0 + ALLOC_(G%dyBu(IsdB:IedB,JsdB:JedB)) ; G%dyBu(:,:) = 0.0 + ALLOC_(G%IdyT(isd:ied,jsd:jed)) ; G%IdyT(:,:) = 0.0 + ALLOC_(G%IdyCu(IsdB:IedB,jsd:jed)) ; G%IdyCu(:,:) = 0.0 + ALLOC_(G%IdyCv(isd:ied,JsdB:JedB)) ; G%IdyCv(:,:) = 0.0 + ALLOC_(G%IdyBu(IsdB:IedB,JsdB:JedB)) ; G%IdyBu(:,:) = 0.0 + + ALLOC_(G%areaT(isd:ied,jsd:jed)) ; G%areaT(:,:) = 0.0 + ALLOC_(G%IareaT(isd:ied,jsd:jed)) ; G%IareaT(:,:) = 0.0 + ALLOC_(G%areaBu(IsdB:IedB,JsdB:JedB)) ; G%areaBu(:,:) = 0.0 + ALLOC_(G%IareaBu(IsdB:IedB,JsdB:JedB)) ; G%IareaBu(:,:) = 0.0 + + ALLOC_(G%mask2dT(isd:ied,jsd:jed)) ; G%mask2dT(:,:) = 0.0 + ALLOC_(G%mask2dCu(IsdB:IedB,jsd:jed)) ; G%mask2dCu(:,:) = 0.0 + ALLOC_(G%mask2dCv(isd:ied,JsdB:JedB)) ; G%mask2dCv(:,:) = 0.0 + ALLOC_(G%mask2dBu(IsdB:IedB,JsdB:JedB)) ; G%mask2dBu(:,:) = 0.0 + ALLOC_(G%geoLatT(isd:ied,jsd:jed)) ; G%geoLatT(:,:) = 0.0 + ALLOC_(G%geoLatCu(IsdB:IedB,jsd:jed)) ; G%geoLatCu(:,:) = 0.0 + ALLOC_(G%geoLatCv(isd:ied,JsdB:JedB)) ; G%geoLatCv(:,:) = 0.0 + ALLOC_(G%geoLatBu(IsdB:IedB,JsdB:JedB)) ; G%geoLatBu(:,:) = 0.0 + ALLOC_(G%geoLonT(isd:ied,jsd:jed)) ; G%geoLonT(:,:) = 0.0 + ALLOC_(G%geoLonCu(IsdB:IedB,jsd:jed)) ; G%geoLonCu(:,:) = 0.0 + ALLOC_(G%geoLonCv(isd:ied,JsdB:JedB)) ; G%geoLonCv(:,:) = 0.0 + ALLOC_(G%geoLonBu(IsdB:IedB,JsdB:JedB)) ; G%geoLonBu(:,:) = 0.0 + + ALLOC_(G%dx_Cv(isd:ied,JsdB:JedB)) ; G%dx_Cv(:,:) = 0.0 + ALLOC_(G%dy_Cu(IsdB:IedB,jsd:jed)) ; G%dy_Cu(:,:) = 0.0 + ALLOC_(G%dx_Cv_obc(isd:ied,JsdB:JedB)) ; G%dx_Cv_obc(:,:) = 0.0 + ALLOC_(G%dy_Cu_obc(IsdB:IedB,jsd:jed)) ; G%dy_Cu_obc(:,:) = 0.0 + + ALLOC_(G%areaCu(IsdB:IedB,jsd:jed)) ; G%areaCu(:,:) = 0.0 + ALLOC_(G%areaCv(isd:ied,JsdB:JedB)) ; G%areaCv(:,:) = 0.0 + ALLOC_(G%IareaCu(IsdB:IedB,jsd:jed)) ; G%IareaCu(:,:) = 0.0 + ALLOC_(G%IareaCv(isd:ied,JsdB:JedB)) ; G%IareaCv(:,:) = 0.0 + + ALLOC_(G%bathyT(isd:ied, jsd:jed)) ; G%bathyT(:,:) = 0.0 + ALLOC_(G%CoriolisBu(IsdB:IedB, JsdB:JedB)) ; G%CoriolisBu(:,:) = 0.0 + ALLOC_(G%dF_dx(isd:ied, jsd:jed)) ; G%dF_dx(:,:) = 0.0 + ALLOC_(G%dF_dy(isd:ied, jsd:jed)) ; G%dF_dy(:,:) = 0.0 + + ALLOC_(G%sin_rot(isd:ied,jsd:jed)) ; G%sin_rot(:,:) = 0.0 + ALLOC_(G%cos_rot(isd:ied,jsd:jed)) ; G%cos_rot(:,:) = 1.0 + + allocate(G%gridLonT(isg:ieg)) ; G%gridLonT(:) = 0.0 + allocate(G%gridLonB(G%IsgB:G%IegB)) ; G%gridLonB(:) = 0.0 + allocate(G%gridLatT(jsg:jeg)) ; G%gridLatT(:) = 0.0 + allocate(G%gridLatB(G%JsgB:G%JegB)) ; G%gridLatB(:) = 0.0 + +end subroutine allocate_metrics + + +end module MOM_transcribe_grid diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index 511ec52888..59632214c7 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -75,7 +75,7 @@ module MOM_grid_initialize use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_grid, only : ocean_grid_type +use MOM_grid, only : ocean_grid_type, set_derived_metrics use MOM_io, only : read_data, slasher, file_exists use MOM_io, only : CORNER, NORTH_FACE, EAST_FACE @@ -154,7 +154,9 @@ subroutine set_grid_metrics(G, param_file) end select ! Calculate derived metrics (i.e. reciprocals and products) - call set_grid_derived_metrics(G, param_file) + call callTree_enter("set_derived_metrics(), MOM_grid_initialize.F90") + call set_derived_metrics(G) + call callTree_leave("set_derived_metrics()") if (debug) call grid_metrics_chksum('MOM_grid_init/set_grid_metrics',G) @@ -163,101 +165,6 @@ end subroutine set_grid_metrics ! ------------------------------------------------------------------------------ - -!> set_grid_derived_metrics is sets additional grid metrics that can be derived -!! from the basic grid lengths and areas. -subroutine set_grid_derived_metrics(G, param_file) - type(ocean_grid_type), intent(inout) :: G !< The horizontal grid structure - type(param_file_type), intent(in) :: param_file !< Parameter file structure - -! Calculate the values of the metric terms that might be used -! and save them in arrays. -! Within this subroutine, the x- and y- grid spacings and their -! inverses and the cell areas centered on T, Bu, Cu, and Cv points are -! calculated, as are the geographic locations of each of these 4 -! sets of points. - character( len = 128) :: warnmesg - integer :: i, j, isd, ied, jsd, jed - integer :: IsdB, IedB, JsdB, JedB - - isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB - - call callTree_enter("set_grid_derived_metrics(), MOM_grid_initialize.F90") - - do j=jsd,jed ; do i=isd,ied - if (G%dxT(i,j) < 0.0) then - write(warnmesg,68) pe_here(),"dxT",i,j,G%dxT(i,j),0.0 - call MOM_mesg(warnmesg, all_print=.true.) - G%dxT(i,j) = 0.0 - endif - if (G%dyT(i,j) < 0.0) then - write(warnmesg,68) pe_here(),"dyT",i,j,G%dyT(i,j),0.0 - call MOM_mesg(warnmesg, all_print=.true.) - G%dyT(i,j) = 0.0 - endif - G%IdxT(i,j) = Adcroft_reciprocal(G%dxT(i,j)) - G%IdyT(i,j) = Adcroft_reciprocal(G%dyT(i,j)) - G%IareaT(i,j) = Adcroft_reciprocal(G%areaT(i,j)) - enddo ; enddo - - do j=jsd,jed ; do I=IsdB,IedB - if (G%dxCu(I,j) < 0.0) then - write(warnmesg,68) pe_here(),"dxCu",I,j,G%dxCu(I,j),0.0 - call MOM_mesg(warnmesg, all_print=.true.) - G%dxCu(I,j) = 0.0 - endif - if (G%dyCu(I,j) < 0.0) then - write(warnmesg,68) pe_here(),"dyCu",I,j,G%dyCu(I,j),0.0 - call MOM_mesg(warnmesg, all_print=.true.) - G%dyCu(I,j) = 0.0 - endif - G%IdxCu(I,j) = Adcroft_reciprocal(G%dxCu(I,j)) - G%IdyCu(I,j) = Adcroft_reciprocal(G%dyCu(I,j)) - enddo ; enddo - - do J=JsdB,JedB ; do i=isd,ied - if (G%dxCv(i,J) < 0.0) then - write(warnmesg,68) pe_here(),"dxCv",i,J,G%dxCv(i,J),0.0 - call MOM_mesg(warnmesg, all_print=.true.) - G%dxCv(i,J) = 0.0 - endif - if (G%dyCv(i,J) < 0.0) then - write(warnmesg,68) pe_here(),"dyCv",i,J,G%dyCv(i,J),0.0 - call MOM_mesg(warnmesg, all_print=.true.) - G%dyCv(i,J) = 0.0 - endif - G%IdxCv(i,J) = Adcroft_reciprocal(G%dxCv(i,J)) - G%IdyCv(i,J) = Adcroft_reciprocal(G%dyCv(i,J)) - enddo ; enddo - - do J=JsdB,JedB ; do I=IsdB,IedB - if (G%dxBu(I,J) < 0.0) then - write(warnmesg,68) pe_here(),"dxBu",I,J,G%dxBu(I,J),0.0 - call MOM_mesg(warnmesg, all_print=.true.) - G%dxBu(I,J) = 0.0 - endif - if (G%dyBu(I,J) < 0.0) then - write(warnmesg,68) pe_here(),"dyBu",I,J,G%dyBu(I,J),0.0 - call MOM_mesg(warnmesg, all_print=.true.) - G%dyBu(I,J) = 0.0 - endif - - G%IdxBu(I,J) = Adcroft_reciprocal(G%dxBu(I,J)) - G%IdyBu(I,J) = Adcroft_reciprocal(G%dyBu(I,J)) - ! areaBu has usually been set to a positive area elsewhere. - if (G%areaBu(I,J) <= 0.0) G%areaBu(I,J) = G%dxBu(I,J) * G%dyBu(I,J) - G%IareaBu(I,J) = Adcroft_reciprocal(G%areaBu(I,J)) - enddo ; enddo - -68 FORMAT ("WARNING: PE ",I4," ",a4,"(",I4,",",I4,") = ",ES12.4, & - " is being changed to ",ES12.4,".") - - call callTree_leave("set_grid_derived_metrics()") -end subroutine set_grid_derived_metrics - -! ------------------------------------------------------------------------------ - !> grid_metrics_chksum performs a set of checksums on metrics on the grid for !! debugging. subroutine grid_metrics_chksum(parent, G)