Merge branch 'dev/gfdl' into more_unit_fixes

src/core/MOM_CoriolisAdv.F90

@@ -25,6 +25,7 @@ module MOM_CoriolisAdv
 
 !> Control structure for mom_coriolisadv
 type, public :: CoriolisAdv_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   integer :: Coriolis_Scheme !< Selects the discretization for the Coriolis terms.
                              !! Valid values are:
                              !! - SADOURNY75_ENERGY - Sadourny, 1975
@@ -245,6 +246,9 @@ subroutine CorAdCalc(u, v, h, uh, vh, CAu, CAv, OBC, AD, G, GV, US, CS)
 !   v(is-1:ie+2,js-1:je+1), u(is-1:ie+1,js-1:je+2), h(is-1:ie+2,js-1:je+2),
 !   uh(is-1,ie,js:je+1) and vh(is:ie+1,js-1:je).
 
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+         "MOM_CoriolisAdv: Module must be initialized before it is used.")
+
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
   Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = GV%ke
   vol_neglect = GV%H_subroundoff * (1e-4 * US%m_to_L)**2
@@ -1123,6 +1127,7 @@ subroutine CoriolisAdv_init(Time, G, GV, US, param_file, diag, AD, CS)
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke
   IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
 
+  CS%initialized = .true.
   CS%diag => diag ; CS%Time => Time
 
   ! Read all relevant parameters and write them to the model log.

src/core/MOM_PressureForce_FV.F90

@@ -34,6 +34,7 @@ module MOM_PressureForce_FV
 
 !> Finite volume pressure gradient control structure
 type, public :: PressureForce_FV_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   logical :: tides          !< If true, apply tidal momentum forcing.
   real    :: Rho0           !< The density used in the Boussinesq
                             !! approximation [R ~> kg m-3].
@@ -163,6 +164,9 @@ subroutine PressureForce_FV_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_
   Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
   EOSdom(1) = Isq - (G%isd-1) ;  EOSdom(2) = G%iec+1 - (G%isd-1)
 
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+       "MOM_PressureForce_FV_nonBouss: Module must be initialized before it is used.")
+
   if (CS%Stanley_T2_det_coeff>=0.) call MOM_error(FATAL, &
        "MOM_PressureForce_FV_nonBouss: The Stanley parameterization is not yet"//&
        "implemented in non-Boussinesq mode.")
@@ -497,6 +501,9 @@ subroutine PressureForce_FV_Bouss(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm
   Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
   EOSdom(1) = Isq - (G%isd-1) ;  EOSdom(2) = G%iec+1 - (G%isd-1)
 
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+       "MOM_PressureForce_FV_Bouss: Module must be initialized before it is used.")
+
   use_p_atm = associated(p_atm)
   use_EOS = associated(tv%eqn_of_state)
   do i=Isq,Ieq+1 ; p0(i) = 0.0 ; enddo
@@ -809,6 +816,7 @@ subroutine PressureForce_FV_init(Time, G, GV, US, param_file, diag, CS, tides_CS
   character(len=40)  :: mdl  ! This module's name.
   logical :: use_ALE
 
+  CS%initialized = .true.
   CS%diag => diag ; CS%Time => Time
   if (present(tides_CSp)) &
     CS%tides_CSp => tides_CSp

src/core/MOM_PressureForce_Montgomery.F90

@@ -30,6 +30,7 @@ module MOM_PressureForce_Mont
 
 !> Control structure for the Montgomery potential form of pressure gradient
 type, public :: PressureForce_Mont_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   logical :: tides          !< If true, apply tidal momentum forcing.
   real    :: Rho0           !< The density used in the Boussinesq
                             !! approximation [R ~> kg m-3].
@@ -137,6 +138,9 @@ subroutine PressureForce_Mont_nonBouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pb
   is_split = present(pbce)
   use_EOS = associated(tv%eqn_of_state)
 
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+      "MOM_PressureForce_Mont: Module must be initialized before it is used.")
+
   if (use_EOS) then
     if (query_compressible(tv%eqn_of_state)) call MOM_error(FATAL, &
       "PressureForce_Mont_nonBouss: The Montgomery form of the pressure force "//&
@@ -422,6 +426,9 @@ subroutine PressureForce_Mont_Bouss(h, tv, PFu, PFv, G, GV, US, CS, p_atm, pbce,
   is_split = present(pbce)
   use_EOS = associated(tv%eqn_of_state)
 
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+       "MOM_PressureForce_Mont: Module must be initialized before it is used.")
+
   if (use_EOS) then
     if (query_compressible(tv%eqn_of_state)) call MOM_error(FATAL, &
       "PressureForce_Mont_Bouss: The Montgomery form of the pressure force "//&
@@ -829,6 +836,7 @@ subroutine PressureForce_Mont_init(Time, G, GV, US, param_file, diag, CS, tides_
 # include "version_variable.h"
   character(len=40)  :: mdl   ! This module's name.
 
+  CS%initialized = .true.
   CS%diag => diag ; CS%Time => Time
   if (present(tides_CSp)) &
     CS%tides_CSp => tides_CSp

src/core/MOM_barotropic.F90

@@ -694,6 +694,9 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, forces, pbce,
   integer :: ioff, joff
   integer :: l_seg
 
+  if (.not.CS%module_is_initialized) call MOM_error(FATAL, &
+      "btstep: Module MOM_barotropic must be initialized before it is used.")
+
   if (.not.CS%split) return
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
   Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
@@ -2764,6 +2767,9 @@ subroutine set_dtbt(G, GV, US, CS, eta, pbce, BT_cont, gtot_est, SSH_add)
   character(len=200) :: mesg
   integer :: i, j, k, is, ie, js, je, nz
 
+  if (.not.CS%module_is_initialized) call MOM_error(FATAL, &
+      "set_dtbt: Module MOM_barotropic must be initialized before it is used.")
+
   if (.not.CS%split) return
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
   MS%isdw = G%isd ; MS%iedw = G%ied ; MS%jsdw = G%jsd ; MS%jedw = G%jed
@@ -3298,6 +3304,9 @@ subroutine btcalc(h, G, GV, CS, h_u, h_v, may_use_default, OBC)
 
 !    This section interpolates thicknesses onto u & v grid points with the
 ! second order accurate estimate h = 2*(h+ * h-)/(h+ + h-).
+  if (.not.CS%module_is_initialized) call MOM_error(FATAL, &
+      "btcalc: Module MOM_barotropic must be initialized before it is used.")
+
   if (.not.CS%split) return
 
   use_default = .false.
@@ -4186,6 +4195,9 @@ subroutine bt_mass_source(h, eta, set_cor, G, GV, CS)
                               ! thicknesses [H ~> m or kg m-2].
   integer :: is, ie, js, je, nz, i, j, k
 
+  if (.not.CS%module_is_initialized) call MOM_error(FATAL, "bt_mass_source: "// &
+        "Module MOM_barotropic must be initialized before it is used.")
+
   if (.not.CS%split) return
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke

src/core/MOM_continuity_PPM.F90

@@ -26,6 +26,7 @@ module MOM_continuity_PPM
 
 !> Control structure for mom_continuity_ppm
 type, public :: continuity_PPM_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   type(diag_ctrl), pointer :: diag !< Diagnostics control structure.
   logical :: upwind_1st      !< If true, use a first-order upwind scheme.
   logical :: monotonic       !< If true, use the Colella & Woodward monotonic
@@ -134,6 +135,9 @@ subroutine continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, uhbt, vh
 
   h_min = GV%Angstrom_H
 
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+         "MOM_continuity_PPM: Module must be initialized before it is used.")
+
   x_first = (MOD(G%first_direction,2) == 0)
 
   if (present(visc_rem_u) .neqv. present(visc_rem_v)) call MOM_error(FATAL, &
@@ -2202,6 +2206,8 @@ subroutine continuity_PPM_init(Time, G, GV, US, param_file, diag, CS)
   real :: tol_eta_m  ! An unscaled version of tol_eta [m].
   character(len=40)  :: mdl = "MOM_continuity_PPM" ! This module's name.
 
+  CS%initialized = .true.
+
 ! Read all relevant parameters and write them to the model log.
   call log_version(param_file, mdl, version, "")
   call get_param(param_file, mdl, "MONOTONIC_CONTINUITY", CS%monotonic, &

src/diagnostics/MOM_diagnostics.F90

@@ -49,6 +49,7 @@ module MOM_diagnostics
 
 !> The control structure for the MOM_diagnostics module
 type, public :: diagnostics_CS ; private
+  logical :: initialized = .false.     !< True if this control structure has been initialized.
   real :: mono_N2_column_fraction = 0. !< The lower fraction of water column over which N2 is limited as
                                        !! monotonic for the purposes of calculating the equivalent
                                        !! barotropic wave speed.
@@ -267,6 +268,9 @@ subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
   if (loc(CS)==0) call MOM_error(FATAL, &
          "calculate_diagnostic_fields: Module must be initialized before used.")
 
+  if (.not. CS%initialized) call MOM_error(FATAL, &
+         "calculate_diagnostic_fields: Module must be initialized before used.")
+
   call calculate_derivs(dt, G, CS)
 
   if (dt > 0.0) then
@@ -1245,6 +1249,9 @@ subroutine register_time_deriv(lb, f_ptr, deriv_ptr, CS)
   integer :: m      !< New index of deriv_ptr in CS%deriv
   integer :: ub(3)  !< Upper index bound of f_ptr, based on shape.
 
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+         "register_time_deriv: Module must be initialized before it is used.")
+
   if (CS%num_time_deriv >= MAX_FIELDS_) then
     call MOM_error(WARNING,"MOM_diagnostics:  Attempted to register more than " // &
                    "MAX_FIELDS_ diagnostic time derivatives via register_time_deriv.")
@@ -1604,6 +1611,8 @@ subroutine MOM_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag
   isd  = G%isd  ; ied  = G%ied  ; jsd  = G%jsd  ; jed  = G%jed ; nz = GV%ke
   IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
 
+  CS%initialized = .true.
+
   CS%diag => diag
   use_temperature = associated(tv%T)
   call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., &

src/diagnostics/MOM_sum_output.F90

@@ -59,6 +59,8 @@ module MOM_sum_output
 
 !> The control structure for the MOM_sum_output module
 type, public :: sum_output_CS ; private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
+
   type(Depth_List)              :: DL !< The sorted depth list.
 
   integer, allocatable, dimension(:) :: lH
@@ -159,6 +161,8 @@ subroutine MOM_sum_output_init(G, GV, US, param_file, directory, ntrnc, &
   endif
   allocate(CS)
 
+  CS%initialized = .true.
+
   ! Read all relevant parameters and write them to the model log.
   call log_version(param_file, mdl, version, "")
   call get_param(param_file, mdl, "CALCULATE_APE", CS%do_APE_calc, &
@@ -488,6 +492,9 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci
   if (.not.associated(CS)) call MOM_error(FATAL, &
          "write_energy: Module must be initialized before it is used.")
 
+  if (.not.CS%initialized) call MOM_error(FATAL, &
+         "write_energy: Module must be initialized before it is used.")
+
   do j=js,je ; do i=is,ie
     areaTm(i,j) = G%mask2dT(i,j)*G%areaT(i,j)
   enddo ; enddo

src/diagnostics/MOM_wave_speed.F90

@@ -26,6 +26,7 @@ module MOM_wave_speed
 
 !> Control structure for MOM_wave_speed
 type, public :: wave_speed_CS ; private
+  logical :: initialized = .false.     !< True if this control structure has been initialized.
   logical :: use_ebt_mode = .false.    !< If true, calculate the equivalent barotropic wave speed instead
                                        !! of the first baroclinic wave speed.
                                        !! This parameter controls the default behavior of wave_speed() which
@@ -146,6 +147,9 @@ subroutine wave_speed(h, tv, G, GV, US, cg1, CS, full_halos, use_ebt_mode, mono_
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
 
+  if (.not. CS%initialized) call MOM_error(FATAL, "MOM_wave_speed: "// &
+           "Module must be initialized before it is used.")
+
   if (present(full_halos)) then ; if (full_halos) then
     is = G%isd ; ie = G%ied ; js = G%jsd ; je = G%jed
   endif ; endif
@@ -722,6 +726,11 @@ subroutine wave_speeds(h, tv, G, GV, US, nmodes, cn, CS, full_halos)
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
 
+  if (present(CS)) then
+    if (.not. CS%initialized) call MOM_error(FATAL, "MOM_wave_speed: "// &
+           "Module must be initialized before it is used.")
+  endif
+
   if (present(full_halos)) then ; if (full_halos) then
     is = G%isd ; ie = G%ied ; js = G%jsd ; je = G%jed
   endif ; endif
@@ -1185,6 +1194,8 @@ subroutine wave_speed_init(CS, use_ebt_mode, mono_N2_column_fraction, mono_N2_de
 # include "version_variable.h"
   character(len=40)  :: mdl = "MOM_wave_speed"  ! This module's name.
 
+  CS%initialized = .true.
+
   ! Write all relevant parameters to the model log.
   call log_version(mdl, version)
 

src/diagnostics/MOM_wave_structure.F90

@@ -36,6 +36,7 @@ module MOM_wave_structure
 
 !> The control structure for the MOM_wave_structure module
 type, public :: wave_structure_CS ; !private
+  logical :: initialized = .false. !< True if this control structure has been initialized.
   type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
                                    !! regulate the timing of diagnostic output.
   real, allocatable, dimension(:,:,:) :: w_strct
@@ -191,6 +192,9 @@ subroutine wave_structure(h, tv, G, GV, US, cn, ModeNum, freq, CS, En, full_halo
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
   I_a_int = 1/a_int
 
+  if (.not. CS%initialized) call MOM_error(FATAL, "MOM_wave_structure: "// &
+         "Module must be initialized before it is used.")
+
   if (present(full_halos)) then ; if (full_halos) then
     is = G%isd ; ie = G%ied ; js = G%jsd ; je = G%jed
   endif ; endif
@@ -728,6 +732,8 @@ subroutine wave_structure_init(Time, G, GV, param_file, diag, CS)
 
   isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke
 
+  CS%initialized = .true.
+
   call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_X", CS%int_tide_source_x, &
                  "X Location of generation site for internal tide", default=1.)
   call get_param(param_file, mdl, "INTERNAL_TIDE_SOURCE_Y", CS%int_tide_source_y, &