+Created allocate_surface_state

  Created a new subroutine, allocate_surface_state, to allocate the variables
that MOM6 shares as its surface values, including the possibility of using an
optional coupler_1d_bc_type argument to allocate the space associated with
surface fluxes during initialization, which will then allow MOM6 to offer
diagnostics associated with these fluxes.  All answers are bitwise identical.

config_src/coupled_driver/ocean_model_MOM.F90

@@ -35,7 +35,7 @@ module ocean_model_mod
 !</OVERVIEW>
 
 use MOM, only : initialize_MOM, step_MOM, MOM_control_struct, MOM_end
-use MOM, only : calculate_surface_state, finish_MOM_initialization
+use MOM, only : calculate_surface_state, allocate_surface_state, finish_MOM_initialization
 use MOM, only : step_offline
 use MOM_constants, only : CELSIUS_KELVIN_OFFSET, hlf
 use MOM_diag_mediator, only : diag_ctrl, enable_averaging, disable_averaging
@@ -240,7 +240,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn)
   character(len=40)  :: mdl = "ocean_model_init"  ! This module's name.
   character(len=48)  :: stagger
   integer :: secs, days
-  integer :: is, ie, js, je
   type(param_file_type) :: param_file !< A structure to parse for run-time parameters
   logical :: offline_tracer_mode
 
@@ -262,16 +261,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn)
   OS%C_p = OS%MOM_CSp%tv%C_p
   OS%fluxes%C_p = OS%MOM_CSp%tv%C_p
 
-  is = OS%grid%isc ; ie = OS%grid%iec ; js = OS%grid%jsc ; je = OS%grid%jec
-  if (present(gas_fields_ocn)) then
-    call coupler_type_spawn(gas_fields_ocn, OS%state%tr_fields, &
-                            (/is,is,ie,ie/), (/js,js,je,je/))
-  elseif (coupler_type_initialized(Ocean_sfc%fields)) then
-    !### I think that this is never used.
-    call coupler_type_spawn(Ocean_sfc%fields, OS%state%tr_fields, &
-                            (/is,is,ie,ie/), (/js,js,je,je/))
-  endif
-
   ! Read all relevant parameters and write them to the model log.
   call log_version(param_file, mdl, version, "")
   call get_param(param_file, mdl, "RESTART_CONTROL", OS%Restart_control, &
@@ -327,23 +316,28 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn)
     call get_param(param_file, mdl, "KV_ICEBERG",  OS%kv_iceberg, &
                  "The viscosity of the icebergs",  units="m2 s-1",default=1.0e10)
     call get_param(param_file, mdl, "DENSITY_ICEBERGS",  OS%density_iceberg, &
-                  "A typical density of icebergs.", units="kg m-3", default=917.0)
+                 "A typical density of icebergs.", units="kg m-3", default=917.0)
     call get_param(param_file, mdl, "LATENT_HEAT_FUSION", OS%latent_heat_fusion, &
                  "The latent heat of fusion.", units="J/kg", default=hlf)
     call get_param(param_file, mdl, "BERG_AREA_THRESHOLD", OS%berg_area_threshold, &
                  "Fraction of grid cell which iceberg must occupy, so that fluxes \n"//&
-                  "below berg are set to zero. Not applied for negative \n"//&
+                 "below berg are set to zero. Not applied for negative \n"//&
                  " values.", units="non-dim", default=-1.0)
   endif
 
   OS%press_to_z = 1.0/(Rho0*G_Earth)
 
+  !   Consider using a run-time flag to determine whether to do the diagnostic
+  ! vertical integrals, since the related 3-d sums are not negligible in cost.
+  call allocate_surface_state(OS%state, OS%grid, OS%MOM_CSp%use_temperature, &
+                              do_integrals=.true., gas_fields_ocn=gas_fields_ocn)
+
   call surface_forcing_init(Time_in, OS%grid, param_file, OS%MOM_CSp%diag, &
                             OS%forcing_CSp, OS%restore_salinity, OS%restore_temp)
 
   if (OS%use_ice_shelf)  then
-     call initialize_ice_shelf(param_file, OS%grid, OS%Time, OS%ice_shelf_CSp, &
-                               OS%MOM_CSp%diag, OS%fluxes)
+    call initialize_ice_shelf(param_file, OS%grid, OS%Time, OS%ice_shelf_CSp, &
+                              OS%MOM_CSp%diag, OS%fluxes)
   endif
   if (OS%icebergs_apply_rigid_boundary)  then
     !call allocate_forcing_type(OS%grid, OS%fluxes, iceberg=.true.)
@@ -891,9 +885,6 @@ subroutine convert_state_to_ocean_type(state, Ocean_sfc, G, use_conT_absS, &
     if (.not.coupler_type_initialized(Ocean_sfc%fields)) then
       call MOM_error(FATAL, "convert_state_to_ocean_type: "//&
                "Ocean_sfc%fields has not been initialized.")
-!      call coupler_type_spawn(state%tr_fields, Ocean_sfc%fields, &
-!                              (/isc_bnd,isc_bnd,iec_bnd,iec_bnd/), &
-!                              (/jsc_bnd,jsc_bnd,jec_bnd,jec_bnd/) )
     endif
     call coupler_type_copy_data(state%tr_fields, Ocean_sfc%fields)
   endif

src/core/MOM.F90

@@ -36,7 +36,7 @@ module MOM
 use MOM_diag_mediator,        only : diag_set_state_ptrs, diag_update_remap_grids
 use MOM_diag_mediator,        only : disable_averaging, post_data, safe_alloc_ptr
 use MOM_diag_mediator,        only : register_diag_field, register_static_field
-use MOM_diag_mediator,        only : register_scalar_field
+use MOM_diag_mediator,        only : register_scalar_field, get_diag_time_end
 use MOM_diag_mediator,        only : set_axes_info, diag_ctrl, diag_masks_set
 use MOM_domains,              only : MOM_domains_init, clone_MOM_domain
 use MOM_domains,              only : sum_across_PEs, pass_var, pass_vector
@@ -62,6 +62,7 @@ module MOM
 use MOM_time_manager,         only : operator(-), operator(>), operator(*), operator(/)
 use MOM_time_manager,         only : increment_date
 use MOM_unit_tests,           only : unit_tests
+use coupler_types_mod,        only : coupler_type_send_data, coupler_1d_bc_type, coupler_type_spawn
 
 ! MOM core modules
 use MOM_ALE,                   only : ALE_init, ALE_end, ALE_main, ALE_CS, adjustGridForIntegrity
@@ -433,6 +434,7 @@ module MOM
 public step_MOM
 public step_offline
 public MOM_end
+public allocate_surface_state
 public calculate_surface_state
 
 integer :: id_clock_ocean
@@ -1929,7 +1931,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo
   endif
 
   if (CS%bulkmixedlayer .or. CS%use_temperature) then
-     allocate(CS%Hml(isd:ied,jsd:jed)) ; CS%Hml(:,:) = 0.0
+    allocate(CS%Hml(isd:ied,jsd:jed)) ; CS%Hml(:,:) = 0.0
   endif
 
   if (CS%bulkmixedlayer) then
@@ -3228,6 +3230,8 @@ subroutine post_surface_diagnostics(CS, G, diag, state)
     call post_data(CS%id_speed, sfc_speed, diag, mask=G%mask2dT)
   endif
 
+  call coupler_type_send_data(state%tr_fields, get_diag_time_end(diag))
+
 end subroutine post_surface_diagnostics
 
 !> Offers the static fields in the ocean grid type
@@ -3434,6 +3438,64 @@ subroutine adjust_ssh_for_p_atm(CS, G, GV, ssh, p_atm)
 
 end subroutine adjust_ssh_for_p_atm
 
+!> This subroutine allocates the fields for the surface (return) properties of
+!! the ocean model.  Unused fields are unallocated.
+subroutine allocate_surface_state(state, G, use_temperature, do_integrals, &
+                                  gas_fields_ocn)
+  type(ocean_grid_type), intent(in)    :: G                !< ocean grid structure
+  type(surface),         intent(inout) :: state            !< ocean surface state type to be allocated.
+  logical,     optional, intent(in)    :: use_temperature  !< If true, allocate the space for thermodynamic variables.
+  logical,     optional, intent(in)    :: do_integrals     !< If true, allocate the space for vertically integrated fields.
+  type(coupler_1d_bc_type), &
+               optional, intent(in)    :: gas_fields_ocn   !< If present, this type describes the ocean
+                                              !! ocean and surface-ice fields that will participate
+                                              !! in the calculation of additional gas or other
+                                              !! tracer fluxes, and can be used to spawn related
+                                              !! internal variables in the ice model.
+
+  logical :: use_temp, alloc_integ
+  integer :: is, ie, js, je, isd, ied, jsd, jed
+  integer :: isdB, iedB, jsdB, jedB
+
+  is  = G%isc ; ie  = G%iec ; js  = G%jsc ; je  = G%jec
+  isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
+  isdB = G%isdB ; iedB = G%iedB; jsdB = G%jsdB ; jedB = G%jedB
+
+  use_temp = .true. ; if (present(use_temperature)) use_temp = use_temperature
+  alloc_integ = .true. ; if (present(do_integrals)) alloc_integ = do_integrals
+
+  if (state%arrays_allocated) return
+
+  if (use_temp) then
+    allocate(state%SST(isd:ied,jsd:jed)) ; state%SST(:,:) = 0.0
+    allocate(state%SSS(isd:ied,jsd:jed)) ; state%SSS(:,:) = 0.0
+  else
+    allocate(state%sfc_density(isd:ied,jsd:jed)) ; state%sfc_density(:,:) = 0.0
+  endif
+  allocate(state%sea_lev(isd:ied,jsd:jed)) ; state%sea_lev(:,:) = 0.0
+  allocate(state%Hml(isd:ied,jsd:jed)) ; state%Hml(:,:) = 0.0
+  allocate(state%u(IsdB:IedB,jsd:jed)) ; state%u(:,:) = 0.0
+  allocate(state%v(isd:ied,JsdB:JedB)) ; state%v(:,:) = 0.0
+
+  if (alloc_integ) then
+    ! Allocate structures for the vertically integrated ocean_mass, ocean_heat,
+    ! and ocean_salt.
+    allocate(state%ocean_mass(isd:ied,jsd:jed)) ; state%ocean_mass(:,:) = 0.0
+    if (use_temp) then
+      allocate(state%ocean_heat(isd:ied,jsd:jed)) ; state%ocean_heat(:,:) = 0.0
+      allocate(state%ocean_salt(isd:ied,jsd:jed)) ; state%ocean_salt(:,:) = 0.0
+    endif
+    allocate(state%salt_deficit(isd:ied,jsd:jed)) ; state%salt_deficit(:,:) = 0.0
+  endif
+
+  if (present(gas_fields_ocn)) &
+    call coupler_type_spawn(gas_fields_ocn, state%tr_fields, &
+                            (/isd,is,ie,ied/), (/jsd,js,je,jed/), as_needed=.true.)
+
+  state%arrays_allocated = .true.
+
+end subroutine allocate_surface_state
+
 !> This subroutine sets the surface (return) properties of the ocean
 !! model by setting the appropriate fields in state.  Unused fields
 !! are set to NULL or are unallocated.
@@ -3468,28 +3530,9 @@ subroutine calculate_surface_state(state, u, v, h, ssh, G, GV, CS)
   isdB = G%isdB ; iedB = G%iedB; jsdB = G%jsdB ; jedB = G%jedB
 
   if (.not.state%arrays_allocated) then
-    if (CS%use_temperature) then
-      allocate(state%SST(isd:ied,jsd:jed)) ; state%SST(:,:) = 0.0
-      allocate(state%SSS(isd:ied,jsd:jed)) ; state%SSS(:,:) = 0.0
-    else
-      allocate(state%sfc_density(isd:ied,jsd:jed)) ; state%sfc_density(:,:) = 0.0
-    endif
-    allocate(state%sea_lev(isd:ied,jsd:jed)) ; state%sea_lev(:,:) = 0.0
-    allocate(state%Hml(isd:ied,jsd:jed)) ; state%Hml(:,:) = 0.0
-    allocate(state%u(IsdB:IedB,jsd:jed)) ; state%u(:,:) = 0.0
-    allocate(state%v(isd:ied,JsdB:JedB)) ; state%v(:,:) = 0.0
-
-  ! Allocate structures for ocean_mass, ocean_heat, and ocean_salt.  This could
-  ! be wrapped in a run-time flag to disable it for economy, since the 3-d
-  ! sums are not negligible.
-    allocate(state%ocean_mass(isd:ied,jsd:jed)) ; state%ocean_mass(:,:) = 0.0
-    if (CS%use_temperature) then
-      allocate(state%ocean_heat(isd:ied,jsd:jed)) ; state%ocean_heat(:,:) = 0.0
-      allocate(state%ocean_salt(isd:ied,jsd:jed)) ; state%ocean_salt(:,:) = 0.0
-    endif
-    allocate(state%salt_deficit(isd:ied,jsd:jed)) ; state%salt_deficit(:,:) = 0.0
-
-    state%arrays_allocated = .true.
+    !  Consider using a run-time flag to determine whether to do the vertical
+    ! integrals, since the 3-d sums are not negligible in cost.
+    call allocate_surface_state(state, G, CS%use_temperature, do_integrals=.true.)
   endif
   state%frazil => CS%tv%frazil
   state%TempxPmE => CS%tv%TempxPmE