Skip to content

Commit

Permalink
Merge branch 'dev/gfdl' into read_salin_in_ppt
Browse files Browse the repository at this point in the history
  • Loading branch information
adcroft authored Dec 15, 2023
2 parents 21f7d41 + 7d334f8 commit eb41d80
Show file tree
Hide file tree
Showing 5 changed files with 979 additions and 547 deletions.
3 changes: 1 addition & 2 deletions ac/Makefile.in
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,6 @@ SRC_DIRS = @SRC_DIRS@

-include Makefile.dep


# Generate Makefile from template
Makefile: @srcdir@/ac/Makefile.in config.status
./config.status
Expand All @@ -33,7 +32,7 @@ rwildcard=$(foreach d,$(wildcard $(1:=/*)),$(call rwildcard,$d,$2) $(filter $(su
.PHONY: depend
depend: Makefile.dep
Makefile.dep: $(MAKEDEP) $(call rwildcard,$(SRC_DIRS),*.h *.c *.inc *.F90)
$(PYTHON) $(MAKEDEP) -o Makefile.dep -e $(SRC_DIRS)
$(PYTHON) $(MAKEDEP) $(DEFS) -o Makefile.dep -e $(SRC_DIRS)


# Delete any files associated with configuration (including the Makefile).
Expand Down
93 changes: 88 additions & 5 deletions ac/makedep
Original file line number Diff line number Diff line change
Expand Up @@ -13,9 +13,16 @@ import sys
# Pre-compile re searches
re_module = re.compile(r"^ *module +([a-z_0-9]+)")
re_use = re.compile(r"^ *use +([a-z_0-9]+)")
re_cpp_define = re.compile(r"^ *# *define +[_a-zA-Z][_a-zA-Z0-9]")
re_cpp_undef = re.compile(r"^ *# *undef +[_a-zA-Z][_a-zA-Z0-9]")
re_cpp_ifdef = re.compile(r"^ *# *ifdef +[_a-zA-Z][_a-zA-Z0-9]")
re_cpp_ifndef = re.compile(r"^ *# *ifndef +[_a-zA-Z][_a-zA-Z0-9]")
re_cpp_if = re.compile(r"^ *# *if +")
re_cpp_else = re.compile(r"^ *# *else")
re_cpp_endif = re.compile(r"^ *# *endif")
re_cpp_include = re.compile(r"^ *# *include *[<\"']([a-zA-Z_0-9\.]+)[>\"']")
re_f90_include = re.compile(r"^ *include +[\"']([a-zA-Z_0-9\.]+)[\"']")
re_program = re.compile(r"^ *[pP][rR][oO][gG][rR][aA][mM] +([a-zA-Z_0-9]+)")
re_program = re.compile(r"^ *program +([a-z_0-9]+)", re.IGNORECASE)
re_end = re.compile(r"^ *end *(module|procedure) ", re.IGNORECASE)
# NOTE: This excludes comments and tokens with substrings containing `function`
# or `subroutine`, but will fail if the keywords appear in other contexts.
Expand All @@ -26,7 +33,7 @@ re_procedure = re.compile(


def create_deps(src_dirs, skip_dirs, makefile, debug, exec_target, fc_rule,
link_externals, script_path):
link_externals, defines):
"""Create "makefile" after scanning "src_dis"."""

# Scan everything Fortran related
Expand Down Expand Up @@ -66,7 +73,7 @@ def create_deps(src_dirs, skip_dirs, makefile, debug, exec_target, fc_rule,
o2mods, o2uses, o2h, o2inc, o2prg, prg2o, mod2o = {}, {}, {}, {}, {}, {}, {}
externals, all_modules = [], []
for f in F90_files:
mods, used, cpp, inc, prg, has_externals = scan_fortran_file(f)
mods, used, cpp, inc, prg, has_externals = scan_fortran_file(f, defines)
# maps object file to modules produced
o2mods[object_file(f)] = mods
# maps module produced to object file
Expand Down Expand Up @@ -272,10 +279,16 @@ def nested_inc(inc_files, f2F):
return inc_files + sorted(set(hlst)), used_mods


def scan_fortran_file(src_file):
def scan_fortran_file(src_file, defines=None):
"""Scan the Fortran file "src_file" and return lists of module defined,
module used, and files included."""
module_decl, used_modules, cpp_includes, f90_includes, programs = [], [], [], [], []

cpp_defines = defines if defines is not None else []

cpp_macros = [define.split('=')[0] for define in cpp_defines]
cpp_group_stack = []

with io.open(src_file, 'r', errors='replace') as file:
lines = file.readlines()

Expand All @@ -285,7 +298,72 @@ def scan_fortran_file(src_file):
file_has_externals = False
# True if the file contains any external objects

cpp_exclude = False
# True if the parser excludes the subsequent lines

cpp_group_stack = []
# Stack of condition group exclusion states

for line in lines:
# Start of #ifdef condition group
match = re_cpp_ifdef.match(line)
if match:
cpp_group_stack.append(cpp_exclude)

# If outer group is excluding or macro is missing, then exclude
macro = line.lstrip()[1:].split()[1]
cpp_exclude = cpp_exclude or macro not in cpp_macros

# Start of #ifndef condition group
match = re_cpp_ifndef.match(line)
if match:
cpp_group_stack.append(cpp_exclude)

# If outer group is excluding or macro is present, then exclude
macro = line.lstrip()[1:].split()[1]
cpp_exclude = cpp_exclude or macro in cpp_macros

# Start of #if condition group
match = re_cpp_if.match(line)
if match:
cpp_group_stack.append(cpp_exclude)

# XXX: Don't attempt to parse #if statements, but store the state.
# if/endif stack. For now, assume that these always fail.
cpp_exclude = False

# Complement #else condition group
match = re_cpp_else.match(line)
if match:
# Reverse the exclude state, if there is no outer exclude state
outer_grp_exclude = cpp_group_stack and cpp_group_stack[-1]
cpp_exclude = not cpp_exclude or outer_grp_exclude

# Restore exclude state when exiting conditional block
match = re_cpp_endif.match(line)
if match:
cpp_exclude = cpp_group_stack.pop()

# Skip lines inside of false condition blocks
if cpp_exclude:
continue

# Activate a new macro (ignoring the value)
match = re_cpp_define.match(line)
if match:
new_macro = line.lstrip()[1:].split()[1]
cpp_macros.append(new_macro)

# Deactivate a macro
match = re_cpp_undef.match(line)
if match:
new_macro = line.lstrip()[1:].split()[1]
try:
cpp_macros.remove(new_macro)
except:
# Ignore missing macros (for now?)
continue

match = re_module.match(line.lower())
if match:
if match.group(1) not in 'procedure': # avoid "module procedure" statements
Expand Down Expand Up @@ -404,8 +482,13 @@ parser.add_argument(
action='append',
help="Skip directory in source code search."
)
parser.add_argument(
'-D', '--define',
action='append',
help="Apply preprocessor define macros (of the form -DMACRO[=value])",
)
args = parser.parse_args()

# Do the thing
create_deps(args.path, args.skip, args.makefile, args.debug, args.exec_target,
args.fc_rule, args.link_externals, sys.argv[0])
args.fc_rule, args.link_externals, args.define)
171 changes: 20 additions & 151 deletions src/core/MOM_continuity.F90
Original file line number Diff line number Diff line change
Expand Up @@ -3,159 +3,28 @@ module MOM_continuity

! This file is part of MOM6. See LICENSE.md for the license.

use MOM_continuity_PPM, only : continuity_PPM, continuity_PPM_init
use MOM_continuity_PPM, only : continuity_PPM_stencil
use MOM_continuity_PPM, only : continuity_PPM_CS
use MOM_diag_mediator, only : time_type, diag_ctrl
use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_pe
use MOM_file_parser, only : get_param, log_version, param_file_type
use MOM_string_functions, only : uppercase
use MOM_grid, only : ocean_grid_type
use MOM_open_boundary, only : ocean_OBC_type
use MOM_unit_scaling, only : unit_scale_type
use MOM_variables, only : BT_cont_type, porous_barrier_type
use MOM_verticalGrid, only : verticalGrid_type
use MOM_continuity_PPM, only : continuity=>continuity_PPM
use MOM_continuity_PPM, only : continuity_stencil=>continuity_PPM_stencil
use MOM_continuity_PPM, only : continuity_init=>continuity_PPM_init
use MOM_continuity_PPM, only : continuity_CS=>continuity_PPM_CS
use MOM_continuity_PPM, only : continuity_fluxes, continuity_adjust_vel
use MOM_continuity_PPM, only : zonal_mass_flux, meridional_mass_flux
use MOM_continuity_PPM, only : zonal_edge_thickness, meridional_edge_thickness
use MOM_continuity_PPM, only : continuity_zonal_convergence, continuity_merdional_convergence
use MOM_continuity_PPM, only : zonal_flux_thickness, meridional_flux_thickness
use MOM_continuity_PPM, only : zonal_BT_mass_flux, meridional_BT_mass_flux
use MOM_continuity_PPM, only : set_continuity_loop_bounds, cont_loop_bounds_type

implicit none ; private

#include <MOM_memory.h>

public continuity, continuity_init, continuity_stencil

!> Control structure for mom_continuity
type, public :: continuity_CS ; private
integer :: continuity_scheme !< Selects the discretization for the continuity solver.
!! Valid values are:
!! - PPM - A directionally split piecewise parabolic reconstruction solver.
!! The default, PPM, seems most appropriate for use with our current
!! time-splitting strategies.
type(continuity_PPM_CS) :: PPM !< Control structure for mom_continuity_ppm
end type continuity_CS

integer, parameter :: PPM_SCHEME = 1 !< Enumerated constant to select PPM
character(len=20), parameter :: PPM_STRING = "PPM" !< String to select PPM

contains

!> Time steps the layer thicknesses, using a monotonically limited, directionally split PPM scheme,
!! based on Lin (1994).
subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, US, CS, OBC, pbv, uhbt, vhbt, &
visc_rem_u, visc_rem_v, u_cor, v_cor, BT_cont)
type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure.
type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure.
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1].
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(in) :: hin !< Initial layer thickness [H ~> m or kg m-2].
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
intent(inout) :: h !< Final layer thickness [H ~> m or kg m-2].
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
intent(out) :: uh !< Volume flux through zonal faces =
!! u*h*dy [H L2 T-1 ~> m3 s-1 or kg s-1].
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
intent(out) :: vh !< Volume flux through meridional faces =
!! v*h*dx [H L2 T-1 ~> m3 s-1 or kg s-1].
real, intent(in) :: dt !< Time increment [T ~> s].
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(continuity_CS), intent(in) :: CS !< Control structure for mom_continuity.
type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
type(porous_barrier_type), intent(in) :: pbv !< porous barrier fractional cell metrics
real, dimension(SZIB_(G),SZJ_(G)), &
optional, intent(in) :: uhbt !< The vertically summed volume
!! flux through zonal faces [H L2 T-1 ~> m3 s-1 or kg s-1].
real, dimension(SZI_(G),SZJB_(G)), &
optional, intent(in) :: vhbt !< The vertically summed volume
!! flux through meridional faces [H L2 T-1 ~> m3 s-1 or kg s-1].
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
optional, intent(in) :: visc_rem_u !< Both the fraction of
!! zonal momentum that remains after a time-step of viscosity, and the fraction of a time-step's
!! worth of a barotropic acceleration that a layer experiences after viscosity is applied [nondim].
!! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
optional, intent(in) :: visc_rem_v !< Both the fraction of
!! meridional momentum that remains after a time-step of viscosity, and the fraction of a time-step's
!! worth of a barotropic acceleration that a layer experiences after viscosity is applied [nondim].
!! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
optional, intent(out) :: u_cor !< The zonal velocities that
!! give uhbt as the depth-integrated transport [L T-1 ~> m s-1].
real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
optional, intent(out) :: v_cor !< The meridional velocities that
!! give vhbt as the depth-integrated transport [L T-1 ~> m s-1].
type(BT_cont_type), &
optional, pointer :: BT_cont !< A structure with elements
!! that describe the effective open face areas as a function of barotropic flow.

if (present(visc_rem_u) .neqv. present(visc_rem_v)) call MOM_error(FATAL, &
"MOM_continuity: Either both visc_rem_u and visc_rem_v or neither"// &
" one must be present in call to continuity.")
if (present(u_cor) .neqv. present(v_cor)) call MOM_error(FATAL, &
"MOM_continuity: Either both u_cor and v_cor or neither"// &
" one must be present in call to continuity.")

if (CS%continuity_scheme == PPM_SCHEME) then
call continuity_PPM(u, v, hin, h, uh, vh, dt, G, GV, US, CS%PPM, OBC, pbv, uhbt, vhbt, &
visc_rem_u, visc_rem_v, u_cor, v_cor, BT_cont=BT_cont)
else
call MOM_error(FATAL, "continuity: Unrecognized value of continuity_scheme")
endif

end subroutine continuity

!> Initializes continuity_cs
subroutine continuity_init(Time, G, GV, US, param_file, diag, CS)
type(time_type), target, intent(in) :: Time !< Current model time.
type(ocean_grid_type), intent(in) :: G !< Ocean grid structure.
type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure.
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
type(param_file_type), intent(in) :: param_file !< Parameter file handles.
type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure.
type(continuity_CS), intent(inout) :: CS !< Control structure for mom_continuity.

! This include declares and sets the variable "version".
# include "version_variable.h"
character(len=40) :: mdl = "MOM_continuity" ! This module's name.
character(len=20) :: tmpstr

! Read all relevant parameters and write them to the model log.
call log_version(param_file, mdl, version, "")
call get_param(param_file, mdl, "CONTINUITY_SCHEME", tmpstr, &
"CONTINUITY_SCHEME selects the discretization for the "//&
"continuity solver. The only valid value currently is: \n"//&
"\t PPM - use a positive-definite (or monotonic) \n"//&
"\t piecewise parabolic reconstruction solver.", &
default=PPM_STRING)

tmpstr = uppercase(tmpstr) ; CS%continuity_scheme = 0
select case (trim(tmpstr))
case (PPM_STRING) ; CS%continuity_scheme = PPM_SCHEME
case default
call MOM_mesg('continuity_init: CONTINUITY_SCHEME ="'//trim(tmpstr)//'"', 0)
call MOM_mesg("continuity_init: The only valid value is currently "// &
trim(PPM_STRING), 0)
call MOM_error(FATAL, "continuity_init: Unrecognized setting "// &
"#define CONTINUITY_SCHEME "//trim(tmpstr)//" found in input file.")
end select

if (CS%continuity_scheme == PPM_SCHEME) then
call continuity_PPM_init(Time, G, GV, US, param_file, diag, CS%PPM)
endif

end subroutine continuity_init


!> continuity_stencil returns the continuity solver stencil size
function continuity_stencil(CS) result(stencil)
type(continuity_CS), intent(in) :: CS !< Module's control structure.
integer :: stencil !< The continuity solver stencil size with the current settings.

stencil = 1

if (CS%continuity_scheme == PPM_SCHEME) then
stencil = continuity_PPM_stencil(CS%PPM)
endif
end function continuity_stencil
! These are direct pass-throughs of routines in continuity_PPM
public continuity, continuity_init, continuity_stencil, continuity_CS
public continuity_fluxes, continuity_adjust_vel
public zonal_mass_flux, meridional_mass_flux
public zonal_edge_thickness, meridional_edge_thickness
public continuity_zonal_convergence, continuity_merdional_convergence
public zonal_flux_thickness, meridional_flux_thickness
public zonal_BT_mass_flux, meridional_BT_mass_flux
public set_continuity_loop_bounds, cont_loop_bounds_type

end module MOM_continuity
Loading

0 comments on commit eb41d80

Please sign in to comment.