Merge pull request mom-ocean#117 from Hallberg-NOAA/Bering_domain_rescaling_fix

+Fix some rescaling issues smoked out by the regional Bering domain

Author: Hallberg-NOAA

src/ALE/MOM_ALE.F90

@@ -718,6 +718,7 @@ subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n_itt, u, v, OBC, Reg, dt, d
   ! we have to keep track of the total dzInterface if for some reason
   ! we're using the old remapping algorithm for u/v
   real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzInterface, dzIntTotal
+  real :: h_neglect, h_neglect_edge ! small thicknesses [H ~> m or kg m-2]
 
   nz = GV%ke
 
@@ -743,7 +744,17 @@ subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n_itt, u, v, OBC, Reg, dt, d
   if (present(dt)) &
     call ALE_update_regrid_weights(dt, CS)
 
+  if (.not. CS%answers_2018) then
+    h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff
+  elseif (GV%Boussinesq) then
+    h_neglect = GV%m_to_H * 1.0e-30 ; h_neglect_edge = GV%m_to_H * 1.0e-10
+  else
+    h_neglect = GV%kg_m2_to_H * 1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H * 1.0e-10
+  endif
+
+
   do itt = 1, n_itt
+
     call do_group_pass(pass_T_S_h, G%domain)
 
     ! generate new grid
@@ -753,8 +764,10 @@ subroutine ALE_regrid_accelerated(CS, G, GV, h, tv, n_itt, u, v, OBC, Reg, dt, d
 
     ! remap from original grid onto new grid
     do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1
-      call remapping_core_h(CS%remapCS, nz, h_orig(i,j,:), tv%S(i,j,:), nz, h(i,j,:), tv_local%S(i,j,:))
-      call remapping_core_h(CS%remapCS, nz, h_orig(i,j,:), tv%T(i,j,:), nz, h(i,j,:), tv_local%T(i,j,:))
+      call remapping_core_h(CS%remapCS, nz, h_orig(i,j,:), tv%S(i,j,:), nz, h(i,j,:), tv_local%S(i,j,:), &
+                            h_neglect, h_neglect_edge)
+      call remapping_core_h(CS%remapCS, nz, h_orig(i,j,:), tv%T(i,j,:), nz, h(i,j,:), tv_local%T(i,j,:), &
+                            h_neglect, h_neglect_edge)
     enddo ; enddo
 
     ! starting grid for next iteration

src/core/MOM_barotropic.F90

@@ -17,7 +17,7 @@ module MOM_barotropic
 use MOM_grid, only : ocean_grid_type
 use MOM_hor_index, only : hor_index_type
 use MOM_io, only : vardesc, var_desc, MOM_read_data, slasher
-use MOM_open_boundary, only : ocean_OBC_type, OBC_SIMPLE, OBC_NONE, open_boundary_query
+use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE, open_boundary_query
 use MOM_open_boundary, only : OBC_DIRECTION_E, OBC_DIRECTION_W
 use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S, OBC_segment_type
 use MOM_restart, only : register_restart_field, register_restart_pair
@@ -1245,13 +1245,13 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
       if (Htot_avg*CS%dy_Cu(I,j) <= 0.0) then
         CS%IDatu(I,j) = 0.0
       elseif (integral_BT_cont) then
-        CS%IDatu(I,j) = CS%dy_Cu(I,j) / (max(find_duhbt_du(ubt(I,j)*dt, BTCL_u(I,j)), &
+        CS%IDatu(I,j) = GV%Z_to_H * CS%dy_Cu(I,j) / (max(find_duhbt_du(ubt(I,j)*dt, BTCL_u(I,j)), &
                                              CS%dy_Cu(I,j)*Htot_avg) )
       elseif (use_BT_cont) then ! Reconsider the max and whether there should be some scaling.
-        CS%IDatu(I,j) = CS%dy_Cu(I,j) / (max(find_duhbt_du(ubt(I,j), BTCL_u(I,j)), &
+        CS%IDatu(I,j) = GV%Z_to_H * CS%dy_Cu(I,j) / (max(find_duhbt_du(ubt(I,j), BTCL_u(I,j)), &
                                              CS%dy_Cu(I,j)*Htot_avg) )
       else
-        CS%IDatu(I,j) = 1.0 / Htot_avg
+        CS%IDatu(I,j) = GV%Z_to_H / Htot_avg
       endif
     endif
 
@@ -1271,13 +1271,13 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
       if (Htot_avg*CS%dx_Cv(i,J) <= 0.0) then
         CS%IDatv(i,J) = 0.0
       elseif (integral_BT_cont) then
-        CS%IDatv(i,J) = CS%dx_Cv(i,J) / (max(find_dvhbt_dv(vbt(i,J)*dt, BTCL_v(i,J)), &
+        CS%IDatv(i,J) = GV%Z_to_H * CS%dx_Cv(i,J) / (max(find_dvhbt_dv(vbt(i,J)*dt, BTCL_v(i,J)), &
                                              CS%dx_Cv(i,J)*Htot_avg) )
       elseif (use_BT_cont) then ! Reconsider the max and whether there should be some scaling.
-        CS%IDatv(i,J) = CS%dx_Cv(i,J) / (max(find_dvhbt_dv(vbt(i,J), BTCL_v(i,J)), &
+        CS%IDatv(i,J) = GV%Z_to_H * CS%dx_Cv(i,J) / (max(find_dvhbt_dv(vbt(i,J), BTCL_v(i,J)), &
                                              CS%dx_Cv(i,J)*Htot_avg) )
       else
-        CS%IDatv(i,J) = 1.0 / Htot_avg
+        CS%IDatv(i,J) = GV%Z_to_H / Htot_avg
       endif
     endif
 

src/core/MOM_open_boundary.F90

@@ -48,6 +48,9 @@ module MOM_open_boundary
 public open_boundary_test_extern_uv
 public open_boundary_test_extern_h
 public open_boundary_zero_normal_flow
+public parse_segment_str
+public parse_segment_manifest_str
+public parse_segment_data_str
 public register_OBC, OBC_registry_init
 public register_file_OBC, file_OBC_end
 public segment_tracer_registry_init
@@ -61,12 +64,10 @@ module MOM_open_boundary
 public rotate_OBC_config
 public rotate_OBC_init
 public initialize_segment_data
+public flood_fill
+public flood_fill2
 
 integer, parameter, public :: OBC_NONE = 0      !< Indicates the use of no open boundary
-integer, parameter, public :: OBC_SIMPLE = 1    !< Indicates the use of a simple inflow open boundary
-integer, parameter, public :: OBC_WALL = 2      !< Indicates the use of a closed wall
-integer, parameter, public :: OBC_FLATHER =  3  !< Indicates the use of a Flather open boundary
-integer, parameter, public :: OBC_RADIATION = 4 !< Indicates the use of a radiation open boundary
 integer, parameter, public :: OBC_DIRECTION_N = 100 !< Indicates the boundary is an effective northern boundary
 integer, parameter, public :: OBC_DIRECTION_S = 200 !< Indicates the boundary is an effective southern boundary
 integer, parameter, public :: OBC_DIRECTION_E = 300 !< Indicates the boundary is an effective eastern boundary
@@ -310,6 +311,9 @@ module MOM_open_boundary
   real :: ramp_value                        !< If ramp is True, where we are on the ramp from
                                             !! zero to one [nondim].
   type(time_type) :: ramp_start_time        !< Time when model was started.
+  logical :: answers_2018   !< If true, use the order of arithmetic and expressions for remapping
+                            !! that recover the answers from the end of 2018.  Otherwise, use more
+                            !! robust and accurate forms of mathematically equivalent expressions.
 end type ocean_OBC_type
 
 !> Control structure for open boundaries that read from files.
@@ -607,15 +611,15 @@ subroutine open_boundary_config(G, US, param_file, OBC)
     call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
                  "This sets the default value for the various _2018_ANSWERS parameters.", &
                  default=.false.)
-    call get_param(param_file, mdl, "REMAPPING_2018_ANSWERS", answers_2018, &
+    call get_param(param_file, mdl, "REMAPPING_2018_ANSWERS", OBC%answers_2018, &
                  "If true, use the order of arithmetic and expressions that recover the "//&
                  "answers from the end of 2018.  Otherwise, use updated and more robust "//&
                  "forms of the same expressions.", default=default_2018_answers)
 
     allocate(OBC%remap_CS)
     call initialize_remapping(OBC%remap_CS, remappingScheme, boundary_extrapolation = .false., &
                check_reconstruction=check_reconstruction, check_remapping=check_remapping, &
-               force_bounds_in_subcell=force_bounds_in_subcell, answers_2018=answers_2018)
+               force_bounds_in_subcell=force_bounds_in_subcell, answers_2018=OBC%answers_2018)
 
   endif ! OBC%number_of_segments > 0
 
@@ -1693,88 +1697,6 @@ subroutine parse_for_tracer_reservoirs(OBC, PF, use_temperature)
 
 end subroutine parse_for_tracer_reservoirs
 
-!> Parse an OBC_SEGMENT_%%%_PARAMS string
-subroutine parse_segment_param_real(segment_str, var, param_value, debug )
-  character(len=*),  intent(in)  :: segment_str !< A string in form of
-                                                !! "VAR1=file:foo1.nc(varnam1),VAR2=file:foo2.nc(varnam2),..."
-  character(len=*),  intent(in)  :: var         !< The name of the variable for which parameters are needed
-  real,              intent(out) :: param_value !< The value of the parameter
-  logical, optional, intent(in)  :: debug       !< If present and true, write verbose debugging messages
-  ! Local variables
-  character(len=128) :: word1, word2, word3, method
-  integer :: lword, nfields, n, m
-  logical :: continue,dbg
-  character(len=32), dimension(MAX_OBC_FIELDS) :: flds
-
-  nfields = 0
-  continue = .true.
-  dbg = .false.
-  if (PRESENT(debug)) dbg = debug
-
-  do while (continue)
-    word1 = extract_word(segment_str,',',nfields+1)
-    if (trim(word1) == '') exit
-    nfields = nfields+1
-    word2 = extract_word(word1,'=',1)
-    flds(nfields) = trim(word2)
-  enddo
-
-  ! if (PRESENT(fields)) then
-  !   do n=1,nfields
-  !     fields(n) = flds(n)
-  !   enddo
-  ! endif
-
-  ! if (PRESENT(num_fields)) then
-  !   num_fields = nfields
-  !   return
-  ! endif
-
-  m=0
-! if (PRESENT(var)) then
-    do n=1,nfields
-      if (trim(var)==trim(flds(n))) then
-        m = n
-        exit
-      endif
-    enddo
-    if (m==0) then
-      error stop
-    endif
-
-    ! Process first word which will start with the fieldname
-    word3 = extract_word(segment_str,',',m)
-!     word1 = extract_word(word3,':',1)
-!     if (trim(word1) == '') exit
-    word2 = extract_word(word1,'=',1)
-    if (trim(word2) == trim(var)) then
-      method=trim(extract_word(word1,'=',2))
-      lword=len_trim(method)
-      read(method(1:lword),*,err=987) param_value
-      ! if (method(lword-3:lword) == 'file') then
-      !   ! raise an error id filename/fieldname not in argument list
-      !   word1 = extract_word(word3,':',2)
-      !   filenam = extract_word(word1,'(',1)
-      !   fieldnam = extract_word(word1,'(',2)
-      !   lword=len_trim(fieldnam)
-      !   fieldnam = fieldnam(1:lword-1)  ! remove trailing parenth
-      !   value=-999.
-      ! elseif (method(lword-4:lword) == 'value') then
-      !   filenam = 'none'
-      !   fieldnam = 'none'
-      !   word1 = extract_word(word3,':',2)
-      !   lword=len_trim(word1)
-      !   read(word1(1:lword),*,end=986,err=987) value
-      ! endif
-    endif
-! endif
-
-  return
-  986 call MOM_error(FATAL,'End of record while parsing segment data specification! '//trim(segment_str))
-  987 call MOM_error(FATAL,'Error while parsing segment parameter specification! '//trim(segment_str))
-
-end subroutine parse_segment_param_real
-
 !> Initialize open boundary control structure and do any necessary rescaling of OBC
 !! fields that have been read from a restart file.
 subroutine open_boundary_init(G, GV, US, param_file, OBC, restart_CS)
@@ -3722,6 +3644,7 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
   real, allocatable :: normal_trans_bt(:,:) ! barotropic transport [H L2 T-1 ~> m3 s-1]
   integer :: turns    ! Number of index quarter turns
   real :: time_delta  ! Time since tidal reference date [T ~> s]
+  real :: h_neglect, h_neglect_edge ! Small thicknesses [H ~> m or kg m-2]
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
@@ -3734,6 +3657,14 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
 
   if (OBC%add_tide_constituents) time_delta = US%s_to_T * time_type_to_real(Time - OBC%time_ref)
 
+  if (.not. OBC%answers_2018) then
+    h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff
+  elseif (GV%Boussinesq) then
+    h_neglect = GV%m_to_H * 1.0e-30 ; h_neglect_edge = GV%m_to_H * 1.0e-10
+  else
+    h_neglect = GV%kg_m2_to_H * 1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H * 1.0e-10
+  endif
+
   do n = 1, OBC%number_of_segments
     segment => OBC%segment(n)
 
@@ -3932,7 +3863,8 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
         ! no dz for tidal variables
         if (segment%field(m)%nk_src > 1 .and.&
             (index(segment%field(m)%name, 'phase') <= 0 .and. index(segment%field(m)%name, 'amp') <= 0)) then
-          call time_interp_external(segment%field(m)%fid_dz,Time, tmp_buffer_in)
+          call time_interp_external(segment%field(m)%fid_dz, Time, tmp_buffer_in)
+          tmp_buffer_in(:,:,:) = tmp_buffer_in(:,:,:) * US%m_to_Z
           if (turns /= 0) then
             ! TODO: This is hardcoded for 90 degrees, and needs to be generalized.
             if (segment%is_E_or_W &
@@ -3998,19 +3930,22 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
                   call remapping_core_h(OBC%remap_CS, &
                        segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), &
                        segment%field(m)%buffer_src(I,J,:), &
-                       GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:))
+                       GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), &
+                       h_neglect, h_neglect_edge)
                 elseif (G%mask2dCu(I,j)>0.) then
                   h_stack(:) = h(i+ishift,j,:)
                   call remapping_core_h(OBC%remap_CS, &
                        segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), &
                        segment%field(m)%buffer_src(I,J,:), &
-                       GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:))
+                       GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), &
+                       h_neglect, h_neglect_edge)
                 elseif (G%mask2dCu(I,j+1)>0.) then
                   h_stack(:) = h(i+ishift,j+1,:)
                   call remapping_core_h(OBC%remap_CS, &
                        segment%field(m)%nk_src,segment%field(m)%dz_src(I,j,:), &
                        segment%field(m)%buffer_src(I,J,:), &
-                       GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:))
+                       GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), &
+                       h_neglect, h_neglect_edge)
                 endif
               enddo
             else
@@ -4025,7 +3960,8 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
                   call remapping_core_h(OBC%remap_CS, &
                        segment%field(m)%nk_src, scl_fac*segment%field(m)%dz_src(I,j,:), &
                        segment%field(m)%buffer_src(I,j,:), &
-                       GV%ke, h(i+ishift,j,:), segment%field(m)%buffer_dst(I,j,:))
+                       GV%ke, h(i+ishift,j,:), segment%field(m)%buffer_dst(I,j,:), &
+                       h_neglect, h_neglect_edge)
                 endif
               enddo
             endif
@@ -4044,19 +3980,22 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
                   call remapping_core_h(OBC%remap_CS, &
                        segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), &
                        segment%field(m)%buffer_src(I,J,:), &
-                       GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:))
+                       GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), &
+                       h_neglect, h_neglect_edge)
                 elseif (G%mask2dCv(i,J)>0.) then
                   h_stack(:) = h(i,j+jshift,:)
                   call remapping_core_h(OBC%remap_CS, &
                        segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), &
                        segment%field(m)%buffer_src(I,J,:), &
-                       GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:))
+                       GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), &
+                       h_neglect, h_neglect_edge)
                 elseif (G%mask2dCv(i+1,J)>0.) then
                   h_stack(:) = h(i+1,j+jshift,:)
                   call remapping_core_h(OBC%remap_CS, &
                        segment%field(m)%nk_src,segment%field(m)%dz_src(I,J,:), &
                        segment%field(m)%buffer_src(I,J,:), &
-                       GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:))
+                       GV%ke, h_stack, segment%field(m)%buffer_dst(I,J,:), &
+                       h_neglect, h_neglect_edge)
                 endif
               enddo
             else
@@ -4071,7 +4010,8 @@ subroutine update_OBC_segment_data(G, GV, US, OBC, tv, h, Time)
                   call remapping_core_h(OBC%remap_CS, &
                        segment%field(m)%nk_src, scl_fac*segment%field(m)%dz_src(i,J,:), &
                        segment%field(m)%buffer_src(i,J,:), &
-                       GV%ke, h(i,j+jshift,:), segment%field(m)%buffer_dst(i,J,:))
+                       GV%ke, h(i,j+jshift,:), segment%field(m)%buffer_dst(i,J,:), &
+                       h_neglect, h_neglect_edge)
                 endif
               enddo
             endif

src/framework/MOM_horizontal_regridding.F90

@@ -856,7 +856,7 @@ subroutine horiz_interp_and_extrap_tracer_fms_id(fms_id,  Time, conversion, G, t
       do k=1,kd
         do j=js,je
           do i=is,ie
-            tr_z(i,j,k)=data_in(i,j,k)
+            tr_z(i,j,k)=data_in(i,j,k) * conversion
             if (.not. ans_2018) mask_z(i,j,k) = 1.
             if (abs(tr_z(i,j,k)-missing_value) < abs(roundoff*missing_value)) mask_z(i,j,k) = 0.
           enddo

src/initialization/MOM_state_initialization.F90

@@ -20,7 +20,7 @@ module MOM_state_initialization
 use MOM_interface_heights, only : find_eta
 use MOM_io, only : file_exists, field_size, MOM_read_data, MOM_read_vector, slasher
 use MOM_open_boundary, only : ocean_OBC_type, open_boundary_init, set_tracer_data
-use MOM_open_boundary, only : OBC_NONE, OBC_SIMPLE
+use MOM_open_boundary, only : OBC_NONE
 use MOM_open_boundary, only : open_boundary_query
 use MOM_open_boundary, only : set_tracer_data, initialize_segment_data
 use MOM_open_boundary, only : open_boundary_test_extern_h