Commit abf8ed35 authored by sparonuz's avatar sparonuz
Browse files
parents 3ee17b6a b8d160b8
vault.pkl
\ No newline at end of file
vault.pkl
Scripts/tmp_processedsources/*
## PyCharm .gitignore
.idea/*
__pycache__/*
venv/*
## PyCharm .gitignore end
# Script that automatically creates different versions of a routine with different precisions
# Import classes and functions from RPE_Utils
from AutoRPE.UtilsRPE.MainFunctions import parse_sources, postprocess_sources, assign_sbits_to_rpe_variables, \
create_read_precisions_module_and_add_to_file, preprocess_sources
from AutoRPE.UtilsRPE.Functions import load_vault, remove_comments
from AutoRPE.UtilsRPE.CurrentBlock import Current_block
# It creates all the possible combinations of the real variables.
#
# If the subroutine has the reals foo and bar, it would generates the routines:
# - subroutine_name_0000
# * foo(dp)
# * bar(dp)
# - subroutine_name_0001
# * foo(sp)
# * bar(dp)
# - subroutine_name_0002
# * foo(dp)
# * bar(sp)
# - subroutine_name_0003
# * foo(sp)
# * bar(sp)
#
# Import classes and functions from UtilsRPE
from AutoRPE.UtilsRPE.CurrentBlock import CurrentBlock
from AutoRPE.UtilsRPE.Inserter import replace_real_declaration
from os import listdir, mkdir
from os.path import join, isdir, isfile
from AutoRPE.UtilsRPE.SourceManager import preprocess_sources, parse_sources, load_vault
from AutoRPE.UtilsRPE.Cleaner import remove_comments
from AutoRPE.UtilsRPE.Finder import find_declaration_line
import AutoRPE.UtilsRPE.RegexPattern as re_pattern
from argparse import ArgumentParser
from os import mkdir
from os.path import join, isdir
from shutil import rmtree
import re
from math import factorial
import itertools
def get_command_line_arguments():
help_text = """
This script takes as arguments the input source (-s/--input-source)
and the name of a routine (-r/--routine) and creates different versions
of this routine accouting for all different precision combinations that can happen.
An additional argument for the vault (-v/--vault) allows to skip the
creation of the vault if this is provided. Optional output (-o/--output-folder) for creating
the files there, if not given the outputs are printed out.
"""
parser = ArgumentParser(usage=help_text)
parser.add_argument("-v", "--vault", dest="vault", required=False, help="Path to the vault")
parser.add_argument("-o", "--output-folder", dest="output_folder", required=False, help="Path to outputs folder")
parser.add_argument("-s", "--input-source", dest="input_source", required=True,
help="Path to the processed sources")
parser.add_argument("-r", "--routine", dest="routine_name", required=True,
help="Name of the routine to create different interfaces")
args = parser.parse_args()
return args.vault, args.input_source, args.routine_name, args.output_folder
# Few functions
def merge_declaration_lines(code_text):
lines = [ l for l in code_text.split("\n")]
lines = [l for l in code_text.split("\n")]
for _index, _line in enumerate(lines):
line = remove_comments(_line).strip()
if line.count("&") and line.count("::"):
......@@ -36,142 +75,120 @@ def merge_declaration_lines(code_text):
def replace_subroutine_name(old_name, new_name, code_text):
pattern = r"(subroutine +\b)%s\b" % old_name
pattern = re_pattern.subroutine_declaration % old_name
pat = re.compile(pattern, re.I)
replacement = r"\1 %s" % new_name
return pat.sub(replacement, code_text)
def find_declaration_line(variable, routine_code):
import re
from AutoRPE.UtilsRPE.Functions import in_which_block, remove_comments
import warnings
# Declaration pattern
# TODO: This might have issues with split lines
pattern = r"::.*\b%s\b" % variable.name
exception_pattern = r"&*\b%s\b" % variable.name
real_pattern = r"real.*::"
lines = [l for l in routine_code.split("\n")]
current_block = Current_block(vault[variable.procedure.module.name].main())
for _index, _line in enumerate(lines):
_line = remove_comments(_line).strip()
if re.search(pattern, _line):
if current_block.procedure == variable.procedure:
return _index
else:
current_block = in_which_block(_line, current_block, vault=vault)
print("Declaration of variable %s with id %i not found in module %s" % (variable.name,variable.id, variable.procedure.module.name))
return None
def obtain_rotine_source(routine, sources_path):
def obtain_routine_source(routine, sources_path):
source_file = join(sources_path, "%s.f90" % routine.module.name)
indices = []
with open(source_file) as f:
pattern = r"routine.*\b%s\b" % routine.name
pattern = re_pattern.subroutine_declaration % routine.name
lines = [l for l in f]
indices = []
for index, line in enumerate(lines):
line = remove_comments(line).strip()
# if it's the start or end of the routine
if re.search(pattern, line, re.I):
indices.append(index)
routine_code = "".join(lines[indices[0]:indices[1]+1])
routine_code = merge_declaration_lines(routine_code)
return routine_code
if __name__ == "__main__":
# Getting Command line options
from optparse import OptionParser
help_text = """
This script takes as arguments the input source (-s/--input-source)
and the name of a routine (-r/--routine) and creates different versions
of this routine accouting for all different precision combinations that can happen.
An additional argument for the vault (-v/--vault) allows to skip the
creation of the vault if this is provided.
"""
parser = OptionParser(usage=help_text)
parser.add_option("-v", "--vault", dest="vault", default=None, metavar="FILE")
parser.add_option("-s", "--input-source", dest="input_source", default=None, metavar="FILE")
parser.add_option("-r", "--routine", dest="routine_name", default=None)
(options, args) = parser.parse_args()
vault_path = options.vault
path_to_input_sources = options.input_source
routine_name = options.routine_name
return routine_code, indices
if path_to_input_sources is None:
print("\nArgument -s/--input-source and -r/--routine are mandatory\n")
parser.print_help()
exit(1)
if routine_name is None:
print("\nArgument -r/--routine is mandatory\n")
parser.print_help()
exit(1)
def replace_routine_name(old_name, new_name, code_text):
pattern = re_pattern.subroutine_declaration % old_name
pat = re.compile(pattern, re.I)
replacement = r"\1 %s" % new_name
return pat.sub(replacement, code_text)
def prepare_vault(vault_path=None):
path_to_preprocessed_sources = "./tmp_processedsources/"
# In case no vault has been provided it will create it from the sources in the input sources folder
if vault_path is None:
# Paths to source folders
# Path_to_input_sources
if not isdir(path_to_preprocessed_sources):
mkdir(path_to_preprocessed_sources)
# List of files that should not be processed, as red from "list_of_files_to_keep_unmodified.txt"
preprocess_blacklist = listdir(path_to_input_sources)
mkdir(path_to_preprocessed_sources)
# Format text + replace real declarations with type(rpe_var)
preprocess_sources(path_to_input_sources, path_to_preprocessed_sources, blacklist=preprocess_blacklist)
preprocess_sources(path_to_input_sources, path_to_preprocessed_sources)
# Obtain source information
vault = parse_sources(path_to_preprocessed_sources, save_vault=True)
# Remove generated folder
rmtree(path_to_preprocessed_sources)
else:
vault = load_vault(vault_path)
routine = None
for procedure in vault.procedures:
if procedure.name == routine_name:
routine = procedure
break
if routine is None:
return vault
def get_routine(vault, routine_name=None):
# Find procedure with the routine name
routine = [p for p in vault.procedures if p.name == routine_name]
# Check if we found it
if len(routine) == 0:
print("Routine %s not found in vault." % routine_name)
exit(0)
routine_code = obtain_rotine_source(routine, path_to_input_sources)
routine_lines = routine_code.split("\n")
sensitive_variables = [v for v in routine.variables if v.is_dummy_argument and v.type == "real"]
# All possible combinations of precisions
precision_combinations = list(itertools.product(["dp","sp"], repeat=len(sensitive_variables)))
return routine[0]
def update_line(version_lines, variable, filepath, vault, routines_indices, variable_precision):
# Find the declaration/original lines
d_lines = find_declaration_line(variable, filepath, vault)
# Search in declaration lines
for d_line in d_lines:
# Update line if it's between routine lines routine
if d_line >= routines_indices[0] and d_line <= routines_indices[1]:
original_line = version_lines[d_line - routines_indices[0]]
new_line = replace_real_declaration(variable, original_line, variable_precision)
version_lines[d_line - routines_indices[0]] = new_line
def generate_routine_versions(precision_combinations, sensitive_variables, vault, routine_lines, routine, filepath, indices):
code_versions = {}
interface = "module procedure "
for c_index, combination in enumerate(precision_combinations):
version_lines = routine_lines[:]
new_routine_name = "%s_%04i" % (routine.name, c_index)
for index, variable in enumerate(sensitive_variables):
variable_precision = combination[index]
dec_line = find_declaration_line(variable, routine_code)
original_line = version_lines[dec_line]
new_line = replace_real_declaration(variable, original_line, variable_precision)
version_lines[dec_line] = new_line
version_text = "\n".join(version_lines).strip()
version_text = replace_routine_name(routine.name, new_routine_name,version_text)
code_versions[c_index] = version_text
if c_index == 0:
interface += "%s" % new_routine_name
else:
interface += ", %s" % new_routine_name
interface = "interface %s\n %s\nend interface %s" % (routine.name, interface, routine.name)
print(interface)
for key, code in code_versions.items():
print(code)
# Deleting temporary folders
if isdir(path_to_preprocessed_sources):
rmtree(path_to_preprocessed_sources)
new_routine_name = "%s_%04i" % (routine.name, c_index)
# Search and update real declaration
for index, variable in enumerate(sensitive_variables):
update_line(routine_lines, variable, filepath, vault, indices, combination[index])
# Join the routine lines and uodate nane
version_text = "\n".join(routine_lines).strip()
version_text = replace_routine_name(routine.name, new_routine_name, version_text)
code_versions[new_routine_name] = version_text
return code_versions
def write_routine_versions(code_versions, output_folder):
if output_folder is None:
# If there's no folder, print the outputs to terminal
for key, code in code_versions.items():
print(code)
else:
# Create folder if it doesn't exists
if not isdir(output_folder):
mkdir(output_folder)
# Create a file for each subroutine interface
for key, code in code_versions.items():
with open(output_folder + key, 'w') as f:
f.write(code)
if __name__ == "__main__":
# Get command line arguments
vault_path, path_to_input_sources, routine_name, output_folder = get_command_line_arguments()
# Get the vault
vault = prepare_vault(vault_path)
# Get routine from vault
routine = get_routine(vault, routine_name)
# Get the routine code
routine_code, indices = obtain_routine_source(routine, path_to_input_sources)
# All possible combinations of precisions x sensitive variables
sensitive_variables = [v for v in routine.variables if v.is_dummy_argument and v.type == "double"]
precision_combinations = list(itertools.product(["dp", "sp"], repeat=len(sensitive_variables)))
# Get filepath
filepath = join(path_to_input_sources, "%s.f90" % routine.module.name)
# Generate the different routine versions
code_versions = generate_routine_versions(precision_combinations, sensitive_variables, vault, routine_code.split("\n"), routine, filepath, indices)
write_routine_versions(code_versions, output_folder)
import os
import unittest
from shutil import rmtree
class TestInterfaceGenerator(unittest.TestCase):
def _changeWorkingDirToCurrent(self):
# Get script absoult path and set it as current work dir
abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)
@classmethod
def setUpClass(cls) -> None:
cls._changeWorkingDirToCurrent(cls)
params = "-s pp_vault_files/ -r dia_hth_dep -o output_folder/"
os.system("python3.7 ../../Scripts/InterfaceGenerator.py " + params)
@classmethod
def tearDownClass(cls) -> None:
if os.path.exists("vault.pkl"):
os.remove("vault.pkl")
if os.path.exists("output_folder"):
rmtree("output_folder")
def test_expect_vault(self):
self.assertEqual(os.path.exists("vault.pkl"), True, msg="The vault.pkl haven't been created")
def test_expect_outputs_folder(self):
self.assertEqual(os.path.exists("output_folder"), True, msg="The output_folder haven't been created")
def test_outputs_contents(self):
from itertools import product
filename="dia_hth_dep_%s"
with open('helper_files/' + filename % ("template"), 'r') as f:
template = f.read()
filename = filename % ("000%i")
precision_combinations = [('dp', 'dp'), ('dp', 'sp'), ('sp', 'dp'), ('sp', 'sp')]
for idx, comb in enumerate(precision_combinations):
with open("output_folder/" + filename % (idx)) as f:
# Set the params from template correctly
txt = template % (idx, comb[0], comb[1], idx)
# split the text by rows
code = f.read().split('\n')
txt = txt.split('\n')
# Check row by row
self.assertEqual([row.strip() for row in code], [row.strip() for row in txt])
if __name__ == '__main__':
unittest.main()
SUBROUTINE dia_hth_dep_000%s( Kmm, ptem, pdept )
!
INTEGER , INTENT(in) :: Kmm ! ocean time level index
REAL(%s), INTENT(in) :: ptem
REAL(%s), DIMENSION(jpi,jpj), INTENT(out) :: pdept
!
INTEGER :: ji, jj, jk, iid
REAL(wp) :: zztmp, zzdep
INTEGER, DIMENSION(jpi,jpj) :: iktem
! --------------------------------------- !
! search deepest level above ptem !
! --------------------------------------- !
iktem(:,:) = 1
DO jk = 1, jpkm1 ; DO jj = 1, jpj ; DO ji = 1, jpi
zztmp = ts(ji,jj,jk,jp_tem,Kmm)
IF( zztmp >= ptem ) iktem(ji,jj) = jk
END DO ; END DO ; END DO
! ------------------------------- !
! Depth of ptem isotherm !
! ------------------------------- !
DO jj = 1, jpj ; DO ji = 1, jpi
!
zzdep = gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ! depth of the ocean bottom
!
iid = iktem(ji,jj)
IF( iid /= 1 ) THEN
zztmp = gdept(ji,jj,iid ,Kmm) & ! linear interpolation
& + ( gdept(ji,jj,iid+1,Kmm) - gdept(ji,jj,iid,Kmm) ) &
& * ( 20.*tmask(ji,jj,iid+1) - ts(ji,jj,iid,jp_tem,Kmm) ) &
& / ( ts(ji,jj,iid+1,jp_tem,Kmm) - ts(ji,jj,iid,jp_tem,Kmm) + (1.-tmask(ji,jj,1)) )
pdept(ji,jj) = MIN( zztmp , zzdep) * tmask(ji,jj,1) ! bound by the ocean depth
ELSE
pdept(ji,jj) = 0._wp
ENDIF
END DO ; END DO
!
END SUBROUTINE dia_hth_dep_000%s
\ No newline at end of file
MODULE diahth
!!======================================================================
!! *** MODULE diahth ***
!! Ocean diagnostics: thermocline and 20 degree depth
!!======================================================================
!! History : OPA ! 1994-09 (J.-P. Boulanger) Original code
!! ! 1996-11 (E. Guilyardi) OPA8
!! ! 1997-08 (G. Madec) optimization
!! ! 1999-07 (E. Guilyardi) hd28 + heat content
!! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module
!! 3.2 ! 2009-07 (S. Masson) hc300 bugfix + cleaning + add new diag
!!----------------------------------------------------------------------
!! dia_hth : Compute varius diagnostics associated with the mixed layer
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers
USE dom_oce ! ocean space and time domain
USE phycst ! physical constants
!
USE in_out_manager ! I/O manager
USE lib_mpp ! MPP library
USE iom ! I/O library
USE timing ! preformance summary
IMPLICIT NONE
PRIVATE
PUBLIC dia_hth ! routine called by step.F90
PUBLIC dia_hth_alloc ! routine called by nemogcm.F90
LOGICAL, SAVE :: l_hth !: thermocline-20d depths flag
! note: following variables should move to local variables once iom_put is always used
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hth !: depth of the max vertical temperature gradient [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hd20 !: depth of 20 C isotherm [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hd26 !: depth of 26 C isotherm [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hd28 !: depth of 28 C isotherm [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: htc3 !: heat content of first 300 m [W]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: htc7 !: heat content of first 700 m [W]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: htc20 !: heat content of first 2000 m [W]
!! * Substitutions
!!----------------------------------------------------------------------
!! *** domzgr_substitute.h90 ***
!!----------------------------------------------------------------------
!! ** purpose : substitute fsdep. and fse.., the vert. depth and scale
!! factors depending on the vertical coord. used, using CPP macro.
!!----------------------------------------------------------------------
!! History : 4.2 ! 2020-02 (S. Techene, G. Madec) star coordinate
!!----------------------------------------------------------------------
!! NEMO/OCE 4.2 , NEMO Consortium (2020)
!! $Id$
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id$
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
FUNCTION dia_hth_alloc()
!!---------------------------------------------------------------------
INTEGER :: dia_hth_alloc
!!---------------------------------------------------------------------
!
ALLOCATE( hth(jpi,jpj), hd20(jpi,jpj), hd26(jpi,jpj), hd28(jpi,jpj), &
& htc3(jpi,jpj), htc7(jpi,jpj), htc20(jpi,jpj), STAT=dia_hth_alloc )
!
CALL mpp_sum ( 'diahth', dia_hth_alloc )
IF(dia_hth_alloc /= 0) CALL ctl_stop( 'STOP', 'dia_hth_alloc: failed to allocate arrays.' )
!
END FUNCTION dia_hth_alloc
SUBROUTINE dia_hth( kt, Kmm )
!!---------------------------------------------------------------------
!! *** ROUTINE dia_hth ***
!!
!! ** Purpose : Computes
!! the mixing layer depth (turbocline): avt = 5.e-4
!! the depth of strongest vertical temperature gradient
!! the mixed layer depth with density criteria: rho = rho(10m or surf) + 0.03(or 0.01)
!! the mixed layer depth with temperature criteria: abs( tn - tn(10m) ) = 0.2
!! the top of the thermochine: tn = tn(10m) - ztem2
!! the pycnocline depth with density criteria equivalent to a temperature variation
!! rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
!! the barrier layer thickness
!! the maximal verical inversion of temperature and its depth max( 0, max of tn - tn(10m) )
!! the depth of the 20 degree isotherm (linear interpolation)
!! the depth of the 28 degree isotherm (linear interpolation)
!! the heat content of first 300 m
!!
!! ** Method :
!!-------------------------------------------------------------------
INTEGER, INTENT( in ) :: kt ! ocean time-step index
INTEGER, INTENT( in ) :: Kmm ! ocean time level index
!!
INTEGER :: ji, jj, jk ! dummy loop arguments
REAL(wp) :: zrho3 = 0.03_wp ! density criterion for mixed layer depth
REAL(wp) :: zrho1 = 0.01_wp ! density criterion for mixed layer depth
REAL(wp) :: ztem2 = 0.2_wp ! temperature criterion for mixed layer depth
REAL(wp) :: zztmp, zzdep ! temporary scalars inside do loop
REAL(wp) :: zu, zv, zw, zut, zvt ! temporary workspace
REAL(wp), DIMENSION(jpi,jpj) :: zabs2 ! MLD: abs( tn - tn(10m) ) = ztem2
REAL(wp), DIMENSION(jpi,jpj) :: ztm2 ! Top of thermocline: tn = tn(10m) - ztem2
REAL(wp), DIMENSION(jpi,jpj) :: zrho10_3 ! MLD: rho = rho10m + zrho3
REAL(wp), DIMENSION(jpi,jpj) :: zpycn ! pycnocline: rho = rho10m + (dr/dT)(T,S,10m)*(-0.2 degC)
REAL(wp), DIMENSION(jpi,jpj) :: ztinv ! max of temperature inversion
REAL(wp), DIMENSION(jpi,jpj) :: zdepinv ! depth of temperature inversion
REAL(wp), DIMENSION(jpi,jpj) :: zrho0_3 ! MLD rho = rho(surf) = 0.03
REAL(wp), DIMENSION(jpi,jpj) :: zrho0_1 ! MLD rho = rho(surf) = 0.01
REAL(wp), DIMENSION(jpi,jpj) :: zmaxdzT ! max of dT/dz
REAL(wp), DIMENSION(jpi,jpj) :: zdelr ! delta rho equivalent to deltaT = 0.2
!!----------------------------------------------------------------------
IF( ln_timing ) CALL timing_start('dia_hth')
IF( kt == nit000 ) THEN
l_hth = .FALSE.
IF( iom_use( 'mlddzt' ) .OR. iom_use( 'mldr0_3' ) .OR. iom_use( 'mldr0_1' ) .OR. &
& iom_use( 'mld_dt02' ) .OR. iom_use( 'topthdep' ) .OR. iom_use( 'mldr10_3' ) .OR. &
& iom_use( '20d' ) .OR. iom_use( '26d' ) .OR. iom_use( '28d' ) .OR. &
& iom_use( 'hc300' ) .OR. iom_use( 'hc700' ) .OR. iom_use( 'hc2000' ) .