Merge branch 'dev/ncar' into ncar_GME_willy_merged

Conflicts:
	src/parameterizations/lateral/MOM_hor_visc.F90
	src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90

config_src/mct_driver/MOM_ocean_model.F90

@@ -559,8 +559,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
   OS%nstep = OS%nstep + 1
 
   call enable_averaging(time_step, OS%Time, OS%diag)
-  call mech_forcing_diags(OS%forces, OS%fluxes, time_step, OS%grid, &
-                          OS%diag, OS%forcing_CSp%handles)
+  call mech_forcing_diags(OS%forces, time_step, OS%grid, OS%diag, OS%forcing_CSp%handles)
   call disable_averaging(OS%diag)
 
   if (OS%fluxes%fluxes_used) then

config_src/mct_driver/MOM_surface_forcing.F90

@@ -152,11 +152,13 @@ module MOM_surface_forcing
 ! MOM-based coupled models.
 type, public :: ice_ocean_boundary_type
   real, pointer, dimension(:,:) :: latent_flux     =>NULL() !< latent flux (W/m2)
-  real, pointer, dimension(:,:) :: rofl_flux       =>NULL() !< liquid runoff (W/m2)
-  real, pointer, dimension(:,:) :: rofi_flux       =>NULL() !< ice runoff (W/m2)
+  real, pointer, dimension(:,:) :: rofl_flux       =>NULL() !< liquid runoff (kg/m2/s)
+  real, pointer, dimension(:,:) :: rofi_flux       =>NULL() !< ice runoff (kg/m2/s)
   real, pointer, dimension(:,:) :: u_flux          =>NULL() !< i-direction wind stress (Pa)
   real, pointer, dimension(:,:) :: v_flux          =>NULL() !< j-direction wind stress (Pa)
   real, pointer, dimension(:,:) :: t_flux          =>NULL() !< sensible heat flux (W/m2)
+  real, pointer, dimension(:,:) :: melth           =>NULL() !< sea ice and snow melt heat flux (W/m2)
+  real, pointer, dimension(:,:) :: meltw           =>NULL() !< water flux due to sea ice and snow melting (kg/m2/s)
   real, pointer, dimension(:,:) :: q_flux          =>NULL() !< specific humidity flux (kg/m2/s)
   real, pointer, dimension(:,:) :: salt_flux       =>NULL() !< salt flux (kg/m2/s)
   real, pointer, dimension(:,:) :: lw_flux         =>NULL() !< long wave radiation (W/m2)
@@ -166,7 +168,6 @@ module MOM_surface_forcing
   real, pointer, dimension(:,:) :: sw_flux_nir_dif =>NULL() !< diffuse Near InfraRed sw radiation (W/m2)
   real, pointer, dimension(:,:) :: lprec           =>NULL() !< mass flux of liquid precip (kg/m2/s)
   real, pointer, dimension(:,:) :: fprec           =>NULL() !< mass flux of frozen precip (kg/m2/s)
-  real, pointer, dimension(:,:) :: runoff          =>NULL() !< mass flux of liquid runoff (kg/m2/s)
   real, pointer, dimension(:,:) :: calving         =>NULL() !< mass flux of frozen runoff (kg/m2/s)
   real, pointer, dimension(:,:) :: ustar_berg      =>NULL() !< frictional velocity beneath icebergs (m/s)
   real, pointer, dimension(:,:) :: area_berg       =>NULL() !< area covered by icebergs(m2/m2)
@@ -443,6 +444,14 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, Time, G, CS, &
     if (associated(fluxes%sens)) &
       fluxes%sens(i,j) = G%mask2dT(i,j) * IOB%t_flux(i-i0,j-j0)
 
+    ! sea ice and snow melt heat flux (W/m2)
+    if (associated(fluxes%melth)) &
+      fluxes%melth(i,j) = G%mask2dT(i,j) * IOB%melth(i-i0,j-j0)
+
+    ! water flux due to sea ice and snow melt (kg/m2/s)
+    if (associated(fluxes%meltw)) &
+      fluxes%meltw(i,j) = G%mask2dT(i,j) * IOB%meltw(i-i0,j-j0)
+
     ! latent heat flux (W/m^2)
     if (associated(fluxes%latent)) &
       fluxes%latent(i,j) = G%mask2dT(i,j) * IOB%latent_flux(i-i0,j-j0)
@@ -474,7 +483,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, Time, G, CS, &
     sign_for_net_FW_bug = 1.
     if (CS%use_net_FW_adjustment_sign_bug) sign_for_net_FW_bug = -1.
     do j=js,je ; do i=is,ie
-      net_FW(i,j) = (((fluxes%lprec(i,j)   + fluxes%fprec(i,j)) + &
+      net_FW(i,j) = (((fluxes%lprec(i,j)   + fluxes%fprec(i,j) + fluxes%meltw(i,j)) + &
           (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j))) + &
           (fluxes%evap(i,j)    + fluxes%vprec(i,j)) ) * G%areaT(i,j)
       !   The following contribution appears to be calculating the volume flux of sea-ice
@@ -483,9 +492,12 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, Time, G, CS, &
       !   Bob thinks this is trying ensure the net fresh-water of the ocean + sea-ice system
       ! is constant.
       !   To do this correctly we will need a sea-ice melt field added to IOB. -AJA
-      if (associated(IOB%salt_flux) .and. (CS%ice_salt_concentration>0.0)) &
-          net_FW(i,j) = net_FW(i,j) + sign_for_net_FW_bug * G%areaT(i,j) * &
-          (IOB%salt_flux(i-i0,j-j0) / CS%ice_salt_concentration)
+      ! GMM: as stated above, the following is wrong. CIME deals with volume/mass and
+      ! heat from sea ice/snow via meltw and melth, respectively.
+      if (associated(fluxes%salt_flux) .and. (CS%ice_salt_concentration>0.0)) &
+          net_FW(i,j) = net_FW(i,j) + G%areaT(i,j) * &
+          (fluxes%salt_flux(i,j) / CS%ice_salt_concentration)
+
       net_FW2(i,j) = net_FW(i,j)/G%areaT(i,j)
     enddo; enddo
 
@@ -599,13 +611,14 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, CS)
         forces%p_surf(i,j) = forces%p_surf_full(i,j)
       endif
 
-      if (CS%use_limited_P_SSH) then
-        forces%p_surf_SSH => forces%p_surf
-      else
-        forces%p_surf_SSH => forces%p_surf_full
-      endif
-    end if
-  end do; end do
+    endif
+  enddo; enddo
+
+  if (CS%use_limited_P_SSH) then
+    forces%p_surf_SSH => forces%p_surf
+  else
+    forces%p_surf_SSH => forces%p_surf_full
+  endif
 
   ! GMM, CIME uses AGRID. All the BGRID_NE code can be cleaned later
   wind_stagger = AGRID
@@ -774,6 +787,8 @@ subroutine IOB_allocate(IOB, isc, iec, jsc, jec)
              IOB% u_flux (isc:iec,jsc:jec),          &
              IOB% v_flux (isc:iec,jsc:jec),          &
              IOB% t_flux (isc:iec,jsc:jec),          &
+             IOB% melth (isc:iec,jsc:jec),          &
+             IOB% meltw (isc:iec,jsc:jec),          &
              IOB% q_flux (isc:iec,jsc:jec),          &
              IOB% salt_flux (isc:iec,jsc:jec),       &
              IOB% lw_flux (isc:iec,jsc:jec),         &
@@ -783,7 +798,6 @@ subroutine IOB_allocate(IOB, isc, iec, jsc, jec)
              IOB% sw_flux_nir_dif (isc:iec,jsc:jec), &
              IOB% lprec (isc:iec,jsc:jec),           &
              IOB% fprec (isc:iec,jsc:jec),           &
-             IOB% runoff (isc:iec,jsc:jec),          &
              IOB% ustar_berg (isc:iec,jsc:jec),      &
              IOB% area_berg (isc:iec,jsc:jec),       &
              IOB% mass_berg (isc:iec,jsc:jec),       &
@@ -799,6 +813,8 @@ subroutine IOB_allocate(IOB, isc, iec, jsc, jec)
   IOB%u_flux          = 0.0
   IOB%v_flux          = 0.0
   IOB%t_flux          = 0.0
+  IOB%melth           = 0.0
+  IOB%meltw           = 0.0
   IOB%q_flux          = 0.0
   IOB%salt_flux       = 0.0
   IOB%lw_flux         = 0.0
@@ -808,7 +824,6 @@ subroutine IOB_allocate(IOB, isc, iec, jsc, jec)
   IOB%sw_flux_nir_dif = 0.0
   IOB%lprec           = 0.0
   IOB%fprec           = 0.0
-  IOB%runoff          = 0.0
   IOB%ustar_berg      = 0.0
   IOB%area_berg       = 0.0
   IOB%mass_berg       = 0.0
@@ -1319,7 +1334,11 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt)
   write(outunit,100) 'iobt%u_flux         ', mpp_chksum( iobt%u_flux         )
   write(outunit,100) 'iobt%v_flux         ', mpp_chksum( iobt%v_flux         )
   write(outunit,100) 'iobt%t_flux         ', mpp_chksum( iobt%t_flux         )
+  write(outunit,100) 'iobt%melth          ', mpp_chksum( iobt%melth          )
+  write(outunit,100) 'iobt%meltw          ', mpp_chksum( iobt%meltw          )
   write(outunit,100) 'iobt%q_flux         ', mpp_chksum( iobt%q_flux         )
+  write(outunit,100) 'iobt%rofl_flux      ', mpp_chksum( iobt%rofl_flux      )
+  write(outunit,100) 'iobt%rofi_flux      ', mpp_chksum( iobt%rofi_flux      )
   write(outunit,100) 'iobt%salt_flux      ', mpp_chksum( iobt%salt_flux      )
   write(outunit,100) 'iobt%lw_flux        ', mpp_chksum( iobt%lw_flux        )
   write(outunit,100) 'iobt%sw_flux_vis_dir', mpp_chksum( iobt%sw_flux_vis_dir)
@@ -1328,7 +1347,6 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt)
   write(outunit,100) 'iobt%sw_flux_nir_dif', mpp_chksum( iobt%sw_flux_nir_dif)
   write(outunit,100) 'iobt%lprec          ', mpp_chksum( iobt%lprec          )
   write(outunit,100) 'iobt%fprec          ', mpp_chksum( iobt%fprec          )
-  write(outunit,100) 'iobt%runoff         ', mpp_chksum( iobt%runoff         )
   write(outunit,100) 'iobt%calving        ', mpp_chksum( iobt%calving        )
   write(outunit,100) 'iobt%p              ', mpp_chksum( iobt%p              )
   if (associated(iobt%ustar_berg)) &

config_src/mct_driver/ocn_cap_methods.F90

@@ -33,7 +33,7 @@ subroutine ocn_import(x2o, ind, grid, ice_ocean_boundary, ocean_public, logunit,
   real(kind=8), optional        , intent(in)    :: c1, c2, c3, c4     !< Coeffs. used in the shortwave decomposition
 
   ! Local variables
-  integer         :: i, j, ig, jg, isc, iec, jsc, jec  ! Grid indices
+  integer         :: i, j, isc, iec, jsc, jec  ! Grid indices
   integer         :: k
   integer         :: day, secs, rc
   type(ESMF_time) :: currTime
@@ -44,16 +44,17 @@ subroutine ocn_import(x2o, ind, grid, ice_ocean_boundary, ocean_public, logunit,
 
   k = 0
   do j = jsc, jec
-    jg = j + grid%jsc - jsc
     do i = isc, iec
-      ig = i + grid%jsc - isc
       k = k + 1 ! Increment position within gindex
 
+      ! rotate taux and tauy from true zonal/meridional to local coordinates
       ! taux
-      ice_ocean_boundary%u_flux(i,j) = x2o(ind%x2o_Foxx_taux,k)
+      ice_ocean_boundary%u_flux(i,j) = GRID%cos_rot(i,j) * x2o(ind%x2o_Foxx_taux,k) &
+                                      + GRID%sin_rot(i,j) * x2o(ind%x2o_Foxx_tauy,k)
 
       ! tauy
-      ice_ocean_boundary%v_flux(i,j) = x2o(ind%x2o_Foxx_tauy,k)
+      ice_ocean_boundary%v_flux(i,j) = GRID%cos_rot(i,j) * x2o(ind%x2o_Foxx_tauy,k) &
+                                      - GRID%sin_rot(i,j) * x2o(ind%x2o_Foxx_taux,k)
 
       ! liquid precipitation (rain)
       ice_ocean_boundary%lprec(i,j) = x2o(ind%x2o_Faxa_rain,k)
@@ -65,41 +66,47 @@ subroutine ocn_import(x2o, ind, grid, ice_ocean_boundary, ocean_public, logunit,
       ice_ocean_boundary%lw_flux(i,j) = (x2o(ind%x2o_Faxa_lwdn,k) + x2o(ind%x2o_Foxx_lwup,k))
 
       ! specific humitidy flux
-      ice_ocean_boundary%q_flux(i,j) = x2o(ind%x2o_Foxx_evap,k) !???TODO: should this be a minus sign
+      ice_ocean_boundary%q_flux(i,j) = x2o(ind%x2o_Foxx_evap,k)
 
       ! sensible heat flux (W/m2)
-      ice_ocean_boundary%t_flux(i,j) = x2o(ind%x2o_Foxx_sen,k)  !???TODO: should this be a minus sign
+      ice_ocean_boundary%t_flux(i,j) = x2o(ind%x2o_Foxx_sen,k)
 
       ! latent heat flux (W/m^2)
-      ice_ocean_boundary%latent_flux(i,j) = x2o(ind%x2o_Foxx_lat,k) !???TODO: should this be a minus sign
+      ice_ocean_boundary%latent_flux(i,j) = x2o(ind%x2o_Foxx_lat,k)
+
+      ! snow&ice melt heat flux  (W/m^2)
+      ice_ocean_boundary%melth(i,j) = x2o(ind%x2o_Fioi_melth,k)
+
+      ! water flux from snow&ice melt (kg/m2/s)
+      ice_ocean_boundary%meltw(i,j) = x2o(ind%x2o_Fioi_meltw,k)
 
       ! liquid runoff
-      ice_ocean_boundary%rofl_flux(i,j) = x2o(ind%x2o_Foxx_rofl,k) * GRID%mask2dT(ig,jg)
+      ice_ocean_boundary%rofl_flux(i,j) = x2o(ind%x2o_Foxx_rofl,k) * GRID%mask2dT(i,j)
 
       ! ice runoff
-      ice_ocean_boundary%rofi_flux(i,j) = x2o(ind%x2o_Foxx_rofi,k) * GRID%mask2dT(ig,jg)
+      ice_ocean_boundary%rofi_flux(i,j) = x2o(ind%x2o_Foxx_rofi,k) * GRID%mask2dT(i,j)
 
       ! surface pressure
-      ice_ocean_boundary%p(i,j) = x2o(ind%x2o_Sa_pslv,k) * GRID%mask2dT(ig,jg)
+      ice_ocean_boundary%p(i,j) = x2o(ind%x2o_Sa_pslv,k) * GRID%mask2dT(i,j)
 
       ! salt flux (minus sign needed here -GMM)
-      ice_ocean_boundary%salt_flux(i,j) = -x2o(ind%x2o_Fioi_salt,k) * GRID%mask2dT(ig,jg)
+      ice_ocean_boundary%salt_flux(i,j) = -x2o(ind%x2o_Fioi_salt,k) * GRID%mask2dT(i,j)
 
       ! 1) visible, direct shortwave  (W/m2)
       ! 2) visible, diffuse shortwave (W/m2)
       ! 3) near-IR, direct shortwave  (W/m2)
       ! 4) near-IR, diffuse shortwave (W/m2)
       if (present(c1) .and. present(c2) .and. present(c3) .and. present(c4)) then
         ! Use runtime coefficients to decompose net short-wave heat flux into 4 components
-        ice_ocean_boundary%sw_flux_vis_dir(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c1 * GRID%mask2dT(ig,jg)
-        ice_ocean_boundary%sw_flux_vis_dif(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c2 * GRID%mask2dT(ig,jg)
-        ice_ocean_boundary%sw_flux_nir_dir(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c3 * GRID%mask2dT(ig,jg)
-        ice_ocean_boundary%sw_flux_nir_dif(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c4 * GRID%mask2dT(ig,jg)
+        ice_ocean_boundary%sw_flux_vis_dir(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c1 * GRID%mask2dT(i,j)
+        ice_ocean_boundary%sw_flux_vis_dif(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c2 * GRID%mask2dT(i,j)
+        ice_ocean_boundary%sw_flux_nir_dir(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c3 * GRID%mask2dT(i,j)
+        ice_ocean_boundary%sw_flux_nir_dif(i,j) = x2o(ind%x2o_Foxx_swnet,k) * c4 * GRID%mask2dT(i,j)
       else
-        ice_ocean_boundary%sw_flux_vis_dir(i,j) = x2o(ind%x2o_Faxa_swvdr,k) * GRID%mask2dT(ig,jg)
-        ice_ocean_boundary%sw_flux_vis_dif(i,j) = x2o(ind%x2o_Faxa_swvdf,k) * GRID%mask2dT(ig,jg)
-        ice_ocean_boundary%sw_flux_nir_dir(i,j) = x2o(ind%x2o_Faxa_swndr,k) * GRID%mask2dT(ig,jg)
-        ice_ocean_boundary%sw_flux_nir_dif(i,j) = x2o(ind%x2o_Faxa_swndf,k) * GRID%mask2dT(ig,jg)
+        ice_ocean_boundary%sw_flux_vis_dir(i,j) = x2o(ind%x2o_Faxa_swvdr,k) * GRID%mask2dT(i,j)
+        ice_ocean_boundary%sw_flux_vis_dif(i,j) = x2o(ind%x2o_Faxa_swvdf,k) * GRID%mask2dT(i,j)
+        ice_ocean_boundary%sw_flux_nir_dir(i,j) = x2o(ind%x2o_Faxa_swndr,k) * GRID%mask2dT(i,j)
+        ice_ocean_boundary%sw_flux_nir_dif(i,j) = x2o(ind%x2o_Faxa_swndf,k) * GRID%mask2dT(i,j)
       endif
     enddo
   enddo
@@ -116,6 +123,8 @@ subroutine ocn_import(x2o, ind, grid, ice_ocean_boundary, ocean_public, logunit,
         write(logunit,F01)'import: day, secs, j, i, lwrad           = ',day,secs,j,i,ice_ocean_boundary%lw_flux(i,j)
         write(logunit,F01)'import: day, secs, j, i, q_flux          = ',day,secs,j,i,ice_ocean_boundary%q_flux(i,j)
         write(logunit,F01)'import: day, secs, j, i, t_flux          = ',day,secs,j,i,ice_ocean_boundary%t_flux(i,j)
+        write(logunit,F01)'import: day, secs, j, i, melth           = ',day,secs,j,i,ice_ocean_boundary%melth(i,j)
+        write(logunit,F01)'import: day, secs, j, i, meltw           = ',day,secs,j,i,ice_ocean_boundary%meltw(i,j)
         write(logunit,F01)'import: day, secs, j, i, latent_flux     = ',&
                           day,secs,j,i,ice_ocean_boundary%latent_flux(i,j)
         write(logunit,F01)'import: day, secs, j, i, runoff          = ',&
@@ -152,6 +161,8 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
 
   ! Local variables
   real, dimension(grid%isd:grid%ied,grid%jsd:grid%jed) :: ssh !< Local copy of sea_lev with updated halo
+  real, dimension(grid%isd:grid%ied,grid%jsd:grid%jed) :: sshx!< Zonal SSH gradient, local coordinate.
+  real, dimension(grid%isd:grid%ied,grid%jsd:grid%jed) :: sshy!< Meridional SSH gradient, local coordinate.
   integer :: i, j, n, ig, jg  !< Grid indices
   real    :: slp_L, slp_R, slp_C, slope, u_min, u_max
   real :: I_time_int  !< The inverse of coupling time interval in s-1.
@@ -174,8 +185,13 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
       ! surface temperature in Kelvin
       o2x(ind%o2x_So_t, n) = ocn_public%t_surf(ig,jg) * grid%mask2dT(i,j)
       o2x(ind%o2x_So_s, n) = ocn_public%s_surf(ig,jg) * grid%mask2dT(i,j)
-      o2x(ind%o2x_So_u, n) = ocn_public%u_surf(ig,jg) * grid%mask2dT(i,j)
-      o2x(ind%o2x_So_v, n) = ocn_public%v_surf(ig,jg) * grid%mask2dT(i,j)
+      ! rotate ocn current from local tripolar grid to true zonal/meridional (inverse transformation)
+      o2x(ind%o2x_So_u, n) = (grid%cos_rot(i,j) * ocn_public%u_surf(ig,jg) - &
+                              grid%sin_rot(i,j) * ocn_public%v_surf(ig,jg)) * grid%mask2dT(i,j)
+      o2x(ind%o2x_So_v, n) = (grid%cos_rot(i,j) * ocn_public%v_surf(ig,jg) + &
+                              grid%sin_rot(i,j) * ocn_public%u_surf(ig,jg)) * grid%mask2dT(i,j)
+
+      ! boundary layer depth (m)
       o2x(ind%o2x_So_bldepth, n) = ocn_public%OBLD(ig,jg) * grid%mask2dT(i,j)
       ! ocean melt and freeze potential (o2x_Fioo_q), W m-2
       if (ocn_public%frazil(ig,jg) > 0.0) then
@@ -197,9 +213,7 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
   call pass_var(ssh, grid%domain)
 
   ! d/dx ssh
-  n = 0
   do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec
-    n = n+1
     ! This is a simple second-order difference
     ! o2x(ind%o2x_So_dhdx, n) = 0.5 * (ssh(i+1,j) - ssh(i-1,j)) * grid%IdxT(i,j) * grid%mask2dT(i,j)
     ! This is a PLM slope which might be less prone to the A-grid null mode
@@ -219,14 +233,12 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
       ! larger extreme values.
       slope = 0.0
     endif
-    o2x(ind%o2x_So_dhdx, n) = slope * grid%IdxT(i,j) * grid%mask2dT(i,j)
-    if (grid%mask2dT(i,j)==0.) o2x(ind%o2x_So_dhdx, n) = 0.0
+    sshx(i,j) = slope * grid%IdxT(i,j) * grid%mask2dT(i,j)
+    if (grid%mask2dT(i,j)==0.) sshx(i,j) = 0.0
   enddo; enddo
 
   ! d/dy ssh
-  n = 0
   do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec
-    n = n+1
     ! This is a simple second-order difference
     ! o2x(ind%o2x_So_dhdy, n) = 0.5 * (ssh(i,j+1) - ssh(i,j-1)) * grid%IdyT(i,j) * grid%mask2dT(i,j)
     ! This is a PLM slope which might be less prone to the A-grid null mode
@@ -237,7 +249,6 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
     if (grid%mask2dCv(i,J+1)==0.) slp_R = 0.
 
     slp_C = 0.5 * (slp_L + slp_R)
-    !write(6,*)'slp_L, slp_R,i,j,slp_L*slp_R', slp_L, slp_R,i,j,slp_L*slp_R
     if ((slp_L * slp_R) > 0.0) then
       ! This limits the slope so that the edge values are bounded by the
       ! two cell averages spanning the edge.
@@ -249,8 +260,16 @@ subroutine ocn_export(ind, ocn_public, grid, o2x, dt_int, ncouple_per_day)
       ! larger extreme values.
       slope = 0.0
     endif
-    o2x(ind%o2x_So_dhdy, n) = slope * grid%IdyT(i,j) * grid%mask2dT(i,j)
-    if (grid%mask2dT(i,j)==0.) o2x(ind%o2x_So_dhdy, n) = 0.0
+    sshy(i,j) = slope * grid%IdyT(i,j) * grid%mask2dT(i,j)
+    if (grid%mask2dT(i,j)==0.) sshy(i,j) = 0.0
+  enddo; enddo
+
+  ! rotate ssh gradients from local coordinates to true zonal/meridional (inverse transformation)
+  n = 0
+  do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec
+    n = n+1
+    o2x(ind%o2x_So_dhdx, n) = grid%cos_rot(i,j) * sshx(i,j) - grid%sin_rot(i,j) * sshy(i,j)
+    o2x(ind%o2x_So_dhdy, n) = grid%cos_rot(i,j) * sshy(i,j) + grid%sin_rot(i,j) * sshx(i,j)
   enddo; enddo
 
 end subroutine ocn_export

config_src/mct_driver/ocn_comp_mct.F90

@@ -773,8 +773,6 @@ end subroutine ocean_model_init_sfc
 !! mi, mass of ice (kg/m2)
 !!
 !! Variables in the coupler that are **NOT** used in MOM6 (i.e., no corresponding field in fluxes):
-!! x2o_Fioi_melth, heat flux from snow & ice melt (W/m2)
-!! x2o_Fioi_meltw, snow melt flux (kg/m2/s)
 !! x2o_Si_ifrac, fractional ice wrt ocean
 !! x2o_So_duu10n, 10m wind speed squared (m^2/s^2)
 !! x2o_Sa_co2prog, bottom atm level prognostic CO2

config_src/mct_driver/ocn_cpl_indices.F90

@@ -41,7 +41,7 @@ module ocn_cpl_indices
     integer :: x2o_Faxa_swndr    !< near-IR, direct shortwave (W/m2)
     integer :: x2o_Faxa_swndf    !< near-IR, direct shortwave (W/m2)
     integer :: x2o_Fioi_melth    !< Heat flux from snow & ice melt (W/m2)
-    integer :: x2o_Fioi_meltw    !< Snow melt flux (kg/m2/s)
+    integer :: x2o_Fioi_meltw    !< Water flux from sea ice and snow melt (kg/m2/s)
     integer :: x2o_Fioi_bcpho    !< Black Carbon hydrophobic release from sea ice component
     integer :: x2o_Fioi_bcphi    !< Black Carbon hydrophilic release from sea ice component
     integer :: x2o_Fioi_flxdst   !< Dust release from sea ice component