Skip to content

Commit

Permalink
Merge pull request mom-ocean#1454 from jiandewang/EMC-IAU-candidate-2…
Browse files Browse the repository at this point in the history 
…0210727

IAU candidate 20210727
  • Loading branch information
@marshallward
marshallward authored Aug 6, 2021
2 parents b629e8d + 3fb3a4c commit e794c41
Show file tree
Hide file tree
Showing 4 changed files with 1,099 additions and 10 deletions.
19 changes: 14 additions & 5 deletions src/core/MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,6 +139,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, init_oda_incupd_diags

! 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
Expand Down Expand Up @@ -369,6 +371,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
Expand Down Expand Up @@ -1667,6 +1671,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"
Expand Down Expand Up @@ -2421,12 +2426,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
Expand Down Expand Up @@ -2468,12 +2473,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

Expand Down Expand Up @@ -2665,7 +2670,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)) &
Expand All @@ -2674,6 +2679,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)
Expand Down
178 changes: 175 additions & 3 deletions src/initialization/MOM_state_initialization.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -94,6 +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_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: calc_oda_increments, output_oda_incupd_inc

implicit none ; private

Expand All @@ -114,7 +117,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
Expand All @@ -140,6 +143,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
Expand All @@ -157,7 +161,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.
Expand Down Expand Up @@ -478,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.

Expand Down Expand Up @@ -615,6 +629,13 @@ 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()')

! 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, restart_CS, Time)
endif


end subroutine MOM_initialize_state

!> Reads the layer thicknesses or interface heights from a file.
Expand Down Expand Up @@ -2022,6 +2043,157 @@ 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, 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(inout) :: h !< Layer thickness [H ~> m or kg m-2] (in)

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
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)
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

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.)
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")
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_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)
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,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(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, restart_CS)
deallocate(hoda)

! 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, 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, oda_incupd_CSp)
deallocate(tmp_tr)
endif

! 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)
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, oda_incupd_CSp)
deallocate(tmp_u,tmp_v)
endif

! 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
if (is_root_pe()) call MOM_error(NOTE,"incupd using full fields ")
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)
Expand Down
Loading

0 comments on commit e794c41

Please sign in to comment.