From fa62eca556281397e3b3588e74bc0cca939e3724 Mon Sep 17 00:00:00 2001 From: abozec Date: Thu, 8 Apr 2021 17:08:46 -0400 Subject: [PATCH 1/9] initial set up for ODA incremental updates This is a basic framework to use NCODA increments on layer space for ODA incremental updates. --- src/core/MOM.F90 | 15 +- .../MOM_state_initialization.F90 | 123 +++- src/ocean_data_assim/MOM_oda_incupd.F90 | 655 ++++++++++++++++++ .../vertical/MOM_diabatic_driver.F90 | 55 +- 4 files changed, 839 insertions(+), 9 deletions(-) create mode 100644 src/ocean_data_assim/MOM_oda_incupd.F90 diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 9b94a96797..4deae909e8 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -132,6 +132,8 @@ module MOM ! ODA modules use MOM_oda_driver_mod, only : ODA_CS, oda, init_oda, oda_end use MOM_oda_driver_mod, only : set_prior_tracer, set_analysis_time, apply_oda_tracer_increments +use MOM_oda_incupd, only : oda_incupd_CS + ! Offline modules use MOM_offline_main, only : offline_transport_CS, offline_transport_init, update_offline_fields use MOM_offline_main, only : insert_offline_main, extract_offline_main, post_offline_convergence_diags @@ -362,6 +364,8 @@ module MOM type(sponge_CS), pointer :: sponge_CSp => NULL() !< Pointer to the layered-mode sponge control structure type(ALE_sponge_CS), pointer :: ALE_sponge_CSp => NULL() + !< Pointer to the oda incremental update control structure + type(oda_incupd_CS), pointer :: oda_incupd_CSp => NULL() !< Pointer to the ALE-mode sponge control structure type(ALE_CS), pointer :: ALE_CSp => NULL() !< Pointer to the Arbitrary Lagrangian Eulerian (ALE) vertical coordinate control structure @@ -1660,6 +1664,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & type(ocean_OBC_type), pointer :: OBC_in => NULL() type(sponge_CS), pointer :: sponge_in_CSp => NULL() type(ALE_sponge_CS), pointer :: ALE_sponge_in_CSp => NULL() + type(oda_incupd_CS),pointer :: oda_incupd_in_CSp => NULL() ! This include declares and sets the variable "version". # include "version_variable.h" @@ -2415,12 +2420,12 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call rotate_array(CS%frac_shelf_h, -turns, frac_shelf_in) call MOM_initialize_state(u_in, v_in, h_in, CS%tv, Time, G_in, GV, US, & param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & - sponge_in_CSp, ALE_sponge_in_CSp, OBC_in, Time_in, & + sponge_in_CSp, ALE_sponge_in_CSp, oda_incupd_in_CSp, OBC_in, Time_in, & frac_shelf_h=frac_shelf_in) else call MOM_initialize_state(u_in, v_in, h_in, CS%tv, Time, G_in, GV, US, & param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & - sponge_in_CSp, ALE_sponge_in_CSp, OBC_in, Time_in) + sponge_in_CSp, ALE_sponge_in_CSp, oda_incupd_in_CSp, OBC_in, Time_in) endif if (use_temperature) then @@ -2462,12 +2467,12 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & call ice_shelf_query(ice_shelf_CSp,G,CS%frac_shelf_h) call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, US, & param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & - CS%sponge_CSp, CS%ALE_sponge_CSp, CS%OBC, Time_in, & + CS%sponge_CSp, CS%ALE_sponge_CSp,CS%oda_incupd_CSp, CS%OBC, Time_in, & frac_shelf_h=CS%frac_shelf_h) else call MOM_initialize_state(CS%u, CS%v, CS%h, CS%tv, Time, G, GV, US, & param_file, dirs, restart_CSp, CS%ALE_CSp, CS%tracer_Reg, & - CS%sponge_CSp, CS%ALE_sponge_CSp, CS%OBC, Time_in) + CS%sponge_CSp, CS%ALE_sponge_CSp, CS%oda_incupd_CSp, CS%OBC, Time_in) endif endif @@ -2659,7 +2664,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & else call diabatic_driver_init(Time, G, GV, US, param_file, CS%use_ALE_algorithm, diag, & CS%ADp, CS%CDp, CS%diabatic_CSp, CS%tracer_flow_CSp, & - CS%sponge_CSp, CS%ALE_sponge_CSp) + CS%sponge_CSp, CS%ALE_sponge_CSp, CS%oda_incupd_CSp) endif if (associated(CS%sponge_CSp)) & diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 5050b6fce3..4650090967 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -94,6 +94,8 @@ module MOM_state_initialization use MOM_remapping, only : remapping_CS, initialize_remapping use MOM_remapping, only : remapping_core_h use MOM_horizontal_regridding, only : horiz_interp_and_extrap_tracer +use MOM_oda_incupd, only: oda_incupd_CS, initialize_oda_incupd +use MOM_oda_incupd, only: set_up_oda_incupd_field, set_up_oda_incupd_vel_field implicit none ; private @@ -114,7 +116,7 @@ module MOM_state_initialization !! conditions or by reading them from a restart (or saves) file. subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & restart_CS, ALE_CSp, tracer_Reg, sponge_CSp, & - ALE_sponge_CSp, OBC, Time_in, frac_shelf_h) + ALE_sponge_CSp, oda_incupd_CSp, OBC, Time_in, frac_shelf_h) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -140,6 +142,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & type(sponge_CS), pointer :: sponge_CSp !< The layerwise sponge control structure. type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< The ALE sponge control structure. type(ocean_OBC_type), pointer :: OBC !< The open boundary condition control structure. + type(oda_incupd_CS), pointer :: oda_incupd_CSp !< The oda_incupd control structure. type(time_type), optional, intent(in) :: Time_in !< Time at the start of the run segment. real, dimension(SZI_(G),SZJ_(G)), & optional, intent(in) :: frac_shelf_h !< The fraction of the grid cell covered @@ -157,7 +160,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & logical :: from_Z_file, useALE logical :: new_sim integer :: write_geom - logical :: use_temperature, use_sponge, use_OBC + logical :: use_temperature, use_sponge, use_OBC, use_oda_incupd logical :: use_EOS ! If true, density is calculated from T & S using an equation of state. logical :: depress_sfc ! If true, remove the mass that would be displaced ! by a large surface pressure by squeezing the column. @@ -615,6 +618,16 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & if (debug_OBC) call open_boundary_test_extern_h(G, GV, OBC, h) call callTree_leave('MOM_initialize_state()') + ! Data Assimilation with incremental update + call get_param(PF, mdl, "ODA_INCUPD", use_oda_incupd, & + "If true, oda incremental updates will be applied "//& + "everywhere in the domain.", default=.false.) + if (use_oda_incupd) then + call initialize_oda_incupd_file(G, GV, US, use_temperature, tv, u, v, & + PF, oda_incupd_CSp) + endif + + end subroutine MOM_initialize_state !> Reads the layer thicknesses or interface heights from a file. @@ -2022,6 +2035,112 @@ end function separate_idamp_for_uv end subroutine initialize_sponges_file +subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, u, v, param_file, oda_incupd_CSp) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + logical, intent(in) :: use_temperature !< If true, T & S are state variables. + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic + !! variables. + real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: u !< The zonal velocity that is being + !! initialized [L T-1 ~> m s-1] + real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + intent(in) :: v !< The meridional velocity that is being + !! initialized [L T-1 ~> m s-1] + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters. + type(oda_incupd_CS), pointer :: oda_incupd_CSp !< A pointer that is set to point to the control + !! structure for this module. + ! Local variables + real, allocatable, dimension(:,:,:) :: p ! The interface depth inc. [H ~> m or kg m-2]. + real, allocatable, dimension(:,:,:) :: tmp_tr ! A temporary array for reading oda fields + real, allocatable, dimension(:,:,:) :: tmp_u,tmp_v ! A temporary array for reading oda fields + + + integer :: i, j, k, is, ie, js, je, nz + integer :: isd, ied, jsd, jed + integer, dimension(4) :: siz + integer :: nz_data ! The size of the sponge source grid + character(len=40) :: tempinc_var, salinc_var, uinc_var, vinc_var, pinc_var + character(len=40) :: mdl = "initialize_oda_incupd_file" + character(len=200) :: inc_file, uv_inc_file ! Strings for filenames + character(len=200) :: filename, inputdir ! Strings for file/path and path. + + logical :: use_ALE ! True if ALE is being used, False if in layered mode + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + + + call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) + + call get_param(param_file, mdl, "ODA_INCUPD_FILE", inc_file, & + "The name of the file with the T,S,h increments.", & + fail_if_missing=.true.) + call get_param(param_file, mdl, "ODA_TEMPINC_VAR", tempinc_var, & + "The name of the potential temperature inc. variable in "//& + "ODA_INCUPD_FILE.", default="ptemp_inc") + call get_param(param_file, mdl, "ODA_SALTINC_VAR", salinc_var, & + "The name of the salinity inc. variable in "//& + "ODA_INCUPD_FILE.", default="sal_inc") + call get_param(param_file, mdl, "ODA_INTDINC_VAR", pinc_var, & + "The name of the interface depth inc. variable in "//& + "ODA_INCUPD_FILE.", default="p_inc") + call get_param(param_file, mdl, "ODA_INCUPD_UV_FILE", uv_inc_file, & + "The name of the file with the U,V increments.", & + default=inc_file) + call get_param(param_file, mdl, "ODA_UINC_VAR", uinc_var, & + "The name of the zonal vel. inc. variable in "//& + "ODA_INCUPD_FILE.", default="u_inc") + call get_param(param_file, mdl, "ODA_VINC_VAR", vinc_var, & + "The name of the meridional vel. inc. variable in "//& + "ODA_INCUPD_FILE.", default="v_inc") + + call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, do_not_log = .true.) + + ! Read in incremental update for tracers + filename = trim(inputdir)//trim(inc_file) + call log_param(param_file, mdl, "INPUTDIR/ODA_INCUPD_FILE", filename) + if (.not.file_exists(filename, G%Domain)) & + call MOM_error(FATAL, " initialize_oda_incupd: Unable to open "//trim(filename)) + + call field_size(filename,pinc_var,siz,no_domain=.true.) + if (siz(1) /= G%ieg-G%isg+1 .or. siz(2) /= G%jeg-G%jsg+1) & + call MOM_error(FATAL,"initialize_oda_incupd_file: Array size mismatch for oda data.") + nz_data = siz(3) + ! get h increments + allocate(p(isd:ied,jsd:jed,nz)) + call MOM_read_data(filename, pinc_var, p(:,:,:), G%Domain, scale=US%m_to_Z) + call initialize_oda_incupd( G, GV, param_file, oda_incupd_CSp, p, nz_data) + deallocate(p) + + ! get T and S increments + if (use_temperature) then + allocate(tmp_tr(isd:ied,jsd:jed,nz_data)) + call MOM_read_data(filename, tempinc_var, tmp_tr(:,:,:), G%Domain) + call set_up_oda_incupd_field(tmp_tr, G, GV, tv%T, oda_incupd_CSp) + call MOM_read_data(filename, salinc_var, tmp_tr(:,:,:), G%Domain) + call set_up_oda_incupd_field(tmp_tr, G, GV, tv%S, oda_incupd_CSp) + deallocate(tmp_tr) + endif + + ! get U and V increments + filename = trim(inputdir)//trim(uv_inc_file) + call log_param(param_file, mdl, "INPUTDIR/ODA_INCUPD_UV_FILE", filename) + if (.not.file_exists(filename, G%Domain)) & + call MOM_error(FATAL, " initialize_oda_incupd_uv: Unable to open "//trim(filename)) + allocate(tmp_u(G%IsdB:G%IedB,jsd:jed,nz_data)) + allocate(tmp_v(isd:ied,G%JsdB:G%JedB,nz_data)) + tmp_u(:,:,:) = 0.0 ; tmp_v(:,:,:) = 0.0 + call MOM_read_vector(filename, uinc_var, vinc_var, tmp_u, tmp_v, G%Domain,scale=US%m_s_to_L_T) + call set_up_oda_incupd_vel_field(tmp_u, tmp_v, G, GV, u, v, oda_incupd_CSp) + deallocate(tmp_u,tmp_v) + + + +end subroutine initialize_oda_incupd_file + !> This subroutine sets the 4 bottom depths at velocity points to be the !! maximum of the adjacent depths. subroutine set_velocity_depth_max(G) diff --git a/src/ocean_data_assim/MOM_oda_incupd.F90 b/src/ocean_data_assim/MOM_oda_incupd.F90 new file mode 100644 index 0000000000..a327b5a869 --- /dev/null +++ b/src/ocean_data_assim/MOM_oda_incupd.F90 @@ -0,0 +1,655 @@ +!> This module contains the routines used to apply incremental updates +!! from data assimilation. +!! +!! Applying incremental updates requires the following: +!! 1. initialize_oda_incupd +!! 2. set_up_oda_incupd_field (tracers) and set_up_oda_incupd_vel_field (vel) +!! 3. apply_oda_incupd +!! 5. oda_incupd_end (not being used for now) + +module MOM_oda_incupd + + +! This file is part of MOM6. See LICENSE.md for the license. +use MOM_array_transform, only: rotate_array +use MOM_coms, only : sum_across_PEs +use MOM_diag_mediator, only : post_data, query_averaging_enabled, register_diag_field +use MOM_diag_mediator, only : diag_ctrl +use MOM_domains, only : pass_var +use MOM_error_handler, only : MOM_error, FATAL, NOTE, WARNING, is_root_pe +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_grid, only : ocean_grid_type +use MOM_remapping, only : remapping_cs, remapping_core_h, initialize_remapping +use MOM_spatial_means, only : global_i_mean +use MOM_time_manager, only : time_type +use MOM_unit_scaling, only : unit_scale_type +use MOM_verticalGrid, only : verticalGrid_type + +use mpp_io_mod, only : mpp_get_axis_length +use mpp_io_mod, only : axistype + +implicit none ; private + +#include + + +! Publicly available functions +public set_up_oda_incupd_field, set_up_oda_incupd_vel_field +public initialize_oda_incupd, apply_oda_incupd, oda_incupd_end +!public rotate_oda_incupd, update_oda_incupd_field + +! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional +! consistency testing. These are noted in comments with units like Z, H, L, and T, along with +! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units +! vary with the Boussinesq approximation, the Boussinesq variant is given first. + +!> A structure for creating arrays of pointers to 3D arrays with extra gridding information +type :: p3d + integer :: id !< id for FMS external time interpolator + integer :: nz_data !< The number of vertical levels in the input field. + real, dimension(:,:,:), pointer :: mask_in => NULL() !< pointer to the data mask. + real, dimension(:,:,:), pointer :: p => NULL() !< pointer to the data. + real, dimension(:,:,:), pointer :: h => NULL() !< pointer to the data grid. +end type p3d + +!> A structure for creating arrays of pointers to 2D arrays with extra gridding information +type :: p2d + integer :: id !< id for FMS external time interpolator + integer :: nz_data !< The number of vertical levels in the input field + real, dimension(:,:), pointer :: mask_in => NULL()!< pointer to the data mask. + real, dimension(:,:), pointer :: p => NULL() !< pointer the data. + real, dimension(:,:), pointer :: h => NULL() !< pointer the data grid. +end type p2d + +!> oda incupd control structure +type, public :: oda_incupd_CS ; private + integer :: nz !< The total number of layers. + integer :: nz_data !< The total number of arbritary layers (used by older code). + integer :: isc !< The starting i-index of the computational domain at h. + integer :: iec !< The ending i-index of the computational domain at h. + integer :: jsc !< The starting j-index of the computational domain at h. + integer :: jec !< The ending j-index of the computational domain at h. + integer :: IscB !< The starting I-index of the computational domain at u/v. + integer :: IecB !< The ending I-index of the computational domain at u/v. + integer :: JscB !< The starting J-index of the computational domain at u/v. + integer :: JecB !< The ending J-index of the computational domain at h. + integer :: isd !< The starting i-index of the data domain at h. + integer :: ied !< The ending i-index of the data domain at h. + integer :: jsd !< The starting j-index of the data domain at h. + integer :: jed !< The ending j-index of the data domain at h. + integer :: IsdB !< The starting I-index of the data domain at u/v. + integer :: IedB !< The ending I-index of the data domain at u/v. + integer :: JsdB !< The starting J-index of the data domain at u/v. + integer :: JedB !< The ending J-index of the data domain at h. + integer :: fldno = 0 !< The number of fields which have already been + !! registered by calls to set_up_oda_incupd_field + + type(p3d) :: var(MAX_FIELDS_) !< Pointers to the fields that are being updated. + type(p3d) :: Ref_val(MAX_FIELDS_) !< The values to which the fields are updated. + type(p3d) :: Ref_oda_u !< The values to which the u-velocities are updated. + type(p3d) :: Ref_oda_v !< The values to which the v-velocities are updated. + type(p3d) :: var_u !< Pointer to the u velocities. that are being updated. + type(p3d) :: var_v !< Pointer to the v velocities. that are being updated. + type(p3d) :: Ref_h !< Grid on which reference data is provided (older code). + type(p3d) :: Ref_p !< grid of layer pressure increment + type(p3d) :: Ref_hu !< u-point grid on which reference data is provided (older code). + type(p3d) :: Ref_hv !< v-point grid on which reference data is provided (older code). + + + integer :: nstep_incupd !< number of time step for full update + integer :: ncount !< increment time step counter + type(remapping_cs) :: remap_cs !< Remapping parameters and work arrays + logical :: remap_answers_2018 !< If true, use the order of arithmetic and expressions that + !! recover the answers for remapping from the end of 2018. + !! Otherwise, use more robust forms of the same expressions. + logical :: hor_regrid_answers_2018 !< If true, use the order of arithmetic for horizonal regridding + !! that recovers the answers from the end of 2018. Otherwise, use + !! rotationally symmetric forms of the same expressions. + + logical :: incupdDataOngrid !< True if the incupd data are on the model horizontal grid + + logical :: reentrant_x !< grid is reentrant in the x direction + logical :: tripolar_N !< grid is folded at its north edge + +end type oda_incupd_CS + +contains + +!> This subroutine defined the number of time step for full update, stores the layer pressure +!! increments and initialize remap structure. +subroutine initialize_oda_incupd( G, GV, param_file, CS, data_p, nz_data) + + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + integer, intent(in) :: nz_data !< The total number of incr. input layers. + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for model parameter values. + type(oda_incupd_CS), pointer :: CS !< A pointer that is set to point to the control + !! structure for this module (in/out). + real, dimension(SZI_(G),SZJ_(G),nz_data), intent(in) :: data_p !< The int. depth incr. + !! [H ~> m or kg m-2]. + + +! This include declares and sets the variable "version". +#include "version_variable.h" + character(len=40) :: mdl = "MOM_oda" ! This module's name. + logical :: use_oda_incupd + logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries + logical :: default_2018_answers + integer :: i, j, k + real :: nhours_incupd, dt_baclin + character(len=256) :: mesg + character(len=10) :: remapScheme + if (associated(CS)) then + call MOM_error(WARNING, "initialize_oda_incupd_fixed called with an associated "// & + "control structure.") + return + endif + +! Set default, read and log parameters + call log_version(param_file, mdl, version, "") + call get_param(param_file, mdl, "ODA_INCUPD", use_oda_incupd, & + "If true, oda incremental updates will be applied "//& + "everywhere in the domain.", default=.false.) + + if (.not.use_oda_incupd) return + + allocate(CS) + + call get_param(param_file, mdl, "ODA_INCUPD_NHOURS", nhours_incupd, & + "Number of hours for full update (0=direct insertion).", & + default=3.0) + + call get_param(param_file, mdl, "DT_THERM", dt_baclin, & + "The thermodynamic and tracer advection time step."// & + "Ideally DT_THERM should be an integer multiple of DT"//& + "and less than the forcing or coupling time-step, unless"//& + "THERMO_SPANS_COUPLING is true, in which case DT_THERM"//& + "can be an integer multiple of the coupling timestep. By"//& + "default DT_THERM is set to DT.",& + default=300.0) + + call get_param(param_file, mdl, "REMAPPING_SCHEME", remapScheme, & + "This sets the reconstruction scheme used "//& + " for vertical remapping for all variables.", & + default="PLM", do_not_log=.true.) + + call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION", bndExtrapolation, & + "When defined, a proper high-order reconstruction "//& + "scheme is used within boundary cells rather "//& + "than PCM. E.g., if PPM is used for remapping, a "//& + "PPM reconstruction will also be used within boundary cells.", & + default=.false., do_not_log=.true.) + call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, & + "This sets the default value for the various _2018_ANSWERS parameters.", & + default=.false.) + call get_param(param_file, mdl, "REMAPPING_2018_ANSWERS", CS%remap_answers_2018, & + "If true, use the order of arithmetic and expressions that recover the "//& + "answers from the end of 2018. Otherwise, use updated and more robust "//& + "forms of the same expressions.", default=default_2018_answers) + call get_param(param_file, mdl, "ODA_INCUPD_DATA_ONGRID", CS%incupdDataOngrid, & + "When defined, the incoming oda_incupd data are "//& + "assumed to be on the model horizontal grid " , & + default=.true.) + call get_param(param_file, mdl, "REENTRANT_X", CS%reentrant_x, & + "If true, the domain is zonally reentrant.", default=.true.) + call get_param(param_file, mdl, "TRIPOLAR_N", CS%tripolar_N, & + "Use tripolar connectivity at the northern edge of the "//& + "domain. With TRIPOLAR_N, NIGLOBAL must be even.", default=.false.) + + CS%nz = GV%ke + CS%isc = G%isc ; CS%iec = G%iec ; CS%jsc = G%jsc ; CS%jec = G%jec + CS%isd = G%isd ; CS%ied = G%ied ; CS%jsd = G%jsd ; CS%jed = G%jed + CS%isdB = G%isdB ; CS%iedB = G%iedB; CS%jsdB = G%jsdB ; CS%jedB = G%jedB + + !! increments on horizontal grid + if (.not. CS%incupdDataOngrid) call MOM_error(FATAL,'initialize_oda_incupd: '// & + 'The oda_incupd code only applies ODA increments on the same horizontal grid. ') + + ! get number of timestep for full update + if (nhours_incupd == 0) then + CS%nstep_incupd = 1 !! direct insertion + else + CS%nstep_incupd = nhours_incupd*3600./dt_baclin + endif + write(mesg,'(i12)') CS%nstep_incupd + if (is_root_pe()) & + call MOM_error(NOTE,"initialize_oda_incupd: Number of Timestep of inc. update:"//& + trim(mesg)) + + ! initialize time counter + if (nhours_incupd == 0) then + CS%ncount=1 !! to assure putting the increment only once + else + CS%ncount=0 + endif + + ! get layer pressure increment p + CS%nz_data = nz_data + allocate(CS%Ref_p%p(G%isd:G%ied,G%jsd:G%jed,CS%nz_data)) + CS%Ref_p%p(:,:,:) = 0.0 ; + do j=G%jsc,G%jec; do i=G%isc,G%iec ; do k=1,CS%nz_data + CS%Ref_p%p(i,j,k) = data_p(i,j,k) + enddo; enddo ; enddo + +! Call the constructor for remapping control structure + call initialize_remapping(CS%remap_cs, remapScheme, boundary_extrapolation=bndExtrapolation, & + answers_2018=CS%remap_answers_2018) + + +end subroutine initialize_oda_incupd + + +!> This subroutine stores theincrements at h points for the variable +!! whose address is given by f_ptr. +subroutine set_up_oda_incupd_field(sp_val, G, GV, f_ptr, CS) + type(ocean_grid_type), intent(in) :: G !< Grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(oda_incupd_CS), pointer :: CS !< oda_incupd control structure (in/out). + real, dimension(SZI_(G),SZJ_(G),CS%nz_data), & + intent(in) :: sp_val !< increment field, it can have an + !! arbitrary number of layers. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + target, intent(in) :: f_ptr !< Pointer to the field to be updated + + integer :: i, j, k + character(len=256) :: mesg ! String for error messages + + if (.not.associated(CS)) return + + CS%fldno = CS%fldno + 1 + if (CS%fldno > MAX_FIELDS_) then + write(mesg,'("Increase MAX_FIELDS_ to at least ",I3," in MOM_memory.h or decrease & + &the number of fields increments in the call to & + &initialize_oda_incupd." )') CS%fldno + call MOM_error(FATAL,"set_up_oda_incupd_field: "//mesg) + endif + + ! stores the reference profile + CS%Ref_val(CS%fldno)%nz_data = CS%nz_data + allocate(CS%Ref_val(CS%fldno)%p(G%isd:G%ied,G%jsd:G%jed,CS%nz_data)) + CS%Ref_val(CS%fldno)%p(:,:,:) = 0.0 + do j=G%jsc,G%jec ; do i=G%isc,G%iec + do k=1,CS%nz_data + CS%Ref_val(CS%fldno)%p(i,j,k) = sp_val(i,j,k) + enddo + enddo ; enddo + + CS%var(CS%fldno)%p => f_ptr + +end subroutine set_up_oda_incupd_field + + +!> This subroutine stores the increments at u and v points for the variable +!! whose address is given by u_ptr and v_ptr. +subroutine set_up_oda_incupd_vel_field(u_val, v_val, G, GV, u_ptr, v_ptr, CS) + type(ocean_grid_type), intent(in) :: G !< Grid structure (in). + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(oda_incupd_CS), pointer :: CS !< oda incupd structure (in/out). + + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: u_val !< u increment, it has arbritary number of layers but + !! not to exceed the total number of model layers + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + intent(in) :: v_val !< v increment, it has arbritary number of layers but + !! not to exceed the number of model layers + real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u_ptr !< u ptr to the field to be assimilated + real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v_ptr !< v ptr to the field to be assimilated + + integer :: i, j, k + + if (.not.associated(CS)) return + + + ! stores the reference profile + allocate(CS%Ref_oda_u%p(G%isdB:G%iedB,G%jsd:G%jed,CS%nz_data)) + CS%Ref_oda_u%p(:,:,:) = 0.0 + do j=G%jsc,G%jec ; do i=G%iscB,G%iecB + do k=1,CS%nz_data + CS%Ref_oda_u%p(i,j,k) = u_val(i,j,k) + enddo + enddo ; enddo + CS%var_u%p => u_ptr + + allocate(CS%Ref_oda_v%p(G%isd:G%ied,G%jsdB:G%jedB,CS%nz_data)) + CS%Ref_oda_v%p(:,:,:) = 0.0 + do j=G%jscB,G%jecB ; do i=G%isc,G%iec + do k=1,CS%nz_data + CS%Ref_oda_v%p(i,j,k) = v_val(i,j,k) + enddo + enddo ; enddo + CS%var_v%p => v_ptr + +end subroutine set_up_oda_incupd_vel_field + + +!> This subroutine applies oda increments to layers thicknesses, temp, salt, U +!! and V everywhere . +subroutine apply_oda_incupd(h, dt, G, GV, US, CS) + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure (in). + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] (in) + real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]. + type(oda_incupd_CS), pointer :: CS !< A pointer to the control structure for this module + !! that is set by a previous call to initialize_oda_incupd (in). + + real :: m_to_Z ! A unit conversion factor from m to Z. + real, allocatable, dimension(:) :: tmp_val2 ! data values on the original grid + real, dimension(SZK_(GV)) :: tmp_val1 ! data values remapped to model grid + real, dimension(SZK_(GV)) :: h_col ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] + real, dimension(SZK_(GV)) :: hu, hv ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] + + real, allocatable, dimension(:,:,:) :: tmp_p ! A temporary array for int. depth field + real, allocatable, dimension(:,:,:) :: tmp_h ! A temporary array for thickness field + real, allocatable, dimension(:,:,:) :: tmp_u !< A temporary array for u inc. + real, allocatable, dimension(:,:,:) :: tmp_v !< A temporary array for v inc. + integer :: m, i, j, k, is, ie, js, je, nz, nz_data + integer :: isB, ieB, jsB, jeB + integer :: ncount ! time step counter + real :: q + real :: h_neglect, h_neglect_edge + character(len=256) :: mesg + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB + if (.not.associated(CS)) return + + ! update counter + if (CS%nstep_incupd /= 1) then + CS%ncount = CS%ncount+1 + if (CS%ncount == 1 .or. CS%ncount == CS%nstep_incupd+1) then + q = 0.5/CS%nstep_incupd !half increment on 1st and CS%step_incupd+1th step + else + q = 1.0/CS%nstep_incupd + endif !ncount + else + CS%ncount = CS%ncount+1 + q = 1.0/CS%nstep_incupd + endif + + ! no assimilation after CS%step_incupd+1 + if (CS%ncount > CS%nstep_incupd+1) then + if (CS%ncount == CS%nstep_incupd+2) then + if (is_root_pe()) call MOM_error(NOTE,"ended updating fields with increments. ") + endif !ncount==CS%nstep_incupd+2 + return + endif !ncount>CS%nstep_incupd+1 + + ! print out increments + write(mesg,'(i8)') CS%ncount + if (is_root_pe()) call MOM_error(NOTE,"updating fields with increments ncount:"//trim(mesg)) + write(mesg,'(f10.8)') q + if (is_root_pe()) call MOM_error(NOTE,"updating fields with weight q:"//trim(mesg)) + + + + if (.not.CS%remap_answers_2018) then + h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff + elseif (GV%Boussinesq) then + h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 + else + h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 + endif + + ! get interface depth from h + allocate(tmp_p(G%isd:G%ied,G%jsd:G%jed,nz+1)) ; tmp_p(:,:,:) = 0.0 + allocate(tmp_h(G%isd:G%ied,G%jsd:G%jed,nz )) ; tmp_h(:,:,:) = 0.0 + do j=js,je ; do i=is,ie + tmp_p(i,j,nz+1) = G%bathyT(i,j) + enddo ; enddo + do k=nz,1,-1 ; do j=js,je ; do i=is,ie + tmp_p(i,j,k) = tmp_p(i,j,k+1) - h(i,j,k) + enddo ; enddo ; enddo + + ! add increments to interface depths + do k=1,nz ; do j=js,je ; do i=is,ie + tmp_p(i,j,k) = tmp_p(i,j,k) + q*CS%Ref_p%p(i,j,k) + enddo ; enddo ; enddo + + ! get the new layer thickness + do k=1,nz ; do j=js,je ; do i=is,ie + tmp_h(i,j,k) = max(0.0, tmp_p(i,j,k+1) - tmp_p(i,j,k)) + enddo ; enddo ; enddo + + + ! add increments to tracers + tmp_val1(:)=0.0;h_col(:)=0.0 + do m=1,CS%fldno + nz_data = CS%Ref_val(m)%nz_data + allocate(tmp_val2(CS%Ref_val(m)%nz_data)) + do j=js,je ; do i=is,ie + tmp_val2(1:nz_data) = CS%Ref_val(m)%p(i,j,1:nz_data) + do k=1,nz + h_col(k)=h(i,j,k) + enddo + call remapping_core_h(CS%remap_cs, nz_data, tmp_h(i,j,1:nz_data), tmp_val2, & + nz, h_col, tmp_val1, h_neglect, h_neglect_edge) +! tmp_val1(1:nz)=tmp_val2(1:nz) + CS%var(m)%p(i,j,1:nz) = CS%var(m)%p(i,j,1:nz) + q*tmp_val1(1:nz) + enddo; enddo + deallocate(tmp_val2) + enddo + +! add increments to u + hu(:) = 0.0; h_col(:)=0 + allocate(tmp_val2(nz_data));tmp_val2(:)=0.0 + do j=js,je ; do i=isB,ieB + tmp_val2(1:nz_data) = 0.5*(CS%Ref_oda_u%p(i ,j,1:nz_data) + & + CS%Ref_oda_u%p(i+1,j,1:nz_data)) + do k=1,nz + hu(k) = 0.5 * (tmp_h(i,j,k) + tmp_h(i+1,j,k)) + h_col(k) = 0.5 * ( h(i,j,k) + h(i+1,j,k)) + enddo + call remapping_core_h(CS%remap_cs, nz_data, hu(1:nz_data), tmp_val2, & + nz, h_col, tmp_val1, h_neglect, h_neglect_edge) +! tmp_val1(1:nz)=tmp_val2(1:nz) + CS%var_u%p(i,j,1:nz) = CS%var_u%p(i,j,1:nz) + q*tmp_val1(1:nz) + enddo; enddo + deallocate(tmp_val2) + +! add increments to v + hv(:) = 0.0; h_col(:)=0 + allocate(tmp_val2(nz_data));tmp_val2(:)=0.0 + do j=jsB,jeB ; do i=is,ie + tmp_val2(1:nz_data) = 0.5*(CS%Ref_oda_v%p(i,j ,1:nz_data) + & + CS%Ref_oda_v%p(i,j+1,1:nz_data)) + do k=1,nz + hv(k) = 0.5 * (tmp_h(i,j,k) + tmp_h(i,j+1,k)) + h_col(k) = 0.5 * ( h(i,j,k) + h(i,j+1,k)) + enddo + call remapping_core_h(CS%remap_cs, nz_data, hv(1:nz_data), tmp_val2, & + nz, h_col, tmp_val1, h_neglect, h_neglect_edge) +! tmp_val1(1:nz)=tmp_val2(1:nz) + CS%var_v%p(i,j,1:nz) = CS%var_v%p(i,j,1:nz) + q*tmp_val1(1:nz) + enddo; enddo + deallocate(tmp_val2) + + +end subroutine apply_oda_incupd + +!> Rotate the incr. fields from the input to the model index map. +!> Not implemented yet +!subroutine rotate_oda_incupd(oda_incupd_in, G_in, oda_incupd, G, GV, turns, param_file) +! type(oda_incupd_CS), intent(in) :: oda_incupd_in !< The control structure for this module with the +! !! original grid rotation +! type(ocean_grid_type), intent(in) :: G_in !< The ocean's grid structure with the original rotation. +! type(oda_incupd_CS), pointer :: oda_incup !< A pointer to the control that will be set up with +! !! the new grid rotation +! type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure with the new rotation. +! type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure +! integer, intent(in) :: turns !< The number of 90-degree turns between grids +! type(param_file_type), intent(in) :: param_file !< A structure indicating the open file +! !! to parse for model parameter values. +! +! ! First part: Index construction +! ! 1. Reconstruct Iresttime(:,:) from sponge_in +! ! 2. rotate Iresttime(:,:) +! ! 3. Call initialize_oda_incupd using new grid and rotated Iresttime(:,:) +! ! All the index adjustment should follow from the Iresttime rotation +! +! real, dimension(:,:), allocatable :: Iresttime_in, Iresttime +! real, dimension(:,:,:), allocatable :: data_h_in, data_h +! real, dimension(:,:,:), allocatable :: sp_val_in, sp_val +! real, dimension(:,:,:), pointer :: sp_ptr => NULL() +! integer :: c, c_i, c_j +! integer :: k, nz_data +! integer :: n +! logical :: fixed_sponge +! +! fixed_sponge = .not. sponge_in%time_varying_sponges +! ! NOTE: nz_data is only conditionally set when fixed_sponge is true. +! +! allocate(Iresttime_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed)) +! allocate(Iresttime(G%isd:G%ied, G%jsd:G%jed)) +! Iresttime_in(:,:) = 0.0 +! +! if (fixed_sponge) then +! nz_data = sponge_in%nz_data +! allocate(data_h_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed, nz_data)) +! allocate(data_h(G%isd:G%ied, G%jsd:G%jed, nz_data)) +! data_h_in(:,:,:) = 0. +! endif +! +! ! Re-populate the 2D Iresttime and data_h arrays on the original grid +! do c=1,sponge_in%num_col +! c_i = sponge_in%col_i(c) +! c_j = sponge_in%col_j(c) +! Iresttime_in(c_i, c_j) = sponge_in%Iresttime_col(c) +! if (fixed_sponge) then +! do k = 1, nz_data +! data_h(c_i, c_j, k) = sponge_in%Ref_h%p(k,c) +! enddo +! endif +! enddo +! +! call rotate_array(Iresttime_in, turns, Iresttime) +! if (fixed_sponge) then +! call rotate_array(data_h_in, turns, data_h) +! call initialize_oda_incupd_fixed(Iresttime, G, GV, param_file, sponge, & +! data_h, nz_data) +! else +! call initialize_oda_incupd_varying(Iresttime, G, GV, param_file, sponge) +! endif +! +! deallocate(Iresttime_in) +! deallocate(Iresttime) +! if (fixed_sponge) then +! deallocate(data_h_in) +! deallocate(data_h) +! endif +! +! ! Second part: Provide rotated fields for which relaxation is applied +! +! sponge%fldno = sponge_in%fldno +! +! if (fixed_sponge) then +! allocate(sp_val_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed, nz_data)) +! allocate(sp_val(G%isd:G%ied, G%jsd:G%jed, nz_data)) +! endif +! +! do n=1,sponge_in%fldno +! ! Assume that tracers are pointers and are remapped in other functions(?) +! sp_ptr => sponge_in%var(n)%p +! if (fixed_sponge) then +! sp_val_in(:,:,:) = 0.0 +! do c = 1, sponge_in%num_col +! c_i = sponge_in%col_i(c) +! c_j = sponge_in%col_j(c) +! do k = 1, nz_data +! sp_val_in(c_i, c_j, k) = sponge_in%Ref_val(n)%p(k,c) +! enddo +! enddo +! +! call rotate_array(sp_val_in, turns, sp_val) +! +! ! NOTE: This points sp_val with the unrotated field. See note below. +! call set_up_oda_incupd_field(sp_val, G, GV, sp_ptr, sponge) +! +! deallocate(sp_val_in) +! else +! ! We don't want to repeat FMS init in set_up_oda_incupd_field_varying() +! ! (time_interp_external_init, init_external_field, etc), so we manually +! ! do a portion of this function below. +! sponge%Ref_val(n)%id = sponge_in%Ref_val(n)%id +! sponge%Ref_val(n)%num_tlevs = sponge_in%Ref_val(n)%num_tlevs +! +! nz_data = sponge_in%Ref_val(n)%nz_data +! sponge%Ref_val(n)%nz_data = nz_data +! +! allocate(sponge%Ref_val(n)%p(nz_data, sponge_in%num_col)) +! allocate(sponge%Ref_val(n)%h(nz_data, sponge_in%num_col)) +! sponge%Ref_val(n)%p(:,:) = 0.0 +! sponge%Ref_val(n)%h(:,:) = 0.0 +! +! ! TODO: There is currently no way to associate a generic field pointer to +! ! its rotated equivalent without introducing a new data structure which +! ! explicitly tracks the pairing. +! ! +! ! As a temporary fix, we store the pointer to the unrotated field in +! ! the rotated sponge, and use this reference to replace the pointer +! ! to the rotated field update_oda_incupd field. +! ! +! ! This makes a lot of unverifiable assumptions, and should not be +! ! considered the final solution. +! sponge%var(n)%p => sp_ptr +! endif +! enddo +! +! ! TODO: var_u and var_v sponge dampling is not yet supported. +! if (associated(sponge_in%var_u%p) .or. associated(sponge_in%var_v%p)) & +! call MOM_error(FATAL, "Rotation of ALE sponge velocities is not yet " & +! // "implemented.") +! +! ! Transfer any existing diag_CS reference pointer +! sponge%diag => sponge_in%diag +! +! ! NOTE: initialize_oda_incupd_* resolves remap_cs +!end subroutine rotate_oda_incupd + + +!> Scan the oda_incupd variables and replace a prescribed pointer to a new value. +!> Not implemented yet +! TODO: This function solely exists to replace field pointers in the sponge +! after rotation. This function is part of a temporary solution until +! something more robust is developed. +!subroutine update_oda_incupd_field(sponge, p_old, G, GV, p_new) +! type(oda_incupd_CS), pointer :: sponge !< A pointer to the control structure for this module +! !! that is set by a previous call to initialize_oda_incupd. +! real, dimension(:,:,:), & +! target, intent(in) :: p_old !< The previous array of target values +! type(ocean_grid_type), intent(in) :: G !< The updated ocean grid structure +! type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure +! real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & +! target, intent(in) :: p_new !< The new array of target values +! +! integer :: n +! +! do n=1,sponge%fldno +! if (associated(sponge%var(n)%p, p_old)) sponge%var(n)%p => p_new +! enddo +! +!end subroutine update_oda_incupd_field + + +! GMM: I could not find where sponge_end is being called, but I am keeping +! oda_incupd_end here so we can add that if needed. +!> This subroutine deallocates any memory associated with the oda_incupd module. +subroutine oda_incupd_end(CS) + type(oda_incupd_CS), pointer :: CS !< A pointer to the control structure that is + !! set by a previous call to initialize_oda_incupd. + + integer :: m + + if (.not.associated(CS)) return + + do m=1,CS%fldno + if (associated(CS%Ref_val(m)%p)) deallocate(CS%Ref_val(m)%p) + enddo + + deallocate(CS) + +end subroutine oda_incupd_end + +end module MOM_oda_incupd diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index c96edb785c..f7bb028aab 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -49,6 +49,7 @@ module MOM_diabatic_driver use MOM_CVMix_KPP, only : KPP_CS, KPP_init, KPP_compute_BLD, KPP_calculate use MOM_CVMix_KPP, only : KPP_end, KPP_get_BLD use MOM_CVMix_KPP, only : KPP_NonLocalTransport_temp, KPP_NonLocalTransport_saln +use MOM_oda_incupd, only : apply_oda_incupd, oda_incupd_CS use MOM_opacity, only : opacity_init, opacity_end, opacity_CS use MOM_opacity, only : absorbRemainingSW, optics_type, optics_nbands use MOM_open_boundary, only : ocean_OBC_type @@ -110,6 +111,8 @@ module MOM_diabatic_driver !! domain. The exact location and properties of !! those sponges are set by calls to !! initialize_sponge and set_up_sponge_field. + logical :: use_oda_incupd !!< If True, DA incremental update is + !! applied everywhere logical :: use_geothermal !< If true, apply geothermal heating. logical :: use_int_tides !< If true, use the code that advances a separate set !! of equations for the internal tide energy density. @@ -233,6 +236,7 @@ module MOM_diabatic_driver type(KPP_CS), pointer :: KPP_CSp => NULL() !< Control structure for a child module type(CVMix_conv_cs), pointer :: CVMix_conv_csp => NULL() !< Control structure for a child module type(diapyc_energy_req_CS), pointer :: diapyc_en_rec_CSp => NULL() !< Control structure for a child module + type(oda_incupd_CS), pointer :: oda_incupd_CSp => NULL() !< Control structure for a child module type(group_pass_type) :: pass_hold_eb_ea !< For group halo pass type(group_pass_type) :: pass_Kv !< For group halo pass @@ -251,7 +255,7 @@ module MOM_diabatic_driver integer :: id_clock_entrain, id_clock_mixedlayer, id_clock_set_diffusivity integer :: id_clock_tracers, id_clock_tridiag, id_clock_pass, id_clock_sponge integer :: id_clock_geothermal, id_clock_differential_diff, id_clock_remap -integer :: id_clock_kpp +integer :: id_clock_kpp, id_clock_oda_incupd !>@} contains @@ -1001,6 +1005,20 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim endif endif ! CS%use_sponge + ! Apply data assimilation incremental update -oda_incupd- + if (CS%use_oda_incupd .and. associated(CS%oda_incupd_CSp)) then + call MOM_mesg("Starting ODA_INCUPD legacy ", 5) + call cpu_clock_begin(id_clock_oda_incupd) + call apply_oda_incupd(h, dt, G, GV, US, CS%oda_incupd_CSp) + call cpu_clock_end(id_clock_oda_incupd) + if (CS%debug) then + call MOM_state_chksum("apply_oda_incupd ", u, v, h, G, GV, US, haloshift=0) + call MOM_thermovar_chksum("apply_oda_incupd ", tv, G) + endif + endif ! CS%use_oda_incupd + + + call disable_averaging(CS%diag) if (showCallTree) call callTree_leave("diabatic_ALE_legacy()") @@ -1484,6 +1502,19 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, endif endif ! CS%use_sponge + ! Apply data assimilation incremental update -oda_incupd- + if (CS%use_oda_incupd .and. associated(CS%oda_incupd_CSp)) then + call MOM_mesg("Starting ODA_INCUPD ", 5) + call cpu_clock_begin(id_clock_oda_incupd) + call apply_oda_incupd(h, dt, G, GV, US, CS%oda_incupd_CSp) + call cpu_clock_end(id_clock_oda_incupd) + if (CS%debug) then + call MOM_state_chksum("apply_oda_incupd ", u, v, h, G, GV, US, haloshift=0) + call MOM_thermovar_chksum("apply_oda_incupd ", tv, G) + endif + endif ! CS%use_oda_incupd + + call cpu_clock_begin(id_clock_pass) ! visc%Kv_slow is not in the group pass because it has larger vertical extent. if (associated(visc%Kv_slow)) & @@ -2267,6 +2298,19 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e endif endif ! CS%use_sponge + ! Apply data assimilation incremental update -oda_incupd- + if (CS%use_oda_incupd .and. associated(CS%oda_incupd_CSp)) then + call cpu_clock_begin(id_clock_oda_incupd) + call apply_oda_incupd(h, dt, G, GV, US, CS%oda_incupd_CSp) + call cpu_clock_end(id_clock_oda_incupd) + if (CS%debug) then + call MOM_state_chksum("apply_oda_incupd ", u, v, h, G, GV, US, haloshift=0) + call MOM_thermovar_chksum("apply_oda_incupd ", tv, G) + endif + endif ! CS%use_oda_incupd + + + ! Save the diapycnal mass fluxes as a diagnostic field. if (associated(CDp%diapyc_vel)) then !$OMP parallel do default(shared) @@ -2768,7 +2812,7 @@ end subroutine adiabatic_driver_init !> This routine initializes the diabatic driver module. subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, diag, & ADp, CDp, CS, tracer_flow_CSp, sponge_CSp, & - ALE_sponge_CSp) + ALE_sponge_CSp, oda_incupd_CSp) type(time_type), target :: Time !< model time type(ocean_grid_type), intent(inout) :: G !< model grid structure type(verticalGrid_type), intent(in) :: GV !< model vertical grid structure @@ -2784,6 +2828,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di !! tracer flow control module type(sponge_CS), pointer :: sponge_CSp !< pointer to the sponge module control structure type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< pointer to the ALE sponge module control structure + type(oda_incupd_CS), pointer :: oda_incupd_CSp !< pointer to the oda incupd module control structure real :: Kd ! A diffusivity used in the default for other tracer diffusivities, in MKS units [m2 s-1] integer :: num_mode @@ -2814,6 +2859,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di if (associated(tracer_flow_CSp)) CS%tracer_flow_CSp => tracer_flow_CSp if (associated(sponge_CSp)) CS%sponge_CSp => sponge_CSp if (associated(ALE_sponge_CSp)) CS%ALE_sponge_CSp => ALE_sponge_CSp + if (associated(oda_incupd_CSp)) CS%oda_incupd_CSp => oda_incupd_CSp CS%useALEalgorithm = useALEalgorithm CS%bulkmixedlayer = (GV%nkml > 0) @@ -2831,6 +2877,9 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di "The exact location and properties of those sponges are "//& "specified via calls to initialize_sponge and possibly "//& "set_up_sponge_field.", default=.false.) + call get_param(param_file, mdl, "ODA_INCUPD", CS%use_oda_incupd, & + "If true, oda incremental updates will be applied "//& + "everywhere in the domain.", default=.false.) call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, & "If true, temperature and salinity are used as state "//& "variables.", default=.true.) @@ -3322,6 +3371,8 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di id_clock_tracers = cpu_clock_id('(Ocean tracer_columns)', grain=CLOCK_MODULE_DRIVER+5) if (CS%use_sponge) & id_clock_sponge = cpu_clock_id('(Ocean sponges)', grain=CLOCK_MODULE) + if (CS%use_oda_incupd) & + id_clock_oda_incupd = cpu_clock_id('(Ocean inc. update data assimilation)', grain=CLOCK_MODULE) id_clock_tridiag = cpu_clock_id('(Ocean diabatic tridiag)', grain=CLOCK_ROUTINE) id_clock_pass = cpu_clock_id('(Ocean diabatic message passing)', grain=CLOCK_ROUTINE) id_clock_differential_diff = -1 ; if (CS%double_diffuse .and. .not.CS%use_CVMix_ddiff) & From 17d327bbd086b99a3d652d4c0089289bcb0b2e3c Mon Sep 17 00:00:00 2001 From: abozec Date: Wed, 28 Apr 2021 15:53:26 -0400 Subject: [PATCH 2/9] add options to incupd P_inc has been replaced by h of analysis state, add options to use full fields instead of increments and add diagnostics to get outputs of increments (for now to check what your increment looks like if full fields). --- src/core/MOM.F90 | 6 +- .../MOM_state_initialization.F90 | 21 +- src/ocean_data_assim/MOM_oda_incupd.F90 | 217 +++++++++++++----- 3 files changed, 175 insertions(+), 69 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 4deae909e8..95a89fdef2 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -132,7 +132,7 @@ module MOM ! ODA modules use MOM_oda_driver_mod, only : ODA_CS, oda, init_oda, oda_end use MOM_oda_driver_mod, only : set_prior_tracer, set_analysis_time, apply_oda_tracer_increments -use MOM_oda_incupd, only : oda_incupd_CS +use MOM_oda_incupd, only : oda_incupd_CS, init_oda_incupd_diags ! Offline modules use MOM_offline_main, only : offline_transport_CS, offline_transport_init, update_offline_fields @@ -2673,6 +2673,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & if (associated(CS%ALE_sponge_CSp)) & call init_ALE_sponge_diags(Time, G, diag, CS%ALE_sponge_CSp, US) + if (associated(CS%oda_incupd_CSp)) & + call init_oda_incupd_diags(Time, G, GV, diag, CS%oda_incupd_CSp, US) + + call tracer_advect_init(Time, G, US, param_file, diag, CS%tracer_adv_CSp) call tracer_hor_diff_init(Time, G, GV, US, param_file, diag, CS%tv%eqn_of_state, CS%diabatic_CSp, & CS%tracer_diff_CSp) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 4650090967..65308b2c25 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -2051,8 +2051,9 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, u, v, para type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters. type(oda_incupd_CS), pointer :: oda_incupd_CSp !< A pointer that is set to point to the control !! structure for this module. + ! Local variables - real, allocatable, dimension(:,:,:) :: p ! The interface depth inc. [H ~> m or kg m-2]. + real, allocatable, dimension(:,:,:) :: hoda ! The layer thk inc. and oda layer thk [H ~> m or kg m-2]. real, allocatable, dimension(:,:,:) :: tmp_tr ! A temporary array for reading oda fields real, allocatable, dimension(:,:,:) :: tmp_u,tmp_v ! A temporary array for reading oda fields @@ -2061,7 +2062,7 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, u, v, para integer :: isd, ied, jsd, jed integer, dimension(4) :: siz integer :: nz_data ! The size of the sponge source grid - character(len=40) :: tempinc_var, salinc_var, uinc_var, vinc_var, pinc_var + character(len=40) :: tempinc_var, salinc_var, uinc_var, vinc_var, h_var character(len=40) :: mdl = "initialize_oda_incupd_file" character(len=200) :: inc_file, uv_inc_file ! Strings for filenames character(len=200) :: filename, inputdir ! Strings for file/path and path. @@ -2084,9 +2085,9 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, u, v, para call get_param(param_file, mdl, "ODA_SALTINC_VAR", salinc_var, & "The name of the salinity inc. variable in "//& "ODA_INCUPD_FILE.", default="sal_inc") - call get_param(param_file, mdl, "ODA_INTDINC_VAR", pinc_var, & - "The name of the interface depth inc. variable in "//& - "ODA_INCUPD_FILE.", default="p_inc") + call get_param(param_file, mdl, "ODA_THK_VAR", h_var, & + "The name of the layer thickness variable in "//& + "ODA_INCUPD_FILE.", default="h") call get_param(param_file, mdl, "ODA_INCUPD_UV_FILE", uv_inc_file, & "The name of the file with the U,V increments.", & default=inc_file) @@ -2105,15 +2106,15 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, u, v, para if (.not.file_exists(filename, G%Domain)) & call MOM_error(FATAL, " initialize_oda_incupd: Unable to open "//trim(filename)) - call field_size(filename,pinc_var,siz,no_domain=.true.) + call field_size(filename,h_var,siz,no_domain=.true.) if (siz(1) /= G%ieg-G%isg+1 .or. siz(2) /= G%jeg-G%jsg+1) & call MOM_error(FATAL,"initialize_oda_incupd_file: Array size mismatch for oda data.") nz_data = siz(3) ! get h increments - allocate(p(isd:ied,jsd:jed,nz)) - call MOM_read_data(filename, pinc_var, p(:,:,:), G%Domain, scale=US%m_to_Z) - call initialize_oda_incupd( G, GV, param_file, oda_incupd_CSp, p, nz_data) - deallocate(p) + allocate(hoda(isd:ied,jsd:jed,nz)) + call MOM_read_data(filename, h_var , hoda(:,:,:), G%Domain, scale=US%m_to_Z) + call initialize_oda_incupd( G, GV, US, param_file, oda_incupd_CSp, hoda, nz_data) + deallocate(hoda) ! get T and S increments if (use_temperature) then diff --git a/src/ocean_data_assim/MOM_oda_incupd.F90 b/src/ocean_data_assim/MOM_oda_incupd.F90 index a327b5a869..8def6446b2 100644 --- a/src/ocean_data_assim/MOM_oda_incupd.F90 +++ b/src/ocean_data_assim/MOM_oda_incupd.F90 @@ -15,7 +15,7 @@ module MOM_oda_incupd use MOM_coms, only : sum_across_PEs use MOM_diag_mediator, only : post_data, query_averaging_enabled, register_diag_field use MOM_diag_mediator, only : diag_ctrl -use MOM_domains, only : pass_var +use MOM_domains, only : pass_var,pass_vector use MOM_error_handler, only : MOM_error, FATAL, NOTE, WARNING, is_root_pe use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type @@ -24,6 +24,7 @@ module MOM_oda_incupd use MOM_time_manager, only : time_type use MOM_unit_scaling, only : unit_scale_type use MOM_verticalGrid, only : verticalGrid_type +use MOM_verticalGrid, only : get_thickness_units use mpp_io_mod, only : mpp_get_axis_length use mpp_io_mod, only : axistype @@ -36,6 +37,7 @@ module MOM_oda_incupd ! Publicly available functions public set_up_oda_incupd_field, set_up_oda_incupd_vel_field public initialize_oda_incupd, apply_oda_incupd, oda_incupd_end +public init_oda_incupd_diags !public rotate_oda_incupd, update_oda_incupd_field ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional @@ -99,6 +101,7 @@ module MOM_oda_incupd integer :: nstep_incupd !< number of time step for full update integer :: ncount !< increment time step counter type(remapping_cs) :: remap_cs !< Remapping parameters and work arrays + logical :: oda_inc !< using increments or full fields logical :: remap_answers_2018 !< If true, use the order of arithmetic and expressions that !! recover the answers for remapping from the end of 2018. !! Otherwise, use more robust forms of the same expressions. @@ -111,22 +114,32 @@ module MOM_oda_incupd logical :: reentrant_x !< grid is reentrant in the x direction logical :: tripolar_N !< grid is folded at its north edge + ! for diagnostics + type(diag_ctrl), pointer :: diag ! This subroutine defined the number of time step for full update, stores the layer pressure !! increments and initialize remap structure. -subroutine initialize_oda_incupd( G, GV, param_file, CS, data_p, nz_data) +subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type integer, intent(in) :: nz_data !< The total number of incr. input layers. type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. type(oda_incupd_CS), pointer :: CS !< A pointer that is set to point to the control !! structure for this module (in/out). - real, dimension(SZI_(G),SZJ_(G),nz_data), intent(in) :: data_p !< The int. depth incr. + real, dimension(SZI_(G),SZJ_(G),nz_data), intent(in) :: data_h !< The ODA h !! [H ~> m or kg m-2]. @@ -159,6 +172,9 @@ subroutine initialize_oda_incupd( G, GV, param_file, CS, data_p, nz_data) call get_param(param_file, mdl, "ODA_INCUPD_NHOURS", nhours_incupd, & "Number of hours for full update (0=direct insertion).", & default=3.0) + call get_param(param_file, mdl, "ODA_INCUPD_INC", CS%oda_inc, & + "INCUPD files are increments and not full fields.", & + default=.true.) call get_param(param_file, mdl, "DT_THERM", dt_baclin, & "The thermodynamic and tracer advection time step."// & @@ -226,10 +242,10 @@ subroutine initialize_oda_incupd( G, GV, param_file, CS, data_p, nz_data) ! get layer pressure increment p CS%nz_data = nz_data - allocate(CS%Ref_p%p(G%isd:G%ied,G%jsd:G%jed,CS%nz_data)) - CS%Ref_p%p(:,:,:) = 0.0 ; + allocate(CS%Ref_h%p(G%isd:G%ied,G%jsd:G%jed,CS%nz_data)) + CS%Ref_h%p(:,:,:) = 0.0 ; do j=G%jsc,G%jec; do i=G%isc,G%iec ; do k=1,CS%nz_data - CS%Ref_p%p(i,j,k) = data_p(i,j,k) + CS%Ref_h%p(i,j,k) = data_h(i,j,k) enddo; enddo ; enddo ! Call the constructor for remapping control structure @@ -341,10 +357,18 @@ subroutine apply_oda_incupd(h, dt, G, GV, US, CS) real, dimension(SZK_(GV)) :: h_col ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] real, dimension(SZK_(GV)) :: hu, hv ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] - real, allocatable, dimension(:,:,:) :: tmp_p ! A temporary array for int. depth field - real, allocatable, dimension(:,:,:) :: tmp_h ! A temporary array for thickness field + !real, allocatable, dimension(:,:,:) :: tmp_p ! A temporary array for int. depth field + !real, allocatable, dimension(:,:,:) :: tmp_h ! A temporary array for thickness field real, allocatable, dimension(:,:,:) :: tmp_u !< A temporary array for u inc. real, allocatable, dimension(:,:,:) :: tmp_v !< A temporary array for v inc. + +! real, allocatable, dimension(:,:,:,:) :: var_inc !< increment fields if input is full fields +! real, allocatable, dimension(:,:,:) :: varu_inc !< increment U field if input is full fields +! real, allocatable, dimension(:,:,:) :: varv_inc !< increment V field if input is full fields + real, allocatable, dimension(:,:,:) :: h_obs !< h of increments + real, allocatable, dimension(:) :: hu_obs,hv_obs ! A column of thicknesses at h points [H ~> m or kg m-2] + + integer :: m, i, j, k, is, ie, js, je, nz, nz_data integer :: isB, ieB, jsB, jeB integer :: ncount ! time step counter @@ -369,6 +393,15 @@ subroutine apply_oda_incupd(h, dt, G, GV, US, CS) q = 1.0/CS%nstep_incupd endif + ! Diagnostics of increments, mostly used for debugging. + if (CS%id_u_oda_inc > 0) call post_data(CS%id_u_oda_inc, CS%Ref_oda_u%p, CS%diag) + if (CS%id_v_oda_inc > 0) call post_data(CS%id_v_oda_inc, CS%Ref_oda_v%p, CS%diag) + if (CS%id_h_oda_inc > 0) call post_data(CS%id_h_oda_inc, CS%Ref_h%p, CS%diag) + if (CS%id_T_oda_inc > 0) call post_data(CS%id_T_oda_inc, CS%Ref_val(1)%p,CS%diag) + if (CS%id_S_oda_inc > 0) call post_data(CS%id_S_oda_inc, CS%Ref_val(2)%p,CS%diag) + + + ! no assimilation after CS%step_incupd+1 if (CS%ncount > CS%nstep_incupd+1) then if (CS%ncount == CS%nstep_incupd+2) then @@ -377,6 +410,69 @@ subroutine apply_oda_incupd(h, dt, G, GV, US, CS) return endif !ncount>CS%nstep_incupd+1 + ! get h_obs + nz_data = CS%Ref_val(1)%nz_data + allocate(h_obs(G%isd:G%ied,G%jsd:G%jed,nz_data)) ; h_obs(:,:,:)=0.0 + do k=1,nz_data ; do j=js,je ; do i=is,ie + h_obs(i,j,k) = CS%Ref_h%p(i,j,k) + enddo ; enddo ; enddo + call pass_var(h_obs,G%Domain) + + + if (CS%oda_inc) then ! input are increments + if (is_root_pe()) call MOM_error(NOTE,"incupd using increments fields ") + + else ! inputs are full fields + if (q == 1.0 .or. CS%ncount == 1) then + if (is_root_pe()) call MOM_error(NOTE,"incupd using full fields ") + + ! remap t,s (on h_init) to h_obs to get increment + tmp_val1(:)=0.0 + do m=1,CS%fldno + allocate(tmp_val2(nz)) + do j=js,je ; do i=is,ie + tmp_val2(1:nz) = CS%var(m)%p(i,j,1:nz) + call remapping_core_h(CS%remap_cs, nz , h(i,j,1:nz ), tmp_val2, & + nz_data, h_obs(i,j,1:nz_data), tmp_val1, & + h_neglect, h_neglect_edge) + CS%Ref_val(m)%p(i,j,1:nz_data) = CS%Ref_val(m)%p(i,j,1:nz_data) - tmp_val1(1:nz_data) + enddo; enddo + deallocate(tmp_val2) + enddo + + ! remap u to h_obs to get increment + allocate(hu_obs(nz_data)) ; hu_obs(:)=0.0 ; hu(:) = 0.0; + allocate(tmp_val2(nz));tmp_val2(:)=0.0 + do j=js,je ; do i=isB,ieB + tmp_val2(1:nz) = CS%var_u%p(i ,j,1:nz) + hu(1:nz) = 0.5*( h(i,j,1:nz )+ h(i+1,j,1:nz )) + hu_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i+1,j,1:nz_data)) + call remapping_core_h(CS%remap_cs, nz , hu(1:nz ), tmp_val2, & + nz_data, hu_obs(1:nz_data), tmp_val1, & + h_neglect, h_neglect_edge) + CS%Ref_oda_u%p(i,j,1:nz_data) = CS%Ref_oda_u%p(i,j,1:nz_data) - tmp_val1(1:nz_data) + enddo; enddo + deallocate(tmp_val2,hu_obs) + + ! remap v to h_obs to get increment + allocate(hv_obs(nz_data)) ; hv_obs(:)=0.0 ; hv(:) = 0.0; + allocate(tmp_val2(nz));tmp_val2(:)=0.0 + do j=jsB,jeB ; do i=is,ie + tmp_val2(1:nz) = CS%var_v%p(i ,j,1:nz) + hv(1:nz) = 0.5*( h(i,j,1:nz )+ h(i,j+1,1:nz )) + hv_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i,j+1,1:nz_data)) + call remapping_core_h(CS%remap_cs, nz , hv(1:nz ), tmp_val2, & + nz_data, hv_obs(1:nz_data), tmp_val1, & + h_neglect, h_neglect_edge) + CS%Ref_oda_v%p(i,j,1:nz_data) = CS%Ref_oda_v%p(i,j,1:nz_data) - tmp_val1(1:nz_data) + enddo; enddo + deallocate(tmp_val2,hv_obs) + + endif ! not q == 1 + + endif ! not oda_inc + + ! print out increments write(mesg,'(i8)') CS%ncount if (is_root_pe()) call MOM_error(NOTE,"updating fields with increments ncount:"//trim(mesg)) @@ -393,82 +489,87 @@ subroutine apply_oda_incupd(h, dt, G, GV, US, CS) h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 endif - ! get interface depth from h - allocate(tmp_p(G%isd:G%ied,G%jsd:G%jed,nz+1)) ; tmp_p(:,:,:) = 0.0 - allocate(tmp_h(G%isd:G%ied,G%jsd:G%jed,nz )) ; tmp_h(:,:,:) = 0.0 - do j=js,je ; do i=is,ie - tmp_p(i,j,nz+1) = G%bathyT(i,j) - enddo ; enddo - do k=nz,1,-1 ; do j=js,je ; do i=is,ie - tmp_p(i,j,k) = tmp_p(i,j,k+1) - h(i,j,k) - enddo ; enddo ; enddo - - ! add increments to interface depths - do k=1,nz ; do j=js,je ; do i=is,ie - tmp_p(i,j,k) = tmp_p(i,j,k) + q*CS%Ref_p%p(i,j,k) - enddo ; enddo ; enddo - - ! get the new layer thickness - do k=1,nz ; do j=js,je ; do i=is,ie - tmp_h(i,j,k) = max(0.0, tmp_p(i,j,k+1) - tmp_p(i,j,k)) - enddo ; enddo ; enddo - ! add increments to tracers - tmp_val1(:)=0.0;h_col(:)=0.0 + tmp_val1(:)=0.0 do m=1,CS%fldno - nz_data = CS%Ref_val(m)%nz_data allocate(tmp_val2(CS%Ref_val(m)%nz_data)) do j=js,je ; do i=is,ie tmp_val2(1:nz_data) = CS%Ref_val(m)%p(i,j,1:nz_data) - do k=1,nz - h_col(k)=h(i,j,k) - enddo - call remapping_core_h(CS%remap_cs, nz_data, tmp_h(i,j,1:nz_data), tmp_val2, & - nz, h_col, tmp_val1, h_neglect, h_neglect_edge) -! tmp_val1(1:nz)=tmp_val2(1:nz) + call remapping_core_h(CS%remap_cs, nz_data, h_obs(i,j,1:nz_data), tmp_val2, & + nz , h(i,j,1:nz ), tmp_val1, & + h_neglect, h_neglect_edge) CS%var(m)%p(i,j,1:nz) = CS%var(m)%p(i,j,1:nz) + q*tmp_val1(1:nz) enddo; enddo deallocate(tmp_val2) enddo + ! bound salinity values ! check if it is correct to do that or if it hides + ! other problems ... + do j=js,je ; do i=is,ie + do k=1,nz + CS%var(2)%p(i,j,k) = max(0.0 , CS%var(2)%p(i,j,k)) + enddo + enddo ; enddo ! add increments to u - hu(:) = 0.0; h_col(:)=0 + allocate(hu_obs(nz_data)) ; hu_obs(:)=0.0 ; hu(:) = 0.0; allocate(tmp_val2(nz_data));tmp_val2(:)=0.0 do j=js,je ; do i=isB,ieB - tmp_val2(1:nz_data) = 0.5*(CS%Ref_oda_u%p(i ,j,1:nz_data) + & - CS%Ref_oda_u%p(i+1,j,1:nz_data)) - do k=1,nz - hu(k) = 0.5 * (tmp_h(i,j,k) + tmp_h(i+1,j,k)) - h_col(k) = 0.5 * ( h(i,j,k) + h(i+1,j,k)) - enddo - call remapping_core_h(CS%remap_cs, nz_data, hu(1:nz_data), tmp_val2, & - nz, h_col, tmp_val1, h_neglect, h_neglect_edge) -! tmp_val1(1:nz)=tmp_val2(1:nz) + tmp_val2(1:nz_data) = CS%Ref_oda_u%p(i,j,1:nz_data) + hu_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i+1,j,1:nz_data)) + hu(1:nz) = 0.5*( h(i,j,1:nz )+ h(i+1,j,1:nz )) + call remapping_core_h(CS%remap_cs, nz_data, hu_obs(1:nz_data), tmp_val2, & + nz , hu(1:nz ), tmp_val1, & + h_neglect, h_neglect_edge) CS%var_u%p(i,j,1:nz) = CS%var_u%p(i,j,1:nz) + q*tmp_val1(1:nz) enddo; enddo - deallocate(tmp_val2) + deallocate(tmp_val2,hu_obs) ! add increments to v - hv(:) = 0.0; h_col(:)=0 + allocate(hv_obs(nz_data)) ; hv_obs(:)=0.0 ; hv(:) = 0.0; allocate(tmp_val2(nz_data));tmp_val2(:)=0.0 do j=jsB,jeB ; do i=is,ie - tmp_val2(1:nz_data) = 0.5*(CS%Ref_oda_v%p(i,j ,1:nz_data) + & - CS%Ref_oda_v%p(i,j+1,1:nz_data)) - do k=1,nz - hv(k) = 0.5 * (tmp_h(i,j,k) + tmp_h(i,j+1,k)) - h_col(k) = 0.5 * ( h(i,j,k) + h(i,j+1,k)) - enddo - call remapping_core_h(CS%remap_cs, nz_data, hv(1:nz_data), tmp_val2, & - nz, h_col, tmp_val1, h_neglect, h_neglect_edge) -! tmp_val1(1:nz)=tmp_val2(1:nz) + tmp_val2(1:nz_data) = CS%Ref_oda_v%p(i,j,1:nz_data) + hv_obs(1:nz_data) = 0.5 * (h_obs(i,j,1:nz_data)+h_obs(i,j+1,1:nz_data)) + hv(1:nz) = 0.5 * ( h(i,j,1:nz )+ h(i,j+1,1:nz )) + call remapping_core_h(CS%remap_cs, nz_data, hv_obs(1:nz_data), tmp_val2, & + nz , hv(1:nz ), tmp_val1, & + h_neglect, h_neglect_edge) CS%var_v%p(i,j,1:nz) = CS%var_v%p(i,j,1:nz) + q*tmp_val1(1:nz) enddo; enddo - deallocate(tmp_val2) + deallocate(tmp_val2,hv_obs) end subroutine apply_oda_incupd +!> Initialize diagnostics for the oda_incupd module. +subroutine init_oda_incupd_diags(Time, G, GV, diag, CS, US) + type(time_type), target, intent(in) :: Time !< The current model time + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic + !! output. + type(oda_incupd_CS), pointer :: CS !< ALE sponge control structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling + + if (.not.associated(CS)) return + + CS%diag => diag + ! These diagnostics of the state variables increments,useful for debugging the + ! ODA code. + CS%id_u_oda_inc = register_diag_field('ocean_model', 'u_oda_inc', diag%axesCuL, Time, & + 'Zonal velocity ODA inc.', 'm s-1') + CS%id_v_oda_inc = register_diag_field('ocean_model', 'v_oda_inc', diag%axesCvL, Time, & + 'Meridional velocity ODA inc.', 'm s-1') + CS%id_h_oda_inc = register_diag_field('ocean_model', 'h_oda_inc', diag%axesTL, Time, & + 'Layer Thickness ODA inc.', get_thickness_units(GV)) + CS%id_T_oda_inc = register_diag_field('ocean_model', 'T_oda_inc', diag%axesTL, Time, & + 'Temperature ODA inc.', 'degC') + CS%id_S_oda_inc = register_diag_field('ocean_model', 'S_oda_inc', diag%axesTL, Time, & + 'Salinity ODA inc.', 'PSU') + +end subroutine init_oda_incupd_diags + !> Rotate the incr. fields from the input to the model index map. !> Not implemented yet !subroutine rotate_oda_incupd(oda_incupd_in, G_in, oda_incupd, G, GV, turns, param_file) From fe39cbf62829d27d6768ee3a0c83ef9875134b98 Mon Sep 17 00:00:00 2001 From: abozec Date: Tue, 4 May 2021 14:01:31 -0400 Subject: [PATCH 3/9] Change increments initialization, update time steps and diagnostics Increments calculation moved at initialization, updates applied over nstep instead of nstep+1 and increments diagnostics has been added --- .../MOM_state_initialization.F90 | 26 ++- src/ocean_data_assim/MOM_oda_incupd.F90 | 219 ++++++++++-------- 2 files changed, 138 insertions(+), 107 deletions(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 65308b2c25..a46304cfe4 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -11,7 +11,7 @@ module MOM_state_initialization use MOM_cpu_clock, only : CLOCK_ROUTINE, CLOCK_LOOP use MOM_domains, only : pass_var, pass_vector, sum_across_PEs, broadcast use MOM_domains, only : root_PE, To_All, SCALAR_PAIR, CGRID_NE, AGRID -use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe +use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, NOTE, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, read_param, log_param, param_file_type use MOM_file_parser, only : log_version @@ -96,6 +96,7 @@ module MOM_state_initialization use MOM_horizontal_regridding, only : horiz_interp_and_extrap_tracer use MOM_oda_incupd, only: oda_incupd_CS, initialize_oda_incupd use MOM_oda_incupd, only: set_up_oda_incupd_field, set_up_oda_incupd_vel_field +use MOM_oda_incupd, only: get_oda_increments implicit none ; private @@ -623,7 +624,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & "If true, oda incremental updates will be applied "//& "everywhere in the domain.", default=.false.) if (use_oda_incupd) then - call initialize_oda_incupd_file(G, GV, US, use_temperature, tv, u, v, & + call initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, & PF, oda_incupd_CSp) endif @@ -2035,13 +2036,16 @@ end function separate_idamp_for_uv end subroutine initialize_sponges_file -subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, u, v, param_file, oda_incupd_CSp) +subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, param_file, oda_incupd_CSp) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type logical, intent(in) :: use_temperature !< If true, T & S are state variables. type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic !! variables. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Layer thickness [H ~> m or kg m-2] (in) + real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & intent(in) :: u !< The zonal velocity that is being !! initialized [L T-1 ~> m s-1] @@ -2057,11 +2061,13 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, u, v, para real, allocatable, dimension(:,:,:) :: tmp_tr ! A temporary array for reading oda fields real, allocatable, dimension(:,:,:) :: tmp_u,tmp_v ! A temporary array for reading oda fields - integer :: i, j, k, is, ie, js, je, nz integer :: isd, ied, jsd, jed + integer, dimension(4) :: siz integer :: nz_data ! The size of the sponge source grid + logical :: oda_inc ! input files are increments (true) or full fields (false) + character(len=40) :: tempinc_var, salinc_var, uinc_var, vinc_var, h_var character(len=40) :: mdl = "initialize_oda_incupd_file" character(len=200) :: inc_file, uv_inc_file ! Strings for filenames @@ -2072,13 +2078,15 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, u, v, para is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") inputdir = slasher(inputdir) call get_param(param_file, mdl, "ODA_INCUPD_FILE", inc_file, & "The name of the file with the T,S,h increments.", & fail_if_missing=.true.) + call get_param(param_file, mdl, "ODA_INCUPD_INC", oda_inc, & + "INCUPD files are increments and not full fields.", & + default=.true.) call get_param(param_file, mdl, "ODA_TEMPINC_VAR", tempinc_var, & "The name of the potential temperature inc. variable in "//& "ODA_INCUPD_FILE.", default="ptemp_inc") @@ -2138,7 +2146,13 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, u, v, para call set_up_oda_incupd_vel_field(tmp_u, tmp_v, G, GV, u, v, oda_incupd_CSp) deallocate(tmp_u,tmp_v) - + ! calculate increments if full fields + if (oda_inc) then ! input are increments + if (is_root_pe()) call MOM_error(NOTE,"incupd using increments fields ") + else ! inputs are full fields + if (is_root_pe()) call MOM_error(NOTE,"incupd using full fields ") + call get_oda_increments(h, G, GV, US, oda_incupd_CSp) + endif ! not oda_inc end subroutine initialize_oda_incupd_file diff --git a/src/ocean_data_assim/MOM_oda_incupd.F90 b/src/ocean_data_assim/MOM_oda_incupd.F90 index 8def6446b2..6549b95bee 100644 --- a/src/ocean_data_assim/MOM_oda_incupd.F90 +++ b/src/ocean_data_assim/MOM_oda_incupd.F90 @@ -37,7 +37,7 @@ module MOM_oda_incupd ! Publicly available functions public set_up_oda_incupd_field, set_up_oda_incupd_vel_field public initialize_oda_incupd, apply_oda_incupd, oda_incupd_end -public init_oda_incupd_diags +public init_oda_incupd_diags,get_oda_increments !public rotate_oda_incupd, update_oda_incupd_field ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional @@ -101,7 +101,6 @@ module MOM_oda_incupd integer :: nstep_incupd !< number of time step for full update integer :: ncount !< increment time step counter type(remapping_cs) :: remap_cs !< Remapping parameters and work arrays - logical :: oda_inc !< using increments or full fields logical :: remap_answers_2018 !< If true, use the order of arithmetic and expressions that !! recover the answers for remapping from the end of 2018. !! Otherwise, use more robust forms of the same expressions. @@ -172,9 +171,6 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data) call get_param(param_file, mdl, "ODA_INCUPD_NHOURS", nhours_incupd, & "Number of hours for full update (0=direct insertion).", & default=3.0) - call get_param(param_file, mdl, "ODA_INCUPD_INC", CS%oda_inc, & - "INCUPD files are increments and not full fields.", & - default=.true.) call get_param(param_file, mdl, "DT_THERM", dt_baclin, & "The thermodynamic and tracer advection time step."// & @@ -234,11 +230,7 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data) trim(mesg)) ! initialize time counter - if (nhours_incupd == 0) then - CS%ncount=1 !! to assure putting the increment only once - else - CS%ncount=0 - endif + CS%ncount=0 ! get layer pressure increment p CS%nz_data = nz_data @@ -248,7 +240,7 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data) CS%Ref_h%p(i,j,k) = data_h(i,j,k) enddo; enddo ; enddo -! Call the constructor for remapping control structure + ! Call the constructor for remapping control structure call initialize_remapping(CS%remap_cs, remapScheme, boundary_extrapolation=bndExtrapolation, & answers_2018=CS%remap_answers_2018) @@ -338,6 +330,92 @@ subroutine set_up_oda_incupd_vel_field(u_val, v_val, G, GV, u_ptr, v_ptr, CS) end subroutine set_up_oda_incupd_vel_field +subroutine get_oda_increments(h, G, GV, US, CS) + + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure (in). + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Layer thickness [H ~> m or kg m-2] (in) + type(oda_incupd_CS), pointer :: CS !< A pointer to the control structure for this module + !! that is set by a previous call to initialize_oda_incupd (in). + + + real, allocatable, dimension(:) :: tmp_val2 ! data values on the original grid + real, dimension(SZK_(GV)) :: tmp_val1 ! data values remapped to model grid + real, allocatable, dimension(:,:,:) :: h_obs !< h of increments + real, allocatable, dimension(:) :: hu_obs,hv_obs ! A column of thicknesses at h points [H ~> m or kg m-2] + real, dimension(SZK_(GV)) :: hu, hv ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] + + + integer :: m, i, j, k, is, ie, js, je, nz, nz_data + integer :: isB, ieB, jsB, jeB + real :: h_neglect, h_neglect_edge + character(len=256) :: mesg + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB + if (.not.associated(CS)) return + + if (.not.CS%remap_answers_2018) then + h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff + elseif (GV%Boussinesq) then + h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 + else + h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 + endif + + ! get h_obs + nz_data = CS%Ref_val(1)%nz_data + allocate(h_obs(G%isd:G%ied,G%jsd:G%jed,nz_data)) ; h_obs(:,:,:)=0.0 + do k=1,nz_data ; do j=js,je ; do i=is,ie + h_obs(i,j,k) = CS%Ref_h%p(i,j,k) + enddo ; enddo ; enddo + call pass_var(h_obs,G%Domain) + + ! remap t,s (on h_init) to h_obs to get increment + tmp_val1(:)=0.0 + do m=1,CS%fldno + allocate(tmp_val2(nz)) + do j=js,je ; do i=is,ie + tmp_val2(1:nz) = CS%var(m)%p(i,j,1:nz) + call remapping_core_h(CS%remap_cs, nz , h(i,j,1:nz ), tmp_val2, & + nz_data, h_obs(i,j,1:nz_data), tmp_val1, & + h_neglect, h_neglect_edge) + CS%Ref_val(m)%p(i,j,1:nz_data) = CS%Ref_val(m)%p(i,j,1:nz_data) - tmp_val1(1:nz_data) + enddo; enddo + deallocate(tmp_val2) + enddo + + ! remap u to h_obs to get increment + allocate(hu_obs(nz_data)) ; hu_obs(:)=0.0 ; hu(:) = 0.0; + allocate(tmp_val2(nz));tmp_val2(:)=0.0 + do j=js,je ; do i=isB,ieB + tmp_val2(1:nz) = CS%var_u%p(i ,j,1:nz) + hu(1:nz) = 0.5*( h(i,j,1:nz )+ h(i+1,j,1:nz )) + hu_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i+1,j,1:nz_data)) + call remapping_core_h(CS%remap_cs, nz , hu(1:nz ), tmp_val2, & + nz_data, hu_obs(1:nz_data), tmp_val1, & + h_neglect, h_neglect_edge) + CS%Ref_oda_u%p(i,j,1:nz_data) = CS%Ref_oda_u%p(i,j,1:nz_data) - tmp_val1(1:nz_data) + enddo; enddo + deallocate(tmp_val2,hu_obs) + + ! remap v to h_obs to get increment + allocate(hv_obs(nz_data)) ; hv_obs(:)=0.0 ; hv(:) = 0.0; + allocate(tmp_val2(nz));tmp_val2(:)=0.0 + do j=jsB,jeB ; do i=is,ie + tmp_val2(1:nz) = CS%var_v%p(i ,j,1:nz) + hv(1:nz) = 0.5*( h(i,j,1:nz )+ h(i,j+1,1:nz )) + hv_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i,j+1,1:nz_data)) + call remapping_core_h(CS%remap_cs, nz , hv(1:nz ), tmp_val2, & + nz_data, hv_obs(1:nz_data), tmp_val1, & + h_neglect, h_neglect_edge) + CS%Ref_oda_v%p(i,j,1:nz_data) = CS%Ref_oda_v%p(i,j,1:nz_data) - tmp_val1(1:nz_data) + enddo; enddo + deallocate(tmp_val2,hv_obs) + +end subroutine !> This subroutine applies oda increments to layers thicknesses, temp, salt, U !! and V everywhere . @@ -357,18 +435,14 @@ subroutine apply_oda_incupd(h, dt, G, GV, US, CS) real, dimension(SZK_(GV)) :: h_col ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] real, dimension(SZK_(GV)) :: hu, hv ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] - !real, allocatable, dimension(:,:,:) :: tmp_p ! A temporary array for int. depth field - !real, allocatable, dimension(:,:,:) :: tmp_h ! A temporary array for thickness field - real, allocatable, dimension(:,:,:) :: tmp_u !< A temporary array for u inc. - real, allocatable, dimension(:,:,:) :: tmp_v !< A temporary array for v inc. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tmp_t !< A temporary array for t inc. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tmp_s !< A temporary array for s inc. + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: tmp_u !< A temporary array for u inc. + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: tmp_v !< A temporary array for v inc. -! real, allocatable, dimension(:,:,:,:) :: var_inc !< increment fields if input is full fields -! real, allocatable, dimension(:,:,:) :: varu_inc !< increment U field if input is full fields -! real, allocatable, dimension(:,:,:) :: varv_inc !< increment V field if input is full fields real, allocatable, dimension(:,:,:) :: h_obs !< h of increments real, allocatable, dimension(:) :: hu_obs,hv_obs ! A column of thicknesses at h points [H ~> m or kg m-2] - integer :: m, i, j, k, is, ie, js, je, nz, nz_data integer :: isB, ieB, jsB, jeB integer :: ncount ! time step counter @@ -381,34 +455,14 @@ subroutine apply_oda_incupd(h, dt, G, GV, US, CS) if (.not.associated(CS)) return ! update counter - if (CS%nstep_incupd /= 1) then - CS%ncount = CS%ncount+1 - if (CS%ncount == 1 .or. CS%ncount == CS%nstep_incupd+1) then - q = 0.5/CS%nstep_incupd !half increment on 1st and CS%step_incupd+1th step - else - q = 1.0/CS%nstep_incupd - endif !ncount - else - CS%ncount = CS%ncount+1 - q = 1.0/CS%nstep_incupd - endif - - ! Diagnostics of increments, mostly used for debugging. - if (CS%id_u_oda_inc > 0) call post_data(CS%id_u_oda_inc, CS%Ref_oda_u%p, CS%diag) - if (CS%id_v_oda_inc > 0) call post_data(CS%id_v_oda_inc, CS%Ref_oda_v%p, CS%diag) - if (CS%id_h_oda_inc > 0) call post_data(CS%id_h_oda_inc, CS%Ref_h%p, CS%diag) - if (CS%id_T_oda_inc > 0) call post_data(CS%id_T_oda_inc, CS%Ref_val(1)%p,CS%diag) - if (CS%id_S_oda_inc > 0) call post_data(CS%id_S_oda_inc, CS%Ref_val(2)%p,CS%diag) - - + CS%ncount = CS%ncount+1 + q = 1.0/CS%nstep_incupd - ! no assimilation after CS%step_incupd+1 - if (CS%ncount > CS%nstep_incupd+1) then - if (CS%ncount == CS%nstep_incupd+2) then - if (is_root_pe()) call MOM_error(NOTE,"ended updating fields with increments. ") - endif !ncount==CS%nstep_incupd+2 + ! no assimilation after CS%step_incupd + if (CS%ncount > CS%nstep_incupd) then + if (is_root_pe()) call MOM_error(NOTE,"ended updating fields with increments. ") return - endif !ncount>CS%nstep_incupd+1 + endif !ncount>CS%nstep_incupd ! get h_obs nz_data = CS%Ref_val(1)%nz_data @@ -419,60 +473,6 @@ subroutine apply_oda_incupd(h, dt, G, GV, US, CS) call pass_var(h_obs,G%Domain) - if (CS%oda_inc) then ! input are increments - if (is_root_pe()) call MOM_error(NOTE,"incupd using increments fields ") - - else ! inputs are full fields - if (q == 1.0 .or. CS%ncount == 1) then - if (is_root_pe()) call MOM_error(NOTE,"incupd using full fields ") - - ! remap t,s (on h_init) to h_obs to get increment - tmp_val1(:)=0.0 - do m=1,CS%fldno - allocate(tmp_val2(nz)) - do j=js,je ; do i=is,ie - tmp_val2(1:nz) = CS%var(m)%p(i,j,1:nz) - call remapping_core_h(CS%remap_cs, nz , h(i,j,1:nz ), tmp_val2, & - nz_data, h_obs(i,j,1:nz_data), tmp_val1, & - h_neglect, h_neglect_edge) - CS%Ref_val(m)%p(i,j,1:nz_data) = CS%Ref_val(m)%p(i,j,1:nz_data) - tmp_val1(1:nz_data) - enddo; enddo - deallocate(tmp_val2) - enddo - - ! remap u to h_obs to get increment - allocate(hu_obs(nz_data)) ; hu_obs(:)=0.0 ; hu(:) = 0.0; - allocate(tmp_val2(nz));tmp_val2(:)=0.0 - do j=js,je ; do i=isB,ieB - tmp_val2(1:nz) = CS%var_u%p(i ,j,1:nz) - hu(1:nz) = 0.5*( h(i,j,1:nz )+ h(i+1,j,1:nz )) - hu_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i+1,j,1:nz_data)) - call remapping_core_h(CS%remap_cs, nz , hu(1:nz ), tmp_val2, & - nz_data, hu_obs(1:nz_data), tmp_val1, & - h_neglect, h_neglect_edge) - CS%Ref_oda_u%p(i,j,1:nz_data) = CS%Ref_oda_u%p(i,j,1:nz_data) - tmp_val1(1:nz_data) - enddo; enddo - deallocate(tmp_val2,hu_obs) - - ! remap v to h_obs to get increment - allocate(hv_obs(nz_data)) ; hv_obs(:)=0.0 ; hv(:) = 0.0; - allocate(tmp_val2(nz));tmp_val2(:)=0.0 - do j=jsB,jeB ; do i=is,ie - tmp_val2(1:nz) = CS%var_v%p(i ,j,1:nz) - hv(1:nz) = 0.5*( h(i,j,1:nz )+ h(i,j+1,1:nz )) - hv_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i,j+1,1:nz_data)) - call remapping_core_h(CS%remap_cs, nz , hv(1:nz ), tmp_val2, & - nz_data, hv_obs(1:nz_data), tmp_val1, & - h_neglect, h_neglect_edge) - CS%Ref_oda_v%p(i,j,1:nz_data) = CS%Ref_oda_v%p(i,j,1:nz_data) - tmp_val1(1:nz_data) - enddo; enddo - deallocate(tmp_val2,hv_obs) - - endif ! not q == 1 - - endif ! not oda_inc - - ! print out increments write(mesg,'(i8)') CS%ncount if (is_root_pe()) call MOM_error(NOTE,"updating fields with increments ncount:"//trim(mesg)) @@ -492,6 +492,7 @@ subroutine apply_oda_incupd(h, dt, G, GV, US, CS) ! add increments to tracers tmp_val1(:)=0.0 + tmp_t(:,:,:) = 0.0 ; tmp_s(:,:,:) = 0.0 ! diagnostics do m=1,CS%fldno allocate(tmp_val2(CS%Ref_val(m)%nz_data)) do j=js,je ; do i=is,ie @@ -500,6 +501,11 @@ subroutine apply_oda_incupd(h, dt, G, GV, US, CS) nz , h(i,j,1:nz ), tmp_val1, & h_neglect, h_neglect_edge) CS%var(m)%p(i,j,1:nz) = CS%var(m)%p(i,j,1:nz) + q*tmp_val1(1:nz) + if (m == 1) then + tmp_t(i,j,1:nz) = tmp_val1(1:nz) + else + tmp_s(i,j,1:nz) = tmp_val1(1:nz) + endif enddo; enddo deallocate(tmp_val2) enddo @@ -511,9 +517,10 @@ subroutine apply_oda_incupd(h, dt, G, GV, US, CS) enddo enddo ; enddo -! add increments to u + ! add increments to u allocate(hu_obs(nz_data)) ; hu_obs(:)=0.0 ; hu(:) = 0.0; allocate(tmp_val2(nz_data));tmp_val2(:)=0.0 + tmp_u(:,:,:) = 0.0 ! diagnostics do j=js,je ; do i=isB,ieB tmp_val2(1:nz_data) = CS%Ref_oda_u%p(i,j,1:nz_data) hu_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i+1,j,1:nz_data)) @@ -522,12 +529,14 @@ subroutine apply_oda_incupd(h, dt, G, GV, US, CS) nz , hu(1:nz ), tmp_val1, & h_neglect, h_neglect_edge) CS%var_u%p(i,j,1:nz) = CS%var_u%p(i,j,1:nz) + q*tmp_val1(1:nz) + tmp_u(i,j,1:nz) = tmp_val1(1:nz) enddo; enddo deallocate(tmp_val2,hu_obs) -! add increments to v + ! add increments to v allocate(hv_obs(nz_data)) ; hv_obs(:)=0.0 ; hv(:) = 0.0; allocate(tmp_val2(nz_data));tmp_val2(:)=0.0 + tmp_v(:,:,:) = 0.0 ! diagnostics do j=jsB,jeB ; do i=is,ie tmp_val2(1:nz_data) = CS%Ref_oda_v%p(i,j,1:nz_data) hv_obs(1:nz_data) = 0.5 * (h_obs(i,j,1:nz_data)+h_obs(i,j+1,1:nz_data)) @@ -536,9 +545,17 @@ subroutine apply_oda_incupd(h, dt, G, GV, US, CS) nz , hv(1:nz ), tmp_val1, & h_neglect, h_neglect_edge) CS%var_v%p(i,j,1:nz) = CS%var_v%p(i,j,1:nz) + q*tmp_val1(1:nz) + tmp_v(i,j,1:nz) = tmp_val1(1:nz) enddo; enddo deallocate(tmp_val2,hv_obs) + ! Diagnostics of increments, mostly used for debugging. + if (CS%id_u_oda_inc > 0) call post_data(CS%id_u_oda_inc, tmp_u, CS%diag) + if (CS%id_v_oda_inc > 0) call post_data(CS%id_v_oda_inc, tmp_v, CS%diag) + if (CS%id_h_oda_inc > 0) call post_data(CS%id_h_oda_inc, h , CS%diag) + if (CS%id_T_oda_inc > 0) call post_data(CS%id_T_oda_inc, tmp_t, CS%diag) + if (CS%id_S_oda_inc > 0) call post_data(CS%id_S_oda_inc, tmp_s, CS%diag) + end subroutine apply_oda_incupd From 3feaf5f991a9ca3d0b50e6d6c6df73468360ad39 Mon Sep 17 00:00:00 2001 From: abozec Date: Wed, 26 May 2021 11:24:17 -0400 Subject: [PATCH 4/9] Add uv inc option, mv get_oda_increment at initialization,add SSH correction and general clean up - Add option to use or not uv increment - The calculation of the increment when using full field has been moved at initialization - A SSH correction has been applied to h_obs to account for the different SSH between the background and target state - The incremental update is now done over nstep-1 - a general clean up has been done with variable name change to make the code easier to follow --- .../MOM_state_initialization.F90 | 29 +- src/ocean_data_assim/MOM_oda_incupd.F90 | 397 ++++++++++-------- 2 files changed, 235 insertions(+), 191 deletions(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index a46304cfe4..48a0a02ef1 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -96,7 +96,7 @@ module MOM_state_initialization use MOM_horizontal_regridding, only : horiz_interp_and_extrap_tracer use MOM_oda_incupd, only: oda_incupd_CS, initialize_oda_incupd use MOM_oda_incupd, only: set_up_oda_incupd_field, set_up_oda_incupd_vel_field -use MOM_oda_incupd, only: get_oda_increments +use MOM_oda_incupd, only: get_oda_increments, apply_oda_incupd implicit none ; private @@ -2067,6 +2067,7 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, p integer, dimension(4) :: siz integer :: nz_data ! The size of the sponge source grid logical :: oda_inc ! input files are increments (true) or full fields (false) + logical :: uv_inc ! use u and v increments character(len=40) :: tempinc_var, salinc_var, uinc_var, vinc_var, h_var character(len=40) :: mdl = "initialize_oda_incupd_file" @@ -2096,6 +2097,9 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, p call get_param(param_file, mdl, "ODA_THK_VAR", h_var, & "The name of the layer thickness variable in "//& "ODA_INCUPD_FILE.", default="h") + call get_param(param_file, mdl, "ODA_INCUPD_UV", uv_inc, & + "use U,V increments.", & + default=.true.) call get_param(param_file, mdl, "ODA_INCUPD_UV_FILE", uv_inc_file, & "The name of the file with the U,V increments.", & default=inc_file) @@ -2135,16 +2139,18 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, p endif ! get U and V increments - filename = trim(inputdir)//trim(uv_inc_file) - call log_param(param_file, mdl, "INPUTDIR/ODA_INCUPD_UV_FILE", filename) - if (.not.file_exists(filename, G%Domain)) & - call MOM_error(FATAL, " initialize_oda_incupd_uv: Unable to open "//trim(filename)) - allocate(tmp_u(G%IsdB:G%IedB,jsd:jed,nz_data)) - allocate(tmp_v(isd:ied,G%JsdB:G%JedB,nz_data)) - tmp_u(:,:,:) = 0.0 ; tmp_v(:,:,:) = 0.0 - call MOM_read_vector(filename, uinc_var, vinc_var, tmp_u, tmp_v, G%Domain,scale=US%m_s_to_L_T) - call set_up_oda_incupd_vel_field(tmp_u, tmp_v, G, GV, u, v, oda_incupd_CSp) - deallocate(tmp_u,tmp_v) + if (uv_inc) then + filename = trim(inputdir)//trim(uv_inc_file) + call log_param(param_file, mdl, "INPUTDIR/ODA_INCUPD_UV_FILE", filename) + if (.not.file_exists(filename, G%Domain)) & + call MOM_error(FATAL, " initialize_oda_incupd_uv: Unable to open "//trim(filename)) + allocate(tmp_u(G%IsdB:G%IedB,jsd:jed,nz_data)) + allocate(tmp_v(isd:ied,G%JsdB:G%JedB,nz_data)) + tmp_u(:,:,:) = 0.0 ; tmp_v(:,:,:) = 0.0 + call MOM_read_vector(filename, uinc_var, vinc_var, tmp_u, tmp_v, G%Domain,scale=US%m_s_to_L_T) + call set_up_oda_incupd_vel_field(tmp_u, tmp_v, G, GV, u, v, oda_incupd_CSp) + deallocate(tmp_u,tmp_v) + endif ! calculate increments if full fields if (oda_inc) then ! input are increments @@ -2152,6 +2158,7 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, p else ! inputs are full fields if (is_root_pe()) call MOM_error(NOTE,"incupd using full fields ") call get_oda_increments(h, G, GV, US, oda_incupd_CSp) + ! call apply_oda_incupd(h, 0.0, G, GV, US, oda_incupd_CSp) endif ! not oda_inc end subroutine initialize_oda_incupd_file diff --git a/src/ocean_data_assim/MOM_oda_incupd.F90 b/src/ocean_data_assim/MOM_oda_incupd.F90 index 6549b95bee..e86c8be2d6 100644 --- a/src/ocean_data_assim/MOM_oda_incupd.F90 +++ b/src/ocean_data_assim/MOM_oda_incupd.F90 @@ -54,15 +54,6 @@ module MOM_oda_incupd real, dimension(:,:,:), pointer :: h => NULL() !< pointer to the data grid. end type p3d -!> A structure for creating arrays of pointers to 2D arrays with extra gridding information -type :: p2d - integer :: id !< id for FMS external time interpolator - integer :: nz_data !< The number of vertical levels in the input field - real, dimension(:,:), pointer :: mask_in => NULL()!< pointer to the data mask. - real, dimension(:,:), pointer :: p => NULL() !< pointer the data. - real, dimension(:,:), pointer :: h => NULL() !< pointer the data grid. -end type p2d - !> oda incupd control structure type, public :: oda_incupd_CS ; private integer :: nz !< The total number of layers. @@ -86,16 +77,13 @@ module MOM_oda_incupd integer :: fldno = 0 !< The number of fields which have already been !! registered by calls to set_up_oda_incupd_field - type(p3d) :: var(MAX_FIELDS_) !< Pointers to the fields that are being updated. - type(p3d) :: Ref_val(MAX_FIELDS_) !< The values to which the fields are updated. - type(p3d) :: Ref_oda_u !< The values to which the u-velocities are updated. - type(p3d) :: Ref_oda_v !< The values to which the v-velocities are updated. - type(p3d) :: var_u !< Pointer to the u velocities. that are being updated. - type(p3d) :: var_v !< Pointer to the v velocities. that are being updated. - type(p3d) :: Ref_h !< Grid on which reference data is provided (older code). - type(p3d) :: Ref_p !< grid of layer pressure increment - type(p3d) :: Ref_hu !< u-point grid on which reference data is provided (older code). - type(p3d) :: Ref_hv !< v-point grid on which reference data is provided (older code). + type(p3d) :: Var(MAX_FIELDS_) !< Pointers to the fields that are being updated. + type(p3d) :: Inc(MAX_FIELDS_) !< The increments to be applied to the field + type(p3d) :: Var_u !< Pointer to the u-velocities. that are being updated. + type(p3d) :: Inc_u !< The increments to be applied to the u-velocities + type(p3d) :: Var_v !< Pointer to the v-velocities. that are being updated. + type(p3d) :: Inc_v !< The increments to be applied to the v-velocities + type(p3d) :: Ref_h !< Vertical grid on which the increments are provided integer :: nstep_incupd !< number of time step for full update @@ -109,9 +97,7 @@ module MOM_oda_incupd !! rotationally symmetric forms of the same expressions. logical :: incupdDataOngrid !< True if the incupd data are on the model horizontal grid - - logical :: reentrant_x !< grid is reentrant in the x direction - logical :: tripolar_N !< grid is folded at its north edge + logical :: uv_inc !< use u and v increments ! for diagnostics type(diag_ctrl), pointer :: diag ! m or kg m-2]. - ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mdl = "MOM_oda" ! This module's name. logical :: use_oda_incupd - logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries + logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries logical :: default_2018_answers integer :: i, j, k - real :: nhours_incupd, dt_baclin + real :: nhours_incupd, dt, dt_therm character(len=256) :: mesg character(len=10) :: remapScheme if (associated(CS)) then @@ -171,16 +156,23 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data) call get_param(param_file, mdl, "ODA_INCUPD_NHOURS", nhours_incupd, & "Number of hours for full update (0=direct insertion).", & default=3.0) - - call get_param(param_file, mdl, "DT_THERM", dt_baclin, & - "The thermodynamic and tracer advection time step."// & - "Ideally DT_THERM should be an integer multiple of DT"//& - "and less than the forcing or coupling time-step, unless"//& - "THERMO_SPANS_COUPLING is true, in which case DT_THERM"//& - "can be an integer multiple of the coupling timestep. By"//& - "default DT_THERM is set to DT.",& - default=300.0) - + call get_param(param_file, mdl, "DT", dt, & + "The (baroclinic) dynamics time step. The time-step that "//& + "is actually used will be an integer fraction of the "//& + "forcing time-step (DT_FORCING in ocean-only mode or the "//& + "coupling timestep in coupled mode.)", units="s", scale=US%s_to_T, & + fail_if_missing=.true.) + call get_param(param_file, mdl, "DT_THERM", dt_therm, & + "The thermodynamic and tracer advection time step. "//& + "Ideally DT_THERM should be an integer multiple of DT "//& + "and less than the forcing or coupling time-step, unless "//& + "THERMO_SPANS_COUPLING is true, in which case DT_THERM "//& + "can be an integer multiple of the coupling timestep. By "//& + "default DT_THERM is set to DT.", & + units="s", scale=US%s_to_T, default=US%T_to_s*dt) + call get_param(param_file, mdl, "ODA_INCUPD_UV", CS%uv_inc, & + "use U,V increments.", & + default=.true.) call get_param(param_file, mdl, "REMAPPING_SCHEME", remapScheme, & "This sets the reconstruction scheme used "//& " for vertical remapping for all variables.", & @@ -203,18 +195,13 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data) "When defined, the incoming oda_incupd data are "//& "assumed to be on the model horizontal grid " , & default=.true.) - call get_param(param_file, mdl, "REENTRANT_X", CS%reentrant_x, & - "If true, the domain is zonally reentrant.", default=.true.) - call get_param(param_file, mdl, "TRIPOLAR_N", CS%tripolar_N, & - "Use tripolar connectivity at the northern edge of the "//& - "domain. With TRIPOLAR_N, NIGLOBAL must be even.", default=.false.) CS%nz = GV%ke CS%isc = G%isc ; CS%iec = G%iec ; CS%jsc = G%jsc ; CS%jec = G%jec CS%isd = G%isd ; CS%ied = G%ied ; CS%jsd = G%jsd ; CS%jed = G%jed CS%isdB = G%isdB ; CS%iedB = G%iedB; CS%jsdB = G%jsdB ; CS%jedB = G%jedB - !! increments on horizontal grid + ! increments on horizontal grid if (.not. CS%incupdDataOngrid) call MOM_error(FATAL,'initialize_oda_incupd: '// & 'The oda_incupd code only applies ODA increments on the same horizontal grid. ') @@ -222,7 +209,7 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data) if (nhours_incupd == 0) then CS%nstep_incupd = 1 !! direct insertion else - CS%nstep_incupd = nhours_incupd*3600./dt_baclin + CS%nstep_incupd = floor(nhours_incupd*3600./dt_therm + 0.001 ) - 1 endif write(mesg,'(i12)') CS%nstep_incupd if (is_root_pe()) & @@ -232,7 +219,7 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data) ! initialize time counter CS%ncount=0 - ! get layer pressure increment p + ! get the vertical grid (h_obs) of the increments CS%nz_data = nz_data allocate(CS%Ref_h%p(G%isd:G%ied,G%jsd:G%jed,CS%nz_data)) CS%Ref_h%p(:,:,:) = 0.0 ; @@ -251,14 +238,14 @@ end subroutine initialize_oda_incupd !> This subroutine stores theincrements at h points for the variable !! whose address is given by f_ptr. subroutine set_up_oda_incupd_field(sp_val, G, GV, f_ptr, CS) - type(ocean_grid_type), intent(in) :: G !< Grid structure - type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure - type(oda_incupd_CS), pointer :: CS !< oda_incupd control structure (in/out). + type(ocean_grid_type), intent(in) :: G !< Grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(oda_incupd_CS), pointer :: CS !< oda_incupd control structure (in/out). real, dimension(SZI_(G),SZJ_(G),CS%nz_data), & intent(in) :: sp_val !< increment field, it can have an - !! arbitrary number of layers. + !! arbitrary number of layers. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - target, intent(in) :: f_ptr !< Pointer to the field to be updated + target, intent(in) :: f_ptr !< Pointer to the field to be updated integer :: i, j, k character(len=256) :: mesg ! String for error messages @@ -273,17 +260,16 @@ subroutine set_up_oda_incupd_field(sp_val, G, GV, f_ptr, CS) call MOM_error(FATAL,"set_up_oda_incupd_field: "//mesg) endif - ! stores the reference profile - CS%Ref_val(CS%fldno)%nz_data = CS%nz_data - allocate(CS%Ref_val(CS%fldno)%p(G%isd:G%ied,G%jsd:G%jed,CS%nz_data)) - CS%Ref_val(CS%fldno)%p(:,:,:) = 0.0 + ! store the increment/full field tracer profiles + CS%Inc(CS%fldno)%nz_data = CS%nz_data + allocate(CS%Inc(CS%fldno)%p(G%isd:G%ied,G%jsd:G%jed,CS%nz_data)) + CS%Inc(CS%fldno)%p(:,:,:) = 0.0 do j=G%jsc,G%jec ; do i=G%isc,G%iec do k=1,CS%nz_data - CS%Ref_val(CS%fldno)%p(i,j,k) = sp_val(i,j,k) + CS%Inc(CS%fldno)%p(i,j,k) = sp_val(i,j,k) enddo enddo ; enddo - - CS%var(CS%fldno)%p => f_ptr + CS%Var(CS%fldno)%p => f_ptr ! pointer to model tracers end subroutine set_up_oda_incupd_field @@ -309,24 +295,25 @@ subroutine set_up_oda_incupd_vel_field(u_val, v_val, G, GV, u_ptr, v_ptr, CS) if (.not.associated(CS)) return - ! stores the reference profile - allocate(CS%Ref_oda_u%p(G%isdB:G%iedB,G%jsd:G%jed,CS%nz_data)) - CS%Ref_oda_u%p(:,:,:) = 0.0 + ! store the increment/full field u profile + allocate(CS%Inc_u%p(G%isdB:G%iedB,G%jsd:G%jed,CS%nz_data)) + CS%Inc_u%p(:,:,:) = 0.0 do j=G%jsc,G%jec ; do i=G%iscB,G%iecB do k=1,CS%nz_data - CS%Ref_oda_u%p(i,j,k) = u_val(i,j,k) + CS%Inc_u%p(i,j,k) = u_val(i,j,k) enddo enddo ; enddo - CS%var_u%p => u_ptr + CS%Var_u%p => u_ptr ! pointer to model u-velocity - allocate(CS%Ref_oda_v%p(G%isd:G%ied,G%jsdB:G%jedB,CS%nz_data)) - CS%Ref_oda_v%p(:,:,:) = 0.0 + ! store the increment/full field v profile + allocate(CS%Inc_v%p(G%isd:G%ied,G%jsdB:G%jedB,CS%nz_data)) + CS%Inc_v%p(:,:,:) = 0.0 do j=G%jscB,G%jecB ; do i=G%isc,G%iec do k=1,CS%nz_data - CS%Ref_oda_v%p(i,j,k) = v_val(i,j,k) + CS%Inc_v%p(i,j,k) = v_val(i,j,k) enddo enddo ; enddo - CS%var_v%p => v_ptr + CS%Var_v%p => v_ptr ! pointer to model v-velocity end subroutine set_up_oda_incupd_vel_field @@ -341,11 +328,12 @@ subroutine get_oda_increments(h, G, GV, US, CS) !! that is set by a previous call to initialize_oda_incupd (in). - real, allocatable, dimension(:) :: tmp_val2 ! data values on the original grid - real, dimension(SZK_(GV)) :: tmp_val1 ! data values remapped to model grid - real, allocatable, dimension(:,:,:) :: h_obs !< h of increments + real, dimension(SZK_(GV)) :: tmp_val1 ! data values on the model grid + real, allocatable, dimension(:) :: tmp_val2 ! data values remapped to increment grid + real, allocatable, dimension(:,:,:) :: h_obs !< h of increments + real, allocatable, dimension(:) :: tmp_h ! temporary array for corrected h_obs real, allocatable, dimension(:) :: hu_obs,hv_obs ! A column of thicknesses at h points [H ~> m or kg m-2] - real, dimension(SZK_(GV)) :: hu, hv ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] + real, dimension(SZK_(GV)) :: hu, hv ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] integer :: m, i, j, k, is, ie, js, je, nz, nz_data @@ -366,81 +354,104 @@ subroutine get_oda_increments(h, G, GV, US, CS) endif ! get h_obs - nz_data = CS%Ref_val(1)%nz_data + nz_data = CS%Inc(1)%nz_data allocate(h_obs(G%isd:G%ied,G%jsd:G%jed,nz_data)) ; h_obs(:,:,:)=0.0 do k=1,nz_data ; do j=js,je ; do i=is,ie h_obs(i,j,k) = CS%Ref_h%p(i,j,k) enddo ; enddo ; enddo call pass_var(h_obs,G%Domain) + ! allocate 1-d arrays + allocate(tmp_h(nz_data)); tmp_h(:)=0.0 + allocate(tmp_val2(nz_data)) ; tmp_val2(:) = 0.0 + allocate(hu_obs(nz_data)) ; hu_obs(:)=0.0 + allocate(hv_obs(nz_data)) ; hv_obs(:)=0.0 + ! remap t,s (on h_init) to h_obs to get increment tmp_val1(:)=0.0 - do m=1,CS%fldno - allocate(tmp_val2(nz)) - do j=js,je ; do i=is,ie - tmp_val2(1:nz) = CS%var(m)%p(i,j,1:nz) - call remapping_core_h(CS%remap_cs, nz , h(i,j,1:nz ), tmp_val2, & - nz_data, h_obs(i,j,1:nz_data), tmp_val1, & + do j=js,je ; do i=is,ie + ! account for the different SSH + tmp_h(1:nz_data)=(sum(h(i,j,1:nz))/sum(h_obs(i,j,1:nz_data)))*h_obs(i,j,1:nz_data) + do m=1,CS%fldno + ! get tracer + tmp_val1(1:nz) = CS%Var(m)%p(i,j,1:nz) + ! remap tracer on h_obs + call remapping_core_h(CS%remap_cs, nz , h(i,j,1:nz ), tmp_val1, & + nz_data, tmp_h( 1:nz_data), tmp_val2, & h_neglect, h_neglect_edge) - CS%Ref_val(m)%p(i,j,1:nz_data) = CS%Ref_val(m)%p(i,j,1:nz_data) - tmp_val1(1:nz_data) - enddo; enddo - deallocate(tmp_val2) - enddo + ! get increment from full field on h_obs + CS%Inc(m)%p(i,j,1:nz_data) = CS%Inc(m)%p(i,j,1:nz_data) - tmp_val2(1:nz_data) + enddo + enddo; enddo ! remap u to h_obs to get increment - allocate(hu_obs(nz_data)) ; hu_obs(:)=0.0 ; hu(:) = 0.0; - allocate(tmp_val2(nz));tmp_val2(:)=0.0 - do j=js,je ; do i=isB,ieB - tmp_val2(1:nz) = CS%var_u%p(i ,j,1:nz) - hu(1:nz) = 0.5*( h(i,j,1:nz )+ h(i+1,j,1:nz )) - hu_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i+1,j,1:nz_data)) - call remapping_core_h(CS%remap_cs, nz , hu(1:nz ), tmp_val2, & - nz_data, hu_obs(1:nz_data), tmp_val1, & - h_neglect, h_neglect_edge) - CS%Ref_oda_u%p(i,j,1:nz_data) = CS%Ref_oda_u%p(i,j,1:nz_data) - tmp_val1(1:nz_data) - enddo; enddo - deallocate(tmp_val2,hu_obs) - - ! remap v to h_obs to get increment - allocate(hv_obs(nz_data)) ; hv_obs(:)=0.0 ; hv(:) = 0.0; - allocate(tmp_val2(nz));tmp_val2(:)=0.0 - do j=jsB,jeB ; do i=is,ie - tmp_val2(1:nz) = CS%var_v%p(i ,j,1:nz) - hv(1:nz) = 0.5*( h(i,j,1:nz )+ h(i,j+1,1:nz )) - hv_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i,j+1,1:nz_data)) - call remapping_core_h(CS%remap_cs, nz , hv(1:nz ), tmp_val2, & - nz_data, hv_obs(1:nz_data), tmp_val1, & - h_neglect, h_neglect_edge) - CS%Ref_oda_v%p(i,j,1:nz_data) = CS%Ref_oda_v%p(i,j,1:nz_data) - tmp_val1(1:nz_data) - enddo; enddo - deallocate(tmp_val2,hv_obs) + if (CS%uv_inc) then + hu(:)=0.0 + do j=js,je ; do i=isB,ieB + ! get u-velocity + tmp_val1(1:nz) = CS%Var_u%p(i ,j,1:nz) + ! get the h and h_obs at u points + hu(1:nz) = 0.5*( h(i,j,1:nz )+ h(i+1,j,1:nz )) + hu_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i+1,j,1:nz_data)) + ! account for the different SSH + hu_obs(1:nz_data) = (sum(hu(1:nz))/sum(hu_obs(1:nz_data))) * hu_obs(1:nz_data) + ! remap model u on hu_obs + call remapping_core_h(CS%remap_cs, nz , hu(1:nz ), tmp_val1, & + nz_data, hu_obs(1:nz_data), tmp_val2, & + h_neglect, h_neglect_edge) + ! get increment from full field on h_obs + CS%Inc_u%p(i,j,1:nz_data) = CS%Inc_u%p(i,j,1:nz_data) - tmp_val2(1:nz_data) + enddo; enddo + + ! remap v to h_obs to get increment + hv(:) = 0.0; + do j=jsB,jeB ; do i=is,ie + ! get v-velocity + tmp_val1(1:nz) = CS%Var_v%p(i ,j,1:nz) + ! get the h and h_obs at v points + hv(1:nz) = 0.5*( h(i,j,1:nz )+ h(i,j+1,1:nz )) + hv_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i,j+1,1:nz_data)) + ! account for the different SSH + hv_obs(1:nz_data) = (sum(hv(1:nz))/sum(hv_obs(1:nz_data))) *hv_obs(1:nz_data) + ! remap model v on hv_obs + call remapping_core_h(CS%remap_cs, nz , hv(1:nz ), tmp_val1, & + nz_data, hv_obs(1:nz_data), tmp_val2, & + h_neglect, h_neglect_edge) + ! get increment from full field on h_obs + CS%Inc_v%p(i,j,1:nz_data) = CS%Inc_v%p(i,j,1:nz_data) - tmp_val2(1:nz_data) + enddo; enddo + endif ! uv_inc + + ! deallocate arrays + deallocate(tmp_h,tmp_val2,hu_obs,hv_obs) + deallocate(h_obs) end subroutine !> This subroutine applies oda increments to layers thicknesses, temp, salt, U !! and V everywhere . subroutine apply_oda_incupd(h, dt, G, GV, US, CS) - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure (in). + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure (in). type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] (in) + intent(in) :: h !< Layer thickness [H ~> m or kg m-2] (in) real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]. type(oda_incupd_CS), pointer :: CS !< A pointer to the control structure for this module !! that is set by a previous call to initialize_oda_incupd (in). real :: m_to_Z ! A unit conversion factor from m to Z. - real, allocatable, dimension(:) :: tmp_val2 ! data values on the original grid + real, allocatable, dimension(:) :: tmp_val2 ! data values on the increment grid real, dimension(SZK_(GV)) :: tmp_val1 ! data values remapped to model grid - real, dimension(SZK_(GV)) :: h_col ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] real, dimension(SZK_(GV)) :: hu, hv ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tmp_t !< A temporary array for t inc. - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tmp_s !< A temporary array for s inc. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tmp_t !< A temporary array for t inc. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tmp_s !< A temporary array for s inc. real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: tmp_u !< A temporary array for u inc. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: tmp_v !< A temporary array for v inc. real, allocatable, dimension(:,:,:) :: h_obs !< h of increments + real, allocatable, dimension(:) :: tmp_h !< temporary array for corrected h_obs real, allocatable, dimension(:) :: hu_obs,hv_obs ! A column of thicknesses at h points [H ~> m or kg m-2] integer :: m, i, j, k, is, ie, js, je, nz, nz_data @@ -464,14 +475,6 @@ subroutine apply_oda_incupd(h, dt, G, GV, US, CS) return endif !ncount>CS%nstep_incupd - ! get h_obs - nz_data = CS%Ref_val(1)%nz_data - allocate(h_obs(G%isd:G%ied,G%jsd:G%jed,nz_data)) ; h_obs(:,:,:)=0.0 - do k=1,nz_data ; do j=js,je ; do i=is,ie - h_obs(i,j,k) = CS%Ref_h%p(i,j,k) - enddo ; enddo ; enddo - call pass_var(h_obs,G%Domain) - ! print out increments write(mesg,'(i8)') CS%ncount @@ -479,8 +482,6 @@ subroutine apply_oda_incupd(h, dt, G, GV, US, CS) write(mesg,'(f10.8)') q if (is_root_pe()) call MOM_error(NOTE,"updating fields with weight q:"//trim(mesg)) - - if (.not.CS%remap_answers_2018) then h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff elseif (GV%Boussinesq) then @@ -488,76 +489,112 @@ subroutine apply_oda_incupd(h, dt, G, GV, US, CS) else h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 endif - + + ! get h_obs + nz_data = CS%Inc(1)%nz_data + allocate(h_obs(G%isd:G%ied,G%jsd:G%jed,nz_data)) ; h_obs(:,:,:)=0.0 + do k=1,nz_data ; do j=js,je ; do i=is,ie + h_obs(i,j,k) = CS%Ref_h%p(i,j,k) + enddo ; enddo ; enddo + call pass_var(h_obs,G%Domain) + + ! allocate 1-d array + allocate(tmp_h(nz_data)); tmp_h(:)=0.0 + allocate(tmp_val2(nz_data)) + allocate(hu_obs(nz_data)) ; hu_obs(:)=0.0 + allocate(hv_obs(nz_data)) ; hv_obs(:)=0.0 ! add increments to tracers - tmp_val1(:)=0.0 + tmp_val1(:)=0.0 tmp_t(:,:,:) = 0.0 ; tmp_s(:,:,:) = 0.0 ! diagnostics - do m=1,CS%fldno - allocate(tmp_val2(CS%Ref_val(m)%nz_data)) - do j=js,je ; do i=is,ie - tmp_val2(1:nz_data) = CS%Ref_val(m)%p(i,j,1:nz_data) - call remapping_core_h(CS%remap_cs, nz_data, h_obs(i,j,1:nz_data), tmp_val2, & - nz , h(i,j,1:nz ), tmp_val1, & + do j=js,je ; do i=is,ie + ! account for the different SSH + tmp_h(1:nz_data)=(sum(h(i,j,1:nz))/sum(h_obs(i,j,1:nz_data)))*h_obs(i,j,1:nz_data) + do m=1,CS%fldno + ! get increment + tmp_val2(1:nz_data) = CS%Inc(m)%p(i,j,1:nz_data) + ! remap increment profile on model h + call remapping_core_h(CS%remap_cs, nz_data, tmp_h( 1:nz_data),tmp_val2, & + nz , h(i,j,1:nz ),tmp_val1, & h_neglect, h_neglect_edge) - CS%var(m)%p(i,j,1:nz) = CS%var(m)%p(i,j,1:nz) + q*tmp_val1(1:nz) + ! add increment to tracer on model h + CS%Var(m)%p(i,j,1:nz) = CS%Var(m)%p(i,j,1:nz) + q*tmp_val1(1:nz) + ! store increment for diagnostics if (m == 1) then tmp_t(i,j,1:nz) = tmp_val1(1:nz) else tmp_s(i,j,1:nz) = tmp_val1(1:nz) endif - enddo; enddo - deallocate(tmp_val2) - enddo - ! bound salinity values ! check if it is correct to do that or if it hides - ! other problems ... - do j=js,je ; do i=is,ie - do k=1,nz - CS%var(2)%p(i,j,k) = max(0.0 , CS%var(2)%p(i,j,k)) enddo - enddo ; enddo - - ! add increments to u - allocate(hu_obs(nz_data)) ; hu_obs(:)=0.0 ; hu(:) = 0.0; - allocate(tmp_val2(nz_data));tmp_val2(:)=0.0 - tmp_u(:,:,:) = 0.0 ! diagnostics - do j=js,je ; do i=isB,ieB - tmp_val2(1:nz_data) = CS%Ref_oda_u%p(i,j,1:nz_data) - hu_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i+1,j,1:nz_data)) - hu(1:nz) = 0.5*( h(i,j,1:nz )+ h(i+1,j,1:nz )) - call remapping_core_h(CS%remap_cs, nz_data, hu_obs(1:nz_data), tmp_val2, & - nz , hu(1:nz ), tmp_val1, & - h_neglect, h_neglect_edge) - CS%var_u%p(i,j,1:nz) = CS%var_u%p(i,j,1:nz) + q*tmp_val1(1:nz) - tmp_u(i,j,1:nz) = tmp_val1(1:nz) - enddo; enddo - deallocate(tmp_val2,hu_obs) - - ! add increments to v - allocate(hv_obs(nz_data)) ; hv_obs(:)=0.0 ; hv(:) = 0.0; - allocate(tmp_val2(nz_data));tmp_val2(:)=0.0 - tmp_v(:,:,:) = 0.0 ! diagnostics - do j=jsB,jeB ; do i=is,ie - tmp_val2(1:nz_data) = CS%Ref_oda_v%p(i,j,1:nz_data) - hv_obs(1:nz_data) = 0.5 * (h_obs(i,j,1:nz_data)+h_obs(i,j+1,1:nz_data)) - hv(1:nz) = 0.5 * ( h(i,j,1:nz )+ h(i,j+1,1:nz )) - call remapping_core_h(CS%remap_cs, nz_data, hv_obs(1:nz_data), tmp_val2, & - nz , hv(1:nz ), tmp_val1, & - h_neglect, h_neglect_edge) - CS%var_v%p(i,j,1:nz) = CS%var_v%p(i,j,1:nz) + q*tmp_val1(1:nz) - tmp_v(i,j,1:nz) = tmp_val1(1:nz) + ! bound salinity values ! check if it is correct to do that or if it hides + ! other problems ... + do k=1,nz + CS%Var(2)%p(i,j,k) = max(0.0 , CS%Var(2)%p(i,j,k)) + enddo enddo; enddo - deallocate(tmp_val2,hv_obs) + + + ! add u and v increments + if (CS%uv_inc) then + + ! add increments to u + hu(:) = 0.0 + tmp_u(:,:,:) = 0.0 ! diagnostics + do j=js,je ; do i=isB,ieB + ! get increment + tmp_val2(1:nz_data) = CS%Inc_u%p(i,j,1:nz_data) + ! get the h and h_obs at u points + hu_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i+1,j,1:nz_data)) + hu(1:nz) = 0.5*( h(i,j,1:nz )+ h(i+1,j,1:nz )) + ! account for different SSH + hu_obs(1:nz_data) = (sum(hu(1:nz))/sum(hu_obs(1:nz_data))) * hu_obs(1:nz_data) + ! remap increment profile on hu + call remapping_core_h(CS%remap_cs, nz_data, hu_obs(1:nz_data), tmp_val2, & + nz , hu(1:nz ), tmp_val1, & + h_neglect, h_neglect_edge) + ! add increment to u-velocity on hu + CS%Var_u%p(i,j,1:nz) = CS%Var_u%p(i,j,1:nz) + q*tmp_val1(1:nz) + ! store increment for diagnostics + tmp_u(i,j,1:nz) = tmp_val1(1:nz) + enddo; enddo + + ! add increments to v + hv(:) = 0.0 + tmp_v(:,:,:) = 0.0 ! diagnostics + do j=jsB,jeB ; do i=is,ie + ! get increment + tmp_val2(1:nz_data) = CS%Inc_v%p(i,j,1:nz_data) + ! get the h and h_obs at v points + hv_obs(1:nz_data) = 0.5 * (h_obs(i,j,1:nz_data)+h_obs(i,j+1,1:nz_data)) + hv(1:nz) = 0.5 * ( h(i,j,1:nz )+ h(i,j+1,1:nz )) + ! account for different SSH + hv_obs(1:nz_data) = (sum(hv(1:nz))/sum(hv_obs(1:nz_data))) * hv_obs(1:nz_data) + ! remap increment profile on hv + call remapping_core_h(CS%remap_cs, nz_data, hv_obs(1:nz_data), tmp_val2, & + nz , hv(1:nz ), tmp_val1, & + h_neglect, h_neglect_edge) + ! add increment to v-velocity on hv + CS%Var_v%p(i,j,1:nz) = CS%Var_v%p(i,j,1:nz) + q*tmp_val1(1:nz) + ! store increment for diagnostics + tmp_v(i,j,1:nz) = tmp_val1(1:nz) + enddo; enddo + + endif ! uv_inc ! Diagnostics of increments, mostly used for debugging. - if (CS%id_u_oda_inc > 0) call post_data(CS%id_u_oda_inc, tmp_u, CS%diag) - if (CS%id_v_oda_inc > 0) call post_data(CS%id_v_oda_inc, tmp_v, CS%diag) + if (CS%uv_inc) then + if (CS%id_u_oda_inc > 0) call post_data(CS%id_u_oda_inc, tmp_u, CS%diag) + if (CS%id_v_oda_inc > 0) call post_data(CS%id_v_oda_inc, tmp_v, CS%diag) + endif if (CS%id_h_oda_inc > 0) call post_data(CS%id_h_oda_inc, h , CS%diag) if (CS%id_T_oda_inc > 0) call post_data(CS%id_T_oda_inc, tmp_t, CS%diag) if (CS%id_S_oda_inc > 0) call post_data(CS%id_S_oda_inc, tmp_s, CS%diag) + ! deallocate arrays + deallocate(tmp_h,tmp_val2,hu_obs,hv_obs) + deallocate(h_obs) -end subroutine apply_oda_incupd + end subroutine apply_oda_incupd !> Initialize diagnostics for the oda_incupd module. subroutine init_oda_incupd_diags(Time, G, GV, diag, CS, US) @@ -763,7 +800,7 @@ subroutine oda_incupd_end(CS) if (.not.associated(CS)) return do m=1,CS%fldno - if (associated(CS%Ref_val(m)%p)) deallocate(CS%Ref_val(m)%p) + if (associated(CS%Inc(m)%p)) deallocate(CS%Inc(m)%p) enddo deallocate(CS) From 97b2ca7d73aa2c630c4a2204b31dd4597a66959a Mon Sep 17 00:00:00 2001 From: abozec Date: Fri, 18 Jun 2021 09:46:47 -0400 Subject: [PATCH 5/9] Updates to ease readability of the code and ensure restart reproducibility MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - remove the u,v,tv pointers from the set_up_oda_incupd subroutines  - explicitly put tv,u,v in arguments of apply_oda_incupd  - change name  get_oda_increments to calc_oda_increments  - explicitly put tv,u,v in arguments of calc_oda_increments  - use loops to calculate sum of h’s - remove the 2018_answers options - added option to output increments when using full fields  - added subroutine output_oda_incupd_inc to write a separate double-precision file for the increments. - add if (mask==1 ) condition to apply/calculate the increments  -add pass_var,pass_vector on increments and u,v,tv at the end of calc_oda_increments and apply_oda_incupd respectively - added CS%ncount to restart file to have reproducibility on restart. - Added option ODA_INCUPD_RESET_NCOUNT (! default= True, If True, reinitialize number of updates already done.) --- .../MOM_state_initialization.F90 | 109 +-- src/ocean_data_assim/MOM_oda_incupd.F90 | 651 +++++++++--------- .../vertical/MOM_diabatic_driver.F90 | 8 +- 3 files changed, 416 insertions(+), 352 deletions(-) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index 48a0a02ef1..322bc8faef 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -94,9 +94,9 @@ module MOM_state_initialization use MOM_remapping, only : remapping_CS, initialize_remapping use MOM_remapping, only : remapping_core_h use MOM_horizontal_regridding, only : horiz_interp_and_extrap_tracer -use MOM_oda_incupd, only: oda_incupd_CS, initialize_oda_incupd +use MOM_oda_incupd, only: oda_incupd_CS, initialize_oda_incupd_fixed, initialize_oda_incupd use MOM_oda_incupd, only: set_up_oda_incupd_field, set_up_oda_incupd_vel_field -use MOM_oda_incupd, only: get_oda_increments, apply_oda_incupd +use MOM_oda_incupd, only: calc_oda_increments, output_oda_incupd_inc implicit none ; private @@ -482,6 +482,16 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & dt=dt, initial=.true.) endif endif + + ! Initialized assimilative incremental update (oda_incupd) structure and + ! register restart. + call get_param(PF, mdl, "ODA_INCUPD", use_oda_incupd, & + "If true, oda incremental updates will be applied "//& + "everywhere in the domain.", default=.false.) + if (use_oda_incupd) then + call initialize_oda_incupd_fixed(G, GV, US, oda_incupd_CSp, restart_CS) + endif + ! This is the end of the block of code that might have initialized fields ! internally at the start of a new run. @@ -619,15 +629,12 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & if (debug_OBC) call open_boundary_test_extern_h(G, GV, OBC, h) call callTree_leave('MOM_initialize_state()') - ! Data Assimilation with incremental update - call get_param(PF, mdl, "ODA_INCUPD", use_oda_incupd, & - "If true, oda incremental updates will be applied "//& - "everywhere in the domain.", default=.false.) + ! Set-up of data Assimilation with incremental update if (use_oda_incupd) then call initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, & - PF, oda_incupd_CSp) + PF, oda_incupd_CSp, restart_CS, Time) endif - + end subroutine MOM_initialize_state @@ -2036,26 +2043,31 @@ end function separate_idamp_for_uv end subroutine initialize_sponges_file -subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, param_file, oda_incupd_CSp) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - logical, intent(in) :: use_temperature !< If true, T & S are state variables. - type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic +subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, param_file, & + oda_incupd_CSp, restart_CS, Time) + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + logical, intent(in) :: use_temperature !< If true, T & S are state variables. + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic !! variables. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: h !< Layer thickness [H ~> m or kg m-2] (in) + intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] (in) - real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: u !< The zonal velocity that is being - !! initialized [L T-1 ~> m s-1] - real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & - intent(in) :: v !< The meridional velocity that is being - !! initialized [L T-1 ~> m s-1] + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: u !< The zonal velocity that is being + !! initialized [L T-1 ~> m s-1] + real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + intent(in) :: v !< The meridional velocity that is being + !! initialized [L T-1 ~> m s-1] type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters. type(oda_incupd_CS), pointer :: oda_incupd_CSp !< A pointer that is set to point to the control !! structure for this module. - + type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control + !! structure. + type(time_type), intent(in) :: Time !< Time at the start of the run segment. Time_in + !! overrides any value set for + !Time. ! Local variables real, allocatable, dimension(:,:,:) :: hoda ! The layer thk inc. and oda layer thk [H ~> m or kg m-2]. real, allocatable, dimension(:,:,:) :: tmp_tr ! A temporary array for reading oda fields @@ -2066,15 +2078,17 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, p integer, dimension(4) :: siz integer :: nz_data ! The size of the sponge source grid - logical :: oda_inc ! input files are increments (true) or full fields (false) - logical :: uv_inc ! use u and v increments + logical :: oda_inc ! input files are increments (true) or full fields (false) + logical :: save_inc ! save increments if using full fields + logical :: uv_inc ! use u and v increments + logical :: reset_ncount ! reset ncount to zero if true character(len=40) :: tempinc_var, salinc_var, uinc_var, vinc_var, h_var character(len=40) :: mdl = "initialize_oda_incupd_file" character(len=200) :: inc_file, uv_inc_file ! Strings for filenames character(len=200) :: filename, inputdir ! Strings for file/path and path. - logical :: use_ALE ! True if ALE is being used, False if in layered mode +! logical :: use_ALE ! True if ALE is being used, False if in layered mode is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -2088,6 +2102,18 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, p call get_param(param_file, mdl, "ODA_INCUPD_INC", oda_inc, & "INCUPD files are increments and not full fields.", & default=.true.) + if (.not.oda_inc) then + call get_param(param_file, mdl, "ODA_INCUPD_SAVE", save_inc, & + "If true, save the increments when using full fields.", & + default=.false.) + endif + call get_param(param_file, mdl, "ODA_INCUPD_RESET_NCOUNT", reset_ncount, & + "If True, reinitialize number of updates already done, ncount.",& + default=.true.) + if (.not.oda_inc .and. .not.reset_ncount) & + call MOM_error(FATAL, " initialize_oda_incupd: restarting during update "// & + "necessitates increments input file") + call get_param(param_file, mdl, "ODA_TEMPINC_VAR", tempinc_var, & "The name of the potential temperature inc. variable in "//& "ODA_INCUPD_FILE.", default="ptemp_inc") @@ -2110,7 +2136,7 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, p "The name of the meridional vel. inc. variable in "//& "ODA_INCUPD_FILE.", default="v_inc") - call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, do_not_log = .true.) +! call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, do_not_log = .true.) ! Read in incremental update for tracers filename = trim(inputdir)//trim(inc_file) @@ -2123,22 +2149,24 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, p call MOM_error(FATAL,"initialize_oda_incupd_file: Array size mismatch for oda data.") nz_data = siz(3) ! get h increments - allocate(hoda(isd:ied,jsd:jed,nz)) + allocate(hoda(isd:ied,jsd:jed,nz_data)) call MOM_read_data(filename, h_var , hoda(:,:,:), G%Domain, scale=US%m_to_Z) - call initialize_oda_incupd( G, GV, US, param_file, oda_incupd_CSp, hoda, nz_data) + call initialize_oda_incupd( G, GV, US, param_file, oda_incupd_CSp, hoda, nz_data, restart_CS) deallocate(hoda) - ! get T and S increments + ! set-up T and S increments arrays if (use_temperature) then allocate(tmp_tr(isd:ied,jsd:jed,nz_data)) + ! temperature inc. in array Inc(1) call MOM_read_data(filename, tempinc_var, tmp_tr(:,:,:), G%Domain) - call set_up_oda_incupd_field(tmp_tr, G, GV, tv%T, oda_incupd_CSp) + call set_up_oda_incupd_field(tmp_tr, G, GV, oda_incupd_CSp) + ! salinity inc. in array Inc(2) call MOM_read_data(filename, salinc_var, tmp_tr(:,:,:), G%Domain) - call set_up_oda_incupd_field(tmp_tr, G, GV, tv%S, oda_incupd_CSp) + call set_up_oda_incupd_field(tmp_tr, G, GV, oda_incupd_CSp) deallocate(tmp_tr) endif - - ! get U and V increments + + ! set-up U and V increments arrays if (uv_inc) then filename = trim(inputdir)//trim(uv_inc_file) call log_param(param_file, mdl, "INPUTDIR/ODA_INCUPD_UV_FILE", filename) @@ -2148,21 +2176,24 @@ subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, p allocate(tmp_v(isd:ied,G%JsdB:G%JedB,nz_data)) tmp_u(:,:,:) = 0.0 ; tmp_v(:,:,:) = 0.0 call MOM_read_vector(filename, uinc_var, vinc_var, tmp_u, tmp_v, G%Domain,scale=US%m_s_to_L_T) - call set_up_oda_incupd_vel_field(tmp_u, tmp_v, G, GV, u, v, oda_incupd_CSp) + call set_up_oda_incupd_vel_field(tmp_u, tmp_v, G, GV, oda_incupd_CSp) deallocate(tmp_u,tmp_v) endif - - ! calculate increments if full fields + + ! calculate increments if input are full fields if (oda_inc) then ! input are increments if (is_root_pe()) call MOM_error(NOTE,"incupd using increments fields ") - else ! inputs are full fields + else ! inputs are full fields if (is_root_pe()) call MOM_error(NOTE,"incupd using full fields ") - call get_oda_increments(h, G, GV, US, oda_incupd_CSp) - ! call apply_oda_incupd(h, 0.0, G, GV, US, oda_incupd_CSp) + call calc_oda_increments(h, tv, u, v, G, GV, US, oda_incupd_CSp) + if (save_inc) then + call output_oda_incupd_inc(Time, G, GV, param_file, oda_incupd_CSp, US) + endif endif ! not oda_inc end subroutine initialize_oda_incupd_file + !> This subroutine sets the 4 bottom depths at velocity points to be the !! maximum of the adjacent depths. subroutine set_velocity_depth_max(G) diff --git a/src/ocean_data_assim/MOM_oda_incupd.F90 b/src/ocean_data_assim/MOM_oda_incupd.F90 index e86c8be2d6..006e9dec03 100644 --- a/src/ocean_data_assim/MOM_oda_incupd.F90 +++ b/src/ocean_data_assim/MOM_oda_incupd.F90 @@ -2,29 +2,37 @@ !! from data assimilation. !! !! Applying incremental updates requires the following: -!! 1. initialize_oda_incupd +!! 1. initialize_oda_incupd_fixed and initialize_oda_incupd !! 2. set_up_oda_incupd_field (tracers) and set_up_oda_incupd_vel_field (vel) -!! 3. apply_oda_incupd -!! 5. oda_incupd_end (not being used for now) +!! 3. calc_oda_increments (if using full fields input) +!! 4. apply_oda_incupd +!! 5. output_oda_incupd_inc (output increment if using full fields input) +!! 6. init_oda_incupd_diags (to output increments in diagnostics) +!! 7. oda_incupd_end (not being used for now) module MOM_oda_incupd ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_array_transform, only: rotate_array -use MOM_coms, only : sum_across_PEs -use MOM_diag_mediator, only : post_data, query_averaging_enabled, register_diag_field -use MOM_diag_mediator, only : diag_ctrl -use MOM_domains, only : pass_var,pass_vector -use MOM_error_handler, only : MOM_error, FATAL, NOTE, WARNING, is_root_pe -use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_grid, only : ocean_grid_type -use MOM_remapping, only : remapping_cs, remapping_core_h, initialize_remapping -use MOM_spatial_means, only : global_i_mean -use MOM_time_manager, only : time_type -use MOM_unit_scaling, only : unit_scale_type -use MOM_verticalGrid, only : verticalGrid_type -use MOM_verticalGrid, only : get_thickness_units +use MOM_array_transform, only : rotate_array +use MOM_coms, only : sum_across_PEs +use MOM_diag_mediator, only : post_data, query_averaging_enabled, register_diag_field +use MOM_diag_mediator, only : diag_ctrl +use MOM_domains, only : pass_var,pass_vector +use MOM_error_handler, only : MOM_error, FATAL, NOTE, WARNING, is_root_pe +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_get_input, only : directories, Get_MOM_input +use MOM_grid, only : ocean_grid_type +use MOM_io, only : vardesc, var_desc +use MOM_remapping, only : remapping_cs, remapping_core_h, initialize_remapping +use MOM_restart, only : register_restart_field, register_restart_pair, MOM_restart_CS +use MOM_restart, only : restart_init, save_restart, query_initialized +use MOM_spatial_means, only : global_i_mean +use MOM_time_manager, only : time_type +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : thermo_var_ptrs +use MOM_verticalGrid, only : verticalGrid_type +use MOM_verticalGrid, only : get_thickness_units use mpp_io_mod, only : mpp_get_axis_length use mpp_io_mod, only : axistype @@ -36,9 +44,8 @@ module MOM_oda_incupd ! Publicly available functions public set_up_oda_incupd_field, set_up_oda_incupd_vel_field -public initialize_oda_incupd, apply_oda_incupd, oda_incupd_end -public init_oda_incupd_diags,get_oda_increments -!public rotate_oda_incupd, update_oda_incupd_field +public initialize_oda_incupd_fixed, initialize_oda_incupd, apply_oda_incupd, oda_incupd_end +public init_oda_incupd_diags,calc_oda_increments,output_oda_incupd_inc ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with @@ -77,27 +84,17 @@ module MOM_oda_incupd integer :: fldno = 0 !< The number of fields which have already been !! registered by calls to set_up_oda_incupd_field - type(p3d) :: Var(MAX_FIELDS_) !< Pointers to the fields that are being updated. type(p3d) :: Inc(MAX_FIELDS_) !< The increments to be applied to the field - type(p3d) :: Var_u !< Pointer to the u-velocities. that are being updated. type(p3d) :: Inc_u !< The increments to be applied to the u-velocities - type(p3d) :: Var_v !< Pointer to the v-velocities. that are being updated. type(p3d) :: Inc_v !< The increments to be applied to the v-velocities - type(p3d) :: Ref_h !< Vertical grid on which the increments are provided + type(p3d) :: Ref_h !< Vertical grid on which the increments are provided integer :: nstep_incupd !< number of time step for full update - integer :: ncount !< increment time step counter + real :: ncount = 0.0 !< increment time step counter type(remapping_cs) :: remap_cs !< Remapping parameters and work arrays - logical :: remap_answers_2018 !< If true, use the order of arithmetic and expressions that - !! recover the answers for remapping from the end of 2018. - !! Otherwise, use more robust forms of the same expressions. - logical :: hor_regrid_answers_2018 !< If true, use the order of arithmetic for horizonal regridding - !! that recovers the answers from the end of 2018. Otherwise, use - !! rotationally symmetric forms of the same expressions. - logical :: incupdDataOngrid !< True if the incupd data are on the model horizontal grid - logical :: uv_inc !< use u and v increments + logical :: uv_inc !< use u and v increments ! for diagnostics type(diag_ctrl), pointer :: diag ! This subroutine defined the number of time step for full update, stores the layer pressure +!> This subroutine defined the control structure of module and register +!the time counter to full update in restart +subroutine initialize_oda_incupd_fixed( G, GV, US, CS, restart_CS) + + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(oda_incupd_CS), pointer :: CS !< A pointer that is set to point to the control + !! structure for this module (in/out). + type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control structure. + + +! This include declares and sets the variable "version". +#include "version_variable.h" + character(len=256) :: mesg + if (associated(CS)) then + call MOM_error(WARNING, "initialize_oda_incupd_fixed called with an associated "// & + "control structure.") + return + endif + allocate(CS) + + ! initialize time counter + CS%ncount = 0.0 + ! register ncount in restart + call register_restart_field(CS%ncount, "oda_incupd_ncount", .false., restart_CS,& + "Number of inc. update already done", "N/A") + + +end subroutine initialize_oda_incupd_fixed + + +!> This subroutine defined the number of time step for full update, stores the layer pressure !! increments and initialize remap structure. -subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data) +subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data, restart_CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure @@ -124,21 +153,25 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data) !! to parse for model parameter values. type(oda_incupd_CS), pointer :: CS !< A pointer that is set to point to the control !! structure for this module (in/out). - real, dimension(SZI_(G),SZJ_(G),nz_data), intent(in) :: data_h !< The ODA h + real, dimension(SZI_(G),SZJ_(G),nz_data), intent(in) :: data_h !< The ODA h !! [H ~> m or kg m-2]. + type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control + !! structure. + ! This include declares and sets the variable "version". #include "version_variable.h" character(len=40) :: mdl = "MOM_oda" ! This module's name. logical :: use_oda_incupd logical :: bndExtrapolation = .true. ! If true, extrapolate boundaries - logical :: default_2018_answers - integer :: i, j, k + logical :: reset_ncount + integer :: i, j, k real :: nhours_incupd, dt, dt_therm - character(len=256) :: mesg + type(vardesc) :: vd + character(len=256) :: mesg character(len=10) :: remapScheme - if (associated(CS)) then - call MOM_error(WARNING, "initialize_oda_incupd_fixed called with an associated "// & + if (.not.associated(CS)) then + call MOM_error(WARNING, "initialize_oda_incupd called without an associated "// & "control structure.") return endif @@ -151,11 +184,12 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data) if (.not.use_oda_incupd) return - allocate(CS) - call get_param(param_file, mdl, "ODA_INCUPD_NHOURS", nhours_incupd, & "Number of hours for full update (0=direct insertion).", & default=3.0) + call get_param(param_file, mdl, "ODA_INCUPD_RESET_NCOUNT", reset_ncount, & + "If True, reinitialize number of updates already done, ncount.", & + default=.true.) call get_param(param_file, mdl, "DT", dt, & "The (baroclinic) dynamics time step. The time-step that "//& "is actually used will be an integer fraction of the "//& @@ -184,13 +218,6 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data) "than PCM. E.g., if PPM is used for remapping, a "//& "PPM reconstruction will also be used within boundary cells.", & default=.false., do_not_log=.true.) - call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, & - "This sets the default value for the various _2018_ANSWERS parameters.", & - default=.false.) - call get_param(param_file, mdl, "REMAPPING_2018_ANSWERS", CS%remap_answers_2018, & - "If true, use the order of arithmetic and expressions that recover the "//& - "answers from the end of 2018. Otherwise, use updated and more robust "//& - "forms of the same expressions.", default=default_2018_answers) call get_param(param_file, mdl, "ODA_INCUPD_DATA_ONGRID", CS%incupdDataOngrid, & "When defined, the incoming oda_incupd data are "//& "assumed to be on the model horizontal grid " , & @@ -206,9 +233,9 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data) 'The oda_incupd code only applies ODA increments on the same horizontal grid. ') ! get number of timestep for full update - if (nhours_incupd == 0) then + if (nhours_incupd == 0) then CS%nstep_incupd = 1 !! direct insertion - else + else CS%nstep_incupd = floor(nhours_incupd*3600./dt_therm + 0.001 ) - 1 endif write(mesg,'(i12)') CS%nstep_incupd @@ -216,38 +243,45 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data) call MOM_error(NOTE,"initialize_oda_incupd: Number of Timestep of inc. update:"//& trim(mesg)) - ! initialize time counter - CS%ncount=0 + ! number of inc. update already done, CS%ncount, either from restart or set to 0.0 + if (query_initialized(CS%ncount, "oda_incupd_ncount", restart_CS) .and. & + .not.reset_ncount) then + CS%ncount = CS%ncount + else + CS%ncount = 0.0 + endif + write(mesg,'(f4.1)') CS%ncount + if (is_root_pe()) & + call MOM_error(NOTE,"initialize_oda_incupd: Inc. update already done:"//& + trim(mesg)) ! get the vertical grid (h_obs) of the increments CS%nz_data = nz_data allocate(CS%Ref_h%p(G%isd:G%ied,G%jsd:G%jed,CS%nz_data)) CS%Ref_h%p(:,:,:) = 0.0 ; do j=G%jsc,G%jec; do i=G%isc,G%iec ; do k=1,CS%nz_data - CS%Ref_h%p(i,j,k) = data_h(i,j,k) + CS%Ref_h%p(i,j,k) = data_h(i,j,k) enddo; enddo ; enddo ! Call the constructor for remapping control structure call initialize_remapping(CS%remap_cs, remapScheme, boundary_extrapolation=bndExtrapolation, & - answers_2018=CS%remap_answers_2018) + answers_2018=.false.) end subroutine initialize_oda_incupd -!> This subroutine stores theincrements at h points for the variable +!> This subroutine stores the increments at h points for the variable !! whose address is given by f_ptr. -subroutine set_up_oda_incupd_field(sp_val, G, GV, f_ptr, CS) +subroutine set_up_oda_incupd_field(sp_val, G, GV, CS) type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(oda_incupd_CS), pointer :: CS !< oda_incupd control structure (in/out). real, dimension(SZI_(G),SZJ_(G),CS%nz_data), & intent(in) :: sp_val !< increment field, it can have an !! arbitrary number of layers. - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - target, intent(in) :: f_ptr !< Pointer to the field to be updated - integer :: i, j, k + integer :: i, j, k character(len=256) :: mesg ! String for error messages if (.not.associated(CS)) return @@ -269,28 +303,24 @@ subroutine set_up_oda_incupd_field(sp_val, G, GV, f_ptr, CS) CS%Inc(CS%fldno)%p(i,j,k) = sp_val(i,j,k) enddo enddo ; enddo - CS%Var(CS%fldno)%p => f_ptr ! pointer to model tracers end subroutine set_up_oda_incupd_field !> This subroutine stores the increments at u and v points for the variable !! whose address is given by u_ptr and v_ptr. -subroutine set_up_oda_incupd_vel_field(u_val, v_val, G, GV, u_ptr, v_ptr, CS) +subroutine set_up_oda_incupd_vel_field(u_val, v_val, G, GV, CS) type(ocean_grid_type), intent(in) :: G !< Grid structure (in). type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(oda_incupd_CS), pointer :: CS !< oda incupd structure (in/out). - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + real, dimension(SZIB_(G),SZJ_(G),CS%nz_data), & intent(in) :: u_val !< u increment, it has arbritary number of layers but !! not to exceed the total number of model layers - real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + real, dimension(SZI_(G),SZJB_(G),CS%nz_data), & intent(in) :: v_val !< v increment, it has arbritary number of layers but !! not to exceed the number of model layers - real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u_ptr !< u ptr to the field to be assimilated - real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v_ptr !< v ptr to the field to be assimilated - - integer :: i, j, k + integer :: i, j, k if (.not.associated(CS)) return @@ -303,7 +333,6 @@ subroutine set_up_oda_incupd_vel_field(u_val, v_val, G, GV, u_ptr, v_ptr, CS) CS%Inc_u%p(i,j,k) = u_val(i,j,k) enddo enddo ; enddo - CS%Var_u%p => u_ptr ! pointer to model u-velocity ! store the increment/full field v profile allocate(CS%Inc_v%p(G%isd:G%ied,G%jsdB:G%jedB,CS%nz_data)) @@ -313,17 +342,25 @@ subroutine set_up_oda_incupd_vel_field(u_val, v_val, G, GV, u_ptr, v_ptr, CS) CS%Inc_v%p(i,j,k) = v_val(i,j,k) enddo enddo ; enddo - CS%Var_v%p => v_ptr ! pointer to model v-velocity end subroutine set_up_oda_incupd_vel_field -subroutine get_oda_increments(h, G, GV, US, CS) +! calculation of the increments if using full fields (ODA_INCUPD_INC=.false.) at initialization +subroutine calc_oda_increments(h, tv, u, v, G, GV, US, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure (in). type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: h !< Layer thickness [H ~> m or kg m-2] (in) + intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] (in) + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables + + real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: u !< The zonal velocity that is being + !! initialized [L T-1 ~> m s-1] + real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + intent(in) :: v !< The meridional velocity that is being + !! initialized [L T-1 ~> m s-1] type(oda_incupd_CS), pointer :: CS !< A pointer to the control structure for this module !! that is set by a previous call to initialize_oda_incupd (in). @@ -336,21 +373,26 @@ subroutine get_oda_increments(h, G, GV, US, CS) real, dimension(SZK_(GV)) :: hu, hv ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] - integer :: m, i, j, k, is, ie, js, je, nz, nz_data + integer :: i, j, k, is, ie, js, je, nz, nz_data integer :: isB, ieB, jsB, jeB real :: h_neglect, h_neglect_edge + real :: sum_h1, sum_h2 !vertical sums of h's character(len=256) :: mesg is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB if (.not.associated(CS)) return - if (.not.CS%remap_answers_2018) then - h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff - elseif (GV%Boussinesq) then + + ! increments calculated on if CS%ncount = 0.0 + if (CS%ncount /= 0.0) call MOM_error(FATAL,'calc_oda_increments: '// & + 'CS%ncount should be 0.0 to get accurate increments.') + + + if (GV%Boussinesq) then h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 else - h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 + h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff endif ! get h_obs @@ -361,6 +403,7 @@ subroutine get_oda_increments(h, G, GV, US, CS) enddo ; enddo ; enddo call pass_var(h_obs,G%Domain) + ! allocate 1-d arrays allocate(tmp_h(nz_data)); tmp_h(:)=0.0 allocate(tmp_val2(nz_data)) ; tmp_val2(:) = 0.0 @@ -371,71 +414,129 @@ subroutine get_oda_increments(h, G, GV, US, CS) tmp_val1(:)=0.0 do j=js,je ; do i=is,ie ! account for the different SSH - tmp_h(1:nz_data)=(sum(h(i,j,1:nz))/sum(h_obs(i,j,1:nz_data)))*h_obs(i,j,1:nz_data) - do m=1,CS%fldno - ! get tracer - tmp_val1(1:nz) = CS%Var(m)%p(i,j,1:nz) - ! remap tracer on h_obs - call remapping_core_h(CS%remap_cs, nz , h(i,j,1:nz ), tmp_val1, & - nz_data, tmp_h( 1:nz_data), tmp_val2, & - h_neglect, h_neglect_edge) - ! get increment from full field on h_obs - CS%Inc(m)%p(i,j,1:nz_data) = CS%Inc(m)%p(i,j,1:nz_data) - tmp_val2(1:nz_data) + sum_h1=0.0 + do k=1,nz + sum_h1 = sum_h1+h(i,j,k) + enddo + sum_h2=0.0 + do k=1,nz_data + sum_h2 = sum_h2+h_obs(i,j,k) + enddo + do k=1,nz_data + tmp_h(k)=(sum_h1/sum_h2)*h_obs(i,j,k) enddo + if (G%mask2dT(i,j) == 1) then + ! get temperature + tmp_val1(1:nz) = tv%T(i,j,1:nz) + ! remap tracer on h_obs + call remapping_core_h(CS%remap_cs, nz , h(i,j,1:nz ), tmp_val1, & + nz_data, tmp_h( 1:nz_data), tmp_val2, & + h_neglect, h_neglect_edge) + ! get increment from full field on h_obs + CS%Inc(1)%p(i,j,1:nz_data) = CS%Inc(1)%p(i,j,1:nz_data) - tmp_val2(1:nz_data) + + ! get salinity + tmp_val1(1:nz) = tv%S(i,j,1:nz) + ! remap tracer on h_obs + call remapping_core_h(CS%remap_cs, nz , h(i,j,1:nz ), tmp_val1, & + nz_data, tmp_h( 1:nz_data), tmp_val2, & + h_neglect, h_neglect_edge) + ! get increment from full field on h_obs + CS%Inc(2)%p(i,j,1:nz_data) = CS%Inc(2)%p(i,j,1:nz_data) - tmp_val2(1:nz_data) + endif enddo; enddo ! remap u to h_obs to get increment - if (CS%uv_inc) then + if (CS%uv_inc) then + call pass_var(h, G%Domain) + hu(:)=0.0 do j=js,je ; do i=isB,ieB + if (G%mask2dCu(i,j) == 1) then ! get u-velocity - tmp_val1(1:nz) = CS%Var_u%p(i ,j,1:nz) + tmp_val1(1:nz) = u(i,j,1:nz) ! get the h and h_obs at u points hu(1:nz) = 0.5*( h(i,j,1:nz )+ h(i+1,j,1:nz )) hu_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i+1,j,1:nz_data)) - ! account for the different SSH - hu_obs(1:nz_data) = (sum(hu(1:nz))/sum(hu_obs(1:nz_data))) * hu_obs(1:nz_data) + ! account for the different SSH + sum_h1=0.0 + do k=1,nz + sum_h1 = sum_h1+hu(k) + enddo + sum_h2=0.0 + do k=1,nz_data + sum_h2 = sum_h2+hu_obs(k) + enddo + do k=1,nz_data + hu_obs(k)=(sum_h1/sum_h2)*hu_obs(k) + enddo ! remap model u on hu_obs call remapping_core_h(CS%remap_cs, nz , hu(1:nz ), tmp_val1, & nz_data, hu_obs(1:nz_data), tmp_val2, & h_neglect, h_neglect_edge) ! get increment from full field on h_obs CS%Inc_u%p(i,j,1:nz_data) = CS%Inc_u%p(i,j,1:nz_data) - tmp_val2(1:nz_data) + endif enddo; enddo - + ! remap v to h_obs to get increment hv(:) = 0.0; do j=jsB,jeB ; do i=is,ie + if (G%mask2dCv(i,j) == 1) then ! get v-velocity - tmp_val1(1:nz) = CS%Var_v%p(i ,j,1:nz) - ! get the h and h_obs at v points + tmp_val1(1:nz) = v(i,j,1:nz) + ! get the h and h_obs at v points hv(1:nz) = 0.5*( h(i,j,1:nz )+ h(i,j+1,1:nz )) hv_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i,j+1,1:nz_data)) ! account for the different SSH - hv_obs(1:nz_data) = (sum(hv(1:nz))/sum(hv_obs(1:nz_data))) *hv_obs(1:nz_data) - ! remap model v on hv_obs + sum_h1=0.0 + do k=1,nz + sum_h1 = sum_h1+hv(k) + enddo + sum_h2=0.0 + do k=1,nz_data + sum_h2 = sum_h2+hv_obs(k) + enddo + do k=1,nz_data + hv_obs(k)=(sum_h1/sum_h2)*hv_obs(k) + enddo + ! remap model v on hv_obs call remapping_core_h(CS%remap_cs, nz , hv(1:nz ), tmp_val1, & nz_data, hv_obs(1:nz_data), tmp_val2, & h_neglect, h_neglect_edge) ! get increment from full field on h_obs CS%Inc_v%p(i,j,1:nz_data) = CS%Inc_v%p(i,j,1:nz_data) - tmp_val2(1:nz_data) + endif enddo; enddo endif ! uv_inc - + + call pass_var(CS%Inc(1)%p, G%Domain) + call pass_var(CS%Inc(2)%p, G%Domain) + call pass_vector(CS%Inc_u%p,CS%Inc_v%p,G%Domain) + ! deallocate arrays deallocate(tmp_h,tmp_val2,hu_obs,hv_obs) deallocate(h_obs) -end subroutine +end subroutine calc_oda_increments !> This subroutine applies oda increments to layers thicknesses, temp, salt, U !! and V everywhere . -subroutine apply_oda_incupd(h, dt, G, GV, US, CS) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure (in). +subroutine apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure (in). type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: h !< Layer thickness [H ~> m or kg m-2] (in) + intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] (in) + type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic variables + + real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), & + intent(inout) :: u !< The zonal velocity that is being + !! initialized [L T-1 ~> m s-1] + real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), & + intent(inout) :: v !< The meridional velocity that is being + !! initialized [L T-1 ~> m s-1] + real, intent(in) :: dt !< The amount of time covered by this call [T ~> s]. type(oda_incupd_CS), pointer :: CS !< A pointer to the control structure for this module !! that is set by a previous call to initialize_oda_incupd (in). @@ -445,49 +546,47 @@ subroutine apply_oda_incupd(h, dt, G, GV, US, CS) real, dimension(SZK_(GV)) :: tmp_val1 ! data values remapped to model grid real, dimension(SZK_(GV)) :: hu, hv ! A column of thicknesses at h, u or v points [H ~> m or kg m-2] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tmp_t !< A temporary array for t inc. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tmp_t !< A temporary array for t inc. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tmp_s !< A temporary array for s inc. - real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: tmp_u !< A temporary array for u inc. + real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: tmp_u !< A temporary array for u inc. real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: tmp_v !< A temporary array for v inc. real, allocatable, dimension(:,:,:) :: h_obs !< h of increments - real, allocatable, dimension(:) :: tmp_h !< temporary array for corrected h_obs + real, allocatable, dimension(:) :: tmp_h !< temporary array for corrected h_obs real, allocatable, dimension(:) :: hu_obs,hv_obs ! A column of thicknesses at h points [H ~> m or kg m-2] - integer :: m, i, j, k, is, ie, js, je, nz, nz_data + integer :: i, j, k, is, ie, js, je, nz, nz_data integer :: isB, ieB, jsB, jeB - integer :: ncount ! time step counter - real :: q +! integer :: ncount ! time step counter + real :: q ! weight of the update real :: h_neglect, h_neglect_edge - character(len=256) :: mesg + real :: sum_h1, sum_h2 !vertical sums of h's + character(len=256) :: mesg is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke isB = G%iscB ; ieB = G%iecB ; jsB = G%jscB ; jeB = G%jecB if (.not.associated(CS)) return - ! update counter - CS%ncount = CS%ncount+1 - q = 1.0/CS%nstep_incupd - ! no assimilation after CS%step_incupd - if (CS%ncount > CS%nstep_incupd) then + if (CS%ncount >= CS%nstep_incupd) then if (is_root_pe()) call MOM_error(NOTE,"ended updating fields with increments. ") return endif !ncount>CS%nstep_incupd + ! update counter + CS%ncount = CS%ncount+1.0 + q = 1.0/CS%nstep_incupd ! print out increments - write(mesg,'(i8)') CS%ncount + write(mesg,'(f10.0)') CS%ncount if (is_root_pe()) call MOM_error(NOTE,"updating fields with increments ncount:"//trim(mesg)) write(mesg,'(f10.8)') q if (is_root_pe()) call MOM_error(NOTE,"updating fields with weight q:"//trim(mesg)) - if (.not.CS%remap_answers_2018) then - h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff - elseif (GV%Boussinesq) then + if (GV%Boussinesq) then h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 else - h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 + h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff endif ! get h_obs @@ -509,77 +608,122 @@ subroutine apply_oda_incupd(h, dt, G, GV, US, CS) tmp_t(:,:,:) = 0.0 ; tmp_s(:,:,:) = 0.0 ! diagnostics do j=js,je ; do i=is,ie ! account for the different SSH - tmp_h(1:nz_data)=(sum(h(i,j,1:nz))/sum(h_obs(i,j,1:nz_data)))*h_obs(i,j,1:nz_data) - do m=1,CS%fldno - ! get increment - tmp_val2(1:nz_data) = CS%Inc(m)%p(i,j,1:nz_data) - ! remap increment profile on model h - call remapping_core_h(CS%remap_cs, nz_data, tmp_h( 1:nz_data),tmp_val2, & - nz , h(i,j,1:nz ),tmp_val1, & - h_neglect, h_neglect_edge) - ! add increment to tracer on model h - CS%Var(m)%p(i,j,1:nz) = CS%Var(m)%p(i,j,1:nz) + q*tmp_val1(1:nz) - ! store increment for diagnostics - if (m == 1) then - tmp_t(i,j,1:nz) = tmp_val1(1:nz) - else - tmp_s(i,j,1:nz) = tmp_val1(1:nz) - endif + sum_h1=0.0 + do k=1,nz + sum_h1 = sum_h1+h(i,j,k) enddo + sum_h2=0.0 + do k=1,nz_data + sum_h2 = sum_h2+h_obs(i,j,k) + enddo + do k=1,nz_data + tmp_h(k)=(sum_h1/sum_h2)*h_obs(i,j,k) + enddo + if (G%mask2dT(i,j) == 1) then + ! get temperature increment + tmp_val2(1:nz_data) = CS%Inc(1)%p(i,j,1:nz_data) + ! remap increment profile on model h + call remapping_core_h(CS%remap_cs, nz_data, tmp_h( 1:nz_data),tmp_val2, & + nz , h(i,j,1:nz ),tmp_val1, & + h_neglect, h_neglect_edge) + ! add increment to tracer on model h + tv%T(i,j,1:nz) = tv%T(i,j,1:nz) + q*tmp_val1(1:nz) + tmp_t(i,j,1:nz) = tmp_val1(1:nz) ! store T increment for diagnostics + + ! get salinity increment + tmp_val2(1:nz_data) = CS%Inc(2)%p(i,j,1:nz_data) + ! remap increment profile on model h + call remapping_core_h(CS%remap_cs, nz_data, tmp_h( 1:nz_data),tmp_val2,& + nz , h(i,j,1:nz ),tmp_val1,& + h_neglect, h_neglect_edge) + ! add increment to tracer on model h + tv%S(i,j,1:nz) = tv%S(i,j,1:nz) + q*tmp_val1(1:nz) + tmp_s(i,j,1:nz) = tmp_val1(1:nz) ! store S increment for diagnostics ! bound salinity values ! check if it is correct to do that or if it hides ! other problems ... do k=1,nz - CS%Var(2)%p(i,j,k) = max(0.0 , CS%Var(2)%p(i,j,k)) + tv%S(i,j,k) = max(0.0 , tv%S(i,j,k)) enddo + endif enddo; enddo ! add u and v increments if (CS%uv_inc) then + call pass_var(h,G%Domain) ! to ensure reproducibility + ! add increments to u hu(:) = 0.0 tmp_u(:,:,:) = 0.0 ! diagnostics do j=js,je ; do i=isB,ieB - ! get increment + if (G%mask2dCu(i,j) == 1) then + ! get u increment tmp_val2(1:nz_data) = CS%Inc_u%p(i,j,1:nz_data) ! get the h and h_obs at u points hu_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i+1,j,1:nz_data)) hu(1:nz) = 0.5*( h(i,j,1:nz )+ h(i+1,j,1:nz )) ! account for different SSH - hu_obs(1:nz_data) = (sum(hu(1:nz))/sum(hu_obs(1:nz_data))) * hu_obs(1:nz_data) + sum_h1=0.0 + do k=1,nz + sum_h1 = sum_h1+hu(k) + enddo + sum_h2=0.0 + do k=1,nz_data + sum_h2 = sum_h2+hu_obs(k) + enddo + do k=1,nz_data + hu_obs(k)=(sum_h1/sum_h2)*hu_obs(k) + enddo ! remap increment profile on hu call remapping_core_h(CS%remap_cs, nz_data, hu_obs(1:nz_data), tmp_val2, & nz , hu(1:nz ), tmp_val1, & h_neglect, h_neglect_edge) ! add increment to u-velocity on hu - CS%Var_u%p(i,j,1:nz) = CS%Var_u%p(i,j,1:nz) + q*tmp_val1(1:nz) + u(i,j,1:nz) = u(i,j,1:nz) + q*tmp_val1(1:nz) ! store increment for diagnostics tmp_u(i,j,1:nz) = tmp_val1(1:nz) + endif enddo; enddo ! add increments to v hv(:) = 0.0 tmp_v(:,:,:) = 0.0 ! diagnostics do j=jsB,jeB ; do i=is,ie - ! get increment + if (G%mask2dCv(i,j) == 1) then + ! get v increment tmp_val2(1:nz_data) = CS%Inc_v%p(i,j,1:nz_data) ! get the h and h_obs at v points hv_obs(1:nz_data) = 0.5 * (h_obs(i,j,1:nz_data)+h_obs(i,j+1,1:nz_data)) hv(1:nz) = 0.5 * ( h(i,j,1:nz )+ h(i,j+1,1:nz )) ! account for different SSH - hv_obs(1:nz_data) = (sum(hv(1:nz))/sum(hv_obs(1:nz_data))) * hv_obs(1:nz_data) + sum_h1=0.0 + do k=1,nz + sum_h1 = sum_h1+hv(k) + enddo + sum_h2=0.0 + do k=1,nz_data + sum_h2 = sum_h2+hv_obs(k) + enddo + do k=1,nz_data + hv_obs(k)=(sum_h1/sum_h2)*hv_obs(k) + enddo ! remap increment profile on hv call remapping_core_h(CS%remap_cs, nz_data, hv_obs(1:nz_data), tmp_val2, & nz , hv(1:nz ), tmp_val1, & h_neglect, h_neglect_edge) ! add increment to v-velocity on hv - CS%Var_v%p(i,j,1:nz) = CS%Var_v%p(i,j,1:nz) + q*tmp_val1(1:nz) + v(i,j,1:nz) = v(i,j,1:nz) + q*tmp_val1(1:nz) ! store increment for diagnostics tmp_v(i,j,1:nz) = tmp_val1(1:nz) + endif enddo; enddo - endif ! uv_inc + endif ! uv_inc + + call pass_var(tv%T, G%Domain) + call pass_var(tv%S, G%Domain) + call pass_vector(u,v,G%Domain) ! Diagnostics of increments, mostly used for debugging. if (CS%uv_inc) then @@ -596,6 +740,61 @@ subroutine apply_oda_incupd(h, dt, G, GV, US, CS) end subroutine apply_oda_incupd +!> Output increment if using full fields for the oda_incupd module. +subroutine output_oda_incupd_inc(Time, G, GV, param_file, CS, US) + type(time_type), target, intent(in) :: Time !< The current model time + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure + type(param_file_type), intent(in) :: param_file !< A structure indicating the open file + !! to parse for + !model parameter + !values. + type(oda_incupd_CS), pointer :: CS !< ODA incupd control structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling + + type(MOM_restart_CS), pointer :: restart_CSp_tmp => NULL() + + type(directories) :: dirs + type(vardesc) :: u_desc, v_desc + + character(len=40) :: mdl = "MOM_oda" ! This module's name. + character(len=200) :: inc_file ! name of the increment file + + if (.not.associated(CS)) return + ! get the output_directory + call Get_MOM_Input(dirs=dirs) + if (is_root_pe()) call MOM_error(NOTE,"output increments in output_directory") + + ! get a restart structure + call restart_init(param_file, restart_CSp_tmp) + + ! register the variables to write + call register_restart_field(CS%Inc(1)%p, "T_inc", .true., restart_CSp_tmp, & + "Pot. T. increment", "degC") + call register_restart_field(CS%Inc(2)%p, "S_inc", .true., restart_CSp_tmp, & + "Salinity increment", "psu") + call register_restart_field(CS%Ref_h%p, "h_obs", .true., restart_CSp_tmp, & + "Observational h", "m") + if (CS%uv_inc) then + u_desc = var_desc("u_inc", "m s-1", "U-vel increment", hor_grid='Cu') + v_desc = var_desc("v_inc", "m s-1", "V-vel increment", hor_grid='Cv') + call register_restart_pair(CS%Inc_u%p, CS%Inc_v%p, u_desc, v_desc, & + .false., restart_CSp_tmp) + endif + + ! get the name of the output file + call get_param(param_file, mdl, "ODA_INCUPD_OUTPUT_FILE", inc_file,& + "The name-root of the output file for the increment if using full fields.", & + default="MOM.inc") + + ! write the increments file + call save_restart(dirs%output_directory, Time, G, restart_CSp_tmp, & + filename=inc_file, GV=GV, write_ic=.true.) + +end subroutine output_oda_incupd_inc + + + !> Initialize diagnostics for the oda_incupd module. subroutine init_oda_incupd_diags(Time, G, GV, diag, CS, US) type(time_type), target, intent(in) :: Time !< The current model time @@ -624,172 +823,6 @@ subroutine init_oda_incupd_diags(Time, G, GV, diag, CS, US) end subroutine init_oda_incupd_diags -!> Rotate the incr. fields from the input to the model index map. -!> Not implemented yet -!subroutine rotate_oda_incupd(oda_incupd_in, G_in, oda_incupd, G, GV, turns, param_file) -! type(oda_incupd_CS), intent(in) :: oda_incupd_in !< The control structure for this module with the -! !! original grid rotation -! type(ocean_grid_type), intent(in) :: G_in !< The ocean's grid structure with the original rotation. -! type(oda_incupd_CS), pointer :: oda_incup !< A pointer to the control that will be set up with -! !! the new grid rotation -! type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure with the new rotation. -! type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure -! integer, intent(in) :: turns !< The number of 90-degree turns between grids -! type(param_file_type), intent(in) :: param_file !< A structure indicating the open file -! !! to parse for model parameter values. -! -! ! First part: Index construction -! ! 1. Reconstruct Iresttime(:,:) from sponge_in -! ! 2. rotate Iresttime(:,:) -! ! 3. Call initialize_oda_incupd using new grid and rotated Iresttime(:,:) -! ! All the index adjustment should follow from the Iresttime rotation -! -! real, dimension(:,:), allocatable :: Iresttime_in, Iresttime -! real, dimension(:,:,:), allocatable :: data_h_in, data_h -! real, dimension(:,:,:), allocatable :: sp_val_in, sp_val -! real, dimension(:,:,:), pointer :: sp_ptr => NULL() -! integer :: c, c_i, c_j -! integer :: k, nz_data -! integer :: n -! logical :: fixed_sponge -! -! fixed_sponge = .not. sponge_in%time_varying_sponges -! ! NOTE: nz_data is only conditionally set when fixed_sponge is true. -! -! allocate(Iresttime_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed)) -! allocate(Iresttime(G%isd:G%ied, G%jsd:G%jed)) -! Iresttime_in(:,:) = 0.0 -! -! if (fixed_sponge) then -! nz_data = sponge_in%nz_data -! allocate(data_h_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed, nz_data)) -! allocate(data_h(G%isd:G%ied, G%jsd:G%jed, nz_data)) -! data_h_in(:,:,:) = 0. -! endif -! -! ! Re-populate the 2D Iresttime and data_h arrays on the original grid -! do c=1,sponge_in%num_col -! c_i = sponge_in%col_i(c) -! c_j = sponge_in%col_j(c) -! Iresttime_in(c_i, c_j) = sponge_in%Iresttime_col(c) -! if (fixed_sponge) then -! do k = 1, nz_data -! data_h(c_i, c_j, k) = sponge_in%Ref_h%p(k,c) -! enddo -! endif -! enddo -! -! call rotate_array(Iresttime_in, turns, Iresttime) -! if (fixed_sponge) then -! call rotate_array(data_h_in, turns, data_h) -! call initialize_oda_incupd_fixed(Iresttime, G, GV, param_file, sponge, & -! data_h, nz_data) -! else -! call initialize_oda_incupd_varying(Iresttime, G, GV, param_file, sponge) -! endif -! -! deallocate(Iresttime_in) -! deallocate(Iresttime) -! if (fixed_sponge) then -! deallocate(data_h_in) -! deallocate(data_h) -! endif -! -! ! Second part: Provide rotated fields for which relaxation is applied -! -! sponge%fldno = sponge_in%fldno -! -! if (fixed_sponge) then -! allocate(sp_val_in(G_in%isd:G_in%ied, G_in%jsd:G_in%jed, nz_data)) -! allocate(sp_val(G%isd:G%ied, G%jsd:G%jed, nz_data)) -! endif -! -! do n=1,sponge_in%fldno -! ! Assume that tracers are pointers and are remapped in other functions(?) -! sp_ptr => sponge_in%var(n)%p -! if (fixed_sponge) then -! sp_val_in(:,:,:) = 0.0 -! do c = 1, sponge_in%num_col -! c_i = sponge_in%col_i(c) -! c_j = sponge_in%col_j(c) -! do k = 1, nz_data -! sp_val_in(c_i, c_j, k) = sponge_in%Ref_val(n)%p(k,c) -! enddo -! enddo -! -! call rotate_array(sp_val_in, turns, sp_val) -! -! ! NOTE: This points sp_val with the unrotated field. See note below. -! call set_up_oda_incupd_field(sp_val, G, GV, sp_ptr, sponge) -! -! deallocate(sp_val_in) -! else -! ! We don't want to repeat FMS init in set_up_oda_incupd_field_varying() -! ! (time_interp_external_init, init_external_field, etc), so we manually -! ! do a portion of this function below. -! sponge%Ref_val(n)%id = sponge_in%Ref_val(n)%id -! sponge%Ref_val(n)%num_tlevs = sponge_in%Ref_val(n)%num_tlevs -! -! nz_data = sponge_in%Ref_val(n)%nz_data -! sponge%Ref_val(n)%nz_data = nz_data -! -! allocate(sponge%Ref_val(n)%p(nz_data, sponge_in%num_col)) -! allocate(sponge%Ref_val(n)%h(nz_data, sponge_in%num_col)) -! sponge%Ref_val(n)%p(:,:) = 0.0 -! sponge%Ref_val(n)%h(:,:) = 0.0 -! -! ! TODO: There is currently no way to associate a generic field pointer to -! ! its rotated equivalent without introducing a new data structure which -! ! explicitly tracks the pairing. -! ! -! ! As a temporary fix, we store the pointer to the unrotated field in -! ! the rotated sponge, and use this reference to replace the pointer -! ! to the rotated field update_oda_incupd field. -! ! -! ! This makes a lot of unverifiable assumptions, and should not be -! ! considered the final solution. -! sponge%var(n)%p => sp_ptr -! endif -! enddo -! -! ! TODO: var_u and var_v sponge dampling is not yet supported. -! if (associated(sponge_in%var_u%p) .or. associated(sponge_in%var_v%p)) & -! call MOM_error(FATAL, "Rotation of ALE sponge velocities is not yet " & -! // "implemented.") -! -! ! Transfer any existing diag_CS reference pointer -! sponge%diag => sponge_in%diag -! -! ! NOTE: initialize_oda_incupd_* resolves remap_cs -!end subroutine rotate_oda_incupd - - -!> Scan the oda_incupd variables and replace a prescribed pointer to a new value. -!> Not implemented yet -! TODO: This function solely exists to replace field pointers in the sponge -! after rotation. This function is part of a temporary solution until -! something more robust is developed. -!subroutine update_oda_incupd_field(sponge, p_old, G, GV, p_new) -! type(oda_incupd_CS), pointer :: sponge !< A pointer to the control structure for this module -! !! that is set by a previous call to initialize_oda_incupd. -! real, dimension(:,:,:), & -! target, intent(in) :: p_old !< The previous array of target values -! type(ocean_grid_type), intent(in) :: G !< The updated ocean grid structure -! type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure -! real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & -! target, intent(in) :: p_new !< The new array of target values -! -! integer :: n -! -! do n=1,sponge%fldno -! if (associated(sponge%var(n)%p, p_old)) sponge%var(n)%p => p_new -! enddo -! -!end subroutine update_oda_incupd_field - - -! GMM: I could not find where sponge_end is being called, but I am keeping -! oda_incupd_end here so we can add that if needed. !> This subroutine deallocates any memory associated with the oda_incupd module. subroutine oda_incupd_end(CS) type(oda_incupd_CS), pointer :: CS !< A pointer to the control structure that is diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index f7bb028aab..9551c5dced 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -1009,7 +1009,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Tim if (CS%use_oda_incupd .and. associated(CS%oda_incupd_CSp)) then call MOM_mesg("Starting ODA_INCUPD legacy ", 5) call cpu_clock_begin(id_clock_oda_incupd) - call apply_oda_incupd(h, dt, G, GV, US, CS%oda_incupd_CSp) + call apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS%oda_incupd_CSp) call cpu_clock_end(id_clock_oda_incupd) if (CS%debug) then call MOM_state_chksum("apply_oda_incupd ", u, v, h, G, GV, US, haloshift=0) @@ -1506,7 +1506,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, if (CS%use_oda_incupd .and. associated(CS%oda_incupd_CSp)) then call MOM_mesg("Starting ODA_INCUPD ", 5) call cpu_clock_begin(id_clock_oda_incupd) - call apply_oda_incupd(h, dt, G, GV, US, CS%oda_incupd_CSp) + call apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS%oda_incupd_CSp) call cpu_clock_end(id_clock_oda_incupd) if (CS%debug) then call MOM_state_chksum("apply_oda_incupd ", u, v, h, G, GV, US, haloshift=0) @@ -2301,7 +2301,7 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e ! Apply data assimilation incremental update -oda_incupd- if (CS%use_oda_incupd .and. associated(CS%oda_incupd_CSp)) then call cpu_clock_begin(id_clock_oda_incupd) - call apply_oda_incupd(h, dt, G, GV, US, CS%oda_incupd_CSp) + call apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS%oda_incupd_CSp) call cpu_clock_end(id_clock_oda_incupd) if (CS%debug) then call MOM_state_chksum("apply_oda_incupd ", u, v, h, G, GV, US, haloshift=0) @@ -2828,7 +2828,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di !! tracer flow control module type(sponge_CS), pointer :: sponge_CSp !< pointer to the sponge module control structure type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< pointer to the ALE sponge module control structure - type(oda_incupd_CS), pointer :: oda_incupd_CSp !< pointer to the oda incupd module control structure + type(oda_incupd_CS), pointer :: oda_incupd_CSp !< pointer to the oda incupd module control structure real :: Kd ! A diffusivity used in the default for other tracer diffusivities, in MKS units [m2 s-1] integer :: num_mode From 4ee64d1bf32eaf741495cc59fe6df5bdce010cab Mon Sep 17 00:00:00 2001 From: pjpegion Date: Thu, 24 Jun 2021 20:02:04 +0000 Subject: [PATCH 6/9] remove write_ic=.true. from save_restart call --- src/ocean_data_assim/MOM_oda_incupd.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ocean_data_assim/MOM_oda_incupd.F90 b/src/ocean_data_assim/MOM_oda_incupd.F90 index 006e9dec03..9f32e70d36 100644 --- a/src/ocean_data_assim/MOM_oda_incupd.F90 +++ b/src/ocean_data_assim/MOM_oda_incupd.F90 @@ -789,7 +789,7 @@ subroutine output_oda_incupd_inc(Time, G, GV, param_file, CS, US) ! write the increments file call save_restart(dirs%output_directory, Time, G, restart_CSp_tmp, & - filename=inc_file, GV=GV, write_ic=.true.) + filename=inc_file, GV=GV) !, write_ic=.true.) end subroutine output_oda_incupd_inc From d0564907d2bec522c390a917cf6a83f37ed7f2e9 Mon Sep 17 00:00:00 2001 From: pjpegion Date: Thu, 1 Jul 2021 14:41:03 +0000 Subject: [PATCH 7/9] code cleanup and add rescaling for ODA_INCUPD_NHOURS --- src/ocean_data_assim/MOM_oda_incupd.F90 | 423 +++++++++++++----------- 1 file changed, 227 insertions(+), 196 deletions(-) diff --git a/src/ocean_data_assim/MOM_oda_incupd.F90 b/src/ocean_data_assim/MOM_oda_incupd.F90 index 9f32e70d36..12d18587bf 100644 --- a/src/ocean_data_assim/MOM_oda_incupd.F90 +++ b/src/ocean_data_assim/MOM_oda_incupd.F90 @@ -186,7 +186,7 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data, res call get_param(param_file, mdl, "ODA_INCUPD_NHOURS", nhours_incupd, & "Number of hours for full update (0=direct insertion).", & - default=3.0) + default=3.0,units="h", scale=US%s_to_T) call get_param(param_file, mdl, "ODA_INCUPD_RESET_NCOUNT", reset_ncount, & "If True, reinitialize number of updates already done, ncount.", & default=.true.) @@ -236,7 +236,7 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data, res if (nhours_incupd == 0) then CS%nstep_incupd = 1 !! direct insertion else - CS%nstep_incupd = floor(nhours_incupd*3600./dt_therm + 0.001 ) - 1 + CS%nstep_incupd = floor( nhours_incupd * 3600. / dt_therm + 0.001 ) - 1 endif write(mesg,'(i12)') CS%nstep_incupd if (is_root_pe()) & @@ -298,11 +298,9 @@ subroutine set_up_oda_incupd_field(sp_val, G, GV, CS) CS%Inc(CS%fldno)%nz_data = CS%nz_data allocate(CS%Inc(CS%fldno)%p(G%isd:G%ied,G%jsd:G%jed,CS%nz_data)) CS%Inc(CS%fldno)%p(:,:,:) = 0.0 - do j=G%jsc,G%jec ; do i=G%isc,G%iec - do k=1,CS%nz_data - CS%Inc(CS%fldno)%p(i,j,k) = sp_val(i,j,k) - enddo - enddo ; enddo + do k=1,CS%nz_data ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + CS%Inc(CS%fldno)%p(i,j,k) = sp_val(i,j,k) + enddo ; enddo ; enddo end subroutine set_up_oda_incupd_field @@ -316,10 +314,10 @@ subroutine set_up_oda_incupd_vel_field(u_val, v_val, G, GV, CS) real, dimension(SZIB_(G),SZJ_(G),CS%nz_data), & intent(in) :: u_val !< u increment, it has arbritary number of layers but - !! not to exceed the total number of model layers + !! not to exceed the total number of model layers [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),CS%nz_data), & intent(in) :: v_val !< v increment, it has arbritary number of layers but - !! not to exceed the number of model layers + !! not to exceed the number of model layers [L T-1 ~> m s-1] integer :: i, j, k if (.not.associated(CS)) return @@ -397,7 +395,7 @@ subroutine calc_oda_increments(h, tv, u, v, G, GV, US, CS) ! get h_obs nz_data = CS%Inc(1)%nz_data - allocate(h_obs(G%isd:G%ied,G%jsd:G%jed,nz_data)) ; h_obs(:,:,:)=0.0 + allocate(h_obs(G%isd:G%ied,G%jsd:G%jed,nz_data)) ; h_obs(:,:,:) = 0.0 do k=1,nz_data ; do j=js,je ; do i=is,ie h_obs(i,j,k) = CS%Ref_h%p(i,j,k) enddo ; enddo ; enddo @@ -405,109 +403,128 @@ subroutine calc_oda_increments(h, tv, u, v, G, GV, US, CS) ! allocate 1-d arrays - allocate(tmp_h(nz_data)); tmp_h(:)=0.0 + allocate(tmp_h(nz_data)); tmp_h(:) = 0.0 allocate(tmp_val2(nz_data)) ; tmp_val2(:) = 0.0 - allocate(hu_obs(nz_data)) ; hu_obs(:)=0.0 - allocate(hv_obs(nz_data)) ; hv_obs(:)=0.0 + allocate(hu_obs(nz_data)) ; hu_obs(:) = 0.0 + allocate(hv_obs(nz_data)) ; hv_obs(:) = 0.0 ! remap t,s (on h_init) to h_obs to get increment - tmp_val1(:)=0.0 + tmp_val1(:) = 0.0 do j=js,je ; do i=is,ie - ! account for the different SSH - sum_h1=0.0 - do k=1,nz - sum_h1 = sum_h1+h(i,j,k) - enddo - sum_h2=0.0 - do k=1,nz_data - sum_h2 = sum_h2+h_obs(i,j,k) - enddo - do k=1,nz_data - tmp_h(k)=(sum_h1/sum_h2)*h_obs(i,j,k) - enddo if (G%mask2dT(i,j) == 1) then - ! get temperature - tmp_val1(1:nz) = tv%T(i,j,1:nz) - ! remap tracer on h_obs - call remapping_core_h(CS%remap_cs, nz , h(i,j,1:nz ), tmp_val1, & - nz_data, tmp_h( 1:nz_data), tmp_val2, & - h_neglect, h_neglect_edge) - ! get increment from full field on h_obs - CS%Inc(1)%p(i,j,1:nz_data) = CS%Inc(1)%p(i,j,1:nz_data) - tmp_val2(1:nz_data) - - ! get salinity - tmp_val1(1:nz) = tv%S(i,j,1:nz) - ! remap tracer on h_obs - call remapping_core_h(CS%remap_cs, nz , h(i,j,1:nz ), tmp_val1, & - nz_data, tmp_h( 1:nz_data), tmp_val2, & - h_neglect, h_neglect_edge) - ! get increment from full field on h_obs - CS%Inc(2)%p(i,j,1:nz_data) = CS%Inc(2)%p(i,j,1:nz_data) - tmp_val2(1:nz_data) + ! account for the different SSH + sum_h1 = 0.0 + sum_h2 = 0.0 + do k=1,nz + sum_h1 = sum_h1+h(i,j,k) + enddo + + do k=1,nz_data + sum_h2 = sum_h2+h_obs(i,j,k) + tmp_h(k)=(sum_h1/sum_h2)*h_obs(i,j,k) + enddo + ! get temperature + do k=1,nz + tmp_val1(k) = tv%T(i,j,k) + enddo + ! remap tracer on h_obs + call remapping_core_h(CS%remap_cs, nz, h(i,j,1:nz), tmp_val1, & + nz_data, tmp_h(1:nz_data), tmp_val2, & + h_neglect, h_neglect_edge) + ! get increment from full field on h_obs + do k=1,nz_data + CS%Inc(1)%p(i,j,k) = CS%Inc(1)%p(i,j,k) - tmp_val2(k) + enddo + + ! get salinity + do k=1,nz + tmp_val1(k) = tv%S(i,j,k) + enddo + ! remap tracer on h_obs + call remapping_core_h(CS%remap_cs, nz, h(i,j,1:nz), tmp_val1, & + nz_data, tmp_h(1:nz_data), tmp_val2, & + h_neglect, h_neglect_edge) + ! get increment from full field on h_obs + do k=1,nz_data + CS%Inc(2)%p(i,j,k) = CS%Inc(2)%p(i,j,k) - tmp_val2(k) + enddo endif enddo; enddo - + ! remap u to h_obs to get increment if (CS%uv_inc) then - call pass_var(h, G%Domain) - - hu(:)=0.0 - do j=js,je ; do i=isB,ieB - if (G%mask2dCu(i,j) == 1) then - ! get u-velocity - tmp_val1(1:nz) = u(i,j,1:nz) - ! get the h and h_obs at u points - hu(1:nz) = 0.5*( h(i,j,1:nz )+ h(i+1,j,1:nz )) - hu_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i+1,j,1:nz_data)) - ! account for the different SSH - sum_h1=0.0 - do k=1,nz - sum_h1 = sum_h1+hu(k) - enddo - sum_h2=0.0 - do k=1,nz_data - sum_h2 = sum_h2+hu_obs(k) - enddo - do k=1,nz_data - hu_obs(k)=(sum_h1/sum_h2)*hu_obs(k) - enddo - ! remap model u on hu_obs - call remapping_core_h(CS%remap_cs, nz , hu(1:nz ), tmp_val1, & - nz_data, hu_obs(1:nz_data), tmp_val2, & - h_neglect, h_neglect_edge) - ! get increment from full field on h_obs - CS%Inc_u%p(i,j,1:nz_data) = CS%Inc_u%p(i,j,1:nz_data) - tmp_val2(1:nz_data) - endif - enddo; enddo - - ! remap v to h_obs to get increment - hv(:) = 0.0; - do j=jsB,jeB ; do i=is,ie - if (G%mask2dCv(i,j) == 1) then - ! get v-velocity - tmp_val1(1:nz) = v(i,j,1:nz) - ! get the h and h_obs at v points - hv(1:nz) = 0.5*( h(i,j,1:nz )+ h(i,j+1,1:nz )) - hv_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i,j+1,1:nz_data)) - ! account for the different SSH - sum_h1=0.0 - do k=1,nz - sum_h1 = sum_h1+hv(k) - enddo - sum_h2=0.0 - do k=1,nz_data - sum_h2 = sum_h2+hv_obs(k) - enddo - do k=1,nz_data - hv_obs(k)=(sum_h1/sum_h2)*hv_obs(k) - enddo - ! remap model v on hv_obs - call remapping_core_h(CS%remap_cs, nz , hv(1:nz ), tmp_val1, & - nz_data, hv_obs(1:nz_data), tmp_val2, & - h_neglect, h_neglect_edge) - ! get increment from full field on h_obs - CS%Inc_v%p(i,j,1:nz_data) = CS%Inc_v%p(i,j,1:nz_data) - tmp_val2(1:nz_data) - endif - enddo; enddo + call pass_var(h, G%Domain) + + hu(:) = 0.0 + do j=js,je ; do i=isB,ieB + if (G%mask2dCu(i,j) == 1) then + ! get u-velocity + do k=1,nz + tmp_val1(k) = u(i,j,k) + ! get the h and h_obs at u points + hu(k) = 0.5*( h(i,j,k)+ h(i+1,j,k)) + enddo + do k=1,nz_data + hu_obs(k) = 0.5*(h_obs(i,j,k)+h_obs(i+1,j,k)) + enddo + ! account for the different SSH + sum_h1 = 0.0 + do k=1,nz + sum_h1 = sum_h1+hu(k) + enddo + sum_h2 = 0.0 + do k=1,nz_data + sum_h2 = sum_h2+hu_obs(k) + enddo + do k=1,nz_data + hu_obs(k)=(sum_h1/sum_h2)*hu_obs(k) + enddo + ! remap model u on hu_obs + call remapping_core_h(CS%remap_cs, nz, hu(1:nz), tmp_val1, & + nz_data, hu_obs(1:nz_data), tmp_val2, & + h_neglect, h_neglect_edge) + ! get increment from full field on h_obs + do k=1,nz_data + CS%Inc_u%p(i,j,k) = CS%Inc_u%p(i,j,k) - tmp_val2(k) + enddo + endif + enddo; enddo + + ! remap v to h_obs to get increment + hv(:) = 0.0; + do j=jsB,jeB ; do i=is,ie + if (G%mask2dCv(i,j) == 1) then + ! get v-velocity + do k=1,nz + tmp_val1(k) = v(i,j,k) + ! get the h and h_obs at v points + hv(k) = 0.5*(h(i,j,k)+h(i,j+1,k)) + enddo + do k=1,nz_data + hv_obs(k) = 0.5*(h_obs(i,j,k)+h_obs(i,j+1,k)) + enddo + ! account for the different SSH + sum_h1 = 0.0 + do k=1,nz + sum_h1 = sum_h1+hv(k) + enddo + sum_h2 = 0.0 + do k=1,nz_data + sum_h2 = sum_h2+hv_obs(k) + enddo + do k=1,nz_data + hv_obs(k)=(sum_h1/sum_h2)*hv_obs(k) + enddo + ! remap model v on hv_obs + call remapping_core_h(CS%remap_cs, nz, hv(1:nz), tmp_val1, & + nz_data, hv_obs(1:nz_data), tmp_val2, & + h_neglect, h_neglect_edge) + ! get increment from full field on h_obs + do k=1,nz_data + CS%Inc_v%p(i,j,k) = CS%Inc_v%p(i,j,k) - tmp_val2(k) + enddo + endif + enddo; enddo endif ! uv_inc call pass_var(CS%Inc(1)%p, G%Domain) @@ -558,7 +575,7 @@ subroutine apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS) integer :: i, j, k, is, ie, js, je, nz, nz_data integer :: isB, ieB, jsB, jeB ! integer :: ncount ! time step counter - real :: q ! weight of the update + real :: inc_wt ! weight of the update for this time-step real :: h_neglect, h_neglect_edge real :: sum_h1, sum_h2 !vertical sums of h's character(len=256) :: mesg @@ -575,13 +592,13 @@ subroutine apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS) ! update counter CS%ncount = CS%ncount+1.0 - q = 1.0/CS%nstep_incupd + inc_wt = 1.0/CS%nstep_incupd ! print out increments write(mesg,'(f10.0)') CS%ncount if (is_root_pe()) call MOM_error(NOTE,"updating fields with increments ncount:"//trim(mesg)) - write(mesg,'(f10.8)') q - if (is_root_pe()) call MOM_error(NOTE,"updating fields with weight q:"//trim(mesg)) + write(mesg,'(f10.8)') inc_wt + if (is_root_pe()) call MOM_error(NOTE,"updating fields with weight inc_wt:"//trim(mesg)) if (GV%Boussinesq) then h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 @@ -591,57 +608,61 @@ subroutine apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS) ! get h_obs nz_data = CS%Inc(1)%nz_data - allocate(h_obs(G%isd:G%ied,G%jsd:G%jed,nz_data)) ; h_obs(:,:,:)=0.0 + allocate(h_obs(G%isd:G%ied,G%jsd:G%jed,nz_data)) ; h_obs(:,:,:) = 0.0 do k=1,nz_data ; do j=js,je ; do i=is,ie h_obs(i,j,k) = CS%Ref_h%p(i,j,k) enddo ; enddo ; enddo call pass_var(h_obs,G%Domain) ! allocate 1-d array - allocate(tmp_h(nz_data)); tmp_h(:)=0.0 + allocate(tmp_h(nz_data)); tmp_h(:) = 0.0 allocate(tmp_val2(nz_data)) - allocate(hu_obs(nz_data)) ; hu_obs(:)=0.0 - allocate(hv_obs(nz_data)) ; hv_obs(:)=0.0 + allocate(hu_obs(nz_data)) ; hu_obs(:) = 0.0 + allocate(hv_obs(nz_data)) ; hv_obs(:) = 0.0 ! add increments to tracers - tmp_val1(:)=0.0 + tmp_val1(:) = 0.0 tmp_t(:,:,:) = 0.0 ; tmp_s(:,:,:) = 0.0 ! diagnostics do j=js,je ; do i=is,ie ! account for the different SSH - sum_h1=0.0 + sum_h1 = 0.0 do k=1,nz sum_h1 = sum_h1+h(i,j,k) enddo - sum_h2=0.0 + sum_h2 = 0.0 do k=1,nz_data sum_h2 = sum_h2+h_obs(i,j,k) enddo do k=1,nz_data - tmp_h(k)=(sum_h1/sum_h2)*h_obs(i,j,k) + tmp_h(k) = ( sum_h1 / sum_h2 ) * h_obs(i,j,k) enddo if (G%mask2dT(i,j) == 1) then ! get temperature increment - tmp_val2(1:nz_data) = CS%Inc(1)%p(i,j,1:nz_data) + do k=1,nz_data + tmp_val2(k) = CS%Inc(1)%p(i,j,k) + enddo ! remap increment profile on model h - call remapping_core_h(CS%remap_cs, nz_data, tmp_h( 1:nz_data),tmp_val2, & - nz , h(i,j,1:nz ),tmp_val1, & - h_neglect, h_neglect_edge) + call remapping_core_h(CS%remap_cs, nz_data, tmp_h(1:nz_data), tmp_val2, & + nz, h(i,j,1:nz),tmp_val1, h_neglect, h_neglect_edge) + do k=1,nz ! add increment to tracer on model h - tv%T(i,j,1:nz) = tv%T(i,j,1:nz) + q*tmp_val1(1:nz) - tmp_t(i,j,1:nz) = tmp_val1(1:nz) ! store T increment for diagnostics + tv%T(i,j,k) = tv%T(i,j,k) + inc_wt * tmp_val1(k) + tmp_t(i,j,k) = tmp_val1(k) ! store T increment for diagnostics + enddo ! get salinity increment - tmp_val2(1:nz_data) = CS%Inc(2)%p(i,j,1:nz_data) + do k=1,nz_data + tmp_val2(k) = CS%Inc(2)%p(i,j,k) + enddo ! remap increment profile on model h - call remapping_core_h(CS%remap_cs, nz_data, tmp_h( 1:nz_data),tmp_val2,& - nz , h(i,j,1:nz ),tmp_val1,& - h_neglect, h_neglect_edge) + call remapping_core_h(CS%remap_cs, nz_data, tmp_h(1:nz_data),tmp_val2,& + nz, h(i,j,1:nz),tmp_val1, h_neglect, h_neglect_edge) ! add increment to tracer on model h - tv%S(i,j,1:nz) = tv%S(i,j,1:nz) + q*tmp_val1(1:nz) - tmp_s(i,j,1:nz) = tmp_val1(1:nz) ! store S increment for diagnostics + do k=1,nz + tv%S(i,j,k) = tv%S(i,j,k) + inc_wt * tmp_val1(k) + tmp_s(i,j,k) = tmp_val1(k) ! store S increment for diagnostics ! bound salinity values ! check if it is correct to do that or if it hides ! other problems ... - do k=1,nz tv%S(i,j,k) = max(0.0 , tv%S(i,j,k)) enddo endif @@ -651,73 +672,83 @@ subroutine apply_oda_incupd(h, tv, u, v, dt, G, GV, US, CS) ! add u and v increments if (CS%uv_inc) then - call pass_var(h,G%Domain) ! to ensure reproducibility - - ! add increments to u - hu(:) = 0.0 - tmp_u(:,:,:) = 0.0 ! diagnostics - do j=js,je ; do i=isB,ieB - if (G%mask2dCu(i,j) == 1) then - ! get u increment - tmp_val2(1:nz_data) = CS%Inc_u%p(i,j,1:nz_data) - ! get the h and h_obs at u points - hu_obs(1:nz_data) = 0.5*(h_obs(i,j,1:nz_data)+h_obs(i+1,j,1:nz_data)) - hu(1:nz) = 0.5*( h(i,j,1:nz )+ h(i+1,j,1:nz )) - ! account for different SSH - sum_h1=0.0 - do k=1,nz - sum_h1 = sum_h1+hu(k) - enddo - sum_h2=0.0 - do k=1,nz_data - sum_h2 = sum_h2+hu_obs(k) - enddo - do k=1,nz_data - hu_obs(k)=(sum_h1/sum_h2)*hu_obs(k) - enddo - ! remap increment profile on hu - call remapping_core_h(CS%remap_cs, nz_data, hu_obs(1:nz_data), tmp_val2, & - nz , hu(1:nz ), tmp_val1, & - h_neglect, h_neglect_edge) - ! add increment to u-velocity on hu - u(i,j,1:nz) = u(i,j,1:nz) + q*tmp_val1(1:nz) - ! store increment for diagnostics - tmp_u(i,j,1:nz) = tmp_val1(1:nz) - endif - enddo; enddo - - ! add increments to v - hv(:) = 0.0 - tmp_v(:,:,:) = 0.0 ! diagnostics - do j=jsB,jeB ; do i=is,ie - if (G%mask2dCv(i,j) == 1) then - ! get v increment - tmp_val2(1:nz_data) = CS%Inc_v%p(i,j,1:nz_data) - ! get the h and h_obs at v points - hv_obs(1:nz_data) = 0.5 * (h_obs(i,j,1:nz_data)+h_obs(i,j+1,1:nz_data)) - hv(1:nz) = 0.5 * ( h(i,j,1:nz )+ h(i,j+1,1:nz )) - ! account for different SSH - sum_h1=0.0 - do k=1,nz - sum_h1 = sum_h1+hv(k) - enddo - sum_h2=0.0 - do k=1,nz_data - sum_h2 = sum_h2+hv_obs(k) - enddo - do k=1,nz_data - hv_obs(k)=(sum_h1/sum_h2)*hv_obs(k) - enddo - ! remap increment profile on hv - call remapping_core_h(CS%remap_cs, nz_data, hv_obs(1:nz_data), tmp_val2, & - nz , hv(1:nz ), tmp_val1, & - h_neglect, h_neglect_edge) - ! add increment to v-velocity on hv - v(i,j,1:nz) = v(i,j,1:nz) + q*tmp_val1(1:nz) - ! store increment for diagnostics - tmp_v(i,j,1:nz) = tmp_val1(1:nz) - endif - enddo; enddo + call pass_var(h,G%Domain) ! to ensure reproducibility + + ! add increments to u + hu(:) = 0.0 + tmp_u(:,:,:) = 0.0 ! diagnostics + do j=js,je ; do i=isB,ieB + if (G%mask2dCu(i,j) == 1) then + do k=1,nz_data + ! get u increment + tmp_val2(k) = CS%Inc_u%p(i,j,k) + ! get the h and h_obs at u points + hu_obs(k) = 0.5 * ( h_obs(i,j,k) + h_obs(i+1,j,k) ) + enddo + do k=1,nz + hu(k) = 0.5 * ( h(i,j,k) + h(i+1,j,k) ) + enddo + ! account for different SSH + sum_h1 = 0.0 + do k=1,nz + sum_h1 = sum_h1 + hu(k) + enddo + sum_h2 = 0.0 + do k=1,nz_data + sum_h2 = sum_h2 + hu_obs(k) + enddo + do k=1,nz_data + hu_obs(k)=( sum_h1 / sum_h2 ) * hu_obs(k) + enddo + ! remap increment profile on hu + call remapping_core_h(CS%remap_cs, nz_data, hu_obs(1:nz_data), tmp_val2, & + nz, hu(1:nz), tmp_val1, h_neglect, h_neglect_edge) + ! add increment to u-velocity on hu + do k=1,nz + u(i,j,k) = u(i,j,k) + inc_wt * tmp_val1(k) + ! store increment for diagnostics + tmp_u(i,j,k) = tmp_val1(k) + enddo + endif + enddo; enddo + + ! add increments to v + hv(:) = 0.0 + tmp_v(:,:,:) = 0.0 ! diagnostics + do j=jsB,jeB ; do i=is,ie + if (G%mask2dCv(i,j) == 1) then + ! get v increment + do k=1,nz_data + tmp_val2(k) = CS%Inc_v%p(i,j,k) + ! get the h and h_obs at v points + hv_obs(k) = 0.5 * ( h_obs(i,j,k) + h_obs(i,j+1,k) ) + enddo + do k=1,nz + hv(k) = 0.5 * (h(i,j,k) + h(i,j+1,k) ) + enddo + ! account for different SSH + sum_h1 = 0.0 + do k=1,nz + sum_h1 = sum_h1 + hv(k) + enddo + sum_h2 = 0.0 + do k=1,nz_data + sum_h2 = sum_h2 + hv_obs(k) + enddo + do k=1,nz_data + hv_obs(k)=( sum_h1 / sum_h2 ) * hv_obs(k) + enddo + ! remap increment profile on hv + call remapping_core_h(CS%remap_cs, nz_data, hv_obs(1:nz_data), tmp_val2, & + nz, hv(1:nz), tmp_val1, h_neglect, h_neglect_edge) + ! add increment to v-velocity on hv + do k=1,nz + v(i,j,k) = v(i,j,k) + inc_wt * tmp_val1(k) + ! store increment for diagnostics + tmp_v(i,j,k) = tmp_val1(k) + enddo + endif + enddo; enddo endif ! uv_inc From a4ca3073c0cd8b7d3a106e4f3fbbe3760c3e3cd2 Mon Sep 17 00:00:00 2001 From: pjpegion Date: Thu, 29 Jul 2021 15:01:26 +0000 Subject: [PATCH 8/9] remove white space and fix comment --- src/ocean_data_assim/MOM_oda_incupd.F90 | 2 +- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ocean_data_assim/MOM_oda_incupd.F90 b/src/ocean_data_assim/MOM_oda_incupd.F90 index 12d18587bf..334597cda7 100644 --- a/src/ocean_data_assim/MOM_oda_incupd.F90 +++ b/src/ocean_data_assim/MOM_oda_incupd.F90 @@ -450,7 +450,7 @@ subroutine calc_oda_increments(h, tv, u, v, G, GV, US, CS) enddo endif enddo; enddo - + ! remap u to h_obs to get increment if (CS%uv_inc) then call pass_var(h, G%Domain) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 86a83c58a9..fb759b13b2 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -111,7 +111,7 @@ module MOM_diabatic_driver !! domain. The exact location and properties of !! those sponges are set by calls to !! initialize_sponge and set_up_sponge_field. - logical :: use_oda_incupd !!< If True, DA incremental update is + logical :: use_oda_incupd !< If True, DA incremental update is !! applied everywhere logical :: use_geothermal !< If true, apply geothermal heating. logical :: use_int_tides !< If true, use the code that advances a separate set From cf46bc90bd04746e72ca3618e3305a2afc2b4e37 Mon Sep 17 00:00:00 2001 From: Phil Pegion <38869668+pjpegion@users.noreply.github.com> Date: Thu, 5 Aug 2021 10:20:00 -0600 Subject: [PATCH 9/9] Update MOM_oda_incupd.F90 remove unused index bounds, and fix sum_h2 loop. --- src/ocean_data_assim/MOM_oda_incupd.F90 | 21 ++------------------- 1 file changed, 2 insertions(+), 19 deletions(-) diff --git a/src/ocean_data_assim/MOM_oda_incupd.F90 b/src/ocean_data_assim/MOM_oda_incupd.F90 index 334597cda7..d3199dcb74 100644 --- a/src/ocean_data_assim/MOM_oda_incupd.F90 +++ b/src/ocean_data_assim/MOM_oda_incupd.F90 @@ -65,22 +65,6 @@ module MOM_oda_incupd type, public :: oda_incupd_CS ; private integer :: nz !< The total number of layers. integer :: nz_data !< The total number of arbritary layers (used by older code). - integer :: isc !< The starting i-index of the computational domain at h. - integer :: iec !< The ending i-index of the computational domain at h. - integer :: jsc !< The starting j-index of the computational domain at h. - integer :: jec !< The ending j-index of the computational domain at h. - integer :: IscB !< The starting I-index of the computational domain at u/v. - integer :: IecB !< The ending I-index of the computational domain at u/v. - integer :: JscB !< The starting J-index of the computational domain at u/v. - integer :: JecB !< The ending J-index of the computational domain at h. - integer :: isd !< The starting i-index of the data domain at h. - integer :: ied !< The ending i-index of the data domain at h. - integer :: jsd !< The starting j-index of the data domain at h. - integer :: jed !< The ending j-index of the data domain at h. - integer :: IsdB !< The starting I-index of the data domain at u/v. - integer :: IedB !< The ending I-index of the data domain at u/v. - integer :: JsdB !< The starting J-index of the data domain at u/v. - integer :: JedB !< The ending J-index of the data domain at h. integer :: fldno = 0 !< The number of fields which have already been !! registered by calls to set_up_oda_incupd_field @@ -224,9 +208,6 @@ subroutine initialize_oda_incupd( G, GV, US, param_file, CS, data_h,nz_data, res default=.true.) CS%nz = GV%ke - CS%isc = G%isc ; CS%iec = G%iec ; CS%jsc = G%jsc ; CS%jec = G%jec - CS%isd = G%isd ; CS%ied = G%ied ; CS%jsd = G%jsd ; CS%jed = G%jed - CS%isdB = G%isdB ; CS%iedB = G%iedB; CS%jsdB = G%jsdB ; CS%jedB = G%jedB ! increments on horizontal grid if (.not. CS%incupdDataOngrid) call MOM_error(FATAL,'initialize_oda_incupd: '// & @@ -421,6 +402,8 @@ subroutine calc_oda_increments(h, tv, u, v, G, GV, US, CS) do k=1,nz_data sum_h2 = sum_h2+h_obs(i,j,k) + enddo + do k=1,nz_data tmp_h(k)=(sum_h1/sum_h2)*h_obs(i,j,k) enddo ! get temperature