Changes following code review

-- cleaned up a lot of comments and whitespace
-- used source argument in more allocate statements, and deallocated more
   arrays
-- 3D diags now have zl:mean in cell_methods attribute
-- marbl_instances%domain%kmt is set once (during initialization)

config_src/drivers/nuopc_cap/mom_cap_methods.F90

@@ -287,7 +287,6 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
   ! nhx deposition
   !----
   if (associated(ice_ocean_boundary%nhx_dep)) then
-    ice_ocean_boundary%nhx_dep(:,:) = 0._ESMF_KIND_R8
     call state_getimport(importState, 'Faxa_ndep',  &
         isc, iec, jsc, jec, ice_ocean_boundary%nhx_dep(:,:), &
         areacor=med2mod_areacor, esmf_ind=1, rc=rc)
@@ -298,7 +297,6 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
     ! noy deposition
     !----
   if (associated(ice_ocean_boundary%noy_dep)) then
-    ice_ocean_boundary%noy_dep(:,:) = 0._ESMF_KIND_R8
     call state_getimport(importState, 'Faxa_ndep',  &
         isc, iec, jsc, jec, ice_ocean_boundary%noy_dep(:,:), &
         areacor=med2mod_areacor, esmf_ind=2, rc=rc)
@@ -311,13 +309,11 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
   ! in which the pointer(s) will not be associated
   !----
   if (associated(ice_ocean_boundary%atm_co2_prog)) then
-    ice_ocean_boundary%atm_co2_prog(:,:) = 0._ESMF_KIND_R8
     call state_getimport(importState, 'Sa_co2prog',  &
         isc, iec, jsc, jec, ice_ocean_boundary%atm_co2_prog(:,:), rc=rc)
   if (ChkErr(rc,__LINE__,u_FILE_u)) return
   endif
   if (associated(ice_ocean_boundary%atm_co2_diag)) then
-    ice_ocean_boundary%atm_co2_diag(:,:) = 0._ESMF_KIND_R8
     call state_getimport(importState, 'Sa_co2diag',  &
         isc, iec, jsc, jec, ice_ocean_boundary%atm_co2_diag(:,:), rc=rc)
     if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
@@ -362,7 +358,6 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,
   ! dust flux from sea ice
   !----
   if (associated(ice_ocean_boundary%seaice_dust_flux)) then
-    ice_ocean_boundary%seaice_dust_flux(:,:) = 0._ESMF_KIND_R8
     call state_getimport(importState, 'Fioi_flxdst',  &
         isc, iec, jsc, jec, ice_ocean_boundary%seaice_dust_flux, &
         areacor=med2mod_areacor, rc=rc)
@@ -921,25 +916,25 @@ subroutine State_GetImport_2d(state, fldname, isc, iec, jsc, jec, output, do_sum
       ! determine output array and apply area correction if present
       n = 0
       do j = jsc,jec
-         do i = isc,iec
-            n = n + 1
-            if (do_sum_loc) then
-               if (present(areacor)) then
-                  output(i,j)  = output(i,j) + dataPtr1d(n) * areacor(n)
-               else
-                  output(i,j)  = output(i,j) + dataPtr1d(n)
-               endif
+        do i = isc,iec
+          n = n + 1
+          if (do_sum_loc) then
+            if (present(areacor)) then
+              output(i,j)  = output(i,j) + dataPtr1d(n) * areacor(n)
             else
-               if (present(areacor)) then
-                  output(i,j)  = dataPtr1d(n) * areacor(n)
-               else
-                  output(i,j)  = dataPtr1d(n)
-               endif
+              output(i,j)  = output(i,j) + dataPtr1d(n)
             endif
-         enddo
+          else
+            if (present(areacor)) then
+              output(i,j)  = dataPtr1d(n) * areacor(n)
+            else
+              output(i,j)  = dataPtr1d(n)
+            endif
+          endif
+        enddo
       enddo
 
-    else if (geomtype == ESMF_GEOMTYPE_GRID) then
+    elseif (geomtype == ESMF_GEOMTYPE_GRID) then
 
       call state_getfldptr(state, trim(fldname), dataptr2d, rc)
       if (ChkErr(rc,__LINE__,u_FILE_u)) return
@@ -948,15 +943,15 @@ subroutine State_GetImport_2d(state, fldname, isc, iec, jsc, jec, output, do_sum
       lbnd2 = lbound(dataPtr2d,2)
 
       do j = jsc, jec
-         j1 = j + lbnd2 - jsc
-         do i = isc, iec
-            i1 = i + lbnd1 - isc
-            if (do_sum_loc) then
-               output(i,j) = output(i,j) + dataPtr2d(i1,j1)
-            else
-               output(i,j) = dataPtr2d(i1,j1)
-            endif
-         enddo
+        j1 = j + lbnd2 - jsc
+        do i = isc, iec
+          i1 = i + lbnd1 - isc
+          if (do_sum_loc) then
+            output(i,j) = output(i,j) + dataPtr2d(i1,j1)
+          else
+            output(i,j) = dataPtr2d(i1,j1)
+          endif
+        enddo
       enddo
 
     endif

config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90

@@ -543,7 +543,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, &
 
   if (OS%fluxes%fluxes_used) then
 
-    ! MNL: Necessary for averaging some MARBL forcing fields
+    ! enable_averages() is necessary to post forcing fields to diagnostics
     call enable_averages(dt_coupling, OS%Time + Ocean_coupling_time_step, OS%diag)
 
     if (do_thermo) &
@@ -789,7 +789,7 @@ subroutine ocean_model_end(Ocean_sfc, Ocean_state, Time, write_restart)
   type(time_type),         intent(in)    :: Time        !< The model time, used for writing restarts.
   logical,                 intent(in)    :: write_restart !< true => write restart file
 
-  if (write_restart)call ocean_model_save_restart(Ocean_state, Time)
+  if (write_restart) call ocean_model_save_restart(Ocean_state, Time)
   call diag_mediator_end(Time, Ocean_state%diag, end_diag_manager=.true.)
   call MOM_end(Ocean_state%MOM_CSp)
   if (Ocean_state%use_ice_shelf) call ice_shelf_end(Ocean_state%Ice_shelf_CSp)

config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90

@@ -39,7 +39,7 @@ module MOM_surface_forcing_nuopc
 use user_revise_forcing,  only : user_revise_forcing_CS
 use iso_fortran_env,      only : int64
 use MARBL_forcing_mod,    only : marbl_forcing_CS, MARBL_forcing_init
-use MARBL_forcing_mod,    only : convert_marbl_IOB_to_forcings
+use MARBL_forcing_mod,    only : convert_driver_fields_to_forcings
 
 implicit none ; private
 
@@ -585,11 +585,11 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
 
   ! Copy MARBL-specific IOB fields into fluxes; also set some MARBL-specific forcings to other values
   ! (constants, values from netCDF, etc)
-  call convert_marbl_IOB_to_forcings(IOB%atm_fine_dust_flux, IOB%atm_coarse_dust_flux, &
-                                     IOB%seaice_dust_flux, IOB%atm_bc_flux, IOB%seaice_bc_flux, &
-                                     IOB%nhx_dep, IOB%noy_dep, IOB%atm_co2_prog, IOB%atm_co2_diag, &
-                                     IOB%afracr, IOB%swnet_afracr, IOB%ifrac_n, IOB%swpen_ifrac_n, &
-                                     Time, G, US, i0, j0, fluxes, CS%marbl_forcing_CSp)
+  call convert_driver_fields_to_forcings(IOB%atm_fine_dust_flux, IOB%atm_coarse_dust_flux, &
+                                         IOB%seaice_dust_flux, IOB%atm_bc_flux, IOB%seaice_bc_flux, &
+                                         IOB%nhx_dep, IOB%noy_dep, IOB%atm_co2_prog, IOB%atm_co2_diag, &
+                                         IOB%afracr, IOB%swnet_afracr, IOB%ifrac_n, IOB%swpen_ifrac_n, &
+                                         Time, G, US, i0, j0, fluxes, CS%marbl_forcing_CSp)
 
   ! wave to ocean coupling
   if ( associated(IOB%lamult)) then

config_src/drivers/solo_driver/MOM_surface_forcing.F90

@@ -57,7 +57,7 @@ module MOM_surface_forcing
 use dumbbell_surface_forcing,    only : dumbbell_surface_forcing_init, dumbbell_surface_forcing_CS
 use dumbbell_surface_forcing, only    : dumbbell_buoyancy_forcing
 use MARBL_forcing_mod,    only : marbl_forcing_CS, MARBL_forcing_init
-use MARBL_forcing_mod,    only : convert_marbl_IOB_to_forcings
+use MARBL_forcing_mod,    only : convert_driver_fields_to_forcings
 
 implicit none ; private
 
@@ -1571,7 +1571,7 @@ subroutine MARBL_forcing_from_data_override(fluxes, day, G, US, CS)
   real, pointer, dimension(:,:) :: noy_dep              =>NULL() !< Nitrogen deposition [kg/m^2/s ~> RZ/T]
   integer :: isc, iec, jsc, jec
 
-  ! Necessary null pointers for arguments to convert_marbl_IOB_to_forcings()
+  ! Necessary null pointers for arguments to convert_driver_fields_to_forcings()
   ! Since they are null, MARBL will not use multiple ice categories
   real, pointer, dimension(:,:)   :: afracr        =>NULL()
   real, pointer, dimension(:,:)   :: swnet_afracr  =>NULL()
@@ -1595,17 +1595,9 @@ subroutine MARBL_forcing_from_data_override(fluxes, day, G, US, CS)
              atm_bc_flux (isc:iec,jsc:jec),          &
              seaice_bc_flux (isc:iec,jsc:jec),       &
              nhx_dep (isc:iec,jsc:jec),              &
-             noy_dep (isc:iec,jsc:jec))
-
-  atm_co2_prog(:,:) = 0.0
-  atm_co2_diag(:,:) = 0.0
-  atm_fine_dust_flux(:,:) = 0.0
-  atm_coarse_dust_flux(:,:) = 0.0
-  atm_bc_flux(:,:) = 0.0
-  seaice_dust_flux(:,:) = 0.0
-  seaice_bc_flux(:,:) = 0.0
-  nhx_dep(:,:) = 0.0
-  noy_dep(:,:) = 0.0
+             noy_dep (isc:iec,jsc:jec),              &
+             source=0.0)
+
 
   ! fluxes used directly as MARBL inputs
   ! (should be scaled)
@@ -1614,7 +1606,7 @@ subroutine MARBL_forcing_from_data_override(fluxes, day, G, US, CS)
 
   ! fluxes used to compute MARBL inputs
   ! These are kept in physical units, and will be scaled appropriately in
-  ! convert_marbl_IOB_to_forcings()
+  ! convert_driver_fields_to_forcings()
   call data_override(G%Domain, 'atm_co2_prog', atm_co2_prog, day)
   call data_override(G%Domain, 'atm_co2_diag', atm_co2_diag, day)
   call data_override(G%Domain, 'atm_fine_dust_flux', atm_fine_dust_flux, day)
@@ -1625,11 +1617,11 @@ subroutine MARBL_forcing_from_data_override(fluxes, day, G, US, CS)
   call data_override(G%Domain, 'nhx_dep', nhx_dep, day)
   call data_override(G%Domain, 'noy_dep', noy_dep, day)
 
-  call convert_marbl_IOB_to_forcings(atm_fine_dust_flux, atm_coarse_dust_flux, &
-                                     seaice_dust_flux, atm_bc_flux, seaice_bc_flux, &
-                                     nhx_dep, noy_dep, atm_co2_prog, atm_co2_diag, &
-                                     afracr, swnet_afracr, ifrac_n, swpen_ifrac_n, &
-                                     day, G, US, 0, 0, fluxes, CS%marbl_forcing_CSp)
+  call convert_driver_fields_to_forcings(atm_fine_dust_flux, atm_coarse_dust_flux, &
+                                         seaice_dust_flux, atm_bc_flux, seaice_bc_flux, &
+                                         nhx_dep, noy_dep, atm_co2_prog, atm_co2_diag, &
+                                         afracr, swnet_afracr, ifrac_n, swpen_ifrac_n, &
+                                         day, G, US, 0, 0, fluxes, CS%marbl_forcing_CSp)
 
   deallocate ( atm_co2_prog,         &
                atm_co2_diag,         &

config_src/external/MARBL/README.md

@@ -1,5 +1,5 @@
-GFDL_ocean_BGC
-==============
+MARBL
+=====
 
 These APIs reflect those for the MARBL library available at https://github.com/marbl-ecosys/MARBL
 

src/core/MOM_forcing_type.F90

@@ -3617,6 +3617,12 @@ subroutine deallocate_forcing_type(fluxes)
   if (associated(fluxes%mass_berg))            deallocate(fluxes%mass_berg)
   if (associated(fluxes%ice_fraction))         deallocate(fluxes%ice_fraction)
   if (associated(fluxes%u10_sqr))              deallocate(fluxes%u10_sqr)
+  if (associated(fluxes%noy_dep))              deallocate(fluxes%noy_dep)
+  if (associated(fluxes%nhx_dep))              deallocate(fluxes%nhx_dep)
+  if (associated(fluxes%atm_co2))              deallocate(fluxes%atm_co2)
+  if (associated(fluxes%atm_alt_co2))          deallocate(fluxes%atm_alt_co2)
+  if (associated(fluxes%dust_flux))            deallocate(fluxes%dust_flux)
+  if (associated(fluxes%iron_flux))            deallocate(fluxes%iron_flux)
   if (associated(fluxes%fracr_cat))            deallocate(fluxes%fracr_cat)
   if (associated(fluxes%qsw_cat))              deallocate(fluxes%qsw_cat)
 

src/tracer/MARBL_forcing_mod.F90

@@ -22,7 +22,7 @@ module MARBL_forcing_mod
 #include <MOM_memory.h>
 
 public :: MARBL_forcing_init
-public :: convert_marbl_IOB_to_forcings
+public :: convert_driver_fields_to_forcings
 
 !> Data type used to store diagnostic index returned from register_diag_field()
 !! For the forcing fields that can be written via post_data()
@@ -161,7 +161,7 @@ subroutine MARBL_forcing_init(G, US, param_file, diag, day, inputdir, use_marbl,
     endif
 
     ! Register diagnostic fields for outputing forcing values
-    ! These fields are posted from convert_marbl_IOB_to_forcings(), and they are received
+    ! These fields are posted from convert_driver_fields_to_forcings(), and they are received
     ! in physical units so no conversion is necessary here.
     CS%diag_ids%atm_fine_dust = register_diag_field("ocean_model", "ATM_FINE_DUST_FLUX_CPL", &
         CS%diag%axesT1, & ! T=> tracer grid? 1 => no vertical grid
@@ -183,11 +183,11 @@ subroutine MARBL_forcing_init(G, US, param_file, diag, day, inputdir, use_marbl,
   end subroutine MARBL_forcing_init
 
   ! Note: ice fraction and u10_sqr are handled in mom_surface_forcing because of CFCs
-  subroutine convert_marbl_IOB_to_forcings(atm_fine_dust_flux, atm_coarse_dust_flux, &
-                                           seaice_dust_flux, atm_bc_flux, seaice_bc_flux, &
-                                           nhx_dep, noy_dep, atm_co2_prog, atm_co2_diag, &
-                                           afracr, swnet_afracr, ifrac_n, &
-                                           swpen_ifrac_n, Time, G, US, i0, j0, fluxes, CS)
+  subroutine convert_driver_fields_to_forcings(atm_fine_dust_flux, atm_coarse_dust_flux, &
+                                               seaice_dust_flux, atm_bc_flux, seaice_bc_flux, &
+                                               nhx_dep, noy_dep, atm_co2_prog, atm_co2_diag, &
+                                               afracr, swnet_afracr, ifrac_n, &
+                                               swpen_ifrac_n, Time, G, US, i0, j0, fluxes, CS)
 
     real, dimension(:,:),   pointer, intent(in)    :: atm_fine_dust_flux   !< atmosphere fine dust flux from IOB
                                                                            !! [kg m-2 s-1]
@@ -223,7 +223,7 @@ subroutine convert_marbl_IOB_to_forcings(atm_fine_dust_flux, atm_coarse_dust_flu
     real :: seaice_fe_bioavail_frac  !< TODO: define this (local) term
     real :: iron_flux_conversion     !< TODO: define this (local) term
     real :: ndep_conversion          !< Combination of unit conversion factors for rescaling
-                                     !! nitrogen deposition [kg(N) m-2 s-1 ~> mol L-2 T-1]
+                                     !! nitrogen deposition [kg(N) m-2 s-1 ~> mol m-3 Z T-1]
 
     if (.not. CS%use_marbl_tracers) return
 
@@ -373,6 +373,6 @@ subroutine convert_marbl_IOB_to_forcings(atm_fine_dust_flux, atm_coarse_dust_flu
       enddo; enddo
     endif
 
-  end subroutine convert_marbl_IOB_to_forcings
+  end subroutine convert_driver_fields_to_forcings
 
-end module MARBL_forcing_mod
+end module MARBL_forcing_mod