Rotation bugfixes

This patch contains several bugfixes associated with the rotational grid
testing.

The following checksums are declared as scalar pairs, rather than
vectors:

* eta_[uv] (in porous barrier)
* por_face_area[UV]
* por_layer_width[UV]
* Ray_[uv] (Rayleigh drag velocity)

Fluxes and surface fields are now permitted to contain tracer fluxes
(tr_fluxes) when rotation is enabled.  The fields are retained in their
unrotated form, since these are accessed and handled outside of MOM6.

The rotated p_surf_SSH pointer in forces now correctly points to either
p_surf or p_surf_full.

read_netcdf_nc() now correctly uses the unrotated horizontal index
struct HI, used to access the contents of the file.  Previously, it was
using the model HI, which may be rotated.

Reading chlorophyll with time_interp_external now uses rotation to
correctly fetch its output.

NOTE: This could be cleaned up so that the rotation details are hidden
from users, but there are some unresolved issues around how to approach
this.

src/core/MOM_forcing_type.F90

@@ -6,6 +6,7 @@ module MOM_forcing_type
 use MOM_array_transform, only : rotate_array, rotate_vector, rotate_array_pair
 use MOM_coupler_types, only : coupler_2d_bc_type, coupler_type_destructor
 use MOM_coupler_types, only : coupler_type_increment_data, coupler_type_initialized
+use MOM_coupler_types, only : coupler_type_copy_data
 use MOM_cpu_clock,     only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE
 use MOM_debugging,     only : hchksum, uvchksum
 use MOM_diag_mediator, only : post_data, register_diag_field, register_scalar_field
@@ -3702,9 +3703,11 @@ subroutine rotate_forcing(fluxes_in, fluxes, turns)
   if (associated(fluxes_in%ustar_tidal)) &
     call rotate_array(fluxes_in%ustar_tidal, turns, fluxes%ustar_tidal)
 
-  ! TODO: tracer flux rotation
-  if (coupler_type_initialized(fluxes%tr_fluxes)) &
-    call MOM_error(FATAL, "Rotation of tracer BC fluxes not yet implemented.")
+  ! NOTE: Tracer fields are handled by FMS, so are left unrotated.  Any
+  ! reads/writes to tr_fields must be appropriately rotated.
+  if (coupler_type_initialized(fluxes%tr_fluxes)) then
+    call coupler_type_copy_data(fluxes_in%tr_fluxes, fluxes%tr_fluxes)
+  endif
 
   ! Scalars and flags
   fluxes%accumulate_p_surf = fluxes_in%accumulate_p_surf
@@ -3757,10 +3760,18 @@ subroutine rotate_mech_forcing(forces_in, turns, forces)
   endif
 
   if (do_press) then
-    ! NOTE: p_surf_SSH either points to p_surf or p_surf_full
     call rotate_array(forces_in%p_surf, turns, forces%p_surf)
     call rotate_array(forces_in%p_surf_full, turns, forces%p_surf_full)
     call rotate_array(forces_in%net_mass_src, turns, forces%net_mass_src)
+
+    ! p_surf_SSH points to either p_surf or p_surf_full
+    if (associated(forces_in%p_surf_SSH, forces_in%p_surf)) then
+      forces%p_surf_SSH => forces%p_surf
+    else if (associated(forces_in%p_surf_SSH, forces_in%p_surf_full)) then
+      forces%p_surf_SSH => forces%p_surf_full
+    else
+      forces%p_surf_SSH => null()
+    endif
   endif
 
   if (do_iceberg) then

src/core/MOM_porous_barriers.F90

@@ -168,9 +168,10 @@ subroutine porous_widths_layer(h, tv, G, GV, US, pbv, CS, eta_bt)
 
   if (CS%debug) then
     call uvchksum("Interface height used by porous barrier for layer weights", &
-                  eta_u, eta_v, G%HI, haloshift=0)
+                  eta_u, eta_v, G%HI, haloshift=0, scalar_pair=.true.)
     call uvchksum("Porous barrier layer-averaged weights: por_face_area[UV]", &
-                  pbv%por_face_areaU, pbv%por_face_areaV, G%HI, haloshift=0)
+                  pbv%por_face_areaU, pbv%por_face_areaV, G%HI, haloshift=0, &
+                  scalar_pair=.true.)
   endif
 
   if (CS%id_por_face_areaU > 0) call post_data(CS%id_por_face_areaU, pbv%por_face_areaU, CS%diag)
@@ -256,9 +257,10 @@ subroutine porous_widths_interface(h, tv, G, GV, US, pbv, CS, eta_bt)
 
   if (CS%debug) then
     call uvchksum("Interface height used by porous barrier for interface weights", &
-                  eta_u, eta_v, G%HI, haloshift=0)
+                  eta_u, eta_v, G%HI, haloshift=0, scalar_pair=.true.)
     call uvchksum("Porous barrier weights at the layer-interface: por_layer_width[UV]", &
-                  pbv%por_layer_widthU, pbv%por_layer_widthV, G%HI, haloshift=0)
+                  pbv%por_layer_widthU, pbv%por_layer_widthV, G%HI, &
+                  haloshift=0, scalar_pair=.true.)
   endif
 
   if (CS%id_por_layer_widthU > 0) call post_data(CS%id_por_layer_widthU, pbv%por_layer_widthU, CS%diag)

src/core/MOM_variables.F90

@@ -6,6 +6,7 @@ module MOM_variables
 use MOM_array_transform, only : rotate_array, rotate_vector
 use MOM_coupler_types, only : coupler_1d_bc_type, coupler_2d_bc_type
 use MOM_coupler_types, only : coupler_type_spawn, coupler_type_destructor, coupler_type_initialized
+use MOM_coupler_types, only : coupler_type_copy_data
 use MOM_debugging,     only : hchksum
 use MOM_domains,       only : MOM_domain_type, get_domain_extent, group_pass_type
 use MOM_EOS,           only : EOS_type
@@ -499,9 +500,11 @@ subroutine rotate_surface_state(sfc_state_in, sfc_state, G, turns)
   sfc_state%T_is_conT = sfc_state_in%T_is_conT
   sfc_state%S_is_absS = sfc_state_in%S_is_absS
 
-  ! TODO: tracer field rotation
-  if (coupler_type_initialized(sfc_state_in%tr_fields)) &
-    call MOM_error(FATAL, "Rotation of surface state tracers is not yet implemented.")
+  ! NOTE: Tracer fields are handled by FMS, so are left unrotated.  Any
+  ! reads/writes to tr_fields must be appropriately rotated.
+  if (coupler_type_initialized(sfc_state_in%tr_fields)) then
+    call coupler_type_copy_data(sfc_state_in%tr_fields, sfc_state%tr_fields)
+  endif
 end subroutine rotate_surface_state
 
 !> Allocates the arrays contained within a BT_cont_type and initializes them to 0.

src/framework/MOM_io_file.F90

@@ -1326,7 +1326,13 @@ subroutine open_file_nc(handle, filename, action, MOM_domain, threading, fileset
 
   if (present(MOM_domain)) then
     handle%domain_decomposed = .true.
-    call hor_index_init(MOM_domain, handle%HI)
+
+    ! Input files use unrotated indexing.
+    if (associated(MOM_domain%domain_in)) then
+      call hor_index_init(MOM_domain%domain_in, handle%HI)
+    else
+      call hor_index_init(MOM_domain, handle%HI)
+    endif
   endif
 
   call handle%axes%init()

src/parameterizations/vertical/MOM_diabatic_aux.F90

@@ -643,7 +643,7 @@ subroutine set_pen_shortwave(optics, fluxes, G, GV, US, CS, opacity, tracer_flow
     if (CS%chl_from_file) then
       ! Only the 2-d surface chlorophyll can be read in from a file.  The
       ! same value is assumed for all layers.
-      call time_interp_external(CS%sbc_chl, CS%Time, chl_2d)
+      call time_interp_external(CS%sbc_chl, CS%Time, chl_2d, turns=G%HI%turns)
       do j=js,je ; do i=is,ie
         if ((G%mask2dT(i,j) > 0.0) .and. (chl_2d(i,j) < 0.0)) then
           write(mesg,'(" Time_interp negative chl of ",(1pe12.4)," at i,j = ",&
@@ -1899,7 +1899,11 @@ subroutine diabatic_aux_init(Time, G, GV, US, param_file, diag, CS, useALEalgori
         call log_param(param_file, mdl, "INPUTDIR/CHL_FILE", chl_filename)
         call get_param(param_file, mdl, "CHL_VARNAME", chl_varname, &
                    "Name of CHL_A variable in CHL_FILE.", default='CHL_A')
-        CS%sbc_chl = init_external_field(chl_filename, trim(chl_varname), MOM_domain=G%Domain)
+        if (modulo(G%Domain%turns, 4) /= 0) then
+          CS%sbc_chl = init_external_field(chl_filename, trim(chl_varname), MOM_domain=G%Domain%domain_in)
+        else
+          CS%sbc_chl = init_external_field(chl_filename, trim(chl_varname), MOM_domain=G%Domain)
+        endif
       endif
 
       CS%id_chl = register_diag_field('ocean_model', 'Chl_opac', diag%axesT1, Time, &

src/parameterizations/vertical/MOM_set_diffusivity.F90

@@ -613,7 +613,8 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i
     endif
 
     if (allocated(visc%Ray_u) .and. allocated(visc%Ray_v)) then
-      call uvchksum("Ray_[uv]", visc%Ray_u, visc%Ray_v, G%HI, 0, symmetric=.true., scale=GV%H_to_m*US%s_to_T)
+      call uvchksum("Ray_[uv]", visc%Ray_u, visc%Ray_v, G%HI, 0, &
+          symmetric=.true., scale=GV%H_to_m*US%s_to_T, scalar_pair=.true.)
     endif
 
   endif

src/parameterizations/vertical/MOM_set_viscosity.F90

@@ -1072,7 +1072,8 @@ subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
 
   if (CS%debug) then
     if (allocated(visc%Ray_u) .and. allocated(visc%Ray_v)) &
-        call uvchksum("Ray [uv]", visc%Ray_u, visc%Ray_v, G%HI, haloshift=0, scale=GV%H_to_m*US%s_to_T)
+        call uvchksum("Ray [uv]", visc%Ray_u, visc%Ray_v, G%HI, haloshift=0, &
+                      scale=GV%H_to_m*US%s_to_T, scalar_pair=.true.)
     if (allocated(visc%kv_bbl_u) .and. allocated(visc%kv_bbl_v)) &
         call uvchksum("kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, G%HI, &
                       haloshift=0, scale=GV%HZ_T_to_m2_s, scalar_pair=.true.)