Ignore:
Timestamp:
10/25/23 11:31:03 (9 months ago)
Author:
omamce
Message:

O.M. : WATER BUDGET

Improved code with pylint analysis

File:
1 edited

Legend:

Unmodified
Added
Removed

  • TOOLS/WATER_BUDGET/ATM_waterbudget.py

    r6651 r6665  
    2323### 
    2424## Import system modules 
    25 import sys, os, shutil#, subprocess, platform 
    26 import configparser, re 
     25import sys 
     26import os 
     27import configparser 
    2728 
    2829## Import needed scientific modules 
    29 import numpy as np, xarray as xr 
    30  
    31 # Check python version 
    32 if sys.version_info < (3, 8, 0) : 
    33     print ( f'Python version : {platform.python_version()}' )  
    34     raise Exception ( "Minimum Python version is 3.8" ) 
     30import numpy as np 
     31import xarray as xr 
    3532 
    3633## Import local modules 
    3734import WaterUtils as wu 
    3835import libIGCM_sys 
    39 import nemo, lmdz 
    40  
    41 from WaterUtils import VarInt, Rho, Ra, Grav, ICE_rho_ice, ICE_rho_sno, OCE_rho_liq, ATM_rho, SRF_rho, RUN_rho, ICE_rho_pnd, YearLength 
     36import lmdz 
     37 
     38from WaterUtils import RA, GRAV, ICE_RHO_ICE, ICE_RHO_SNO, \ 
     39                       OCE_RHO_LIQ, ATM_RHO, SRF_RHO, RUN_RHO, ICE_RHO_PND, YEAR_LENGTH 
    4240 
    4341## Creates parser for reading .ini input file 
    … …  
    4846## Experiment parameters 
    4947## --------------------- 
    50 ATM=None ; ATM_HIS='latlon' ; SRF_HIS='latlon' ; RUN_HIS='latlon' ; ORCA=None ; NEMO=None ; OCE_relax=False 
    51 OCE_icb=False ; Coupled=False ; Routing=None ; TestInterp=None 
    52 TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=None 
     48JobName=None ; TagName=None ; ExperimentName=None ; SpaceName=None 
     49SRF=True ; RUN=True 
     50ATM=None ; ATM_HIS='latlon' ; SRF_HIS='latlon' ; RUN_HIS='latlon' ; ORCA=None 
     51NEMO=None ; OCE_relax=False 
     52OCE_icb=False ; Coupled=False ; Rsouting=None ; TestInterp=None 
     53TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None 
     54Period=None ; Title=None 
    5355YearBegin=None ; YearEnd=None ; DateBegin=None ; DateEnd=None 
    54  
     56Freq=None ; libIGCM=None ; User=None; Group=None ; Routing=None 
     57PackFrequency = None ; FreqDir=None 
     58 
     59Timer=False ; Debug=False 
    5560## 
    56 ARCHIVE=None ; STORAGE=None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None ; TmpDir=None 
    57 FileDir=None ; FileOut=None 
     61ARCHIVE=None ; STORAGE=None ; SCRATCHDIR=None ; R_IN=None 
     62TmpDir=None ; FileDir=None ; FileOut=None ; R_OUT=None ; R_FIG=None 
     63R_FIGR=None ; R_BUF=None ; R_BUFR=None ; R_SAVE=None 
     64R_BUF_KSH=None ; POST_DIR=None ; L_EXP=None 
     65 
    5866dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None 
    5967FileCommon=None ; file_ATM_his=None ; file_SRF_his=None ; file_RUN_his=None 
    6068file_OCE_his=None ;  file_ICE_his=None ; file_OCE_sca=None 
    61 tar_restart_beg=None ; tar_restart_end=None ; file_ATM_beg=None ; file_ATM_end=None ; file_DYN_beg=None 
     69tar_restart_beg=None ; tar_restart_end=None ; file_ATM_beg=None 
     70file_ATM_end=None ; file_DYN_beg=None 
    6271file_DYN_end=None ; file_SRF_beg=None ; file_SRF_end=None 
    6372file_RUN_beg=None ; file_RUN_end=None ; file_RUN_end=None ; file_OCE_beg=None 
    6473file_ICE_beg=None ; file_OCE_beg=None 
    6574file_OCE_end=None ; file_ICE_beg=None ; file_OCE_end=None ; file_ICE_end=None 
     75TarRestartDate_beg=None ; TarRestartDate_end=None 
     76file_DYN_aire=None 
    6677tar_restart_beg_ATM=None ; tar_restart_beg_DYN=None ; tar_restart_beg_SRF=None 
    6778tar_restart_beg_RUN=None ; tar_restart_beg_OCE=None ; tar_restart_beg_ICE=None 
    6879tar_restart_end_ATM=None ; tar_restart_end_DYN=None ; tar_restart_end_SRF=None 
    6980tar_restart_end_RUN=None ; tar_restart_end_OCE=None ; tar_restart_end_ICE=None 
    70 ContinueOnError=False ; ErrorCount=0 # ; SortIco = False 
     81ContinueOnError=False ; ErrorCount=0 ; SortIco = False 
     82 
     83FileDirOCE=None ; FileDirATM=None ; FileDirICE=None ; FileDirSRF=None ; FileDirRUN=None 
    7184 
    7285## 
    … …  
    7487## --------------------------------- 
    7588# Default is float (full precision). Degrade the precision by using np.float32 
    76 # Restart file are always read at the full precision 
     89# Restart files are always read at the full precision 
    7790readPrec=float 
    7891 
    … …  
    88101if 'full' in IniFile : FullIniFile = IniFile 
    89102else                 : FullIniFile = 'full_' + IniFile 
    90      
     103 
    91104print ("Input file : ", IniFile ) 
    92105config.read (IniFile) 
    … …  
    109122        MyReader = configparser.ConfigParser (interpolation=configparser.ExtendedInterpolation() ) 
    110123        MyReader.optionxform = str # To keep capitals 
    111          
     124 
    112125        MyReader.read (ConfigCard) 
    113126 
    … …  
    120133                exec  ( f'wu.{VarName} = {VarName}' ) 
    121134                print ( f'    {VarName:21} set to : {locals()[VarName]:}' ) 
    122                  
     135 
    123136        for VarName in ['PackFrequency'] : 
    124137            if VarName in MyReader['Post'].keys() : 
    … …  
    131144    else : 
    132145        raise FileExistsError ( f"File not found : {ConfigCard = }" ) 
    133                  
     146 
    134147## Reading config file 
    135148## ------------------- 
    136 for Section in ['Config', 'Experiment', 'libIGCM', 'Files', 'Physics' ] : 
    137    if Section in config.keys () :  
    138         print ( f'\nReading [{Section}]' ) 
    139         for VarName in config[Section].keys() : 
    140             locals()[VarName] = config[Section][VarName] 
    141             exec  ( f'{VarName} = wu.setBool ({VarName})' ) 
    142             exec  ( f'{VarName} = wu.setNum  ({VarName})' ) 
    143             exec  ( f'{VarName} = wu.setNone ({VarName})' ) 
    144             exec  ( f'wu.{VarName} = {VarName}' ) 
    145             print ( f'    {VarName:21} set to : {locals()[VarName]}' ) 
    146             #exec ( f'del {VarName}' ) 
     149# Each entry in the .ini file will create a Python variable with the same name 
     150for Section in config.keys () :  
     151    print ( f'\nReading [{Section}]' ) 
     152    for VarName in config[Section].keys() : 
     153        locals()[VarName] = config[Section][VarName] 
     154        exec  ( f'{VarName} = wu.setBool ({VarName})' ) 
     155        exec  ( f'{VarName} = wu.setNum  ({VarName})' ) 
     156        exec  ( f'{VarName} = wu.setNone ({VarName})' ) 
     157        exec  ( f'wu.{VarName} = {VarName}' ) 
     158        print ( f'    {VarName:21} set to : {locals()[VarName]}' ) 
    147159 
    148160print ( f'\nConfig file readed : {IniFile} ' ) 
    … …  
    150162## 
    151163## Reading prec 
    152 if wu.unDefined ( 'readPrec' ) : 
     164if readPrec : 
    153165    readPrec = np.float64 
    154166else : 
    … …  
    156168    if readPrec in [         "float32", "r4", "single", "<class 'numpy.float32'>" ] : readPrec = np.float32 
    157169    if readPrec in [         "float16", "r2", "half"  , "<class 'numpy.float16'>" ] : readPrec = np.float16 
    158      
     170 
    159171## Some physical constants 
    160172## ======================= 
    161 if wu.unDefined ( 'Ra' )           : Ra          = wu.Ra           #-- Earth Radius (m) 
    162 if wu.unDefined ( 'Grav' )         : Grav        = wu.Grav         #-- Gravity (m^2/s 
    163 if wu.unDefined ( 'ICE_rho_ice' )  : ICE_rho_ice = wu.ICE_rho_ice  #-- Ice volumic mass (kg/m3) in LIM3 
    164 if wu.unDefined ( 'ICE_rho_sno')   : ICE_rho_sno = wu.ICE_rho_sno  #-- Snow volumic mass (kg/m3) in LIM3 
    165 if wu.unDefined ( 'OCE_rho_liq' )  : OCE_rho_liq = wu.OCE_rho_liq  #-- Ocean water volumic mass (kg/m3) in NEMO 
    166 if wu.unDefined ( 'ATM_rho' )      : ATM_rho     = wu.ATM_rho      #-- Water volumic mass in atmosphere (kg/m^3) 
    167 if wu.unDefined ( 'SRF_rho' )      : SRF_rho     = wu.SRF_rho      #-- Water volumic mass in surface reservoir (kg/m^3) 
    168 if wu.unDefined ( 'RUN_rho' )      : RUN_rho     = wu.RUN_rho      #-- Water volumic mass of rivers (kg/m^3) 
    169 if wu.unDefined ( 'ICE_rho_pnd' )  : ICE_rho_pnd = wu.ICE_rho_pnd  #-- Water volumic mass in ice ponds (kg/m^3) 
    170 if wu.unDefined ( 'YearLength' )   : YearLength  = wu.YearLength   #-- Year length (s) 
    171      
     173if not RA           : RA          = wu.RA           #-- Earth Radius (m) 
     174if not GRAV         : GRAV        = wu.GRAV         #-- Gravity (m^2/s 
     175if not ICE_RHO_ICE  : ICE_RHO_ICE = wu.ICE_RHO_ICE  #-- Ice volumic mass (kg/m3) in LIM3 
     176if not ICE_RHO_SNO  : ICE_RHO_SNO = wu.ICE_RHO_SNO  #-- Snow volumic mass (kg/m3) in LIM3 
     177if not OCE_RHO_LIQ  : OCE_RHO_LIQ = wu.OCE_RHO_LIQ  #-- Ocean water volumic mass (kg/m3) in NEMO 
     178if not ATM_RHO      : ATM_RHO     = wu.ATM_RHO      #-- Water volumic mass in atmosphere (kg/m^3) 
     179if not SRF_RHO      : SRF_RHO     = wu.SRF_RHO      #-- Water volumic mass in surface reservoir (kg/m^3) 
     180if not RUN_RHO      : RUN_RHO     = wu.RUN_RHO      #-- Water volumic mass of rivers (kg/m^3) 
     181if not ICE_RHO_PND  : ICE_RHO_PND = wu.ICE_RHO_PND  #-- Water volumic mass in ice ponds (kg/m^3) 
     182if not YEAR_LENGTH  : YEAR_LENGTH = wu.YEAR_LENGTH  #-- Year length (s) 
     183 
    172184## Set libIGCM and machine dependant values 
    173185## ---------------------------------------- 
    174 if not 'Files' in config.keys () : config['Files'] = {} 
    175  
    176 config['Physics'] = { 'Ra':str(Ra), 'Grav':str(Grav), 'ICE_rho_ice':str(ICE_rho_ice), 'ICE_rho_sno':str(ICE_rho_sno), 
    177                       'OCE_rho_liq':str(OCE_rho_liq), 'ATM_rho':str(ATM_rho), 'SRF_rho':str(SRF_rho), 'RUN_rho':str(RUN_rho)} 
    178  
    179 config['Config'] = { 'ContinueOnError':str(ContinueOnError), 'TestInterp':str(TestInterp), 'readPrec':str(readPrec) } 
     186if 'Files'   not in config.keys () : config['Files']   = {} 
     187if 'Physics' not in config.keys () : config['Physics'] = {} 
     188if 'Config'  not in config.keys () : config['Physics'] = {} 
     189 
     190 
     191config['Physics'].update ( { 'RA':str(RA), 'GRAV':str(GRAV), 'ICE_RHO_ICE':str(ICE_RHO_ICE), 'ICE_RHO_SNO':str(ICE_RHO_SNO), 
     192                             'OCE_RHO_LIQ':str(OCE_RHO_LIQ), 'ATM_RHO':str(ATM_RHO), 'SRF_RHO':str(SRF_RHO), 'RUN_RHO':str(RUN_RHO), 
     193                             'YEAR_LENGTH':str(YEAR_LENGTH)} ) 
     194 
     195config['Config'].update ( { 'ContinueOnError':str(ContinueOnError), 'TestInterp':str(TestInterp), 'readPrec':str(readPrec), 
     196                            'Debug':str(Debug), 'Timer':str(Timer)    } ) 
    180197 
    181198## -------------------------- 
    182 ICO  = ( 'ICO' in wu.ATM ) 
    183 LMDZ = ( 'LMD' in wu.ATM ) 
     199ICO  = 'ICO' in ATM 
     200LMDZ = 'LMD' in ATM 
    184201 
    185202mm = libIGCM_sys.config ( TagName=TagName, SpaceName=SpaceName, ExperimentName=ExperimentName, JobName=JobName, User=User, Group=Group, 
    186              ARCHIVE=None, SCRATCHDIR=None, STORAGE=None, R_IN=None, R_OUT=None, R_FIG=None, rebuild=None, TmpDir=None,  
    187              R_SAVE=None, R_FIGR=None, R_BUFR=None, R_BUF_KSH=None, REBUILD_DIR=None, POST_DIR=None ) 
     203             ARCHIVE=ARCHIVE, SCRATCHDIR=SCRATCHDIR, STORAGE=STORAGE, R_IN=R_IN, R_OUT=R_OUT, R_FIG=R_FIG, TmpDir=TmpDir, 
     204             R_SAVE=R_SAVE, R_FIGR=R_FIGR, R_BUFR=R_BUFR, R_BUF_KSH=R_BUF_KSH, POST_DIR=POST_DIR, L_EXP=L_EXP ) 
    188205globals().update(mm) 
    189206 
    190207config['Files']['TmpDir'] = TmpDir 
    191 config['libIGCM'] = { 'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'TmpDir':TmpDir, 'R_IN':R_IN, 'rebuild':rebuild } 
     208 
     209if 'libIGCM' not in config.keys () : config['libIGCM'] = {} 
     210config['libIGCM'].update ( { 'ARCHIVE':str(ARCHIVE), 'STORAGE':str(STORAGE), 'TmpDir':str(TmpDir), 'R_IN':str(R_IN) } ) 
     211 
     212## Debuging and timer 
     213Timer = wu.functools.partial (wu.Timer, debug=Debug, timer=Timer) 
    192214 
    193215## Defines begining and end of experiment 
    194216## -------------------------------------- 
    195 if wu.unDefined ( 'DateBegin' ) : 
     217if not DateBegin : 
    196218    DateBegin = f'{YearBegin}0101' 
    197219    config['Experiment']['DateBegin'] = str(DateBegin) 
    … …  
    201223    config['Experiment']['YearBegin'] = str(YearBegin) 
    202224 
    203 if wu.unDefined ( 'DateEnd' )   : 
     225if not DateEnd  : 
    204226    DateEnd   = f'{YearEnd}1231' 
    205227    config['Experiment']['DateEnd'] = str(DateEnd) 
    … …  
    209231    config['Experiment']['DateEnd'] = str(DateEnd) 
    210232 
    211 if wu.unDefined ( 'PackFrequency' ) : 
     233if not PackFrequency : 
    212234    PackFrequency = YearEnd - YearBegin + 1 
    213235    config['Experiment']['PackFrequency']   = f'{PackFrequency}' 
    214236 
    215 if type ( PackFrequency ) == str :  
    216     if 'Y' in PackFrequency : PackFrequency = PackFrequency.replace ( 'Y', '')  
     237if isinstance ( PackFrequency, str) : 
     238    if 'Y' in PackFrequency : PackFrequency = PackFrequency.replace ( 'Y', '') 
    217239    if 'M' in PackFrequency : PackFrequency = PackFrequency.replace ( 'M', '') 
    218240    PackFrequency = int ( PackFrequency ) 
    219      
     241 
    220242## Output file with water budget diagnostics 
    221243## ----------------------------------------- 
    222 if wu.unDefined ( 'FileOut' ) : 
     244if not FileOut : 
    223245    FileOut = f'ATM_waterbudget_{JobName}_{YearBegin}_{YearEnd}' 
    224     if ICO :  
     246    if ICO : 
    225247        if ATM_HIS == 'latlon' : FileOut = f'{FileOut}_LATLON' 
    226248        if ATM_HIS == 'ico'    : FileOut = f'{FileOut}_ICO' 
    … …  
    230252config['Files']['FileOut'] = FileOut 
    231253 
    232 f_out = open ( FileOut, mode = 'w' ) 
     254f_out = open ( FileOut, mode = 'w', encoding="utf-8" ) 
    233255 
    234256## Useful functions 
    235257## ---------------- 
    236 if readPrec == float : 
    237     def rprec (tab) : return tab 
     258if repr(readPrec) == "<class 'numpy.float64'>" or readPrec == float : 
     259    def rprec (ptab) : 
     260        '''This version does nothing 
     261 
     262        rprec may be used to reduce floating precision when reading history files 
     263        ''' 
     264        return ptab 
    238265else                 : 
    239     def rprec (tab) : return tab.astype(readPrec).astype(float) 
    240      
    241 def kg2Sv    (val, rho=ATM_rho) : 
     266    def rprec (ptab) : 
     267        '''Returns float with a different precision''' 
     268        return ptab.astype(readPrec).astype(float) 
     269 
     270def kg2Sv    (pval, rho=ATM_RHO) : 
    242271    '''From kg to Sverdrup''' 
    243     return val/dtime_sec*1.0e-6/rho 
    244  
    245 def kg2myear (val, rho=ATM_rho) : 
     272    return pval/dtime_sec*1.0e-6/rho 
     273 
     274def kg2myear (pval, rho=ATM_RHO) : 
    246275    '''From kg to m/year''' 
    247     return val/ATM_aire_sea_tot/rho/NbYear 
    248  
    249 def var2prt (var, small=False, rho=ATM_rho) : 
    250     if small :  return var , kg2Sv(var, rho=rho)*1000., kg2myear(var, rho=rho)*1000 
    251     else     :  return var , kg2Sv(var, rho=rho)      , kg2myear(var, rho=rho) 
    252  
    253 def prtFlux (Desc, var, Form='F', small=False, rho=ATM_rho, width=15) : 
     276    return pval/ATM_aire_sea_tot/rho/NbYear 
     277 
     278def var2prt (pvar, small=False, rho=ATM_RHO) : 
     279    '''Formats value for printing''' 
     280    if small :  return pvar, kg2Sv (pvar, rho=rho)*1000., kg2myear (pvar, rho=rho)*1000 
     281    else     :  return pvar, kg2Sv (pvar, rho=rho)      , kg2myear (pvar, rho=rho) 
     282 
     283def prtFlux (Desc, pvar, Form='F', small=False, rho=ATM_RHO, width=15) : 
     284    '''Pretty print of formattd value''' 
    254285    if small : 
    255286        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} mSv | {:12.4f} mm/year " 
    … …  
    258289        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} Sv  | {:12.4f} m/year  " 
    259290        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} Sv  | {:12.4e} m/year  " 
    260     echo ( (' {:>{width}} = ' +ff).format (Desc, *var2prt(var, small=small, rho=rho), width=width ) ) 
    261     return None 
     291    echo ( (' {:>{width}} = ' +ff).format (Desc, *var2prt (pvar, small=small, rho=rho), width=width ) ) 
    262292 
    263293def echo (string, end='\n') : 
    … …  
    267297    f_out.write ( str(string) + end ) 
    268298    f_out.flush () 
    269     return None 
    270299 
    271300echo ( f'{ContinueOnError = }' ) 
    … …  
    274303echo ( f'{JobName         = }' ) 
    275304echo ( f'{ConfigCard      = }' ) 
    276 echo ( f'{libIGCM         = }' )      
    277 echo ( f'{User            = }' )        
    278 echo ( f'{Group           = }' )        
    279 echo ( f'{Freq            = }' )        
    280 echo ( f'{YearBegin       = }' )      
    281 echo ( f'{YearEnd         = }' )      
     305echo ( f'{libIGCM         = }' ) 
     306echo ( f'{User            = }' ) 
     307echo ( f'{Group           = }' ) 
     308echo ( f'{Freq            = }' ) 
     309echo ( f'{YearBegin       = }' ) 
     310echo ( f'{YearEnd         = }' ) 
    282311echo ( f'{DateBegin       = }' ) 
    283312echo ( f'{DateEnd         = }' ) 
    284 echo ( f'{PackFrequency   = }' )  
    285 echo ( f'{ATM             = }' )        
    286 echo ( f'{Routing         = }' )        
    287 echo ( f'{ORCA            = }' )       
    288 echo ( f'{NEMO            = }' )       
    289 echo ( f'{Coupled         = }' )       
    290 echo ( f'{ATM_HIS         = }' )       
    291 echo ( f'{SRF_HIS         = }' )       
    292 echo ( f'{RUN_HIS         = }' )       
     313echo ( f'{PackFrequency   = }' ) 
     314echo ( f'{ATM             = }' ) 
     315echo ( f'{Routing         = }' ) 
     316echo ( f'{ORCA            = }' ) 
     317echo ( f'{NEMO            = }' ) 
     318echo ( f'{Coupled         = }' ) 
     319echo ( f'{ATM_HIS         = }' ) 
     320echo ( f'{SRF_HIS         = }' ) 
     321echo ( f'{RUN_HIS         = }' ) 
    293322 
    294323## Set libIGCM directories 
    295324## ----------------------- 
    296 if wu.unDefined ('R_OUT'      ) : R_OUT       = os.path.join ( ARCHIVE   , 'IGCM_OUT' ) 
    297 if wu.unDefined ('R_BUF'      ) : R_BUF       = os.path.join ( SCRATCHDIR, 'IGCM_OUT' ) 
    298 if wu.unDefined ('L_EXP'      ) : L_EXP       = os.path.join (TagName, SpaceName, ExperimentName, JobName) 
    299 if wu.unDefined ('R_SAVE'     ) : R_SAVE      = os.path.join ( R_OUT, L_EXP ) 
    300 if wu.unDefined ('R_BUFR'     ) : R_BUFR      = os.path.join ( R_BUF, L_EXP ) 
    301 if wu.unDefined ('POST_DIR'   ) : POST_DIR    = os.path.join ( R_BUFR, 'Out' ) 
    302 if wu.unDefined ('REBUILD_DIR') : REBUILD_DIR = os.path.join ( R_BUFR, 'REBUILD' ) 
    303 if wu.unDefined ('R_BUF_KSH'  ) : R_BUF_KSH   = os.path.join ( R_BUFR, 'Out' ) 
    304 if wu.unDefined ('R_FIGR'     ) : R_FIGR      = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP ) 
    305  
    306 config['libIGCM'].update ( { 'R_OUT':R_OUT, 'R_BUF':R_BUF, 'L_EXP':L_EXP, 'R_BUFR':R_BUFR, 'POST_DIR':POST_DIR, 
    307                       'REBUILD_DIR':REBUILD_DIR, 'R_BUF_KSH':R_BUF_KSH, 'R_FIGR':R_FIGR, 'rebuild':rebuild } ) 
     325if not R_OUT       : R_OUT       = os.path.join ( ARCHIVE   , 'IGCM_OUT' ) 
     326if not R_BUF       : R_BUF       = os.path.join ( SCRATCHDIR, 'IGCM_OUT' ) 
     327if not L_EXP       : 
     328    if TagName and SpaceName and ExperimentName and JobName : L_EXP = os.path.join (TagName, SpaceName, ExperimentName, JobName) 
     329    else : L_EXP = '/' 
     330if not R_SAVE      : R_SAVE      = os.path.join ( R_OUT, L_EXP ) 
     331if not R_BUFR      : R_BUFR      = os.path.join ( R_BUF, L_EXP ) 
     332if not POST_DIR    : POST_DIR    = os.path.join ( R_BUFR, 'Out' ) 
     333if not R_BUF_KSH   : R_BUF_KSH   = os.path.join ( R_BUFR, 'Out' ) 
     334if not R_FIGR      : R_FIGR      = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP ) 
     335 
     336config['libIGCM'].update ( { 'R_OUT':R_OUT, 'R_BUF':R_BUF, 'L_EXP':L_EXP, 'R_BUFR':R_BUFR, 'POST_DIR':POST_DIR, 'R_SAVE':R_SAVE, 
     337                       'R_BUF_KSH':R_BUF_KSH, 'R_FIGR':R_BUFR } ) 
    308338 
    309339## Set directory to extract files 
    310340## ------------------------------ 
    311 if wu.unDefined ( 'FileDir' ) : FileDir = os.path.join ( TmpDir, f'WATER_{JobName}' ) 
     341if not FileDir : FileDir = os.path.join ( TmpDir, f'WATER_{JobName}' ) 
    312342config['Files']['FileDir'] = FileDir 
    313343 
    314344if not os.path.isdir ( FileDir ) : os.makedirs ( FileDir ) 
    315  
    316 ##- Set directories to rebuild ocean and ice restart files 
    317 if wu.unDefined ( 'FileDirOCE' ) : FileDirOCE = os.path.join ( FileDir, 'OCE' ) 
    318 if wu.unDefined ( 'FileDirICE' ) : FileDirICE = os.path.join ( FileDir, 'ICE' ) 
    319 if not os.path.exists ( FileDirOCE ) : os.mkdir ( FileDirOCE ) 
    320 if not os.path.exists ( FileDirICE ) : os.mkdir ( FileDirICE ) 
    321345 
    322346echo (' ') 
    … …  
    325349echo ( f'TmpDir      : {TmpDir}'     ) 
    326350echo ( f'FileDir     : {FileDir}'    ) 
    327 echo ( f'FileDirOCE  : {FileDirOCE}' ) 
    328 echo ( f'FileDirICE  : {FileDirICE}' ) 
    329351 
    330352echo ( f'\nDealing with {L_EXP}'  ) 
    … …  
    334356if Freq == "MO" : FreqDir = os.path.join ('Output' , 'MO' ) 
    335357if Freq == "SE" : FreqDir = os.path.join ('Analyse', 'SE' ) 
    336 if wu.unDefined ('dir_ATM_his' ) : 
     358if not dir_ATM_his : 
    337359    dir_ATM_his = os.path.join ( R_SAVE, "ATM", FreqDir ) 
    338360    config['Files']['dir_ATM_his'] = dir_ATM_his 
    339 if wu.unDefined ( 'dir_SRF_his' ) :  
     361if not dir_SRF_his : 
    340362    dir_SRF_his = os.path.join ( R_SAVE, "SRF", FreqDir ) 
    341363    config['Files']['dir_SRF_his'] = dir_SRF_his 
    342      
    343 echo ( f'The analysis relies on files from the following model output directories : ' ) 
     364 
     365echo (  'The analysis relies on files from the following model output directories : ' ) 
    344366echo ( f'{dir_ATM_his = }' ) 
    345367echo ( f'{dir_SRF_his = }' ) 
    346368 
    347369##-- Creates files names 
    348 if wu.unDefined ( 'Period' ) : 
     370if not Period : 
    349371    if Freq == 'MO' : Period = f'{DateBegin}_{DateEnd}_1M' 
    350372    if Freq == 'SE' : Period = f'SE_{DateBegin}_{DateEnd}_1M' 
    … …  
    353375config['Files']['DateBegin'] = DateBegin 
    354376config['Files']['DateBegin'] = DateEnd 
    355         
     377 
    356378echo ( f'Period   : {Period}' ) 
    357      
    358 if wu.unDefined ( 'FileCommon' ) : 
     379 
     380if not FileCommon : 
    359381    FileCommon = f'{JobName}_{Period}' 
    360382    config['Files']['FileCommon'] = FileCommon 
    361383 
    362 if wu.unDefined ( 'Title' ) : 
     384if not Title : 
    363385    Title = f'{JobName} : {Freq} : {DateBegin} - {DateEnd}' 
    364386    config['Files']['Title'] = Title 
    365       
     387 
    366388echo ('\nOpen history files' ) 
    367 if wu.unDefined ( 'file_ATM_his' ) : 
     389if not file_ATM_his : 
    368390    if ATM_HIS == 'latlon' : 
    369391        file_ATM_his = os.path.join ( dir_ATM_his, f'{FileCommon}_histmth.nc' ) 
    … …  
    371393        file_ATM_his = os.path.join ( dir_ATM_his, f'{FileCommon}_histmth_ico.nc' ) 
    372394    config['Files']['file_ATM_his'] = file_ATM_his 
    373 if wu.unDefined ( 'file_SRF_his' ) : 
     395if not file_SRF_his : 
    374396    if ATM_HIS == 'latlon' : 
    375397        file_SRF_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 
    … …  
    378400    config['Files']['file_SRF_his'] = file_SRF_his 
    379401 
    380 if Routing == 'SIMPLE' : 
    381     if file_RUN_his == None : 
     402if SRF and Routing == 'SIMPLE' : 
     403    if file_RUN_his is None : 
    382404        if ATM_HIS == 'latlon' : 
    383405            file_RUN_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 
    … …  
    387409 
    388410d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 
    389 d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 
    390 if Routing == 'SECHIBA' : d_RUN_his = d_SRF_his 
    391 if Routing == 'SIMPLE'  : d_RUN_his = xr.open_dataset ( file_RUN_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 
     411if SRF : 
     412    d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 
     413    if Routing == 'SECHIBA' : d_RUN_his = d_SRF_his 
     414    if Routing == 'SIMPLE'  : d_RUN_his = xr.open_dataset ( file_RUN_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 
    392415 
    393416echo ( f'{file_ATM_his = }' ) 
    394 echo ( f'{file_SRF_his = }' ) 
    395 if Routing == 'SIMPLE' : echo ( f'{file_RUN_his = }' ) 
     417if SRF : 
     418    echo ( f'{file_SRF_his = }' ) 
     419    if Routing == 'SIMPLE' : echo ( f'{file_RUN_his = }' ) 
    396420 
    397421## Compute run length 
    … …  
    409433 
    410434##-- Number of years (approximative) 
    411 NbYear = dtime_sec / YearLength 
     435NbYear = dtime_sec / YEAR_LENGTH 
    412436 
    413437##-- Extract restart files from tar 
    414438 
    415 if wu.unDefined ('TarRestartDate_beg' ) : TarRestartDate_beg = wu.DateMinusOneDay ( DateBegin ) 
    416 if wu.unDefined ('TarRestartDate_end' ) : TarRestartDate_end = wu.FormatToGregorian ( DateEnd ) 
    417  
    418 if wu.unDefined ( 'TarRestartPeriod_beg' ) : 
     439if not TarRestartDate_beg : TarRestartDate_beg = wu.DateMinusOneDay ( DateBegin ) 
     440if not TarRestartDate_end : TarRestartDate_end = wu.FormatToGregorian ( DateEnd ) 
     441 
     442if not TarRestartPeriod_beg : 
    419443 
    420444    TarRestartPeriod_beg_DateEnd = TarRestartDate_beg 
    421445    TarRestartPeriod_beg_DateBeg = wu.DateAddYear ( TarRestartPeriod_beg_DateEnd, -PackFrequency ) 
    422446    TarRestartPeriod_beg_DateBeg = wu.DatePlusOneDay ( TarRestartPeriod_beg_DateBeg ) 
    423      
     447 
    424448    TarRestartPeriod_beg = f'{TarRestartPeriod_beg_DateBeg}_{TarRestartPeriod_beg_DateEnd}' 
    425449    echo (f'Tar period for initial restart : {TarRestartPeriod_beg}') 
    426450    config['Files']['TarRestartPeriod_beg'] = TarRestartPeriod_beg 
    427451 
    428 if wu.unDefined ( 'TarRestartPeriod_end' ) : 
     452if not TarRestartPeriod_end : 
    429453 
    430454    TarRestartPeriod_end_DateEnd = TarRestartDate_end 
    431455    TarRestartPeriod_end_DateBeg = wu.DateAddYear ( TarRestartPeriod_end_DateEnd, -PackFrequency ) 
    432456    TarRestartPeriod_end_DateBeg = wu.DatePlusOneDay ( TarRestartPeriod_end_DateBeg ) 
    433       
     457 
    434458    TarRestartPeriod_end = f'{TarRestartPeriod_end_DateBeg}_{TarRestartPeriod_end_DateEnd}' 
    435459    echo (f'Tar period for final restart : {TarRestartPeriod_end}') 
    … …  
    438462echo (f'Restart dates - Start : {TarRestartPeriod_beg}  /  End : {TarRestartPeriod_end}') 
    439463 
    440 if wu.unDefined ( 'tar_restart_beg' ) : 
     464if not tar_restart_beg : 
    441465    tar_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_beg}_restart.tar' ) 
    442466    config['Files']['tar_restart_beg'] = tar_restart_beg 
    443 if wu.unDefined ( 'tar_restart_end' ) : 
     467if not tar_restart_end : 
    444468    tar_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_end}_restart.tar' ) 
    445469    config['Files']['tar_restart_end'] = tar_restart_end 
    446      
     470 
    447471echo ( f'{tar_restart_beg = }' ) 
    448472echo ( f'{tar_restart_end = }' ) 
    449473 
    450474##-- Names of tar files with restarts 
    451 if wu.unDefined ( 'SRF_HIS' ) : SRF_HIS = ATM_HIS 
    452  
    453 if wu.unDefined ( 'tar_restart_beg_ATM' ) : tar_restart_beg_ATM = tar_restart_beg 
    454 if wu.unDefined ( 'tar_restart_beg_DYN' ) : tar_restart_beg_DYN = tar_restart_beg 
    455 if wu.unDefined ( 'tar_restart_beg_SRF' ) : tar_restart_beg_SRF = tar_restart_beg 
    456 if wu.unDefined ( 'tar_restart_beg_RUN' ) : tar_restart_beg_RUN = tar_restart_beg 
    457 if wu.unDefined ( 'tar_restart_beg_OCE' ) : tar_restart_beg_OCE = tar_restart_beg 
    458 if wu.unDefined ( 'tar_restart_beg_ICE' ) : tar_restart_beg_ICE = tar_restart_beg 
    459  
    460 if wu.unDefined ( 'tar_restart_end_ATM' ) : tar_restart_end_ATM = tar_restart_end 
    461 if wu.unDefined ( 'tar_restart_end_DYN' ) : tar_restart_end_DYN = tar_restart_end 
    462 if wu.unDefined ( 'tar_restart_end_SRF' ) : tar_restart_end_SRF = tar_restart_end 
    463 if wu.unDefined ( 'tar_restart_end_RUN' ) : tar_restart_end_RUN = tar_restart_end 
    464 if wu.unDefined ( 'tar_restart_end_OCE' ) : tar_restart_end_OCE = tar_restart_end 
    465 if wu.unDefined ( 'tar_restart_end_ICE' ) : tar_restart_end_ICE = tar_restart_end 
    466  
    467 if wu.unDefined ( 'file_ATM_beg' ) :  
     475 
     476if not tar_restart_beg_ATM : tar_restart_beg_ATM = tar_restart_beg 
     477if not tar_restart_beg_DYN : tar_restart_beg_DYN = tar_restart_beg 
     478if not tar_restart_beg_RUN : tar_restart_beg_RUN = tar_restart_beg 
     479if not tar_restart_beg_OCE : tar_restart_beg_OCE = tar_restart_beg 
     480if not tar_restart_beg_ICE : tar_restart_beg_ICE = tar_restart_beg 
     481 
     482if not tar_restart_end_ATM : tar_restart_end_ATM = tar_restart_end 
     483if not tar_restart_end_DYN : tar_restart_end_DYN = tar_restart_end 
     484if not tar_restart_end_RUN : tar_restart_end_RUN = tar_restart_end 
     485if not tar_restart_end_OCE : tar_restart_end_OCE = tar_restart_end 
     486if not tar_restart_end_ICE : tar_restart_end_ICE = tar_restart_end 
     487 
     488if SRF : 
     489    if not SRF_HIS : SRF_HIS = ATM_HIS 
     490    if not tar_restart_beg_SRF : tar_restart_beg_SRF = tar_restart_beg 
     491    if not tar_restart_end_SRF : tar_restart_end_SRF = tar_restart_end 
     492 
     493if not file_ATM_beg : 
    468494    file_ATM_beg = f'{FileDir}/ATM_{JobName}_{TarRestartDate_beg}_restartphy.nc' 
    469495    config['Files']['file_ATM_beg'] = file_ATM_beg 
    470 if wu.unDefined ( 'file_ATM_end' ) : 
     496if not file_ATM_end : 
    471497    file_ATM_end = f'{FileDir}/ATM_{JobName}_{TarRestartDate_end}_restartphy.nc' 
    472498    config['Files']['file_ATM_end'] = file_ATM_end 
    … …  
    475501liste_end = [file_ATM_end, ] 
    476502 
    477 if wu.unDefined ( 'file_DYN_beg' ) : 
     503if not file_DYN_beg : 
    478504    if LMDZ : file_DYN_beg = f'{FileDir}/ATM_{JobName}_{TarRestartDate_beg}_restart.nc' 
    479505    if ICO  : file_DYN_beg = f'{FileDir}/ICO_{JobName}_{TarRestartDate_beg}_restart.nc' 
    480506    liste_beg.append (file_DYN_beg) 
    481507    config['Files']['file_DYN_beg'] = file_DYN_beg 
    482      
    483 if wu.unDefined ( 'file_DYN_end' ) :  
     508 
     509if not file_DYN_end : 
    484510    if LMDZ : file_DYN_end = f'{FileDir}/ATM_{JobName}_{TarRestartDate_end}_restart.nc' 
    485511    if ICO  : file_DYN_end = f'{FileDir}/ICO_{JobName}_{TarRestartDate_end}_restart.nc' 
    … …  
    487513    config['Files']['file_DYN_end'] = file_DYN_end 
    488514 
    489 if wu.unDefined ( 'file_SRF_beg' ) : 
    490     file_SRF_beg = f'{FileDir}/SRF_{JobName}_{TarRestartDate_beg}_sechiba_rest.nc' 
    491     config['Files']['file_SRF_beg'] = file_SRF_beg 
    492 if wu.unDefined ( 'file_SRF_end' ) : 
    493     file_SRF_end = f'{FileDir}/SRF_{JobName}_{TarRestartDate_end}_sechiba_rest.nc' 
    494     config['Files']['file_SRF_end'] = file_SRF_end 
    495      
     515if SRF : 
     516    if not file_SRF_beg : 
     517        file_SRF_beg = f'{FileDir}/SRF_{JobName}_{TarRestartDate_beg}_sechiba_rest.nc' 
     518        config['Files']['file_SRF_beg'] = file_SRF_beg 
     519    if not file_SRF_end : 
     520        file_SRF_end = f'{FileDir}/SRF_{JobName}_{TarRestartDate_end}_sechiba_rest.nc' 
     521        config['Files']['file_SRF_end'] = file_SRF_end 
     522 
    496523liste_beg.append ( file_SRF_beg ) 
    497524liste_end.append ( file_SRF_end ) 
    … …  
    501528echo ( f'{file_DYN_beg = }') 
    502529echo ( f'{file_DYN_end = }') 
    503 echo ( f'{file_SRF_beg = }') 
    504 echo ( f'{file_SRF_end = }') 
     530if SRF : 
     531    echo ( f'{file_SRF_beg = }') 
     532    echo ( f'{file_SRF_end = }') 
    505533 
    506534if ICO : 
    507     if wu.unDefined ('file_DYN_aire') : file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ATM+'_grid.nc' ) 
     535    if not file_DYN_aire : file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ATM+'_grid.nc' ) 
    508536    config['Files']['file_DYN_aire'] = file_DYN_aire 
    509537 
    510 if Routing == 'SIMPLE' : 
    511     if wu.unDefined ( 'file_RUN_beg' ) :  
     538if SRF and Routing == 'SIMPLE' : 
     539    if not file_RUN_beg : 
    512540        file_RUN_beg = f'{FileDir}/SRF_{JobName}_{TarRestartDate_beg}_routing_restart.nc' 
    513541        config['Files']['file_RUN_beg'] = file_RUN_beg 
    514     if wu.unDefined ( 'file_RUN_end' ) :  
     542    if not file_RUN_end : 
    515543        file_RUN_end = f'{FileDir}/SRF_{JobName}_{TarRestartDate_end}_routing_restart.nc' 
    516544        config['Files']['file_RUN_end'] = file_RUN_end 
    517          
     545 
    518546    liste_beg.append ( file_RUN_beg ) 
    519547    liste_end.append ( file_RUN_end ) 
    … …  
    521549    echo ( f'{file_RUN_end = }' ) 
    522550 
    523 echo ('\nExtract restart files from tar : ATM, ICO and SRF') 
    524  
    525 def extract_and_rebuild ( file_name=file_ATM_beg, tar_restart=tar_restart_end, FileDirComp=FileDir, ErrorCount=ErrorCount ) : 
     551echo ('\nExtract restart files from tar : ATM, ICO', end='') 
     552if SRF : echo ( ' and SRF') 
     553else   : echo (' ') 
     554 
     555    
     556@Timer 
     557def extract ( file_name=file_ATM_beg, tar_restart=tar_restart_end, file_dir_comp=FileDir, error_count=ErrorCount ) : 
     558    ''' 
     559    Extract restart files from a tar 
     560    ''' 
    526561    echo ( f'----------') 
    527562    echo ( f'file to extract : {file_name = }' ) 
    … …  
    532567        base_resFile = os.path.basename (file_name) 
    533568        if os.path.exists ( tar_restart ) : 
    534             command =  f'cd {FileDir} ; tar xf {tar_restart} {base_resFile}' 
     569            command =  f'cd {file_dir_comp} ; tar xf {tar_restart} {base_resFile}' 
    535570            echo ( f'{command = }' ) 
    536571            try : os.system ( command ) 
    537572            except : 
    538573                if ContinueOnError : 
    539                     ErrorCount += 1 
    540                     echo ( f'****** Command failed : {command}' ) 
    541                     echo ( f'****** Trying to continue' ) 
     574                    error_count += 1 
     575                    echo ( '****** Command failed : {command}' ) 
     576                    echo ( '****** Trying to continue' ) 
    542577                    echo ( ' ') 
    543578                else : 
    544                     raise Exception ( f'**** command failed : {command} - Stopping' ) 
     579                    raise OSError ( f'**** command failed : {command} - Stopping' ) 
    545580            else : 
    546581                echo ( f'tar done : {base_resFile}' ) 
    … …  
    548583            echo ( f'****** Tar restart file {tar_restart = } not found ' ) 
    549584            if ContinueOnError : 
    550                 ErrorCount += 1 
    551             else :  
    552                 raise Exception ( f'****** tar file not found {tar_restart = } - Stopping' ) 
    553     return ErrorCount 
    554   
    555 ErrorCount += extract_and_rebuild ( file_name=file_ATM_beg, tar_restart=tar_restart_beg_ATM, FileDirComp=FileDir ) 
    556 ErrorCount += extract_and_rebuild ( file_name=file_DYN_beg, tar_restart=tar_restart_beg_DYN, FileDirComp=FileDir ) 
    557 ErrorCount += extract_and_rebuild ( file_name=file_SRF_beg, tar_restart=tar_restart_beg_SRF, FileDirComp=FileDir ) 
    558  
    559 ErrorCount += extract_and_rebuild ( file_name=file_ATM_end, tar_restart=tar_restart_end_ATM, FileDirComp=FileDir ) 
    560 ErrorCount += extract_and_rebuild ( file_name=file_DYN_end, tar_restart=tar_restart_end_DYN, FileDirComp=FileDir ) 
    561 ErrorCount += extract_and_rebuild ( file_name=file_SRF_end, tar_restart=tar_restart_end_SRF, FileDirComp=FileDir ) 
    562  
    563 if Routing == 'SIMPLE' : 
    564     ErrorCount += extract_and_rebuild ( file_name=file_RUN_beg, tar_restart=tar_restart_beg_RUN, FileDirComp=FileDir ) 
    565     ErrorCount += extract_and_rebuild ( file_name=file_RUN_end, tar_restart=tar_restart_end_RUN, FileDirComp=FileDir ) 
     585                error_count += 1 
     586            else : 
     587                raise OSError ( f'****** tar file not found {tar_restart = } - Stopping' ) 
     588    return error_count 
     589 
     590ErrorCount += extract ( file_name=file_ATM_beg, tar_restart=tar_restart_beg_ATM, file_dir_comp=FileDir, error_count=ErrorCount ) 
     591ErrorCount += extract ( file_name=file_DYN_beg, tar_restart=tar_restart_beg_DYN, file_dir_comp=FileDir, error_count=ErrorCount ) 
     592 
     593ErrorCount += extract ( file_name=file_ATM_end, tar_restart=tar_restart_end_ATM, file_dir_comp=FileDir, error_count=ErrorCount ) 
     594ErrorCount += extract ( file_name=file_DYN_end, tar_restart=tar_restart_end_DYN, file_dir_comp=FileDir, error_count=ErrorCount ) 
     595 
     596if SRF : 
     597    ErrorCount += extract ( file_name=file_SRF_beg, tar_restart=tar_restart_beg_SRF, file_dir_comp=FileDir, error_count=ErrorCount ) 
     598    ErrorCount += extract ( file_name=file_SRF_end, tar_restart=tar_restart_end_SRF, file_dir_comp=FileDir, error_count=ErrorCount ) 
     599 
     600    if Routing == 'SIMPLE' : 
     601        ErrorCount += extract ( file_name=file_RUN_beg, tar_restart=tar_restart_beg_RUN, file_dir_comp=FileDir, error_count=ErrorCount ) 
     602        ErrorCount += extract ( file_name=file_RUN_end, tar_restart=tar_restart_end_RUN, file_dir_comp=FileDir, error_count=ErrorCount ) 
    566603 
    567604##-- Exit in case of error in the opening file phase 
    568605if ErrorCount > 0 : 
    569606    echo ( '  ' ) 
    570     raise Exception ( f'**** Some files missing - Stopping - {ErrorCount = }' ) 
    571  
    572 ##  
     607    raise RuntimeError ( f'**** Some files missing - Stopping - {ErrorCount = }' ) 
     608 
     609## 
    573610echo ('\nOpening ATM SRF and ICO restart files') 
    574611d_ATM_beg = xr.open_dataset ( os.path.join (FileDir, file_ATM_beg), decode_times=False, decode_cf=True ).squeeze() 
    575612d_ATM_end = xr.open_dataset ( os.path.join (FileDir, file_ATM_end), decode_times=False, decode_cf=True ).squeeze() 
    576 d_SRF_beg = xr.open_dataset ( os.path.join (FileDir, file_SRF_beg), decode_times=False, decode_cf=True ).squeeze() 
    577 d_SRF_end = xr.open_dataset ( os.path.join (FileDir, file_SRF_end), decode_times=False, decode_cf=True ).squeeze() 
     613if SRF : 
     614    d_SRF_beg = xr.open_dataset ( os.path.join (FileDir, file_SRF_beg), decode_times=False, decode_cf=True ).squeeze() 
     615    d_SRF_end = xr.open_dataset ( os.path.join (FileDir, file_SRF_end), decode_times=False, decode_cf=True ).squeeze() 
    578616d_DYN_beg = xr.open_dataset ( os.path.join (FileDir, file_DYN_beg), decode_times=False, decode_cf=True ).squeeze() 
    579617d_DYN_end = xr.open_dataset ( os.path.join (FileDir, file_DYN_end), decode_times=False, decode_cf=True ).squeeze() 
    580618 
    581 for var in d_SRF_beg.variables : 
    582     d_SRF_beg[var] = d_SRF_beg[var].where ( d_SRF_beg[var]<1.e20, 0.) 
    583     d_SRF_end[var] = d_SRF_end[var].where ( d_SRF_end[var]<1.e20, 0.) 
    584  
    585 if Routing == 'SIMPLE' : 
    586     d_RUN_beg = xr.open_dataset ( os.path.join (FileDir, file_RUN_beg), decode_times=False, decode_cf=True ).squeeze() 
    587     d_RUN_end = xr.open_dataset ( os.path.join (FileDir, file_RUN_end), decode_times=False, decode_cf=True ).squeeze() 
    588  
     619if SRF : 
     620    for var in d_SRF_beg.variables : 
     621        d_SRF_beg[var] = d_SRF_beg[var].where ( d_SRF_beg[var]<1.e20, 0.) 
     622        d_SRF_end[var] = d_SRF_end[var].where ( d_SRF_end[var]<1.e20, 0.) 
     623 
     624    if Routing == 'SIMPLE' : 
     625        d_RUN_beg = xr.open_dataset ( os.path.join (FileDir, file_RUN_beg), decode_times=False, decode_cf=True ).squeeze() 
     626        d_RUN_end = xr.open_dataset ( os.path.join (FileDir, file_RUN_end), decode_times=False, decode_cf=True ).squeeze() 
     627 
     628@Timer 
    589629def to_cell ( dd, newname='cell' ) : 
    590     '''Set space dimension to 'cell' ''' 
     630    '''Set space dimension to newname 
     631    ''' 
    591632    for oldname in [ 'cell_mesh', 'y', 'points_physiques' ] : 
    592         try    : dd = dd.rename ( {oldname : newname} ) 
    593         except : pass 
     633        if oldname in dd.dims and oldname != newname : 
     634            dd = dd.rename ( {oldname : newname} ) 
    594635    return dd 
    595636 
    596637d_ATM_beg = to_cell ( d_ATM_beg ) 
    597638d_ATM_end = to_cell ( d_ATM_end ) 
    598 d_SRF_beg = to_cell ( d_SRF_beg ) 
    599 d_SRF_end = to_cell ( d_SRF_end ) 
     639if SRF : 
     640    d_SRF_beg = to_cell ( d_SRF_beg ) 
     641    d_SRF_end = to_cell ( d_SRF_end ) 
    600642d_DYN_beg = to_cell ( d_DYN_beg ) 
    601643d_DYN_end = to_cell ( d_DYN_end ) 
    602644 
    603 if Routing == 'SIMPLE' : 
     645if SRF and Routing == 'SIMPLE' : 
    604646    d_RUN_beg = to_cell ( d_RUN_beg ) 
    605647    d_RUN_end = to_cell ( d_RUN_end ) 
    606648 
    607649d_ATM_his = to_cell ( d_ATM_his ) 
    608 d_SRF_his = to_cell ( d_SRF_his ) 
    609      
     650if SRF : d_SRF_his = to_cell ( d_SRF_his ) 
     651 
    610652echo ( f'{file_ATM_beg = }' ) 
    611653echo ( f'{file_ATM_end = }' ) 
    612654echo ( f'{file_DYN_beg = }' ) 
    613655echo ( f'{file_DYN_end = }' ) 
    614 echo ( f'{file_SRF_beg = }' ) 
    615 echo ( f'{file_SRF_end = }' ) 
    616 if Routing == 'SIMPLE' : 
    617     echo ( f'{file_RUN_beg = }' ) 
    618     echo ( f'{file_RUN_end = }' ) 
     656if SRF : 
     657    echo ( f'{file_SRF_beg = }' ) 
     658    echo ( f'{file_SRF_end = }' ) 
     659    if Routing == 'SIMPLE' : 
     660        echo ( f'{file_RUN_beg = }' ) 
     661        echo ( f'{file_RUN_end = }' ) 
    619662 
    620663## Write the full configuration 
    621 config_out = open (FullIniFile, 'w') 
     664config_out = open (FullIniFile, 'w', encoding = 'utf-8') 
    622665config.write ( config_out ) 
    623666config_out.close () 
    … …  
    626669if LMDZ : 
    627670    echo ('ATM grid with cell surfaces : LMDZ') 
    628     ATM_lat       = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1D='cell' ) 
    629     ATM_lon       = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1D='cell' ) 
    630     ATM_aire      = lmdz.geo2point ( rprec (d_ATM_his ['aire']     [0]), cumulPoles=True, dim1D='cell' ) 
    631     ATM_fter      = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1D='cell' ) 
    632     ATM_foce      = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1D='cell' ) 
    633     ATM_fsic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1D='cell' ) 
    634     ATM_flic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1D='cell' ) 
    635     SRF_lat       = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1D='cell' ) 
    636     SRF_lon       = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1D='cell' ) 
    637     SRF_aire      = lmdz.geo2point ( rprec (d_SRF_his ['Areas']) * rprec (d_SRF_his ['Contfrac']), dim1D='cell', cumulPoles=True ) 
    638     SRF_areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas'])  ,  dim1D='cell', cumulPoles=True ) 
    639     SRF_contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac']), dim1D='cell' ) 
     671    ATM_lat       = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1d='cell' ) 
     672    ATM_lon       = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1d='cell' ) 
     673    ATM_aire      = lmdz.geo2point ( rprec (d_ATM_his ['aire']     [0]), cumul_poles=True, dim1d='cell' ) 
     674    ATM_fter      = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1d='cell' ) 
     675    ATM_foce      = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1d='cell' ) 
     676    ATM_fsic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1d='cell' ) 
     677    ATM_flic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1d='cell' ) 
     678    if SRF : 
     679        SRF_lat       = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1d='cell' ) 
     680        SRF_lon       = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1d='cell' ) 
     681        SRF_aire      = lmdz.geo2point ( rprec (d_SRF_his ['Areas']) * rprec (d_SRF_his ['Contfrac']), dim1d='cell', cumul_poles=True ) 
     682        SRF_areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas'])  ,  dim1d='cell', cumul_poles=True ) 
     683        SRF_contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac']), dim1d='cell' ) 
    640684 
    641685if ICO : 
    … …  
    643687        echo ( 'ATM areas and fractions on latlon grid' ) 
    644688        if 'lat_dom_out' in d_ATM_his.variables : 
    645             ATM_lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat_dom_out'])+0*rprec (d_ATM_his ['lon_dom_out']), dim1D='cell' ) 
    646             ATM_lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat_dom_out'])+  rprec (d_ATM_his ['lon_dom_out']), dim1D='cell' ) 
     689            ATM_lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat_dom_out'])+0*rprec (d_ATM_his ['lon_dom_out']), dim1d='cell' ) 
     690            ATM_lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat_dom_out'])+  rprec (d_ATM_his ['lon_dom_out']), dim1d='cell' ) 
    647691        else : 
    648             ATM_lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1D='cell' ) 
    649             ATM_lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1D='cell' ) 
    650         ATM_aire = lmdz.geo2point ( rprec (d_ATM_his ['aire'][0]).squeeze(), cumulPoles=True, dim1D='cell' ) 
    651         ATM_fter = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1D='cell' ) 
    652         ATM_foce = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1D='cell' ) 
    653         ATM_fsic = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1D='cell' ) 
    654         ATM_flic = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1D='cell' ) 
     692            ATM_lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1d='cell' ) 
     693            ATM_lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1d='cell' ) 
     694        ATM_aire = lmdz.geo2point ( rprec (d_ATM_his ['aire'][0]).squeeze(), cumul_poles=True, dim1d='cell' ) 
     695        ATM_fter = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1d='cell' ) 
     696        ATM_foce = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1d='cell' ) 
     697        ATM_fsic = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1d='cell' ) 
     698        ATM_flic = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1d='cell' ) 
    655699 
    656700    if ATM_HIS == 'ico' : 
    … …  
    664708        ATM_flic =  rprec (d_ATM_his ['fract_lic'][0]) 
    665709 
    666     if SRF_HIS == 'latlon' : 
    667         echo ( 'SRF areas and fractions on latlon grid' ) 
    668         if 'lat_domain_landpoints_out' in d_SRF_his  :  
    669             SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat_domain_landpoints_out'])+0*rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1D='cell' ) 
    670             SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_domain_landpoints_out'])+  rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1D='cell' ) 
    671         else :  
    672             if 'lat_domain_landpoints_out' in d_SRF_his : 
    673                 SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat_dom_out'])+0*rprec (d_SRF_his ['lon_dom_out']), dim1D='cell' ) 
    674                 SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_dom_out'])+  rprec (d_SRF_his ['lon_dom_out']), dim1D='cell' ) 
     710    if SRF : 
     711        if SRF_HIS == 'latlon' : 
     712            echo ( 'SRF areas and fractions on latlon grid' ) 
     713            if 'lat_domain_landpoints_out' in d_SRF_his  : 
     714                SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat_domain_landpoints_out'])+0*rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1d='cell' ) 
     715                SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_domain_landpoints_out'])+  rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1d='cell' ) 
    675716            else : 
    676                 SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1D='cell' ) 
    677                 SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1D='cell' ) 
    678                  
    679         SRF_areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas']   )   , dim1D='cell', cumulPoles=True ) 
    680         SRF_areafrac  = lmdz.geo2point ( rprec (d_SRF_his ['AreaFrac'])   , dim1D='cell', cumulPoles=True ) 
    681         SRF_contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac'])   , dim1D='cell', cumulPoles=True ) 
    682         SRF_aire = SRF_areafrac 
    683   
    684     if SRF_HIS == 'ico' : 
    685         echo ( 'SRF areas and fractions on latlon grid' ) 
    686         SRF_lat       =  rprec (d_SRF_his ['lat']     ) 
    687         SRF_lon       =  rprec (d_SRF_his ['lon']     ) 
    688         SRF_areas     =  rprec (d_SRF_his ['Areas']   )  
    689         SRF_contfrac  =  rprec (d_SRF_his ['Contfrac']) 
    690         SRF_aire      =  SRF_areas * SRF_contfrac 
     717                if 'lat_domain_landpoints_out' in d_SRF_his : 
     718                    SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat_dom_out'])+0*rprec (d_SRF_his ['lon_dom_out']), dim1d='cell' ) 
     719                    SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_dom_out'])+  rprec (d_SRF_his ['lon_dom_out']), dim1d='cell' ) 
     720                else : 
     721                    SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1d='cell' ) 
     722                    SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1d='cell' ) 
     723 
     724            SRF_areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas']   )   , dim1d='cell', cumul_poles=True ) 
     725            SRF_areafrac  = lmdz.geo2point ( rprec (d_SRF_his ['AreaFrac'])   , dim1d='cell', cumul_poles=True ) 
     726            SRF_contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac'])   , dim1d='cell', cumul_poles=True ) 
     727            SRF_aire = SRF_areafrac 
     728 
     729        if SRF_HIS == 'ico' : 
     730            echo ( 'SRF areas and fractions on latlon grid' ) 
     731            SRF_lat       =  rprec (d_SRF_his ['lat']     ) 
     732            SRF_lon       =  rprec (d_SRF_his ['lon']     ) 
     733            SRF_areas     =  rprec (d_SRF_his ['Areas']   ) 
     734            SRF_contfrac  =  rprec (d_SRF_his ['Contfrac']) 
     735            SRF_aire      =  SRF_areas * SRF_contfrac 
    691736 
    692737ATM_fsea      = ATM_foce + ATM_fsic 
    … …  
    702747 
    703748## Write the full configuration 
    704 config_out = open (FullIniFile, 'w') 
     749config_out = open (FullIniFile, 'w', encoding='utf-8') 
    705750config.write ( config_out ) 
    706751config_out.close () 
    707              
     752 
    708753if ICO : 
    709754    # Area on icosahedron grid 
    710755    d_DYN_aire = xr.open_dataset ( file_DYN_aire, decode_times=False ).squeeze() 
    711756 
    712     if SortIco :  
     757    if SortIco : 
    713758        # Creation d'une clef de tri pour le fichier aire 
    714759        DYN_aire_keysort = np.lexsort ( (d_DYN_aire['lat'], d_DYN_aire['lon']) ) 
    715     else :  
     760    else : 
    716761        DYN_aire_keysort = np.arange ( len ( d_DYN_aire['lat'] ) ) 
    717              
     762 
    718763    DYN_lat = d_DYN_aire['lat'] 
    719764    DYN_lon = d_DYN_aire['lon'] 
    … …  
    725770    DYN_flic = d_ATM_beg['FLIC'] 
    726771    DYN_aire_fter = DYN_aire * DYN_fter 
    727      
     772 
    728773if LMDZ : 
    729774    # Area on lon/lat grid 
    … …  
    736781 
    737782# Functions computing integrals and sum 
     783@Timer 
    738784def DYN_stock_int (stock) : 
    739785    '''Integrate (* surface) stock on atmosphere grid''' 
    740     DYN_stock_int  = wu.Psum ( (stock * DYN_aire).to_masked_array().ravel() )  
    741     return DYN_stock_int 
    742  
     786    return wu.Psum ( (stock * DYN_aire).to_masked_array().ravel() ) 
     787 
     788@Timer 
    743789def ATM_flux_int (flux) : 
    744790    '''Integrate (* time * surface) flux on atmosphere grid''' 
    745     ATM_flux_int  = wu.Psum ( (flux * dtime_per_sec * ATM_aire).to_masked_array().ravel() ) 
    746     return ATM_flux_int 
    747  
     791    return wu.Psum ( (flux * dtime_per_sec * ATM_aire).to_masked_array().ravel() ) 
     792 
     793@Timer 
    748794def LIC_flux_int (flux) : 
    749795    '''Integrate (* time * surface) flux on land ice grid''' 
    750     LIC_flux_int  = wu.Psum ( (flux * dtime_per_sec * ATM_aire_flic).to_masked_array().ravel() ) 
    751     return LIC_flux_int 
    752  
    753 def SRF_stock_int (stock) : 
    754     '''Integrate (* surface) stock on land grid''' 
    755     SRF_stock_int  = wu.Ksum (  ( (stock * DYN_aire_fter).to_masked_array().ravel()) ) 
    756     return SRF_stock_int 
    757  
    758 def SRF_flux_int (flux) : 
    759     '''Integrate (* time * surface) flux on land grid''' 
    760     SRF_flux_int  = wu.Psum (  (flux * dtime_per_sec * SRF_aire).to_masked_array().ravel() ) 
    761     return SRF_flux_int 
    762  
     796    return wu.Psum ( (flux * dtime_per_sec * ATM_aire_flic).to_masked_array().ravel() ) 
     797 
     798if SRF : 
     799    @Timer 
     800    def SRF_stock_int (stock) : 
     801        '''Integrate (* surface) stock on land grid''' 
     802        return wu.Ksum (  ( (stock * DYN_aire_fter).to_masked_array().ravel()) ) 
     803 
     804    @Timer 
     805    def SRF_flux_int (flux) : 
     806        '''Integrate (* time * surface) flux on land grid''' 
     807        return wu.Psum (  (flux * dtime_per_sec * SRF_aire).to_masked_array().ravel() ) 
     808 
     809@Timer 
    763810def ONE_stock_int (stock) : 
    764811    '''Sum stock ''' 
    765     ONE_stock_int  = wu.Psum ( stock.to_masked_array().ravel() ) 
    766     return ONE_stock_int 
    767  
     812    return wu.Psum ( stock.to_masked_array().ravel() ) 
     813 
     814@Timer 
    768815def ONE_flux_int (flux) : 
    769816    '''Integrate (* time) flux on area=1 grid ''' 
    770     ONE_flux_int  = wu.Psum ( (flux * dtime_per_sec ).to_masked_array().ravel() ) 
    771     return ONE_flux_int 
    772  
     817    return wu.Psum ( (flux * dtime_per_sec ).to_masked_array().ravel() ) 
     818 
     819@Timer 
    773820def Sprec ( tlist ) : 
    774821    '''Accurate sum of list of scalar elements''' 
    … …  
    783830 
    784831DYN_aire_tot     = ONE_stock_int ( DYN_aire ) 
    785 SRF_aire_tot     = ONE_stock_int ( SRF_aire ) 
     832if SRF : SRF_aire_tot     = ONE_stock_int ( SRF_aire ) 
    786833 
    787834echo ('') 
    … …  
    790837echo ( f'ATM HIS : Area of ter in atmosphere : {ATM_aire_ter_tot:18.9e}' ) 
    791838echo ( f'ATM HIS : Area of lic in atmosphere : {ATM_aire_lic_tot:18.9e}' ) 
    792 echo ( f'ATM SRF : Area of atmosphere        : {SRF_aire_tot:18.9e}' ) 
     839if SRF : 
     840    echo ( f'ATM SRF : Area of atmosphere        : {SRF_aire_tot:18.9e}' ) 
    793841echo ('') 
    794 echo ( 'ATM DYN : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(DYN_aire_tot/(Ra*Ra*4*np.pi) ) ) 
    795 echo ( 'ATM HIS : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(ATM_aire_tot/(Ra*Ra*4*np.pi) ) ) 
    796 echo ( 'ATM HIS : Area of ter in atmosphere/(4pi R^2) : {:18.9f}'.format(ATM_aire_ter_tot/(Ra*Ra*4*np.pi) ) ) 
    797 echo ( 'ATM SRF : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(SRF_aire_tot/(Ra*Ra*4*np.pi) ) ) 
    798 echo ('') 
    799 echo ( f'ATM SRF : Area of atmosphere (no contfrac): {ONE_stock_int (SRF_areas):18.9e}' ) 
    800  
    801  
    802 if (  np.abs (ATM_aire_tot/(Ra*Ra*4*np.pi) - 1.0) > 0.01 ) : 
    803     raise Exception ('Error of atmosphere surface interpolated on lon/lat grid') 
     842echo ( 'ATM DYN : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(DYN_aire_tot/(RA*RA*4*np.pi) ) ) 
     843echo ( 'ATM HIS : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(ATM_aire_tot/(RA*RA*4*np.pi) ) ) 
     844echo ( 'ATM HIS : Area of ter in atmosphere/(4pi R^2) : {:18.9f}'.format(ATM_aire_ter_tot/(RA*RA*4*np.pi) ) ) 
     845if SRF : 
     846    echo ( 'ATM SRF : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(SRF_aire_tot/(RA*RA*4*np.pi) ) ) 
     847    echo ('') 
     848    echo ( f'ATM SRF : Area of atmosphere (no contfrac): {ONE_stock_int (SRF_areas):18.9e}' ) 
     849 
     850 
     851if np.abs (ATM_aire_tot/(RA*RA*4*np.pi) - 1.0) > 0.01 : 
     852    raise RuntimeError ('Error of atmosphere surface interpolated on lon/lat grid') 
    804853 
    805854echo ( '\n====================================================================================' ) 
    … …  
    817866    DYN_psol_beg = d_DYN_beg['ps'] 
    818867    DYN_psol_end = d_DYN_end['ps'] 
    819 if LMDZ :  
    820     DYN_psol_beg = lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)), dim1D='cell' ) 
    821     DYN_psol_end = lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)), dim1D='cell' ) 
    822      
     868if LMDZ : 
     869    DYN_psol_beg = lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)), dim1d='cell' ) 
     870    DYN_psol_end = lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)), dim1d='cell' ) 
     871 
    823872echo ( '3D Pressure at the interface layers (not scalar points)' ) 
    824873DYN_pres_beg = ATM_Ahyb + ATM_Bhyb * DYN_psol_beg 
    … …  
    848897    echo ( f'(DYN_pres_end[-1]).min().item()              = {ind[6]}' ) 
    849898    echo ( f'(DYN_pres_end[-1]).max().item()              = {ind[7]}' ) 
    850     raise Exception 
    851      
     899    raise RuntimeError 
     900 
    852901klevp1 = ATM_Bhyb.shape[-1] 
    853902cell   = DYN_psol_beg.shape[-1] 
    … …  
    857906DYN_mass_beg = xr.DataArray ( np.empty( (klev, cell)), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) ) 
    858907DYN_mass_end = xr.DataArray ( np.empty( (klev, cell)), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) ) 
    859      
     908 
    860909for k in np.arange (klev) : 
    861     DYN_mass_beg[k,:] = ( DYN_pres_beg[k,:] - DYN_pres_beg[k+1,:] ) / Grav 
    862     DYN_mass_end[k,:] = ( DYN_pres_end[k,:] - DYN_pres_end[k+1,:] ) / Grav 
     910    DYN_mass_beg[k,:] = ( DYN_pres_beg[k,:] - DYN_pres_beg[k+1,:] ) / GRAV 
     911    DYN_mass_end[k,:] = ( DYN_pres_end[k,:] - DYN_pres_end[k+1,:] ) / GRAV 
    863912 
    864913DYN_mass_beg_2D = DYN_mass_beg.sum (dim='sigs') 
    … …  
    872921    if 'H2Ov' in d_DYN_beg.variables : 
    873922        echo ('reading LATLON : H2Ov, H2Ol, H2Oi' ) 
    874         DYN_wat_beg = lmdz.geo3point ( d_DYN_beg['H2Ov'] + d_DYN_beg['H2Ol'] + d_DYN_beg['H2Oi'].isel(rlonv=slice(0,-1) ), dim1D='cell' ) 
    875         DYN_wat_end = lmdz.geo3point ( d_DYN_end['H2Ov'] + d_DYN_end['H2Ol'] + d_DYN_end['H2Oi'].isel(rlonv=slice(0,-1) ), dim1D='cell' ) 
     923        DYN_wat_beg = lmdz.geo3point ( d_DYN_beg['H2Ov'] + d_DYN_beg['H2Ol'] + d_DYN_beg['H2Oi'].isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
     924        DYN_wat_end = lmdz.geo3point ( d_DYN_end['H2Ov'] + d_DYN_end['H2Ol'] + d_DYN_end['H2Oi'].isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
    876925    if 'H2Ov_g' in d_DYN_beg.variables : 
    877926        echo ('reading LATLON : H2O_g, H2O_l, H2O_s' ) 
    878         DYN_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s']).isel(rlonv=slice(0,-1) ), dim1D='cell' ) 
    879         DYN_wat_end = lmdz.geo3point ( (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s']).isel(rlonv=slice(0,-1) ), dim1D='cell' ) 
     927        DYN_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s']).isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
     928        DYN_wat_end = lmdz.geo3point ( (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s']).isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
    880929if ICO : 
    881930    if 'H2Ov_g' in d_DYN_beg.variables : 
    882931        echo ('reading ICO : H2O_g, H2O_l, H2O_s' ) 
    883         DYN_wat_beg = (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s']) 
    884         DYN_wat_end = (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s']) 
     932        DYN_wat_beg = d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s'] 
     933        DYN_wat_end = d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s'] 
    885934    elif 'H2O_g' in d_DYN_beg.variables : 
    886935        echo ('reading ICO : H2O_g, H2O_l, H2O_s' ) 
    887         DYN_wat_beg = (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s']) 
    888         DYN_wat_end = (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s']) 
     936        DYN_wat_beg = d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s'] 
     937        DYN_wat_end = d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s'] 
    889938    elif 'q' in d_DYN_beg.variables : 
    890939        echo ('reading ICO : q' ) 
    891         DYN_wat_beg = (d_DYN_beg['q'].isel(nq=0) + d_DYN_beg['q'].isel(nq=1) + d_DYN_beg['q'].isel(nq=2) ) 
    892         DYN_wat_end = (d_DYN_end['q'].isel(nq=0) + d_DYN_end['q'].isel(nq=1) + d_DYN_end['q'].isel(nq=2) ) 
    893  
    894 if 'lev' in DYN_wat_beg.dims :  
     940        DYN_wat_beg = d_DYN_beg['q'].isel(nq=0) + d_DYN_beg['q'].isel(nq=1) + d_DYN_beg['q'].isel(nq=2) 
     941        DYN_wat_end = d_DYN_end['q'].isel(nq=0) + d_DYN_end['q'].isel(nq=1) + d_DYN_end['q'].isel(nq=2) 
     942 
     943if 'lev' in DYN_wat_beg.dims : 
    895944    DYN_wat_beg = DYN_wat_beg.rename ( {'lev':'sigs'} ) 
    896945    DYN_wat_end = DYN_wat_end.rename ( {'lev':'sigs'} ) 
    897      
     946 
    898947echo ( 'Compute water content : vertical and horizontal integral' ) 
    899948 
    … …  
    940989ATM_qs02_end     = d_ATM_end['QS02'] * d_ATM_end['FLIC'] 
    941990ATM_qs03_end     = d_ATM_end['QS03'] * d_ATM_end['FOCE'] 
    942 ATM_qs04_end     = d_ATM_end['QS04'] * d_ATM_end['FSIC']   
    943  
    944 if ICO : 
    945      ATM_sno_beg     = ATM_sno_beg     
    946      ATM_sno_end     = ATM_sno_end     
    947      ATM_qs_beg      = ATM_qs_beg      
    948      ATM_qs_end      = ATM_qs_end      
    949      ATM_qsol_beg    = ATM_qsol_beg    
    950      ATM_qs01_beg    = ATM_qs01_beg    
    951      ATM_qs02_beg    = ATM_qs02_beg    
    952      ATM_qs03_beg    = ATM_qs03_beg    
    953      ATM_qs04_beg    = ATM_qs04_beg    
    954      ATM_qsol_end    = ATM_qsol_end    
    955      ATM_qs01_end    = ATM_qs01_end    
    956      ATM_qs02_end    = ATM_qs02_end    
    957      ATM_qs03_end    = ATM_qs03_end    
    958      ATM_qs04_end    = ATM_qs04_end    
    959      LIC_sno_beg     = LIC_sno_beg     
    960      LIC_sno_end     = LIC_sno_end     
    961      LIC_runlic0_beg = LIC_runlic0_beg 
    962      LIC_runlic0_end = LIC_runlic0_end 
    963     
     991ATM_qs04_end     = d_ATM_end['QS04'] * d_ATM_end['FSIC'] 
     992 
    964993LIC_qs_beg = ATM_qs02_beg 
    965994LIC_qs_end = ATM_qs02_end 
    … …  
    10201049prtFlux ( 'dMass (eau + neige atm) ', dDYN_mas_wat + dATM_mas_sno , 'e', True) 
    10211050 
    1022 echo ( '\n====================================================================================' ) 
    1023 echo ( f'-- SRF changes  -- {Title} ' ) 
    1024  
    1025 if Routing == 'SIMPLE' : 
    1026     RUN_mas_wat_fast_beg   = ONE_stock_int ( d_RUN_beg ['fast_reservoir']   ) 
    1027     RUN_mas_wat_slow_beg   = ONE_stock_int ( d_RUN_beg ['slow_reservoir']   ) 
    1028     RUN_mas_wat_stream_beg = ONE_stock_int ( d_RUN_beg ['stream_reservoir'] ) 
    1029      
    1030     RUN_mas_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  ) 
    1031     RUN_mas_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   ) 
    1032     RUN_mas_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   ) 
    1033      
    1034     RUN_mas_wat_fast_end   = ONE_stock_int ( d_RUN_end ['fast_reservoir']   ) 
    1035     RUN_mas_wat_slow_end   = ONE_stock_int ( d_RUN_end ['slow_reservoir']   ) 
    1036     RUN_mas_wat_stream_end = ONE_stock_int ( d_RUN_end ['stream_reservoir'] ) 
    1037      
    1038     RUN_mas_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  ) 
    1039     RUN_mas_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   ) 
    1040     RUN_mas_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   ) 
    1041  
    1042 if Routing == 'SECHIBA' : 
    1043     RUN_mas_wat_fast_beg   = ONE_stock_int ( d_SRF_beg ['fastres']   ) 
    1044     RUN_mas_wat_slow_beg   = ONE_stock_int ( d_SRF_beg ['slowres']   ) 
    1045     RUN_mas_wat_stream_beg = ONE_stock_int ( d_SRF_beg ['streamres'] ) 
    1046     RUN_mas_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  ) 
    1047     RUN_mas_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   ) 
    1048     RUN_mas_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   ) 
    1049      
    1050     RUN_mas_wat_fast_end   = ONE_stock_int ( d_SRF_end ['fastres']   ) 
    1051     RUN_mas_wat_slow_end   = ONE_stock_int ( d_SRF_end ['slowres']   ) 
    1052     RUN_mas_wat_stream_end = ONE_stock_int ( d_SRF_end ['streamres'] ) 
    1053     RUN_mas_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  ) 
    1054     RUN_mas_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   ) 
    1055     RUN_mas_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   ) 
    1056  
    1057 RUN_mas_wat_beg = Sprec ( [RUN_mas_wat_fast_beg , RUN_mas_wat_slow_beg, RUN_mas_wat_stream_beg, 
    1058                            RUN_mas_wat_flood_beg, RUN_mas_wat_lake_beg, RUN_mas_wat_pond_beg] ) 
    1059                    
    1060 RUN_mas_wat_end = Sprec ( [RUN_mas_wat_fast_end  , RUN_mas_wat_slow_end , RUN_mas_wat_stream_end, 
    1061                            RUN_mas_wat_flood_end , RUN_mas_wat_lake_end , RUN_mas_wat_pond_end] ) 
    1062  
    1063 dRUN_mas_wat_fast   = RUN_mas_wat_fast_end   - RUN_mas_wat_fast_beg 
    1064 dRUN_mas_wat_slow   = RUN_mas_wat_slow_end   - RUN_mas_wat_slow_beg 
    1065 dRUN_mas_wat_stream = RUN_mas_wat_stream_end - RUN_mas_wat_stream_beg 
    1066 dRUN_mas_wat_flood  = RUN_mas_wat_flood_end  - RUN_mas_wat_flood_beg 
    1067 dRUN_mas_wat_lake   = RUN_mas_wat_lake_end   - RUN_mas_wat_lake_beg 
    1068 dRUN_mas_wat_pond   = RUN_mas_wat_pond_end   - RUN_mas_wat_pond_beg 
    1069  
    1070 dRUN_mas_wat        = RUN_mas_wat_end        - RUN_mas_wat_beg 
    1071  
    1072 echo ( f'\nRunoff reservoirs -- {Title} ') 
    1073 echo ( f'------------------------------------------------------------------------------------' ) 
    1074 echo ( f'RUN_mas_wat_fast_beg   = {RUN_mas_wat_fast_beg  :12.6e} kg | RUN_mas_wat_fast_end   = {RUN_mas_wat_fast_end  :12.6e} kg ' ) 
    1075 echo ( f'RUN_mas_wat_slow_beg   = {RUN_mas_wat_slow_beg  :12.6e} kg | RUN_mas_wat_slow_end   = {RUN_mas_wat_slow_end  :12.6e} kg ' ) 
    1076 echo ( f'RUN_mas_wat_stream_beg = {RUN_mas_wat_stream_beg:12.6e} kg | RUN_mas_wat_stream_end = {RUN_mas_wat_stream_end:12.6e} kg ' ) 
    1077 echo ( f'RUN_mas_wat_flood_beg  = {RUN_mas_wat_flood_beg :12.6e} kg | RUN_mas_wat_flood_end  = {RUN_mas_wat_flood_end :12.6e} kg ' ) 
    1078 echo ( f'RUN_mas_wat_lake_beg   = {RUN_mas_wat_lake_beg  :12.6e} kg | RUN_mas_wat_lake_end   = {RUN_mas_wat_lake_end  :12.6e} kg ' ) 
    1079 echo ( f'RUN_mas_wat_pond_beg   = {RUN_mas_wat_pond_beg  :12.6e} kg | RUN_mas_wat_pond_end   = {RUN_mas_wat_pond_end  :12.6e} kg ' ) 
    1080 echo ( f'RUN_mas_wat_beg        = {RUN_mas_wat_beg       :12.6e} kg | RUN_mas_wat_end        = {RUN_mas_wat_end       :12.6e} kg ' ) 
    1081  
    1082 echo ( '------------------------------------------------------------------------------------' ) 
    1083  
    1084 prtFlux ( 'dMass (fast)   ', dRUN_mas_wat_fast  , 'e', True ) 
    1085 prtFlux ( 'dMass (slow)   ', dRUN_mas_wat_slow  , 'e', True ) 
    1086 prtFlux ( 'dMass (stream) ', dRUN_mas_wat_stream, 'e', True ) 
    1087 prtFlux ( 'dMass (flood)  ', dRUN_mas_wat_flood , 'e', True ) 
    1088 prtFlux ( 'dMass (lake)   ', dRUN_mas_wat_lake  , 'e', True ) 
    1089 prtFlux ( 'dMass (pond)   ', dRUN_mas_wat_pond  , 'e', True ) 
    1090 prtFlux ( 'dMass (all)    ', dRUN_mas_wat       , 'e', True ) 
    1091  
    1092 echo ( f'\nWater content in routing  -- {Title} ' ) 
    1093 echo ( '------------------------------------------------------------------------------------' ) 
    1094 echo ( f'RUN_mas_wat_beg = {RUN_mas_wat_end:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg' ) 
    1095 prtFlux ( 'dMass (routing) ', dRUN_mas_wat , 'e', True   ) 
    1096  
    1097 echo ( '\n====================================================================================' ) 
    1098 print (f'Reading SRF restart') 
    1099 SRF_tot_watveg_beg  = d_SRF_beg['tot_watveg_beg']  ; SRF_tot_watveg_beg  = SRF_tot_watveg_beg .where (SRF_tot_watveg_beg  < 1E15, 0.) 
    1100 SRF_tot_watsoil_beg = d_SRF_beg['tot_watsoil_beg'] ; SRF_tot_watsoil_beg = SRF_tot_watsoil_beg.where (SRF_tot_watsoil_beg < 1E15, 0.) 
    1101 SRF_snow_beg        = d_SRF_beg['snow_beg']        ; SRF_snow_beg        = SRF_snow_beg       .where (SRF_snow_beg        < 1E15, 0.) 
    1102 SRF_lakeres_beg     = d_SRF_beg['lakeres']         ; SRF_lakeres_beg     = SRF_lakeres_beg    .where (SRF_lakeres_beg     < 1E15, 0.) 
    1103  
    1104 SRF_tot_watveg_end  = d_SRF_end['tot_watveg_beg']  ; SRF_tot_watveg_end  = SRF_tot_watveg_end .where (SRF_tot_watveg_end  < 1E15, 0.) 
    1105 SRF_tot_watsoil_end = d_SRF_end['tot_watsoil_beg'] ; SRF_tot_watsoil_end = SRF_tot_watsoil_end.where (SRF_tot_watsoil_end < 1E15, 0.) 
    1106 SRF_snow_end        = d_SRF_end['snow_beg']        ; SRF_snow_end        = SRF_snow_end       .where (SRF_snow_end        < 1E15, 0.) 
    1107 SRF_lakeres_end     = d_SRF_end['lakeres']         ; SRF_lakeres_end     = SRF_lakeres_end    .where (SRF_lakeres_end     < 1E15, 0.) 
    1108  
    1109 if LMDZ : 
    1110     SRF_tot_watveg_beg  = lmdz.geo2point (SRF_tot_watveg_beg , dim1D='cell') 
    1111     SRF_tot_watsoil_beg = lmdz.geo2point (SRF_tot_watsoil_beg, dim1D='cell') 
    1112     SRF_snow_beg        = lmdz.geo2point (SRF_snow_beg       , dim1D='cell') 
    1113     SRF_lakeres_beg     = lmdz.geo2point (SRF_lakeres_beg    , dim1D='cell') 
    1114     SRF_tot_watveg_end  = lmdz.geo2point (SRF_tot_watveg_end , dim1D='cell') 
    1115     SRF_tot_watsoil_end = lmdz.geo2point (SRF_tot_watsoil_end, dim1D='cell') 
    1116     SRF_snow_end        = lmdz.geo2point (SRF_snow_end       , dim1D='cell') 
    1117     SRF_lakeres_end     = lmdz.geo2point (SRF_lakeres_end    , dim1D='cell') 
    1118  
    1119  
    1120 # Stock dSoilHum dInterce dSWE dStream dFastR dSlowR dLake dPond dFlood 
    1121  
    1122 SRF_wat_beg = SRF_tot_watveg_beg + SRF_tot_watsoil_beg + SRF_snow_beg 
    1123 SRF_wat_end = SRF_tot_watveg_end + SRF_tot_watsoil_end + SRF_snow_end 
    1124  
    1125 echo ( '\n====================================================================================' ) 
    1126 print ('Computing integrals') 
    1127  
    1128 print ( ' 1/8', end='' ) ; sys.stdout.flush () 
    1129 SRF_mas_watveg_beg   = SRF_stock_int ( SRF_tot_watveg_beg  ) 
    1130 print ( ' 2/8', end='' ) ; sys.stdout.flush () 
    1131 SRF_mas_watsoil_beg  = SRF_stock_int ( SRF_tot_watsoil_beg ) 
    1132 print ( ' 3/8', end='' ) ; sys.stdout.flush () 
    1133 SRF_mas_snow_beg     = SRF_stock_int ( SRF_snow_beg        ) 
    1134 print ( ' 4/8', end='' ) ; sys.stdout.flush () 
    1135 SRF_mas_lake_beg     = ONE_stock_int ( SRF_lakeres_beg     ) 
    1136 print ( ' 5/8', end='' ) ; sys.stdout.flush () 
    1137  
    1138 SRF_mas_watveg_end   = SRF_stock_int ( SRF_tot_watveg_end  ) 
    1139 print ( ' 6/8', end='' ) ; sys.stdout.flush () 
    1140 SRF_mas_watsoil_end  = SRF_stock_int ( SRF_tot_watsoil_end ) 
    1141 print ( ' 7/8', end='' ) ; sys.stdout.flush () 
    1142 SRF_mas_snow_end     = SRF_stock_int ( SRF_snow_end        ) 
    1143 print ( ' 8/8', end='' ) ; sys.stdout.flush () 
    1144 SRF_mas_lake_end     = ONE_stock_int ( SRF_lakeres_end     ) 
    1145  
    1146 print (' -- ') ; sys.stdout.flush () 
    1147  
    1148 dSRF_mas_watveg   = Sprec ( [SRF_mas_watveg_end , -SRF_mas_watveg_beg] ) 
    1149 dSRF_mas_watsoil  = Sprec ( [SRF_mas_watsoil_end, -SRF_mas_watsoil_beg] ) 
    1150 dSRF_mas_snow     = Sprec ( [SRF_mas_snow_end   , -SRF_mas_snow_beg] ) 
    1151 dSRF_mas_lake     = Sprec ( [SRF_mas_lake_end   , -SRF_mas_lake_beg] ) 
    1152  
    1153 echo ( '------------------------------------------------------------------------------------' ) 
    1154 echo ( f'\nSurface reservoirs  -- {Title} ') 
    1155 echo ( f'SRF_mas_watveg_beg   = {SRF_mas_watveg_beg :12.6e} kg | SRF_mas_watveg_end   = {SRF_mas_watveg_end :12.6e} kg ' ) 
    1156 echo ( f'SRF_mas_watsoil_beg  = {SRF_mas_watsoil_beg:12.6e} kg | SRF_mas_watsoil_end  = {SRF_mas_watsoil_end:12.6e} kg ' ) 
    1157 echo ( f'SRF_mas_snow_beg     = {SRF_mas_snow_beg   :12.6e} kg | SRF_mas_snow_end     = {SRF_mas_snow_end   :12.6e} kg ' ) 
    1158 echo ( f'SRF_mas_lake_beg     = {SRF_mas_lake_beg   :12.6e} kg | SRF_mas_lake_end     = {SRF_mas_lake_end   :12.6e} kg ' ) 
    1159  
    1160 prtFlux ( 'dMass (watveg) ', dSRF_mas_watveg , 'e' , True ) 
    1161 prtFlux ( 'dMass (watsoil)', dSRF_mas_watsoil, 'e' , True ) 
    1162 prtFlux ( 'dMass (snow)   ', dSRF_mas_snow   , 'e' , True ) 
    1163 prtFlux ( 'dMass (lake)   ', dSRF_mas_lake   , 'e' , True ) 
    1164  
    1165 SRF_mas_wat_beg = Sprec ( [SRF_mas_watveg_beg , SRF_mas_watsoil_beg, SRF_mas_snow_beg] ) 
    1166 SRF_mas_wat_end = Sprec ( [SRF_mas_watveg_end , SRF_mas_watsoil_end, SRF_mas_snow_end] ) 
    1167  
    1168 dSRF_mas_wat    = Sprec ( [+SRF_mas_watveg_end , +SRF_mas_watsoil_end, +SRF_mas_snow_end,  
    1169                            -SRF_mas_watveg_beg , -SRF_mas_watsoil_beg, -SRF_mas_snow_beg]  ) 
    1170  
    1171 echo ( '------------------------------------------------------------------------------------' ) 
    1172 echo ( f'Water content in surface  -- {Title} ' ) 
    1173 echo ( f'SRF_mas_wat_beg   = {SRF_mas_wat_beg:12.6e} kg | SRF_mas_wat_end  = {SRF_mas_wat_end:12.6e} kg ') 
    1174 prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', True ) 
    1175  
    1176 echo ( '------------------------------------------------------------------------------------' ) 
    1177 echo ( 'Water content in  ATM + SRF + RUN + LAKE' ) 
    1178 echo ( 'mas_wat_beg = {:12.6e} kg | mas_wat_end = {:12.6e} kg '. 
    1179            format (DYN_mas_wat_beg + ATM_mas_sno_beg + RUN_mas_wat_beg + SRF_mas_wat_beg + SRF_mas_lake_beg , 
    1180                    DYN_mas_wat_end + ATM_mas_sno_end + RUN_mas_wat_end + SRF_mas_wat_end + SRF_mas_lake_end ) ) 
    1181 prtFlux ( 'dMass (water atm+srf+run+lake)', dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat + dSRF_mas_lake, 'e', True) 
     1051if SRF : 
     1052    echo ( '\n====================================================================================' ) 
     1053    echo ( f'-- SRF changes  -- {Title} ' ) 
     1054 
     1055    if Routing == 'SIMPLE' : 
     1056        RUN_mas_wat_fast_beg   = ONE_stock_int ( d_RUN_beg ['fast_reservoir']   ) 
     1057        RUN_mas_wat_slow_beg   = ONE_stock_int ( d_RUN_beg ['slow_reservoir']  ) 
     1058        RUN_mas_wat_stream_beg = ONE_stock_int ( d_RUN_beg ['stream_reservoir'] ) 
     1059        RUN_mas_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  ) 
     1060        RUN_mas_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   ) 
     1061        RUN_mas_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   ) 
     1062 
     1063        RUN_mas_wat_fast_end   = ONE_stock_int ( d_RUN_end ['fast_reservoir']   ) 
     1064        RUN_mas_wat_slow_end   = ONE_stock_int ( d_RUN_end ['slow_reservoir']   ) 
     1065        RUN_mas_wat_stream_end = ONE_stock_int ( d_RUN_end ['stream_reservoir'] ) 
     1066 
     1067        RUN_mas_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  ) 
     1068        RUN_mas_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   ) 
     1069        RUN_mas_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   ) 
     1070 
     1071    if Routing == 'SECHIBA' : 
     1072        RUN_mas_wat_fast_beg   = ONE_stock_int ( d_SRF_beg ['fastres']   ) 
     1073        RUN_mas_wat_slow_beg   = ONE_stock_int ( d_SRF_beg ['slowres']   ) 
     1074        RUN_mas_wat_stream_beg = ONE_stock_int ( d_SRF_beg ['streamres'] ) 
     1075        RUN_mas_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  ) 
     1076        RUN_mas_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   ) 
     1077        RUN_mas_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   ) 
     1078 
     1079        RUN_mas_wat_fast_end   = ONE_stock_int ( d_SRF_end ['fastres']   ) 
     1080        RUN_mas_wat_slow_end   = ONE_stock_int ( d_SRF_end ['slowres']   ) 
     1081        RUN_mas_wat_stream_end = ONE_stock_int ( d_SRF_end ['streamres'] ) 
     1082        RUN_mas_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  ) 
     1083        RUN_mas_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   ) 
     1084        RUN_mas_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   ) 
     1085 
     1086    RUN_mas_wat_beg = Sprec ( [RUN_mas_wat_fast_beg , RUN_mas_wat_slow_beg, RUN_mas_wat_stream_beg, 
     1087                            RUN_mas_wat_flood_beg, RUN_mas_wat_lake_beg, RUN_mas_wat_pond_beg] ) 
     1088 
     1089    RUN_mas_wat_end = Sprec ( [RUN_mas_wat_fast_end  , RUN_mas_wat_slow_end , RUN_mas_wat_stream_end, 
     1090                            RUN_mas_wat_flood_end , RUN_mas_wat_lake_end , RUN_mas_wat_pond_end] ) 
     1091 
     1092    dRUN_mas_wat_fast   = RUN_mas_wat_fast_end   - RUN_mas_wat_fast_beg 
     1093    dRUN_mas_wat_slow   = RUN_mas_wat_slow_end   - RUN_mas_wat_slow_beg 
     1094    dRUN_mas_wat_stream = RUN_mas_wat_stream_end - RUN_mas_wat_stream_beg 
     1095    dRUN_mas_wat_flood  = RUN_mas_wat_flood_end  - RUN_mas_wat_flood_beg 
     1096    dRUN_mas_wat_lake   = RUN_mas_wat_lake_end   - RUN_mas_wat_lake_beg 
     1097    dRUN_mas_wat_pond   = RUN_mas_wat_pond_end   - RUN_mas_wat_pond_beg 
     1098 
     1099    dRUN_mas_wat        = RUN_mas_wat_end        - RUN_mas_wat_beg 
     1100 
     1101    echo ( f'\nRunoff reservoirs -- {Title} ') 
     1102    echo (  '------------------------------------------------------------------------------------' ) 
     1103    echo ( f'RUN_mas_wat_fast_beg   = {RUN_mas_wat_fast_beg  :12.6e} kg | RUN_mas_wat_fast_end   = {RUN_mas_wat_fast_end  :12.6e} kg ' ) 
     1104    echo ( f'RUN_mas_wat_slow_beg   = {RUN_mas_wat_slow_beg  :12.6e} kg | RUN_mas_wat_slow_end   = {RUN_mas_wat_slow_end  :12.6e} kg ' ) 
     1105    echo ( f'RUN_mas_wat_stream_beg = {RUN_mas_wat_stream_beg:12.6e} kg | RUN_mas_wat_stream_end = {RUN_mas_wat_stream_end:12.6e} kg ' ) 
     1106    echo ( f'RUN_mas_wat_flood_beg  = {RUN_mas_wat_flood_beg :12.6e} kg | RUN_mas_wat_flood_end  = {RUN_mas_wat_flood_end :12.6e} kg ' ) 
     1107    echo ( f'RUN_mas_wat_lake_beg   = {RUN_mas_wat_lake_beg  :12.6e} kg | RUN_mas_wat_lake_end   = {RUN_mas_wat_lake_end  :12.6e} kg ' ) 
     1108    echo ( f'RUN_mas_wat_pond_beg   = {RUN_mas_wat_pond_beg  :12.6e} kg | RUN_mas_wat_pond_end   = {RUN_mas_wat_pond_end  :12.6e} kg ' ) 
     1109    echo ( f'RUN_mas_wat_beg        = {RUN_mas_wat_beg       :12.6e} kg | RUN_mas_wat_end        = {RUN_mas_wat_end       :12.6e} kg ' ) 
     1110 
     1111    echo ( '------------------------------------------------------------------------------------' ) 
     1112 
     1113    prtFlux ( 'dMass (fast)   ', dRUN_mas_wat_fast  , 'e', True ) 
     1114    prtFlux ( 'dMass (slow)   ', dRUN_mas_wat_slow  , 'e', True ) 
     1115    prtFlux ( 'dMass (stream) ', dRUN_mas_wat_stream, 'e', True ) 
     1116    prtFlux ( 'dMass (flood)  ', dRUN_mas_wat_flood , 'e', True ) 
     1117    prtFlux ( 'dMass (lake)   ', dRUN_mas_wat_lake  , 'e', True ) 
     1118    prtFlux ( 'dMass (pond)   ', dRUN_mas_wat_pond  , 'e', True ) 
     1119    prtFlux ( 'dMass (all)    ', dRUN_mas_wat       , 'e', True ) 
     1120 
     1121    echo ( f'\nWater content in routing  -- {Title} ' ) 
     1122    echo ( '------------------------------------------------------------------------------------' ) 
     1123    echo ( f'RUN_mas_wat_beg = {RUN_mas_wat_end:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg' ) 
     1124    prtFlux ( 'dMass (routing) ', dRUN_mas_wat , 'e', True   ) 
     1125 
     1126    echo ( '\n====================================================================================' ) 
     1127    print ( 'Reading SRF restart') 
     1128    SRF_tot_watveg_beg  = d_SRF_beg['tot_watveg_beg']  ; SRF_tot_watveg_beg  = SRF_tot_watveg_beg .where (SRF_tot_watveg_beg  < 1E15, 0.) 
     1129    SRF_tot_watsoil_beg = d_SRF_beg['tot_watsoil_beg'] ; SRF_tot_watsoil_beg = SRF_tot_watsoil_beg.where (SRF_tot_watsoil_beg < 1E15, 0.) 
     1130    SRF_snow_beg        = d_SRF_beg['snow_beg']        ; SRF_snow_beg        = SRF_snow_beg       .where (SRF_snow_beg        < 1E15, 0.) 
     1131    SRF_lakeres_beg     = d_SRF_beg['lakeres']         ; SRF_lakeres_beg     = SRF_lakeres_beg    .where (SRF_lakeres_beg     < 1E15, 0.) 
     1132 
     1133    SRF_tot_watveg_end  = d_SRF_end['tot_watveg_beg']  ; SRF_tot_watveg_end  = SRF_tot_watveg_end .where (SRF_tot_watveg_end  < 1E15, 0.) 
     1134    SRF_tot_watsoil_end = d_SRF_end['tot_watsoil_beg'] ; SRF_tot_watsoil_end = SRF_tot_watsoil_end.where (SRF_tot_watsoil_end < 1E15, 0.) 
     1135    SRF_snow_end        = d_SRF_end['snow_beg']        ; SRF_snow_end        = SRF_snow_end       .where (SRF_snow_end        < 1E15, 0.) 
     1136    SRF_lakeres_end     = d_SRF_end['lakeres']         ; SRF_lakeres_end     = SRF_lakeres_end    .where (SRF_lakeres_end     < 1E15, 0.) 
     1137 
     1138    if LMDZ : 
     1139        SRF_tot_watveg_beg  = lmdz.geo2point (SRF_tot_watveg_beg , dim1d='cell') 
     1140        SRF_tot_watsoil_beg = lmdz.geo2point (SRF_tot_watsoil_beg, dim1d='cell') 
     1141        SRF_snow_beg        = lmdz.geo2point (SRF_snow_beg       , dim1d='cell') 
     1142        SRF_lakeres_beg     = lmdz.geo2point (SRF_lakeres_beg    , dim1d='cell') 
     1143        SRF_tot_watveg_end  = lmdz.geo2point (SRF_tot_watveg_end , dim1d='cell') 
     1144        SRF_tot_watsoil_end = lmdz.geo2point (SRF_tot_watsoil_end, dim1d='cell') 
     1145        SRF_snow_end        = lmdz.geo2point (SRF_snow_end       , dim1d='cell') 
     1146        SRF_lakeres_end     = lmdz.geo2point (SRF_lakeres_end    , dim1d='cell') 
     1147 
     1148 
     1149    # Stock dSoilHum dInterce dSWE dStream dFastR dSlowR dLake dPond dFlood 
     1150 
     1151    SRF_wat_beg = SRF_tot_watveg_beg + SRF_tot_watsoil_beg + SRF_snow_beg 
     1152    SRF_wat_end = SRF_tot_watveg_end + SRF_tot_watsoil_end + SRF_snow_end 
     1153 
     1154    echo ( '\n====================================================================================' ) 
     1155    print ('Computing integrals') 
     1156 
     1157    print ( ' 1/8', end='' ) ; sys.stdout.flush () 
     1158    SRF_mas_watveg_beg   = SRF_stock_int ( SRF_tot_watveg_beg  ) 
     1159    print ( ' 2/8', end='' ) ; sys.stdout.flush () 
     1160    SRF_mas_watsoil_beg  = SRF_stock_int ( SRF_tot_watsoil_beg ) 
     1161    print ( ' 3/8', end='' ) ; sys.stdout.flush () 
     1162    SRF_mas_snow_beg     = SRF_stock_int ( SRF_snow_beg        ) 
     1163    print ( ' 4/8', end='' ) ; sys.stdout.flush () 
     1164    SRF_mas_lake_beg     = ONE_stock_int ( SRF_lakeres_beg     ) 
     1165    print ( ' 5/8', end='' ) ; sys.stdout.flush () 
     1166 
     1167    SRF_mas_watveg_end   = SRF_stock_int ( SRF_tot_watveg_end  ) 
     1168    print ( ' 6/8', end='' ) ; sys.stdout.flush () 
     1169    SRF_mas_watsoil_end  = SRF_stock_int ( SRF_tot_watsoil_end ) 
     1170    print ( ' 7/8', end='' ) ; sys.stdout.flush () 
     1171    SRF_mas_snow_end     = SRF_stock_int ( SRF_snow_end        ) 
     1172    print ( ' 8/8', end='' ) ; sys.stdout.flush () 
     1173    SRF_mas_lake_end     = ONE_stock_int ( SRF_lakeres_end     ) 
     1174 
     1175    print (' -- ') ; sys.stdout.flush () 
     1176 
     1177    dSRF_mas_watveg   = Sprec ( [SRF_mas_watveg_end , -SRF_mas_watveg_beg] ) 
     1178    dSRF_mas_watsoil  = Sprec ( [SRF_mas_watsoil_end, -SRF_mas_watsoil_beg] ) 
     1179    dSRF_mas_snow     = Sprec ( [SRF_mas_snow_end   , -SRF_mas_snow_beg] ) 
     1180    dSRF_mas_lake     = Sprec ( [SRF_mas_lake_end   , -SRF_mas_lake_beg] ) 
     1181 
     1182    echo ( '------------------------------------------------------------------------------------' ) 
     1183    echo ( f'\nSurface reservoirs  -- {Title} ') 
     1184    echo ( f'SRF_mas_watveg_beg   = {SRF_mas_watveg_beg :12.6e} kg | SRF_mas_watveg_end   = {SRF_mas_watveg_end :12.6e} kg ' ) 
     1185    echo ( f'SRF_mas_watsoil_beg  = {SRF_mas_watsoil_beg:12.6e} kg | SRF_mas_watsoil_end  = {SRF_mas_watsoil_end:12.6e} kg ' ) 
     1186    echo ( f'SRF_mas_snow_beg     = {SRF_mas_snow_beg   :12.6e} kg | SRF_mas_snow_end     = {SRF_mas_snow_end   :12.6e} kg ' ) 
     1187    echo ( f'SRF_mas_lake_beg     = {SRF_mas_lake_beg   :12.6e} kg | SRF_mas_lake_end     = {SRF_mas_lake_end   :12.6e} kg ' ) 
     1188 
     1189    prtFlux ( 'dMass (watveg) ', dSRF_mas_watveg , 'e' , True ) 
     1190    prtFlux ( 'dMass (watsoil)', dSRF_mas_watsoil, 'e' , True ) 
     1191    prtFlux ( 'dMass (snow)   ', dSRF_mas_snow   , 'e' , True ) 
     1192    prtFlux ( 'dMass (lake)   ', dSRF_mas_lake   , 'e' , True ) 
     1193 
     1194    SRF_mas_wat_beg = Sprec ( [SRF_mas_watveg_beg , SRF_mas_watsoil_beg, SRF_mas_snow_beg] ) 
     1195    SRF_mas_wat_end = Sprec ( [SRF_mas_watveg_end , SRF_mas_watsoil_end, SRF_mas_snow_end] ) 
     1196 
     1197    dSRF_mas_wat    = Sprec ( [+SRF_mas_watveg_end , +SRF_mas_watsoil_end, +SRF_mas_snow_end, 
     1198                               -SRF_mas_watveg_beg , -SRF_mas_watsoil_beg, -SRF_mas_snow_beg]  ) 
     1199 
     1200    echo ( '------------------------------------------------------------------------------------' ) 
     1201    echo ( f'Water content in surface  -- {Title} ' ) 
     1202    echo ( f'SRF_mas_wat_beg   = {SRF_mas_wat_beg:12.6e} kg | SRF_mas_wat_end  = {SRF_mas_wat_end:12.6e} kg ') 
     1203    prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', True ) 
     1204 
     1205    echo ( '------------------------------------------------------------------------------------' ) 
     1206    echo ( 'Water content in  ATM + SRF + RUN + LAKE' ) 
     1207    echo ( 'mas_wat_beg = {:12.6e} kg | mas_wat_end = {:12.6e} kg '. 
     1208            format (DYN_mas_wat_beg + ATM_mas_sno_beg + RUN_mas_wat_beg + SRF_mas_wat_beg + SRF_mas_lake_beg , 
     1209                    DYN_mas_wat_end + ATM_mas_sno_end + RUN_mas_wat_end + SRF_mas_wat_end + SRF_mas_lake_end ) ) 
     1210    prtFlux ( 'dMass (water atm+srf+run+lake)', dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat + dSRF_mas_lake, 'e', True) 
    11821211 
    11831212echo ( '\n====================================================================================' ) 
    … …  
    11861215if ATM_HIS == 'latlon' : 
    11871216    echo ( ' latlon case' ) 
    1188     ATM_wbilo_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_oce']), dim1D='cell' ) 
    1189     ATM_wbilo_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_sic']), dim1D='cell' ) 
    1190     ATM_wbilo_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_ter']), dim1D='cell' ) 
    1191     ATM_wbilo_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_lic']), dim1D='cell' ) 
    1192     ATM_runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1D='cell' ) 
    1193     ATM_fqcalving   = lmdz.geo2point ( rprec (d_ATM_his ['fqcalving']), dim1D='cell' ) 
    1194     ATM_fqfonte     = lmdz.geo2point ( rprec (d_ATM_his ['fqfonte']  ), dim1D='cell' ) 
    1195     ATM_precip      = lmdz.geo2point ( rprec (d_ATM_his ['precip']   ), dim1D='cell' ) 
    1196     ATM_snowf       = lmdz.geo2point ( rprec (d_ATM_his ['snow']     ), dim1D='cell' ) 
    1197     ATM_evap        = lmdz.geo2point ( rprec (d_ATM_his ['evap']     ), dim1D='cell' ) 
    1198     ATM_wevap_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_ter']), dim1D='cell' ) 
    1199     ATM_wevap_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_oce']), dim1D='cell' ) 
    1200     ATM_wevap_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_lic']), dim1D='cell' ) 
    1201     ATM_wevap_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_sic']), dim1D='cell' ) 
    1202     ATM_wrain_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_ter']), dim1D='cell' ) 
    1203     ATM_wrain_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_oce']), dim1D='cell' ) 
    1204     ATM_wrain_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_lic']), dim1D='cell' ) 
    1205     ATM_wrain_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_sic']), dim1D='cell' ) 
    1206     ATM_wsnow_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_ter']), dim1D='cell' ) 
    1207     ATM_wsnow_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_oce']), dim1D='cell' ) 
    1208     ATM_wsnow_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_lic']), dim1D='cell' ) 
    1209     ATM_wsnow_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_sic']), dim1D='cell' ) 
    1210     ATM_runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1D='cell' ) 
    1211     echo ( f'End of LATLON case') 
    1212      
     1217    ATM_wbilo_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_oce']), dim1d='cell' ) 
     1218    ATM_wbilo_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_sic']), dim1d='cell' ) 
     1219    ATM_wbilo_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_ter']), dim1d='cell' ) 
     1220    ATM_wbilo_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_lic']), dim1d='cell' ) 
     1221    ATM_runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1d='cell' ) 
     1222    ATM_fqcalving   = lmdz.geo2point ( rprec (d_ATM_his ['fqcalving']), dim1d='cell' ) 
     1223    ATM_fqfonte     = lmdz.geo2point ( rprec (d_ATM_his ['fqfonte']  ), dim1d='cell' ) 
     1224    ATM_precip      = lmdz.geo2point ( rprec (d_ATM_his ['precip']   ), dim1d='cell' ) 
     1225    ATM_snowf       = lmdz.geo2point ( rprec (d_ATM_his ['snow']     ), dim1d='cell' ) 
     1226    ATM_evap        = lmdz.geo2point ( rprec (d_ATM_his ['evap']     ), dim1d='cell' ) 
     1227    ATM_wevap_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_ter']), dim1d='cell' ) 
     1228    ATM_wevap_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_oce']), dim1d='cell' ) 
     1229    ATM_wevap_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_lic']), dim1d='cell' ) 
     1230    ATM_wevap_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_sic']), dim1d='cell' ) 
     1231    ATM_wrain_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_ter']), dim1d='cell' ) 
     1232    ATM_wrain_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_oce']), dim1d='cell' ) 
     1233    ATM_wrain_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_lic']), dim1d='cell' ) 
     1234    ATM_wrain_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_sic']), dim1d='cell' ) 
     1235    ATM_wsnow_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_ter']), dim1d='cell' ) 
     1236    ATM_wsnow_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_oce']), dim1d='cell' ) 
     1237    ATM_wsnow_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_lic']), dim1d='cell' ) 
     1238    ATM_wsnow_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_sic']), dim1d='cell' ) 
     1239    ATM_runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1d='cell' ) 
     1240    echo ( 'End of LATLON case') 
     1241 
    12131242if ATM_HIS == 'ico' : 
    12141243    echo (' ico case') 
    … …  
    12401269    ATM_wsnow_lic   = rprec (d_ATM_his ['wsnow_lic']) 
    12411270    ATM_wsnow_sic   = rprec (d_ATM_his ['wsnow_sic']) 
    1242     echo ( f'End of ico case ') 
     1271    echo ( 'End of ico case ') 
    12431272 
    12441273echo ( 'ATM wprecip_oce' ) 
    … …  
    12681297ATM_wemp_sea    = ATM_wevap_sic - ATM_wprecip_oce 
    12691298 
    1270 if RUN_HIS == 'latlon' : 
    1271     echo ( f'RUN costalflow Grille LATLON' ) 
    1272     if TestInterp : 
    1273          echo ( f'RUN runoff TestInterp' ) 
    1274          RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff_contfrac_interp']  )   , dim1D='cell' ) 
    1275          RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage_contfrac_interp'])   , dim1D='cell' ) 
    1276     else :  
    1277         echo ( f'RUN runoff' ) 
    1278         RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff']         ), dim1D='cell' ) 
    1279         RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage']       ), dim1D='cell' ) 
    1280          
    1281     RUN_coastalflow     = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow']    ), dim1D='cell' )  
    1282     RUN_riverflow       = lmdz.geo2point ( rprec (d_RUN_his ['riverflow']      ), dim1D='cell' ) 
    1283     RUN_riversret       = lmdz.geo2point ( rprec (d_RUN_his ['riversret']      ), dim1D='cell' ) 
    1284     RUN_coastalflow_cpl = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow_cpl']), dim1D='cell' )  
    1285     RUN_riverflow_cpl   = lmdz.geo2point ( rprec (d_RUN_his ['riverflow_cpl']  ), dim1D='cell' ) 
    1286  
    1287 if RUN_HIS == 'ico' : 
    1288     echo ( f'RUN costalflow Grille ICO' ) 
    1289     RUN_coastalflow =  rprec (d_RUN_his ['coastalflow']) 
    1290     RUN_riverflow   =  rprec (d_RUN_his ['riverflow']  ) 
    1291     RUN_runoff      =  rprec (d_RUN_his ['runoff']     ) 
    1292     RUN_drainage    =  rprec (d_RUN_his ['drainage']   ) 
    1293     RUN_riversret   =  rprec (d_RUN_his ['riversret']  ) 
    1294      
    1295     RUN_coastalflow_cpl = rprec (d_RUN_his ['coastalflow_cpl']) 
    1296     RUN_riverflow_cpl   = rprec (d_RUN_his ['riverflow_cpl']  ) 
    1297  
    1298 Step = 0 
    1299  
    1300 if SRF_HIS == 'latlon' : 
    1301     if TestInterp : 
    1302          echo ( f'SRF rain TestInterp' ) 
    1303          SRF_rain     = lmdz.geo2point ( rprec (d_SRF_his ['rain_contfrac_interp'] ), dim1D='cell') 
    1304          SRF_evap     = lmdz.geo2point ( rprec (d_SRF_his ['evap_contfrac_interp'] ), dim1D='cell') 
    1305          SRF_snowf    = lmdz.geo2point ( rprec (d_SRF_his ['snow_contfrac_interp'] ), dim1D='cell') 
    1306          SRF_subli    = lmdz.geo2point ( rprec (d_SRF_his ['subli_contfrac_interp']), dim1D='cell') 
    1307          SRF_transpir = lmdz.geo2point ( rprec (d_SRF_his ['transpir_contfrac_interp']).sum(dim='veget'), dim1D='cell' ) 
    1308          #SRF_rain.attrs.update     ( d_SRF_his ['rain_contfrac_interp'].attrs ) 
    1309          #SRF_evap.attrs.update     ( d_SRF_his ['evap_contfrac_interp'].attrs ) 
    1310          #SRF_snowf.attrs.update    ( d_SRF_his ['snow_contfrac_interp'].attrs ) 
    1311          #SRF_subli.attrs.update    ( d_SRF_his ['subli_contfrac_interp'].attrs ) 
    1312          #SRF_transpir.attrs.update ( d_SRF_his ['transpir_contfrac_interp'].attrs ) 
    1313     else :  
    1314         echo ( f'SRF rain' ) 
    1315         SRF_rain     = lmdz.geo2point ( rprec (d_SRF_his ['rain'] ) , dim1D='cell') 
    1316         SRF_evap     = lmdz.geo2point ( rprec (d_SRF_his ['evap'] ) , dim1D='cell') 
    1317         SRF_snowf    = lmdz.geo2point ( rprec (d_SRF_his ['snowf']) , dim1D='cell') 
    1318         SRF_subli    = lmdz.geo2point ( rprec (d_SRF_his ['subli']) , dim1D='cell') 
    1319         SRF_transpir = lmdz.geo2point ( rprec (d_SRF_his ['transpir']).sum(dim='veget'), dim1D='cell' ) 
    1320  
    1321 if SRF_HIS == 'ico' : 
    1322     echo ( f'SRF rain') 
    1323     SRF_rain     = rprec (d_SRF_his ['rain'] ) 
    1324     SRF_evap     = rprec (d_SRF_his ['evap'] ) 
    1325     SRF_snowf    = rprec (d_SRF_his ['snowf']) 
    1326     SRF_subli    = rprec (d_SRF_his ['subli']) 
    1327     SRF_transpir = rprec (d_SRF_his ['transpir']).sum(dim='veget') 
    1328  
    1329 echo ( f'SRF emp' ) 
    1330 SRF_transpir.attrs['units'] = d_SRF_his ['transpir'].attrs['units'] 
    1331 SRF_emp      = SRF_evap - SRF_rain - SRF_snowf ; SRF_emp.attrs['units'] = SRF_rain.attrs['units'] 
    1332      
    1333 ## Correcting units of SECHIBA variables 
    1334 def mmd2SI ( Var ) : 
    1335     '''Change unit from mm/d or m^3/s to kg/s if needed''' 
    1336     if 'units' in VarT.attrs :  
    1337         if VarT.attrs['units'] in ['m^3/s', 'm3/s', 'm3.s-1',] : 
    1338             VarT.values = VarT.values  * ATM_rho                 ;  VarT.attrs['units'] = 'kg/s' 
    1339         if VarT.attrs['units'] == 'mm/d' : 
    1340             VarT.values = VarT.values  * ATM_rho * (1e-3/86400.) ;  VarT.attrs['units'] = 'kg/s' 
    1341         if VarT.attrs['units'] in ['m^3', 'm3', ] : 
    1342             VarT.values = VarT.values  * ATM_rho                 ;  VarT.attrs['units'] = 'kg' 
    1343              
    1344 for var in [ 'runoff', 'drainage', 'riversret', 'coastalflow', 'riverflow', 'coastalflow_cpl', 'riverflow_cpl' ] : 
    1345     VarT = locals()['RUN_' + var] 
    1346     mmd2SI (VarT) 
    1347  
    1348 for var in ['evap', 'snowf', 'subli', 'transpir', 'rain', 'emp' ] : 
    1349     VarT = locals()['SRF_' + var] 
    1350     mmd2SI (VarT) 
    1351  
    1352 echo ( f'RUN input' ) 
    1353 RUN_input  = RUN_runoff      + RUN_drainage 
    1354 RUN_output = RUN_coastalflow + RUN_riverflow 
    1355  
    1356 echo ( f'ATM flw_wbilo' ) 
     1299if SRF : 
     1300    if RUN_HIS == 'latlon' : 
     1301        echo ( 'RUN costalflow Grille LATLON' ) 
     1302        if TestInterp : 
     1303            echo ( 'RUN runoff TestInterp' ) 
     1304            RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff_contfrac_interp']  )   , dim1d='cell' ) 
     1305            RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage_contfrac_interp'])   , dim1d='cell' ) 
     1306        else : 
     1307            echo ( 'RUN runoff' ) 
     1308            RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff']         ), dim1d='cell' ) 
     1309            RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage']       ), dim1d='cell' ) 
     1310 
     1311        RUN_coastalflow     = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow']    ), dim1d='cell' ) 
     1312        RUN_riverflow       = lmdz.geo2point ( rprec (d_RUN_his ['riverflow']      ), dim1d='cell' ) 
     1313        RUN_riversret       = lmdz.geo2point ( rprec (d_RUN_his ['riversret']      ), dim1d='cell' ) 
     1314        RUN_coastalflow_cpl = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow_cpl']), dim1d='cell' ) 
     1315        RUN_riverflow_cpl   = lmdz.geo2point ( rprec (d_RUN_his ['riverflow_cpl']  ), dim1d='cell' ) 
     1316 
     1317    if RUN_HIS == 'ico' : 
     1318        echo ( 'RUN costalflow Grille ICO' ) 
     1319        RUN_coastalflow =  rprec (d_RUN_his ['coastalflow']) 
     1320        RUN_riverflow   =  rprec (d_RUN_his ['riverflow']  ) 
     1321        RUN_runoff      =  rprec (d_RUN_his ['runoff']     ) 
     1322        RUN_drainage    =  rprec (d_RUN_his ['drainage']   ) 
     1323        RUN_riversret   =  rprec (d_RUN_his ['riversret']  ) 
     1324 
     1325        RUN_coastalflow_cpl = rprec (d_RUN_his ['coastalflow_cpl']) 
     1326        RUN_riverflow_cpl   = rprec (d_RUN_his ['riverflow_cpl']  ) 
     1327 
     1328    Step = 0 
     1329 
     1330    if SRF_HIS == 'latlon' : 
     1331        if TestInterp : 
     1332            echo ( 'SRF rain TestInterp' ) 
     1333            SRF_rain     = lmdz.geo2point ( rprec (d_SRF_his ['rain_contfrac_interp'] ), dim1d='cell') 
     1334            SRF_evap     = lmdz.geo2point ( rprec (d_SRF_his ['evap_contfrac_interp'] ), dim1d='cell') 
     1335            SRF_snowf    = lmdz.geo2point ( rprec (d_SRF_his ['snow_contfrac_interp'] ), dim1d='cell') 
     1336            SRF_subli    = lmdz.geo2point ( rprec (d_SRF_his ['subli_contfrac_interp']), dim1d='cell') 
     1337            SRF_transpir = lmdz.geo2point ( rprec (d_SRF_his ['transpir_contfrac_interp']).sum(dim='veget'), dim1d='cell' ) 
     1338            #SRF_rain.attrs.update     ( d_SRF_his ['rain_contfrac_interp'].attrs ) 
     1339            #SRF_evap.attrs.update     ( d_SRF_his ['evap_contfrac_interp'].attrs ) 
     1340            #SRF_snowf.attrs.update    ( d_SRF_his ['snow_contfrac_interp'].attrs ) 
     1341            #SRF_subli.attrs.update    ( d_SRF_his ['subli_contfrac_interp'].attrs ) 
     1342            #SRF_transpir.attrs.update ( d_SRF_his ['transpir_contfrac_interp'].attrs ) 
     1343        else : 
     1344            echo ( 'SRF rain' ) 
     1345            SRF_rain     = lmdz.geo2point ( rprec (d_SRF_his ['rain'] ) , dim1d='cell') 
     1346            SRF_evap     = lmdz.geo2point ( rprec (d_SRF_his ['evap'] ) , dim1d='cell') 
     1347            SRF_snowf    = lmdz.geo2point ( rprec (d_SRF_his ['snowf']) , dim1d='cell') 
     1348            SRF_subli    = lmdz.geo2point ( rprec (d_SRF_his ['subli']) , dim1d='cell') 
     1349            SRF_transpir = lmdz.geo2point ( rprec (d_SRF_his ['transpir']).sum(dim='veget'), dim1d='cell' ) 
     1350 
     1351    if SRF_HIS == 'ico' : 
     1352        echo ( 'SRF rain') 
     1353        SRF_rain     = rprec (d_SRF_his ['rain'] ) 
     1354        SRF_evap     = rprec (d_SRF_his ['evap'] ) 
     1355        SRF_snowf    = rprec (d_SRF_his ['snowf']) 
     1356        SRF_subli    = rprec (d_SRF_his ['subli']) 
     1357        SRF_transpir = rprec (d_SRF_his ['transpir']).sum(dim='veget') 
     1358 
     1359    echo ( 'SRF emp' ) 
     1360    SRF_transpir.attrs['units'] = d_SRF_his ['transpir'].attrs['units'] 
     1361    SRF_emp      = SRF_evap - SRF_rain - SRF_snowf ; SRF_emp.attrs['units'] = SRF_rain.attrs['units'] 
     1362 
     1363    ## Correcting units of SECHIBA variables 
     1364    def mmd2si ( pvar ) : 
     1365        '''Change unit from mm/d or m^3/s to kg/s if needed''' 
     1366        if 'units' in pvar.attrs : 
     1367            if pvar.attrs['units'] in ['m^3/s', 'm3/s', 'm3.s-1',] : 
     1368                pvar.values = pvar.values  * ATM_RHO                    ;  pvar.attrs['units'] = 'kg/s' 
     1369            if pvar.attrs['units'] == 'mm/d' : 
     1370                pvar.values = pvar.values  * ATM_RHO * (1e-3/lmdz.RDAY) ;  pvar.attrs['units'] = 'kg/s' 
     1371            if pvar.attrs['units'] in ['m^3', 'm3', ] : 
     1372                pvar.values = pvar.values  * ATM_RHO                    ;  pvar.attrs['units'] = 'kg' 
     1373 
     1374    for var in [ 'runoff', 'drainage', 'riversret', 'coastalflow', 'riverflow', 'coastalflow_cpl', 'riverflow_cpl' ] : 
     1375        zvar = locals()['RUN_' + var] 
     1376        mmd2si (zvar) 
     1377 
     1378    for var in ['evap', 'snowf', 'subli', 'transpir', 'rain', 'emp' ] : 
     1379        zvar = locals()['SRF_' + var] 
     1380        mmd2si (zvar) 
     1381 
     1382    echo ( 'RUN input' ) 
     1383    RUN_input  = RUN_runoff      + RUN_drainage 
     1384    RUN_output = RUN_coastalflow + RUN_riverflow 
     1385 
     1386echo ( 'ATM flw_wbilo' ) 
    13571387ATM_flx_wbilo       = ATM_flux_int ( ATM_wbilo      ) 
    13581388ATM_flx_wevap       = ATM_flux_int ( ATM_wevap      ) 
    … …  
    13671397ATM_flx_wbilo_sic   = ATM_flux_int ( ATM_wbilo_sic  ) 
    13681398ATM_flx_wbilo_ter   = ATM_flux_int ( ATM_wbilo_ter  ) 
     1399# Type d'integration a verifier 
    13691400ATM_flx_calving     = ATM_flux_int ( ATM_fqcalving  ) 
    13701401ATM_flx_fqfonte     = ATM_flux_int ( ATM_fqfonte    ) 
    … …  
    13731404LIC_flx_fqfonte     = LIC_flux_int ( ATM_fqfonte    ) 
    13741405 
    1375 echo ( f'ATM flx precip' ) 
     1406echo ( 'ATM flx precip' ) 
    13761407ATM_flx_precip      = ATM_flux_int ( ATM_precip     ) 
    13771408ATM_flx_snowf       = ATM_flux_int ( ATM_snowf      ) 
    … …  
    13841415LIC_flx_runlic      = LIC_flux_int ( ATM_runofflic  ) 
    13851416 
    1386 echo ( f'ATM flx_wrain_ter' ) 
     1417echo ( 'ATM flx_wrain_ter' ) 
    13871418ATM_flx_wrain_ter    = ATM_flux_int ( ATM_wrain_ter ) 
    13881419ATM_flx_wrain_oce    = ATM_flux_int ( ATM_wrain_oce ) 
    … …  
    13971428ATM_flx_wsnow_sea    = ATM_flux_int ( ATM_wsnow_sea ) 
    13981429 
    1399 echo ( f'ATM flx_evap_ter' )  
     1430echo ( 'ATM flx_evap_ter' ) 
    14001431ATM_flx_wevap_ter    = ATM_flux_int ( ATM_wevap_ter ) 
    14011432ATM_flx_wevap_oce    = ATM_flux_int ( ATM_wevap_oce ) 
    … …  
    14161447ATM_flx_emp          = ATM_flux_int ( ATM_emp ) 
    14171448 
    1418 echo ( f'RUN flx_coastal' ) 
    1419 RUN_flx_coastal     = ONE_flux_int ( RUN_coastalflow) 
    1420 echo ( f'RUN flx_river' ) 
    1421 RUN_flx_river       = ONE_flux_int ( RUN_riverflow  ) 
    1422 echo ( f'RUN flx_coastal_cpl' ) 
    1423 RUN_flx_coastal_cpl = ONE_flux_int ( RUN_coastalflow_cpl) 
    1424 echo ( f'RUN flx_river_cpl' ) 
    1425 RUN_flx_river_cpl   = ONE_flux_int ( RUN_riverflow_cpl  ) 
    1426 echo ( f'RUN flx_drainage' ) 
    1427 RUN_flx_drainage    = SRF_flux_int ( RUN_drainage   ) 
    1428 echo ( f'RUN flx_riversset' ) 
    1429 RUN_flx_riversret   = SRF_flux_int ( RUN_riversret  ) 
    1430 echo ( f'RUN flx_runoff' ) 
    1431 RUN_flx_runoff      = SRF_flux_int ( RUN_runoff     ) 
    1432 echo ( f'RUN flx_input' ) 
    1433 RUN_flx_input       = SRF_flux_int ( RUN_input      ) 
    1434 echo ( f'RUN flx_output' ) 
    1435 RUN_flx_output      = ONE_flux_int ( RUN_output     ) 
    1436  
    1437 echo ( f'RUN flx_bil' ) ; Step += 1 
    1438 #RUN_flx_bil    = RUN_flx_input   - RUN_flx_output 
    1439 #RUN_flx_rivcoa = RUN_flx_coastal + RUN_flx_river 
    1440  
    1441 RUN_flx_bil    = ONE_flux_int ( RUN_input       - RUN_output) 
    1442 RUN_flx_rivcoa = ONE_flux_int ( RUN_coastalflow + RUN_riverflow) 
    1443  
    1444 prtFlux ('wbilo_oce            ', ATM_flx_wbilo_oce     , 'f' )           
    1445 prtFlux ('wbilo_sic            ', ATM_flx_wbilo_sic     , 'f' )           
    1446 prtFlux ('wbilo_sic+oce        ', ATM_flx_wbilo_sea     , 'f' )           
    1447 prtFlux ('wbilo_ter            ', ATM_flx_wbilo_ter     , 'f' )           
    1448 prtFlux ('wbilo_lic            ', ATM_flx_wbilo_lic     , 'f' )           
    1449 prtFlux ('Sum wbilo_*          ', ATM_flx_wbilo         , 'f', True)   
    1450 prtFlux ('E-P                  ', ATM_flx_emp           , 'f', True)   
    1451 prtFlux ('calving              ', ATM_flx_calving       , 'f' )  
    1452 prtFlux ('fqfonte              ', ATM_flx_fqfonte       , 'f' )        
    1453 prtFlux ('precip               ', ATM_flx_precip        , 'f' )        
    1454 prtFlux ('snowf                ', ATM_flx_snowf         , 'f' )         
     1449if SRF : 
     1450    echo ( 'RUN flx_coastal' ) 
     1451    RUN_flx_coastal     = ONE_flux_int ( RUN_coastalflow) 
     1452    echo ( 'RUN flx_river' ) 
     1453    RUN_flx_river       = ONE_flux_int ( RUN_riverflow  ) 
     1454    echo ( 'RUN flx_coastal_cpl' ) 
     1455    RUN_flx_coastal_cpl = ONE_flux_int ( RUN_coastalflow_cpl) 
     1456    echo ( 'RUN flx_river_cpl' ) 
     1457    RUN_flx_river_cpl   = ONE_flux_int ( RUN_riverflow_cpl  ) 
     1458    echo ( 'RUN flx_drainage' ) 
     1459    RUN_flx_drainage    = SRF_flux_int ( RUN_drainage   ) 
     1460    echo ( 'RUN flx_riversset' ) 
     1461    RUN_flx_riversret   = SRF_flux_int ( RUN_riversret  ) 
     1462    echo ( 'RUN flx_runoff' ) 
     1463    RUN_flx_runoff      = SRF_flux_int ( RUN_runoff     ) 
     1464    echo ( 'RUN flx_input' ) 
     1465    RUN_flx_input       = SRF_flux_int ( RUN_input      ) 
     1466    echo ( 'RUN flx_output' ) 
     1467    RUN_flx_output      = ONE_flux_int ( RUN_output     ) 
     1468 
     1469    echo ( 'RUN flx_bil' ) ; Step += 1 
     1470    #RUN_flx_bil    = RUN_flx_input   - RUN_flx_output 
     1471    #RUN_flx_rivcoa = RUN_flx_coastal + RUN_flx_river 
     1472 
     1473    RUN_flx_bil    = ONE_flux_int ( RUN_input       - RUN_output) 
     1474    RUN_flx_rivcoa = ONE_flux_int ( RUN_coastalflow + RUN_riverflow) 
     1475 
     1476prtFlux ('wbilo_oce            ', ATM_flx_wbilo_oce     , 'f' ) 
     1477prtFlux ('wbilo_sic            ', ATM_flx_wbilo_sic     , 'f' ) 
     1478prtFlux ('wbilo_sic+oce        ', ATM_flx_wbilo_sea     , 'f' ) 
     1479prtFlux ('wbilo_ter            ', ATM_flx_wbilo_ter     , 'f' ) 
     1480prtFlux ('wbilo_lic            ', ATM_flx_wbilo_lic     , 'f' ) 
     1481prtFlux ('Sum wbilo_*          ', ATM_flx_wbilo         , 'f', True) 
     1482prtFlux ('E-P                  ', ATM_flx_emp           , 'f', True) 
     1483prtFlux ('calving              ', ATM_flx_calving       , 'f' ) 
     1484prtFlux ('fqfonte              ', ATM_flx_fqfonte       , 'f' ) 
     1485prtFlux ('precip               ', ATM_flx_precip        , 'f' ) 
     1486prtFlux ('snowf                ', ATM_flx_snowf         , 'f' ) 
    14551487prtFlux ('evap                 ', ATM_flx_evap          , 'f' ) 
    14561488prtFlux ('runoff lic           ', ATM_flx_runlic        , 'f' ) 
    … …  
    14671499prtFlux ('ERROR emp            ', ATM_flx_wemp    - ATM_flx_emp   , 'e', True ) 
    14681500 
    1469  
    1470 echo ( '\n====================================================================================' ) 
    1471 echo ( f'-- RUNOFF Fluxes  -- {Title} ' ) 
    1472 prtFlux ('coastalflow   ', RUN_flx_coastal    , 'f' )  
    1473 prtFlux ('riverflow     ', RUN_flx_river      , 'f' )         
    1474 prtFlux ('coastal_cpl   ', RUN_flx_coastal_cpl, 'f' )   
    1475 prtFlux ('riverf_cpl    ', RUN_flx_river_cpl  , 'f' )     
    1476 prtFlux ('river+coastal ', RUN_flx_rivcoa     , 'f' )    
    1477 prtFlux ('drainage      ', RUN_flx_drainage   , 'f' )    
    1478 prtFlux ('riversret     ', RUN_flx_riversret  , 'f' )    
    1479 prtFlux ('runoff        ', RUN_flx_runoff     , 'f' )    
    1480 prtFlux ('river in      ', RUN_flx_input      , 'f' )    
    1481 prtFlux ('river out     ', RUN_flx_output     , 'f' )    
    1482 prtFlux ('river bil     ', RUN_flx_bil        , 'f' )    
     1501if SRF : 
     1502    echo ( '\n====================================================================================' ) 
     1503    echo ( f'-- RUNOFF Fluxes  -- {Title} ' ) 
     1504    prtFlux ('coastalflow   ', RUN_flx_coastal    , 'f' ) 
     1505    prtFlux ('riverflow     ', RUN_flx_river      , 'f' ) 
     1506    prtFlux ('coastal_cpl   ', RUN_flx_coastal_cpl, 'f' ) 
     1507    prtFlux ('riverf_cpl    ', RUN_flx_river_cpl  , 'f' ) 
     1508    prtFlux ('river+coastal ', RUN_flx_rivcoa     , 'f' ) 
     1509    prtFlux ('drainage      ', RUN_flx_drainage   , 'f' ) 
     1510    prtFlux ('riversret     ', RUN_flx_riversret  , 'f' ) 
     1511    prtFlux ('runoff        ', RUN_flx_runoff     , 'f' ) 
     1512    prtFlux ('river in      ', RUN_flx_input      , 'f' ) 
     1513    prtFlux ('river out     ', RUN_flx_output     , 'f' ) 
     1514    prtFlux ('river bil     ', RUN_flx_bil        , 'f' ) 
    14831515 
    14841516ATM_flx_budget = -ATM_flx_wbilo + ATM_flx_calving + ATM_flx_runlic #+ ATM_flx_fqfonte + RUN_flx_river 
    … …  
    15421574echo ( 'LIC error (-wbilo_lic - runofflic*frac_lic)                   = {:12.4e} (rel) '.format ( (LIC_flx_budget3-dLIC_mas_wat)/dLIC_mas_wat) ) 
    15431575 
    1544 echo ( '\n====================================================================================' ) 
    1545 echo ( f'-- SECHIBA fluxes  -- {Title} ' ) 
    1546  
    1547 SRF_flx_rain     = SRF_flux_int ( SRF_rain     ) 
    1548 SRF_flx_evap     = SRF_flux_int ( SRF_evap     ) 
    1549 SRF_flx_snowf    = SRF_flux_int ( SRF_snowf    ) 
    1550 SRF_flx_subli    = SRF_flux_int ( SRF_subli    ) 
    1551 SRF_flx_transpir = SRF_flux_int ( SRF_transpir ) 
    1552 SRF_flx_emp      = SRF_flux_int ( SRF_emp      ) 
    1553  
    1554 RUN_flx_torouting   = SRF_flux_int  ( RUN_runoff       + RUN_drainage) 
    1555 RUN_flx_fromrouting = ONE_flux_int  ( RUN_coastalflow + RUN_riverflow ) 
    1556  
    1557 SRF_flx_all =  SRF_flux_int  ( SRF_rain + SRF_snowf - SRF_evap - RUN_runoff - RUN_drainage ) 
    1558  
    1559 prtFlux ('rain         ', SRF_flx_rain       , 'f' ) 
    1560 prtFlux ('evap         ', SRF_flx_evap       , 'f' ) 
    1561 prtFlux ('snowf        ', SRF_flx_snowf      , 'f' ) 
    1562 prtFlux ('E-P          ', SRF_flx_emp        , 'f' ) 
    1563 prtFlux ('subli        ', SRF_flx_subli      , 'f' ) 
    1564 prtFlux ('transpir     ', SRF_flx_transpir   , 'f' ) 
    1565 prtFlux ('to routing   ', RUN_flx_torouting  , 'f' ) 
    1566 prtFlux ('budget       ', SRF_flx_all        , 'f', small=True ) 
    1567  
    1568 echo ( '\n------------------------------------------------------------------------------------' ) 
    1569 echo ( 'Water content in surface ' ) 
    1570 echo ( f'SRF_mas_wat_beg  = {SRF_mas_wat_beg:12.6e} kg | SRF_mas_wat_end  = {SRF_mas_wat_end:12.6e} kg ' ) 
    1571 prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', small=True) 
    1572 prtFlux ( 'Error            ',  SRF_flx_all-dSRF_mas_wat, 'e', small=True ) 
    1573 echo ( 'dMass (water srf) = {:12.4e} (rel)   '.format ( (SRF_flx_all-dSRF_mas_wat)/dSRF_mas_wat) ) 
    1574  
    1575 echo ( '\n====================================================================================' ) 
    1576 echo ( f'-- Check ATM vs. SRF -- {Title} ' ) 
    1577 prtFlux ('E-P ATM       ', ATM_flx_wemp_ter   , 'f' ) 
    1578 prtFlux ('wbilo ter     ', ATM_flx_wbilo_ter  , 'f' ) 
    1579 prtFlux ('E-P SRF       ', SRF_flx_emp        , 'f' ) 
    1580 prtFlux ('SRF/ATM error ', ATM_flx_wbilo_ter - SRF_flx_emp, 'e', True) 
    1581 echo ( 'SRF/ATM error {:12.3e} (rel)  '.format ( (ATM_flx_wbilo_ter - SRF_flx_emp)/SRF_flx_emp ) ) 
    1582  
    1583 echo ('') 
    1584 echo ( '\n====================================================================================' ) 
    1585 echo ( f'-- RUNOFF fluxes  -- {Title} ' ) 
    1586 RUN_flx_all = RUN_flx_torouting - RUN_flx_river - RUN_flx_coastal 
    1587 prtFlux ('runoff    ', RUN_flx_runoff      , 'f' )  
    1588 prtFlux ('drainage  ', RUN_flx_drainage    , 'f' )  
    1589 prtFlux ('run+drain ', RUN_flx_torouting   , 'f' )  
    1590 prtFlux ('river     ', RUN_flx_river       , 'f' )  
    1591 prtFlux ('coastal   ', RUN_flx_coastal     , 'f' )  
    1592 prtFlux ('riv+coa   ', RUN_flx_fromrouting , 'f' )  
    1593 prtFlux ('budget    ', RUN_flx_all         , 'f' , small=True) 
    1594  
    1595 echo ( '\n------------------------------------------------------------------------------------' ) 
    1596 echo ( f'Water content in routing+lake  -- {Title} ' ) 
    1597 echo ( f'RUN_mas_wat_beg  = {RUN_mas_wat_beg:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg ' ) 
    1598 prtFlux ( 'dMass (routing)  ', dRUN_mas_wat+dSRF_mas_lake, 'f', small=True) 
    1599 prtFlux ( 'Routing error    ', RUN_flx_all+dSRF_mas_lake-dRUN_mas_wat, 'e', small=True ) 
    1600 echo ( 'Routing error : {:12.3e} (rel)'.format ( (RUN_flx_all-dSRF_mas_lake-dRUN_mas_wat)/(dRUN_mas_wat+dSRF_mas_lake) ) ) 
    1601  
    1602 echo ( '\n------------------------------------------------------------------------------------' ) 
    1603 echo ( f'Water content in routing  -- {Title} ' ) 
    1604 echo ( f'RUN_mas_wat_beg  = {RUN_mas_wat_beg:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg ' ) 
    1605 prtFlux ( 'dMass (routing) ', dRUN_mas_wat, 'f', small=True  ) 
    1606 prtFlux ( 'Routing error   ', RUN_flx_all-dRUN_mas_wat, 'e', small=True) 
    1607 echo ( 'Routing error : {:12.3e} (rel)'.format ( (RUN_flx_all-dRUN_mas_wat)/dRUN_mas_wat ) ) 
     1576if SRF : 
     1577    echo ( '\n====================================================================================' ) 
     1578    echo ( f'-- SECHIBA fluxes  -- {Title} ' ) 
     1579 
     1580    SRF_flx_rain     = SRF_flux_int ( SRF_rain     ) 
     1581    SRF_flx_evap     = SRF_flux_int ( SRF_evap     ) 
     1582    SRF_flx_snowf    = SRF_flux_int ( SRF_snowf    ) 
     1583    SRF_flx_subli    = SRF_flux_int ( SRF_subli    ) 
     1584    SRF_flx_transpir = SRF_flux_int ( SRF_transpir ) 
     1585    SRF_flx_emp      = SRF_flux_int ( SRF_emp      ) 
     1586 
     1587    RUN_flx_torouting   = SRF_flux_int  ( RUN_runoff       + RUN_drainage) 
     1588    RUN_flx_fromrouting = ONE_flux_int  ( RUN_coastalflow + RUN_riverflow ) 
     1589 
     1590    SRF_flx_all =  SRF_flux_int  ( SRF_rain + SRF_snowf - SRF_evap - RUN_runoff - RUN_drainage ) 
     1591 
     1592    prtFlux ('rain         ', SRF_flx_rain       , 'f' ) 
     1593    prtFlux ('evap         ', SRF_flx_evap       , 'f' ) 
     1594    prtFlux ('snowf        ', SRF_flx_snowf      , 'f' ) 
     1595    prtFlux ('E-P          ', SRF_flx_emp        , 'f' ) 
     1596    prtFlux ('subli        ', SRF_flx_subli      , 'f' ) 
     1597    prtFlux ('transpir     ', SRF_flx_transpir   , 'f' ) 
     1598    prtFlux ('to routing   ', RUN_flx_torouting  , 'f' ) 
     1599    prtFlux ('budget       ', SRF_flx_all        , 'f', small=True ) 
     1600 
     1601    echo ( '\n------------------------------------------------------------------------------------' ) 
     1602    echo ( 'Water content in surface ' ) 
     1603    echo ( f'SRF_mas_wat_beg  = {SRF_mas_wat_beg:12.6e} kg | SRF_mas_wat_end  = {SRF_mas_wat_end:12.6e} kg ' ) 
     1604    prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', small=True) 
     1605    prtFlux ( 'Error            ',  SRF_flx_all-dSRF_mas_wat, 'e', small=True ) 
     1606    echo ( 'dMass (water srf) = {:12.4e} (rel)   '.format ( (SRF_flx_all-dSRF_mas_wat)/dSRF_mas_wat) ) 
     1607 
     1608    echo ( '\n====================================================================================' ) 
     1609    echo ( f'-- Check ATM vs. SRF -- {Title} ' ) 
     1610    prtFlux ('E-P ATM       ', ATM_flx_wemp_ter   , 'f' ) 
     1611    prtFlux ('wbilo ter     ', ATM_flx_wbilo_ter  , 'f' ) 
     1612    prtFlux ('E-P SRF       ', SRF_flx_emp        , 'f' ) 
     1613    prtFlux ('SRF/ATM error ', ATM_flx_wbilo_ter - SRF_flx_emp, 'e', True) 
     1614    echo ( 'SRF/ATM error {:12.3e} (rel)  '.format ( (ATM_flx_wbilo_ter - SRF_flx_emp)/SRF_flx_emp ) ) 
     1615 
     1616    echo ('') 
     1617    echo ( '\n====================================================================================' ) 
     1618    echo ( f'-- RUNOFF fluxes  -- {Title} ' ) 
     1619    RUN_flx_all = RUN_flx_torouting - RUN_flx_river - RUN_flx_coastal 
     1620    prtFlux ('runoff    ', RUN_flx_runoff      , 'f' ) 
     1621    prtFlux ('drainage  ', RUN_flx_drainage    , 'f' ) 
     1622    prtFlux ('run+drain ', RUN_flx_torouting   , 'f' ) 
     1623    prtFlux ('river     ', RUN_flx_river       , 'f' ) 
     1624    prtFlux ('coastal   ', RUN_flx_coastal     , 'f' ) 
     1625    prtFlux ('riv+coa   ', RUN_flx_fromrouting , 'f' ) 
     1626    prtFlux ('budget    ', RUN_flx_all         , 'f' , small=True) 
     1627 
     1628    echo ( '\n------------------------------------------------------------------------------------' ) 
     1629    echo ( f'Water content in routing+lake  -- {Title} ' ) 
     1630    echo ( f'RUN_mas_wat_beg  = {RUN_mas_wat_beg:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg ' ) 
     1631    prtFlux ( 'dMass (routing)  ', dRUN_mas_wat+dSRF_mas_lake, 'f', small=True) 
     1632    prtFlux ( 'Routing error    ', RUN_flx_all+dSRF_mas_lake-dRUN_mas_wat, 'e', small=True ) 
     1633    echo ( 'Routing error : {:12.3e} (rel)'.format ( (RUN_flx_all-dSRF_mas_lake-dRUN_mas_wat)/(dRUN_mas_wat+dSRF_mas_lake) ) ) 
     1634 
     1635    echo ( '\n------------------------------------------------------------------------------------' ) 
     1636    echo ( f'Water content in routing  -- {Title} ' ) 
     1637    echo ( f'RUN_mas_wat_beg  = {RUN_mas_wat_beg:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg ' ) 
     1638    prtFlux ( 'dMass (routing) ', dRUN_mas_wat, 'f', small=True  ) 
     1639    prtFlux ( 'Routing error   ', RUN_flx_all-dRUN_mas_wat, 'e', small=True) 
     1640    echo ( 'Routing error : {:12.3e} (rel)'.format ( (RUN_flx_all-dRUN_mas_wat)/dRUN_mas_wat ) ) 
    16081641 
    16091642echo ( ' ' ) 
    … …  
    16111644 
    16121645echo ( 'SVN Information' ) 
    1613 for clef in SVN.keys () : echo ( SVN[clef] ) 
     1646for kk in SVN.keys(): 
     1647    print ( SVN[kk] ) 
Note: See TracChangeset for help on using the changeset viewer.