Changeset 6277 for TOOLS


Ignore:
Timestamp:
12/08/22 10:24:05 (17 months ago)
Author:
omamce
Message:

O.M. : WATER_BUDGET

Switch to .ini file to read configuratio
Correct SECHIBA budget

Location:
TOOLS/WATER_BUDGET
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • TOOLS/WATER_BUDGET/ATM_waterbudget.py

    r6271 r6277  
    2121### 
    2222## Import system modules 
    23 import sys, os, shutil, subprocess 
     23import sys, os, shutil, subprocess, platform 
    2424import numpy as np 
    2525import configparser, re 
    … …  
    2929config.optionxform = str # To keep capitals 
    3030 
    31 config['Files']  = {} 
    32 config['System'] = {} 
    33  
    34 ##-- Some physical constants 
    35 #-- Earth Radius 
    36 Ra = 6366197.7236758135 
    37 #-- Gravity 
    38 Grav = 9.81 
    39 #-- Ice volumic mass (kg/m3) in LIM3 
    40 ICE_rho_ice = 917.0 
    41 #-- Snow volumic mass (kg/m3) in LIM3 
    42 ICE_rho_sno = 330.0 
    43 #-- Ocean water volumic mass (kg/m3) in NEMO 
    44 OCE_rho_liq = 1026. 
    45 #-- Water volumic mass in atmosphere 
    46 ATM_rho = 1.0e3 
    47 #-- Water volumic mass in surface reservoirs 
    48 SRF_rho = 1.0e3 
    49 #-- Water volumic mass of rivers 
    50 RUN_rho = 1.0e3 
    51  
    5231## Read experiment parameters 
    53 ATM = None ; ORCA = None ; NEMO = None ; OCE_relax = False ; OCE_icb = False ; Coupled = False ; Routing = None 
     32ATM=None ; ORCA=None ; NEMO=None ; OCE_relax=False ; OCE_icb=False ; Coupled=False ; Routing=None ;  TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=None 
     33 
     34## 
     35ARCHIVE =None ; STORAGE = None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None 
     36TmpDir=None ; FileOut=None ; dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None ; FileCommon=None 
     37file_ATM_his=None ; file_SRF_his=None ; file_RUN_his=None ; file_OCE_his=None ;  file_ICE_his=None ; file_OCE_sca=None 
     38file_restart_beg=None ; file_restart_end=None ; file_ATM_beg=None ; file_ATM_end=None ; file_DYN_beg=None ; file_DYN_end=None ; file_SRF_beg=None ; file_SRF_end=None 
     39file_RUN_beg=None ; file_RUN_end=None ; file_RUN_end=None ; file_OCE_beg=None ; file_ICE_beg=None ; file_OCE_beg=None ; file_OCE_end=None 
     40 
     41d_ATM_his=None ; d_SRF_his=None ; d_RUN_his=None ; d_OCE_his=None ;  d_ICE_his=None ; d_OCE_sca=None 
     42d_restart_beg=None ; d_restart_end=None ; d_ATM_beg=None ; d_ATM_end=None ; d_DYN_beg=None ; d_DYN_end=None ; d_SRF_beg=None ; d_SRF_end=None 
     43d_RUN_beg=None ; d_RUN_end=None ; d_RUN_end=None ; d_OCE_beg=None ; d_ICE_beg=None ; d_OCE_beg=None ; d_OCE_end=None 
    5444 
    5545# Arguments passed 
    5646print ( "Name of Python script:", sys.argv[0] ) 
    57 IniFile =  sys.argv[1] 
     47IniFile = sys.argv[1] 
    5848print ("Input file : ", IniFile ) 
    5949config.read (IniFile) 
     50FullIniFile = 'full_' + IniFile 
    6051 
    6152def setBool (chars) : 
    … …  
    7970    return setNum 
    8071 
    81 print ('[Experiment]') 
    82 for VarName in config['Experiment'].keys() : 
    83     locals()[VarName] = config['Experiment'][VarName] 
    84     locals()[VarName] = setBool (locals()[VarName]) 
    85     locals()[VarName] = setNum  (locals()[VarName]) 
    86     print ( '{:25} set to : {:}'.format (VarName, locals()[VarName]) ) 
    87  
    88 # ### 
     72def setNone (chars) : 
     73    '''Convert specific char string to None if possible''' 
     74    if type (chars) == str : 
     75        if chars in ['None', 'NONE', 'none'] : setNone = None 
     76        else : setNone = chars 
     77    else : setNone = chars 
     78    return setNone 
     79 
     80## Reading config 
     81for Section in ['Experiment', 'libIGCM', 'Files', 'Physics' ] : 
     82    if Section in config.keys() :  
     83        print ( f'[{Section}]' ) 
     84        for VarName in config[Section].keys() : 
     85            locals()[VarName] = config[Section][VarName] 
     86            locals()[VarName] = setBool (locals()[VarName]) 
     87            locals()[VarName] = setNum  (locals()[VarName]) 
     88            print ( '{:25} set to : {:}'.format (VarName, locals()[VarName]) ) 
     89 
     90if not 'Files' in config.keys() : config['Files'] = {} 
     91 
     92def unDefined (char) : 
     93    if char in globals () : 
     94        if char == None : return True 
     95        else : return False 
     96    else : return True 
     97 
     98##-- Some physical constants 
     99#-- Earth Radius 
     100if not 'Ra'            in locals () : Ra = 6366197.7236758135 
     101#-- Gravity 
     102if not 'Grav'          in locals () : Grav = 9.81 
     103#-- Ice volumic mass (kg/m3) in LIM3 
     104if not 'ICE_rho_ice'   in locals () : ICE_rho_ice = 917.0 
     105#-- Snow volumic mass (kg/m3) in LIM3 
     106if not 'ICE_rho_sno'   in locals () : ICE_rho_sno = 330.0 
     107#-- Ocean water volumic mass (kg/m3) in NEMO 
     108if not 'OCE_rho_liq'   in locals () : OCE_rho_liq = 1026. 
     109#-- Water volumic mass in atmosphere 
     110if not 'ATM_rho'       in locals () : ATM_rho = 1000. 
     111#-- Water volumic mass in surface reservoirs 
     112if not 'SRF_rho'       in locals () : SRF_rho = 1000. 
     113#-- Water volumic mass of rivers 
     114if not 'RUN_rho'       in locals () : RUN_rho = 1000. 
     115#-- Year length 
     116if not 'YearLength'    in locals () : YearLength = 365.25 * 24. * 60. * 60. 
     117 
     118config['Physics'] = { 'Ra':Ra, 'Grav':Grav, 'ICE_rho_ice':ICE_rho_ice, 'ICE_rho_sno':ICE_rho_sno, 'OCE_rho_liq':OCE_rho_liq, 'ATM_rho':ATM_rho, 'SRF_rho':SRF_rho, 'RUN_rho':RUN_rho} 
     119 
     120## -- 
    89121ICO  = ( 'ICO' in ATM ) 
    90122LMDZ = ( 'LMD' in ATM ) 
    91123     
    92 ### 
    93 ## Import system modules 
    94 import sys, os, shutil, subprocess 
    95 import configparser, re 
    96  
    97 config = configparser.ConfigParser() 
    98 config['Files'] = {} 
    99124 
    100125# Where do we run ? 
    101 SysName, NodeName, Release, Version, Machine = os.uname() 
     126SysName, NodeName, Release, Version, Machine = os.uname () 
    102127TGCC  = ( 'irene'   in NodeName ) 
    103128IDRIS = ( 'jeanzay' in NodeName ) 
    … …  
    108133    if "Intel(R) Xeon(R) Platinum" in CPU : Machine = 'irene' 
    109134    if "AMD"                       in CPU : Machine = 'irene-amd' 
    110          
    111     ARCHIVE     = subprocess.getoutput ( f'ccc_home --cccstore   -d {Project} -u {User}' ) 
    112     STORAGE     = subprocess.getoutput ( f'ccc_home --cccwork    -d {Project} -u {User}' ) 
    113     SCRATCHDIR  = subprocess.getoutput ( f'ccc_home --cccscratch -d {Project} -u {User}' ) 
    114     R_IN        = os.path.join ( subprocess.getoutput ( f'ccc_home --cccwork -d igcmg -u igcmg' ), 'IGCM') 
    115     rebuild     = os.path.join ( subprocess.getoutput ( f'ccc_home --ccchome -d igcmg -u igcmg' ), 'Tools', Machine, 'rebuild_nemo', 'bin', 'rebuild_nemo' ) 
     135 
     136    if libIGCM :  
     137        if ARCHIVE    == None : ARCHIVE     = subprocess.getoutput ( f'ccc_home --cccstore   -d {Project} -u {User}' ) 
     138        if STORAGE    == None : STORAGE     = subprocess.getoutput ( f'ccc_home --cccwork    -d {Project} -u {User}' ) 
     139        if SCRATCHDIR == None : SCRATCHDIR  = subprocess.getoutput ( f'ccc_home --cccscratch -d {Project} -u {User}' ) 
     140        if R_IN       == None : R_IN        = os.path.join ( subprocess.getoutput ( f'ccc_home --cccwork -d igcmg -u igcmg' ), 'IGCM') 
     141        if rebuild    == None : rebuild     = os.path.join ( subprocess.getoutput ( f'ccc_home --ccchome -d igcmg -u igcmg' ), 'Tools', Machine, 'rebuild_nemo', 'bin', 'rebuild_nemo' ) 
    116142 
    117143    ## Specific to run at TGCC. 
    … …  
    121147         
    122148    ## Creates output directory 
    123     #TmpDir = os.path.join ( subprocess.getoutput ( 'ccc_home --cccscratch' ), f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) 
    124     TmpDir = os.path.join ( '/ccc/scratch/cont003/drf/p86mart', f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) 
    125  
     149    if TmpDir == None : 
     150        TmpDir = os.path.join ( subprocess.getoutput ( 'ccc_home --cccscratch' ), f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) 
     151        config['Files']['TmpDir'] = TmpDir 
     152         
    126153if IDRIS : 
    127154    raise Exception ("Pour IDRIS : repertoires et chemins a definir")  
     155 
     156config['System'] = {'SysName':SysName, 'NodeName':NodeName, 'Release':Release, 'Version':Version,'Machine':Machine, 'TGCC':TGCC,'IDRIS':IDRIS, 'CPU':CPU, 
     157                    'Program'  : "Generated by : " + str(sys.argv),  
     158                    'HOSTNAME' : platform.node (), 'LOGNAME'  : os.getlogin (), 
     159                    'Python'   : f'{platform.python_version ()}', 
     160                    'OS'       : f'{platform.system ()} release : {platform.release ()} hardware : {platform.machine ()}', 
     161                    'SVN_Author'   : "$Author$",  
     162                    'SVN_Date'     : "$Date$", 
     163                    'SVN_Revision' : "$Revision$", 
     164                    'SVN_Id'       : "$Id$", 
     165                    'SVN_HeadURL'  : "$HeadURL$"} 
     166 
     167if libIGCM : 
     168    config['libIGCM'] = {'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'SCRATCHDIR':SCRATCHDIR, 'R_IN':R_IN, 'rebuild':rebuild } 
    128169 
    129170## Import specific module 
    … …  
    132173import xarray as xr 
    133174 
    134 config['Files'][TmpDir] = TmpDir 
    135  
    136175# Output file 
    137 FileOut = f'ATM_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out' 
     176if FileOut == None :  
     177    FileOut = f'ATM_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out' 
     178    config['Files']['FileOut'] = FileOut 
     179 
    138180f_out = open ( FileOut, mode = 'w' ) 
    139181 
    140182# Function to print to stdout *and* output file 
    141183def echo (string, end='\n') : 
    142     print ( string, end=end  ) 
     184    print ( str(string), end=end  ) 
    143185    sys.stdout.flush () 
    144     f_out.write ( string + end ) 
     186    f_out.write ( str(string) + end ) 
    145187    f_out.flush () 
    146188    return None 
    … …  
    165207if not os.path.exists (TmpDirICE) : os.mkdir (TmpDirICE ) 
    166208 
     209echo (' ') 
     210echo ( f'JobName : {JobName}' ) 
     211echo (Comment) 
    167212echo ( f'Working in TMPDIR : {TmpDir}' ) 
    168213 
    … …  
    172217if Freq == "MO" : FreqDir = os.path.join ('Output' , 'MO' ) 
    173218if Freq == "SE" : FreqDir = os.path.join ('Analyse', 'SE' ) 
    174 dir_ATM_his = os.path.join ( R_SAVE, "ATM", FreqDir ) 
    175 dir_SRF_his = os.path.join ( R_SAVE, "SRF", FreqDir ) 
    176  
     219if dir_ATM_his == None : 
     220    dir_ATM_his = os.path.join ( R_SAVE, "ATM", FreqDir ) 
     221    config['Files']['dir_ATM_his'] = dir_ATM_his 
     222if dir_SRF_his == None :  
     223    dir_SRF_his = os.path.join ( R_SAVE, "SRF", FreqDir ) 
     224    config['Files']['dir_SRF_his'] = dir_SRF_his 
     225     
    177226echo ( f'The analysis relies on files from the following model output directories : ' ) 
    178227echo ( f'{dir_ATM_his}' ) 
    … …  
    180229 
    181230#-- Files Names 
    182 if Freq == 'MO' : FileCommon = f'{JobName}_{YearBegin}0101_{YearEnd}1231_1M' 
    183 if Freq == 'SE' : FileCommon = f'{JobName}_SE_{YearBegin}0101_{YearEnd}1231_1M' 
    184  
     231if Period == None : 
     232    if Freq == 'MO' : Period = f'{YearBegin}0101_{YearEnd}1231_1M' 
     233    if Freq == 'SE' : Period = f'SE_{YearBegin}0101_{YearEnd}1231_1M' 
     234    config['Files']['Period'] = Period 
     235if FileCommon == None : 
     236    FileCommon = f'{JobName}_{Period}' 
     237    config['Files']['FileCommon'] = FileCommon 
     238 
     239if Title == None : 
     240    Title = f'{JobName} : {Freq} : {YearBegin}-01-01 - {YearEnd}-12-31' 
     241    config['Files']['Title'] = Title 
     242      
    185243echo ('\nOpen history files' ) 
    186 file_ATM_his = os.path.join (  dir_ATM_his, f'{FileCommon}_histmth.nc' ) 
    187 file_SRF_his = os.path.join (  dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 
    188 if Routing == 'ORCHIDEE'  : 
    189     file_RUN_his = os.path.join (  dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 
     244if file_ATM_his == None :  
     245    file_ATM_his = os.path.join (  dir_ATM_his, f'{FileCommon}_histmth.nc' ) 
     246    config['Files']['file_ATM_his'] = file_ATM_his 
     247if file_SRF_his == None : 
     248    file_SRF_his = os.path.join (  dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 
     249    config['Files']['file_SRF_his'] = file_SRF_his 
     250#if Routing == 'SECHIBA'  : 
     251#     file_RUN_his = os.path.join (  dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 
    190252if Routing == 'SIMPLE' : 
    191     file_RUN_his = os.path.join (  dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 
     253    if file_RUN_his == None :  
     254        file_RUN_his = os.path.join (  dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 
     255        config['Files']['file_RUN_his'] = file_RUN_his 
    192256 
    193257d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 
    194258d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 
    195 if Routing == 'ORCHIDEE' : 
     259if Routing == 'SECHIBA' : 
    196260    d_RUN_his = d_SRF_his 
    197261if Routing == 'SIMPLE' :  
    … …  
    200264echo ( file_ATM_his ) 
    201265echo ( file_SRF_his ) 
    202 echo ( file_RUN_his ) 
    203  
    204  
    205 config['Files']['file_ATM_his'] = file_ATM_his 
    206 config['Files']['file_SRF_his'] = file_SRF_his 
    207 config['Files']['file_RUN_his'] = file_SRF_his 
     266if Routing == 'SIMPLE' :  
     267    echo ( file_RUN_his ) 
    208268 
    209269## Compute run length 
    … …  
    218278dtime_per_sec = xr.DataArray (dtime_per_sec, dims=["time_counter", ], coords=[d_ATM_his.time_counter,] ) 
    219279 
     280## Number of years 
     281NbYear = dtime_sec / YearLength 
     282 
    220283#-- Open restart files 
     284 
    221285YearRes = YearBegin - 1              # Year of the restart of beginning of simulation 
    222286YearPre = YearBegin - PackFrequency  # Year to find the tarfile of the restart of beginning of simulation 
    223287 
    224 echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ') 
    225  
    226 file_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{YearPre}0101_{YearRes}1231_restart.tar' ) 
    227 file_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{YearBegin}0101_{YearEnd}1231_restart.tar' ) 
    228  
     288if TarRestartPeriod_beg == None :  
     289    echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ') 
     290    TarRestartPeriod_beg = f'{YearPre}0101_{YearRes}1231' 
     291    config['Files']['TarRestartPeriod_beg'] = TarRestartPeriod_beg 
     292 
     293if TarRestartPeriod_end == None :  
     294    YearPre = YearBegin - PackFrequency  # Year to find the tarfile of the restart of beginning of simulation 
     295    echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ') 
     296    TarRestartPeriod_end = f'{YearBegin}0101_{YearEnd}1231' 
     297    config['Files']['TarRestartPeriod_end'] = TarRestartPeriod_end 
     298 
     299if file_restart_beg == None : 
     300    file_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_beg}_restart.tar' ) 
     301    config['Files']['file_restart_beg'] = file_restart_beg 
     302if file_restart_end == None : 
     303    file_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_end}_restart.tar' ) 
     304    config['Files']['file_restart_end'] = file_restart_end 
     305     
    229306echo ( f'{file_restart_beg}' ) 
    230307echo ( f'{file_restart_end}' ) 
    231308 
    232 file_ATM_beg = f'ATM_{JobName}_{YearRes}1231_restartphy.nc' 
    233 file_ATM_end = f'ATM_{JobName}_{YearEnd}1231_restartphy.nc' 
    234 if LMDZ : 
    235     file_DYN_beg = f'ATM_{JobName}_{YearRes}1231_restart.nc' 
    236     file_DYN_end = f'ATM_{JobName}_{YearEnd}1231_restart.nc' 
    237 if ICO : 
    238     file_DYN_beg = f'ICO_{JobName}_{YearRes}1231_restart.nc' 
    239     file_DYN_end = f'ICO_{JobName}_{YearEnd}1231_restart.nc' 
    240 file_SRF_beg = f'SRF_{JobName}_{YearRes}1231_sechiba_rest.nc' 
    241 file_SRF_end = f'SRF_{JobName}_{YearEnd}1231_sechiba_rest.nc' 
    242 liste_beg = [file_ATM_beg, file_DYN_beg, file_SRF_beg, ] 
    243 liste_end = [file_ATM_end, file_DYN_end, file_SRF_end, ] 
     309if file_ATM_beg == None :  
     310    file_ATM_beg = f'{TmpDir}/ATM_{JobName}_{YearRes}1231_restartphy.nc' 
     311    config['Files']['file_ATM_beg'] = file_ATM_beg 
     312if file_ATM_end == None : 
     313    file_ATM_end = f'{TmpDir}/ATM_{JobName}_{YearEnd}1231_restartphy.nc' 
     314    config['Files']['file_ATM_end'] = file_ATM_end 
     315 
     316liste_beg = [file_ATM_beg, ] 
     317liste_end = [file_ATM_end, ] 
     318     
     319 
     320if file_DYN_beg == None : 
     321    if LMDZ : 
     322        file_DYN_beg = f'{TmpDir}/ATM_{JobName}_{YearRes}1231_restart.nc' 
     323    if ICO  : 
     324        file_DYN_beg = f'{TmpDir}/ICO_{JobName}_{YearRes}1231_restart.nc' 
     325    liste_beg.append (file_DYN_beg) 
     326    config['Files']['file_DYN_beg'] = file_DYN_beg 
     327     
     328if file_DYN_end == None :  
     329    if LMDZ : 
     330        file_DYN_end = f'{TmpDir}/ATM_{JobName}_{YearEnd}1231_restart.nc' 
     331    if ICO  : 
     332        file_DYN_end = f'{TmpDir}/ICO_{JobName}_{YearEnd}1231_restart.nc' 
     333    liste_end.append ( file_DYN_end ) 
     334    config['Files']['file_DYN_end'] = file_DYN_end 
     335 
     336if file_SRF_beg == None :  
     337    file_SRF_beg = f'{TmpDir}/SRF_{JobName}_{YearRes}1231_sechiba_rest.nc' 
     338    config['Files']['file_SRF_beg'] = file_SRF_beg 
     339if file_SRF_end == None :  
     340    file_SRF_end = f'{TmpDir}/SRF_{JobName}_{YearEnd}1231_sechiba_rest.nc' 
     341    config['Files']['file_SRF_end'] = file_SRF_end 
     342     
     343liste_beg.append ( file_SRF_beg ) 
     344liste_end.append ( file_SRF_end ) 
     345 
    244346echo ( f'{file_ATM_beg}') 
    245347echo ( f'{file_ATM_end}') 
     348echo ( f'{file_DYN_beg}') 
     349echo ( f'{file_DYN_end}') 
    246350echo ( f'{file_SRF_beg}') 
    247351echo ( f'{file_SRF_end}') 
    248352 
    249 if Routing == 'SIMPLE' :  
    250     file_RUN_beg = f'SRF_{JobName}_{YearRes}1231_routing_restart.nc' 
    251     file_RUN_end = f'SRF_{JobName}_{YearEnd}1231_routing_restart.nc' 
     353if Routing == 'SIMPLE' : 
     354    if file_RUN_beg == None :  
     355        file_RUN_beg = f'{TmpDir}/SRF_{JobName}_{YearRes}1231_routing_restart.nc' 
     356        config['Files']['file_RUN_beg'] = file_RUN_beg 
     357    if file_RUN_end == None :  
     358        file_RUN_end = f'{TmpDir}/SRF_{JobName}_{YearEnd}1231_routing_restart.nc' 
     359        config['Files']['file_RUN_end'] = file_RUN_end 
     360         
    252361    liste_beg.append ( file_RUN_beg ) 
    253362    liste_end.append ( file_RUN_end ) 
    … …  
    257366echo ('\nExtract restart files from tar : ATM, ICO and SRF') 
    258367for resFile in liste_beg : 
    259     if not os.path.exists ( os.path.join (TmpDir, resFile) ) : 
    260         command =  f'cd {TmpDir} ; tar xf {file_restart_beg} {resFile}' 
     368    if os.path.exists ( os.path.join (TmpDir, resFile) ) : 
     369        echo ( f'file found : {resFile}' ) 
     370    else : 
     371        base_resFile = os.path.basename (resFile) 
     372        command =  f'cd {TmpDir} ; tar xf {file_restart_beg} {base_resFile}' 
    261373        echo ( command ) 
    262374        os.system ( command ) 
    263375         
    264376for resFile in liste_end : 
    265     if not os.path.exists ( os.path.join (TmpDir, resFile) ) : 
    266         command =  f'cd {TmpDir} ; tar xf {file_restart_end} {resFile}' 
     377    if os.path.exists ( os.path.join (TmpDir, resFile) ) : 
     378        echo ( f'file found : {resFile}' ) 
     379    else : 
     380        base_resFile = os.path.basename (resFile) 
     381        command =  f'cd {TmpDir} ; tar xf {file_restart_end} {base_resFile}' 
    267382        echo ( command ) 
    268383        os.system ( command ) 
    269  
    270 config['Files']['file_ATM_beg'] = file_ATM_beg 
    271 config['Files']['file_ATM_end'] = file_ATM_end 
    272 config['Files']['file_SRF_beg'] = file_SRF_beg 
    273 config['Files']['file_SRF_end'] = file_SRF_end 
    274 if Routing == 'SIMPLE' :  
    275     config['Files']['file_RUN_beg'] = file_RUN_beg 
    276     config['Files']['file_RUN_end'] = file_RUN_end 
    277 config['Files']['file_DYN_beg'] = file_DYN_beg 
    278 config['Files']['file_DYN_end'] = file_DYN_end 
    279          
     384        
    280385echo ('\nOpening ATM SRF and ICO restart files') 
    281386d_ATM_beg = xr.open_dataset ( os.path.join (TmpDir, file_ATM_beg), decode_times=False, decode_cf=True).squeeze() 
    … …  
    287392 
    288393for var in d_SRF_beg.variables : 
    289     d_SRF_beg[var] = d_SRF_beg[var].where (  d_SRF_beg[var]<1.e20, 0.) 
    290     d_SRF_end[var] = d_SRF_end[var].where (  d_SRF_end[var]<1.e20, 0.) 
    291  
    292 if ICO : 
     394    d_SRF_beg[var] = d_SRF_beg[var].where ( d_SRF_beg[var]<1.e20, 0.) 
     395    d_SRF_end[var] = d_SRF_end[var].where ( d_SRF_end[var]<1.e20, 0.) 
     396 
     397if Routing == 'SIMPLE' : 
    293398    d_RUN_beg = xr.open_dataset ( os.path.join (TmpDir, file_RUN_beg), decode_times=False, decode_cf=True).squeeze() 
    294399    d_RUN_end = xr.open_dataset ( os.path.join (TmpDir, file_RUN_end), decode_times=False, decode_cf=True).squeeze() 
    … …  
    304409    echo ( file_RUN_end ) 
    305410 
     411## 
     412config_out = open (FullIniFile, 'w') 
     413config.write (config_out ) 
     414config_out.close () 
     415 
     416def Ksum (tab) : 
     417    ''' 
     418    Kahan summation algorithm, also known as compensated summation 
     419    https://en.wikipedia.org/wiki/Kahan_summation_algorithm 
     420    ''' 
     421    Ksum = 0.0                   # Prepare the accumulator. 
     422    comp = 0.0                   # A running compensation for lost low-order bits. 
     423    
     424    for xx in np.where ( np.isnan(tab), 0., tab ) :  
     425        yy = xx - comp           # comp is zero the first time around. 
     426        tt = Ksum + yy           # Alas, sum is big, y small, so low-order digits of y are lost. 
     427        comp = (tt - Ksum) - yy  # (tt - Ksum) cancels the high-order part of y; subtracting y recovers negative (low part of yy) 
     428        Ksum = tt                # Algebraically, comp should always be zero. Beware overly-aggressive optimizing compilers! 
     429                                 # Next time around, the lost low part will be added to y in a fresh attempt. 
     430    return Ksum 
     431 
     432def Ssum (tab) : 
     433    ''' 
     434    Precision summation by sorting by absolute value 
     435    ''' 
     436    Ssum = np.sum ( tab[np.argsort(np.abs(tab))] ) 
     437    return Ssum 
     438 
     439def KSsum (tab) : 
     440    ''' 
     441    Precision summation by sorting by absolute value, and applying Kahan summation 
     442    ''' 
     443    KSsum = Ksum ( tab[np.argsort(np.abs(tab))] ) 
     444    return KSsum 
     445 
     446# Choosing algorithm 
     447Psum = KSsum 
     448 
    306449def ATM_stock_int (stock) : 
    307450    '''Integrate stock on atmosphere grid''' 
    308     ATM_stock_int  = np.sum (  np.sort ( (stock * DYN_aire).to_masked_array().ravel()) ) 
     451    ATM_stock_int  = Psum ( (stock * DYN_aire).to_masked_array().ravel() )  
    309452    return ATM_stock_int 
    310453 
    311454def ATM_flux_int (flux) : 
    312455    '''Integrate flux on atmosphere grid''' 
    313     ATM_stock_int  = np.sum (  np.sort ( (flux * dtime_per_sec * ATM_aire).to_masked_array().ravel()) ) 
     456    ATM_flux_int  = Psum ( (flux * dtime_per_sec * ATM_aire).to_masked_array().ravel() ) 
     457    return ATM_flux_int 
     458 
     459def SRF_stock_int (stock) : 
     460    '''Integrate stock on land grid''' 
     461    ATM_stock_int  = Ksum (  ( (stock * DYN_aire_fter).to_masked_array().ravel()) ) 
    314462    return ATM_stock_int 
     463 
     464def SRF_flux_int (flux) : 
     465    '''Integrate flux on land grid''' 
     466    SRF_flux_int  = Psum (  (flux * dtime_per_sec * SRF_aire).to_masked_array().ravel() ) 
     467    return SRF_flux_int 
    315468 
    316469def ONE_stock_int (stock) : 
    317470    '''Sum stock ''' 
    318     ONE_stock_int  = np.sum (  np.sort ( (stock).to_masked_array().ravel()) ) 
     471    ONE_stock_int  = Psum ( stock.to_masked_array().ravel() ) 
    319472    return ONE_stock_int 
    320473 
    321474def ONE_flux_int (flux) : 
    322475    '''Sum flux ''' 
    323     ONE_flux_int  = np.sum (  np.sort ( (flux * dtime_per_sec ).to_masked_array().ravel()) ) 
     476    ONE_flux_int  = Psum ( (flux * dtime_per_sec ).to_masked_array().ravel() ) 
    324477    return ONE_flux_int 
     478 
     479def kg2Sv  (val, rho=ATM_rho) : 
     480    '''From kg to Sverdrup''' 
     481    return val/dtime_sec*1.0e-6/rho 
     482 
     483def kg2myear (val, rho=ATM_rho) : 
     484    '''From kg to m/year''' 
     485    return val/ATM_aire_sea_tot/rho/NbYear 
    325486     
    326487# ATM grid with cell surfaces 
    … …  
    329490    file_ATM_aire = os.path.join ( R_IN, 'ATM', 'GRID', f'aire_{ATM}_to_{jpia}x{jpja}.nc' ) 
    330491    config['Files']['file_ATM_aire'] = file_ATM_aire 
    331     echo ('Aire sur grille reguliere :', file_ATM_aire) 
     492    echo ( f'Aire sur grille reguliere : {file_ATM_aire}' ) 
    332493    d_ATM_aire = xr.open_dataset ( file_ATM_aire, decode_times=False ).squeeze() 
    333494    ATM_aire = lmdz.geo2point ( d_ATM_aire ['aire'].squeeze(), cumulPoles=True ) 
    334     ATM_fsea = lmdz.geo2point ( d_ATM_his ['fract_oce'][0] + d_ATM_his ['fract_sic'][0] ) 
    335     ATM_flnd = lmdz.geo2point ( d_ATM_his ['fract_ter'][0] + d_ATM_his ['fract_lic'][0] ) 
     495    ATM_fter = lmdz.geo2point ( d_ATM_his ['fract_ter'][0] ) 
     496    ATM_foce = lmdz.geo2point ( d_ATM_his ['fract_oce'][0] ) 
     497    ATM_fsic = lmdz.geo2point ( d_ATM_his ['fract_sic'][0] ) 
     498    ATM_flic = lmdz.geo2point ( d_ATM_his ['fract_lic'][0] ) 
     499    ATM_fsea = ATM_foce + ATM_fsic 
     500    ATM_flnd = ATM_fter + ATM_flic 
    336501    SRF_aire = lmdz.geo2point ( d_SRF_his ['Areas'] * d_SRF_his ['Contfrac'] ) 
     502    #SRF_aire = ATM_aire * lmdz.geo2point (d_SRF_his ['Contfrac'] ) 
     503    #SRF_aire = ATM_aire * ATM_fter 
    337504     
    338505if LMDZ : 
    … …  
    340507    ATM_fsea = lmdz.geo2point ( d_ATM_his ['fract_oce'][0] + d_ATM_his ['fract_sic'][0] ) 
    341508    ATM_flnd = lmdz.geo2point ( d_ATM_his ['fract_ter'][0] ) 
    342     SRF_aire = lmdz.geo2point ( d_SRF_his['Areas'] * d_SRF_his['Contfrac'] ) 
     509    #SRF_aire = lmdz.geo2point ( d_SRF_his['Areas'] * d_SRF_his['Contfrac'] ) 
     510    SRF_aire = ATM_aire * lmdz.geo2point ( d_SRF_his['Contfrac'] ) 
    343511 
    344512SRF_aire = SRF_aire.where ( SRF_aire < 1E15, 0.) 
    … …  
    353521    DYN_fsea = d_DYN_aire ['fract_oce'] + d_DYN_aire ['fract_sic'] 
    354522    DYN_flnd = 1.0 - DYN_fsea 
     523    DYN_fter = d_ATM_beg['FTER'].rename({'points_physiques':'cell_mesh'}) 
     524    DYN_flic = d_ATM_beg['FTER'].rename({'points_physiques':'cell_mesh'}) 
    355525     
    356526if LMDZ : 
    … …  
    358528    DYN_fsea = ATM_fsea 
    359529    DYN_flnd = ATM_flnd 
     530    DYN_fter = d_ATM_beg['FTER'] 
     531    DYN_flic = d_ATM_beg['FTER'] 
     532 
     533     
     534DYN_aire_fter = DYN_aire * DYN_fter 
    360535 
    361536#if LMDZ : 
    … …  
    374549    raise Exception ('Erreur surface interpolee sur grille reguliere') 
    375550 
    376 echo ( '\n------------------------------------------------------------------------------------' ) 
    377 echo ( '-- ATM changes in stores ' ) 
     551echo ( '\n====================================================================================' ) 
     552echo ( f'-- ATM changes in stores  -- {Title} ' ) 
    378553 
    379554#-- Change in precipitable water from the atmosphere daily and monthly files 
    … …  
    386561# Surface pressure 
    387562if ICO :  
    388     DYN_ps_beg = d_DYN_beg['ps'] 
    389     DYN_ps_end = d_DYN_end['ps'] 
    390      
     563    DYN_psol_beg = d_DYN_beg['ps'] 
     564    DYN_psol_end = d_DYN_end['ps'] 
    391565if LMDZ :  
    392     DYN_ps_beg = lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)) ) 
    393     DYN_ps_end = lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)) ) 
    394      
    395 # 3D Pressure 
    396 DYN_p_beg = ATM_Ahyb + ATM_Bhyb * DYN_ps_beg 
    397 DYN_p_end = ATM_Ahyb + ATM_Bhyb * DYN_ps_end 
    398  
    399 if ICO  : klevp1, cell_mesh        = DYN_p_beg.shape 
    400 if LMDZ : klevp1, points_physiques = DYN_p_beg.shape 
     566    DYN_psol_beg = lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)) ) 
     567    DYN_psol_end = lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)) ) 
     568     
     569# 3D Pressure at the interface layers (not scalar points) 
     570DYN_pres_beg = ATM_Ahyb + ATM_Bhyb * DYN_psol_beg 
     571DYN_pres_end = ATM_Ahyb + ATM_Bhyb * DYN_psol_end 
     572 
     573klevp1 = ATM_Bhyb.shape[-1] 
     574if ICO  : cell_mesh        = DYN_psol_beg.shape[-1] 
     575if LMDZ : points_physiques = DYN_psol_beg.shape[-1] 
    401576klev = klevp1 - 1 
    402577 
    403 # Layer thickness 
     578# Layer thickness (pressure) 
    404579if ICO :  
    405     DYN_dsigma_beg = xr.DataArray ( np.empty( (klev, cell_mesh       )), dims=('sigs', 'cell_mesh'       ), coords=(np.arange(klev), np.arange(cell_mesh) ) ) 
    406     DYN_dsigma_end = xr.DataArray ( np.empty( (klev, cell_mesh       )), dims=('sigs', 'cell_mesh'       ), coords=(np.arange(klev), np.arange(cell_mesh) ) ) 
     580    DYN_dpres_beg = xr.DataArray ( np.empty( (klev, cell_mesh       )), dims=('sigs', 'cell_mesh'       ), coords=(np.arange(klev), np.arange(cell_mesh)        ) ) 
     581    DYN_dpres_end = xr.DataArray ( np.empty( (klev, cell_mesh       )), dims=('sigs', 'cell_mesh'       ), coords=(np.arange(klev), np.arange(cell_mesh)        ) ) 
    407582if LMDZ :  
    408     DYN_dsigma_beg = xr.DataArray ( np.empty( (klev, points_physiques)), dims=('sigs', 'points_physiques'), coords=(np.arange(klev), np.arange(points_physiques) ) ) 
    409     DYN_dsigma_end = xr.DataArray ( np.empty( (klev, points_physiques)), dims=('sigs', 'points_physiques'), coords=(np.arange(klev), np.arange(points_physiques) ) ) 
     583    DYN_dpres_beg = xr.DataArray ( np.empty( (klev, points_physiques)), dims=('sigs', 'points_physiques'), coords=(np.arange(klev), np.arange(points_physiques) ) ) 
     584    DYN_dpres_end = xr.DataArray ( np.empty( (klev, points_physiques)), dims=('sigs', 'points_physiques'), coords=(np.arange(klev), np.arange(points_physiques) ) ) 
    410585 
    411586for k in np.arange (klevp1-1) : 
    412     DYN_dsigma_beg[k,:] = DYN_p_beg[k,:] - DYN_p_beg[k+1,:] 
    413     DYN_dsigma_end[k,:] = DYN_p_end[k,:] - DYN_p_end[k+1,:] 
     587    DYN_dpres_beg[k,:] = DYN_pres_beg[k,:] - DYN_pres_beg[k+1,:] 
     588    DYN_dpres_end[k,:] = DYN_pres_end[k,:] - DYN_pres_end[k+1,:] 
    414589 
    415590##-- Vertical and horizontal integral, and sum of liquid, solid and vapor water phases 
    … …  
    429604        DYN_wat_end = (d_DYN_end['q'].isel(nq=0) + d_DYN_end['q'].isel(nq=1) + d_DYN_end['q'].isel(nq=2) ).rename ( {'lev':'sigs'} ) 
    430605     
    431 # Integral 
    432 DYN_mas_wat_beg = ATM_stock_int (DYN_dsigma_beg * DYN_wat_beg) / Grav 
    433 DYN_mas_wat_end = ATM_stock_int (DYN_dsigma_end * DYN_wat_end) / Grav 
    434  
     606# Compute water content : vertical and horizontal integral 
     607DYN_mas_wat_beg = ATM_stock_int (DYN_dpres_beg * DYN_wat_beg) / Grav 
     608DYN_mas_wat_end = ATM_stock_int (DYN_dpres_end * DYN_wat_end) / Grav 
     609 
     610# Variation of water content 
    435611dDYN_mas_wat = DYN_mas_wat_end - DYN_mas_wat_beg 
    436612 
    437 echo ( '\nVariation du contenu en eau atmosphere (dynamique) ' ) 
     613# \([a-z,A-Z,_]*\)/dtime_sec\*1e-9  kg2Sv(\1) 
     614# \([a-z,A-Z,_]*\)\/ATM_aire_sea_tot\/ATM_rho\/NbYear  kg2myear(\1) 
     615 
     616def var2prt (var, small=False) : 
     617    if small :  return var , kg2Sv(var)*1000., kg2myear(var)*1000. 
     618    else     :  return var , kg2Sv(var)      , kg2myear(var) 
     619 
     620def prtFlux (Desc, var, Form='F', small=False) : 
     621    if small : 
     622        if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} mSv | {:12.4f} mm/year " 
     623        if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} mSv | {:12.4e} mm/year " 
     624    else : 
     625        if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} Sv  | {:12.4f} m/year  " 
     626        if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} Sv  | {:12.4e} m/year  " 
     627    echo ( (' {:>15} = ' +ff).format (Desc, *var2prt(var, small) ) ) 
     628    return None 
     629 
     630echo ( f'\nVariation du contenu en eau atmosphere (dynamique)  -- {Title} ' ) 
     631echo ( '------------------------------------------------------------------------------------' ) 
    438632echo ( 'DYN_mas_beg = {:12.6e} kg | DYN_mas_end = {:12.6e} kg'.format (DYN_mas_wat_beg, DYN_mas_wat_end) ) 
    439 echo ( 'dMass (atm)   = {:12.3e} kg '.format (dDYN_mas_wat) ) 
    440 echo ( 'dMass (atm)   = {:12.3e} Sv '.format (dDYN_mas_wat/dtime_sec*1.e-6/ATM_rho) ) 
    441 echo ( 'dMass (atm)   = {:12.3e} m  '.format (dDYN_mas_wat/ATM_aire_sea_tot/ATM_rho) ) 
     633prtFlux ( 'dMass (atm)  ', dDYN_mas_wat, 'e', True ) 
    442634 
    443635ATM_sno_beg = d_ATM_beg['SNOW01']*d_ATM_beg['FTER']+d_ATM_beg['SNOW02']*d_ATM_beg['FLIC']+d_ATM_beg['SNOW03']*d_ATM_beg['FOCE']+d_ATM_beg['SNOW04']*d_ATM_beg['FSIC'] 
    … …  
    499691dATM_mas_qs04 = ATM_mas_qs04_end - ATM_mas_qs04_beg 
    500692 
    501 echo ( '\nVariation du contenu en neige atmosphere (calottes)' ) 
     693dATM_mas_sno_2 = ATM_stock_int ( ATM_sno_end - ATM_sno_beg ) 
     694 
     695echo ( f'\nVariation du contenu en neige atmosphere (calottes)  -- {Title} ' ) 
     696echo ( '------------------------------------------------------------------------------------' ) 
    502697echo ( 'ATM_mas_sno_beg  = {:12.6e} kg | ATM_mas_sno_end  = {:12.6e} kg'.format (ATM_mas_sno_beg, ATM_mas_sno_end) ) 
    503 echo ( 'dMass (neige atm) = {:12.3e} kg '.format (dATM_mas_sno ) ) 
    504 echo ( 'dMass (neige atm) = {:12.3e} Sv '.format (dATM_mas_sno/dtime_sec*1e-6/ICE_rho_ice) ) 
    505 echo ( 'dMass (neige atm) = {:12.3e} m  '.format (dATM_mas_sno/ATM_aire_sea_tot/ATM_rho) ) 
    506  
    507 echo ( '\nVariation du contenu humidite du sol' ) 
     698prtFlux ( 'dMass (neige atm) ', dATM_mas_sno  , 'e', True ) 
     699prtFlux ( 'dMass (neige atm) ', dATM_mas_sno_2, 'e', True ) 
     700 
     701echo ( f'\nVariation du contenu humidite du sol  -- {Title} ' ) 
     702echo ( '------------------------------------------------------------------------------------' ) 
    508703echo ( 'ATM_mas_qs_beg  = {:12.6e} kg | ATM_mas_qs_end  = {:12.6e} kg'.format (ATM_mas_qs_beg, ATM_mas_qs_end) ) 
    509 echo ( 'dMass (neige atm) = {:12.3e} kg '.format (dATM_mas_qs ) ) 
    510 echo ( 'dMass (neige atm) = {:12.3e} Sv '.format (dATM_mas_qs/dtime_sec*1e-6/ATM_rho) ) 
    511 echo ( 'dMass (neige atm) = {:12.3e} m  '.format (dATM_mas_qs/ATM_aire_sea_tot/ATM_rho) ) 
    512  
    513 echo ( '\nVariation du contenu en eau+neige atmosphere ' ) 
    514 echo ( 'dMass (eau + neige atm) = {:12.3e} kg '.format (  dDYN_mas_wat + dATM_mas_sno) ) 
    515 echo ( 'dMass (eau + neige atm) = {:12.3e} Sv '.format ( (dDYN_mas_wat + dATM_mas_sno)/dtime_sec*1E-6/ATM_rho) ) 
    516 echo ( 'dMass (eau + neige atm) = {:12.3e} m  '.format ( (dDYN_mas_wat + dATM_mas_sno)/ATM_aire_sea_tot/ATM_rho) ) 
    517  
    518 echo ( '\n------------------------------------------------------------------------------------' ) 
    519 echo ( '-- SRF changes ' ) 
    520  
    521 # dSoilHum_in_Sv=`cdo outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=(maxvegetfrac*DelSoilMoist_daily)*Contfrac' ${file} -gridarea ${file}` 
    522 # dInterce_in_Sv=`cdo outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=DelIntercept_daily*Contfrac' ${file} -gridarea ${file}` 
    523 # dSWE_in_Sv=`cdo outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=DelSWE_daily*Contfrac' ${file} -gridarea ${file}` 
    524 # dStream_in_Sv=`cdo outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=delstreamr_daily*Contfrac' ${file} -gridarea ${file}` 
    525 # dFastR_in_Sv=`cdo outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=delfastr_daily*Contfrac' ${file} -gridarea ${file}` 
    526 # dSlowR_in_Sv=`cdo outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=delslowr_daily*Contfrac' ${file} -gridarea ${file}` 
    527 # dFlood_in_Sv=`cdo outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=delfloodr_daily*Contfrac' ${file} -gridarea ${file}` 
    528 # dPond_in_Sv=`cdo outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=delpondr_daily*Contfrac' ${file} -gridarea ${file}` 
    529 # dLake_in_Sv=`cdo outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=dellakevol_daily*Contfrac' ${file} -gridarea ${file}` 
    530 #echo 'dSTOCK (Sv) Soil Intercept SWE Stream FastR SlowR Lake Pond Flood=' $dSoilHum_in_Sv, $dInterce_in_Sv, $dSWE_in_Sv, $dStream_in_Sv, $dFastR_in_Sv, $dSlowR_in_Sv, $dLake_in_Sv, $dPond_in_Sv, $dFlood_in_Sv >> ${fileout} 
    531 #dSRF_tot=`python -c "print $dSoilHum_in_Sv+$dInterce_in_Sv+$dSWE_in_Sv+$dStream_in_Sv+$dFastR_in_Sv+$dSlowR_in_Sv+$dLake_in_Sv+$dPond_in_Sv+$dFlood_in_Sv"` 
    532 #echo 'dSTOCK (Sv) total='${dSRF_tot} >> ${fileout} 
    533  
     704prtFlux ( 'dMass (neige atm) ', dATM_mas_qs, 'e', True ) 
     705 
     706echo ( f'\nVariation du contenu en eau+neige atmosphere  -- {Title}  ' ) 
     707echo ( '------------------------------------------------------------------------------------' ) 
     708prtFlux ( 'dMass (eau + neige atm) ', dDYN_mas_wat + dATM_mas_sno , 'e', True) 
     709 
     710echo ( '\n====================================================================================' ) 
     711echo ( f'-- SRF changes  -- {Title} ' ) 
    534712 
    535713if Routing == 'SIMPLE' : 
    536     RUN_mas_wat_beg = ONE_stock_int ( d_RUN_beg ['fast_reservoir'] +  d_RUN_beg ['slow_reservoir'] +  d_RUN_beg ['stream_reservoir'] ) 
    537     RUN_mas_wat_end = ONE_stock_int ( d_RUN_end ['fast_reservoir'] +  d_RUN_end ['slow_reservoir'] +  d_RUN_end ['stream_reservoir'] ) 
    538  
    539 if Routing == 'ORCHIDEE' :  
    540     RUN_mas_wat_beg = ONE_stock_int ( d_SRF_beg['fastres']  + d_SRF_beg['slowres'] + d_SRF_beg['streamres'] \ 
    541                                     + d_SRF_beg['floodres'] + d_SRF_beg['lakeres'] + d_SRF_beg['pondres']  ) 
    542     RUN_mas_wat_end = ONE_stock_int ( d_SRF_end['fastres']  + d_SRF_end['slowres'] + d_SRF_end['streamres'] \ 
    543                                     + d_SRF_end['floodres'] + d_SRF_end['lakeres'] + d_SRF_end['pondres']  ) 
    544  
    545 dRUN_mas_wat = RUN_mas_wat_end - RUN_mas_wat_beg 
    546      
    547 echo ( '\nWater content in routing ' ) 
     714    RUN_mas_wat_fast_beg   = ONE_stock_int ( d_RUN_beg ['fast_reservoir']   ) 
     715    RUN_mas_wat_slow_beg   = ONE_stock_int ( d_RUN_beg ['slow_reservoir']   ) 
     716    RUN_mas_wat_stream_beg = ONE_stock_int ( d_RUN_beg ['stream_reservoir'] ) 
     717    RUN_mas_wat_flood_beg  = 0.0 
     718    RUN_mas_wat_lake_beg   = 0.0 
     719    RUN_mas_wat_pond_beg   = 0.0 
     720     
     721    RUN_mas_wat_fast_end   = ONE_stock_int ( d_RUN_end ['fast_reservoir']   ) 
     722    RUN_mas_wat_slow_end   = ONE_stock_int ( d_RUN_end ['slow_reservoir']   ) 
     723    RUN_mas_wat_stream_end = ONE_stock_int ( d_RUN_end ['stream_reservoir'] ) 
     724    RUN_mas_wat_flood_end  = 0.0 
     725    RUN_mas_wat_lake_end   = 0.0 
     726    RUN_mas_wat_pond_end   = 0.0 
     727 
     728if Routing == 'SECHIBA' : 
     729    RUN_mas_wat_fast_beg   = ONE_stock_int ( d_SRF_beg ['fastres']   ) 
     730    RUN_mas_wat_slow_beg   = ONE_stock_int ( d_SRF_beg ['slowres']   ) 
     731    RUN_mas_wat_stream_beg = ONE_stock_int ( d_SRF_beg ['streamres'] ) 
     732    RUN_mas_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  ) 
     733    RUN_mas_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   ) 
     734    RUN_mas_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   ) 
     735 
     736     
     737    RUN_mas_wat_fast_end   = ONE_stock_int ( d_SRF_end ['fastres']   ) 
     738    RUN_mas_wat_slow_end   = ONE_stock_int ( d_SRF_end ['slowres']   ) 
     739    RUN_mas_wat_stream_end = ONE_stock_int ( d_SRF_end ['streamres'] ) 
     740    RUN_mas_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  ) 
     741    RUN_mas_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   ) 
     742    RUN_mas_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   ) 
     743 
     744RUN_mas_wat_beg = RUN_mas_wat_fast_beg + RUN_mas_wat_slow_beg + RUN_mas_wat_stream_beg + RUN_mas_wat_flood_beg + RUN_mas_wat_lake_beg + RUN_mas_wat_pond_beg 
     745RUN_mas_wat_end = RUN_mas_wat_fast_end + RUN_mas_wat_slow_end + RUN_mas_wat_stream_end + RUN_mas_wat_flood_end + RUN_mas_wat_lake_end + RUN_mas_wat_pond_end 
     746 
     747dRUN_mas_wat_fast   = RUN_mas_wat_fast_end   - RUN_mas_wat_fast_beg 
     748dRUN_mas_wat_slow   = RUN_mas_wat_slow_end   - RUN_mas_wat_slow_beg 
     749dRUN_mas_wat_stream = RUN_mas_wat_stream_end - RUN_mas_wat_stream_beg 
     750dRUN_mas_wat_flood  = RUN_mas_wat_flood_end  - RUN_mas_wat_flood_beg 
     751dRUN_mas_wat_lake   = RUN_mas_wat_lake_end   - RUN_mas_wat_lake_beg 
     752dRUN_mas_wat_pond   = RUN_mas_wat_pond_end   - RUN_mas_wat_pond_beg 
     753 
     754dRUN_mas_wat        = RUN_mas_wat_end        - RUN_mas_wat_beg 
     755 
     756echo ( f'\nLes differents reservoirs  -- {Title} ') 
     757echo ( '------------------------------------------------------------------------------------' ) 
     758echo ( 'RUN_mas_wat_fast_beg   = {:12.6e} kg | RUN_mas_wat_fast_end   = {:12.6e} kg '.format (RUN_mas_wat_fast_beg  , RUN_mas_wat_fast_end   ) ) 
     759echo ( 'RUN_mas_wat_slow_beg   = {:12.6e} kg | RUN_mas_wat_slow_end   = {:12.6e} kg '.format (RUN_mas_wat_slow_beg  , RUN_mas_wat_slow_end   ) ) 
     760echo ( 'RUN_mas_wat_stream_beg = {:12.6e} kg | RUN_mas_wat_stream_end = {:12.6e} kg '.format (RUN_mas_wat_stream_beg, RUN_mas_wat_stream_end ) ) 
     761echo ( 'RUN_mas_wat_flood_beg  = {:12.6e} kg | RUN_mas_wat_flood_end  = {:12.6e} kg '.format (RUN_mas_wat_flood_beg , RUN_mas_wat_flood_end  ) ) 
     762echo ( 'RUN_mas_wat_lake_beg   = {:12.6e} kg | RUN_mas_wat_lake_end   = {:12.6e} kg '.format (RUN_mas_wat_lake_beg  , RUN_mas_wat_lake_end   ) ) 
     763echo ( 'RUN_mas_wat_pond_beg   = {:12.6e} kg | RUN_mas_wat_pond_end   = {:12.6e} kg '.format (RUN_mas_wat_pond_beg  , RUN_mas_wat_pond_end   ) ) 
     764 
     765echo ( '------------------------------------------------------------------------------------' ) 
     766 
     767prtFlux ( 'dMass (fast)   ', dRUN_mas_wat_fast  , 'e', True ) 
     768prtFlux ( 'dMass (slow)   ', dRUN_mas_wat_slow  , 'e', True ) 
     769prtFlux ( 'dMass (stream) ', dRUN_mas_wat_stream, 'e', True ) 
     770prtFlux ( 'dMass (flood)  ', dRUN_mas_wat_flood , 'e', True ) 
     771prtFlux ( 'dMass (lake)   ', dRUN_mas_wat_lake  , 'e', True ) 
     772prtFlux ( 'dMass (pond)   ', dRUN_mas_wat_pond  , 'e', True ) 
     773prtFlux ( 'dMass (all)    ', dRUN_mas_wat       , 'e', True ) 
     774 
     775echo ( f'\nWater content in routing  -- {Title} ' ) 
     776echo ( '------------------------------------------------------------------------------------' ) 
    548777echo ( 'RUN_mas_wat_beg = {:12.6e} kg | RUN_mas_wat_end = {:12.6e} kg '.format (RUN_mas_wat_end, RUN_mas_wat_end) ) 
    549 echo ( 'dMass (routing)  = {:12.3e} kg '.format(dRUN_mas_wat) ) 
    550 echo ( 'dMass (routing)  = {:12.3e} Sv '.format(dRUN_mas_wat/dtime_sec*1E-9) ) 
    551 echo ( 'dMass (routing)  = {:12.3e} m  '.format(dRUN_mas_wat/ATM_aire_sea_tot*1E-3) ) 
    552  
    553 print ('Reading SRF restart') 
    554 tot_watveg_beg  = d_SRF_beg['tot_watveg_beg']  ; tot_watveg_beg  = tot_watveg_beg .where (tot_watveg_beg < 1E10, 0.) 
    555 tot_watsoil_beg = d_SRF_beg['tot_watsoil_beg'] ; tot_watsoil_beg = tot_watsoil_beg.where (tot_watsoil_beg< 1E10, 0.) 
    556 snow_beg        = d_SRF_beg['snow_beg']        ; snow_beg        = snow_beg       .where (snow_beg       < 1E10, 0.) 
    557  
    558 tot_watveg_end  = d_SRF_end['tot_watveg_beg']  ; tot_watveg_end  = tot_watveg_end .where (tot_watveg_end < 1E10, 0.) 
    559 tot_watsoil_end = d_SRF_end['tot_watsoil_beg'] ; tot_watsoil_end = tot_watsoil_end.where (tot_watsoil_end< 1E10, 0.) 
    560 snow_end        = d_SRF_end['snow_beg']        ; snow_end        = snow_end       .where (snow_end       < 1E10, 0.) 
     778prtFlux ( 'dMass (routing) ', dRUN_mas_wat       ) 
     779 
     780echo ( '\n====================================================================================' ) 
     781print (f'Reading SRF restart') 
     782SRF_tot_watveg_beg  = d_SRF_beg['tot_watveg_beg']  ; SRF_tot_watveg_beg  = SRF_tot_watveg_beg .where (SRF_tot_watveg_beg < 1E15, 0.) 
     783SRF_tot_watsoil_beg = d_SRF_beg['tot_watsoil_beg'] ; SRF_tot_watsoil_beg = SRF_tot_watsoil_beg.where (SRF_tot_watsoil_beg< 1E15, 0.) 
     784SRF_snow_beg        = d_SRF_beg['snow_beg']        ; SRF_snow_beg        = SRF_snow_beg       .where (SRF_snow_beg       < 1E15, 0.) 
     785SRF_lakeres_beg     = d_SRF_beg['lakeres']         ; SRF_lakeres_beg     = SRF_lakeres_beg    .where (SRF_lakeres_beg    < 1E15, 0.) 
     786 
     787SRF_tot_watveg_end  = d_SRF_end['tot_watveg_beg']  ; SRF_tot_watveg_end  = SRF_tot_watveg_end .where (SRF_tot_watveg_end < 1E15, 0.) 
     788SRF_tot_watsoil_end = d_SRF_end['tot_watsoil_beg'] ; SRF_tot_watsoil_end = SRF_tot_watsoil_end.where (SRF_tot_watsoil_end< 1E15, 0.) 
     789SRF_snow_end        = d_SRF_end['snow_beg']        ; SRF_snow_end        = SRF_snow_end       .where (SRF_snow_end       < 1E15, 0.) 
     790SRF_lakeres_end     = d_SRF_end['lakeres']         ; SRF_lakeres_end     = SRF_lakeres_end    .where (SRF_lakeres_end    < 1E15, 0.) 
    561791 
    562792if LMDZ : 
    563     tot_watveg_beg  = lmdz.geo2point (tot_watveg_beg) 
    564     tot_watsoil_beg = lmdz.geo2point (tot_watsoil_beg) 
    565     snow_beg        = lmdz.geo2point (snow_beg) 
    566     tot_watveg_end  = lmdz.geo2point (tot_watveg_end) 
    567     tot_watsoil_end = lmdz.geo2point (tot_watsoil_end) 
    568     snow_end        = lmdz.geo2point (snow_end) 
     793    SRF_tot_watveg_beg  = lmdz.geo2point (SRF_tot_watveg_beg ) 
     794    SRF_tot_watsoil_beg = lmdz.geo2point (SRF_tot_watsoil_beg) 
     795    SRF_snow_beg        = lmdz.geo2point (SRF_snow_beg       ) 
     796    SRF_lakeres_beg     = lmdz.geo2point (SRF_lakeres_beg    ) 
     797    SRF_tot_watveg_end  = lmdz.geo2point (SRF_tot_watveg_end ) 
     798    SRF_tot_watsoil_end = lmdz.geo2point (SRF_tot_watsoil_end) 
     799    SRF_snow_end        = lmdz.geo2point (SRF_snow_end       ) 
     800    SRF_lakeres_end     = lmdz.geo2point (SRF_lakeres_end    ) 
     801   
     802if ICO : 
     803    SRF_tot_watveg_beg  = SRF_tot_watveg_beg .rename ( {'y':'cell_mesh'} ) 
     804    SRF_tot_watsoil_beg = SRF_tot_watsoil_beg.rename ( {'y':'cell_mesh'} ) 
     805    SRF_snow_beg        = SRF_snow_beg       .rename ( {'y':'cell_mesh'} ) 
     806    SRF_lakeres_beg     = SRF_lakeres_beg    .rename ( {'y':'cell_mesh'} ) 
     807    SRF_tot_watveg_end  = SRF_tot_watveg_end .rename ( {'y':'cell_mesh'} ) 
     808    SRF_tot_watsoil_end = SRF_tot_watsoil_end.rename ( {'y':'cell_mesh'} ) 
     809    SRF_snow_end        = SRF_snow_end       .rename ( {'y':'cell_mesh'} ) 
     810    SRF_lakeres_end     = SRF_lakeres_end    .rename ( {'y':'cell_mesh'} ) 
    569811 
    570812# Stock dSoilHum dInterce dSWE dStream dFastR dSlowR dLake dPond dFlood 
    571813 
    572 SRF_wat_beg = tot_watveg_beg + tot_watsoil_beg + snow_beg 
    573 SRF_wat_end = tot_watveg_end + tot_watsoil_end + snow_end 
    574      
    575 if ICO : 
    576     tot_watveg_beg  = tot_watveg_beg .rename ( {'y':'cell_mesh'} ) 
    577     tot_watsoil_beg = tot_watsoil_beg.rename ( {'y':'cell_mesh'} ) 
    578     snow_beg        = snow_beg       .rename ( {'y':'cell_mesh'} ) 
    579     tot_watveg_end  = tot_watveg_end .rename ( {'y':'cell_mesh'} ) 
    580     tot_watsoil_end = tot_watsoil_end.rename ( {'y':'cell_mesh'} ) 
    581     snow_end        = snow_end       .rename ( {'y':'cell_mesh'} ) 
    582     SRF_wat_beg     = SRF_wat_beg    .rename ( {'y':'cell_mesh'} ) 
    583     SRF_wat_end     = SRF_wat_end    .rename ( {'y':'cell_mesh'} ) 
    584  
     814SRF_wat_beg = SRF_tot_watveg_beg + SRF_tot_watsoil_beg + SRF_snow_beg  
     815SRF_wat_end = SRF_tot_watveg_end + SRF_tot_watsoil_end + SRF_snow_end 
     816 
     817echo ( '\n====================================================================================' ) 
    585818print ('Computing integrals') 
    586819 
    587 print ( ' 1/6', end='' ) ; sys.stdout.flush () 
    588 SRF_mas_watveg_beg   = ATM_stock_int ( tot_watveg_beg  ) 
    589 print ( ' 2/6', end='' ) ; sys.stdout.flush () 
    590 SRF_mas_watsoil_beg  = ATM_stock_int ( tot_watsoil_beg ) 
    591 print ( ' 3/6', end='' ) ; sys.stdout.flush () 
    592 SRF_mas_snow_beg     = ATM_stock_int ( snow_beg        ) 
    593 print ( ' 4/6', end='' ) ; sys.stdout.flush () 
    594 SRF_mas_watveg_end   = ATM_stock_int ( tot_watveg_end  ) 
    595 print ( ' 5/6', end='' ) ; sys.stdout.flush () 
    596 SRF_mas_watsoil_end  = ATM_stock_int ( tot_watsoil_end ) 
    597 print ( ' 6/6', end='' ) ; sys.stdout.flush () 
    598 SRF_mas_snow_end     = ATM_stock_int ( snow_end        ) 
     820print ( ' 1/8', end='' ) ; sys.stdout.flush () 
     821SRF_mas_watveg_beg   = SRF_stock_int ( SRF_tot_watveg_beg  ) 
     822print ( ' 2/8', end='' ) ; sys.stdout.flush () 
     823SRF_mas_watsoil_beg  = SRF_stock_int ( SRF_tot_watsoil_beg ) 
     824print ( ' 3/8', end='' ) ; sys.stdout.flush () 
     825SRF_mas_snow_beg     = SRF_stock_int ( SRF_snow_beg        ) 
     826print ( ' 4/8', end='' ) ; sys.stdout.flush () 
     827SRF_mas_lake_beg     = ONE_stock_int ( SRF_lakeres_beg     ) 
     828print ( ' 5/8', end='' ) ; sys.stdout.flush () 
     829 
     830SRF_mas_watveg_end   = SRF_stock_int ( SRF_tot_watveg_end  ) 
     831print ( ' 6/8', end='' ) ; sys.stdout.flush () 
     832SRF_mas_watsoil_end  = SRF_stock_int ( SRF_tot_watsoil_end ) 
     833print ( ' 7/8', end='' ) ; sys.stdout.flush () 
     834SRF_mas_snow_end     = SRF_stock_int ( SRF_snow_end        ) 
     835print ( ' 8/8', end='' ) ; sys.stdout.flush () 
     836SRF_mas_lake_end     = ONE_stock_int ( SRF_lakeres_end     ) 
     837 
    599838print (' -- ') ; sys.stdout.flush () 
    600839 
    601 dSRF_mas_watveg   = SRF_mas_watveg_end   - SRF_mas_watveg_beg 
    602 dSRF_mas_watsoil  = SRF_mas_watsoil_end  - SRF_mas_watsoil_beg 
    603 dSRF_mas_snow     = SRF_mas_snow_end     - SRF_mas_snow_beg 
    604  
    605 echo ('\nLes differents reservoirs') 
     840dSRF_mas_watveg   = SRF_mas_watveg_end  - SRF_mas_watveg_beg 
     841dSRF_mas_watsoil  = SRF_mas_watsoil_end - SRF_mas_watsoil_beg 
     842dSRF_mas_snow     = SRF_mas_snow_end    - SRF_mas_snow_beg 
     843dSRF_mas_lake     = SRF_mas_lake_end    - SRF_mas_lake_beg 
     844 
     845 
     846 
     847echo ( '------------------------------------------------------------------------------------' ) 
     848echo ( f'\nLes differents reservoirs  -- {Title} ') 
    606849echo ( 'SRF_mas_watveg_beg   = {:12.6e} kg | SRF_mas_watveg_end   = {:12.6e} kg '.format (SRF_mas_watveg_beg  , SRF_mas_watveg_end   ) ) 
    607850echo ( 'SRF_mas_watsoil_beg  = {:12.6e} kg | SRF_mas_watsoil_end  = {:12.6e} kg '.format (SRF_mas_watsoil_beg , SRF_mas_watsoil_end  ) ) 
    608851echo ( 'SRF_mas_snow_beg     = {:12.6e} kg | SRF_mas_snow_end     = {:12.6e} kg '.format (SRF_mas_snow_beg    , SRF_mas_snow_end     ) ) 
    609  
    610 echo ( 'dMass (watveg)    = {:12.3e} kg | {:12.2e} Sv | {:12.2e} m '.format (dSRF_mas_watveg  , dSRF_mas_watveg  /dtime_sec*1E-9, dSRF_mas_watveg  /ATM_aire_sea_tot*1E-3) ) 
    611 echo ( 'dMass (watsoil)   = {:12.3e} kg | {:12.2e} Sv | {:12.2e} m '.format (dSRF_mas_watsoil , dSRF_mas_watsoil /dtime_sec*1E-9, dSRF_mas_watsoil /ATM_aire_sea_tot*1E-3) ) 
    612 echo ( 'dMass (sno)       = {:12.3e} kg | {:12.2e} Sv | {:12.2e} m '.format (dSRF_mas_snow    , dSRF_mas_snow    /dtime_sec*1E-9, dSRF_mas_snow    /ATM_aire_sea_tot*1E-3) ) 
     852echo ( 'SRF_mas_lake_beg     = {:12.6e} kg | SRF_mas_lake_end     = {:12.6e} kg '.format (SRF_mas_lake_beg    , SRF_mas_lake_end     ) ) 
     853 
     854prtFlux ('dMass (watveg) ', dSRF_mas_watveg , 'e' , True) 
     855prtFlux ('dMass (watsoil)', dSRF_mas_watsoil, 'e' , True) 
     856prtFlux ('dMass (snow)   ', dSRF_mas_snow   , 'e' , True) 
     857prtFlux ('dMass (lake)   ', dSRF_mas_lake   , 'e' , True) 
    613858 
    614859SRF_mas_wat_beg = SRF_mas_watveg_beg + SRF_mas_watsoil_beg + SRF_mas_snow_beg  
    615860SRF_mas_wat_end = SRF_mas_watveg_end + SRF_mas_watsoil_end + SRF_mas_snow_end  
    616 dSRF_mas_wat = SRF_mas_wat_end - SRF_mas_wat_beg 
    617  
    618 echo ( '\nWater content in surface ' ) 
    619 echo ( 'SRF_mas_wat_beg  = {:12.6e} kg | SRF_mas_wat_end  = {:12.6e} kg '.format (SRF_mas_wat_beg, SRF_mas_wat_end) ) 
    620 echo ( 'dMass (water srf) = {:12.3e} kg '.format (dSRF_mas_wat) ) 
    621 echo ( 'dMass (water srf) = {:12.3e} Sv '.format (dSRF_mas_wat/dtime_sec*1E-6/ATM_rho) ) 
    622 echo ( 'dMass (water srf) = {:12.3e} m  '.format (dSRF_mas_wat/ATM_aire_sea_tot/ATM_rho) ) 
    623  
    624 echo ( '\nWater content in  ATM + SRF + RUN ' ) 
     861dSRF_mas_wat    = SRF_mas_wat_end - SRF_mas_wat_beg 
     862 
     863echo ( '------------------------------------------------------------------------------------' ) 
     864echo ( f'Water content in surface  -- {Title} ' ) 
     865echo ( 'SRF_mas_wat_beg   = {:12.6e} kg | SRF_mas_wat_end  = {:12.6e} kg '.format (SRF_mas_wat_beg, SRF_mas_wat_end) ) 
     866prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', True ) 
     867 
     868echo ( '------------------------------------------------------------------------------------' ) 
     869echo ( 'Water content in  ATM + SRF + RUN + LAKE' ) 
    625870echo ( 'mas_wat_beg = {:12.6e} kg | mas_wat_end = {:12.6e} kg '. 
    626            format (DYN_mas_wat_beg + ATM_mas_sno_beg + RUN_mas_wat_beg + SRF_mas_wat_beg, 
    627                    DYN_mas_wat_end + ATM_mas_sno_end + RUN_mas_wat_end + SRF_mas_wat_end) ) 
    628 echo ( 'dMass (water atm+srf+run) = {:12.6e} kg '.format ( dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat) ) 
    629 echo ( 'dMass (water atm+srf+run) = {:12.3e} Sv '.format ((dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat)/dtime_sec*1E-6/ATM_rho) ) 
    630 echo ( 'dMass (water atm+srf+run) = {:12.3e} m  '.format ((dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat)/ATM_aire_sea_tot/ATM_rho) ) 
    631  
    632 echo ( '\n------------------------------------------------------------------------------------' ) 
    633 echo ( '-- ATM Fluxes ' ) 
     871           format (DYN_mas_wat_beg + ATM_mas_sno_beg + RUN_mas_wat_beg + SRF_mas_wat_beg + SRF_mas_lake_beg , 
     872                   DYN_mas_wat_end + ATM_mas_sno_end + RUN_mas_wat_end + SRF_mas_wat_end + SRF_mas_lake_end ) ) 
     873prtFlux ( 'dMass (water atm+srf+run+lake)', dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat + dSRF_mas_lake, 'e', True) 
     874 
     875echo ( '\n====================================================================================' ) 
     876echo ( f'-- ATM Fluxes  -- {Title} ' ) 
    634877 
    635878ATM_wbilo_oce   = lmdz.geo2point ( d_ATM_his ['wbilo_oce']   ) 
    … …  
    643886ATM_snowf       = lmdz.geo2point ( d_ATM_his ['snow']        ) 
    644887ATM_evap        = lmdz.geo2point ( d_ATM_his ['evap']        ) 
    645 ATM_bilo = ATM_wbilo_oce + ATM_wbilo_ter + ATM_wbilo_lic + ATM_wbilo_sic 
     888ATM_wbilo       = ATM_wbilo_oce + ATM_wbilo_sic + ATM_wbilo_ter + ATM_wbilo_lic 
     889ATM_emp         = ATM_evap - ATM_precip 
    646890 
    647891RUN_coastalflow = lmdz.geo2point ( d_RUN_his ['coastalflow'] )  
    … …  
    651895RUN_riversret   = lmdz.geo2point ( d_RUN_his ['riversret']   ) 
    652896 
    653 SRF_evap     = lmdz.geo2point ( d_SRF_his ['evap'] ) 
     897RUN_coastalflow_cpl = lmdz.geo2point ( d_RUN_his ['coastalflow_cpl'] )  
     898RUN_riverflow_cpl   = lmdz.geo2point ( d_RUN_his ['riverflow_cpl']   ) 
     899 
     900SRF_rain     = lmdz.geo2point ( d_SRF_his ['rain']  ) 
     901SRF_evap     = lmdz.geo2point ( d_SRF_his ['evap']  ) 
    654902SRF_snowf    = lmdz.geo2point ( d_SRF_his ['snowf'] ) 
    655 SRF_TWS      = lmdz.geo2point ( d_SRF_his ['TWS'] ) 
    656903SRF_subli    = lmdz.geo2point ( d_SRF_his ['subli'] ) 
    657904SRF_transpir = lmdz.geo2point ( np.sum(d_SRF_his ['transpir'], axis=1) ) ; SRF_transpir.attrs['units'] = d_SRF_his ['transpir'].attrs['units'] 
    658  
    659 def mmd2SI ( Var) : 
     905SRF_emp      = SRF_evap - SRF_rain - SRF_snowf ; SRF_emp.attrs['units'] = SRF_rain.attrs['units'] 
     906 
     907## Correcting units of SECHIBA variables 
     908def mmd2SI ( Var ) : 
    660909    '''Change unit from mm/d or m^3/s to kg/s if needed''' 
    661910    if 'units' in VarT.attrs :  
    … …  
    665914            VarT.values = VarT.values  * ATM_rho * (1e-3/86400.) ;  VarT.attrs['units'] = 'kg/s' 
    666915 
    667 for var in ['runoff', 'drainage', 'riversret', 'coastalflow', 'riverflow'] : 
     916for var in [ 'runoff', 'drainage', 'riversret', 'coastalflow', 'riverflow', 'coastalflow_cpl', 'riverflow_cpl' ] : 
    668917    VarT = locals()['RUN_' + var] 
    669918    mmd2SI (VarT) 
    670919 
    671 for var in ['evap', 'snowf', 'TWS', 'subli', 'transpir'] : 
     920for var in ['evap', 'snowf', 'subli', 'transpir', 'rain', 'emp' ] : 
    672921    VarT = locals()['SRF_' + var] 
    673922    mmd2SI (VarT) 
    674923 
     924 
    675925RUN_input  = RUN_runoff      + RUN_drainage 
    676926RUN_output = RUN_coastalflow + RUN_riverflow 
    … …  
    678928ATM_wbilo_sea = ATM_wbilo_oce + ATM_wbilo_sic 
    679929 
    680 ATM_flx_oce       = ATM_flux_int ( ATM_wbilo_oce  ) 
    681 ATM_flx_sic       = ATM_flux_int ( ATM_wbilo_sic  ) 
    682 ATM_flx_sea       = ATM_flux_int ( ATM_wbilo_sea  ) 
    683 ATM_flx_ter       = ATM_flux_int ( ATM_wbilo_ter  ) 
    684 ATM_flx_lic       = ATM_flux_int ( ATM_wbilo_lic  ) 
    685 ATM_flx_calving   = ATM_flux_int ( ATM_fqcalving  ) 
    686 ATM_flx_qfonte    = ATM_flux_int ( ATM_fqfonte    ) 
    687 ATM_flx_precip    = ATM_flux_int ( ATM_precip     ) 
    688 ATM_flx_snowf     = ATM_flux_int ( ATM_snowf      ) 
    689 ATM_flx_evap      = ATM_flux_int ( ATM_evap       ) 
    690 ATM_flx_runlic    = ATM_flux_int ( ATM_runofflic  ) 
    691  
    692 RUN_flx_coastal   = ONE_flux_int ( RUN_coastalflow) 
    693 RUN_flx_river     = ONE_flux_int ( RUN_riverflow  ) 
    694 RUN_flx_drainage  = ATM_flux_int ( RUN_drainage   ) 
    695 RUN_flx_riversret = ATM_flux_int ( RUN_riversret  ) 
    696 RUN_flx_runoff    = ATM_flux_int ( RUN_runoff     ) 
    697 RUN_flx_input     = ATM_flux_int ( RUN_input      ) 
    698 RUN_flx_output    = ONE_flux_int ( RUN_output     ) 
    699  
    700 ATM_flx_emp   = ATM_flx_evap - ATM_flx_precip 
    701 ATM_flx_bilo  = ATM_flx_oce + ATM_flx_sic + ATM_flx_ter + ATM_flx_lic 
     930ATM_flx_oce         = ATM_flux_int ( ATM_wbilo_oce  ) 
     931ATM_flx_sic         = ATM_flux_int ( ATM_wbilo_sic  ) 
     932ATM_flx_sea         = ATM_flux_int ( ATM_wbilo_sea  ) 
     933ATM_flx_ter         = ATM_flux_int ( ATM_wbilo_ter  ) 
     934ATM_flx_lic         = ATM_flux_int ( ATM_wbilo_lic  ) 
     935ATM_flx_calving     = ATM_flux_int ( ATM_fqcalving  ) 
     936ATM_flx_qfonte      = ATM_flux_int ( ATM_fqfonte    ) 
     937ATM_flx_precip      = ATM_flux_int ( ATM_precip     ) 
     938ATM_flx_snowf       = ATM_flux_int ( ATM_snowf      ) 
     939ATM_flx_evap        = ATM_flux_int ( ATM_evap       ) 
     940ATM_flx_runlic      = ATM_flux_int ( ATM_runofflic  ) 
     941 
     942RUN_flx_coastal     = ONE_flux_int ( RUN_coastalflow) 
     943RUN_flx_river       = ONE_flux_int ( RUN_riverflow  ) 
     944RUN_flx_coastal_cpl = ONE_flux_int ( RUN_coastalflow_cpl) 
     945RUN_flx_river_cpl   = ONE_flux_int ( RUN_riverflow_cpl  ) 
     946RUN_flx_drainage    = SRF_flux_int ( RUN_drainage   ) 
     947RUN_flx_riversret   = SRF_flux_int ( RUN_riversret  ) 
     948RUN_flx_runoff      = SRF_flux_int ( RUN_runoff     ) 
     949RUN_flx_input       = SRF_flux_int ( RUN_input      ) 
     950RUN_flx_output      = ONE_flux_int ( RUN_output     ) 
     951 
     952ATM_flx_emp   = ATM_flux_int (ATM_emp) 
     953ATM_flx_wbilo = ATM_flux_int (ATM_wbilo) 
     954 
    702955RUN_flx_bil   = RUN_flx_input - RUN_flx_output 
    703  
    704956RUN_flx_rivcoa = RUN_flx_coastal + RUN_flx_river 
    705957 
    706 ATM_flx_bilo2 = ATM_flux_int (ATM_bilo) 
    707  
    708  
    709 echo ('  wbilo_oce     {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_oce      , ATM_flx_oce      / dtime_sec*1E-6/ATM_rho, ATM_flx_oce      /ATM_aire_sea_tot/ATM_rho )) 
    710 echo ('  wbilo_sic     {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_sic      , ATM_flx_sic      / dtime_sec*1E-6/ATM_rho, ATM_flx_sic      /ATM_aire_sea_tot/ATM_rho )) 
    711 echo ('  wbilo_sic+oce {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_sea      , ATM_flx_sea      / dtime_sec*1E-6/ATM_rho, ATM_flx_sea      /ATM_aire_sea_tot/ATM_rho )) 
    712 echo ('  wbilo_ter     {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_ter      , ATM_flx_ter      / dtime_sec*1E-6/ATM_rho, ATM_flx_ter      /ATM_aire_sea_tot/ATM_rho )) 
    713 echo ('  wbilo_lic     {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_lic      , ATM_flx_lic      / dtime_sec*1E-6/ATM_rho, ATM_flx_lic      /ATM_aire_sea_tot/ATM_rho )) 
    714 echo ('  Sum wbilo_*   {:12.3e} (kg) | {:12.4e} (Sv) | {:12.4f} m '.format ( ATM_flx_bilo     , ATM_flx_bilo     / dtime_sec*1E-6/ATM_rho, ATM_flx_bilo     /ATM_aire_sea_tot/ATM_rho )) 
    715 echo ('  E-P           {:12.3e} (kg) | {:12.4e} (Sv) | {:12.4f} m '.format ( ATM_flx_emp      , ATM_flx_emp      / dtime_sec*1E-6/ATM_rho, ATM_flx_emp      /ATM_aire_sea_tot/ATM_rho )) 
    716 echo ('  calving       {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_calving  , ATM_flx_calving  / dtime_sec*1E-6/ATM_rho, ATM_flx_calving  /ATM_aire_sea_tot/ATM_rho )) 
    717 echo ('  fqfonte       {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_qfonte   , ATM_flx_qfonte   / dtime_sec*1E-6/ATM_rho, ATM_flx_qfonte   /ATM_aire_sea_tot/ATM_rho )) 
    718 echo ('  precip        {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_precip   , ATM_flx_precip   / dtime_sec*1E-6/ATM_rho, ATM_flx_precip   /ATM_aire_sea_tot/ATM_rho )) 
    719 echo ('  snowf         {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_snowf    , ATM_flx_snowf    / dtime_sec*1E-6/ATM_rho, ATM_flx_snowf    /ATM_aire_sea_tot/ATM_rho )) 
    720 echo ('  evap          {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_evap     , ATM_flx_evap     / dtime_sec*1E-6/ATM_rho, ATM_flx_evap     /ATM_aire_sea_tot/ATM_rho )) 
    721 echo ('  coastalflow   {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( RUN_flx_coastal  , RUN_flx_coastal  / dtime_sec*1E-6/ATM_rho, RUN_flx_coastal  /ATM_aire_sea_tot/ATM_rho )) 
    722 echo ('  riverflow     {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( RUN_flx_river    , RUN_flx_river    / dtime_sec*1E-6/ATM_rho, RUN_flx_river    /ATM_aire_sea_tot/ATM_rho )) 
    723 echo ('  river+coastal {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( RUN_flx_rivcoa   , RUN_flx_rivcoa   / dtime_sec*1E-6/ATM_rho, RUN_flx_rivcoa   /ATM_aire_sea_tot/ATM_rho )) 
    724 echo ('  drainage      {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( RUN_flx_drainage , RUN_flx_drainage / dtime_sec*1E-6/ATM_rho, RUN_flx_drainage /ATM_aire_sea_tot/ATM_rho )) 
    725 echo ('  riversret     {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( RUN_flx_riversret, RUN_flx_riversret/ dtime_sec*1E-6/ATM_rho, RUN_flx_riversret/ATM_aire_sea_tot/ATM_rho )) 
    726 echo ('  runoff        {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( RUN_flx_runoff   , RUN_flx_runoff   / dtime_sec*1E-6/ATM_rho, RUN_flx_runoff   /ATM_aire_sea_tot/ATM_rho )) 
    727 echo ('  river in      {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( RUN_flx_input    , RUN_flx_input    / dtime_sec*1E-6/ATM_rho, RUN_flx_input    /ATM_aire_sea_tot/ATM_rho )) 
    728 echo ('  river out     {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( RUN_flx_output   , RUN_flx_output   / dtime_sec*1E-6/ATM_rho, RUN_flx_output   /ATM_aire_sea_tot/ATM_rho )) 
    729 echo ('  river bil     {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( RUN_flx_bil      , RUN_flx_bil      / dtime_sec*1E-6/ATM_rho, RUN_flx_bil      /ATM_aire_sea_tot/ATM_rho )) 
    730 echo ('  runoff lic    {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_runlic   , ATM_flx_runlic   / dtime_sec*1E-6/ATM_rho, ATM_flx_runlic   /ATM_aire_sea_tot/ATM_rho )) 
    731  
    732 ATM_flx_budget = -ATM_flx_sea + ATM_flx_calving + ATM_flx_runlic #+ ATM_flx_qfonte + RUN_flx_river  
    733 echo ('\n  Global      {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( ATM_flx_budget , ATM_flx_budget / dtime_sec*1E-9, ATM_flx_budget /ATM_aire_sea_tot/ATM_rho )) 
    734  
    735 #echo ('  E-P-R         {:12.3e} (kg) | {:12.4e} (Sv) | {:12.4f} m '.format ( ATM_flx_emp      , ATM_flx_emp      / dtime_sec*1E-6/ATM_rho, ATM_flx_emp      /ATM_aire_sea_tot/ATM_rho )) 
    736  
    737  
    738 echo ( '------------------------------------------------------------------------------------' ) 
    739 echo ( '-- SECHIBA fluxes' ) 
    740  
    741  
    742 SRF_flx_evap     = ATM_flux_int ( SRF_evap     ) 
    743 SRF_flx_snowf    = ATM_flux_int ( SRF_snowf    ) 
    744 SRF_flx_subli    = ATM_flux_int ( SRF_subli    ) 
    745 SRF_flx_transpir = ATM_flux_int ( SRF_transpir ) 
    746  
    747 echo ('  evap     {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( SRF_flx_evap      , SRF_flx_evap      / dtime_sec*1E-6/ATM_rho, SRF_flx_evap      /ATM_aire_sea_tot/ATM_rho )) 
    748 echo ('  snowf    {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( SRF_flx_snowf     , SRF_flx_snowf     / dtime_sec*1E-6/ATM_rho, SRF_flx_snowf     /ATM_aire_sea_tot/ATM_rho )) 
    749 echo ('  subli    {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( SRF_flx_subli     , SRF_flx_subli     / dtime_sec*1E-6/ATM_rho, SRF_flx_subli     /ATM_aire_sea_tot/ATM_rho )) 
    750 echo ('  transpir {:12.3e} (kg) | {:12.4f} (Sv) | {:12.4f} m '.format ( SRF_flx_transpir  , SRF_flx_transpir  / dtime_sec*1E-6/ATM_rho, SRF_flx_transpir  /ATM_aire_sea_tot/ATM_rho )) 
     958prtFlux ('wbilo_oce     ', ATM_flx_oce, 'f' ) 
     959prtFlux ('wbilo_oce     ', ATM_flx_oce, 'f' )           
     960prtFlux ('wbilo_sic     ', ATM_flx_sic, 'f' )           
     961prtFlux ('wbilo_sic+oce ', ATM_flx_sea, 'f' )           
     962prtFlux ('wbilo_ter     ', ATM_flx_ter, 'f' )           
     963prtFlux ('wbilo_lic     ', ATM_flx_lic, 'f' )           
     964prtFlux ('Sum wbilo_*   ', ATM_flx_wbilo      , 'f', True)   
     965prtFlux ('E-P           ', ATM_flx_emp        , 'f', True)   
     966prtFlux ('calving       ', ATM_flx_calving    , 'f' )  
     967prtFlux ('fqfonte       ', ATM_flx_qfonte     , 'f' )        
     968prtFlux ('precip        ', ATM_flx_precip     , 'f' )        
     969prtFlux ('snowf         ', ATM_flx_snowf      , 'f' )         
     970prtFlux ('evap          ', ATM_flx_evap       , 'f' )          
     971prtFlux ('coastalflow   ', RUN_flx_coastal    , 'f' )  
     972prtFlux ('riverflow     ', RUN_flx_river      , 'f' )         
     973prtFlux ('coastal_cpl   ', RUN_flx_coastal_cpl, 'f' )   
     974prtFlux ('riverf_cpl    ', RUN_flx_river_cpl  , 'f' )     
     975prtFlux ('river+coastal ', RUN_flx_rivcoa     , 'f' )    
     976prtFlux ('drainage      ', RUN_flx_drainage   , 'f' )    
     977prtFlux ('riversret     ', RUN_flx_riversret  , 'f' )    
     978prtFlux ('runoff        ', RUN_flx_runoff     , 'f' )    
     979prtFlux ('river in      ', RUN_flx_input      , 'f' )    
     980prtFlux ('river out     ', RUN_flx_output     , 'f' )    
     981prtFlux ('river bil     ', RUN_flx_bil        , 'f' )    
     982prtFlux ('runoff lic    ', ATM_flx_runlic     , 'f' )    
     983 
     984ATM_flx_budget = -ATM_flx_sea + ATM_flx_calving + ATM_flx_runlic #+ ATM_flx_qfonte + RUN_flx_river 
     985echo ('') 
     986echo ('  Global        {:12.3e} kg | {:12.4f} Sv | {:12.4f} m '.format ( ATM_flx_budget , ATM_flx_budget / dtime_sec*1E-9, ATM_flx_budget /ATM_aire_sea_tot/ATM_rho )) 
     987 
     988#echo ('  E-P-R         {:12.3e} kg | {:12.4e} Sv | {:12.4f} m '.format ( ATM_flx_emp      , ATM_flx_emp      / dtime_sec*1E-6/ATM_rho, ATM_flx_emp      /ATM_aire_sea_tot/ATM_rho )) 
     989 
     990echo (' ') 
     991echo ( '\n====================================================================================' ) 
     992echo ( f'--  Atmosphere  -- {Title} ' ) 
     993echo ( 'Mass begin = {:12.6e} kg | Mass end = {:12.6e} kg'.format (DYN_mas_wat_beg, DYN_mas_wat_end) ) 
     994echo ( 'dmass (atm)  = {:12.3e} kg | {:12.4e} Sv | {:12.4f} mm/year '.format ( dDYN_mas_wat , kg2Sv(dDYN_mas_wat) , kg2myear(dDYN_mas_wat )*1000. )) 
     995echo ( 'Sum wbilo_*  = {:12.3e} kg | {:12.4e} Sv | {:12.4f} mm/year '.format ( ATM_flx_wbilo, kg2Sv(ATM_flx_wbilo), kg2myear(ATM_flx_wbilo )*1000. )) 
     996echo ( 'E-P          = {:12.3e} kg | {:12.4e} Sv | {:12.4f} mm/year '.format ( ATM_flx_emp  , kg2Sv(ATM_flx_emp ) , kg2myear(ATM_flx_emp  )*1000. )) 
     997echo ( ' ' ) 
     998prtFlux ( 'Perte eau atm ', ATM_flx_wbilo - dDYN_mas_wat, 'e', True ) 
     999echo ( 'Perte eau atm = {:12.3e} (rel)  '.format ( (ATM_flx_wbilo - dDYN_mas_wat)/dDYN_mas_wat  ) ) 
     1000 
     1001echo (' ') 
     1002prtFlux ( 'Perte eau atm', ATM_flx_emp  - dDYN_mas_wat , 'e', True ) 
     1003echo ( 'Perte eau atm = {:12.3e} (rel)  '.format ( (ATM_flx_emp-dDYN_mas_wat)/dDYN_mas_wat  ) ) 
     1004echo (' ') 
     1005 
     1006echo ( '\n====================================================================================' ) 
     1007echo ( f'-- SECHIBA fluxes  -- {Title} ' ) 
     1008 
     1009SRF_flx_rain     = SRF_flux_int ( SRF_rain     ) 
     1010SRF_flx_evap     = SRF_flux_int ( SRF_evap     ) 
     1011SRF_flx_snowf    = SRF_flux_int ( SRF_snowf    ) 
     1012SRF_flx_subli    = SRF_flux_int ( SRF_subli    ) 
     1013SRF_flx_transpir = SRF_flux_int ( SRF_transpir ) 
     1014SRF_flx_emp      = SRF_flux_int ( SRF_emp      ) 
     1015 
     1016RUN_flx_torouting  =  RUN_flx_runoff  + RUN_flx_drainage 
     1017RUN_flx_fromrouting = RUN_flx_coastal + RUN_flx_river 
     1018 
     1019SRF_flx_all = SRF_flx_rain + SRF_flx_snowf - SRF_flx_evap - RUN_flx_runoff - RUN_flx_drainage 
     1020 
     1021prtFlux ('rain       ', SRF_flx_rain     , 'f' ) 
     1022prtFlux ('evap       ', SRF_flx_evap     , 'f' ) 
     1023prtFlux ('snowf      ', SRF_flx_snowf    , 'f' ) 
     1024prtFlux ('E-P        ', SRF_flx_emp      , 'f' ) 
     1025prtFlux ('subli      ', SRF_flx_subli    , 'f' ) 
     1026prtFlux ('transpir   ', SRF_flx_transpir , 'f' ) 
     1027prtFlux ('to routing ', RUN_flx_torouting, 'f' ) 
     1028prtFlux ('budget     ', SRF_flx_all      , 'f', True ) 
     1029 
     1030echo ( '\n------------------------------------------------------------------------------------' ) 
     1031echo ( 'Water content in surface ' ) 
     1032echo ( 'SRF_mas_wat_beg  = {:12.6e} kg | SRF_mas_wat_end  = {:12.6e} kg '.format (SRF_mas_wat_beg, SRF_mas_wat_end) ) 
     1033prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', True) 
     1034echo ( 'dMass (water srf) = {:12.4e} (rel)   '.format ( (SRF_flx_all-dSRF_mas_wat)/dSRF_mas_wat) ) 
     1035 
     1036echo ('') 
     1037echo ( '\n====================================================================================' ) 
     1038echo ( f'-- RUNOFF fluxes  -- {Title} ' ) 
     1039RUN_flx_all = RUN_flx_torouting - RUN_flx_river - RUN_flx_coastal 
     1040prtFlux ('runoff    ', RUN_flx_runoff      , 'f' )  
     1041prtFlux ('drainage  ', RUN_flx_drainage    , 'f' )  
     1042prtFlux ('run+drain ', RUN_flx_torouting   , 'f' )  
     1043prtFlux ('river     ', RUN_flx_river       , 'f' )  
     1044prtFlux ('coastal   ', RUN_flx_coastal     , 'f' )  
     1045prtFlux ('riv+coa   ', RUN_flx_fromrouting , 'f' )  
     1046prtFlux ('budget    ', RUN_flx_all         , 'f' , True) 
     1047 
     1048echo ( '\n------------------------------------------------------------------------------------' ) 
     1049echo ( f'Water content in routing  -- {Title} ' ) 
     1050echo ( 'RUN_mas_wat_beg  = {:12.6e} kg | RUN_mas_wat_end = {:12.6e} kg '.format (RUN_mas_wat_beg, RUN_mas_wat_end) ) 
     1051prtFlux ( 'dMass (routing)  ', dRUN_mas_wat+dSRF_mas_lake, 'f', True) 
     1052 
     1053echo ( '\n------------------------------------------------------------------------------------' ) 
     1054echo ( f'Water content in routing  -- {Title} ' ) 
     1055echo ( 'RUN_mas_wat_beg  = {:12.6e} kg | RUN_mas_wat_end = {:12.6e} kg '.format (RUN_mas_wat_beg, RUN_mas_wat_end) ) 
     1056prtFlux ( 'dMass (routing)  ', dRUN_mas_wat, 'f', True  ) 
     1057 
     1058print ( 'Routing error : {:12.3e} (rel)'.format ( (RUN_flx_all-dRUN_mas_wat)/dRUN_mas_wat ) ) 
     1059 
     1060# PRINT *,"routing water Balance ;  before : ", water_balance_before," ; after : ",water_balance_after,  &a 
     1061# 1150              " ; delta : ", 100*(water_balance_after-water_balance_before)/(0.5*(water_balance_after+water_balance_before)),"%" 

  • TOOLS/WATER_BUDGET/CPL_waterbudget.py

    r6271 r6277  
    692692    RUN_mas_wat_end = ONE_stock_int ( d_RUN_end ['fast_reservoir'] +  d_RUN_end ['slow_reservoir'] +  d_RUN_end ['stream_reservoir']) 
    693693 
    694 if Routing == 'ORCHIDEE' :  
     694if Routing == 'SECHIBA' :  
    695695    RUN_mas_wat_beg = ONE_stock_int ( d_SRF_beg['fastres']  + d_SRF_beg['slowres'] + d_SRF_beg['streamres'] \ 
    696696                                    + d_SRF_beg['floodres'] + d_SRF_beg['lakeres'] + d_SRF_beg['pondres'] ) 

  • TOOLS/WATER_BUDGET/OCE_waterbudget.py

    r6271 r6277  
    2121### 
    2222## Import system modules 
    23 import sys, os, shutil, subprocess 
     23import sys, os, shutil, subprocess, platform 
    2424import numpy as np 
    2525import configparser, re 
    … …  
    2929config.optionxform = str # To keep capitals 
    3030 
    31 config['Files']  = {} 
    32 config['System'] = {} 
    33  
    34 ##-- Some physical constants 
    35 #-- Earth Radius 
    36 Ra = 40.e6 / (2. * np.pi) 
    37 #-- Gravity 
    38 g = 9.81 
    39 #-- Ice density (kg/m3) in LIM3 
    40 ICE_rho_ice = 917.0 
    41 #-- Snow density (kg/m3) in LIM3 
    42 ICE_rho_sno = 330.0 
    43 #-- Ocean water density (kg/m3) in NEMO 
    44 OCE_rho_liq = 1026. 
    45 #-- Water density in ice pounds in SI3 
    46 ICE_rho_pnd = 1000. 
    47  
    4831## Read experiment parameters 
    49 ATM = None ; ORCA = None ; NEMO = None ; OCE_relax = False ; OCE_icb = False ; Coupled = False ; Routing = None 
     32ATM=None ; ORCA=None ; NEMO=None ; OCE_relax=False ; OCE_icb=False ; Coupled=False ; Routing=None ;  TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=None 
     33 
     34## 
     35ARCHIVE=None ; STORAGE = None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None 
     36TmpDir=None ; FileOut=None ; dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None ; FileCommon=None 
     37file_ATM_his=None ; file_SRF_his=None ; file_RUN_his=None ; file_OCE_his=None ;  file_ICE_his=None ; file_OCE_sca=None  
     38file_restart_beg=None ; file_restart_end=None ; file_ATM_beg=None ; file_ATM_end=None ; file_DYN_beg=None ; file_DYN_end=None ; file_SRF_beg=None ; file_SRF_end=None 
     39file_RUN_beg=None ; file_RUN_end=None ; file_OCE_beg=None ; file_ICE_beg=None ; file_OCE_end=None ; file_ICE_end=None 
    5040 
    5141# Arguments passed 
    5242print ( "Name of Python script:", sys.argv[0] ) 
    53 IniFile =  sys.argv[1] 
     43IniFile = sys.argv[1] 
    5444print ("Input file : ", IniFile ) 
    5545config.read (IniFile) 
     46if 'full' in IniFile : FullIniFile = IniFile 
     47else                 : FullIniFile = 'full_' + IniFile 
    5648 
    5749def setBool (chars) : 
    … …  
    7567    return setNum 
    7668 
    77 print ('[Experiment]') 
    78 for VarName in config['Experiment'].keys() : 
    79     locals()[VarName] = config['Experiment'][VarName] 
    80     locals()[VarName] = setBool (locals()[VarName]) 
    81     locals()[VarName] = setNum  (locals()[VarName]) 
    82     print ( '{:25} set to : {:}'.format (VarName, locals()[VarName]) ) 
     69def setNone (chars) : 
     70    '''Convert specific char string to None if possible''' 
     71    if type (chars) == str : 
     72        if chars in ['None', 'NONE', 'none'] : setNone = None 
     73        else : setNone = chars 
     74    else : setNone = chars 
     75    return setNone 
     76 
     77## Reading config 
     78for Section in ['Experiment', 'libIGCM', 'Files', 'Physics' ] : 
     79    if Section in config.keys() :  
     80        print ( f'[{Section}]' ) 
     81        for VarName in config[Section].keys() : 
     82            locals()[VarName] = config[Section][VarName] 
     83            locals()[VarName] = setBool (locals()[VarName]) 
     84            locals()[VarName] = setNum  (locals()[VarName]) 
     85            print ( '{:25} set to : {:}'.format (VarName, locals()[VarName]) ) 
     86 
     87if not 'Files' in config.keys() : config['Files'] = {} 
     88 
     89def unDefined (char) : 
     90    if char in globals () : 
     91        if char == None : return True 
     92        else : return False 
     93    else : return True 
     94 
     95##-- Some physical constants 
     96#-- Earth Radius 
     97if not 'Ra'            in locals () : Ra = 6366197.7236758135 
     98#-- Gravity 
     99if not 'Grav'          in locals () : Grav = 9.81 
     100#-- Ice volumic mass (kg/m3) in LIM3 
     101if not 'ICE_rho_ice'   in locals () : ICE_rho_ice = 917.0 
     102#-- Snow volumic mass (kg/m3) in LIM3 
     103if not 'ICE_rho_sno'   in locals () : ICE_rho_sno = 330.0 
     104    #-- Water density in ice pounds in SI3 
     105if not 'ICE_rho_pnd'   in locals () : ICE_rho_pnd = 1000. 
     106#-- Ocean water volumic mass (kg/m3) in NEMO 
     107if not 'OCE_rho_liq'   in locals () : OCE_rho_liq = 1026. 
     108#-- Water volumic mass in atmosphere 
     109if not 'ATM_rho'       in locals () : ATM_rho = 1000. 
     110#-- Water volumic mass in surface reservoirs 
     111if not 'SRF_rho'       in locals () : SRF_rho = 1000. 
     112#-- Water volumic mass of rivers 
     113if not 'RUN_rho'       in locals () : RUN_rho = 1000. 
     114#-- Year length 
     115if not 'YearLength'    in locals () : YearLength = 365.25 * 24. * 60. * 60. 
     116 
     117config['Physics'] = { 'Ra':Ra, 'Grav':Grav, 'ICE_rho_ice':ICE_rho_ice, 'ICE_rho_sno':ICE_rho_sno, 'OCE_rho_liq':OCE_rho_liq, 'ICE_rho_pnd':ICE_rho_pnd, 
     118                          'ATM_rho':ATM_rho, 'SRF_rho':SRF_rho, 'RUN_rho':RUN_rho} 
     119         
     120 
    83121 
    84122# Where do we run ? 
    … …  
    109147if IDRIS : 
    110148    raise Exception ("Pour IDRIS : repertoires et chemins a definir")  
    111      
     149 
     150config['System'] = {'SysName':SysName, 'NodeName':NodeName, 'Release':Release, 'Version':Version,'Machine':Machine, 'TGCC':TGCC,'IDRIS':IDRIS, 'CPU':CPU, 
     151                    'Program'  : "Generated by : " + str(sys.argv),  
     152                    'HOSTNAME' : platform.node (), 'LOGNAME'  : os.getlogin (), 
     153                    'Python'   : f'{platform.python_version ()}', 
     154                    'OS'       : f'{platform.system ()} release : {platform.release ()} hardware : {platform.machine ()}', 
     155                    'SVN_Author'   : "$Author$",  
     156                    'SVN_Date'     : "$Date$", 
     157                    'SVN_Revision' : "$Revision$", 
     158                    'SVN_Id'       : "$Id$", 
     159                    'SVN_HeadURL'  : "$HeadURL$"} 
     160 
     161if libIGCM : 
     162    config['libIGCM'] = {'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'SCRATCHDIR':SCRATCHDIR, 'R_IN':R_IN, 'rebuild':rebuild } 
    112163 
    113164## Import specific module 
    … …  
    117168     
    118169# Output file 
    119 FileOut = f'OCE_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out' 
     170if FileOut == None :  
     171    FileOut = f'OCE_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out' 
     172    config['Files']['FileOut'] = FileOut 
     173 
    120174f_out = open ( FileOut, mode = 'w' ) 
    121  
     175     
    122176# Function to print to stdout *and* output file 
    123 def echo (string) : 
    124     print ( string  ) 
    125     f_out.write ( string + '\n' ) 
     177def echo (string, end='\n') : 
     178    print ( str(string), end=end  ) 
     179    sys.stdout.flush () 
     180    f_out.write ( str(string) + end ) 
    126181    f_out.flush () 
    127182    return None 
    … …  
    146201if not os.path.exists (TmpDirICE) : os.mkdir (TmpDirICE ) 
    147202 
     203echo (' ') 
     204echo ( f'JobName : {JobName}' ) 
     205echo (Comment) 
    148206echo ( f'Working in TMPDIR : {TmpDir}' ) 
    149207 
    … …  
    153211if Freq == "MO" : FreqDir =  os.path.join ('Output' , 'MO' ) 
    154212if Freq == "SE" : FreqDir =  os.path.join ('Analyse', 'SE' ) 
    155 dir_OCE_his = os.path.join ( R_SAVE, "OCE", FreqDir ) 
    156 dir_ICE_his = os.path.join ( R_SAVE, "ICE", FreqDir ) 
    157  
     213if dir_OCE_his == None : 
     214    dir_OCE_his = os.path.join ( R_SAVE, "OCE", FreqDir ) 
     215    config['Files']['dir_OCE_his'] = dir_OCE_his 
     216if dir_ICE_his == None : 
     217    dir_ICE_his = os.path.join ( R_SAVE, "ICE", FreqDir ) 
     218    config['Files']['dir_OCE_his'] = dir_OCE_his 
     219  
    158220echo ( f'The analysis relies on files from the following model output directories : ' ) 
    159221echo ( f'{dir_OCE_his}' ) 
    … …  
    161223 
    162224#-- Files Names 
    163 if Freq == 'MO' : FileCommon = f'{JobName}_{YearBegin}0101_{YearEnd}1231_1M' 
    164 if Freq == 'SE' : FileCommon = f'{JobName}_SE_{YearBegin}0101_{YearEnd}1231_1M' 
    165  
     225if Period == None : 
     226    if Freq == 'MO' : Period = f'{YearBegin}0101_{YearEnd}1231_1M' 
     227    if Freq == 'SE' : Period = f'SE_{YearBegin}0101_{YearEnd}1231_1M' 
     228    config['Files']['Period'] = Period 
     229if FileCommon == None : 
     230    FileCommon = f'{JobName}_{Period}' 
     231    config['Files']['FileCommon'] = FileCommon 
     232 
     233if Title == None : 
     234    Title = f'{JobName} : {Freq} : {YearBegin}-01-01 - {YearEnd}-12-31' 
     235    config['Files']['Title'] = Title 
     236         
     237         
    166238echo ('\nOpen history files' ) 
    167 file_OCE_his = os.path.join ( dir_OCE_his, f'{FileCommon}_grid_T.nc' ) 
    168 file_OCE_sca = os.path.join ( dir_OCE_his, f'{FileCommon}_scalar.nc' ) 
    169 file_ICE_his = os.path.join ( dir_ICE_his, f'{FileCommon}_icemod.nc' ) 
     239if file_OCE_his == None : 
     240    file_OCE_his = os.path.join ( dir_OCE_his, f'{FileCommon}_grid_T.nc' ) 
     241    file_OCE_his = file_OCE_his 
     242if file_OCE_sca == None :     
     243    file_OCE_sca = os.path.join ( dir_OCE_his, f'{FileCommon}_scalar.nc' ) 
     244    config['Files']['file_OCE_sca'] = file_OCE_sca 
     245if file_ICE_his == None :  
     246    file_ICE_his = os.path.join ( dir_ICE_his, f'{FileCommon}_icemod.nc' ) 
     247    config['Files']['file_ICE_his'] = file_ICE_his 
    170248 
    171249d_OCE_his = xr.open_dataset ( file_OCE_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 
    … …  
    191269dtime_per_sec.attrs['unit'] = 's' 
    192270 
     271## Number of years 
     272NbYear = dtime_sec / YearLength 
    193273#-- Open restart files 
    194 YearRes = YearBegin - 1        # Year of the restart of beginning of simulation 
     274YearRes = YearBegin - 1              # Year of the restart of beginning of simulation 
    195275YearPre = YearBegin - PackFrequency  # Year to find the tarfile of the restart of beginning of simulation 
    196276 
    197277echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ') 
    198278 
    199 file_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{YearPre}0101_{YearRes}1231_restart.tar' ) 
    200 file_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{YearBegin}0101_{YearEnd}1231_restart.tar' ) 
     279if TarRestartPeriod_beg == None :  
     280    echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ') 
     281    TarRestartPeriod_beg = f'{YearPre}0101_{YearRes}1231' 
     282    config['Files']['TarRestartPeriod_beg'] = TarRestartPeriod_beg 
     283 
     284if TarRestartPeriod_end == None :  
     285    YearPre = YearBegin - PackFrequency  # Year to find the tarfile of the restart of beginning of simulation 
     286    echo (f'Restart dates - Start : {YearRes}-12-31  /  End : {YearEnd}-12-31 ') 
     287    TarRestartPeriod_end = f'{YearBegin}0101_{YearEnd}1231' 
     288    config['Files']['TarRestartPeriod_end'] = TarRestartPeriod_end 
     289 
     290if file_restart_beg == None : 
     291    file_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_beg}_restart.tar' ) 
     292    config['Files']['file_restart_beg'] = file_restart_beg 
     293if file_restart_end == None : 
     294    file_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_end}_restart.tar' ) 
     295    config['Files']['file_restart_end'] = file_restart_end 
    201296 
    202297echo ( f'{file_restart_beg}' ) 
    203298echo ( f'{file_restart_end}' ) 
    204299 
    205 file_OCE_beg = f'OCE_{JobName}_{YearRes}1231_restart.nc' 
    206 file_OCE_end = f'OCE_{JobName}_{YearEnd}1231_restart.nc' 
    207 file_ICE_beg = f'ICE_{JobName}_{YearRes}1231_restart_icemod.nc' 
    208 file_ICE_end = f'ICE_{JobName}_{YearEnd}1231_restart_icemod.nc' 
     300if file_OCE_beg == None :  
     301    file_OCE_beg = f'{TmpDir}/OCE_{JobName}_{YearRes}1231_restart.nc' 
     302    config['Files']['file_OCE_beg'] = file_OCE_beg 
     303if file_OCE_end == None : 
     304    file_OCE_end = f'{TmpDir}/OCE_{JobName}_{YearEnd}1231_restart.nc' 
     305    config['Files']['file_OCE_end'] = file_OCE_end 
     306if file_ICE_beg == None : 
     307    file_ICE_beg = f'{TmpDir}/ICE_{JobName}_{YearRes}1231_restart_icemod.nc' 
     308    config['Files']['file_ICE_beg'] = file_ICE_beg 
     309if file_ICE_end == None :  
     310    file_ICE_end = f'{TmpDir}/ICE_{JobName}_{YearEnd}1231_restart_icemod.nc' 
     311    config['Files']['file_ICE_end'] = file_ICE_end 
    209312 
    210313echo ( f'{file_OCE_beg}' ) 
    … …  
    221324    return int (ndomain_opa) 
    222325 
    223 if not os.path.exists ( os.path.join (TmpDir, file_OCE_beg) ) : 
    224     echo ( 'Extracting '+ os.path.join (TmpDir, file_OCE_beg) ) 
     326if not os.path.exists ( file_OCE_beg) : 
     327    echo ( f'Extracting {file_OCE_beg}' ) 
     328    base_file_OCE_beg = os.path.basename (file_OCE_beg) 
    225329    if not os.path.exists ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearRes}1231_restart_0000.nc') ) : 
    226330        command =  f'cd {TmpDirOCE} ; tar xf {file_restart_beg}  OCE_{JobName}_{YearRes}1231_restart_*.nc' 
    … …  
    229333    echo ('extract ndomain' ) 
    230334    ndomain_opa = get_ndomain ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearRes}1231_restart') ) 
    231     command = f'cd {TmpDirOCE} ; {rebuild} OCE_{JobName}_{YearRes}1231_restart {ndomain_opa} ; mv {file_OCE_beg} {TmpDir}' 
     335    command = f'cd {TmpDirOCE} ; {rebuild} OCE_{JobName}_{YearRes}1231_restart {ndomain_opa} ; mv {base_file_OCE_beg} {TmpDir}' 
    232336    echo ( command ) 
    233337    os.system ( command ) 
    234338    echo ( f'Rebuild done : {file_OCE_beg}') 
    235339     
    236 if not os.path.exists ( os.path.join (TmpDir, file_OCE_end) ) : 
     340if not os.path.exists ( file_OCE_end) : 
     341    echo ( f'Extracting {file_OCE_end}' ) 
     342    base_file_OCE_end = os.path.basename (file_OCE_end) 
    237343    if not os.path.exists ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearEnd}1231_restart_0000.nc') ): 
    238344        command =  f'cd {TmpDirOCE} ; tar xf {file_restart_end}  OCE_{JobName}_{YearEnd}1231_restart_*.nc' 
    … …  
    241347    echo ('extract ndomain' ) 
    242348    ndomain_opa = get_ndomain ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearEnd}1231_restart') ) 
    243     command = f'cd {TmpDirOCE} ; {rebuild} OCE_{JobName}_{YearEnd}1231_restart {ndomain_opa} ; mv {file_OCE_end} {TmpDir}' 
     349    command = f'cd {TmpDirOCE} ; {rebuild} OCE_{JobName}_{YearEnd}1231_restart {ndomain_opa} ; mv {base_file_OCE_end} {TmpDir}' 
    244350    echo ( command ) 
    245351    os.system ( command ) 
    246352    echo ( f'Rebuild done : {file_OCE_end}') 
    247353 
    248 if not os.path.exists ( os.path.join (TmpDir, file_ICE_beg) ) : 
     354if not os.path.exists ( file_ICE_beg) : 
     355    echo ( f'Extracting {file_ICE_beg}' ) 
     356    base_file_ICE_beg = os.path.basename (file_ICE_beg) 
    249357    if not os.path.exists ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearRes}1231_restart_icemod_0000.nc') ): 
    250358        command =  f'cd {TmpDirICE} ; tar xf {file_restart_beg}  ICE_{JobName}_{YearRes}1231_restart_icemod_*.nc' 
    … …  
    253361    echo ('extract ndomain' ) 
    254362    ndomain_opa = get_ndomain ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearRes}1231_restart_icemod') ) 
    255     command = f'cd {TmpDirICE} ; {rebuild} ICE_{JobName}_{YearRes}1231_restart_icemod {ndomain_opa} ; mv {file_ICE_beg} {TmpDir} ' 
     363    command = f'cd {TmpDirICE} ; {rebuild} ICE_{JobName}_{YearRes}1231_restart_icemod {ndomain_opa} ; mv {base_file_ICE_beg} {TmpDir} ' 
    256364    echo ( command ) 
    257365    os.system ( command ) 
    258366    echo ( f'Rebuild done : {file_OCE_beg}') 
    259367     
    260 if not os.path.exists ( os.path.join (TmpDir, file_ICE_end) ) : 
     368if not os.path.exists ( file_ICE_end ) : 
     369    echo ( f'Extracting {file_ICE_end}' ) 
     370    base_file_ICE_end = os.path.basename (file_ICE_end) 
    261371    if not os.path.exists ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearEnd}1231_restart_icemod_0000.nc') ): 
    262372        command =  f'cd {TmpDirICE} ; tar xf {file_restart_end} ICE_{JobName}_{YearEnd}1231_restart_icemod_*.nc' 
    … …  
    265375    echo ('extract ndomain' ) 
    266376    ndomain_opa = get_ndomain ( os.path.join (TmpDirICE, f'ICE_{JobName}_{YearRes}1231_restart_icemod') ) 
    267     command = f'cd {TmpDirICE} ; {rebuild} ICE_{JobName}_{YearEnd}1231_restart_icemod {ndomain_opa} ; mv {file_ICE_end} {TmpDir}' 
     377    command = f'cd {TmpDirICE} ; {rebuild} ICE_{JobName}_{YearEnd}1231_restart_icemod {ndomain_opa} ; mv {base_file_ICE_end} {TmpDir}' 
    268378    echo ( command ) 
    269379    os.system ( command ) 
    … …  
    282392    d_ICE_end = xr.open_dataset ( os.path.join (TmpDir, file_ICE_end), decode_times=False, decode_cf=True, drop_variables=['y', 'x']).squeeze() 
    283393 
     394## 
     395config_out = open (FullIniFile, 'w') 
     396config.write (config_out ) 
     397config_out.close () 
     398 
     399     
     400def kg2Sv  (val, rho=OCE_rho_liq) : 
     401    '''From kg to Sverdrup''' 
     402    return val/dtime_sec*1.0e-6/rho 
     403 
     404def kg2myear (val, rho=OCE_rho_liq) : 
     405    '''From kg to m/year''' 
     406    return val/OCE_aire_tot/rho/NbYear 
     407 
     408def var2prt (var, small=False) : 
     409    if small :  return var , kg2Sv(var)*1000., kg2myear(var)*1000. 
     410    else     :  return var , kg2Sv(var)      , kg2myear(var) 
     411 
     412def prtFlux (Desc, var, Form='F', small=False) : 
     413    if small : 
     414        if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} mSv | {:12.4f} mm/year " 
     415        if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} mSv | {:12.4e} mm/year " 
     416    else : 
     417        if Form in ['f', 'F'] : ff=" {:12.4e} kg | {:12.4f} Sv  | {:12.4f} m/year  " 
     418        if Form in ['e', 'E'] : ff=" {:12.4e} kg | {:12.4e} Sv  | {:12.4e} m/year  " 
     419    echo ( (' {:>15} = ' +ff).format (Desc, *var2prt(var, small) ) ) 
     420    return None 
     421 
    284422# Get mask and surfaces 
    285423sos = d_OCE_his ['sos'][0].squeeze() 
    … …  
    300438def ONE_stock_int (stock) : 
    301439    '''Sum stock''' 
    302     ONE_stock_int = np.sum (  np.sort ( (stock * OCE_aire).to_masked_array().ravel()) ) 
     440    ONE_stock_int = np.sum (  np.sort ( (stock ).to_masked_array().ravel()) ) 
    303441    return ONE_stock_int 
    304442 
    … …  
    306444    '''Integrate flux on oce grid''' 
    307445    OCE_flux_int = np.sum (  np.sort ( (flux * OCE_aire * dtime_per_sec).to_masked_array().ravel()) ) 
     446    return OCE_flux_int 
     447 
     448def ONE_flux_int (flux) : 
     449    '''Integrate flux on oce grid''' 
     450    OCE_flux_int = np.sum (  np.sort ( (flux * dtime_per_sec).to_masked_array().ravel()) ) 
    308451    return OCE_flux_int 
    309452 
    … …  
    343486echo ( 'dOCE vol    = {:12.3e} m^3'.format ( dOCE_vol_wat) ) 
    344487echo ( 'dOCE ssh    = {:12.3e} m  '.format ( dOCE_vol_wat/OCE_aire_tot) ) 
    345 echo ( 'dOCE mass   = {:12.3e} kg '.format ( dOCE_mas_wat) ) 
    346 echo ( 'dOCE mass   = {:12.3e} Sv '.format ( dOCE_mas_wat/dtime_sec*1E-9) ) 
     488prtFlux ( 'dOCE mass ', dOCE_mas_wat, 'e' ) 
    347489 
    348490if NEMO == 3.6 : 
    349491    echo ( 'dOCE e3tn   vol    = {:12.3e} m^3'.format ( dOCE_e3tn_vol) ) 
    350     echo ( 'dOCE e3tn   mass   = {:12.3e} kg '.format ( dOCE_e3tn_mas) ) 
     492    prtFlux ( 'dOCE e3tn   mass', dOCE_e3tn_mas, 'e' ) 
    351493 
    352494## Glace et neige 
    353 if NEMO == 3.6 : 
     495if NEMO == 3.6 :         
    354496    ICE_ice_beg = d_ICE_beg['v_i_htc1']+d_ICE_beg['v_i_htc2']+d_ICE_beg['v_i_htc3']+d_ICE_beg['v_i_htc4']+d_ICE_beg['v_i_htc5'] 
    355497    ICE_ice_end = d_ICE_end['v_i_htc1']+d_ICE_end['v_i_htc2']+d_ICE_end['v_i_htc3']+d_ICE_end['v_i_htc4']+d_ICE_end['v_i_htc5'] 
    … …  
    358500    ICE_sno_end = d_ICE_end['v_s_htc1']+d_ICE_end['v_s_htc2']+d_ICE_end['v_s_htc3']+d_ICE_end['v_s_htc4']+d_ICE_end['v_s_htc5'] 
    359501     
    360     ICE_pnd_beg = 0.0 
    361     ICE_pnd_end = 0.0 
    362      
    363     ICE_fzl_beg = 0.0 
    364     ICE_fzl_end = 0.0 
     502    ICE_pnd_beg = 0.0 ; ICE_pnd_end = 0.0 ; ICE_fzl_beg = 0.0 ; ICE_fzl_end = 0.0 
    365503 
    366504    ICE_mas_wat_beg = ICE_ice_beg*ICE_rho_ice + ICE_sno_beg*ICE_rho_sno 
    … …  
    424562echo ( '\n------------------------------------------------------------') 
    425563echo ( 'Variation du contenu en eau ocean + glace ' ) 
    426 echo ( 'dMass (ocean)   = {:12.6e} kg '.format(dSEA_mas_wat) ) 
    427 echo ( 'dVol  (ocean)   = {:12.3e} Sv '.format(dSEA_mas_wat/dtime_sec*1E-9) ) 
    428 echo ( 'dVol  (ocean)   = {:12.3e} m  '.format(dSEA_mas_wat*1E-3/OCE_aire_tot) ) 
     564prtFlux ( 'dMass (ocean)', dSEA_mas_wat, 'e', True ) 
    429565 
    430566 
    … …  
    504640 
    505641echo ('\n-- Fields:' )  
    506 echo ('OCE+ICE budget {:12.3e} (kg) {:12.3e} (Sv) '.format ( OCE_mas_all       , OCE_mas_all        / dtime_sec*1E-9 )) 
    507 echo ('  EMPMR        {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_empmr     , OCE_mas_empmr      / dtime_sec*1E-9 )) 
    508 echo ('  WFOB         {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfob      , OCE_mas_wfob       / dtime_sec*1E-9 )) 
    509 echo ('  EMP_OCE      {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_emp_oce   , OCE_mas_emp_oce    / dtime_sec*1E-9 )) 
    510 echo ('  EMP_ICE      {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_emp_ice   , OCE_mas_emp_ice    / dtime_sec*1E-9 )) 
    511 echo ('  EMP          {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_emp       , OCE_mas_emp        / dtime_sec*1E-9 )) 
    512 echo ('  ICEBERG      {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_iceberg   , OCE_mas_iceberg    / dtime_sec*1E-9 )) 
    513 echo ('  ICESHELF     {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_iceshelf  , OCE_mas_iceshelf   / dtime_sec*1E-9 )) 
    514 echo ('  CALVING      {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_calving   , OCE_mas_calving    / dtime_sec*1E-9 )) 
    515 echo ('  FRIVER       {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_friver    , OCE_mas_friver     / dtime_sec*1E-9 ))   
    516 echo ('  RUNOFFS      {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_runoffs   , OCE_mas_runoffs    / dtime_sec*1E-9 )) 
    517 echo ('  WFXICE       {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfxice    , OCE_mas_wfxice     / dtime_sec*1E-9 )) 
    518 echo ('  WFXSNW       {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfxsnw    , OCE_mas_wfxsnw     / dtime_sec*1E-9 )) 
    519 echo ('  WFXSUB       {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfxsub    , OCE_mas_wfxsub     / dtime_sec*1E-9 )) 
    520 echo ('  WFXPND       {:12.3e} (kg) {:12.3f} (Sv) '.format ( ICE_mas_wfxpnd    , ICE_mas_wfxpnd     / dtime_sec*1E-9 )) 
    521 echo ('  WFXSNW_SUB   {:12.3e} (kg) {:12.3f} (Sv) '.format ( ICE_mas_wfxsnw_sub, ICE_mas_wfxsnw_sub / dtime_sec*1E-9 )) 
    522 echo ('  WFXSNW_PRE   {:12.3e} (kg) {:12.3f} (Sv) '.format ( ICE_mas_wfxsnw_pre, ICE_mas_wfxsnw_pre / dtime_sec*1E-9 )) 
    523 echo ('  WFXSUB_ERR   {:12.3e} (kg) {:12.3e} (Sv) '.format ( ICE_mas_wfxsub_err, ICE_mas_wfxsub_err / dtime_sec*1E-9 )) 
    524 echo ('  EVAP_OCE     {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_evap_oce  , OCE_mas_evap_oce   / dtime_sec*1E-9 )) 
    525 echo ('  EVAP_ICE     {:12.3e} (kg) {:12.3f} (Sv) '.format ( ICE_mas_evap_ice  , ICE_mas_evap_ice   / dtime_sec*1E-9 )) 
    526 echo ('  SNOW_OCE     {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_snow_oce  , OCE_mas_snow_oce   / dtime_sec*1E-9 )) 
    527 echo ('  SNOW_ICE     {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_snow_ice  , OCE_mas_snow_ice   / dtime_sec*1E-9 )) 
    528 echo ('  RAIN         {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_rain      , OCE_mas_rain       / dtime_sec*1E-9 )) 
    529 echo ('  FWB          {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_fwb       , OCE_mas_fwb        / dtime_sec*1E-9 )) 
    530 echo ('  SSR          {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_ssr       , OCE_mas_ssr        / dtime_sec*1E-9 )) 
    531 echo ('  WFCORR       {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfcorr    , OCE_mas_wfcorr     / dtime_sec*1E-9 )) 
    532 echo ('  BERG_ICB     {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_berg_icb  , OCE_mas_berg_icb   / dtime_sec*1E-9 )) 
    533 echo ('  CALV_ICB     {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_calv_icb  , OCE_mas_calv_icb   / dtime_sec*1E-9 )) 
     642prtFlux ('OCE+ICE budget ', OCE_mas_all       , 'e', True) 
     643prtFlux ('  EMPMR        ', OCE_mas_empmr     , 'e', True) 
     644prtFlux ('  WFOB         ', OCE_mas_wfob      , 'e', True) 
     645prtFlux ('  EMP_OCE      ', OCE_mas_emp_oce   , 'e', True) 
     646prtFlux ('  EMP_ICE      ', OCE_mas_emp_ice   , 'e', True) 
     647prtFlux ('  EMP          ', OCE_mas_emp       , 'e', True) 
     648prtFlux ('  ICEBERG      ', OCE_mas_iceberg   , 'e', True) 
     649prtFlux ('  ICESHELF     ', OCE_mas_iceshelf  , 'e', True) 
     650prtFlux ('  CALVING      ', OCE_mas_calving   , 'e', True) 
     651prtFlux ('  FRIVER       ', OCE_mas_friver    , 'e', True)   
     652prtFlux ('  RUNOFFS      ', OCE_mas_runoffs   , 'e', True) 
     653prtFlux ('  WFXICE       ', OCE_mas_wfxice    , 'e', True) 
     654prtFlux ('  WFXSNW       ', OCE_mas_wfxsnw    , 'e', True) 
     655prtFlux ('  WFXSUB       ', OCE_mas_wfxsub    , 'e', True) 
     656prtFlux ('  WFXPND       ', ICE_mas_wfxpnd    , 'e', True) 
     657prtFlux ('  WFXSNW_SUB   ', ICE_mas_wfxsnw_sub, 'e', True) 
     658prtFlux ('  WFXSNW_PRE   ', ICE_mas_wfxsnw_pre, 'e', True) 
     659prtFlux ('  WFXSUB_ERR   ', ICE_mas_wfxsub_err, 'e', True) 
     660prtFlux ('  EVAP_OCE     ', OCE_mas_evap_oce  , 'e' ) 
     661prtFlux ('  EVAP_ICE     ', ICE_mas_evap_ice  , 'e', True) 
     662prtFlux ('  SNOW_OCE     ', OCE_mas_snow_oce  , 'e', True) 
     663prtFlux ('  SNOW_ICE     ', OCE_mas_snow_ice  , 'e', True) 
     664prtFlux ('  RAIN         ', OCE_mas_rain      , 'e' ) 
     665prtFlux ('  FWB          ', OCE_mas_fwb       , 'e', True) 
     666prtFlux ('  SSR          ', OCE_mas_ssr       , 'e', True) 
     667prtFlux ('  WFCORR       ', OCE_mas_wfcorr    , 'e', True) 
     668prtFlux ('  BERG_ICB     ', OCE_mas_berg_icb  , 'e', True) 
     669prtFlux ('  CALV_ICB     ', OCE_mas_calv_icb  , 'e', True) 
    534670 
    535671echo     ('\n===>  Comparing ===>' )  
Note: See TracChangeset for help on using the changeset viewer.