Changeset 6620 for TOOLS

Timestamp:

09/14/23 15:08:35 (8 months ago)

Author:

omamce

Message:

WATER_BUDET (by OM) : OK for SECHIBA on LMD grid. Still need improvment for ICO.

Location:

TOOLS/WATER_BUDGET

Files:

• ATM_waterbudget.py (modified) (37 diffs)
• OCE_waterbudget.py (modified) (21 diffs)
• WaterUtils.py (modified) (5 diffs)
• libIGCM.py (added)

TOOLS/WATER_BUDGET/ATM_waterbudget.py

@@ -23 +23 @@
 ###
 ## Import system modules
-import sys, os, shutil, subprocess, platform
-import numpy as np
+import sys, os, shutil#, subprocess, platform
 import configparser, re
+
+## Import needed scientific modules
+import numpy as np, xarray as xr
 
 # Check python version
@@ -34 +36 @@
 ## Import local module
 import WaterUtils as wu
+import libIGCM_sys
+import nemo, lmdz
 
 ## Creates parser for reading .ini input file
+## -------------------------------------------
 config = configparser.ConfigParser ( interpolation=configparser.ExtendedInterpolation() )
 config.optionxform = str # To keep capitals
 
 ## Experiment parameters
+## ---------------------
 ATM=None ; ATM_HIS='latlon' ; SRF_HIS='latlon' ; RUN_HIS='latlon' ; ORCA=None ; NEMO=None ; OCE_relax=False
 OCE_icb=False ; Coupled=False ; Routing=None
@@ -45 +51 @@
 
 ##
-ARCHIVE=None ; STORAGE=None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None
-TmpDir=None ; RunDir=None ; FileOut=None
+ARCHIVE=None ; STORAGE=None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None ; TmpDir=None
+FileDir=None ; FileOut=None
 dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None
 FileCommon=None ; file_ATM_his=None ; file_SRF_his=None ; file_RUN_his=None
@@ -61 +67 @@
 ContinueOnError=False ; ErrorCount=0 ; SortIco = False
 
-# Read command line arguments
+## Read command line arguments
+## ---------------------------
 print ( "Name of Python script:", sys.argv[0] )
 IniFile = sys.argv[1]
-# Text existence of IniFile
+
+# Test existence of IniFile
 if not os.path.exists (IniFile ) :
     raise FileExistsError ( f"File not found : {IniFile = }" )
@@ -75 +83 @@
 FullIniFile = 'full_' + IniFile
 
+## Reading config.card if possible
+## -------------------------------
+ConfigCard = None
+
+if 'Experiment' in config.keys ()  : ## Read Experiment on Config file if possible
+    if 'ConfigCard' in config['Experiment'].keys () :
+        ConfigCard = config['Experiment']['ConfigCard']
+        print ( f'{ConfigCard=}' )
+
+if ConfigCard : ## Read config card if it exists
+    # Text existence of ConfigCard
+    if os.path.exists ( ConfigCard ) :
+        print ( f'Reading Config Card : {ConfigCard}' )
+        ## Creates parser for reading .ini input file
+        MyReader = configparser.ConfigParser (interpolation=configparser.ExtendedInterpolation() )
+        MyReader.optionxform = str # To keep capitals
+        
+        MyReader.read (ConfigCard)
+
+        for VarName in ['JobName', 'ExperimentName', 'SpaceName', 'LongName', 'ModelName', 'TagName'] :
+            if VarName in MyReader['UserChoices'].keys() :
+                locals()[VarName] = MyReader['UserChoices'][VarName]
+                exec  ( f'{VarName} = wu.setBool ({VarName})' )
+                exec  ( f'{VarName} = wu.setNum  ({VarName})' )
+                exec  ( f'{VarName} = wu.setNone ({VarName})' )
+                exec  ( f'wu.{VarName} = {VarName}' )
+                print ( f'    {VarName:21} set to : {locals()[VarName]:}' )
+                
+        for VarName in ['PackFrequency'] :
+            if VarName in MyReader['Post'].keys() :
+                locals()[VarName] = MyReader['Post'][VarName]
+                exec  ( f'{VarName} = wu.setBool ({VarName})' )
+                exec  ( f'{VarName} = wu.setNum  ({VarName})' )
+                exec  ( f'{VarName} = wu.setNone ({VarName})' )
+                exec  ( f'wu.{VarName} = {VarName}' )
+                print ( f'    {VarName:21} set to : {locals()[VarName]:}' )
+    else :
+        raise FileExistsError ( f"File not found : {ConfigCard = }" )
+                
 ## Reading config file
+## -------------------
 for Section in ['Config', 'Experiment', 'libIGCM', 'Files', 'Physics' ] :
    if Section in config.keys () : 
@@ -85 +133 @@
             exec  ( f'{VarName} = wu.setNone ({VarName})' )
             exec  ( f'wu.{VarName} = {VarName}' )
-            print ( f'    {VarName:21} set to : {locals()[VarName]:}' )
+            print ( f'    {VarName:21} set to : {locals()[VarName]}' )
             #exec ( f'del {VarName}' )
 
 print ( f'\nConfig file readed : {IniFile} ' )
-            
-##-- Some physical constants
-if wu.unDefined ( 'Ra' )           : Ra          = 6366197.7236758135 #-- Earth Radius (m)
-if wu.unDefined ( 'Grav' )         : Grav        = 9.81               #-- Gravity (m^2/s)
-if wu.unDefined ( 'ICE_rho_ice' )  : ICE_rho_ice = 917.0              #-- Ice volumic mass (kg/m3) in LIM3 or SI3
-if wu.unDefined ( 'ICE_rho_sno')   : ICE_rho_sno = 330.0              #-- Snow volumic mass (kg/m3) in LIM3 or SI3
-if wu.unDefined ( 'OCE_rho_liq' )  : OCE_rho_liq = 1026.0             #-- Ocean water volumic mass (kg/m3) in NEMO
-if wu.unDefined ( 'ATM_rho' )      : ATM_rho     = 1000.0             #-- Water volumic mass in atmosphere (kg/m^3)
-if wu.unDefined ( 'SRF_rho' )      : SRF_rho     = 1000.0             #-- Water volumic mass in surface reservoir (kg/m^3)
-if wu.unDefined ( 'RUN_rho' )      : RUN_rho     = 1000.0             #-- Water volumic mass of rivers (kg/m^3)
-if wu.unDefined ( 'ICE_rho_pnd' )  : ICE_rho_pnd = 1000.              #-- Water volumic mass in ice ponds (kg/m^3)
-if wu.unDefined ( 'YearLength' )   : YearLength  = 365.25 * 86400.    #-- Year length (s) - Use only to compute drif in approximate unit 
-            
-##-- Set libIGCM and machine dependant values
+
+## Some physical constants
+## =======================
+if wu.unDefined ( 'Ra' )           : Ra          = wu.Ra           #-- Earth Radius (m)
+if wu.unDefined ( 'Grav' )         : Grav        = wu.Grav         #-- Gravity (m^2/s
+if wu.unDefined ( 'ICE_rho_ice' )  : ICE_rho_ice = wu.ICE_rho_ice  #-- Ice volumic mass (kg/m3) in LIM3
+if wu.unDefined ( 'ICE_rho_sno')   : ICE_rho_sno = wu.ICE_rho_sno  #-- Snow volumic mass (kg/m3) in LIM3
+if wu.unDefined ( 'OCE_rho_liq' )  : OCE_rho_liq = wu.OCE_rho_liq  #-- Ocean water volumic mass (kg/m3) in NEMO
+if wu.unDefined ( 'ATM_rho' )      : ATM_rho     = wu.ATM_rho      #-- Water volumic mass in atmosphere (kg/m^3)
+if wu.unDefined ( 'SRF_rho' )      : SRF_rho     = wu.SRF_rho      #-- Water volumic mass in surface reservoir (kg/m^3)
+if wu.unDefined ( 'RUN_rho' )      : RUN_rho     = wu.RUN_rho      #-- Water volumic mass of rivers (kg/m^3)
+if wu.unDefined ( 'ICE_rho_pnd' )  : ICE_rho_pnd = wu.ICE_rho_pnd  #-- Water volumic mass in ice ponds (kg/m^3)
+if wu.unDefined ( 'YearLength' )   : YearLength  = wu.YearLength   #-- Year length (s)
+    
+## Set libIGCM and machine dependant values
+## ----------------------------------------
 if not 'Files' in config.keys () : config['Files'] = {}
 
@@ -113 +163 @@
 ICO  = ( 'ICO' in wu.ATM )
 LMDZ = ( 'LMD' in wu.ATM )
-    
-with open ('SetLibIGCM.py') as f: exec ( f.read() )
+
+mm = libIGCM_sys.config ( TagName=TagName, SpaceName=SpaceName, ExperimentName=ExperimentName, JobName=JobName, User=User, Group=Group,
+             ARCHIVE=None, SCRATCHDIR=None, STORAGE=None, R_IN=None, R_OUT=None, R_FIG=None, rebuild=None, TmpDir=None, 
+             R_SAVE=None, R_FIGR=None, R_BUFR=None, R_BUF_KSH=None, REBUILD_DIR=None, POST_DIR=None )
+globals().update(mm)
+
 config['Files']['TmpDir'] = TmpDir
-
-if libIGCM :
-    config['libIGCM'] = { 'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'SCRATCHDIR':SCRATCHDIR, 'R_IN':R_IN, 'rebuild':rebuild }
-
-##-- Import specific module
-import nemo, lmdz
-##-- Now import needed scientific modules
-import xarray as xr
-
-##- Output file with water budget diagnostics
-if wu.unDefined ( 'FileOut' ) : FileOut = f'ATM_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out'
+config['libIGCM'] = { 'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'TmpDir':TmpDir, 'R_IN':R_IN, 'rebuild':rebuild }
+
+## Defines begining and end of experiment
+## --------------------------------------
+if wu.unDefined ( 'DateBegin' ) :
+    DateBegin = f'{YearBegin}0101'
+    config['Experiment']['DateBegin'] = DateBegin
+else :
+    YearBegin, MonthBegin, DayBegin = wu.DateSplit ( DateBegin )
+    DateBegin = wu.FormatToGregorian ( DateBegin)
+    config['Experiment']['YearBegin'] = YearBegin
+
+if wu.unDefined ( 'DateEnd' )   :
+    DateEnd   = f'{YearEnd}1231'
+    config['Experiment']['DateEnd']   = DateEnd
+else :
+    YearEnd, MonthEnd, DayEnd = wu.DateSplit ( DateEnd )
+    DateEnd   = wu.FormatToGregorian ( DateEnd)
+
+if wu.unDefined ( 'PackFrequency' ) :
+    PackFrequency = YearEnd - YearBegin + 1
+    config['Experiment']['PackFrequency']   = f'{PackFrequency}'
+
+if type ( PackFrequency ) == str : 
+    if 'Y' in PackFrequency : PackFrequency = PackFrequency.replace ( 'Y', '') 
+    if 'M' in PackFrequency : PackFrequency = PackFrequency.replace ( 'M', '')
+    PackFrequency = int ( PackFrequency )
+    
+## Output file with water budget diagnostics
+## -----------------------------------------
+if wu.unDefined ( 'FileOut' ) : FileOut = f'ATM_waterbudget_{JobName}_{DateBegin}_{DateEnd}.out'
 config['Files']['FileOut'] = FileOut
 
 f_out = open ( FileOut, mode = 'w' )
 
-##- Useful functions
+## Useful functions
+## ----------------
 def kg2Sv    (val, rho=ATM_rho) :
     '''From kg to Sverdrup'''
@@ -141 +216 @@
 
 def var2prt (var, small=False, rho=ATM_rho) :
-    if small :  return var , kg2Sv(var, rho=rho)*1000., kg2myear(var, rho=rho)*1000.
+    if small :  return var , kg2Sv(var, rho=rho)*1000., kg2myear(var, rho=rho)*100#
     else     :  return var , kg2Sv(var, rho=rho)      , kg2myear(var, rho=rho)
 
 def prtFlux (Desc, var, Form='F', small=False, rho=ATM_rho, width=15) :
     if small :
-        if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} mSv | {:12.4f} mm/year "
-        if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} mSv | {:12.4e} mm/year "
+        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} mSv | {:12.4f} mm/year "
+        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} mSv | {:12.4e} mm/year "
     else :
-        if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} Sv  | {:12.4f} m/year  "
-        if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} Sv  | {:12.4e} m/year  "
-    echo ( (' {:>{width}} = ' +ff).format (Desc, *var2prt(var, small, rho=rho), width=width ) )
+        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} Sv  | {:12.4f} m/year  "
+        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} Sv  | {:12.4e} m/year  "
+    echo ( (' {:>{width}} = ' +ff).format (Desc, *var2prt(var, small=small, rho=rho), width=width ) )
     return None
 
@@ -162 +237 @@
     return None
     
-##- Set libIGCM directories
+## Set libIGCM directories
+## -----------------------
 if wu.unDefined ('R_OUT'      ) : R_OUT       = os.path.join ( ARCHIVE   , 'IGCM_OUT' )
 if wu.unDefined ('R_BUF'      ) : R_BUF       = os.path.join ( SCRATCHDIR, 'IGCM_OUT' )
@@ -173 +249 @@
 if wu.unDefined ('R_FIGR'     ) : R_FIGR      = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP )
 
-config['libIGCM'] = { 'R_OUT':R_OUT, 'R_BUF':R_BUF, 'L_EXP':L_EXP, 'R_BUFR':R_BUFR, 'POST_DIR':POST_DIR,
-                      'REBUILD_DIR':REBUILD_DIR, 'R_BUF_KSH':R_BUF_KSH, 'R_FIGR':R_FIGR }
-
-##- Set directory to extract files
-if wu.unDefined ( 'RunDir' ) : RunDir = os.path.join ( TmpDir, f'WATER_{JobName}_{YearBegin}_{YearEnd}' )
-config['Files']['RunDir'] = RunDir
-
-if not os.path.isdir ( RunDir ) : os.makedirs ( RunDir )
+config['libIGCM'].update ( { 'R_OUT':R_OUT, 'R_BUF':R_BUF, 'L_EXP':L_EXP, 'R_BUFR':R_BUFR, 'POST_DIR':POST_DIR,
+                      'REBUILD_DIR':REBUILD_DIR, 'R_BUF_KSH':R_BUF_KSH, 'R_FIGR':R_FIGR, 'rebuild':rebuild } )
+
+## Set directory to extract files
+## ------------------------------
+if wu.unDefined ( 'FileDir' ) : FileDir = os.path.join ( TmpDir, f'WATER_{JobName}' )
+config['Files']['FileDir'] = FileDir
+
+if not os.path.isdir ( FileDir ) : os.makedirs ( FileDir )
 
 ##- Set directories to rebuild ocean and ice restart files
-if wu.unDefined ( 'RunDirOCE' ) : RunDirOCE = os.path.join ( RunDir, 'OCE' )
-if wu.unDefined ( 'RunDirICE' ) : RunDirICE = os.path.join ( RunDir, 'ICE' )
-if not os.path.exists ( RunDirOCE ) : os.mkdir ( RunDirOCE )
-if not os.path.exists ( RunDirICE ) : os.mkdir ( RunDirICE )
+if wu.unDefined ( 'FileDirOCE' ) : FileDirOCE = os.path.join ( FileDir, 'OCE' )
+if wu.unDefined ( 'FileDirICE' ) : FileDirICE = os.path.join ( FileDir, 'ICE' )
+if not os.path.exists ( FileDirOCE ) : os.mkdir ( FileDirOCE )
+if not os.path.exists ( FileDirICE ) : os.mkdir ( FileDirICE )
 
 echo (' ')
-echo ( f'JobName   : {JobName}'   )
-echo ( f'Comment   : {Comment}'   )
-echo ( f'TmpDir    : {TmpDir}'    )
-echo ( f'RunDir    : {RunDir}'    )
-echo ( f'RunDirOCE : {RunDirOCE}' )
-echo ( f'RunDirICE : {RunDirICE}' )
+echo ( f'JobName    : {JobName}'   )
+echo ( f'Comment    : {Comment}'   )
+echo ( f'TmpDir     : {TmpDir}'    )
+echo ( f'FileDir     : {FileDir}'    )
+echo ( f'FileDirOCE  : {FileDirOCE}' )
+echo ( f'FileDirICE  : {FileDirICE}' )
 
 echo ( f'\nDealing with {L_EXP}'  )
 
-##-- Creates model output directory names
+## Creates model output directory names
+## ------------------------------------
 if Freq == "MO" : FreqDir = os.path.join ('Output' , 'MO' )
 if Freq == "SE" : FreqDir = os.path.join ('Analyse', 'SE' )
@@ -214 +292 @@
 ##-- Creates files names
 if wu.unDefined ( 'Period' ) :
-    if Freq == 'MO' : Period = f'{YearBegin}0101_{YearEnd}1231_1M'
-    if Freq == 'SE' : Period = f'SE_{YearBegin}0101_{YearEnd}1231_1M'
+    if Freq == 'MO' : Period = f'{DateBegin}_{DateEnd}_1M'
+    if Freq == 'SE' : Period = f'SE_{DateBegin}_{DateEnd}_1M'
     config['Files']['Period'] = Period
+
+config['Files']['DateBegin'] = DateBegin
+config['Files']['DateBegin'] = DateEnd
+       
+echo ( f'Period   : {Period}' )
+    
 if wu.unDefined ( 'FileCommon' ) :
     FileCommon = f'{JobName}_{Period}'
@@ -222 +306 @@
 
 if wu.unDefined ( 'Title' ) :
-    Title = f'{JobName} : {Freq} : {YearBegin}-01-01 - {YearEnd}-12-31'
+    Title = f'{JobName} : {Freq} : {DateBegin} - {DateEnd}'
     config['Files']['Title'] = Title
      
@@ -256 +340 @@
 if Routing == 'SIMPLE' : echo ( f'{file_RUN_his = }' )
 
-##-- Compute run length
+## Compute run length
+## ------------------
 dtime = ( d_ATM_his.time_counter_bounds.max() - d_ATM_his.time_counter_bounds.min() )
 echo ('\nRun length : {:8.2f} days'.format ( (dtime/np.timedelta64(1, "D")).values ) )
@@ -272 +357 @@
 
 ##-- Extract restart files from tar
-YearRes = YearBegin - 1              # Year of the restart of beginning of simulation
-YearPre = YearBegin - PackFrequency  # Year to find the tarfile of the restart of beginning of simulation
-
-config['Files']['YearPre'] = f'{YearBegin}'
-
-if wu.unDefined ( 'TarRestartPeriod_beg' ) : 
-    echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ')
-    TarRestartPeriod_beg = f'{YearPre}0101_{YearRes}1231'
+
+if wu.unDefined ('TarRestartDate_beg' ) : TarRestartDate_beg = wu.DateMinusOneDay ( DateBegin )
+if wu.unDefined ('TarRestartDate_end' ) : TarRestartDate_end = wu.FormatToGregorian ( DateEnd )
+
+if wu.unDefined ( 'TarRestartPeriod_beg' ) :
+
+    TarRestartPeriod_beg_DateEnd = TarRestartDate_beg
+    TarRestartPeriod_beg_DateBeg = wu.DateAddYear ( TarRestartPeriod_beg_DateEnd, -PackFrequency )
+    TarRestartPeriod_beg_DateBeg = wu.DatePlusOneDay ( TarRestartPeriod_beg_DateBeg )
+    
+    TarRestartPeriod_beg = f'{TarRestartPeriod_beg_DateBeg}_{TarRestartPeriod_beg_DateEnd}'
+    echo (f'Tar period for initial restart : {TarRestartPeriod_beg}')
     config['Files']['TarRestartPeriod_beg'] = TarRestartPeriod_beg
 
-if wu.unDefined ( 'TarRestartPeriod_end ' ) : 
-    YearPre = YearBegin - PackFrequency  # Year to find the tarfile of the restart of beginning of simulation
-    echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ')
-    TarRestartPeriod_end = f'{YearBegin}0101_{YearEnd}1231'
+if wu.unDefined ( 'TarRestartPeriod_end' ) :
+
+    TarRestartPeriod_end_DateEnd = TarRestartDate_end
+    TarRestartPeriod_end_DateBeg = wu.DateAddYear ( TarRestartPeriod_end_DateEnd, -PackFrequency )
+    TarRestartPeriod_end_DateBeg = wu.DatePlusOneDay ( TarRestartPeriod_end_DateBeg )
+     
+    TarRestartPeriod_end = f'{TarRestartPeriod_end_DateBeg}_{TarRestartPeriod_end_DateEnd}'
+    echo (f'Tar period for final restart : {TarRestartPeriod_end}')
     config['Files']['TarRestartPeriod_end'] = TarRestartPeriod_end
+
+echo (f'Restart dates - Start : {TarRestartPeriod_beg}  /  End : {TarRestartPeriod_end}')
 
 if wu.unDefined ( 'tar_restart_beg' ) :
@@ -316 +411 @@
 
 if wu.unDefined ( 'file_ATM_beg' ) : 
-    file_ATM_beg = f'{RunDir}/ATM_{JobName}_{YearRes}1231_restartphy.nc'
+    file_ATM_beg = f'{FileDir}/ATM_{JobName}_{TarRestartDate_beg}_restartphy.nc'
     config['Files']['file_ATM_beg'] = file_ATM_beg
 if wu.unDefined ( 'file_ATM_end' ) :
-    file_ATM_end = f'{RunDir}/ATM_{JobName}_{YearEnd}1231_restartphy.nc'
+    file_ATM_end = f'{FileDir}/ATM_{JobName}_{TarRestartDate_end}_restartphy.nc'
     config['Files']['file_ATM_end'] = file_ATM_end
 
@@ -326 +421 @@
 
 if wu.unDefined ( 'file_DYN_beg' ) :
-    if LMDZ : file_DYN_beg = f'{RunDir}/ATM_{JobName}_{YearRes}1231_restart.nc'
-    if ICO  : file_DYN_beg = f'{RunDir}/ICO_{JobName}_{YearRes}1231_restart.nc'
+    if LMDZ : file_DYN_beg = f'{FileDir}/ATM_{JobName}_{TarRestartDate_beg}_restart.nc'
+    if ICO  : file_DYN_beg = f'{FileDir}/ICO_{JobName}_{TarRestartDate_beg}_restart.nc'
     liste_beg.append (file_DYN_beg)
     config['Files']['file_DYN_beg'] = file_DYN_beg
     
 if wu.unDefined ( 'file_DYN_end' ) : 
-    if LMDZ : file_DYN_end = f'{RunDir}/ATM_{JobName}_{YearEnd}1231_restart.nc'
-    if ICO  : file_DYN_end = f'{RunDir}/ICO_{JobName}_{YearEnd}1231_restart.nc'
+    if LMDZ : file_DYN_end = f'{FileDir}/ATM_{JobName}_{TarRestartDate_end}_restart.nc'
+    if ICO  : file_DYN_end = f'{FileDir}/ICO_{JobName}_{TarRestartDate_end}_restart.nc'
     liste_end.append ( file_DYN_end )
     config['Files']['file_DYN_end'] = file_DYN_end
 
 if wu.unDefined ( 'file_SRF_beg' ) :
-    file_SRF_beg = f'{RunDir}/SRF_{JobName}_{YearRes}1231_sechiba_rest.nc'
+    file_SRF_beg = f'{FileDir}/SRF_{JobName}_{TarRestartDate_beg}_sechiba_rest.nc'
     config['Files']['file_SRF_beg'] = file_SRF_beg
 if wu.unDefined ( 'file_SRF_end' ) :
-    file_SRF_end = f'{RunDir}/SRF_{JobName}_{YearEnd}1231_sechiba_rest.nc'
+    file_SRF_end = f'{FileDir}/SRF_{JobName}_{TarRestartDate_end}_sechiba_rest.nc'
     config['Files']['file_SRF_end'] = file_SRF_end
     
@@ -356 +451 @@
 if ICO :
     if wu.unDefined ('file_DYN_aire') : file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ATM+'_grid.nc' )
-    config['Files'] ['file_DYN_aire'] = file_DYN_aire
+    config['Files']['file_DYN_aire'] = file_DYN_aire
 
 if Routing == 'SIMPLE' :
     if wu.unDefined ( 'file_RUN_beg' ) : 
-        file_RUN_beg = f'{RunDir}/SRF_{JobName}_{YearRes}1231_routing_restart.nc'
+        file_RUN_beg = f'{FileDir}/SRF_{JobName}_{TarRestartDate_beg}_routing_restart.nc'
         config['Files']['file_RUN_beg'] = file_RUN_beg
     if wu.unDefined ( 'file_RUN_end' ) : 
-        file_RUN_end = f'{RunDir}/SRF_{JobName}_{YearEnd}1231_routing_restart.nc'
+        file_RUN_end = f'{FileDir}/SRF_{JobName}_{TarRestartDate_end}_routing_restart.nc'
         config['Files']['file_RUN_end'] = file_RUN_end
         
@@ -373 +468 @@
 echo ('\nExtract restart files from tar : ATM, ICO and SRF')
 
-def extract_and_rebuild ( file_name=file_ATM_beg, tar_restart=tar_restart_end, RunDirComp=RunDir, ErrorCount=ErrorCount ) :
+def extract_and_rebuild ( file_name=file_ATM_beg, tar_restart=tar_restart_end, FileDirComp=FileDir, ErrorCount=ErrorCount ) :
     echo ( f'----------')
     echo ( f'file to extract : {file_name = }' )
-    if os.path.exists ( os.path.join (RunDir, file_name) ) :
+    if os.path.exists ( os.path.join (FileDir, file_name) ) :
         echo ( f'file found : {file_name = }' )
     else :
@@ -382 +477 @@
         base_resFile = os.path.basename (file_name)
         if os.path.exists ( tar_restart ) :
-            command =  f'cd {RunDir} ; tar xf {tar_restart} {base_resFile}'
+            command =  f'cd {FileDir} ; tar xf {tar_restart} {base_resFile}'
             echo ( f'{command = }' )
-            try : 
-                os.system ( command )
+            try : os.system ( command )
             except :
                 if ContinueOnError :
@@ -404 +498 @@
     return ErrorCount
  
-ErrorCount += extract_and_rebuild ( file_name=file_ATM_beg, tar_restart=tar_restart_beg_ATM, RunDirComp=RunDir )
-ErrorCount += extract_and_rebuild ( file_name=file_DYN_beg, tar_restart=tar_restart_beg_DYN, RunDirComp=RunDir )
-ErrorCount += extract_and_rebuild ( file_name=file_SRF_beg, tar_restart=tar_restart_beg_SRF, RunDirComp=RunDir )
-
-ErrorCount += extract_and_rebuild ( file_name=file_ATM_end, tar_restart=tar_restart_end_ATM, RunDirComp=RunDir )
-ErrorCount += extract_and_rebuild ( file_name=file_DYN_end, tar_restart=tar_restart_end_DYN, RunDirComp=RunDir )
-ErrorCount += extract_and_rebuild ( file_name=file_SRF_end, tar_restart=tar_restart_end_SRF, RunDirComp=RunDir )
+ErrorCount += extract_and_rebuild ( file_name=file_ATM_beg, tar_restart=tar_restart_beg_ATM, FileDirComp=FileDir )
+ErrorCount += extract_and_rebuild ( file_name=file_DYN_beg, tar_restart=tar_restart_beg_DYN, FileDirComp=FileDir )
+ErrorCount += extract_and_rebuild ( file_name=file_SRF_beg, tar_restart=tar_restart_beg_SRF, FileDirComp=FileDir )
+
+ErrorCount += extract_and_rebuild ( file_name=file_ATM_end, tar_restart=tar_restart_end_ATM, FileDirComp=FileDir )
+ErrorCount += extract_and_rebuild ( file_name=file_DYN_end, tar_restart=tar_restart_end_DYN, FileDirComp=FileDir )
+ErrorCount += extract_and_rebuild ( file_name=file_SRF_end, tar_restart=tar_restart_end_SRF, FileDirComp=FileDir )
 
 if Routing == 'SIMPLE' :
-    ErrorCount += extract_and_rebuild ( file_name=file_RUN_beg, tar_restart=tar_restart_beg_RUN, RunDirComp=RunDir )
-    ErrorCount += extract_and_rebuild ( file_name=file_RUN_end, tar_restart=tar_restart_end_RUN, RunDirComp=RunDir )
+    ErrorCount += extract_and_rebuild ( file_name=file_RUN_beg, tar_restart=tar_restart_beg_RUN, FileDirComp=FileDir )
+    ErrorCount += extract_and_rebuild ( file_name=file_RUN_end, tar_restart=tar_restart_end_RUN, FileDirComp=FileDir )
 
 ##-- Exit in case of error in the opening file phase
@@ -423 +517 @@
 ## 
 echo ('\nOpening ATM SRF and ICO restart files')
-d_ATM_beg = xr.open_dataset ( os.path.join (RunDir, file_ATM_beg), decode_times=False, decode_cf=True ).squeeze()
-d_ATM_end = xr.open_dataset ( os.path.join (RunDir, file_ATM_end), decode_times=False, decode_cf=True ).squeeze()
-d_SRF_beg = xr.open_dataset ( os.path.join (RunDir, file_SRF_beg), decode_times=False, decode_cf=True ).squeeze()
-d_SRF_end = xr.open_dataset ( os.path.join (RunDir, file_SRF_end), decode_times=False, decode_cf=True ).squeeze()
-d_DYN_beg = xr.open_dataset ( os.path.join (RunDir, file_DYN_beg), decode_times=False, decode_cf=True ).squeeze()
-d_DYN_end = xr.open_dataset ( os.path.join (RunDir, file_DYN_end), decode_times=False, decode_cf=True ).squeeze()
+d_ATM_beg = xr.open_dataset ( os.path.join (FileDir, file_ATM_beg), decode_times=False, decode_cf=True ).squeeze()
+d_ATM_end = xr.open_dataset ( os.path.join (FileDir, file_ATM_end), decode_times=False, decode_cf=True ).squeeze()
+d_SRF_beg = xr.open_dataset ( os.path.join (FileDir, file_SRF_beg), decode_times=False, decode_cf=True ).squeeze()
+d_SRF_end = xr.open_dataset ( os.path.join (FileDir, file_SRF_end), decode_times=False, decode_cf=True ).squeeze()
+d_DYN_beg = xr.open_dataset ( os.path.join (FileDir, file_DYN_beg), decode_times=False, decode_cf=True ).squeeze()
+d_DYN_end = xr.open_dataset ( os.path.join (FileDir, file_DYN_end), decode_times=False, decode_cf=True ).squeeze()
 
 for var in d_SRF_beg.variables :
@@ -435 +529 @@
 
 if Routing == 'SIMPLE' :
-    d_RUN_beg = xr.open_dataset ( os.path.join (RunDir, file_RUN_beg), decode_times=False, decode_cf=True ).squeeze()
-    d_RUN_end = xr.open_dataset ( os.path.join (RunDir, file_RUN_end), decode_times=False, decode_cf=True ).squeeze()
+    d_RUN_beg = xr.open_dataset ( os.path.join (FileDir, file_RUN_beg), decode_times=False, decode_cf=True ).squeeze()
+    d_RUN_end = xr.open_dataset ( os.path.join (FileDir, file_RUN_end), decode_times=False, decode_cf=True ).squeeze()
 
 def to_cell ( dd, newname='cell' ) :
@@ -470 +564 @@
 
 ## Write the full configuration
-#config_out = open (FullIniFile, 'w')
-#config.write ( config_out )
-#config_out.close ()
+config_out = open (FullIniFile, 'w')
+config.write ( config_out )
+config_out.close ()
 
 if ICO :
@@ -554 +648 @@
         SRF_areas    = xr.DataArray ( SRF_areas   , coords=d_SRF_his ['Contfrac'].coords, dims=d_SRF_his ['Contfrac'].dims)
         SRF_areafrac = xr.DataArray ( SRF_areafrac, coords=d_SRF_his ['Contfrac'].coords, dims=d_SRF_his ['Contfrac'].dims)
-        SRF_contfrac = SRF_areafrac / SRF_areas
-        #SRF_aire      = lmdz.geo2point (SRF_areas * SRF_contfrac, dim1D='cell_latlon', cumulPoles=True )
-        #SRF_aire      = lmdz.geo2point (SRF_areas, dim1D='cell_latlon', cumulPoles=True )
-        #SRF_aire      = lmdz.geo2point (d_ATM_aire ['aire'].squeeze().values * SRF_contfrac, dim1D='cell_latlon', cumulPoles=True )
+
         SRF_areas     = lmdz.geo2point ( SRF_areas   , dim1D='cell_latlon', cumulPoles=True )
         SRF_areafrac  = lmdz.geo2point ( SRF_areafrac, dim1D='cell_latlon', cumulPoles=True )
-        SRF_contfrac  = lmdz.geo2point ( SRF_contfrac, dim1D='cell_latlon' )
+        SRF_contfrac  = SRF_areafrac / SRF_areas
         SRF_aire      = SRF_areafrac
 
@@ -574 +665 @@
         
     if SRF_HIS == 'ico' : 
-        SRF_lat       =  d_SRF_his ['lat'][SRF_his_keysort]
-        SRF_lon       =  d_SRF_his ['lon'][SRF_his_keysort]
-        SRF_aire      =  d_SRF_his ['Areas'][SRF_his_keysort] * d_SRF_his ['Contfrac'][SRF_his_keysort]
-        SRF_areas     =  d_SRF_his ['Areas'][SRF_his_keysort] 
+        SRF_lat       =  d_SRF_his ['lat']      [SRF_his_keysort]
+        SRF_lon       =  d_SRF_his ['lon']      [SRF_his_keysort]
+        SRF_areas     =  d_SRF_his ['Areas']    [SRF_his_keysort] 
+        #SRF_areafrac =  d_SRF_his ['AreaFrac'] [SRF_his_keysort] 
         SRF_contfrac  =  d_SRF_his ['Contfrac'][SRF_his_keysort]
+        SRF_aire      =  SRF_areas * SRF_contfrac
 
 ATM_fsea      = ATM_foce + ATM_fsic
@@ -589 +681 @@
 ATM_aire_fsea = ATM_aire * ATM_fsea
 
-SRF_aire = SRF_aire.where ( np.abs (SRF_aire) < 1E15, 0.)
+SRF_aire = SRF_aire.where ( np.abs (SRF_aire) < 1E15, 0. )
 
 ## Write the full configuration
@@ -1023 +1115 @@
 echo ( f'-- ATM Fluxes  -- {Title} ' )
 
-print ( 'Step 1', end='' )
+Step=1
+print ( 'Step {Step}', end='' ) ; Step += 1
 if ATM_HIS == 'latlon' :
     ATM_wbilo_oce   = lmdz.geo2point ( d_ATM_his ['wbilo_oce'], dim1D='cell_latlon' )
@@ -1039 +1132 @@
     ATM_wevap_lic   = lmdz.geo2point ( d_ATM_his ['wevap_lic'], dim1D='cell_latlon' )
     ATM_wevap_sic   = lmdz.geo2point ( d_ATM_his ['wevap_sic'], dim1D='cell_latlon' )
+
+    ATM_wrain_ter   = lmdz.geo2point ( d_ATM_his ['wrain_ter'], dim1D='cell_latlon' )
+    ATM_wrain_oce   = lmdz.geo2point ( d_ATM_his ['wrain_oce'], dim1D='cell_latlon' )
+    ATM_wrain_lic   = lmdz.geo2point ( d_ATM_his ['wrain_lic'], dim1D='cell_latlon' )
+    ATM_wrain_sic   = lmdz.geo2point ( d_ATM_his ['wrain_sic'], dim1D='cell_latlon' )
+
+    ATM_wsnow_ter   = lmdz.geo2point ( d_ATM_his ['wsnow_ter'], dim1D='cell_latlon' )
+    ATM_wsnow_oce   = lmdz.geo2point ( d_ATM_his ['wsnow_oce'], dim1D='cell_latlon' )
+    ATM_wsnow_lic   = lmdz.geo2point ( d_ATM_his ['wsnow_lic'], dim1D='cell_latlon' )
+    ATM_wsnow_sic   = lmdz.geo2point ( d_ATM_his ['wsnow_sic'], dim1D='cell_latlon' )
+
+    ATM_wevap_ter   = lmdz.geo2point ( d_ATM_his ['wevap_ter'], dim1D='cell_latlon' )
+    ATM_wevap_oce   = lmdz.geo2point ( d_ATM_his ['wevap_oce'], dim1D='cell_latlon' )
+    ATM_wevap_lic   = lmdz.geo2point ( d_ATM_his ['wevap_lic'], dim1D='cell_latlon' )
+    ATM_wevap_sic   = lmdz.geo2point ( d_ATM_his ['wevap_sic'], dim1D='cell_latlon' )
+    
     ATM_runofflic   = lmdz.geo2point ( d_ATM_his ['runofflic'], dim1D='cell_latlon' )
     
-print ( ' - Step 2', end='' )
+print ( 'Step {Step}', end='' ) ; Step += 1
 if ATM_HIS == 'ico' :
     ATM_wbilo_oce   = d_ATM_his ['wbilo_oce'][..., ATM_his_keysort]
@@ -1059 +1168 @@
     ATM_runofflic   = d_ATM_his ['runofflic'][..., ATM_his_keysort]
 
-print ( ' - Step 3', end='' )
+    ATM_wevap_ter   = d_ATM_his ['wevap_ter'][..., ATM_his_keysort]
+    ATM_wevap_oce   = d_ATM_his ['wevap_oce'][..., ATM_his_keysort]
+    ATM_wevap_lic   = d_ATM_his ['wevap_lic'][..., ATM_his_keysort]
+    ATM_wevap_sic   = d_ATM_his ['wevap_sic'][..., ATM_his_keysort]
+
+    ATM_wrain_ter   = d_ATM_his ['wrain_ter'][..., ATM_his_keysort]
+    ATM_wrain_oce   = d_ATM_his ['wrain_oce'][..., ATM_his_keysort]
+    ATM_wrain_lic   = d_ATM_his ['wrain_lic'][..., ATM_his_keysort]
+    ATM_wrain_sic   = d_ATM_his ['wrain_sic'][..., ATM_his_keysort]
+
+    ATM_wsnow_ter   = d_ATM_his ['wsnow_ter'][..., ATM_his_keysort]
+    ATM_wsnow_oce   = d_ATM_his ['wsnow_oce'][..., ATM_his_keysort]
+    ATM_wsnow_lic   = d_ATM_his ['wsnow_lic'][..., ATM_his_keysort]
+    ATM_wsnow_sic   = d_ATM_his ['wsnow_sic'][..., ATM_his_keysort]
+
+print ( 'Step {Step}', end='' ) ; Step += 1
+
+ATM_wprecip_oce = ATM_wrain_oce + ATM_wsnow_oce
+ATM_wprecip_ter = ATM_wrain_ter + ATM_wsnow_ter
+ATM_wprecip_sic = ATM_wrain_sic + ATM_wsnow_sic
+ATM_wprecip_lic = ATM_wrain_lic + ATM_wsnow_lic
+
 ATM_wbilo       = ATM_wbilo_oce + ATM_wbilo_sic + ATM_wbilo_ter + ATM_wbilo_lic
 ATM_emp         = ATM_evap - ATM_precip
-ATM_wbilo_sea   = ATM_wbilo_oce + ATM_wbilo_sic
-ATM_wevap_sea   = ATM_wevap_sic + ATM_wevap_oce
-
-ATM_wemp_ter = ATM_wevap_ter - ATM_precip * ATM_fter
-ATM_wemp_oce = ATM_wevap_oce - ATM_precip * ATM_foce
-ATM_wemp_sic = ATM_wevap_sic - ATM_precip * ATM_fsic
-ATM_wemp_lic = ATM_wevap_lic - ATM_precip * ATM_flic
-ATM_wemp_sea = ATM_wevap_sic + ATM_wevap_oce - ATM_precip * ATM_fsea
-
-print ( ' - Step 4', end='' )
+ATM_wprecip_sea = ATM_wprecip_oce + ATM_wprecip_sic
+ATM_wsnow_sea   = ATM_wsnow_oce   + ATM_wsnow_sic
+ATM_wrain_sea   = ATM_wrain_oce   + ATM_wrain_sic
+ATM_wbilo_sea   = ATM_wbilo_oce   + ATM_wbilo_sic
+ATM_wevap_sea   = ATM_wevap_sic   + ATM_wevap_oce
+
+#ATM_wemp_ter = ATM_wevap_ter - ATM_precip * ATM_fter
+#ATM_wemp_oce = ATM_wevap_oce - ATM_precip * ATM_foce
+#ATM_wemp_sic = ATM_wevap_sic - ATM_precip * ATM_fsic
+#ATM_wemp_lic = ATM_wevap_lic - ATM_precip * ATM_flic
+#ATM_wemp_sea = ATM_wevap_sic + ATM_wevap_oce - ATM_precip * ATM_fsea
+
+ATM_wemp_ter = ATM_wevap_ter - ATM_wprecip_ter 
+ATM_wemp_oce = ATM_wevap_oce - ATM_wprecip_oce
+ATM_wemp_sic = ATM_wevap_sic - ATM_wprecip_sic
+ATM_wemp_lic = ATM_wevap_lic - ATM_wprecip_lic
+ATM_wemp_sea = ATM_wevap_sic - ATM_wprecip_oce
+
+print ( 'Step {Step}', end='' ) ; Step += 1
 if RUN_HIS == 'latlon' : 
     RUN_coastalflow = lmdz.geo2point ( d_RUN_his ['coastalflow'], dim1D='cell_latlon' ) 
@@ -1082 +1221 @@
     RUN_riverflow_cpl   = lmdz.geo2point ( d_RUN_his ['riverflow_cpl']  , dim1D='cell_latlon' )
 
-print ( ' - Step 5', end='' )
+print ( 'Step {Step}', end='' ) ; Step += 1
 if RUN_HIS == 'ico' :
     RUN_coastalflow =  d_RUN_his ['coastalflow'][..., RUN_his_keysort]
@@ -1093 +1232 @@
     RUN_riverflow_cpl   = d_RUN_his ['riverflow_cpl']  [..., RUN_his_keysort]
 
-print ( ' - Step 6', end='' )
+print ( 'Step {Step}', end='' ) ; Step += 1
 if SRF_HIS == 'latlon' :
     SRF_rain     = lmdz.geo2point ( d_SRF_his ['rain']  , dim1D='cell_latlon')
@@ -1101 +1240 @@
     SRF_transpir = lmdz.geo2point ( np.sum (d_SRF_his ['transpir'], axis=1), dim1D='cell_latlon' )
 
-print ( '- Step 7', end='' )
+print ( 'Step {Step}', end='' ) ; Step += 1
 if SRF_HIS == 'ico' :
     SRF_rain     = d_SRF_his ['rain'] [..., SRF_his_keysort]
@@ -1109 +1248 @@
     SRF_transpir = np.sum (d_SRF_his ['transpir'], axis=1)[..., SRF_his_keysort]
 
-print ( '- Step 8', end='' )
+print ( 'Step {Step}', end='' ) ; Step += 1
 SRF_transpir.attrs['units'] = d_SRF_his ['transpir'].attrs['units']
 SRF_emp      = SRF_evap - SRF_rain - SRF_snowf ; SRF_emp.attrs['units'] = SRF_rain.attrs['units']
@@ -1132 +1271 @@
     mmd2SI (VarT)
 
-print ( '- Step 9', end='' )
+print ( 'Step {Step}', end='' ) ; Step += 1
 RUN_input  = RUN_runoff      + RUN_drainage
 RUN_output = RUN_coastalflow + RUN_riverflow
 
-print ( '- Step 10', end='' )
+print ( 'Step {Step}', end='' ) ; Step += 1
 ATM_flx_wbilo       = ATM_flux_int ( ATM_wbilo      )
+ATM_flx_wbilo_lic   = ATM_flux_int ( ATM_wbilo_lic  )
 ATM_flx_wbilo_oce   = ATM_flux_int ( ATM_wbilo_oce  )
+ATM_flx_wbilo_sea   = ATM_flux_int ( ATM_wbilo_sea  )
 ATM_flx_wbilo_sic   = ATM_flux_int ( ATM_wbilo_sic  )
-ATM_flx_wbilo_sea   = ATM_flux_int ( ATM_wbilo_sea  )
 ATM_flx_wbilo_ter   = ATM_flux_int ( ATM_wbilo_ter  )
-ATM_flx_wbilo_lic   = ATM_flux_int ( ATM_wbilo_lic  )
 ATM_flx_calving     = ATM_flux_int ( ATM_fqcalving  )
 ATM_flx_fqfonte     = ATM_flux_int ( ATM_fqfonte    )
 
-print ( '- Step 11', end='' )
+print ( 'Step {Step}', end='' ) ; Step += 1
 ATM_flx_precip      = ATM_flux_int ( ATM_precip     )
 ATM_flx_snowf       = ATM_flux_int ( ATM_snowf      )
@@ -1152 +1291 @@
 ATM_flx_runlic      = ATM_flux_int ( ATM_runofflic  )
 
-print ( '- Step 12', end='' )
-ATM_flx_wevap_ter   = ATM_flux_int ( ATM_wevap_ter )
-ATM_flx_wevap_sea   = ATM_flux_int ( ATM_wevap_sea )
-ATM_flx_precip_ter  = ATM_flux_int ( ATM_precip * ATM_fter )
-ATM_flx_precip_sea  = ATM_flux_int ( ATM_precip * ATM_fsea )
-ATM_flx_emp         = ATM_flux_int ( ATM_emp)
-ATM_flx_wemp_ter    = ATM_flux_int ( ATM_wemp_ter )
-ATM_flx_wemp_oce    = ATM_flux_int ( ATM_wemp_oce )
-ATM_flx_wemp_lic    = ATM_flux_int ( ATM_wemp_lic )
-ATM_flx_wemp_sea    = ATM_flux_int ( ATM_wemp_sea )
-
-print ( '- Step 13', end='' )
+print ( 'Step {Step}', end='' ) ; Step += 1
+ATM_flx_wrain_ter    = ATM_flux_int ( ATM_wrain_ter )
+ATM_flx_wrain_oce    = ATM_flux_int ( ATM_wrain_oce )
+ATM_flx_wrain_lic    = ATM_flux_int ( ATM_wrain_lic )
+ATM_flx_wrain_sic    = ATM_flux_int ( ATM_wrain_sic )
+ATM_flx_wrain_sea    = ATM_flux_int ( ATM_wrain_sea )
+
+ATM_flx_wsnow_ter    = ATM_flux_int ( ATM_wsnow_ter )
+ATM_flx_wsnow_oce    = ATM_flux_int ( ATM_wsnow_oce )
+ATM_flx_wsnow_lic    = ATM_flux_int ( ATM_wsnow_lic )
+ATM_flx_wsnow_sic    = ATM_flux_int ( ATM_wsnow_sic )
+ATM_flx_wsnow_sea    = ATM_flux_int ( ATM_wsnow_sea )
+print ( 'Step {Step}', end='' ) ; Step += 1
+ATM_flx_wrain_ter    = ATM_flux_int ( ATM_wrain_ter )
+ATM_flx_wrain_oce    = ATM_flux_int ( ATM_wrain_oce )
+ATM_flx_wrain_lic    = ATM_flux_int ( ATM_wrain_lic )
+ATM_flx_wrain_sic    = ATM_flux_int ( ATM_wrain_sic )
+ATM_flx_wrain_sea    = ATM_flux_int ( ATM_wrain_sea )
+ATM_flx_wevap_ter    = ATM_flux_int ( ATM_wevap_ter )
+ATM_flx_wevap_oce    = ATM_flux_int ( ATM_wevap_oce )
+ATM_flx_wevap_lic    = ATM_flux_int ( ATM_wevap_lic )
+ATM_flx_wevap_sic    = ATM_flux_int ( ATM_wevap_sic )
+ATM_flx_wevap_sea    = ATM_flux_int ( ATM_wevap_sea )
+ATM_flx_wprecip_lic  = ATM_flux_int ( ATM_wprecip_lic )
+ATM_flx_wprecip_oce  = ATM_flux_int ( ATM_wprecip_oce )
+ATM_flx_wprecip_sic  = ATM_flux_int ( ATM_wprecip_sic )
+ATM_flx_wprecip_ter  = ATM_flux_int ( ATM_wprecip_ter )
+ATM_flx_wprecip_sea  = ATM_flux_int ( ATM_wprecip_sea )
+ATM_flx_wemp_lic     = ATM_flux_int ( ATM_wemp_lic )
+ATM_flx_wemp_oce     = ATM_flux_int ( ATM_wemp_oce )
+ATM_flx_wemp_sic     = ATM_flux_int ( ATM_wemp_sic )
+ATM_flx_wemp_ter     = ATM_flux_int ( ATM_wemp_ter )
+ATM_flx_wemp_sea     = ATM_flux_int ( ATM_wemp_sea )
+
+ATM_flx_emp          = ATM_flux_int ( ATM_emp)
+
+print ( 'Step {Step}', end='' ) ; Step += 1
 RUN_flx_coastal     = ONE_flux_int ( RUN_coastalflow)
 RUN_flx_river       = ONE_flux_int ( RUN_riverflow  )
@@ -1174 +1338 @@
 RUN_flx_output      = ONE_flux_int ( RUN_output     )
 
-print ( '- Step 14' )
+print ( 'Step {Step}', end='' ) ; Step += 1
 RUN_flx_bil    = RUN_flx_input   - RUN_flx_output
 RUN_flx_rivcoa = RUN_flx_coastal + RUN_flx_river
@@ -1190 +1354 @@
 prtFlux ('snowf         ', ATM_flx_snowf      , 'f' )        
 prtFlux ('evap          ', ATM_flx_evap       , 'f' )
-prtFlux ('runoff lic    ', ATM_flx_runlic     , 'f' )   
+prtFlux ('runoff lic    ', ATM_flx_runlic     , 'f' )
+
+prtFlux ('ATM_flx_wemp_lic     ', ATM_flx_wemp_lic      , 'f' )
+prtFlux ('ATM_flx_wemp_oce     ', ATM_flx_wemp_oce      , 'f' )
+prtFlux ('ATM_flx_wemp_sic     ', ATM_flx_wemp_sic      , 'f' )
+prtFlux ('ATM_flx_wemp_ter     ', ATM_flx_wemp_ter      , 'f' )
+
+prtFlux ('ATM_flx_wprecip_lic  ', ATM_flx_wprecip_lic   , 'f' )
+prtFlux ('ATM_flx_wprecip_oce  ', ATM_flx_wprecip_oce   , 'f' )
+prtFlux ('ATM_flx_wprecip_sic  ', ATM_flx_wprecip_sic   , 'f' )
+prtFlux ('ATM_flx_wprecip_ter  ', ATM_flx_wprecip_ter   , 'f' )
+
+prtFlux ('ATM_flx_wrain_lic    ', ATM_flx_wrain_lic     , 'f' )
+prtFlux ('ATM_flx_wrain_oce    ', ATM_flx_wrain_oce     , 'f' )
+prtFlux ('ATM_flx_wrain_sic    ', ATM_flx_wrain_sic     , 'f' )
+prtFlux ('ATM_flx_wrain_ter    ', ATM_flx_wrain_ter     , 'f' )
+
+prtFlux ('ATM_flx_wsnow_lic    ', ATM_flx_wsnow_lic     , 'f' )
+prtFlux ('ATM_flx_wsnow_oce    ', ATM_flx_wsnow_oce     , 'f' )
+prtFlux ('ATM_flx_wsnow_sic    ', ATM_flx_wsnow_sic     , 'f' )
+prtFlux ('ATM_flx_wsnow_ter    ', ATM_flx_wsnow_ter     , 'f' )
+
+prtFlux ('ATM_flx_wevap_lic    ', ATM_flx_wevap_lic     , 'f' )
+prtFlux ('ATM_flx_wevap_oce    ', ATM_flx_wevap_oce     , 'f' )
+prtFlux ('ATM_flx_wevap_sic    ', ATM_flx_wevap_sic     , 'f' )
+prtFlux ('ATM_flx_wevap_ter    ', ATM_flx_wevap_ter     , 'f' )
+
+prtFlux ('ATM_flx_wevap_lic    ', ATM_flx_wevap_lic     , 'f' )
+prtFlux ('ATM_flx_wevap_oce    ', ATM_flx_wevap_oce     , 'f' )
+prtFlux ('ATM_flx_wevap_sic    ', ATM_flx_wevap_sic     , 'f' )
+prtFlux ('ATM_flx_wevap_ter    ', ATM_flx_wevap_ter     , 'f' )
 
 echo ( '\n====================================================================================' )

TOOLS/WATER_BUDGET/OCE_waterbudget.py

@@ -22 +22 @@
 ###
 ## Import system modules
-import sys, os, shutil, subprocess, platform
-import numpy as np
+import sys, os, shutil#, subprocess, platform
 import configparser, re
 from pathlib import Path
 
-## Import local module
-import WaterUtils as wu
+## Import needed scientific modules
+import numpy as np, xarray as xr
 
 # Check python version
@@ -35 +34 @@
     raise Exception ( "Minimum Python version is 3.8" )
 
+## Import local module
+import WaterUtils as wu
+import libIGCM_sys
+import nemo, lmdz
+
 ## Creates parser for reading .ini input file
-config = configparser.ConfigParser()
+## -------------------------------------------
+config = configparser.ConfigParser ( interpolation=configparser.ExtendedInterpolation() )
 config.optionxform = str # To keep capitals
 
-## Read experiment parameters
-ATM=None ; ATM_HIS='latlon' ; SRF_HIS='latlon' ; RUN_HIS='latlon' ; ORCA=None ; NEMO=None ; OCE_relax=False ; OCE_icb=False ; Coupled=False ; Routing=None ; TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=None
+## Experiment parameters
+## ---------------------
+ATM=None ; ATM_HIS='latlon' ; SRF_HIS='latlon' ; RUN_HIS='latlon' ; ORCA=None ; NEMO=None ; OCE_relax=False
+OCE_icb=False ; Coupled=False ; Routing=None
+TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=None
+
 ##
-ARCHIVE=None ; STORAGE=None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None
-TmpDir=None ; RunDir=None ; FileOut=None ; dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None ; FileCommon=None
-file_ATM_his=None ; file_SRF_his=None ; file_RUN_his=None ; file_OCE_his=None ;  file_ICE_his=None ; file_OCE_sca=None
-tar_restart_beg=None ; tar_restart_end=None ; file_ATM_beg=None ; file_ATM_end=None ; file_DYN_beg=None ; file_DYN_end=None ; file_SRF_beg=None ; file_SRF_end=None
-file_RUN_beg=None ; file_RUN_end=None ; file_RUN_end=None ; file_OCE_beg=None ; file_ICE_beg=None ; file_OCE_beg=None ; file_OCE_end=None ; file_ICE_beg=None ; file_OCE_end=None ; file_ICE_end=None ; ContinueOnError=False ; ErrorCount=0
-tar_restart_beg_ATM=None ; tar_restart_beg_DYN=None ; tar_restart_beg_SRF=None ; tar_restart_beg_RUN=None ; tar_restart_beg_OCE=None ; tar_restart_beg_ICE=None
-tar_restart_end_ATM=None ; tar_restart_end_DYN=None ; tar_restart_end_SRF=None ; tar_restart_end_RUN=None ; tar_restart_end_OCE=None ; tar_restart_end_ICE=None
-
-# Read command line arguments
+ARCHIVE=None ; STORAGE=None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None ; TmpDir=None
+FileDir=None ; FileOut=None
+dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None
+FileCommon=None ; file_ATM_his=None ; file_SRF_his=None ; file_RUN_his=None
+file_OCE_his=None ;  file_ICE_his=None ; file_OCE_sca=None
+tar_restart_beg=None ; tar_restart_end=None ; file_ATM_beg=None ; file_ATM_end=None ; file_DYN_beg=None
+file_DYN_end=None ; file_SRF_beg=None ; file_SRF_end=None
+file_RUN_beg=None ; file_RUN_end=None ; file_RUN_end=None ; file_OCE_beg=None
+file_ICE_beg=None ; file_OCE_beg=None
+file_OCE_end=None ; file_ICE_beg=None ; file_OCE_end=None ; file_ICE_end=None
+tar_restart_beg_ATM=None ; tar_restart_beg_DYN=None ; tar_restart_beg_SRF=None
+tar_restart_beg_RUN=None ; tar_restart_beg_OCE=None ; tar_restart_beg_ICE=None
+tar_restart_end_ATM=None ; tar_restart_end_DYN=None ; tar_restart_end_SRF=None
+tar_restart_end_RUN=None ; tar_restart_end_OCE=None ; tar_restart_end_ICE=None
+ContinueOnError=False ; ErrorCount=0 ; SortIco = False
+
+## Read command line arguments
+## ---------------------------
 print ( "Name of Python script:", sys.argv[0] )
 IniFile = sys.argv[1]
-# Text existence of IniFile
+
+# Test existence of IniFile
 if not os.path.exists (IniFile ) :
     raise FileExistsError ( f"File not found : {IniFile = }" )
@@ -63 +82 @@
 config.read (IniFile)
 FullIniFile = 'full_' + IniFile
-    
+
+# Reading config.card if possible
+## -------------------------------
+ConfigCard = None
+
+if 'Experiment' in config.keys ()  : ## Read Experiment on Config file if possible
+    if 'ConfigCard' in config['Experiment'].keys () :
+        ConfigCard = config['Experiment']['ConfigCard']
+        print ( f'{ConfigCard=}' )
+
+if ConfigCard : ## Read config card if it exists
+    # Text existence of ConfigCard
+    if os.path.exists ( ConfigCard ) :
+        print ( f'Reading Config Card : {ConfigCard}' )
+        ## Creates parser for reading .ini input file
+        MyReader = configparser.ConfigParser (interpolation=configparser.ExtendedInterpolation() )
+        MyReader.optionxform = str # To keep capitals
+        
+        MyReader.read (ConfigCard)
+
+        for VarName in ['JobName', 'ExperimentName', 'SpaceName', 'LongName', 'ModelName', 'TagName'] :
+            if VarName in MyReader['UserChoices'].keys() :
+                locals()[VarName] = MyReader['UserChoices'][VarName]
+                exec  ( f'{VarName} = wu.setBool ({VarName})' )
+                exec  ( f'{VarName} = wu.setNum  ({VarName})' )
+                exec  ( f'{VarName} = wu.setNone ({VarName})' )
+                exec  ( f'wu.{VarName} = {VarName}' )
+                print ( f'    {VarName:21} set to : {locals()[VarName]:}' )
+                
+        for VarName in ['PackFrequency'] :
+            if VarName in MyReader['Post'].keys() :
+                locals()[VarName] = MyReader['Post'][VarName]
+                exec  ( f'{VarName} = wu.setBool ({VarName})' )
+                exec  ( f'{VarName} = wu.setNum  ({VarName})' )
+                exec  ( f'{VarName} = wu.setNone ({VarName})' )
+                exec  ( f'wu.{VarName} = {VarName}' )
+                print ( f'    {VarName:21} set to : {locals()[VarName]:}' )
+    else :
+        raise FileExistsError ( f"File not found : {ConfigCard = }" )
+
 ## Reading config file
-for Section in ['Experiment', 'libIGCM', 'Files', 'Physics' ] :
+## -------------------
+for Section in ['Config', 'Experiment', 'libIGCM', 'Files', 'Physics' ] :
    if Section in config.keys () : 
-        print ( f'Reading [{Section}]' )
+        print ( f'\nReading [{Section}]' )
         for VarName in config[Section].keys() :
             locals()[VarName] = config[Section][VarName]
@@ -74 +133 @@
             exec ( f'{VarName} = wu.setNone ({VarName})' )
             exec ( f'wu.{VarName} = {VarName}' )
-            print ( f'{VarName:25} set to : {locals()[VarName])}' )
+            print ( f'    {VarName:21} set to : {locals()[VarName]}' )
             #exec ( f'del {VarName}' )
 
 print ( f'\nConfig file readed : {IniFile} ' )
             
-#-- Some physical constants
-if wu.unDefined ( 'Ra' )           : Ra          = 6366197.7236758135 #-- Earth Radius (m)
-if wu.unDefined ( 'Grav' )         : Grav        = 9.81               #-- Gravity (m^2/s
-if wu.unDefined ( 'ICE_rho_ice' )  : ICE_rho_ice = 917.0              #-- Ice volumic mass (kg/m3) in LIM3
-if wu.unDefined ( 'ICE_rho_sno')   : ICE_rho_sno = 330.0              #-- Snow volumic mass (kg/m3) in LIM3
-if wu.unDefined ( 'OCE_rho_liq' )  : OCE_rho_liq = 1026.0             #-- Ocean water volumic mass (kg/m3) in NEMO
-if wu.unDefined ( 'ATM_rho' )      : ATM_rho     = 1000.0             #-- Water volumic mass in atmosphere (kg/m^3)
-if wu.unDefined ( 'SRF_rho' )      : SRF_rho     = 1000.0             #-- Water volumic mass in surface reservoir (kg/m^3)
-if wu.unDefined ( 'RUN_rho' )      : RUN_rho     = 1000.0             #-- Water volumic mass of rivers (kg/m^3)
-if wu.unDefined ( 'ICE_rho_pnd' )  : ICE_rho_pnd = 1000.              #-- Water volumic mass in ice ponds (kg/m^3)
-if wu.unDefined ( 'YearLength' )   : YearLength  = 365.25 * 24. * 60. * 60. #-- Year length (s)
-
-# Set libIGCM and machine dependant values
-if not 'Files' in config.keys() : config['Files'] = {}
-    
-config['Physics'] = { 'Ra':Ra, 'Grav':Grav, 'ICE_rho_ice':ICE_rho_ice, 'ICE_rho_sno':ICE_rho_sno, 'OCE_rho_liq':OCE_rho_liq, 'ATM_rho':ATM_rho, 'SRF_rho':SRF_rho, 'RUN_rho':RUN_rho}
-
-config['Config'] = { 'ContinueOnError':ContinueOnError, 'SortIco':SortIco}
-
-with open ('SetLibIGCM.py') as f: exec ( f.read() )
+## Some physical constants
+## =======================
+if wu.unDefined ( 'Ra' )           : Ra          = wu.Ra           #-- Earth Radius (m)
+if wu.unDefined ( 'Grav' )         : Grav        = wu.Grav         #-- Gravity (m^2/s
+if wu.unDefined ( 'ICE_rho_ice' )  : ICE_rho_ice = wu.ICE_rho_ice  #-- Ice volumic mass (kg/m3) in LIM3
+if wu.unDefined ( 'ICE_rho_sno')   : ICE_rho_sno = wu.ICE_rho_sno  #-- Snow volumic mass (kg/m3) in LIM3
+if wu.unDefined ( 'OCE_rho_liq' )  : OCE_rho_liq = wu.OCE_rho_liq  #-- Ocean water volumic mass (kg/m3) in NEMO
+if wu.unDefined ( 'ATM_rho' )      : ATM_rho     = wu.ATM_rho      #-- Water volumic mass in atmosphere (kg/m^3)
+if wu.unDefined ( 'SRF_rho' )      : SRF_rho     = wu.SRF_rho      #-- Water volumic mass in surface reservoir (kg/m^3)
+if wu.unDefined ( 'RUN_rho' )      : RUN_rho     = wu.RUN_rho      #-- Water volumic mass of rivers (kg/m^3)
+if wu.unDefined ( 'ICE_rho_pnd' )  : ICE_rho_pnd = wu.ICE_rho_pnd  #-- Water volumic mass in ice ponds (kg/m^3)
+if wu.unDefined ( 'YearLength' )   : YearLength  = wu.YearLength   #-- Year length (s)
+
+## Set libIGCM and machine dependant values
+## ----------------------------------------
+if not 'Files' in config.keys () : config['Files'] = {}
+
+config['Physics'] = { 'Ra':Ra, 'Grav':Grav, 'ICE_rho_ice':ICE_rho_ice, 'ICE_rho_sno':ICE_rho_sno,
+                      'OCE_rho_liq':OCE_rho_liq, 'ATM_rho':ATM_rho, 'SRF_rho':SRF_rho, 'RUN_rho':RUN_rho}
+
+config['Config'] = { 'ContinueOnError':ContinueOnError, 'SortIco':SortIco }
+
+## --------------------------
+ICO  = ( 'ICO' in wu.ATM )
+LMDZ = ( 'LMD' in wu.ATM )
+
+mm = libIGCM_sys.config ( TagName=TagName, SpaceName=SpaceName, ExperimentName=ExperimentName, JobName=JobName, User=User, Group=Group,
+             ARCHIVE=None, SCRATCHDIR=None, STORAGE=None, R_IN=None, R_OUT=None, R_FIG=None, rebuild=None, TmpDir=None, 
+             R_SAVE=None, R_FIGR=None, R_BUFR=None, R_BUF_KSH=None, REBUILD_DIR=None, POST_DIR=None )
+globals().update(mm)
+
 config['Files']['TmpDir'] = TmpDir
-
-if libIGCM :
-    config['libIGCM'] = {'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'SCRATCHDIR':SCRATCHDIR, 'R_IN':R_IN, 'rebuild':rebuild }
-
-## Import specific module
-import nemo
-## Now import needed scientific modules
-import xarray as xr
-    
-# Output file with water budget diagnostics
-if wu.unDefined ( 'FileOut' ) : FileOut = f'OCE_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out'
+config['libIGCM'] = { 'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'TmpDir':TmpDir, 'R_IN':R_IN, 'rebuild':rebuild }
+
+## Defines begining and end of experiment
+## --------------------------------------
+if wu.unDefined ( 'DateBegin' ) :
+    DateBegin = f'{YearBegin}0101'
+    config['Experiment']['DateBegin'] = DateBegin
+else :
+    YearBegin, MonthBegin, DayBegin = wu.DateSplit ( DateBegin )
+    DateBegin = wu.FormatToGregorian ( DateBegin)
+    config['Experiment']['YearBegin'] = YearBegin
+
+if wu.unDefined ( 'DateEnd' )   :
+    DateEnd   = f'{YearEnd}1231'
+    config['Experiment']['DateEnd']   = DateEnd
+else :
+    YearEnd, MonthEnd, DayEnd = wu.DateSplit ( DateEnd )
+    DateEnd   = wu.FormatToGregorian ( DateEnd)
+
+if wu.unDefined ( 'PackFrequency' ) :
+    PackFrequency = YearEnd - YearBegin +1
+    config['Experiment']['PackFrequency']   = f'{PackFrequency}'
+
+if type ( PackFrequency ) == str : 
+    if 'Y' in PackFrequency : PackFrequency = PackFrequency.replace ( 'Y', '') 
+    if 'M' in PackFrequency : PackFrequency = PackFrequency.replace ( 'M', '')
+    PackFrequency = int ( PackFrequency )
+
+print ( f'{YearBegin=} {DateBegin=}' )
+print ( f'{YearEnd  =} {DateEnd  =}' )
+print ( f'{PackFrequency=}' )
+    
+## Output file with water budget diagnostics
+## -----------------------------------------
+if wu.unDefined ( 'FileOut' ) : FileOut = f'ATM_waterbudget_{JobName}_{DateBegin}_{DateEnd}.out'
 config['Files']['FileOut'] = FileOut
 
 f_out = open ( FileOut, mode = 'w' )
 
-# Useful functions
-def kg2Sv    (val, rho=ATM_rho) :
+## Useful functions
+## ----------------
+def kg2Sv  (val, rho=OCE_rho_liq) :
     '''From kg to Sverdrup'''
     return val/dtime_sec*1.0e-6/rho
 
-def kg2myear (val, rho=ATM_rho) :
+def kg2myear (val, rho=OCE_rho_liq) :
     '''From kg to m/year'''
-    return val/ATM_aire_sea_tot/rho/NbYear
-
-def var2prt (var, small=False, rho=ATM_rho) :
+    return val/OCE_aire_tot/rho/NbYear
+
+def var2prt (var, small=False, rho=OCE_rho_liq) :
     if small :  return var , kg2Sv(var, rho=rho)*1000., kg2myear(var, rho=rho)*1000.
     else     :  return var , kg2Sv(var, rho=rho)      , kg2myear(var, rho=rho)
 
-def prtFlux (Desc, var, Form='F', small=False, rho=ATM_rho, width=15) :
+def prtFlux (Desc, var, Form='F', small=False, rho=OCE_rho_liq, width=15) :
     if small :
-        if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} mSv | {:12.4f} mm/year "
-        if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} mSv | {:12.4e} mm/year "
+        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} mSv | {:12.4f} mm/year "
+        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} mSv | {:12.4e} mm/year "
     else :
-        if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} Sv  | {:12.4f} m/year  "
-        if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} Sv  | {:12.4e} m/year  "
+        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} Sv  | {:12.4f} m/year  "
+        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} Sv  | {:12.4e} m/year  "
     echo ( (' {:>{width}} = ' +ff).format (Desc, *var2prt(var, small, rho=rho), width=width ) )
     return None
@@ -146 +241 @@
     return None
     
-### Set libIGCM directories
+   
+## Set libIGCM directories
+## -----------------------
 if wu.unDefined ('R_OUT'      ) : R_OUT       = os.path.join ( ARCHIVE   , 'IGCM_OUT' )
 if wu.unDefined ('R_BUF'      ) : R_BUF       = os.path.join ( SCRATCHDIR, 'IGCM_OUT' )
@@ -157 +254 @@
 if wu.unDefined ('R_FIGR'     ) : R_FIGR      = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP )
 
-config['libIGCM'] = { 'R_OUT':R_OUT, 'R_BUF':R_BUF, 'L_EXP':L_EXP, 'R_BUFR':R_BUFR, 'POST_DIR':POST_DIR, 'REBUILD_DIR':REBUILD_DIR, 'R_BUF_KSH':R_BUF_KSH, 'R_FIGR':R_FIGR }
-
-# Set directory to extract files
-if wu.unDefined ( 'RunDir' ) : RunDir = os.path.join ( TmpDir, f'WATER_{JobName}_{YearBegin}_{YearEnd}' )
-config['Files']['RunDir'] = RunDir
-
-if not os.path.isdir (RunDir) : os.makedirs (RunDir)
-
-# Set directories to rebuild ocean and ice restart files
-RunDirOCE = os.path.join (RunDir, 'OCE')
-RunDirICE = os.path.join (RunDir, 'ICE')
-if not os.path.exists (RunDirOCE) : os.mkdir (RunDirOCE )
-if not os.path.exists (RunDirICE) : os.mkdir (RunDirICE )
+config['libIGCM'].update ( { 'R_OUT':R_OUT, 'R_BUF':R_BUF, 'L_EXP':L_EXP, 'R_BUFR':R_BUFR, 'POST_DIR':POST_DIR,
+                      'REBUILD_DIR':REBUILD_DIR, 'R_BUF_KSH':R_BUF_KSH, 'R_FIGR':R_FIGR, 'rebuild':rebuild } )
+
+## Set directory to extract files
+## ------------------------------
+if wu.unDefined ( 'FileDir' ) : FileDir = os.path.join ( TmpDir, f'WATER_{JobName}' )
+config['Files']['FileDir'] = FileDir
+
+if not os.path.isdir ( FileDir ) : os.makedirs ( FileDir )
+
+##- Set directories to rebuild ocean and ice restart files
+if wu.unDefined ( 'FileDirOCE' ) : FileDirOCE = os.path.join ( FileDir, 'OCE' )
+if wu.unDefined ( 'FileDirICE' ) : FileDirICE = os.path.join ( FileDir, 'ICE' )
+if not os.path.exists ( FileDirOCE ) : os.mkdir ( FileDirOCE )
+if not os.path.exists ( FileDirICE ) : os.mkdir ( FileDirICE )
 
 echo (' ')
-echo ( f'JobName   : {JobName}'   )
-echo ( f'Comment   : {Comment}'   )
-echo ( f'TmpDir    : {TmpDir}'    )
-echo ( f'RunDirOCE : {RunDirOCE}' )
-echo ( f'RunDirICE : {RunDirICE}' )
+echo ( f'JobName    : {JobName}'   )
+echo ( f'Comment    : {Comment}'   )
+echo ( f'TmpDir     : {TmpDir}'    )
+echo ( f'FileDir    : {FileDir}'    )
+echo ( f'FileDirOCE : {FileDirOCE}' )
+echo ( f'FileDirICE : {FileDirICE}' )
 
 echo ( f'\nDealing with {L_EXP}'  )
 
-#-- Creates model output directory names
+## Creates model output directory names
+## ------------------------------------
 if Freq == "MO" : FreqDir = os.path.join ( 'Output' , 'MO' )
 if Freq == "SE" : FreqDir = os.path.join ( 'Analyse', 'SE' )
@@ -196 +297 @@
 #-- Creates files names
 if wu.unDefined ( 'Period ' ) :
-    if Freq == 'MO' : Period = f'{YearBegin}0101_{YearEnd}1231_1M'
-    if Freq == 'SE' : Period = f'SE_{YearBegin}0101_{YearEnd}1231_1M'
+    if Freq == 'MO' : Period = f'{DateBegin}_{DateEnd}_1M'
+    if Freq == 'SE' : Period = f'SE_{DateBegin}_{DateEnd}_1M'
     config['Files']['Period'] = Period
+
+config['Files']['DateBegin'] = DateBegin
+config['Files']['DateBegin'] = DateEnd
+       
+echo ( f'Period   : {Period}' )
+
 if wu.unDefined ( 'FileCommon' ) :
     FileCommon = f'{JobName}_{Period}'
@@ -204 +311 @@
 
 if wu.unDefined ( 'Title' ) :
-    Title = f'{JobName} : {Freq} : {YearBegin}-01-01 - {YearEnd}-12-31'
+    Title = f'{JobName} : {Freq} : {DateBegin} - {DateEnd}'
     config['Files']['Title'] = Title
         
@@ -228 +335 @@
 
 ## Compute run length
+## ------------------
 dtime = ( d_OCE_his.time_counter_bounds.max() - d_OCE_his.time_counter_bounds.min() )
 echo ( f'\nRun length : {(dtime/np.timedelta64(1, "D")).values:8.2f} days' )
@@ -233 +341 @@
 
 ## Compute length of each period
+## -----------------------------
 dtime_per = ( d_OCE_his.time_counter_bounds[:,-1] - d_OCE_his.time_counter_bounds[:,0] )
 echo ( f'\nPeriods lengths (days) : ')
@@ -243 +352 @@
 NbYear = dtime_sec / YearLength
 
-#-- Open restart files
-YearRes = YearBegin - 1              # Year of the restart of beginning of simulation
-YearPre = YearBegin - PackFrequency  # Year to find the tarfile of the restart of beginning of simulation
-
-echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ')
-
-if wu.unDefined ( 'TarRestartPeriod_beg' ) : 
-    echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ')
-    TarRestartPeriod_beg = f'{YearPre}0101_{YearRes}1231'
+##-- Extract restart files from tar
+
+if wu.unDefined ('TarRestartDate_beg' ) : TarRestartDate_beg = wu.DateMinusOneDay ( DateBegin )
+if wu.unDefined ('TarRestartDate_end' ) : TarRestartDate_end = wu.FormatToGregorian ( DateEnd )
+
+if wu.unDefined ( 'TarRestartPeriod_beg' ) :
+
+    TarRestartPeriod_beg_DateEnd = TarRestartDate_beg
+    TarRestartPeriod_beg_DateBeg = wu.DateAddYear ( TarRestartPeriod_beg_DateEnd, -PackFrequency )
+    TarRestartPeriod_beg_DateBeg = wu.DatePlusOneDay ( TarRestartPeriod_beg_DateBeg )
+    
+    TarRestartPeriod_beg = f'{TarRestartPeriod_beg_DateBeg}_{TarRestartPeriod_beg_DateEnd}'
+    echo (f'Tar period for initial restart : {TarRestartPeriod_beg}')
     config['Files']['TarRestartPeriod_beg'] = TarRestartPeriod_beg
 
-if wu.unDefined ( 'TarRestartPeriod_end' ) : 
-    YearPre = YearBegin - PackFrequency  # Year to find the tarfile of the restart of beginning of simulation
-    echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ')
-    TarRestartPeriod_end = f'{YearBegin}0101_{YearEnd}1231'
+if wu.unDefined ( 'TarRestartPeriod_end' ) :
+
+    TarRestartPeriod_end_DateEnd = TarRestartDate_end
+    TarRestartPeriod_end_DateBeg = wu.DateAddYear ( TarRestartPeriod_end_DateEnd, -PackFrequency )
+    TarRestartPeriod_end_DateBeg = wu.DatePlusOneDay ( TarRestartPeriod_end_DateBeg )
+     
+    TarRestartPeriod_end = f'{TarRestartPeriod_end_DateBeg}_{TarRestartPeriod_end_DateEnd}'
+    echo (f'Tar period for final restart : {TarRestartPeriod_end}')
     config['Files']['TarRestartPeriod_end'] = TarRestartPeriod_end
+
+echo (f'Restart dates - Start : {TarRestartPeriod_beg}  /  End : {TarRestartPeriod_end}')
 
 if wu.unDefined ( 'tar_restart_beg' ) :
@@ -267 +386 @@
     config['Files']['tar_restart_end'] = tar_restart_end
 
-echo ( f'{tar_restart_beg}' )
-echo ( f'{tar_restart_end}' )
-
+echo ( f'{tar_restart_beg = }' )
+echo ( f'{tar_restart_end = }' )
+    
 #
 if wu.unDefined ('tar_restart_beg_ATM' ) : tar_restart_beg_ATM = tar_restart_beg
@@ -286 +405 @@
 
 if wu.unDefined ( 'file_OCE_beg' ) : 
-    file_OCE_beg = f'{RunDir}/OCE_{JobName}_{YearRes}1231_restart.nc'
+    file_OCE_beg = f'{FileDir}/OCE_{JobName}_{TarRestartDate_beg}_restart.nc'
     config['Files']['file_OCE_beg'] = file_OCE_beg
 if wu.unDefined ( 'file_OCE_end' ) :
-    file_OCE_end = f'{RunDir}/OCE_{JobName}_{YearEnd}1231_restart.nc'
+    file_OCE_end = f'{FileDir}/OCE_{JobName}_{TarRestartDate_end}_restart.nc'
     config['Files']['file_OCE_end'] = file_OCE_end
 if wu.unDefined ( 'file_ICE_beg' ) :
-    file_ICE_beg = f'{RunDir}/ICE_{JobName}_{YearRes}1231_restart_icemod.nc'
+    file_ICE_beg = f'{FileDir}/ICE_{JobName}_{TarRestartDate_beg}_restart_icemod.nc'
     config['Files']['file_ICE_beg'] = file_ICE_beg
 if wu.unDefined ( 'file_ICE_end' ) : 
-    file_ICE_end = f'{RunDir}/ICE_{JobName}_{YearEnd}1231_restart_icemod.nc'
+    file_ICE_end = f'{FileDir}/ICE_{JobName}_{TarRestartDate_end}_restart_icemod.nc'
     config['Files']['file_ICE_end'] = file_ICE_end
 
@@ -302 +421 @@
 echo ( f'{file_ICE_beg}' )
 echo ( f'{file_ICE_end}' )
+
+liste_beg = [file_OCE_beg, file_ICE_beg ]
+liste_end = [file_ATM_end, file_ICE_end ]
 
 echo ('\nExtract and rebuild OCE and ICE restarts')
@@ -311 +433 @@
     return int (ndomain_opa)
 
-def extract_and_rebuild ( file_name=file_OCE_beg, tar_restart=tar_restart_end, RunDirComp=RunDirOCE, ErrorCount=ErrorCount ) :
+def extract_and_rebuild ( file_name=file_OCE_beg, tar_restart=tar_restart_end, FileDirComp=FileDirOCE, ErrorCount=ErrorCount ) :
     '''Extract restart file from tar. Rebuild ocean files if needed'''
     echo ( f'----------')
@@ -321 +443 @@
         # Try to extract the rebuilded file
         if os.path.exists ( tar_restart ) :
-            command =  f'cd {RunDirComp} ; tar xf {tar_restart} {base_resFile}.nc'
+            command =  f'cd {FileDirComp} ; tar xf {tar_restart} {base_resFile}.nc'
             echo ( command )
             try : 
                 os.system ( command )
             except :
-                if not os.path.exists ( os.path.join (RunDir, f'{base_resFile}_0000.nc') ):
-                    command =  f'cd {RunDirComp} ; tar xf {tar_restart_end} {base_file}_*.nc'
+                if not os.path.exists ( os.path.join (FileDir, f'{base_resFile}_0000.nc') ):
+                    command =  f'cd {FileDirComp} ; tar xf {tar_restart_end} {base_file}_*.nc'
                     echo ( command )
                     ierr = os.system ( command )
@@ -333 +455 @@
                     else         : raise Exception ( f'command failed : {command}' )
                     echo ( f'extract ndomain' )
-                ndomain_opa = get_ndomain ( os.path.join (RunDir, f'{base_file}') )
-                command = f'cd {RunDirComp} ; {rebuild} {base_resFile} {ndomain_opa} ; mv {base_resFile}.nc {RunDir}'
+                ndomain_opa = get_ndomain ( os.path.join (FileDir, f'{base_file}') )
+                command = f'cd {FileDirComp} ; {rebuild} {base_resFile} {ndomain_opa} ; mv {base_resFile}.nc {FileDir}'
                 echo ( command )
                 ierr = os.system ( command )
@@ -341 +463 @@
             else : 
                 echo ( f'tar done : {base_resFile}')
-                command = f'cd {RunDirComp} ; mv {base_resFile}.nc {RunDir}'
+                command = f'cd {FileDirComp} ; mv {base_resFile}.nc {FileDir}'
                 ierr = os.system ( command )
                 if ierr == 0 : echo ( f'command done : {command}' )
@@ -354 +476 @@
 
             
-ErrorCount += extract_and_rebuild ( file_name=file_OCE_beg, tar_restart=tar_restart_beg, RunDirComp=RunDirOCE )
-ErrorCount += extract_and_rebuild ( file_name=file_OCE_end, tar_restart=tar_restart_end, RunDirComp=RunDirOCE )
-ErrorCount += extract_and_rebuild ( file_name=file_ICE_beg, tar_restart=tar_restart_beg, RunDirComp=RunDirICE )
-ErrorCount += extract_and_rebuild ( file_name=file_ICE_end, tar_restart=tar_restart_end, RunDirComp=RunDirICE )
+ErrorCount += extract_and_rebuild ( file_name=file_OCE_beg, tar_restart=tar_restart_beg, FileDirComp=FileDirOCE )
+ErrorCount += extract_and_rebuild ( file_name=file_OCE_end, tar_restart=tar_restart_end, FileDirComp=FileDirOCE )
+ErrorCount += extract_and_rebuild ( file_name=file_ICE_beg, tar_restart=tar_restart_beg, FileDirComp=FileDirICE )
+ErrorCount += extract_and_rebuild ( file_name=file_ICE_end, tar_restart=tar_restart_end, FileDirComp=FileDirICE )
 
 echo ('Opening OCE and ICE restart files')
 if NEMO == 3.6 : 
-    d_OCE_beg = xr.open_dataset ( os.path.join (RunDir, file_OCE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
-    d_OCE_end = xr.open_dataset ( os.path.join (RunDir, file_OCE_end), decode_times=False, decode_cf=True).squeeze()
-    d_ICE_beg = xr.open_dataset ( os.path.join (RunDir, file_ICE_beg), decode_times=False, decode_cf=True).squeeze()
-    d_ICE_end = xr.open_dataset ( os.path.join (RunDir, file_ICE_end), decode_times=False, decode_cf=True).squeeze()
+    d_OCE_beg = xr.open_dataset ( os.path.join (FileDir, file_OCE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
+    d_OCE_end = xr.open_dataset ( os.path.join (FileDir, file_OCE_end), decode_times=False, decode_cf=True).squeeze()
+    d_ICE_beg = xr.open_dataset ( os.path.join (FileDir, file_ICE_beg), decode_times=False, decode_cf=True).squeeze()
+    d_ICE_end = xr.open_dataset ( os.path.join (FileDir, file_ICE_end), decode_times=False, decode_cf=True).squeeze()
 if NEMO == 4.0 or NEMO == 4.2 : 
-    d_OCE_beg = xr.open_dataset ( os.path.join (RunDir, file_OCE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
-    d_OCE_end = xr.open_dataset ( os.path.join (RunDir, file_OCE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
-    d_ICE_beg = xr.open_dataset ( os.path.join (RunDir, file_ICE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
-    d_ICE_end = xr.open_dataset ( os.path.join (RunDir, file_ICE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
+    d_OCE_beg = xr.open_dataset ( os.path.join (FileDir, file_OCE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
+    d_OCE_end = xr.open_dataset ( os.path.join (FileDir, file_OCE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
+    d_ICE_beg = xr.open_dataset ( os.path.join (FileDir, file_ICE_beg), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
+    d_ICE_end = xr.open_dataset ( os.path.join (FileDir, file_ICE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze()
 
 ## Write the full configuration
@@ -377 +499 @@
 
     
-def kg2Sv  (val, rho=OCE_rho_liq) :
-    '''From kg to Sverdrup'''
-    return val/dtime_sec*1.0e-6/rho
-
-def kg2myear (val, rho=OCE_rho_liq) :
-    '''From kg to m/year'''
-    return val/OCE_aire_tot/rho/NbYear
-
-def var2prt (var, small=False) :
-    if small :  return var , kg2Sv(var)*1000., kg2myear(var)*1000.
-    else     :  return var , kg2Sv(var)      , kg2myear(var)
-
-def prtFlux (Desc, var, Form='F', small=False) :
-    if small :
-        if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} mSv | {:12.4f} mm/year "
-        if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} mSv | {:12.4e} mm/year "
-    else :
-        if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} Sv  | {:12.4f} m/year  "
-        if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} Sv  | {:12.4e} m/year  "
-    echo ( (' {:>15} = ' +ff).format (Desc, *var2prt(var, small) ) )
-    return None
-
 # Get mask and surfaces
 sos = d_OCE_his ['sos'][0].squeeze()
@@ -670 +770 @@
 echo ('                        = (dSEA_mas_wat + emp_oce + emp_ice - runoffs - iceshelf)*/abs(dSEA_mas_wat)      = {:12.1e} (-)  '.format ( (dSEA_mas_wat + OCE_mas_emp_oce + OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf)/abs(dSEA_mas_wat) ) )
 
+prtFlux ('   Leak OCE             ', ( dOCE_mas_wat + OCE_mas_empmr - OCE_mas_iceshelf) , 'e' )
+prtFlux ('   Leak ICE+SNW+PND     ', ( dICE_mas_wat + OCE_mas_emp_ice + OCE_mas_wfxice + OCE_mas_wfxsnw + ICE_mas_wfxpnd + ICE_mas_wfxsub_err ) , 'e' )
+prtFlux ('   Leak OCE+ICE+SNW+PND ',  ( dSEA_mas_wat + OCE_mas_emp_oce +OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf ) , 'e' )
+
+
+
 
 # check if emp     = emp_ice + emp_oce - calving

TOOLS/WATER_BUDGET/WaterUtils.py

@@ -1 +1 @@
 #!/usr/bin/env python3
-###
-### Script to check water conservation in the IPSL coupled model
-###
-## Here, some common utilities to all scripts
-##
-##  Warning, to install, configure, run, use any of included software or
-##  to read the associated documentation you'll need at least one (1)
-##  brain in a reasonably working order. Lack of this implement will
-##  void any warranties (either express or implied).  Authors assumes
-##  no responsability for errors, omissions, data loss, or any other
-##  consequences caused directly or indirectly by the usage of his
-##  software by incorrectly or partially configured personal
-##
-##
-## SVN information
-#  $Author$
-#  $Date$
-#  $Revision$
-#  $Id: ATM_waterbudget.py 6277 2022-12-08 09:24:05Z omamce $
-#  $HeadURL$
-
-## Import system modules
-#import sys, os, shutil, subprocess, platform
+'''
+  Script to check water conservation in the IPSL coupled model
+ 
+  Here, some common utilities to all scripts
+   
+  Warning, to install, configure, run, use any of included software or
+  to read the associated documentation you'll need at least one (1)
+  brain in a reasonably working order. Lack of this implement will
+  void any warranties (either express or implied).  Authors assumes
+  no responsability for errors, omissions, data loss, or any other
+  consequences caused directly or indirectly by the usage of his
+  software by incorrectly or partially configured personal
+
+
+ SVN information
+ $Author$
+ $Date$
+ $Revision$
+ $Id: ATM_waterbudget.py 6277 2022-12-08 09:24:05Z omamce $
+ $HeadURL$
+'''
+
 import numpy as np
 import configparser, re
 
-VarInt = 10
-Ra     = None
-
-def ftest () :
-    print ( toto )
+VarInt      = 10
+Rho         = 1000.0
+Ra          = 6366197.7236758135 #-- Earth Radius (m)
+Grav        = 9.81               #-- Gravity (m^2/s
+ICE_rho_ice = 917.0              #-- Ice volumic mass (kg/m3) in LIM3
+ICE_rho_sno = 330.0              #-- Snow volumic mass (kg/m3) in LIM3
+OCE_rho_liq = 1026.0             #-- Ocean water volumic mass (kg/m3) in NEMO
+ATM_rho     = 1000.0             #-- Water volumic mass in atmosphere (kg/m^3)
+SRF_rho     = 1000.0             #-- Water volumic mass in surface reservoir (kg/m^3)
+RUN_rho     = 1000.0             #-- Water volumic mass of rivers (kg/m^3)
+ICE_rho_pnd = 1000.              #-- Water volumic mass in ice ponds (kg/m^3)
+YearLength  = 365.25 * 24. * 60. * 60. #-- Year length (s)
 
 def setBool (chars) :
@@ -43 +49 @@
     '''Convert specific char string in integer or real if possible'''
     if type (chars) == str :
-        realnum = re.compile ("^[-+]?[0-9]*\.?[0-9]+(e[-+]?[0-9]+)?$")
+        realnum  = re.compile ("^[-+]?[0-9]*\.?[0-9]+(e[-+]?[0-9]+)?$")
         isNumber = realnum.match(chars.strip()) != None
         isInt    = chars.strip().isdigit()
@@ -63 +69 @@
 
 def unDefined (char) :
-    if char in globals () :
-        if globals()[char] == None : unDefined = True
+    '''True if a variable is not defined, or set to None'''
+    if char in locals () :
+        if locals()[char] == None : unDefined = True
         else            : unDefined = False
     else                : unDefined = True
     return unDefined
+
+def Defined (char) :
+    '''True if a variable is defined and not equal to None'''
+    if char in globals () :
+        if globals()[char] == None : Defined = False
+        else            : Defined = True
+    else                : Defined = False
+    return Defined
 
 def Ksum (tab) :
@@ -87 +102 @@
 def Ssum (tab) :
     '''
-    Precision summation by sorting by absolute value
+    Precision summation by sorting by absolute values
     '''
     Ssum = np.sum ( tab[np.argsort(np.abs(tab))] )
@@ -102 +117 @@
 Psum = KSsum
 
-
-
-
+def IsLeapYear ( Year, CalendarType="Gregorian" ) :
+    ''' True is Year is a leap year ''' 
+
+    # What is the calendar :
+    if CalendarType in [ '360d', '360_day', 'noleap', '365_day'] : 
+        IsLeapYear = False
+        return IsLeapYear
+
+    if CalendarType in [ 'all_leap', '366_day' ] : 
+        IsLeapYear = True
+        return IsLeapYear
+    
+    Year = int ( Year )
+  
+    # a year is a leap year if it is even divisible by 4
+    # but not evenly divisible by 100
+    # unless it is evenly divisible by 400
+
+    # if it is evenly divisible by 400 it must be a leap year
+    if np.mod ( Year, 400 ) == 0 : 
+        IsLeapYear = True
+        return IsLeapYear
+        
+    # if it is evenly divisible by 100 it must not be a leap year
+    if np.mod ( Year, 100 ) == 0 :
+        IsLeapYear = False
+        return IsLeapYear
+
+    # if it is evenly divisible by 4 it must be a leap year
+    if np.mod ( Year, 4 ) == 0 :
+        IsLeapYear = True
+        return IsLeapYear
+
+    IsLeapYear = False
+    return IsLeapYear
+
+def DateFormat ( Date ) :
+    '''
+    Get date format :
+      [yy]yymmdd   is Gregorian
+      [yy]yy-mm-dd is Human
+    '''
+    if type (Date) == str :
+        if '-' in Date : 
+            DateFormat = 'Human'
+        else : 
+            DateFormat = 'Gregorian'
+    if type (Date) == int : DateFormat = 'Gregorian'
+    return DateFormat
+
+def PrintDate ( YE, MO, DA, Format ) :
+    '''Return a date in the requested format'''
+    if Format == 'Human'     : PrintDate = f'{YE:04d}-{MO:02d}-{DA:02d}'
+    if Format == 'Gregorian' : PrintDate = f'{YE:04d}{MO:02d}{DA:02d}'
+    return PrintDate
+        
+def FormatToGregorian ( Date ) :
+    ''' from a yyyy-mm-dd or yyymmdd date format returns
+        a yyymmdd date format
+    '''
+    YE, MO, DA = SplitDate ( Date )
+    FormatToGregorian = f'{YE:04d}{MO:02d}{DA:02d}'
+    return FormatToGregorian
+  
+def FormatToHuman ( Date ) :
+    ''' From a yyyymmdd or yyymmdd date format returns
+        a yyy-mm-dd date format
+    '''
+    YE, MO, DA = SplitDate ( Date )
+    FormatToHuman =  f'{YE_new:04d}-{MO:02d}-{DA:02d}'
+    return FormatToHuman
+
+def SplitDate  ( Date, Out='int' ) :
+    ''' Split Date in format [yy]yymmdd or [yy]yy-mm-dd
+        to yy, mm, dd '''
+    if type (Date) == str :
+        if '-' in Date : 
+            Dash = True
+            YE, MO, DA = Date.split ('-')
+        else : 
+            Dash = False
+            YE = Date[:-4] ; MO = Date[-4:-2] ; DA = Date[-2:]
+            YearDigits = len (YE)
+    if type (Date) == int :
+        DA = np.mod ( Date, 100) 
+        MO = np.mod ( Date//100, 100)
+        YE = Date // 10000
+    
+    YE = int(YE) ; MO = int(MO) ; DA=int(DA)
+    return YE, MO, DA
+
+def DateAddYear ( Date, YearInc=1 ) :
+    ''' Add on year to date in format [yy]yymmdd or [yy]yy-mm-dd'''
+    Format = DateFormat ( Date )
+    YE, MO, DA = SplitDate ( Date )
+    YE_new = YE + YearInc
+    DateAddYear = PrintDate ( YE_new, MO, DA, Format)
+    return DateAddYear 
+
+def DateMinusOneDay ( Date ) :
+    ''' Substracts one day to date in format [yy]yymmdd or [yy]yy-mm-dd'''
+    Format = DateFormat ( Date )
+    mth_length = np.array ( [31, 28, 31,  30,  31,  30,  31,  31,  30,  31,  30,  31] )
+    YE, MO, DA = SplitDate ( Date )
+    if IsLeapYear ( YE) : mth_length[1] += 1
+    
+    YE = int(YE) ; MO = int(MO) ; DA=int(DA)
+    if DA ==  1 : 
+        if MO == 1 : DA_new = mth_length[-1]   ; MO_new = 12     ; YE_new = YE - 1
+        else       : DA_new = mth_length[MO-2] ; MO_new = MO - 1 ; YE_new = YE
+    else        :    DA_new = DA - 1           ; MO_new = MO     ; YE_new = YE
+
+    DateMinusOneDay = PrintDate ( YE_new, MO_new, DA_new, Format)
+    return DateMinusOneDay
+
+def DatePlusOneDay ( Date ) :
+    ''' Add one day to date in format [yy]yymmdd or [yy]yy-mm-dd'''
+    Format = DateFormat ( Date )
+    mth_length = np.array ( [31, 28, 31,  30,  31,  30,  31,  31,  30,  31,  30,  31] )
+    YE, MO, DA = SplitDate ( Date )
+    if IsLeapYear ( YE ) : mth_length[1] += 1
+    
+    YE_new = YE ; MO_new=MO ; DA_new=DA+1
+    if DA_new > mth_length [MO_new-1] :
+        DA_new = 1 ; MO_new = MO_new + 1
+        if MO_new == 13 : 
+            MO_new = 1 ; YE_new = YE_new+1
+
+    DatePlusOneDay = PrintDate ( YE_new, MO_new, DA_new, Format )
+    return DatePlusOneDay