Merge branch 'dev/gfdl' into fix_diags_intwaves

config_src/infra/FMS1/MOM_interp_infra.F90

@@ -231,7 +231,7 @@ end subroutine time_interp_extern_3d
 
 !> initialize an external field
 integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose, &
-                                   threading, ierr, ignore_axis_atts )
+                                   threading, ierr, ignore_axis_atts, correct_leap_year_inconsistency )
 
   character(len=*),         intent(in)  :: file  !< The name of the file to read
   character(len=*),         intent(in)  :: fieldname !< The name of the field in the file
@@ -246,13 +246,20 @@ integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose,
   logical,        optional, intent(in)  :: ignore_axis_atts !< If present and true, do not issue a
                                                  !! fatal error if the axis Cartesian attribute is
                                                  !! not set to a recognized value.
+  logical,        optional, intent(in)  :: correct_leap_year_inconsistency !< If present and true,
+                                                 !! then if, (1) a calendar containing leap years
+                                                 !! is in use, and (2) the modulo time period of the
+                                                 !! data is an integer number of years, then map
+                                                 !! a model date of Feb 29. onto a common year on Feb. 28.
 
   if (present(MOM_Domain)) then
     init_extern_field = init_external_field(file, fieldname, domain=MOM_domain%mpp_domain, &
-             verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts)
+             verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, &
+             correct_leap_year_inconsistency=correct_leap_year_inconsistency)
   else
     init_extern_field = init_external_field(file, fieldname, domain=domain, &
-             verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts)
+             verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, &
+             correct_leap_year_inconsistency=correct_leap_year_inconsistency)
   endif
 
 end function init_extern_field

config_src/infra/FMS2/MOM_interp_infra.F90

@@ -231,7 +231,7 @@ end subroutine time_interp_extern_3d
 
 !> initialize an external field
 integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose, &
-                                   threading, ierr, ignore_axis_atts )
+                                   threading, ierr, ignore_axis_atts, correct_leap_year_inconsistency )
 
   character(len=*),         intent(in)  :: file  !< The name of the file to read
   character(len=*),         intent(in)  :: fieldname !< The name of the field in the file
@@ -246,13 +246,22 @@ integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose,
   logical,        optional, intent(in)  :: ignore_axis_atts !< If present and true, do not issue a
                                                  !! fatal error if the axis Cartesian attribute is
                                                  !! not set to a recognized value.
+  logical,        optional, intent(in)  :: correct_leap_year_inconsistency !< If present and true,
+                                                 !! then if, (1) a calendar containing leap years
+                                                 !! is in use, and (2) the modulo time period of the
+                                                 !! data is an integer number of years, then map
+                                                 !! a model date of Feb 29. onto a common year on Feb. 28.
+
+
 
   if (present(MOM_Domain)) then
     init_extern_field = init_external_field(file, fieldname, domain=MOM_domain%mpp_domain, &
-             verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts)
+             verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, &
+             correct_leap_year_inconsistency=correct_leap_year_inconsistency)
   else
     init_extern_field = init_external_field(file, fieldname, domain=domain, &
-             verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts)
+             verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, &
+             correct_leap_year_inconsistency=correct_leap_year_inconsistency)
   endif
 
 end function init_extern_field

src/core/MOM.F90

@@ -1276,7 +1276,13 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
   call enable_averages(dtdia, Time_end_thermo, CS%diag)
 
   if (associated(CS%odaCS)) then
-    call apply_oda_tracer_increments(US%T_to_s*dtdia, G, GV, tv, h, CS%odaCS)
+    if (CS%debug) then
+      call MOM_thermo_chksum("Pre-oda ", tv, G, US, haloshift=0)
+    endif
+    call apply_oda_tracer_increments(US%T_to_s*dtdia, Time_end_thermo, G, GV, tv, h, CS%odaCS)
+    if (CS%debug) then
+      call MOM_thermo_chksum("Post-oda ", tv, G, US, haloshift=0)
+    endif
   endif
 
   if (associated(fluxes%p_surf) .or. associated(fluxes%p_surf_full)) then
@@ -2807,7 +2813,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
                       (LEN_TRIM(dirs%input_filename) == 1))
 
   if (CS%ensemble_ocean) then
-    call init_oda(Time, G, GV, CS%odaCS)
+    call init_oda(Time, G, GV, CS%diag, CS%odaCS)
   endif
 
   !### This could perhaps go here instead of in finish_MOM_initialization?

src/core/MOM_barotropic.F90

@@ -2685,6 +2685,16 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
       ADp%diag_hv(i,J,k) = BT_cont%h_v(i,J,k)
     enddo ; enddo ; enddo
   endif
+  if (present(ADp) .and. (associated(ADp%visc_rem_u))) then
+    do k=1,nz ; do j=js,je ; do I=is-1,ie
+      ADp%visc_rem_u(I,j,k) = visc_rem_u(I,j,k)
+    enddo ; enddo ; enddo
+  endif
+  if (present(ADp) .and. (associated(ADp%visc_rem_u))) then
+    do k=1,nz ; do J=js-1,je ; do i=is,ie
+      ADp%visc_rem_v(i,J,k) = visc_rem_v(i,J,k)
+    enddo ; enddo ; enddo
+  endif
 
   if (G%nonblocking_updates) then
     if (find_etaav) call complete_group_pass(CS%pass_etaav, G%Domain)

src/core/MOM_dynamics_split_RK2.F90

@@ -169,10 +169,12 @@ module MOM_dynamics_split_RK2
   integer :: id_h_PFu    = -1, id_h_PFv    = -1
   integer :: id_hf_PFu_2d = -1, id_hf_PFv_2d = -1
   integer :: id_intz_PFu_2d = -1, id_intz_PFv_2d = -1
+  integer :: id_PFu_visc_rem = -1, id_PFv_visc_rem = -1
   ! integer :: id_hf_CAu    = -1, id_hf_CAv    = -1
   integer :: id_h_CAu    = -1, id_h_CAv    = -1
   integer :: id_hf_CAu_2d = -1, id_hf_CAv_2d = -1
   integer :: id_intz_CAu_2d = -1, id_intz_CAv_2d = -1
+  integer :: id_CAu_visc_rem = -1, id_CAv_visc_rem = -1
 
   ! Split scheme only.
   integer :: id_uav        = -1, id_vav        = -1
@@ -181,6 +183,7 @@ module MOM_dynamics_split_RK2
   integer :: id_h_u_BT_accel    = -1, id_h_v_BT_accel    = -1
   integer :: id_hf_u_BT_accel_2d = -1, id_hf_v_BT_accel_2d = -1
   integer :: id_intz_u_BT_accel_2d = -1, id_intz_v_BT_accel_2d = -1
+  integer :: id_u_BT_accel_visc_rem    = -1, id_v_BT_accel_visc_rem    = -1
   !>@}
 
   type(diag_ctrl), pointer       :: diag !< A structure that is used to regulate the
@@ -360,6 +363,12 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
   real, dimension(SZI_(G),SZJB_(G)) :: &
     intz_PFv_2d, intz_CAv_2d, intz_v_BT_accel_2d ! [H L T-2 ~> m2 s-2].
 
+  ! Diagnostics for momentum budget terms multiplied by visc_rem_[uv],
+  real, allocatable, dimension(:,:,:) :: &
+    PFu_visc_rem, PFv_visc_rem, & ! Pressure force accel. x visc_rem_[uv] [L T-2 ~> m s-2].
+    CAu_visc_rem, CAv_visc_rem, & ! Coriolis force accel. x visc_rem_[uv] [L T-2 ~> m s-2].
+    u_BT_accel_visc_rem, v_BT_accel_visc_rem ! barotropic correction accel. x visc_rem_[uv] [L T-2 ~> m s-2].
+
   real :: dt_pred   ! The time step for the predictor part of the baroclinic time stepping [T ~> s].
 
   logical :: dyn_p_surf
@@ -1108,6 +1117,61 @@ subroutine step_MOM_dyn_split_RK2(u, v, h, tv, visc, Time_local, dt, forces, p_s
     deallocate(h_v_BT_accel)
   endif
 
+  if (CS%id_PFu_visc_rem > 0) then
+    allocate(PFu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
+    PFu_visc_rem(:,:,:) = 0.0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      PFu_visc_rem(I,j,k) = CS%PFu(I,j,k) * CS%ADp%visc_rem_u(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_PFu_visc_rem, PFu_visc_rem, CS%diag)
+    deallocate(PFu_visc_rem)
+  endif
+  if (CS%id_PFv_visc_rem > 0) then
+    allocate(PFv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
+    PFv_visc_rem(:,:,:) = 0.0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      PFv_visc_rem(i,J,k) = CS%PFv(i,J,k) * CS%ADp%visc_rem_v(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_PFv_visc_rem, PFv_visc_rem, CS%diag)
+    deallocate(PFv_visc_rem)
+  endif
+  if (CS%id_CAu_visc_rem > 0) then
+    allocate(CAu_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
+    CAu_visc_rem(:,:,:) = 0.0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      CAu_visc_rem(I,j,k) = CS%CAu(I,j,k) * CS%ADp%visc_rem_u(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_CAu_visc_rem, CAu_visc_rem, CS%diag)
+    deallocate(CAu_visc_rem)
+  endif
+  if (CS%id_CAv_visc_rem > 0) then
+    allocate(CAv_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
+    CAv_visc_rem(:,:,:) = 0.0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      CAv_visc_rem(i,J,k) = CS%CAv(i,J,k) * CS%ADp%visc_rem_v(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_CAv_visc_rem, CAv_visc_rem, CS%diag)
+    deallocate(CAv_visc_rem)
+  endif
+  if (CS%id_u_BT_accel_visc_rem > 0) then
+    allocate(u_BT_accel_visc_rem(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke))
+    u_BT_accel_visc_rem(:,:,:) = 0.0
+    do k=1,nz ; do j=js,je ; do I=Isq,Ieq
+      u_BT_accel_visc_rem(I,j,k) = CS%u_accel_bt(I,j,k) * CS%ADp%visc_rem_u(I,j,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_u_BT_accel_visc_rem, u_BT_accel_visc_rem, CS%diag)
+    deallocate(u_BT_accel_visc_rem)
+  endif
+  if (CS%id_v_BT_accel_visc_rem > 0) then
+    allocate(v_BT_accel_visc_rem(G%isd:G%ied,G%JsdB:G%JedB,GV%ke))
+    v_BT_accel_visc_rem(:,:,:) = 0.0
+    do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
+      v_BT_accel_visc_rem(i,J,k) = CS%v_accel_bt(i,J,k) * CS%ADp%visc_rem_v(i,J,k)
+    enddo ; enddo ; enddo
+    call post_data(CS%id_v_BT_accel_visc_rem, v_BT_accel_visc_rem, CS%diag)
+    deallocate(v_BT_accel_visc_rem)
+  endif
+
   if (CS%debug) then
     call MOM_state_chksum("Corrector ", u, v, h, uh, vh, G, GV, US, symmetric=sym)
     call uvchksum("Corrector avg [uv]", u_av, v_av, G%HI, haloshift=1, symmetric=sym, scale=US%L_T_to_m_s)
@@ -1612,6 +1676,33 @@ subroutine initialize_dyn_split_RK2(u, v, h, uh, vh, eta, Time, G, GV, US, param
       'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
   if (CS%id_intz_v_BT_accel_2d > 0) call safe_alloc_ptr(CS%ADp%diag_hv,isd,ied,JsdB,JedB,nz)
 
+  CS%id_PFu_visc_rem = register_diag_field('ocean_model', 'PFu_visc_rem', diag%axesCuL, Time, &
+      'Zonal Pressure Force Acceleration multiplied by the viscous remnant', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if (CS%id_PFu_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
+  CS%id_PFv_visc_rem = register_diag_field('ocean_model', 'PFv_visc_rem', diag%axesCvL, Time, &
+      'Meridional Pressure Force Acceleration multiplied by the viscous remnant', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if(CS%id_PFv_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
+
+  CS%id_CAu_visc_rem = register_diag_field('ocean_model', 'CAu_visc_rem', diag%axesCuL, Time, &
+      'Zonal Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if (CS%id_CAu_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
+  CS%id_CAv_visc_rem = register_diag_field('ocean_model', 'CAv_visc_rem', diag%axesCvL, Time, &
+      'Meridional Coriolis and Advective Acceleration multiplied by the viscous remnant', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if(CS%id_CAv_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
+
+  CS%id_u_BT_accel_visc_rem = register_diag_field('ocean_model', 'u_BT_accel_visc_rem', diag%axesCuL, Time, &
+      'Barotropic Anomaly Zonal Acceleration multiplied by the viscous remnant', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if (CS%id_u_BT_accel_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
+  CS%id_v_BT_accel_visc_rem = register_diag_field('ocean_model', 'v_BT_accel_visc_rem', diag%axesCvL, Time, &
+      'Barotropic Anomaly Meridional Acceleration multiplied by the viscous remnant', 'm s-2', &
+      conversion=US%L_T2_to_m_s2)
+  if(CS%id_v_BT_accel_visc_rem > 0) call safe_alloc_ptr(CS%ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
+
   id_clock_Cor        = cpu_clock_id('(Ocean Coriolis & mom advection)', grain=CLOCK_MODULE)
   id_clock_continuity = cpu_clock_id('(Ocean continuity equation)',      grain=CLOCK_MODULE)
   id_clock_pres       = cpu_clock_id('(Ocean pressure force)',           grain=CLOCK_MODULE)

src/core/MOM_variables.F90

@@ -195,6 +195,9 @@ module MOM_variables
   real, pointer :: diag_hu(:,:,:) => NULL() !< layer thickness at u points
   real, pointer :: diag_hv(:,:,:) => NULL() !< layer thickness at v points
 
+  real, pointer :: visc_rem_u(:,:,:) => NULL() !< viscous remnant at u points
+  real, pointer :: visc_rem_v(:,:,:) => NULL() !< viscous remnant at v points
+
 end type accel_diag_ptrs
 
 !> Pointers to arrays with transports, which can later be used for derived diagnostics, like energy balances.

src/diagnostics/MOM_diagnostics.F90

@@ -2341,7 +2341,6 @@ subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, GV, CS)
     call safe_alloc_ptr(ADp%gradKEu,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(ADp%gradKEv,isd,ied,JsdB,JedB,nz)
   endif
-
   if (associated(CS%KE_visc)) then
     call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
     call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
@@ -2395,7 +2394,7 @@ subroutine MOM_diagnostics_end(CS, ADp, CDp)
   if (associated(CS%vhGM_Rlay))  deallocate(CS%vhGM_Rlay)
 
   if (associated(ADp%gradKEu))    deallocate(ADp%gradKEu)
-  if (associated(ADp%gradKEu))    deallocate(ADp%gradKEu)
+  if (associated(ADp%gradKEv))    deallocate(ADp%gradKEv)
   if (associated(ADp%du_dt_visc)) deallocate(ADp%du_dt_visc)
   if (associated(ADp%dv_dt_visc)) deallocate(ADp%dv_dt_visc)
   if (associated(ADp%du_dt_str))  deallocate(ADp%du_dt_str)