diff --git a/config_src/coupled_driver/MOM_surface_forcing.F90 b/config_src/coupled_driver/MOM_surface_forcing.F90 index 02b54daefe..57eb9cfcbc 100644 --- a/config_src/coupled_driver/MOM_surface_forcing.F90 +++ b/config_src/coupled_driver/MOM_surface_forcing.F90 @@ -50,107 +50,106 @@ module MOM_surface_forcing public ice_ocn_bnd_type_chksum public forcing_save_restart - -! surface_forcing_CS is a structure containing pointers to the forcing fields -! which may be used to drive MOM. All fluxes are positive downward. +!> surface_forcing_CS is a structure containing pointers to the forcing fields +!! which may be used to drive MOM. All fluxes are positive downward. type, public :: surface_forcing_CS ; private - integer :: wind_stagger ! AGRID, BGRID_NE, or CGRID_NE (integer values - ! from MOM_domains) to indicate the staggering of - ! the winds that are being provided in calls to - ! update_ocean_model. - logical :: use_temperature ! If true, temp and saln used as state variables + integer :: wind_stagger !< AGRID, BGRID_NE, or CGRID_NE (integer values + !! from MOM_domains) to indicate the staggering of + !! the winds that are being provided in calls to + !! update_ocean_model. + logical :: use_temperature !< If true, temp and saln used as state variables real :: wind_stress_multiplier !< A multiplier applied to incoming wind stress (nondim). - real :: Rho0 ! Boussinesq reference density (kg/m^3) - real :: area_surf = -1.0 ! total ocean surface area (m^2) - real :: latent_heat_fusion ! latent heat of fusion (J/kg) - real :: latent_heat_vapor ! latent heat of vaporization (J/kg) - - real :: max_p_surf ! maximum surface pressure that can be - ! exerted by the atmosphere and floating sea-ice, - ! in Pa. This is needed because the FMS coupling - ! structure does not limit the water that can be - ! frozen out of the ocean and the ice-ocean heat - ! fluxes are treated explicitly. - logical :: use_limited_P_SSH ! If true, return the sea surface height with - ! the correction for the atmospheric (and sea-ice) - ! pressure limited by max_p_surf instead of the - ! full atmospheric pressure. The default is true. - - real :: gust_const ! constant unresolved background gustiness for ustar (Pa) - logical :: read_gust_2d ! If true, use a 2-dimensional gustiness supplied - ! from an input file. + real :: Rho0 !< Boussinesq reference density (kg/m^3) + real :: area_surf = -1.0 !< Total ocean surface area (m^2) + real :: latent_heat_fusion !< Latent heat of fusion (J/kg) + real :: latent_heat_vapor !< Latent heat of vaporization (J/kg) + + real :: max_p_surf !< The maximum surface pressure that can be + !! exerted by the atmosphere and floating sea-ice, + !! in Pa. This is needed because the FMS coupling + !! structure does not limit the water that can be + !! frozen out of the ocean and the ice-ocean heat + !! fluxes are treated explicitly. + logical :: use_limited_P_SSH !< If true, return the sea surface height with + !! the correction for the atmospheric (and sea-ice) + !! pressure limited by max_p_surf instead of the + !! full atmospheric pressure. The default is true. + logical :: approx_net_mass_src !< If true, use the net mass sources from the ice-ocean boundary + !! type without any further adjustments to drive the ocean dynamics. + !! The actual net mass source may differ due to corrections. + + real :: gust_const !< Constant unresolved background gustiness for ustar (Pa) + logical :: read_gust_2d !< If true, use a 2-dimensional gustiness supplied from an input file. + real, pointer, dimension(:,:) :: & + TKE_tidal => NULL() !< Turbulent kinetic energy introduced to the bottom boundary layer + !! by drag on the tidal flows, in W m-2. + real, pointer, dimension(:,:) :: & + gust => NULL() !< A spatially varying unresolved background gustiness that + !! contributes to ustar (Pa). gust is used when read_gust_2d is true. real, pointer, dimension(:,:) :: & - TKE_tidal => NULL(), & ! turbulent kinetic energy introduced to the - ! bottom boundary layer by drag on the tidal flows, - ! in W m-2. - gust => NULL(), & ! spatially varying unresolved background - ! gustiness that contributes to ustar (Pa). - ! gust is used when read_gust_2d is true. - ustar_tidal => NULL() ! tidal contribution to the bottom friction velocity (m/s) - real :: cd_tides ! drag coefficient that applies to the tides (nondimensional) - real :: utide ! constant tidal velocity to use if read_tideamp - ! is false, in m s-1. - logical :: read_tideamp ! If true, spatially varying tidal amplitude read from a file. - - logical :: rigid_sea_ice ! If true, sea-ice exerts a rigidity that acts - ! to damp surface deflections (especially surface - ! gravity waves). The default is false. - real :: Kv_sea_ice ! viscosity in sea-ice that resists sheared vertical motions (m^2/s) - real :: density_sea_ice ! typical density of sea-ice (kg/m^3). The value is - ! only used to convert the ice pressure into - ! appropriate units for use with Kv_sea_ice. - real :: rigid_sea_ice_mass ! A mass per unit area of sea-ice beyond which - ! sea-ice viscosity becomes effective, in kg m-2, - ! typically of order 1000 kg m-2. - logical :: allow_flux_adjustments ! If true, use data_override to obtain flux adjustments - - real :: Flux_const ! piston velocity for surface restoring (m/s) - logical :: salt_restore_as_sflux ! If true, SSS restore as salt flux instead of water flux - logical :: adjust_net_srestore_to_zero ! adjust srestore to zero (for both salt_flux or vprec) - logical :: adjust_net_srestore_by_scaling ! adjust srestore w/o moving zero contour - logical :: adjust_net_fresh_water_to_zero ! adjust net surface fresh-water (w/ restoring) to zero - logical :: use_net_FW_adjustment_sign_bug ! use the wrong sign when adjusting net FW - logical :: adjust_net_fresh_water_by_scaling ! adjust net surface fresh-water w/o moving zero contour - logical :: mask_srestore_under_ice ! If true, use an ice mask defined by frazil - ! criteria for salinity restoring. - real :: ice_salt_concentration ! salt concentration for sea ice (kg/kg) - logical :: mask_srestore_marginal_seas ! if true, then mask SSS restoring in marginal seas - real :: max_delta_srestore ! maximum delta salinity used for restoring - real :: max_delta_trestore ! maximum delta sst used for restoring - real, pointer, dimension(:,:) :: basin_mask => NULL() ! mask for SSS restoring by basin - - type(diag_ctrl), pointer :: diag ! structure to regulate diagnostic output timing - character(len=200) :: inputdir ! directory where NetCDF input files are - character(len=200) :: salt_restore_file ! filename for salt restoring data - character(len=30) :: salt_restore_var_name ! name of surface salinity in salt_restore_file - logical :: mask_srestore ! if true, apply a 2-dimensional mask to the surface - ! salinity restoring fluxes. The masking file should be - ! in inputdir/salt_restore_mask.nc and the field should - ! be named 'mask' - real, pointer, dimension(:,:) :: srestore_mask => NULL() ! mask for SSS restoring - character(len=200) :: temp_restore_file ! filename for sst restoring data - character(len=30) :: temp_restore_var_name ! name of surface temperature in temp_restore_file - logical :: mask_trestore ! if true, apply a 2-dimensional mask to the surface - ! temperature restoring fluxes. The masking file should be - ! in inputdir/temp_restore_mask.nc and the field should - ! be named 'mask' - real, pointer, dimension(:,:) :: trestore_mask => NULL() ! mask for SST restoring - integer :: id_srestore = -1 ! id number for time_interp_external. - integer :: id_trestore = -1 ! id number for time_interp_external. - - ! Diagnostics handles - type(forcing_diags), public :: handles + ustar_tidal => NULL() !< Tidal contribution to the bottom friction velocity (m/s) + real :: cd_tides !< Drag coefficient that applies to the tides (nondimensional) + real :: utide !< Constant tidal velocity to use if read_tideamp is false, in m s-1. + logical :: read_tideamp !< If true, spatially varying tidal amplitude read from a file. + + logical :: rigid_sea_ice !< If true, sea-ice exerts a rigidity that acts to damp surface + !! deflections (especially surface gravity waves). The default is false. + real :: Kv_sea_ice !< Viscosity in sea-ice that resists sheared vertical motions (m^2/s) + real :: density_sea_ice !< Typical density of sea-ice (kg/m^3). The value is only used to convert + !! the ice pressure into appropriate units for use with Kv_sea_ice. + real :: rigid_sea_ice_mass !< A mass per unit area of sea-ice beyond which sea-ice viscosity + !! becomes effective, in kg m-2, typically of order 1000 kg m-2. + logical :: allow_flux_adjustments !< If true, use data_override to obtain flux adjustments + + logical :: restore_salt !< If true, the coupled MOM driver adds a term to restore surface + !! salinity to a specified value. + logical :: restore_temp !< If true, the coupled MOM driver adds a term to restore sea + !! surface temperature to a specified value. + real :: Flux_const !< Piston velocity for surface restoring (m/s) + logical :: salt_restore_as_sflux !< If true, SSS restore as salt flux instead of water flux + logical :: adjust_net_srestore_to_zero !< Adjust srestore to zero (for both salt_flux or vprec) + logical :: adjust_net_srestore_by_scaling !< Adjust srestore w/o moving zero contour + logical :: adjust_net_fresh_water_to_zero !< Adjust net surface fresh-water (with restoring) to zero + logical :: use_net_FW_adjustment_sign_bug !< Use the wrong sign when adjusting net FW + logical :: adjust_net_fresh_water_by_scaling !< Adjust net surface fresh-water w/o moving zero contour + logical :: mask_srestore_under_ice !< If true, use an ice mask defined by frazil criteria + !! for salinity restoring. + real :: ice_salt_concentration !< Salt concentration for sea ice (kg/kg) + logical :: mask_srestore_marginal_seas !< If true, then mask SSS restoring in marginal seas + real :: max_delta_srestore !< Maximum delta salinity used for restoring + real :: max_delta_trestore !< Maximum delta sst used for restoring + real, pointer, dimension(:,:) :: basin_mask => NULL() !< Mask for surface salinity restoring by basin + + type(diag_ctrl), pointer :: diag => NULL() !< Structure to regulate diagnostic output timing + character(len=200) :: inputdir !< Directory where NetCDF input files are + character(len=200) :: salt_restore_file !< Filename for salt restoring data + character(len=30) :: salt_restore_var_name !< Name of surface salinity in salt_restore_file + logical :: mask_srestore !< If true, apply a 2-dimensional mask to the surface + !! salinity restoring fluxes. The masking file should be + !! in inputdir/salt_restore_mask.nc and the field should + !! be named 'mask' + real, pointer, dimension(:,:) :: srestore_mask => NULL() !< mask for SSS restoring + character(len=200) :: temp_restore_file !< Filename for sst restoring data + character(len=30) :: temp_restore_var_name !< Name of surface temperature in temp_restore_file + logical :: mask_trestore !< If true, apply a 2-dimensional mask to the surface + !! temperature restoring fluxes. The masking file should be + !! in inputdir/temp_restore_mask.nc and the field should + !! be named 'mask' + real, pointer, dimension(:,:) :: trestore_mask => NULL() !< Mask for SST restoring + integer :: id_srestore = -1 !< An id number for time_interp_external. + integer :: id_trestore = -1 !< An id number for time_interp_external. + + type(forcing_diags), public :: handles !< Diagnostics handles !### type(ctrl_forcing_CS), pointer :: ctrl_forcing_CSp => NULL() - type(MOM_restart_CS), pointer :: restart_CSp => NULL() - type(user_revise_forcing_CS), pointer :: urf_CS => NULL() + type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< A pointer to the restart control structure + type(user_revise_forcing_CS), pointer :: urf_CS => NULL() !< A control structure for user forcing revisions end type surface_forcing_CS -! ice_ocean_boundary_type is a structure corresponding to forcing, but with -! the elements, units, and conventions that exactly conform to the use for -! MOM-based coupled models. +!> ice_ocean_boundary_type is a structure corresponding to forcing, but with the elements, units, +!! and conventions that exactly conform to the use for MOM6-based coupled models. type, public :: ice_ocean_boundary_type real, pointer, dimension(:,:) :: u_flux =>NULL() !< i-direction wind stress (Pa) real, pointer, dimension(:,:) :: v_flux =>NULL() !< j-direction wind stress (Pa) @@ -179,25 +178,23 @@ module MOM_surface_forcing !! ice-shelves, expressed as a coefficient !! for divergence damping, as determined !! outside of the ocean model in (m3/s) - integer :: xtype !< The type of the exchange - REGRID, REDIST or DIRECT - type(coupler_2d_bc_type) :: fluxes !< A structure that may contain an array of - !! named fields used for passive tracer fluxes. - integer :: wind_stagger = -999 !< A flag indicating the spatial discretization of - !! wind stresses. This flag may be set by the - !! flux-exchange code, based on what the sea-ice - !! model is providing. Otherwise, the value from - !! the surface_forcing_CS is used. + integer :: xtype !< The type of the exchange - REGRID, REDIST or DIRECT + type(coupler_2d_bc_type) :: fluxes !< A structure that may contain an array of named fields + !! used for passive tracer fluxes. + integer :: wind_stagger = -999 !< A flag indicating the spatial discretization of wind stresses. + !! This flag may be set by the flux-exchange code, based on what + !! the sea-ice model is providing. Otherwise, the value from + !! the surface_forcing_CS is used. end type ice_ocean_boundary_type -integer :: id_clock_forcing +integer :: id_clock_forcing !< A CPU time clock contains !> This subroutine translates the Ice_ocean_boundary_type into a MOM !! thermodynamic forcing type, including changes of units, sign conventions, !! and putting the fields into arrays with MOM-standard halos. -subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, & - sfc_state, restore_salt, restore_temp) +subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, sfc_state) type(ice_ocean_boundary_type), & target, intent(in) :: IOB !< An ice-ocean boundary type with fluxes to drive !! the ocean in a coupled model @@ -212,9 +209,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, & !! previous call to surface_forcing_init. type(surface), intent(in) :: sfc_state !< A structure containing fields that describe the !! surface state of the ocean. - logical, optional, intent(in) :: restore_salt !< If true, salinity is restored to a target value. - logical, optional, intent(in) :: restore_temp !< If true, temperature is restored to a target value. - real, dimension(SZI_(G),SZJ_(G)) :: & data_restore, & ! The surface value toward which to restore (g/kg or degC) @@ -234,10 +228,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, & integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isr, ier, jsr, jer integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd - logical :: restore_salinity ! local copy of the argument restore_salt, if it - ! is present, or false (no restoring) otherwise. - logical :: restore_sst ! local copy of the argument restore_temp, if it - ! is present, or false (no restoring) otherwise. real :: delta_sss ! temporary storage for sss diff from restoring value real :: delta_sst ! temporary storage for sst diff from restoring value @@ -264,11 +254,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, & fluxes%netFWGlobalAdj = 0.0 fluxes%netFWGlobalScl = 0.0 - restore_salinity = .false. - if (present(restore_salt)) restore_salinity = restore_salt - restore_sst = .false. - if (present(restore_temp)) restore_sst = restore_temp - ! allocation and initialization if this is the first time that this ! flux type has been used. if (fluxes%dt_buoy_accum < 0) then @@ -305,7 +290,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, & fluxes%ustar_tidal(i,j) = CS%ustar_tidal(i,j) enddo ; enddo - if (restore_temp) call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed) + if (CS%restore_temp) call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed) fluxes%dt_buoy_accum = 0.0 endif ! endif for allocation and initialization @@ -324,7 +309,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, & ! ocean model, rather than using haloless arrays, in which case the last line ! would be: ( (/isd,is,ie,ied/), (/jsd,js,je,jed/)) - if (CS%allow_flux_adjustments) then fluxes%heat_added(:,:)=0.0 fluxes%salt_flux_added(:,:)=0.0 @@ -344,7 +328,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, & enddo ; enddo ! Salinity restoring logic - if (restore_salinity) then + if (CS%restore_salt) then call time_interp_external(CS%id_srestore,Time,data_restore) ! open_ocn_mask indicates where to restore salinity (1 means restore, 0 does not) open_ocn_mask(:,:) = 1.0 @@ -397,7 +381,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, & endif ! SST restoring logic - if (restore_sst) then + if (CS%restore_temp) then call time_interp_external(CS%id_trestore,Time,data_restore) do j=js,je ; do i=is,ie delta_sst = data_restore(i,j)- sfc_state%SST(i,j) @@ -545,6 +529,16 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, & endif + ! Set the wind stresses and ustar. + if (associated(fluxes%ustar) .and. associated(fluxes%ustar_gustless)) then + call extract_IOB_stresses(IOB, index_bounds, Time, G, CS, ustar=fluxes%ustar, & + gustless_ustar=fluxes%ustar_gustless) + elseif (associated(fluxes%ustar)) then + call extract_IOB_stresses(IOB, index_bounds, Time, G, CS, ustar=fluxes%ustar) + elseif (associated(fluxes%ustar_gustless)) then + call extract_IOB_stresses(IOB, index_bounds, Time, G, CS, gustless_ustar=fluxes%ustar_gustless) + endif + if (coupler_type_initialized(fluxes%tr_fluxes) .and. & coupler_type_initialized(IOB%fluxes)) & call coupler_type_copy_data(IOB%fluxes, fluxes%tr_fluxes) @@ -564,7 +558,7 @@ end subroutine convert_IOB_to_fluxes !> This subroutine translates the Ice_ocean_boundary_type into a MOM !! mechanical forcing type, including changes of units, sign conventions, !! and putting the fields into arrays with MOM-standard halos. -subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, CS) +subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, CS, dt_forcing, reset_avg) type(ice_ocean_boundary_type), & target, intent(in) :: IOB !< An ice-ocean boundary type with fluxes to drive !! the ocean in a coupled model @@ -575,27 +569,24 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, CS) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure type(surface_forcing_CS),pointer :: CS !< A pointer to the control structure returned by a !! previous call to surface_forcing_init. + real, optional, intent(in) :: dt_forcing !< A time interval over which to apply the + !! current value of ustar as a weighted running + !! average, in s, or if 0 do not average ustar. + !! Missing is equivalent to 0. + logical, optional, intent(in) :: reset_avg !< If true, reset the time average. - - real, dimension(SZIB_(G),SZJB_(G)) :: & - taux_at_q, & ! Zonal wind stresses at q points (Pa) - tauy_at_q ! Meridional wind stresses at q points (Pa) - + ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: & - rigidity_at_h, & ! Ice rigidity at tracer points (m3 s-1) - taux_at_h, & ! Zonal wind stresses at h points (Pa) - tauy_at_h ! Meridional wind stresses at h points (Pa) + rigidity_at_h, & ! Ice rigidity at tracer points (m3 s-1) + net_mass_src, & ! A temporary of net mass sources, in kg m-2 s-1. + ustar_tmp ! A temporary array of ustar values, in m s-1. - real :: gustiness ! unresolved gustiness that contributes to ustar (Pa) - real :: Irho0 ! inverse of the mean density in (m^3/kg) - real :: taux2, tauy2 ! squared wind stresses (Pa^2) - real :: tau_mag ! magnitude of the wind stress (Pa) real :: I_GEarth ! 1.0 / G%G_Earth (s^2/m) real :: Kv_rho_ice ! (CS%kv_sea_ice / CS%density_sea_ice) ( m^5/(s*kg) ) real :: mass_ice ! mass of sea ice at a face (kg/m^2) real :: mass_eff ! effective mass of sea ice for rigidity (kg/m^2) + real :: wt1, wt2 ! Relative weights of previous and current values of ustar, ND. - integer :: wind_stagger ! AGRID, BGRID_NE, or CGRID_NE (integers from MOM_domains) integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isr, ier, jsr, jer integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd @@ -611,8 +602,6 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, CS) isr = is-isd+1 ; ier = ie-isd+1 ; jsr = js-jsd+1 ; jer = je-jsd+1 i0 = is - isc_bnd ; j0 = js - jsc_bnd - Irho0 = 1.0/CS%Rho0 - ! allocation and initialization if this is the first time that this ! mechanical forcing type has been used. if (.not.forces%initialized) then @@ -638,6 +627,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, CS) if ( (associated(IOB%area_berg) .and. (.not. associated(forces%area_berg))) .or. & (associated(IOB%mass_berg) .and. (.not. associated(forces%mass_berg))) ) & call allocate_mech_forcing(G, forces, iceberg=.true.) + if (associated(IOB%ice_rigidity)) then rigidity_at_h(:,:) = 0.0 call safe_alloc_ptr(forces%rigidity_ice_u,IsdB,IedB,jsd,jed) @@ -648,6 +638,23 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, CS) if (associated(forces%rigidity_ice_u)) forces%rigidity_ice_u(:,:) = 0.0 if (associated(forces%rigidity_ice_v)) forces%rigidity_ice_v(:,:) = 0.0 + ! Set the weights for forcing fields that use running time averages. + if (present(reset_avg)) then ; if (reset_avg) forces%dt_force_accum = 0.0 ; endif + wt1 = 0.0 ; wt2 = 1.0 + if (present(dt_forcing)) then + if ((forces%dt_force_accum > 0.0) .and. (dt_forcing > 0.0)) then + wt1 = forces%dt_force_accum / (forces%dt_force_accum + dt_forcing) + wt2 = 1.0 - wt1 + endif + if (dt_forcing > 0.0) then + forces%dt_force_accum = max(forces%dt_force_accum, 0.0) + dt_forcing + else + forces%dt_force_accum = 0.0 ! Reset the averaging time interval. + endif + else + forces%dt_force_accum = 0.0 ! Reset the averaging time interval. + endif + ! applied surface pressure from atmosphere and cryosphere if (associated(IOB%p)) then if (CS%max_p_surf >= 0.0) then @@ -669,136 +676,62 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, CS) endif forces%accumulate_p_surf = .true. ! Multiple components may contribute to surface pressure. - wind_stagger = CS%wind_stagger - if ((IOB%wind_stagger == AGRID) .or. (IOB%wind_stagger == BGRID_NE) .or. & - (IOB%wind_stagger == CGRID_NE)) wind_stagger = IOB%wind_stagger - if (wind_stagger == BGRID_NE) then - ! This is necessary to fill in the halo points. - taux_at_q(:,:) = 0.0 ; tauy_at_q(:,:) = 0.0 - endif - if (wind_stagger == AGRID) then - ! This is necessary to fill in the halo points. - taux_at_h(:,:) = 0.0 ; tauy_at_h(:,:) = 0.0 - endif - - ! obtain fluxes from IOB; note the staggering of indices - do j=js,je ; do i=is,ie - if (associated(IOB%area_berg)) & - forces%area_berg(i,j) = IOB%area_berg(i-i0,j-j0) * G%mask2dT(i,j) - - if (associated(IOB%mass_berg)) & - forces%mass_berg(i,j) = IOB%mass_berg(i-i0,j-j0) * G%mask2dT(i,j) - - if (associated(IOB%ice_rigidity)) & - rigidity_at_h(i,j) = IOB%ice_rigidity(i-i0,j-j0) * G%mask2dT(i,j) - - if (wind_stagger == BGRID_NE) then - if (associated(IOB%u_flux)) taux_at_q(I,J) = IOB%u_flux(i-i0,j-j0) * CS%wind_stress_multiplier - if (associated(IOB%v_flux)) tauy_at_q(I,J) = IOB%v_flux(i-i0,j-j0) * CS%wind_stress_multiplier - elseif (wind_stagger == AGRID) then - if (associated(IOB%u_flux)) taux_at_h(i,j) = IOB%u_flux(i-i0,j-j0) * CS%wind_stress_multiplier - if (associated(IOB%v_flux)) tauy_at_h(i,j) = IOB%v_flux(i-i0,j-j0) * CS%wind_stress_multiplier - else ! C-grid wind stresses. - if (associated(IOB%u_flux)) forces%taux(I,j) = IOB%u_flux(i-i0,j-j0) * CS%wind_stress_multiplier - if (associated(IOB%v_flux)) forces%tauy(i,J) = IOB%v_flux(i-i0,j-j0) * CS%wind_stress_multiplier - endif - - enddo ; enddo - - ! surface momentum stress related fields as function of staggering - if (wind_stagger == BGRID_NE) then - if (G%symmetric) & - call fill_symmetric_edges(taux_at_q, tauy_at_q, G%Domain, stagger=BGRID_NE) - call pass_vector(taux_at_q, tauy_at_q, G%Domain, stagger=BGRID_NE, halo=1) - - do j=js,je ; do I=Isq,Ieq - forces%taux(I,j) = 0.0 - if ((G%mask2dBu(I,J) + G%mask2dBu(I,J-1)) > 0) & - forces%taux(I,j) = (G%mask2dBu(I,J)*taux_at_q(I,J) + & - G%mask2dBu(I,J-1)*taux_at_q(I,J-1)) / & - (G%mask2dBu(I,J) + G%mask2dBu(I,J-1)) - enddo ; enddo - - do J=Jsq,Jeq ; do i=is,ie - forces%tauy(i,J) = 0.0 - if ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J)) > 0) & - forces%tauy(i,J) = (G%mask2dBu(I,J)*tauy_at_q(I,J) + & - G%mask2dBu(I-1,J)*tauy_at_q(I-1,J)) / & - (G%mask2dBu(I,J) + G%mask2dBu(I-1,J)) - enddo ; enddo - - ! ustar is required for the bulk mixed layer formulation. The background value - ! of 0.02 Pa is a relatively small value intended to give reasonable behavior - ! in regions of very weak winds. - + ! Set the wind stresses and ustar. + if (wt1 <= 0.0) then + call extract_IOB_stresses(IOB, index_bounds, Time, G, CS, taux=forces%taux, tauy=forces%tauy, & + ustar=forces%ustar, tau_halo=1) + else + call extract_IOB_stresses(IOB, index_bounds, Time, G, CS, taux=forces%taux, tauy=forces%tauy, & + ustar=ustar_tmp, tau_halo=1) do j=js,je ; do i=is,ie - tau_mag = 0.0 ; gustiness = CS%gust_const - if (((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + & - (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0) then - tau_mag = sqrt(((G%mask2dBu(I,J)*(taux_at_q(I,J)**2 + tauy_at_q(I,J)**2) + & - G%mask2dBu(I-1,J-1)*(taux_at_q(I-1,J-1)**2 + tauy_at_q(I-1,J-1)**2)) + & - (G%mask2dBu(I,J-1)*(taux_at_q(I,J-1)**2 + tauy_at_q(I,J-1)**2) + & - G%mask2dBu(I-1,J)*(taux_at_q(I-1,J)**2 + tauy_at_q(I-1,J)**2)) ) / & - ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) ) - if (CS%read_gust_2d) gustiness = CS%gust(i,j) - endif - forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0*tau_mag) - enddo ; enddo - - elseif (wind_stagger == AGRID) then - call pass_vector(taux_at_h, tauy_at_h, G%Domain, To_All+Omit_Corners, & - stagger=AGRID, halo=1) - - do j=js,je ; do I=Isq,Ieq - forces%taux(I,j) = 0.0 - if ((G%mask2dT(i,j) + G%mask2dT(i+1,j)) > 0) & - forces%taux(I,j) = (G%mask2dT(i,j)*taux_at_h(i,j) + & - G%mask2dT(i+1,j)*taux_at_h(i+1,j)) / & - (G%mask2dT(i,j) + G%mask2dT(i+1,j)) + forces%ustar(i,j) = wt1*forces%ustar(i,j) + wt2*ustar_tmp(i,j) enddo ; enddo + endif - do J=Jsq,Jeq ; do i=is,ie - forces%tauy(i,J) = 0.0 - if ((G%mask2dT(i,j) + G%mask2dT(i,j+1)) > 0) & - forces%tauy(i,J) = (G%mask2dT(i,j)*tauy_at_h(i,j) + & - G%mask2dT(i,J+1)*tauy_at_h(i,j+1)) / & - (G%mask2dT(i,j) + G%mask2dT(i,j+1)) - enddo ; enddo + ! Find the net mass source in the input forcing without other adjustments. + if (CS%approx_net_mass_src .and. associated(forces%net_mass_src)) then + net_mass_src(:,:) = 0.0 + i0 = is - isc_bnd ; j0 = js - jsc_bnd + do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then + if (associated(IOB%lprec)) & + net_mass_src(i,j) = net_mass_src(i,j) + IOB%lprec(i-i0,j-j0) + if (associated(IOB%fprec)) & + net_mass_src(i,j) = net_mass_src(i,j) + IOB%fprec(i-i0,j-j0) + if (associated(IOB%runoff)) & + net_mass_src(i,j) = net_mass_src(i,j) + IOB%runoff(i-i0,j-j0) + if (associated(IOB%calving)) & + net_mass_src(i,j) = net_mass_src(i,j) + IOB%calving(i-i0,j-j0) + if (associated(IOB%q_flux)) & + net_mass_src(i,j) = net_mass_src(i,j) - IOB%q_flux(i-i0,j-j0) + endif ; enddo ; enddo + if (wt1 <= 0.0) then + do j=js,je ; do i=is,ie + forces%net_mass_src(i,j) = wt2*net_mass_src(i,j) + enddo ; enddo + else + do j=js,je ; do i=is,ie + forces%net_mass_src(i,j) = wt1*forces%net_mass_src(i,j) + wt2*net_mass_src(i,j) + enddo ; enddo + endif + forces%net_mass_src_set = .true. + else + forces%net_mass_src_set = .false. + endif - do j=js,je ; do i=is,ie - gustiness = CS%gust_const - if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0)) gustiness = CS%gust(i,j) - forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * G%mask2dT(i,j) * & - sqrt(taux_at_h(i,j)**2 + tauy_at_h(i,j)**2)) - enddo ; enddo + ! Obtain optional ice-berg related fluxes from the IOB type: + if (associated(IOB%area_berg)) then ; do j=js,je ; do i=is,ie + forces%area_berg(i,j) = IOB%area_berg(i-i0,j-j0) * G%mask2dT(i,j) + enddo ; enddo ; endif - else ! C-grid wind stresses. - if (G%symmetric) & - call fill_symmetric_edges(forces%taux, forces%tauy, G%Domain) - call pass_vector(forces%taux, forces%tauy, G%Domain, halo=1) + if (associated(IOB%mass_berg)) then ; do j=js,je ; do i=is,ie + forces%mass_berg(i,j) = IOB%mass_berg(i-i0,j-j0) * G%mask2dT(i,j) + enddo ; enddo ; endif + ! Obtain sea ice related dynamic fields + if (associated(IOB%ice_rigidity)) then do j=js,je ; do i=is,ie - taux2 = 0.0 - if ((G%mask2dCu(I-1,j) + G%mask2dCu(I,j)) > 0) & - taux2 = (G%mask2dCu(I-1,j)*forces%taux(I-1,j)**2 + & - G%mask2dCu(I,j)*forces%taux(I,j)**2) / (G%mask2dCu(I-1,j) + G%mask2dCu(I,j)) - - tauy2 = 0.0 - if ((G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) > 0) & - tauy2 = (G%mask2dCv(i,J-1)*forces%tauy(i,J-1)**2 + & - G%mask2dCv(i,J)*forces%tauy(i,J)**2) / (G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) - - if (CS%read_gust_2d) then - forces%ustar(i,j) = sqrt(CS%gust(i,j)*Irho0 + Irho0*sqrt(taux2 + tauy2)) - else - forces%ustar(i,j) = sqrt(CS%gust_const*Irho0 + Irho0*sqrt(taux2 + tauy2)) - endif + rigidity_at_h(i,j) = IOB%ice_rigidity(i-i0,j-j0) * G%mask2dT(i,j) enddo ; enddo - - endif ! endif for wind related fields - - ! sea ice related dynamic fields - if (associated(IOB%ice_rigidity)) then call pass_var(rigidity_at_h, G%Domain, halo=1) do I=is-1,ie ; do j=js,je forces%rigidity_ice_u(I,j) = forces%rigidity_ice_u(I,j) + & @@ -845,6 +778,226 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, CS) call cpu_clock_end(id_clock_forcing) end subroutine convert_IOB_to_forces + +!> This subroutine extracts the wind stresses and related fields like ustar from an +!! Ice_ocean_boundary_type into optional argument arrays, including changes of units, sign +!! conventions, and putting the fields into arrays with MOM-standard sized halos. +subroutine extract_IOB_stresses(IOB, index_bounds, Time, G, CS, taux, tauy, ustar, & + gustless_ustar, tau_halo) + type(ice_ocean_boundary_type), & + target, intent(in) :: IOB !< An ice-ocean boundary type with fluxes to drive + !! the ocean in a coupled model + integer, dimension(4), intent(in) :: index_bounds !< The i- and j- size of the arrays in IOB. + type(time_type), intent(in) :: Time !< The time of the fluxes, used for interpolating the + !! salinity to the right time, when it is being restored. + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + type(surface_forcing_CS),pointer :: CS !< A pointer to the control structure returned by a + !! previous call to surface_forcing_init. + real, dimension(SZIB_(G),SZJ_(G)), & + optional, intent(inout) :: taux !< The zonal wind stresses on a C-grid, in Pa. + real, dimension(SZI_(G),SZJB_(G)), & + optional, intent(inout) :: tauy !< The meridional wind stresses on a C-grid, in Pa. + real, dimension(SZI_(G),SZJ_(G)), & + optional, intent(inout) :: ustar !< The surface friction velocity, in m s-1. + real, dimension(SZI_(G),SZJ_(G)), & + optional, intent(out) :: gustless_ustar !< The surface friction velocity without + !! any contributions from gustiness, in m s-1. + integer, optional, intent(in) :: tau_halo !< The halo size of wind stresses to set, 0 by default. + + ! Local variables + real, dimension(SZI_(G),SZJ_(G)) :: taux_in_A ! Zonal wind stresses (in Pa) at h points + real, dimension(SZI_(G),SZJ_(G)) :: tauy_in_A ! Meridional wind stresses (in Pa) at h points + real, dimension(SZIB_(G),SZJ_(G)) :: taux_in_C ! Zonal wind stresses (in Pa) at u points + real, dimension(SZI_(G),SZJB_(G)) :: tauy_in_C ! Meridional wind stresses (in Pa) at v points + real, dimension(SZIB_(G),SZJB_(G)) :: taux_in_B ! Zonal wind stresses (in Pa) at q points + real, dimension(SZIB_(G),SZJB_(G)) :: tauy_in_B ! Meridional wind stresses (in Pa) at q points + + real :: gustiness ! unresolved gustiness that contributes to ustar (Pa) + real :: Irho0 ! inverse of the mean density in (m^3/kg) + real :: taux2, tauy2 ! squared wind stresses (Pa^2) + real :: tau_mag ! magnitude of the wind stress (Pa) + + logical :: do_ustar, do_gustless + integer :: wind_stagger ! AGRID, BGRID_NE, or CGRID_NE (integers from MOM_domains) + integer :: i, j, is, ie, js, je, ish, ieh, jsh, jeh, Isqh, Ieqh, Jsqh, Jeqh, i0, j0, halo + + halo = 0 ; if (present(tau_halo)) halo = tau_halo + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + ish = G%isc-halo ; ieh = G%iec+halo ; jsh = G%jsc-halo ; jeh = G%jec+halo + Isqh = G%IscB-halo ; Ieqh = G%IecB+halo ; Jsqh = G%JscB-halo ; Jeqh = G%JecB+halo + i0 = is - index_bounds(1) ; j0 = js - index_bounds(3) + + Irho0 = 1.0/CS%Rho0 + + do_ustar = present(ustar) ; do_gustless = present(gustless_ustar) + + wind_stagger = CS%wind_stagger + if ((IOB%wind_stagger == AGRID) .or. (IOB%wind_stagger == BGRID_NE) .or. & + (IOB%wind_stagger == CGRID_NE)) wind_stagger = IOB%wind_stagger + + if (associated(IOB%u_flux).neqv.associated(IOB%v_flux)) call MOM_error(FATAL,"extract_IOB_stresses: "//& + "associated(IOB%u_flux) /= associated(IOB%v_flux !!!") + if (present(taux).neqv.present(tauy)) call MOM_error(FATAL,"extract_IOB_stresses: "//& + "present(taux) /= present(tauy) !!!") + + ! Set surface momentum stress related fields as a function of staggering. + if (present(taux) .or. present(tauy) .or. & + ((do_ustar.or.do_gustless) .and. .not.associated(IOB%stress_mag)) ) then + + if (wind_stagger == BGRID_NE) then + taux_in_B(:,:) = 0.0 ; tauy_in_B(:,:) = 0.0 + if (associated(IOB%u_flux).and.associated(IOB%v_flux)) then + do J=js,je ; do I=is,ie + taux_in_B(I,J) = IOB%u_flux(i-i0,j-j0) * CS%wind_stress_multiplier + tauy_in_B(I,J) = IOB%v_flux(i-i0,j-j0) * CS%wind_stress_multiplier + enddo ; enddo + endif + + if (G%symmetric) call fill_symmetric_edges(taux_in_B, tauy_in_B, G%Domain, stagger=BGRID_NE) + call pass_vector(taux_in_B, tauy_in_B, G%Domain, stagger=BGRID_NE, halo=max(1,halo)) + + if (present(taux).and.present(tauy)) then + do j=jsh,jeh ; do I=Isqh,Ieqh + taux(I,j) = 0.0 + if ((G%mask2dBu(I,J) + G%mask2dBu(I,J-1)) > 0) & + taux(I,j) = (G%mask2dBu(I,J)*taux_in_B(I,J) + G%mask2dBu(I,J-1)*taux_in_B(I,J-1)) / & + (G%mask2dBu(I,J) + G%mask2dBu(I,J-1)) + enddo ; enddo + do J=Jsqh,Jeqh ; do i=ish,ieh + tauy(i,J) = 0.0 + if ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J)) > 0) & + tauy(i,J) = (G%mask2dBu(I,J)*tauy_in_B(I,J) + G%mask2dBu(I-1,J)*tauy_in_B(I-1,J)) / & + (G%mask2dBu(I,J) + G%mask2dBu(I-1,J)) + enddo ; enddo + endif + elseif (wind_stagger == AGRID) then + taux_in_A(:,:) = 0.0 ; tauy_in_A(:,:) = 0.0 + if (associated(IOB%u_flux).and.associated(IOB%v_flux)) then + do j=js,je ; do i=is,ie + taux_in_A(i,j) = IOB%u_flux(i-i0,j-j0) * CS%wind_stress_multiplier + tauy_in_A(i,j) = IOB%v_flux(i-i0,j-j0) * CS%wind_stress_multiplier + enddo ; enddo + endif + + if (halo == 0) then + call pass_vector(taux_in_A, tauy_in_A, G%Domain, To_All+Omit_Corners, stagger=AGRID, halo=1) + else + call pass_vector(taux_in_A, tauy_in_A, G%Domain, stagger=AGRID, halo=max(1,halo)) + endif + + if (present(taux)) then ; do j=jsh,jeh ; do I=Isqh,Ieqh + taux(I,j) = 0.0 + if ((G%mask2dT(i,j) + G%mask2dT(i+1,j)) > 0) & + taux(I,j) = (G%mask2dT(i,j)*taux_in_A(i,j) + G%mask2dT(i+1,j)*taux_in_A(i+1,j)) / & + (G%mask2dT(i,j) + G%mask2dT(i+1,j)) + enddo ; enddo ; endif + + if (present(tauy)) then ; do J=Jsqh,Jeqh ; do i=ish,ieh + tauy(i,J) = 0.0 + if ((G%mask2dT(i,j) + G%mask2dT(i,j+1)) > 0) & + tauy(i,J) = (G%mask2dT(i,j)*tauy_in_A(i,j) + G%mask2dT(i,J+1)*tauy_in_A(i,j+1)) / & + (G%mask2dT(i,j) + G%mask2dT(i,j+1)) + enddo ; enddo ; endif + + else ! C-grid wind stresses. + taux_in_C(:,:) = 0.0 ; tauy_in_C(:,:) = 0.0 + if (associated(IOB%u_flux).and.associated(IOB%v_flux)) then + do j=js,je ; do i=is,ie + taux_in_C(I,j) = IOB%u_flux(i-i0,j-j0) * CS%wind_stress_multiplier + tauy_in_C(i,J) = IOB%v_flux(i-i0,j-j0) * CS%wind_stress_multiplier + enddo ; enddo + endif + + if (G%symmetric) call fill_symmetric_edges(taux_in_C, tauy_in_C, G%Domain) + call pass_vector(taux_in_C, tauy_in_C, G%Domain, halo=max(1,halo)) + + if (present(taux).and.present(tauy)) then + do j=jsh,jeh ; do I=Isqh,Ieqh + taux(I,j) = G%mask2dCu(I,j)*taux_in_C(I,j) + enddo ; enddo + do J=Jsqh,Jeqh ; do i=ish,ieh + tauy(i,J) = G%mask2dCv(i,J)*tauy_in_C(i,J) + enddo ; enddo + endif + endif ! endif for extracting wind stress fields with various staggerings + endif + + if (do_ustar .or. do_gustless) then + ! Set surface friction velocity directly or as a function of staggering. + ! ustar is required for the bulk mixed layer formulation and other turbulent mixing + ! parametizations. The background gustiness (for example with a relatively small value + ! of 0.02 Pa) is intended to give reasonable behavior in regions of very weak winds. + if (associated(IOB%stress_mag)) then + if (do_ustar) then ; do j=js,je ; do i=is,ie + gustiness = CS%gust_const + !### SIMPLIFY THE TREATMENT OF GUSTINESS! + if (CS%read_gust_2d) then + if ((wind_stagger == CGRID_NE) .or. & + ((wind_stagger == AGRID) .and. (G%mask2dT(i,j) > 0)) .or. & + ((wind_stagger == BGRID_NE) .and. & + (((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + & + (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0)) ) & + gustiness = CS%gust(i,j) + endif + ustar(i,j) = sqrt(gustiness*Irho0 + Irho0*IOB%stress_mag(i-i0,j-j0)) + enddo ; enddo ; endif + if (do_gustless) then ; do j=js,je ; do i=is,ie + gustless_ustar(i,j) = sqrt(IOB%stress_mag(i-i0,j-j0) / CS%Rho0) +!### Change to: +! gustless_ustar(i,j) = sqrt(Irho0 * IOB%stress_mag(i-i0,j-j0)) + enddo ; enddo ; endif + elseif (wind_stagger == BGRID_NE) then + do j=js,je ; do i=is,ie + tau_mag = 0.0 ; gustiness = CS%gust_const + if (((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + & + (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0) then + tau_mag = sqrt(((G%mask2dBu(I,J)*(taux_in_B(I,J)**2 + tauy_in_B(I,J)**2) + & + G%mask2dBu(I-1,J-1)*(taux_in_B(I-1,J-1)**2 + tauy_in_B(I-1,J-1)**2)) + & + (G%mask2dBu(I,J-1)*(taux_in_B(I,J-1)**2 + tauy_in_B(I,J-1)**2) + & + G%mask2dBu(I-1,J)*(taux_in_B(I-1,J)**2 + tauy_in_B(I-1,J)**2)) ) / & + ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) ) + if (CS%read_gust_2d) gustiness = CS%gust(i,j) + endif + if (do_ustar) ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * tau_mag) + if (do_gustless) gustless_ustar(i,j) = sqrt(tau_mag / CS%Rho0) +!### Change to: +! if (do_gustless) gustless_ustar(i,j) = sqrt(Irho0 * tau_mag) + enddo ; enddo + elseif (wind_stagger == AGRID) then + do j=js,je ; do i=is,ie + tau_mag = G%mask2dT(i,j) * sqrt(taux_in_A(i,j)**2 + tauy_in_A(i,j)**2) + gustiness = CS%gust_const + if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0)) gustiness = CS%gust(i,j) + if (do_ustar) ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * tau_mag) + if (do_gustless) gustless_ustar(i,j) = sqrt(tau_mag / CS%Rho0) +!### Change to: +! if (do_gustless) gustless_ustar(i,j) = sqrt(Irho0 * tau_mag) + enddo ; enddo + else ! C-grid wind stresses. + do j=js,je ; do i=is,ie + taux2 = 0.0 ; tauy2 = 0.0 + if ((G%mask2dCu(I-1,j) + G%mask2dCu(I,j)) > 0) & + taux2 = (G%mask2dCu(I-1,j)*taux_in_C(I-1,j)**2 + & + G%mask2dCu(I,j)*taux_in_C(I,j)**2) / (G%mask2dCu(I-1,j) + G%mask2dCu(I,j)) + if ((G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) > 0) & + tauy2 = (G%mask2dCv(i,J-1)*tauy_in_C(i,J-1)**2 + & + G%mask2dCv(i,J)*tauy_in_C(i,J)**2) / (G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) + tau_mag = sqrt(taux2 + tauy2) + + gustiness = CS%gust_const + if (CS%read_gust_2d) gustiness = CS%gust(i,j) + + if (do_ustar) ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * tau_mag) + if (do_gustless) gustless_ustar(i,j) = sqrt(tau_mag / CS%Rho0) +!### Change to: +! if (do_gustless) gustless_ustar(i,j) = sqrt(Irho0 * tau_mag) + enddo ; enddo + endif ! endif for wind friction velocity fields + endif + +end subroutine extract_IOB_stresses + + !> Adds thermodynamic flux adjustments obtained via data_override !! Component name is 'OCN' !! Available adjustments are: @@ -969,7 +1122,7 @@ subroutine forcing_save_restart(CS, G, Time, directory, time_stamped, & end subroutine forcing_save_restart !> Initialize the surface forcing, including setting parameters and allocating permanent memory. -subroutine surface_forcing_init(Time, G, param_file, diag, CS, restore_salt, restore_temp) +subroutine surface_forcing_init(Time, G, param_file, diag, CS) type(time_type), intent(in) :: Time !< The current model time type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters @@ -977,10 +1130,6 @@ subroutine surface_forcing_init(Time, G, param_file, diag, CS, restore_salt, res !! diagnostic output type(surface_forcing_CS), pointer :: CS !< A pointer that is set to point to the control !! structure for this module - logical, optional, intent(in) :: restore_salt !< If present and true surface salinity - !! restoring will be applied in this model. - logical, optional, intent(in) :: restore_temp !< If present and true surface temperature - !! restoring will be applied in this model. ! Local variables real :: utide ! The RMS tidal velocity, in m s-1. @@ -1039,11 +1188,19 @@ subroutine surface_forcing_init(Time, G, param_file, diag, CS, restore_salt, res "the ice-ocean heat fluxes are treated explicitly. No \n"//& "limit is applied if a negative value is used.", units="Pa", & default=-1.0) + call get_param(param_file, mdl, "RESTORE_SALINITY", CS%restore_salt, & + "If true, the coupled driver will add a globally-balanced \n"//& + "fresh-water flux that drives sea-surface salinity \n"//& + "toward specified values.", default=.false.) + call get_param(param_file, mdl, "RESTORE_TEMPERATURE", CS%restore_temp, & + "If true, the coupled driver will add a \n"//& + "heat flux that drives sea-surface temperauture \n"//& + "toward specified values.", default=.false.) call get_param(param_file, mdl, "ADJUST_NET_SRESTORE_TO_ZERO", & CS%adjust_net_srestore_to_zero, & "If true, adjusts the salinity restoring seen to zero\n"//& "whether restoring is via a salt flux or virtual precip.",& - default=restore_salt) + default=CS%restore_salt) call get_param(param_file, mdl, "ADJUST_NET_SRESTORE_BY_SCALING", & CS%adjust_net_srestore_by_scaling, & "If true, adjustments to salt restoring to achieve zero net are\n"//& @@ -1073,6 +1230,11 @@ subroutine surface_forcing_init(Time, G, param_file, diag, CS, restore_salt, res "correction for the atmospheric (and sea-ice) pressure \n"//& "limited by max_p_surf instead of the full atmospheric \n"//& "pressure.", default=.true.) + call get_param(param_file, mdl, "APPROX_NET_MASS_SRC", CS%approx_net_mass_src, & + "If true, use the net mass sources from the ice-ocean \n"//& + "boundary type without any further adjustments to drive \n"//& + "the ocean dynamics. The actual net mass source may differ \n"//& + "due to internal corrections.", default=.false.) call get_param(param_file, mdl, "WIND_STAGGER", stagger, & "A case-insensitive character string to indicate the \n"//& @@ -1088,7 +1250,7 @@ subroutine surface_forcing_init(Time, G, param_file, diag, CS, restore_salt, res "coupler. This is used for testing and should be =1.0 for any\n"//& "production runs.", default=1.0) - if (restore_salt) then + if (CS%restore_salt) then call get_param(param_file, mdl, "FLUXCONST", CS%Flux_const, & "The constant that relates the restoring surface fluxes \n"//& "to the relative surface anomalies (akin to a piston \n"//& @@ -1136,7 +1298,7 @@ subroutine surface_forcing_init(Time, G, param_file, diag, CS, restore_salt, res "a mask for SSS restoring.", default=.false.) endif - if (restore_temp) then + if (CS%restore_temp) then call get_param(param_file, mdl, "FLUXCONST", CS%Flux_const, & "The constant that relates the restoring surface fluxes \n"//& "to the relative surface anomalies (akin to a piston \n"//& @@ -1255,7 +1417,7 @@ subroutine surface_forcing_init(Time, G, param_file, diag, CS, restore_salt, res call data_override_init(Ocean_domain_in=G%Domain%mpp_domain) endif - if (present(restore_salt)) then ; if (restore_salt) then + if (CS%restore_salt) then salt_file = trim(CS%inputdir) // trim(CS%salt_restore_file) CS%id_srestore = init_external_field(salt_file, CS%salt_restore_var_name, domain=G%Domain%mpp_domain) call safe_alloc_ptr(CS%srestore_mask,isd,ied,jsd,jed); CS%srestore_mask(:,:) = 1.0 @@ -1263,9 +1425,9 @@ subroutine surface_forcing_init(Time, G, param_file, diag, CS, restore_salt, res flnam = trim(CS%inputdir) // 'salt_restore_mask.nc' call MOM_read_data(flnam,'mask', CS%srestore_mask, G%domain, timelevel=1) endif - endif ; endif + endif - if (present(restore_temp)) then ; if (restore_temp) then + if (CS%restore_temp) then temp_file = trim(CS%inputdir) // trim(CS%temp_restore_file) CS%id_trestore = init_external_field(temp_file, CS%temp_restore_var_name, domain=G%Domain%mpp_domain) call safe_alloc_ptr(CS%trestore_mask,isd,ied,jsd,jed); CS%trestore_mask(:,:) = 1.0 @@ -1273,7 +1435,7 @@ subroutine surface_forcing_init(Time, G, param_file, diag, CS, restore_salt, res flnam = trim(CS%inputdir) // 'temp_restore_mask.nc' call MOM_read_data(flnam, 'mask', CS%trestore_mask, G%domain, timelevel=1) endif - endif ; endif + endif ! Set up any restart fields associated with the forcing. call restart_init(param_file, CS%restart_CSp, "MOM_forcing.res") diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index a09a5bfe29..70437d0e4c 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -20,14 +20,12 @@ module ocean_model_mod use MOM_diag_mediator, only : diag_mediator_close_registration, diag_mediator_end use MOM_domains, only : pass_var, pass_vector, AGRID, BGRID_NE, CGRID_NE use MOM_domains, only : TO_ALL, Omit_Corners -use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe +use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave use MOM_file_parser, only : get_param, log_version, close_param_file, param_file_type -use MOM_forcing_type, only : allocate_forcing_type -use MOM_forcing_type, only : forcing, mech_forcing -use MOM_forcing_type, only : forcing_accumulate, copy_common_forcing_fields -use MOM_forcing_type, only : copy_back_forcing_fields, set_net_mass_forcing -use MOM_forcing_type, only : set_derived_forcing_fields +use MOM_forcing_type, only : forcing, mech_forcing, allocate_forcing_type +use MOM_forcing_type, only : fluxes_accumulate, get_net_mass_forcing +use MOM_forcing_type, only : copy_back_forcing_fields use MOM_forcing_type, only : forcing_diagnostics, mech_forcing_diags use MOM_get_input, only : Get_MOM_Input, directories use MOM_grid, only : ocean_grid_type @@ -39,10 +37,10 @@ module ocean_model_mod use MOM_surface_forcing, only : convert_IOB_to_forces, ice_ocn_bnd_type_chksum use MOM_surface_forcing, only : ice_ocean_boundary_type, surface_forcing_CS use MOM_surface_forcing, only : forcing_save_restart -use MOM_time_manager, only : time_type, get_time, set_time, operator(>) -use MOM_time_manager, only : operator(+), operator(-), operator(*), operator(/) -use MOM_time_manager, only : operator(/=), operator(<=), operator(>=) -use MOM_time_manager, only : operator(<), real_to_time_type, time_type_to_real +use MOM_time_manager, only : time_type, operator(>), operator(+), operator(-) +use MOM_time_manager, only : operator(*), operator(/), operator(/=) +use MOM_time_manager, only : operator(<=), operator(>=), operator(<) +use MOM_time_manager, only : real_to_time, time_type_to_real use MOM_tracer_flow_control, only : call_tracer_register, tracer_flow_control_init use MOM_tracer_flow_control, only : call_tracer_flux_init use MOM_variables, only : surface @@ -138,6 +136,8 @@ module ocean_model_mod ! This type is private, and can therefore vary between different ocean models. logical :: is_ocean_PE = .false. !< True if this is an ocean PE. type(time_type) :: Time !< The ocean model's time and master clock. + type(time_type) :: Time_dyn !< The ocean model's time for the dynamics. Time and Time_dyn + !! should be the same after a full time step. integer :: Restart_control !< An integer that is bit-tested to determine whether !! incremental restart files are saved and whether they !! have a time stamped name. +1 (bit 0) for generic @@ -145,16 +145,13 @@ module ocean_model_mod !! restart file is saved at the end of a run segment !! unless Restart_control is negative. - integer :: nstep = 0 !< The number of calls to update_ocean. + integer :: nstep = 0 !< The number of calls to update_ocean that update the dynamics. + integer :: nstep_thermo = 0 !< The number of calls to update_ocean that update the thermodynamics. logical :: use_ice_shelf !< If true, the ice shelf model is enabled. logical :: use_waves !< If true use wave coupling. logical :: icebergs_alter_ocean !< If true, the icebergs can change ocean the !! ocean dynamics and forcing fluxes. - logical :: restore_salinity !< If true, the coupled MOM driver adds a term to - !! restore salinity to a specified value. - logical :: restore_temp !< If true, the coupled MOM driver adds a term to - !! restore sst to a specified value. real :: press_to_z !< A conversion factor between pressure and ocean !! depth in m, usually 1/(rho_0*g), in m Pa-1. real :: C_p !< The heat capacity of seawater, in J K-1 kg-1. @@ -243,11 +240,10 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn) ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mdl = "ocean_model_init" ! This module's name. - character(len=48) :: stagger - integer :: secs, days + character(len=48) :: stagger ! A string indicating the staggering locations for the + ! surface velocities returned to the coupler. type(param_file_type) :: param_file !< A structure to parse for run-time parameters - logical :: use_temperature - type(time_type) :: dt_geometric, dt_savedays, dt_from_base + logical :: use_temperature ! If true, temperature and salinity are state variables. call callTree_enter("ocean_model_init(), ocean_model_MOM.F90") if (associated(OS)) then @@ -260,7 +256,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn) OS%is_ocean_pe = Ocean_sfc%is_ocean_pe if (.not.OS%is_ocean_pe) return - OS%Time = Time_in + OS%Time = Time_in ; OS%Time_dyn = Time_in call initialize_MOM(OS%Time, Time_init, param_file, OS%dirs, OS%MOM_CSp, & OS%restart_CSp, Time_in, offline_tracer_mode=OS%offline_tracer_mode, & diag_ptr=OS%diag, count_calls=.true.) @@ -314,14 +310,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn) else ; call MOM_error(FATAL,"ocean_model_init: OCEAN_SURFACE_STAGGER = "// & trim(stagger)//" is invalid.") ; endif - call get_param(param_file, mdl, "RESTORE_SALINITY",OS%restore_salinity, & - "If true, the coupled driver will add a globally-balanced \n"//& - "fresh-water flux that drives sea-surface salinity \n"//& - "toward specified values.", default=.false.) - call get_param(param_file, mdl, "RESTORE_TEMPERATURE",OS%restore_temp, & - "If true, the coupled driver will add a \n"//& - "heat flux that drives sea-surface temperauture \n"//& - "toward specified values.", default=.false.) call get_param(param_file, mdl, "RHO_0", Rho0, & "The mean ocean density used with BOUSSINESQ true to \n"//& "calculate accelerations and the mass for conservation \n"//& @@ -346,7 +334,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn) do_integrals=.true., gas_fields_ocn=gas_fields_ocn) call surface_forcing_init(Time_in, OS%grid, param_file, OS%diag, & - OS%forcing_CSp, OS%restore_salinity, OS%restore_temp) + OS%forcing_CSp) if (OS%use_ice_shelf) then call initialize_ice_shelf(param_file, OS%grid, OS%Time, OS%ice_shelf_CSp, & @@ -390,8 +378,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn) call close_param_file(param_file) call diag_mediator_close_registration(OS%diag) - if (is_root_pe()) & - write(*,'(/12x,a/)') '======== COMPLETED MOM INITIALIZATION ========' + call MOM_mesg('==== Completed MOM6 Coupled Initialization ====', 2) call callTree_leave("ocean_model_init(") end subroutine ocean_model_init @@ -401,23 +388,22 @@ end subroutine ocean_model_init !! time time_start_update) for a time interval of Ocean_coupling_time_step, !! returning the publicly visible ocean surface properties in Ocean_sfc and !! storing the new ocean properties in Ocean_state. -subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & - time_start_update, Ocean_coupling_time_step, & - update_dyn, update_thermo, Ocn_fluxes_used) +subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_update, & + Ocean_coupling_time_step, update_dyn, update_thermo, & + Ocn_fluxes_used, start_cycle, end_cycle, cycle_length) type(ice_ocean_boundary_type), & - intent(in) :: Ice_ocean_boundary !< A structure containing the - !! various forcing fields coming from the ice. + intent(in) :: Ice_ocean_boundary !< A structure containing the various + !! forcing fields coming from the ice and atmosphere. type(ocean_state_type), & - pointer :: OS !< A pointer to a private structure containing - !! the internal ocean state. + pointer :: OS !< A pointer to a private structure containing the + !! internal ocean state. type(ocean_public_type), & - intent(inout) :: Ocean_sfc !< A structure containing all the - !! publicly visible ocean surface fields after - !! a coupling time step. The data in this type is - !! intent out. + intent(inout) :: Ocean_sfc !< A structure containing all the publicly visible + !! ocean surface fields after a coupling time step. + !! The data in this type is intent out. type(time_type), intent(in) :: time_start_update !< The time at the beginning of the update step. - type(time_type), intent(in) :: Ocean_coupling_time_step !< The amount of time over - !! which to advance the ocean. + type(time_type), intent(in) :: Ocean_coupling_time_step !< The amount of time over which to + !! advance the ocean. logical, optional, intent(in) :: update_dyn !< If present and false, do not do updates !! due to the ocean dynamics. logical, optional, intent(in) :: update_thermo !< If present and false, do not do updates @@ -425,39 +411,38 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & logical, optional, intent(in) :: Ocn_fluxes_used !< If present, this indicates whether the !! cumulative thermodynamic fluxes from the ocean, !! like frazil, have been used and should be reset. + logical, optional, intent(in) :: start_cycle !< This indicates whether this call is to be + !! treated as the first call to step_MOM in a + !! time-stepping cycle; missing is like true. + logical, optional, intent(in) :: end_cycle !< This indicates whether this call is to be + !! treated as the last call to step_MOM in a + !! time-stepping cycle; missing is like true. + real, optional, intent(in) :: cycle_length !< The duration of a coupled time stepping cycle, in s. + ! Local variables - type(time_type) :: Master_time ! This allows step_MOM to temporarily change - ! the time that is seen by internal modules. - type(time_type) :: Time1 ! The value of the ocean model's time at the - ! start of a call to step_MOM. - integer :: index_bnds(4) ! The computational domain index bounds in the - ! ice-ocean boundary type. - real :: weight ! Flux accumulation weight - real :: dt_coupling ! The coupling time step in seconds. - integer :: nts ! The number of baroclinic dynamics time steps - ! within dt_coupling. - real :: dt_therm ! A limited and quantized version of OS%dt_therm (sec) - real :: dt_dyn ! The dynamics time step in sec. - real :: dtdia ! The diabatic time step in sec. - real :: t_elapsed_seg ! The elapsed time in this update segment, in s. - integer :: n, n_max, n_last_thermo - type(time_type) :: Time2 ! A temporary time. - logical :: thermo_does_span_coupling ! If true, thermodynamic forcing spans - ! multiple dynamic timesteps. - logical :: do_dyn ! If true, step the ocean dynamics and transport. - logical :: do_thermo ! If true, step the ocean thermodynamics. - logical :: step_thermo ! If true, take a thermodynamic step. - integer :: secs, days + type(time_type) :: Time_seg_start ! Stores the ocean model time at the start of this call to allow + ! step_MOM to temporarily change the time as seen by internal modules. + type(time_type) :: Time1 ! The value of the ocean model's time at the start of a call to step_MOM. + integer :: index_bnds(4) ! The computational domain index bounds in the ice-ocean boundary type. + real :: weight ! Flux accumulation weight of the current fluxes. + real :: dt_coupling ! The coupling time step in seconds. + real :: dt_therm ! A limited and quantized version of OS%dt_therm (sec) + real :: dt_dyn ! The dynamics time step in sec. + real :: dtdia ! The diabatic time step in sec. + real :: t_elapsed_seg ! The elapsed time in this update segment, in s. + integer :: n ! The internal iteration counter. + integer :: nts ! The number of baroclinic dynamics time steps in a thermodynamic step. + integer :: n_max ! The number of calls to step_MOM dynamics in this call to update_ocean_model. + integer :: n_last_thermo ! The iteration number the last time thermodynamics were updated. + logical :: thermo_does_span_coupling ! If true, thermodynamic forcing spans multiple dynamic timesteps. + logical :: do_dyn ! If true, step the ocean dynamics and transport. + logical :: do_thermo ! If true, step the ocean thermodynamics. + logical :: step_thermo ! If true, take a thermodynamic step. integer :: is, ie, js, je call callTree_enter("update_ocean_model(), ocean_model_MOM.F90") - call get_time(Ocean_coupling_time_step, secs, days) - dt_coupling = 86400.0*real(days) + real(secs) + dt_coupling = time_type_to_real(Ocean_coupling_time_step) - if (time_start_update /= OS%Time) then - call MOM_error(WARNING, "update_ocean_model: internal clock does not "//& - "agree with time_start_update argument.") - endif if (.not.associated(OS)) then call MOM_error(FATAL, "update_ocean_model called with an unassociated "// & "ocean_state_type structure. ocean_model_init must be "// & @@ -468,113 +453,111 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & do_dyn = .true. ; if (present(update_dyn)) do_dyn = update_dyn do_thermo = .true. ; if (present(update_thermo)) do_thermo = update_thermo + if (do_thermo .and. (time_start_update /= OS%Time)) & + call MOM_error(WARNING, "update_ocean_model: internal clock does not "//& + "agree with time_start_update argument.") + if (do_dyn .and. (time_start_update /= OS%Time_dyn)) & + call MOM_error(WARNING, "update_ocean_model: internal dynamics clock does not "//& + "agree with time_start_update argument.") + + if (.not.(do_dyn .or. do_thermo)) call MOM_error(FATAL, & + "update_ocean_model called without updating either dynamics or thermodynamics.") + if (do_dyn .and. do_thermo .and. (OS%Time /= OS%Time_dyn)) call MOM_error(FATAL, & + "update_ocean_model called to update both dynamics and thermodynamics with inconsistent clocks.") + ! This is benign but not necessary if ocean_model_init_sfc was called or if ! OS%sfc_state%tr_fields was spawned in ocean_model_init. Consider removing it. is = OS%grid%isc ; ie = OS%grid%iec ; js = OS%grid%jsc ; je = OS%grid%jec call coupler_type_spawn(Ocean_sfc%fields, OS%sfc_state%tr_fields, & (/is,is,ie,ie/), (/js,js,je,je/), as_needed=.true.) - ! Translate Ice_ocean_boundary into fluxes. + ! Translate Ice_ocean_boundary into fluxes and forces. call mpp_get_compute_domain(Ocean_sfc%Domain, index_bnds(1), index_bnds(2), & index_bnds(3), index_bnds(4)) - weight = 1.0 - - call convert_IOB_to_forces(Ice_ocean_boundary, OS%forces, index_bnds, OS%Time, & - OS%grid, OS%forcing_CSp) + if (do_dyn) then + call convert_IOB_to_forces(Ice_ocean_boundary, OS%forces, index_bnds, OS%Time_dyn, OS%grid, & + OS%forcing_CSp, dt_forcing=dt_coupling, reset_avg=OS%fluxes%fluxes_used) + if (OS%use_ice_shelf) & + call add_shelf_forces(OS%grid, OS%Ice_shelf_CSp, OS%forces) + if (OS%icebergs_alter_ocean) & + call iceberg_forces(OS%grid, OS%forces, OS%use_ice_shelf, & + OS%sfc_state, dt_coupling, OS%marine_ice_CSp) + endif - if (OS%fluxes%fluxes_used) then - if (do_thermo) & + if (do_thermo) then + if (OS%fluxes%fluxes_used) then call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%fluxes, index_bnds, OS%Time, & - OS%grid, OS%forcing_CSp, OS%sfc_state, & - OS%restore_salinity, OS%restore_temp) + OS%grid, OS%forcing_CSp, OS%sfc_state) - ! Add ice shelf fluxes - if (OS%use_ice_shelf) then - if (do_thermo) & + ! Add ice shelf fluxes + if (OS%use_ice_shelf) & call shelf_calc_flux(OS%sfc_state, OS%fluxes, OS%Time, dt_coupling, OS%Ice_shelf_CSp) - if (do_dyn) & - call add_shelf_forces(OS%grid, OS%Ice_shelf_CSp, OS%forces) - endif - if (OS%icebergs_alter_ocean) then - if (do_dyn) & - call iceberg_forces(OS%grid, OS%forces, OS%use_ice_shelf, & - OS%sfc_state, dt_coupling, OS%marine_ice_CSp) - if (do_thermo) & + if (OS%icebergs_alter_ocean) & call iceberg_fluxes(OS%grid, OS%fluxes, OS%use_ice_shelf, & - OS%sfc_state, dt_coupling, OS%marine_ice_CSp) - endif - - ! Fields that exist in both the forcing and mech_forcing types must be copied. - call copy_common_forcing_fields(OS%forces, OS%fluxes, OS%grid, skip_pres=.true.) + OS%sfc_state, dt_coupling, OS%marine_ice_CSp) #ifdef _USE_GENERIC_TRACER - call enable_averaging(dt_coupling, OS%Time + Ocean_coupling_time_step, OS%diag) !Is this needed? - call MOM_generic_tracer_fluxes_accumulate(OS%fluxes, weight) !here weight=1, just saving the current fluxes + call enable_averaging(dt_coupling, OS%Time + Ocean_coupling_time_step, OS%diag) !Is this needed? + call MOM_generic_tracer_fluxes_accumulate(OS%fluxes, 1.0) ! Here weight=1, so just store the current fluxes + call disable_averaging(OS%diag) #endif - ! Indicate that there are new unused fluxes. - OS%fluxes%fluxes_used = .false. - OS%fluxes%dt_buoy_accum = dt_coupling - else - OS%flux_tmp%C_p = OS%fluxes%C_p - if (do_thermo) & + ! Indicate that there are new unused fluxes. + OS%fluxes%fluxes_used = .false. + OS%fluxes%dt_buoy_accum = dt_coupling + else + ! The previous fluxes have not been used yet, so translate the input fluxes + ! into a temporary type and then accumulate them in about 20 lines. + OS%flux_tmp%C_p = OS%fluxes%C_p call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%flux_tmp, index_bnds, OS%Time, & - OS%grid, OS%forcing_CSp, OS%sfc_state, OS%restore_salinity,OS%restore_temp) + OS%grid, OS%forcing_CSp, OS%sfc_state) - if (OS%use_ice_shelf) then - if (do_thermo) & + if (OS%use_ice_shelf) & call shelf_calc_flux(OS%sfc_state, OS%flux_tmp, OS%Time, dt_coupling, OS%Ice_shelf_CSp) - if (do_dyn) & - call add_shelf_forces(OS%grid, OS%Ice_shelf_CSp, OS%forces) - endif - if (OS%icebergs_alter_ocean) then - if (do_dyn) & - call iceberg_forces(OS%grid, OS%forces, OS%use_ice_shelf, & - OS%sfc_state, dt_coupling, OS%marine_ice_CSp) - if (do_thermo) & + if (OS%icebergs_alter_ocean) & call iceberg_fluxes(OS%grid, OS%flux_tmp, OS%use_ice_shelf, & - OS%sfc_state, dt_coupling, OS%marine_ice_CSp) - endif - - call forcing_accumulate(OS%flux_tmp, OS%forces, OS%fluxes, dt_coupling, OS%grid, weight) - ! Some of the fields that exist in both the forcing and mech_forcing types - ! (e.g., ustar) are time-averages must be copied back to the forces type. - call copy_back_forcing_fields(OS%fluxes, OS%forces, OS%grid) + OS%sfc_state, dt_coupling, OS%marine_ice_CSp) + call fluxes_accumulate(OS%flux_tmp, OS%fluxes, dt_coupling, OS%grid, weight) #ifdef _USE_GENERIC_TRACER - call MOM_generic_tracer_fluxes_accumulate(OS%flux_tmp, weight) !weight of the current flux in the running average + ! Incorporate the current tracer fluxes into the running averages + call MOM_generic_tracer_fluxes_accumulate(OS%flux_tmp, weight) #endif + endif endif - call set_derived_forcing_fields(OS%forces, OS%fluxes, OS%grid, OS%GV%Rho0) - call set_net_mass_forcing(OS%fluxes, OS%forces, OS%grid) - if (OS%use_waves) then + ! The net mass forcing is not currently used in the MOM6 dynamics solvers, so this is may be unnecessary. + if (do_dyn .and. associated(OS%forces%net_mass_src) .and. .not.OS%forces%net_mass_src_set) & + call get_net_mass_forcing(OS%fluxes, OS%grid, OS%forces%net_mass_src) + + if (OS%use_waves .and. do_thermo) then + ! For now, the waves are only updated on the thermodynamics steps, because that is where + ! the wave intensities are actually used to drive mixing. At some point, the wave updates + ! might also need to become a part of the ocean dynamics, according to B. Reichl. call Update_Surface_Waves(OS%grid, OS%GV, OS%time, ocean_coupling_time_step, OS%waves) endif - if (OS%nstep==0) then + if ((OS%nstep==0) .and. (OS%nstep_thermo==0)) then ! This is the first call to update_ocean_model. call finish_MOM_initialization(OS%Time, OS%dirs, OS%MOM_CSp, OS%restart_CSp) endif - call disable_averaging(OS%diag) - Master_time = OS%Time ; Time1 = OS%Time + Time_seg_start = OS%Time ; if (do_dyn) Time_seg_start = OS%Time_dyn + Time1 = Time_seg_start - if (OS%offline_tracer_mode) then + if (OS%offline_tracer_mode .and. do_thermo) then call step_offline(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp) elseif ((.not.do_thermo) .or. (.not.do_dyn)) then ! The call sequence is being orchestrated from outside of update_ocean_model. call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp, & - Waves=OS%Waves, do_dynamics=do_thermo, do_thermodynamics=do_dyn, & + Waves=OS%Waves, do_dynamics=do_dyn, do_thermodynamics=do_thermo, & + start_cycle=start_cycle, end_cycle=end_cycle, cycle_length=cycle_length, & reset_therm=Ocn_fluxes_used) - !### What to do with these? , start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling) - elseif (OS%single_step_call) then call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_coupling, OS%MOM_CSp, Waves=OS%Waves) else n_max = 1 ; if (dt_coupling > OS%dt) n_max = ceiling(dt_coupling/OS%dt - 0.001) dt_dyn = dt_coupling / real(n_max) - thermo_does_span_coupling = (OS%thermo_spans_coupling .and. & - (OS%dt_therm > 1.5*dt_coupling)) + thermo_does_span_coupling = (OS%thermo_spans_coupling .and. (OS%dt_therm > 1.5*dt_coupling)) if (thermo_does_span_coupling) then dt_therm = dt_coupling * floor(OS%dt_therm / dt_coupling + 0.001) @@ -584,7 +567,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & n_last_thermo = 0 endif - Time2 = Time1 ; t_elapsed_seg = 0.0 + Time1 = Time_seg_start ; t_elapsed_seg = 0.0 do n=1,n_max if (OS%diabatic_first) then if (thermo_does_span_coupling) call MOM_error(FATAL, & @@ -592,16 +575,16 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & "THERMO_SPANS_COUPLING and DIABATIC_FIRST.") if (modulo(n-1,nts)==0) then dtdia = dt_dyn*min(nts,n_max-(n-1)) - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dtdia, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dtdia, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=.false., do_thermodynamics=.true., & start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling) endif - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dt_dyn, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_dyn, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=.true., do_thermodynamics=.false., & start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling) else - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dt_dyn, OS%MOM_CSp, & + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dt_dyn, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=.true., do_thermodynamics=.false., & start_cycle=(n==1), end_cycle=.false., cycle_length=dt_coupling) @@ -616,28 +599,32 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & endif if (step_thermo) then - ! Back up Time2 to the start of the thermodynamic segment. - Time2 = Time2 - set_time(int(floor((dtdia - dt_dyn) + 0.5))) - call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time2, dtdia, OS%MOM_CSp, & + ! Back up Time1 to the start of the thermodynamic segment. + Time1 = Time1 - real_to_time(dtdia - dt_dyn) + call step_MOM(OS%forces, OS%fluxes, OS%sfc_state, Time1, dtdia, OS%MOM_CSp, & Waves=OS%Waves, do_dynamics=.false., do_thermodynamics=.true., & start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_coupling) endif endif t_elapsed_seg = t_elapsed_seg + dt_dyn - Time2 = Time1 + set_time(int(floor(t_elapsed_seg + 0.5))) + Time1 = Time_seg_start + real_to_time(t_elapsed_seg) enddo endif - OS%Time = Master_time + Ocean_coupling_time_step - OS%nstep = OS%nstep + 1 + if (do_dyn) OS%Time_dyn = Time_seg_start + Ocean_coupling_time_step + if (do_dyn) OS%nstep = OS%nstep + 1 + OS%Time = Time_seg_start ! Reset the clock to compensate for shared pointers. + if (do_thermo) OS%Time = OS%Time + Ocean_coupling_time_step + if (do_thermo) OS%nstep_thermo = OS%nstep_thermo + 1 - call enable_averaging(dt_coupling, OS%Time, OS%diag) - call mech_forcing_diags(OS%forces, OS%fluxes, dt_coupling, OS%grid, & - OS%diag, OS%forcing_CSp%handles) - call disable_averaging(OS%diag) + if (do_dyn) then + call enable_averaging(dt_coupling, OS%Time_dyn, OS%diag) + call mech_forcing_diags(OS%forces, dt_coupling, OS%grid, OS%diag, OS%forcing_CSp%handles) + call disable_averaging(OS%diag) + endif - if (OS%fluxes%fluxes_used) then + if (OS%fluxes%fluxes_used .and. do_thermo) then call enable_averaging(OS%fluxes%dt_buoy_accum, OS%Time, OS%diag) call forcing_diagnostics(OS%fluxes, OS%sfc_state, OS%fluxes%dt_buoy_accum, & OS%grid, OS%diag, OS%forcing_CSp%handles) @@ -648,7 +635,8 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & ! call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid, & ! Ice_ocean_boundary%p, OS%press_to_z) call convert_state_to_ocean_type(OS%sfc_state, Ocean_sfc, OS%grid) - call coupler_type_send_data(Ocean_sfc%fields, OS%Time) + Time1 = OS%Time ; if (do_dyn) Time1 = OS%Time_dyn + call coupler_type_send_data(Ocean_sfc%fields, Time1) call callTree_leave("update_ocean_model()") end subroutine update_ocean_model diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index d294c29656..d2de157a49 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -1790,8 +1790,7 @@ subroutine update_ocean_model(OS, Ocean_sfc, time_start_update, & 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 diff --git a/config_src/solo_driver/MESO_surface_forcing.F90 b/config_src/solo_driver/MESO_surface_forcing.F90 index eaa11da6c1..68852f89d9 100644 --- a/config_src/solo_driver/MESO_surface_forcing.F90 +++ b/config_src/solo_driver/MESO_surface_forcing.F90 @@ -12,14 +12,14 @@ module MESO_surface_forcing use MOM_forcing_type, only : allocate_forcing_type, allocate_mech_forcing use MOM_grid, only : ocean_grid_type use MOM_io, only : file_exists, MOM_read_data, slasher -use MOM_time_manager, only : time_type, operator(+), operator(/), get_time +use MOM_time_manager, only : time_type, operator(+), operator(/) use MOM_tracer_flow_control, only : call_tracer_set_forcing use MOM_tracer_flow_control, only : tracer_flow_control_CS use MOM_variables, only : surface implicit none ; private -public MESO_wind_forcing, MESO_buoyancy_forcing, MESO_surface_forcing_init +public MESO_buoyancy_forcing, MESO_surface_forcing_init !> This control structure is used to store parameters associated with the MESO forcing. type, public :: MESO_surface_forcing_CS ; private @@ -52,71 +52,6 @@ module MESO_surface_forcing contains -!### This subroutine sets zero surface wind stresses, but it is not even -!### used by the MESO experimeents. This subroutine can be deleted. -RWH -subroutine MESO_wind_forcing(sfc_state, forces, day, G, CS) - type(surface), intent(inout) :: sfc_state !< A structure containing fields that - !! describe the surface state of the ocean. - type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces - type(time_type), intent(in) :: day !< The time of the fluxes - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure - type(MESO_surface_forcing_CS), pointer :: CS !< A pointer to the control structure returned by a previous - !! call to MESO_surface_forcing_init - -! This subroutine sets the surface wind stresses, forces%taux and forces%tauy. -! These are the stresses in the direction of the model grid (i.e. the same -! direction as the u- and v- velocities.) They are both in Pa. -! In addition, this subroutine can be used to set the surface friction -! velocity, forces%ustar, in m s-1. This is needed with a bulk mixed layer. -! -! Arguments: state - A structure containing fields that describe the -! surface state of the ocean. -! (out) fluxes - A structure containing pointers to any possible -! forcing fields. Unused fields have NULL ptrs. -! (in) day - Time of the fluxes. -! (in) G - The ocean's grid structure. -! (in) CS - A pointer to the control structure returned by a previous -! call to MESO_surface_forcing_init - - integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq - integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB - - ! When modifying the code, comment out this error message. It is here - ! so that the original (unmodified) version is not accidentally used. - call MOM_error(FATAL, "MESO_wind_surface_forcing: " // & - "User forcing routine called without modification." ) - - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec - Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB - - ! Allocate the forcing arrays, if necessary. - call allocate_mech_forcing(G, forces, stress=.true., ustar=.true.) - - ! Set the surface wind stresses, in units of Pa. A positive taux - ! accelerates the ocean to the (pseudo-)east. - - ! The i-loop extends to is-1 so that taux can be used later in the - ! calculation of ustar - otherwise the lower bound would be Isq. - do j=js,je ; do I=is-1,Ieq - forces%taux(I,j) = G%mask2dCu(I,j) * 0.0 ! Change this to the desired expression. - enddo ; enddo - do J=js-1,Jeq ; do i=is,ie - forces%tauy(i,J) = G%mask2dCv(i,J) * 0.0 ! Change this to the desired expression. - enddo ; enddo - - ! Set the surface friction velocity, in units of m s-1. ustar - ! is always positive. - if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie - ! This expression can be changed if desired, but need not be. - forces%ustar(i,j) = G%mask2dT(i,j) * sqrt(CS%gust_const/CS%Rho0 + & - sqrt(0.5*(forces%taux(I-1,j)**2 + forces%taux(I,j)**2) + & - 0.5*(forces%tauy(i,J-1)**2 + forces%tauy(i,J)**2))/CS%Rho0) - enddo ; enddo ; endif - -end subroutine MESO_wind_forcing - !> This subroutine sets up the MESO buoyancy forcing, which uses control-theory style !! specification restorative buoyancy fluxes at large scales. subroutine MESO_buoyancy_forcing(sfc_state, fluxes, day, dt, G, CS) @@ -130,10 +65,6 @@ subroutine MESO_buoyancy_forcing(sfc_state, fluxes, day, dt, G, CS) type(MESO_surface_forcing_CS), pointer :: CS !< A pointer to the control structure returned by !! a previous call to MESO_surface_forcing_init -! This subroutine specifies the current surface fluxes of buoyancy or -! temperature and fresh water. It may also be modified to add -! surface fluxes of user provided tracers. - ! When temperature is used, there are long list of fluxes that need to be ! set - essentially the same as for a full coupled model, but most of these ! can be simply set to zero. The net fresh water flux should probably be @@ -144,17 +75,6 @@ subroutine MESO_buoyancy_forcing(sfc_state, fluxes, day, dt, G, CS) ! are in W m-2 and positive for heat going into the ocean. All fresh water ! fluxes are in kg m-2 s-1 and positive for water moving into the ocean. -! Arguments: state - A structure containing fields that describe the -! surface state of the ocean. -! (out) fluxes - A structure containing pointers to any possible -! forcing fields. Unused fields have NULL ptrs. -! (in) day_start - Start time of the fluxes. -! (in) day_interval - Length of time over which these fluxes -! will be applied. -! (in) G - The ocean's grid structure. -! (in) CS - A pointer to the control structure returned by a previous -! call to MESO_surface_forcing_init - real :: Temp_restore ! The temperature that is being restored toward, in C. real :: Salin_restore ! The salinity that is being restored toward, in PSU. real :: density_restore ! The potential density that is being restored @@ -293,14 +213,6 @@ subroutine MESO_surface_forcing_init(Time, G, param_file, diag, CS) type(MESO_surface_forcing_CS), pointer :: CS !< A pointer that is set to point to the !! control structure for this module -! Arguments: Time - The current model time. -! (in) G - The ocean's grid structure. -! (in) param_file - A structure indicating the open file to parse for -! model parameter values. -! (in) diag - A structure that is used to regulate diagnostic output. -! (in/out) CS - A pointer that is set to point to the control structure -! for this module - ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mdl = "MESO_surface_forcing" ! This module's name. @@ -383,9 +295,6 @@ end subroutine MESO_surface_forcing_init !! it is probably a good idea to read the forcing from input files !! using "file" for WIND_CONFIG and BUOY_CONFIG. !! -!! MESO_wind_forcing should set the surface wind stresses (taux and -!! tauy) perhaps along with the surface friction velocity (ustar). -!! !! MESO_buoyancy forcing is used to set the surface buoyancy !! forcing, which may include a number of fresh water flux fields !! (evap, liq_precip, froz_precip, liq_runoff, froz_runoff, and diff --git a/config_src/solo_driver/MOM_driver.F90 b/config_src/solo_driver/MOM_driver.F90 index 61c3f4a509..4933f29182 100644 --- a/config_src/solo_driver/MOM_driver.F90 +++ b/config_src/solo_driver/MOM_driver.F90 @@ -48,7 +48,8 @@ program MOM_main use MOM_string_functions,only : uppercase use MOM_surface_forcing, only : set_forcing, forcing_save_restart use MOM_surface_forcing, only : surface_forcing_init, surface_forcing_CS - use MOM_time_manager, only : time_type, set_date, set_time, get_date, time_type_to_real + use MOM_time_manager, only : time_type, set_date, get_date + use MOM_time_manager, only : real_to_time, time_type_to_real use MOM_time_manager, only : operator(+), operator(-), operator(*), operator(/) use MOM_time_manager, only : operator(>), operator(<), operator(>=) use MOM_time_manager, only : increment_date, set_calendar_type, month_name @@ -137,7 +138,7 @@ program MOM_main real :: dt_dyn, dtdia, t_elapsed_seg integer :: n, n_max, nts, n_last_thermo logical :: diabatic_first, single_step_call - type(time_type) :: Time2 + type(time_type) :: Time2, time_chg integer :: Restart_control ! An integer that is bit-tested to determine whether ! incremental restart files are saved and whether they @@ -290,7 +291,7 @@ program MOM_main Start_time = set_date(date_init(1),date_init(2), date_init(3), & date_init(4),date_init(5),date_init(6)) else - Start_time = set_time(0,days=0) + Start_time = real_to_time(0.0) endif call time_interp_external_init @@ -356,7 +357,7 @@ program MOM_main endif ntstep = MAX(1,ceiling(dt_forcing/dt - 0.001)) - Time_step_ocean = set_time(int(floor(dt_forcing+0.5))) + Time_step_ocean = real_to_time(dt_forcing) elapsed_time_master = (abs(dt_forcing - time_type_to_real(Time_step_ocean)) > 1.0e-12*dt_forcing) if (elapsed_time_master) & call MOM_mesg("Using real elapsed time for the master clock.", 2) @@ -415,7 +416,7 @@ program MOM_main call get_param(param_file, mod_name, "RESTINT", restint, & "The interval between saves of the restart file in units \n"//& "of TIMEUNIT. Use 0 (the default) to not save \n"//& - "incremental restart files at all.", default=set_time(0), & + "incremental restart files at all.", default=real_to_time(0.0), & timeunit=Time_unit) call get_param(param_file, mod_name, "WRITE_CPU_STEPS", cpu_steps, & "The number of coupled timesteps between writing the cpu \n"//& @@ -454,7 +455,7 @@ program MOM_main if (((.not.BTEST(Restart_control,1)) .and. (.not.BTEST(Restart_control,0))) & .or. (Restart_control < 0)) permit_incr_restart = .false. - if (restint > set_time(0)) then + if (restint > real_to_time(0.0)) then ! restart_time is the next integral multiple of restint. restart_time = Start_time + restint * & (1 + ((Time + Time_step_ocean) - Start_time) / restint) @@ -532,7 +533,7 @@ program MOM_main dtdia = dt_dyn*(n - n_last_thermo) ! Back up Time2 to the start of the thermodynamic segment. if (n > n_last_thermo+1) & - Time2 = Time2 - set_time(int(floor((dtdia - dt_dyn) + 0.5))) + Time2 = Time2 - real_to_time(dtdia - dt_dyn) call step_MOM(forces, fluxes, sfc_state, Time2, dtdia, MOM_CSp, & do_dynamics=.false., do_thermodynamics=.true., & start_cycle=.false., end_cycle=(n==n_max), cycle_length=dt_forcing) @@ -541,7 +542,7 @@ program MOM_main endif t_elapsed_seg = t_elapsed_seg + dt_dyn - Time2 = Time1 + set_time(int(floor(t_elapsed_seg + 0.5))) + Time2 = Time1 + real_to_time(t_elapsed_seg) enddo endif @@ -549,17 +550,17 @@ program MOM_main ! This is here to enable fractional-second time steps. elapsed_time = elapsed_time + dt_forcing if (elapsed_time > 2e9) then - ! This is here to ensure that the conversion from a real to an integer - ! can be accurately represented in long runs (longer than ~63 years). - ! It will also ensure that elapsed time does not lose resolution of order - ! the timetype's resolution, provided that the timestep and tick are - ! larger than 10-5 seconds. If a clock with a finer resolution is used, - ! a smaller value would be required. - segment_start_time = segment_start_time + set_time(int(floor(elapsed_time))) - elapsed_time = elapsed_time - floor(elapsed_time) + ! This is here to ensure that the conversion from a real to an integer can be accurately + ! represented in long runs (longer than ~63 years). It will also ensure that elapsed time + ! does not lose resolution of order the timetype's resolution, provided that the timestep and + ! tick are larger than 10-5 seconds. If a clock with a finer resolution is used, a smaller + ! value would be required. + time_chg = real_to_time(elapsed_time) + segment_start_time = segment_start_time + time_chg + elapsed_time = elapsed_time - time_type_to_real(time_chg) endif if (elapsed_time_master) then - Master_Time = segment_start_time + set_time(int(floor(elapsed_time+0.5))) + Master_Time = segment_start_time + real_to_time(elapsed_time) else Master_Time = Master_Time + Time_step_ocean endif @@ -570,8 +571,7 @@ program MOM_main endif ; endif call enable_averaging(dt_forcing, Time, diag) - call mech_forcing_diags(forces, fluxes, dt_forcing, grid, diag, & - surface_forcing_CSp%handles) + call mech_forcing_diags(forces, dt_forcing, grid, diag, surface_forcing_CSp%handles) call disable_averaging(diag) if (.not. offline_tracer_mode) then diff --git a/config_src/solo_driver/MOM_surface_forcing.F90 b/config_src/solo_driver/MOM_surface_forcing.F90 index 351b149830..a3a9a12204 100644 --- a/config_src/solo_driver/MOM_surface_forcing.F90 +++ b/config_src/solo_driver/MOM_surface_forcing.F90 @@ -32,11 +32,11 @@ module MOM_surface_forcing use MOM_io, only : EAST_FACE, NORTH_FACE, num_timelevels use MOM_restart, only : register_restart_field, restart_init, MOM_restart_CS use MOM_restart, only : restart_init_end, save_restart, restore_state -use MOM_time_manager, only : time_type, operator(+), operator(/), get_time, set_time +use MOM_time_manager, only : time_type, operator(+), operator(/), get_time, time_type_to_real use MOM_tracer_flow_control, only : call_tracer_set_forcing use MOM_tracer_flow_control, only : tracer_flow_control_CS use MOM_variables, only : surface -use MESO_surface_forcing, only : MESO_wind_forcing, MESO_buoyancy_forcing +use MESO_surface_forcing, only : MESO_buoyancy_forcing use MESO_surface_forcing, only : MESO_surface_forcing_init, MESO_surface_forcing_CS use Neverland_surface_forcing, only : Neverland_wind_forcing, Neverland_buoyancy_forcing use Neverland_surface_forcing, only : Neverland_surface_forcing_init, Neverland_surface_forcing_CS @@ -226,7 +226,6 @@ subroutine set_forcing(sfc_state, forces, fluxes, day_start, day_interval, G, CS ! Local variables real :: dt ! length of time in seconds over which fluxes applied type(time_type) :: day_center ! central time of the fluxes. - integer :: intdt integer :: isd, ied, jsd, jed isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -234,8 +233,7 @@ subroutine set_forcing(sfc_state, forces, fluxes, day_start, day_interval, G, CS call callTree_enter("set_forcing, MOM_surface_forcing.F90") day_center = day_start + day_interval/2 - call get_time(day_interval, intdt) - dt = real(intdt) + dt = time_type_to_real(day_interval) if (CS%first_call_set_forcing) then ! Allocate memory for the mechanical and thermodyanmic forcing fields. @@ -275,8 +273,6 @@ subroutine set_forcing(sfc_state, forces, fluxes, day_start, day_interval, G, CS call wind_forcing_const(sfc_state, forces, 0., 0., day_center, G, CS) elseif (trim(CS%wind_config) == "const") then call wind_forcing_const(sfc_state, forces, CS%tau_x0, CS%tau_y0, day_center, G, CS) - elseif (trim(CS%wind_config) == "MESO") then - call MESO_wind_forcing(sfc_state, forces, day_center, G, CS%MESO_forcing_CSp) elseif (trim(CS%wind_config) == "Neverland") then call Neverland_wind_forcing(sfc_state, forces, day_center, G, CS%Neverland_forcing_CSp) elseif (trim(CS%wind_config) == "SCM_ideal_hurr") then @@ -369,13 +365,10 @@ subroutine wind_forcing_const(sfc_state, forces, tau_x0, tau_y0, day, G, CS) ! Local variables real :: mag_tau integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq - integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB call callTree_enter("wind_forcing_const, MOM_surface_forcing.F90") is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB !set steady surface wind stresses, in units of Pa. mag_tau = sqrt( tau_x0**2 + tau_y0**2) @@ -414,13 +407,10 @@ subroutine wind_forcing_2gyre(sfc_state, forces, day, G, CS) ! Local variables real :: PI integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq - integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB call callTree_enter("wind_forcing_2gyre, MOM_surface_forcing.F90") is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB !set the steady surface wind stresses, in units of Pa. PI = 4.0*atan(1.0) @@ -450,13 +440,10 @@ subroutine wind_forcing_1gyre(sfc_state, forces, day, G, CS) ! Local variables real :: PI integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq - integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB call callTree_enter("wind_forcing_1gyre, MOM_surface_forcing.F90") is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB ! set the steady surface wind stresses, in units of Pa. PI = 4.0*atan(1.0) @@ -484,25 +471,22 @@ subroutine wind_forcing_gyres(sfc_state, forces, day, G, CS) ! Local variables real :: PI, y integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq - integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB call callTree_enter("wind_forcing_gyres, MOM_surface_forcing.F90") is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB ! steady surface wind stresses (Pa) PI = 4.0*atan(1.0) - do j=jsd,jed ; do I=is-1,IedB + do j=js-1,je+1 ; do I=is-1,Ieq y = (G%geoLatCu(I,j)-CS%South_lat)/CS%len_lat forces%taux(I,j) = CS%gyres_taux_const + & ( CS%gyres_taux_sin_amp*sin(CS%gyres_taux_n_pis*PI*y) & + CS%gyres_taux_cos_amp*cos(CS%gyres_taux_n_pis*PI*y) ) enddo ; enddo - do J=js-1,JedB ; do i=isd,ied + do J=js-1,Jeq ; do i=is-1,ie+1 forces%tauy(i,J) = 0.0 enddo ; enddo @@ -535,16 +519,13 @@ subroutine wind_forcing_from_file(sfc_state, forces, day, G, CS) integer :: time_lev ! The time level that is used for a field. integer :: days, seconds integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq - integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB logical :: read_Ustar call callTree_enter("wind_forcing_from_file, MOM_surface_forcing.F90") is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB - call get_time(day,seconds,days) + call get_time(day, seconds, days) time_lev_daily = days - 365*floor(real(days) / 365.0) if (time_lev_daily < 31) then ; time_lev_monthly = 0 @@ -774,7 +755,7 @@ subroutine buoyancy_forcing_from_files(sfc_state, fluxes, day, dt, G, CS) Irho0 = 1.0/CS%Rho0 ! Read the buoyancy forcing file - call get_time(day,seconds,days) + call get_time(day, seconds, days) time_lev_daily = days - 365*floor(real(days) / 365.0) diff --git a/config_src/solo_driver/Neverland_surface_forcing.F90 b/config_src/solo_driver/Neverland_surface_forcing.F90 index 65a5ca1339..e6111b2a19 100644 --- a/config_src/solo_driver/Neverland_surface_forcing.F90 +++ b/config_src/solo_driver/Neverland_surface_forcing.F90 @@ -12,7 +12,7 @@ module Neverland_surface_forcing use MOM_forcing_type, only : allocate_forcing_type, allocate_mech_forcing use MOM_grid, only : ocean_grid_type use MOM_io, only : file_exists, read_data, slasher -use MOM_time_manager, only : time_type, operator(+), operator(/), get_time +use MOM_time_manager, only : time_type, operator(+), operator(/) use MOM_variables, only : surface implicit none ; private @@ -48,15 +48,15 @@ module Neverland_surface_forcing !! Neverland forcing configuration. subroutine Neverland_wind_forcing(sfc_state, forces, day, G, CS) type(surface), intent(inout) :: sfc_state !< A structure containing fields that - !! describe the surface state of the ocean. + !! describe the surface state of the ocean. type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces - type(time_type), intent(in) :: day !< Time used for determining the fluxes. - type(ocean_grid_type), intent(inout) :: G !< Grid structure. - type(Neverland_surface_forcing_CS), pointer :: CS !< Control structure for this module. - ! Local variable + type(time_type), intent(in) :: day !< Time used for determining the fluxes. + type(ocean_grid_type), intent(inout) :: G !< Grid structure. + type(Neverland_surface_forcing_CS), pointer :: CS !< Control structure for this module. + + ! Local variables integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB - real :: x, y real :: PI real :: tau_max, off @@ -110,26 +110,26 @@ subroutine Neverland_wind_forcing(sfc_state, forces, day, G, CS) end subroutine Neverland_wind_forcing !> Returns the value of a cosine-bell function evaluated at x/L - real function cosbell(x,L) +real function cosbell(x,L) - real , intent(in) :: x !< non-dimensional position - real , intent(in) :: L !< non-dimensional width - real :: PI !< 3.1415926... calculated as 4*atan(1) + real , intent(in) :: x !< non-dimensional position + real , intent(in) :: L !< non-dimensional width + real :: PI !< 3.1415926... calculated as 4*atan(1) - PI = 4.0*atan(1.0) - cosbell = 0.5 * (1 + cos(PI*MIN(ABS(x/L),1.0))) - end function cosbell + PI = 4.0*atan(1.0) + cosbell = 0.5 * (1 + cos(PI*MIN(ABS(x/L),1.0))) +end function cosbell !> Returns the value of a sin-spike function evaluated at x/L - real function spike(x,L) +real function spike(x,L) - real , intent(in) :: x !< non-dimensional position - real , intent(in) :: L !< non-dimensional width - real :: PI !< 3.1415926... calculated as 4*atan(1) + real , intent(in) :: x !< non-dimensional position + real , intent(in) :: L !< non-dimensional width + real :: PI !< 3.1415926... calculated as 4*atan(1) - PI = 4.0*atan(1.0) - spike = (1 - sin(PI*MIN(ABS(x/L),0.5))) - end function spike + PI = 4.0*atan(1.0) + spike = (1 - sin(PI*MIN(ABS(x/L),0.5))) +end function spike !> Surface fluxes of buoyancy for the Neverland configurations. diff --git a/config_src/solo_driver/user_surface_forcing.F90 b/config_src/solo_driver/user_surface_forcing.F90 index e0136abf0f..7a27c75e18 100644 --- a/config_src/solo_driver/user_surface_forcing.F90 +++ b/config_src/solo_driver/user_surface_forcing.F90 @@ -12,7 +12,7 @@ module user_surface_forcing use MOM_forcing_type, only : allocate_forcing_type, allocate_mech_forcing use MOM_grid, only : ocean_grid_type use MOM_io, only : file_exists, read_data -use MOM_time_manager, only : time_type, operator(+), operator(/), get_time +use MOM_time_manager, only : time_type, operator(+), operator(/) use MOM_tracer_flow_control, only : call_tracer_set_forcing use MOM_tracer_flow_control, only : tracer_flow_control_CS use MOM_variables, only : surface @@ -49,30 +49,15 @@ module user_surface_forcing !! direction as the u- and v- velocities.) They are both in Pa. subroutine USER_wind_forcing(sfc_state, forces, day, G, CS) type(surface), intent(inout) :: sfc_state !< A structure containing fields that - !! describe the surface state of the ocean. + !! describe the surface state of the ocean. type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces type(time_type), intent(in) :: day !< The time of the fluxes type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure type(user_surface_forcing_CS), pointer :: CS !< A pointer to the control structure returned !! by a previous call to user_surface_forcing_init -! This subroutine sets the surface wind stresses, forces%taux and forces%tauy. -! These are the stresses in the direction of the model grid (i.e. the same -! direction as the u- and v- velocities.) They are both in Pa. -! In addition, this subroutine can be used to set the surface friction -! velocity, forces%ustar, in m s-1. This is needed with a bulk mixed layer. -! -! Arguments: state - A structure containing fields that describe the -! surface state of the ocean. -! (out) fluxes - A structure containing pointers to any possible -! forcing fields. Unused fields have NULL ptrs. -! (in) day - Time of the fluxes. -! (in) G - The ocean's grid structure. -! (in) CS - A pointer to the control structure returned by a previous -! call to user_surface_forcing_init - + ! Local variables integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq - integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB ! When modifying the code, comment out this error message. It is here ! so that the original (unmodified) version is not accidentally used. @@ -81,8 +66,6 @@ subroutine USER_wind_forcing(sfc_state, forces, day, G, CS) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB ! Allocate the forcing arrays, if necessary. call allocate_mech_forcing(G, forces, stress=.true., ustar=.true.) @@ -138,22 +121,12 @@ subroutine USER_buoyancy_forcing(sfc_state, fluxes, day, dt, G, CS) ! are in W m-2 and positive for heat going into the ocean. All fresh water ! fluxes are in kg m-2 s-1 and positive for water moving into the ocean. -! Arguments: state - A structure containing fields that describe the -! surface state of the ocean. -! (out) fluxes - A structure containing pointers to any possible -! forcing fields. Unused fields have NULL ptrs. -! (in) day_start - Start time of the fluxes. -! (in) day_interval - Length of time over which these fluxes -! will be applied. -! (in) G - The ocean's grid structure. -! (in) CS - A pointer to the control structure returned by a previous -! call to user_surface_forcing_init - + ! Local variables real :: Temp_restore ! The temperature that is being restored toward, in C. real :: Salin_restore ! The salinity that is being restored toward, in PSU. real :: density_restore ! The potential density that is being restored ! toward, in kg m-3. - real :: rhoXcp ! The mean density times the heat capacity, in J m-3 K-1. + real :: rhoXcp ! The mean density times the heat capacity, in J m-3 K-1. real :: buoy_rest_const ! A constant relating density anomalies to the ! restoring buoyancy flux, in m5 s-3 kg-1. @@ -266,14 +239,6 @@ subroutine USER_surface_forcing_init(Time, G, param_file, diag, CS) type(user_surface_forcing_CS), pointer :: CS !< A pointer that is set to point to !! the control structure for this module -! Arguments: Time - The current model time. -! (in) G - The ocean's grid structure. -! (in) param_file - A structure indicating the open file to parse for -! model parameter values. -! (in) diag - A structure that is used to regulate diagnostic output. -! (in/out) CS - A pointer that is set to point to the control structure -! for this module - ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mdl = "user_surface_forcing" ! This module's name. diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index a71dfb557c..7e2885fd6f 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -760,14 +760,12 @@ subroutine remap_all_state_vars(CS_remapping, CS_ALE, G, GV, h_old, h_new, Reg, nz = GV%ke ppt2mks = 0.001 - if (associated(Reg)) then - ntr = Reg%ntr - else - ntr = 0 - endif + ntr = 0 ; if (associated(Reg)) ntr = Reg%ntr if (present(dt)) then Idt = 1.0/dt + work_conc(:,:,:) = 0.0 + work_cont(:,:,:) = 0.0 endif ! Remap tracer @@ -801,22 +799,23 @@ subroutine remap_all_state_vars(CS_remapping, CS_ALE, G, GV, h_old, h_new, Reg, endif ; enddo ; enddo ! tendency diagnostics. - if (Tr%id_remap_conc > 0) then - call post_data(Tr%id_remap_conc, work_conc, CS_ALE%diag) - endif - if (Tr%id_remap_cont > 0) then - call post_data(Tr%id_remap_cont, work_cont, CS_ALE%diag) - endif - if (Tr%id_remap_cont_2d > 0) then - do j = G%jsc,G%jec ; do i = G%isc,G%iec - work_2d(i,j) = 0.0 - do k = 1,GV%ke - work_2d(i,j) = work_2d(i,j) + work_cont(i,j,k) - enddo - enddo ; enddo - call post_data(Tr%id_remap_cont_2d, work_2d, CS_ALE%diag) + if (present(dt)) then + if (Tr%id_remap_conc > 0) then + call post_data(Tr%id_remap_conc, work_conc, CS_ALE%diag) + endif + if (Tr%id_remap_cont > 0) then + call post_data(Tr%id_remap_cont, work_cont, CS_ALE%diag) + endif + if (Tr%id_remap_cont_2d > 0) then + do j = G%jsc,G%jec ; do i = G%isc,G%iec + work_2d(i,j) = 0.0 + do k = 1,GV%ke + work_2d(i,j) = work_2d(i,j) + work_cont(i,j,k) + enddo + enddo ; enddo + call post_data(Tr%id_remap_cont_2d, work_2d, CS_ALE%diag) + endif endif - enddo ! m=1,ntr endif ! endif for ntr > 0 @@ -866,7 +865,7 @@ subroutine remap_all_state_vars(CS_remapping, CS_ALE, G, GV, h_old, h_new, Reg, endif if (CS_ALE%id_vert_remap_h > 0) call post_data(CS_ALE%id_vert_remap_h, h_old, CS_ALE%diag) - if (CS_ALE%id_vert_remap_h_tendency > 0) then + if ((CS_ALE%id_vert_remap_h_tendency > 0) .and. present(dt)) then do k = 1, nz ; do j = G%jsc,G%jec ; do i = G%isc,G%iec work_cont(i,j,k) = (h_new(i,j,k) - h_old(i,j,k))*Idt enddo ; enddo ; enddo diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index c554e4f92e..54b5197ad9 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -39,7 +39,7 @@ module MOM use MOM_restart, only : register_restart_field, query_initialized, save_restart use MOM_restart, only : restart_init, is_new_run, MOM_restart_CS use MOM_spatial_means, only : global_mass_integral -use MOM_time_manager, only : time_type, set_time, time_type_to_real, operator(+) +use MOM_time_manager, only : time_type, real_to_time, time_type_to_real, operator(+) use MOM_time_manager, only : operator(-), operator(>), operator(*), operator(/) use MOM_time_manager, only : operator(>=), increment_date use MOM_unit_tests, only : unit_tests @@ -556,7 +556,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & do j=js,je ; do i=is,ie ; CS%ssh_rint(i,j) = 0.0 ; enddo ; enddo if (associated(CS%VarMix)) then - call enable_averaging(cycle_time, Time_start+set_time(int(cycle_time)), & + call enable_averaging(cycle_time, Time_start + real_to_time(cycle_time), & CS%diag) call calc_resoln_function(h, CS%tv, G, GV, CS%VarMix) call disable_averaging(CS%diag) @@ -582,7 +582,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & if (CS%UseWaves) then ! Update wave information, which is presently kept static over each call to step_mom - call enable_averaging(time_interval, Time_start + set_time(int(floor(time_interval+0.5))), CS%diag) + call enable_averaging(time_interval, Time_start + real_to_time(time_interval), CS%diag) call Update_Stokes_Drift(G, GV, Waves, h, forces%ustar) call disable_averaging(CS%diag) endif @@ -604,9 +604,9 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & do n=1,n_max rel_time = rel_time + dt ! The relative time at the end of the step. ! Set the universally visible time to the middle of the time step. - CS%Time = Time_start + set_time(int(floor(0.5 + rel_time - 0.5*dt))) + CS%Time = Time_start + real_to_time(rel_time - 0.5*dt) ! Set the local time to the end of the time step. - Time_local = Time_start + set_time(int(floor(rel_time+0.5))) + Time_local = Time_start + real_to_time(rel_time) if (showCallTree) call callTree_enter("DT cycles (step_MOM) n=",n) @@ -633,10 +633,10 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & if (dtdia > dt) then ! If necessary, temporarily reset CS%Time to the center of the period covered ! by the call to step_MOM_thermo, noting that they begin at the same time. - CS%Time = CS%Time + set_time(int(floor(0.5*(dtdia-dt) + 0.5))) + CS%Time = CS%Time + real_to_time(0.5*(dtdia-dt)) ! The end-time of the diagnostic interval needs to be set ahead if there ! are multiple dynamic time steps worth of thermodynamics applied here. - end_time_thermo = Time_local + set_time(int(floor(dtdia-dt+0.5))) + end_time_thermo = Time_local + real_to_time(dtdia-dt) endif ! Apply diabatic forcing, do mixing, and regrid. @@ -649,7 +649,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & if (showCallTree) call callTree_waypoint("finished diabatic_first (step_MOM)") if (dtdia > dt) & ! Reset CS%Time to its previous value. - CS%Time = Time_start + set_time(int(floor(0.5 + rel_time - 0.5*dt))) + CS%Time = Time_start + real_to_time(rel_time - 0.5*dt) endif ! end of block "(CS%diabatic_first .and. (CS%t_dyn_rel_adv==0.0))" if (do_dyn) then @@ -731,7 +731,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & ! If necessary, temporarily reset CS%Time to the center of the period covered ! by the call to step_MOM_thermo, noting that they end at the same time. - if (dtdia > dt) CS%Time = CS%Time - set_time(int(floor(0.5*(dtdia-dt) + 0.5))) + if (dtdia > dt) CS%Time = CS%Time - real_to_time(0.5*(dtdia-dt)) ! Apply diabatic forcing, do mixing, and regrid. call step_MOM_thermo(CS, G, GV, u, v, h, CS%tv, fluxes, dtdia, & @@ -740,7 +740,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & CS%t_dyn_rel_thermo = 0.0 if (dtdia > dt) & ! Reset CS%Time to its previous value. - CS%Time = Time_start + set_time(int(floor(0.5 + rel_time - 0.5*dt))) + CS%Time = Time_start + real_to_time(rel_time - 0.5*dt) endif if (do_dyn) then @@ -774,7 +774,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & CS%t_dyn_rel_diag = 0.0 call cpu_clock_begin(id_clock_Z_diag) - if (Time_local + set_time(int(0.5*dt_therm)) > CS%Z_diag_time) then + if (Time_local + real_to_time(0.5*dt_therm) > CS%Z_diag_time) then call enable_averaging(real(time_type_to_real(CS%Z_diag_interval)), & CS%Z_diag_time, CS%diag) !### This is the one place where fluxes might used if do_thermo=.false. Is this correct? @@ -852,7 +852,7 @@ subroutine step_MOM(forces, fluxes, sfc_state, Time_start, time_interval, CS, & if (MOM_state_is_synchronized(CS)) & call write_energy(CS%u, CS%v, CS%h, CS%tv, Time_local, CS%nstep_tot, & G, GV, CS%sum_output_CSp, CS%tracer_flow_CSp, & - dt_forcing=set_time(int(floor(time_interval+0.5))) ) + dt_forcing=real_to_time(time_interval) ) call cpu_clock_end(id_clock_other) @@ -912,7 +912,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & if ((CS%t_dyn_rel_adv == 0.0) .and. CS%thickness_diffuse .and. CS%thickness_diffuse_first) then - call enable_averaging(dt_thermo, Time_local+set_time(int(floor(dt_thermo-dt+0.5))), CS%diag) + call enable_averaging(dt_thermo, Time_local+real_to_time(dt_thermo-dt), CS%diag) call cpu_clock_begin(id_clock_thick_diff) if (associated(CS%VarMix)) & call calc_slope_functions(h, CS%tv, dt, G, GV, CS%VarMix) @@ -931,7 +931,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & ! The bottom boundary layer properties need to be recalculated. if (bbl_time_int > 0.0) then call enable_averaging(bbl_time_int, & - Time_local + set_time(int(bbl_time_int-dt+0.5)), CS%diag) + Time_local + real_to_time(bbl_time_int-dt), CS%diag) ! Calculate the BBL properties and store them inside visc (u,h). call cpu_clock_begin(id_clock_BBL_visc) call set_viscous_BBL(CS%u, CS%v, CS%h, CS%tv, CS%visc, G, GV, & @@ -2135,8 +2135,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & ! Set a few remaining fields that are specific to the ocean grid type. call set_first_direction(G, first_direction) ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes. - if (CS%debug .or. G%symmetric) & + if (CS%debug .or. G%symmetric) then call clone_MOM_domain(G%Domain, G%Domain_aux, symmetric=.false.) + else ; G%Domain_aux => G%Domain ; endif ! Copy common variables from the vertical grid to the horizontal grid. ! Consider removing this later? G%ke = GV%ke ; G%g_Earth = GV%g_Earth @@ -2165,8 +2166,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call MOM_grid_end(G) ; deallocate(G) G => CS%G - if (CS%debug .or. CS%G%symmetric) & + if (CS%debug .or. CS%G%symmetric) then call clone_MOM_domain(CS%G%Domain, CS%G%Domain_aux, symmetric=.false.) + else ; CS%G%Domain_aux => CS%G%Domain ;endif G%ke = GV%ke ; G%g_Earth = GV%g_Earth endif @@ -2286,7 +2288,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & CS%OBC, CS%update_OBC_CSp, CS%ALE_CSp, CS%set_visc_CSp, & CS%visc, dirs, CS%ntrunc, calc_dtbt=calc_dtbt) if (CS%dtbt_reset_period > 0.0) then - CS%dtbt_reset_interval = set_time(int(floor(CS%dtbt_reset_period))) + CS%dtbt_reset_interval = real_to_time(CS%dtbt_reset_period) ! Set dtbt_reset_time to be the next even multiple of dtbt_reset_interval. CS%dtbt_reset_time = Time_init + CS%dtbt_reset_interval * & ((Time - Time_init) / CS%dtbt_reset_interval) @@ -2325,11 +2327,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & param_file, diag, CS%diagnostics_CSp, CS%tv) call diag_copy_diag_to_storage(CS%diag_pre_sync, CS%h, CS%diag) - CS%Z_diag_interval = set_time(int((CS%dt_therm) * & - max(1,floor(0.01 + Z_diag_int/(CS%dt_therm))))) + CS%Z_diag_interval = real_to_time(CS%dt_therm * max(1,floor(0.01 + Z_diag_int/CS%dt_therm))) call MOM_diag_to_Z_init(Time, G, GV, param_file, diag, CS%diag_to_Z_CSp) CS%Z_diag_time = Start_time + CS%Z_diag_interval * (1 + & - ((Time + set_time(int(CS%dt_therm))) - Start_time) / CS%Z_diag_interval) + ((Time + real_to_time(CS%dt_therm)) - Start_time) / CS%Z_diag_interval) if (associated(CS%sponge_CSp)) & call init_sponge_diags(Time, G, diag, CS%sponge_CSp) diff --git a/src/core/MOM_barotropic.F90 b/src/core/MOM_barotropic.F90 index 940c99b8be..674f6f1bff 100644 --- a/src/core/MOM_barotropic.F90 +++ b/src/core/MOM_barotropic.F90 @@ -22,7 +22,7 @@ module MOM_barotropic use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S, OBC_segment_type use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS use MOM_tidal_forcing, only : tidal_forcing_sensitivity, tidal_forcing_CS -use MOM_time_manager, only : time_type, set_time, operator(+), operator(-) +use MOM_time_manager, only : time_type, real_to_time, operator(+), operator(-) use MOM_variables, only : BT_cont_type, alloc_bt_cont_type use MOM_verticalGrid, only : verticalGrid_type @@ -723,7 +723,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, & (CS%id_uhbt_hifreq > 0) .or. (CS%id_vhbt_hifreq > 0)) then do_hifreq_output = query_averaging_enabled(CS%diag, time_int_in, time_end_in) if (do_hifreq_output) & - time_bt_start = time_end_in - set_time(int(floor(dt+0.5))) + time_bt_start = time_end_in - real_to_time(dt) endif !--- begin setup for group halo update @@ -2008,7 +2008,7 @@ subroutine btstep(U_in, V_in, eta_in, dt, bc_accel_u, bc_accel_v, & enddo ; enddo if (do_hifreq_output) then - time_step_end = time_bt_start + set_time(int(floor(n*dtbt+0.5))) + time_step_end = time_bt_start + real_to_time(n*dtbt) call enable_averaging(dtbt, time_step_end, CS%diag) if (CS%id_ubt_hifreq > 0) call post_data(CS%id_ubt_hifreq, ubt(IsdB:IedB,jsd:jed), CS%diag) if (CS%id_vbt_hifreq > 0) call post_data(CS%id_vbt_hifreq, vbt(isd:ied,JsdB:JedB), CS%diag) diff --git a/src/core/MOM_dynamics_split_RK2.F90 b/src/core/MOM_dynamics_split_RK2.F90 index d02285148a..0f4bd88111 100644 --- a/src/core/MOM_dynamics_split_RK2.F90 +++ b/src/core/MOM_dynamics_split_RK2.F90 @@ -30,7 +30,7 @@ module MOM_dynamics_split_RK2 use MOM_io, only : MOM_io_init, vardesc, var_desc use MOM_restart, only : register_restart_field, query_initialized, save_restart use MOM_restart, only : restart_init, is_new_run, MOM_restart_CS -use MOM_time_manager, only : time_type, set_time, time_type_to_real, operator(+) +use MOM_time_manager, only : time_type, time_type_to_real, operator(+) use MOM_time_manager, only : operator(-), operator(>), operator(*), operator(/) use MOM_ALE, only : ALE_CS diff --git a/src/core/MOM_dynamics_unsplit.F90 b/src/core/MOM_dynamics_unsplit.F90 index 47d3510c5a..3965758510 100644 --- a/src/core/MOM_dynamics_unsplit.F90 +++ b/src/core/MOM_dynamics_unsplit.F90 @@ -72,7 +72,7 @@ module MOM_dynamics_unsplit use MOM_io, only : MOM_io_init use MOM_restart, only : register_restart_field, query_initialized, save_restart use MOM_restart, only : restart_init, MOM_restart_CS -use MOM_time_manager, only : time_type, set_time, time_type_to_real, operator(+) +use MOM_time_manager, only : time_type, real_to_time, operator(+) use MOM_time_manager, only : operator(-), operator(>), operator(*), operator(/) use MOM_ALE, only : ALE_CS @@ -267,7 +267,7 @@ subroutine step_MOM_dyn_unsplit(u, v, h, tv, visc, Time_local, dt, forces, & call pass_var(hp, G%Domain, clock=id_clock_pass) call pass_vector(uh, vh, G%Domain, clock=id_clock_pass) - call enable_averaging(0.5*dt,Time_local-set_time(int(0.5*dt)), CS%diag) + call enable_averaging(0.5*dt,Time_local-real_to_time(0.5*dt), CS%diag) ! Here the first half of the thickness fluxes are offered for averaging. if (CS%id_uh > 0) call post_data(CS%id_uh, uh, CS%diag) if (CS%id_vh > 0) call post_data(CS%id_vh, vh, CS%diag) diff --git a/src/core/MOM_dynamics_unsplit_RK2.F90 b/src/core/MOM_dynamics_unsplit_RK2.F90 index a1615ad413..0f6d61905e 100644 --- a/src/core/MOM_dynamics_unsplit_RK2.F90 +++ b/src/core/MOM_dynamics_unsplit_RK2.F90 @@ -70,7 +70,7 @@ module MOM_dynamics_unsplit_RK2 use MOM_io, only : MOM_io_init use MOM_restart, only : register_restart_field, query_initialized, save_restart use MOM_restart, only : restart_init, MOM_restart_CS -use MOM_time_manager, only : time_type, set_time, time_type_to_real, operator(+) +use MOM_time_manager, only : time_type, time_type_to_real, operator(+) use MOM_time_manager, only : operator(-), operator(>), operator(*), operator(/) use MOM_ALE, only : ALE_CS diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 9486967b40..9ac616dac0 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -26,11 +26,13 @@ module MOM_forcing_type public extractFluxes1d, extractFluxes2d, optics_type public MOM_forcing_chksum, MOM_mech_forcing_chksum -public calculateBuoyancyFlux1d, calculateBuoyancyFlux2d, forcing_accumulate +public calculateBuoyancyFlux1d, calculateBuoyancyFlux2d +public forcing_accumulate, fluxes_accumulate public forcing_SinglePointPrint, mech_forcing_diags, forcing_diagnostics public register_forcing_type_diags, allocate_forcing_type, deallocate_forcing_type public copy_common_forcing_fields, allocate_mech_forcing, deallocate_mech_forcing -public set_derived_forcing_fields, copy_back_forcing_fields, set_net_mass_forcing +public set_derived_forcing_fields, copy_back_forcing_fields +public set_net_mass_forcing, get_net_mass_forcing !> Structure that contains pointers to the boundary forcing used to drive the !! liquid ocean simulated by MOM. @@ -40,8 +42,6 @@ module MOM_forcing_type !! MESO_surface_forcing.F90, which is a special case of solo_driver/MOM_surface_forcing.F90. type, public :: forcing - ! Pointers in this module should be initialized to NULL. - ! surface stress components and turbulent velocity scale real, pointer, dimension(:,:) :: & ustar => NULL(), & !< surface friction velocity scale (m/s) @@ -152,11 +152,10 @@ module MOM_forcing_type logical :: fluxes_used = .true. !< If true, all of the heat, salt, and mass !! fluxes have been applied to the ocean. - real :: dt_buoy_accum = -1.0 !< The amount of time over which the buoyancy fluxes + real :: dt_buoy_accum = -1.0 !< The amount of time over which the buoyancy fluxes !! should be applied, in s. If negative, this forcing !! type variable has not yet been inialized. - ! heat capacity real :: C_p !< heat capacity of seawater ( J/(K kg) ). !! C_p is is the same value as in thermovar_ptrs_type. @@ -167,7 +166,7 @@ module MOM_forcing_type !! This is not a convenient convention, but imposed on MOM6 by the coupler. ! For internal error tracking - integer :: num_msg = 0 !< Number of messages issues about excessive SW penetration + integer :: num_msg = 0 !< Number of messages issued about excessive SW penetration integer :: max_msg = 2 !< Maximum number of messages to issue about excessive SW penetration end type forcing @@ -211,6 +210,9 @@ module MOM_forcing_type real, pointer, dimension(:,:) :: & rigidity_ice_u => NULL(), & !< Depth-integrated lateral viscosity of ice shelves or sea ice at u-points (m3/s) rigidity_ice_v => NULL() !< Depth-integrated lateral viscosity of ice shelves or sea ice at v-points (m3/s) + real :: dt_force_accum = -1.0 !< The amount of time over which the mechanical forcing fluxes + !! have been averaged, in s. + logical :: net_mass_src_set = .false. !< If true, an estimate of net_mass_src has been provided. logical :: accumulate_p_surf = .false. !< If true, the surface pressure due to the atmosphere !! and various types of ice needs to be accumulated, and the !! surface pressure explicitly reset to zero at the driver level @@ -296,7 +298,7 @@ module MOM_forcing_type integer :: id_netFWGlobalAdj = -1 integer :: id_netFWGlobalScl = -1 - ! momentum flux amd forcing diagnostic handles + ! momentum flux and forcing diagnostic handles integer :: id_taux = -1 integer :: id_tauy = -1 integer :: id_ustar = -1 @@ -1036,6 +1038,11 @@ subroutine MOM_mech_forcing_chksum(mesg, forces, G, haloshift) haloshift=hshift, symmetric=.true.) if (associated(forces%p_surf)) & call hchksum(forces%p_surf, mesg//" forces%p_surf",G%HI,haloshift=hshift) + if (associated(forces%ustar)) & + call hchksum(forces%ustar, mesg//" forces%ustar",G%HI,haloshift=hshift) + if (associated(forces%rigidity_ice_u) .and. associated(forces%rigidity_ice_v)) & + call uvchksum(mesg//" forces%rigidity_ice_[uv]", forces%rigidity_ice_u, & + forces%rigidity_ice_v, G%HI, haloshift=hshift, symmetric=.true.) end subroutine MOM_mech_forcing_chksum @@ -1741,7 +1748,8 @@ subroutine register_forcing_type_diags(Time, diag, use_temperature, handles, use end subroutine register_forcing_type_diags -!> Accumulate the forcing over time steps +!> Accumulate the forcing over time steps, taking input from a mechanical forcing type +!! and a temporary forcing-flux type. subroutine forcing_accumulate(flux_tmp, forces, fluxes, dt, G, wt2) type(forcing), intent(in) :: flux_tmp !< A temporary structure with current !!thermodynamic forcing fields @@ -1752,6 +1760,26 @@ subroutine forcing_accumulate(flux_tmp, forces, fluxes, dt, G, wt2) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure real, intent(out) :: wt2 !< The relative weight of the new fluxes + ! This subroutine copies mechancal forcing from flux_tmp to fluxes and + ! stores the time-weighted averages of the various buoyancy fluxes in fluxes, + ! and increments the amount of time over which the buoyancy forcing should be + ! applied, all via a call to fluxes accumulate. + + call fluxes_accumulate(flux_tmp, fluxes, dt, G, wt2, forces) + +end subroutine forcing_accumulate + +!> Accumulate the thermodynamic fluxes over time steps +subroutine fluxes_accumulate(flux_tmp, fluxes, dt, G, wt2, forces) + type(forcing), intent(in) :: flux_tmp !< A temporary structure with current + !! thermodynamic forcing fields + type(forcing), intent(inout) :: fluxes !< A structure containing time-averaged + !! thermodynamic forcing fields + real, intent(in) :: dt !< The elapsed time since the last call to this subroutine, in s + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + real, intent(out) :: wt2 !< The relative weight of the new fluxes + type(mech_forcing), optional, intent(in) :: forces !< A structure with the driving mechanical forces + ! This subroutine copies mechancal forcing from flux_tmp to fluxes and ! stores the time-weighted averages of the various buoyancy fluxes in fluxes, ! and increments the amount of time over which the buoyancy forcing should be @@ -1774,15 +1802,29 @@ subroutine forcing_accumulate(flux_tmp, forces, fluxes, dt, G, wt2) wt2 = 1.0 - wt1 ! = dt / (fluxes%dt_buoy_accum + dt) fluxes%dt_buoy_accum = fluxes%dt_buoy_accum + dt - ! Copy over the pressure fields. - do j=js,je ; do i=is,ie - fluxes%p_surf(i,j) = forces%p_surf(i,j) - fluxes%p_surf_full(i,j) = forces%p_surf_full(i,j) - enddo ; enddo + ! Copy over the pressure fields and accumulate averages of ustar, either from the forcing + ! type or from the temporary fluxes type. + if (present(forces)) then + do j=js,je ; do i=is,ie + fluxes%p_surf(i,j) = forces%p_surf(i,j) + fluxes%p_surf_full(i,j) = forces%p_surf_full(i,j) + + fluxes%ustar(i,j) = wt1*fluxes%ustar(i,j) + wt2*forces%ustar(i,j) + enddo ; enddo + else + do j=js,je ; do i=is,ie + fluxes%p_surf(i,j) = flux_tmp%p_surf(i,j) + fluxes%p_surf_full(i,j) = flux_tmp%p_surf_full(i,j) + + fluxes%ustar(i,j) = wt1*fluxes%ustar(i,j) + wt2*flux_tmp%ustar(i,j) + enddo ; enddo + endif ! Average the water, heat, and salt fluxes, and ustar. do j=js,je ; do i=is,ie - fluxes%ustar(i,j) = wt1*fluxes%ustar(i,j) + wt2*forces%ustar(i,j) +!### Replace the expression for ustar_gustless with this one... +! fluxes%ustar_gustless(i,j) = wt1*fluxes%ustar_gustless(i,j) + wt2*flux_tmp%ustar_gustless(i,j) + fluxes%ustar_gustless(i,j) = flux_tmp%ustar_gustless(i,j) fluxes%evap(i,j) = wt1*fluxes%evap(i,j) + wt2*flux_tmp%evap(i,j) fluxes%lprec(i,j) = wt1*fluxes%lprec(i,j) + wt2*flux_tmp%lprec(i,j) @@ -1866,7 +1908,7 @@ subroutine forcing_accumulate(flux_tmp, forces, fluxes, dt, G, wt2) call coupler_type_increment_data(flux_tmp%tr_fluxes, fluxes%tr_fluxes, & scale_factor=wt2, scale_prev=wt1) -end subroutine forcing_accumulate +end subroutine fluxes_accumulate !> This subroutine copies the computational domains of common forcing fields !! from a mech_forcing type to a (thermodynamic) forcing type. @@ -1922,9 +1964,12 @@ subroutine set_derived_forcing_fields(forces, fluxes, G, Rho0) !! as used to calculate ustar. real :: taux2, tauy2 ! Squared wind stress components, in Pa^2. + real :: Irho0 ! Inverse of the mean density in (m^3/kg) integer :: i, j, is, ie, js, je is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + Irho0 = 1.0/Rho0 + if (associated(forces%taux) .and. associated(forces%tauy) .and. & associated(fluxes%ustar_gustless)) then do j=js,je ; do i=is,ie @@ -1940,45 +1985,58 @@ subroutine set_derived_forcing_fields(forces, fluxes, G, Rho0) (G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) fluxes%ustar_gustless(i,j) = sqrt(sqrt(taux2 + tauy2) / Rho0) +!### Change to: +! fluxes%ustar_gustless(i,j) = sqrt(sqrt(taux2 + tauy2) * Irho0) enddo ; enddo endif end subroutine set_derived_forcing_fields -!> This subroutine calculates determines the net mass source to th eocean from +!> This subroutine calculates determines the net mass source to the ocean from !! a (thermodynamic) forcing type and stores it in a mech_forcing type. subroutine set_net_mass_forcing(fluxes, forces, G) type(forcing), intent(in) :: fluxes !< A structure containing thermodynamic forcing fields type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces - type(ocean_grid_type), intent(in) :: G !< grid type + type(ocean_grid_type), intent(in) :: G !< The ocean grid type + + if (associated(forces%net_mass_src)) & + call get_net_mass_forcing(fluxes, G, forces%net_mass_src) + +end subroutine set_net_mass_forcing + +!> This subroutine calculates determines the net mass source to the ocean from +!! a (thermodynamic) forcing type and stores it in a provided array. +subroutine get_net_mass_forcing(fluxes, G, net_mass_src) + type(forcing), intent(in) :: fluxes !< A structure containing thermodynamic forcing fields + type(ocean_grid_type), intent(in) :: G !< The ocean grid type + real, dimension(SZI_(G),SZJ_(G)), intent(out) :: net_mass_src !< The net mass flux of water into the ocean + !! in kg m-2 s-1. integer :: i, j, is, ie, js, je is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec - if (associated(forces%net_mass_src)) then - forces%net_mass_src(:,:) = 0.0 - if (associated(fluxes%lprec)) then ; do j=js,je ; do i=is,ie - forces%net_mass_src(i,j) = forces%net_mass_src(i,j) + fluxes%lprec(i,j) - enddo ; enddo ; endif - if (associated(fluxes%fprec)) then ; do j=js,je ; do i=is,ie - forces%net_mass_src(i,j) = forces%net_mass_src(i,j) + fluxes%fprec(i,j) - enddo ; enddo ; endif - if (associated(fluxes%vprec)) then ; do j=js,je ; do i=is,ie - forces%net_mass_src(i,j) = forces%net_mass_src(i,j) + fluxes%vprec(i,j) - enddo ; enddo ; endif - if (associated(fluxes%lrunoff)) then ; do j=js,je ; do i=is,ie - forces%net_mass_src(i,j) = forces%net_mass_src(i,j) + fluxes%lrunoff(i,j) - enddo ; enddo ; endif - if (associated(fluxes%frunoff)) then ; do j=js,je ; do i=is,ie - forces%net_mass_src(i,j) = forces%net_mass_src(i,j) + fluxes%frunoff(i,j) - enddo ; enddo ; endif - if (associated(fluxes%evap)) then ; do j=js,je ; do i=is,ie - forces%net_mass_src(i,j) = forces%net_mass_src(i,j) + fluxes%evap(i,j) - enddo ; enddo ; endif - endif - -end subroutine set_net_mass_forcing + net_mass_src(:,:) = 0.0 + if (associated(fluxes%lprec)) then ; do j=js,je ; do i=is,ie + net_mass_src(i,j) = net_mass_src(i,j) + fluxes%lprec(i,j) + enddo ; enddo ; endif + if (associated(fluxes%fprec)) then ; do j=js,je ; do i=is,ie + net_mass_src(i,j) = net_mass_src(i,j) + fluxes%fprec(i,j) + enddo ; enddo ; endif + if (associated(fluxes%vprec)) then ; do j=js,je ; do i=is,ie + net_mass_src(i,j) = net_mass_src(i,j) + fluxes%vprec(i,j) + enddo ; enddo ; endif + if (associated(fluxes%lrunoff)) then ; do j=js,je ; do i=is,ie + net_mass_src(i,j) = net_mass_src(i,j) + fluxes%lrunoff(i,j) + enddo ; enddo ; endif + if (associated(fluxes%frunoff)) then ; do j=js,je ; do i=is,ie + net_mass_src(i,j) = net_mass_src(i,j) + fluxes%frunoff(i,j) + enddo ; enddo ; endif + if (associated(fluxes%evap)) then ; do j=js,je ; do i=is,ie + net_mass_src(i,j) = net_mass_src(i,j) + fluxes%evap(i,j) + enddo ; enddo ; endif + +end subroutine get_net_mass_forcing !> This subroutine copies the computational domains of common forcing fields !! from a mech_forcing type to a (thermodynamic) forcing type. @@ -2001,9 +2059,8 @@ end subroutine copy_back_forcing_fields !> Offer mechanical forcing fields for diagnostics for those !! fields registered as part of register_forcing_type_diags. -subroutine mech_forcing_diags(forces, fluxes, dt, G, diag, handles) +subroutine mech_forcing_diags(forces, dt, G, diag, handles) type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces - type(forcing), intent(in) :: fluxes !< A structure containing thermodynamic forcing fields real, intent(in) :: dt !< time step type(ocean_grid_type), intent(in) :: G !< grid type type(diag_ctrl), intent(in) :: diag !< diagnostic type @@ -2018,20 +2075,15 @@ subroutine mech_forcing_diags(forces, fluxes, dt, G, diag, handles) if ((handles%id_taux > 0) .and. associated(forces%taux)) & call post_data(handles%id_taux, forces%taux, diag) + if ((handles%id_tauy > 0) .and. associated(forces%tauy)) & call post_data(handles%id_tauy, forces%tauy, diag) - if ((handles%id_ustar > 0) .and. associated(fluxes%ustar)) & - call post_data(handles%id_ustar, fluxes%ustar, diag) - if (handles%id_ustar_berg > 0) & - call post_data(handles%id_ustar_berg, fluxes%ustar_berg, diag) - if (handles%id_area_berg > 0) & - call post_data(handles%id_area_berg, fluxes%area_berg, diag) - if (handles%id_mass_berg > 0) & - call post_data(handles%id_mass_berg, fluxes%mass_berg, diag) - if (handles%id_frac_ice_cover > 0) & - call post_data(handles%id_frac_ice_cover, fluxes%frac_shelf_h, diag) - if (handles%id_ustar_ice_cover > 0) & - call post_data(handles%id_ustar_ice_cover, fluxes%ustar_shelf, diag) + + if ((handles%id_mass_berg > 0) .and. associated(forces%mass_berg)) & + call post_data(handles%id_mass_berg, forces%mass_berg, diag) + + if ((handles%id_area_berg > 0) .and. associated(forces%area_berg)) & + call post_data(handles%id_area_berg, forces%area_berg, diag) endif @@ -2522,8 +2574,19 @@ subroutine forcing_diagnostics(fluxes, sfc_state, dt, G, diag, handles) if ((handles%id_buoy > 0) .and. associated(fluxes%buoy)) & call post_data(handles%id_buoy, fluxes%buoy, diag) + if ((handles%id_ustar > 0) .and. associated(fluxes%ustar)) & + call post_data(handles%id_ustar, fluxes%ustar, diag) - endif + if ((handles%id_ustar_berg > 0) .and. associated(fluxes%ustar_berg)) & + call post_data(handles%id_ustar_berg, fluxes%ustar_berg, diag) + + if ((handles%id_frac_ice_cover > 0) .and. associated(fluxes%frac_shelf_h)) & + call post_data(handles%id_frac_ice_cover, fluxes%frac_shelf_h, diag) + + if ((handles%id_ustar_ice_cover > 0) .and. associated(fluxes%ustar_shelf)) & + call post_data(handles%id_ustar_ice_cover, fluxes%ustar_shelf, diag) + + endif ! query_averaging_enabled call cpu_clock_end(handles%id_clock_forcing) end subroutine forcing_diagnostics diff --git a/src/core/MOM_open_boundary.F90 b/src/core/MOM_open_boundary.F90 index 259714e984..de8c2fe174 100644 --- a/src/core/MOM_open_boundary.F90 +++ b/src/core/MOM_open_boundary.F90 @@ -112,6 +112,10 @@ module MOM_open_boundary logical :: radiation_grad !< If true, 1D Orlanksi radiation boundary conditions are applied to !! dudv and dvdx. logical :: oblique !< Oblique waves supported at radiation boundary. + logical :: oblique_tan !< If true, 2D radiation boundary conditions are applied to + !! tangential flows. + logical :: oblique_grad !< If true, 2D radiation boundary conditions are applied to + !! dudv and dvdx. logical :: nudged !< Optional supplement to radiation boundary. logical :: nudged_tan !< Optional supplement to nudge tangential velocity. logical :: nudged_grad !< Optional supplement to nudge normal gradient of tangential velocity. @@ -151,7 +155,11 @@ module MOM_open_boundary !! the OB segment (m s-1). real, pointer, dimension(:,:) :: eta=>NULL() !< The sea-surface elevation along the segment (m). real, pointer, dimension(:,:,:) :: grad_normal=>NULL() !< The gradient of the normal flow along the - !! segment (m s-1) + !! segment (s-1) + real, pointer, dimension(:,:,:) :: grad_tan=>NULL() !< The gradient of the tangential flow along the + !! segment (s-1) + real, pointer, dimension(:,:,:) :: grad_gradient=>NULL() !< The gradient of the gradient of tangential flow along the + !! segment (m-1 s-1) real, pointer, dimension(:,:,:) :: rx_normal=>NULL() !< The rx_old_u value for radiation coeff !! for normal velocity real, pointer, dimension(:,:,:) :: ry_normal=>NULL() !< The tangential value for radiation coeff @@ -397,6 +405,8 @@ subroutine open_boundary_config(G, param_file, OBC) OBC%segment(l)%radiation_tan = .false. OBC%segment(l)%radiation_grad = .false. OBC%segment(l)%oblique = .false. + OBC%segment(l)%oblique_tan = .false. + OBC%segment(l)%oblique_grad = .false. OBC%segment(l)%nudged = .false. OBC%segment(l)%nudged_tan = .false. OBC%segment(l)%nudged_grad = .false. @@ -818,6 +828,13 @@ subroutine setup_u_point_obc(OBC, G, segment_str, l_seg, PF, reentrant_y) OBC%segment(l_seg)%open = .true. OBC%oblique_BCs_exist_globally = .true. OBC%open_u_BCs_exist_globally = .true. + elseif (trim(action_str(a_loop)) == 'OBLIQUE_TAN') then + OBC%segment(l_seg)%oblique = .true. + OBC%segment(l_seg)%oblique_tan = .true. + OBC%oblique_BCs_exist_globally = .true. + elseif (trim(action_str(a_loop)) == 'OBLIQUE_GRAD') then + OBC%segment(l_seg)%oblique = .true. + OBC%segment(l_seg)%oblique_grad = .true. elseif (trim(action_str(a_loop)) == 'NUDGED') then OBC%segment(l_seg)%nudged = .true. OBC%nudged_u_BCs_exist_globally = .true. @@ -871,6 +888,10 @@ subroutine setup_u_point_obc(OBC, G, segment_str, l_seg, PF, reentrant_y) OBC%segment(l_seg)%Je_obc = Je_obc call allocate_OBC_segment_data(OBC, OBC%segment(l_seg)) + if (OBC%segment(l_seg)%oblique .and. OBC%segment(l_seg)%radiation) & + call MOM_error(FATAL, "MOM_open_boundary.F90, setup_u_point_obc:\n"//& + "Orlanski and Oblique OBC options cannot be used together on one segment.") + end subroutine setup_u_point_obc !> Parse an OBC_SEGMENT_%%% string starting with "J=" and configure placement and type of OBC accordingly @@ -931,6 +952,13 @@ subroutine setup_v_point_obc(OBC, G, segment_str, l_seg, PF, reentrant_x) OBC%segment(l_seg)%open = .true. OBC%oblique_BCs_exist_globally = .true. OBC%open_v_BCs_exist_globally = .true. + elseif (trim(action_str(a_loop)) == 'OBLIQUE_TAN') then + OBC%segment(l_seg)%oblique = .true. + OBC%segment(l_seg)%oblique_tan = .true. + OBC%oblique_BCs_exist_globally = .true. + elseif (trim(action_str(a_loop)) == 'OBLIQUE_GRAD') then + OBC%segment(l_seg)%oblique = .true. + OBC%segment(l_seg)%oblique_grad = .true. elseif (trim(action_str(a_loop)) == 'NUDGED') then OBC%segment(l_seg)%nudged = .true. OBC%nudged_v_BCs_exist_globally = .true. @@ -984,6 +1012,10 @@ subroutine setup_v_point_obc(OBC, G, segment_str, l_seg, PF, reentrant_x) OBC%segment(l_seg)%Je_obc = J_obc call allocate_OBC_segment_data(OBC, OBC%segment(l_seg)) + if (OBC%segment(l_seg)%oblique .and. OBC%segment(l_seg)%radiation) & + call MOM_error(FATAL, "MOM_open_boundary.F90, setup_v_point_obc:\n"//& + "Orlanski and Oblique OBC options cannot be used together on one segment.") + end subroutine setup_v_point_obc !> Parse an OBC_SEGMENT_%%% string @@ -1505,6 +1537,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) real :: cff_new, cff_avg ! denominator in oblique real, pointer, dimension(:,:,:) :: rx_tangential=>NULL() real, pointer, dimension(:,:,:) :: ry_tangential=>NULL() + real, pointer, dimension(:,:,:) :: cff_tangential=>NULL() real, parameter :: eps = 1.0e-20 type(OBC_segment_type), pointer :: segment => NULL() integer :: i, j, k, is, ie, js, je, nz, n @@ -1606,6 +1639,8 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) segment%normal_vel(I,j,k) = ((cff_avg*u_new(I,j,k) + rx_avg*u_new(I-1,j,k)) - & (max(ry_avg,0.0)*segment%grad_normal(J-1,2,k) + min(ry_avg,0.0)*segment%grad_normal(J,2,k))) / & (cff_avg + rx_avg) + ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues + ! implemented as a work-around to limitations in restart capability OBC%rx_normal(I,j,k) = segment%rx_normal(I,j,k) OBC%ry_normal(i,J,k) = segment%ry_normal(i,J,k) OBC%cff_normal(I,j,k) = segment%cff_normal(I,j,k) @@ -1642,7 +1677,8 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) endif if (segment%nudged_tan) then do k=1,nz ; do J=segment%HI%JsdB,segment%HI%JedB - if (rx_tangential(I,J,k) < 0.0) then + ! dhdt gets set to 0 on inflow in oblique case + if (rx_tangential(I,J,k) <= 0.0) then tau = segment%Velocity_nudging_timescale_in else tau = segment%Velocity_nudging_timescale_out @@ -1672,7 +1708,77 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) endif if (segment%nudged_grad) then do k=1,nz ; do J=segment%HI%JsdB,segment%HI%JedB - if (rx_tangential(I,J,k) < 0.0) then + ! dhdt gets set to 0 on inflow in oblique case + if (rx_tangential(I,J,k) <= 0.0) then + tau = segment%Velocity_nudging_timescale_in + else + tau = segment%Velocity_nudging_timescale_out + endif + gamma_2 = dt / (tau + dt) + segment%tangential_grad(I,J,k) = (1 - gamma_2) * segment%tangential_grad(I,J,k) + & + gamma_2 * segment%nudged_tangential_grad(I,J,k) + enddo ; enddo + endif + deallocate(rx_tangential) + endif + if (segment%oblique_tan .or. segment%oblique_grad) then + I=segment%HI%IsdB + allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz)) + allocate(ry_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz)) + allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz)) + do k=1,nz + rx_tangential(I,segment%HI%JsdB,k) = segment%rx_normal(I,segment%HI%jsd,k) + rx_tangential(I,segment%HI%JedB,k) = segment%rx_normal(I,segment%HI%jed,k) + ry_tangential(I,segment%HI%JsdB,k) = segment%ry_normal(I,segment%HI%jsd,k) + ry_tangential(I,segment%HI%JedB,k) = segment%ry_normal(I,segment%HI%jed,k) + cff_tangential(I,segment%HI%JsdB,k) = segment%cff_normal(I,segment%HI%jsd,k) + cff_tangential(I,segment%HI%JedB,k) = segment%cff_normal(I,segment%HI%jed,k) + do J=segment%HI%JsdB+1,segment%HI%JedB-1 + rx_tangential(I,J,k) = 0.5*(segment%rx_normal(I,j,k) + segment%rx_normal(I,j+1,k)) + ry_tangential(I,J,k) = 0.5*(segment%ry_normal(I,j,k) + segment%ry_normal(I,j+1,k)) + cff_tangential(I,J,k) = 0.5*(segment%cff_normal(I,j,k) + segment%cff_normal(I,j+1,k)) + enddo + enddo + if (segment%oblique_tan) then + do k=1,nz ; do J=segment%HI%JsdB+1,segment%HI%JedB-1 + rx_avg = rx_tangential(I,J,k) + ry_avg = ry_tangential(I,J,k) + cff_avg = cff_tangential(I,J,k) + segment%tangential_vel(I,J,k) = ((cff_avg*v_new(i,J,k) + rx_avg*v_new(i-1,J,k)) - & + (max(ry_avg,0.0)*segment%grad_tan(j,2,k) + min(ry_avg,0.0)*segment%grad_tan(j+1,2,k))) / & + (cff_avg + rx_avg) + enddo ; enddo + endif + if (segment%nudged_tan) then + do k=1,nz ; do J=segment%HI%JsdB,segment%HI%JedB + ! dhdt gets set to 0 on inflow in oblique case + if (rx_tangential(I,J,k) <= 0.0) then + tau = segment%Velocity_nudging_timescale_in + else + tau = segment%Velocity_nudging_timescale_out + endif + gamma_2 = dt / (tau + dt) + segment%tangential_vel(I,J,k) = (1 - gamma_2) * segment%tangential_vel(I,J,k) + & + gamma_2 * segment%nudged_tangential_vel(I,J,k) + enddo ; enddo + endif + if (segment%oblique_grad) then + Js_obc = max(segment%HI%JsdB,G%jsd+1) + Je_obc = min(segment%HI%JedB,G%jed-1) + do k=1,nz ; do J=segment%HI%JsdB+1,segment%HI%JedB-1 + rx_avg = rx_tangential(I,J,k) + ry_avg = ry_tangential(I,J,k) + cff_avg = cff_tangential(I,J,k) + segment%tangential_grad(I,J,k) = ((cff_avg*(v_new(i,J,k) - v_new(i-1,J,k))*G%IdxBu(I-1,J) & + + rx_avg*(v_new(i-1,J,k) - v_new(i-2,J,k))*G%IdxBu(I-2,J)) - & + (max(ry_avg,0.0)*segment%grad_gradient(J,2,k) + min(ry_avg,0.0)*segment%grad_gradient(J+1,2,k))) / & + (cff_avg + rx_avg) + enddo ; enddo + endif + if (segment%nudged_grad) then + do k=1,nz ; do J=segment%HI%JsdB,segment%HI%JedB + ! dhdt gets set to 0 on inflow in oblique case + if (rx_tangential(I,J,k) <= 0.0) then tau = segment%Velocity_nudging_timescale_in else tau = segment%Velocity_nudging_timescale_out @@ -1683,6 +1789,8 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) enddo ; enddo endif deallocate(rx_tangential) + deallocate(ry_tangential) + deallocate(cff_tangential) endif endif @@ -1727,6 +1835,8 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) segment%normal_vel(I,j,k) = ((cff_avg*u_new(I,j,k) + rx_avg*u_new(I+1,j,k)) - & (max(ry_avg,0.0)*segment%grad_normal(J-1,2,k) + min(ry_avg,0.0)*segment%grad_normal(J,2,k))) / & (cff_avg + rx_avg) + ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues + ! implemented as a work-around to limitations in restart capability OBC%rx_normal(I,j,k) = segment%rx_normal(I,j,k) OBC%ry_normal(i,J,k) = segment%ry_normal(i,J,k) OBC%cff_normal(I,j,k) = segment%cff_normal(I,j,k) @@ -1758,12 +1868,13 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) if (segment%radiation_tan) then do k=1,nz ; do J=segment%HI%JsdB,segment%HI%JedB rx_avg = rx_tangential(I,J,k) - segment%tangential_vel(I,J,k) = (v_new(I+1,J,k) + rx_avg*v_new(I+2,J,k)) / (1.0+rx_avg) + segment%tangential_vel(I,J,k) = (v_new(i+1,J,k) + rx_avg*v_new(i+2,J,k)) / (1.0+rx_avg) enddo ; enddo endif if (segment%nudged_tan) then do k=1,nz ; do J=segment%HI%JsdB,segment%HI%JedB - if (rx_tangential(I,J,k) < 0.0) then + ! dhdt gets set to 0 on inflow in oblique case + if (rx_tangential(I,J,k) <= 0.0) then tau = segment%Velocity_nudging_timescale_in else tau = segment%Velocity_nudging_timescale_out @@ -1793,7 +1904,8 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) endif if (segment%nudged_grad) then do k=1,nz ; do J=segment%HI%JsdB,segment%HI%JedB - if (rx_tangential(I,J,k) < 0.0) then + ! dhdt gets set to 0 on inflow in oblique case + if (rx_tangential(I,J,k) <= 0.0) then tau = segment%Velocity_nudging_timescale_in else tau = segment%Velocity_nudging_timescale_out @@ -1805,6 +1917,77 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) endif deallocate(rx_tangential) endif + if (segment%oblique_tan .or. segment%oblique_grad) then + I=segment%HI%IsdB + allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz)) + allocate(ry_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz)) + allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz)) + do k=1,nz + rx_tangential(I,segment%HI%JsdB,k) = segment%rx_normal(I,segment%HI%jsd,k) + rx_tangential(I,segment%HI%JedB,k) = segment%rx_normal(I,segment%HI%jed,k) + ry_tangential(I,segment%HI%JsdB,k) = segment%ry_normal(I,segment%HI%jsd,k) + ry_tangential(I,segment%HI%JedB,k) = segment%ry_normal(I,segment%HI%jed,k) + cff_tangential(I,segment%HI%JsdB,k) = segment%cff_normal(I,segment%HI%jsd,k) + cff_tangential(I,segment%HI%JedB,k) = segment%cff_normal(I,segment%HI%jed,k) + do J=segment%HI%JsdB+1,segment%HI%JedB-1 + rx_tangential(I,J,k) = 0.5*(segment%rx_normal(I,j,k) + segment%rx_normal(I,j+1,k)) + ry_tangential(I,J,k) = 0.5*(segment%ry_normal(I,j,k) + segment%ry_normal(I,j+1,k)) + cff_tangential(I,J,k) = 0.5*(segment%cff_normal(I,j,k) + segment%cff_normal(I,j+1,k)) + enddo + enddo + if (segment%oblique_tan) then + do k=1,nz ; do J=segment%HI%JsdB+1,segment%HI%JedB-1 + rx_avg = rx_tangential(I,J,k) + ry_avg = ry_tangential(I,J,k) + cff_avg = cff_tangential(I,J,k) + segment%tangential_vel(I,J,k) = ((cff_avg*v_new(i+1,J,k) + rx_avg*v_new(i+2,J,k)) - & + (max(ry_avg,0.0)*segment%grad_tan(j,2,k) + min(ry_avg,0.0)*segment%grad_tan(j+1,2,k))) / & + (cff_avg + rx_avg) + enddo ; enddo + endif + if (segment%nudged_tan) then + do k=1,nz ; do J=segment%HI%JsdB,segment%HI%JedB + ! dhdt gets set to 0 on inflow in oblique case + if (rx_tangential(I,J,k) <= 0.0) then + tau = segment%Velocity_nudging_timescale_in + else + tau = segment%Velocity_nudging_timescale_out + endif + gamma_2 = dt / (tau + dt) + segment%tangential_vel(I,J,k) = (1 - gamma_2) * segment%tangential_vel(I,J,k) + & + gamma_2 * segment%nudged_tangential_vel(I,J,k) + enddo ; enddo + endif + if (segment%oblique_grad) then + Js_obc = max(segment%HI%JsdB,G%jsd+1) + Je_obc = min(segment%HI%JedB,G%jed-1) + do k=1,nz ; do J=segment%HI%JsdB+1,segment%HI%JedB-1 + rx_avg = rx_tangential(I,J,k) + ry_avg = ry_tangential(I,J,k) + cff_avg = cff_tangential(I,J,k) + segment%tangential_grad(I,J,k) = ((cff_avg*(v_new(i+2,J,k) - v_new(i+1,J,k))*G%IdxBu(I+1,J) & + + rx_avg*(v_new(i+3,J,k) - v_new(i+2,J,k))*G%IdxBu(I+2,J)) - & + (max(ry_avg,0.0)*segment%grad_gradient(J,2,k) + min(ry_avg,0.0)*segment%grad_gradient(J+1,2,k))) / & + (cff_avg + rx_avg) + enddo ; enddo + endif + if (segment%nudged_grad) then + do k=1,nz ; do J=segment%HI%JsdB,segment%HI%JedB + ! dhdt gets set to 0 on inflow in oblique case + if (rx_tangential(I,J,k) <= 0.0) then + tau = segment%Velocity_nudging_timescale_in + else + tau = segment%Velocity_nudging_timescale_out + endif + gamma_2 = dt / (tau + dt) + segment%tangential_grad(I,J,k) = (1 - gamma_2) * segment%tangential_grad(I,J,k) + & + gamma_2 * segment%nudged_tangential_grad(I,J,k) + enddo ; enddo + endif + deallocate(rx_tangential) + deallocate(ry_tangential) + deallocate(cff_tangential) + endif endif if (segment%direction == OBC_DIRECTION_N) then @@ -1849,6 +2032,8 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) segment%normal_vel(i,J,k) = ((cff_avg*v_new(i,J,k) + ry_avg*v_new(i,J-1,k)) - & (max(rx_avg,0.0)*segment%grad_normal(I-1,2,k) + min(rx_avg,0.0)*segment%grad_normal(I,2,k))) / & (cff_avg + ry_avg) + ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues + ! implemented as a work-around to limitations in restart capability OBC%rx_normal(I,j,k) = segment%rx_normal(I,j,k) OBC%ry_normal(i,J,k) = segment%ry_normal(i,J,k) OBC%cff_normal(i,J,k) = segment%cff_normal(i,J,k) @@ -1857,7 +2042,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) endif if ((segment%radiation .or. segment%oblique) .and. segment%nudged) then ! dhdt gets set to 0 on inflow in oblique case - if (dhdt*dhdx <= 0.0) then + if (dhdt*dhdy <= 0.0) then tau = segment%Velocity_nudging_timescale_in else tau = segment%Velocity_nudging_timescale_out @@ -1885,7 +2070,8 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) endif if (segment%nudged_tan) then do k=1,nz ; do I=segment%HI%IsdB,segment%HI%IedB - if (rx_tangential(I,J,k) < 0.0) then + ! dhdt gets set to 0 on inflow in oblique case + if (rx_tangential(I,J,k) <= 0.0) then tau = segment%Velocity_nudging_timescale_in else tau = segment%Velocity_nudging_timescale_out @@ -1909,13 +2095,14 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) ! else ! rx_avg = 0.0 ! endif - segment%tangential_grad(I,J,k) = ((u_new(I,j,k) - u_new(I-1,j,k))*G%IdyBu(I,J-1) + & + segment%tangential_grad(I,J,k) = ((u_new(I,j,k) - u_new(I,j-1,k))*G%IdyBu(I,J-1) + & rx_avg*(u_new(I,j-1,k) - u_new(I,j-2,k))*G%IdyBu(I,J-2)) / (1.0+rx_avg) enddo ; enddo endif if (segment%nudged_grad) then do k=1,nz ; do I=segment%HI%IsdB,segment%HI%IedB - if (rx_tangential(I,J,k) < 0.0) then + ! dhdt gets set to 0 on inflow in oblique case + if (ry_tangential(I,J,k) <= 0.0) then tau = segment%Velocity_nudging_timescale_in else tau = segment%Velocity_nudging_timescale_out @@ -1927,9 +2114,79 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) endif deallocate(rx_tangential) endif + if (segment%oblique_tan .or. segment%oblique_grad) then + J=segment%HI%JsdB + allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz)) + allocate(ry_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz)) + allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz)) + do k=1,nz + rx_tangential(segment%HI%IsdB,J,k) = segment%ry_normal(segment%HI%isd,J,k) + rx_tangential(segment%HI%IedB,J,k) = segment%ry_normal(segment%HI%ied,J,k) + ry_tangential(segment%HI%IsdB,J,k) = segment%rx_normal(segment%HI%isd,J,k) + ry_tangential(segment%HI%IedB,J,k) = segment%rx_normal(segment%HI%ied,J,k) + cff_tangential(segment%HI%IsdB,J,k) = segment%cff_normal(segment%HI%isd,J,k) + cff_tangential(segment%HI%IedB,J,k) = segment%cff_normal(segment%HI%ied,J,k) + do I=segment%HI%IsdB+1,segment%HI%IedB-1 + rx_tangential(I,J,k) = 0.5*(segment%ry_normal(i,J,k) + segment%ry_normal(i+1,J,k)) + ry_tangential(I,J,k) = 0.5*(segment%rx_normal(i,J,k) + segment%rx_normal(i+1,J,k)) + cff_tangential(I,J,k) = 0.5*(segment%cff_normal(i,J,k) + segment%cff_normal(i+1,J,k)) + enddo + enddo + if (segment%oblique_tan) then + do k=1,nz ; do I=segment%HI%IsdB+1,segment%HI%IedB-1 + rx_avg = rx_tangential(I,J,k) + ry_avg = ry_tangential(I,J,k) + cff_avg = cff_tangential(I,J,k) + segment%tangential_vel(I,J,k) = ((cff_avg*u_new(I,j,k) + rx_avg*u_new(I,j-1,k)) - & + (max(ry_avg,0.0)*segment%grad_tan(i,2,k) + min(ry_avg,0.0)*segment%grad_tan(i+1,2,k))) / & + (cff_avg + rx_avg) + enddo ; enddo + endif + if (segment%nudged_tan) then + do k=1,nz ; do I=segment%HI%IsdB,segment%HI%IedB + ! dhdt gets set to 0 on inflow in oblique case + if (ry_tangential(I,J,k) <= 0.0) then + tau = segment%Velocity_nudging_timescale_in + else + tau = segment%Velocity_nudging_timescale_out + endif + gamma_2 = dt / (tau + dt) + segment%tangential_vel(I,J,k) = (1 - gamma_2) * segment%tangential_vel(I,J,k) + & + gamma_2 * segment%nudged_tangential_vel(I,J,k) + enddo ; enddo + endif + if (segment%oblique_grad) then + Is_obc = max(segment%HI%IsdB,G%isd+1) + Ie_obc = min(segment%HI%IedB,G%ied-1) + do k=1,nz ; do I=segment%HI%IsdB+1,segment%HI%IedB-1 + rx_avg = rx_tangential(I,J,k) + ry_avg = ry_tangential(I,J,k) + cff_avg = cff_tangential(I,J,k) + segment%tangential_grad(I,J,k) = ((cff_avg*(u_new(I,j,k) - u_new(I,j-1,k))*G%IdyBu(I,J-1) & + + rx_avg*(u_new(I,j-1,k) - u_new(I,j-2,k))*G%IdyBu(I,J-2)) - & + (max(ry_avg,0.0)*segment%grad_gradient(I,2,k) + min(ry_avg,0.0)*segment%grad_gradient(I+1,2,k))) / & + (cff_avg + rx_avg) + enddo ; enddo + endif + if (segment%nudged_grad) then + do k=1,nz ; do I=segment%HI%IsdB,segment%HI%IedB + ! dhdt gets set to 0 on inflow in oblique case + if (ry_tangential(I,J,k) <= 0.0) then + tau = segment%Velocity_nudging_timescale_in + else + tau = segment%Velocity_nudging_timescale_out + endif + gamma_2 = dt / (tau + dt) + segment%tangential_grad(I,J,k) = (1 - gamma_2) * segment%tangential_grad(I,J,k) + & + gamma_2 * segment%nudged_tangential_grad(I,J,k) + enddo ; enddo + endif + deallocate(rx_tangential) + deallocate(ry_tangential) + deallocate(cff_tangential) + endif endif - if (segment%direction == OBC_DIRECTION_S) then J=segment%HI%JsdB if (J>G%HI%JecB) cycle @@ -1971,6 +2228,8 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) segment%normal_vel(i,J,k) = ((cff_avg*v_new(i,J,k) + ry_avg*v_new(i,J+1,k)) - & (max(rx_avg,0.0)*segment%grad_normal(I-1,2,k) + min(rx_avg,0.0)*segment%grad_normal(I,2,k))) / & (cff_avg + ry_avg) + ! Copy restart fields into 3-d arrays. This is an inefficient and temporary issues + ! implemented as a work-around to limitations in restart capability OBC%rx_normal(I,j,k) = segment%rx_normal(I,j,k) OBC%ry_normal(i,J,k) = segment%ry_normal(i,J,k) OBC%cff_normal(i,J,k) = segment%cff_normal(i,J,k) @@ -1979,7 +2238,7 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) endif if ((segment%radiation .or. segment%oblique) .and. segment%nudged) then ! dhdt gets set to 0 on inflow in oblique case - if (dhdt*dhdx <= 0.0) then + if (dhdt*dhdy <= 0.0) then tau = segment%Velocity_nudging_timescale_in else tau = segment%Velocity_nudging_timescale_out @@ -2007,7 +2266,8 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) endif if (segment%nudged_tan) then do k=1,nz ; do I=segment%HI%IsdB,segment%HI%IedB - if (rx_tangential(I,J,k) < 0.0) then + ! dhdt gets set to 0 on inflow in oblique case + if (rx_tangential(I,J,k) <= 0.0) then tau = segment%Velocity_nudging_timescale_in else tau = segment%Velocity_nudging_timescale_out @@ -2037,7 +2297,77 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) endif if (segment%nudged_grad) then do k=1,nz ; do J=segment%HI%JsdB,segment%HI%JedB - if (rx_tangential(I,J,k) < 0.0) then + ! dhdt gets set to 0 on inflow in oblique case + if (ry_tangential(I,J,k) <= 0.0) then + tau = segment%Velocity_nudging_timescale_in + else + tau = segment%Velocity_nudging_timescale_out + endif + gamma_2 = dt / (tau + dt) + segment%tangential_grad(I,J,k) = (1 - gamma_2) * segment%tangential_grad(I,J,k) + & + gamma_2 * segment%nudged_tangential_grad(I,J,k) + enddo ; enddo + endif + deallocate(rx_tangential) + endif + if (segment%oblique_tan .or. segment%oblique_grad) then + J=segment%HI%JsdB + allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz)) + allocate(ry_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz)) + allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz)) + do k=1,nz + rx_tangential(segment%HI%IsdB,J,k) = segment%ry_normal(segment%HI%isd,J,k) + rx_tangential(segment%HI%IedB,J,k) = segment%ry_normal(segment%HI%ied,J,k) + ry_tangential(segment%HI%IsdB,J,k) = segment%rx_normal(segment%HI%isd,J,k) + ry_tangential(segment%HI%IedB,J,k) = segment%rx_normal(segment%HI%ied,J,k) + cff_tangential(segment%HI%IsdB,J,k) = segment%cff_normal(segment%HI%isd,J,k) + cff_tangential(segment%HI%IedB,J,k) = segment%cff_normal(segment%HI%ied,J,k) + do I=segment%HI%IsdB+1,segment%HI%IedB-1 + rx_tangential(I,J,k) = 0.5*(segment%ry_normal(i,J,k) + segment%ry_normal(i+1,J,k)) + ry_tangential(I,J,k) = 0.5*(segment%rx_normal(i,J,k) + segment%rx_normal(i+1,J,k)) + cff_tangential(I,J,k) = 0.5*(segment%cff_normal(i,J,k) + segment%cff_normal(i+1,J,k)) + enddo + enddo + if (segment%oblique_tan) then + do k=1,nz ; do I=segment%HI%IsdB+1,segment%HI%IedB-1 + rx_avg = rx_tangential(I,J,k) + ry_avg = ry_tangential(I,J,k) + cff_avg = cff_tangential(I,J,k) + segment%tangential_vel(I,J,k) = ((cff_avg*u_new(I,j+1,k) + rx_avg*u_new(I,j+2,k)) - & + (max(ry_avg,0.0)*segment%grad_tan(i,2,k) + min(ry_avg,0.0)*segment%grad_tan(i+1,2,k))) / & + (cff_avg + rx_avg) + enddo ; enddo + endif + if (segment%nudged_tan) then + do k=1,nz ; do I=segment%HI%IsdB,segment%HI%IedB + ! dhdt gets set to 0 on inflow in oblique case + if (ry_tangential(I,J,k) <= 0.0) then + tau = segment%Velocity_nudging_timescale_in + else + tau = segment%Velocity_nudging_timescale_out + endif + gamma_2 = dt / (tau + dt) + segment%tangential_vel(I,J,k) = (1 - gamma_2) * segment%tangential_vel(I,J,k) + & + gamma_2 * segment%nudged_tangential_vel(I,J,k) + enddo ; enddo + endif + if (segment%oblique_grad) then + Is_obc = max(segment%HI%IsdB,G%isd+1) + Ie_obc = min(segment%HI%IedB,G%ied-1) + do k=1,nz ; do I=segment%HI%IsdB+1,segment%HI%IedB-1 + rx_avg = rx_tangential(I,J,k) + ry_avg = ry_tangential(I,J,k) + cff_avg = cff_tangential(I,J,k) + segment%tangential_grad(I,J,k) = ((cff_avg*(u_new(I,j+2,k) - u_new(I,j+1,k))*G%IdyBu(I,J+1) & + + rx_avg*(u_new(I,j+3,k) - u_new(I,j+2,k))*G%IdyBu(I,J+2)) - & + (max(ry_avg,0.0)*segment%grad_gradient(i,2,k) + min(ry_avg,0.0)*segment%grad_gradient(i+1,2,k))) / & + (cff_avg + rx_avg) + enddo ; enddo + endif + if (segment%nudged_grad) then + do k=1,nz ; do J=segment%HI%JsdB,segment%HI%JedB + ! dhdt gets set to 0 on inflow in oblique case + if (ry_tangential(I,J,k) <= 0.0) then tau = segment%Velocity_nudging_timescale_in else tau = segment%Velocity_nudging_timescale_out @@ -2048,6 +2378,8 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, dt) enddo ; enddo endif deallocate(rx_tangential) + deallocate(ry_tangential) + deallocate(cff_tangential) endif endif enddo @@ -2144,6 +2476,24 @@ subroutine gradient_at_q_points(G, segment, uvel, vvel) segment%grad_normal(J,2,k) = (uvel(I,j+1,k)-uvel(I,j,k)) * G%mask2dBu(I,J) enddo enddo + if (segment%oblique_tan) then + do k=1,G%ke + do J=max(segment%HI%jsd, G%HI%jsd+1),min(segment%HI%jed, G%HI%jed-1) + segment%grad_tan(j,1,k) = (vvel(i-1,J,k)-vvel(i-1,J-1,k)) * G%mask2dT(i-1,j) + segment%grad_tan(j,2,k) = (vvel(i,J,k)-vvel(i,J-1,k)) * G%mask2dT(i,j) + enddo + enddo + endif + if (segment%oblique_grad) then + do k=1,G%ke + do J=max(segment%HI%jsd, G%HI%jsd+1),min(segment%HI%jed, G%HI%jed-1) + segment%grad_gradient(j,1,k) = (((vvel(i-1,J,k) - vvel(i-2,J,k))*G%IdxBu(I-2,J)) - & + (vvel(i-1,J-1,k) - vvel(i-2,J-1,k))*G%IdxBu(I-2,J-1)) * G%mask2dCu(I-1,j) + segment%grad_gradient(j,2,k) = (((vvel(i,J,k) - vvel(i-1,J,k))*G%IdxBu(I-1,J)) - & + (vvel(i,J-1,k) - vvel(i-1,J-1,k))*G%IdxBu(I-1,J-1)) * G%mask2dCu(I,j) + enddo + enddo + endif else ! western segment I=segment%HI%isdB do k=1,G%ke @@ -2152,6 +2502,24 @@ subroutine gradient_at_q_points(G, segment, uvel, vvel) segment%grad_normal(J,2,k) = (uvel(I,j+1,k)-uvel(I,j,k)) * G%mask2dBu(I,J) enddo enddo + if (segment%oblique_tan) then + do k=1,G%ke + do J=max(segment%HI%jsd, G%HI%jsd+1),min(segment%HI%jed, G%HI%jed-1) + segment%grad_tan(j,1,k) = (vvel(i+2,J,k)-vvel(i+2,J-1,k)) * G%mask2dT(i+2,j) + segment%grad_tan(j,2,k) = (vvel(i+1,J,k)-vvel(i+1,J-1,k)) * G%mask2dT(i+1,j) + enddo + enddo + endif + if (segment%oblique_grad) then + do k=1,G%ke + do J=max(segment%HI%jsd, G%HI%jsd+1),min(segment%HI%jed, G%HI%jed-1) + segment%grad_gradient(j,1,k) = (((vvel(i+3,J,k) - vvel(i+2,J,k))*G%IdxBu(I+2,J)) - & + (vvel(i+3,J-1,k) - vvel(i+2,J-1,k))*G%IdxBu(I+2,J-1)) * G%mask2dCu(I+2,j) + segment%grad_gradient(j,2,k) = (((vvel(i+2,J,k) - vvel(i+1,J,k))*G%IdxBu(I+1,J)) - & + (vvel(i+2,J-1,k) - vvel(i+1,J-1,k))*G%IdxBu(I+1,J-1)) * G%mask2dCu(I+1,j) + enddo + enddo + endif endif elseif (segment%is_N_or_S) then if (segment%direction == OBC_DIRECTION_N) then @@ -2162,6 +2530,24 @@ subroutine gradient_at_q_points(G, segment, uvel, vvel) segment%grad_normal(I,2,k) = (vvel(i+1,J,k)-vvel(i,J,k)) * G%mask2dBu(I,J) enddo enddo + if (segment%oblique_tan) then + do k=1,G%ke + do I=max(segment%HI%isd, G%HI%isd+1),min(segment%HI%ied, G%HI%ied-1) + segment%grad_tan(i,1,k) = (uvel(I,j-1,k)-uvel(I-1,j-1,k)) * G%mask2dT(i,j-1) + segment%grad_tan(i,2,k) = (uvel(I,j,k)-uvel(I-1,j,k)) * G%mask2dT(i,j) + enddo + enddo + endif + if (segment%oblique_grad) then + do k=1,G%ke + do I=max(segment%HI%isd, G%HI%isd+1),min(segment%HI%ied, G%HI%ied-1) + segment%grad_gradient(i,1,k) = (((uvel(I,j-1,k) - uvel(I,j-2,k))*G%IdxBu(I,J-2)) - & + (uvel(I-1,j-1,k) - uvel(I-1,j-2,k))*G%IdxBu(I-1,J-2)) * G%mask2dCv(I,j-1) + segment%grad_gradient(i,2,k) = (((uvel(I,j,k) - uvel(I,j-1,k))*G%IdyBu(I,J-1)) - & + (uvel(I-1,j,k) - uvel(I-1,j-1,k))*G%IdyBu(I-1,J-1)) * G%mask2dCv(i,J) + enddo + enddo + endif else ! south segment J=segment%HI%jsdB do k=1,G%ke @@ -2170,6 +2556,24 @@ subroutine gradient_at_q_points(G, segment, uvel, vvel) segment%grad_normal(I,2,k) = (vvel(i+1,J,k)-vvel(i,J,k)) * G%mask2dBu(I,J) enddo enddo + if (segment%oblique_tan) then + do k=1,G%ke + do I=max(segment%HI%isd, G%HI%isd+1),min(segment%HI%ied, G%HI%ied-1) + segment%grad_tan(i,1,k) = (uvel(I,j+2,k)-uvel(I-1,j+2,k)) * G%mask2dT(i,j+2) + segment%grad_tan(i,2,k) = (uvel(I,j+1,k)-uvel(I-1,j+1,k)) * G%mask2dT(i,j+1) + enddo + enddo + endif + if (segment%oblique_grad) then + do k=1,G%ke + do I=max(segment%HI%isd, G%HI%isd+1),min(segment%HI%ied, G%HI%ied-1) + segment%grad_gradient(i,1,k) = (((uvel(I,j+3,k) - uvel(I,j+2,k))*G%IdxBu(I,J+2)) - & + (uvel(I-1,j+3,k) - uvel(I-1,j+2,k))*G%IdyBu(I-1,J+2)) * G%mask2dCv(i,J+2) + segment%grad_gradient(i,2,k) = (((uvel(I,j+2,k) - uvel(I,j+1,k))*G%IdxBu(I,J+1)) - & + (uvel(I-1,j+2,k) - uvel(I-1,j+1,k))*G%IdyBu(I-1,J+1)) * G%mask2dCv(i,J+1) + enddo + enddo + endif endif endif @@ -2329,7 +2733,8 @@ subroutine allocate_OBC_segment_data(OBC, segment) if (segment%nudged_grad) then allocate(segment%nudged_tangential_grad(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%nudged_tangential_grad(:,:,:)=0.0 endif - if (OBC%specified_vorticity .or. OBC%specified_strain .or. segment%radiation_grad) then + if (OBC%specified_vorticity .or. OBC%specified_strain .or. segment%radiation_grad .or. & + segment%oblique_grad) then allocate(segment%tangential_grad(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%tangential_grad(:,:,:)=0.0 endif if (segment%oblique) then @@ -2338,6 +2743,12 @@ subroutine allocate_OBC_segment_data(OBC, segment) allocate(segment%ry_normal(IsdB:IedB,jsd:jed,OBC%ke)); segment%ry_normal(:,:,:)=0.0 allocate(segment%cff_normal(IsdB:IedB,jsd:jed,OBC%ke)); segment%cff_normal(:,:,:)=0.0 endif + if (segment%oblique_tan) then + allocate(segment%grad_tan(jsd:jed,2,OBC%ke)); segment%grad_tan(:,:,:) = 0.0 + endif + if (segment%oblique_grad) then + allocate(segment%grad_gradient(jsd:jed,2,OBC%ke)); segment%grad_gradient(:,:,:) = 0.0 + endif endif if (segment%is_N_or_S) then @@ -2365,7 +2776,8 @@ subroutine allocate_OBC_segment_data(OBC, segment) if (segment%nudged_grad) then allocate(segment%nudged_tangential_grad(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%nudged_tangential_grad(:,:,:)=0.0 endif - if (OBC%specified_vorticity .or. OBC%specified_strain .or. segment%radiation_grad) then + if (OBC%specified_vorticity .or. OBC%specified_strain .or. segment%radiation_grad .or. & + segment%oblique_grad) then allocate(segment%tangential_grad(IsdB:IedB,JsdB:JedB,OBC%ke)); segment%tangential_grad(:,:,:)=0.0 endif if (segment%oblique) then @@ -2374,6 +2786,12 @@ subroutine allocate_OBC_segment_data(OBC, segment) allocate(segment%ry_normal(isd:ied,JsdB:JedB,OBC%ke)); segment%ry_normal(:,:,:)=0.0 allocate(segment%cff_normal(isd:ied,JsdB:JedB,OBC%ke)); segment%cff_normal(:,:,:)=0.0 endif + if (segment%oblique_tan) then + allocate(segment%grad_tan(isd:ied,2,OBC%ke)); segment%grad_tan(:,:,:) = 0.0 + endif + if (segment%oblique_grad) then + allocate(segment%grad_gradient(isd:ied,2,OBC%ke)); segment%grad_gradient(:,:,:) = 0.0 + endif endif end subroutine allocate_OBC_segment_data @@ -2394,12 +2812,16 @@ subroutine deallocate_OBC_segment_data(OBC, segment) if (associated (segment%rx_normal)) deallocate(segment%rx_normal) if (associated (segment%ry_normal)) deallocate(segment%ry_normal) if (associated (segment%cff_normal)) deallocate(segment%cff_normal) + if (associated (segment%grad_normal)) deallocate(segment%grad_normal) + if (associated (segment%grad_tan)) deallocate(segment%grad_tan) + if (associated (segment%grad_gradient)) deallocate(segment%grad_gradient) if (associated (segment%normal_vel)) deallocate(segment%normal_vel) if (associated (segment%normal_vel_bt)) deallocate(segment%normal_vel_bt) if (associated (segment%normal_trans)) deallocate(segment%normal_trans) if (associated (segment%nudged_normal_vel)) deallocate(segment%nudged_normal_vel) if (associated (segment%tangential_vel)) deallocate(segment%tangential_vel) if (associated (segment%nudged_tangential_vel)) deallocate(segment%nudged_tangential_vel) + if (associated (segment%nudged_tangential_grad)) deallocate(segment%nudged_tangential_grad) if (associated (segment%tangential_grad)) deallocate(segment%tangential_grad) if (associated (segment%tr_Reg)) call segment_tracer_registry_end(segment%tr_Reg) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 4165fb0e11..4a2dbbea54 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -338,7 +338,7 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & if (present(gas_fields_ocn)) & call coupler_type_spawn(gas_fields_ocn, sfc_state%tr_fields, & - (/isd,is,ie,ied/), (/jsd,js,je,jed/), as_needed=.true.) + (/is,is,ie,ie/), (/js,js,je,je/), as_needed=.true.) sfc_state%arrays_allocated = .true. diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index f200a15bed..8e18ed5a01 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/src/diagnostics/MOM_diagnostics.F90 @@ -1952,6 +1952,14 @@ subroutine write_static_fields(G, GV, tv, diag) 'Open zonal grid spacing at v points (meter)', 'm', interp_method='none') if (id > 0) call post_data(id, G%dx_Cv, diag, .true.) + id = register_static_field('ocean_model', 'sin_rot', diag%axesT1, & + 'sine of the clockwise angle of the ocean grid north to true north', 'none') + if (id > 0) call post_data(id, G%sin_rot, diag, .true.) + + id = register_static_field('ocean_model', 'cos_rot', diag%axesT1, & + 'cosine of the clockwise angle of the ocean grid north to true north', 'none') + if (id > 0) call post_data(id, G%cos_rot, diag, .true.) + ! This static diagnostic is from CF 1.8, and is the fraction of a cell ! covered by ocean, given as a percentage (poorly named). diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index 103b328aa1..a38facf79a 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -118,7 +118,6 @@ module MOM_domains !! domain in the i-direction in a define_domain call. integer :: Y_FLAGS !< Flag that specifies the properties of the !! domain in the j-direction in a define_domain call. - logical :: use_io_layout !< True if an I/O layout is available. logical, pointer :: maskmap(:,:) => NULL() !< A pointer to an array indicating !! which logical processors are actually used for !! the ocean code. The other logical processors @@ -1401,11 +1400,11 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & "STATIC_MEMORY_ this is set in "//trim(inc_nm)//" at compile time.",& layoutParam=.true.) call log_param(param_file, mdl, trim(njproc_nm), layout(2), & - "The number of processors in the x-direction. With \n"//& !### FIX THIS COMMENT + "The number of processors in the y-direction. With \n"//& "STATIC_MEMORY_ this is set in "//trim(inc_nm)//" at compile time.",& layoutParam=.true.) call log_param(param_file, mdl, trim(layout_nm), layout, & - "The processor layout that was acutally used.",& + "The processor layout that was actually used.",& layoutParam=.true.) ! Idiot check that fewer PEs than columns have been requested @@ -1490,7 +1489,6 @@ subroutine MOM_domains_init(MOM_dom, param_file, symmetric, static_memory, & MOM_dom%Y_FLAGS = Y_FLAGS MOM_dom%layout = layout MOM_dom%io_layout = io_layout - MOM_dom%use_io_layout = (io_layout(1) + io_layout(2) > 0) if (is_static) then ! A requirement of equal sized compute domains is necessary when STATIC_MEMORY_ @@ -1554,7 +1552,6 @@ subroutine clone_MD_to_MD(MD_in, MOM_dom, min_halo, halo_size, symmetric, & MOM_dom%X_FLAGS = MD_in%X_FLAGS ; MOM_dom%Y_FLAGS = MD_in%Y_FLAGS MOM_dom%layout(:) = MD_in%layout(:) ; MOM_dom%io_layout(:) = MD_in%io_layout(:) - MOM_dom%use_io_layout = (MOM_dom%io_layout(1) + MOM_dom%io_layout(2) > 0) if (associated(MD_in%maskmap)) then mask_table_exists = .true. diff --git a/src/framework/MOM_file_parser.F90 b/src/framework/MOM_file_parser.F90 index ae876b16dd..72944c4f7a 100644 --- a/src/framework/MOM_file_parser.F90 +++ b/src/framework/MOM_file_parser.F90 @@ -6,8 +6,8 @@ module MOM_file_parser use MOM_coms, only : root_PE, broadcast use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg use MOM_error_handler, only : is_root_pe, stdlog, stdout -use MOM_time_manager, only : set_time, get_time, time_type, get_ticks_per_second -use MOM_time_manager, only : set_date, get_date +use MOM_time_manager, only : get_time, time_type, get_ticks_per_second +use MOM_time_manager, only : set_date, get_date, real_to_time, operator(-), set_time use MOM_document, only : doc_param, doc_module, doc_init, doc_end, doc_type use MOM_document, only : doc_openBlock, doc_closeBlock use MOM_string_functions, only : left_int, left_ints, slasher @@ -821,7 +821,7 @@ subroutine read_param_time(CS, varname, value, timeunit, fail_if_missing, date_f character(len=240) :: err_msg logical :: found, defined real :: real_time, time_unit - integer :: days, secs, vals(7) + integer :: vals(7) if (present(date_format)) date_format = .false. @@ -854,10 +854,8 @@ subroutine read_param_time(CS, varname, value, timeunit, fail_if_missing, date_f else time_unit = 1.0 ; if (present(timeunit)) time_unit = timeunit read( value_string(1), *) real_time - days = int(real_time*(time_unit/86400.0)) - secs = int(floor((real_time*(time_unit/86400.0)-days)*86400.0 + 0.5)) - value = set_time(secs, days) - endif + value = real_to_time(real_time*time_unit) + endif else if (present(fail_if_missing)) then ; if (fail_if_missing) then if (.not.found) then diff --git a/src/framework/MOM_horizontal_regridding.F90 b/src/framework/MOM_horizontal_regridding.F90 index 19b73ee07f..c7befad3b3 100644 --- a/src/framework/MOM_horizontal_regridding.F90 +++ b/src/framework/MOM_horizontal_regridding.F90 @@ -19,7 +19,7 @@ module MOM_horizontal_regridding use MOM_io, only : open_file, read_data, read_axis_data, SINGLE_FILE, MULTIPLE use MOM_io, only : slasher, vardesc, write_field use MOM_string_functions, only : uppercase -use MOM_time_manager, only : time_type, set_time, get_external_field_size +use MOM_time_manager, only : time_type, get_external_field_size use MOM_time_manager, only : init_external_field, time_interp_external use MOM_time_manager, only : get_external_field_axes, get_external_field_missing use MOM_variables, only : thermo_var_ptrs diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index 21d42ea436..e523270802 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -30,7 +30,7 @@ module MOM_io use mpp_io_mod, only : MPP_APPEND, MPP_MULTI, MPP_OVERWR, MPP_NETCDF, MPP_RDONLY use mpp_io_mod, only : get_file_info=>mpp_get_info, get_file_atts=>mpp_get_atts use mpp_io_mod, only : get_file_fields=>mpp_get_fields, get_file_times=>mpp_get_times -use mpp_io_mod, only : read_field=>mpp_read, io_infra_init=>mpp_io_init +use mpp_io_mod, only : io_infra_init=>mpp_io_init use netcdf @@ -38,7 +38,7 @@ module MOM_io public :: close_file, create_file, field_exists, field_size, fieldtype, get_filename_appendix public :: file_exists, flush_file, get_file_info, get_file_atts, get_file_fields -public :: get_file_times, open_file, read_axis_data, read_data, read_field +public :: get_file_times, open_file, read_axis_data, read_data public :: num_timelevels, MOM_read_data, MOM_read_vector, ensembler public :: reopen_file, slasher, write_field, write_version_number, MOM_io_init public :: open_namelist_file, check_nml_error, io_infra_init, io_infra_end @@ -154,9 +154,7 @@ subroutine create_file(unit, filename, vars, novars, fields, threading, timeunit endif one_file = .true. - if (domain_set) then - one_file = ((thread == SINGLE_FILE) .or. .not.Domain%use_io_layout) - endif + if (domain_set) one_file = (thread == SINGLE_FILE) if (one_file) then call open_file(unit, filename, MPP_OVERWR, MPP_NETCDF, threading=thread) @@ -398,9 +396,7 @@ subroutine reopen_file(unit, filename, vars, novars, fields, threading, timeunit endif one_file = .true. - if (domain_set) then - one_file = ((thread == SINGLE_FILE) .or. .not.Domain%use_io_layout) - endif + if (domain_set) one_file = (thread == SINGLE_FILE) if (one_file) then call open_file(unit, filename, MPP_APPEND, MPP_NETCDF, threading=thread) @@ -1012,7 +1008,7 @@ end subroutine MOM_io_init !! !! * write_field: write a field to an open file. !! * write_time: write a value of the time axis to an open file. -!! * read_field: read a field from an open file. +!! * read_data: read a variable from an open file. !! * read_time: read a time from an open file. !! !! * name_output_file: provide a name for an output file based on a diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index bf40da4897..81407d783c 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -9,15 +9,15 @@ module MOM_restart use MOM_string_functions, only : lowercase use MOM_grid, only : ocean_grid_type use MOM_io, only : create_file, fieldtype, file_exists, open_file, close_file -use MOM_io, only : read_field, write_field, MOM_read_data, read_data, get_filename_appendix +use MOM_io, only : write_field, MOM_read_data, read_data, get_filename_appendix use MOM_io, only : get_file_info, get_file_atts, get_file_fields, get_file_times use MOM_io, only : vardesc, var_desc, query_vardesc, modify_vardesc use MOM_io, only : MULTIPLE, NETCDF_FILE, READONLY_FILE, SINGLE_FILE use MOM_io, only : CENTER, CORNER, NORTH_FACE, EAST_FACE -use MOM_time_manager, only : time_type, get_time, get_date, set_date, set_time -use MOM_time_manager, only : days_in_month +use MOM_time_manager, only : time_type, time_type_to_real, real_to_time +use MOM_time_manager, only : days_in_month, get_date, set_date use MOM_verticalGrid, only : verticalGrid_type -use mpp_mod, only: mpp_chksum +use mpp_mod, only: mpp_chksum,mpp_pe use mpp_io_mod, only: mpp_attribute_exist, mpp_get_atts implicit none ; private @@ -801,15 +801,14 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV) ! With parallel read & write, it is possible to disable the following... -! jgj: this was set to 4294967292, changed to 4294967295 (see mpp_parameter.F90) - if (CS%large_file_support) max_file_size = 4294967295_8 + ! The maximum file size is 4294967292, according to the NetCDF documentation. + if (CS%large_file_support) max_file_size = 4294967292_8 num_files = 0 next_var = 0 nz = 1 ; if (present(GV)) nz = GV%ke - call get_time(time,seconds,days) - restart_time = real(days) + real(seconds)/86400.0 + restart_time = time_type_to_real(time) / 86400.0 restartname = trim(CS%restartfile) if (present(filename)) restartname = trim(filename) @@ -918,7 +917,7 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV) elseif (associated(CS%var_ptr1d(m)%p)) then check_val(m-start_var+1,1) = mpp_chksum(CS%var_ptr1d(m)%p) elseif (associated(CS%var_ptr0d(m)%p)) then - check_val(m-start_var+1,1) = mpp_chksum(CS%var_ptr0d(m)%p) + check_val(m-start_var+1,1) = mpp_chksum(CS%var_ptr0d(m)%p,pelist=(/mpp_pe()/)) endif enddo @@ -982,7 +981,7 @@ subroutine restore_state(filename, directory, day, G, CS) character(len=80) :: varname ! A variable's name. integer :: num_file ! The number of files (restart files and others ! explicitly in filename) that are open. - integer :: i, n, m, start_of_day, num_days, missing_fields + integer :: i, n, m, missing_fields integer :: isL, ieL, jsL, jeL, is0, js0 integer :: sizes(7) integer :: ndim, nvar, natt, ntime, pos @@ -1028,9 +1027,7 @@ subroutine restore_state(filename, directory, day, G, CS) t1 = time_vals(1) deallocate(time_vals) - start_of_day = INT((t1 - INT(t1)) *86400) ! Number of seconds. - num_days = INT(t1) - day = set_time(start_of_day, num_days) + day = real_to_time(t1*86400.0) exit enddo @@ -1096,117 +1093,46 @@ subroutine restore_state(filename, directory, day, G, CS) call mpp_get_atts(fields(i),checksum=checksum_file) is_there_a_checksum = .true. endif - if (.NOT. CS%checksum_required ) is_there_a_checksum = .false. ! Do not need to do data checksumming. + if (.NOT. CS%checksum_required) is_there_a_checksum = .false. ! Do not need to do data checksumming. if (associated(CS%var_ptr1d(m)%p)) then ! Read a 1d array, which should be invariant to domain decomposition. call read_data(unit_path(n), varname, CS%var_ptr1d(m)%p, & - no_domain=.true., timelevel=1) - if ( is_there_a_checksum ) checksum_data = mpp_chksum(CS%var_ptr1d(m)%p) + G%Domain%mpp_domain, timelevel=1) + if (is_there_a_checksum) checksum_data = mpp_chksum(CS%var_ptr1d(m)%p) elseif (associated(CS%var_ptr0d(m)%p)) then ! Read a scalar... call read_data(unit_path(n), varname, CS%var_ptr0d(m)%p, & - no_domain=.true., timelevel=1) - if ( is_there_a_checksum ) checksum_data = mpp_chksum(CS%var_ptr0d(m)%p) - elseif ((pos == 0) .and. associated(CS%var_ptr2d(m)%p)) then ! Read a non-decomposed 2d array. - ! Probably should query the field type to make sure that the sizes are right. - call read_data(unit_path(n), varname, CS%var_ptr2d(m)%p, & - no_domain=.true., timelevel=1) - if ( is_there_a_checksum ) checksum_data = mpp_chksum(CS%var_ptr2d(m)%p(isL:ieL,jsL:jeL)) - elseif ((pos == 0) .and. associated(CS%var_ptr3d(m)%p)) then ! Read a non-decomposed 3d array. - ! Probably should query the field type to make sure that the sizes are right. - call read_data(unit_path(n), varname, CS%var_ptr3d(m)%p, & - no_domain=.true., timelevel=1) - if ( is_there_a_checksum ) checksum_data = mpp_chksum(CS%var_ptr3d(m)%p(isL:ieL,jsL:jeL,:)) - elseif ((pos == 0) .and. associated(CS%var_ptr4d(m)%p)) then ! Read a non-decomposed 4d array. - ! Probably should query the field type to make sure that the sizes are right. - call read_data(unit_path(n), varname, CS%var_ptr4d(m)%p, & - no_domain=.true., timelevel=1) - if ( is_there_a_checksum ) checksum_data = mpp_chksum(CS%var_ptr4d(m)%p(isL:ieL,jsL:jeL,:,:)) - elseif (unit_is_global(n) .or. G%Domain%use_io_layout) then - if (associated(CS%var_ptr3d(m)%p)) then - ! Read 3d array... Time level 1 is always used. - call MOM_read_data(unit_path(n), varname, CS%var_ptr3d(m)%p, & - G%Domain, 1, position=pos) - if ( is_there_a_checksum ) checksum_data = mpp_chksum(CS%var_ptr3d(m)%p(isL:ieL,jsL:jeL,:)) - elseif (associated(CS%var_ptr2d(m)%p)) then ! Read 2d array... + G%Domain%mpp_domain, timelevel=1) + if (is_there_a_checksum) checksum_data = mpp_chksum(CS%var_ptr0d(m)%p,pelist=(/mpp_pe()/)) + elseif (associated(CS%var_ptr2d(m)%p)) then ! Read a 2d array. + if (pos /= 0) then call MOM_read_data(unit_path(n), varname, CS%var_ptr2d(m)%p, & - G%Domain, 1, position=pos) - if ( is_there_a_checksum ) checksum_data = mpp_chksum(CS%var_ptr2d(m)%p(isL:ieL,jsL:jeL)) - elseif (associated(CS%var_ptr4d(m)%p)) then ! Read 4d array... - call MOM_read_data(unit_path(n), varname, CS%var_ptr4d(m)%p, & - G%Domain, 1, position=pos) - if ( is_there_a_checksum ) checksum_data = mpp_chksum(CS%var_ptr4d(m)%p(isL:ieL,jsL:jeL,:,:)) - else - call MOM_error(FATAL, "MOM_restart restore_state: "//& - "No pointers set for "//trim(varname)) + G%Domain, timelevel=1, position=pos) + else ! This array is not domain-decomposed. This variant may be under-tested. + call read_data(unit_path(n), varname, CS%var_ptr2d(m)%p, & + no_domain=.true., timelevel=1) endif - else ! Do not use an io_layout. !### GET RID OF THIS BRANCH ONCE read_data_4d_new IS AVAILABLE. - ! This file is decomposed onto the current processors. We need - ! to check whether the sizes look right, and abort if not. - call get_file_atts(fields(i),ndim=ndim,siz=sizes) - - ! NOTE: The index ranges f var_ptrs always start with 1, so with - ! symmetric memory the staggering is swapped from NE to SW! - is0 = 1-G%isd - if ((pos == EAST_FACE) .or. (pos == CORNER)) is0 = 1-G%IsdB - if (sizes(1) == G%iec-G%isc+1) then - isL = G%isc+is0 ; ieL = G%iec+is0 - elseif (sizes(1) == G%IecB-G%IscB+1) then - isL = G%IscB+is0 ; ieL = G%IecB+is0 - elseif (((pos == EAST_FACE) .or. (pos == CORNER)) .and. & - (G%IscB == G%isc) .and. (sizes(1) == G%iec-G%isc+2)) then - ! This is reading a symmetric file in a non-symmetric model. - isL = G%isc-1+is0 ; ieL = G%iec+is0 - else - call MOM_error(WARNING, "MOM_restart restore_state, "//trim(varname)//& - " has the wrong i-size in "//trim(unit_path(n))) - exit - endif - - js0 = 1-G%jsd - if ((pos == NORTH_FACE) .or. (pos == CORNER)) js0 = 1-G%JsdB - if (sizes(2) == G%jec-G%jsc+1) then - jsL = G%jsc+js0 ; jeL = G%jec+js0 - elseif (sizes(2) == G%jecB-G%jscB+1) then - jsL = G%jscB+js0 ; jeL = G%jecB+js0 - elseif (((pos == NORTH_FACE) .or. (pos == CORNER)) .and. & - (G%JscB == G%jsc) .and. (sizes(2) == G%jec-G%jsc+2)) then - ! This is reading a symmetric file in a non-symmetric model. - jsL = G%jsc-1+js0 ; jeL = G%jec+js0 - else - call MOM_error(WARNING, "MOM_restart restore_state, "//trim(varname)//& - " has the wrong j-size in "//trim(unit_path(n))) - exit + if (is_there_a_checksum) checksum_data = mpp_chksum(CS%var_ptr2d(m)%p(isL:ieL,jsL:jeL)) + elseif (associated(CS%var_ptr3d(m)%p)) then ! Read a 3d array. + if (pos /= 0) then + call MOM_read_data(unit_path(n), varname, CS%var_ptr3d(m)%p, & + G%Domain, timelevel=1, position=pos) + else ! This array is not domain-decomposed. This variant may be under-tested. + call read_data(unit_path(n), varname, CS%var_ptr3d(m)%p, & + no_domain=.true., timelevel=1) endif - - if (associated(CS%var_ptr3d(m)%p)) then - if (ntime == 0) then - call read_field(unit(n), fields(i), & - CS%var_ptr3d(m)%p(isL:ieL,jsL:jeL,:)) - else - call read_field(unit(n), fields(i), & - CS%var_ptr3d(m)%p(isL:ieL,jsL:jeL,:), 1) - endif - elseif (associated(CS%var_ptr2d(m)%p)) then - if (ntime == 0) then - call read_field(unit(n), fields(i), & - CS%var_ptr2d(m)%p(isL:ieL,jsL:jeL)) - else - call read_field(unit(n), fields(i), & - CS%var_ptr2d(m)%p(isL:ieL,jsL:jeL), 1) - endif - elseif (associated(CS%var_ptr4d(m)%p)) then - if (ntime == 0) then - call read_field(unit(n), fields(i), & - CS%var_ptr4d(m)%p(isL:ieL,jsL:jeL,:,:)) - else - call read_field(unit(n), fields(i), & - CS%var_ptr4d(m)%p(isL:ieL,jsL:jeL,:,:), 1) - endif - else - call MOM_error(FATAL, "MOM_restart restore_state: "//& - "No pointers set for "//trim(varname)) + if (is_there_a_checksum) checksum_data = mpp_chksum(CS%var_ptr3d(m)%p(isL:ieL,jsL:jeL,:)) + elseif (associated(CS%var_ptr4d(m)%p)) then ! Read a 4d array. + if (pos /= 0) then + call MOM_read_data(unit_path(n), varname, CS%var_ptr4d(m)%p, & + G%Domain, timelevel=1, position=pos) + else ! This array is not domain-decomposed. This variant may be under-tested. + call read_data(unit_path(n), varname, CS%var_ptr4d(m)%p, & + no_domain=.true., timelevel=1) endif + if (is_there_a_checksum) checksum_data = mpp_chksum(CS%var_ptr4d(m)%p(isL:ieL,jsL:jeL,:,:)) + else + call MOM_error(FATAL, "MOM_restart restore_state: No pointers set for "//trim(varname)) endif if (is_root_pe() .and. is_there_a_checksum .and. (checksum_file(1) /= checksum_data)) then @@ -1371,8 +1297,12 @@ function open_restart_units(filename, directory, G, CS, units, file_paths, & enddo fname = filename(start_char:m-1) start_char = m - do while ((start_char <= len_trim(filename)) .and. (filename(start_char:start_char) == ' ')) - start_char = start_char + 1 + do while (start_char <= len_trim(filename)) + if (filename(start_char:start_char) == ' ') then + start_char = start_char + 1 + else + exit + endif enddo if ((fname(1:1)=='r') .and. ( len_trim(fname) == 1)) then @@ -1411,24 +1341,11 @@ function open_restart_units(filename, directory, G, CS, units, file_paths, & threading = MULTIPLE, fileset = SINGLE_FILE) if (present(global_files)) global_files(n) = .true. elseif (CS%parallel_restartfiles) then - if (G%Domain%use_io_layout) then - ! Look for decomposed files using the I/O Layout. - fexists = file_exists(filepath, G%Domain) - if (fexists .and. (present(units))) & - call open_file(units(n), trim(filepath), READONLY_FILE, NETCDF_FILE, & - domain=G%Domain%mpp_domain) - else - ! Look for any PE-specific files of the form NAME.nc.####. - if (num_PEs()>10000) then - write(filepath, '(a,i6.6)' ) trim(filepath)//'.', pe_here() - else - write(filepath, '(a,i4.4)' ) trim(filepath)//'.', pe_here() - endif - inquire(file=filepath, exist=fexists) - if (fexists .and. (present(units))) & - call open_file(units(n), trim(filepath), READONLY_FILE, NETCDF_FILE, & - threading = MULTIPLE, fileset = SINGLE_FILE) - endif + ! Look for decomposed files using the I/O Layout. + fexists = file_exists(filepath, G%Domain) + if (fexists .and. (present(units))) & + call open_file(units(n), trim(filepath), READONLY_FILE, NETCDF_FILE, & + domain=G%Domain%mpp_domain) if (fexists .and. present(global_files)) global_files(n) = .false. endif diff --git a/src/framework/MOM_time_manager.F90 b/src/framework/MOM_time_manager.F90 index 25c367c1ef..229c3ded3a 100644 --- a/src/framework/MOM_time_manager.F90 +++ b/src/framework/MOM_time_manager.F90 @@ -20,8 +20,9 @@ module MOM_time_manager implicit none ; private -public :: time_type, get_time, set_time, time_type_to_real, real_to_time_type -public :: set_ticks_per_second , get_ticks_per_second +public :: time_type, get_time, set_time +public :: time_type_to_real, real_to_time_type, real_to_time +public :: set_ticks_per_second, get_ticks_per_second public :: operator(+), operator(-), operator(*), operator(/) public :: operator(>), operator(<), operator(>=), operator(<=) public :: operator(==), operator(/=), operator(//) @@ -35,4 +36,29 @@ module MOM_time_manager public :: get_external_field_axes public :: get_external_field_missing +contains + +!> This is an alternate implementation of the FMS function real_to_time_type that is accurate over +!! a larger range of input values. With 32 bit signed integers, this version should work over the +!! entire valid range (2^31 days or ~5.8835 million years) of time_types, whereas the standard +!! version in the FMS time_manager stops working for conversions of times greater than 2^31 seconds, +!! or ~68.1 years. +function real_to_time(x, err_msg) + type(time_type) :: real_to_time !< The output time as a time_type + real, intent(in) :: x !< The input time in real seconds. + character(len=*), intent(out), optional :: err_msg !< An optional returned error message. + + ! Local variables + integer :: seconds, days, ticks + real :: real_subsecond_remainder + + days = floor(x/86400.) + seconds = floor(x - 86400.*days) + real_subsecond_remainder = x - (days*86400. + seconds) + ticks = nint(real_subsecond_remainder * get_ticks_per_second()) + + real_to_time = set_time(seconds=seconds, days=days, ticks=ticks, err_msg=err_msg) +end function real_to_time + + end module MOM_time_manager diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index e1a61f355c..e6989caa54 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -25,7 +25,7 @@ module MOM_ice_shelf use MOM_io, only : write_field, close_file, SINGLE_FILE, MULTIPLE use MOM_restart, only : register_restart_field, query_initialized, save_restart use MOM_restart, only : restart_init, restore_state, MOM_restart_CS -use MOM_time_manager, only : time_type, set_time, time_type_to_real +use MOM_time_manager, only : time_type, time_type_to_real, time_type_to_real, real_to_time use MOM_transcribe_grid, only : copy_dyngrid_to_MOM_grid, copy_MOM_grid_to_dyngrid use MOM_variables, only : surface use MOM_forcing_type, only : forcing, allocate_forcing_type, MOM_forcing_chksum @@ -47,7 +47,7 @@ module MOM_ice_shelf use MOM_checksums, only : hchksum, qchksum, chksum, uchksum, vchksum, uvchksum use time_interp_external_mod, only : init_external_field, time_interp_external use time_interp_external_mod, only : time_interp_external_init -use time_manager_mod, only : print_time, time_type_to_real, real_to_time_type +use time_manager_mod, only : print_time implicit none ; private #include @@ -979,7 +979,7 @@ subroutine add_shelf_flux(G, CS, state, fluxes) ! just compute changes in mass after first time step if (t0>0.0) then - Time0 = real_to_time_type(t0) + Time0 = real_to_time(t0) last_hmask(:,:) = ISS%hmask(:,:) ; last_area_shelf_h(:,:) = ISS%area_shelf_h(:,:) call time_interp_external(CS%id_read_mass, Time0, last_mass_shelf) last_h_shelf = last_mass_shelf/CS%density_ice diff --git a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 index e504db90c7..9d25d8c8a3 100644 --- a/src/ice_shelf/MOM_ice_shelf_dynamics.F90 +++ b/src/ice_shelf/MOM_ice_shelf_dynamics.F90 @@ -17,7 +17,7 @@ module MOM_ice_shelf_dynamics use MOM_io, only : file_exists, slasher, MOM_read_data use MOM_restart, only : register_restart_field, query_initialized use MOM_restart, only : MOM_restart_CS -use MOM_time_manager, only : time_type, set_time, time_type_to_real +use MOM_time_manager, only : time_type, set_time !MJH use MOM_ice_shelf_initialize, only : initialize_ice_shelf_boundary use MOM_ice_shelf_state, only : ice_shelf_state use MOM_coms, only : reproducing_sum, sum_across_PEs, max_across_PEs, min_across_PEs @@ -523,13 +523,13 @@ subroutine initialize_diagnostic_fields(CS, ISS, G, Time) type(ocean_grid_type), intent(inout) :: G !< The grid structure used by the ice shelf. type(time_type), intent(in) :: Time !< The current model time - integer :: i, j, iters, isd, ied, jsd, jed - real :: rhoi, rhow, OD - type(time_type) :: dummy_time + integer :: i, j, iters, isd, ied, jsd, jed + real :: rhoi, rhow, OD + type(time_type) :: dummy_time rhoi = CS%density_ice rhow = CS%density_ocean_avg - dummy_time = set_time (0,0) + dummy_time = set_time(0,0) isd=G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed do j=jsd,jed diff --git a/src/initialization/MOM_grid_initialize.F90 b/src/initialization/MOM_grid_initialize.F90 index be82ffc33f..9f7c5dcc28 100644 --- a/src/initialization/MOM_grid_initialize.F90 +++ b/src/initialization/MOM_grid_initialize.F90 @@ -220,7 +220,6 @@ subroutine set_grid_metrics_from_mosaic(G, param_file) SGdom%niglobal = 2*G%domain%niglobal SGdom%njglobal = 2*G%domain%njglobal SGdom%layout(:) = G%domain%layout(:) - SGdom%use_io_layout = G%domain%use_io_layout SGdom%io_layout(:) = G%domain%io_layout(:) global_indices(1) = 1+SGdom%nihalo global_indices(2) = SGdom%niglobal+SGdom%nihalo @@ -241,8 +240,7 @@ subroutine set_grid_metrics_from_mosaic(G, param_file) symmetry=.true., name="MOM_MOSAIC") endif - if (SGdom%use_io_layout) & - call MOM_define_IO_domain(SGdom%mpp_domain, SGdom%io_layout) + call MOM_define_IO_domain(SGdom%mpp_domain, SGdom%io_layout) deallocate(exni) deallocate(exnj) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 0f9c17022d..57820accc0 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -34,7 +34,7 @@ module MOM_state_initialization use MOM_ALE_sponge, only : set_up_ALE_sponge_field, initialize_ALE_sponge use MOM_ALE_sponge, only : ALE_sponge_CS use MOM_string_functions, only : uppercase, lowercase -use MOM_time_manager, only : time_type, set_time +use MOM_time_manager, only : time_type use MOM_tracer_registry, only : tracer_registry_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : setVerticalGridAxes, verticalGrid_type diff --git a/src/initialization/MOM_tracer_initialization_from_Z.F90 b/src/initialization/MOM_tracer_initialization_from_Z.F90 index 67cf7bbd24..07be1ee340 100644 --- a/src/initialization/MOM_tracer_initialization_from_Z.F90 +++ b/src/initialization/MOM_tracer_initialization_from_Z.F90 @@ -16,7 +16,6 @@ module MOM_tracer_initialization_from_Z use MOM_get_input, only : directories use MOM_grid, only : ocean_grid_type, isPointInCell use MOM_string_functions, only : uppercase -use MOM_time_manager, only : time_type, set_time use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : setVerticalGridAxes use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type diff --git a/src/ocean_data_assim/MOM_oda_driver.F90 b/src/ocean_data_assim/MOM_oda_driver.F90 index 17cc300bd2..f9dae9b246 100644 --- a/src/ocean_data_assim/MOM_oda_driver.F90 +++ b/src/ocean_data_assim/MOM_oda_driver.F90 @@ -19,7 +19,7 @@ module MOM_oda_driver_mod use ensemble_manager_mod, only : get_ensemble_id, get_ensemble_size use ensemble_manager_mod, only : get_ensemble_pelist, get_ensemble_filter_pelist use time_manager_mod, only : time_type, decrement_time, increment_time -use time_manager_mod, only : get_date, get_time, operator(>=),operator(/=),operator(==),operator(<) +use time_manager_mod, only : get_date, operator(>=),operator(/=),operator(==),operator(<) use constants_mod, only : radius, epsln ! ODA Modules use ocean_da_types_mod, only : grid_type, ocean_profile_type, ocean_control_struct diff --git a/src/parameterizations/lateral/MOM_internal_tides.F90 b/src/parameterizations/lateral/MOM_internal_tides.F90 index 3205f81b02..822c11470e 100644 --- a/src/parameterizations/lateral/MOM_internal_tides.F90 +++ b/src/parameterizations/lateral/MOM_internal_tides.F90 @@ -18,9 +18,7 @@ module MOM_internal_tides use MOM_io, only : slasher, vardesc, MOM_read_data use MOM_restart, only : register_restart_field, MOM_restart_CS, restart_init, save_restart use MOM_spatial_means, only : global_area_mean -use MOM_time_manager, only : time_type, operator(+), operator(/), operator(-) -use MOM_time_manager, only : get_time, get_date, set_time, set_date -use MOM_time_manager, only : time_type_to_real +use MOM_time_manager, only : time_type, time_type_to_real, operator(+), operator(/), operator(-) use MOM_variables, only : surface, thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type use MOM_wave_structure, only: wave_structure_init, wave_structure, wave_structure_CS @@ -592,13 +590,8 @@ subroutine sum_En(G, CS, En, label) integer :: m,fr,a real :: En_sum, tmpForSumming, En_sum_diff, En_sum_pdiff character(len=160) :: mesg ! The text of an error message - integer :: seconds - real :: Isecs_per_day = 1.0 / 86400.0 real :: days - call get_time(CS%Time, seconds) - days = real(seconds) * Isecs_per_day - En_sum = 0.0 tmpForSumming = 0.0 do a=1,CS%nAngle @@ -614,6 +607,7 @@ subroutine sum_En(G, CS, En, label) CS%En_sum = En_sum !! Print to screen !if (is_root_pe()) then + ! days = time_type_to_real(CS%Time) / 86400.0 ! write(mesg,*) trim(label)//': days =', days, ', En_sum=', En_sum, & ! ', En_sum_diff=', En_sum_diff, ', Percent change=', En_sum_pdiff, '%' ! call MOM_mesg(mesg) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 188ba9c8f3..e3806fd684 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -61,8 +61,7 @@ module MOM_diabatic_driver use MOM_shortwave_abs, only : absorbRemainingSW, optics_type use MOM_sponge, only : apply_sponge, sponge_CS use MOM_ALE_sponge, only : apply_ALE_sponge, ALE_sponge_CS -use MOM_time_manager, only : operator(-), set_time -use MOM_time_manager, only : operator(<=), time_type ! for testing itides (BDM) +use MOM_time_manager, only : time_type, real_to_time, operator(-), operator(<=) use MOM_tracer_flow_control, only : call_tracer_column_fns, tracer_flow_control_CS use MOM_tracer_diabatic, only : tracer_vertdiff use MOM_variables, only : thermo_var_ptrs, vertvisc_type, accel_diag_ptrs @@ -440,7 +439,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! the end of the diabatic processes. if (associated(tv%T) .AND. associated(tv%frazil)) then ! For frazil diagnostic, the first call covers the first half of the time step - call enable_averaging(0.5*dt, Time_end - set_time(int(floor(0.5*dt+0.5))), CS%diag) + call enable_averaging(0.5*dt, Time_end - real_to_time(0.5*dt), CS%diag) if (CS%frazil_tendency_diag) then do k=1,nz ; do j=js,je ; do i=is,ie temp_diag(i,j,k) = tv%T(i,j,k) @@ -1316,7 +1315,7 @@ subroutine legacy_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_en ! the end of the diabatic processes. if (associated(tv%T) .AND. associated(tv%frazil)) then ! For frazil diagnostic, the first call covers the first half of the time step - call enable_averaging(0.5*dt, Time_end - set_time(int(floor(0.5*dt+0.5))), CS%diag) + call enable_averaging(0.5*dt, Time_end - real_to_time(0.5*dt), CS%diag) if (CS%frazil_tendency_diag) then do k=1,nz ; do j=js,je ; do i=is,ie temp_diag(i,j,k) = tv%T(i,j,k) diff --git a/src/parameterizations/vertical/MOM_opacity.F90 b/src/parameterizations/vertical/MOM_opacity.F90 index 26a23a0f0d..ca2afdc655 100644 --- a/src/parameterizations/vertical/MOM_opacity.F90 +++ b/src/parameterizations/vertical/MOM_opacity.F90 @@ -5,7 +5,6 @@ module MOM_opacity use MOM_diag_mediator, only : time_type, diag_ctrl, safe_alloc_ptr, post_data use MOM_diag_mediator, only : query_averaging_enabled, register_diag_field -use MOM_time_manager, only : get_time use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_string_functions, only : uppercase @@ -225,7 +224,6 @@ subroutine opacity_from_chl(optics, fluxes, G, CS, chl_in) ! radiation, in W m-2. type(time_type) :: day character(len=128) :: mesg - integer :: days, seconds integer :: i, j, k, n, is, ie, js, je, nz, nbands logical :: multiband_vis_input, multiband_nir_input @@ -271,7 +269,6 @@ subroutine opacity_from_chl(optics, fluxes, G, CS, chl_in) else ! Only the 2-d surface chlorophyll can be read in from a file. The ! same value is assumed for all layers. - call get_time(CS%Time,seconds,days) call time_interp_external(CS%sbc_chl, CS%Time, chl_data) do j=js,je ; do i=is,ie if ((G%mask2dT(i,j) > 0.5) .and. (chl_data(i,j) < 0.0)) then diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 88da20bb4d..6b5fcb3202 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -217,14 +217,12 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, CS, & Idt = 1.0 / dt !Check if Stokes mixing allowed if requested (present and associated) + DoStokesMixing=.false. if (CS%StokesMixing) then - DoStokesMixing=(present(Waves) .and. associated(Waves)) - if (.not.DoStokesMixing) then + if (present(Waves)) DoStokesMixing = associated(Waves) + if (.not. DoStokesMixing) & call MOM_error(FATAL,"Stokes Mixing called without allocated"//& - "Waves Control Structure") - endif - else - DoStokesMixing=.false. + "Waves Control Structure") endif do k=1,nz ; do i=Isq,Ieq ; Ray(i,k) = 0.0 ; enddo ; enddo @@ -232,17 +230,17 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, CS, & ! Update the zonal velocity component using a modification of a standard ! tridagonal solver. - ! When mixing down Eulerian current + Stokes drift add before calling solver - if (DoStokesMixing) then ; do k=1,nz ; do j=G%jsc,G%jec ; do I=Isq,Ieq - if (G%mask2dCu(I,j) > 0) u(I,j,k) = u(I,j,k) + Waves%Us_x(I,j,k) - enddo ; enddo ; enddo ; endif - !$OMP parallel do default(shared) firstprivate(Ray) & !$OMP private(do_i,surface_stress,zDS,stress,h_a,hfr, & !$OMP b_denom_1,b1,d1,c1) do j=G%jsc,G%jec do I=Isq,Ieq ; do_i(I) = (G%mask2dCu(I,j) > 0) ; enddo + ! When mixing down Eulerian current + Stokes drift add before calling solver + if (DoStokesMixing) then ; do k=1,nz ; do I=Isq,Ieq + if (do_i(I)) u(I,j,k) = u(I,j,k) + Waves%Us_x(I,j,k) + enddo ; enddo ; endif + if (associated(ADp%du_dt_visc)) then ; do k=1,nz ; do I=Isq,Ieq ADp%du_dt_visc(I,j,k) = u(I,j,k) enddo ; enddo ; endif @@ -330,19 +328,15 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, CS, & taux_bot(I,j) = taux_bot(I,j) + Rho0 * (Ray(I,k)*u(I,j,k)) enddo ; enddo ; endif endif - enddo ! end u-component j loop - ! When mixing down Eulerian current + Stokes drift subtract after calling solver - if (DoStokesMixing) then ; do k=1,nz ; do j=G%jsc,G%jec ; do I=Isq,Ieq - if (G%mask2dCu(I,j) > 0) u(I,j,k) = u(I,j,k) - Waves%Us_x(I,j,k) - enddo ; enddo ; enddo ; endif + ! When mixing down Eulerian current + Stokes drift subtract after calling solver + if (DoStokesMixing) then ; do k=1,nz ; do I=Isq,Ieq + if (do_i(I)) u(I,j,k) = u(I,j,k) - Waves%Us_x(I,j,k) + enddo ; enddo ; endif + + enddo ! end u-component j loop ! Now work on the meridional velocity component. - ! When mixing down Eulerian current + Stokes drift add before calling solver - if (DoStokesMixing) then ; do k=1,nz ; do j=Jsq,Jeq ; do I=Is,Ie - if (G%mask2dCv(I,j) > 0) & - v(i,j,k) = v(i,j,k) + Waves%Us_y(i,j,k) - enddo ; enddo ; enddo ; endif !$OMP parallel do default(shared) firstprivate(Ray) & !$OMP private(do_i,surface_stress,zDS,stress,h_a,hfr, & @@ -350,6 +344,11 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, CS, & do J=Jsq,Jeq do i=is,ie ; do_i(i) = (G%mask2dCv(i,J) > 0) ; enddo + ! When mixing down Eulerian current + Stokes drift add before calling solver + if (DoStokesMixing) then ; do k=1,nz ; do i=is,ie + if (do_i(i)) v(i,j,k) = v(i,j,k) + Waves%Us_y(i,j,k) + enddo ; enddo ; endif + if (associated(ADp%dv_dt_visc)) then ; do k=1,nz ; do i=is,ie ADp%dv_dt_visc(i,J,k) = v(i,J,k) enddo ; enddo ; endif @@ -411,12 +410,13 @@ subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, CS, & tauy_bot(i,J) = tauy_bot(i,J) + Rho0 * (Ray(i,k)*v(i,J,k)) enddo ; enddo ; endif endif - enddo ! end of v-component J loop - ! When mixing down Eulerian current + Stokes drift subtract after calling solver - if (DoStokesMixing) then ; do k=1,nz ; do J=Jsq,Jeq ; do i=Is,Ie - if (G%mask2dCv(i,J) > 0) v(i,J,k) = v(i,J,k) - Waves%Us_y(i,J,k) - enddo ; enddo ; enddo ; endif + ! When mixing down Eulerian current + Stokes drift subtract after calling solver + if (DoStokesMixing) then ; do k=1,nz ; do i=is,ie + if (do_i(i)) v(i,J,k) = v(i,J,k) - Waves%Us_y(i,J,k) + enddo ; enddo ; endif + + enddo ! end of v-component J loop call vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, CS) diff --git a/src/tracer/DOME_tracer.F90 b/src/tracer/DOME_tracer.F90 index 91b156751f..0354f90a51 100644 --- a/src/tracer/DOME_tracer.F90 +++ b/src/tracer/DOME_tracer.F90 @@ -15,7 +15,7 @@ module DOME_tracer use MOM_open_boundary, only : OBC_segment_type, register_segment_tracer use MOM_restart, only : MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS -use MOM_time_manager, only : time_type, get_time +use MOM_time_manager, only : time_type use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut use MOM_variables, only : surface diff --git a/src/tracer/ISOMIP_tracer.F90 b/src/tracer/ISOMIP_tracer.F90 index 40e8ef6db5..0707b54fb3 100644 --- a/src/tracer/ISOMIP_tracer.F90 +++ b/src/tracer/ISOMIP_tracer.F90 @@ -20,7 +20,7 @@ module ISOMIP_tracer use MOM_io, only : file_exists, MOM_read_data, slasher, vardesc, var_desc, query_vardesc use MOM_restart, only : MOM_restart_CS use MOM_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS -use MOM_time_manager, only : time_type, get_time +use MOM_time_manager, only : time_type use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut use MOM_variables, only : surface diff --git a/src/tracer/MOM_OCMIP2_CFC.F90 b/src/tracer/MOM_OCMIP2_CFC.F90 index e8c3387cea..ebff38508c 100644 --- a/src/tracer/MOM_OCMIP2_CFC.F90 +++ b/src/tracer/MOM_OCMIP2_CFC.F90 @@ -14,7 +14,7 @@ module MOM_OCMIP2_CFC use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : query_initialized, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS -use MOM_time_manager, only : time_type, get_time +use MOM_time_manager, only : time_type use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut use MOM_tracer_Z_init, only : tracer_Z_init diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90 index d06ffe0e2c..b373fc064a 100644 --- a/src/tracer/MOM_generic_tracer.F90 +++ b/src/tracer/MOM_generic_tracer.F90 @@ -35,7 +35,7 @@ module MOM_generic_tracer use MOM_spatial_means, only : global_area_mean use MOM_sponge, only : set_up_sponge_field, sponge_CS use MOM_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS - use MOM_time_manager, only : time_type, get_time, set_time + use MOM_time_manager, only : time_type, set_time use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_Z_init, only : tracer_Z_init @@ -498,7 +498,7 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, ! g_tracer=>CS%g_tracer_list do - if (_allocated(g_tracer%trunoff)) then + if (_ALLOCATED(g_tracer%trunoff)) then call g_tracer_get_alias(g_tracer,g_tracer_name) call g_tracer_get_pointer(g_tracer,g_tracer_name,'stf', stf_array) call g_tracer_get_pointer(g_tracer,g_tracer_name,'trunoff',trunoff_array) diff --git a/src/tracer/advection_test_tracer.F90 b/src/tracer/advection_test_tracer.F90 index 4ed395bac8..aeb1b3aae9 100644 --- a/src/tracer/advection_test_tracer.F90 +++ b/src/tracer/advection_test_tracer.F90 @@ -14,7 +14,7 @@ module advection_test_tracer use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : query_initialized, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS -use MOM_time_manager, only : time_type, get_time +use MOM_time_manager, only : time_type use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut use MOM_variables, only : surface diff --git a/src/tracer/boundary_impulse_tracer.F90 b/src/tracer/boundary_impulse_tracer.F90 index 7995b712e3..9b785fe41d 100644 --- a/src/tracer/boundary_impulse_tracer.F90 +++ b/src/tracer/boundary_impulse_tracer.F90 @@ -14,7 +14,7 @@ module boundary_impulse_tracer use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS -use MOM_time_manager, only : time_type, get_time +use MOM_time_manager, only : time_type use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut use MOM_tracer_Z_init, only : tracer_Z_init diff --git a/src/tracer/dye_example.F90 b/src/tracer/dye_example.F90 index a597b1fc8c..0e1b9a06b9 100644 --- a/src/tracer/dye_example.F90 +++ b/src/tracer/dye_example.F90 @@ -14,7 +14,7 @@ module regional_dyes use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : query_initialized, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS -use MOM_time_manager, only : time_type, get_time +use MOM_time_manager, only : time_type use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut use MOM_tracer_Z_init, only : tracer_Z_init diff --git a/src/tracer/dyed_obc_tracer.F90 b/src/tracer/dyed_obc_tracer.F90 index 2102f1cc71..af69a39c52 100644 --- a/src/tracer/dyed_obc_tracer.F90 +++ b/src/tracer/dyed_obc_tracer.F90 @@ -13,7 +13,7 @@ module dyed_obc_tracer use MOM_io, only : file_exists, MOM_read_data, slasher, vardesc, var_desc, query_vardesc use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : MOM_restart_CS -use MOM_time_manager, only : time_type, get_time +use MOM_time_manager, only : time_type use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut use MOM_variables, only : surface diff --git a/src/tracer/ideal_age_example.F90 b/src/tracer/ideal_age_example.F90 index 1f77bd639e..d7fcb53324 100644 --- a/src/tracer/ideal_age_example.F90 +++ b/src/tracer/ideal_age_example.F90 @@ -14,7 +14,7 @@ module ideal_age_example use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : query_initialized, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS -use MOM_time_manager, only : time_type, get_time +use MOM_time_manager, only : time_type, time_type_to_real use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut use MOM_tracer_Z_init, only : tracer_Z_init @@ -317,7 +317,6 @@ subroutine ideal_age_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, real :: sfc_val ! The surface value for the tracers. real :: Isecs_per_year ! The number of seconds in a year. real :: year ! The time in years. - integer :: secs, days ! Integer components of the time type. integer :: i, j, k, is, ie, js, je, nz, m is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -342,8 +341,7 @@ subroutine ideal_age_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, Isecs_per_year = 1.0 / (365.0*86400.0) ! Set the surface value of tracer 1 to increase exponentially ! with a 30 year time scale. - call get_time(CS%Time, secs, days) - year = (86400.0*days + real(secs)) * Isecs_per_year + year = time_type_to_real(CS%Time) * Isecs_per_year do m=1,CS%ntr if (CS%sfc_growth_rate(m) == 0.0) then diff --git a/src/tracer/oil_tracer.F90 b/src/tracer/oil_tracer.F90 index fd794aff0b..3b98c19a73 100644 --- a/src/tracer/oil_tracer.F90 +++ b/src/tracer/oil_tracer.F90 @@ -14,7 +14,7 @@ module oil_tracer use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : query_initialized, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS -use MOM_time_manager, only : time_type, get_time +use MOM_time_manager, only : time_type, time_type_to_real use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut use MOM_tracer_Z_init, only : tracer_Z_init @@ -334,7 +334,6 @@ subroutine oil_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified real :: Isecs_per_year = 1.0 / (365.0*86400.0) real :: year, h_total, ldecay - integer :: secs, days integer :: i, j, k, is, ie, js, je, nz, m, k_max is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -356,10 +355,7 @@ subroutine oil_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS enddo endif - ! Set the surface value of tracer 1 to increase exponentially - ! with a 30 year time scale. - call get_time(CS%Time, secs, days) - year = (86400.0*days + real(secs)) * Isecs_per_year + year = time_type_to_real(CS%Time) * Isecs_per_year ! Decay tracer (limit decay rate to 1/dt - just in case) do m=2,CS%ntr diff --git a/src/tracer/pseudo_salt_tracer.F90 b/src/tracer/pseudo_salt_tracer.F90 index fb0d38d86a..d9f4d3f682 100644 --- a/src/tracer/pseudo_salt_tracer.F90 +++ b/src/tracer/pseudo_salt_tracer.F90 @@ -16,7 +16,7 @@ module pseudo_salt_tracer use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : query_initialized, MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS -use MOM_time_manager, only : time_type, get_time +use MOM_time_manager, only : time_type use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut use MOM_tracer_Z_init, only : tracer_Z_init diff --git a/src/tracer/tracer_example.F90 b/src/tracer/tracer_example.F90 index 966fa07410..bf6b504658 100644 --- a/src/tracer/tracer_example.F90 +++ b/src/tracer/tracer_example.F90 @@ -14,7 +14,7 @@ module USER_tracer_example use MOM_open_boundary, only : ocean_OBC_type use MOM_restart, only : MOM_restart_CS use MOM_sponge, only : set_up_sponge_field, sponge_CS -use MOM_time_manager, only : time_type, get_time +use MOM_time_manager, only : time_type use MOM_tracer_registry, only : register_tracer, tracer_registry_type use MOM_variables, only : surface use MOM_verticalGrid, only : verticalGrid_type diff --git a/src/user/BFB_surface_forcing.F90 b/src/user/BFB_surface_forcing.F90 index 27168618be..edcdb002cf 100644 --- a/src/user/BFB_surface_forcing.F90 +++ b/src/user/BFB_surface_forcing.F90 @@ -12,7 +12,7 @@ module BFB_surface_forcing use MOM_grid, only : ocean_grid_type use MOM_io, only : file_exists, read_data use MOM_safe_alloc, only : safe_alloc_ptr -use MOM_time_manager, only : time_type, operator(+), operator(/), get_time +use MOM_time_manager, only : time_type, operator(+), operator(/) use MOM_tracer_flow_control, only : call_tracer_set_forcing use MOM_tracer_flow_control, only : tracer_flow_control_CS use MOM_variables, only : surface diff --git a/src/user/Kelvin_initialization.F90 b/src/user/Kelvin_initialization.F90 index eeda2e267f..8cf56a42ac 100644 --- a/src/user/Kelvin_initialization.F90 +++ b/src/user/Kelvin_initialization.F90 @@ -17,7 +17,7 @@ module Kelvin_initialization use MOM_open_boundary, only : OBC_DIRECTION_S, OBC_DIRECTION_W use MOM_open_boundary, only : OBC_registry_type use MOM_verticalGrid, only : verticalGrid_type -use MOM_time_manager, only : time_type, set_time, time_type_to_real +use MOM_time_manager, only : time_type, time_type_to_real implicit none ; private diff --git a/src/user/MOM_controlled_forcing.F90 b/src/user/MOM_controlled_forcing.F90 index c361a37176..2034a16bb4 100644 --- a/src/user/MOM_controlled_forcing.F90 +++ b/src/user/MOM_controlled_forcing.F90 @@ -18,8 +18,8 @@ module MOM_controlled_forcing use MOM_io, only : vardesc, var_desc use MOM_restart, only : register_restart_field, MOM_restart_CS use MOM_time_manager, only : time_type, operator(+), operator(/), operator(-) -use MOM_time_manager, only : get_time, get_date, set_time, set_date -use MOM_time_manager, only : time_type_to_real +use MOM_time_manager, only : get_date, set_date +use MOM_time_manager, only : time_type_to_real, real_to_time use MOM_variables, only : surface implicit none ; private @@ -121,7 +121,7 @@ subroutine apply_ctrl_forcing(SST_anom, SSS_anom, SSS_mean, virt_heat, virt_prec if (.not.associated(CS)) return if ((CS%num_cycle <= 0) .and. (.not.CS%do_integrated)) return - day_end = day_start + set_time(floor(dt+0.5)) + day_end = day_start + real_to_time(dt) do j=js,je ; do i=is,ie virt_heat(i,j) = 0.0 ; virt_precip(i,j) = 0.0 diff --git a/src/user/MOM_wave_interface.F90 b/src/user/MOM_wave_interface.F90 index 950fe4729d..c8ce37ad55 100644 --- a/src/user/MOM_wave_interface.F90 +++ b/src/user/MOM_wave_interface.F90 @@ -12,8 +12,7 @@ module MOM_wave_interface use MOM_grid, only : ocean_grid_type use MOM_verticalgrid, only : verticalGrid_type use MOM_safe_alloc, only : safe_alloc_ptr -use MOM_time_manager, only : time_type, operator(+), operator(/), get_time,& - time_type_to_real,real_to_time_type +use MOM_time_manager, only : time_type, operator(+), operator(/) use MOM_variables, only : thermo_var_ptrs, surface use data_override_mod, only : data_override_init, data_override diff --git a/src/user/SCM_CVmix_tests.F90 b/src/user/SCM_CVmix_tests.F90 index 2f2026c848..fca5ffa1d2 100644 --- a/src/user/SCM_CVmix_tests.F90 +++ b/src/user/SCM_CVmix_tests.F90 @@ -8,11 +8,10 @@ module SCM_CVMix_tests use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_forcing_type, only : forcing, mech_forcing use MOM_grid, only : ocean_grid_type -use MOM_verticalgrid, only: verticalGrid_type +use MOM_verticalgrid, only : verticalGrid_type use MOM_safe_alloc, only : safe_alloc_ptr -use MOM_time_manager, only : time_type, operator(+), operator(/), get_time,& - time_type_to_real -use MOM_variables, only : thermo_var_ptrs, surface +use MOM_time_manager, only : time_type, operator(+), operator(/), time_type_to_real +use MOM_variables, only : thermo_var_ptrs, surface implicit none ; private #include diff --git a/src/user/SCM_idealized_hurricane.F90 b/src/user/SCM_idealized_hurricane.F90 index f688c40ec6..2bb04b30f9 100644 --- a/src/user/SCM_idealized_hurricane.F90 +++ b/src/user/SCM_idealized_hurricane.F90 @@ -10,8 +10,7 @@ module SCM_idealized_hurricane use MOM_forcing_type, only : allocate_forcing_type, allocate_mech_forcing use MOM_grid, only : ocean_grid_type use MOM_safe_alloc, only : safe_alloc_ptr -use MOM_time_manager, only : time_type, operator(+), operator(/), get_time,& - time_type_to_real +use MOM_time_manager, only : time_type, operator(+), operator(/), time_type_to_real use MOM_variables, only : thermo_var_ptrs, surface use MOM_verticalGrid, only : verticalGrid_type diff --git a/src/user/dumbbell_surface_forcing.F90 b/src/user/dumbbell_surface_forcing.F90 index 0718ceb09c..d206914e2a 100644 --- a/src/user/dumbbell_surface_forcing.F90 +++ b/src/user/dumbbell_surface_forcing.F90 @@ -162,6 +162,7 @@ subroutine dumbbell_dynamic_forcing(state, fluxes, day, dt, G, CS) call get_time(day,isecs,idays) rdays = real(idays) + real(isecs)/8.64e4 + ! This could be: rdays = time_type_to_real(day)/8.64e4 ! Allocate and zero out the forcing arrays, as necessary. call safe_alloc_ptr(fluxes%p_surf, isd, ied, jsd, jed) diff --git a/src/user/dyed_channel_initialization.F90 b/src/user/dyed_channel_initialization.F90 index 133b5388cb..cb1b9a6b2f 100644 --- a/src/user/dyed_channel_initialization.F90 +++ b/src/user/dyed_channel_initialization.F90 @@ -11,7 +11,7 @@ module dyed_channel_initialization use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE, OBC_SIMPLE use MOM_open_boundary, only : OBC_segment_type, register_segment_tracer use MOM_open_boundary, only : OBC_registry_type, register_OBC -use MOM_time_manager, only : time_type, set_time, time_type_to_real +use MOM_time_manager, only : time_type, time_type_to_real use MOM_tracer_registry, only : tracer_registry_type, tracer_name_lookup use MOM_tracer_registry, only : tracer_type use MOM_variables, only : thermo_var_ptrs diff --git a/src/user/shelfwave_initialization.F90 b/src/user/shelfwave_initialization.F90 index 1640c9ec5a..9207830032 100644 --- a/src/user/shelfwave_initialization.F90 +++ b/src/user/shelfwave_initialization.F90 @@ -11,7 +11,7 @@ module shelfwave_initialization use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE, OBC_DIRECTION_W use MOM_open_boundary, only : OBC_segment_type, register_OBC use MOM_open_boundary, only : OBC_registry_type -use MOM_time_manager, only : time_type, set_time, time_type_to_real +use MOM_time_manager, only : time_type, time_type_to_real implicit none ; private diff --git a/src/user/soliton_initialization.F90 b/src/user/soliton_initialization.F90 index c9e7eec40e..e258b87bf1 100644 --- a/src/user/soliton_initialization.F90 +++ b/src/user/soliton_initialization.F90 @@ -53,7 +53,7 @@ subroutine soliton_initialize_thickness(h, G, GV) y = G%geoLatT(i,j)-y0 val3 = exp(-val1*x) val4 = val2*((2.0*val3/(1.0+(val3*val3)))**2) - h(i,j,k) = GV%m_to_H * (0.25*val4 * (6.0*y*y+3.0) * exp(-0.5*y*y)) + h(i,j,k) = GV%m_to_H * (0.25*val4 * (6.0*y*y+3.0) * exp(-0.5*y*y) + G%bathyT(i,j)) enddo enddo ; enddo diff --git a/src/user/supercritical_initialization.F90 b/src/user/supercritical_initialization.F90 index 6b10664d57..f12378c3d9 100644 --- a/src/user/supercritical_initialization.F90 +++ b/src/user/supercritical_initialization.F90 @@ -9,7 +9,7 @@ module supercritical_initialization use MOM_grid, only : ocean_grid_type use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE, OBC_SIMPLE, OBC_segment_type use MOM_verticalGrid, only : verticalGrid_type -use MOM_time_manager, only : time_type, set_time, time_type_to_real +use MOM_time_manager, only : time_type, time_type_to_real implicit none ; private diff --git a/src/user/tidal_bay_initialization.F90 b/src/user/tidal_bay_initialization.F90 index 7726dbf171..161ad25c11 100644 --- a/src/user/tidal_bay_initialization.F90 +++ b/src/user/tidal_bay_initialization.F90 @@ -13,7 +13,7 @@ module tidal_bay_initialization use MOM_open_boundary, only : OBC_segment_type, register_OBC use MOM_open_boundary, only : OBC_registry_type use MOM_verticalGrid, only : verticalGrid_type -use MOM_time_manager, only : time_type, set_time, time_type_to_real +use MOM_time_manager, only : time_type, time_type_to_real implicit none ; private diff --git a/src/user/user_revise_forcing.F90 b/src/user/user_revise_forcing.F90 index f2e381cc4a..d1be729734 100644 --- a/src/user/user_revise_forcing.F90 +++ b/src/user/user_revise_forcing.F90 @@ -10,7 +10,7 @@ module user_revise_forcing use MOM_grid, only : ocean_grid_type use MOM_io, only : file_exists, read_data use MOM_restart, only : register_restart_field, MOM_restart_CS -use MOM_time_manager, only : time_type, operator(+), operator(/), get_time +use MOM_time_manager, only : time_type, operator(+), operator(/) use MOM_tracer_flow_control, only : call_tracer_set_forcing use MOM_tracer_flow_control, only : tracer_flow_control_CS use MOM_variables, only : surface