diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 9b94a96797..8b88b58e0b 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1265,7 +1265,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, & call enable_averages(dtdia, Time_end_thermo, CS%diag) if (associated(CS%odaCS)) then - call apply_oda_tracer_increments(US%T_to_s*dtdia, G, GV, tv, h, CS%odaCS) + call apply_oda_tracer_increments(US%T_to_s*dtdia, Time_end_thermo, G, GV, tv, h, CS%odaCS) endif if (associated(fluxes%p_surf) .or. associated(fluxes%p_surf_full)) then @@ -2796,7 +2796,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & (LEN_TRIM(dirs%input_filename) == 1)) if (CS%ensemble_ocean) then - call init_oda(Time, G, GV, CS%odaCS) + call init_oda(Time, G, GV, US, CS%diag, CS%odaCS) endif !### This could perhaps go here instead of in finish_MOM_initialization? diff --git a/src/framework/MOM_interp_infra.F90 b/src/framework/MOM_interp_infra.F90 index ca5b2b8516..de5cac7879 100644 --- a/src/framework/MOM_interp_infra.F90 +++ b/src/framework/MOM_interp_infra.F90 @@ -222,7 +222,7 @@ end subroutine time_interp_extern_3d !> initialize an external field integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose, & - threading, ierr, ignore_axis_atts ) + threading, ierr, ignore_axis_atts, correct_leap_year_inconsistency ) character(len=*), intent(in) :: file !< The name of the file to read character(len=*), intent(in) :: fieldname !< The name of the field in the file @@ -237,13 +237,18 @@ integer function init_extern_field(file, fieldname, MOM_domain, domain, verbose, logical, optional, intent(in) :: ignore_axis_atts !< If present and true, do not issue a !! fatal error if the axis Cartesian attribute is !! not set to a recognized value. + logical, optional, intent(in) :: correct_leap_year_inconsistency !< If present and true, + !! then allow for leap year inconsistency + if (present(MOM_Domain)) then init_extern_field = init_external_field(file, fieldname, domain=MOM_domain%mpp_domain, & - verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts) + verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & + correct_leap_year_inconsistency=correct_leap_year_inconsistency) else init_extern_field = init_external_field(file, fieldname, domain=domain, & - verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts) + verbose=verbose, threading=threading, ierr=ierr, ignore_axis_atts=ignore_axis_atts, & + correct_leap_year_inconsistency=correct_leap_year_inconsistency) endif end function init_extern_field diff --git a/src/ocean_data_assim/MOM_oda_driver.F90 b/src/ocean_data_assim/MOM_oda_driver.F90 index 8057234cdc..7a0ca291f3 100644 --- a/src/ocean_data_assim/MOM_oda_driver.F90 +++ b/src/ocean_data_assim/MOM_oda_driver.F90 @@ -9,14 +9,18 @@ module MOM_oda_driver_mod use MOM_domains, only : domain2d, global_field, get_domain_extent use MOM_domains, only : pass_var, redistribute_array, broadcast_domain use MOM_diag_mediator, only : register_diag_field, diag_axis_init, post_data +use MOM_diag_mediator, only : enable_averaging use MOM_ensemble_manager, only : get_ensemble_id, get_ensemble_size use MOM_ensemble_manager, only : get_ensemble_pelist, get_ensemble_filter_pelist use MOM_error_handler, only : stdout, stdlog, MOM_error use MOM_io, only : SINGLE_FILE +use MOM_interp_infra, only : init_extern_field, get_external_field_info +use MOM_interp_infra, only : time_interp_extern use MOM_time_manager, only : time_type, real_to_time, get_date use MOM_time_manager, only : operator(+), operator(>=), operator(/=) use MOM_time_manager, only : operator(==), operator(<) - +use MOM_cpu_clock, only : cpu_clock_begin, cpu_clock_end, cpu_clock_id +use MOM_horizontal_regridding, only : horiz_interp_and_extrap_tracer ! ODA Modules use ocean_da_types_mod, only : grid_type, ocean_profile_type, ocean_control_struct use ocean_da_core_mod, only : ocean_da_core_init, get_profiles @@ -54,6 +58,15 @@ module MOM_oda_driver_mod public :: init_oda, oda_end, set_prior_tracer, get_posterior_tracer public :: set_analysis_time, oda, save_obs_diff, apply_oda_tracer_increments +!>@{ CPU time clock ID +integer :: id_clock_oda_init +integer :: id_clock_oda_filter +integer :: id_clock_bias_adjustment +integer :: id_clock_apply_increments +integer :: id_clock_oda_prior +integer :: id_clock_oda_posterior +!>@} + #include !> A structure with a pointer to a domain2d, to allow for the creation of arrays of pointers. @@ -61,13 +74,23 @@ module MOM_oda_driver_mod type(domain2d), pointer :: mpp_domain => NULL() !< pointer to a domain2d end type ptr_mpp_domain +!> A structure containing integer handles for bias adjustment of tracers +type :: INC_CS + integer :: fldno = 0 !< The number of tracers + integer :: T_id !< The integer handle for the temperature file + integer :: S_id !< The integer handle for the salinity file +end type INC_CS + !> Control structure that contains a transpose of the ocean state across ensemble members. type, public :: ODA_CS ; private type(ocean_control_struct), pointer :: Ocean_prior=> NULL() !< ensemble ocean prior states in DA space type(ocean_control_struct), pointer :: Ocean_posterior=> NULL() !< ensemble ocean posterior states !! or increments to prior in DA space + type(ocean_control_struct), pointer :: Ocean_increment=> NULL() !< A separate structure for + !! increment diagnostics integer :: nk !< number of vertical layers used for DA type(ocean_grid_type), pointer :: Grid => NULL() !< MOM6 grid type and decomposition for the DA + type(ocean_grid_type), pointer :: G => NULL() !< MOM6 grid type and decomposition for the model type(MOM_domain_type), pointer, dimension(:) :: domains => NULL() !< Pointer to mpp_domain objects !! for ensemble members type(verticalGrid_type), pointer :: GV => NULL() !< vertical grid for DA @@ -78,12 +101,17 @@ module MOM_oda_driver_mod type(grid_type), pointer :: oda_grid !< local tracer grid real, pointer, dimension(:,:,:) :: h => NULL() ! m or kg m-2] for DA type(thermo_var_ptrs), pointer :: tv => NULL() !< pointer to thermodynamic variables + type(thermo_var_ptrs), pointer :: tv_bc => NULL() !< pointer to thermodynamic bias correction integer :: ni !< global i-direction grid size integer :: nj !< global j-direction grid size logical :: reentrant_x !< grid is reentrant in the x direction logical :: reentrant_y !< grid is reentrant in the y direction logical :: tripolar_N !< grid is folded at its north edge logical :: symmetric !< Values at C-grid locations are symmetric + logical :: use_basin_mask !< If true, use a basin file to delineate weakly coupled ocean basins + logical :: do_bias_adjustment !< If true, use spatio-temporally varying climatological tendency + !! adjustment for Temperature and Salinity + real :: bias_adjustment_multiplier !< A scaling for the bias adjustment integer :: assim_method !< Method: NO_ASSIM,EAKF_ASSIM or OI_ASSIM integer :: ensemble_size !< Size of the ensemble integer :: ensemble_id = 0 !< id of the current ensemble member @@ -99,7 +127,10 @@ module MOM_oda_driver_mod type(regridding_CS) :: regridCS !< ALE control structure for regridding type(remapping_CS) :: remapCS !< ALE control structure for remapping type(time_type) :: Time !< Current Analysis time - type(diag_ctrl) :: diag_cs ! NULL() ! initialize First_guess (prior) and Analysis grid !! information for all ensemble members -subroutine init_oda(Time, G, GV, CS) +subroutine init_oda(Time, G, GV, US, diag_CS, CS) type(time_type), intent(in) :: Time !< The current model time. type(ocean_grid_type), pointer :: G !< domain and grid information for ocean model type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), pointer :: US + type(diag_ctrl), target, intent(inout) :: diag_CS !< A pointer to a diagnostic control structure type(ODA_CS), pointer, intent(inout) :: CS !< The DA control structure ! Local variables @@ -133,6 +166,7 @@ subroutine init_oda(Time, G, GV, CS) integer :: isg,ieg,jsg,jeg integer :: idg_offset, jdg_offset integer :: stdout_unit + integer, dimension(4) :: fld_sz character(len=32) :: assim_method integer :: npes_pm, ens_info(6), ni, nj character(len=128) :: mesg @@ -140,14 +174,22 @@ subroutine init_oda(Time, G, GV, CS) character(len=30) :: coord_mode character(len=200) :: inputdir, basin_file logical :: reentrant_x, reentrant_y, tripolar_N, symmetric + character(len=80) :: bias_correction_file, inc_file if (associated(CS)) call MOM_error(FATAL, 'Calling oda_init with associated control structure') allocate(CS) + + id_clock_oda_init=cpu_clock_id('(ODA initialization)') + id_clock_oda_prior=cpu_clock_id('(ODA setting prior)') + id_clock_oda_filter=cpu_clock_id('(ODA filter computation)') + id_clock_oda_posterior=cpu_clock_id('(ODA getting posterior)') + call cpu_clock_begin(id_clock_oda_init) + ! Use ens1 parameters , this could be changed at a later time ! if it were desirable to have alternate parameters, e.g. for the grid ! for the analysis call get_MOM_input(PF,dirs,ensemble_num=0) - + CS%US=>US call unit_scaling_init(PF, CS%US) call get_param(PF, "MOM", "ASSIM_METHOD", assim_method, & @@ -167,6 +209,19 @@ subroutine init_oda(Time, G, GV, CS) "Use tripolar connectivity at the northern edge of the "//& "domain. With TRIPOLAR_N, NIGLOBAL must be even.", & default=.false.) + call get_param(PF,"MOM", "USE_TEMP_SALT_BIAS_ADJUSTMENT", CS%do_bias_adjustment, & + "If true, add a spatio-temporally varying climatological adjustment "//& + "to temperature and salinity.", & + default=.false.) + if (CS%do_bias_adjustment) then + call get_param(PF,"MOM", "BIAS_ADJUSTMENT_FACTOR", CS%bias_adjustment_multiplier, & + "A multiplicative scaling factor for the climatological tracer tendency adjustment ", & + default=1.0) + endif + call get_param(PF,"MOM", "USE_BASIN_MASK", CS%use_basin_mask, & + "If true, add a basin mask to delineate weakly connected "//& + "ocean basins for the purpose of data assimilation.", & + default=.false.) call get_param(PF,"MOM", "NIGLOBAL", CS%ni, & "The total number of thickness grid points in the "//& "x-direction in the physical domain.") @@ -200,19 +255,19 @@ subroutine init_oda(Time, G, GV, CS) call set_PElist(CS%filter_pelist) allocate(CS%domains(CS%ensemble_size)) - CS%domains(CS%ensemble_id)%mpp_domain => G%Domain%mpp_domain + CS%domains(CS%ensemble_id)%mpp_domain => G%Domain%mpp_domain ! this should go away do n=1,CS%ensemble_size if (.not. associated(CS%domains(n)%mpp_domain)) allocate(CS%domains(n)%mpp_domain) - call set_rootPE(CS%ensemble_pelist(n,1)) + call set_rootPE(CS%ensemble_pelist(n,1)) ! this line is not in Feiyu's version (needed?) call broadcast_domain(CS%domains(n)%mpp_domain) enddo - call set_rootPE(CS%filter_pelist(1)) + call set_rootPE(CS%filter_pelist(1)) ! this line is not in Feiyu's version (needed?) + CS%G => G allocate(CS%Grid) ! params NIHALO_ODA, NJHALO_ODA set the DA halo size call MOM_domains_init(CS%Grid%Domain,PF,param_suffix='_ODA') allocate(HI) - call hor_index_init(CS%Grid%Domain, HI, PF, & - local_indexing=.false.) ! Use global indexing for DA + call hor_index_init(CS%Grid%Domain, HI, PF) call verticalGridInit( PF, CS%GV, CS%US ) allocate(dG) call create_dyn_horgrid(dG, HI) @@ -233,7 +288,9 @@ subroutine init_oda(Time, G, GV, CS) call init_ocean_ensemble(CS%Ocean_prior,CS%Grid,CS%GV,CS%ensemble_size) allocate(CS%Ocean_posterior) call init_ocean_ensemble(CS%Ocean_posterior,CS%Grid,CS%GV,CS%ensemble_size) - allocate(CS%tv) + allocate(CS%Ocean_increment) + call init_ocean_ensemble(CS%Ocean_increment,CS%Grid,CS%GV,CS%ensemble_size) + call get_param(PF, 'oda_driver', "REGRIDDING_COORDINATE_MODE", coord_mode, & "Coordinate mode for vertical regridding.", & @@ -241,76 +298,75 @@ subroutine init_oda(Time, G, GV, CS) call initialize_regridding(CS%regridCS, CS%GV, CS%US, dG%max_depth,PF,'oda_driver',coord_mode,'','') call initialize_remapping(CS%remapCS,'PLM') call set_regrid_params(CS%regridCS, min_thickness=0.) + isd = G%isd; ied = G%ied; jsd = G%jsd; jed = G%jed + ! breaking with the MOM6 convention and using global indices - call get_domain_extent(G%Domain,is,ie,js,je,isd,ied,jsd,jed,& - isg,ieg,jsg,jeg,idg_offset,jdg_offset,symmetric) - isd=isd+idg_offset; ied=ied+idg_offset - jsd=jsd+jdg_offset; jed=jed+jdg_offset - !call mpp_get_data_domain(G%Domain%mpp_domain,isd,ied,jsd,jed) + !call get_domain_extent(G%Domain,is,ie,js,je,isd,ied,jsd,jed,& + ! isg,ieg,jsg,jeg,idg_offset,jdg_offset,symmetric) + !isd=isd+idg_offset; ied=ied+idg_offset ! using global indexing within the DA module + !jsd=jsd+jdg_offset; jed=jed+jdg_offset ! TODO: switch to local indexing? (mjh) + if (.not. associated(CS%h)) then - allocate(CS%h(isd:ied,jsd:jed,CS%GV%ke)); CS%h(:,:,:)=0.0 + allocate(CS%h(isd:ied,jsd:jed,CS%GV%ke)); CS%h(:,:,:)=CS%GV%Angstrom_m ! assign thicknesses call ALE_initThicknessToCoord(CS%ALE_CS,G,CS%GV,CS%h) endif + allocate(CS%tv) allocate(CS%tv%T(isd:ied,jsd:jed,CS%GV%ke)); CS%tv%T(:,:,:)=0.0 allocate(CS%tv%S(isd:ied,jsd:jed,CS%GV%ke)); CS%tv%S(:,:,:)=0.0 - call set_axes_info(CS%Grid, CS%GV, CS%US, PF, CS%diag_cs, set_vertical=.true.) - ! get domain extents for the analysis grid and use global indexing - !call get_domain_extent(CS%Grid%Domain,is,ie,js,je,isd,ied,jsd,jed,& - ! isg,ieg,jsg,jeg,idg_offset,jdg_offset,symmetric) - !isd=isd+idg_offset; ied=ied+idg_offset - !jsd=jsd+jdg_offset; jed=jed+jdg_offset - !call mpp_get_data_domain(CS%mpp_domain,isd,ied,jsd,jed) + call set_axes_info(CS%Grid, CS%GV, CS%US, PF, CS%diag_cs, set_vertical=.true.) ! missing in Feiyu's fork allocate(CS%oda_grid) CS%oda_grid%x => CS%Grid%geolonT CS%oda_grid%y => CS%Grid%geolatT - call get_param(PF, 'oda_driver', "BASIN_FILE", basin_file, & + + if (CS%use_basin_mask) then + call get_param(PF, 'oda_driver', "BASIN_FILE", basin_file, & "A file in which to find the basin masks, in variable 'basin'.", & default="basin.nc") - basin_file = trim(inputdir) // trim(basin_file) - allocate(CS%oda_grid%basin_mask(isd:ied,jsd:jed)) - CS%oda_grid%basin_mask(:,:) = 0.0 - call MOM_read_data(basin_file,'basin',CS%oda_grid%basin_mask,CS%Grid%domain, timelevel=1) - -! get global grid information from ocean_model - allocate(T_grid) - allocate(T_grid%x(CS%ni,CS%nj)) - allocate(T_grid%y(CS%ni,CS%nj)) - allocate(T_grid%basin_mask(CS%ni,CS%nj)) - call global_field(CS%mpp_domain, CS%Grid%geolonT, T_grid%x) - call global_field(CS%mpp_domain, CS%Grid%geolatT, T_grid%y) - call global_field(CS%mpp_domain, CS%oda_grid%basin_mask, T_grid%basin_mask) - T_grid%ni = CS%ni - T_grid%nj = CS%nj - T_grid%nk = CS%nk - allocate(T_grid%mask(CS%ni,CS%nj,CS%nk)) - allocate(T_grid%z(CS%ni,CS%nj,CS%nk)) - allocate(global2D(CS%ni,CS%nj)) - allocate(global2D_old(CS%ni,CS%nj)) - T_grid%mask(:,:,:) = 0.0 - T_grid%z(:,:,:) = 0.0 - - do k = 1, CS%nk - call global_field(G%Domain%mpp_domain, CS%h(:,:,k), global2D) - do i=1,CS%ni ; do j=1,CS%nj - if ( global2D(i,j) > 1 ) then - T_grid%mask(i,j,k) = 1.0 - endif - enddo ; enddo - if (k == 1) then - T_grid%z(:,:,k) = global2D/2 - else - T_grid%z(:,:,k) = T_grid%z(:,:,k-1) + (global2D + global2D_old)/2 - endif - global2D_old = global2D - enddo + basin_file = trim(inputdir) // trim(basin_file) + allocate(CS%oda_grid%basin_mask(isd:ied,jsd:jed)) + CS%oda_grid%basin_mask(:,:) = 0.0 + call MOM_read_data(basin_file,'basin',CS%oda_grid%basin_mask,CS%Grid%domain, timelevel=1) + endif - call ocean_da_core_init(CS%mpp_domain, T_grid, CS%Profiles, Time) + ! set up diag variables for analysis increments + CS%diag_CS => diag_CS + CS%id_inc_t=register_diag_field('ocean_model','temp_increment',diag_CS%axesTL,& + Time,'ocean potential temperature increments','degC') + CS%id_inc_s=register_diag_field('ocean_model','salt_increment',diag_CS%axesTL,& + Time,'ocean salinity increments','psu') + !! get global grid information from ocean model needed for ODA initialization + call set_up_global_tgrid(T_grid, CS, G) + + call ocean_da_core_init(CS%mpp_domain, T_grid, CS%Profiles, Time) + deallocate(T_grid) CS%Time=Time !! switch back to ensemble member pelist call set_PElist(CS%ensemble_pelist(CS%ensemble_id,:)) + + if (CS%do_bias_adjustment) then + inc_file = trim(inputdir) // trim(bias_correction_file) + CS%INC_CS%T_id = init_extern_field(inc_file, "temp_increment", & + correct_leap_year_inconsistency=.true.,verbose=.true.,domain=G%Domain%mpp_domain) + CS%INC_CS%S_id = init_extern_field(inc_file, "salt_increment", & + correct_leap_year_inconsistency=.true.,verbose=.true.,domain=G%Domain%mpp_domain) + call get_external_field_info(CS%INC_CS%T_id,size=fld_sz) + CS%INC_CS%fldno = 2 + if (CS%nk .ne. fld_sz(3)) call MOM_error(FATAL,'Increment levels /= ODA levels') + allocate(CS%tv_bc) ! storage for increment + allocate(CS%tv_bc%T(G%isd:G%ied,G%jsd:G%jed,CS%GV%ke)); CS%tv_bc%T(:,:,:)=0.0 + allocate(CS%tv_bc%S(G%isd:G%ied,G%jsd:G%jed,CS%GV%ke)); CS%tv_bc%S(:,:,:)=0.0 + endif + + call cpu_clock_end(id_clock_oda_init) + +! if (CS%write_obs) then +! temp_fid = open_profile_file("temp_"//trim(obs_file)) +! salt_fid = open_profile_file("salt_"//trim(obs_file)) +! end if + end subroutine init_oda !> Copy ensemble member tracers to ensemble vector. @@ -340,7 +396,8 @@ subroutine set_prior_tracer(Time, G, GV, h, tv, CS) !! switch to global pelist call set_PElist(CS%filter_pelist) - call MOM_mesg('Setting prior') + !call MOM_mesg('Setting prior') + call cpu_clock_begin(id_clock_oda_prior) ! computational domain for the analysis grid isc=CS%Grid%isc;iec=CS%Grid%iec;jsc=CS%Grid%jsc;jec=CS%Grid%jec @@ -367,6 +424,7 @@ subroutine set_prior_tracer(Time, G, GV, h, tv, CS) call pass_var(CS%Ocean_prior%S(:,:,:,m),CS%Grid%domain) enddo + call cpu_clock_end(id_clock_oda_prior) !! switch back to ensemble member pelist call set_PElist(CS%ensemble_pelist(CS%ensemble_id,:)) @@ -379,28 +437,31 @@ end subroutine set_prior_tracer subroutine get_posterior_tracer(Time, CS, h, tv, increment) type(time_type), intent(in) :: Time !< the current model time type(ODA_CS), pointer :: CS !< ocean DA control structure - real, dimension(:,:,:), pointer :: h !< Layer thicknesses [H ~> m or kg m-2] - type(thermo_var_ptrs), pointer :: tv !< A structure pointing to various thermodynamic variables + real, dimension(:,:,:), pointer, optional :: h !< Layer thicknesses [H ~> m or kg m-2] + type(thermo_var_ptrs), pointer, optional :: tv !< A structure pointing to various thermodynamic variables logical, optional, intent(in) :: increment !< True if returning increment only type(ocean_control_struct), pointer :: Ocean_increment=>NULL() integer :: i, j, m logical :: used, get_inc + integer :: seconds_per_hour = 3600. ! return if not analysis time (retain pointers for h and tv) - if (Time < CS%Time) return + if (Time < CS%Time .or. CS%assim_method .eq. NO_ASSIM) return !! switch to global pelist call set_PElist(CS%filter_pelist) call MOM_mesg('Getting posterior') - + call cpu_clock_begin(id_clock_oda_posterior) + if (present(h)) h => CS%h ! get analysis thickness + !! Calculate and redistribute increments to CS%tv right after assimilation + !! Retain CS%tv to calculate increments for IAU updates CS%tv_inc otherwise get_inc = .true. if (present(increment)) get_inc = increment if (get_inc) then allocate(Ocean_increment) - call init_ocean_ensemble(Ocean_increment,CS%Grid,CS%GV,CS%ensemble_size) Ocean_increment%T = CS%Ocean_posterior%T - CS%Ocean_prior%T Ocean_increment%S = CS%Ocean_posterior%S - CS%Ocean_prior%S endif @@ -418,17 +479,28 @@ subroutine get_posterior_tracer(Time, CS, h, tv, increment) endif enddo - tv => CS%tv - h => CS%h + if (present(tv)) tv => CS%tv + if (present(h)) h => CS%h + + call cpu_clock_end(id_clock_oda_posterior) + !! switch back to ensemble member pelist call set_PElist(CS%ensemble_pelist(CS%ensemble_id,:)) + call pass_var(CS%tv%T,CS%domains(CS%ensemble_id)) + call pass_var(CS%tv%S,CS%domains(CS%ensemble_id)) + + !convert to a tendency (degC or PSU per second) + CS%tv%T = CS%tv%T / (CS%assim_frequency * seconds_per_hour) + CS%tv%S = CS%tv%S / (CS%assim_frequency * seconds_per_hour) + + end subroutine get_posterior_tracer !> Gather observations and call ODA routines subroutine oda(Time, CS) type(time_type), intent(in) :: Time !< the current model time - type(oda_CS), intent(inout) :: CS !< the ocean DA control structure + type(oda_CS), pointer :: CS !< A pointer the ocean DA control structure integer :: i, j integer :: m @@ -438,20 +510,60 @@ subroutine oda(Time, CS) !! switch to global pelist call set_PElist(CS%filter_pelist) - + call cpu_clock_begin(id_clock_oda_filter) call get_profiles(Time, CS%Profiles, CS%CProfiles) -#ifdef ENABLE_ECDA call ensemble_filter(CS%Ocean_prior, CS%Ocean_posterior, CS%CProfiles, CS%kdroot, CS%mpp_domain, CS%oda_grid) -#endif - + call cpu_clock_end(id_clock_oda_filter) + !if (CS%write_obs) call save_obs_diff(CS%CProfiles) ! not fully implemented !! switch back to ensemble member pelist call set_PElist(CS%ensemble_pelist(CS%ensemble_id,:)) + call get_posterior_tracer(Time, CS, increment=.true.) + if (CS%do_bias_adjustment) call get_bias_correction_tracer(Time, CS) endif return end subroutine oda +subroutine get_bias_correction_tracer(Time, CS) + type(time_type), intent(in) :: Time !< the current model time + type(ODA_CS), pointer :: CS !< ocean DA control structure + + integer :: i,j,k + real, allocatable, dimension(:,:,:) :: T_bias, S_bias + real, allocatable, dimension(:,:,:) :: mask_z + real, allocatable, dimension(:), target :: z_in, z_edges_in + real :: missing_value + integer,dimension(3) :: fld_sz + + call cpu_clock_begin(id_clock_bias_adjustment) + call horiz_interp_and_extrap_tracer(CS%INC_CS%T_id,Time,1.0,CS%G,T_bias,& + mask_z,z_in,z_edges_in,missing_value,.true.,.false.,.false.,.true.) + call horiz_interp_and_extrap_tracer(CS%INC_CS%S_id,Time,1.0,CS%G,S_bias,& + mask_z,z_in,z_edges_in,missing_value,.true.,.false.,.false.,.true.) + + ! This should be replaced to use mask_z instead of the following lines + ! which are intended to zero land values using an arbitrary limit. + fld_sz=shape(T_bias) + do i=1,fld_sz(1) + do j=1,fld_sz(2) + do k=1,fld_sz(3) + if (T_bias(i,j,k) .gt. 1.0E-3) T_bias(i,j,k) = 0.0 + if (S_bias(i,j,k) .gt. 1.0E-3) S_bias(i,j,k) = 0.0 + enddo + enddo + enddo + + CS%tv_bc%T = T_bias * CS%bias_adjustment_multiplier + CS%tv_bc%S = S_bias * CS%bias_adjustment_multiplier + + call pass_var(CS%tv_bc%T, CS%domains(CS%ensemble_id)) + call pass_var(CS%tv_bc%S, CS%domains(CS%ensemble_id)) + + call cpu_clock_end(id_clock_bias_adjustment) + + end subroutine get_bias_correction_tracer + !> Finalize DA module subroutine oda_end(CS) type(ODA_CS), intent(inout) :: CS !< the ocean DA control structure @@ -541,17 +653,123 @@ end subroutine save_obs_diff !> Apply increments to tracers -subroutine apply_oda_tracer_increments(dt, G, GV, tv, h, CS) +subroutine apply_oda_tracer_increments(dt, Time_end, G, GV, tv, h, CS) real, intent(in) :: dt !< The tracer timestep [s] + type(time_type), intent(in) :: Time_end !< Time at the end of the interval type(ocean_grid_type), intent(in) :: G !< ocean grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(thermo_var_ptrs), intent(inout) :: 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] - type(ODA_CS), intent(inout) :: CS !< the data assimilation structure + type(ODA_CS), pointer :: CS !< the data assimilation structure + + !! local variables + integer :: yr, mon, day, hr, min, sec + integer :: i, j, k + integer :: isc, iec, jsc, jec + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: T_inc !< an adjustment to the temperature + !! tendency [degC T-1 -> degC s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: S_inc !< an adjustment to the salinity + !! tendency [g kg-1 T-1 -> g kg-1 s-1] + real, dimension(SZI_(G),SZJ_(G),SZK_(CS%Grid)) :: T !< The updated temperature [degC] + real, dimension(SZI_(G),SZJ_(G),SZK_(CS%Grid)) :: S !< The updated salinity [g kg-1] + real :: missing_value + + if (.not. associated(CS)) return + if (CS%assim_method .eq. NO_ASSIM .and. (.not. CS%do_bias_adjustment)) return + + call cpu_clock_begin(id_clock_apply_increments) + + T_inc(:,:,:) = 0.0; S_inc(:,:,:) = 0.0; T(:,:,:) = 0.0; S(:,:,:) = 0.0 + if (CS%assim_method > 0 ) then + T = T + CS%tv%T + S = S + CS%tv%S + endif + if (CS%do_bias_adjustment ) then + T = T + CS%tv_bc%T + S = S + CS%tv_bc%S + endif + + isc=G%isc; iec=G%iec; jsc=G%jsc; jec=G%jec + do j=jsc,jec; do i=isc,iec + call remapping_core_h(CS%remapCS, CS%nk, CS%h(i,j,:), T(i,j,:), & + G%ke, h(i,j,:), T_inc(i,j,:)) + call remapping_core_h(CS%remapCS, CS%nk, CS%h(i,j,:), S(i,j,:), & + G%ke, h(i,j,:), S_inc(i,j,:)) + enddo; enddo + + + call pass_var(T_inc, G%Domain) + call pass_var(S_inc, G%Domain) + + tv%T(isc:iec,jsc:jec,:)=tv%T(isc:iec,jsc:jec,:)+T_inc(isc:iec,jsc:jec,:)*dt + tv%S(isc:iec,jsc:jec,:)=tv%S(isc:iec,jsc:jec,:)+S_inc(isc:iec,jsc:jec,:)*dt + + call pass_var(tv%T, G%Domain) + call pass_var(tv%S, G%Domain) + + call enable_averaging(dt, Time_end, CS%diag_CS) + if (CS%id_inc_t > 0) call post_data(CS%id_inc_t, T_inc, CS%diag_CS) + if (CS%id_inc_s > 0) call post_data(CS%id_inc_s, S_inc, CS%diag_CS) + call disable_averaging(CS%diag_CS) + + call diag_update_remap_grids(CS%diag_CS) + call cpu_clock_end(id_clock_apply_increments) + + return end subroutine apply_oda_tracer_increments + subroutine set_up_global_tgrid(T_grid, CS, G) + type(grid_type), pointer :: T_grid !< global tracer grid + type(ODA_CS), pointer, intent(in) :: CS + type(ocean_grid_type), pointer :: G !< domain and grid information for ocean model + + ! local variables + real, dimension(:,:), allocatable :: global2D, global2D_old + integer :: i, j, k + + ! get global grid information from ocean_model + if (associated(T_grid)) call MOM_error(FATAL,'MOM_oda_driver:set_up_global_tgrid called with associated T_grid') + + allocate(T_grid) + allocate(T_grid%x(CS%ni,CS%nj)) + allocate(T_grid%y(CS%ni,CS%nj)) + if (CS%do_bias_adjustment) then + allocate(T_grid%basin_mask(CS%ni,CS%nj)) + endif + allocate(T_grid%bathyT(CS%ni,CS%nj)) + call global_field(CS%mpp_domain, CS%Grid%geolonT, T_grid%x) + call global_field(CS%mpp_domain, CS%Grid%geolatT, T_grid%y) + call global_field(CS%mpp_domain, CS%oda_grid%basin_mask, T_grid%basin_mask) + call global_field(CS%domains(CS%ensemble_id)%mpp_domain, G%bathyT, T_grid%bathyT) + + allocate(T_grid%mask(CS%ni,CS%nj,CS%nk)) + allocate(T_grid%z(CS%ni,CS%nj,CS%nk)) + allocate(global2D(CS%ni,CS%nj)) + allocate(global2D_old(CS%ni,CS%nj)) + T_grid%mask(:,:,:) = 0.0 + T_grid%z(:,:,:) = 0.0 + + do k = 1, CS%nk + call global_field(G%Domain%mpp_domain, CS%h(:,:,k), global2D) + do i=1,CS%ni ; do j=1,CS%nj + if ( global2D(i,j) > 1 ) then + T_grid%mask(i,j,k) = 1.0 + endif + enddo; enddo + if (k == 1) then + T_grid%z(:,:,k) = global2D/2 + else + T_grid%z(:,:,k) = T_grid%z(:,:,k-1) + (global2D + global2D_old)/2 + endif + global2D_old = global2D + enddo + + deallocate(global2D) + deallocate(global2D_old) + end subroutine set_up_global_tgrid + !> \namespace MOM_oda_driver_mod !! !! \section section_ODA The Ocean data assimilation (DA) and Ensemble Framework