From 67c6668a310ef0d56bf071924cb04ff75a5750be Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 24 Aug 2017 09:19:22 -0600 Subject: [PATCH 01/20] Fixes inequality sign to avoid division by zero --- src/equation_of_state/MOM_EOS.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/equation_of_state/MOM_EOS.F90 b/src/equation_of_state/MOM_EOS.F90 index d3d6d7dc33..d4604b42e3 100644 --- a/src/equation_of_state/MOM_EOS.F90 +++ b/src/equation_of_state/MOM_EOS.F90 @@ -1393,12 +1393,12 @@ subroutine find_depth_of_pressure_in_cell(T_t, T_b, S_t, S_b, z_t, z_b, P_t, P_t P_b = P_t + dp ! Anomalous pressure at bottom of cell - if (P_tgt < P_t ) then + if (P_tgt <= P_t ) then z_out = z_t return endif - if (P_tgt > P_b) then + if (P_tgt >= P_b) then z_out = z_b return endif From 44bd09b40606b7e35d579de6356412e7cb0b4e97 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 24 Aug 2017 09:21:46 -0600 Subject: [PATCH 02/20] Corrects axis types for outputs in the ice_shelf module --- src/ice_shelf/MOM_ice_shelf.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index fe45471dbd..985ab33a9f 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -1905,9 +1905,9 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti 'salinity in the boundary layer minus salinity at the ice-ocean interface.', 'PPT') CS%id_Sbdry = register_diag_field('ocean_model', 'sbdry', CS%diag%axesT1, CS%Time, & 'salinity at the ice-ocean interface.', 'PPT') - CS%id_u_ml = register_diag_field('ocean_model', 'u_ml', CS%diag%axesT1, CS%Time, & + CS%id_u_ml = register_diag_field('ocean_model', 'u_ml', CS%diag%axesCu1, CS%Time, & 'Eastward vel. in the boundary layer (used to compute ustar)', 'meter second-1') - CS%id_v_ml = register_diag_field('ocean_model', 'v_ml', CS%diag%axesT1, CS%Time, & + CS%id_v_ml = register_diag_field('ocean_model', 'v_ml', CS%diag%axesCv1, CS%Time, & 'Northward vel. in the boundary layer (used to compute ustar)', 'meter second-1') CS%id_exch_vel_s = register_diag_field('ocean_model', 'exch_vel_s', CS%diag%axesT1, CS%Time, & 'Sub-shelf salinity exchange velocity', 'meter second-1') @@ -1921,13 +1921,13 @@ subroutine initialize_ice_shelf(param_file, ocn_grid, Time, CS, diag, fluxes, Ti 'Fric vel under shelf', 'm/s') if (CS%shelf_mass_is_dynamic .and. .not.CS%override_shelf_movement) then - CS%id_u_shelf = register_diag_field('ocean_model','u_shelf',CS%diag%axesB1,CS%Time, & + CS%id_u_shelf = register_diag_field('ocean_model','u_shelf',CS%diag%axesCu1,CS%Time, & 'x-velocity of ice', 'm year') - CS%id_v_shelf = register_diag_field('ocean_model','v_shelf',CS%diag%axesB1,CS%Time, & + CS%id_v_shelf = register_diag_field('ocean_model','v_shelf',CS%diag%axesCv1,CS%Time, & 'y-velocity of ice', 'm year') - CS%id_u_mask = register_diag_field('ocean_model','u_mask',CS%diag%axesB1,CS%Time, & + CS%id_u_mask = register_diag_field('ocean_model','u_mask',CS%diag%axesCu1,CS%Time, & 'mask for u-nodes', 'none') - CS%id_v_mask = register_diag_field('ocean_model','v_mask',CS%diag%axesB1,CS%Time, & + CS%id_v_mask = register_diag_field('ocean_model','v_mask',CS%diag%axesCv1,CS%Time, & 'mask for v-nodes', 'none') CS%id_h_mask = register_diag_field('ocean_model','h_mask',CS%diag%axesT1,CS%Time, & 'ice shelf thickness', 'none') From a96c4ddc37155ca613b2104a62dab8e7ca53bcac Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 22 Aug 2017 16:36:31 -0400 Subject: [PATCH 03/20] Doxygenized MOM_get_input.F90 - Completed doxygenization of MOM_get_input.F90 - No answer changes. --- src/framework/MOM_get_input.F90 | 53 ++++++++++++++++----------------- 1 file changed, 25 insertions(+), 28 deletions(-) diff --git a/src/framework/MOM_get_input.F90 b/src/framework/MOM_get_input.F90 index 1aed45038a..08bdd8623d 100644 --- a/src/framework/MOM_get_input.F90 +++ b/src/framework/MOM_get_input.F90 @@ -1,17 +1,12 @@ +!> \brief Reads the only Fortran name list needed to boot-strap the model. +!! +!! The name list parameters indicate which directories to use for +!! certain types of input and output, and which files to look in for +!! the full parsable input parameter file(s). module MOM_get_input ! This file is part of MOM6. See LICENSE.md for the license. -!********+*********+*********+*********+*********+*********+*********+** -!* * -!* By Robert Hallberg, April 2013 * -!* * -!* The subroutine in this file reads the MOM6 namelist input, which * -!* indicates which directories to use for certain types of input and * -!* output, and where to look for the full parsable input file(s). * -!* * -!********+*********+*********+*********+*********+*********+*********+** - use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe use MOM_file_parser, only : open_param_file, param_file_type use MOM_io, only : file_exists, close_file, slasher, ensembler @@ -19,30 +14,28 @@ module MOM_get_input implicit none ; private -public Get_MOM_Input - -! This structure is to simplify communication with the calling code. +public get_MOM_input +!> Container for paths and parameter file names. type, public :: directories character(len=240) :: & - restart_input_dir = ' ',& ! The directory to read restart and input files. - restart_output_dir = ' ',&! The directory into which to write restart files. - output_directory = ' ', & ! The directory to use to write the model output. - input_filename = ' ' ! A string that indicates the input files or how - ! the run segment should be started. + restart_input_dir = ' ',& !< The directory to read restart and input files. + restart_output_dir = ' ',&!< The directory into which to write restart files. + output_directory = ' ', & !< The directory to use to write the model output. + input_filename = ' ' !< A string that indicates the input files or how + !! the run segment should be started. end type directories contains -subroutine Get_MOM_Input(param_file, dirs, check_params) - type(param_file_type), optional, intent(out) :: param_file !< A structure to parse for run-time parameters - type(directories), optional, intent(out) :: dirs - logical, optional, intent(in) :: check_params - -! See if the run is to be started from saved conditions, and get ! -! the names of the I/O directories and initialization file. This ! -! subroutine also calls the subroutine that allows run-time changes ! -! in parameters. ! +!> Get the names of the I/O directories and initialization file. +!! Also calls the subroutine that opens run-time parameter files. +subroutine get_MOM_input(param_file, dirs, check_params) + type(param_file_type), optional, intent(out) :: param_file !< A structure to parse for run-time parameters. + type(directories), optional, intent(out) :: dirs !< Container for paths and parameter file names. + logical, optional, intent(in) :: check_params !< If present and False will stop error checking for + !! run-time parameters. + ! Local variables integer, parameter :: npf = 5 ! Maximum number of parameter files character(len=240) :: & parameter_filename(npf) = ' ', & ! List of files containing parameters. @@ -57,18 +50,21 @@ subroutine Get_MOM_Input(param_file, dirs, check_params) namelist /MOM_input_nml/ output_directory, input_filename, parameter_filename, & restart_input_dir, restart_output_dir + ! Open namelist if (file_exists('input.nml')) then unit = open_namelist_file(file='input.nml') else call MOM_error(FATAL,'Required namelist file input.nml does not exist.') endif + ! Read namelist parameters ierr=1 ; do while (ierr /= 0) read(unit, nml=MOM_input_nml, iostat=io, end=10) ierr = check_nml_error(io, 'MOM_input_nml') enddo 10 call close_file(unit) + ! Store parameters in container if (present(dirs)) then dirs%output_directory = slasher(ensembler(output_directory)) dirs%restart_output_dir = slasher(ensembler(restart_output_dir)) @@ -76,6 +72,7 @@ subroutine Get_MOM_Input(param_file, dirs, check_params) dirs%input_filename = ensembler(input_filename) endif + ! Open run-time parameter file(s) if (present(param_file)) then output_dir = slasher(ensembler(output_directory)) valid_param_files = 0 @@ -90,6 +87,6 @@ subroutine Get_MOM_Input(param_file, dirs, check_params) "least 1 valid entry in input_filename in MOM_input_nml in input.nml.") endif -end subroutine Get_MOM_Input +end subroutine get_MOM_input end module MOM_get_input From 07d0f36fc9a3a63f03fe5a393db0cf217548a2e6 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 22 Aug 2017 17:28:10 -0400 Subject: [PATCH 04/20] Fix initialization of namelist parameters - Explicitly initialize variables before reading from namelist. - Initialization of variables in a declaration statement only does what you expect the first time the subroutine/function is called. Subsequent calls do not initialize the variables. - get_MOM_input() is called several times in some configurations. --- src/framework/MOM_get_input.F90 | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/src/framework/MOM_get_input.F90 b/src/framework/MOM_get_input.F90 index 08bdd8623d..c6c709a704 100644 --- a/src/framework/MOM_get_input.F90 +++ b/src/framework/MOM_get_input.F90 @@ -38,18 +38,25 @@ subroutine get_MOM_input(param_file, dirs, check_params) ! Local variables integer, parameter :: npf = 5 ! Maximum number of parameter files character(len=240) :: & - parameter_filename(npf) = ' ', & ! List of files containing parameters. - output_directory = ' ', & ! Directory to use to write the model output. - restart_input_dir = ' ', & ! Directory for reading restart and input files. - restart_output_dir = ' ', & ! Directory into which to write restart files. - input_filename = ' ' ! A string that indicates the input files or how - ! the run segment should be started. + parameter_filename(npf), & ! List of files containing parameters. + output_directory, & ! Directory to use to write the model output. + restart_input_dir, & ! Directory for reading restart and input files. + restart_output_dir, & ! Directory into which to write restart files. + input_filename ! A string that indicates the input files or how + ! the run segment should be started. character(len=240) :: output_dir integer :: unit, io, ierr, valid_param_files namelist /MOM_input_nml/ output_directory, input_filename, parameter_filename, & restart_input_dir, restart_output_dir + ! Default values in case parameter is not set in file input.nml + parameter_filename(:) = ' ' + output_directory = ' ' + restart_input_dir = ' ' + restart_output_dir = ' ' + input_filename = ' ' + ! Open namelist if (file_exists('input.nml')) then unit = open_namelist_file(file='input.nml') From aa6280893138f5a1f9a6b9316eaa188cfec7b711 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Thu, 24 Aug 2017 11:17:49 -0400 Subject: [PATCH 05/20] Allows coupled driver to specify restart file - Some of the changes in this commit are temporary; once the mct_driver is no longer using ocean_model_MOM.F90, some of the optional arguments can be removed. - NCAR's MCT driver specifics the absolute path of a restart file with a generated name (dates). This commit allows the driver to override the namelist parameter input_filename. - To use, set restart_input_dir='/' in input.nml since leaving restart_input_dir blank automatically replaces it with './' - No answer changes. --- config_src/coupled_driver/ocean_model_MOM.F90 | 7 ++++--- src/core/MOM.F90 | 7 +++++-- src/framework/MOM_get_input.F90 | 6 +++++- src/framework/MOM_restart.F90 | 2 +- 4 files changed, 15 insertions(+), 7 deletions(-) diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index a30f7a7a81..8a4809b04e 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -28,7 +28,7 @@ module ocean_model_mod use MOM_error_handler, only : callTree_enter, callTree_leave use MOM_file_parser, only : get_param, log_version, close_param_file, param_file_type use MOM_forcing_type, only : forcing, forcing_diagnostics, mech_forcing_diags, forcing_accumulate -use MOM_get_input, only : Get_MOM_Input, directories +use MOM_get_input, only : directories use MOM_grid, only : ocean_grid_type use MOM_io, only : close_file, file_exists, read_data, write_version_number use MOM_restart, only : save_restart @@ -189,7 +189,7 @@ module ocean_model_mod !> ocean_model_init initializes the ocean model, including registering fields !! for restarts and reading restart files if appropriate. -subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn) +subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, input_restart_file) type(ocean_public_type), target, & intent(inout) :: Ocean_sfc !< A structure containing various !! publicly visible ocean surface properties after initialization, @@ -205,6 +205,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn) !! in the calculation of additional gas or other !! tracer fluxes, and can be used to spawn related !! internal variables in the ice model. + character(len=*), optional, intent(in) :: input_restart_file !< If present, name of restart file to read ! This subroutine initializes both the ocean state and the ocean surface type. ! Because of the way that indicies and domains are handled, Ocean_sfc must have @@ -241,7 +242,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn) OS%Time = Time_in call initialize_MOM(OS%Time, param_file, OS%dirs, OS%MOM_CSp, Time_in, & - offline_tracer_mode=offline_tracer_mode) + offline_tracer_mode=offline_tracer_mode, input_restart_file=input_restart_file) OS%grid => OS%MOM_CSp%G ; OS%GV => OS%MOM_CSp%GV OS%C_p = OS%MOM_CSp%tv%C_p OS%fluxes%C_p = OS%MOM_CSp%tv%C_p diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index d192f0faad..85581bfe1b 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1475,7 +1475,7 @@ subroutine step_offline(fluxes, state, Time_start, time_interval, CS) end subroutine step_offline !> This subroutine initializes MOM. -subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mode) +subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mode, input_restart_file) type(time_type), target, intent(inout) :: Time !< model time, set in this routine type(param_file_type), intent(out) :: param_file !< structure indicating paramater file to parse type(directories), intent(out) :: dirs !< structure with directory paths @@ -1483,6 +1483,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo type(time_type), optional, intent(in) :: Time_in !< time passed to MOM_initialize_state when !! model is not being started from a restart file logical, optional, intent(out) :: offline_tracer_mode !< True if tracers are being run offline + character(len=*),optional, intent(in) :: input_restart_file !< If present, name of restart file to read ! local type(ocean_grid_type), pointer :: G => NULL() ! A pointer to a structure with metrics and related @@ -1552,7 +1553,9 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in, offline_tracer_mo Start_time = Time ; if (present(Time_in)) Start_time = Time_in - call Get_MOM_Input(param_file, dirs) + ! Read paths and filenames from namelist and store in "dirs". + ! Also open the parsed input parameter file(s) and setup param_file. + call get_MOM_input(param_file, dirs, default_input_filename=input_restart_file) verbosity = 2 ; call read_param(param_file, "VERBOSITY", verbosity) call MOM_set_verbosity(verbosity) diff --git a/src/framework/MOM_get_input.F90 b/src/framework/MOM_get_input.F90 index c6c709a704..2ee3e93bbd 100644 --- a/src/framework/MOM_get_input.F90 +++ b/src/framework/MOM_get_input.F90 @@ -30,11 +30,14 @@ module MOM_get_input !> Get the names of the I/O directories and initialization file. !! Also calls the subroutine that opens run-time parameter files. -subroutine get_MOM_input(param_file, dirs, check_params) +subroutine get_MOM_input(param_file, dirs, check_params, default_input_filename) type(param_file_type), optional, intent(out) :: param_file !< A structure to parse for run-time parameters. type(directories), optional, intent(out) :: dirs !< Container for paths and parameter file names. logical, optional, intent(in) :: check_params !< If present and False will stop error checking for !! run-time parameters. + character(len=*), optional, intent(in) :: default_input_filename !< If present, is the value assumed for + !! input_filename if input_filename is not listed + !! in the namelist MOM_input_nml. ! Local variables integer, parameter :: npf = 5 ! Maximum number of parameter files character(len=240) :: & @@ -56,6 +59,7 @@ subroutine get_MOM_input(param_file, dirs, check_params) restart_input_dir = ' ' restart_output_dir = ' ' input_filename = ' ' + if (present(default_input_filename)) input_filename = trim(default_input_filename) ! Open namelist if (file_exists('input.nml')) then diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index d57c894545..f538d3f3e4 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -889,7 +889,7 @@ subroutine restore_state(filename, directory, day, G, CS) ! restart_init. character(len=200) :: filepath ! The path (dir/file) to the file being opened. - character(len=80) :: fname ! The name of the current file. + character(len=200) :: fname ! The name of the current file. character(len=8) :: suffix ! A suffix (like "_2") that is added to any ! additional restart files. character(len=256) :: mesg ! A message for warnings. From 09fa6c307781b94ed4b8879482ea8e387d0d6ccc Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 28 Aug 2017 11:20:32 -0600 Subject: [PATCH 06/20] Deletes get_state_pointers from ocean_model_MOM, its no longer needed there --- config_src/coupled_driver/ocean_model_MOM.F90 | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 8a4809b04e..4495641fb7 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -1080,15 +1080,4 @@ subroutine ocean_public_type_chksum(id, timestep, ocn) end subroutine ocean_public_type_chksum -!> Returns pointers to objects within ocean_state_type -subroutine get_state_pointers(OS, grid, surf) - type(ocean_state_type), pointer :: OS !< Ocean state type - type(ocean_grid_type), optional, pointer :: grid !< Ocean grid - type(surface), optional, pointer :: surf !< Ocean surface state - - if (present(grid)) grid => OS%grid - if (present(surf)) surf=> OS%state - -end subroutine get_state_pointers - end module ocean_model_mod From d043f87ad3f0b3b98055f2d33de03f6732317aa4 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 28 Aug 2017 13:24:34 -0600 Subject: [PATCH 07/20] Deletes public declaration for get_state_pointers --- config_src/coupled_driver/ocean_model_MOM.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 4495641fb7..341fdfcf7c 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -76,7 +76,6 @@ module ocean_model_mod public ice_ocn_bnd_type_chksum public ocean_public_type_chksum public ocean_model_data_get -public get_state_pointers interface ocean_model_data_get module procedure ocean_model_data1D_get From 69f5ac70d77d7a03277bb19cc0b09aa87a170d1d Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 28 Aug 2017 14:23:06 -0600 Subject: [PATCH 08/20] Separates ocn_comp_mct from ocean_model_mod This PR is the first attempt to separates ocn_comp_mct from ocean_model_mod. --- config_src/mct_driver/ocn_comp_mct.F90 | 712 ++++++++++++++++++++++++- 1 file changed, 690 insertions(+), 22 deletions(-) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index 847c3daba9..b750264589 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -14,7 +14,7 @@ module ocn_comp_mct mct_aVect_nRattr use mct_mod, only: mct_gGrid, mct_gGrid_init, mct_gGrid_importRAttr, & mct_gGrid_importIAttr -use mct_mod, only : mct_avect_indexra, mct_aVect_clean +use mct_mod, only: mct_avect_indexra, mct_aVect_clean use seq_flds_mod, only: seq_flds_x2o_fields, seq_flds_o2x_fields, seq_flds_dom_coord, & seq_flds_dom_other use seq_infodata_mod, only: seq_infodata_type, seq_infodata_GetData, & @@ -26,17 +26,44 @@ module ocn_comp_mct use shr_kind_mod, only: shr_kind_r8 ! MOM6 modules -use ocean_model_mod, only: ocean_state_type, ocean_public_type, ocean_model_init_sfc -use ocean_model_mod, only: ocean_model_init, get_state_pointers, ocean_model_restart -use ocean_model_mod, only: ice_ocean_boundary_type, update_ocean_model -use MOM_domains, only: MOM_infra_init, num_pes, root_pe, pe_here -use MOM_domains, only : pass_var, AGRID -use MOM_grid, only: ocean_grid_type, get_global_grid_size -use MOM_variables, only: surface -use MOM_error_handler, only: MOM_error, FATAL, is_root_pe -use MOM_time_manager, only: time_type, set_date, set_time, set_calendar_type, NOLEAP -use MOM_file_parser, only: get_param, log_version, param_file_type -use MOM_get_input, only: Get_MOM_Input, directories +use MOM, only: initialize_MOM, step_MOM, MOM_control_struct, MOM_end +use MOM, only: calculate_surface_state, allocate_surface_state +use MOM, only: finish_MOM_initialization, step_offline +use MOM_surface_forcing, only: ice_ocean_boundary_type, surface_forcing_CS +use MOM_forcing_type, only: forcing, forcing_diagnostics, mech_forcing_diags, forcing_accumulate +use MOM_forcing_type, only: allocate_forcing_type +use MOM_surface_forcing, only: surface_forcing_init, convert_IOB_to_fluxes +use MOM_restart, only: save_restart +use MOM_domains, only: MOM_infra_init, num_pes, root_pe, pe_here +use MOM_domains, only: pass_vector, BGRID_NE, CGRID_NE +use MOM_domains, only: pass_var, AGRID +use MOM_grid, only: ocean_grid_type, get_global_grid_size +use MOM_verticalGrid, only: verticalGrid_type +use MOM_variables, only: surface +use MOM_error_handler, only: MOM_error, FATAL, is_root_pe, WARNING +use MOM_error_handler, only: callTree_enter, callTree_leave +use MOM_time_manager, only: time_type, set_date, set_time, set_calendar_type, NOLEAP, get_date +use MOM_time_manager, only : operator(+), operator(-), operator(*), operator(/) +use MOM_time_manager, only : operator(/=), operator(>), get_time +use MOM_file_parser, only: get_param, log_version, param_file_type, close_param_file +use MOM_get_input, only: Get_MOM_Input, directories +use MOM_diag_mediator, only : diag_ctrl, enable_averaging, disable_averaging +use MOM_diag_mediator, only : diag_mediator_close_registration, diag_mediator_end +use MOM_ice_shelf, only: initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS +use MOM_sum_output, only: MOM_sum_output_init, sum_output_CS +use MOM_sum_output, only: write_energy, accumulate_net_input +use MOM_string_functions, only : uppercase +use MOM_constants, only : CELSIUS_KELVIN_OFFSET, hlf +use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct + +! FMS +use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain +use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain + +! GFDL coupler +use coupler_types_mod, only : coupler_1d_bc_type, coupler_2d_bc_type +use coupler_types_mod, only : coupler_type_spawn +use coupler_types_mod, only : coupler_type_initialized, coupler_type_copy_data ! By default make data private implicit none; private @@ -118,6 +145,99 @@ module ocn_comp_mct end type cpl_indices + !> This type is used for communication with other components via the FMS coupler. + ! The element names and types can be changed only with great deliberation, hence + ! the persistnce of things like the cutsy element name "avg_kount". + type, public :: ocean_public_type + type(domain2d) :: Domain ! The domain for the surface fields. + logical :: is_ocean_pe ! .true. on processors that run the ocean model. + character(len=32) :: instance_name = '' ! A name that can be used to identify + ! this instance of an ocean model, for example + ! in ensembles when writing messages. + integer, pointer, dimension(:) :: pelist => NULL() ! The list of ocean PEs. + logical, pointer, dimension(:,:) :: maskmap =>NULL() ! A pointer to an array + ! indicating which logical processors are actually used for + ! the ocean code. The other logical processors would be all + ! land points and are not assigned to actual processors. + ! This need not be assigned if all logical processors are used. + + integer :: stagger = -999 ! The staggering relative to the tracer points + ! points of the two velocity components. Valid entries + ! include AGRID, BGRID_NE, CGRID_NE, BGRID_SW, and CGRID_SW, + ! corresponding to the community-standard Arakawa notation. + ! (These are named integers taken from mpp_parameter_mod.) + ! Following MOM, this is BGRID_NE by default when the ocean + ! is initialized, but here it is set to -999 so that a + ! global max across ocean and non-ocean processors can be + ! used to determine its value. + real, pointer, dimension(:,:) :: & + t_surf => NULL(), & ! SST on t-cell (degrees Kelvin) + s_surf => NULL(), & ! SSS on t-cell (psu) + u_surf => NULL(), & ! i-velocity at the locations indicated by stagger, m/s. + v_surf => NULL(), & ! j-velocity at the locations indicated by stagger, m/s. + sea_lev => NULL(), & ! Sea level in m after correction for surface pressure, + ! i.e. dzt(1) + eta_t + patm/rho0/grav (m) + frazil =>NULL(), & ! Accumulated heating (in Joules/m^2) from frazil + ! formation in the ocean. + area => NULL() ! cell area of the ocean surface, in m2. + type(coupler_2d_bc_type) :: fields ! A structure that may contain an + ! array of named tracer-related fields. + integer :: avg_kount ! Used for accumulating averages of this type. + integer, dimension(2) :: axes = 0 ! Axis numbers that are available + ! for I/O using this surface data. + end type ocean_public_type + + !> Contains information about the ocean state, although it is not necessary that + !! this is implemented with all models. This type is private, and can therefore vary + !! between different ocean models. + type, public :: ocean_state_type ; private + logical :: is_ocean_PE = .false. !< True if this is an ocean PE. + type(time_type) :: Time !< The ocean model's time and master clock. + integer :: Restart_control !< An integer that is bit-tested to determine whether + !! incremental restart files are saved and whether they + !! have a time stamped name. +1 (bit 0) for generic + !! files and +2 (bit 1) for time-stamped files. A + !! restart file is saved at the end of a run segment + !! unless Restart_control is negative. + type(time_type) :: energysavedays ! The interval between writing the energies + ! and other integral quantities of the run. + type(time_type) :: write_energy_time ! The next time to write to the energy file. + + integer :: nstep = 0 ! The number of calls to update_ocean. + logical :: use_ice_shelf ! If true, the ice shelf model is enabled. + logical :: icebergs_apply_rigid_boundary ! If true, the icebergs can change ocean bd condition. + real :: kv_iceberg ! The viscosity of the icebergs in m2/s (for ice rigidity) + real :: berg_area_threshold ! Fraction of grid cell which iceberg must occupy + !so that fluxes below are set to zero. (0.5 is a + !good value to use. Not applied for negative values. + real :: latent_heat_fusion ! Latent heat of fusion + real :: density_iceberg ! A typical density of icebergs in kg/m3 (for ice rigidity) + type(ice_shelf_CS), pointer :: Ice_shelf_CSp => NULL() + logical :: restore_salinity ! If true, the coupled MOM driver adds a term to + ! restore salinity to a specified value. + logical :: restore_temp ! If true, the coupled MOM driver adds a term to + ! restore sst to a specified value. + real :: press_to_z ! A conversion factor between pressure and ocean + ! depth in m, usually 1/(rho_0*g), in m Pa-1. + real :: C_p ! The heat capacity of seawater, in J K-1 kg-1. + + type(directories) :: dirs ! A structure containing several relevant directory paths. + type(forcing) :: fluxes ! A structure containing pointers to + ! the ocean forcing fields. + type(forcing) :: flux_tmp ! A secondary structure containing pointers to the + ! ocean forcing fields for when multiple coupled + ! timesteps are taken per thermodynamic step. + type(surface) :: state ! A structure containing pointers to + ! the ocean surface state fields. + type(ocean_grid_type), pointer :: grid => NULL() ! A pointer to a grid structure + ! containing metrics and related information. + type(verticalGrid_type), pointer :: GV => NULL() ! A pointer to a vertical grid + ! structure containing metrics and related information. + type(MOM_control_struct), pointer :: MOM_CSp => NULL() + type(surface_forcing_CS), pointer :: forcing_CSp => NULL() + type(sum_output_CS), pointer :: sum_output_CSp => NULL() + end type ocean_state_type + !> Control structure for this module type MCT_MOM_Data @@ -156,6 +276,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) integer :: year, month, day, hour, minute, seconds, seconds_n, seconds_d, rc character(len=384) :: runid !< Run ID character(len=384) :: runtype !< Run type + character(len=384) :: restartfile !< Path/Name of restart file character(len=32) :: starttype !< infodata start type integer :: mpicom_ocn !< MPI ocn communicator integer :: npes, pe0 !< # of processors and current processor @@ -303,8 +424,22 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) glb%c1 = 0.0; glb%c2 = 0.0; glb%c3 = 0.0; glb%c4 = 0.0 endif + restartfile = '/glade/scratch/gmarques/mom_test/run/RESTART/MOM.res_Y0001_D002_S00000.nc' ! Initialize the MOM6 model - call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in) + if (runtype == "initial") then ! startup (new run) - 'n' is needed below since we don't + ! specify input_filename in input.nml + call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in, input_restart_file = 'n') + else ! hybrid or branch or continuos runs + !call seq_infodata_GetData(glb%infodata, restart_file=restartfile) + call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in, input_restart_file = restartfile) + endif + + if (debug .and. root_pe().eq.pe_here()) then + write(6,*)'runtype, restartfile',runtype, restartfile + endif + + ! Initialize the MOM6 model + !call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in) ! Initialize ocn_state%state out of sight call ocean_model_init_sfc(glb%ocn_state, glb%ocn_public) @@ -528,6 +663,362 @@ subroutine coupler_indices_init(ind) end subroutine coupler_indices_init +!> Initializes the ocean model, including registering fields +!! for restarts and reading restart files if appropriate. +subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, input_restart_file) + type(ocean_public_type), target, & + intent(inout) :: Ocean_sfc !< A structure containing various + !! publicly visible ocean surface properties after initialization, + !! the data in this type is intent(out). + type(ocean_state_type), pointer :: OS !< A structure whose internal + !! contents are private to ocean_model_mod that may be used to + !! contain all information about the ocean's interior state. + type(time_type), intent(in) :: Time_init !< The start time for the coupled model's calendar + type(time_type), intent(in) :: Time_in !< The time at which to initialize the ocean model. + type(coupler_1d_bc_type), & + optional, intent(in) :: gas_fields_ocn !< If present, this type describes the + !! ocean and surface-ice fields that will participate + !! in the calculation of additional gas or other + !! tracer fluxes, and can be used to spawn related + !! internal variables in the ice model. + character(len=*), optional, intent(in) :: input_restart_file !< If present, name of restart file to read + +! This subroutine initializes both the ocean state and the ocean surface type. +! Because of the way that indicies and domains are handled, Ocean_sfc must have +! been used in a previous call to initialize_ocean_type. + +! Arguments: Ocean_sfc - A structure containing various publicly visible ocean +! surface properties after initialization, this is intent(out). +! (out,private) OS - A structure whose internal contents are private +! to ocean_model_mod that may be used to contain all +! information about the ocean's interior state. +! (in) Time_init - The start time for the coupled model's calendar. +! (in) Time_in - The time at which to initialize the ocean model. + real :: Time_unit ! The time unit in seconds for ENERGYSAVEDAYS. + real :: Rho0 ! The Boussinesq ocean density, in kg m-3. + real :: G_Earth ! The gravitational acceleration in m s-2. +! This include declares and sets the variable "version". +#include "version_variable.h" + character(len=40) :: mdl = "ocean_model_init" ! This module's name. + character(len=48) :: stagger + integer :: secs, days + type(param_file_type) :: param_file !< A structure to parse for run-time parameters + logical :: offline_tracer_mode + + call callTree_enter("ocean_model_init(), ocn_comp_mct.F90") + if (associated(OS)) then + call MOM_error(WARNING, "ocean_model_init called with an associated "// & + "ocean_state_type structure. Model is already initialized.") + return + endif + allocate(OS) + + OS%is_ocean_pe = Ocean_sfc%is_ocean_pe + if (.not.OS%is_ocean_pe) return + + OS%Time = Time_in + call initialize_MOM(OS%Time, param_file, OS%dirs, OS%MOM_CSp, Time_in, & + offline_tracer_mode=offline_tracer_mode, input_restart_file=input_restart_file) + OS%grid => OS%MOM_CSp%G ; OS%GV => OS%MOM_CSp%GV + OS%C_p = OS%MOM_CSp%tv%C_p + OS%fluxes%C_p = OS%MOM_CSp%tv%C_p + + ! Read all relevant parameters and write them to the model log. + call log_version(param_file, mdl, version, "") + call get_param(param_file, mdl, "RESTART_CONTROL", OS%Restart_control, & + "An integer whose bits encode which restart files are \n"//& + "written. Add 2 (bit 1) for a time-stamped file, and odd \n"//& + "(bit 0) for a non-time-stamped file. A restart file \n"//& + "will be saved at the end of the run segment for any \n"//& + "non-negative value.", default=1) + call get_param(param_file, mdl, "TIMEUNIT", Time_unit, & + "The time unit for ENERGYSAVEDAYS.", & + units="s", default=86400.0) + call get_param(param_file, mdl, "ENERGYSAVEDAYS",OS%energysavedays, & + "The interval in units of TIMEUNIT between saves of the \n"//& + "energies of the run and other globally summed diagnostics.", & + default=set_time(0,days=1), timeunit=Time_unit) + + call get_param(param_file, mdl, "OCEAN_SURFACE_STAGGER", stagger, & + "A case-insensitive character string to indicate the \n"//& + "staggering of the surface velocity field that is \n"//& + "returned to the coupler. Valid values include \n"//& + "'A', 'B', or 'C'.", default="C") + if (uppercase(stagger(1:1)) == 'A') then ; Ocean_sfc%stagger = AGRID + elseif (uppercase(stagger(1:1)) == 'B') then ; Ocean_sfc%stagger = BGRID_NE + elseif (uppercase(stagger(1:1)) == 'C') then ; Ocean_sfc%stagger = CGRID_NE + else ; call MOM_error(FATAL,"ocean_model_init: OCEAN_SURFACE_STAGGER = "// & + trim(stagger)//" is invalid.") ; endif + + call get_param(param_file, mdl, "RESTORE_SALINITY",OS%restore_salinity, & + "If true, the coupled driver will add a globally-balanced \n"//& + "fresh-water flux that drives sea-surface salinity \n"//& + "toward specified values.", default=.false.) + call get_param(param_file, mdl, "RESTORE_TEMPERATURE",OS%restore_temp, & + "If true, the coupled driver will add a \n"//& + "heat flux that drives sea-surface temperauture \n"//& + "toward specified values.", default=.false.) + call get_param(param_file, mdl, "RHO_0", Rho0, & + "The mean ocean density used with BOUSSINESQ true to \n"//& + "calculate accelerations and the mass for conservation \n"//& + "properties, or with BOUSSINSEQ false to convert some \n"//& + "parameters from vertical units of m to kg m-2.", & + units="kg m-3", default=1035.0) + call get_param(param_file, mdl, "G_EARTH", G_Earth, & + "The gravitational acceleration of the Earth.", & + units="m s-2", default = 9.80) + + call get_param(param_file, mdl, "ICE_SHELF", OS%use_ice_shelf, & + "If true, enables the ice shelf model.", default=.false.) + + call get_param(param_file, mdl, "ICEBERGS_APPLY_RIGID_BOUNDARY", OS%icebergs_apply_rigid_boundary, & + "If true, allows icebergs to change boundary condition felt by ocean", default=.false.) + + if (OS%icebergs_apply_rigid_boundary) then + call get_param(param_file, mdl, "KV_ICEBERG", OS%kv_iceberg, & + "The viscosity of the icebergs", units="m2 s-1",default=1.0e10) + call get_param(param_file, mdl, "DENSITY_ICEBERGS", OS%density_iceberg, & + "A typical density of icebergs.", units="kg m-3", default=917.0) + call get_param(param_file, mdl, "LATENT_HEAT_FUSION", OS%latent_heat_fusion, & + "The latent heat of fusion.", units="J/kg", default=hlf) + call get_param(param_file, mdl, "BERG_AREA_THRESHOLD", OS%berg_area_threshold, & + "Fraction of grid cell which iceberg must occupy, so that fluxes \n"//& + "below berg are set to zero. Not applied for negative \n"//& + " values.", units="non-dim", default=-1.0) + endif + + OS%press_to_z = 1.0/(Rho0*G_Earth) + + ! Consider using a run-time flag to determine whether to do the diagnostic + ! vertical integrals, since the related 3-d sums are not negligible in cost. + call allocate_surface_state(OS%state, OS%grid, OS%MOM_CSp%use_temperature, & + do_integrals=.true., gas_fields_ocn=gas_fields_ocn) + + call surface_forcing_init(Time_in, OS%grid, param_file, OS%MOM_CSp%diag, & + OS%forcing_CSp, OS%restore_salinity, OS%restore_temp) + + if (OS%use_ice_shelf) then + call initialize_ice_shelf(param_file, OS%grid, OS%Time, OS%ice_shelf_CSp, & + OS%MOM_CSp%diag, OS%fluxes) + endif + if (OS%icebergs_apply_rigid_boundary) then + !call allocate_forcing_type(OS%grid, OS%fluxes, iceberg=.true.) + !This assumes that the iceshelf and ocean are on the same grid. I hope this is true + if (.not. OS%use_ice_shelf) call allocate_forcing_type(OS%grid, OS%fluxes, ustar=.true., shelf=.true.) + endif + + call MOM_sum_output_init(OS%grid, param_file, OS%dirs%output_directory, & + OS%MOM_CSp%ntrunc, Time_init, OS%sum_output_CSp) + + ! This call has been moved into the first call to update_ocean_model. +! call write_energy(OS%MOM_CSp%u, OS%MOM_CSp%v, OS%MOM_CSp%h, OS%MOM_CSp%tv, & +! OS%Time, 0, OS%grid, OS%GV, OS%sum_output_CSp, OS%MOM_CSp%tracer_flow_CSp) + + ! write_energy_time is the next integral multiple of energysavedays. + OS%write_energy_time = Time_init + OS%energysavedays * & + (1 + (OS%Time - Time_init) / OS%energysavedays) + + if (ASSOCIATED(OS%grid%Domain%maskmap)) then + call initialize_ocean_public_type(OS%grid%Domain%mpp_domain, Ocean_sfc, & + OS%MOM_CSp%diag, maskmap=OS%grid%Domain%maskmap, & + gas_fields_ocn=gas_fields_ocn) + else + call initialize_ocean_public_type(OS%grid%Domain%mpp_domain, Ocean_sfc, & + OS%MOM_CSp%diag, gas_fields_ocn=gas_fields_ocn) + endif + + ! This call can only occur here if the coupler_bc_type variables have been + ! initialized already using the information from gas_fields_ocn. + if (present(gas_fields_ocn)) then + call calculate_surface_state(OS%state, OS%MOM_CSp%u, & + OS%MOM_CSp%v, OS%MOM_CSp%h, OS%MOM_CSp%ave_ssh,& + OS%grid, OS%GV, OS%MOM_CSp) + + call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid, & + OS%MOM_CSp%use_conT_absS) + endif + + call close_param_file(param_file) + call diag_mediator_close_registration(OS%MOM_CSp%diag) + + if (is_root_pe()) & + write(*,'(/12x,a/)') '======== COMPLETED MOM INITIALIZATION ========' + + call callTree_leave("ocean_model_init(") +end subroutine ocean_model_init + +!> This subroutine extracts the surface properties from the ocean's internal +!! state and stores them in the ocean type returned to the calling ice model. +!! It has to be separate from the ocean_initialization call because the coupler +!! module allocates the space for some of these variables. +subroutine ocean_model_init_sfc(OS, Ocean_sfc) + type(ocean_state_type), pointer :: OS + type(ocean_public_type), intent(inout) :: Ocean_sfc + + integer :: is, ie, js, je + + is = OS%grid%isc ; ie = OS%grid%iec ; js = OS%grid%jsc ; je = OS%grid%jec + call coupler_type_spawn(Ocean_sfc%fields, OS%state%tr_fields, & + (/is,is,ie,ie/), (/js,js,je,je/), as_needed=.true.) + + call calculate_surface_state(OS%state, OS%MOM_CSp%u, & + OS%MOM_CSp%v, OS%MOM_CSp%h, OS%MOM_CSp%ave_ssh,& + OS%grid, OS%GV, OS%MOM_CSp) + + call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid, & + OS%MOM_CSp%use_conT_absS) + +end subroutine ocean_model_init_sfc + +subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, maskmap, & + gas_fields_ocn) + type(domain2D), intent(in) :: input_domain + type(ocean_public_type), intent(inout) :: Ocean_sfc + type(diag_ctrl), intent(in) :: diag + logical, intent(in), optional :: maskmap(:,:) + type(coupler_1d_bc_type), & + optional, intent(in) :: gas_fields_ocn !< If present, this type describes the + !! ocean and surface-ice fields that will participate + !! in the calculation of additional gas or other + !! tracer fluxes. + + integer :: xsz, ysz, layout(2) + ! ice-ocean-boundary fields are always allocated using absolute indicies + ! and have no halos. + integer :: isc, iec, jsc, jec + + call mpp_get_layout(input_domain,layout) + call mpp_get_global_domain(input_domain, xsize=xsz, ysize=ysz) + if(PRESENT(maskmap)) then + call mpp_define_domains((/1,xsz,1,ysz/),layout,Ocean_sfc%Domain, maskmap=maskmap) + else + call mpp_define_domains((/1,xsz,1,ysz/),layout,Ocean_sfc%Domain) + endif + call mpp_get_compute_domain(Ocean_sfc%Domain, isc, iec, jsc, jec) + + allocate ( Ocean_sfc%t_surf (isc:iec,jsc:jec), & + Ocean_sfc%s_surf (isc:iec,jsc:jec), & + Ocean_sfc%u_surf (isc:iec,jsc:jec), & + Ocean_sfc%v_surf (isc:iec,jsc:jec), & + Ocean_sfc%sea_lev(isc:iec,jsc:jec), & + Ocean_sfc%area (isc:iec,jsc:jec), & + Ocean_sfc%frazil (isc:iec,jsc:jec)) + + Ocean_sfc%t_surf = 0.0 ! time averaged sst (Kelvin) passed to atmosphere/ice model + Ocean_sfc%s_surf = 0.0 ! time averaged sss (psu) passed to atmosphere/ice models + Ocean_sfc%u_surf = 0.0 ! time averaged u-current (m/sec) passed to atmosphere/ice models + Ocean_sfc%v_surf = 0.0 ! time averaged v-current (m/sec) passed to atmosphere/ice models + Ocean_sfc%sea_lev = 0.0 ! time averaged thickness of top model grid cell (m) plus patm/rho0/grav + Ocean_sfc%frazil = 0.0 ! time accumulated frazil (J/m^2) passed to ice model + Ocean_sfc%area = 0.0 + Ocean_sfc%axes = diag%axesT1%handles !diag axes to be used by coupler tracer flux diagnostics + + if (present(gas_fields_ocn)) then + call coupler_type_spawn(gas_fields_ocn, Ocean_sfc%fields, (/isc,isc,iec,iec/), & + (/jsc,jsc,jec,jec/), suffix = '_ocn', as_needed=.true.) + endif + +end subroutine initialize_ocean_public_type + +!> Translates the coupler's ocean_data_type into MOM6's surface state variable. +!! This may eventually be folded into the MOM6's code that calculates the +!! surface state in the first place. Note the offset in the arrays because +!! the ocean_data_type has no halo points in its arrays and always uses +!! absolute indicies. +subroutine convert_state_to_ocean_type(state, Ocean_sfc, G, use_conT_absS, & + patm, press_to_z) + type(surface), intent(inout) :: state + type(ocean_public_type), target, intent(inout) :: Ocean_sfc + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + logical, intent(in) :: use_conT_absS + real, optional, intent(in) :: patm(:,:) + real, optional, intent(in) :: press_to_z + + ! local variables + real :: IgR0 + character(len=48) :: val_str + integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd + integer :: i, j, i0, j0, is, ie, js, je + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + call pass_vector(state%u,state%v,G%Domain) + + call mpp_get_compute_domain(Ocean_sfc%Domain, isc_bnd, iec_bnd, & + jsc_bnd, jec_bnd) + if (present(patm)) then + ! Check that the inidicies in patm are (isc_bnd:iec_bnd,jsc_bnd:jec_bnd). + if (.not.present(press_to_z)) call MOM_error(FATAL, & + 'convert_state_to_ocean_type: press_to_z must be present if patm is.') + endif + + i0 = is - isc_bnd ; j0 = js - jsc_bnd + !If directed convert the surface T&S + !from conservative T to potential T and + !from absolute (reference) salinity to practical salinity + ! + if(use_conT_absS) then + do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd + Ocean_sfc%s_surf(i,j) = gsw_sp_from_sr(state%SSS(i+i0,j+j0)) + Ocean_sfc%t_surf(i,j) = gsw_pt_from_ct(state%SSS(i+i0,j+j0),state%SST(i+i0,j+j0)) + CELSIUS_KELVIN_OFFSET + enddo ; enddo + else + do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd + Ocean_sfc%t_surf(i,j) = state%SST(i+i0,j+j0) + CELSIUS_KELVIN_OFFSET + Ocean_sfc%s_surf(i,j) = state%SSS(i+i0,j+j0) + enddo ; enddo + endif + + do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd + Ocean_sfc%sea_lev(i,j) = state%sea_lev(i+i0,j+j0) + if (present(patm)) & + Ocean_sfc%sea_lev(i,j) = Ocean_sfc%sea_lev(i,j) + patm(i,j) * press_to_z + if (associated(state%frazil)) & + Ocean_sfc%frazil(i,j) = state%frazil(i+i0,j+j0) + Ocean_sfc%area(i,j) = G%areaT(i+i0,j+j0) + enddo ; enddo + + if (Ocean_sfc%stagger == AGRID) then + do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd + Ocean_sfc%u_surf(i,j) = G%mask2dT(i+i0,j+j0)*0.5*(state%u(I+i0,j+j0)+state%u(I-1+i0,j+j0)) + Ocean_sfc%v_surf(i,j) = G%mask2dT(i+i0,j+j0)*0.5*(state%v(i+i0,J+j0)+state%v(i+i0,J-1+j0)) + enddo ; enddo + elseif (Ocean_sfc%stagger == BGRID_NE) then + do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd + Ocean_sfc%u_surf(i,j) = G%mask2dBu(I+i0,J+j0)*0.5*(state%u(I+i0,j+j0)+state%u(I+i0,j+j0+1)) + Ocean_sfc%v_surf(i,j) = G%mask2dBu(I+i0,J+j0)*0.5*(state%v(i+i0,J+j0)+state%v(i+i0+1,J+j0)) + enddo ; enddo + elseif (Ocean_sfc%stagger == CGRID_NE) then + do j=jsc_bnd,jec_bnd ; do i=isc_bnd,iec_bnd + Ocean_sfc%u_surf(i,j) = G%mask2dCu(I+i0,j+j0)*state%u(I+i0,j+j0) + Ocean_sfc%v_surf(i,j) = G%mask2dCv(i+i0,J+j0)*state%v(i+i0,J+j0) + enddo ; enddo + else + write(val_str, '(I8)') Ocean_sfc%stagger + call MOM_error(FATAL, "convert_state_to_ocean_type: "//& + "Ocean_sfc%stagger has the unrecognized value of "//trim(val_str)) + endif + + if (coupler_type_initialized(state%tr_fields)) then + if (.not.coupler_type_initialized(Ocean_sfc%fields)) then + call MOM_error(FATAL, "convert_state_to_ocean_type: "//& + "Ocean_sfc%fields has not been initialized.") + endif + call coupler_type_copy_data(state%tr_fields, Ocean_sfc%fields) + endif + +end subroutine convert_state_to_ocean_type + +!> Returns pointers to objects within ocean_state_type +subroutine get_state_pointers(OS, grid, surf) + type(ocean_state_type), pointer :: OS !< Ocean state type + type(ocean_grid_type), optional, pointer :: grid !< Ocean grid + type(surface), optional, pointer :: surf !< Ocean surface state + + if (present(grid)) grid => OS%grid + if (present(surf)) surf=> OS%state + +end subroutine get_state_pointers + !> Maps outgoing ocean data to MCT buffer subroutine ocn_export(ind, ocn_public, grid, o2x) type(cpl_indices), intent(inout) :: ind !< Structure with coupler @@ -741,12 +1232,13 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) type(ESMF_time) :: time_start_ESMF !< Time at the start of the coupling interval type(ESMF_timeInterval) :: ocn_cpl_interval !< The length of one ocean coupling interval integer :: year, month, day, hour, minute, seconds, seconds_n, seconds_d, rc - logical :: write_restart_at_eod + logical :: write_restart_at_eod !< Controls if restart files must be written type(time_type) :: time_start !< Start of coupled time interval to pass to MOM6 type(time_type) :: coupling_timestep !< Coupled time interval to pass to MOM6 - character(len=128) :: err_msg + character(len=128) :: err_msg !< Error message character(len=32) :: timestamp !< Name of intermediate restart file - + character(len=80) :: restartname !< The restart file name (no dir) + character(len=384) :: runid !< Run ID ! Compute the time at the start of this coupling interval call ESMF_ClockGet(EClock, PrevTime=time_start_ESMF, rc=rc) call ESMF_TimeGet(time_start_ESMF, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) @@ -783,7 +1275,6 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) ! fill ice ocean boundary call fill_ice_ocean_bnd(glb%ice_ocean_boundary, glb%grid, x2o_o%rattr, glb%ind, glb%sw_decomp, & glb%c1, glb%c2, glb%c3, glb%c4) - if (debug .and. is_root_pe()) write(6,*) 'fill_ice_ocean_bnd' call update_ocean_model(glb%ice_ocean_boundary, glb%ocn_state, glb%ocn_public, & time_start, coupling_timestep) @@ -793,15 +1284,192 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) write_restart_at_eod = seq_timemgr_RestartAlarmIsOn(EClock) if (debug .and. is_root_pe()) write(6,*) 'ocn_run_mct, write_restart_at_eod=', write_restart_at_eod - if (write_restart_at_eod) then - !timestamp = date_to_string(EClock) - ! \todo add time stamp to ocean_model_restart + !!if (write_restart_at_eod) then + ! case name + !!call seq_infodata_GetData( glb%infodata, case_name=runid ) + ! add time stamp to restart filename + ! \todo use MCT call to get year, month etc + !!call ESMF_TimeGet(current_time, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) + !call get_date(glb%ocn_state%Time,year,month,days,hour,minute,seconds) + !!seconds = seconds + hour*3600 + minute*60 + !!write(restartname,'(A,".mom6.",I4.4,"-",I2.2,"-",I2.2,"-",I5.5,".nc")') runid, year, month, day, seconds !call ocean_model_restart(glb%ocn_state, timestamp) - call ocean_model_restart(glb%ocn_state) - endif + !!call save_restart(glb%ocn_state%dirs%restart_output_dir, glb%ocn_state%Time, glb%grid, & + !! glb%ocn_state%MOM_CSp%restart_CSp, filename=restartname,GV=glb%ocn_state%GV) + + !! GM call via ocean_model_MOM.F90 call ocean_model_restart(glb%ocn_state) + + ! Is this needed? + !call forcing_save_restart(OS%forcing_CSp, OS%grid, OS%Time, & + ! OS%dirs%restart_output_dir, .true.) + ! Once we start using the ice shelf module, the following will be needed + !if (glb%ocn_state%use_ice_shelf) then + ! call ice_shelf_save_restart(OS%Ice_shelf_CSp, OS%Time, OS%dirs%restart_output_dir, .true.) + !endif + + !!endif end subroutine ocn_run_mct +!> Updates the ocean model fields. This code wraps the call to step_MOM with MOM6's call. +!! It uses the forcing in Ice_ocean_boundary to advance the ocean model's state from the +!! input value of Ocean_state (which must be for time time_start_update) for a time interval +!! of Ocean_coupling_time_step, eturning the publicly visible ocean surface properties in +!! Ocean_sfc and storing the new ocean properties in Ocean_state. +subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & + time_start_update, Ocean_coupling_time_step) + type(ice_ocean_boundary_type), intent(inout) :: Ice_ocean_boundary + type(ocean_state_type), pointer :: OS + type(ocean_public_type), intent(inout) :: Ocean_sfc + type(time_type), intent(in) :: time_start_update + type(time_type), intent(in) :: Ocean_coupling_time_step + +! Arguments: Ice_ocean_boundary - A structure containing the various forcing +! fields coming from the ice. It is intent in. +! (inout) Ocean_state - A structure containing the internal ocean state. +! (out) Ocean_sfc - A structure containing all the publicly visible ocean +! surface fields after a coupling time step. +! (in) time_start_update - The time at the beginning of the update step. +! (in) Ocean_coupling_time_step - The amount of time over which to advance +! the ocean. + +! Note: although several types are declared intent(inout), this is to allow for +! the possibility of halo updates and to keep previously allocated memory. +! In practice, Ice_ocean_boundary is intent in, Ocean_state is private to +! this module and intent inout, and Ocean_sfc is intent out. + type(time_type) :: Master_time ! This allows step_MOM to temporarily change + ! the time that is seen by internal modules. + type(time_type) :: Time1 ! The value of the ocean model's time at the + ! start of a call to step_MOM. + integer :: index_bnds(4) ! The computational domain index bounds in the + ! ice-ocean boundary type. + real :: weight ! Flux accumulation weight + real :: time_step ! The time step of a call to step_MOM in seconds. + integer :: secs, days + integer :: is, ie, js, je + + call callTree_enter("update_ocean_model(), ocn_comp_mct.F90") + call get_time(Ocean_coupling_time_step, secs, days) + time_step = 86400.0*real(days) + real(secs) + + if (time_start_update /= OS%Time) then + call MOM_error(WARNING, "update_ocean_model: internal clock does not "//& + "agree with time_start_update argument.") + endif + + if (.not.associated(OS)) then + call MOM_error(FATAL, "update_ocean_model called with an unassociated "// & + "ocean_state_type structure. ocean_model_init must be "// & + "called first to allocate this structure.") + return + endif + + ! This is benign but not necessary if ocean_model_init_sfc was called or if + ! OS%state%tr_fields was spawnded in ocean_model_init. Consider removing it. + is = OS%grid%isc ; ie = OS%grid%iec ; js = OS%grid%jsc ; je = OS%grid%jec + call coupler_type_spawn(Ocean_sfc%fields, OS%state%tr_fields, & + (/is,is,ie,ie/), (/js,js,je,je/), as_needed=.true.) + +! Translate Ice_ocean_boundary into fluxes. + call mpp_get_compute_domain(Ocean_sfc%Domain, index_bnds(1), index_bnds(2), & + index_bnds(3), index_bnds(4)) + + weight = 1.0 + + if (OS%fluxes%fluxes_used) then + call enable_averaging(time_step, OS%Time + Ocean_coupling_time_step, OS%MOM_CSp%diag) ! Needed to allow diagnostics in convert_IOB + call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%fluxes, index_bnds, OS%Time, & + OS%grid, OS%forcing_CSp, OS%state, OS%restore_salinity,OS%restore_temp) +#ifdef _USE_GENERIC_TRACER + call MOM_generic_tracer_fluxes_accumulate(OS%fluxes, weight) !here weight=1, just saving the current fluxes +#endif + + ! Add ice shelf fluxes + if (OS%use_ice_shelf) then + call shelf_calc_flux(OS%State, OS%fluxes, OS%Time, time_step, OS%Ice_shelf_CSp) + endif + + ! GMM, check ocean_model_MOM.F90 to enable the following option + !if (OS%icebergs_apply_rigid_boundary) then + !This assumes that the iceshelf and ocean are on the same grid. I hope this is true + ! call add_berg_flux_to_shelf(OS%grid, OS%fluxes,OS%use_ice_shelf,OS%density_iceberg,OS%kv_iceberg, OS%latent_heat_fusion, OS%State, time_step, OS%berg_area_threshold) + !endif + + ! Indicate that there are new unused fluxes. + OS%fluxes%fluxes_used = .false. + OS%fluxes%dt_buoy_accum = time_step + else + OS%flux_tmp%C_p = OS%fluxes%C_p + call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%flux_tmp, index_bnds, OS%Time, & + OS%grid, OS%forcing_CSp, OS%state, OS%restore_salinity,OS%restore_temp) + if (OS%use_ice_shelf) then + call shelf_calc_flux(OS%State, OS%flux_tmp, OS%Time, time_step, OS%Ice_shelf_CSp) + endif + + ! GMM, check ocean_model_MOM.F90 to enable the following option + !if (OS%icebergs_apply_rigid_boundary) then + !This assumes that the iceshelf and ocean are on the same grid. I hope this is true + ! call add_berg_flux_to_shelf(OS%grid, OS%flux_tmp, OS%use_ice_shelf,OS%density_iceberg,OS%kv_iceberg, OS%latent_heat_fusion, OS%State, time_step, OS%berg_area_threshold) + !endif + + call forcing_accumulate(OS%flux_tmp, OS%fluxes, time_step, OS%grid, weight) +#ifdef _USE_GENERIC_TRACER + call MOM_generic_tracer_fluxes_accumulate(OS%flux_tmp, weight) !weight of the current flux in the running average +#endif + endif + + if (OS%nstep==0) then + call finish_MOM_initialization(OS%Time, OS%dirs, OS%MOM_CSp, OS%fluxes) + + call write_energy(OS%MOM_CSp%u, OS%MOM_CSp%v, OS%MOM_CSp%h, OS%MOM_CSp%tv, & + OS%Time, 0, OS%grid, OS%GV, OS%sum_output_CSp, & + OS%MOM_CSp%tracer_flow_CSp) + endif + + call disable_averaging(OS%MOM_CSp%diag) + Master_time = OS%Time ; Time1 = OS%Time + + if(OS%MOM_Csp%offline_tracer_mode) then + call step_offline(OS%fluxes, OS%state, Time1, time_step, OS%MOM_CSp) + else + call step_MOM(OS%fluxes, OS%state, Time1, time_step, OS%MOM_CSp) + endif + + OS%Time = Master_time + Ocean_coupling_time_step + OS%nstep = OS%nstep + 1 + + call enable_averaging(time_step, OS%Time, OS%MOM_CSp%diag) + call mech_forcing_diags(OS%fluxes, time_step, OS%grid, & + OS%MOM_CSp%diag, OS%forcing_CSp%handles) + call disable_averaging(OS%MOM_CSp%diag) + + if (OS%fluxes%fluxes_used) then + call enable_averaging(OS%fluxes%dt_buoy_accum, OS%Time, OS%MOM_CSp%diag) + call forcing_diagnostics(OS%fluxes, OS%state, OS%fluxes%dt_buoy_accum, & + OS%grid, OS%MOM_CSp%diag, OS%forcing_CSp%handles) + call accumulate_net_input(OS%fluxes, OS%state, OS%fluxes%dt_buoy_accum, & + OS%grid, OS%sum_output_CSp) + call disable_averaging(OS%MOM_CSp%diag) + endif + +! See if it is time to write out the energy. + if ((OS%Time + ((Ocean_coupling_time_step)/2) > OS%write_energy_time) .and. & + (OS%MOM_CSp%t_dyn_rel_adv==0.0)) then + call write_energy(OS%MOM_CSp%u, OS%MOM_CSp%v, OS%MOM_CSp%h, OS%MOM_CSp%tv, & + OS%Time, OS%nstep, OS%grid, OS%GV, OS%sum_output_CSp, & + OS%MOM_CSp%tracer_flow_CSp) + OS%write_energy_time = OS%write_energy_time + OS%energysavedays + endif + +! Translate state into Ocean. +! call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid, & +! Ice_ocean_boundary%p, OS%press_to_z) + call convert_state_to_ocean_type(OS%state, Ocean_sfc, OS%grid, & + OS%MOM_CSp%use_conT_absS) + + call callTree_leave("update_ocean_model()") +end subroutine update_ocean_model + !> Finalizes MOM6 !! !! \todo This needs to be done here. From 5ecf3bc1b13018794821d86ab0237c277c8d0cd7 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 28 Aug 2017 16:16:40 -0600 Subject: [PATCH 09/20] Writes restart files using CESM convention This PR allows MOM6 to write restart files using the CESM convention ($CASE.mom6.r.yyyy-mm-dd-sssss.nc). TODO: * Check if forcing_save_restart is needed; * Apply restart name convention to ice_shelf_save_restart --- config_src/mct_driver/ocn_comp_mct.F90 | 55 +++++++++++++------------- 1 file changed, 27 insertions(+), 28 deletions(-) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index b750264589..cbdf846178 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -29,10 +29,11 @@ module ocn_comp_mct use MOM, only: initialize_MOM, step_MOM, MOM_control_struct, MOM_end use MOM, only: calculate_surface_state, allocate_surface_state use MOM, only: finish_MOM_initialization, step_offline -use MOM_surface_forcing, only: ice_ocean_boundary_type, surface_forcing_CS use MOM_forcing_type, only: forcing, forcing_diagnostics, mech_forcing_diags, forcing_accumulate use MOM_forcing_type, only: allocate_forcing_type use MOM_surface_forcing, only: surface_forcing_init, convert_IOB_to_fluxes +use MOM_surface_forcing, only: ice_ocean_boundary_type, surface_forcing_CS +use MOM_surface_forcing, only: forcing_save_restart use MOM_restart, only: save_restart use MOM_domains, only: MOM_infra_init, num_pes, root_pe, pe_here use MOM_domains, only: pass_vector, BGRID_NE, CGRID_NE @@ -43,18 +44,19 @@ module ocn_comp_mct use MOM_error_handler, only: MOM_error, FATAL, is_root_pe, WARNING use MOM_error_handler, only: callTree_enter, callTree_leave use MOM_time_manager, only: time_type, set_date, set_time, set_calendar_type, NOLEAP, get_date -use MOM_time_manager, only : operator(+), operator(-), operator(*), operator(/) -use MOM_time_manager, only : operator(/=), operator(>), get_time +use MOM_time_manager, only: operator(+), operator(-), operator(*), operator(/) +use MOM_time_manager, only: operator(/=), operator(>), get_time use MOM_file_parser, only: get_param, log_version, param_file_type, close_param_file use MOM_get_input, only: Get_MOM_Input, directories -use MOM_diag_mediator, only : diag_ctrl, enable_averaging, disable_averaging -use MOM_diag_mediator, only : diag_mediator_close_registration, diag_mediator_end +use MOM_diag_mediator, only: diag_ctrl, enable_averaging, disable_averaging +use MOM_diag_mediator, only: diag_mediator_close_registration, diag_mediator_end use MOM_ice_shelf, only: initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS +use MOM_ice_shelf, only: ice_shelf_end, ice_shelf_save_restart use MOM_sum_output, only: MOM_sum_output_init, sum_output_CS use MOM_sum_output, only: write_energy, accumulate_net_input -use MOM_string_functions, only : uppercase -use MOM_constants, only : CELSIUS_KELVIN_OFFSET, hlf -use MOM_EOS, only : gsw_sp_from_sr, gsw_pt_from_ct +use MOM_string_functions, only: uppercase +use MOM_constants, only: CELSIUS_KELVIN_OFFSET, hlf +use MOM_EOS, only: gsw_sp_from_sr, gsw_pt_from_ct ! FMS use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain @@ -1237,7 +1239,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) type(time_type) :: coupling_timestep !< Coupled time interval to pass to MOM6 character(len=128) :: err_msg !< Error message character(len=32) :: timestamp !< Name of intermediate restart file - character(len=80) :: restartname !< The restart file name (no dir) + character(len=384) :: restartname !< The restart file name (no dir) character(len=384) :: runid !< Run ID ! Compute the time at the start of this coupling interval call ESMF_ClockGet(EClock, PrevTime=time_start_ESMF, rc=rc) @@ -1284,30 +1286,27 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) write_restart_at_eod = seq_timemgr_RestartAlarmIsOn(EClock) if (debug .and. is_root_pe()) write(6,*) 'ocn_run_mct, write_restart_at_eod=', write_restart_at_eod - !!if (write_restart_at_eod) then + if (write_restart_at_eod) then ! case name - !!call seq_infodata_GetData( glb%infodata, case_name=runid ) - ! add time stamp to restart filename - ! \todo use MCT call to get year, month etc - !!call ESMF_TimeGet(current_time, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) - !call get_date(glb%ocn_state%Time,year,month,days,hour,minute,seconds) - !!seconds = seconds + hour*3600 + minute*60 - !!write(restartname,'(A,".mom6.",I4.4,"-",I2.2,"-",I2.2,"-",I5.5,".nc")') runid, year, month, day, seconds - !call ocean_model_restart(glb%ocn_state, timestamp) - !!call save_restart(glb%ocn_state%dirs%restart_output_dir, glb%ocn_state%Time, glb%grid, & - !! glb%ocn_state%MOM_CSp%restart_CSp, filename=restartname,GV=glb%ocn_state%GV) - - !! GM call via ocean_model_MOM.F90 call ocean_model_restart(glb%ocn_state) + call seq_infodata_GetData( glb%infodata, case_name=runid ) + ! add time stamp to the restart filename + call ESMF_ClockGet(EClock, CurrTime=current_time, rc=rc) + call ESMF_TimeGet(current_time, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) + seconds = seconds + hour*3600 + minute*60 + write(restartname,'(A,".mom6.r.",I4.4,"-",I2.2,"-",I2.2,"-",I5.5)') trim(runid), year, month, day, seconds + + call save_restart(glb%ocn_state%dirs%restart_output_dir, glb%ocn_state%Time, glb%grid, & + glb%ocn_state%MOM_CSp%restart_CSp, .false., filename=restartname,GV=glb%ocn_state%GV) ! Is this needed? - !call forcing_save_restart(OS%forcing_CSp, OS%grid, OS%Time, & - ! OS%dirs%restart_output_dir, .true.) + call forcing_save_restart(glb%ocn_state%forcing_CSp, glb%grid, glb%ocn_state%Time, & + glb%ocn_state%dirs%restart_output_dir, .true.) ! Once we start using the ice shelf module, the following will be needed - !if (glb%ocn_state%use_ice_shelf) then - ! call ice_shelf_save_restart(OS%Ice_shelf_CSp, OS%Time, OS%dirs%restart_output_dir, .true.) - !endif + if (glb%ocn_state%use_ice_shelf) then + call ice_shelf_save_restart(glb%ocn_state%Ice_shelf_CSp, glb%ocn_state%Time, glb%ocn_state%dirs%restart_output_dir, .true.) + endif - !!endif + endif end subroutine ocn_run_mct From 2733fd29ce600fd53ad6fe4b68ffa7d170aae2a9 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 29 Aug 2017 11:49:14 -0600 Subject: [PATCH 10/20] writes name of restart file in ocn rpointer file --- config_src/mct_driver/ocn_comp_mct.F90 | 52 +++++++++++++++++++++++--- 1 file changed, 47 insertions(+), 5 deletions(-) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index cbdf846178..214e9b57aa 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -24,6 +24,7 @@ module ocn_comp_mct use seq_timemgr_mod, only: seq_timemgr_EClockGetData, seq_timemgr_RestartAlarmIsOn use perf_mod, only: t_startf, t_stopf use shr_kind_mod, only: shr_kind_r8 +use shr_file_mod, only: shr_file_getUnit, shr_file_freeUnit ! MOM6 modules use MOM, only: initialize_MOM, step_MOM, MOM_control_struct, MOM_end @@ -252,8 +253,9 @@ module ocn_comp_mct type(cpl_indices), public :: ind !< Variable IDs ! runtime params logical :: sw_decomp !< Controls whether shortwave is decomposed into four components - real :: c1, c2, c3, c4 !< Coeffs. used in the shortwave decomposition - + real :: c1, c2, c3, c4 !< Coeffs. used in the shortwave decomposition + character(len=384) :: pointer_filename !< Name of the ascii file that contains the path + !! and filename of the latest restart file. end type MCT_MOM_Data type(MCT_MOM_Data) :: glb !< global structure @@ -402,6 +404,9 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) ! read useful runtime params call get_MOM_Input(param_file, dirs_tmp, check_params=.false.) !call log_version(param_file, mdl, version, "") + call get_param(param_file, mdl, "POINTER_FILENAME", glb%pointer_filename, & + "Name of the ascii file that contains the path and filename of" // & + " the latest restart file.", default='rpointer.ocn') call get_param(param_file, mdl, "SW_DECOMP", glb%sw_decomp, & "If True, read coeffs c1, c2, c3 and c4 and decompose" // & "the net shortwave radiation (SW) into four components:\n" // & @@ -426,13 +431,18 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) glb%c1 = 0.0; glb%c2 = 0.0; glb%c3 = 0.0; glb%c4 = 0.0 endif - restartfile = '/glade/scratch/gmarques/mom_test/run/RESTART/MOM.res_Y0001_D002_S00000.nc' ! Initialize the MOM6 model if (runtype == "initial") then ! startup (new run) - 'n' is needed below since we don't ! specify input_filename in input.nml call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in, input_restart_file = 'n') else ! hybrid or branch or continuos runs - !call seq_infodata_GetData(glb%infodata, restart_file=restartfile) + ! read pointer_filename + !call seq_infodata_GetData( glb%infodata, case_name=runid ) + !call ESMF_ClockGet(EClock, CurrTime=current_time, rc=rc) + !call ESMF_TimeGet(current_time, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) + !seconds = seconds + hour*3600 + minute*60 + !write(restartfile,'(A,".mom6.r.",I4.4,"-",I2.2,"-",I2.2,"-",I5.5,".nc")') trim(runid), year, month, day, seconds + call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in, input_restart_file = restartfile) endif @@ -1240,7 +1250,9 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) character(len=128) :: err_msg !< Error message character(len=32) :: timestamp !< Name of intermediate restart file character(len=384) :: restartname !< The restart file name (no dir) - character(len=384) :: runid !< Run ID + character(len=384) :: restart_pointer_file !< File name for restart pointer file + character(len=384) :: runid !< Run ID + integer :: nu !< i/o unit to write pointer file ! Compute the time at the start of this coupling interval call ESMF_ClockGet(EClock, PrevTime=time_start_ESMF, rc=rc) call ESMF_TimeGet(time_start_ESMF, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) @@ -1298,6 +1310,17 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) call save_restart(glb%ocn_state%dirs%restart_output_dir, glb%ocn_state%Time, glb%grid, & glb%ocn_state%MOM_CSp%restart_CSp, .false., filename=restartname,GV=glb%ocn_state%GV) + ! write name of restart file in the rpointer file + nu = shr_file_getUnit() + if (is_root_pe()) then + restart_pointer_file = trim(glb%pointer_filename) + open(nu, file=restart_pointer_file, form='formatted', status='unknown') + write(nu,'(a)') trim(restartname) + close(nu) + write(6,*) 'ocn restart pointer file written: ',trim(restartname) + endif + call shr_file_freeUnit(nu) + ! Is this needed? call forcing_save_restart(glb%ocn_state%forcing_CSp, glb%grid, glb%ocn_state%Time, & glb%ocn_state%dirs%restart_output_dir, .true.) @@ -1478,8 +1501,27 @@ subroutine ocn_final_mct( EClock, cdata_o, x2o_o, o2x_o) type(mct_aVect) , intent(inout) :: x2o_o !< Fluxes from coupler to ocean, computed by ocean type(mct_aVect) , intent(inout) :: o2x_o !< Fluxes from ocean to coupler, computed by ocean + call ocean_model_end(glb%ocn_public, glb%ocn_state, glb%ocn_state%Time) + end subroutine ocn_final_mct +!> Terminates the model run, saving the ocean state in a +!! restart file and deallocating any data associated with the ocean. +subroutine ocean_model_end(Ocean_sfc, Ocean_state, Time) + type(ocean_public_type), intent(inout) :: Ocean_sfc !< An ocean_public_type structure that is to be + !! deallocated upon termination. + type(ocean_state_type), pointer :: Ocean_state!< pointer to the structure containing the internal + ! !! ocean state to be deallocated upon termination. + type(time_type), intent(in) :: Time !< The model time, used for writing restarts. + + if (debug .and. is_root_pe()) write(6,*)'Here 1' + !GMM call save_restart(Ocean_state, Time) + call diag_mediator_end(Time, Ocean_state%MOM_CSp%diag) + call MOM_end(Ocean_state%MOM_CSp) + if (Ocean_state%use_ice_shelf) call ice_shelf_end(Ocean_state%Ice_shelf_CSp) + if (debug .and. is_root_pe()) write(6,*)'Here 2' + +end subroutine ocean_model_end !> Sets mct global segment maps for the MOM decomposition. !! From f14d6640ff935bac9ac7573659ab4c95790fb05d Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 1 Sep 2017 13:10:25 -0600 Subject: [PATCH 11/20] Reads a restart file using the ocn rpointer file --- config_src/mct_driver/ocn_comp_mct.F90 | 38 ++++++++++++++++++++------ 1 file changed, 29 insertions(+), 9 deletions(-) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index 214e9b57aa..c94302432a 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -278,9 +278,12 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) type(ESMF_time) :: time_in_ESMF !< Time after first ocean coupling interval (of type ESMF_time) type(ESMF_timeInterval) :: ocn_cpl_interval !< Ocean coupling interval integer :: year, month, day, hour, minute, seconds, seconds_n, seconds_d, rc - character(len=384) :: runid !< Run ID - character(len=384) :: runtype !< Run type - character(len=384) :: restartfile !< Path/Name of restart file + character(len=240) :: runid !< Run ID + character(len=240) :: runtype !< Run type + character(len=240) :: restartfile !< Path/Name of restart file + integer :: nu !< i/o unit to read pointer file + character(len=240) :: restart_pointer_file !< File name for restart pointer file + character(len=240) :: restartpath !< Path of the restart file character(len=32) :: starttype !< infodata start type integer :: mpicom_ocn !< MPI ocn communicator integer :: npes, pe0 !< # of processors and current processor @@ -436,19 +439,36 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) ! specify input_filename in input.nml call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in, input_restart_file = 'n') else ! hybrid or branch or continuos runs + ! output path root + call seq_infodata_GetData( glb%infodata, outPathRoot=restartpath ) ! read pointer_filename + ! write name of restart file in the rpointer file + nu = shr_file_getUnit() + !if (is_root_pe()) then + restart_pointer_file = trim(glb%pointer_filename) + write(6,*) 'Reading ocn pointer file: ',restart_pointer_file + open(nu, file=restart_pointer_file, form='formatted', status='unknown') + read(nu,'(a)') restartfile + close(nu) + !restartfile = trim(restartpath) // trim(restartfile) + write(6,*) 'Reading ocn pointer file: ',trim(restartfile) + !endif + call shr_file_freeUnit(nu) + !restartfile = '/glade/scratch/gmarques/mom_test/run/mom_test.mom6.r.0001-01-01-21600.nc' !call seq_infodata_GetData( glb%infodata, case_name=runid ) !call ESMF_ClockGet(EClock, CurrTime=current_time, rc=rc) !call ESMF_TimeGet(current_time, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) !seconds = seconds + hour*3600 + minute*60 !write(restartfile,'(A,".mom6.r.",I4.4,"-",I2.2,"-",I2.2,"-",I5.5,".nc")') trim(runid), year, month, day, seconds - - call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in, input_restart_file = restartfile) + if (debug .and. root_pe().eq.pe_here()) then + write(6,*)'runtype, restartfile,time_init,time_in',runtype, restartfile !, time_init, time_in + endif + call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in, input_restart_file=trim(restartfile)) endif - if (debug .and. root_pe().eq.pe_here()) then - write(6,*)'runtype, restartfile',runtype, restartfile - endif + !if (debug .and. root_pe().eq.pe_here()) then + ! write(6,*)'runtype, restartfile,time_init,time_in',runtype, restartfile,time_init,time_in + !endif ! Initialize the MOM6 model !call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in) @@ -1315,7 +1335,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) if (is_root_pe()) then restart_pointer_file = trim(glb%pointer_filename) open(nu, file=restart_pointer_file, form='formatted', status='unknown') - write(nu,'(a)') trim(restartname) + write(nu,'(a)') trim(restartname) //'.nc' close(nu) write(6,*) 'ocn restart pointer file written: ',trim(restartname) endif From e734271204b8b52a3a24c67a2faeff951fe73eca Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 7 Sep 2017 12:41:55 -0600 Subject: [PATCH 12/20] Fixed initial time for cont'd runs. Implemented ocn.log file --- config_src/mct_driver/ocn_comp_mct.F90 | 112 +++++++++++++++++++------ 1 file changed, 86 insertions(+), 26 deletions(-) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index c94302432a..bf1a58aeaf 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -24,7 +24,9 @@ module ocn_comp_mct use seq_timemgr_mod, only: seq_timemgr_EClockGetData, seq_timemgr_RestartAlarmIsOn use perf_mod, only: t_startf, t_stopf use shr_kind_mod, only: shr_kind_r8 -use shr_file_mod, only: shr_file_getUnit, shr_file_freeUnit +use shr_file_mod, only: shr_file_getUnit, shr_file_freeUnit, shr_file_setIO, & + shr_file_getLogUnit, shr_file_getLogLevel, & + shr_file_setLogUnit, shr_file_setLogLevel ! MOM6 modules use MOM, only: initialize_MOM, step_MOM, MOM_control_struct, MOM_end @@ -254,8 +256,10 @@ module ocn_comp_mct ! runtime params logical :: sw_decomp !< Controls whether shortwave is decomposed into four components real :: c1, c2, c3, c4 !< Coeffs. used in the shortwave decomposition + ! i/o character(len=384) :: pointer_filename !< Name of the ascii file that contains the path !! and filename of the latest restart file. + integer :: stdout !< standard output unit. (by default, it should point to ocn.log.* file) end type MCT_MOM_Data type(MCT_MOM_Data) :: glb !< global structure @@ -275,7 +279,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) type(time_type) :: time_init !< Start time of coupled model's calendar type(time_type) :: time_in !< Time at the beginning of the first ocn coupling interval type(ESMF_time) :: current_time !< Current time - type(ESMF_time) :: time_in_ESMF !< Time after first ocean coupling interval (of type ESMF_time) + type(ESMF_time) :: time_in_ESMF !< Initial time for ocean type(ESMF_timeInterval) :: ocn_cpl_interval !< Ocean coupling interval integer :: year, month, day, hour, minute, seconds, seconds_n, seconds_d, rc character(len=240) :: runid !< Run ID @@ -313,6 +317,10 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) seconds_in_day = 86400.0d0, & minutes_in_hour = 60.0d0 + ! ocn model input namelist filename + character(len=99) :: ocn_modelio_name + integer :: shrlogunit ! original log file unit + integer :: shrloglev ! original log level ! instance control vars (these are local for now) integer(kind=4) :: inst_index @@ -346,7 +354,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) else if (trim(starttype) == trim(seq_infodata_start_type_brnch)) then runtype = "branch" else - write(*,*) 'ocn_comp_mct ERROR: unknown starttype' + write(glb%stdout,*) 'ocn_comp_mct ERROR: unknown starttype' call exit(0) end if @@ -360,6 +368,32 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) ! Initialize MOM6 comm call MOM_infra_init(mpicom_ocn) + ! initialize ocn log file + if (is_root_pe()) then + + ! get original log file properties + call shr_file_getLogUnit (shrlogunit) + call shr_file_getLogLevel(shrloglev) + + glb%stdout = shr_file_getUnit() ! get an unused unit number + + ! open the ocn_modelio.nml file and then open a log file associated with stdout + ocn_modelio_name = 'ocn_modelio.nml' // trim(inst_suffix) + call shr_file_setIO(ocn_modelio_name,glb%stdout) + + !if (debug) write(*,*) "stdout-1" + !if (debug) write(6,*) "stdout-2" + !if (debug) write(glb%stdout,*) "stdout-3" + + ! set the shr log io unit number + call shr_file_setLogUnit(glb%stdout) + + !if (debug) write(*,*) "stdout-4" + !if (debug) write(6,*) "stdout-5" + !if (debug) write(glb%stdout,*) "stdout-6" + + end if + call set_calendar_type(NOLEAP) !TODO: confirm this ! Get the ESMF clock instance (assigned by CESM for MOM6) @@ -370,28 +404,32 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) time_init = set_date(year, month, day, hour, minute, seconds, err_msg=err_msg) ! Compute time_in: time at the beginning of the first ocn coupling interval - ! (In CESM, ocean component is lagged by one ocean coupling interval, so: - ! time_in = time_init + ocn_cpl_interval ) call ESMF_ClockGet(EClock, TimeStep=ocn_cpl_interval, rc=rc) - time_in_ESMF = ESMF_TimeInc(current_time, ocn_cpl_interval) + if (runtype /= "continue") then + ! In startup runs, take the one ocn coupling interval lag into account to + ! compute the initial ocn time. (time_in = time_init + ocn_cpl_interval) + time_in_ESMF = ESMF_TimeInc(current_time, ocn_cpl_interval) + else + time_in_ESMF = current_time + endif call ESMF_TimeGet(time_in_ESMF, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) time_in = set_date(year, month, day, hour, minute, seconds, err_msg=err_msg) ! Debugging clocks if (debug .and. is_root_pe()) then - write(6,*) 'ocn_init_mct, current time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds + write(glb%stdout,*) 'ocn_init_mct, current time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds call ESMF_ClockGet(EClock, StartTime=current_time, rc=rc) call ESMF_TimeGet(current_time, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) - write(6,*) 'ocn_init_mct, start time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds + write(glb%stdout,*) 'ocn_init_mct, start time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds call ESMF_ClockGet(EClock, StopTime=current_time, rc=rc) call ESMF_TimeGet(current_time, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) - write(6,*) 'ocn_init_mct, stop time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds + write(glb%stdout,*) 'ocn_init_mct, stop time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds call ESMF_ClockGet(EClock, PrevTime=current_time, rc=rc) call ESMF_TimeGet(current_time, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) - write(6,*) 'ocn_init_mct, previous time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds + write(glb%stdout,*) 'ocn_init_mct, previous time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds call ESMF_ClockGet(EClock, TimeStep=ocn_cpl_interval, rc=rc) call ESMF_TimeIntervalGet(ocn_cpl_interval, yy=year, mm=month, d=day, s=seconds, sn=seconds_n, sd=seconds_d, rc=rc) - write(6,*) 'ocn_init_mct, time step: y,m,d-',year,month,day,'s,sn,sd=',seconds,seconds_n,seconds_d + write(glb%stdout,*) 'ocn_init_mct, time step: y,m,d-',year,month,day,'s,sn,sd=',seconds,seconds_n,seconds_d endif npes = num_pes() @@ -446,12 +484,12 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) nu = shr_file_getUnit() !if (is_root_pe()) then restart_pointer_file = trim(glb%pointer_filename) - write(6,*) 'Reading ocn pointer file: ',restart_pointer_file + write(glb%stdout,*) 'Reading ocn pointer file: ',restart_pointer_file open(nu, file=restart_pointer_file, form='formatted', status='unknown') read(nu,'(a)') restartfile close(nu) !restartfile = trim(restartpath) // trim(restartfile) - write(6,*) 'Reading ocn pointer file: ',trim(restartfile) + write(glb%stdout,*) 'Reading ocn pointer file: ',trim(restartfile) !endif call shr_file_freeUnit(nu) !restartfile = '/glade/scratch/gmarques/mom_test/run/mom_test.mom6.r.0001-01-01-21600.nc' @@ -461,13 +499,13 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) !seconds = seconds + hour*3600 + minute*60 !write(restartfile,'(A,".mom6.r.",I4.4,"-",I2.2,"-",I2.2,"-",I5.5,".nc")') trim(runid), year, month, day, seconds if (debug .and. root_pe().eq.pe_here()) then - write(6,*)'runtype, restartfile,time_init,time_in',runtype, restartfile !, time_init, time_in + write(glb%stdout,*)'runtype, restartfile,time_init,time_in',runtype, restartfile !, time_init, time_in endif call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in, input_restart_file=trim(restartfile)) endif !if (debug .and. root_pe().eq.pe_here()) then - ! write(6,*)'runtype, restartfile,time_init,time_in',runtype, restartfile,time_init,time_in + ! write(glb%stdout,*)'runtype, restartfile,time_init,time_in',runtype, restartfile,time_init,time_in !endif ! Initialize the MOM6 model @@ -527,7 +565,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) ncouple_per_day = seconds_in_day / ocn_cpl_dt mom_cpl_dt = seconds_in_day / ncouple_per_day if (mom_cpl_dt /= ocn_cpl_dt) then - write(*,*) 'ERROR mom_cpl_dt and ocn_cpl_dt must be identical' + write(glb%stdout,*) 'ERROR mom_cpl_dt and ocn_cpl_dt must be identical' call exit(0) end if @@ -584,6 +622,12 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) if (debug .and. root_pe().eq.pe_here()) print *, "leaving ocean_init_mct" + ! Reset shr logging to original values + if (is_root_pe()) then + call shr_file_setLogUnit (shrlogunit) + call shr_file_setLogLevel(shrloglev) + end if + end subroutine ocn_init_mct !> Determines attribute vector indices @@ -874,7 +918,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i call diag_mediator_close_registration(OS%MOM_CSp%diag) if (is_root_pe()) & - write(*,'(/12x,a/)') '======== COMPLETED MOM INITIALIZATION ========' + write(glb%stdout,'(/12x,a/)') '======== COMPLETED MOM INITIALIZATION ========' call callTree_leave("ocean_model_init(") end subroutine ocean_model_init @@ -1273,6 +1317,16 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) character(len=384) :: restart_pointer_file !< File name for restart pointer file character(len=384) :: runid !< Run ID integer :: nu !< i/o unit to write pointer file + integer :: shrlogunit ! original log file unit + integer :: shrloglev ! original log level + + ! reset shr logging to ocn log file: + if (is_root_pe()) then + call shr_file_getLogUnit(shrlogunit) + call shr_file_getLogLevel(shrloglev) + call shr_file_setLogUnit(glb%stdout) + endif + ! Compute the time at the start of this coupling interval call ESMF_ClockGet(EClock, PrevTime=time_start_ESMF, rc=rc) call ESMF_TimeGet(time_start_ESMF, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) @@ -1282,19 +1336,19 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) if (debug .and. is_root_pe()) then call ESMF_ClockGet(EClock, CurrTime=current_time, rc=rc) call ESMF_TimeGet(current_time, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) - write(6,*) 'ocn_run_mct, current time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds + write(glb%stdout,*) 'ocn_run_mct, current time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds call ESMF_ClockGet(EClock, StartTime=current_time, rc=rc) call ESMF_TimeGet(current_time, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) - write(6,*) 'ocn_run_mct, start time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds + write(glb%stdout,*) 'ocn_run_mct, start time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds call ESMF_ClockGet(EClock, StopTime=current_time, rc=rc) call ESMF_TimeGet(current_time, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) - write(6,*) 'ocn_run_mct, stop time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds + write(glb%stdout,*) 'ocn_run_mct, stop time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds call ESMF_ClockGet(EClock, PrevTime=current_time, rc=rc) call ESMF_TimeGet(current_time, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) - write(6,*) 'ocn_run_mct, previous time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds + write(glb%stdout,*) 'ocn_run_mct, previous time: y,m,d-',year,month,day,'h,m,s=',hour,minute,seconds call ESMF_ClockGet(EClock, TimeStep=ocn_cpl_interval, rc=rc) call ESMF_TimeIntervalGet(ocn_cpl_interval, yy=year, mm=month, d=day, s=seconds, sn=seconds_n, sd=seconds_d, rc=rc) - write(6,*) 'ocn_init_mct, time step: y,m,d-',year,month,day,'s,sn,sd=',seconds,seconds_n,seconds_d + write(glb%stdout,*) 'ocn_init_mct, time step: y,m,d-',year,month,day,'s,sn,sd=',seconds,seconds_n,seconds_d endif ! Translate the coupling time interval @@ -1316,7 +1370,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) !--- write out intermediate restart file when needed. ! Check alarms for flag to write restart at end of day write_restart_at_eod = seq_timemgr_RestartAlarmIsOn(EClock) - if (debug .and. is_root_pe()) write(6,*) 'ocn_run_mct, write_restart_at_eod=', write_restart_at_eod + if (debug .and. is_root_pe()) write(glb%stdout,*) 'ocn_run_mct, write_restart_at_eod=', write_restart_at_eod if (write_restart_at_eod) then ! case name @@ -1337,7 +1391,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) open(nu, file=restart_pointer_file, form='formatted', status='unknown') write(nu,'(a)') trim(restartname) //'.nc' close(nu) - write(6,*) 'ocn restart pointer file written: ',trim(restartname) + write(glb%stdout,*) 'ocn restart pointer file written: ',trim(restartname) endif call shr_file_freeUnit(nu) @@ -1351,6 +1405,12 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) endif + ! reset shr logging to original values + if (is_root_pe()) then + call shr_file_setLogUnit(shrlogunit) + call shr_file_setLogLevel(shrloglev) + endif + end subroutine ocn_run_mct !> Updates the ocean model fields. This code wraps the call to step_MOM with MOM6's call. @@ -1534,12 +1594,12 @@ subroutine ocean_model_end(Ocean_sfc, Ocean_state, Time) ! !! ocean state to be deallocated upon termination. type(time_type), intent(in) :: Time !< The model time, used for writing restarts. - if (debug .and. is_root_pe()) write(6,*)'Here 1' + if (debug .and. is_root_pe()) write(glb%stdout,*)'Here 1' !GMM call save_restart(Ocean_state, Time) call diag_mediator_end(Time, Ocean_state%MOM_CSp%diag) call MOM_end(Ocean_state%MOM_CSp) if (Ocean_state%use_ice_shelf) call ice_shelf_end(Ocean_state%Ice_shelf_CSp) - if (debug .and. is_root_pe()) write(6,*)'Here 2' + if (debug .and. is_root_pe()) write(glb%stdout,*)'Here 2' end subroutine ocean_model_end From 8c672046daabceeff985072d35bfb99984bbf5e8 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 7 Sep 2017 13:42:20 -0600 Subject: [PATCH 13/20] Adds option to pass latent heat flux via IOB Before, latent heat flux was computed by the ocean model using evap,lprec etc. This PR adds the option to pass latent, computed by the coupler, to the ocean via ice ocean boundary type (IOB). If IOB%latent has not been allocated, MOM6 will then compute latent heat flux as previously. --- .../coupled_driver/MOM_surface_forcing.F90 | 27 +++++++++++-------- config_src/mct_driver/ocn_comp_mct.F90 | 3 +++ 2 files changed, 19 insertions(+), 11 deletions(-) diff --git a/config_src/coupled_driver/MOM_surface_forcing.F90 b/config_src/coupled_driver/MOM_surface_forcing.F90 index 178a122087..46fdb2232b 100644 --- a/config_src/coupled_driver/MOM_surface_forcing.F90 +++ b/config_src/coupled_driver/MOM_surface_forcing.F90 @@ -152,6 +152,7 @@ module MOM_surface_forcing real, pointer, dimension(:,:) :: sw_flux_vis_dif =>NULL() ! diffuse visible sw radiation (W/m2) real, pointer, dimension(:,:) :: sw_flux_nir_dir =>NULL() ! direct Near InfraRed sw radiation (W/m2) real, pointer, dimension(:,:) :: sw_flux_nir_dif =>NULL() ! diffuse Near InfraRed sw radiation (W/m2) + real, pointer, dimension(:,:) :: latent =>NULL() ! latent heat flux (W/m2) real, pointer, dimension(:,:) :: lprec =>NULL() ! mass flux of liquid precip (kg/m2/s) real, pointer, dimension(:,:) :: fprec =>NULL() ! mass flux of frozen precip (kg/m2/s) real, pointer, dimension(:,:) :: runoff =>NULL() ! mass flux of liquid runoff (kg/m2/s) @@ -480,17 +481,21 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, state, fluxes%sens(i,j) = - IOB%t_flux(i-i0,j-j0) * G%mask2dT(i,j) fluxes%latent(i,j) = 0.0 - if (ASSOCIATED(IOB%fprec)) then - fluxes%latent(i,j) = fluxes%latent(i,j) - IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion - fluxes%latent_fprec_diag(i,j) = -G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion - endif - if (ASSOCIATED(IOB%calving)) then - fluxes%latent(i,j) = fluxes%latent(i,j) - IOB%calving(i-i0,j-j0)*CS%latent_heat_fusion - fluxes%latent_frunoff_diag(i,j) = -G%mask2dT(i,j) * IOB%calving(i-i0,j-j0)*CS%latent_heat_fusion - endif - if (ASSOCIATED(IOB%q_flux)) then - fluxes%latent(i,j) = fluxes%latent(i,j) - IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor - fluxes%latent_evap_diag(i,j) = -G%mask2dT(i,j) * IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor + if (ASSOCIATED(IOB%latent)) then + fluxes%latent(i,j) = IOB%latent(i-i0,j-j0) + else + if (ASSOCIATED(IOB%fprec)) then + fluxes%latent(i,j) = fluxes%latent(i,j) - IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion + fluxes%latent_fprec_diag(i,j) = -G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion + endif + if (ASSOCIATED(IOB%calving)) then + fluxes%latent(i,j) = fluxes%latent(i,j) - IOB%calving(i-i0,j-j0)*CS%latent_heat_fusion + fluxes%latent_frunoff_diag(i,j) = -G%mask2dT(i,j) * IOB%calving(i-i0,j-j0)*CS%latent_heat_fusion + endif + if (ASSOCIATED(IOB%q_flux)) then + fluxes%latent(i,j) = fluxes%latent(i,j) - IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor + fluxes%latent_evap_diag(i,j) = -G%mask2dT(i,j) * IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor + endif endif fluxes%latent(i,j) = G%mask2dT(i,j) * fluxes%latent(i,j) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index c94302432a..db82e95269 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -558,6 +558,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) allocate(glb%ice_ocean_boundary%v_flux(isc:iec,jsc:jec)); glb%ice_ocean_boundary%v_flux(:,:) = 0.0 allocate(glb%ice_ocean_boundary%t_flux(isc:iec,jsc:jec)); glb%ice_ocean_boundary%t_flux(:,:) = 0.0 allocate(glb%ice_ocean_boundary%q_flux(isc:iec,jsc:jec)); glb%ice_ocean_boundary%q_flux(:,:) = 0.0 + allocate(glb%ice_ocean_boundary%latent(isc:iec,jsc:jec)); glb%ice_ocean_boundary%latent(:,:) = 0.0 allocate(glb%ice_ocean_boundary%salt_flux(isc:iec,jsc:jec)); glb%ice_ocean_boundary%salt_flux(:,:) = 0.0 allocate(glb%ice_ocean_boundary%lw_flux(isc:iec,jsc:jec)); glb%ice_ocean_boundary%lw_flux(:,:) = 0.0 allocate(glb%ice_ocean_boundary%sw_flux_vis_dir(isc:iec,jsc:jec)); glb%ice_ocean_boundary%sw_flux_vis_dir(:,:) = 0.0 @@ -1214,6 +1215,8 @@ subroutine fill_ice_ocean_bnd(ice_ocean_boundary, grid, x2o_o, ind, & ice_ocean_boundary%v_flux(ig,jg) = x2o_o(ind%x2o_Foxx_tauy,k) ! sensible heat flux (W/m2) ice_ocean_boundary%t_flux(ig,jg) = -x2o_o(ind%x2o_Foxx_sen,k) + ! latent heat flux (W/m^2) + ice_ocean_boundary%latent(ig,jg) = x2o_o(ind%x2o_Foxx_lat,k) ! salt flux ice_ocean_boundary%salt_flux(ig,jg) = -x2o_o(ind%x2o_Fioi_salt,k) ! heat content from frozen runoff From a748f968dba7ee6c3c33ac7a76a88d8662b42a74 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 11 Sep 2017 11:08:08 -0600 Subject: [PATCH 14/20] Commented SSH slopes using PLM This PLM code leads to a segfault. I am commenting it out for now and using a simple second-order difference to get SSH gradients. --- config_src/mct_driver/ocn_comp_mct.F90 | 57 ++++++++++++++------------ 1 file changed, 31 insertions(+), 26 deletions(-) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index db82e95269..6297785322 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -1087,54 +1087,55 @@ subroutine ocn_export(ind, ocn_public, grid, o2x) ! Update halo of ssh so we can calculate gradients call pass_var(ssh, grid%domain) - ! d/dx ssh n = 0 do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec n = n+1 ! This is a simple second-order difference - ! o2x(ind%o2x_So_dhdx, n) = 0.5 * (ssh(i+1,j) + ssh(i-1,j)) * grid%IdxT(i,j) * grid%mask2dT(i,j) + o2x(ind%o2x_So_dhdx, n) = 0.5 * (ssh(i+1,j) + ssh(i-1,j)) * grid%IdxT(i,j) * grid%mask2dT(i,j) ! This is a PLM slope which might be less prone to the A-grid null mode - slp_L = (ssh(i,j) - ssh(i-1,j)) * grid%mask2dCu(I-1,j) +! slp_L = (ssh(I,j) - ssh(I-1,j)) * grid%mask2dCu(I-1,j) !if (grid%mask2dCu(I-1,j)==0.) slp_L = 0. - slp_R = (ssh(i+1,j) - ssh(i,j)) * grid%mask2dCu(I,j) +! slp_R = (ssh(I+1,j) - ssh(I,j)) * grid%mask2dCu(I,j) !if (grid%mask2dCu(I,j)==0.) slp_R = 0. - slp_C = 0.5 * (slp_L + slp_R) - if ( (slp_L * slp_R) > 0.0 ) then +! slp_C = 0.5 * (slp_L + slp_R) +! if ( (slp_L * slp_R) > 0.0 ) then ! This limits the slope so that the edge values are bounded by the ! two cell averages spanning the edge. - u_min = min( ssh(i-1,j), ssh(i,j), ssh(i+1,j) ) - u_max = max( ssh(i-1,j), ssh(i,j), ssh(i+1,j) ) - slope = sign( min( abs(slp_C), 2.*min( ssh(i,j) - u_min, u_max - ssh(i,j) ) ), slp_C ) - else +! u_min = min( ssh(i-1,j), ssh(i,j), ssh(i+1,j) ) +! u_max = max( ssh(i-1,j), ssh(i,j), ssh(i+1,j) ) +! slope = sign( min( abs(slp_C), 2.*min( ssh(i,j) - u_min, u_max - ssh(i,j) ) ), slp_C ) +! else ! Extrema in the mean values require a PCM reconstruction avoid generating ! larger extreme values. - slope = 0.0 - end if - o2x(ind%o2x_So_dhdx, n) = slope * grid%IdxT(i,j) * grid%mask2dT(i,j) +! slope = 0.0 +! end if +! o2x(ind%o2x_So_dhdx, n) = slope * grid%IdxT(i,j) * grid%mask2dT(i,j) end do; end do ! d/dy ssh + n = 0 do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec + n = n+1 ! This is a simple second-order difference - ! o2x(ind%o2x_So_dhdy, n) = 0.5 * (ssh(i,j+1) + ssh(i,j-1)) * grid%IdyT(i,j) * grid%mask2dT(i,j) + o2x(ind%o2x_So_dhdy, n) = 0.5 * (ssh(i,j+1) + ssh(i,j-1)) * grid%IdyT(i,j) * grid%mask2dT(i,j) ! This is a PLM slope which might be less prone to the A-grid null mode - slp_L = ssh(i,j) - ssh(i,j-1) - slp_R = ssh(i,j+1) - ssh(i,j) - slp_L=0.; slp_R=0. - slp_C = 0.5 * (slp_L + slp_R) - if ((slp_L * slp_R) > 0.0) then +! slp_L = ssh(i,J) - ssh(i,J-1) * grid%mask2dCv(i,J-1) +! slp_R = ssh(i,J+1) - ssh(i,J) * grid%mask2dCv(i,J) +! slp_C = 0.5 * (slp_L + slp_R) +! write(6,*)'slp_L, slp_R,i,j,slp_L*slp_R', slp_L, slp_R,i,j,slp_L*slp_R +! if ((slp_L * slp_R) > 0.0) then ! This limits the slope so that the edge values are bounded by the ! two cell averages spanning the edge. - u_min = min( ssh(i,j-1), ssh(i,j), ssh(i,j+1) ) - u_max = max( ssh(i,j-1), ssh(i,j), ssh(i,j+1) ) - slope = sign( min( abs(slp_C), 2.*min( ssh(i,j) - u_min, u_max - ssh(i,j) ) ), slp_C ) - else +! u_min = min( ssh(i,j-1), ssh(i,j), ssh(i,j+1) ) +! u_max = max( ssh(i,j-1), ssh(i,j), ssh(i,j+1) ) +! slope = sign( min( abs(slp_C), 2.*min( ssh(i,j) - u_min, u_max - ssh(i,j) ) ), slp_C ) +! else ! Extrema in the mean values require a PCM reconstruction avoid generating ! larger extreme values. - slope = 0.0 - end if - o2x(ind%o2x_So_dhdy, n) = slope * grid%IdyT(i,j) * grid%mask2dT(i,j) +! slope = 0.0 +! end if +! o2x(ind%o2x_So_dhdy, n) = slope * grid%IdyT(i,j) * grid%mask2dT(i,j) end do; end do end subroutine ocn_export @@ -1268,6 +1269,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) type(ESMF_timeInterval) :: ocn_cpl_interval !< The length of one ocean coupling interval integer :: year, month, day, hour, minute, seconds, seconds_n, seconds_d, rc logical :: write_restart_at_eod !< Controls if restart files must be written + logical :: debug=.false. type(time_type) :: time_start !< Start of coupled time interval to pass to MOM6 type(time_type) :: coupling_timestep !< Coupled time interval to pass to MOM6 character(len=128) :: err_msg !< Error message @@ -1316,6 +1318,9 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) call update_ocean_model(glb%ice_ocean_boundary, glb%ocn_state, glb%ocn_public, & time_start, coupling_timestep) + ! return export state to driver + call ocn_export(glb%ind, glb%ocn_public, glb%grid, o2x_o%rattr) + !--- write out intermediate restart file when needed. ! Check alarms for flag to write restart at end of day write_restart_at_eod = seq_timemgr_RestartAlarmIsOn(EClock) From eda526e3cfac0032ff849c65ace17ff92f94bf2d Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 11 Sep 2017 11:29:15 -0600 Subject: [PATCH 15/20] Revert temporary changes in ocean_model_MOM.F90 The option argument input_restart_file was a temporary change to allow restart file specification via mct_driver. Since the latter no longer depends on ocean_model_MOM.F90, this change can be reverted. PS: MOM-examples fails if input_restart_file is kept as an argument in ocean_model_init. --- config_src/coupled_driver/ocean_model_MOM.F90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 341fdfcf7c..e19057e77c 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -188,7 +188,7 @@ module ocean_model_mod !> ocean_model_init initializes the ocean model, including registering fields !! for restarts and reading restart files if appropriate. -subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, input_restart_file) +subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn) type(ocean_public_type), target, & intent(inout) :: Ocean_sfc !< A structure containing various !! publicly visible ocean surface properties after initialization, @@ -204,7 +204,6 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i !! in the calculation of additional gas or other !! tracer fluxes, and can be used to spawn related !! internal variables in the ice model. - character(len=*), optional, intent(in) :: input_restart_file !< If present, name of restart file to read ! This subroutine initializes both the ocean state and the ocean surface type. ! Because of the way that indicies and domains are handled, Ocean_sfc must have @@ -241,7 +240,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i OS%Time = Time_in call initialize_MOM(OS%Time, param_file, OS%dirs, OS%MOM_CSp, Time_in, & - offline_tracer_mode=offline_tracer_mode, input_restart_file=input_restart_file) + offline_tracer_mode=offline_tracer_mode) OS%grid => OS%MOM_CSp%G ; OS%GV => OS%MOM_CSp%GV OS%C_p = OS%MOM_CSp%tv%C_p OS%fluxes%C_p = OS%MOM_CSp%tv%C_p From 725df7d333c40a95ce9ca575ff826d8369052ebd Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 25 Sep 2017 11:34:36 -0600 Subject: [PATCH 16/20] Fix calculation of SSH slopes This PR fixes the calculation of SSH slopes using PLM. A bug (wrong sign) in the second-order scheme, which is not used at the moment, was also fixed. This PR closes https://github.com/NCAR/MOM6/issues/29 --- config_src/mct_driver/ocn_comp_mct.F90 | 85 +++++++++++++------------- 1 file changed, 41 insertions(+), 44 deletions(-) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index 9e77b183d6..c693a1d90e 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -479,17 +479,17 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) else ! hybrid or branch or continuos runs ! output path root call seq_infodata_GetData( glb%infodata, outPathRoot=restartpath ) - ! read pointer_filename - ! write name of restart file in the rpointer file + ! read name of restart file in the pointer file nu = shr_file_getUnit() - !if (is_root_pe()) then - restart_pointer_file = trim(glb%pointer_filename) - write(glb%stdout,*) 'Reading ocn pointer file: ',restart_pointer_file - open(nu, file=restart_pointer_file, form='formatted', status='unknown') - read(nu,'(a)') restartfile - close(nu) - !restartfile = trim(restartpath) // trim(restartfile) - write(glb%stdout,*) 'Reading ocn pointer file: ',trim(restartfile) + restart_pointer_file = trim(glb%pointer_filename) + write(glb%stdout,*) 'Reading ocn pointer file: ',restart_pointer_file + open(nu, file=restart_pointer_file, form='formatted', status='unknown') + read(nu,'(a)') restartfile + close(nu) + restartfile = trim(restartpath) // trim(restartfile) + write(6,*) 'Reading ocn pointer file: ', restartfile + write(6,*)'runtype, restartfile',runtype, restartfile + write(glb%stdout,*) 'Reading ocn pointer file: ',trim(restartfile) !endif call shr_file_freeUnit(nu) !restartfile = '/glade/scratch/gmarques/mom_test/run/mom_test.mom6.r.0001-01-01-21600.nc' @@ -504,13 +504,6 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in, input_restart_file=trim(restartfile)) endif - !if (debug .and. root_pe().eq.pe_here()) then - ! write(glb%stdout,*)'runtype, restartfile,time_init,time_in',runtype, restartfile,time_init,time_in - !endif - - ! Initialize the MOM6 model - !call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in) - ! Initialize ocn_state%state out of sight call ocean_model_init_sfc(glb%ocn_state, glb%ocn_public) @@ -1136,25 +1129,26 @@ subroutine ocn_export(ind, ocn_public, grid, o2x) do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec n = n+1 ! This is a simple second-order difference - o2x(ind%o2x_So_dhdx, n) = 0.5 * (ssh(i+1,j) + ssh(i-1,j)) * grid%IdxT(i,j) * grid%mask2dT(i,j) + !o2x(ind%o2x_So_dhdx, n) = 0.5 * (ssh(i+1,j) - ssh(i-1,j)) * grid%IdxT(i,j) * grid%mask2dT(i,j) ! This is a PLM slope which might be less prone to the A-grid null mode -! slp_L = (ssh(I,j) - ssh(I-1,j)) * grid%mask2dCu(I-1,j) - !if (grid%mask2dCu(I-1,j)==0.) slp_L = 0. -! slp_R = (ssh(I+1,j) - ssh(I,j)) * grid%mask2dCu(I,j) - !if (grid%mask2dCu(I,j)==0.) slp_R = 0. -! slp_C = 0.5 * (slp_L + slp_R) -! if ( (slp_L * slp_R) > 0.0 ) then + slp_L = (ssh(I,j) - ssh(I-1,j)) * grid%mask2dCu(I-1,j) + if (grid%mask2dCu(I-1,j)==0.) slp_L = 0. + slp_R = (ssh(I+1,j) - ssh(I,j)) * grid%mask2dCu(I,j) + if (grid%mask2dCu(I+1,j)==0.) slp_R = 0. + slp_C = 0.5 * (slp_L + slp_R) + if ( (slp_L * slp_R) > 0.0 ) then ! This limits the slope so that the edge values are bounded by the ! two cell averages spanning the edge. -! u_min = min( ssh(i-1,j), ssh(i,j), ssh(i+1,j) ) -! u_max = max( ssh(i-1,j), ssh(i,j), ssh(i+1,j) ) -! slope = sign( min( abs(slp_C), 2.*min( ssh(i,j) - u_min, u_max - ssh(i,j) ) ), slp_C ) -! else + u_min = min( ssh(i-1,j), ssh(i,j), ssh(i+1,j) ) + u_max = max( ssh(i-1,j), ssh(i,j), ssh(i+1,j) ) + slope = sign( min( abs(slp_C), 2.*min( ssh(i,j) - u_min, u_max - ssh(i,j) ) ), slp_C ) + else ! Extrema in the mean values require a PCM reconstruction avoid generating ! larger extreme values. -! slope = 0.0 -! end if -! o2x(ind%o2x_So_dhdx, n) = slope * grid%IdxT(i,j) * grid%mask2dT(i,j) + slope = 0.0 + end if + o2x(ind%o2x_So_dhdx, n) = slope * grid%IdxT(i,j) * grid%mask2dT(i,j) + if (grid%mask2dT(i,j)==0.) o2x(ind%o2x_So_dhdx, n) = 0.0 end do; end do ! d/dy ssh @@ -1162,24 +1156,27 @@ subroutine ocn_export(ind, ocn_public, grid, o2x) do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec n = n+1 ! This is a simple second-order difference - o2x(ind%o2x_So_dhdy, n) = 0.5 * (ssh(i,j+1) + ssh(i,j-1)) * grid%IdyT(i,j) * grid%mask2dT(i,j) + !o2x(ind%o2x_So_dhdy, n) = 0.5 * (ssh(i,j+1) - ssh(i,j-1)) * grid%IdyT(i,j) * grid%mask2dT(i,j) ! This is a PLM slope which might be less prone to the A-grid null mode -! slp_L = ssh(i,J) - ssh(i,J-1) * grid%mask2dCv(i,J-1) -! slp_R = ssh(i,J+1) - ssh(i,J) * grid%mask2dCv(i,J) -! slp_C = 0.5 * (slp_L + slp_R) -! write(6,*)'slp_L, slp_R,i,j,slp_L*slp_R', slp_L, slp_R,i,j,slp_L*slp_R -! if ((slp_L * slp_R) > 0.0) then + slp_L = ssh(i,J) - ssh(i,J-1) * grid%mask2dCv(i,J-1) + if (grid%mask2dCv(i,J-1)==0.) slp_L = 0. + slp_R = ssh(i,J+1) - ssh(i,J) * grid%mask2dCv(i,J) + if (grid%mask2dCv(i,J+1)==0.) slp_R = 0. + slp_C = 0.5 * (slp_L + slp_R) + !write(6,*)'slp_L, slp_R,i,j,slp_L*slp_R', slp_L, slp_R,i,j,slp_L*slp_R + if ((slp_L * slp_R) > 0.0) then ! This limits the slope so that the edge values are bounded by the ! two cell averages spanning the edge. -! u_min = min( ssh(i,j-1), ssh(i,j), ssh(i,j+1) ) -! u_max = max( ssh(i,j-1), ssh(i,j), ssh(i,j+1) ) -! slope = sign( min( abs(slp_C), 2.*min( ssh(i,j) - u_min, u_max - ssh(i,j) ) ), slp_C ) -! else + u_min = min( ssh(i,j-1), ssh(i,j), ssh(i,j+1) ) + u_max = max( ssh(i,j-1), ssh(i,j), ssh(i,j+1) ) + slope = sign( min( abs(slp_C), 2.*min( ssh(i,j) - u_min, u_max - ssh(i,j) ) ), slp_C ) + else ! Extrema in the mean values require a PCM reconstruction avoid generating ! larger extreme values. -! slope = 0.0 -! end if -! o2x(ind%o2x_So_dhdy, n) = slope * grid%IdyT(i,j) * grid%mask2dT(i,j) + slope = 0.0 + end if + o2x(ind%o2x_So_dhdy, n) = slope * grid%IdyT(i,j) * grid%mask2dT(i,j) + if (grid%mask2dT(i,j)==0.) o2x(ind%o2x_So_dhdy, n) = 0.0 end do; end do end subroutine ocn_export From 01bcdecb0f6d8a17b580fa7745975aad2dcf9f58 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 27 Sep 2017 09:57:37 -0600 Subject: [PATCH 17/20] Changes the path of restart files This commit eliminates the need to specify a restartpath. In all cases, cold or restart, the user must set restart_output_dir = './' in input.nml. --- config_src/mct_driver/ocn_comp_mct.F90 | 26 +++----------------------- 1 file changed, 3 insertions(+), 23 deletions(-) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index c693a1d90e..5cd929c901 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -381,17 +381,8 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) ocn_modelio_name = 'ocn_modelio.nml' // trim(inst_suffix) call shr_file_setIO(ocn_modelio_name,glb%stdout) - !if (debug) write(*,*) "stdout-1" - !if (debug) write(6,*) "stdout-2" - !if (debug) write(glb%stdout,*) "stdout-3" - ! set the shr log io unit number call shr_file_setLogUnit(glb%stdout) - - !if (debug) write(*,*) "stdout-4" - !if (debug) write(6,*) "stdout-5" - !if (debug) write(glb%stdout,*) "stdout-6" - end if call set_calendar_type(NOLEAP) !TODO: confirm this @@ -482,25 +473,14 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) ! read name of restart file in the pointer file nu = shr_file_getUnit() restart_pointer_file = trim(glb%pointer_filename) - write(glb%stdout,*) 'Reading ocn pointer file: ',restart_pointer_file + if (is_root_pe()) write(glb%stdout,*) 'Reading ocn pointer file: ',restart_pointer_file open(nu, file=restart_pointer_file, form='formatted', status='unknown') read(nu,'(a)') restartfile close(nu) - restartfile = trim(restartpath) // trim(restartfile) - write(6,*) 'Reading ocn pointer file: ', restartfile - write(6,*)'runtype, restartfile',runtype, restartfile - write(glb%stdout,*) 'Reading ocn pointer file: ',trim(restartfile) + !restartfile = trim(restartpath) // trim(restartfile) + if (is_root_pe()) write(glb%stdout,*) 'Reading restart file: ',trim(restartfile) !endif call shr_file_freeUnit(nu) - !restartfile = '/glade/scratch/gmarques/mom_test/run/mom_test.mom6.r.0001-01-01-21600.nc' - !call seq_infodata_GetData( glb%infodata, case_name=runid ) - !call ESMF_ClockGet(EClock, CurrTime=current_time, rc=rc) - !call ESMF_TimeGet(current_time, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) - !seconds = seconds + hour*3600 + minute*60 - !write(restartfile,'(A,".mom6.r.",I4.4,"-",I2.2,"-",I2.2,"-",I5.5,".nc")') trim(runid), year, month, day, seconds - if (debug .and. root_pe().eq.pe_here()) then - write(glb%stdout,*)'runtype, restartfile,time_init,time_in',runtype, restartfile !, time_init, time_in - endif call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in, input_restart_file=trim(restartfile)) endif From 8efe390e4abf57b6537e5c0f9fe88eb48cc22faa Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 2 Oct 2017 16:18:25 -0600 Subject: [PATCH 18/20] Makes ocn_comp_mct independent of MOM_surface_forcing This commit is the first step towards making ocn_comp_mct independent of MOM_surface_forcing. All subroutines from MOM_surface_forcing have been copied to ocn_comp_mct. Next steps: 1) get rid of ice_ocean_boundary_type 2) clean the module (remove unused code) --- config_src/mct_driver/ocn_comp_mct.F90 | 1126 +++++++++++++++++++++++- 1 file changed, 1118 insertions(+), 8 deletions(-) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index 5cd929c901..10b28cc462 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -29,18 +29,19 @@ module ocn_comp_mct shr_file_setLogUnit, shr_file_setLogLevel ! MOM6 modules +use MOM_coms, only : reproducing_sum +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end +use MOM_cpu_clock, only : CLOCK_SUBCOMPONENT use MOM, only: initialize_MOM, step_MOM, MOM_control_struct, MOM_end use MOM, only: calculate_surface_state, allocate_surface_state use MOM, only: finish_MOM_initialization, step_offline -use MOM_forcing_type, only: forcing, forcing_diagnostics, mech_forcing_diags, forcing_accumulate -use MOM_forcing_type, only: allocate_forcing_type -use MOM_surface_forcing, only: surface_forcing_init, convert_IOB_to_fluxes -use MOM_surface_forcing, only: ice_ocean_boundary_type, surface_forcing_CS -use MOM_surface_forcing, only: forcing_save_restart +use MOM_forcing_type, only: forcing, forcing_diags, register_forcing_type_diags +use MOM_forcing_type, only: allocate_forcing_type, deallocate_forcing_type +use MOM_forcing_type, only: mech_forcing_diags, forcing_accumulate, forcing_diagnostics use MOM_restart, only: save_restart use MOM_domains, only: MOM_infra_init, num_pes, root_pe, pe_here -use MOM_domains, only: pass_vector, BGRID_NE, CGRID_NE -use MOM_domains, only: pass_var, AGRID +use MOM_domains, only: pass_vector, BGRID_NE, CGRID_NE, To_All +use MOM_domains, only: pass_var, AGRID, fill_symmetric_edges use MOM_grid, only: ocean_grid_type, get_global_grid_size use MOM_verticalGrid, only: verticalGrid_type use MOM_variables, only: surface @@ -53,17 +54,28 @@ module ocn_comp_mct use MOM_get_input, only: Get_MOM_Input, directories use MOM_diag_mediator, only: diag_ctrl, enable_averaging, disable_averaging use MOM_diag_mediator, only: diag_mediator_close_registration, diag_mediator_end +use MOM_diag_mediator, only: safe_alloc_ptr use MOM_ice_shelf, only: initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS use MOM_ice_shelf, only: ice_shelf_end, ice_shelf_save_restart use MOM_sum_output, only: MOM_sum_output_init, sum_output_CS use MOM_sum_output, only: write_energy, accumulate_net_input use MOM_string_functions, only: uppercase -use MOM_constants, only: CELSIUS_KELVIN_OFFSET, hlf +use MOM_constants, only: CELSIUS_KELVIN_OFFSET, hlf, hlv use MOM_EOS, only: gsw_sp_from_sr, gsw_pt_from_ct +use user_revise_forcing, only : user_alter_forcing, user_revise_forcing_init +use user_revise_forcing, only : user_revise_forcing_CS +use MOM_restart, only : restart_init, MOM_restart_CS +use MOM_restart, only : restart_init_end, save_restart, restore_state +use data_override_mod, only : data_override_init, data_override +use MOM_io, only : slasher, write_version_number +use MOM_spatial_means, only : adjust_area_mean_to_zero ! FMS use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain +use time_interp_external_mod, only : init_external_field, time_interp_external +use time_interp_external_mod, only : time_interp_external_init +use fms_mod, only : read_data ! GFDL coupler use coupler_types_mod, only : coupler_1d_bc_type, coupler_2d_bc_type @@ -72,6 +84,9 @@ module ocn_comp_mct ! By default make data private implicit none; private + +#include + ! Public member functions public :: ocn_init_mct public :: ocn_run_mct @@ -192,6 +207,95 @@ module ocn_comp_mct ! for I/O using this surface data. end type ocean_public_type + !> This type contains pointers to the forcing fields which may be used to drive MOM. + !! All fluxes are positive downward. + type, public :: surface_forcing_CS ; private + integer :: wind_stagger ! AGRID, BGRID_NE, or CGRID_NE (integer values + ! from MOM_domains) to indicate the staggering of + ! the winds that are being provided in calls to + ! update_ocean_model. + logical :: use_temperature ! If true, temp and saln used as state variables + real :: wind_stress_multiplier !< A multiplier applied to incoming wind stress (nondim). + + ! smg: remove when have A=B code reconciled + logical :: bulkmixedlayer ! If true, model based on bulk mixed layer code + + real :: Rho0 ! Boussinesq reference density (kg/m^3) + real :: area_surf = -1.0 ! total ocean surface area (m^2) + real :: latent_heat_fusion ! latent heat of fusion (J/kg) + real :: latent_heat_vapor ! latent heat of vaporization (J/kg) + + real :: max_p_surf ! maximum surface pressure that can be + ! exerted by the atmosphere and floating sea-ice, + ! in Pa. This is needed because the FMS coupling + ! structure does not limit the water that can be + ! frozen out of the ocean and the ice-ocean heat + ! fluxes are treated explicitly. + logical :: use_limited_P_SSH ! If true, return the sea surface height with + ! the correction for the atmospheric (and sea-ice) + ! pressure limited by max_p_surf instead of the + ! full atmospheric pressure. The default is true. + + real :: gust_const ! constant unresolved background gustiness for ustar (Pa) + logical :: read_gust_2d ! If true, use a 2-dimensional gustiness supplied + ! from an input file. + real, pointer, dimension(:,:) :: & + TKE_tidal => NULL(), & ! turbulent kinetic energy introduced to the + ! bottom boundary layer by drag on the tidal flows, + ! in W m-2. + gust => NULL(), & ! spatially varying unresolved background + ! gustiness that contributes to ustar (Pa). + ! gust is used when read_gust_2d is true. + ustar_tidal => NULL() ! tidal contribution to the bottom friction velocity (m/s) + real :: cd_tides ! drag coefficient that applies to the tides (nondimensional) + real :: utide ! constant tidal velocity to use if read_tideamp + ! is false, in m s-1. + logical :: read_tideamp ! If true, spatially varying tidal amplitude read from a file. + + logical :: rigid_sea_ice ! If true, sea-ice exerts a rigidity that acts + ! to damp surface deflections (especially surface + ! gravity waves). The default is false. + real :: Kv_sea_ice ! viscosity in sea-ice that resists sheared vertical motions (m^2/s) + real :: density_sea_ice ! typical density of sea-ice (kg/m^3). The value is + ! only used to convert the ice pressure into + ! appropriate units for use with Kv_sea_ice. + real :: rigid_sea_ice_mass ! A mass per unit area of sea-ice beyond which + ! sea-ice viscosity becomes effective, in kg m-2, + ! typically of order 1000 kg m-2. + logical :: allow_flux_adjustments ! If true, use data_override to obtain flux adjustments + + real :: Flux_const ! piston velocity for surface restoring (m/s) + logical :: salt_restore_as_sflux ! If true, SSS restore as salt flux instead of water flux + logical :: adjust_net_srestore_to_zero ! adjust srestore to zero (for both salt_flux or vprec) + logical :: adjust_net_srestore_by_scaling ! adjust srestore w/o moving zero contour + logical :: adjust_net_fresh_water_to_zero ! adjust net surface fresh-water (w/ restoring) to zero + logical :: adjust_net_fresh_water_by_scaling ! adjust net surface fresh-water w/o moving zero contour + logical :: mask_srestore_under_ice ! If true, use an ice mask defined by frazil + ! criteria for salinity restoring. + real :: ice_salt_concentration ! salt concentration for sea ice (kg/kg) + logical :: mask_srestore_marginal_seas ! if true, then mask SSS restoring in marginal seas + real :: max_delta_srestore ! maximum delta salinity used for restoring + real :: max_delta_trestore ! maximum delta sst used for restoring + real, pointer, dimension(:,:) :: basin_mask => NULL() ! mask for SSS restoring + + type(diag_ctrl), pointer :: diag ! structure to regulate diagnostic output timing + character(len=200) :: inputdir ! directory where NetCDF input files are + character(len=200) :: salt_restore_file ! filename for salt restoring data + character(len=30) :: salt_restore_var_name ! name of surface salinity in salt_restore_file + character(len=200) :: temp_restore_file ! filename for sst restoring data + character(len=30) :: temp_restore_var_name ! name of surface temperature in temp_restore_file + + integer :: id_srestore = -1 ! id number for time_interp_external. + integer :: id_trestore = -1 ! id number for time_interp_external. + + ! Diagnostics handles + type(forcing_diags), public :: handles + + !### type(ctrl_forcing_CS), pointer :: ctrl_forcing_CSp => NULL() + type(MOM_restart_CS), pointer :: restart_CSp => NULL() + type(user_revise_forcing_CS), pointer :: urf_CS => NULL() + end type surface_forcing_CS + !> Contains information about the ocean state, although it is not necessary that !! this is implemented with all models. This type is private, and can therefore vary !! between different ocean models. @@ -243,6 +347,44 @@ module ocn_comp_mct type(sum_output_CS), pointer :: sum_output_CSp => NULL() end type ocean_state_type + !> Control structure corresponding to forcing, but with + !! the elements, units, and conventions that exactly conform to the use for + !! MOM-based coupled models. + type, public :: ice_ocean_boundary_type + real, pointer, dimension(:,:) :: u_flux =>NULL() ! i-direction wind stress (Pa) + real, pointer, dimension(:,:) :: v_flux =>NULL() ! j-direction wind stress (Pa) + real, pointer, dimension(:,:) :: t_flux =>NULL() ! sensible heat flux (W/m2) + real, pointer, dimension(:,:) :: q_flux =>NULL() ! specific humidity flux (kg/m2/s) + real, pointer, dimension(:,:) :: salt_flux =>NULL() ! salt flux (kg/m2/s) + real, pointer, dimension(:,:) :: lw_flux =>NULL() ! long wave radiation (W/m2) + real, pointer, dimension(:,:) :: sw_flux_vis_dir =>NULL() ! direct visible sw radiation (W/m2) + real, pointer, dimension(:,:) :: sw_flux_vis_dif =>NULL() ! diffuse visible sw radiation (W/m2) + real, pointer, dimension(:,:) :: sw_flux_nir_dir =>NULL() ! direct Near InfraRed sw radiation (W/m2) + real, pointer, dimension(:,:) :: sw_flux_nir_dif =>NULL() ! diffuse Near InfraRed sw radiation (W/m2) + real, pointer, dimension(:,:) :: latent =>NULL() ! latent heat flux (W/m2) + real, pointer, dimension(:,:) :: lprec =>NULL() ! mass flux of liquid precip (kg/m2/s) + real, pointer, dimension(:,:) :: fprec =>NULL() ! mass flux of frozen precip (kg/m2/s) + real, pointer, dimension(:,:) :: runoff =>NULL() ! mass flux of liquid runoff (kg/m2/s) + real, pointer, dimension(:,:) :: calving =>NULL() ! mass flux of frozen runoff (kg/m2/s) + real, pointer, dimension(:,:) :: ustar_berg =>NULL() ! frictional velocity beneath icebergs (m/s) + real, pointer, dimension(:,:) :: area_berg =>NULL() ! area covered by icebergs(m2/m2) + real, pointer, dimension(:,:) :: mass_berg =>NULL() ! mass of icebergs(kg/m2) + real, pointer, dimension(:,:) :: runoff_hflx =>NULL() ! heat content of liquid runoff (W/m2) + real, pointer, dimension(:,:) :: calving_hflx =>NULL() ! heat content of frozen runoff (W/m2) + real, pointer, dimension(:,:) :: p =>NULL() ! pressure of overlying ice and atmosphere + ! on ocean surface (Pa) + real, pointer, dimension(:,:) :: mi =>NULL() ! mass of ice (kg/m2) + integer :: xtype ! REGRID, REDIST or DIRECT + type(coupler_2d_bc_type) :: fluxes ! A structure that may contain an + ! array of named fields used for + ! passive tracer fluxes. + integer :: wind_stagger = -999 ! A flag indicating the spatial discretization of + ! wind stresses. This flag may be set by the + ! flux-exchange code, based on what the sea-ice + ! model is providing. Otherwise, the value from + ! the surface_forcing_CS is used. + end type ice_ocean_boundary_type + !> Control structure for this module type MCT_MOM_Data @@ -264,6 +406,8 @@ module ocn_comp_mct type(MCT_MOM_Data) :: glb !< global structure + integer :: id_clock_forcing + contains !> Initializes MOM6 @@ -920,6 +1064,319 @@ subroutine ocean_model_init_sfc(OS, Ocean_sfc) end subroutine ocean_model_init_sfc +!> Initialize surface forcing: get relevant parameters and allocate arrays +subroutine surface_forcing_init(Time, G, param_file, diag, CS, restore_salt, restore_temp) + type(time_type), intent(in) :: Time !< The current model time + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic output + type(surface_forcing_CS), pointer :: CS !< A pointer that is set to point to the + !! control structure for this module + logical, optional, intent(in) :: restore_salt, restore_temp !< If present and true, + !! temp/salt restoring will be applied + + ! local variables + real :: utide !< The RMS tidal velocity, in m s-1. + type(directories) :: dirs + logical :: new_sim, iceberg_flux_diags + type(time_type) :: Time_frc + character(len=200) :: TideAmp_file, gust_file, salt_file, temp_file ! Input file names. +! This include declares and sets the variable "version". +#include "version_variable.h" + character(len=40) :: mdl = "ocn_comp_mct" ! This module's name. + character(len=48) :: stagger + character(len=240) :: basin_file + integer :: i, j, isd, ied, jsd, jed + + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + + if (associated(CS)) then + call MOM_error(WARNING, "surface_forcing_init called with an associated "// & + "control structure.") + return + endif + allocate(CS) + + id_clock_forcing=cpu_clock_id('Ocean surface forcing', grain=CLOCK_SUBCOMPONENT) + call cpu_clock_begin(id_clock_forcing) + + CS%diag => diag + + call write_version_number (version) + ! Read all relevant parameters and write them to the model log. + call log_version(param_file, mdl, version, "") + + call get_param(param_file, mdl, "INPUTDIR", CS%inputdir, & + "The directory in which all input files are found.", & + default=".") + CS%inputdir = slasher(CS%inputdir) + call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", CS%use_temperature, & + "If true, Temperature and salinity are used as state \n"//& + "variables.", default=.true.) + call get_param(param_file, mdl, "RHO_0", CS%Rho0, & + "The mean ocean density used with BOUSSINESQ true to \n"//& + "calculate accelerations and the mass for conservation \n"//& + "properties, or with BOUSSINSEQ false to convert some \n"//& + "parameters from vertical units of m to kg m-2.", & + units="kg m-3", default=1035.0) + call get_param(param_file, mdl, "LATENT_HEAT_FUSION", CS%latent_heat_fusion, & + "The latent heat of fusion.", units="J/kg", default=hlf) + call get_param(param_file, mdl, "LATENT_HEAT_VAPORIZATION", CS%latent_heat_vapor, & + "The latent heat of fusion.", units="J/kg", default=hlv) + call get_param(param_file, mdl, "MAX_P_SURF", CS%max_p_surf, & + "The maximum surface pressure that can be exerted by the \n"//& + "atmosphere and floating sea-ice or ice shelves. This is \n"//& + "needed because the FMS coupling structure does not \n"//& + "limit the water that can be frozen out of the ocean and \n"//& + "the ice-ocean heat fluxes are treated explicitly. No \n"//& + "limit is applied if a negative value is used.", units="Pa", & + default=-1.0) + call get_param(param_file, mdl, "ADJUST_NET_SRESTORE_TO_ZERO", & + CS%adjust_net_srestore_to_zero, & + "If true, adjusts the salinity restoring seen to zero\n"//& + "whether restoring is via a salt flux or virtual precip.",& + default=restore_salt) + call get_param(param_file, mdl, "ADJUST_NET_SRESTORE_BY_SCALING", & + CS%adjust_net_srestore_by_scaling, & + "If true, adjustments to salt restoring to achieve zero net are\n"//& + "made by scaling values without moving the zero contour.",& + default=.false.) + call get_param(param_file, mdl, "ADJUST_NET_FRESH_WATER_TO_ZERO", & + CS%adjust_net_fresh_water_to_zero, & + "If true, adjusts the net fresh-water forcing seen \n"//& + "by the ocean (including restoring) to zero.", default=.false.) + call get_param(param_file, mdl, "ADJUST_NET_FRESH_WATER_BY_SCALING", & + CS%adjust_net_fresh_water_by_scaling, & + "If true, adjustments to net fresh water to achieve zero net are\n"//& + "made by scaling values without moving the zero contour.",& + default=.false.) + call get_param(param_file, mdl, "ICE_SALT_CONCENTRATION", & + CS%ice_salt_concentration, & + "The assumed sea-ice salinity needed to reverse engineer the \n"//& + "melt flux (or ice-ocean fresh-water flux).", & + units="kg/kg", default=0.005) + call get_param(param_file, mdl, "USE_LIMITED_PATM_SSH", CS%use_limited_P_SSH, & + "If true, return the sea surface height with the \n"//& + "correction for the atmospheric (and sea-ice) pressure \n"//& + "limited by max_p_surf instead of the full atmospheric \n"//& + "pressure.", default=.true.) + +! smg: should get_param call should be removed when have A=B code reconciled. +! this param is used to distinguish how to diagnose surface heat content from water. + call get_param(param_file, mdl, "BULKMIXEDLAYER", CS%bulkmixedlayer, & + default=CS%use_temperature,do_not_log=.true.) + + call get_param(param_file, mdl, "WIND_STAGGER", stagger, & + "A case-insensitive character string to indicate the \n"//& + "staggering of the input wind stress field. Valid \n"//& + "values are 'A', 'B', or 'C'.", default="C") + if (uppercase(stagger(1:1)) == 'A') then ; CS%wind_stagger = AGRID + elseif (uppercase(stagger(1:1)) == 'B') then ; CS%wind_stagger = BGRID_NE + elseif (uppercase(stagger(1:1)) == 'C') then ; CS%wind_stagger = CGRID_NE + else ; call MOM_error(FATAL,"surface_forcing_init: WIND_STAGGER = "// & + trim(stagger)//" is invalid.") ; endif + call get_param(param_file, mdl, "WIND_STRESS_MULTIPLIER", CS%wind_stress_multiplier, & + "A factor multiplying the wind-stress given to the ocean by the\n"//& + "coupler. This is used for testing and should be =1.0 for any\n"//& + "production runs.", default=1.0) + + if (restore_salt) then + call get_param(param_file, mdl, "FLUXCONST", CS%Flux_const, & + "The constant that relates the restoring surface fluxes \n"//& + "to the relative surface anomalies (akin to a piston \n"//& + "velocity). Note the non-MKS units.", units="m day-1", & + fail_if_missing=.true.) + call get_param(param_file, mdl, "SALT_RESTORE_FILE", CS%salt_restore_file, & + "A file in which to find the surface salinity to use for restoring.", & + default="salt_restore.nc") + call get_param(param_file, mdl, "SALT_RESTORE_VARIABLE", CS%salt_restore_var_name, & + "The name of the surface salinity variable to read from "//& + "SALT_RESTORE_FILE for restoring salinity.", & + default="salt") +! Convert CS%Flux_const from m day-1 to m s-1. + CS%Flux_const = CS%Flux_const / 86400.0 + + call get_param(param_file, mdl, "SRESTORE_AS_SFLUX", CS%salt_restore_as_sflux, & + "If true, the restoring of salinity is applied as a salt \n"//& + "flux instead of as a freshwater flux.", default=.false.) + call get_param(param_file, mdl, "MAX_DELTA_SRESTORE", CS%max_delta_srestore, & + "The maximum salinity difference used in restoring terms.", & + units="PSU or g kg-1", default=999.0) + call get_param(param_file, mdl, "MASK_SRESTORE_UNDER_ICE", & + CS%mask_srestore_under_ice, & + "If true, disables SSS restoring under sea-ice based on a frazil\n"//& + "criteria (SST<=Tf). Only used when RESTORE_SALINITY is True.", & + default=.false.) + call get_param(param_file, mdl, "MASK_SRESTORE_MARGINAL_SEAS", & + CS%mask_srestore_marginal_seas, & + "If true, disable SSS restoring in marginal seas. Only used when\n"//& + "RESTORE_SALINITY is True.", default=.false.) + call get_param(param_file, mdl, "BASIN_FILE", basin_file, & + "A file in which to find the basin masks, in variable 'basin'.", & + default="basin.nc") + basin_file = trim(CS%inputdir) // trim(basin_file) + call safe_alloc_ptr(CS%basin_mask,isd,ied,jsd,jed) ; CS%basin_mask(:,:) = 1.0 + if (CS%mask_srestore_marginal_seas) then + call read_data(basin_file,'basin',CS%basin_mask,domain=G%domain%mpp_domain,timelevel=1) + do j=jsd,jed ; do i=isd,ied + if (CS%basin_mask(i,j) >= 6.0) then ; CS%basin_mask(i,j) = 0.0 + else ; CS%basin_mask(i,j) = 1.0 ; endif + enddo ; enddo + endif + endif + + if (restore_temp) then + call get_param(param_file, mdl, "FLUXCONST", CS%Flux_const, & + "The constant that relates the restoring surface fluxes \n"//& + "to the relative surface anomalies (akin to a piston \n"//& + "velocity). Note the non-MKS units.", units="m day-1", & + fail_if_missing=.true.) + call get_param(param_file, mdl, "SST_RESTORE_FILE", CS%temp_restore_file, & + "A file in which to find the surface temperature to use for restoring.", & + default="temp_restore.nc") + call get_param(param_file, mdl, "SST_RESTORE_VARIABLE", CS%temp_restore_var_name, & + "The name of the surface temperature variable to read from "//& + "SST_RESTORE_FILE for restoring sst.", & + default="temp") +! Convert CS%Flux_const from m day-1 to m s-1. + CS%Flux_const = CS%Flux_const / 86400.0 + + call get_param(param_file, mdl, "MAX_DELTA_TRESTORE", CS%max_delta_trestore, & + "The maximum sst difference used in restoring terms.", & + units="degC ", default=999.0) + + endif + +! Optionally read tidal amplitude from input file (m s-1) on model grid. +! Otherwise use default tidal amplitude for bottom frictionally-generated +! dissipation. Default cd_tides is chosen to yield approx 1 TWatt of +! work done against tides globally using OSU tidal amplitude. + call get_param(param_file, mdl, "CD_TIDES", CS%cd_tides, & + "The drag coefficient that applies to the tides.", & + units="nondim", default=1.0e-4) + call get_param(param_file, mdl, "READ_TIDEAMP", CS%read_TIDEAMP, & + "If true, read a file (given by TIDEAMP_FILE) containing \n"//& + "the tidal amplitude with INT_TIDE_DISSIPATION.", default=.false.) + if (CS%read_TIDEAMP) then + call get_param(param_file, mdl, "TIDEAMP_FILE", TideAmp_file, & + "The path to the file containing the spatially varying \n"//& + "tidal amplitudes with INT_TIDE_DISSIPATION.", & + default="tideamp.nc") + CS%utide=0.0 + else + call get_param(param_file, mdl, "UTIDE", CS%utide, & + "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", & + units="m s-1", default=0.0) + endif + + call safe_alloc_ptr(CS%TKE_tidal,isd,ied,jsd,jed) + call safe_alloc_ptr(CS%ustar_tidal,isd,ied,jsd,jed) + + if (CS%read_TIDEAMP) then + TideAmp_file = trim(CS%inputdir) // trim(TideAmp_file) + call read_data(TideAmp_file,'tideamp',CS%TKE_tidal,domain=G%domain%mpp_domain,timelevel=1) + do j=jsd, jed; do i=isd, ied + utide = CS%TKE_tidal(i,j) + CS%TKE_tidal(i,j) = G%mask2dT(i,j)*CS%Rho0*CS%cd_tides*(utide*utide*utide) + CS%ustar_tidal(i,j)=sqrt(CS%cd_tides)*utide + enddo ; enddo + else + do j=jsd,jed; do i=isd,ied + utide=CS%utide + CS%TKE_tidal(i,j) = CS%Rho0*CS%cd_tides*(utide*utide*utide) + CS%ustar_tidal(i,j)=sqrt(CS%cd_tides)*utide + enddo ; enddo + endif + + call time_interp_external_init + +! Optionally read a x-y gustiness field in place of a global +! constant. + + call get_param(param_file, mdl, "READ_GUST_2D", CS%read_gust_2d, & + "If true, use a 2-dimensional gustiness supplied from \n"//& + "an input file", default=.false.) + call get_param(param_file, mdl, "GUST_CONST", CS%gust_const, & + "The background gustiness in the winds.", units="Pa", & + default=0.02) + if (CS%read_gust_2d) then + call get_param(param_file, mdl, "GUST_2D_FILE", gust_file, & + "The file in which the wind gustiness is found in \n"//& + "variable gustiness.") + + call safe_alloc_ptr(CS%gust,isd,ied,jsd,jed) + gust_file = trim(CS%inputdir) // trim(gust_file) + call read_data(gust_file,'gustiness',CS%gust,domain=G%domain%mpp_domain, & + timelevel=1) ! units should be Pa + endif + +! See whether sufficiently thick sea ice should be treated as rigid. + call get_param(param_file, mdl, "USE_RIGID_SEA_ICE", CS%rigid_sea_ice, & + "If true, sea-ice is rigid enough to exert a \n"//& + "nonhydrostatic pressure that resist vertical motion.", & + default=.false.) + if (CS%rigid_sea_ice) then + call get_param(param_file, mdl, "SEA_ICE_MEAN_DENSITY", CS%density_sea_ice, & + "A typical density of sea ice, used with the kinematic \n"//& + "viscosity, when USE_RIGID_SEA_ICE is true.", units="kg m-3", & + default=900.0) + call get_param(param_file, mdl, "SEA_ICE_VISCOSITY", CS%Kv_sea_ice, & + "The kinematic viscosity of sufficiently thick sea ice \n"//& + "for use in calculating the rigidity of sea ice.", & + units="m2 s-1", default=1.0e9) + call get_param(param_file, mdl, "SEA_ICE_RIGID_MASS", CS%rigid_sea_ice_mass, & + "The mass of sea-ice per unit area at which the sea-ice \n"//& + "starts to exhibit rigidity", units="kg m-2", default=1000.0) + endif + + call get_param(param_file, mdl, "ALLOW_ICEBERG_FLUX_DIAGNOSTICS", iceberg_flux_diags, & + "If true, makes available diagnostics of fluxes from icebergs\n"//& + "as seen by MOM6.", default=.false.) + call register_forcing_type_diags(Time, diag, CS%use_temperature, CS%handles, & + use_berg_fluxes=iceberg_flux_diags) + + call get_param(param_file, mdl, "ALLOW_FLUX_ADJUSTMENTS", CS%allow_flux_adjustments, & + "If true, allows flux adjustments to specified via the \n"//& + "data_table using the component name 'OCN'.", default=.false.) + if (CS%allow_flux_adjustments) then + call data_override_init(Ocean_domain_in=G%Domain%mpp_domain) + endif + + if (present(restore_salt)) then ; if (restore_salt) then + salt_file = trim(CS%inputdir) // trim(CS%salt_restore_file) + CS%id_srestore = init_external_field(salt_file, CS%salt_restore_var_name, domain=G%Domain%mpp_domain) + endif ; endif + + if (present(restore_temp)) then ; if (restore_temp) then + temp_file = trim(CS%inputdir) // trim(CS%temp_restore_file) + CS%id_trestore = init_external_field(temp_file, CS%temp_restore_var_name, domain=G%Domain%mpp_domain) + endif ; endif + + ! Set up any restart fields associated with the forcing. + call restart_init(param_file, CS%restart_CSp, "MOM_forcing.res") +!### call register_ctrl_forcing_restarts(G, param_file, CS%ctrl_forcing_CSp, & +!### CS%restart_CSp) + call restart_init_end(CS%restart_CSp) + + if (associated(CS%restart_CSp)) then + call Get_MOM_Input(dirs=dirs) + + new_sim = .false. + if ((dirs%input_filename(1:1) == 'n') .and. & + (LEN_TRIM(dirs%input_filename) == 1)) new_sim = .true. + if (.not.new_sim) then + call restore_state(dirs%input_filename, dirs%restart_input_dir, Time_frc, & + G, CS%restart_CSp) + endif + endif + +!### call controlled_forcing_init(Time, G, param_file, diag, CS%ctrl_forcing_CSp) + + call user_revise_forcing_init(param_file, CS%urf_CS) + + call cpu_clock_end(id_clock_forcing) +end subroutine surface_forcing_init + subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, maskmap, & gas_fields_ocn) type(domain2D), intent(in) :: input_domain @@ -1398,6 +1855,28 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) end subroutine ocn_run_mct +!> Saves restart fields associated with the forcing +subroutine forcing_save_restart(CS, G, Time, directory, time_stamped, & + filename_suffix) + type(surface_forcing_CS), pointer :: CS !< pointer to the control structure + !! returned by a previous call to + !! surface_forcing_init + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + type(time_type), intent(in) :: Time !< model time at this call + character(len=*), intent(in) :: directory !< optional directory into which + !! to write these restart files + logical, optional, intent(in) :: time_stamped !< If true, the restart file + !! names include a unique time + !! stamp + character(len=*), optional, intent(in) :: filename_suffix !< optional suffix + !! (e.g., a time-stamp) to append to the + !! restart file names + if (.not.associated(CS)) return + if (.not.associated(CS%restart_CSp)) return + call save_restart(directory, Time, G, CS%restart_CSp, time_stamped) + +end subroutine forcing_save_restart + !> Updates the ocean model fields. This code wraps the call to step_MOM with MOM6's call. !! It uses the forcing in Ice_ocean_boundary to advance the ocean model's state from the !! input value of Ocean_state (which must be for time time_start_update) for a time interval @@ -1557,6 +2036,637 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & call callTree_leave("update_ocean_model()") end subroutine update_ocean_model +subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, state, restore_salt, restore_temp) + type(ice_ocean_boundary_type), intent(in), target :: IOB + type(forcing), intent(inout) :: fluxes + integer, dimension(4), intent(in) :: index_bounds + type(time_type), intent(in) :: Time + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + type(surface_forcing_CS), pointer :: CS + type(surface), intent(in) :: state + logical, optional, intent(in) :: restore_salt, restore_temp + +! This subroutine translates the Ice_ocean_boundary_type into a +! MOM forcing type, including changes of units, sign conventions, +! and puting the fields into arrays with MOM-standard halos. + +! Arguments: +! IOB ice-ocean boundary type w/ fluxes to drive ocean in a coupled model +! (out) fluxes - A structure containing pointers to any possible +! forcing fields. Unused fields have NULL ptrs. +! (in) index_bounds - the i- and j- size of the arrays in IOB. +! (in) Time - The time of the fluxes, used for interpolating the salinity +! to the right time, when it is being restored. +! (in) G - The ocean's grid structure. +! (in) CS - A pointer to the control structure returned by a previous +! call to surface_forcing_init. +! (in) state - A structure containing fields that describe the +! surface state of the ocean. +! (in) restore_salt - if true, salinity is restored to a target value. +! (in) restore_temp - if true, temperature is restored to a target value. + + real, dimension(SZIB_(G),SZJB_(G)) :: & + taux_at_q, & ! Zonal wind stresses at q points (Pa) + tauy_at_q ! Meridional wind stresses at q points (Pa) + + real, dimension(SZI_(G),SZJ_(G)) :: & + taux_at_h, & ! Zonal wind stresses at h points (Pa) + tauy_at_h, & ! Meridional wind stresses at h points (Pa) + data_restore, & ! The surface value toward which to restore (g/kg or degC) + SST_anom, & ! Instantaneous sea surface temperature anomalies from a target value (deg C) + SSS_anom, & ! Instantaneous sea surface salinity anomalies from a target value (g/kg) + SSS_mean, & ! A (mean?) salinity about which to normalize local salinity + ! anomalies when calculating restorative precipitation anomalies (g/kg) + PmE_adj, & ! The adjustment to PminusE that will cause the salinity + ! to be restored toward its target value (kg/(m^2 * s)) + net_FW, & ! The area integrated net freshwater flux into the ocean (kg/s) + net_FW2, & ! The area integrated net freshwater flux into the ocean (kg/s) + work_sum, & ! A 2-d array that is used as the work space for a global + ! sum, used with units of m2 or (kg/s) + open_ocn_mask ! a binary field indicating where ice is present based on frazil criteria + + real :: gustiness ! unresolved gustiness that contributes to ustar (Pa) + real :: Irho0 ! inverse of the mean density in (m^3/kg) + real :: taux2, tauy2 ! squared wind stresses (Pa^2) + real :: tau_mag ! magnitude of the wind stress (Pa) + real :: I_GEarth ! 1.0 / G%G_Earth (s^2/m) + real :: Kv_rho_ice ! (CS%kv_sea_ice / CS%density_sea_ice) ( m^5/(s*kg) ) + real :: mass_ice ! mass of sea ice at a face (kg/m^2) + real :: mass_eff ! effective mass of sea ice for rigidity (kg/m^2) + + integer :: wind_stagger ! AGRID, BGRID_NE, or CGRID_NE (integers from MOM_domains) + integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0 + integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isr, ier, jsr, jer + integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd + + logical :: restore_salinity ! local copy of the argument restore_salt, if it + ! is present, or false (no restoring) otherwise. + logical :: restore_sst ! local copy of the argument restore_temp, if it + ! is present, or false (no restoring) otherwise. + real :: delta_sss ! temporary storage for sss diff from restoring value + real :: delta_sst ! temporary storage for sst diff from restoring value + + real :: C_p ! heat capacity of seawater ( J/(K kg) ) + + call cpu_clock_begin(id_clock_forcing) + + isc_bnd = index_bounds(1) ; iec_bnd = index_bounds(2) + jsc_bnd = index_bounds(3) ; jec_bnd = index_bounds(4) + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + isr = is-isd+1 ; ier = ie-isd+1 ; jsr = js-jsd+1 ; jer = je-jsd+1 + + C_p = fluxes%C_p + Irho0 = 1.0/CS%Rho0 + open_ocn_mask(:,:) = 1.0 + pme_adj(:,:) = 0.0 + fluxes%vPrecGlobalAdj = 0.0 + fluxes%vPrecGlobalScl = 0.0 + fluxes%saltFluxGlobalAdj = 0.0 + fluxes%saltFluxGlobalScl = 0.0 + fluxes%netFWGlobalAdj = 0.0 + fluxes%netFWGlobalScl = 0.0 + + restore_salinity = .false. + if (present(restore_salt)) restore_salinity = restore_salt + restore_sst = .false. + if (present(restore_temp)) restore_sst = restore_temp + + ! allocation and initialization if this is the first time that this + ! flux type has been used. + if (fluxes%dt_buoy_accum < 0) then + call allocate_forcing_type(G, fluxes, stress=.true., ustar=.true., & + water=.true., heat=.true.) + + call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%sw_nir_dir,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%sw_nir_dif,isd,ied,jsd,jed) + + call safe_alloc_ptr(fluxes%p_surf,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%p_surf_full,isd,ied,jsd,jed) + + call safe_alloc_ptr(fluxes%salt_flux,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%salt_flux_in,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%salt_flux_added,isd,ied,jsd,jed) + + call safe_alloc_ptr(fluxes%TKE_tidal,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%ustar_tidal,isd,ied,jsd,jed) + + if (CS%allow_flux_adjustments) then + call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%salt_flux_added,isd,ied,jsd,jed) + endif + + do j=js-2,je+2 ; do i=is-2,ie+2 + fluxes%TKE_tidal(i,j) = CS%TKE_tidal(i,j) + fluxes%ustar_tidal(i,j) = CS%ustar_tidal(i,j) + enddo; enddo + + if (CS%rigid_sea_ice) then + call safe_alloc_ptr(fluxes%rigidity_ice_u,IsdB,IedB,jsd,jed) + call safe_alloc_ptr(fluxes%rigidity_ice_v,isd,ied,JsdB,JedB) + endif + + if (restore_temp) call safe_alloc_ptr(fluxes%heat_added,isd,ied,jsd,jed) + + fluxes%dt_buoy_accum = 0.0 + endif ! endif for allocation and initialization + + if ((.not.coupler_type_initialized(fluxes%tr_fluxes)) .and. & + coupler_type_initialized(IOB%fluxes)) & + call coupler_type_spawn(IOB%fluxes, fluxes%tr_fluxes, & + (/is,is,ie,ie/), (/js,js,je,je/)) + ! It might prove valuable to use the same array extents as the rest of the + ! ocean model, rather than using haloless arrays, in which case the last line + ! would be: ( (/isd,is,ie,ied/), (/jsd,js,je,jed/)) + + + if (CS%allow_flux_adjustments) then + fluxes%heat_added(:,:)=0.0 + fluxes%salt_flux_added(:,:)=0.0 + endif + + ! allocation and initialization on first call to this routine + if (CS%area_surf < 0.0) then + do j=js,je ; do i=is,ie + work_sum(i,j) = G%areaT(i,j) * G%mask2dT(i,j) + enddo ; enddo + CS%area_surf = reproducing_sum(work_sum, isr, ier, jsr, jer) + endif ! endif for allocation and initialization + + do j=js,je ; do i=is,ie + fluxes%salt_flux(i,j) = 0.0 + fluxes%vprec(i,j) = 0.0 + enddo ; enddo + + ! Salinity restoring logic + if (restore_salinity) then + call time_interp_external(CS%id_srestore,Time,data_restore) + ! open_ocn_mask indicates where to restore salinity (1 means restore, 0 does not) + open_ocn_mask(:,:) = 1.0 + if (CS%mask_srestore_under_ice) then ! Do not restore under sea-ice + do j=js,je ; do i=is,ie + if (state%SST(i,j) .le. -0.0539*state%SSS(i,j)) open_ocn_mask(i,j)=0.0 + enddo; enddo + endif + if (CS%salt_restore_as_sflux) then + do j=js,je ; do i=is,ie + delta_sss = data_restore(i,j)- state%SSS(i,j) + delta_sss = sign(1.0,delta_sss)*min(abs(delta_sss),CS%max_delta_srestore) + fluxes%salt_flux(i,j) = 1.e-3*G%mask2dT(i,j) * (CS%Rho0*CS%Flux_const)* & + (CS%basin_mask(i,j)*open_ocn_mask(i,j)) *delta_sss ! kg Salt m-2 s-1 + enddo; enddo + if (CS%adjust_net_srestore_to_zero) then + if (CS%adjust_net_srestore_by_scaling) then + call adjust_area_mean_to_zero(fluxes%salt_flux, G, fluxes%saltFluxGlobalScl) + fluxes%saltFluxGlobalAdj = 0. + else + work_sum(is:ie,js:je) = G%areaT(is:ie,js:je)*fluxes%salt_flux(is:ie,js:je) + fluxes%saltFluxGlobalAdj = reproducing_sum(work_sum(:,:), isr,ier, jsr,jer)/CS%area_surf + fluxes%salt_flux(is:ie,js:je) = fluxes%salt_flux(is:ie,js:je) - fluxes%saltFluxGlobalAdj + endif + endif + fluxes%salt_flux_added(is:ie,js:je) = fluxes%salt_flux(is:ie,js:je) ! Diagnostic + else + do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) > 0.5) then + delta_sss = state%SSS(i,j) - data_restore(i,j) + delta_sss = sign(1.0,delta_sss)*min(abs(delta_sss),CS%max_delta_srestore) + fluxes%vprec(i,j) = (CS%basin_mask(i,j)*open_ocn_mask(i,j))* & + (CS%Rho0*CS%Flux_const) * & + delta_sss / (0.5*(state%SSS(i,j) + data_restore(i,j))) + endif + enddo; enddo + if (CS%adjust_net_srestore_to_zero) then + if (CS%adjust_net_srestore_by_scaling) then + call adjust_area_mean_to_zero(fluxes%vprec, G, fluxes%vPrecGlobalScl) + fluxes%vPrecGlobalAdj = 0. + else + work_sum(is:ie,js:je) = G%areaT(is:ie,js:je)*fluxes%vprec(is:ie,js:je) + fluxes%vPrecGlobalAdj = reproducing_sum(work_sum(:,:), isr, ier, jsr, jer) / CS%area_surf + do j=js,je ; do i=is,ie + fluxes%vprec(i,j) = ( fluxes%vprec(i,j) - fluxes%vPrecGlobalAdj ) * G%mask2dT(i,j) + enddo; enddo + endif + endif + endif + endif + + ! SST restoring logic + if (restore_sst) then + call time_interp_external(CS%id_trestore,Time,data_restore) + do j=js,je ; do i=is,ie + delta_sst = data_restore(i,j)- state%SST(i,j) + delta_sst = sign(1.0,delta_sst)*min(abs(delta_sst),CS%max_delta_trestore) + fluxes%heat_added(i,j) = G%mask2dT(i,j) * (CS%Rho0*fluxes%C_p) * delta_sst * CS%Flux_const ! W m-2 + enddo; enddo + endif + + wind_stagger = CS%wind_stagger + if ((IOB%wind_stagger == AGRID) .or. (IOB%wind_stagger == BGRID_NE) .or. & + (IOB%wind_stagger == CGRID_NE)) wind_stagger = IOB%wind_stagger + if (wind_stagger == BGRID_NE) then + ! This is necessary to fill in the halo points. + taux_at_q(:,:) = 0.0 ; tauy_at_q(:,:) = 0.0 + endif + if (wind_stagger == AGRID) then + ! This is necessary to fill in the halo points. + taux_at_h(:,:) = 0.0 ; tauy_at_h(:,:) = 0.0 + endif + + + ! obtain fluxes from IOB; note the staggering of indices + i0 = is - isc_bnd ; j0 = js - jsc_bnd + do j=js,je ; do i=is,ie + + if (wind_stagger == BGRID_NE) then + if (ASSOCIATED(IOB%u_flux)) taux_at_q(I,J) = IOB%u_flux(i-i0,j-j0) * CS%wind_stress_multiplier + if (ASSOCIATED(IOB%v_flux)) tauy_at_q(I,J) = IOB%v_flux(i-i0,j-j0) * CS%wind_stress_multiplier + elseif (wind_stagger == AGRID) then + if (ASSOCIATED(IOB%u_flux)) taux_at_h(i,j) = IOB%u_flux(i-i0,j-j0) * CS%wind_stress_multiplier + if (ASSOCIATED(IOB%v_flux)) tauy_at_h(i,j) = IOB%v_flux(i-i0,j-j0) * CS%wind_stress_multiplier + else ! C-grid wind stresses. + if (ASSOCIATED(IOB%u_flux)) fluxes%taux(I,j) = IOB%u_flux(i-i0,j-j0) * CS%wind_stress_multiplier + if (ASSOCIATED(IOB%v_flux)) fluxes%tauy(i,J) = IOB%v_flux(i-i0,j-j0) * CS%wind_stress_multiplier + endif + + if (ASSOCIATED(IOB%lprec)) & + fluxes%lprec(i,j) = IOB%lprec(i-i0,j-j0) * G%mask2dT(i,j) + + if (ASSOCIATED(IOB%fprec)) & + fluxes%fprec(i,j) = IOB%fprec(i-i0,j-j0) * G%mask2dT(i,j) + + if (ASSOCIATED(IOB%q_flux)) & + fluxes%evap(i,j) = - IOB%q_flux(i-i0,j-j0) * G%mask2dT(i,j) + + if (ASSOCIATED(IOB%runoff)) & + fluxes%lrunoff(i,j) = IOB%runoff(i-i0,j-j0) * G%mask2dT(i,j) + + if (ASSOCIATED(IOB%calving)) & + fluxes%frunoff(i,j) = IOB%calving(i-i0,j-j0) * G%mask2dT(i,j) + + if (((ASSOCIATED(IOB%ustar_berg) .and. (.not. ASSOCIATED(fluxes%ustar_berg))) & + .or. (ASSOCIATED(IOB%area_berg) .and. (.not. ASSOCIATED(fluxes%area_berg)))) & + .or. (ASSOCIATED(IOB%mass_berg) .and. (.not. ASSOCIATED(fluxes%mass_berg)))) & + call allocate_forcing_type(G, fluxes, iceberg=.true.) + + if (ASSOCIATED(IOB%ustar_berg)) & + fluxes%ustar_berg(i,j) = IOB%ustar_berg(i-i0,j-j0) * G%mask2dT(i,j) + + if (ASSOCIATED(IOB%area_berg)) & + fluxes%area_berg(i,j) = IOB%area_berg(i-i0,j-j0) * G%mask2dT(i,j) + + if (ASSOCIATED(IOB%mass_berg)) & + fluxes%mass_berg(i,j) = IOB%mass_berg(i-i0,j-j0) * G%mask2dT(i,j) + + if (ASSOCIATED(IOB%runoff_hflx)) & + fluxes%heat_content_lrunoff(i,j) = IOB%runoff_hflx(i-i0,j-j0) * G%mask2dT(i,j) + + if (ASSOCIATED(IOB%calving_hflx)) & + fluxes%heat_content_frunoff(i,j) = IOB%calving_hflx(i-i0,j-j0) * G%mask2dT(i,j) + + if (ASSOCIATED(IOB%lw_flux)) & + fluxes%LW(i,j) = IOB%lw_flux(i-i0,j-j0) * G%mask2dT(i,j) + + if (ASSOCIATED(IOB%t_flux)) & + fluxes%sens(i,j) = - IOB%t_flux(i-i0,j-j0) * G%mask2dT(i,j) + + fluxes%latent(i,j) = 0.0 + if (ASSOCIATED(IOB%latent)) then + fluxes%latent(i,j) = IOB%latent(i-i0,j-j0) + else + if (ASSOCIATED(IOB%fprec)) then + fluxes%latent(i,j) = fluxes%latent(i,j) - IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion + fluxes%latent_fprec_diag(i,j) = -G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion + endif + if (ASSOCIATED(IOB%calving)) then + fluxes%latent(i,j) = fluxes%latent(i,j) - IOB%calving(i-i0,j-j0)*CS%latent_heat_fusion + fluxes%latent_frunoff_diag(i,j) = -G%mask2dT(i,j) * IOB%calving(i-i0,j-j0)*CS%latent_heat_fusion + endif + if (ASSOCIATED(IOB%q_flux)) then + fluxes%latent(i,j) = fluxes%latent(i,j) - IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor + fluxes%latent_evap_diag(i,j) = -G%mask2dT(i,j) * IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor + endif + endif + + fluxes%latent(i,j) = G%mask2dT(i,j) * fluxes%latent(i,j) + + if (ASSOCIATED(IOB%sw_flux_vis_dir)) & + fluxes%sw_vis_dir(i,j) = G%mask2dT(i,j) * IOB%sw_flux_vis_dir(i-i0,j-j0) + if (ASSOCIATED(IOB%sw_flux_vis_dif)) & + fluxes%sw_vis_dif(i,j) = G%mask2dT(i,j) * IOB%sw_flux_vis_dif(i-i0,j-j0) + if (ASSOCIATED(IOB%sw_flux_nir_dir)) & + fluxes%sw_nir_dir(i,j) = G%mask2dT(i,j) * IOB%sw_flux_nir_dir(i-i0,j-j0) + if (ASSOCIATED(IOB%sw_flux_nir_dif)) & + fluxes%sw_nir_dif(i,j) = G%mask2dT(i,j) * IOB%sw_flux_nir_dif(i-i0,j-j0) + fluxes%sw(i,j) = fluxes%sw_vis_dir(i,j) + fluxes%sw_vis_dif(i,j) + & + fluxes%sw_nir_dir(i,j) + fluxes%sw_nir_dif(i,j) + + enddo ; enddo + + ! more salt restoring logic + if (ASSOCIATED(IOB%salt_flux)) then + do j=js,je ; do i=is,ie + fluxes%salt_flux(i,j) = G%mask2dT(i,j)*(fluxes%salt_flux(i,j) - IOB%salt_flux(i-i0,j-j0)) + fluxes%salt_flux_in(i,j) = G%mask2dT(i,j)*( -IOB%salt_flux(i-i0,j-j0) ) + enddo ; enddo + endif + +!### if (associated(CS%ctrl_forcing_CSp)) then +!### do j=js,je ; do i=is,ie +!### SST_anom(i,j) = state%SST(i,j) - CS%T_Restore(i,j) +!### SSS_anom(i,j) = state%SSS(i,j) - CS%S_Restore(i,j) +!### SSS_mean(i,j) = 0.5*(state%SSS(i,j) + CS%S_Restore(i,j)) +!### enddo ; enddo +!### call apply_ctrl_forcing(SST_anom, SSS_anom, SSS_mean, fluxes%heat_restore, & +!### fluxes%vprec, day, dt, G, CS%ctrl_forcing_CSp) +!### endif + + ! adjust the NET fresh-water flux to zero, if flagged + if (CS%adjust_net_fresh_water_to_zero) then + do j=js,je ; do i=is,ie + net_FW(i,j) = (((fluxes%lprec(i,j) + fluxes%fprec(i,j)) + & + (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j))) + & + (fluxes%evap(i,j) + fluxes%vprec(i,j)) ) * G%areaT(i,j) + ! The following contribution appears to be calculating the volume flux of sea-ice + ! melt. This calculation is clearly WRONG if either sea-ice has variable + ! salinity or the sea-ice is completely fresh. + ! Bob thinks this is trying ensure the net fresh-water of the ocean + sea-ice system + ! is constant. + ! To do this correctly we will need a sea-ice melt field added to IOB. -AJA + if (ASSOCIATED(IOB%salt_flux) .and. (CS%ice_salt_concentration>0.0)) & + net_FW(i,j) = net_FW(i,j) - G%areaT(i,j) * & + (IOB%salt_flux(i-i0,j-j0) / CS%ice_salt_concentration) + net_FW2(i,j) = net_FW(i,j) + enddo ; enddo + + if (CS%adjust_net_fresh_water_by_scaling) then + call adjust_area_mean_to_zero(net_FW2, G, fluxes%netFWGlobalScl) + do j=js,je ; do i=is,ie + fluxes%vprec(i,j) = fluxes%vprec(i,j) + (net_FW2(i,j) - net_FW(i,j)) * G%mask2dT(i,j) + enddo; enddo + else + fluxes%netFWGlobalAdj = reproducing_sum(net_FW(:,:), isr, ier, jsr, jer) / CS%area_surf + do j=js,je ; do i=is,ie + fluxes%vprec(i,j) = ( fluxes%vprec(i,j) - fluxes%netFWGlobalAdj ) * G%mask2dT(i,j) + enddo; enddo + endif + + endif + + ! applied surface pressure from atmosphere and cryosphere + if (ASSOCIATED(IOB%p)) then + if (CS%max_p_surf >= 0.0) then + do j=js,je ; do i=is,ie + fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0) + fluxes%p_surf(i,j) = MIN(fluxes%p_surf_full(i,j),CS%max_p_surf) + enddo ; enddo + else + do j=js,je ; do i=is,ie + fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0) + fluxes%p_surf(i,j) = fluxes%p_surf_full(i,j) + enddo ; enddo + endif + if (CS%use_limited_P_SSH) then + fluxes%p_surf_SSH => fluxes%p_surf + else + fluxes%p_surf_SSH => fluxes%p_surf_full + endif + endif + + ! surface momentum stress related fields as function of staggering + if (wind_stagger == BGRID_NE) then + if (G%symmetric) & + call fill_symmetric_edges(taux_at_q, tauy_at_q, G%Domain, stagger=BGRID_NE) + call pass_vector(taux_at_q, tauy_at_q, G%Domain, stagger=BGRID_NE) + + do j=js,je ; do I=Isq,Ieq + fluxes%taux(I,j) = 0.0 + if ((G%mask2dBu(I,J) + G%mask2dBu(I,J-1)) > 0) & + fluxes%taux(I,j) = (G%mask2dBu(I,J)*taux_at_q(I,J) + & + G%mask2dBu(I,J-1)*taux_at_q(I,J-1)) / & + (G%mask2dBu(I,J) + G%mask2dBu(I,J-1)) + enddo ; enddo + + do J=Jsq,Jeq ; do i=is,ie + fluxes%tauy(i,J) = 0.0 + if ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J)) > 0) & + fluxes%tauy(i,J) = (G%mask2dBu(I,J)*tauy_at_q(I,J) + & + G%mask2dBu(I-1,J)*tauy_at_q(I-1,J)) / & + (G%mask2dBu(I,J) + G%mask2dBu(I-1,J)) + enddo ; enddo + + ! ustar is required for the bulk mixed layer formulation. The background value + ! of 0.02 Pa is a relatively small value intended to give reasonable behavior + ! in regions of very weak winds. + + do j=js,je ; do i=is,ie + tau_mag = 0.0 ; gustiness = CS%gust_const + if (((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + & + (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) > 0) then + tau_mag = sqrt(((G%mask2dBu(I,J)*(taux_at_q(I,J)**2 + tauy_at_q(I,J)**2) + & + G%mask2dBu(I-1,J-1)*(taux_at_q(I-1,J-1)**2 + tauy_at_q(I-1,J-1)**2)) + & + (G%mask2dBu(I,J-1)*(taux_at_q(I,J-1)**2 + tauy_at_q(I,J-1)**2) + & + G%mask2dBu(I-1,J)*(taux_at_q(I-1,J)**2 + tauy_at_q(I-1,J)**2)) ) / & + ((G%mask2dBu(I,J) + G%mask2dBu(I-1,J-1)) + (G%mask2dBu(I,J-1) + G%mask2dBu(I-1,J))) ) + if (CS%read_gust_2d) gustiness = CS%gust(i,j) + endif + fluxes%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0*tau_mag) + enddo ; enddo + + elseif (wind_stagger == AGRID) then + call pass_vector(taux_at_h, tauy_at_h, G%Domain,stagger=AGRID) + + do j=js,je ; do I=Isq,Ieq + fluxes%taux(I,j) = 0.0 + if ((G%mask2dT(i,j) + G%mask2dT(i+1,j)) > 0) & + fluxes%taux(I,j) = (G%mask2dT(i,j)*taux_at_h(i,j) + & + G%mask2dT(i+1,j)*taux_at_h(i+1,j)) / & + (G%mask2dT(i,j) + G%mask2dT(i+1,j)) + enddo ; enddo + + do J=Jsq,Jeq ; do i=is,ie + fluxes%tauy(i,J) = 0.0 + if ((G%mask2dT(i,j) + G%mask2dT(i,j+1)) > 0) & + fluxes%tauy(i,J) = (G%mask2dT(i,j)*tauy_at_h(i,j) + & + G%mask2dT(i,J+1)*tauy_at_h(i,j+1)) / & + (G%mask2dT(i,j) + G%mask2dT(i,j+1)) + enddo ; enddo + + do j=js,je ; do i=is,ie + gustiness = CS%gust_const + if (CS%read_gust_2d .and. (G%mask2dT(i,j) > 0)) gustiness = CS%gust(i,j) + fluxes%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * G%mask2dT(i,j) * & + sqrt(taux_at_h(i,j)**2 + tauy_at_h(i,j)**2)) + enddo ; enddo + + else ! C-grid wind stresses. + if (G%symmetric) & + call fill_symmetric_edges(fluxes%taux, fluxes%tauy, G%Domain) + call pass_vector(fluxes%taux, fluxes%tauy, G%Domain) + + do j=js,je ; do i=is,ie + taux2 = 0.0 + if ((G%mask2dCu(I-1,j) + G%mask2dCu(I,j)) > 0) & + taux2 = (G%mask2dCu(I-1,j)*fluxes%taux(I-1,j)**2 + & + G%mask2dCu(I,j)*fluxes%taux(I,j)**2) / (G%mask2dCu(I-1,j) + G%mask2dCu(I,j)) + + tauy2 = 0.0 + if ((G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) > 0) & + tauy2 = (G%mask2dCv(i,J-1)*fluxes%tauy(i,J-1)**2 + & + G%mask2dCv(i,J)*fluxes%tauy(i,J)**2) / (G%mask2dCv(i,J-1) + G%mask2dCv(i,J)) + + if (CS%read_gust_2d) then + fluxes%ustar(i,j) = sqrt(CS%gust(i,j)*Irho0 + Irho0*sqrt(taux2 + tauy2)) + else + fluxes%ustar(i,j) = sqrt(CS%gust_const*Irho0 + Irho0*sqrt(taux2 + tauy2)) + endif + enddo ; enddo + + endif ! endif for wind related fields + + + ! sea ice related fields + if (CS%rigid_sea_ice) then + ! The commented out code here and in the following lines is the correct + ! version, but the incorrect version is being retained temporarily to avoid + ! changing answers. + call pass_var(fluxes%p_surf_full, G%Domain) + I_GEarth = 1.0 / G%G_Earth + Kv_rho_ice = (CS%kv_sea_ice / CS%density_sea_ice) + do I=isd,ied-1 ; do j=jsd,jed + mass_ice = min(fluxes%p_surf_full(i,j), fluxes%p_surf_full(i+1,j)) * I_GEarth + mass_eff = 0.0 + if (mass_ice > CS%rigid_sea_ice_mass) then + mass_eff = (mass_ice - CS%rigid_sea_ice_mass) **2 / & + (mass_ice + CS%rigid_sea_ice_mass) + endif + ! CAUTION: with both rigid_sea_ice and ice shelves, we will need to make this + ! a maximum for the second call. + fluxes%rigidity_ice_u(I,j) = Kv_rho_ice * mass_eff + enddo ; enddo + do i=isd,ied ; do J=jsd,jed-1 + mass_ice = min(fluxes%p_surf_full(i,j), fluxes%p_surf_full(i,j+1)) * I_GEarth + mass_eff = 0.0 + if (mass_ice > CS%rigid_sea_ice_mass) then + mass_eff = (mass_ice - CS%rigid_sea_ice_mass) **2 / & + (mass_ice + CS%rigid_sea_ice_mass) + endif + fluxes%rigidity_ice_v(i,J) = Kv_rho_ice * mass_eff + enddo ; enddo + endif + + if (coupler_type_initialized(fluxes%tr_fluxes) .and. & + coupler_type_initialized(IOB%fluxes)) & + call coupler_type_copy_data(IOB%fluxes, fluxes%tr_fluxes) + + if (CS%allow_flux_adjustments) then + ! Apply adjustments to fluxes + call apply_flux_adjustments(G, CS, Time, fluxes) + endif + + ! Allow for user-written code to alter fluxes after all the above + call user_alter_forcing(state, fluxes, Time, G, CS%urf_CS) + + call cpu_clock_end(id_clock_forcing) +end subroutine convert_IOB_to_fluxes + +!> Adds flux adjustments obtained via data_override +!! Component name is 'OCN' +!! Available adjustments are: +!! - taux_adj (Zonal wind stress delta, positive to the east, in Pa) +!! - tauy_adj (Meridional wind stress delta, positive to the north, in Pa) +subroutine apply_flux_adjustments(G, CS, Time, fluxes) + type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure + type(surface_forcing_CS), pointer :: CS !< Surface forcing control structure + type(time_type), intent(in) :: Time !< Model time structure + type(forcing), optional, intent(inout) :: fluxes !< Surface fluxes structure + ! Local variables + real, dimension(SZI_(G),SZJ_(G)) :: tempx_at_h ! Delta to zonal wind stress at h points (Pa) + real, dimension(SZI_(G),SZJ_(G)) :: tempy_at_h ! Delta to meridional wind stress at h points (Pa) + real, dimension(SZI_(G),SZJ_(G)) :: temp_at_h ! Fluxes at h points (W m-2 or kg m-2 s-1) + + integer :: isc, iec, jsc, jec, i, j + real :: dLonDx, dLonDy, rDlon, cosA, sinA, zonal_tau, merid_tau + logical :: overrode_x, overrode_y, overrode_h + + isc = G%isc; iec = G%iec + jsc = G%jsc; jec = G%jec + + overrode_h = .false. + call data_override('OCN', 'hflx_adj', temp_at_h(isc:iec,jsc:jec), Time, override=overrode_h) + + if (overrode_h) then + do j=G%jsc,G%jec ; do i=G%isc,G%iec + fluxes%heat_added(i,j) = fluxes%heat_added(i,j) + temp_at_h(i,j)* G%mask2dT(i,j) + enddo; enddo + endif + + call pass_var(fluxes%heat_added, G%Domain) + + overrode_h = .false. + call data_override('OCN', 'sflx_adj', temp_at_h(isc:iec,jsc:jec), Time, override=overrode_h) + + if (overrode_h) then + do j=G%jsc,G%jec ; do i=G%isc,G%iec + fluxes%salt_flux_added(i,j) = fluxes%salt_flux_added(i,j) + temp_at_h(i,j)* G%mask2dT(i,j) + enddo; enddo + endif + + call pass_var(fluxes%salt_flux_added, G%Domain) + overrode_h = .false. + + call data_override('OCN', 'prcme_adj', temp_at_h(isc:iec,jsc:jec), Time, override=overrode_h) + + if (overrode_h) then + do j=G%jsc,G%jec ; do i=G%isc,G%iec + fluxes%vprec(i,j) = fluxes%vprec(i,j) + temp_at_h(i,j)* G%mask2dT(i,j) + enddo; enddo + endif + + call pass_var(fluxes%vprec, G%Domain) + + + tempx_at_h(:,:) = 0.0 ; tempy_at_h(:,:) = 0.0 + ! Either reads data or leaves contents unchanged + overrode_x = .false. ; overrode_y = .false. + call data_override('OCN', 'taux_adj', tempx_at_h(isc:iec,jsc:jec), Time, override=overrode_x) + call data_override('OCN', 'tauy_adj', tempy_at_h(isc:iec,jsc:jec), Time, override=overrode_y) + + if (overrode_x .or. overrode_y) then + if (.not. (overrode_x .and. overrode_y)) call MOM_error(FATAL,"apply_flux_adjustments: "//& + "Both taux_adj and tauy_adj must be specified, or neither, in data_table") + + ! Rotate winds + call pass_vector(tempx_at_h, tempy_at_h, G%Domain, To_All, AGRID) + do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1 + dLonDx = G%geoLonCu(I,j) - G%geoLonCu(I-1,j) + dLonDy = G%geoLonCv(i,J) - G%geoLonCv(i,J-1) + rDlon = sqrt( dLonDx * dLonDx + dLonDy * dLonDy ) + if (rDlon > 0.) rDlon = 1. / rDlon + cosA = dLonDx * rDlon + sinA = dLonDy * rDlon + zonal_tau = tempx_at_h(i,j) + merid_tau = tempy_at_h(i,j) + tempx_at_h(i,j) = cosA * zonal_tau - sinA * merid_tau + tempy_at_h(i,j) = sinA * zonal_tau + cosA * merid_tau + enddo ; enddo + + ! Average to C-grid locations + do j=G%jsc,G%jec ; do I=G%isc-1,G%iec + fluxes%taux(I,j) = fluxes%taux(I,j) + 0.5 * ( tempx_at_h(i,j) + tempx_at_h(i+1,j) ) + enddo ; enddo + + do J=G%jsc-1,G%jec ; do i=G%isc,G%iec + fluxes%tauy(i,J) = fluxes%tauy(i,J) + 0.5 * ( tempy_at_h(i,j) + tempy_at_h(i,j+1) ) + enddo ; enddo + endif ! overrode_x .or. overrode_y + +end subroutine apply_flux_adjustments + !> Finalizes MOM6 !! !! \todo This needs to be done here. From 63c77868f143846f1702eb6770dc6a464ee65bb8 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 16 Oct 2017 13:19:39 -0600 Subject: [PATCH 19/20] Removes dependencies from modules in coupled_driver This commit removes all the dependencies in ocn_comp_mct.F90 related to modules located inside coupled_driver. The subroutines that were used by ocn_comp_mct.F90 are now located inside this module. I also deleted the ice_ocean_boundary type since it is no longer needed. Fluxes are now imported directly to the type fluxes. The module ocn_comp_mct is now fully Doxygenized. --- config_src/mct_driver/ocn_comp_mct.F90 | 1329 +++++++++++------------- 1 file changed, 578 insertions(+), 751 deletions(-) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index 10b28cc462..99e2feeabe 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -29,36 +29,36 @@ module ocn_comp_mct shr_file_setLogUnit, shr_file_setLogLevel ! MOM6 modules -use MOM_coms, only : reproducing_sum -use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end -use MOM_cpu_clock, only : CLOCK_SUBCOMPONENT -use MOM, only: initialize_MOM, step_MOM, MOM_control_struct, MOM_end -use MOM, only: calculate_surface_state, allocate_surface_state -use MOM, only: finish_MOM_initialization, step_offline -use MOM_forcing_type, only: forcing, forcing_diags, register_forcing_type_diags -use MOM_forcing_type, only: allocate_forcing_type, deallocate_forcing_type -use MOM_forcing_type, only: mech_forcing_diags, forcing_accumulate, forcing_diagnostics -use MOM_restart, only: save_restart -use MOM_domains, only: MOM_infra_init, num_pes, root_pe, pe_here -use MOM_domains, only: pass_vector, BGRID_NE, CGRID_NE, To_All -use MOM_domains, only: pass_var, AGRID, fill_symmetric_edges -use MOM_grid, only: ocean_grid_type, get_global_grid_size -use MOM_verticalGrid, only: verticalGrid_type -use MOM_variables, only: surface -use MOM_error_handler, only: MOM_error, FATAL, is_root_pe, WARNING -use MOM_error_handler, only: callTree_enter, callTree_leave -use MOM_time_manager, only: time_type, set_date, set_time, set_calendar_type, NOLEAP, get_date -use MOM_time_manager, only: operator(+), operator(-), operator(*), operator(/) -use MOM_time_manager, only: operator(/=), operator(>), get_time -use MOM_file_parser, only: get_param, log_version, param_file_type, close_param_file -use MOM_get_input, only: Get_MOM_Input, directories -use MOM_diag_mediator, only: diag_ctrl, enable_averaging, disable_averaging -use MOM_diag_mediator, only: diag_mediator_close_registration, diag_mediator_end -use MOM_diag_mediator, only: safe_alloc_ptr -use MOM_ice_shelf, only: initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS -use MOM_ice_shelf, only: ice_shelf_end, ice_shelf_save_restart -use MOM_sum_output, only: MOM_sum_output_init, sum_output_CS -use MOM_sum_output, only: write_energy, accumulate_net_input +use MOM_coms, only : reproducing_sum +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end +use MOM_cpu_clock, only : CLOCK_SUBCOMPONENT +use MOM, only: initialize_MOM, step_MOM, MOM_control_struct, MOM_end +use MOM, only: calculate_surface_state, allocate_surface_state +use MOM, only: finish_MOM_initialization, step_offline +use MOM_forcing_type, only: forcing, forcing_diags, register_forcing_type_diags +use MOM_forcing_type, only: allocate_forcing_type, deallocate_forcing_type +use MOM_forcing_type, only: mech_forcing_diags, forcing_accumulate, forcing_diagnostics +use MOM_restart, only: save_restart +use MOM_domains, only: MOM_infra_init, num_pes, root_pe, pe_here +use MOM_domains, only: pass_vector, BGRID_NE, CGRID_NE, To_All +use MOM_domains, only: pass_var, AGRID, fill_symmetric_edges +use MOM_grid, only: ocean_grid_type, get_global_grid_size +use MOM_verticalGrid, only: verticalGrid_type +use MOM_variables, only: surface +use MOM_error_handler, only: MOM_error, FATAL, is_root_pe, WARNING +use MOM_error_handler, only: callTree_enter, callTree_leave +use MOM_time_manager, only: time_type, set_date, set_time, set_calendar_type, NOLEAP, get_date +use MOM_time_manager, only: operator(+), operator(-), operator(*), operator(/) +use MOM_time_manager, only: operator(/=), operator(>), get_time +use MOM_file_parser, only: get_param, log_version, param_file_type, close_param_file +use MOM_get_input, only: Get_MOM_Input, directories +use MOM_diag_mediator, only: diag_ctrl, enable_averaging, disable_averaging +use MOM_diag_mediator, only: diag_mediator_close_registration, diag_mediator_end +use MOM_diag_mediator, only: safe_alloc_ptr +use MOM_ice_shelf, only: initialize_ice_shelf, shelf_calc_flux, ice_shelf_CS +use MOM_ice_shelf, only: ice_shelf_end, ice_shelf_save_restart +use MOM_sum_output, only: MOM_sum_output_init, sum_output_CS +use MOM_sum_output, only: write_energy, accumulate_net_input use MOM_string_functions, only: uppercase use MOM_constants, only: CELSIUS_KELVIN_OFFSET, hlf, hlv use MOM_EOS, only: gsw_sp_from_sr, gsw_pt_from_ct @@ -70,347 +70,297 @@ module ocn_comp_mct use MOM_io, only : slasher, write_version_number use MOM_spatial_means, only : adjust_area_mean_to_zero -! FMS -use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain -use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain +! FMS modules +use mpp_domains_mod, only : domain2d, mpp_get_layout, mpp_get_global_domain +use mpp_domains_mod, only : mpp_define_domains, mpp_get_compute_domain use time_interp_external_mod, only : init_external_field, time_interp_external use time_interp_external_mod, only : time_interp_external_init use fms_mod, only : read_data -! GFDL coupler +! GFDL coupler modules use coupler_types_mod, only : coupler_1d_bc_type, coupler_2d_bc_type -use coupler_types_mod, only : coupler_type_spawn -use coupler_types_mod, only : coupler_type_initialized, coupler_type_copy_data +use coupler_types_mod, only : coupler_type_spawn +use coupler_types_mod, only : coupler_type_initialized, coupler_type_copy_data ! By default make data private implicit none; private #include - ! Public member functions - public :: ocn_init_mct - public :: ocn_run_mct - public :: ocn_final_mct - ! Flag for debugging - logical, parameter :: debug=.true. - - !> Structure with MCT attribute vectors and indices - type cpl_indices - - ! ocean to coupler - integer :: o2x_So_t !< Surface potential temperature (deg C) - integer :: o2x_So_u !< Surface zonal velocity (m/s) - integer :: o2x_So_v !< Surface meridional velocity (m/s) - integer :: o2x_So_s !< Surface salinity (PSU) - integer :: o2x_So_dhdx !< Zonal slope in the sea surface height - integer :: o2x_So_dhdy !< Meridional lope in the sea surface height - integer :: o2x_So_bldepth !< Boundary layer depth (m) - integer :: o2x_Fioo_q !< Heat flux? - integer :: o2x_Faoo_fco2_ocn!< CO2 flux - integer :: o2x_Faoo_fdms_ocn!< DMS flux - - ! coupler to ocean - integer :: x2o_Si_ifrac !< Fractional ice wrt ocean - integer :: x2o_So_duu10n !< 10m wind speed squared (m^2/s^2) - integer :: x2o_Sa_pslv !< Sea-level pressure (Pa) - integer :: x2o_Sa_co2prog !< Bottom atm level prognostic CO2 - integer :: x2o_Sa_co2diag !< Bottom atm level diagnostic CO2 - integer :: x2o_Sw_lamult !< Wave model langmuir multiplier - integer :: x2o_Sw_ustokes !< Surface Stokes drift, x-component - integer :: x2o_Sw_vstokes !< Surface Stokes drift, y-component - integer :: x2o_Foxx_taux !< Zonal wind stress (W/m2) - integer :: x2o_Foxx_tauy !< Meridonal wind stress (W/m2) - integer :: x2o_Foxx_swnet !< Net short-wave heat flux (W/m2) - integer :: x2o_Foxx_sen !< Sensible heat flux (W/m2) - integer :: x2o_Foxx_lat !< Latent heat flux (W/m2) - integer :: x2o_Foxx_lwup !< Longwave radiation, up (W/m2) - integer :: x2o_Faxa_lwdn !< Longwave radiation, down (W/m2) - integer :: x2o_Fioi_melth !< Heat flux from snow & ice melt (W/m2) - integer :: x2o_Fioi_meltw !< Snow melt flux (kg/m2/s) - integer :: x2o_Fioi_bcpho !< Black Carbon hydrophobic release - !! from sea ice component - integer :: x2o_Fioi_bcphi !< Black Carbon hydrophilic release from - !! sea ice component - integer :: x2o_Fioi_flxdst !< Dust release from sea ice component - integer :: x2o_Fioi_salt !< Salt flux (kg(salt)/m2/s) - integer :: x2o_Foxx_evap !< Evaporation flux (kg/m2/s) - integer :: x2o_Faxa_prec !< Total precipitation flux (kg/m2/s) - integer :: x2o_Faxa_snow !< Water flux due to snow (kg/m2/s) - integer :: x2o_Faxa_rain !< Water flux due to rain (kg/m2/s) - integer :: x2o_Faxa_bcphidry !< Black Carbon hydrophilic dry deposition - integer :: x2o_Faxa_bcphodry !< Black Carbon hydrophobic dry deposition - integer :: x2o_Faxa_bcphiwet !< Black Carbon hydrophilic wet deposition - integer :: x2o_Faxa_ocphidry !< Organic Carbon hydrophilic dry deposition - integer :: x2o_Faxa_ocphodry !< Organic Carbon hydrophobic dry deposition - integer :: x2o_Faxa_ocphiwet !< Organic Carbon hydrophilic dry deposition - integer :: x2o_Faxa_dstwet1 !< Size 1 dust -- wet deposition - integer :: x2o_Faxa_dstwet2 !< Size 2 dust -- wet deposition - integer :: x2o_Faxa_dstwet3 !< Size 3 dust -- wet deposition - integer :: x2o_Faxa_dstwet4 !< Size 4 dust -- wet deposition - integer :: x2o_Faxa_dstdry1 !< Size 1 dust -- dry deposition - integer :: x2o_Faxa_dstdry2 !< Size 2 dust -- dry deposition - integer :: x2o_Faxa_dstdry3 !< Size 3 dust -- dry deposition - integer :: x2o_Faxa_dstdry4 !< Size 4 dust -- dry deposition - integer :: x2o_Foxx_rofl !< River runoff flux (kg/m2/s) - integer :: x2o_Foxx_rofi !< Ice runoff flux (kg/m2/s) - - ! optional per thickness category fields - - integer, dimension(:), allocatable :: x2o_frac_col !< Fraction of ocean cell, - !! per column - integer, dimension(:), allocatable :: x2o_fracr_col!< Fraction of ocean cell used - !! in radiation computations, - !! per column - integer, dimension(:), allocatable :: x2o_qsw_fracr_col !< qsw * fracr, per column +! Public member functions +public :: ocn_init_mct +public :: ocn_run_mct +public :: ocn_final_mct +! Flag for debugging +logical, parameter :: debug=.true. +!> Structure with MCT attribute vectors and indices +type cpl_indices + + ! ocean to coupler + integer :: o2x_So_t !< Surface potential temperature (deg C) + integer :: o2x_So_u !< Surface zonal velocity (m/s) + integer :: o2x_So_v !< Surface meridional velocity (m/s) + integer :: o2x_So_s !< Surface salinity (PSU) + integer :: o2x_So_dhdx !< Zonal slope in the sea surface height + integer :: o2x_So_dhdy !< Meridional lope in the sea surface height + integer :: o2x_So_bldepth !< Boundary layer depth (m) + integer :: o2x_Fioo_q !< Heat flux? + integer :: o2x_Faoo_fco2_ocn!< CO2 flux + integer :: o2x_Faoo_fdms_ocn!< DMS flux + + ! coupler to ocean + integer :: x2o_Si_ifrac !< Fractional ice wrt ocean + integer :: x2o_So_duu10n !< 10m wind speed squared (m^2/s^2) + integer :: x2o_Sa_pslv !< Sea-level pressure (Pa) + integer :: x2o_Sa_co2prog !< Bottom atm level prognostic CO2 + integer :: x2o_Sa_co2diag !< Bottom atm level diagnostic CO2 + integer :: x2o_Sw_lamult !< Wave model langmuir multiplier + integer :: x2o_Sw_ustokes !< Surface Stokes drift, x-component + integer :: x2o_Sw_vstokes !< Surface Stokes drift, y-component + integer :: x2o_Foxx_taux !< Zonal wind stress (W/m2) + integer :: x2o_Foxx_tauy !< Meridonal wind stress (W/m2) + integer :: x2o_Foxx_swnet !< Net short-wave heat flux (W/m2) + integer :: x2o_Foxx_sen !< Sensible heat flux (W/m2) + integer :: x2o_Foxx_lat !< Latent heat flux (W/m2) + integer :: x2o_Foxx_lwup !< Longwave radiation, up (W/m2) + integer :: x2o_Faxa_lwdn !< Longwave radiation, down (W/m2) + integer :: x2o_Fioi_melth !< Heat flux from snow & ice melt (W/m2) + integer :: x2o_Fioi_meltw !< Snow melt flux (kg/m2/s) + integer :: x2o_Fioi_bcpho !< Black Carbon hydrophobic release + !! from sea ice component + integer :: x2o_Fioi_bcphi !< Black Carbon hydrophilic release from + !! sea ice component + integer :: x2o_Fioi_flxdst !< Dust release from sea ice component + integer :: x2o_Fioi_salt !< Salt flux (kg(salt)/m2/s) + integer :: x2o_Foxx_evap !< Evaporation flux (kg/m2/s) + integer :: x2o_Faxa_prec !< Total precipitation flux (kg/m2/s) + integer :: x2o_Faxa_snow !< Water flux due to snow (kg/m2/s) + integer :: x2o_Faxa_rain !< Water flux due to rain (kg/m2/s) + integer :: x2o_Faxa_bcphidry !< Black Carbon hydrophilic dry deposition + integer :: x2o_Faxa_bcphodry !< Black Carbon hydrophobic dry deposition + integer :: x2o_Faxa_bcphiwet !< Black Carbon hydrophilic wet deposition + integer :: x2o_Faxa_ocphidry !< Organic Carbon hydrophilic dry deposition + integer :: x2o_Faxa_ocphodry !< Organic Carbon hydrophobic dry deposition + integer :: x2o_Faxa_ocphiwet !< Organic Carbon hydrophilic dry deposition + integer :: x2o_Faxa_dstwet1 !< Size 1 dust -- wet deposition + integer :: x2o_Faxa_dstwet2 !< Size 2 dust -- wet deposition + integer :: x2o_Faxa_dstwet3 !< Size 3 dust -- wet deposition + integer :: x2o_Faxa_dstwet4 !< Size 4 dust -- wet deposition + integer :: x2o_Faxa_dstdry1 !< Size 1 dust -- dry deposition + integer :: x2o_Faxa_dstdry2 !< Size 2 dust -- dry deposition + integer :: x2o_Faxa_dstdry3 !< Size 3 dust -- dry deposition + integer :: x2o_Faxa_dstdry4 !< Size 4 dust -- dry deposition + integer :: x2o_Foxx_rofl !< River runoff flux (kg/m2/s) + integer :: x2o_Foxx_rofi !< Ice runoff flux (kg/m2/s) + + ! optional per thickness category fields + integer, dimension(:), allocatable :: x2o_frac_col !< Fraction of ocean cell, + !! per column + integer, dimension(:), allocatable :: x2o_fracr_col!< Fraction of ocean cell used + !! in radiation computations, + !! per column + integer, dimension(:), allocatable :: x2o_qsw_fracr_col !< qsw * fracr, per column end type cpl_indices - !> This type is used for communication with other components via the FMS coupler. - ! The element names and types can be changed only with great deliberation, hence - ! the persistnce of things like the cutsy element name "avg_kount". - type, public :: ocean_public_type - type(domain2d) :: Domain ! The domain for the surface fields. - logical :: is_ocean_pe ! .true. on processors that run the ocean model. - character(len=32) :: instance_name = '' ! A name that can be used to identify - ! this instance of an ocean model, for example - ! in ensembles when writing messages. - integer, pointer, dimension(:) :: pelist => NULL() ! The list of ocean PEs. - logical, pointer, dimension(:,:) :: maskmap =>NULL() ! A pointer to an array - ! indicating which logical processors are actually used for - ! the ocean code. The other logical processors would be all - ! land points and are not assigned to actual processors. - ! This need not be assigned if all logical processors are used. - - integer :: stagger = -999 ! The staggering relative to the tracer points - ! points of the two velocity components. Valid entries - ! include AGRID, BGRID_NE, CGRID_NE, BGRID_SW, and CGRID_SW, - ! corresponding to the community-standard Arakawa notation. - ! (These are named integers taken from mpp_parameter_mod.) - ! Following MOM, this is BGRID_NE by default when the ocean - ! is initialized, but here it is set to -999 so that a - ! global max across ocean and non-ocean processors can be - ! used to determine its value. - real, pointer, dimension(:,:) :: & - t_surf => NULL(), & ! SST on t-cell (degrees Kelvin) - s_surf => NULL(), & ! SSS on t-cell (psu) - u_surf => NULL(), & ! i-velocity at the locations indicated by stagger, m/s. - v_surf => NULL(), & ! j-velocity at the locations indicated by stagger, m/s. - sea_lev => NULL(), & ! Sea level in m after correction for surface pressure, - ! i.e. dzt(1) + eta_t + patm/rho0/grav (m) - frazil =>NULL(), & ! Accumulated heating (in Joules/m^2) from frazil - ! formation in the ocean. - area => NULL() ! cell area of the ocean surface, in m2. - type(coupler_2d_bc_type) :: fields ! A structure that may contain an - ! array of named tracer-related fields. - integer :: avg_kount ! Used for accumulating averages of this type. - integer, dimension(2) :: axes = 0 ! Axis numbers that are available +!> This type is used for communication with other components via the FMS coupler. +! The element names and types can be changed only with great deliberation, hence +! the persistnce of things like the cutsy element name "avg_kount". +type, public :: ocean_public_type + type(domain2d) :: Domain !< The domain for the surface fields. + logical :: is_ocean_pe !! .true. on processors that run the ocean model. + character(len=32) :: instance_name = '' !< A name that can be used to identify + !! this instance of an ocean model, for example + !! in ensembles when writing messages. + integer, pointer, dimension(:) :: pelist => NULL() !< The list of ocean PEs. + logical, pointer, dimension(:,:) :: maskmap =>NULL() !< A pointer to an array + !! indicating which logical processors are actually + !! used for the ocean code. The other logical + !! processors would be all land points and are not + !! assigned to actual processors. This need not be + !! assigned if all logical processors are used. + + integer :: stagger = -999 !< The staggering relative to the tracer points + !! of the two velocity components. Valid entries + !! include AGRID, BGRID_NE, CGRID_NE, BGRID_SW, and CGRID_SW, + !! corresponding to the community-standard Arakawa notation. + !! (These are named integers taken from mpp_parameter_mod.) + !! Following MOM, this is BGRID_NE by default when the ocean + !! is initialized, but here it is set to -999 so that a + !! global max across ocean and non-ocean processors can be + !! used to determine its value. + real, pointer, dimension(:,:) :: & + t_surf => NULL(), & !< SST on t-cell (degrees Kelvin) + s_surf => NULL(), & !< SSS on t-cell (psu) + u_surf => NULL(), & !< i-velocity at the locations indicated by stagger, m/s. + v_surf => NULL(), & !< j-velocity at the locations indicated by stagger, m/s. + sea_lev => NULL(), & !< Sea level in m after correction for surface pressure, + !! i.e. dzt(1) + eta_t + patm/rho0/grav (m) + frazil =>NULL(), & !< Accumulated heating (in Joules/m^2) from frazil + !! formation in the ocean. + area => NULL() !< cell area of the ocean surface, in m2. + type(coupler_2d_bc_type) :: fields !< A structure that may contain an + !! array of named tracer-related fields. + integer :: avg_kount !< Used for accumulating averages of this type. + integer, dimension(2) :: axes = 0 !< Axis numbers that are available ! for I/O using this surface data. - end type ocean_public_type - - !> This type contains pointers to the forcing fields which may be used to drive MOM. - !! All fluxes are positive downward. - type, public :: surface_forcing_CS ; private - integer :: wind_stagger ! AGRID, BGRID_NE, or CGRID_NE (integer values - ! from MOM_domains) to indicate the staggering of - ! the winds that are being provided in calls to - ! update_ocean_model. - logical :: use_temperature ! If true, temp and saln used as state variables - real :: wind_stress_multiplier !< A multiplier applied to incoming wind stress (nondim). - - ! smg: remove when have A=B code reconciled - logical :: bulkmixedlayer ! If true, model based on bulk mixed layer code - - real :: Rho0 ! Boussinesq reference density (kg/m^3) - real :: area_surf = -1.0 ! total ocean surface area (m^2) - real :: latent_heat_fusion ! latent heat of fusion (J/kg) - real :: latent_heat_vapor ! latent heat of vaporization (J/kg) - - real :: max_p_surf ! maximum surface pressure that can be - ! exerted by the atmosphere and floating sea-ice, - ! in Pa. This is needed because the FMS coupling - ! structure does not limit the water that can be - ! frozen out of the ocean and the ice-ocean heat - ! fluxes are treated explicitly. - logical :: use_limited_P_SSH ! If true, return the sea surface height with - ! the correction for the atmospheric (and sea-ice) - ! pressure limited by max_p_surf instead of the - ! full atmospheric pressure. The default is true. - - real :: gust_const ! constant unresolved background gustiness for ustar (Pa) - logical :: read_gust_2d ! If true, use a 2-dimensional gustiness supplied - ! from an input file. - real, pointer, dimension(:,:) :: & - TKE_tidal => NULL(), & ! turbulent kinetic energy introduced to the - ! bottom boundary layer by drag on the tidal flows, - ! in W m-2. - gust => NULL(), & ! spatially varying unresolved background - ! gustiness that contributes to ustar (Pa). - ! gust is used when read_gust_2d is true. - ustar_tidal => NULL() ! tidal contribution to the bottom friction velocity (m/s) - real :: cd_tides ! drag coefficient that applies to the tides (nondimensional) - real :: utide ! constant tidal velocity to use if read_tideamp - ! is false, in m s-1. - logical :: read_tideamp ! If true, spatially varying tidal amplitude read from a file. - - logical :: rigid_sea_ice ! If true, sea-ice exerts a rigidity that acts - ! to damp surface deflections (especially surface - ! gravity waves). The default is false. - real :: Kv_sea_ice ! viscosity in sea-ice that resists sheared vertical motions (m^2/s) - real :: density_sea_ice ! typical density of sea-ice (kg/m^3). The value is - ! only used to convert the ice pressure into - ! appropriate units for use with Kv_sea_ice. - real :: rigid_sea_ice_mass ! A mass per unit area of sea-ice beyond which - ! sea-ice viscosity becomes effective, in kg m-2, - ! typically of order 1000 kg m-2. - logical :: allow_flux_adjustments ! If true, use data_override to obtain flux adjustments - - real :: Flux_const ! piston velocity for surface restoring (m/s) - logical :: salt_restore_as_sflux ! If true, SSS restore as salt flux instead of water flux - logical :: adjust_net_srestore_to_zero ! adjust srestore to zero (for both salt_flux or vprec) - logical :: adjust_net_srestore_by_scaling ! adjust srestore w/o moving zero contour - logical :: adjust_net_fresh_water_to_zero ! adjust net surface fresh-water (w/ restoring) to zero - logical :: adjust_net_fresh_water_by_scaling ! adjust net surface fresh-water w/o moving zero contour - logical :: mask_srestore_under_ice ! If true, use an ice mask defined by frazil - ! criteria for salinity restoring. - real :: ice_salt_concentration ! salt concentration for sea ice (kg/kg) - logical :: mask_srestore_marginal_seas ! if true, then mask SSS restoring in marginal seas - real :: max_delta_srestore ! maximum delta salinity used for restoring - real :: max_delta_trestore ! maximum delta sst used for restoring - real, pointer, dimension(:,:) :: basin_mask => NULL() ! mask for SSS restoring - - type(diag_ctrl), pointer :: diag ! structure to regulate diagnostic output timing - character(len=200) :: inputdir ! directory where NetCDF input files are - character(len=200) :: salt_restore_file ! filename for salt restoring data - character(len=30) :: salt_restore_var_name ! name of surface salinity in salt_restore_file - character(len=200) :: temp_restore_file ! filename for sst restoring data - character(len=30) :: temp_restore_var_name ! name of surface temperature in temp_restore_file - - integer :: id_srestore = -1 ! id number for time_interp_external. - integer :: id_trestore = -1 ! id number for time_interp_external. - - ! Diagnostics handles - type(forcing_diags), public :: handles - - !### type(ctrl_forcing_CS), pointer :: ctrl_forcing_CSp => NULL() - type(MOM_restart_CS), pointer :: restart_CSp => NULL() - type(user_revise_forcing_CS), pointer :: urf_CS => NULL() - end type surface_forcing_CS - - !> Contains information about the ocean state, although it is not necessary that - !! this is implemented with all models. This type is private, and can therefore vary - !! between different ocean models. - type, public :: ocean_state_type ; private - logical :: is_ocean_PE = .false. !< True if this is an ocean PE. - type(time_type) :: Time !< The ocean model's time and master clock. - integer :: Restart_control !< An integer that is bit-tested to determine whether - !! incremental restart files are saved and whether they - !! have a time stamped name. +1 (bit 0) for generic - !! files and +2 (bit 1) for time-stamped files. A - !! restart file is saved at the end of a run segment - !! unless Restart_control is negative. - type(time_type) :: energysavedays ! The interval between writing the energies - ! and other integral quantities of the run. - type(time_type) :: write_energy_time ! The next time to write to the energy file. - - integer :: nstep = 0 ! The number of calls to update_ocean. - logical :: use_ice_shelf ! If true, the ice shelf model is enabled. - logical :: icebergs_apply_rigid_boundary ! If true, the icebergs can change ocean bd condition. - real :: kv_iceberg ! The viscosity of the icebergs in m2/s (for ice rigidity) - real :: berg_area_threshold ! Fraction of grid cell which iceberg must occupy - !so that fluxes below are set to zero. (0.5 is a - !good value to use. Not applied for negative values. - real :: latent_heat_fusion ! Latent heat of fusion - real :: density_iceberg ! A typical density of icebergs in kg/m3 (for ice rigidity) - type(ice_shelf_CS), pointer :: Ice_shelf_CSp => NULL() - logical :: restore_salinity ! If true, the coupled MOM driver adds a term to - ! restore salinity to a specified value. - logical :: restore_temp ! If true, the coupled MOM driver adds a term to - ! restore sst to a specified value. - real :: press_to_z ! A conversion factor between pressure and ocean - ! depth in m, usually 1/(rho_0*g), in m Pa-1. - real :: C_p ! The heat capacity of seawater, in J K-1 kg-1. - - type(directories) :: dirs ! A structure containing several relevant directory paths. - type(forcing) :: fluxes ! A structure containing pointers to - ! the ocean forcing fields. - type(forcing) :: flux_tmp ! A secondary structure containing pointers to the - ! ocean forcing fields for when multiple coupled - ! timesteps are taken per thermodynamic step. - type(surface) :: state ! A structure containing pointers to - ! the ocean surface state fields. - type(ocean_grid_type), pointer :: grid => NULL() ! A pointer to a grid structure - ! containing metrics and related information. - type(verticalGrid_type), pointer :: GV => NULL() ! A pointer to a vertical grid - ! structure containing metrics and related information. - type(MOM_control_struct), pointer :: MOM_CSp => NULL() - type(surface_forcing_CS), pointer :: forcing_CSp => NULL() - type(sum_output_CS), pointer :: sum_output_CSp => NULL() - end type ocean_state_type - - !> Control structure corresponding to forcing, but with - !! the elements, units, and conventions that exactly conform to the use for - !! MOM-based coupled models. - type, public :: ice_ocean_boundary_type - real, pointer, dimension(:,:) :: u_flux =>NULL() ! i-direction wind stress (Pa) - real, pointer, dimension(:,:) :: v_flux =>NULL() ! j-direction wind stress (Pa) - real, pointer, dimension(:,:) :: t_flux =>NULL() ! sensible heat flux (W/m2) - real, pointer, dimension(:,:) :: q_flux =>NULL() ! specific humidity flux (kg/m2/s) - real, pointer, dimension(:,:) :: salt_flux =>NULL() ! salt flux (kg/m2/s) - real, pointer, dimension(:,:) :: lw_flux =>NULL() ! long wave radiation (W/m2) - real, pointer, dimension(:,:) :: sw_flux_vis_dir =>NULL() ! direct visible sw radiation (W/m2) - real, pointer, dimension(:,:) :: sw_flux_vis_dif =>NULL() ! diffuse visible sw radiation (W/m2) - real, pointer, dimension(:,:) :: sw_flux_nir_dir =>NULL() ! direct Near InfraRed sw radiation (W/m2) - real, pointer, dimension(:,:) :: sw_flux_nir_dif =>NULL() ! diffuse Near InfraRed sw radiation (W/m2) - real, pointer, dimension(:,:) :: latent =>NULL() ! latent heat flux (W/m2) - real, pointer, dimension(:,:) :: lprec =>NULL() ! mass flux of liquid precip (kg/m2/s) - real, pointer, dimension(:,:) :: fprec =>NULL() ! mass flux of frozen precip (kg/m2/s) - real, pointer, dimension(:,:) :: runoff =>NULL() ! mass flux of liquid runoff (kg/m2/s) - real, pointer, dimension(:,:) :: calving =>NULL() ! mass flux of frozen runoff (kg/m2/s) - real, pointer, dimension(:,:) :: ustar_berg =>NULL() ! frictional velocity beneath icebergs (m/s) - real, pointer, dimension(:,:) :: area_berg =>NULL() ! area covered by icebergs(m2/m2) - real, pointer, dimension(:,:) :: mass_berg =>NULL() ! mass of icebergs(kg/m2) - real, pointer, dimension(:,:) :: runoff_hflx =>NULL() ! heat content of liquid runoff (W/m2) - real, pointer, dimension(:,:) :: calving_hflx =>NULL() ! heat content of frozen runoff (W/m2) - real, pointer, dimension(:,:) :: p =>NULL() ! pressure of overlying ice and atmosphere - ! on ocean surface (Pa) - real, pointer, dimension(:,:) :: mi =>NULL() ! mass of ice (kg/m2) - integer :: xtype ! REGRID, REDIST or DIRECT - type(coupler_2d_bc_type) :: fluxes ! A structure that may contain an - ! array of named fields used for - ! passive tracer fluxes. - integer :: wind_stagger = -999 ! A flag indicating the spatial discretization of - ! wind stresses. This flag may be set by the - ! flux-exchange code, based on what the sea-ice - ! model is providing. Otherwise, the value from - ! the surface_forcing_CS is used. - end type ice_ocean_boundary_type - - !> Control structure for this module - type MCT_MOM_Data - - type(ocean_state_type), pointer :: ocn_state => NULL() !< The private state of ocean - type(ocean_public_type), pointer :: ocn_public => NULL() !< The public state of ocean - type(ocean_grid_type), pointer :: grid => NULL() !< The grid structure - type(surface), pointer :: ocn_surface => NULL() !< The ocean surface state - type(ice_ocean_boundary_type) :: ice_ocean_boundary !< The ice ocean boundary type - type(seq_infodata_type), pointer :: infodata !< The input info type - type(cpl_indices), public :: ind !< Variable IDs - ! runtime params - logical :: sw_decomp !< Controls whether shortwave is decomposed into four components - real :: c1, c2, c3, c4 !< Coeffs. used in the shortwave decomposition - ! i/o - character(len=384) :: pointer_filename !< Name of the ascii file that contains the path - !! and filename of the latest restart file. - integer :: stdout !< standard output unit. (by default, it should point to ocn.log.* file) - end type MCT_MOM_Data - - type(MCT_MOM_Data) :: glb !< global structure - - integer :: id_clock_forcing +end type ocean_public_type + +!> Contains pointers to the forcing fields which may be used to drive MOM. +!! All fluxes are positive downward. +type, public :: surface_forcing_CS ; private + integer :: wind_stagger !< AGRID, BGRID_NE, or CGRID_NE (integer values + !! from MOM_domains) to indicate the staggering of + !! the winds that are being provided in calls to + !! update_ocean_model. CIME uses AGRID, so this option + !! is being hard coded for now. + logical :: use_temperature !< If true, temp and saln used as state variables + real :: wind_stress_multiplier!< A multiplier applied to incoming wind stress (nondim). + ! smg: remove when have A=B code reconciled + logical :: bulkmixedlayer !< If true, model based on bulk mixed layer code + real :: Rho0 !< Boussinesq reference density (kg/m^3) + real :: area_surf = -1.0 !< total ocean surface area (m^2) + real :: latent_heat_fusion !< latent heat of fusion (J/kg) + real :: latent_heat_vapor !< latent heat of vaporization (J/kg) + real :: max_p_surf !< maximum surface pressure that can be + !! exerted by the atmosphere and floating sea-ice, + !! in Pa. This is needed because the FMS coupling + !! structure does not limit the water that can be + !! frozen out of the ocean and the ice-ocean heat + !! fluxes are treated explicitly. + logical :: use_limited_P_SSH !< If true, return the sea surface height with + !! the correction for the atmospheric (and sea-ice) + !! pressure limited by max_p_surf instead of the + !! full atmospheric pressure. The default is true. + real :: gust_const !< constant unresolved background gustiness for ustar (Pa) + logical :: read_gust_2d !< If true, use a 2-dimensional gustiness supplied + !! from an input file. + real, pointer, dimension(:,:) :: & + TKE_tidal => NULL(), & !< turbulent kinetic energy introduced to the + !! bottom boundary layer by drag on the tidal flows, + !! in W m-2. + gust => NULL(), & !< spatially varying unresolved background + !! gustiness that contributes to ustar (Pa). + !! gust is used when read_gust_2d is true. + ustar_tidal => NULL() !< tidal contribution to the bottom friction velocity (m/s) + real :: cd_tides !< drag coefficient that applies to the tides (nondimensional) + real :: utide !< constant tidal velocity to use if read_tideamp + !! is false, in m s-1. + logical :: read_tideamp !< If true, spatially varying tidal amplitude read from a file. + logical :: rigid_sea_ice !< If true, sea-ice exerts a rigidity that acts + !! to damp surface deflections (especially surface + !! gravity waves). The default is false. + real :: Kv_sea_ice !< viscosity in sea-ice that resists sheared vertical motions (m^2/s) + real :: density_sea_ice !< typical density of sea-ice (kg/m^3). The value is + !! only used to convert the ice pressure into + !! appropriate units for use with Kv_sea_ice. + real :: rigid_sea_ice_mass !< A mass per unit area of sea-ice beyond which + !! sea-ice viscosity becomes effective, in kg m-2, + !! typically of order 1000 kg m-2. + logical :: allow_flux_adjustments !< If true, use data_override to obtain flux adjustments + real :: Flux_const !< piston velocity for surface restoring (m/s) + logical :: salt_restore_as_sflux !< If true, SSS restore as salt flux instead of water flux + logical :: adjust_net_srestore_to_zero !< adjust srestore to zero (for both salt_flux or vprec) + logical :: adjust_net_srestore_by_scaling !< adjust srestore w/o moving zero contour + logical :: adjust_net_fresh_water_to_zero !< adjust net surface fresh-water (w/ restoring) to zero + logical :: adjust_net_fresh_water_by_scaling !< adjust net surface fresh-water w/o moving zero contour + logical :: mask_srestore_under_ice !< If true, use an ice mask defined by frazil + !! criteria for salinity restoring. + real :: ice_salt_concentration !< salt concentration for sea ice (kg/kg) + logical :: mask_srestore_marginal_seas !< if true, then mask SSS restoring in marginal seas + real :: max_delta_srestore !< maximum delta salinity used for restoring + real :: max_delta_trestore !< maximum delta sst used for restoring + real, pointer, dimension(:,:) :: basin_mask => NULL() !< mask for SSS restoring + type(diag_ctrl), pointer :: diag !< structure to regulate diagnostic output timing + character(len=200) :: inputdir !< directory where NetCDF input files are + character(len=200) :: salt_restore_file !< filename for salt restoring data + character(len=30) :: salt_restore_var_name !< name of surface salinity in salt_restore_file + character(len=200) :: temp_restore_file !< filename for sst restoring data + character(len=30) :: temp_restore_var_name !< name of surface temperature in temp_restore_file + integer :: id_srestore = -1 !< id number for time_interp_external. + integer :: id_trestore = -1 !< id number for time_interp_external. + type(forcing_diags), public :: handles !< diagnostics handles + !### type(ctrl_forcing_CS), pointer :: ctrl_forcing_CSp => NULL() + type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< restart pointer + type(user_revise_forcing_CS), pointer :: urf_CS => NULL()!< user revise pointer +end type surface_forcing_CS + +!> Contains information about the ocean state, although it is not necessary that +!! this is implemented with all models. This type is private, and can therefore vary +!! between different ocean models. +type, public :: ocean_state_type ; private + logical :: is_ocean_PE = .false. !< True if this is an ocean PE. + type(time_type) :: Time !< The ocean model's time and master clock. + integer :: Restart_control !< An integer that is bit-tested to determine whether + !! incremental restart files are saved and whether they + !! have a time stamped name. +1 (bit 0) for generic + !! files and +2 (bit 1) for time-stamped files. A + !! restart file is saved at the end of a run segment + !! unless Restart_control is negative. + type(time_type) :: energysavedays !< The interval between writing the energies + !! and other integral quantities of the run. + type(time_type) :: write_energy_time !< The next time to write to the energy file. + integer :: nstep = 0 !< The number of calls to update_ocean. + logical :: use_ice_shelf !< If true, the ice shelf model is enabled. + logical :: icebergs_apply_rigid_boundary !< If true, the icebergs can change ocean bd condition. + real :: kv_iceberg !< The viscosity of the icebergs in m2/s (for ice rigidity) + real :: berg_area_threshold !< Fraction of grid cell which iceberg must occupy + !! so that fluxes below are set to zero. (0.5 is a + !! good value to use. Not applied for negative values. + real :: latent_heat_fusion !< Latent heat of fusion + real :: density_iceberg !< A typical density of icebergs in kg/m3 (for ice rigidity) + type(ice_shelf_CS), pointer :: Ice_shelf_CSp => NULL() !< ice shelf structure. + logical :: restore_salinity !< If true, the coupled MOM driver adds a term to + !! restore salinity to a specified value. + logical :: restore_temp !< If true, the coupled MOM driver adds a term to + !! restore sst to a specified value. + real :: press_to_z !< A conversion factor between pressure and ocean + !! depth in m, usually 1/(rho_0*g), in m Pa-1. + real :: C_p !< The heat capacity of seawater, in J K-1 kg-1. + type(directories) :: dirs !< A structure containing several relevant directory paths. + type(forcing) :: fluxes !< A structure containing pointers to + !! the ocean forcing fields. + type(forcing) :: flux_tmp !< A secondary structure containing pointers to the + !! ocean forcing fields for when multiple coupled + !! timesteps are taken per thermodynamic step. + type(surface) :: state !< A structure containing pointers to + !! the ocean surface state fields. + type(ocean_grid_type), pointer :: grid => NULL() !< A pointer to a grid structure + !! containing metrics and related information. + type(verticalGrid_type), pointer :: GV => NULL() !< A pointer to a vertical grid + !! structure containing metrics and related information. + type(MOM_control_struct), pointer :: MOM_CSp => NULL() + type(surface_forcing_CS), pointer :: forcing_CSp => NULL() + type(sum_output_CS), pointer :: sum_output_CSp => NULL() +end type ocean_state_type + +!> Control structure for this module +type MCT_MOM_Data + + type(ocean_state_type), pointer :: ocn_state => NULL() !< The private state of ocean + type(ocean_public_type), pointer :: ocn_public => NULL() !< The public state of ocean + type(ocean_grid_type), pointer :: grid => NULL() !< The grid structure + type(surface), pointer :: ocn_surface => NULL() !< The ocean surface state + type(forcing) :: fluxes !< Structure that contains pointers to the + !! boundary forcing used to drive the liquid + !! ocean simulated by MOM. + type(seq_infodata_type), pointer :: infodata !< The input info type + type(cpl_indices), public :: ind !< Variable IDs + ! runtime params + logical :: sw_decomp !< Controls whether shortwave is decomposed into four components + real :: c1, c2, c3, c4 !< Coeffs. used in the shortwave decomposition + ! i/o + character(len=384) :: pointer_filename !< Name of the ascii file that contains the path + !! and filename of the latest restart file. + integer :: stdout !< standard output unit. (by default, it should point to ocn.log.* file) +end type MCT_MOM_Data + +type(MCT_MOM_Data) :: glb !< global structure +integer :: id_clock_forcing contains -!> Initializes MOM6 +!> This subroutine initializes MOM6. subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) type(ESMF_Clock), intent(inout) :: EClock !< Time and time step ? \todo Why must this !! be intent(inout)? @@ -461,13 +411,11 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) seconds_in_day = 86400.0d0, & minutes_in_hour = 60.0d0 - ! ocn model input namelist filename - character(len=99) :: ocn_modelio_name - integer :: shrlogunit ! original log file unit - integer :: shrloglev ! original log level + character(len=99) :: ocn_modelio_name !< ocn model input namelist filename + integer :: shrlogunit !< original log file unit + integer :: shrloglev !< original log level - ! instance control vars (these are local for now) - integer(kind=4) :: inst_index + integer(kind=4) :: inst_index !< instance control vars (these are local for now) character(len=16) :: inst_name character(len=16) :: inst_suffix @@ -704,32 +652,6 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) ! Size of global domain call get_global_grid_size(glb%grid, ni, nj) - ! Allocate ice_ocean_boundary using global indexing without halos - isc = glb%grid%isc + glb%grid%idg_offset - iec = glb%grid%iec + glb%grid%idg_offset - jsc = glb%grid%jsc + glb%grid%jdg_offset - jec = glb%grid%jec + glb%grid%jdg_offset - allocate(glb%ice_ocean_boundary%u_flux(isc:iec,jsc:jec)); glb%ice_ocean_boundary%u_flux(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%v_flux(isc:iec,jsc:jec)); glb%ice_ocean_boundary%v_flux(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%t_flux(isc:iec,jsc:jec)); glb%ice_ocean_boundary%t_flux(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%q_flux(isc:iec,jsc:jec)); glb%ice_ocean_boundary%q_flux(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%latent(isc:iec,jsc:jec)); glb%ice_ocean_boundary%latent(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%salt_flux(isc:iec,jsc:jec)); glb%ice_ocean_boundary%salt_flux(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%lw_flux(isc:iec,jsc:jec)); glb%ice_ocean_boundary%lw_flux(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%sw_flux_vis_dir(isc:iec,jsc:jec)); glb%ice_ocean_boundary%sw_flux_vis_dir(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%sw_flux_vis_dif(isc:iec,jsc:jec)); glb%ice_ocean_boundary%sw_flux_vis_dif(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%sw_flux_nir_dir(isc:iec,jsc:jec)); glb%ice_ocean_boundary%sw_flux_nir_dir(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%sw_flux_nir_dif(isc:iec,jsc:jec)); glb%ice_ocean_boundary%sw_flux_nir_dif(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%lprec(isc:iec,jsc:jec)); glb%ice_ocean_boundary%lprec(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%fprec(isc:iec,jsc:jec)); glb%ice_ocean_boundary%fprec(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%runoff(isc:iec,jsc:jec)); glb%ice_ocean_boundary%runoff(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%calving(isc:iec,jsc:jec)); glb%ice_ocean_boundary%calving(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%runoff_hflx(isc:iec,jsc:jec)); glb%ice_ocean_boundary%runoff_hflx(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%calving_hflx(isc:iec,jsc:jec)); glb%ice_ocean_boundary%calving_hflx(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%p(isc:iec,jsc:jec)); glb%ice_ocean_boundary%p(:,:) = 0.0 - allocate(glb%ice_ocean_boundary%mi(isc:iec,jsc:jec)); glb%ice_ocean_boundary%mi(:,:) = 0.0 - - if (debug .and. root_pe().eq.pe_here()) print *, "calling seq_infodata_putdata" call seq_infodata_PutData( glb%infodata, & @@ -737,7 +659,6 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) call seq_infodata_PutData( glb%infodata, & ocn_prognostic=.true., ocnrof_prognostic=.true.) - if (debug .and. root_pe().eq.pe_here()) print *, "leaving ocean_init_mct" ! Reset shr logging to original values @@ -877,23 +798,16 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i !! internal variables in the ice model. character(len=*), optional, intent(in) :: input_restart_file !< If present, name of restart file to read -! This subroutine initializes both the ocean state and the ocean surface type. +! This subroutine initializes both the ocean state and the ocean surface type. ! Because of the way that indicies and domains are handled, Ocean_sfc must have ! been used in a previous call to initialize_ocean_type. -! Arguments: Ocean_sfc - A structure containing various publicly visible ocean -! surface properties after initialization, this is intent(out). -! (out,private) OS - A structure whose internal contents are private -! to ocean_model_mod that may be used to contain all -! information about the ocean's interior state. -! (in) Time_init - The start time for the coupled model's calendar. -! (in) Time_in - The time at which to initialize the ocean model. - real :: Time_unit ! The time unit in seconds for ENERGYSAVEDAYS. - real :: Rho0 ! The Boussinesq ocean density, in kg m-3. - real :: G_Earth ! The gravitational acceleration in m s-2. -! This include declares and sets the variable "version". + real :: Time_unit !< The time unit in seconds for ENERGYSAVEDAYS. + real :: Rho0 !< The Boussinesq ocean density, in kg m-3. + real :: G_Earth !< The gravitational acceleration in m s-2. + !! This include declares and sets the variable "version". #include "version_variable.h" - character(len=40) :: mdl = "ocean_model_init" ! This module's name. + character(len=40) :: mdl = "ocean_model_init" !< This module's name. character(len=48) :: stagger integer :: secs, days type(param_file_type) :: param_file !< A structure to parse for run-time parameters @@ -1041,7 +955,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i call callTree_leave("ocean_model_init(") end subroutine ocean_model_init -!> This subroutine extracts the surface properties from the ocean's internal +!> Extracts the surface properties from the ocean's internal !! state and stores them in the ocean type returned to the calling ice model. !! It has to be separate from the ocean_initialization call because the coupler !! module allocates the space for some of these variables. @@ -1064,7 +978,7 @@ subroutine ocean_model_init_sfc(OS, Ocean_sfc) end subroutine ocean_model_init_sfc -!> Initialize surface forcing: get relevant parameters and allocate arrays +!> Initializes surface forcing: get relevant parameters and allocate arrays. subroutine surface_forcing_init(Time, G, param_file, diag, CS, restore_salt, restore_temp) type(time_type), intent(in) :: Time !< The current model time type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure @@ -1377,21 +1291,21 @@ subroutine surface_forcing_init(Time, G, param_file, diag, CS, restore_salt, res call cpu_clock_end(id_clock_forcing) end subroutine surface_forcing_init +!> Initializes domain and state variables contained in the ocean public type. subroutine initialize_ocean_public_type(input_domain, Ocean_sfc, diag, maskmap, & gas_fields_ocn) - type(domain2D), intent(in) :: input_domain - type(ocean_public_type), intent(inout) :: Ocean_sfc - type(diag_ctrl), intent(in) :: diag - logical, intent(in), optional :: maskmap(:,:) + type(domain2D), intent(in) :: input_domain !< The FMS domain for the input structure + type(ocean_public_type), intent(inout) :: Ocean_sfc !< Ocean surface state + type(diag_ctrl), intent(in) :: diag !< A structure used to control diagnostics. + logical, intent(in), optional :: maskmap(:,:) !< A pointer to an array indicating which + !! logical processors are actually used for the ocean code. type(coupler_1d_bc_type), & optional, intent(in) :: gas_fields_ocn !< If present, this type describes the - !! ocean and surface-ice fields that will participate - !! in the calculation of additional gas or other - !! tracer fluxes. - + !! ocean and surface-ice fields that will participate + !! in the calculation of additional gas or other + !! tracer fluxes. + ! local variables integer :: xsz, ysz, layout(2) - ! ice-ocean-boundary fields are always allocated using absolute indicies - ! and have no halos. integer :: isc, iec, jsc, jec call mpp_get_layout(input_domain,layout) @@ -1429,17 +1343,17 @@ end subroutine initialize_ocean_public_type !> Translates the coupler's ocean_data_type into MOM6's surface state variable. !! This may eventually be folded into the MOM6's code that calculates the -!! surface state in the first place. Note the offset in the arrays because -!! the ocean_data_type has no halo points in its arrays and always uses -!! absolute indicies. +!! surface state in the first place. subroutine convert_state_to_ocean_type(state, Ocean_sfc, G, use_conT_absS, & patm, press_to_z) type(surface), intent(inout) :: state - type(ocean_public_type), target, intent(inout) :: Ocean_sfc - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure - logical, intent(in) :: use_conT_absS - real, optional, intent(in) :: patm(:,:) - real, optional, intent(in) :: press_to_z + type(ocean_public_type), target, intent(inout) :: Ocean_sfc !< Ocean surface state + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + logical, intent(in) :: use_conT_absS !< If true, , the prognostics + !! T&S are the conservative temperature + real, optional, intent(in) :: patm(:,:) !< Atmospheric pressure. + real, optional, intent(in) :: press_to_z !< Factor to tranform atmospheric + !! pressure to z? ! local variables real :: IgR0 @@ -1526,7 +1440,9 @@ subroutine get_state_pointers(OS, grid, surf) end subroutine get_state_pointers -!> Maps outgoing ocean data to MCT buffer +!> Maps outgoing ocean data to MCT buffer. +!! See \ref section_ocn_export for a summary of the data +!! that is transferred from MOM6 to MCT. subroutine ocn_export(ind, ocn_public, grid, o2x) type(cpl_indices), intent(inout) :: ind !< Structure with coupler !! indices and vectors @@ -1618,123 +1534,6 @@ subroutine ocn_export(ind, ocn_public, grid, o2x) end subroutine ocn_export -!> Fills the ice ocean boundary type -subroutine fill_ice_ocean_bnd(ice_ocean_boundary, grid, x2o_o, ind, & - sw_decomp, c1, c2, c3, c4) - - type(ice_ocean_boundary_type), intent(inout) :: ice_ocean_boundary !< A type for the ice ocean boundary - type(ocean_grid_type), intent(in) :: grid !< - !type(mct_aVect), intent(in) :: x2o_o - real(kind=8), intent(in) :: x2o_o(:,:) - type(cpl_indices), intent(inout) :: ind - logical, intent(in) :: sw_decomp !< controls if shortwave is - !!decomposed into four components - real(kind=8), intent(in), optional :: c1, c2, c3, c4 !< Coeffs. used in the shortwave decomposition - - ! local variables - integer :: i, j, k, ig, jg !< grid indices - - ! /TODO The following comments summarizes the mismatches between MCT and MOM6 in terms - ! of ice ocean fluxes. - - ! Redundancies: - ! x2o_Foxx_lat, x2o_Faxa_prec, x2o_Faxa_evap - **latent from MOM6 and coupler are different** - ! x2o_Faxa_prec = x2o_Faxa_rain + x2o_Faxa_snow - - ! Variables that **could not** be verified so far: - ! x2o_Foxx_rofl - ! x2o_Foxx_rof - - ! Variables in IOB that are **NOT** filled by the coupler: - ! ustar_berg, frictional velocity beneath icebergs (m/s) - ! area_berg, area covered by icebergs(m2/m2) - ! mass_berg, mass of icebergs(kg/m2) - ! runoff_hflx, heat content of liquid runoff (W/m2) - ! calving_hflx, heat content of frozen runoff (W/m2) - ! mi, mass of ice (kg/m2) - - ! Variables in the coupler that are **NOT** used in MOM6 (i.e., no corresponding field in IOB): - ! x2o_Fioi_melth, heat flux from snow & ice melt (W/m2) - ! x2o_Fioi_meltw, snow melt flux (kg/m2/s) - ! x2o_Si_ifrac, fractional ice wrt ocean - ! x2o_So_duu10n, 10m wind speed squared (m^2/s^2) - ! x2o_Sa_co2prog, bottom atm level prognostic CO2 - ! x2o_Sa_co2diag, bottom atm level diagnostic CO2 - ! x2o_Foxx_lat, latent heat flux (W/m2) - - ! Langmuir related fields: - ! surface Stokes drift, x-comp. (x2o_Sw_ustokes) - ! surface Stokes drift, y-comp. (x2o_Sw_vstokes) - ! wave model langmuir multiplier (x2o_Sw_lamult) - - ! Biogeochemistry: - ! x2o_Fioi_bcpho, Black Carbon hydrophobic release from sea ice component - ! x2o_Fioi_bcphi, Black Carbon hydrophilic release from sea ice component - ! x2o_Fioi_flxdst, Dust release from sea ice component - ! x2o_Faxa_bcphidry, Black Carbon hydrophilic dry deposition - ! x2o_Faxa_bcphodry, Black Carbon hydrophobic dry deposition - ! x2o_Faxa_bcphiwet, Black Carbon hydrophobic wet deposition - ! x2o_Faxa_ocphidry, Organic Carbon hydrophilic dry deposition - ! x2o_Faxa_ocphodry, Organic Carbon hydrophobic dry deposition - ! x2o_Faxa_ocphiwet, Organic Carbon hydrophilic dry deposition - ! x2o_Faxa_dstwet, Sizes 1 to 4 dust - wet deposition - ! x2o_Faxa_dstdry, Sizes 1 to 4 dust - dry depositi - - k = 0 - do j = grid%jsc, grid%jec - jg = j + grid%jdg_offset - do i = grid%isc, grid%iec - k = k + 1 ! Increment position within gindex - ig = i + grid%idg_offset - ! sea-level pressure (Pa) - ice_ocean_boundary%p(ig,jg) = x2o_o(ind%x2o_Sa_pslv,k) - ! zonal wind stress (taux) - ice_ocean_boundary%u_flux(ig,jg) = x2o_o(ind%x2o_Foxx_taux,k) - ! meridional wind stress (tauy) - ice_ocean_boundary%v_flux(ig,jg) = x2o_o(ind%x2o_Foxx_tauy,k) - ! sensible heat flux (W/m2) - ice_ocean_boundary%t_flux(ig,jg) = -x2o_o(ind%x2o_Foxx_sen,k) - ! latent heat flux (W/m^2) - ice_ocean_boundary%latent(ig,jg) = x2o_o(ind%x2o_Foxx_lat,k) - ! salt flux - ice_ocean_boundary%salt_flux(ig,jg) = -x2o_o(ind%x2o_Fioi_salt,k) - ! heat content from frozen runoff - ice_ocean_boundary%calving_hflx(ig,jg) = 0.0 - ! snow melt flux - !ice_ocean_boundary%fprec(ig,jg) = x2o_o(ind%x2o_Fioi_meltw,k) - ! river runoff flux - ice_ocean_boundary%runoff(ig,jg) = x2o_o(ind%x2o_Foxx_rofl,k) - ! ice runoff flux - ice_ocean_boundary%calving(ig,jg) = x2o_o(ind%x2o_Foxx_rofi,k) - ! liquid precipitation (rain) - ice_ocean_boundary%lprec(ig,jg) = x2o_o(ind%x2o_Faxa_rain,k) - ! frozen precipitation (snow) - ice_ocean_boundary%fprec(ig,jg) = x2o_o(ind%x2o_Faxa_snow,k) - ! evaporation flux, MOM6 calls q_flux specific humidity (kg/m2/s) - ice_ocean_boundary%q_flux(ig,jg) = -x2o_o(ind%x2o_Foxx_evap,k) - ! longwave radiation, sum up and down (W/m2) - ice_ocean_boundary%lw_flux(ig,jg) = x2o_o(ind%x2o_Faxa_lwdn,k) + x2o_o(ind%x2o_Foxx_lwup,k) - if (sw_decomp) then - ! Use runtime coefficients to decompose net short-wave heat flux into 4 components - ! 1) visible, direct shortwave (W/m2) - ice_ocean_boundary%sw_flux_vis_dir(ig,jg) = x2o_o(ind%x2o_Foxx_swnet,k)*c1 - ! 2) visible, diffuse shortwave (W/m2) - ice_ocean_boundary%sw_flux_vis_dif(ig,jg) = x2o_o(ind%x2o_Foxx_swnet,k)*c2 - ! 3) near-IR, direct shortwave (W/m2) - ice_ocean_boundary%sw_flux_nir_dir(ig,jg) = x2o_o(ind%x2o_Foxx_swnet,k)*c3 - ! 4) near-IR, diffuse shortwave (W/m2) - ice_ocean_boundary%sw_flux_nir_dif(ig,jg) = x2o_o(ind%x2o_Foxx_swnet,k)*c4 - else - call MOM_error(FATAL,"fill_ice_ocean_bnd: this option has not been implemented yet."// & - "Shortwave must be decomposed using coeffs. c1, c2, c3, c4."); - endif - enddo - enddo - - ice_ocean_boundary%wind_stagger = AGRID - -end subroutine fill_ice_ocean_bnd - !> Step forward ocean model for coupling interval subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) type(ESMF_Clock), intent(inout) :: EClock !< Time and time step ? \todo Why must this be intent(inout)? @@ -1797,14 +1596,11 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o) ! set (actually, get from mct) the cdata pointers: ! \todo this was done in _init_, is it needed again. Does this infodata need to be in glb%? + ! GMM, check if this is needed! call seq_cdata_setptrs(cdata_o, infodata=glb%infodata) - ! fill ice ocean boundary - call fill_ice_ocean_bnd(glb%ice_ocean_boundary, glb%grid, x2o_o%rattr, glb%ind, glb%sw_decomp, & - glb%c1, glb%c2, glb%c3, glb%c4) - - call update_ocean_model(glb%ice_ocean_boundary, glb%ocn_state, glb%ocn_public, & - time_start, coupling_timestep) + call update_ocean_model(glb%ocn_state, glb%ocn_public, time_start, coupling_timestep, & + x2o_o%rattr, glb%ind, glb%sw_decomp, glb%c1, glb%c2, glb%c3, glb%c4) ! return export state to driver call ocn_export(glb%ind, glb%ocn_public, glb%grid, o2x_o%rattr) @@ -1878,39 +1674,32 @@ subroutine forcing_save_restart(CS, G, Time, directory, time_stamped, & end subroutine forcing_save_restart !> Updates the ocean model fields. This code wraps the call to step_MOM with MOM6's call. -!! It uses the forcing in Ice_ocean_boundary to advance the ocean model's state from the +!! It uses the forcing to advance the ocean model's state from the !! input value of Ocean_state (which must be for time time_start_update) for a time interval -!! of Ocean_coupling_time_step, eturning the publicly visible ocean surface properties in +!! of Ocean_coupling_time_step, returning the publicly visible ocean surface properties in !! Ocean_sfc and storing the new ocean properties in Ocean_state. -subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & - time_start_update, Ocean_coupling_time_step) - type(ice_ocean_boundary_type), intent(inout) :: Ice_ocean_boundary - type(ocean_state_type), pointer :: OS - type(ocean_public_type), intent(inout) :: Ocean_sfc - type(time_type), intent(in) :: time_start_update - type(time_type), intent(in) :: Ocean_coupling_time_step - -! Arguments: Ice_ocean_boundary - A structure containing the various forcing -! fields coming from the ice. It is intent in. -! (inout) Ocean_state - A structure containing the internal ocean state. -! (out) Ocean_sfc - A structure containing all the publicly visible ocean -! surface fields after a coupling time step. -! (in) time_start_update - The time at the beginning of the update step. -! (in) Ocean_coupling_time_step - The amount of time over which to advance -! the ocean. - -! Note: although several types are declared intent(inout), this is to allow for -! the possibility of halo updates and to keep previously allocated memory. -! In practice, Ice_ocean_boundary is intent in, Ocean_state is private to -! this module and intent inout, and Ocean_sfc is intent out. - type(time_type) :: Master_time ! This allows step_MOM to temporarily change - ! the time that is seen by internal modules. - type(time_type) :: Time1 ! The value of the ocean model's time at the - ! start of a call to step_MOM. - integer :: index_bnds(4) ! The computational domain index bounds in the - ! ice-ocean boundary type. - real :: weight ! Flux accumulation weight - real :: time_step ! The time step of a call to step_MOM in seconds. +subroutine update_ocean_model(OS, Ocean_sfc, time_start_update, & + Ocean_coupling_time_step, x2o_o, ind, sw_decomp, & + c1, c2, c3, c4) + type(ocean_state_type), pointer :: OS !< Structure containing the internal ocean state + type(ocean_public_type), intent(inout) :: Ocean_sfc !< Structure containing all the publicly + !! visible ocean surface fields after a coupling time step + type(time_type), intent(in) :: time_start_update !< Time at the beginning of the update step + type(time_type), intent(in) :: Ocean_coupling_time_step !< Amount of time over which to + !! advance the ocean + real(kind=8), intent(in) :: x2o_o(:,:) !< Fluxes from coupler to ocean, computed by ocean + type(cpl_indices), intent(inout) :: ind !< Structure with MCT attribute vectors and indices + logical, intent(in) :: sw_decomp !< controls if shortwave is + !!decomposed into four components + real(kind=8), intent(in), optional :: c1, c2, c3, c4 !< Coeffs. used in the shortwave decomposition + + ! local variables + type(time_type) :: Master_time !< This allows step_MOM to temporarily change + !! the time that is seen by internal modules. + type(time_type) :: Time1 !< The value of the ocean model's time at the + !! start of a call to step_MOM. + real :: weight !< Flux accumulation weight + real :: time_step !< The time step of a call to step_MOM in seconds. integer :: secs, days integer :: is, ie, js, je @@ -1936,16 +1725,13 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & call coupler_type_spawn(Ocean_sfc%fields, OS%state%tr_fields, & (/is,is,ie,ie/), (/js,js,je,je/), as_needed=.true.) -! Translate Ice_ocean_boundary into fluxes. - call mpp_get_compute_domain(Ocean_sfc%Domain, index_bnds(1), index_bnds(2), & - index_bnds(3), index_bnds(4)) - weight = 1.0 if (OS%fluxes%fluxes_used) then - call enable_averaging(time_step, OS%Time + Ocean_coupling_time_step, OS%MOM_CSp%diag) ! Needed to allow diagnostics in convert_IOB - call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%fluxes, index_bnds, OS%Time, & - OS%grid, OS%forcing_CSp, OS%state, OS%restore_salinity,OS%restore_temp) + ! GMM, is enable_averaging needed now? + call enable_averaging(time_step, OS%Time + Ocean_coupling_time_step, OS%MOM_CSp%diag) + call ocn_import(OS%fluxes, OS%Time, OS%grid, OS%forcing_CSp, OS%state, x2o_o, ind, sw_decomp, & + c1, c2, c3, c4, OS%restore_salinity,OS%restore_temp) #ifdef _USE_GENERIC_TRACER call MOM_generic_tracer_fluxes_accumulate(OS%fluxes, weight) !here weight=1, just saving the current fluxes #endif @@ -1966,8 +1752,11 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & OS%fluxes%dt_buoy_accum = time_step else OS%flux_tmp%C_p = OS%fluxes%C_p - call convert_IOB_to_fluxes(Ice_ocean_boundary, OS%flux_tmp, index_bnds, OS%Time, & - OS%grid, OS%forcing_CSp, OS%state, OS%restore_salinity,OS%restore_temp) + ! Import fluxes from coupler to ocean. Also, perform do SST and SSS restoring, if needed. + call ocn_import(OS%flux_tmp, OS%Time, OS%grid, OS%forcing_CSp, & + OS%state, x2o_o, ind, sw_decomp, c1, c2, c3, c4, & + OS%restore_salinity,OS%restore_temp) + if (OS%use_ice_shelf) then call shelf_calc_flux(OS%State, OS%flux_tmp, OS%Time, time_step, OS%Ice_shelf_CSp) endif @@ -1978,6 +1767,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & ! call add_berg_flux_to_shelf(OS%grid, OS%flux_tmp, OS%use_ice_shelf,OS%density_iceberg,OS%kv_iceberg, OS%latent_heat_fusion, OS%State, time_step, OS%berg_area_threshold) !endif + ! Accumulate the forcing over time steps call forcing_accumulate(OS%flux_tmp, OS%fluxes, time_step, OS%grid, weight) #ifdef _USE_GENERIC_TRACER call MOM_generic_tracer_fluxes_accumulate(OS%flux_tmp, weight) !weight of the current flux in the running average @@ -2036,35 +1826,30 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, & call callTree_leave("update_ocean_model()") end subroutine update_ocean_model -subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, state, restore_salt, restore_temp) - type(ice_ocean_boundary_type), intent(in), target :: IOB - type(forcing), intent(inout) :: fluxes - integer, dimension(4), intent(in) :: index_bounds - type(time_type), intent(in) :: Time - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure - type(surface_forcing_CS), pointer :: CS - type(surface), intent(in) :: state - logical, optional, intent(in) :: restore_salt, restore_temp - -! This subroutine translates the Ice_ocean_boundary_type into a -! MOM forcing type, including changes of units, sign conventions, -! and puting the fields into arrays with MOM-standard halos. - -! Arguments: -! IOB ice-ocean boundary type w/ fluxes to drive ocean in a coupled model -! (out) fluxes - A structure containing pointers to any possible -! forcing fields. Unused fields have NULL ptrs. -! (in) index_bounds - the i- and j- size of the arrays in IOB. -! (in) Time - The time of the fluxes, used for interpolating the salinity -! to the right time, when it is being restored. -! (in) G - The ocean's grid structure. -! (in) CS - A pointer to the control structure returned by a previous -! call to surface_forcing_init. -! (in) state - A structure containing fields that describe the -! surface state of the ocean. -! (in) restore_salt - if true, salinity is restored to a target value. -! (in) restore_temp - if true, temperature is restored to a target value. +!> This function has a few purposes: 1) it allocates and initializes the data +!! in the fluxes structure; 2) it imports surface fluxes using data from +!! the coupler; and 3) it can apply restoring in SST and SSS. +!! See \ref section_ocn_import for a summary of the surface fluxes that are +!! passed from MCT to MOM6, including fluxes that need to be included in +!! the future. +subroutine ocn_import(fluxes, Time, G, CS, state, x2o_o, ind, sw_decomp, & + c1, c2, c3, c4, restore_salt, restore_temp) + type(forcing), intent(inout) :: fluxes !< Surface fluxes structure + type(time_type), intent(in) :: Time !< Model time structure + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + type(surface_forcing_CS), pointer :: CS !< control structure returned by + !! a previous call to surface_forcing_init + type(surface), intent(in) :: state !< control structure to ocean + !! surface state fields. + real(kind=8), intent(in) :: x2o_o(:,:)!< Fluxes from coupler to ocean, computed by ocean + type(cpl_indices), intent(inout) :: ind !< Structure with MCT attribute vectors and indices + logical, intent(in) :: sw_decomp !< controls if shortwave is + !!decomposed into four components + real(kind=8), intent(in), optional :: c1, c2, c3, c4 !< Coeffs. used in the shortwave decomposition + logical, optional, intent(in) :: restore_salt, restore_temp !< Controls if salt and temp are + !! restored + ! local variables real, dimension(SZIB_(G),SZJB_(G)) :: & taux_at_q, & ! Zonal wind stresses at q points (Pa) tauy_at_q ! Meridional wind stresses at q points (Pa) @@ -2095,7 +1880,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, state, real :: mass_eff ! effective mass of sea ice for rigidity (kg/m^2) integer :: wind_stagger ! AGRID, BGRID_NE, or CGRID_NE (integers from MOM_domains) - integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0 + integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isr, ier, jsr, jer integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd @@ -2110,8 +1895,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, state, call cpu_clock_begin(id_clock_forcing) - isc_bnd = index_bounds(1) ; iec_bnd = index_bounds(2) - jsc_bnd = index_bounds(3) ; jec_bnd = index_bounds(4) is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -2134,8 +1917,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, state, restore_sst = .false. if (present(restore_temp)) restore_sst = restore_temp - ! allocation and initialization if this is the first time that this - ! flux type has been used. + ! if true, allocation and initialization if (fluxes%dt_buoy_accum < 0) then call allocate_forcing_type(G, fluxes, stress=.true., ustar=.true., & water=.true., heat=.true.) @@ -2147,6 +1929,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, state, call safe_alloc_ptr(fluxes%p_surf,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%p_surf_full,isd,ied,jsd,jed) + call safe_alloc_ptr(fluxes%p_surf_SSH,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%salt_flux,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%salt_flux_in,isd,ied,jsd,jed) @@ -2175,21 +1958,11 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, state, fluxes%dt_buoy_accum = 0.0 endif ! endif for allocation and initialization - if ((.not.coupler_type_initialized(fluxes%tr_fluxes)) .and. & - coupler_type_initialized(IOB%fluxes)) & - call coupler_type_spawn(IOB%fluxes, fluxes%tr_fluxes, & - (/is,is,ie,ie/), (/js,js,je,je/)) - ! It might prove valuable to use the same array extents as the rest of the - ! ocean model, rather than using haloless arrays, in which case the last line - ! would be: ( (/isd,is,ie,ied/), (/jsd,js,je,jed/)) - - if (CS%allow_flux_adjustments) then fluxes%heat_added(:,:)=0.0 fluxes%salt_flux_added(:,:)=0.0 endif - ! allocation and initialization on first call to this routine if (CS%area_surf < 0.0) then do j=js,je ; do i=is,ie work_sum(i,j) = G%areaT(i,j) * G%mask2dT(i,j) @@ -2265,9 +2038,9 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, state, enddo; enddo endif - wind_stagger = CS%wind_stagger - if ((IOB%wind_stagger == AGRID) .or. (IOB%wind_stagger == BGRID_NE) .or. & - (IOB%wind_stagger == CGRID_NE)) wind_stagger = IOB%wind_stagger + ! GMM, CIME uses AGRID. All the BGRID_NE code can be cleaned later + wind_stagger = AGRID + if (wind_stagger == BGRID_NE) then ! This is necessary to fill in the halo points. taux_at_q(:,:) = 0.0 ; tauy_at_q(:,:) = 0.0 @@ -2277,113 +2050,124 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, state, taux_at_h(:,:) = 0.0 ; tauy_at_h(:,:) = 0.0 endif - - ! obtain fluxes from IOB; note the staggering of indices - i0 = is - isc_bnd ; j0 = js - jsc_bnd + k = 0 do j=js,je ; do i=is,ie + k = k + 1 ! Increment position within gindex if (wind_stagger == BGRID_NE) then - if (ASSOCIATED(IOB%u_flux)) taux_at_q(I,J) = IOB%u_flux(i-i0,j-j0) * CS%wind_stress_multiplier - if (ASSOCIATED(IOB%v_flux)) tauy_at_q(I,J) = IOB%v_flux(i-i0,j-j0) * CS%wind_stress_multiplier + taux_at_q(I,J) = x2o_o(ind%x2o_Foxx_taux,k) * CS%wind_stress_multiplier + tauy_at_q(I,J) = x2o_o(ind%x2o_Foxx_tauy,k) * CS%wind_stress_multiplier + ! GMM, cime uses AGRID elseif (wind_stagger == AGRID) then - if (ASSOCIATED(IOB%u_flux)) taux_at_h(i,j) = IOB%u_flux(i-i0,j-j0) * CS%wind_stress_multiplier - if (ASSOCIATED(IOB%v_flux)) tauy_at_h(i,j) = IOB%v_flux(i-i0,j-j0) * CS%wind_stress_multiplier + taux_at_h(i,j) = x2o_o(ind%x2o_Foxx_taux,k) * CS%wind_stress_multiplier + tauy_at_h(i,j) = x2o_o(ind%x2o_Foxx_tauy,k) * CS%wind_stress_multiplier else ! C-grid wind stresses. - if (ASSOCIATED(IOB%u_flux)) fluxes%taux(I,j) = IOB%u_flux(i-i0,j-j0) * CS%wind_stress_multiplier - if (ASSOCIATED(IOB%v_flux)) fluxes%tauy(i,J) = IOB%v_flux(i-i0,j-j0) * CS%wind_stress_multiplier + fluxes%taux(I,j) = x2o_o(ind%x2o_Foxx_taux,k) * CS%wind_stress_multiplier + fluxes%tauy(i,J) = x2o_o(ind%x2o_Foxx_tauy,k) * CS%wind_stress_multiplier endif - if (ASSOCIATED(IOB%lprec)) & - fluxes%lprec(i,j) = IOB%lprec(i-i0,j-j0) * G%mask2dT(i,j) - - if (ASSOCIATED(IOB%fprec)) & - fluxes%fprec(i,j) = IOB%fprec(i-i0,j-j0) * G%mask2dT(i,j) - - if (ASSOCIATED(IOB%q_flux)) & - fluxes%evap(i,j) = - IOB%q_flux(i-i0,j-j0) * G%mask2dT(i,j) - - if (ASSOCIATED(IOB%runoff)) & - fluxes%lrunoff(i,j) = IOB%runoff(i-i0,j-j0) * G%mask2dT(i,j) - - if (ASSOCIATED(IOB%calving)) & - fluxes%frunoff(i,j) = IOB%calving(i-i0,j-j0) * G%mask2dT(i,j) - - if (((ASSOCIATED(IOB%ustar_berg) .and. (.not. ASSOCIATED(fluxes%ustar_berg))) & - .or. (ASSOCIATED(IOB%area_berg) .and. (.not. ASSOCIATED(fluxes%area_berg)))) & - .or. (ASSOCIATED(IOB%mass_berg) .and. (.not. ASSOCIATED(fluxes%mass_berg)))) & - call allocate_forcing_type(G, fluxes, iceberg=.true.) - - if (ASSOCIATED(IOB%ustar_berg)) & - fluxes%ustar_berg(i,j) = IOB%ustar_berg(i-i0,j-j0) * G%mask2dT(i,j) - - if (ASSOCIATED(IOB%area_berg)) & - fluxes%area_berg(i,j) = IOB%area_berg(i-i0,j-j0) * G%mask2dT(i,j) - - if (ASSOCIATED(IOB%mass_berg)) & - fluxes%mass_berg(i,j) = IOB%mass_berg(i-i0,j-j0) * G%mask2dT(i,j) - - if (ASSOCIATED(IOB%runoff_hflx)) & - fluxes%heat_content_lrunoff(i,j) = IOB%runoff_hflx(i-i0,j-j0) * G%mask2dT(i,j) - - if (ASSOCIATED(IOB%calving_hflx)) & - fluxes%heat_content_frunoff(i,j) = IOB%calving_hflx(i-i0,j-j0) * G%mask2dT(i,j) - - if (ASSOCIATED(IOB%lw_flux)) & - fluxes%LW(i,j) = IOB%lw_flux(i-i0,j-j0) * G%mask2dT(i,j) - - if (ASSOCIATED(IOB%t_flux)) & - fluxes%sens(i,j) = - IOB%t_flux(i-i0,j-j0) * G%mask2dT(i,j) - - fluxes%latent(i,j) = 0.0 - if (ASSOCIATED(IOB%latent)) then - fluxes%latent(i,j) = IOB%latent(i-i0,j-j0) + ! liquid precipitation (rain) + if (ASSOCIATED(fluxes%lprec)) & + fluxes%lprec(i,j) = x2o_o(ind%x2o_Faxa_rain,k) * G%mask2dT(i,j) + + ! frozen precipitation (snow) + if (ASSOCIATED(fluxes%fprec)) & + fluxes%fprec(i,j) = x2o_o(ind%x2o_Faxa_snow,k) * G%mask2dT(i,j) + + ! evaporation + if (ASSOCIATED(fluxes%evap)) & + fluxes%evap(i,j) = x2o_o(ind%x2o_Foxx_evap,k) * G%mask2dT(i,j) + + ! river runoff flux + if (ASSOCIATED(fluxes%lrunoff)) & + fluxes%lrunoff(i,j) = x2o_o(ind%x2o_Foxx_rofl,k) * G%mask2dT(i,j) + + ! ice runoff flux + if (ASSOCIATED(fluxes%frunoff)) & + fluxes%frunoff(i,j) = x2o_o(ind%x2o_Foxx_rofi,k) * G%mask2dT(i,j) + + ! GMM, we don't have an icebergs yet so the following is not needed + !if (((ASSOCIATED(IOB%ustar_berg) .and. (.not. ASSOCIATED(fluxes%ustar_berg))) & + ! .or. (ASSOCIATED(IOB%area_berg) .and. (.not. ASSOCIATED(fluxes%area_berg)))) & + ! .or. (ASSOCIATED(IOB%mass_berg) .and. (.not. ASSOCIATED(fluxes%mass_berg)))) & + ! call allocate_forcing_type(G, fluxes, iceberg=.true.) + !if (ASSOCIATED(IOB%ustar_berg)) & + ! fluxes%ustar_berg(i,j) = IOB%ustar_berg(i-i0,j-j0) * G%mask2dT(i,j) + !if (ASSOCIATED(IOB%area_berg)) & + ! fluxes%area_berg(i,j) = IOB%area_berg(i-i0,j-j0) * G%mask2dT(i,j) + !if (ASSOCIATED(IOB%mass_berg)) & + ! fluxes%mass_berg(i,j) = IOB%mass_berg(i-i0,j-j0) * G%mask2dT(i,j) + + ! GMM, cime does not not have an equivalent for heat_content_lrunoff and + ! heat_content_frunoff. I am seeting these to zero for now. + if (ASSOCIATED(fluxes%heat_content_lrunoff)) & + fluxes%heat_content_lrunoff(i,j) = 0.0 * G%mask2dT(i,j) + + if (ASSOCIATED(fluxes%heat_content_frunoff)) & + fluxes%heat_content_frunoff(i,j) = 0.0 * G%mask2dT(i,j) + + ! longwave radiation, sum up and down (W/m2) + if (ASSOCIATED(fluxes%LW)) & + fluxes%LW(i,j) = (x2o_o(ind%x2o_Faxa_lwdn,k) + x2o_o(ind%x2o_Foxx_lwup,k)) * G%mask2dT(i,j) + + ! sensible heat flux (W/m2) + if (ASSOCIATED(fluxes%sens)) & + fluxes%sens(i,j) = x2o_o(ind%x2o_Foxx_sen,k) * G%mask2dT(i,j) + + ! latent heat flux (W/m^2) + if (ASSOCIATED(fluxes%latent)) & + fluxes%latent(i,j) = x2o_o(ind%x2o_Foxx_lat,k) * G%mask2dT(i,j) + + if (sw_decomp) then + ! Use runtime coefficients to decompose net short-wave heat flux into 4 components + ! 1) visible, direct shortwave (W/m2) + if (ASSOCIATED(fluxes%sw_vis_dir)) & + fluxes%sw_vis_dir(i,j) = G%mask2dT(i,j) * x2o_o(ind%x2o_Foxx_swnet,k)*c1 + ! 2) visible, diffuse shortwave (W/m2) + if (ASSOCIATED(fluxes%sw_vis_dif)) & + fluxes%sw_vis_dif(i,j) = G%mask2dT(i,j) * x2o_o(ind%x2o_Foxx_swnet,k)*c2 + ! 3) near-IR, direct shortwave (W/m2) + if (ASSOCIATED(fluxes%sw_nir_dir)) & + fluxes%sw_nir_dir(i,j) = G%mask2dT(i,j) * x2o_o(ind%x2o_Foxx_swnet,k)*c3 + ! 4) near-IR, diffuse shortwave (W/m2) + if (ASSOCIATED(fluxes%sw_nir_dif)) & + fluxes%sw_nir_dif(i,j) = G%mask2dT(i,j) * x2o_o(ind%x2o_Foxx_swnet,k)*c4 + + fluxes%sw(i,j) = fluxes%sw_vis_dir(i,j) + fluxes%sw_vis_dif(i,j) + & + fluxes%sw_nir_dir(i,j) + fluxes%sw_nir_dif(i,j) else - if (ASSOCIATED(IOB%fprec)) then - fluxes%latent(i,j) = fluxes%latent(i,j) - IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion - fluxes%latent_fprec_diag(i,j) = -G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion - endif - if (ASSOCIATED(IOB%calving)) then - fluxes%latent(i,j) = fluxes%latent(i,j) - IOB%calving(i-i0,j-j0)*CS%latent_heat_fusion - fluxes%latent_frunoff_diag(i,j) = -G%mask2dT(i,j) * IOB%calving(i-i0,j-j0)*CS%latent_heat_fusion + call MOM_error(FATAL,"fill_ice_ocean_bnd: this option has not been implemented yet."// & + "Shortwave must be decomposed using coeffs. c1, c2, c3, c4."); + endif + + ! applied surface pressure from atmosphere and cryosphere + ! sea-level pressure (Pa) + if (ASSOCIATED(fluxes%p_surf_full) .and. ASSOCIATED(fluxes%p_surf)) then + fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * x2o_o(ind%x2o_Sa_pslv,k) + if (CS%max_p_surf >= 0.0) then + fluxes%p_surf(i,j) = MIN(fluxes%p_surf_full(i,j),CS%max_p_surf) + else + fluxes%p_surf(i,j) = fluxes%p_surf_full(i,j) endif - if (ASSOCIATED(IOB%q_flux)) then - fluxes%latent(i,j) = fluxes%latent(i,j) - IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor - fluxes%latent_evap_diag(i,j) = -G%mask2dT(i,j) * IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor + + if (CS%use_limited_P_SSH) then + fluxes%p_surf_SSH(i,j) = fluxes%p_surf(i,j) + else + fluxes%p_surf_SSH(i,j) = fluxes%p_surf_full(i,j) endif - endif - fluxes%latent(i,j) = G%mask2dT(i,j) * fluxes%latent(i,j) - - if (ASSOCIATED(IOB%sw_flux_vis_dir)) & - fluxes%sw_vis_dir(i,j) = G%mask2dT(i,j) * IOB%sw_flux_vis_dir(i-i0,j-j0) - if (ASSOCIATED(IOB%sw_flux_vis_dif)) & - fluxes%sw_vis_dif(i,j) = G%mask2dT(i,j) * IOB%sw_flux_vis_dif(i-i0,j-j0) - if (ASSOCIATED(IOB%sw_flux_nir_dir)) & - fluxes%sw_nir_dir(i,j) = G%mask2dT(i,j) * IOB%sw_flux_nir_dir(i-i0,j-j0) - if (ASSOCIATED(IOB%sw_flux_nir_dif)) & - fluxes%sw_nir_dif(i,j) = G%mask2dT(i,j) * IOB%sw_flux_nir_dif(i-i0,j-j0) - fluxes%sw(i,j) = fluxes%sw_vis_dir(i,j) + fluxes%sw_vis_dif(i,j) + & - fluxes%sw_nir_dir(i,j) + fluxes%sw_nir_dif(i,j) + endif - enddo ; enddo + ! salt flux + ! more salt restoring logic + if (ASSOCIATED(fluxes%salt_flux)) & + fluxes%salt_flux(i,j) = G%mask2dT(i,j)*(x2o_o(ind%x2o_Fioi_salt,k) + fluxes%salt_flux(i,j)) - ! more salt restoring logic - if (ASSOCIATED(IOB%salt_flux)) then - do j=js,je ; do i=is,ie - fluxes%salt_flux(i,j) = G%mask2dT(i,j)*(fluxes%salt_flux(i,j) - IOB%salt_flux(i-i0,j-j0)) - fluxes%salt_flux_in(i,j) = G%mask2dT(i,j)*( -IOB%salt_flux(i-i0,j-j0) ) - enddo ; enddo - endif + if (ASSOCIATED(fluxes%salt_flux_in)) & + fluxes%salt_flux_in(i,j) = G%mask2dT(i,j)*x2o_o(ind%x2o_Fioi_salt,k) -!### if (associated(CS%ctrl_forcing_CSp)) then -!### do j=js,je ; do i=is,ie -!### SST_anom(i,j) = state%SST(i,j) - CS%T_Restore(i,j) -!### SSS_anom(i,j) = state%SSS(i,j) - CS%S_Restore(i,j) -!### SSS_mean(i,j) = 0.5*(state%SSS(i,j) + CS%S_Restore(i,j)) -!### enddo ; enddo -!### call apply_ctrl_forcing(SST_anom, SSS_anom, SSS_mean, fluxes%heat_restore, & -!### fluxes%vprec, day, dt, G, CS%ctrl_forcing_CSp) -!### endif + enddo ; enddo + ! ############################ END OF MCT to MOM ############################## ! adjust the NET fresh-water flux to zero, if flagged if (CS%adjust_net_fresh_water_to_zero) then @@ -2397,9 +2181,9 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, state, ! Bob thinks this is trying ensure the net fresh-water of the ocean + sea-ice system ! is constant. ! To do this correctly we will need a sea-ice melt field added to IOB. -AJA - if (ASSOCIATED(IOB%salt_flux) .and. (CS%ice_salt_concentration>0.0)) & + if (ASSOCIATED(fluxes%salt_flux) .and. (CS%ice_salt_concentration>0.0)) & net_FW(i,j) = net_FW(i,j) - G%areaT(i,j) * & - (IOB%salt_flux(i-i0,j-j0) / CS%ice_salt_concentration) + (fluxes%salt_flux(i,j) / CS%ice_salt_concentration) net_FW2(i,j) = net_FW(i,j) enddo ; enddo @@ -2417,26 +2201,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, state, endif - ! applied surface pressure from atmosphere and cryosphere - if (ASSOCIATED(IOB%p)) then - if (CS%max_p_surf >= 0.0) then - do j=js,je ; do i=is,ie - fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0) - fluxes%p_surf(i,j) = MIN(fluxes%p_surf_full(i,j),CS%max_p_surf) - enddo ; enddo - else - do j=js,je ; do i=is,ie - fluxes%p_surf_full(i,j) = G%mask2dT(i,j) * IOB%p(i-i0,j-j0) - fluxes%p_surf(i,j) = fluxes%p_surf_full(i,j) - enddo ; enddo - endif - if (CS%use_limited_P_SSH) then - fluxes%p_surf_SSH => fluxes%p_surf - else - fluxes%p_surf_SSH => fluxes%p_surf_full - endif - endif - ! surface momentum stress related fields as function of staggering if (wind_stagger == BGRID_NE) then if (G%symmetric) & @@ -2559,10 +2323,6 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, state, enddo ; enddo endif - if (coupler_type_initialized(fluxes%tr_fluxes) .and. & - coupler_type_initialized(IOB%fluxes)) & - call coupler_type_copy_data(IOB%fluxes, fluxes%tr_fluxes) - if (CS%allow_flux_adjustments) then ! Apply adjustments to fluxes call apply_flux_adjustments(G, CS, Time, fluxes) @@ -2572,7 +2332,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, state, call user_alter_forcing(state, fluxes, Time, G, CS%urf_CS) call cpu_clock_end(id_clock_forcing) -end subroutine convert_IOB_to_fluxes + +end subroutine ocn_import !> Adds flux adjustments obtained via data_override !! Component name is 'OCN' @@ -2584,6 +2345,7 @@ subroutine apply_flux_adjustments(G, CS, Time, fluxes) type(surface_forcing_CS), pointer :: CS !< Surface forcing control structure type(time_type), intent(in) :: Time !< Model time structure type(forcing), optional, intent(inout) :: fluxes !< Surface fluxes structure + ! Local variables real, dimension(SZI_(G),SZJ_(G)) :: tempx_at_h ! Delta to zonal wind stress at h points (Pa) real, dimension(SZI_(G),SZJ_(G)) :: tempy_at_h ! Delta to meridional wind stress at h points (Pa) @@ -2689,12 +2451,12 @@ subroutine ocean_model_end(Ocean_sfc, Ocean_state, Time) ! !! ocean state to be deallocated upon termination. type(time_type), intent(in) :: Time !< The model time, used for writing restarts. - if (debug .and. is_root_pe()) write(glb%stdout,*)'Here 1' + !if (debug .and. is_root_pe()) write(glb%stdout,*)'Here 1' !GMM call save_restart(Ocean_state, Time) call diag_mediator_end(Time, Ocean_state%MOM_CSp%diag) call MOM_end(Ocean_state%MOM_CSp) if (Ocean_state%use_ice_shelf) call ice_shelf_end(Ocean_state%Ice_shelf_CSp) - if (debug .and. is_root_pe()) write(glb%stdout,*)'Here 2' + !if (debug .and. is_root_pe()) write(glb%stdout,*)'Here 2' end subroutine ocean_model_end @@ -2826,4 +2588,69 @@ subroutine ocn_domain_mct( lsize, gsMap_ocn, dom_ocn) end subroutine ocn_domain_mct +!> \namespace ocn_comp_mct +!! +!! \section section_ocn_import Fluxes imported from the coupler (MCT) to MOM6 +!! The following summarizes the mismatches between MCT and MOM6 in terms +!! of ice ocean fluxes. +!! +!! Redundancies: +!! x2o_Faxa_prec = x2o_Faxa_rain + x2o_Faxa_snow +!! +!! Variables whose units and sign **could not** be verified so far: +!! x2o_Foxx_rofl +!! x2o_Foxx_rof +!! +!! Variables in MOM6 fluxes that are **NOT** filled by the coupler: +!! ustar_berg, frictional velocity beneath icebergs (m/s) +!! area_berg, area covered by icebergs(m2/m2) +!! mass_berg, mass of icebergs(kg/m2) +!! runoff_hflx, heat content of liquid runoff (W/m2) +!! calving_hflx, heat content of frozen runoff (W/m2) +!! mi, mass of ice (kg/m2) +!! +!! Variables in the coupler that are **NOT** used in MOM6 (i.e., no corresponding field in fluxes): +!! x2o_Fioi_melth, heat flux from snow & ice melt (W/m2) +!! x2o_Fioi_meltw, snow melt flux (kg/m2/s) +!! x2o_Si_ifrac, fractional ice wrt ocean +!! x2o_So_duu10n, 10m wind speed squared (m^2/s^2) +!! x2o_Sa_co2prog, bottom atm level prognostic CO2 +!! x2o_Sa_co2diag, bottom atm level diagnostic CO2 +!! +!! \TODO Langmuir related fields: +!! surface Stokes drift, x-comp. (x2o_Sw_ustokes) +!! surface Stokes drift, y-comp. (x2o_Sw_vstokes) +!! wave model langmuir multiplier (x2o_Sw_lamult) +!! +!! \TODO Biogeochemistry: +!! x2o_Fioi_bcpho, Black Carbon hydrophobic release from sea ice component +!! x2o_Fioi_bcphi, Black Carbon hydrophilic release from sea ice component +!! x2o_Fioi_flxdst, Dust release from sea ice component +!! x2o_Faxa_bcphidry, Black Carbon hydrophilic dry deposition +!! x2o_Faxa_bcphodry, Black Carbon hydrophobic dry deposition +!! x2o_Faxa_bcphiwet, Black Carbon hydrophobic wet deposition +!! x2o_Faxa_ocphidry, Organic Carbon hydrophilic dry deposition +!! x2o_Faxa_ocphodry, Organic Carbon hydrophobic dry deposition +!! x2o_Faxa_ocphiwet, Organic Carbon hydrophilic dry deposition +!! x2o_Faxa_dstwet, Sizes 1 to 4 dust - wet deposition +!! x2o_Faxa_dstdry, Sizes 1 to 4 dust - dry deposition +!! +!! \section section_ocn_export Fluxes exported from MOM6 to the coupler (MCT) +!! +!! Variables that are currently being exported: +!! +!! Surface temperature (Kelvin) +!! Surface salinity (psu) +!! Surface eastward velocity (m/s) +!! Surface northward velocity (m/s) +!! Zonal slope in the sea surface height +!! Meridional slope in the sea surface height +!! +!! \TODO Variables that **are not** currently being exported: +!! +!! Boundary layer depth +!! CO2 +!! DMS +!! o2x_Fioo_q !< Heat flux? + end module ocn_comp_mct From 0e92dfa405b35906114c63d2398a4d2d1e88018b Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 17 Oct 2017 16:10:58 -0600 Subject: [PATCH 20/20] Delete unecessary changes I had made a few small changes in MOM_surface_forcing and MOM_restart files that are no longer necessary. To avoid confusion when merging with dev/master, these changes were deleted. These files are now identical to the latest version in dev/master. --- .../coupled_driver/MOM_surface_forcing.F90 | 27 ++++++++----------- src/framework/MOM_restart.F90 | 2 +- 2 files changed, 12 insertions(+), 17 deletions(-) diff --git a/config_src/coupled_driver/MOM_surface_forcing.F90 b/config_src/coupled_driver/MOM_surface_forcing.F90 index 46fdb2232b..178a122087 100644 --- a/config_src/coupled_driver/MOM_surface_forcing.F90 +++ b/config_src/coupled_driver/MOM_surface_forcing.F90 @@ -152,7 +152,6 @@ module MOM_surface_forcing real, pointer, dimension(:,:) :: sw_flux_vis_dif =>NULL() ! diffuse visible sw radiation (W/m2) real, pointer, dimension(:,:) :: sw_flux_nir_dir =>NULL() ! direct Near InfraRed sw radiation (W/m2) real, pointer, dimension(:,:) :: sw_flux_nir_dif =>NULL() ! diffuse Near InfraRed sw radiation (W/m2) - real, pointer, dimension(:,:) :: latent =>NULL() ! latent heat flux (W/m2) real, pointer, dimension(:,:) :: lprec =>NULL() ! mass flux of liquid precip (kg/m2/s) real, pointer, dimension(:,:) :: fprec =>NULL() ! mass flux of frozen precip (kg/m2/s) real, pointer, dimension(:,:) :: runoff =>NULL() ! mass flux of liquid runoff (kg/m2/s) @@ -481,21 +480,17 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, G, CS, state, fluxes%sens(i,j) = - IOB%t_flux(i-i0,j-j0) * G%mask2dT(i,j) fluxes%latent(i,j) = 0.0 - if (ASSOCIATED(IOB%latent)) then - fluxes%latent(i,j) = IOB%latent(i-i0,j-j0) - else - if (ASSOCIATED(IOB%fprec)) then - fluxes%latent(i,j) = fluxes%latent(i,j) - IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion - fluxes%latent_fprec_diag(i,j) = -G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion - endif - if (ASSOCIATED(IOB%calving)) then - fluxes%latent(i,j) = fluxes%latent(i,j) - IOB%calving(i-i0,j-j0)*CS%latent_heat_fusion - fluxes%latent_frunoff_diag(i,j) = -G%mask2dT(i,j) * IOB%calving(i-i0,j-j0)*CS%latent_heat_fusion - endif - if (ASSOCIATED(IOB%q_flux)) then - fluxes%latent(i,j) = fluxes%latent(i,j) - IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor - fluxes%latent_evap_diag(i,j) = -G%mask2dT(i,j) * IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor - endif + if (ASSOCIATED(IOB%fprec)) then + fluxes%latent(i,j) = fluxes%latent(i,j) - IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion + fluxes%latent_fprec_diag(i,j) = -G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*CS%latent_heat_fusion + endif + if (ASSOCIATED(IOB%calving)) then + fluxes%latent(i,j) = fluxes%latent(i,j) - IOB%calving(i-i0,j-j0)*CS%latent_heat_fusion + fluxes%latent_frunoff_diag(i,j) = -G%mask2dT(i,j) * IOB%calving(i-i0,j-j0)*CS%latent_heat_fusion + endif + if (ASSOCIATED(IOB%q_flux)) then + fluxes%latent(i,j) = fluxes%latent(i,j) - IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor + fluxes%latent_evap_diag(i,j) = -G%mask2dT(i,j) * IOB%q_flux(i-i0,j-j0)*CS%latent_heat_vapor endif fluxes%latent(i,j) = G%mask2dT(i,j) * fluxes%latent(i,j) diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index f538d3f3e4..42c41a59a4 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -889,7 +889,7 @@ subroutine restore_state(filename, directory, day, G, CS) ! restart_init. character(len=200) :: filepath ! The path (dir/file) to the file being opened. - character(len=200) :: fname ! The name of the current file. + character(len=80) :: fname ! The name of the current file. character(len=8) :: suffix ! A suffix (like "_2") that is added to any ! additional restart files. character(len=256) :: mesg ! A message for warnings.