Merge pull request #3 from NOAA-GFDL/dev/gfdl

Latest MOM6 dev/gfdl updates

config_src/coupled_driver/MOM_surface_forcing.F90

@@ -1413,9 +1413,8 @@ subroutine surface_forcing_init(Time, G, param_file, diag, CS)
   call get_param(param_file, mdl, "ALLOW_FLUX_ADJUSTMENTS", CS%allow_flux_adjustments, &
                  "If true, allows flux adjustments to specified via the \n"//&
                  "data_table using the component name 'OCN'.", default=.false.)
-  if (CS%allow_flux_adjustments) then
-    call data_override_init(Ocean_domain_in=G%Domain%mpp_domain)
-  endif
+
+  call data_override_init(Ocean_domain_in=G%Domain%mpp_domain)
 
   if (CS%restore_salt) then
     salt_file = trim(CS%inputdir) // trim(CS%salt_restore_file)

src/core/MOM_forcing_type.F90

@@ -2133,7 +2133,7 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, diag, handles)
         if (associated(fluxes%frunoff))     res(i,j) = res(i,j)+fluxes%frunoff(i,j)
         if (associated(fluxes%vprec))       res(i,j) = res(i,j)+fluxes%vprec(i,j)
       enddo ; enddo
-      call post_data(handles%id_prcme, res, diag)
+      if (handles%id_prcme > 0) call post_data(handles%id_prcme, res, diag)
       if (handles%id_total_prcme > 0) then
         total_transport = global_area_integral(res,G)
         call post_data(handles%id_total_prcme, total_transport, diag)
@@ -2151,7 +2151,7 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, diag, handles)
         if (fluxes%vprec(i,j) < 0.0) res(i,j) = res(i,j) + fluxes%vprec(i,j)
         if (fluxes%evap(i,j)  < 0.0) res(i,j) = res(i,j) + fluxes%evap(i,j)
       enddo ; enddo
-      call post_data(handles%id_net_massout, res, diag)
+      if (handles%id_net_massout > 0) call post_data(handles%id_net_massout, res, diag)
       if (handles%id_total_net_massout > 0) then
         total_transport = global_area_integral(res,G)
         call post_data(handles%id_total_net_massout, total_transport, diag)
@@ -2168,7 +2168,7 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, diag, handles)
         ! fluxes%cond is not needed because it is derived from %evap > 0
         if (fluxes%evap(i,j)  > 0.0) res(i,j) = res(i,j) + fluxes%evap(i,j)
       enddo ; enddo
-      call post_data(handles%id_net_massin, res, diag)
+      if (handles%id_net_massin > 0) call post_data(handles%id_net_massin, res, diag)
       if (handles%id_total_net_massin > 0) then
         total_transport = global_area_integral(res,G)
         call post_data(handles%id_total_net_massin, total_transport, diag)
@@ -2322,7 +2322,7 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, diag, handles)
       if (associated(fluxes%sens))       res(i,j) = res(i,j) + fluxes%sens(i,j)
       if (associated(fluxes%SW))         res(i,j) = res(i,j) + fluxes%SW(i,j)
       enddo ; enddo
-      call post_data(handles%id_net_heat_coupler, res, diag)
+      if (handles%id_net_heat_coupler > 0) call post_data(handles%id_net_heat_coupler, res, diag)
       if (handles%id_total_net_heat_coupler > 0) then
         total_transport = global_area_integral(res,G)
         call post_data(handles%id_total_net_heat_coupler, total_transport, diag)
@@ -2382,7 +2382,7 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, diag, handles)
           if (associated(fluxes%heat_content_massout)) res(i,j) = res(i,j) + fluxes%heat_content_massout(i,j)
       ! endif
       enddo ; enddo
-      call post_data(handles%id_heat_content_surfwater, res, diag)
+      if (handles%id_heat_content_surfwater > 0) call post_data(handles%id_heat_content_surfwater, res, diag)
       if (handles%id_total_heat_content_surfwater > 0) then
         total_transport = global_area_integral(res,G)
         call post_data(handles%id_total_heat_content_surfwater, total_transport, diag)

src/core/MOM_transcribe_grid.F90

@@ -143,7 +143,7 @@ subroutine copy_dyngrid_to_MOM_grid(dG, oG)
   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)
+  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)
@@ -287,7 +287,7 @@ subroutine copy_MOM_grid_to_dyngrid(oG, dG)
   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)
+  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)

src/framework/MOM_domains.F90

@@ -27,18 +27,18 @@ module MOM_domains
 use mpp_domains_mod, only : compute_block_extent => mpp_compute_block_extent
 use mpp_parameter_mod, only : AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR, BITWISE_EXACT_SUM, CORNER
 use mpp_parameter_mod, only : To_East => WUPDATE, To_West => EUPDATE, Omit_Corners => EDGEUPDATE
-use mpp_parameter_mod, only : To_North => SUPDATE, To_South => NUPDATE
+use mpp_parameter_mod, only : To_North => SUPDATE, To_South => NUPDATE, CENTER
 use fms_io_mod,        only : file_exist, parse_mask_table
 
 implicit none ; private
 
 public :: MOM_domains_init, MOM_infra_init, MOM_infra_end, get_domain_extent
 public :: MOM_define_domain, MOM_define_io_domain, clone_MOM_domain
-public :: pass_var, pass_vector, broadcast, PE_here, root_PE, num_PEs
-public :: pass_var_start, pass_var_complete, fill_symmetric_edges
+public :: pass_var, pass_vector, PE_here, root_PE, num_PEs
+public :: pass_var_start, pass_var_complete, fill_symmetric_edges, broadcast
 public :: pass_vector_start, pass_vector_complete
 public :: global_field_sum, sum_across_PEs, min_across_PEs, max_across_PEs
-public :: AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR, BITWISE_EXACT_SUM, CORNER
+public :: AGRID, BGRID_NE, CGRID_NE, SCALAR_PAIR, BITWISE_EXACT_SUM, CORNER, CENTER
 public :: To_East, To_West, To_North, To_South, To_All, Omit_Corners
 public :: create_group_pass, do_group_pass, group_pass_type
 public :: start_group_pass, complete_group_pass
@@ -178,8 +178,7 @@ subroutine pass_var_3d(array, MOM_dom, sideflag, complete, position, halo, &
 end subroutine pass_var_3d
 
 !> pass_var_2d does a halo update for a two-dimensional array.
-subroutine pass_var_2d(array, MOM_dom, sideflag, complete, position, halo, &
-                       clock)
+subroutine pass_var_2d(array, MOM_dom, sideflag, complete, position, halo, inner_halo, clock)
   real, dimension(:,:),  intent(inout) :: array    !< The array which is having its halos points
                                                    !! exchanged.
   type(MOM_domain_type), intent(inout) :: MOM_dom  !< The MOM_domain_type containing the mpp_domain
@@ -197,18 +196,34 @@ subroutine pass_var_2d(array, MOM_dom, sideflag, complete, position, halo, &
                                                    !! by default.
   integer,     optional, intent(in)    :: halo     !< The size of the halo to update - the full halo
                                                    !! by default.
-  integer,      optional, intent(in)   :: clock    !< The handle for a cpu time clock that should be
+  integer,     optional, intent(in)    :: inner_halo !< The size of an inner halo to avoid updating,
+                                                   !! or 0 to avoid updating symmetric memory
+                                                   !! computational domain points.  Setting this >=0
+                                                   !! also enforces that complete=.true.
+  integer,     optional, intent(in)    :: clock    !< The handle for a cpu time clock that should be
                                                    !! started then stopped to time this routine.
 
+  ! Local variables
+  real, allocatable, dimension(:,:) :: tmp
+  integer :: pos, i_halo, j_halo
+  integer :: isc, iec, jsc, jec, isd, ied, jsd, jed, IscB, IecB, JscB, JecB
+  integer :: inner, i, j, isfw, iefw, isfe, iefe, jsfs, jefs, jsfn, jefn
   integer :: dirflag
   logical :: block_til_complete
 
   if (present(clock)) then ; if (clock>0) call cpu_clock_begin(clock) ; endif
 
   dirflag = To_All ! 60
   if (present(sideflag)) then ; if (sideflag > 0) dirflag = sideflag ; endif
-  block_til_complete = .true.
-  if (present(complete)) block_til_complete = complete
+  block_til_complete = .true. ; if (present(complete)) block_til_complete = complete
+  pos = CENTER ; if (present(position)) pos = position
+
+  if (present(inner_halo)) then ; if (inner_halo >= 0) then
+    ! Store the original values.
+    allocate(tmp(size(array,1), size(array,2)))
+    tmp(:,:) = array(:,:)
+    block_til_complete = .true.
+  endif ; endif
 
   if (present(halo) .and. MOM_dom%thin_halo_updates) then
     call mpp_update_domains(array, MOM_dom%mpp_domain, flags=dirflag, &
@@ -219,6 +234,46 @@ subroutine pass_var_2d(array, MOM_dom, sideflag, complete, position, halo, &
                         complete=block_til_complete, position=position)
   endif
 
+  if (present(inner_halo)) then ; if (inner_halo >= 0) then
+    call mpp_get_compute_domain(MOM_dom%mpp_domain, isc, iec, jsc, jec)
+    call mpp_get_data_domain(MOM_dom%mpp_domain, isd, ied, jsd, jed)
+    ! Convert to local indices for arrays starting at 1.
+    isc = isc - (isd-1) ; iec = iec - (isd-1) ; ied = ied - (isd-1) ; isd = 1
+    jsc = jsc - (jsd-1) ; jec = jec - (jsd-1) ; jed = jed - (jsd-1) ; jsd = 1
+    i_halo = min(inner_halo, isc-1) ; j_halo = min(inner_halo, jsc-1)
+
+    ! Figure out the array index extents of the eastern, western, northern and southern regions to copy.
+    if (pos == CENTER) then
+      if (size(array,1) == ied) then
+        isfw = isc - i_halo ; iefw = isc ; isfe = iec ; iefe = iec + i_halo
+      else ; call MOM_error(FATAL, "pass_var_2d: wrong i-size for CENTER array.") ; endif
+      if (size(array,2) == jed) then
+        isfw = isc - i_halo ; iefw = isc ; isfe = iec ; iefe = iec + i_halo
+      else ; call MOM_error(FATAL, "pass_var_2d: wrong j-size for CENTER array.") ; endif
+    elseif (pos == CORNER) then
+      if (size(array,1) == ied) then
+        isfw = max(isc - (i_halo+1), 1) ; iefw = isc ; isfe = iec ; iefe = iec + i_halo
+      elseif (size(array,1) == ied+1) then
+        isfw = isc - i_halo ; iefw = isc+1 ; isfe = iec+1 ; iefe = min(iec + 1 + i_halo, ied+1)
+      else ; call MOM_error(FATAL, "pass_var_2d: wrong i-size for CORNER array.") ; endif
+      if (size(array,2) == jed) then
+        jsfs = max(jsc - (j_halo+1), 1) ; jefs = jsc ; jsfn = jec ; jefn = jec + j_halo
+      elseif (size(array,2) == jed+1) then
+        jsfs = jsc - j_halo ; jefs = jsc+1 ; jsfn = jec+1 ; jefn = min(jec + 1 + j_halo, jed+1)
+      else ; call MOM_error(FATAL, "pass_var_2d: wrong j-size for CORNER array.") ; endif
+    else
+      call MOM_error(FATAL, "pass_var_2d: Unrecognized position")
+    endif
+
+    ! Copy back the stored inner halo points
+    do j=jsfs,jefn ; do i=isfw,iefw ; array(i,j) = tmp(i,j) ; enddo ; enddo
+    do j=jsfs,jefn ; do i=isfe,iefe ; array(i,j) = tmp(i,j) ; enddo ; enddo
+    do j=jsfs,jefs ; do i=isfw,iefe ; array(i,j) = tmp(i,j) ; enddo ; enddo
+    do j=jsfn,jefn ; do i=isfw,iefe ; array(i,j) = tmp(i,j) ; enddo ; enddo
+
+    deallocate(tmp)
+  endif ; endif
+
   if (present(clock)) then ; if (clock>0) call cpu_clock_end(clock) ; endif
 
 end subroutine pass_var_2d

src/initialization/MOM_grid_initialize.F90

@@ -183,6 +183,7 @@ subroutine set_grid_metrics_from_mosaic(G, param_file)
   character(len=64)  :: mdl = "MOM_grid_init set_grid_metrics_from_mosaic"
   integer :: err=0, ni, nj, global_indices(4)
   type(MOM_domain_type) :: SGdom ! Supergrid domain
+  logical :: lon_bug  ! If true use an older buggy answer in the tripolar longitude.
   integer :: i, j, i2, j2
   integer :: npei,npej
   integer, dimension(:), allocatable :: exni,exnj
@@ -193,6 +194,10 @@ subroutine set_grid_metrics_from_mosaic(G, param_file)
   call get_param(param_file, mdl, "GRID_FILE", grid_file, &
                  "Name of the file from which to read horizontal grid data.", &
                  fail_if_missing=.true.)
+  call get_param(param_file, mdl, "USE_TRIPOLAR_GEOLONB_BUG", lon_bug, &
+                 "If true, use older code that incorrectly sets the longitude \n"//&
+                 "in some points along the tripolar fold to be off by 360 degrees.", &
+                 default=.true.)
   call get_param(param_file,  mdl, "INPUTDIR", inputdir, default=".")
   inputdir = slasher(inputdir)
   filename = trim(adjustl(inputdir)) // trim(adjustl(grid_file))
@@ -248,7 +253,11 @@ subroutine set_grid_metrics_from_mosaic(G, param_file)
   tmpZ(:,:) = 999.
   call MOM_read_data(filename, 'x', tmpZ, SGdom, position=CORNER)
 
-  call pass_var(tmpZ, SGdom, position=CORNER)
+  if (lon_bug) then
+    call pass_var(tmpZ, SGdom, position=CORNER)
+  else
+    call pass_var(tmpZ, SGdom, position=CORNER, inner_halo=0)
+  endif
   call extrapolate_metric(tmpZ, 2*(G%jsc-G%jsd)+2, missing=999.)
   do j=G%jsd,G%jed ; do i=G%isd,G%ied ; i2 = 2*i ; j2 = 2*j
     G%geoLonT(i,j) = tmpZ(i2-1,j2-1)