Changeset 6271


Ignore:
Timestamp:
11/25/22 10:02:09 (23 months ago)
Author:
omamce
Message:

O.M. : Water Bugget

Various improvments
OCE budget is OK
Use .ini file as input

Location:
TOOLS/WATER_BUDGET
Files:
1 added
2 edited

Legend:

Unmodified
Added
Removed

  • TOOLS/WATER_BUDGET/ATM_waterbudget.py

    r6265 r6271  
    1919#  $HeadURL$ 
    2020 
    21  
    22 ## Define Experiment 
    23 if False :  
    24     JobName="TEST-CM72-SIMPLE-ROUTING.12" ; TagName="IPSLCM7" ; SpaceName="DEVT" ; ExperimentName="piControl" ; User="p86caub" ; Project="gencmip6"  
    25     Freq = 'MO' ; YearBegin = 1970 ; YearEnd = 1979 ; PackFrequency = 10 
    26     ATM = 'ICO40' ; Routing = 'SIMPLE' 
    27  
    28 if False : 
    29     JobName="TEST-CM72-SIMPLE-ROUTING.10" ; TagName="IPSLCM7" ; SpaceName="DEVT" ; ExperimentName="piControl" ; User="p86caub" ; Project="gencmip6"  
    30     Freq = 'MO' ; YearBegin = 1860 ; YearEnd = 1869 ; PackFrequency = 10 
    31     ATM = 'ICO40' ; Routing = 'SIMPLE' 
    32  
    33 if False :  
    34     JobName="VALID-CM622-LR.01" ; TagName="IPSLCM6" ; SpaceName="DEVT" ; ExperimentName="piControl" ; User="p86caub" ; Project="gencmip6"  
    35     Freq = 'MO' ; YearBegin = 1890 ; YearEnd = 1899 ; PackFrequency = 10 
    36     ATM = 'LMD144142' ; Routing = 'ORCHIDEE' 
    37      
    38 if True : 
    39     JobName="CM65v420-LR-SKL-pi-05" ; TagName="IPSLCM6" ; SpaceName="DEVT" ; ExperimentName="piControl" ; User="p48ethe" ; Project="gencmip6"  
    40     Freq = 'MO' ;  YearBegin = 1890 ; YearEnd = 1899 ; PackFrequency = 10 
    41     ORCA = 'eORCA1.2'  ; ATM = 'LMD144142' ; Routing = 'ORCHIDEE' ; NEMO=4.2 ; OCE_relax = False ; OCE_icb = False ; Coupled = True 
    42      
    43 Coupled = True 
    44  
    45 import numpy as np 
    46 ##-- Some physical constants 
    47 #-- Earth Radius 
    48 Ra = 40.e6 / (2. * np.pi) 
    49 #-- Gravity 
    50 g = 9.81 
    51 #-- Ice density (kg/m3) in LIM3 
    52 ICE_rho_ice = 917.0 
    53 #-- Snow density (kg/m3) in LIM3 
    54 ICE_rho_sno = 330.0 
    55 #-- Ocean water density (kg/m3) in NEMO 
    56 OCE_rho_liq = 1026. 
    57  
    58 ### 
    59 ICO  = False 
    60 if 'ICO' in ATM : ICO  = True 
    61 LMDZ = False 
    62 if 'LMD' in ATM : LMDZ = True 
    63  
    6421### 
    6522## Import system modules 
    6623import sys, os, shutil, subprocess 
     24import numpy as np 
     25import configparser, re 
     26 
     27## Creates parser 
     28config = configparser.ConfigParser() 
     29config.optionxform = str # To keep capitals 
     30 
     31config['Files']  = {} 
     32config['System'] = {} 
     33 
     34##-- Some physical constants 
     35#-- Earth Radius 
     36Ra = 6366197.7236758135 
     37#-- Gravity 
     38Grav = 9.81 
     39#-- Ice volumic mass (kg/m3) in LIM3 
     40ICE_rho_ice = 917.0 
     41#-- Snow volumic mass (kg/m3) in LIM3 
     42ICE_rho_sno = 330.0 
     43#-- Ocean water volumic mass (kg/m3) in NEMO 
     44OCE_rho_liq = 1026. 
     45#-- Water volumic mass in atmosphere 
     46ATM_rho = 1.0e3 
     47#-- Water volumic mass in surface reservoirs 
     48SRF_rho = 1.0e3 
     49#-- Water volumic mass of rivers 
     50RUN_rho = 1.0e3 
     51 
     52## Read experiment parameters 
     53ATM = None ; ORCA = None ; NEMO = None ; OCE_relax = False ; OCE_icb = False ; Coupled = False ; Routing = None 
     54 
     55# Arguments passed 
     56print ( "Name of Python script:", sys.argv[0] ) 
     57IniFile =  sys.argv[1] 
     58print ("Input file : ", IniFile ) 
     59config.read (IniFile) 
     60 
     61def setBool (chars) : 
     62    '''Convert specific char string in boolean if possible''' 
     63    setBool = chars 
     64    for key in configparser.ConfigParser.BOOLEAN_STATES.keys () : 
     65        if chars.lower() == key : setBool = configparser.ConfigParser.BOOLEAN_STATES[key] 
     66    return setBool 
     67 
     68def setNum (chars) : 
     69    '''Convert specific char string in integer or real if possible''' 
     70    if type (chars) == str : 
     71        realnum = re.compile ("^[-+]?[0-9]*\.?[0-9]+(e[-+]?[0-9]+)?$") 
     72        isReal = realnum.match(chars.strip()) != None 
     73        isInt  = chars.strip().isdigit() 
     74        if isReal : 
     75            if isInt : setNum = int   (chars) 
     76            else     : setNum = float (chars) 
     77        else : setNum = chars 
     78    else : setNum = chars 
     79    return setNum 
     80 
     81print ('[Experiment]') 
     82for 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# ### 
     89ICO  = ( 'ICO' in ATM ) 
     90LMDZ = ( 'LMD' in ATM ) 
     91     
     92### 
     93## Import system modules 
     94import sys, os, shutil, subprocess 
     95import configparser, re 
     96 
     97config = configparser.ConfigParser() 
     98config['Files'] = {} 
    6799 
    68100# Where do we run ? 
    69 TGCC = False ; IDRIS = False 
    70101SysName, NodeName, Release, Version, Machine = os.uname() 
    71 if 'irene'   in NodeName : TGCC  = True 
    72 if 'jeanzay' in NodeName : IDRIS = True 
     102TGCC  = ( 'irene'   in NodeName ) 
     103IDRIS = ( 'jeanzay' in NodeName ) 
    73104 
    74105## Set site specific libIGCM directories, and other specific stuff 
    … …  
    90121         
    91122    ## Creates output directory 
    92     TmpDir = os.path.join ( subprocess.getoutput ( 'ccc_home --cccscratch' ), f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) 
     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}' ) 
    93125 
    94126if IDRIS : 
    95     raise Exception("Pour IDRIS : repertoires et chemins a definir")  
     127    raise Exception ("Pour IDRIS : repertoires et chemins a definir")  
    96128 
    97129## Import specific module 
    … …  
    99131## Now import needed scientific modules 
    100132import xarray as xr 
    101      
     133 
     134config['Files'][TmpDir] = TmpDir 
     135 
    102136# Output file 
    103137FileOut = f'ATM_waterbudget_{JobName}_{YearBegin}_{YearEnd}.out' 
    … …  
    105139 
    106140# Function to print to stdout *and* output file 
    107 def echo (string) : 
    108     print ( string  ) 
    109     f_out.write ( string + '\n' ) 
     141def echo (string, end='\n') : 
     142    print ( string, end=end  ) 
     143    sys.stdout.flush () 
     144    f_out.write ( string + end ) 
    110145    f_out.flush () 
    111146    return None 
    … …  
    135170 
    136171#-- Model output directories 
    137 if Freq == "MO" : FreqDir =  os.path.join ('Output' , 'MO' ) 
    138 if Freq == "SE" : FreqDir =  os.path.join ('Analyse', 'SE' ) 
     172if Freq == "MO" : FreqDir = os.path.join ('Output' , 'MO' ) 
     173if Freq == "SE" : FreqDir = os.path.join ('Analyse', 'SE' ) 
    139174dir_ATM_his = os.path.join ( R_SAVE, "ATM", FreqDir ) 
    140175dir_SRF_his = os.path.join ( R_SAVE, "SRF", FreqDir ) 
    … …  
    145180 
    146181#-- Files Names 
    147 if Freq == 'MO' : File = f'{JobName}_{YearBegin}0101_{YearEnd}1231_1M' 
    148 if Freq == 'SE' : File = f'{JobName}_SE_{YearBegin}0101_{YearEnd}1231_1M' 
     182if Freq == 'MO' : FileCommon = f'{JobName}_{YearBegin}0101_{YearEnd}1231_1M' 
     183if Freq == 'SE' : FileCommon = f'{JobName}_SE_{YearBegin}0101_{YearEnd}1231_1M' 
    149184 
    150185echo ('\nOpen history files' ) 
    151 file_ATM_his = os.path.join (  dir_ATM_his, f'{File}_histmth.nc' ) 
    152 file_SRF_his = os.path.join (  dir_SRF_his, f'{File}_sechiba_history.nc' ) 
     186file_ATM_his = os.path.join (  dir_ATM_his, f'{FileCommon}_histmth.nc' ) 
     187file_SRF_his = os.path.join (  dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 
     188if Routing == 'ORCHIDEE'  : 
     189    file_RUN_his = os.path.join (  dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 
     190if Routing == 'SIMPLE' : 
     191    file_RUN_his = os.path.join (  dir_SRF_his, f'{FileCommon}_sechiba_history.nc' ) 
    153192 
    154193d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 
    155194d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 
     195if Routing == 'ORCHIDEE' : 
     196    d_RUN_his = d_SRF_his 
     197if Routing == 'SIMPLE' :  
     198    d_RUN_his = xr.open_dataset ( file_RUN_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 
    156199 
    157200echo ( file_ATM_his ) 
    158201echo ( file_SRF_his ) 
     202echo ( file_RUN_his ) 
     203 
     204 
     205config['Files']['file_ATM_his'] = file_ATM_his 
     206config['Files']['file_SRF_his'] = file_SRF_his 
     207config['Files']['file_RUN_his'] = file_SRF_his 
    159208 
    160209## Compute run length 
    … …  
    170219 
    171220#-- Open restart files 
    172 YearRes = YearBegin - 1        # Year of the restart of beginning of simulation 
     221YearRes = YearBegin - 1              # Year of the restart of beginning of simulation 
    173222YearPre = YearBegin - PackFrequency  # Year to find the tarfile of the restart of beginning of simulation 
    174223 
    … …  
    219268        os.system ( command ) 
    220269 
     270config['Files']['file_ATM_beg'] = file_ATM_beg 
     271config['Files']['file_ATM_end'] = file_ATM_end 
     272config['Files']['file_SRF_beg'] = file_SRF_beg 
     273config['Files']['file_SRF_end'] = file_SRF_end 
     274if Routing == 'SIMPLE' :  
     275    config['Files']['file_RUN_beg'] = file_RUN_beg 
     276    config['Files']['file_RUN_end'] = file_RUN_end 
     277config['Files']['file_DYN_beg'] = file_DYN_beg 
     278config['Files']['file_DYN_end'] = file_DYN_end 
     279         
    221280echo ('\nOpening ATM SRF and ICO restart files') 
    222281d_ATM_beg = xr.open_dataset ( os.path.join (TmpDir, file_ATM_beg), decode_times=False, decode_cf=True).squeeze() 
    … …  
    228287 
    229288for var in d_SRF_beg.variables : 
    230     #d_SRF_beg[var].attrs['_FillValue'] = 1.e20 
    231     #d_SRF_end[var].attrs['_FillValue'] = 1.e20 
    232289    d_SRF_beg[var] = d_SRF_beg[var].where (  d_SRF_beg[var]<1.e20, 0.) 
    233290    d_SRF_end[var] = d_SRF_end[var].where (  d_SRF_end[var]<1.e20, 0.) 
    … …  
    247304    echo ( file_RUN_end ) 
    248305 
     306def ATM_stock_int (stock) : 
     307    '''Integrate stock on atmosphere grid''' 
     308    ATM_stock_int  = np.sum (  np.sort ( (stock * DYN_aire).to_masked_array().ravel()) ) 
     309    return ATM_stock_int 
     310 
     311def ATM_flux_int (flux) : 
     312    '''Integrate flux on atmosphere grid''' 
     313    ATM_stock_int  = np.sum (  np.sort ( (flux * dtime_per_sec * ATM_aire).to_masked_array().ravel()) ) 
     314    return ATM_stock_int 
     315 
     316def ONE_stock_int (stock) : 
     317    '''Sum stock ''' 
     318    ONE_stock_int  = np.sum (  np.sort ( (stock).to_masked_array().ravel()) ) 
     319    return ONE_stock_int 
     320 
     321def ONE_flux_int (flux) : 
     322    '''Sum flux ''' 
     323    ONE_flux_int  = np.sum (  np.sort ( (flux * dtime_per_sec ).to_masked_array().ravel()) ) 
     324    return ONE_flux_int 
     325     
    249326# ATM grid with cell surfaces 
    250 if ICO :  
    251     ATM_aire = d_ATM_his ['aire'][0] 
    252     ATM_fsea = d_ATM_his ['fract_ter'][Ø] + d_ATM_his ['fract_sic'][Ø] 
     327if ICO : 
     328    jpja, jpia = d_ATM_his['aire'][0].shape    
     329    file_ATM_aire = os.path.join ( R_IN, 'ATM', 'GRID', f'aire_{ATM}_to_{jpia}x{jpja}.nc' ) 
     330    config['Files']['file_ATM_aire'] = file_ATM_aire 
     331    echo ('Aire sur grille reguliere :', file_ATM_aire) 
     332    d_ATM_aire = xr.open_dataset ( file_ATM_aire, decode_times=False ).squeeze() 
     333    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] ) 
     336    SRF_aire = lmdz.geo2point ( d_SRF_his ['Areas'] * d_SRF_his ['Contfrac'] ) 
     337     
    253338if LMDZ : 
    254     ATM_aire = lmdz.geo2point ( d_ATM_his ['aire'][0] ) 
    255     ATM_fsea = lmdz.geo2point ( d_ATM_his ['fract_ter'][0] + d_ATM_his ['fract_sic'][0] ) 
    256     ATM_aire[0]  = np.sum ( d_ATM_his ['aire'][0, 0] ) 
    257     ATM_aire[-1] = np.sum ( d_ATM_his ['aire'][0,-1] ) 
     339    ATM_aire = lmdz.geo2point ( d_ATM_his ['aire'][0], cumulPoles=True ) 
     340    ATM_fsea = lmdz.geo2point ( d_ATM_his ['fract_oce'][0] + d_ATM_his ['fract_sic'][0] ) 
     341    ATM_flnd = lmdz.geo2point ( d_ATM_his ['fract_ter'][0] ) 
     342    SRF_aire = lmdz.geo2point ( d_SRF_his['Areas'] * d_SRF_his['Contfrac'] ) 
     343 
     344SRF_aire = SRF_aire.where ( SRF_aire < 1E15, 0.) 
    258345 
    259346if ICO : 
     347    # Area on icosahedron grid 
    260348    file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ATM+'_grid.nc' ) 
    261349    d_DYN_aire = xr.open_dataset ( file_DYN_aire, decode_times=False).squeeze() 
    262     d_DYN_aire = d_DYN_aire.rename( {'cell':'cell_mesh'}) 
     350    d_DYN_aire = d_DYN_aire.rename ( {'cell':'cell_mesh'} ) 
    263351    DYN_aire   = d_DYN_aire['aire'] 
     352 
     353    DYN_fsea = d_DYN_aire ['fract_oce'] + d_DYN_aire ['fract_sic'] 
     354    DYN_flnd = 1.0 - DYN_fsea 
     355     
    264356if LMDZ : 
    265357    DYN_aire = ATM_aire 
     358    DYN_fsea = ATM_fsea 
     359    DYN_flnd = ATM_flnd 
    266360 
    267361#if LMDZ : 
    268362#    d_ATM_beg = d_ATM_beg.assign ( coords={'lon':d_ATM_beg.lon*180./np.pi} ) 
    269363 
    270 ATM_aire_tot = np.sum (ATM_aire).values.item() 
    271 ATM_aire_sea_tot = np.sum (ATM_aire*ATM_fsea).values.item() 
     364ATM_aire_sea     = ATM_aire * ATM_fsea 
     365ATM_aire_sea_tot = ONE_stock_int ( ATM_aire_sea ) 
     366 
     367ATM_aire_tot     = ONE_stock_int (ATM_aire) 
     368ATM_aire_sea_tot = ONE_stock_int (ATM_aire*ATM_fsea) 
     369 
     370 
     371echo ( 'Aire atmosphere/4pi R^2 : {:12.5f}'.format(ATM_aire_tot/(Ra*Ra*4*np.pi) ) ) 
     372 
     373if (  np.abs (ATM_aire_tot/(Ra*Ra*4*np.pi) - 1.0) > 0.01 ) : 
     374    raise Exception ('Erreur surface interpolee sur grille reguliere') 
    272375 
    273376echo ( '\n------------------------------------------------------------------------------------' ) 
    274 echo ( '-- LMDZ changes in stores (for the records)' ) 
     377echo ( '-- ATM changes in stores ' ) 
    275378 
    276379#-- Change in precipitable water from the atmosphere daily and monthly files 
    … …  
    280383ATM_Ahyb = d_ATM_his['Ahyb'].squeeze() 
    281384ATM_Bhyb = d_ATM_his['Bhyb'].squeeze() 
    282 klevp1 = ATM_Ahyb.shape[0] 
    283385 
    284386# Surface pressure 
    285387if ICO :  
    286     ATM_ps_beg = d_DYN_beg['ps'] 
    287     ATM_ps_end = d_DYN_end['ps'] 
     388    DYN_ps_beg = d_DYN_beg['ps'] 
     389    DYN_ps_end = d_DYN_end['ps'] 
    288390     
    289391if LMDZ :  
    290     ATM_ps_beg = lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)) ) 
    291     ATM_ps_end = lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)) ) 
     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)) ) 
    292394     
    293395# 3D Pressure 
    294 ATM_p_beg = ATM_Ahyb + ATM_Bhyb * ATM_ps_beg 
    295 ATM_p_end = ATM_Ahyb + ATM_Bhyb * ATM_ps_end 
     396DYN_p_beg = ATM_Ahyb + ATM_Bhyb * DYN_ps_beg 
     397DYN_p_end = ATM_Ahyb + ATM_Bhyb * DYN_ps_end 
     398 
     399if ICO  : klevp1, cell_mesh        = DYN_p_beg.shape 
     400if LMDZ : klevp1, points_physiques = DYN_p_beg.shape 
     401klev = klevp1 - 1 
    296402 
    297403# Layer thickness 
    298 ATM_sigma_beg = ATM_p_beg[0:-1]*0. 
    299 ATM_sigma_end = ATM_p_end[0:-1]*0.  
     404if 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) ) ) 
     407if 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) ) ) 
    300410 
    301411for k in np.arange (klevp1-1) : 
    302     ATM_sigma_beg[k,:] = (ATM_p_beg[k,:] - ATM_p_beg[k+1,:]) / g 
    303     ATM_sigma_end[k,:] = (ATM_p_end[k,:] - ATM_p_end[k+1,:]) / g 
    304  
    305 ATM_sigma_beg = ATM_sigma_beg.rename ( {'klevp1':'sigs'} ) 
    306 ATM_sigma_end = ATM_sigma_end.rename ( {'klevp1':'sigs'} ) 
     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,:] 
    307414 
    308415##-- Vertical and horizontal integral, and sum of liquid, solid and vapor water phases 
    309416if LMDZ : 
    310417    try :  
    311         ATM_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2Ov'] + d_DYN_beg['H2Ol'] + d_DYN_beg['H2Oi']).isel(rlonv=slice(0,-1) ) ) 
    312         ATM_wat_end = lmdz.geo3point ( (d_DYN_end['H2Ov'] + d_DYN_end['H2Ol'] + d_DYN_end['H2Oi']).isel(rlonv=slice(0,-1) ) ) 
     418        DYN_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2Ov']  + d_DYN_beg['H2Ol']  + d_DYN_beg['H2Oi'] ).isel(rlonv=slice(0,-1) ) ) 
     419        DYN_wat_end = lmdz.geo3point ( (d_DYN_end['H2Ov']  + d_DYN_end['H2Ol']  + d_DYN_end['H2Oi'] ).isel(rlonv=slice(0,-1) ) ) 
    313420    except : 
    314         ATM_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s']).isel(rlonv=slice(0,-1) ) ) 
    315         ATM_wat_end = lmdz.geo3point ( (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s']).isel(rlonv=slice(0,-1) ) ) 
     421        DYN_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s']).isel(rlonv=slice(0,-1) ) ) 
     422        DYN_wat_end = lmdz.geo3point ( (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s']).isel(rlonv=slice(0,-1) ) ) 
    316423if ICO : 
    317     ATM_wat_beg = ( d_DYN_beg['q'][0:4] ).sum (dim = ['nq',] ).rename ( {'lev':'sigs'} ) 
    318     ATM_wat_end = ( d_DYN_end['q'][0:4] ).sum (dim = ['nq',] ).rename ( {'lev':'sigs'} ) 
    319      
    320 ATM_mas_wat_beg = np.sum (ATM_sigma_beg * ATM_wat_beg * ATM_aire).values.item() 
    321 ATM_mas_wat_end = np.sum (ATM_sigma_end * ATM_wat_end * ATM_aire).values.item() 
    322  
    323 dATM_mas_wat = ATM_mas_wat_end - ATM_mas_wat_beg 
    324  
    325 echo ( 'Variation du contenu en eau atmosphere ' ) 
    326 echo ( 'ATM_mass_beg = {:12.6e} kg - ATM_mass_end = {:12.6e} kg'.format (ATM_mas_wat_beg, ATM_mas_wat_end) ) 
    327 echo ( 'dMass(atm)   = {:12.3e} kg '.format (dATM_mas_wat) ) 
    328 echo ( 'dMass(atm)   = {:12.3e} Sv '.format (dATM_mas_wat/dtime_sec*1.E-9) ) 
    329 echo ( 'dMass(atm)   = {:12.3e}m   '.format (dATM_mas_wat/ATM_aire_sea_tot*1E-3) ) 
     424    try : 
     425        DYN_wat_beg = (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s']).rename ( {'lev':'sigs'} ) 
     426        DYN_wat_end = (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s']).rename ( {'lev':'sigs'} ) 
     427    except :  
     428        DYN_wat_beg = (d_DYN_beg['q'].isel(nq=0) + d_DYN_beg['q'].isel(nq=1) + d_DYN_beg['q'].isel(nq=2) ).rename ( {'lev':'sigs'} ) 
     429        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'} ) 
     430     
     431# Integral 
     432DYN_mas_wat_beg = ATM_stock_int (DYN_dsigma_beg * DYN_wat_beg) / Grav 
     433DYN_mas_wat_end = ATM_stock_int (DYN_dsigma_end * DYN_wat_end) / Grav 
     434 
     435dDYN_mas_wat = DYN_mas_wat_end - DYN_mas_wat_beg 
     436 
     437echo ( '\nVariation du contenu en eau atmosphere (dynamique) ' ) 
     438echo ( 'DYN_mas_beg = {:12.6e} kg | DYN_mas_end = {:12.6e} kg'.format (DYN_mas_wat_beg, DYN_mas_wat_end) ) 
     439echo ( 'dMass (atm)   = {:12.3e} kg '.format (dDYN_mas_wat) ) 
     440echo ( 'dMass (atm)   = {:12.3e} Sv '.format (dDYN_mas_wat/dtime_sec*1.e-6/ATM_rho) ) 
     441echo ( 'dMass (atm)   = {:12.3e} m  '.format (dDYN_mas_wat/ATM_aire_sea_tot/ATM_rho) ) 
    330442 
    331443ATM_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'] 
    332444ATM_sno_end = d_ATM_end['SNOW01']*d_ATM_end['FTER']+d_ATM_end['SNOW02']*d_ATM_end['FLIC']+d_ATM_end['SNOW03']*d_ATM_end['FOCE']+d_ATM_end['SNOW04']*d_ATM_end['FSIC'] 
     445 
     446ATM_qs_beg  = d_ATM_beg['QS01']*d_ATM_beg['FTER']+d_ATM_beg['QS02']*d_ATM_beg['FLIC']+d_ATM_beg['QS03']*d_ATM_beg['FOCE']+d_ATM_beg['QS04']*d_ATM_beg['FSIC'] 
     447ATM_qs_end  = d_ATM_end['QS01']*d_ATM_end['FTER']+d_ATM_end['QS02']*d_ATM_end['FLIC']+d_ATM_end['QS03']*d_ATM_end['FOCE']+d_ATM_end['QS04']*d_ATM_end['FSIC'] 
     448 
     449ATM_qsol_beg = d_ATM_beg['QSOL'] 
     450ATM_qsol_end = d_ATM_end['QSOL'] 
     451 
     452ATM_qs01_beg  = d_ATM_beg['QS01'] * d_ATM_beg['FTER'] 
     453ATM_qs01_end  = d_ATM_end['QS01'] * d_ATM_end['FTER'] 
     454ATM_qs02_beg  = d_ATM_beg['QS02'] * d_ATM_beg['FLIC'] 
     455ATM_qs02_end  = d_ATM_end['QS02'] * d_ATM_end['FLIC'] 
     456ATM_qs03_beg  = d_ATM_beg['QS03'] * d_ATM_beg['FOCE'] 
     457ATM_qs03_end  = d_ATM_end['QS03'] * d_ATM_end['FOCE'] 
     458ATM_qs04_beg  = d_ATM_beg['QS04'] * d_ATM_beg['FSIC'] 
     459ATM_qs04_end  = d_ATM_end['QS04'] * d_ATM_end['FSIC']   
     460 
    333461if ICO : 
    334    ATM_sno_beg = ATM_sno_beg.rename ( {'points_physiques':'cell_mesh'} ) 
    335    ATM_sno_end = ATM_sno_end.rename ( {'points_physiques':'cell_mesh'} ) 
    336  
    337 ATM_mas_sno_beg = np.sum ( ATM_sno_beg * DYN_aire ).values.item() 
    338 ATM_mas_sno_end = np.sum ( ATM_sno_end * DYN_aire ).values.item() 
    339  
    340 dATM_mas_sno = ATM_mas_sno_end - ATM_mas_sno_beg 
    341  
    342 echo ( 'Variation du contenu en neige atmosphere ' ) 
    343 echo ( 'ATM_mas_sno_beg  = {:12.6e} kg  - ATM_mas_sno_end  = {:12.6e} kg'.format (ATM_mas_sno_beg, ATM_mas_sno_end) ) 
    344 echo ( 'dMass(neige atm) = {:12.3e} kg '.format (dATM_mas_sno) ) 
    345 echo ( 'dMass(neige atm) = {:12.3e} Sv '.format (dATM_mas_sno/dtime_sec*1E-6/ICE_rho_ice) ) 
    346 echo ( 'dMass(neige atm) = {:12.3e} m  '.format (dATM_mas_sno/ATM_aire_sea_tot*1E-3) ) 
     462   ATM_sno_beg   = ATM_sno_beg .rename ( {'points_physiques':'cell_mesh'} ) 
     463   ATM_sno_end   = ATM_sno_end .rename ( {'points_physiques':'cell_mesh'} ) 
     464   ATM_qs_beg    = ATM_qs_beg  .rename ( {'points_physiques':'cell_mesh'} ) 
     465   ATM_qs_end    = ATM_qs_end  .rename ( {'points_physiques':'cell_mesh'} ) 
     466   ATM_qsol_beg  = ATM_qsol_beg.rename ( {'points_physiques':'cell_mesh'} ) 
     467   ATM_qsol_end  = ATM_qsol_end.rename ( {'points_physiques':'cell_mesh'} ) 
     468   ATM_qs01_beg  = ATM_qs01_beg.rename ( {'points_physiques':'cell_mesh'} ) 
     469   ATM_qs01_end  = ATM_qs01_end.rename ( {'points_physiques':'cell_mesh'} ) 
     470   ATM_qs02_beg  = ATM_qs02_beg.rename ( {'points_physiques':'cell_mesh'} ) 
     471   ATM_qs02_end  = ATM_qs02_end.rename ( {'points_physiques':'cell_mesh'} ) 
     472   ATM_qs03_beg  = ATM_qs03_beg.rename ( {'points_physiques':'cell_mesh'} ) 
     473   ATM_qs03_end  = ATM_qs03_end.rename ( {'points_physiques':'cell_mesh'} ) 
     474   ATM_qs04_beg  = ATM_qs04_beg.rename ( {'points_physiques':'cell_mesh'} ) 
     475   ATM_qs04_end  = ATM_qs04_end.rename ( {'points_physiques':'cell_mesh'} )  
     476 
     477ATM_mas_sno_beg   = ATM_stock_int ( ATM_sno_beg  ) 
     478ATM_mas_sno_end   = ATM_stock_int ( ATM_sno_end  ) 
     479ATM_mas_qs_beg    = ATM_stock_int ( ATM_qs_beg   ) 
     480ATM_mas_qs_end    = ATM_stock_int ( ATM_qs_end   ) 
     481ATM_mas_qsol_beg  = ATM_stock_int ( ATM_qsol_beg ) 
     482ATM_mas_qsol_end  = ATM_stock_int ( ATM_qsol_end ) 
     483ATM_mas_qs01_beg  = ATM_stock_int ( ATM_qs01_beg ) 
     484ATM_mas_qs01_end  = ATM_stock_int ( ATM_qs01_end ) 
     485ATM_mas_qs02_beg  = ATM_stock_int ( ATM_qs02_beg ) 
     486ATM_mas_qs02_end  = ATM_stock_int ( ATM_qs02_end ) 
     487ATM_mas_qs03_beg  = ATM_stock_int ( ATM_qs03_beg ) 
     488ATM_mas_qs03_end  = ATM_stock_int ( ATM_qs03_end ) 
     489ATM_mas_qs04_beg  = ATM_stock_int ( ATM_qs04_beg ) 
     490ATM_mas_qs04_end  = ATM_stock_int ( ATM_qs04_end ) 
     491 
     492dATM_mas_sno  = ATM_mas_sno_end  - ATM_mas_sno_beg 
     493dATM_mas_qs   = ATM_mas_qs_end   - ATM_mas_qs_beg 
     494dATM_mas_qsol = ATM_mas_qsol_end - ATM_mas_qsol_beg 
     495 
     496dATM_mas_qs01 = ATM_mas_qs01_end - ATM_mas_qs01_beg 
     497dATM_mas_qs02 = ATM_mas_qs02_end - ATM_mas_qs02_beg 
     498dATM_mas_qs03 = ATM_mas_qs03_end - ATM_mas_qs03_beg 
     499dATM_mas_qs04 = ATM_mas_qs04_end - ATM_mas_qs04_beg 
     500 
     501echo ( '\nVariation du contenu en neige atmosphere (calottes)' ) 
     502echo ( 'ATM_mas_sno_beg  = {:12.6e} kg | ATM_mas_sno_end  = {:12.6e} kg'.format (ATM_mas_sno_beg, ATM_mas_sno_end) ) 
     503echo ( 'dMass (neige atm) = {:12.3e} kg '.format (dATM_mas_sno ) ) 
     504echo ( 'dMass (neige atm) = {:12.3e} Sv '.format (dATM_mas_sno/dtime_sec*1e-6/ICE_rho_ice) ) 
     505echo ( 'dMass (neige atm) = {:12.3e} m  '.format (dATM_mas_sno/ATM_aire_sea_tot/ATM_rho) ) 
     506 
     507echo ( '\nVariation du contenu humidite du sol' ) 
     508echo ( 'ATM_mas_qs_beg  = {:12.6e} kg | ATM_mas_qs_end  = {:12.6e} kg'.format (ATM_mas_qs_beg, ATM_mas_qs_end) ) 
     509echo ( 'dMass (neige atm) = {:12.3e} kg '.format (dATM_mas_qs ) ) 
     510echo ( 'dMass (neige atm) = {:12.3e} Sv '.format (dATM_mas_qs/dtime_sec*1e-6/ATM_rho) ) 
     511echo ( 'dMass (neige atm) = {:12.3e} m  '.format (dATM_mas_qs/ATM_aire_sea_tot/ATM_rho) ) 
    347512 
    348513echo ( '\nVariation du contenu en eau+neige atmosphere ' ) 
    349 echo ( 'dMass(eau + neige atm) = {:12.3e} kg '.format (  dATM_mas_wat + dATM_mas_sno) ) 
    350 echo ( 'dMass(eau + neige atm) = {:12.3e} Sv '.format ( (dATM_mas_wat + dATM_mas_sno)/dtime_sec*1E-9) ) 
    351 echo ( 'dMass(eau + neige atm) = {:12.3e} m  '.format ( (dATM_mas_wat + dATM_mas_sno)/ATM_aire_sea_tot*1E-3) ) 
     514echo ( 'dMass (eau + neige atm) = {:12.3e} kg '.format (  dDYN_mas_wat + dATM_mas_sno) ) 
     515echo ( 'dMass (eau + neige atm) = {:12.3e} Sv '.format ( (dDYN_mas_wat + dATM_mas_sno)/dtime_sec*1E-6/ATM_rho) ) 
     516echo ( 'dMass (eau + neige atm) = {:12.3e} m  '.format ( (dDYN_mas_wat + dATM_mas_sno)/ATM_aire_sea_tot/ATM_rho) ) 
    352517 
    353518echo ( '\n------------------------------------------------------------------------------------' ) 
    354 echo ( '-- SRF changes in routing reservoirs' ) 
     519echo ( '-- 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 
    355534 
    356535if Routing == 'SIMPLE' : 
    357     RUN_mas_wat_beg = np.sum ( d_RUN_beg ['fast_reservoir'] +  d_RUN_beg ['slow_reservoir'] +  d_RUN_beg ['stream_reservoir']).values.item() 
    358     RUN_mas_wat_end = np.sum ( d_RUN_end ['fast_reservoir'] +  d_RUN_end ['slow_reservoir'] +  d_RUN_end ['stream_reservoir']).values.item() 
     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'] ) 
    359538 
    360539if Routing == 'ORCHIDEE' :  
    361     RUN_mas_wat_beg = np.sum ( d_SRF_beg['fastres']  + d_SRF_beg['slowres'] + d_SRF_beg['streamres'] \ 
    362                              + d_SRF_beg['floodres'] + d_SRF_beg['lakeres'] + d_SRF_beg['pondres'] ) .values.item() 
    363     RUN_mas_wat_end = np.sum ( d_SRF_end['fastres']  + d_SRF_end['slowres'] + d_SRF_end['streamres'] \ 
    364                              + d_SRF_end['floodres'] + d_SRF_end['lakeres'] + d_SRF_end['pondres'] ) .values.item() 
     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']  ) 
    365544 
    366545dRUN_mas_wat = RUN_mas_wat_end - RUN_mas_wat_beg 
    367546     
    368547echo ( '\nWater content in routing ' ) 
    369 echo ( 'RUN_mas_wat_beg = {:12.6e} kg '.format (RUN_mas_wat_beg) ) 
    370 echo ( 'RUN_mas_wat_end = {:12.6e} kg '.format (RUN_mas_wat_end) ) 
    371 echo ( 'dMass(routing)  = {:12.3e} kg '.format(dRUN_mas_wat) ) 
    372 echo ( 'dMass(routing)  = {:12.3e} Sv '.format(dRUN_mas_wat/dtime_sec*1E-9) ) 
    373 echo ( 'dMass(routing)  = {:12.3e} m  '.format(dRUN_mas_wat/ATM_aire_sea_tot*1E-3) ) 
    374      
    375 SRF_mas_wat_beg = d_SRF_beg['tot_watveg_beg']+d_SRF_beg['tot_watsoil_beg'] + d_SRF_beg['snow_beg'] 
    376 SRF_mas_wat_end = d_SRF_end['tot_watveg_beg']+d_SRF_end['tot_watsoil_beg'] + d_SRF_end['snow_beg'] 
    377  
    378 if LMDZ :  
    379     SRF_mas_wat_beg = SRF_mas_wat_beg.rename ( {'y':'points_phyiques'} ) 
    380     SRF_mas_wat_end = SRF_mas_wat_end.rename ( {'y':'points_phyiques'} ) 
    381 if ICO :  
    382     SRF_mas_wat_beg = SRF_mas_wat_beg.rename ( {'y':'cell_mesh'} ) 
    383     SRF_mas_wat_end = SRF_mas_wat_end.rename ( {'y':'cell_mesh'} )   
    384  
    385 SRF_mas_wat_beg = np.sum ( SRF_mas_wat_beg * SRF_aire ).values.item() 
    386 SRF_mas_wat_end = np.sum ( SRF_mas_wat_end * SRF_aire ).values.item() 
    387  
     548echo ( 'RUN_mas_wat_beg = {:12.6e} kg | RUN_mas_wat_end = {:12.6e} kg '.format (RUN_mas_wat_end, RUN_mas_wat_end) ) 
     549echo ( 'dMass (routing)  = {:12.3e} kg '.format(dRUN_mas_wat) ) 
     550echo ( 'dMass (routing)  = {:12.3e} Sv '.format(dRUN_mas_wat/dtime_sec*1E-9) ) 
     551echo ( 'dMass (routing)  = {:12.3e} m  '.format(dRUN_mas_wat/ATM_aire_sea_tot*1E-3) ) 
     552 
     553print ('Reading SRF restart') 
     554tot_watveg_beg  = d_SRF_beg['tot_watveg_beg']  ; tot_watveg_beg  = tot_watveg_beg .where (tot_watveg_beg < 1E10, 0.) 
     555tot_watsoil_beg = d_SRF_beg['tot_watsoil_beg'] ; tot_watsoil_beg = tot_watsoil_beg.where (tot_watsoil_beg< 1E10, 0.) 
     556snow_beg        = d_SRF_beg['snow_beg']        ; snow_beg        = snow_beg       .where (snow_beg       < 1E10, 0.) 
     557 
     558tot_watveg_end  = d_SRF_end['tot_watveg_beg']  ; tot_watveg_end  = tot_watveg_end .where (tot_watveg_end < 1E10, 0.) 
     559tot_watsoil_end = d_SRF_end['tot_watsoil_beg'] ; tot_watsoil_end = tot_watsoil_end.where (tot_watsoil_end< 1E10, 0.) 
     560snow_end        = d_SRF_end['snow_beg']        ; snow_end        = snow_end       .where (snow_end       < 1E10, 0.) 
     561 
     562if 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) 
     569 
     570# Stock dSoilHum dInterce dSWE dStream dFastR dSlowR dLake dPond dFlood 
     571 
     572SRF_wat_beg = tot_watveg_beg + tot_watsoil_beg + snow_beg 
     573SRF_wat_end = tot_watveg_end + tot_watsoil_end + snow_end 
     574     
     575if 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 
     585print ('Computing integrals') 
     586 
     587print ( ' 1/6', end='' ) ; sys.stdout.flush () 
     588SRF_mas_watveg_beg   = ATM_stock_int ( tot_watveg_beg  ) 
     589print ( ' 2/6', end='' ) ; sys.stdout.flush () 
     590SRF_mas_watsoil_beg  = ATM_stock_int ( tot_watsoil_beg ) 
     591print ( ' 3/6', end='' ) ; sys.stdout.flush () 
     592SRF_mas_snow_beg     = ATM_stock_int ( snow_beg        ) 
     593print ( ' 4/6', end='' ) ; sys.stdout.flush () 
     594SRF_mas_watveg_end   = ATM_stock_int ( tot_watveg_end  ) 
     595print ( ' 5/6', end='' ) ; sys.stdout.flush () 
     596SRF_mas_watsoil_end  = ATM_stock_int ( tot_watsoil_end ) 
     597print ( ' 6/6', end='' ) ; sys.stdout.flush () 
     598SRF_mas_snow_end     = ATM_stock_int ( snow_end        ) 
     599print (' -- ') ; sys.stdout.flush () 
     600 
     601dSRF_mas_watveg   = SRF_mas_watveg_end   - SRF_mas_watveg_beg 
     602dSRF_mas_watsoil  = SRF_mas_watsoil_end  - SRF_mas_watsoil_beg 
     603dSRF_mas_snow     = SRF_mas_snow_end     - SRF_mas_snow_beg 
     604 
     605echo ('\nLes differents reservoirs') 
     606echo ( 'SRF_mas_watveg_beg   = {:12.6e} kg | SRF_mas_watveg_end   = {:12.6e} kg '.format (SRF_mas_watveg_beg  , SRF_mas_watveg_end   ) ) 
     607echo ( 'SRF_mas_watsoil_beg  = {:12.6e} kg | SRF_mas_watsoil_end  = {:12.6e} kg '.format (SRF_mas_watsoil_beg , SRF_mas_watsoil_end  ) ) 
     608echo ( 'SRF_mas_snow_beg     = {:12.6e} kg | SRF_mas_snow_end     = {:12.6e} kg '.format (SRF_mas_snow_beg    , SRF_mas_snow_end     ) ) 
     609 
     610echo ( '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) ) 
     611echo ( '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) ) 
     612echo ( '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) ) 
     613 
     614SRF_mas_wat_beg = SRF_mas_watveg_beg + SRF_mas_watsoil_beg + SRF_mas_snow_beg  
     615SRF_mas_wat_end = SRF_mas_watveg_end + SRF_mas_watsoil_end + SRF_mas_snow_end  
    388616dSRF_mas_wat = SRF_mas_wat_end - SRF_mas_wat_beg 
    389617 
    390 echo ( '\nWater content in and surface ' ) 
    391 echo ( 'SRF_mas_wat_beg  = {:12.6e} kg - SRF_mas_wat_end  = {:12.6e} kg '.format (SRF_mas_wat_beg, SRF_mas_wat_end) ) 
    392 echo ( 'dMass(water srf) = {:12.3e} kg '.format (dSRF_mas_wat) ) 
    393 echo ( 'dMass(water srf) = {:12.3e} Sv '.format (dSRF_mas_wat/dtime_sec*1E-9) ) 
    394 echo ( 'dMass(water srf) = {:12.3e} m  '.format (dSRF_mas_wat/ATM_aire_sea_tot*1E-3) ) 
     618echo ( '\nWater content in surface ' ) 
     619echo ( 'SRF_mas_wat_beg  = {:12.6e} kg | SRF_mas_wat_end  = {:12.6e} kg '.format (SRF_mas_wat_beg, SRF_mas_wat_end) ) 
     620echo ( 'dMass (water srf) = {:12.3e} kg '.format (dSRF_mas_wat) ) 
     621echo ( 'dMass (water srf) = {:12.3e} Sv '.format (dSRF_mas_wat/dtime_sec*1E-6/ATM_rho) ) 
     622echo ( 'dMass (water srf) = {:12.3e} m  '.format (dSRF_mas_wat/ATM_aire_sea_tot/ATM_rho) ) 
     623 
     624echo ( '\nWater content in  ATM + SRF + RUN ' ) 
     625echo ( '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) ) 
     628echo ( 'dMass (water atm+srf+run) = {:12.6e} kg '.format ( dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat) ) 
     629echo ( '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) ) 
     630echo ( '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) ) 
    395631 
    396632echo ( '\n------------------------------------------------------------------------------------' ) 
    397 echo ( '\nWater content in  ATM + SRF + RUN ' ) 
    398 echo ( 'mas_wat_beg              = {:12.6e} kg '.format (ATM_mas_wat_beg + ATM_mas_sno_beg + RUN_mas_wat_beg + SRF_mas_wat_beg) ) 
    399 echo ( 'mas_wat_end              = {:12.6e} kg '.format (ATM_mas_wat_end + ATM_mas_sno_end + RUN_mas_wat_end + SRF_mas_wat_end) ) 
    400 echo ( 'dMass(water atm+srf+run) = {:12.6e} kg '.format ( dATM_mas_wat   + dATM_mas_sno    + dRUN_mas_wat    + dSRF_mas_wat) ) 
    401 echo ( 'dMass(water atm+srf+run) = {:12.3e} Sv '.format ((dATM_mas_wat   + dATM_mas_sno    + dRUN_mas_wat    + dSRF_mas_wat)/dtime_sec*1E-9) ) 
    402 echo ( 'dMass(water atm+srf+run) = {:12.3e} m  '.format ((dATM_mas_wat   + dATM_mas_sno    + dRUN_mas_wat    + dSRF_mas_wat)/ATM_aire_sea_tot*1E-3) ) 
    403  
    404  
    405 echo ( "\n-----------------" ) 
    406 echo ( " Atm -> Oce fluxes" ) 
    407  
    408 ATM_wbilo_sea = np.sum ( lmdz.geo2point (d_ATM_his['wbilo_oce'] + d_ATM_his['wbilo_sic'])*dtime_per_sec*ATM_aire ).values.item() 
    409 ATM_runoff    = np.sum ( lmdz.geo2point (d_ATM_his['runofflic'])*ATM_aire*dtime_per_sec ).values.item() 
    410 ATM_calving   = np.sum ( lmdz.geo2point (d_ATM_his['fqcalving'])*ATM_aire*dtime_per_sec ).values.item() 
    411  
    412 echo (' wbilo oce+sic         = {:12.5e} (kg) '.format ( ATM_wbilo_sea ) ) 
    413 echo (' runoff liq            = {:12.5e} (kg) '.format ( ATM_runoff    ) ) 
    414 echo (' calving               = {:12.5e} (kg) '.format ( ATM_calving   ) ) 
    415 echo (' total                 = {:12.5e} (kg) '.format ( ATM_wbilo_sea - ATM_runoff - ATM_calving ) ) 
     633echo ( '-- ATM Fluxes ' ) 
     634 
     635ATM_wbilo_oce   = lmdz.geo2point ( d_ATM_his ['wbilo_oce']   ) 
     636ATM_wbilo_sic   = lmdz.geo2point ( d_ATM_his ['wbilo_sic']   ) 
     637ATM_wbilo_ter   = lmdz.geo2point ( d_ATM_his ['wbilo_ter']   ) 
     638ATM_wbilo_lic   = lmdz.geo2point ( d_ATM_his ['wbilo_lic']   ) 
     639ATM_runofflic   = lmdz.geo2point ( d_ATM_his ['runofflic']   ) 
     640ATM_fqcalving   = lmdz.geo2point ( d_ATM_his ['fqcalving']   ) 
     641ATM_fqfonte     = lmdz.geo2point ( d_ATM_his ['fqfonte']     ) 
     642ATM_precip      = lmdz.geo2point ( d_ATM_his ['precip']      ) 
     643ATM_snowf       = lmdz.geo2point ( d_ATM_his ['snow']        ) 
     644ATM_evap        = lmdz.geo2point ( d_ATM_his ['evap']        ) 
     645ATM_bilo = ATM_wbilo_oce + ATM_wbilo_ter + ATM_wbilo_lic + ATM_wbilo_sic 
     646 
     647RUN_coastalflow = lmdz.geo2point ( d_RUN_his ['coastalflow'] )  
     648RUN_riverflow   = lmdz.geo2point ( d_RUN_his ['riverflow']   ) 
     649RUN_runoff      = lmdz.geo2point ( d_RUN_his ['runoff']      ) 
     650RUN_drainage    = lmdz.geo2point ( d_RUN_his ['drainage']    ) 
     651RUN_riversret   = lmdz.geo2point ( d_RUN_his ['riversret']   ) 
     652 
     653SRF_evap     = lmdz.geo2point ( d_SRF_his ['evap'] ) 
     654SRF_snowf    = lmdz.geo2point ( d_SRF_his ['snowf'] ) 
     655SRF_TWS      = lmdz.geo2point ( d_SRF_his ['TWS'] ) 
     656SRF_subli    = lmdz.geo2point ( d_SRF_his ['subli'] ) 
     657SRF_transpir = lmdz.geo2point ( np.sum(d_SRF_his ['transpir'], axis=1) ) ; SRF_transpir.attrs['units'] = d_SRF_his ['transpir'].attrs['units'] 
     658 
     659def mmd2SI ( Var) : 
     660    '''Change unit from mm/d or m^3/s to kg/s if needed''' 
     661    if 'units' in VarT.attrs :  
     662        if VarT.attrs['units'] in ['m^3/s', 'm3/s', 'm3.s-3'] : 
     663            VarT.values = VarT.values  * ATM_rho                 ;  VarT.attrs['units'] = 'kg/s' 
     664        if VarT.attrs['units'] == 'mm/d' : 
     665            VarT.values = VarT.values  * ATM_rho * (1e-3/86400.) ;  VarT.attrs['units'] = 'kg/s' 
     666 
     667for var in ['runoff', 'drainage', 'riversret', 'coastalflow', 'riverflow'] : 
     668    VarT = locals()['RUN_' + var] 
     669    mmd2SI (VarT) 
     670 
     671for var in ['evap', 'snowf', 'TWS', 'subli', 'transpir'] : 
     672    VarT = locals()['SRF_' + var] 
     673    mmd2SI (VarT) 
     674 
     675RUN_input  = RUN_runoff      + RUN_drainage 
     676RUN_output = RUN_coastalflow + RUN_riverflow 
     677 
     678ATM_wbilo_sea = ATM_wbilo_oce + ATM_wbilo_sic 
     679 
     680ATM_flx_oce       = ATM_flux_int ( ATM_wbilo_oce  ) 
     681ATM_flx_sic       = ATM_flux_int ( ATM_wbilo_sic  ) 
     682ATM_flx_sea       = ATM_flux_int ( ATM_wbilo_sea  ) 
     683ATM_flx_ter       = ATM_flux_int ( ATM_wbilo_ter  ) 
     684ATM_flx_lic       = ATM_flux_int ( ATM_wbilo_lic  ) 
     685ATM_flx_calving   = ATM_flux_int ( ATM_fqcalving  ) 
     686ATM_flx_qfonte    = ATM_flux_int ( ATM_fqfonte    ) 
     687ATM_flx_precip    = ATM_flux_int ( ATM_precip     ) 
     688ATM_flx_snowf     = ATM_flux_int ( ATM_snowf      ) 
     689ATM_flx_evap      = ATM_flux_int ( ATM_evap       ) 
     690ATM_flx_runlic    = ATM_flux_int ( ATM_runofflic  ) 
     691 
     692RUN_flx_coastal   = ONE_flux_int ( RUN_coastalflow) 
     693RUN_flx_river     = ONE_flux_int ( RUN_riverflow  ) 
     694RUN_flx_drainage  = ATM_flux_int ( RUN_drainage   ) 
     695RUN_flx_riversret = ATM_flux_int ( RUN_riversret  ) 
     696RUN_flx_runoff    = ATM_flux_int ( RUN_runoff     ) 
     697RUN_flx_input     = ATM_flux_int ( RUN_input      ) 
     698RUN_flx_output    = ONE_flux_int ( RUN_output     ) 
     699 
     700ATM_flx_emp   = ATM_flx_evap - ATM_flx_precip 
     701ATM_flx_bilo  = ATM_flx_oce + ATM_flx_sic + ATM_flx_ter + ATM_flx_lic 
     702RUN_flx_bil   = RUN_flx_input - RUN_flx_output 
     703 
     704RUN_flx_rivcoa = RUN_flx_coastal + RUN_flx_river 
     705 
     706ATM_flx_bilo2 = ATM_flux_int (ATM_bilo) 
     707 
     708 
     709echo ('  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 )) 
     710echo ('  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 )) 
     711echo ('  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 )) 
     712echo ('  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 )) 
     713echo ('  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 )) 
     714echo ('  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 )) 
     715echo ('  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 )) 
     716echo ('  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 )) 
     717echo ('  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 )) 
     718echo ('  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 )) 
     719echo ('  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 )) 
     720echo ('  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 )) 
     721echo ('  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 )) 
     722echo ('  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 )) 
     723echo ('  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 )) 
     724echo ('  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 )) 
     725echo ('  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 )) 
     726echo ('  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 )) 
     727echo ('  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 )) 
     728echo ('  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 )) 
     729echo ('  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 )) 
     730echo ('  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 
     732ATM_flx_budget = -ATM_flx_sea + ATM_flx_calving + ATM_flx_runlic #+ ATM_flx_qfonte + RUN_flx_river  
     733echo ('\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 
     738echo ( '------------------------------------------------------------------------------------' ) 
     739echo ( '-- SECHIBA fluxes' ) 
     740 
     741 
     742SRF_flx_evap     = ATM_flux_int ( SRF_evap     ) 
     743SRF_flx_snowf    = ATM_flux_int ( SRF_snowf    ) 
     744SRF_flx_subli    = ATM_flux_int ( SRF_subli    ) 
     745SRF_flx_transpir = ATM_flux_int ( SRF_transpir ) 
     746 
     747echo ('  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 )) 
     748echo ('  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 )) 
     749echo ('  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 )) 
     750echo ('  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 )) 

  • TOOLS/WATER_BUDGET/OCE_waterbudget.py

    r6265 r6271  
    1919#  $HeadURL$ 
    2020 
    21  
    22 ## Define Experiment 
    23 ## ----------------- 
    24  
    25 if False :  
    26     JobName="TEST-CM72-SIMPLE-ROUTING.12" ; TagName="IPSLCM7" ; SpaceName="DEVT" ; ExperimentName="piControl" ; User="p86caub" ; Project="gencmip6"  
    27     Freq = 'MO' ; YearBegin = 1970 ; YearEnd = 1979 ; YearInc = 10 ; PackFrequency=10 
    28     #Freq = 'MO' ; YearBegin = 1852 ; YearEnd = 1852 ; PackFrequency = 1 
    29     ORCA = 'eORCA1.2' ; NEMO=3.6 ; OCE_relax = False ; OCE_icb = False ; Coupled = True 
    30  
    31 if False :  
    32     JobName="TEST-CM72-SIMPLE-ROUTING.10" ; TagName="IPSLCM7" ; SpaceName="DEVT" ; ExperimentName="piControl" ; User="p86caub" ; Project="gencmip6"  
    33     Freq = 'MO' ; YearBegin = 1860 ; YearEnd = 1869 ; PackFrequency = 10 
    34     ORCA = 'eORCA1.2' ; NEMO=3.6 ; OCE_relax = False ; OCE_icb = False ; Coupled = True 
    35  
    36 if False :  
    37     JobName="VALID-CM622-LR.01" ; TagName="IPSLCM6" ; SpaceName="DEVT" ; ExperimentName="piControl" ; User="p86caub" ; Project="gencmip6"  
    38     Freq = 'MO' ; YearBegin = 1890 ; YearEnd = 1899 ; PackFrequency = 10 
    39     ORCA = 'eORCA1.2' ; NEMO=3.6 ; OCE_relax = False ; OCE_icb = False ; Coupled = True 
    40  
    41 if False :  
    42     JobName="VALID-CM622-SIMPLE-ROUTING" ; TagName="IPSLCM6" ; SpaceName="DEVT" ; ExperimentName="piControl" ; User="p86caub" ; Project="gencmip6"  
    43     Freq = 'MO' ; YearBegin = 1898 ; YearEnd = 1898 ; PackFrequency = 1 
    44     ORCA = 'eORCA1.2' ; NEMO=3.6 ; OCE_relax = False ; OCE_icb = False ; Coupled = True 
    45      
    46 if False :  
    47     JobName="VALID-CM622-NORESTART" ; TagName="IPSLCM6" ; SpaceName="DEVT" ; ExperimentName="piControl" ; User="p86caub" ; Project="gencmip6"  
    48     Freq = 'MO' ; YearBegin = 1890 ; YearEnd = 1899 ; PackFrequency = 10 
    49     ORCA = 'eORCA1.2' ; NEMO=3.6 ; OCE_relax = False ; OCE_icb = False ; Coupled = True 
    50  
    51 if True : 
    52     JobName="CM65v420-LR-SKL-pi-05" ; TagName="IPSLCM6" ; SpaceName="DEVT" ; ExperimentName="piControl" ; User="p48ethe" ; Project="gencmip6"  
    53     Freq = 'MO' ; YearBegin = 1890 ; YearEnd = 1899 ; PackFrequency = 10 
    54     ORCA = 'eORCA1.2' ; NEMO=4.2 ; OCE_relax = False ; OCE_icb = False ; Coupled = True 
    55     #/ccc/store/cont003/gencmip6/p48ethe/IGCM_OUT/IPSLCM6/DEVT/piControl/CM65v420-LR-SKL-pi-05 
    56  
     21### 
     22## Import system modules 
     23import sys, os, shutil, subprocess 
    5724import numpy as np 
     25import configparser, re 
     26 
     27## Creates parser 
     28config = configparser.ConfigParser() 
     29config.optionxform = str # To keep capitals 
     30 
     31config['Files']  = {} 
     32config['System'] = {} 
     33 
    5834##-- Some physical constants 
    5935#-- Earth Radius 
    … …  
    7046ICE_rho_pnd = 1000. 
    7147 
    72 ## 
    73 ## Import system modules 
    74 import sys, os, shutil, subprocess 
     48## Read experiment parameters 
     49ATM = None ; ORCA = None ; NEMO = None ; OCE_relax = False ; OCE_icb = False ; Coupled = False ; Routing = None 
     50 
     51# Arguments passed 
     52print ( "Name of Python script:", sys.argv[0] ) 
     53IniFile =  sys.argv[1] 
     54print ("Input file : ", IniFile ) 
     55config.read (IniFile) 
     56 
     57def setBool (chars) : 
     58    '''Convert specific char string in boolean if possible''' 
     59    setBool = chars 
     60    for key in configparser.ConfigParser.BOOLEAN_STATES.keys () : 
     61        if chars.lower() == key : setBool = configparser.ConfigParser.BOOLEAN_STATES[key] 
     62    return setBool 
     63 
     64def setNum (chars) : 
     65    '''Convert specific char string in integer or real if possible''' 
     66    if type (chars) == str : 
     67        realnum = re.compile ("^[-+]?[0-9]*\.?[0-9]+(e[-+]?[0-9]+)?$") 
     68        isReal = realnum.match(chars.strip()) != None 
     69        isInt  = chars.strip().isdigit() 
     70        if isReal : 
     71            if isInt : setNum = int   (chars) 
     72            else     : setNum = float (chars) 
     73        else : setNum = chars 
     74    else : setNum = chars 
     75    return setNum 
     76 
     77print ('[Experiment]') 
     78for 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]) ) 
    7583 
    7684# Where do we run ? 
    77 TGCC = False ; IDRIS = False 
    78 SysName, NodeName, Release, Version, Machine = os.uname () 
    79 if 'irene'   in NodeName : TGCC  = True 
    80 if 'jeanzay' in NodeName : IDRIS = True 
     85SysName, NodeName, Release, Version, Machine = os.uname() 
     86TGCC  = ( 'irene'   in NodeName ) 
     87IDRIS = ( 'jeanzay' in NodeName ) 
    8188 
    8289## Set site specific libIGCM directories, and other specific stuff 
    … …  
    154161 
    155162#-- Files Names 
    156 if Freq == 'MO' : File = f'{JobName}_{YearBegin}0101_{YearEnd}1231_1M' 
    157 if Freq == 'SE' : File = f'{JobName}_SE_{YearBegin}0101_{YearEnd}1231_1M' 
     163if Freq == 'MO' : FileCommon = f'{JobName}_{YearBegin}0101_{YearEnd}1231_1M' 
     164if Freq == 'SE' : FileCommon = f'{JobName}_SE_{YearBegin}0101_{YearEnd}1231_1M' 
    158165 
    159166echo ('\nOpen history files' ) 
    160 file_OCE_his = os.path.join ( dir_OCE_his, f'{File}_grid_T.nc' ) 
    161 file_OCE_sca = os.path.join ( dir_OCE_his, f'{File}_scalar.nc' ) 
    162 file_ICE_his = os.path.join ( dir_ICE_his, f'{File}_icemod.nc' ) 
     167file_OCE_his = os.path.join ( dir_OCE_his, f'{FileCommon}_grid_T.nc' ) 
     168file_OCE_sca = os.path.join ( dir_OCE_his, f'{FileCommon}_scalar.nc' ) 
     169file_ICE_his = os.path.join ( dir_ICE_his, f'{FileCommon}_icemod.nc' ) 
    163170 
    164171d_OCE_his = xr.open_dataset ( file_OCE_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 
    … …  
    177184 
    178185## Compute length of each period 
    179 dtime_per = (d_OCE_his.time_counter_bounds[:,-1] - d_OCE_his.time_counter_bounds[:,0] ) 
     186dtime_per = ( d_OCE_his.time_counter_bounds[:,-1] - d_OCE_his.time_counter_bounds[:,0] ) 
    180187echo ('\nPeriods lengths (days) : ') 
    181188echo (' {:}'.format ( (dtime_per/np.timedelta64(1, "D")).values ) ) 
    … …  
    220227        echo ( command ) 
    221228        os.system ( command ) 
    222     print ('1.2') 
    223229    echo ('extract ndomain' ) 
    224230    ndomain_opa = get_ndomain ( os.path.join (TmpDirOCE, f'OCE_{JobName}_{YearRes}1231_restart') ) 
    … …  
    287293ICE_aire = OCE_aire 
    288294 
    289 OCE_aire_tot = np.sum (OCE_aire).values.item() 
    290 ICE_aire_tot = np.sum (ICE_aire).values.item() 
    291  
    292  
     295def OCE_stock_int (stock) : 
     296    '''Integrate stock on ocean grid''' 
     297    OCE_stock_int = np.sum (  np.sort ( (stock * OCE_aire ).to_masked_array().ravel()) ) 
     298    return OCE_stock_int 
     299 
     300def ONE_stock_int (stock) : 
     301    '''Sum stock''' 
     302    ONE_stock_int = np.sum (  np.sort ( (stock * OCE_aire).to_masked_array().ravel()) ) 
     303    return ONE_stock_int 
     304 
     305def OCE_flux_int (flux) : 
     306    '''Integrate flux on oce grid''' 
     307    OCE_flux_int = np.sum (  np.sort ( (flux * OCE_aire * dtime_per_sec).to_masked_array().ravel()) ) 
     308    return OCE_flux_int 
     309 
     310OCE_aire_tot = ONE_stock_int ( OCE_aire ) 
     311ICE_aire_tot = ONE_stock_int ( ICE_aire ) 
     312     
    293313echo ( '\n------------------------------------------------------------------------------------' ) 
    294314echo ( '-- NEMO change in stores (for the records)' ) 
    … …  
    300320OCE_ssh_beg = d_OCE_beg['sshn'] 
    301321OCE_ssh_end = d_OCE_end['sshn'] 
    302 OCE_sum_ssh_beg = np.sum ( OCE_ssh_beg * OCE_aire ).values.item() 
    303 OCE_sum_ssh_end = np.sum ( OCE_ssh_end * OCE_aire ).values.item() 
     322 
     323OCE_sum_ssh_beg = OCE_stock_int ( OCE_ssh_beg ) 
     324OCE_sum_ssh_end = OCE_stock_int ( OCE_ssh_end ) 
    304325 
    305326if NEMO == 3.6 : 
    306327    OCE_e3tn_beg = d_OCE_beg['fse3t_n'] 
    307328    OCE_e3tn_end = d_OCE_end['fse3t_n'] 
    308     OCE_sum_e3tn_beg = np.sum ( OCE_e3tn_beg * OCE_aire * OCE_msk3).values.item() 
    309     OCE_sum_e3tn_end = np.sum ( OCE_e3tn_end * OCE_aire * OCE_msk3).values.item() 
     329    OCE_sum_e3tn_beg = OCE_stock_int ( OCE_e3tn_beg * OCE_msk3) 
     330    OCE_sum_e3tn_end = OCE_stock_int ( OCE_e3tn_end * OCE_msk3) 
    310331 
    311332echo ( 'OCE_sum_ssh_beg = {:12.6e} m^3  - OCE_sum_ssh_end = {:12.6e} m^3'.format (OCE_sum_ssh_beg, OCE_sum_ssh_end) ) 
    … …  
    329350    echo ( 'dOCE e3tn   mass   = {:12.3e} kg '.format ( dOCE_e3tn_mas) ) 
    330351 
    331  
    332352## Glace et neige 
    333353if NEMO == 3.6 : 
    … …  
    348368     
    349369if NEMO == 4.0 or NEMO == 4.2 : 
    350     ICE_ice_beg = d_ICE_beg['v_i'] 
    351     ICE_ice_end = d_ICE_end['v_i'] 
    352  
    353     ICE_sno_beg = d_ICE_beg['v_s'] 
    354     ICE_sno_end = d_ICE_end['v_s'] 
    355  
    356     ICE_pnd_beg = d_ICE_beg['v_ip'] # Ice ponds 
    357     ICE_pnd_end = d_ICE_end['v_ip'] 
    358  
    359     ICE_fzl_beg = d_ICE_beg['v_il'] # Frozen Ice Ponds 
    360     ICE_fzl_end = d_ICE_end['v_il'] 
    361      
    362     ICE_mas_wat_beg =  np.sum ( d_ICE_beg['snwice_mass']*ICE_aire ).values.item() 
    363     ICE_mas_wat_end =  np.sum ( d_ICE_end['snwice_mass']*ICE_aire ).values.item() 
    364      
    365      
    366 ICE_vol_ice_beg = np.sum ( ICE_ice_beg*ICE_aire ).values.item() 
    367 ICE_vol_ice_end = np.sum ( ICE_ice_end*ICE_aire ).values.item() 
    368  
    369 ICE_vol_sno_beg = np.sum ( ICE_sno_beg*ICE_aire ).values.item() 
    370 ICE_vol_sno_end = np.sum ( ICE_sno_end*ICE_aire ).values.item() 
    371  
    372 ICE_vol_pnd_beg = np.sum ( ICE_pnd_beg*ICE_aire ).values.item() 
    373 ICE_vol_pnd_end = np.sum ( ICE_pnd_end*ICE_aire ).values.item() 
    374  
    375 ICE_vol_fzl_beg = np.sum ( ICE_fzl_beg*ICE_aire ).values.item() 
    376 ICE_vol_fzl_end = np.sum ( ICE_fzl_end*ICE_aire ).values.item() 
     370    ICE_ice_beg = d_ICE_beg['v_i']  ; ICE_ice_end = d_ICE_end['v_i'] 
     371    ICE_sno_beg = d_ICE_beg['v_s']  ; ICE_sno_end = d_ICE_end['v_s']  
     372    ICE_pnd_beg = d_ICE_beg['v_ip'] ; ICE_pnd_end = d_ICE_end['v_ip'] # Ice ponds 
     373    ICE_fzl_beg = d_ICE_beg['v_il'] ; ICE_fzl_end = d_ICE_end['v_il'] # Frozen Ice Ponds 
     374     
     375    ICE_mas_wat_beg =  OCE_stock_int ( d_ICE_beg['snwice_mass'] ) 
     376    ICE_mas_wat_end =  OCE_stock_int ( d_ICE_end['snwice_mass'] ) 
     377     
     378ICE_vol_ice_beg = OCE_stock_int ( ICE_ice_beg ) 
     379ICE_vol_ice_end = OCE_stock_int ( ICE_ice_end ) 
     380 
     381ICE_vol_sno_beg = OCE_stock_int ( ICE_sno_beg ) 
     382ICE_vol_sno_end = OCE_stock_int ( ICE_sno_end ) 
     383 
     384ICE_vol_pnd_beg = OCE_stock_int ( ICE_pnd_beg ) 
     385ICE_vol_pnd_end = OCE_stock_int ( ICE_pnd_end ) 
     386 
     387ICE_vol_fzl_beg = OCE_stock_int ( ICE_fzl_beg ) 
     388ICE_vol_fzl_end = OCE_stock_int ( ICE_fzl_end ) 
    377389 
    378390#-- Converting to freswater volume 
    … …  
    397409    dSEA_mas_wat = dOCE_mas_wat + dICE_mas_wat 
    398410 
    399 echo ( 'ICE_vol_ice_beg = {:12.6e} m^3 - ICE_vol_ice_end = {:12.6e} m^3'.format (ICE_vol_ice_beg, ICE_vol_ice_end) ) 
    400 echo ( 'ICE_vol_sno_beg = {:12.6e} m^3 - ICE_vol_sno_end = {:12.6e} m^3'.format (ICE_vol_sno_beg, ICE_vol_sno_end) ) 
    401 echo ( 'ICE_vol_pnd_beg = {:12.6e} m^3 - ICE_vol_pnd_end = {:12.6e} m^3'.format (ICE_vol_pnd_beg, ICE_vol_pnd_end) ) 
    402 echo ( 'ICE_vol_fzl_beg = {:12.6e} m^3 - ICE_vol_fzl_end = {:12.6e} m^3'.format (ICE_vol_fzl_beg, ICE_vol_fzl_end) ) 
     411echo ( 'ICE_vol_ice_beg = {:12.6e} m^3 | ICE_vol_ice_end = {:12.6e} m^3'.format (ICE_vol_ice_beg, ICE_vol_ice_end) ) 
     412echo ( 'ICE_vol_sno_beg = {:12.6e} m^3 | ICE_vol_sno_end = {:12.6e} m^3'.format (ICE_vol_sno_beg, ICE_vol_sno_end) ) 
     413echo ( 'ICE_vol_pnd_beg = {:12.6e} m^3 | ICE_vol_pnd_end = {:12.6e} m^3'.format (ICE_vol_pnd_beg, ICE_vol_pnd_end) ) 
     414echo ( 'ICE_vol_fzl_beg = {:12.6e} m^3 | ICE_vol_fzl_end = {:12.6e} m^3'.format (ICE_vol_fzl_beg, ICE_vol_fzl_end) ) 
    403415 
    404416echo ( 'dICE_vol_ice   = {:12.3e} m^3'.format (dICE_vol_ice) ) 
    … …  
    410422echo ( 'dICE_mas_fzl   = {:12.3e} m^3'.format (dICE_mas_fzl) )   
    411423 
    412  
    413424echo ( '\n------------------------------------------------------------') 
    414425echo ( 'Variation du contenu en eau ocean + glace ' ) 
    415 echo ( 'dMass(ocean)   = {:12.6e} kg '.format(dSEA_mas_wat) ) 
    416 echo ( 'dVol (ocean)   = {:12.3e} Sv '.format(dSEA_mas_wat/dtime_sec*1E-9) ) 
    417 echo ( 'dVol (ocean)   = {:12.3e} m  '.format(dSEA_mas_wat*1E-3/OCE_aire_tot) ) 
     426echo ( 'dMass (ocean)   = {:12.6e} kg '.format(dSEA_mas_wat) ) 
     427echo ( 'dVol  (ocean)   = {:12.3e} Sv '.format(dSEA_mas_wat/dtime_sec*1E-9) ) 
     428echo ( 'dVol  (ocean)   = {:12.3e} m  '.format(dSEA_mas_wat*1E-3/OCE_aire_tot) ) 
    418429 
    419430 
    … …  
    438449echo ( '\n------------------------------------------------------------------------------------' ) 
    439450echo ( '-- Checks in NEMO - from budget_modipsl.sh (Clement Rousset)' ) 
    440  
    441 def OCE_flux_int (flux, aire=OCE_aire, dt=dtime_per_sec) : 
    442     '''Integrate flux on ocean grid''' 
    443     OCE_flux_int  = np.sum ( flux * aire * dt).values.item() 
    444     #OCE_flux_int  = np.sum ( flux.to_masked_array()[:,:] * aire.to_masked_array()[np.newaxis,:,:] * dt.to_masked_array()[:,np.newaxis,np.newaxis]) 
    445     return OCE_flux_int 
    446451 
    447452# Read variable and computes integral over space and time 
    … …  
    496501OCE_mas_emp = OCE_mas_emp_oce - OCE_mas_wfxice - OCE_mas_wfxsnw - ICE_mas_wfxpnd - ICE_mas_wfxsub_err 
    497502 
     503OCE_mas_all = OCE_mas_emp_oce +OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf 
    498504 
    499505echo ('\n-- Fields:' )  
    500 echo ('  EMPMR      {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_empmr     , OCE_mas_empmr      / dtime_sec*1E-9 )) 
    501 echo ('  WFOB       {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfob      , OCE_mas_wfob       / dtime_sec*1E-9 )) 
    502 echo ('  EMP_OCE    {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_emp_oce   , OCE_mas_emp_oce    / dtime_sec*1E-9 )) 
    503 echo ('  EMP_ICE    {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_emp_ice   , OCE_mas_emp_ice    / dtime_sec*1E-9 )) 
    504 echo ('  ICEBERG    {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_iceberg   , OCE_mas_iceberg    / dtime_sec*1E-9 )) 
    505 echo ('  ICESHELF   {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_iceshelf  , OCE_mas_iceshelf   / dtime_sec*1E-9 )) 
    506 echo ('  CALVING    {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_calving   , OCE_mas_calving    / dtime_sec*1E-9 )) 
    507 echo ('  FRIVER     {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_friver    , OCE_mas_friver     / dtime_sec*1E-9 ))   
    508 echo ('  RUNOFFS    {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_runoffs   , OCE_mas_runoffs    / dtime_sec*1E-9 )) 
    509 echo ('  WFXICE     {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfxice    , OCE_mas_wfxice     / dtime_sec*1E-9 )) 
    510 echo ('  WFXSNW     {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfxsnw    , OCE_mas_wfxsnw     / dtime_sec*1E-9 )) 
    511 echo ('  WFXSUB     {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfxsub    , OCE_mas_wfxsub     / dtime_sec*1E-9 )) 
    512 echo ('  WFXPND     {:12.3e} (kg) {:12.3f} (Sv) '.format ( ICE_mas_wfxpnd    , ICE_mas_wfxpnd     / dtime_sec*1E-9 )) 
    513 echo ('  WFXSNW_SUB {:12.3e} (kg) {:12.3f} (Sv) '.format ( ICE_mas_wfxsnw_sub, ICE_mas_wfxsnw_sub / dtime_sec*1E-9 )) 
    514 echo ('  WFXSNW_PRE {:12.3e} (kg) {:12.3f} (Sv) '.format ( ICE_mas_wfxsnw_pre, ICE_mas_wfxsnw_pre / dtime_sec*1E-9 )) 
    515 echo ('  WFXSUB_ERR {:12.3e} (kg) {:12.3e} (Sv) '.format ( ICE_mas_wfxsub_err, ICE_mas_wfxsub_err / dtime_sec*1E-9 )) 
    516 echo ('  EVAP_OCE   {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_evap_oce  , OCE_mas_evap_oce   / dtime_sec*1E-9 )) 
    517 echo ('  EVAP_ICE   {:12.3e} (kg) {:12.3f} (Sv) '.format ( ICE_mas_evap_ice  , ICE_mas_evap_ice   / dtime_sec*1E-9 )) 
    518 echo ('  SNOW_OCE   {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_snow_oce  , OCE_mas_snow_oce   / dtime_sec*1E-9 )) 
    519 echo ('  SNOW_ICE   {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_snow_ice  , OCE_mas_snow_ice   / dtime_sec*1E-9 )) 
    520 echo ('  RAIN       {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_rain      , OCE_mas_rain       / dtime_sec*1E-9 )) 
    521 echo ('  FWB        {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_fwb       , OCE_mas_fwb        / dtime_sec*1E-9 )) 
    522 echo ('  SSR        {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_ssr       , OCE_mas_ssr        / dtime_sec*1E-9 )) 
    523 echo ('  WFCORR     {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfcorr    , OCE_mas_wfcorr     / dtime_sec*1E-9 )) 
    524 echo ('  BERG_ICB   {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_berg_icb  , OCE_mas_berg_icb   / dtime_sec*1E-9 )) 
    525 echo ('  CALV_ICB   {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_calv_icb  , OCE_mas_calv_icb   / dtime_sec*1E-9 )) 
     506echo ('OCE+ICE budget {:12.3e} (kg) {:12.3e} (Sv) '.format ( OCE_mas_all       , OCE_mas_all        / dtime_sec*1E-9 )) 
     507echo ('  EMPMR        {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_empmr     , OCE_mas_empmr      / dtime_sec*1E-9 )) 
     508echo ('  WFOB         {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfob      , OCE_mas_wfob       / dtime_sec*1E-9 )) 
     509echo ('  EMP_OCE      {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_emp_oce   , OCE_mas_emp_oce    / dtime_sec*1E-9 )) 
     510echo ('  EMP_ICE      {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_emp_ice   , OCE_mas_emp_ice    / dtime_sec*1E-9 )) 
     511echo ('  EMP          {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_emp       , OCE_mas_emp        / dtime_sec*1E-9 )) 
     512echo ('  ICEBERG      {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_iceberg   , OCE_mas_iceberg    / dtime_sec*1E-9 )) 
     513echo ('  ICESHELF     {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_iceshelf  , OCE_mas_iceshelf   / dtime_sec*1E-9 )) 
     514echo ('  CALVING      {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_calving   , OCE_mas_calving    / dtime_sec*1E-9 )) 
     515echo ('  FRIVER       {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_friver    , OCE_mas_friver     / dtime_sec*1E-9 ))   
     516echo ('  RUNOFFS      {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_runoffs   , OCE_mas_runoffs    / dtime_sec*1E-9 )) 
     517echo ('  WFXICE       {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfxice    , OCE_mas_wfxice     / dtime_sec*1E-9 )) 
     518echo ('  WFXSNW       {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfxsnw    , OCE_mas_wfxsnw     / dtime_sec*1E-9 )) 
     519echo ('  WFXSUB       {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfxsub    , OCE_mas_wfxsub     / dtime_sec*1E-9 )) 
     520echo ('  WFXPND       {:12.3e} (kg) {:12.3f} (Sv) '.format ( ICE_mas_wfxpnd    , ICE_mas_wfxpnd     / dtime_sec*1E-9 )) 
     521echo ('  WFXSNW_SUB   {:12.3e} (kg) {:12.3f} (Sv) '.format ( ICE_mas_wfxsnw_sub, ICE_mas_wfxsnw_sub / dtime_sec*1E-9 )) 
     522echo ('  WFXSNW_PRE   {:12.3e} (kg) {:12.3f} (Sv) '.format ( ICE_mas_wfxsnw_pre, ICE_mas_wfxsnw_pre / dtime_sec*1E-9 )) 
     523echo ('  WFXSUB_ERR   {:12.3e} (kg) {:12.3e} (Sv) '.format ( ICE_mas_wfxsub_err, ICE_mas_wfxsub_err / dtime_sec*1E-9 )) 
     524echo ('  EVAP_OCE     {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_evap_oce  , OCE_mas_evap_oce   / dtime_sec*1E-9 )) 
     525echo ('  EVAP_ICE     {:12.3e} (kg) {:12.3f} (Sv) '.format ( ICE_mas_evap_ice  , ICE_mas_evap_ice   / dtime_sec*1E-9 )) 
     526echo ('  SNOW_OCE     {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_snow_oce  , OCE_mas_snow_oce   / dtime_sec*1E-9 )) 
     527echo ('  SNOW_ICE     {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_snow_ice  , OCE_mas_snow_ice   / dtime_sec*1E-9 )) 
     528echo ('  RAIN         {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_rain      , OCE_mas_rain       / dtime_sec*1E-9 )) 
     529echo ('  FWB          {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_fwb       , OCE_mas_fwb        / dtime_sec*1E-9 )) 
     530echo ('  SSR          {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_ssr       , OCE_mas_ssr        / dtime_sec*1E-9 )) 
     531echo ('  WFCORR       {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_wfcorr    , OCE_mas_wfcorr     / dtime_sec*1E-9 )) 
     532echo ('  BERG_ICB     {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_berg_icb  , OCE_mas_berg_icb   / dtime_sec*1E-9 )) 
     533echo ('  CALV_ICB     {:12.3e} (kg) {:12.3f} (Sv) '.format ( OCE_mas_calv_icb  , OCE_mas_calv_icb   / dtime_sec*1E-9 )) 
    526534 
    527535echo     ('\n===>  Comparing ===>' )  
Note: See TracChangeset for help on using the changeset viewer.