Merge branch 'user/ksh/open_bc' into user/ksh/open_bc_uv

src/core/MOM.F90

@@ -1424,29 +1424,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
 
   call find_obsolete_params(param_file)
 
-#ifdef SYMMETRIC_MEMORY_
-  symmetric = .true.
-#else
-  symmetric = .false.
-#endif
-#ifdef STATIC_MEMORY_
-  call MOM_domains_init(G%domain, param_file, symmetric=symmetric, &
-            static_memory=.true., NIHALO=NIHALO_, NJHALO=NJHALO_, &
-            NIGLOBAL=NIGLOBAL_, NJGLOBAL=NJGLOBAL_, NIPROC=NIPROC_, &
-            NJPROC=NJPROC_)
-#else
-  call MOM_domains_init(G%domain, param_file, symmetric=symmetric)
-#endif
-  call callTree_waypoint("domains initialized (initialize_MOM)")
-
-  call MOM_checksums_init(param_file)
-
-  call diag_mediator_infrastructure_init()
-  call MOM_io_init(param_file)
-  call MOM_grid_init(G, param_file)
-  call verticalGridInit( param_file, CS%GV )
-  GV => CS%GV
-
   ! Read relevant parameters and write them to the model log.
   call log_version(param_file, "MOM", version, "")
   call get_param(param_file, "MOM", "VERBOSITY", verbosity,  &
@@ -1656,21 +1633,53 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
   if (CS%adiabatic .and. CS%bulkmixedlayer) call MOM_error(FATAL, &
     "MOM: ADIABATIC and BULKMIXEDLAYER can not both be defined.")
 
+  ! Set up the model domain and grids.
+#ifdef SYMMETRIC_MEMORY_
+  symmetric = .true.
+#else
+  symmetric = .false.
+#endif
+#ifdef STATIC_MEMORY_
+  call MOM_domains_init(G%domain, param_file, symmetric=symmetric, &
+            static_memory=.true., NIHALO=NIHALO_, NJHALO=NJHALO_, &
+            NIGLOBAL=NIGLOBAL_, NJGLOBAL=NJGLOBAL_, NIPROC=NIPROC_, &
+            NJPROC=NJPROC_)
+#else
+  call MOM_domains_init(G%domain, param_file, symmetric=symmetric)
+#endif
+  call callTree_waypoint("domains initialized (initialize_MOM)")
+
+  call MOM_checksums_init(param_file)
+
+  call diag_mediator_infrastructure_init()
+  call MOM_io_init(param_file)
+  call MOM_grid_init(G, param_file)
+
+  call create_dyn_horgrid(dG, G%HI)
+  dG%first_direction = G%first_direction
+  dG%bathymetry_at_vel = G%bathymetry_at_vel
+  call clone_MOM_domain(G%Domain, dG%Domain)
+
+  call verticalGridInit( param_file, CS%GV )
+  GV => CS%GV
+  dG%g_Earth = GV%g_Earth
+
+
   ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes.
-  if (CS%debug .or. G%symmetric) &
-    call clone_MOM_domain(G%Domain, G%Domain_aux, symmetric=.false.)
+  if (CS%debug .or. dG%symmetric) &
+    call clone_MOM_domain(dG%Domain, dG%Domain_aux, symmetric=.false.)
 
   call MOM_timing_init(CS)
 
   call tracer_registry_init(param_file, CS%tracer_Reg)
 
   ! Copy a common variable from the vertical grid to the horizontal grid.
   ! Consider removing this later?
-  G%ke = GV%ke
+!  G%ke = GV%ke
 
-  is   = G%isc   ; ie   = G%iec  ; js   = G%jsc  ; je   = G%jec ; nz = GV%ke
-  isd  = G%isd   ; ied  = G%ied  ; jsd  = G%jsd  ; jed  = G%jed
-  IsdB = G%IsdB  ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
+  is   = dG%isc   ; ie   = dG%iec  ; js   = dG%jsc  ; je   = dG%jec ; nz = GV%ke
+  isd  = dG%isd   ; ied  = dG%ied  ; jsd  = dG%jsd  ; jed  = dG%jed
+  IsdB = dG%IsdB  ; IedB = dG%IedB ; JsdB = dG%JsdB ; JedB = dG%JedB
 
   ! Allocate and initialize space for primary MOM variables.
   ALLOC_(CS%u(IsdB:IedB,jsd:jed,nz))   ; CS%u(:,:,:) = 0.0
@@ -1785,17 +1794,28 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
   call callTree_waypoint("restart registration complete (initialize_MOM)")
 
   call cpu_clock_begin(id_clock_MOM_init)
-  call MOM_initialize_fixed(G, CS%OBC, param_file, write_geom_files, dirs%output_directory)
+  call MOM_initialize_fixed(dG, CS%OBC, param_file, write_geom_files, dirs%output_directory)
   call callTree_waypoint("returned from MOM_initialize_fixed() (initialize_MOM)")
   call MOM_initialize_coord(GV, param_file, write_geom_files, &
-                            dirs%output_directory, CS%tv, G%max_depth)
+                            dirs%output_directory, CS%tv, dG%max_depth)
   call callTree_waypoint("returned from MOM_initialize_coord() (initialize_MOM)")
 
   if (CS%use_ALE_algorithm) then
-    call ALE_init(param_file, GV, G%max_depth, CS%ALE_CSp)
+    call ALE_init(param_file, GV, dG%max_depth, CS%ALE_CSp)
     call callTree_waypoint("returned from ALE_init() (initialize_MOM)")
   endif
 
+  !   Shift from using the temporary dynamic grid type to using the final (potentially
+  ! static) ocean grid type.
+  ! call clone_MOM_domain(dG%Domain, CS%G%Domain)
+  ! call MOM_grid_init(CS%G, param_file)
+
+  call copy_dyngrid_to_MOM_grid(dg, G)
+  ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes.
+  if (CS%debug .or. G%symmetric) &
+    call clone_MOM_domain(G%Domain, G%Domain_aux, symmetric=.false.)
+  G%ke = GV%ke
+
   call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, param_file, &
                             dirs, CS%restart_CSp, CS%ALE_CSp, CS%tracer_Reg, &
                             CS%sponge_CSp, CS%ALE_sponge_CSp, CS%OBC, Time_in)
@@ -1805,20 +1825,6 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
   ! From this point, there may be pointers being set, so the final grid type
   ! that will persist through the run has to be used.
 
-  !   Shift from using the temporary dynamic grid type to using the final (potentially
-  ! static) ocean grid type.
-  ! call clone_MOM_domain(dG%Domain, CS%G%Domain)
-  ! call MOM_grid_init(CS%G, param_file)
-  ! call copy_dyngrid_to_MOM_grid(dg, CS%G)
-  ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes.
-  ! if (CS%debug .or. CS%G%symmetric) &
-  !  call clone_MOM_domain(CS%G%Domain, CS%G%Domain_aux, symmetric=.false.)
-
-  ! ! Copy a common variable from the vertical grid to the horizontal grid.
-  ! ! Consider removing this later?
-  ! CS%G%ke = GV%ke
-  ! G => CS%G
-
   if (test_grid_copy) then
     !  Copy the data from the temporary grid to the dyn_hor_grid to CS%G.
     call create_dyn_horgrid(dG, G%HI)
@@ -1885,12 +1891,12 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
 
    diag    => CS%diag
   ! Initialize the diag mediator.
-  call diag_mediator_init(G, param_file, diag, doc_file_dir=dirs%output_directory)
+  call diag_mediator_init(G, GV%ke, param_file, diag, doc_file_dir=dirs%output_directory)
 
   ! Initialize the diagnostics mask arrays.
   ! This step has to be done after call to MOM_initialize_state
   ! and before MOM_diagnostics_init
-  call diag_masks_set(G, CS%missing, diag)
+  call diag_masks_set(G, GV%ke, CS%missing, diag)
 
   ! Set up a pointers h within diag mediator control structure,
   ! this needs to occur _after_ CS%h has been allocated.

src/core/MOM_CoriolisAdv.F90

@@ -875,7 +875,7 @@ subroutine CoriolisAdv_init(Time, G, param_file, diag, AD, CS)
   CS%diag => diag ; CS%Time => Time
 
   ! Read all relevant parameters and write them to the model log.
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   call get_param(param_file, mod, "NOSLIP", CS%no_slip, &
                  "If true, no slip boundary conditions are used; otherwise \n"//&
                  "free slip boundary conditions are assumed. The \n"//&

src/core/MOM_PressureForce.F90

@@ -144,7 +144,7 @@ subroutine PressureForce_init(Time, G, GV, param_file, diag, CS, tides_CSp)
   else ; allocate(CS) ; endif
 
   ! Read all relevant parameters and write them to the model log.
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   call get_param(param_file, mod, "ANALYTIC_FV_PGF", CS%Analytic_FV_PGF, &
                  "If true the pressure gradient forces are calculated \n"//&
                  "with a finite volume form that analytically integrates \n"//&

src/core/MOM_PressureForce_Montgomery.F90

@@ -967,7 +967,7 @@ subroutine PressureForce_Mont_init(Time, G, GV, param_file, diag, CS, tides_CSp)
   endif
 
   mod = "MOM_PressureForce_Mont"
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   call get_param(param_file, mod, "RHO_0", CS%Rho0, &
                  "The mean ocean density used with BOUSSINESQ true to \n"//&
                  "calculate accelerations and the mass for conservation \n"//&

src/core/MOM_PressureForce_analytic_FV.F90

@@ -871,7 +871,7 @@ subroutine PressureForce_AFV_init(Time, G, GV, param_file, diag, CS, tides_CSp)
   endif
 
   mod = "MOM_PressureForce_AFV"
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   call get_param(param_file, mod, "RHO_0", CS%Rho0, &
                  "The mean ocean density used with BOUSSINESQ true to \n"//&
                  "calculate accelerations and the mass for conservation \n"//&

src/core/MOM_continuity.F90

@@ -190,7 +190,7 @@ subroutine continuity_init(Time, G, GV, param_file, diag, CS)
   allocate(CS)
 
   ! Read all relevant parameters and write them to the model log.
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   call get_param(param_file, mod, "CONTINUITY_SCHEME", tmpstr, &
                  "CONTINUITY_SCHEME selects the discretization for the \n"//&
                  "continuity solver. The only valid value currently is: \n"//&

src/core/MOM_continuity_PPM.F90

@@ -192,7 +192,7 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC,
 
   apply_OBC_u_flather_east = .false. ; apply_OBC_u_flather_west = .false.
   apply_OBC_v_flather_north = .false. ; apply_OBC_v_flather_south = .false.
-  if (present(OBC)) then ; if (associated(OBC)) then
+  if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%this_pe) then
     apply_OBC_u_flather_east = OBC%apply_OBC_u_flather_east
     apply_OBC_u_flather_west = OBC%apply_OBC_u_flather_west
     apply_OBC_v_flather_north = OBC%apply_OBC_v_flather_north
@@ -202,7 +202,7 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, CS, uhbt, vhbt, OBC,
     if (apply_OBC_u_flather_east .or. apply_OBC_u_flather_west .or. &
         apply_OBC_v_flather_north .or. apply_OBC_v_flather_south) &
       h_input(:,:,:) = hin(:,:,:)
-  endif ; endif
+  endif ; endif ; endif
 
   if (present(visc_rem_u) .neqv. present(visc_rem_v)) call MOM_error(FATAL, &
       "MOM_continuity_PPM: Either both visc_rem_u and visc_rem_v or neither"// &
@@ -424,9 +424,9 @@ subroutine zonal_mass_flux(u, h_in, uh, dt, G, GV, CS, LB, uhbt, OBC, &
   use_visc_rem = present(visc_rem_u)
   apply_OBC_u = .false. ; set_BT_cont = .false.
   if (present(BT_cont)) set_BT_cont = (associated(BT_cont))
-  if (present(OBC)) then ; if (associated(OBC)) then
+  if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%this_pe) then
     apply_OBC_u = OBC%apply_OBC_u
-  endif ; endif
+  endif ; endif ; endif
   ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = G%ke
 
   CFL_dt = CS%CFL_limit_adjust / dt
@@ -1181,9 +1181,9 @@ subroutine meridional_mass_flux(v, h_in, vh, dt, G, GV, CS, LB, vhbt, OBC, &
   use_visc_rem = present(visc_rem_v)
   apply_OBC_v = .false. ; set_BT_cont = .false.
   if (present(BT_cont)) set_BT_cont = (associated(BT_cont))
-  if (present(OBC)) then ; if (associated(OBC)) then
+  if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%this_pe) then
     apply_OBC_v = OBC%apply_OBC_v
-  endif ; endif
+  endif ; endif ; endif
   ish = LB%ish ; ieh = LB%ieh ; jsh = LB%jsh ; jeh = LB%jeh ; nz = G%ke
 
   CFL_dt = CS%CFL_limit_adjust / dt
@@ -2163,7 +2163,7 @@ subroutine continuity_PPM_init(Time, G, GV, param_file, diag, CS)
   allocate(CS)
 
 ! Read all relevant parameters and write them to the model log.
-  call log_version(param_file, mod, version)
+  call log_version(param_file, mod, version, "")
   call get_param(param_file, mod, "MONOTONIC_CONTINUITY", CS%monotonic, &
                  "If true, CONTINUITY_PPM uses the Colella and Woodward \n"//&
                  "monotonic limiter.  The default (false) is to use a \n"//&

src/core/MOM_open_boundary.F90

@@ -1,5 +1,5 @@
 ! This file is part of MOM6. See LICENSE.md for the license.
-!> Controls where open boundary conditions are applied 
+!> Controls where open boundary conditions are applied
 module MOM_open_boundary
 
 ! This file is part of MOM6. See LICENSE.md for the license.
@@ -11,6 +11,7 @@ module MOM_open_boundary
 use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING
 use MOM_file_parser, only : get_param, log_version, param_file_type, log_param
 use MOM_grid, only : ocean_grid_type
+use MOM_dyn_horgrid, only : dyn_horgrid_type
 use MOM_io, only : EAST_FACE, NORTH_FACE
 use MOM_io, only : slasher, read_data
 use MOM_tracer_registry, only : add_tracer_OBC_values, tracer_registry_type
@@ -100,6 +101,7 @@ module MOM_open_boundary
   real :: rx_max   !< The maximum magnitude of the baroclinic radiation
                    !! velocity (or speed of characteristics), in m s-1.  The
                    !! default value is 10 m s-1.
+  logical :: this_pe !< Is there an open boundary on this tile?
 end type ocean_OBC_type
 
 integer :: id_clock_pass
@@ -112,11 +114,9 @@ module MOM_open_boundary
 
 !> Enables OBC module and reads configuration parameters
 subroutine open_boundary_config(G, param_file, OBC)
-  type(ocean_grid_type),   intent(in)    :: G !< Ocean grid structure
+  type(dyn_horgrid_type),  intent(in)    :: G !< Ocean grid structure
   type(param_file_type),   intent(in)    :: param_file !< Parameter file handle
   type(ocean_OBC_type),    pointer       :: OBC !< Open boundary control structure
-  ! Local variables
-  logical :: flather_east, flather_west, flather_north, flather_south
 
   allocate(OBC)
 
@@ -164,9 +164,9 @@ end subroutine open_boundary_config
 
 !> Initialize open boundary control structure
 subroutine open_boundary_init(G, param_file, OBC)
-  type(ocean_grid_type),   intent(in)    :: G !< Ocean grid structure
-  type(param_file_type),   intent(in)    :: param_file !< Parameter file handle
-  type(ocean_OBC_type),    pointer       :: OBC !< Open boundary control structure
+  type(ocean_grid_type), intent(in)    :: G !< Ocean grid structure
+  type(param_file_type), intent(in)    :: param_file !< Parameter file handle
+  type(ocean_OBC_type),  pointer       :: OBC !< Open boundary control structure
   ! Local variables
 
   if (.not.associated(OBC)) return
@@ -244,7 +244,7 @@ end subroutine open_boundary_end
 !> Sets the slope of bathymetry normal to an open bounndary to zero.
 subroutine open_boundary_impose_normal_slope(OBC, G, depth)
   type(ocean_OBC_type),             pointer       :: OBC !< Open boundary control structure
-  type(ocean_grid_type),            intent(in)    :: G !< Ocean grid structure
+  type(dyn_horgrid_type),           intent(in)    :: G !< Ocean grid structure
   real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: depth !< Bathymetry at h-points
   ! Local variables
   integer :: i, j
@@ -270,7 +270,7 @@ end subroutine open_boundary_impose_normal_slope
 !> Reconcile masks and open boundaries, deallocate OBC on PEs where it is not needed
 subroutine open_boundary_impose_land_mask(OBC, G)
   type(ocean_OBC_type),              pointer       :: OBC !< Open boundary control structure
-  type(ocean_grid_type),             intent(in) :: G !< Ocean grid structure
+  type(dyn_horgrid_type),            intent(in) :: G !< Ocean grid structure
   ! Local variables
   integer :: i, j
   logical :: any_U, any_V
@@ -306,22 +306,24 @@ subroutine open_boundary_impose_land_mask(OBC, G)
       ! bathymetry inside the boundary was do shallow and flagged as land.
       if (OBC%OBC_mask_u(I,j)) any_U = .true.
     enddo ; enddo
-    if (.not. any_U) then
-      deallocate(OBC%OBC_mask_u)
-    endif
+!    if (.not. any_U) then
+!      deallocate(OBC%OBC_mask_u)
+!    endif
   endif
 
   any_V = .false.
   if (associated(OBC%OBC_mask_v)) then
     do J=G%JsdB,G%JedB ; do i=G%isd,G%ied
       if (OBC%OBC_mask_v(i,J)) any_V = .true.
     enddo ; enddo
-    if (.not. any_V) then
-      deallocate(OBC%OBC_mask_v)
-    endif
+!    if (.not. any_V) then
+!      deallocate(OBC%OBC_mask_v)
+!    endif
   endif
 
-  if (.not.(any_U .or. any_V)) call open_boundary_dealloc(OBC)
+!  if (.not.(any_U .or. any_V)) call open_boundary_dealloc(OBC)
+  OBC%this_pe = .true.
+  if (.not.(any_U .or. any_V)) OBC%this_pe = .false.
 
 end subroutine open_boundary_impose_land_mask
 
@@ -442,7 +444,7 @@ end subroutine Radiation_Open_Bdry_Conds
 !> Sets the domain boundaries as Flather open boundaries using the original
 !! Flather run-time logicals
 subroutine set_Flather_positions(G, OBC)
-  type(ocean_grid_type),                  intent(inout) :: G
+  type(dyn_horgrid_type),                 intent(inout) :: G
   type(ocean_OBC_type),                   pointer    :: OBC
   ! Local variables
   integer :: east_boundary, west_boundary, north_boundary, south_boundary
@@ -566,8 +568,6 @@ subroutine set_Flather_positions(G, OBC)
     enddo ; enddo
   endif
 
-  !   If there are no OBC points on this PE, there is no reason to keep the OBC
-  ! type, and it could be deallocated.
 end subroutine set_Flather_positions
 
 !> Sets the initial definitions of the characteristic open boundary conditions.
@@ -595,8 +595,8 @@ subroutine set_Flather_data(OBC, tv, h, G, PF, tracer_Reg)
   real, pointer, dimension(:,:,:) :: &
     OBC_T_u => NULL(), &    ! These arrays should be allocated and set to
     OBC_T_v => NULL(), &    ! specify the values of T and S that should come
-    OBC_S_u => NULL(), & 
-    OBC_S_v => NULL()     
+    OBC_S_u => NULL(), &
+    OBC_S_v => NULL()
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
@@ -639,7 +639,7 @@ subroutine set_Flather_data(OBC, tv, h, G, PF, tracer_Reg)
   if (.not.associated(OBC%eta_outer_v)) then
     allocate(OBC%eta_outer_v(isd:ied,JsdB:JedB)) ; OBC%eta_outer_v(:,:) = 0.0
   endif
-  
+
   if (read_OBC_uv) then
     call read_data(filename, 'ubt', OBC%ubt_outer, &
                    domain=G%Domain%mpp_domain, position=EAST_FACE)