diff --git a/generic_tracers/MOM_generic_tracer.F90 b/generic_tracers/MOM_generic_tracer.F90 new file mode 100644 index 00000000..8e81f095 --- /dev/null +++ b/generic_tracers/MOM_generic_tracer.F90 @@ -0,0 +1,1160 @@ +!> Drives the generic version of tracers TOPAZ and CFC and other GFDL BGC components +module MOM_generic_tracer + +! This file is part of MOM6. See LICENSE.md for the license. + +#include + +! The following macro is usually defined in but since MOM6 should not directly +! include files from FMS we replicate the macro lines here: +#ifdef NO_F2000 +#define _ALLOCATED associated +#else +#define _ALLOCATED allocated +#endif + +! ### These imports should not reach into FMS directly ### +use field_manager_mod, only: fm_string_len + +use generic_tracer, only: generic_tracer_register, generic_tracer_get_diag_list +use generic_tracer, only: generic_tracer_init, generic_tracer_source, generic_tracer_register_diag +use generic_tracer, only: generic_tracer_coupler_get, generic_tracer_coupler_set +use generic_tracer, only: generic_tracer_end, generic_tracer_get_list, do_generic_tracer +use generic_tracer, only: generic_tracer_update_from_bottom,generic_tracer_vertdiff_G +use generic_tracer, only: generic_tracer_coupler_accumulate + +use g_tracer_utils, only: g_tracer_get_name,g_tracer_set_values,g_tracer_set_common,g_tracer_get_common +use g_tracer_utils, only: g_tracer_get_next,g_tracer_type,g_tracer_is_prog,g_tracer_flux_init +use g_tracer_utils, only: g_tracer_send_diag,g_tracer_get_values +use g_tracer_utils, only: g_tracer_get_pointer,g_tracer_get_alias,g_tracer_set_csdiag +use g_tracer_utils, only: g_tracer_get_obc_segment_props + +use MOM_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS +use MOM_ALE_sponge, only : ALE_sponge_CS, initialize_ALE_sponge +use MOM_coms, only : EFP_type, max_across_PEs, min_across_PEs, PE_here +use MOM_diagnose_mld, only : diagnoseMLDbyDensityDifference, diagnoseMLDbyEnergy +use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr +use MOM_diag_mediator, only : diag_ctrl, get_diag_time_end +use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, is_root_pe +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_forcing_type, only : forcing, optics_type +use MOM_grid, only : ocean_grid_type +use MOM_hor_index, only : hor_index_type +use MOM_interface_heights, only : thickness_to_dz +use MOM_io, only : file_exists, MOM_read_data, slasher +use MOM_open_boundary, only : ocean_OBC_type +use MOM_open_boundary, only : register_obgc_segments, fill_obgc_segments +use MOM_open_boundary, only : set_obgc_segments_props +use MOM_restart, only : register_restart_field, query_initialized, set_initialized, MOM_restart_CS +use MOM_spatial_means, only : global_area_mean, global_mass_int_EFP, array_global_min_max +use MOM_sponge, only : set_up_sponge_field, sponge_CS +use MOM_time_manager, only : time_type, set_time +use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut +use MOM_tracer_registry, only : register_tracer, tracer_registry_type +use MOM_tracer_Z_init, only : tracer_Z_init +use MOM_tracer_initialization_from_Z, only : MOM_initialize_tracer_from_Z +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : surface, thermo_var_ptrs +use MOM_verticalGrid, only : verticalGrid_type + + +implicit none ; private + +!> A state hidden in module data that is very much not allowed in MOM6 +! ### This needs to be fixed +logical :: g_registered = .false. + +public register_MOM_generic_tracer, initialize_MOM_generic_tracer +public MOM_generic_tracer_column_physics, MOM_generic_tracer_surface_state +public end_MOM_generic_tracer, MOM_generic_tracer_get +public MOM_generic_tracer_stock +public MOM_generic_flux_init +public MOM_generic_tracer_min_max +public MOM_generic_tracer_fluxes_accumulate +public register_MOM_generic_tracer_segments + +!> Control structure for generic tracers +type, public :: MOM_generic_tracer_CS ; private + character(len = 200) :: IC_file !< The file in which the generic tracer initial values can + !! be found, or an empty string for internal initialization. + logical :: Z_IC_file !< If true, the generic_tracer IC_file is in Z-space. The default is false. + real :: tracer_IC_val = 0.0 !< The initial value assigned to tracers, in + !! concentration units [conc] + real :: tracer_land_val = -1.0 !< The values of tracers used where land is masked out, in + !! concentration units [conc] + logical :: tracers_may_reinit !< If true, tracers may go through the + !! initialization code if they are not found in the restart files. + + logical :: mld_pha_calc = .False. !< If true, use a fixed value for photoacclimation MLD + real :: mld_pha_val = 0.0 !< The value of fixed photoacclimation MLD + logical :: mld_pha_use_delta_rho = .False. !< If true, use a density diference to find the MLD + real :: mld_pha_href = 0.0 !< The reference depth for density difference based MLD + real :: mld_pha_drho = 0.03 !< The density thershold for a density difference based MLD + logical :: mld_pha_use_delta_eng = .False. !< If true, use an energy diference to find the MLD + real :: mld_pha_deng = 25.0 !< The energy threshold for an energy d ifference based MLD + + type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to + !! regulate the timing of diagnostic output. + type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< Restart control structure + type(ocean_OBC_type), pointer :: OBC => NULL() ! Pointer to the first element of the linked list of generic tracers. + type(g_tracer_type), pointer :: g_tracer_list => NULL() + +end type MOM_generic_tracer_CS + +contains + +!> Initializes the generic tracer packages and adds their tracers to the list +!! Adds the tracers in the list of generic tracers to the set of MOM tracers (i.e., MOM-register them) +!! Register these tracers for restart +function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS) + type(hor_index_type), intent(in) :: HI !< Horizontal index ranges + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module + type(tracer_registry_type), pointer :: tr_Reg !< Pointer to the control structure for the tracer + !! advection and diffusion module. + type(MOM_restart_CS), target, intent(inout) :: restart_CS !< MOM restart control struct + + ! Local variables + logical :: register_MOM_generic_tracer + logical :: obc_has + ! This include declares and sets the variable "version". +# include "version_variable.h" + + character(len=128), parameter :: sub_name = 'register_MOM_generic_tracer' + character(len=200) :: inputdir ! The directory where NetCDF input files are. + ! These can be overridden later in via the field manager? + + integer :: ntau, axes(3) + type(g_tracer_type), pointer :: g_tracer, g_tracer_next + character(len=fm_string_len) :: g_tracer_name, longname,units + character(len=fm_string_len) :: obc_src_file_name, obc_src_field_name + real :: lfac_in ! Multiplicative factor used in setting the tracer-specific inverse length + ! scales associated with inflowing tracer reservoirs at OBCs [nondim] + real :: lfac_out ! Multiplicative factor used in setting the tracer-specific inverse length + ! scales associated with outflowing tracer reservoirs at OBCs [nondim] + real, dimension(:,:,:,:), pointer :: tr_field ! A pointer to a generic tracer field, in concentration units [conc] + real, dimension(:,:,:), pointer :: tr_ptr ! A pointer to a generic tracer field, in concentration units [conc] + real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)) :: grid_tmask ! A 3-d copy of G%mask2dT [nondim] + integer, dimension(SZI_(HI),SZJ_(HI)) :: grid_kmt ! A 2-d array of nk + + register_MOM_generic_tracer = .false. + if (associated(CS)) then + call MOM_error(FATAL, "register_MOM_generic_tracer called with an "// & + "associated control structure.") + endif + allocate(CS) + + + !Register all the generic tracers used and create the list of them. + !This can be called by ALL PE's. No array fields allocated. + if (.not. g_registered) then + call generic_tracer_register() + g_registered = .true. + endif + + +! Read all relevant parameters and write them to the model log. + call log_version(param_file, sub_name, version, "") + call get_param(param_file, sub_name, "GENERIC_TRACER_IC_FILE", CS%IC_file, & + "The file in which the generic tracer initial values can "//& + "be found, or an empty string for internal initialization.", & + default=" ") + if ((len_trim(CS%IC_file) > 0) .and. (scan(CS%IC_file,'/') == 0)) then + ! Add the directory if CS%IC_file is not already a complete path. + call get_param(param_file, sub_name, "INPUTDIR", inputdir, default=".") + CS%IC_file = trim(slasher(inputdir))//trim(CS%IC_file) + call log_param(param_file, sub_name, "INPUTDIR/GENERIC_TRACER_IC_FILE", CS%IC_file) + endif + call get_param(param_file, sub_name, "GENERIC_TRACER_IC_FILE_IS_Z", CS%Z_IC_file, & + "If true, GENERIC_TRACER_IC_FILE is in depth space, not "//& + "layer space.",default=.false.) + call get_param(param_file, sub_name, "TRACERS_MAY_REINIT", CS%tracers_may_reinit, & + "If true, tracers may go through the initialization code "//& + "if they are not found in the restart files. Otherwise "//& + "it is a fatal error if tracers are not found in the "//& + "restart files of a restarted run.", default=.false.) + + CS%restart_CSp => restart_CS + + ntau=1 ! MOM needs the fields at only one time step + + + ! At this point G%mask2dT and CS%diag%axesTL are not allocated. + ! postpone diag_registeration to initialize_MOM_generic_tracer + + !Fields cannot be diag registered as they are allocated and have to registered later. + grid_tmask(:,:,:) = 0.0 + grid_kmt(:,:) = 0 + axes(:) = -1 + + ! + ! Initialize all generic tracers + ! + call generic_tracer_init(HI%isc,HI%iec,HI%jsc,HI%jec,HI%isd,HI%ied,HI%jsd,HI%jed,& + GV%ke,ntau,axes,grid_tmask,grid_kmt,set_time(0,0)) + + + ! + ! MOM-register the generic tracers + ! + + !Get the tracer list + call generic_tracer_get_list(CS%g_tracer_list) + if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//& + ": No tracer in the list.") + ! For each tracer name get its T_prog index and get its fields + + g_tracer=>CS%g_tracer_list + do + call g_tracer_get_alias(g_tracer,g_tracer_name) + + call g_tracer_get_pointer(g_tracer,g_tracer_name,'field',tr_field) + call g_tracer_get_values(g_tracer,g_tracer_name,'longname', longname) + call g_tracer_get_values(g_tracer,g_tracer_name,'units',units ) + + !!nnz: MOM field is 3D. Does this affect performance? Need it be override field? + tr_ptr => tr_field(:,:,:,1) + ! Register prognostic tracer for horizontal advection, diffusion, and restarts. + if (g_tracer_is_prog(g_tracer)) then + call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, & + name=g_tracer_name, longname=longname, units=units, & + registry_diags=.true., & !### CHANGE TO TRUE? + restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit) + else + call register_restart_field(tr_ptr, g_tracer_name, .not.CS%tracers_may_reinit, & + restart_CS, longname=longname, units=units) + endif + + !traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer=>g_tracer_next + + enddo + + register_MOM_generic_tracer = .true. +end function register_MOM_generic_tracer + +!> Register OBC segments for generic tracers +subroutine register_MOM_generic_tracer_segments(CS, GV, OBC, tr_Reg, param_file) + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies whether, + !! where, and what open boundary conditions are used. + type(tracer_registry_type), pointer :: tr_Reg !< Pointer to the control structure for the tracer + !! advection and diffusion module. + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + + ! Local variables + logical :: obc_has + ! This include declares and sets the variable "version". +# include "version_variable.h" + character(len=128), parameter :: sub_name = 'register_MOM_generic_tracer_segments' + type(g_tracer_type), pointer :: g_tracer,g_tracer_next + character(len=fm_string_len) :: g_tracer_name + character(len=fm_string_len) :: obc_src_file_name, obc_src_field_name + real :: lfac_in ! Multiplicative factor used in setting the tracer-specific inverse length + ! scales associated with inflowing tracer reservoirs at OBCs [nondim] + real :: lfac_out ! Multiplicative factor used in setting the tracer-specific inverse length + ! scales associated with outflowing tracer reservoirs at OBCs [nondim] + + if (.NOT. associated(OBC)) return + !Get the tracer list + call generic_tracer_get_list(CS%g_tracer_list) + if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//& + ": No tracer in the list.") + + g_tracer=>CS%g_tracer_list + do + call g_tracer_get_alias(g_tracer,g_tracer_name) + if (g_tracer_is_prog(g_tracer)) then + call g_tracer_get_obc_segment_props(g_tracer,g_tracer_name,obc_has ,& + obc_src_file_name,obc_src_field_name,lfac_in,lfac_out) + if (obc_has) then + call set_obgc_segments_props(OBC,g_tracer_name,obc_src_file_name,obc_src_field_name,lfac_in,lfac_out) + call register_obgc_segments(GV, OBC, tr_Reg, param_file, g_tracer_name) + endif + endif + + !traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer=>g_tracer_next + + enddo + +end subroutine register_MOM_generic_tracer_segments + +!> Initialize phase II: Initialize required variables for generic tracers +!! There are some steps of initialization that cannot be done in register_MOM_generic_tracer +!! This is the place and time to do them: +!! Set the grid mask and initial time for all generic tracers. +!! Diag_register them. +!! Z_diag_register them. +!! +!! This subroutine initializes the NTR tracer fields in tr(:,:,:,:) +!! and it sets up the tracer output. +subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, tv, param_file, diag, OBC, & + CS, sponge_CSp, ALE_sponge_CSp) + logical, intent(in) :: restart !< .true. if the fields have already been + !! read from a restart file. + type(time_type), target, intent(in) :: day !< Time of the start of the run. + type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic + !! variables + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + type(diag_ctrl), target, intent(in) :: diag !< Regulates diagnostic output. + type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type specifies whether, + !! where, and what open boundary conditions are used. + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + type(sponge_CS), pointer :: sponge_CSp !< Pointer to the control structure for the sponges. + type(ALE_sponge_CS), pointer :: ALE_sponge_CSp !< Pointer to the control structure for the + !! ALE sponges. + + character(len=128), parameter :: sub_name = 'initialize_MOM_generic_tracer' + logical :: OK,obc_has + logical :: do_use_GT_sponge + integer :: i, j, k, isc, iec, jsc, jec, nk + type(g_tracer_type), pointer :: g_tracer,g_tracer_next + character(len=fm_string_len) :: g_tracer_name + real, dimension(:,:,:,:), pointer :: tr_field ! A pointer to a generic tracer field, in concentration units [conc] + real, dimension(:,:,:), pointer :: tr_ptr ! A pointer to a generic tracer field, in concentration units [conc] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Layer vertical extent [Z ~> m] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: grid_tmask ! A 3-d copy of G%mask2dT [nondim] + integer, dimension(SZI_(G),SZJ_(G)) :: grid_kmt ! A 2-d array of nk + + !! 2010/02/04 Add code to re-initialize Generic Tracers if needed during a model simulation + !! By default, restart cpio should not contain a Generic Tracer IC file and step below will be skipped. + !! Ideally, the generic tracer IC file should have the tracers on Z levels. + + isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke + + CS%diag=>diag + !Get the tracer list + if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//& + ": No tracer in the list.") + !For each tracer name get its fields + g_tracer=>CS%g_tracer_list + + call thickness_to_dz(h, tv, dz, G, GV, US) + + do + if (INDEX(CS%IC_file, '_NULL_') /= 0) then + call MOM_error(WARNING, "The name of the IC_file "//trim(CS%IC_file)//& + " indicates no MOM initialization was asked for the generic tracers."//& + "Bypassing the MOM initialization of ALL generic tracers!") + exit + endif + call g_tracer_get_alias(g_tracer,g_tracer_name) + call g_tracer_get_pointer(g_tracer,g_tracer_name,'field',tr_field) + tr_ptr => tr_field(:,:,:,1) + + if (.not.restart .or. (CS%tracers_may_reinit .and. & + .not.query_initialized(tr_ptr, g_tracer_name, CS%restart_CSp))) then + + if (g_tracer%requires_src_info ) then + call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//& + "initializing generic tracer "//trim(g_tracer_name)//& + " using MOM_initialize_tracer_from_Z ") + + call MOM_initialize_tracer_from_Z(dz, tr_ptr, G, GV, US, param_file, & + src_file=g_tracer%src_file, src_var_nam=g_tracer%src_var_name, & + src_var_unit_conversion=g_tracer%src_var_unit_conversion, & + src_var_record=g_tracer%src_var_record, src_var_gridspec=g_tracer%src_var_gridspec, & + h_in_Z_units=.true.) + + !Check/apply the bounds for each g_tracer + do k=1,nk ; do j=jsc,jec ; do i=isc,iec + if (tr_ptr(i,j,k) /= CS%tracer_land_val) then + if (tr_ptr(i,j,k) < g_tracer%src_var_valid_min) tr_ptr(i,j,k) = g_tracer%src_var_valid_min + !Jasmin does not want to apply the maximum for now + !if (tr_ptr(i,j,k) > g_tracer%src_var_valid_max) tr_ptr(i,j,k) = g_tracer%src_var_valid_max + endif + enddo ; enddo ; enddo + + !jgj: Reset CASED to 0 below K=1 + if ( (trim(g_tracer_name) == 'cased') .or. (trim(g_tracer_name) == 'ca13csed') ) then + do k=2,nk ; do j=jsc,jec ; do i=isc,iec + if (tr_ptr(i,j,k) /= CS%tracer_land_val) then + tr_ptr(i,j,k) = 0.0 + endif + enddo ; enddo ; enddo + endif + elseif(.not. g_tracer%requires_restart) then + !Do nothing for this tracer, it is initialized by the tracer package + call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//& + "skip initialization of generic tracer "//trim(g_tracer_name)) + else !Do it old way if the tracer is not registered to start from a specific source file. + !This path should be deprecated if all generic tracers are required to start from specified sources. + if (len_trim(CS%IC_file) > 0) then + ! Read the tracer concentrations from a netcdf file. + if (.not.file_exists(CS%IC_file)) call MOM_error(FATAL, & + "initialize_MOM_Generic_tracer: Unable to open "//CS%IC_file) + if (CS%Z_IC_file) then + OK = tracer_Z_init(tr_ptr, h, CS%IC_file, g_tracer_name, G, GV, US) + if (.not.OK) then + OK = tracer_Z_init(tr_ptr, h, CS%IC_file, trim(g_tracer_name), G, GV, US) + if (.not.OK) call MOM_error(FATAL,"initialize_MOM_Generic_tracer: "//& + "Unable to read "//trim(g_tracer_name)//" from "//& + trim(CS%IC_file)//".") + endif + call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//& + "initialized generic tracer "//trim(g_tracer_name)//& + " using Generic Tracer File on Z: "//CS%IC_file) + else + ! native grid + call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//& + "Using Generic Tracer IC file on native grid "//trim(CS%IC_file)//& + " for tracer "//trim(g_tracer_name)) + call MOM_read_data(CS%IC_file, trim(g_tracer_name), tr_ptr, G%Domain) + endif + else + call MOM_error(FATAL,"initialize_MOM_generic_tracer: "//& + "check Generic Tracer IC filename "//trim(CS%IC_file)//& + " for tracer "//trim(g_tracer_name)) + endif + + endif + + call set_initialized(tr_ptr, g_tracer_name, CS%restart_CSp) + endif + + call g_tracer_get_obc_segment_props(g_tracer,g_tracer_name,obc_has ) + if(obc_has .and. g_tracer_is_prog(g_tracer) .and. .not.restart) & + call fill_obgc_segments(G, GV, OBC, tr_ptr, g_tracer_name) + !traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer=>g_tracer_next + enddo + !! end section to re-initialize generic tracers + + call get_param(param_file, "MOM", "PHA_MLD_CALC", CS%mld_pha_calc, & + "If false, use a fixed value for the photoacclimation mixed layer depth within the "//& + "generic tracer update. This MLD is only used for photoacclimation. This variable should "//& + "be set to true if using COBALTv3 for the BGC.", default=.false.) + if (CS%mld_pha_calc) then + call get_param(param_file, "MOM", "PHA_MLD_USE_DELTA_RHO", CS%mld_pha_use_delta_rho, & + "If true, use a density difference to find the photoacclimation mixed layer depth "//& + "within the generic tracer update. This MLD is only used for photoacclimation.", default=.false.) + call get_param(param_file, "MOM", "PHA_MLD_USE_DELTA_ENG", CS%mld_pha_use_delta_eng, & + "If true, use an energy difference to find the photoacclimation mixed layer depth "//& + "with the generic tracer update. This MLD is only used for photoacclimation.", default=.false.) + if (CS%mld_pha_use_delta_rho .and. CS%mld_pha_use_delta_eng) then + call MOM_error(FATAL, "PHA_MLD_CALC is set to true and PHA_MLD_USE_DELTA_RHO and PHA_MLD_USE_DELTA_ENG "// & + "are both true. Choose only one option for the calculated photoacclimation MLD!") + elseif ((.not.CS%mld_pha_use_delta_rho) .and. (.not.CS%mld_pha_use_delta_eng)) then + call MOM_error(FATAL, "PHA_MLD_CALC is set to true but PHA_MLD_USE_DELTA_RHO and PHA_MLD_USE_DELTA_ENG "// & + "are both false. Choose an option for the calculated photoacclimation MLD!") + endif + if (CS%mld_pha_use_delta_rho) then + call get_param(param_file, "MOM", "PHA_MLD_HREF", CS%mld_pha_href, & + "The reference depth for a density difference based photoacclimation MLD [m].", & + units='m', default=0.0, scale=US%m_to_Z, do_not_log=.not.CS%mld_pha_use_delta_rho) + call get_param(param_file, "MOM", "PHA_MLD_DRHO", CS%mld_pha_drho, & + "The density difference for a density difference based photoacclimation MLD [kg m-3].", & + units='kg/m3', default=0.03, scale=US%kg_m3_to_R, do_not_log=.not.CS%mld_pha_use_delta_rho) + elseif (CS%mld_pha_use_delta_eng) then + call get_param(param_file, "MOM", "PHA_MLD_DENG", CS%mld_pha_deng, & + "The energy for an energy difference based photoacclimation MLD.", default=25.0, & + units='J/m2',scale=US%W_m2_to_RZ3_T3*US%s_to_T, do_not_log=.not.CS%mld_pha_use_delta_eng) + endif + else + call get_param(param_file, "MOM", "PHA_MLD_VAL", CS%mld_pha_val, & + "The depth of photoacclimation if fixed depth is used [m].", & + units='m', default=0.0, scale=US%m_to_Z) + endif + + + !Now we can reset the grid mask, axes and time to their true values + !Note that grid_tmask must be set correctly on the data domain boundary + !so that coast mask can be deduced from it. + grid_tmask(:,:,:) = 0.0 + grid_kmt(:,:) = 0 + do j = G%jsd, G%jed ; do i = G%isd, G%ied + if (G%mask2dT(i,j) > 0.0) then + grid_tmask(i,j,:) = 1.0 + grid_kmt(i,j) = GV%ke ! Tell the code that a layer thicker than 1m is the bottom layer. + endif + enddo ; enddo + call g_tracer_set_common(G%isc,G%iec,G%jsc,G%jec,G%isd,G%ied,G%jsd,G%jed,& + GV%ke,1,CS%diag%axesTL%handles,grid_tmask,grid_kmt,day) + + call get_param(param_file, "initialize_sponges_file", "DO_SPONGE_GENERIC_TRACER", do_use_gt_sponge, & + "If true, then some generic tracers may be nudged.", default=.false.) + if (do_use_GT_sponge) call g_tracer_initialize_sponges(G, GV, US, CS, param_file, sponge_CSp, ALE_sponge_CSp, day) + + ! Register generic tracer modules diagnostics + +#ifdef _USE_MOM6_DIAG + call g_tracer_set_csdiag(CS%diag) +#endif + call generic_tracer_register_diag() +#ifdef _USE_MOM6_DIAG + call g_tracer_set_csdiag(CS%diag) +#endif + +end subroutine initialize_MOM_generic_tracer + +subroutine g_tracer_initialize_sponges(G, GV, US, CS, param_file, Layer_CSp, ALE_CSp, Time) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G)) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters. + type(sponge_CS), pointer :: Layer_CSp !< A pointer that is set to point to the control + !! structure for this module (in layered mode). + type(ALE_sponge_CS), pointer :: ALE_CSp !< A pointer that is set to point to the control + !! structure for this module (in ALE mode). + type(time_type), intent(in) :: Time !< Time at the start of the run segment. Time_in + !! overrides any value set for Time. + ! Local variables + real, allocatable, dimension(:,:,:) :: eta ! The target interface heights [Z ~> m]. + real, allocatable, dimension(:,:,:) :: dz ! The target interface thicknesses in height units [Z ~> m] + real, allocatable, dimension(:,:,:) :: h ! The target interface thicknesses [H ~> m or kg m-2]. + + real, allocatable, dimension(:,:,:) :: tmp_GT ! A temporary array for reading sponge target temperatures + ! on the vertical grid of the input file [C ~> degC] + + real :: Idamp(SZI_(G),SZJ_(G)) ! The sponge damping rate [T-1 ~> s-1] + + integer :: i, j, k, is, ie, js, je, nz + integer :: isd, ied, jsd, jed + integer, dimension(4) :: siz + integer :: nz_data ! The size of the sponge source grid + character(len=40) :: tmp_var, Idamp_var, eta_var + character(len=40) :: mdl = "initialize_sponges_file" + character(len=200) :: damping_file, state_file, state_uv_file ! Strings for filenames + character(len=200) :: filename, inputdir ! Strings for file/path and path. + + character(len=200) :: g_tracer_name + type(g_tracer_type), pointer :: g_tracer,g_tracer_next + logical :: do_sponge_gt + real, dimension(:,:,:), pointer :: g_tracer_ptr + + logical :: use_ALE ! True if ALE is being used, False if in layered mode + logical :: time_space_interp_sponge ! If true use sponge data that need to be interpolated in both + ! the horizontal dimension and in time prior to vertical remapping. + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed + + Idamp(:,:) = 0.0 + + call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") + inputdir = slasher(inputdir) + call get_param(param_file, mdl, "SPONGE_GT_DAMPING_FILE", damping_file, & + "The name of the file with the sponge damping rates.") !, & + call get_param(param_file, mdl, "SPONGE_IDAMP_VAR", Idamp_var, & + "The name of the inverse damping rate variable in "//& + "SPONGE_DAMPING_FILE.", default="Idamp") + call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, default=.false., do_not_log=.true.) + call get_param(param_file, mdl, "INTERPOLATE_SPONGE_TIME_SPACE", time_space_interp_sponge, & + "If True, perform on-the-fly regridding in lat-lon-time of sponge restoring data.", & + default=.false.) + + ! Read in sponge damping rate for tracers + filename = trim(inputdir)//trim(damping_file) + call log_param(param_file, mdl, "INPUTDIR/SPONGE_DAMPING_FILE", filename) + if (.not.file_exists(filename, G%Domain)) & + call MOM_error(FATAL, " initialize_sponges: Unable to open "//trim(filename)) + + call MOM_read_data(filename, Idamp_var, Idamp(:,:), G%Domain, scale=US%T_to_s) + + ! Now register all of the fields which are damped in the sponge. + ! By default, momentum is advected vertically within the sponge, but + ! momentum is typically not damped within the sponge. + + g_tracer=>CS%g_tracer_list + do ! check for each generic tracer if it is nudged + + call g_tracer_get_alias(g_tracer,g_tracer_name) + + call get_param(param_file, mdl, "DO_SPONGE_GT_"//trim(g_tracer_name), do_sponge_gt, & + "If true, then generic tracer "//trim(g_tracer_name)//" is nudged.", & + default=.false.) + + if (do_sponge_gt) then + call get_param(param_file, mdl, "SPONGE_GT_"//trim(g_tracer_name)//"_FILE", state_file, & + "The name of the file with the state to damp the generic ", & + "tracer "//trim(g_tracer_name)//" toward.", default="sponge_"//trim(g_tracer_name)//".nc.") + call get_param(param_file, mdl, "SPONGE_GT_"//trim(g_tracer_name)//"_VAR", tmp_var, & + "The name of the variable to use in the sponge_GT file for generic ", & + "tracer "//trim(g_tracer_name)//".", default=trim(g_tracer_name)) + call get_param(param_file, mdl, "SPONGE_GT_"//trim(g_tracer_name)//"_ETA_VAR", eta_var, & + "The name of the interface height variable in "//& + "SPONGE_GT_"//trim(g_tracer_name)//"_FILE.", default="ETA") + + filename = trim(inputdir)//trim(state_file) + call log_param(param_file, mdl, "INPUTDIR/SPONGE_GT"//trim(g_tracer_name)//"_FILE", filename) + if (.not.file_exists(filename, G%Domain)) & + call MOM_error(FATAL, " initialize_sponges: Unable to open "//trim(filename)) + + ! get the pointer for this tracer + call g_tracer_get_pointer(g_tracer,g_tracer_name,'field',g_tracer_ptr) + + if (use_ALE) then ! ALE mode + if (.not. time_space_interp_sponge) then + !call field_size(filename,eta_var,siz,no_domain=.true.) + if (siz(1) /= G%ieg-G%isg+1 .or. siz(2) /= G%jeg-G%jsg+1) & + call MOM_error(FATAL,"initialize_sponge_file: Array size mismatch for sponge data.") + nz_data = siz(3)-1 + allocate(eta(isd:ied,jsd:jed,nz_data+1)) + allocate(dz(isd:ied,jsd:jed,nz_data)) + call MOM_read_data(filename, eta_var, eta(:,:,:), G%Domain, scale=US%m_to_Z) + do j=js,je ; do i=is,ie + eta(i,j,nz_data+1) = -depth_tot(i,j) + enddo ; enddo + do k=nz_data,1,-1 ; do j=js,je ; do i=is,ie + if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_Z)) & + eta(i,j,K) = eta(i,j,K+1) + GV%Angstrom_Z + enddo ; enddo ; enddo + do k=1,nz_data ; do j=js,je ; do i=is,ie + dz(i,j,k) = eta(i,j,k)-eta(i,j,k+1) + enddo; enddo ; enddo + deallocate(eta) + + allocate(tmp_GT(isd:ied,jsd:jed,nz_data)) + call MOM_read_data(filename, tmp_var, tmp_GT(:,:,:), G%Domain) !, scale=US%degC_to_C) + + call set_up_ALE_sponge_field(tmp_GT, G, GV, g_tracer_ptr, & + ALE_CSp, trim(g_tracer_name)) + deallocate(tmp_GT) + deallocate(dz) + + else + call set_up_ALE_sponge_field(filename, tmp_var, Time, G, GV, US, g_tracer_ptr, & + ALE_CSp, trim(g_tracer_name)) + endif + endif + endif + + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer=>g_tracer_next + enddo + +end subroutine g_tracer_initialize_sponges + +!> Column physics for generic tracers. +!! Get the coupler values for generic tracers that exchange with atmosphere +!! Update generic tracer concentration fields from sources and sinks. +!! Vertically diffuse generic tracer concentration fields. +!! Update generic tracers from bottom and their bottom reservoir. +!! +!! This subroutine applies diapycnal diffusion and any other column +!! tracer physics or chemistry to the tracers from this file. +!! CFCs are relatively simple, as they are passive tracers. with only a surface +!! flux as a source. +subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, US, CS, tv, optics, & + evap_CFL_limit, minimum_forcing_depth) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: ea !< The amount of fluid entrained from the layer + !! above during this call [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: eb !< The amount of fluid entrained from the layer + !! below during this call [H ~> m or kg m-2]. + type(forcing), intent(in) :: fluxes !< A structure containing pointers to thermodynamic + !! and tracer forcing fields. + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: Hml !< Mixed layer depth [Z ~> m] + real, intent(in) :: dt !< The amount of time covered by this call [T ~> s] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables + type(optics_type), intent(in) :: optics !< The structure containing optical properties. + real, optional, intent(in) :: evap_CFL_limit !< Limit on the fraction of the water that can + !! be fluxed out of the top layer in a timestep [nondim] + ! Stored previously in diabatic CS. + real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which fluxes + !! can be applied [H ~> m or kg m-2] + ! Stored previously in diabatic CS. + ! The arguments to this subroutine are redundant in that + ! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1) + + ! Local variables + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_column_physics' + + type(g_tracer_type), pointer :: g_tracer, g_tracer_next + character(len=fm_string_len) :: g_tracer_name + real, dimension(:,:), pointer :: stf_array ! The surface flux of the tracer [conc kg m-2 s-1] + real, dimension(:,:), pointer :: trunoff_array ! The tracer concentration in the river runoff [conc] + real, dimension(:,:), pointer :: runoff_tracer_flux_array ! The runoff tracer flux [conc kg m-2 s-1] + + real :: surface_field(SZI_(G),SZJ_(G)) ! The surface value of some field, here only used for salinity [S ~> ppt] + real :: dz_ml(SZI_(G),SZJ_(G)) ! The mixed layer depth in the MKS units used for generic tracers [m] + real :: sosga ! The global mean surface salinity [ppt] + + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: rho_dzt ! Layer mass per unit area [kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dzt ! Layer vertical extents [m] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! A work array of thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G)) :: mld_pha ! The mixed layer depth calculated for photoacclimation + ! that is used in COBALTv3 + integer :: i, j, k, isc, iec, jsc, jec, nk + + isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke + + !Get the tracer list + if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL,& + trim(sub_name)//": No tracer in the list.") + +#ifdef _USE_MOM6_DIAG + call g_tracer_set_csdiag(CS%diag) +#endif + + ! + !Extract the tracer surface fields from coupler and update tracer fields from sources + ! + !call generic_tracer_coupler_get(fluxes%tr_fluxes) + !Niki: This is moved out to ocean_model_MOM.F90 because if dt_therm>dt_cpld we need to average + ! the fluxes without coming into this subroutine. + ! MOM5 has to modified to conform. + + ! + !Add contribution of river to surface flux + ! + g_tracer=>CS%g_tracer_list + do + if (_ALLOCATED(g_tracer%trunoff) .and. (.NOT. g_tracer%runoff_added_to_stf)) then + call g_tracer_get_alias(g_tracer,g_tracer_name) + call g_tracer_get_pointer(g_tracer,g_tracer_name,'stf', stf_array) + call g_tracer_get_pointer(g_tracer,g_tracer_name,'trunoff',trunoff_array) + call g_tracer_get_pointer(g_tracer,g_tracer_name,'runoff_tracer_flux',runoff_tracer_flux_array) + !nnz: Why is fluxes%river = 0? + runoff_tracer_flux_array(:,:) = trunoff_array(:,:) * & + US%RZ_T_to_kg_m2s*fluxes%lrunoff(:,:) + stf_array = stf_array + runoff_tracer_flux_array + g_tracer%runoff_added_to_stf = .true. + endif + + !traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer => g_tracer_next + + enddo + + ! + !Prepare input arrays for source update + ! + + rho_dzt(:,:,:) = GV%H_to_kg_m2 * GV%Angstrom_H + do k=1,nk ; do j=jsc,jec ; do i=isc,iec + rho_dzt(i,j,k) = GV%H_to_kg_m2 * h_old(i,j,k) + enddo ; enddo ; enddo + + dzt(:,:,:) = 1.0 + call thickness_to_dz(h_old, tv, dzt, G, GV, US) + do k=1,nk ; do j=jsc,jec ; do i=isc,iec + dzt(i,j,k) = US%Z_to_m * dzt(i,j,k) + enddo ; enddo ; enddo + dz_ml(:,:) = 0.0 + do j=jsc,jec ; do i=isc,iec + surface_field(i,j) = tv%S(i,j,1) + dz_ml(i,j) = US%Z_to_m * Hml(i,j) + enddo ; enddo + sosga = global_area_mean(surface_field, G, unscale=US%S_to_ppt) + + mld_pha(:,:) = 0.0 + if (.not.CS%mld_pha_calc) then + mld_pha(:,:) = CS%mld_pha_val + else + if (CS%mld_pha_use_delta_rho) then + call diagnoseMLDbyDensityDifference(-1, h_old, tv, CS%mld_pha_drho, G, GV, US, CS%diag, & + CS%mld_pha_href, id_ref_z=-1, id_ref_rho=-1, MLD_out=mld_pha) + elseif (CS%mld_pha_use_delta_eng) then + call diagnoseMLDbyEnergy((/-1, -1, -1/), h_old, tv, G, GV, US, (/CS%mld_pha_deng, & + CS%mld_pha_deng, CS%mld_pha_deng/), CS%diag, MLD_out=mld_pha) + endif + endif + + ! + !Calculate tendencies (i.e., field changes at dt) from the sources / sinks + ! + if ((G%US%L_to_m == 1.0) .and. (G%US%s_to_T == 1.0) .and. (G%US%Z_to_m == 1.0) .and. & + (G%US%Q_to_J_kg == 1.0) .and. (G%US%RZ_to_kg_m2 == 1.0) .and. & + (US%C_to_degC == 1.0) .and. (US%S_to_ppt == 1.0)) then + ! Avoid unnecessary copies when no unit conversion is needed. + call generic_tracer_source(tv%T, tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, & + G%areaT, get_diag_time_end(CS%diag), & + optics%nbands, optics%max_wavelength_band, optics%sw_pen_band, optics%opacity_band, & + internal_heat=tv%internal_heat, frunoff=fluxes%frunoff, sosga=sosga, geolat=G%geolatT, & + photo_acc_dpth=mld_pha) + else + ! tv%internal_heat is a null pointer unless DO_GEOTHERMAL = True, + ! so we have to check and only do the scaling if it is associated. + if(associated(tv%internal_heat)) then + call generic_tracer_source(US%C_to_degC*tv%T, US%S_to_ppt*tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, & + G%US%L_to_m**2*G%areaT(:,:), get_diag_time_end(CS%diag), & + optics%nbands, optics%max_wavelength_band, & + sw_pen_band=G%US%QRZ_T_to_W_m2*optics%sw_pen_band(:,:,:), & + opacity_band=G%US%m_to_Z*optics%opacity_band(:,:,:,:), & + internal_heat=G%US%RZ_to_kg_m2*US%C_to_degC*tv%internal_heat(:,:), & + frunoff=G%US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga, geolat=G%geolatT, & + photo_acc_dpth=mld_pha*US%Z_to_m) + else + call generic_tracer_source(US%C_to_degC*tv%T, US%S_to_ppt*tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, & + G%US%L_to_m**2*G%areaT(:,:), get_diag_time_end(CS%diag), & + optics%nbands, optics%max_wavelength_band, & + sw_pen_band=G%US%QRZ_T_to_W_m2*optics%sw_pen_band(:,:,:), & + opacity_band=G%US%m_to_Z*optics%opacity_band(:,:,:,:), & + frunoff=G%US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga, geolat=G%geolatT, & + photo_acc_dpth=mld_pha*US%Z_to_m) + endif + endif + + ! This uses applyTracerBoundaryFluxesInOut to handle the change in tracer due to freshwater fluxes + ! usually in ALE mode + if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then + g_tracer=>CS%g_tracer_list + do + if (g_tracer_is_prog(g_tracer)) then + do k=1,nk ;do j=jsc,jec ; do i=isc,iec + h_work(i,j,k) = h_old(i,j,k) + if (g_tracer%diag_id_boundary_forcing_tend .gt. 0) then + g_tracer%boundary_forcing_tend(i,j,k) = g_tracer%field(i,j,k,1) + endif + enddo ; enddo ; enddo + call applyTracerBoundaryFluxesInOut(G, GV, g_tracer%field(:,:,:,1), dt, & + fluxes, h_work, evap_CFL_limit, minimum_forcing_depth) + if (g_tracer%diag_id_boundary_forcing_tend .gt. 0) then + do k=1,nk ;do j=jsc,jec ; do i=isc,iec + g_tracer%boundary_forcing_tend(i,j,k)=G%mask2dT(i,j)*(g_tracer%field(i,j,k,1) & + - g_tracer%boundary_forcing_tend(i,j,k))/dt + enddo ; enddo ; enddo + endif + endif + + !traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer=>g_tracer_next + enddo + endif + + ! + !Update Tr(n)%field from explicit vertical diffusion + ! + ! Use a tridiagonal solver to determine the concentrations after the + ! surface source is applied and diapycnal advection and diffusion occurs. + if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then + ! Last arg is tau which is always 1 for MOM6 + call generic_tracer_vertdiff_G(h_work, ea, eb, US%T_to_s*dt, GV%kg_m2_to_H, GV%m_to_H, 1) + else + ! Last arg is tau which is always 1 for MOM6 + call generic_tracer_vertdiff_G(h_old, ea, eb, US%T_to_s*dt, GV%kg_m2_to_H, GV%m_to_H, 1) + endif + + ! Update bottom fields after vertical processes + + ! Second arg is tau which is always 1 for MOM6 + call generic_tracer_update_from_bottom(US%T_to_s*dt, 1, get_diag_time_end(CS%diag)) + + !Output diagnostics via diag_manager for all generic tracers and their fluxes + call g_tracer_send_diag(CS%g_tracer_list, get_diag_time_end(CS%diag), tau=1) +#ifdef _USE_MOM6_DIAG + call g_tracer_set_csdiag(CS%diag) +#endif + +end subroutine MOM_generic_tracer_column_physics + +!> This subroutine calculates mass-weighted integral on the PE either +!! of all available tracer concentrations, or of a tracer that is +!! being requested specifically, returning the number of stocks it has +!! calculated. If the stock_index is present, only the stock corresponding +!! to that coded index is returned. +function MOM_generic_tracer_stock(h, stocks, G, GV, CS, names, units, stock_index) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + type(EFP_type), dimension(:), intent(out) :: stocks !< The mass-weighted integrated amount of each + !! tracer, in kg times concentration units [kg conc] + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + character(len=*), dimension(:), intent(out) :: names !< The names of the stocks calculated. + character(len=*), dimension(:), intent(out) :: units !< The units of the stocks calculated. + integer, optional, intent(in) :: stock_index !< The coded index of a specific stock + !! being sought. + integer :: MOM_generic_tracer_stock !< Return value, the + !! number of stocks calculated here. + + ! Local variables + type(g_tracer_type), pointer :: g_tracer, g_tracer_next + real, dimension(:,:,:,:), pointer :: tr_field ! A pointer to a generic tracer field, in concentration units [conc] + real, dimension(:,:,:), pointer :: tr_ptr ! A pointer to a generic tracer field, in concentration units [conc] + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_stock' + + integer :: m + + MOM_generic_tracer_stock = 0 + if (.not.associated(CS)) return + + if (present(stock_index)) then ; if (stock_index > 0) then + ! Check whether this stock is available from this routine. + + ! No stocks from this routine are being checked yet. Return 0. + return + endif ; endif + + if (.NOT. associated(CS%g_tracer_list)) return ! No stocks. + + m=1 ; g_tracer=>CS%g_tracer_list + do + call g_tracer_get_alias(g_tracer,names(m)) + call g_tracer_get_values(g_tracer,names(m),'units',units(m)) + units(m) = trim(units(m))//" kg" + call g_tracer_get_pointer(g_tracer,names(m),'field',tr_field) + + tr_ptr => tr_field(:,:,:,1) + stocks(m) = global_mass_int_EFP(h, G, GV, tr_ptr, on_PE_only=.true.) + + !traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer=>g_tracer_next + m = m+1 + enddo + + MOM_generic_tracer_stock = m + +end function MOM_generic_tracer_stock + +!> This subroutine finds the global min and max of either of all available +!! tracer concentrations, or of a tracer that is being requested specifically, +!! returning the number of tracers it has evaluated. +!! It also optionally returns the locations of the extrema. +function MOM_generic_tracer_min_max(ind_start, got_minmax, gmin, gmax, G, CS, names, units, & + xgmin, ygmin, zgmin, xgmax, ygmax, zgmax) + integer, intent(in) :: ind_start !< The index of the tracer to start with + logical, dimension(:), intent(out) :: got_minmax !< Indicates whether the global min and + !! max are found for each tracer + real, dimension(:), intent(out) :: gmin !< Global minimum of each tracer [conc] + real, dimension(:), intent(out) :: gmax !< Global maximum of each tracer [conc] + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + character(len=*), dimension(:), intent(out) :: names !< The names of the stocks calculated. + character(len=*), dimension(:), intent(out) :: units !< The units of the stocks calculated. + real, dimension(:), optional, intent(out) :: xgmin !< The x-position of the global minimum in the + !! units of G%geoLonT, often [degrees_E] or [km] or [m] + real, dimension(:), optional, intent(out) :: ygmin !< The y-position of the global minimum in the + !! units of G%geoLatT, often [degrees_N] or [km] or [m] + real, dimension(:), optional, intent(out) :: zgmin !< The z-position of the global minimum [layer] + real, dimension(:), optional, intent(out) :: xgmax !< The x-position of the global maximum in the + !! units of G%geoLonT, often [degrees_E] or [km] or [m] + real, dimension(:), optional, intent(out) :: ygmax !< The y-position of the global maximum in the + !! units of G%geoLatT, often [degrees_N] or [km] or [m] + real, dimension(:), optional, intent(out) :: zgmax !< The z-position of the global maximum [layer] + integer :: MOM_generic_tracer_min_max !< Return value, the + !! number of tracers done here. + + ! Local variables + type(g_tracer_type), pointer :: g_tracer, g_tracer_next + real, dimension(:,:,:,:), pointer :: tr_field ! The tracer array whose extrema are being sought [conc] + real, dimension(:,:,:), pointer :: tr_ptr ! The tracer array whose extrema are being sought [conc] + real :: x_min ! The x-position of the global minimum in the units of G%geoLonT, often [degrees_E] or [km] or [m] + real :: y_min ! The y-position of the global minimum in the units of G%geoLatT, often [degrees_N] or [km] or [m] + real :: z_min ! The z-position of the global minimum [layer] + real :: x_max ! The x-position of the global maximum in the units of G%geoLonT, often [degrees_E] or [km] or [m] + real :: y_max ! The y-position of the global maximum in the units of G%geoLatT, often [degrees_N] or [km] or [m] + real :: z_max ! The z-position of the global maximum [layer] + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_min_max' + + logical :: find_location + integer :: isc, iec, jsc, jec, isd, ied, jsd, jed, nk, ntau + integer :: k, is, ie, js, je, m + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + + MOM_generic_tracer_min_max = 0 + if (.not.associated(CS)) return + + if (.NOT. associated(CS%g_tracer_list)) return ! No stocks. + + call g_tracer_get_common(isc, iec, jsc, jec, isd, ied, jsd, jed, nk, ntau) + find_location = present(xgmin) .or. present(ygmin) .or. present(zgmin) .or. & + present(xgmax) .or. present(ygmax) .or. present(zgmax) + + m=ind_start ; g_tracer=>CS%g_tracer_list + do + call g_tracer_get_alias(g_tracer,names(m)) + call g_tracer_get_values(g_tracer,names(m),'units',units(m)) + call g_tracer_get_pointer(g_tracer,names(m),'field',tr_field) + + gmin(m) = -1.0 + gmax(m) = -1.0 + + tr_ptr => tr_field(:,:,:,1) + + if (find_location) then + call array_global_min_max(tr_ptr, G, nk, gmin(m), gmax(m), & + x_min, y_min, z_min, x_max, y_max, z_max) + if (present(xgmin)) xgmin(m) = x_min + if (present(ygmin)) ygmin(m) = y_min + if (present(zgmin)) zgmin(m) = z_min + if (present(xgmax)) xgmax(m) = x_max + if (present(ygmax)) ygmax(m) = y_max + if (present(zgmax)) zgmax(m) = z_max + else + call array_global_min_max(tr_ptr, G, nk, gmin(m), gmax(m)) + endif + + got_minmax(m) = .true. + + !traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer=>g_tracer_next + m = m+1 + enddo + + MOM_generic_tracer_min_max = m + +end function MOM_generic_tracer_min_max + +!> This subroutine calculates the surface state and sets coupler values for +!! those generic tracers that have flux exchange with atmosphere. +!! +!! This subroutine sets up the fields that the coupler needs to calculate the +!! CFC fluxes between the ocean and atmosphere. +subroutine MOM_generic_tracer_surface_state(sfc_state, h, G, GV, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(surface), intent(inout) :: sfc_state !< A structure containing fields that + !! describe the surface state of the ocean. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + + ! Local variables + real :: sosga ! The global mean surface salinity [ppt] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV),1) :: rho0 ! An unused array of densities [kg m-3] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dzt ! Layer vertical extents [m] + + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_surface_state' + + !Set coupler values + !nnz: fake rho0 + rho0(:,:,:,:) = 1.0 + + dzt(:,:,:) = GV%H_to_m * h(:,:,:) + + sosga = global_area_mean(sfc_state%SSS, G, unscale=G%US%S_to_ppt) + + if ((G%US%C_to_degC == 1.0) .and. (G%US%S_to_ppt == 1.0)) then + call generic_tracer_coupler_set(sfc_state%tr_fields, & + ST=sfc_state%SST, SS=sfc_state%SSS, & + rho=rho0, & !nnz: required for MOM5 and previous versions. + ilb=G%isd, jlb=G%jsd, & + dzt=dzt,& !This is needed for the Mocsy method of carbonate system vars + tau=1, sosga=sosga, model_time=get_diag_time_end(CS%diag)) + else + call generic_tracer_coupler_set(sfc_state%tr_fields, & + ST=G%US%C_to_degC*sfc_state%SST, SS=G%US%S_to_ppt*sfc_state%SSS, & + rho=rho0, & !nnz: required for MOM5 and previous versions. + ilb=G%isd, jlb=G%jsd, & + dzt=dzt,& !This is needed for the Mocsy method of carbonate system vars + tau=1, sosga=sosga, model_time=get_diag_time_end(CS%diag)) + endif + + !Output diagnostics via diag_manager for all tracers in this module +! if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//& +! "No tracer in the list.") +! call g_tracer_send_diag(CS%g_tracer_list, get_diag_time_end(CS%diag), tau=1) + !Niki: The problem with calling diagnostic outputs here is that this subroutine is called every dt_cpld + ! hence if dt_therm > dt_cpld we get output (and contribution to the mean) at times that tracers + ! had not been updated. + ! Moving this to the end of column physics subroutine fixes this issue. + +end subroutine MOM_generic_tracer_surface_state + +!ALL PE subroutine on Ocean! Due to otpm design the fluxes should be initialized like this on ALL PE's! +subroutine MOM_generic_flux_init(verbosity) + integer, optional, intent(in) :: verbosity !< A 0-9 integer indicating a level of verbosity. + + character(len=128), parameter :: sub_name = 'MOM_generic_flux_init' + type(g_tracer_type), pointer :: g_tracer_list,g_tracer,g_tracer_next + + if (.not. g_registered) then + call generic_tracer_register() + g_registered = .true. + endif + + call generic_tracer_get_list(g_tracer_list) + if (.NOT. associated(g_tracer_list)) then + call MOM_error(WARNING, trim(sub_name)// ": No generic tracer in the list.") + return + endif + + g_tracer=>g_tracer_list + do + + call g_tracer_flux_init(g_tracer, verbosity=verbosity) + + ! traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if (.NOT. associated(g_tracer_next)) exit + g_tracer=>g_tracer_next + + enddo + +end subroutine MOM_generic_flux_init + +subroutine MOM_generic_tracer_fluxes_accumulate(flux_tmp, weight) + type(forcing), intent(in) :: flux_tmp !< A structure containing pointers to + !! thermodynamic and tracer forcing fields. + real, intent(in) :: weight !< A weight for accumulating this flux [nondim] + + call generic_tracer_coupler_accumulate(flux_tmp%tr_fluxes, weight) + +end subroutine MOM_generic_tracer_fluxes_accumulate + +!> Copy the requested tracer into an array. +subroutine MOM_generic_tracer_get(name,member,array, CS) + character(len=*), intent(in) :: name !< Name of requested tracer. + character(len=*), intent(in) :: member !< The tracer element to return. + real, dimension(:,:,:), intent(out) :: array !< Array filled by this routine, in arbitrary units [A] + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + + ! Local variables + real, dimension(:,:,:), pointer :: array_ptr ! The tracer in the generic tracer structures, in + ! arbitrary units [A] + character(len=128), parameter :: sub_name = 'MOM_generic_tracer_get' + + call g_tracer_get_pointer(CS%g_tracer_list,name,member,array_ptr) + array(:,:,:) = array_ptr(:,:,:) + +end subroutine MOM_generic_tracer_get + +!> This subroutine deallocates the memory owned by this module. +subroutine end_MOM_generic_tracer(CS) + type(MOM_generic_tracer_CS), pointer :: CS !< Pointer to the control structure for this module. + + call generic_tracer_end() + + if (associated(CS)) then + deallocate(CS) + endif +end subroutine end_MOM_generic_tracer + +!---------------------------------------------------------------- +! Niki Zadeh +! +! +! William Cooke +! +! +! +! This module drives the generic version of tracers TOPAZ and CFC +! +!---------------------------------------------------------------- + +end module MOM_generic_tracer