Ignore:
Timestamp:
12/06/23 10:13:00 (7 months ago)
Author:
omamce
Message:

O.M. :

WATER_BUDGET : put function common to ATM and OCE in WAterUtils library

use of dictionary for all input parameters

File:
1 edited

Legend:

Unmodified
Added
Removed
  • TOOLS/WATER_BUDGET/ATM_waterbudget.py

    r6676 r6688  
    2525import sys 
    2626import os 
     27import types 
    2728import configparser 
    2829 
     
    5152print ("Output file : ", FullIniFile ) 
    5253 
    53 ## Experiment parameters 
    54 ## -------------------- 
     54## Experiment parameters : read in .ini file 
     55## ----------------------------------------- 
    5556dpar = wu.ReadConfig ( IniFile ) 
    5657 
     
    9495def kg2myear (pval, rho=ATM_RHO) : 
    9596    '''From kg to m/year''' 
    96     return pval/ATM_aire_sea_tot/rho/NbYear 
     97    return pval/ATM.aire_sea_tot/rho/NbYear 
    9798 
    9899def var2prt (pvar, small=False, rho=ATM_RHO) : 
     
    121122## ------------------ 
    122123d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 
    123 if SRF : 
     124if SECHIBA : 
    124125    d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 
    125126    if Routing == 'SECHIBA' : d_RUN_his = d_SRF_his 
     
    127128 
    128129echo ( f'{file_ATM_his = }' ) 
    129 if SRF : 
     130if SECHIBA : 
    130131    echo ( f'{file_SRF_his = }' ) 
    131132     
     
    135136echo ( f'\nRun length : {(dtime/np.timedelta64(1, "D")).values:8.2f} days' ) 
    136137dtime_sec = (dtime/np.timedelta64(1, "s")).values.item() # Convert in seconds 
     138dpar['Experiment']['dtime_sec'] = dtime_sec 
    137139 
    138140##-- Compute length of each period 
     
    142144dtime_per_sec = xr.DataArray (dtime_per_sec, dims=["time_counter", ], coords=[d_ATM_his.time_counter,] ) 
    143145dtime_per_sec.attrs['unit'] = 's' 
     146dpar['Experiment']['dtime_per_sec'] = dtime_per_sec 
    144147 
    145148##-- Number of years (approximative) 
    146149NbYear = dtime_sec / YEAR_LENGTH 
     150dpar['Experiment']['NbYear'] = NbYear 
    147151 
    148152## Define restart periods and file names 
    149153## ------------------------------------- 
    150 dpar = wu.SetRestartNames ( dpar, f_out)  
     154dpar = wu.SetRestartNames ( dpar, f_out )  
    151155 
    152156## Put dpar values in local namespace 
     
    159163         
    160164## Extract restart files from tar 
    161 ## ---------------------------------- 
    162  
     165## ------------------------------ 
    163166liste_beg = [file_ATM_beg, ] 
    164167liste_end = [file_ATM_end, ] 
     
    173176    dpar['Files']['file_DYN_aire'] = file_DYN_aire 
    174177 
    175 if SRF and Routing == 'SIMPLE' : 
     178if SECHIBA and Routing == 'SIMPLE' : 
    176179    liste_beg.append ( file_RUN_beg ) 
    177180    liste_end.append ( file_RUN_end ) 
    178181 
    179182echo ( '\nExtract restart files from tar : ATM, ICO', end='') 
    180 if SRF : echo ( ' and SRF') 
     183if SECHIBA : echo ( ' and SECHIBA') 
    181184else   : echo (' ') 
    182185 
     
    228231ErrorCount = 0 
    229232 
     233echo ( f'Extract {file_ATM_beg = }' ) 
    230234ErrorCount += extract ( file_name=file_ATM_beg, tar_restart=tar_restart_beg_ATM, file_dir_comp=FileDir ) 
     235echo ( f'Extract {file_DYN_beg = }' ) 
    231236ErrorCount += extract ( file_name=file_DYN_beg, tar_restart=tar_restart_beg_DYN, file_dir_comp=FileDir ) 
    232237 
     238echo ( f'Extract {file_ATM_end = }' ) 
    233239ErrorCount += extract ( file_name=file_ATM_end, tar_restart=tar_restart_end_ATM, file_dir_comp=FileDir ) 
     240echo ( f'Extract {file_DYN_end = }' ) 
    234241ErrorCount += extract ( file_name=file_DYN_end, tar_restart=tar_restart_end_DYN, file_dir_comp=FileDir ) 
    235242 
    236 if SRF : 
     243if SECHIBA : 
     244    echo ( f'Extract {file_SRF_beg = }' ) 
    237245    ErrorCount += extract ( file_name=file_SRF_beg, tar_restart=tar_restart_beg_SRF, file_dir_comp=FileDir ) 
     246    echo ( f'Extract {file_SRF_end = }' ) 
    238247    ErrorCount += extract ( file_name=file_SRF_end, tar_restart=tar_restart_end_SRF, file_dir_comp=FileDir ) 
    239248 
    240249    if Routing == 'SIMPLE' : 
     250        echo ( f'Extract {file_RUN_beg = }' ) 
    241251        ErrorCount += extract ( file_name=file_RUN_beg, tar_restart=tar_restart_beg_RUN, file_dir_comp=FileDir ) 
     252        echo ( f'Extract {file_RUN_end = }' ) 
    242253        ErrorCount += extract ( file_name=file_RUN_end, tar_restart=tar_restart_end_RUN, file_dir_comp=FileDir ) 
    243254 
    244 ##-- Exit in case of error in the opening file phase 
     255##-- Exit in case of error in the extracting file phase 
    245256if ErrorCount > 0 : 
    246257    echo ( '  ' ) 
     
    248259 
    249260## 
    250 echo ('\nOpening ATM SRF and ICO restart files') 
     261echo ('\nOpening ATM SECHIBA and ICO restart files') 
    251262d_ATM_beg = xr.open_dataset ( os.path.join (FileDir, file_ATM_beg), decode_times=False, decode_cf=True ).squeeze() 
    252263d_ATM_end = xr.open_dataset ( os.path.join (FileDir, file_ATM_end), decode_times=False, decode_cf=True ).squeeze() 
    253 if SRF : 
     264if SECHIBA : 
    254265    d_SRF_beg = xr.open_dataset ( os.path.join (FileDir, file_SRF_beg), decode_times=False, decode_cf=True ).squeeze() 
    255266    d_SRF_end = xr.open_dataset ( os.path.join (FileDir, file_SRF_end), decode_times=False, decode_cf=True ).squeeze() 
     
    257268d_DYN_end = xr.open_dataset ( os.path.join (FileDir, file_DYN_end), decode_times=False, decode_cf=True ).squeeze() 
    258269 
    259 if SRF : 
     270if SECHIBA : 
    260271    for var in d_SRF_beg.variables : 
    261272        d_SRF_beg[var] = d_SRF_beg[var].where ( d_SRF_beg[var]<1.e20, 0.) 
     
    277288d_ATM_beg = to_cell ( d_ATM_beg ) 
    278289d_ATM_end = to_cell ( d_ATM_end ) 
    279 if SRF : 
     290if SECHIBA : 
    280291    d_SRF_beg = to_cell ( d_SRF_beg ) 
    281292    d_SRF_end = to_cell ( d_SRF_end ) 
     
    283294d_DYN_end = to_cell ( d_DYN_end ) 
    284295 
    285 if SRF and Routing == 'SIMPLE' : 
     296if SECHIBA and Routing == 'SIMPLE' : 
    286297    d_RUN_beg = to_cell ( d_RUN_beg ) 
    287298    d_RUN_end = to_cell ( d_RUN_end ) 
    288299 
    289300d_ATM_his = to_cell ( d_ATM_his ) 
    290 if SRF : d_SRF_his = to_cell ( d_SRF_his ) 
     301if SECHIBA : d_SRF_his = to_cell ( d_SRF_his ) 
    291302 
    292303echo ( f'{file_ATM_beg = }' ) 
     
    294305echo ( f'{file_DYN_beg = }' ) 
    295306echo ( f'{file_DYN_end = }' ) 
    296 if SRF : 
     307if SECHIBA : 
    297308    echo ( f'{file_SRF_beg = }' ) 
    298309    echo ( f'{file_SRF_end = }' ) 
     
    301312        echo ( f'{file_RUN_end = }' ) 
    302313 
    303 # ATM grid with cell surfaces 
    304 if LMDZ : 
    305     echo ('ATM grid with cell surfaces : LMDZ') 
    306     ATM_lat       = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1d='cell' ) 
    307     ATM_lon       = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1d='cell' ) 
    308     ATM_aire      = lmdz.geo2point ( rprec (d_ATM_his ['aire']     [0]), cumul_poles=True, dim1d='cell' ) 
    309     ATM_fter      = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1d='cell' ) 
    310     ATM_foce      = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1d='cell' ) 
    311     ATM_fsic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1d='cell' ) 
    312     ATM_flic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1d='cell' ) 
    313     if SRF : 
    314         SRF_lat       = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1d='cell' ) 
    315         SRF_lon       = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1d='cell' ) 
    316         SRF_aire      = lmdz.geo2point ( rprec (d_SRF_his ['Areas']) * rprec (d_SRF_his ['Contfrac']), dim1d='cell', cumul_poles=True ) 
    317         SRF_areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas'])  ,  dim1d='cell', cumul_poles=True ) 
    318         SRF_contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac']), dim1d='cell' ) 
    319  
    320 if ICO : 
    321     if ATM_HIS == 'latlon' : 
    322         echo ( 'ATM areas and fractions on LATLON grid' ) 
    323         if 'lat_dom_out' in d_ATM_his.variables : 
    324             ATM_lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat_dom_out'])+0*rprec (d_ATM_his ['lon_dom_out']), dim1d='cell' ) 
    325             ATM_lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat_dom_out'])+  rprec (d_ATM_his ['lon_dom_out']), dim1d='cell' ) 
    326         else : 
    327             ATM_lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1d='cell' ) 
    328             ATM_lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1d='cell' ) 
    329         ATM_aire = lmdz.geo2point ( rprec (d_ATM_his ['aire'][0]).squeeze(), cumul_poles=True, dim1d='cell' ) 
    330         ATM_fter = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1d='cell' ) 
    331         ATM_foce = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1d='cell' ) 
    332         ATM_fsic = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1d='cell' ) 
    333         ATM_flic = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1d='cell' ) 
    334  
    335     if ATM_HIS == 'ico' : 
    336         echo ( 'ATM areas and fractions on ICO grid' ) 
    337         ATM_aire =  rprec (d_ATM_his ['aire']     [0]).squeeze() 
    338         ATM_lat  =  rprec (d_ATM_his ['lat']         ) 
    339         ATM_lon  =  rprec (d_ATM_his ['lon']         ) 
    340         ATM_fter =  rprec (d_ATM_his ['fract_ter'][0]) 
    341         ATM_foce =  rprec (d_ATM_his ['fract_oce'][0]) 
    342         ATM_fsic =  rprec (d_ATM_his ['fract_sic'][0]) 
    343         ATM_flic =  rprec (d_ATM_his ['fract_lic'][0]) 
    344  
    345     if SRF : 
    346         if SRF_HIS == 'latlon' : 
    347             echo ( 'SRF areas and fractions on LATLON grid' ) 
    348             if 'lat_domain_landpoints_out' in d_SRF_his  : 
    349                 SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat_domain_landpoints_out'])+0*rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1d='cell' ) 
    350                 SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_domain_landpoints_out'])+  rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1d='cell' ) 
    351             else : 
    352                 if 'lat_domain_landpoints_out' in d_SRF_his : 
    353                     SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat_dom_out'])+0*rprec (d_SRF_his ['lon_dom_out']), dim1d='cell' ) 
    354                     SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_dom_out'])+  rprec (d_SRF_his ['lon_dom_out']), dim1d='cell' ) 
    355                 else : 
    356                     SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1d='cell' ) 
    357                     SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1d='cell' ) 
    358  
    359             SRF_areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas']   )   , dim1d='cell', cumul_poles=True ) 
    360             SRF_areafrac  = lmdz.geo2point ( rprec (d_SRF_his ['AreaFrac'])   , dim1d='cell', cumul_poles=True ) 
    361             SRF_contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac'])   , dim1d='cell', cumul_poles=True ) 
    362             SRF_aire = SRF_areafrac 
    363  
    364         if SRF_HIS == 'ico' : 
    365             echo ( 'SRF areas and fractions on ICO grid' ) 
    366             SRF_lat       =  rprec (d_SRF_his ['lat']     ) 
    367             SRF_lon       =  rprec (d_SRF_his ['lon']     ) 
    368             SRF_areas     =  rprec (d_SRF_his ['Areas']   ) 
    369             SRF_contfrac  =  rprec (d_SRF_his ['Contfrac']) 
    370             SRF_aire      =  SRF_areas * SRF_contfrac 
    371  
    372 ATM_fsea      = ATM_foce + ATM_fsic 
    373 ATM_flnd      = ATM_fter + ATM_flic 
    374 ATM_aire_fter = ATM_aire * ATM_fter 
    375 ATM_aire_flic = ATM_aire * ATM_flic 
    376 ATM_aire_fsic = ATM_aire * ATM_fsic 
    377 ATM_aire_foce = ATM_aire * ATM_foce 
    378 ATM_aire_flnd = ATM_aire * ATM_flnd 
    379 ATM_aire_fsea = ATM_aire * ATM_fsea 
    380  
    381 #SRF_aire = SRF_aire.where ( np.abs (SRF_aire) < 1E15, 0. ) 
    382  
     314## Compute aire, fractions, etc ... 
     315## -------------------------------- 
    383316if ICO : 
    384317    # Area on icosahedron grid 
    385     echo ( f'{file_DYN_aire = }' ) 
     318    #echo ( f'{file_DYN_aire = }' ) 
    386319    d_DYN_aire = xr.open_dataset ( file_DYN_aire, decode_times=False ).squeeze() 
    387      
    388     if SortIco : 
    389             # Creation d'une clef de tri pour le fichier aire 
    390             DYN_aire_keysort = np.lexsort ( (d_DYN_aire['lat'], d_DYN_aire['lon']) ) 
    391     else : 
    392             DYN_aire_keysort = np.arange ( len ( d_DYN_aire['lat'] ) ) 
    393              
    394     DYN_lat = d_DYN_aire['lat'] 
    395     DYN_lon = d_DYN_aire['lon'] 
     320else : 
     321    d_DYN_aire = None 
    396322         
    397     DYN_aire = d_DYN_aire['aire'] 
    398     DYN_fsea = d_DYN_aire['fract_oce'] + d_DYN_aire['fract_sic'] 
    399          
    400     DYN_flnd = 1.0 - DYN_fsea 
    401     DYN_fter = d_ATM_beg['FTER'] 
    402     DYN_flic = d_ATM_beg['FLIC'] 
    403     DYN_foce = d_ATM_beg['FOCE'] 
    404     DYN_aire_fter = DYN_aire * DYN_fter 
    405  
    406 #if ATM_HIS == 'ico' : 
    407 #    ATM_aire      = DYN_aire 
    408 #    ATM_aire_fter = DYN_aire * ATM_fter 
    409 #    ATM_aire_flic = DYN_aire * ATM_flic 
    410 #    ATM_aire_fsic = DYN_aire * ATM_fsic 
    411 #    ATM_aire_foce = DYN_aire * ATM_foce 
    412 #    ATM_aire_flnd = DYN_aire * ATM_flnd 
    413 #    ATM_aire_fsea = DYN_aire * ATM_fsea 
    414  
    415 if ATM_HIS == 'ico' : 
    416     DYN_aire      = ATM_aire 
    417     DYN_foce      = ATM_foce 
    418     DYN_fsic      = ATM_fsic 
    419     DYN_flic      = ATM_flic 
    420     DYN_fter      = ATM_fter 
    421     DYN_fsea      = ATM_fsea 
    422     DYN_flnd      = ATM_flnd 
    423     DYN_aire_fter = ATM_aire_fter 
    424     DYN_aire_flic = ATM_aire_flic 
    425     DYN_aire_fsic = ATM_aire_fsic 
    426     DYN_aire_foce = ATM_aire_foce 
    427     DYN_aire_flnd = ATM_aire_flnd 
    428     DYN_aire_fsea = ATM_aire_fsea 
    429  
    430 if LMDZ : 
    431     # Area on lon/lat grid 
    432     DYN_aire = ATM_aire 
    433     DYN_fsea = ATM_fsea 
    434     DYN_flnd = ATM_flnd 
    435     DYN_fter = rprec (d_ATM_beg['FTER']) 
    436     DYN_flic = rprec (d_ATM_beg['FLIC']) 
    437     DYN_aire_fter = DYN_aire * DYN_fter 
     323dpar, ATM = wu.ComputeGridATM ( dpar, d_ATM_his )  
     324dpar, SRF = wu.ComputeGridSRF ( dpar, d_SRF_his ) 
     325dpar, DYN = wu.ComputeGridDYN ( dpar, ATM, d_DYN_aire, d_ATM_beg ) 
    438326 
    439327if ICO and ATM_HIS == 'ico' : 
    440328    # Comparaison des aires sur ATM et DYN 
    441     aire_diff = ATM_aire - DYN_aire 
    442     echo ( 'f{Difference Aire hist file vs. grid file {aire_diff.mean()=} {aire_diff.min()=}  {aire_diff.max()=} ' ) 
    443      
     329    aire_diff = ATM.aire - DYN.aire 
     330    echo ( f'Difference Aire hist file vs. grid file mean:{aire_diff.mean().values} min:{aire_diff.min().values} max:{aire_diff.max().values} ' ) 
    444331 
    445332# Functions computing integrals and sum 
     
    447334def DYN_stock_int (stock) : 
    448335    '''Integrate (* surface) stock on atmosphere grid''' 
    449     return wu.Psum ( (stock * DYN_aire).to_masked_array().ravel() ) 
     336    return wu.P1sum ( stock * DYN.aire ) 
    450337 
    451338@Timer 
    452339def ATM_flux_int (flux) : 
    453340    '''Integrate (* time * surface) flux on atmosphere grid''' 
    454     return wu.Psum ( (flux * dtime_per_sec * ATM_aire).to_masked_array().ravel() ) 
     341    return wu.P1sum ( flux * dtime_per_sec * ATM.aire ) 
    455342 
    456343@Timer 
    457344def LIC_flux_int (flux) : 
    458345    '''Integrate (* time * surface) flux on land ice grid''' 
    459     return wu.Psum ( (flux * dtime_per_sec * ATM_aire_flic).to_masked_array().ravel() ) 
    460  
    461 if SRF : 
     346    return wu.P1sum ( flux * dtime_per_sec * ATM.aire_flic ) 
     347 
     348if SECHIBA : 
    462349    @Timer 
    463350    def SRF_stock_int (stock) : 
    464351        '''Integrate (* surface) stock on land grid''' 
    465         return wu.Ksum (  ( (stock * DYN_aire_fter).to_masked_array().ravel()) ) 
     352        return wu.P1sum ( stock * DYN.aire_fter ) 
    466353 
    467354    @Timer 
    468355    def SRF_flux_int (flux) : 
    469356        '''Integrate (* time * surface) flux on land grid''' 
    470         return wu.Psum (  (flux * dtime_per_sec * SRF_aire).to_masked_array().ravel() ) 
     357        return wu.P1sum ( flux * dtime_per_sec * SRF.aire ) 
    471358 
    472359@Timer 
    473360def ONE_stock_int (stock) : 
    474361    '''Sum stock''' 
    475     return wu.Psum ( stock.to_masked_array().ravel() ) 
     362    return wu.P1sum ( stock )  
    476363 
    477364@Timer 
     
    485372    return wu.Psum ( np.array ( tlist) ) 
    486373 
    487 ATM_aire_sea     = ATM_aire * ATM_fsea 
    488  
    489 ATM_aire_tot     = ONE_stock_int (ATM_aire) 
    490 ATM_aire_sea_tot = ONE_stock_int (ATM_aire_fsea) 
    491 ATM_aire_ter_tot = ONE_stock_int (ATM_aire_fter) 
    492 ATM_aire_lic_tot = ONE_stock_int (ATM_aire_flic) 
    493  
    494 DYN_aire_tot     = ONE_stock_int ( DYN_aire ) 
    495 if SRF : SRF_aire_tot     = ONE_stock_int ( SRF_aire ) 
    496  
    497374echo ('') 
    498 echo ( f'ATM DYN : Area of atmosphere        : {DYN_aire_tot:18.9e}' ) 
    499 echo ( f'ATM HIS : Area of atmosphere        : {ATM_aire_tot:18.9e}' ) 
    500 echo ( f'ATM HIS : Area of ter in atmosphere : {ATM_aire_ter_tot:18.9e}' ) 
    501 echo ( f'ATM HIS : Area of lic in atmosphere : {ATM_aire_lic_tot:18.9e}' ) 
    502 if SRF : 
    503     echo ( f'ATM SRF : Area of atmosphere        : {SRF_aire_tot:18.9e}' ) 
     375echo ( f'ATM DYN : Area of atmosphere        : {DYN.aire_tot:18.9e}' ) 
     376echo ( f'ATM HIS : Area of atmosphere        : {ATM.aire_tot:18.9e}' ) 
     377echo ( f'ATM HIS : Area of ter in atmosphere : {ATM.aire_ter_tot:18.9e}' ) 
     378echo ( f'ATM HIS : Area of lic in atmosphere : {ATM.aire_lic_tot:18.9e}' ) 
     379if SECHIBA : 
     380    echo ( f'ATM SECHIBA : Area of atmosphere        : {SRF.aire_tot:18.9e}' ) 
    504381echo ('') 
    505 echo ( 'ATM DYN : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(DYN_aire_tot/(RA*RA*4*np.pi) ) ) 
    506 echo ( 'ATM HIS : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(ATM_aire_tot/(RA*RA*4*np.pi) ) ) 
    507 echo ( 'ATM HIS : Area of ter in atmosphere/(4pi R^2) : {:18.9f}'.format(ATM_aire_ter_tot/(RA*RA*4*np.pi) ) ) 
    508 if SRF : 
    509     echo ( 'ATM SRF : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(SRF_aire_tot/(RA*RA*4*np.pi) ) ) 
     382echo ( 'ATM DYN : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(DYN.aire_tot/(RA*RA*4*np.pi) ) ) 
     383echo ( 'ATM HIS : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(ATM.aire_tot/(RA*RA*4*np.pi) ) ) 
     384echo ( 'ATM HIS : Area of ter in atmosphere/(4pi R^2) : {:18.9f}'.format(ATM.aire_ter_tot/(RA*RA*4*np.pi) ) ) 
     385if SECHIBA : 
     386    echo ( 'ATM SECHIBA : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(SRF.aire_tot/(RA*RA*4*np.pi) ) ) 
    510387    echo ('') 
    511     echo ( f'ATM SRF : Area of atmosphere (no contfrac): {ONE_stock_int (SRF_areas):18.9e}' ) 
    512  
    513  
    514 if np.abs (ATM_aire_tot/(RA*RA*4*np.pi) - 1.0) > 0.01 : 
     388    echo ( f'ATM SECHIBA : Area of atmosphere (no contfrac): {ONE_stock_int (SRF.areas):18.9e}' ) 
     389 
     390 
     391if np.abs (ATM.aire_tot/(RA*RA*4*np.pi) - 1.0) > 0.01 : 
    515392    raise RuntimeError ('Error of atmosphere surface interpolated on lon/lat grid') 
    516393 
     
    522399 
    523400echo ( 'ATM vertical grid' ) 
    524 ATM_Ahyb = d_ATM_his['Ahyb'].squeeze() 
    525 ATM_Bhyb = d_ATM_his['Bhyb'].squeeze() 
     401ATM.Ahyb = d_ATM_his['Ahyb'].squeeze() 
     402ATM.Bhyb = d_ATM_his['Bhyb'].squeeze() 
    526403 
    527404echo ( 'Surface pressure' ) 
    528405if ICO : 
    529     DYN_psol_beg = d_DYN_beg['ps'] 
    530     DYN_psol_end = d_DYN_end['ps'] 
     406    DYN.psol_beg = d_DYN_beg['ps'] 
     407    DYN.psol_end = d_DYN_end['ps'] 
    531408if LMDZ : 
    532     DYN_psol_beg = lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)), dim1d='cell' ) 
    533     DYN_psol_end = lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)), dim1d='cell' ) 
    534  
    535 echo ( '3D Pressure at the interface layers (not scalar points)' ) 
    536 DYN_pres_beg = ATM_Ahyb + ATM_Bhyb * DYN_psol_beg 
    537 DYN_pres_end = ATM_Ahyb + ATM_Bhyb * DYN_psol_end 
     409    DYN.psol_beg = lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)), dim1d='cell' ) 
     410    DYN.psol_end = lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)), dim1d='cell' ) 
     411 
     412echo ( '3D Pressure at the interface layers (not at scalar points)' ) 
     413DYN.pres_beg = ATM.Ahyb + ATM.Bhyb * DYN.psol_beg 
     414DYN.pres_end = ATM.Ahyb + ATM.Bhyb * DYN.psol_end 
    538415 
    539416echo ( 'Check computation of pressure levels' ) 
     
    541418ind = np.empty (8) 
    542419 
    543 ind[0] = (DYN_pres_beg[ 0]-DYN_psol_beg).min().item() 
    544 ind[1] = (DYN_pres_beg[ 0]-DYN_psol_beg).max().item() 
    545 ind[2] = (DYN_pres_beg[-1]).min().item() 
    546 ind[3] = (DYN_pres_beg[-1]).max().item() 
    547 ind[4] = (DYN_pres_end[ 0]-DYN_psol_end).min().item() 
    548 ind[5] = (DYN_pres_end[ 0]-DYN_psol_end).max().item() 
    549 ind[6] = (DYN_pres_end[-1]).min().item() 
    550 ind[7] = (DYN_pres_end[-1]).max().item() 
     420ind[0] = (DYN.pres_beg[ 0]-DYN.psol_beg).min().item() 
     421ind[1] = (DYN.pres_beg[ 0]-DYN.psol_beg).max().item() 
     422ind[2] = (DYN.pres_beg[-1]).min().item() 
     423ind[3] = (DYN.pres_beg[-1]).max().item() 
     424ind[4] = (DYN.pres_end[ 0]-DYN.psol_end).min().item() 
     425ind[5] = (DYN.pres_end[ 0]-DYN.psol_end).max().item() 
     426ind[6] = (DYN.pres_end[-1]).min().item() 
     427ind[7] = (DYN.pres_end[-1]).max().item() 
    551428 
    552429if any ( ind != 0) : 
    553430    echo ( 'All values should be zero' ) 
    554     echo ( f'(DYN_pres_beg[ 0]-DYN_psol_beg).min().item() = {ind[0]}' ) 
    555     echo ( f'(DYN_pres_beg[ 0]-DYN_psol_beg).max().item() = {ind[1]}' ) 
    556     echo ( f'(DYN_pres_beg[-1]).min().item()              = {ind[2]}' ) 
    557     echo ( f'(DYN_pres_beg[-1]).max().item()              = {ind[3]}' ) 
    558     echo ( f'(DYN_pres_end[ 0]-DYN_psol_end).min().item() = {ind[4]}' ) 
    559     echo ( f'(DYN_pres_end[ 0]-DYN_psol_end).max().item() = {ind[5]}' ) 
    560     echo ( f'(DYN_pres_end[-1]).min().item()              = {ind[6]}' ) 
    561     echo ( f'(DYN_pres_end[-1]).max().item()              = {ind[7]}' ) 
     431    echo ( f'(DYN.pres_beg[ 0]-DYN.psol_beg).min().item() = {ind[0]}' ) 
     432    echo ( f'(DYN.pres_beg[ 0]-DYN.psol_beg).max().item() = {ind[1]}' ) 
     433    echo ( f'(DYN.pres_beg[-1]).min().item()              = {ind[2]}' ) 
     434    echo ( f'(DYN.pres_beg[-1]).max().item()              = {ind[3]}' ) 
     435    echo ( f'(DYN.pres_end[ 0]-DYN.psol_end).min().item() = {ind[4]}' ) 
     436    echo ( f'(DYN.pres_end[ 0]-DYN.psol_end).max().item() = {ind[5]}' ) 
     437    echo ( f'(DYN.pres_end[-1]).min().item()              = {ind[6]}' ) 
     438    echo ( f'(DYN.pres_end[-1]).max().item()              = {ind[7]}' ) 
    562439    raise RuntimeError 
    563440 
    564 klevp1 = ATM_Bhyb.shape[-1] 
    565 cell   = DYN_psol_beg.shape[-1] 
     441klevp1 = ATM.Bhyb.shape[-1] 
     442cell   = DYN.psol_beg.shape[-1] 
    566443klev   = klevp1 - 1 
    567444 
    568 echo ( 'Layer thickness (pressure)' ) 
    569 DYN_mass_beg = xr.DataArray ( np.empty( (klev, cell)), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) ) 
    570 DYN_mass_end = xr.DataArray ( np.empty( (klev, cell)), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) ) 
    571  
     445echo ( 'Layer thickness (in pressure)' ) 
     446DYN.mass_beg = xr.DataArray ( np.empty( (klev, cell)), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) ) 
     447DYN.mass_end = xr.DataArray ( np.empty( (klev, cell)), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) ) 
     448 
     449echo ( 'Layer mass (kg/m2)' ) 
    572450for k in np.arange (klev) : 
    573     DYN_mass_beg[k,:] = ( DYN_pres_beg[k,:] - DYN_pres_beg[k+1,:] ) / GRAV 
    574     DYN_mass_end[k,:] = ( DYN_pres_end[k,:] - DYN_pres_end[k+1,:] ) / GRAV 
    575  
    576 DYN_mass_beg_2D = DYN_mass_beg.sum (dim='sigs') 
    577 DYN_mass_end_2D = DYN_mass_end.sum (dim='sigs') 
    578  
    579 DYN_mas_air_beg = DYN_stock_int ( DYN_mass_beg_2D ) 
    580 DYN_mas_air_end = DYN_stock_int ( DYN_mass_end_2D ) 
     451    DYN.mass_beg[k,:] = ( DYN.pres_beg[k,:] - DYN.pres_beg[k+1,:] ) / GRAV 
     452    DYN.mass_end[k,:] = ( DYN.pres_end[k,:] - DYN.pres_end[k+1,:] ) / GRAV 
     453 
     454echo ( 'Column mass (kg/m2)' ) 
     455DYN.mass_beg_2D = DYN.mass_beg.sum (dim='sigs') 
     456DYN.mass_end_2D = DYN.mass_end.sum (dim='sigs') 
     457 
     458echo ( 'Atmosphere mass (kg)' ) 
     459DYN.mass_air_beg = DYN_stock_int ( DYN.mass_beg_2D ) 
     460DYN.mass_air_end = DYN_stock_int ( DYN.mass_end_2D ) 
     461 
     462echo ( 'Atmosphere mass change (kg)' ) 
     463DYN.diff_mass_air = DYN.mass_air_end - DYN.mass_air_beg 
     464 
     465echo ( f'\nChange of atmosphere mass (dynamics)  -- {Title} ' ) 
     466echo ( '------------------------------------------------------------------------------------' ) 
     467echo ( 'DYN.mass_air_beg = {:12.6e} kg | DYN.mass_air_end = {:12.6e} kg'.format (DYN.mass_air_beg, DYN.mass_air_end) ) 
     468echo ( f'dMass (atm) : {DYN.diff_mass_air:9.4e} kg' ) 
     469echo ( f'dMass (atm) : {(DYN.diff_mass_air/DYN.mass_air_beg):9.4e} (-)' ) 
    581470 
    582471echo ( 'Vertical and horizontal integral, and sum of liquid, solid and vapor water phases' ) 
     
    584473    if 'H2Ov' in d_DYN_beg.variables : 
    585474        echo ('reading LATLON : H2Ov, H2Ol, H2Oi' ) 
    586         DYN_wat_beg = lmdz.geo3point ( d_DYN_beg['H2Ov'] + d_DYN_beg['H2Ol'] + d_DYN_beg['H2Oi'].isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
    587         DYN_wat_end = lmdz.geo3point ( d_DYN_end['H2Ov'] + d_DYN_end['H2Ol'] + d_DYN_end['H2Oi'].isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
    588     if 'H2Ov_g' in d_DYN_beg.variables : 
     475        DYN.wat_beg = lmdz.geo3point ( (d_DYN_beg['H2Ov']   + d_DYN_beg['H2Ol']  + d_DYN_beg['H2Oi']).isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
     476        DYN.wat_end = lmdz.geo3point ( (d_DYN_end['H2Ov']   + d_DYN_end['H2Ol']  + d_DYN_end['H2Oi']).isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
     477    if 'H2O_g' in d_DYN_beg.variables : 
    589478        echo ('reading LATLON : H2O_g, H2O_l, H2O_s' ) 
    590         DYN_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s']).isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
    591         DYN_wat_end = lmdz.geo3point ( (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s']).isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
     479        DYN.wat_beg = lmdz.geo3point ( (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s']).isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
     480        DYN.wat_end = lmdz.geo3point ( (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s']).isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
    592481if ICO : 
    593     if 'H2Ov_g' in d_DYN_beg.variables : 
     482    if   'H2Ov_g' in d_DYN_beg.variables : 
    594483        echo ('reading ICO : H2Ov_g, H2Ov_l, H2Ov_s' ) 
    595         DYN_wat_beg = d_DYN_beg['H2Ov_g'] + d_DYN_beg['H2Ov_l'] + d_DYN_beg['H2Ov_s'] 
    596         DYN_wat_end = d_DYN_end['H2Ov_g'] + d_DYN_end['H2Ov_l'] + d_DYN_end['H2Ov_s'] 
    597     elif 'H2O_g' in d_DYN_beg.variables : 
     484        DYN.wat_beg = d_DYN_beg['H2Ov_g'] + d_DYN_beg['H2Ov_l'] + d_DYN_beg['H2Ov_s'] 
     485        DYN.wat_end = d_DYN_end['H2Ov_g'] + d_DYN_end['H2Ov_l'] + d_DYN_end['H2Ov_s'] 
     486    if 'H2O_g' in d_DYN_beg.variables : 
    598487        echo ('reading ICO : H2O_g, H2O_l, H2O_s' ) 
    599         DYN_wat_beg = d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l']   + d_DYN_beg['H2O_s'] 
    600         DYN_wat_end = d_DYN_end['H2O_g'] + d_DYN_end['H2O_l']   + d_DYN_end['H2O_s'] 
    601     elif 'q' in d_DYN_beg.variables : 
     488        DYN.wat_beg = d_DYN_beg['H2O_g']  + d_DYN_beg['H2O_l']  + d_DYN_beg['H2O_s'] 
     489        DYN.wat_end = d_DYN_end['H2O_g']  + d_DYN_end['H2O_l']  + d_DYN_end['H2O_s'] 
     490    if 'q' in d_DYN_beg.variables : 
    602491        echo ('reading ICO : q' ) 
    603         DYN_wat_beg = d_DYN_beg['q'].isel(nq=0) + d_DYN_beg['q'].isel(nq=1) + d_DYN_beg['q'].isel(nq=2) 
    604         DYN_wat_end = d_DYN_end['q'].isel(nq=0) + d_DYN_end['q'].isel(nq=1) + d_DYN_end['q'].isel(nq=2) 
    605  
    606 if 'lev' in DYN_wat_beg.dims : 
    607     DYN_wat_beg = DYN_wat_beg.rename ( {'lev':'sigs'} ) 
    608     DYN_wat_end = DYN_wat_end.rename ( {'lev':'sigs'} ) 
     492        DYN.wat_beg = d_DYN_beg['q'].isel(nq=0) + d_DYN_beg['q'].isel(nq=1) + d_DYN_beg['q'].isel(nq=2) 
     493        DYN.wat_end = d_DYN_end['q'].isel(nq=0) + d_DYN_end['q'].isel(nq=1) + d_DYN_end['q'].isel(nq=2) 
     494 
     495if 'lev' in DYN.wat_beg.dims : 
     496    DYN.wat_beg = DYN.wat_beg.rename ( {'lev':'sigs'} ) 
     497    DYN.wat_end = DYN.wat_end.rename ( {'lev':'sigs'} ) 
    609498 
    610499echo ( 'Compute water content : vertical and horizontal integral' ) 
    611500 
    612 DYN_wat_beg_2D =  (DYN_mass_beg * DYN_wat_beg).sum (dim='sigs') 
    613 DYN_wat_end_2D =  (DYN_mass_end * DYN_wat_end).sum (dim='sigs') 
    614  
    615 DYN_mas_wat_beg = DYN_stock_int ( DYN_wat_beg_2D ) 
    616 DYN_mas_wat_end = DYN_stock_int ( DYN_wat_end_2D ) 
     501DYN.wat_beg_2D =  (DYN.mass_beg * DYN.wat_beg).sum (dim='sigs') 
     502DYN.wat_end_2D =  (DYN.mass_end * DYN.wat_end).sum (dim='sigs') 
     503 
     504DYN.mass_wat_beg = DYN_stock_int ( DYN.wat_beg_2D ) 
     505DYN.mass_wat_end = DYN_stock_int ( DYN.wat_end_2D ) 
    617506 
    618507echo ( 'Variation of water content' ) 
    619 dDYN_mas_wat = DYN_mas_wat_end - DYN_mas_wat_beg 
     508DYN.diff_mass_wat = DYN.mass_wat_end - DYN.mass_wat_beg 
    620509 
    621510echo ( f'\nChange of atmosphere water content (dynamics)  -- {Title} ' ) 
    622511echo ( '------------------------------------------------------------------------------------' ) 
    623 echo ( 'DYN_mas_air_beg = {:12.6e} kg | DYN_mas_air_end = {:12.6e} kg'.format (DYN_mas_air_beg, DYN_mas_air_end) ) 
    624 echo ( 'DYN_mas_wat_beg = {:12.6e} kg | DYN_mas_wat_end = {:12.6e} kg'.format (DYN_mas_wat_beg, DYN_mas_wat_end) ) 
    625 prtFlux ( 'dMass (atm)  ', dDYN_mas_wat, 'e', True ) 
    626  
    627 ATM_sno_beg = d_ATM_beg['SNOW01']*d_ATM_beg['FTER'] + d_ATM_beg['SNOW02']*d_ATM_beg['FLIC'] + \ 
     512echo ( 'DYN.mass_wat_beg = {:12.6e} kg | DYN.mass_wat_end = {:12.6e} kg'.format (DYN.mass_wat_beg, DYN.mass_wat_end) ) 
     513prtFlux ( 'dMass water (atm)', DYN.diff_mass_wat, 'e', True ) 
     514echo (   f'dMass water (atm) : {(DYN.diff_mass_wat/DYN.mass_wat_beg):9.4e} (-)' ) 
     515 
     516 
     517ATM.sno_beg = d_ATM_beg['SNOW01']*d_ATM_beg['FTER'] + d_ATM_beg['SNOW02']*d_ATM_beg['FLIC'] + \ 
    628518              d_ATM_beg['SNOW03']*d_ATM_beg['FOCE'] + d_ATM_beg['SNOW04']*d_ATM_beg['FSIC'] 
    629 ATM_sno_end = d_ATM_end['SNOW01']*d_ATM_end['FTER'] + d_ATM_end['SNOW02']*d_ATM_end['FLIC'] + \ 
     519ATM.sno_end = d_ATM_end['SNOW01']*d_ATM_end['FTER'] + d_ATM_end['SNOW02']*d_ATM_end['FLIC'] + \ 
    630520              d_ATM_end['SNOW03']*d_ATM_end['FOCE'] + d_ATM_end['SNOW04']*d_ATM_end['FSIC'] 
    631521 
    632 ATM_qs_beg  = d_ATM_beg['QS01']  *d_ATM_beg['FTER'] + d_ATM_beg['QS02']  *d_ATM_beg['FLIC'] + \ 
     522ATM.qs_beg  = d_ATM_beg['QS01']  *d_ATM_beg['FTER'] + d_ATM_beg['QS02']  *d_ATM_beg['FLIC'] + \ 
    633523              d_ATM_beg['QS03']  *d_ATM_beg['FOCE'] + d_ATM_beg['QS04']  *d_ATM_beg['FSIC'] 
    634 ATM_qs_end  = d_ATM_end['QS01']  *d_ATM_end['FTER'] + d_ATM_end['QS02']  *d_ATM_end['FLIC'] + \ 
     524ATM.qs_end  = d_ATM_end['QS01']  *d_ATM_end['FTER'] + d_ATM_end['QS02']  *d_ATM_end['FLIC'] + \ 
    635525              d_ATM_end['QS03']  *d_ATM_end['FOCE'] + d_ATM_end['QS04']  *d_ATM_end['FSIC'] 
    636526 
    637 ATM_qsol_beg     = d_ATM_beg['QSOL'] 
    638 ATM_qsol_end     = d_ATM_end['QSOL'] 
    639  
    640 LIC_sno_beg      = d_ATM_beg['SNOW02']*d_ATM_beg['FLIC'] 
    641 LIC_sno_end      = d_ATM_end['SNOW02']*d_ATM_end['FLIC'] 
    642  
    643 LIC_runlic0_beg  = d_ATM_beg['RUNOFFLIC0'] 
    644 LIC_runlic0_end  = d_ATM_end['RUNOFFLIC0'] 
    645  
    646 ATM_qs01_beg     = d_ATM_beg['QS01'] * d_ATM_beg['FTER'] 
    647 ATM_qs02_beg     = d_ATM_beg['QS02'] * d_ATM_beg['FLIC'] 
    648 ATM_qs03_beg     = d_ATM_beg['QS03'] * d_ATM_beg['FOCE'] 
    649 ATM_qs04_beg     = d_ATM_beg['QS04'] * d_ATM_beg['FSIC'] 
    650  
    651 ATM_qs01_end     = d_ATM_end['QS01'] * d_ATM_end['FTER'] 
    652 ATM_qs02_end     = d_ATM_end['QS02'] * d_ATM_end['FLIC'] 
    653 ATM_qs03_end     = d_ATM_end['QS03'] * d_ATM_end['FOCE'] 
    654 ATM_qs04_end     = d_ATM_end['QS04'] * d_ATM_end['FSIC'] 
    655  
    656 LIC_qs_beg = ATM_qs02_beg 
    657 LIC_qs_end = ATM_qs02_end 
    658  
    659 ATM_mas_sno_beg     = DYN_stock_int ( ATM_sno_beg  ) 
    660 ATM_mas_sno_end     = DYN_stock_int ( ATM_sno_end  ) 
    661  
    662 ATM_mas_qs_beg      = DYN_stock_int ( ATM_qs_beg   ) 
    663 ATM_mas_qs_end      = DYN_stock_int ( ATM_qs_end   ) 
    664 ATM_mas_qsol_beg    = DYN_stock_int ( ATM_qsol_beg ) 
    665 ATM_mas_qs01_beg    = DYN_stock_int ( ATM_qs01_beg ) 
    666 ATM_mas_qs02_beg    = DYN_stock_int ( ATM_qs02_beg ) 
    667 ATM_mas_qs03_beg    = DYN_stock_int ( ATM_qs03_beg ) 
    668 ATM_mas_qs04_beg    = DYN_stock_int ( ATM_qs04_beg ) 
    669 ATM_mas_qsol_end    = DYN_stock_int ( ATM_qsol_end ) 
    670 ATM_mas_qs01_end    = DYN_stock_int ( ATM_qs01_end ) 
    671 ATM_mas_qs02_end    = DYN_stock_int ( ATM_qs02_end ) 
    672 ATM_mas_qs03_end    = DYN_stock_int ( ATM_qs03_end ) 
    673 ATM_mas_qs04_end    = DYN_stock_int ( ATM_qs04_end ) 
    674  
    675 LIC_mas_sno_beg     = DYN_stock_int ( LIC_sno_beg     ) 
    676 LIC_mas_sno_end     = DYN_stock_int ( LIC_sno_end     ) 
    677 LIC_mas_runlic0_beg = DYN_stock_int ( LIC_runlic0_beg ) 
    678 LIC_mas_runlic0_end = DYN_stock_int ( LIC_runlic0_end ) 
    679  
    680 LIC_mas_qs_beg  = ATM_mas_qs02_beg 
    681 LIC_mas_qs_end  = ATM_mas_qs02_end 
    682 LIC_mas_wat_beg = LIC_mas_qs_beg + LIC_mas_sno_beg 
    683 LIC_mas_wat_end = LIC_mas_qs_end + LIC_mas_sno_end 
    684  
    685 dATM_mas_sno  = ATM_mas_sno_end  - ATM_mas_sno_beg 
    686 dATM_mas_qs   = ATM_mas_qs_end   - ATM_mas_qs_beg 
    687 dATM_mas_qsol = ATM_mas_qsol_end - ATM_mas_qsol_beg 
    688  
    689 dATM_mas_qs01 = ATM_mas_qs01_end - ATM_mas_qs01_beg 
    690 dATM_mas_qs02 = ATM_mas_qs02_end - ATM_mas_qs02_beg 
    691 dATM_mas_qs03 = ATM_mas_qs03_end - ATM_mas_qs03_beg 
    692 dATM_mas_qs04 = ATM_mas_qs04_end - ATM_mas_qs04_beg 
    693  
    694 dLIC_mas_qs      = LIC_mas_qs_end      - LIC_mas_qs_beg 
    695 dLIC_mas_sno     = LIC_mas_sno_end     - LIC_mas_sno_beg 
    696 dLIC_mas_runlic0 = LIC_mas_runlic0_end - LIC_mas_runlic0_beg 
    697  
    698 dLIC_mas_wat     = dLIC_mas_qs + dLIC_mas_sno # + dLIC_mas_runlic0 
     527ATM.qsol_beg     = d_ATM_beg['QSOL'] 
     528ATM.qsol_end     = d_ATM_end['QSOL'] 
     529 
     530ATM.LIC_sno_beg      = d_ATM_beg['SNOW02']*d_ATM_beg['FLIC'] 
     531ATM.LIC_sno_end      = d_ATM_end['SNOW02']*d_ATM_end['FLIC'] 
     532 
     533ATM.LIC_runlic0_beg  = d_ATM_beg['RUNOFFLIC0'] 
     534ATM.LIC_runlic0_end  = d_ATM_end['RUNOFFLIC0'] 
     535 
     536ATM.qs01_beg     = d_ATM_beg['QS01'] * d_ATM_beg['FTER'] 
     537ATM.qs02_beg     = d_ATM_beg['QS02'] * d_ATM_beg['FLIC'] 
     538ATM.qs03_beg     = d_ATM_beg['QS03'] * d_ATM_beg['FOCE'] 
     539ATM.qs04_beg     = d_ATM_beg['QS04'] * d_ATM_beg['FSIC'] 
     540 
     541ATM.qs01_end     = d_ATM_end['QS01'] * d_ATM_end['FTER'] 
     542ATM.qs02_end     = d_ATM_end['QS02'] * d_ATM_end['FLIC'] 
     543ATM.qs03_end     = d_ATM_end['QS03'] * d_ATM_end['FOCE'] 
     544ATM.qs04_end     = d_ATM_end['QS04'] * d_ATM_end['FSIC'] 
     545 
     546ATM.LIC_qs_beg = ATM.qs02_beg 
     547ATM.LIC_qs_end = ATM.qs02_end 
     548 
     549ATM.mass_sno_beg     = DYN_stock_int ( ATM.sno_beg  ) 
     550ATM.mass_sno_end     = DYN_stock_int ( ATM.sno_end  ) 
     551 
     552ATM.mass_qs_beg      = DYN_stock_int ( ATM.qs_beg   ) 
     553ATM.mass_qs_end      = DYN_stock_int ( ATM.qs_end   ) 
     554ATM.mass_qsol_beg    = DYN_stock_int ( ATM.qsol_beg ) 
     555ATM.mass_qs01_beg    = DYN_stock_int ( ATM.qs01_beg ) 
     556ATM.mass_qs02_beg    = DYN_stock_int ( ATM.qs02_beg ) 
     557ATM.mass_qs03_beg    = DYN_stock_int ( ATM.qs03_beg ) 
     558ATM.mass_qs04_beg    = DYN_stock_int ( ATM.qs04_beg ) 
     559ATM.mass_qsol_end    = DYN_stock_int ( ATM.qsol_end ) 
     560ATM.mass_qs01_end    = DYN_stock_int ( ATM.qs01_end ) 
     561ATM.mass_qs02_end    = DYN_stock_int ( ATM.qs02_end ) 
     562ATM.mass_qs03_end    = DYN_stock_int ( ATM.qs03_end ) 
     563ATM.mass_qs04_end    = DYN_stock_int ( ATM.qs04_end ) 
     564 
     565ATM.LIC_mass_sno_beg     = DYN_stock_int ( ATM.LIC_sno_beg     ) 
     566ATM.LIC_mass_sno_end     = DYN_stock_int ( ATM.LIC_sno_end     ) 
     567ATM.LIC_mass_runlic0_beg = DYN_stock_int ( ATM.LIC_runlic0_beg ) 
     568ATM.LIC_mass_runlic0_end = DYN_stock_int ( ATM.LIC_runlic0_end ) 
     569 
     570ATM.LIC_mass_qs_beg  = ATM.mass_qs02_beg 
     571ATM.LIC_mass_qs_end  = ATM.mass_qs02_end 
     572ATM.LIC_mass_wat_beg = ATM.LIC_mass_qs_beg + ATM.LIC_mass_sno_beg 
     573ATM.LIC_mass_wat_end = ATM.LIC_mass_qs_end + ATM.LIC_mass_sno_end 
     574 
     575ATM.diff_mass_sno  = ATM.mass_sno_end  - ATM.mass_sno_beg 
     576ATM.diff_mass_qs   = ATM.mass_qs_end   - ATM.mass_qs_beg 
     577ATM.diff_mass_qsol = ATM.mass_qsol_end - ATM.mass_qsol_beg 
     578 
     579ATM.diff_mass_qs01 = ATM.mass_qs01_end - ATM.mass_qs01_beg 
     580ATM.diff_mass_qs02 = ATM.mass_qs02_end - ATM.mass_qs02_beg 
     581ATM.diff_mass_qs03 = ATM.mass_qs03_end - ATM.mass_qs03_beg 
     582ATM.diff_mass_qs04 = ATM.mass_qs04_end - ATM.mass_qs04_beg 
     583 
     584ATM.LIC_diff_mass_qs      = ATM.LIC_mass_qs_end      - ATM.LIC_mass_qs_beg 
     585ATM.LIC_diff_mass_sno     = ATM.LIC_mass_sno_end     - ATM.LIC_mass_sno_beg 
     586ATM.LIC_diff_mass_runlic0 = ATM.LIC_mass_runlic0_end - ATM.LIC_mass_runlic0_beg 
     587 
     588ATM.LIC_diff_mass_wat     = ATM.LIC_diff_mass_qs + ATM.LIC_diff_mass_sno # + ATM.LIC_diff_mass_runlic0 
    699589 
    700590echo ( f'\nChange of atmosphere snow content (Land ice) -- {Title} ' ) 
    701591echo ( '------------------------------------------------------------------------------------' ) 
    702 echo ( 'ATM_mas_sno_beg  = {:12.6e} kg | ATM_mas_sno_end  = {:12.6e} kg'.format (ATM_mas_sno_beg, ATM_mas_sno_end) ) 
    703 prtFlux ( 'dMass (neige atm) ', dATM_mas_sno  , 'e', True  ) 
     592echo ( 'ATM.mass_sno_beg  = {:12.6e} kg | ATM.mass_sno_end  = {:12.6e} kg'.format (ATM.mass_sno_beg, ATM.mass_sno_end) ) 
     593prtFlux ( 'dMass (neige atm) ', ATM.diff_mass_sno  , 'e', True  ) 
    704594 
    705595echo ( f'\nChange of soil humidity content -- {Title} ' ) 
    706596echo ( '------------------------------------------------------------------------------------' ) 
    707 echo ( 'ATM_mas_qs_beg  = {:12.6e} kg | ATM_mas_qs_end  = {:12.6e} kg'.format (ATM_mas_qs_beg, ATM_mas_qs_end) ) 
    708 prtFlux ( 'dMass (neige atm) ', dATM_mas_qs, 'e', True ) 
     597echo ( 'ATM.mass_qs_beg  = {:12.6e} kg | ATM.mass_qs_end  = {:12.6e} kg'.format (ATM.mass_qs_beg, ATM.mass_qs_end) ) 
     598prtFlux ( 'dMass (neige atm) ', ATM.diff_mass_qs, 'e', True ) 
    709599 
    710600echo ( f'\nChange of atmosphere water+snow content -- {Title}  ' ) 
    711601echo ( '------------------------------------------------------------------------------------' ) 
    712 prtFlux ( 'dMass (eau + neige atm) ', dDYN_mas_wat + dATM_mas_sno , 'e', True) 
    713  
    714 if SRF : 
     602prtFlux ( 'dMass (eau + neige atm) ', DYN.diff_mass_wat + ATM.diff_mass_sno , 'e', True) 
     603 
     604if SECHIBA : 
    715605    echo ( '\n====================================================================================' ) 
    716     echo ( f'-- SRF changes  -- {Title} ' ) 
     606    echo ( f'-- SECHIBA changes  -- {Title} ' ) 
    717607 
    718608    if Routing == 'SIMPLE' : 
    719         RUN_mas_wat_fast_beg   = ONE_stock_int ( d_RUN_beg ['fast_reservoir']   ) 
    720         RUN_mas_wat_slow_beg   = ONE_stock_int ( d_RUN_beg ['slow_reservoir']   ) 
    721         RUN_mas_wat_stream_beg = ONE_stock_int ( d_RUN_beg ['stream_reservoir'] ) 
    722         RUN_mas_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  ) 
    723         RUN_mas_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   ) 
    724         RUN_mas_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   ) 
    725  
    726         RUN_mas_wat_fast_end   = ONE_stock_int ( d_RUN_end ['fast_reservoir']   ) 
    727         RUN_mas_wat_slow_end   = ONE_stock_int ( d_RUN_end ['slow_reservoir']   ) 
    728         RUN_mas_wat_stream_end = ONE_stock_int ( d_RUN_end ['stream_reservoir'] ) 
    729  
    730         RUN_mas_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  ) 
    731         RUN_mas_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   ) 
    732         RUN_mas_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   ) 
     609        SRF.RUN_mass_wat_fast_beg   = ONE_stock_int ( d_RUN_beg ['fast_reservoir']   ) 
     610        SRF.RUN_mass_wat_slow_beg   = ONE_stock_int ( d_RUN_beg ['slow_reservoir']   ) 
     611        SRF.RUN_mass_wat_stream_beg = ONE_stock_int ( d_RUN_beg ['stream_reservoir'] ) 
     612        SRF.RUN_mass_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  ) 
     613        SRF.RUN_mass_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   ) 
     614        SRF.RUN_mass_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   ) 
     615 
     616        SRF.RUN_mass_wat_fast_end   = ONE_stock_int ( d_RUN_end ['fast_reservoir']   ) 
     617        SRF.RUN_mass_wat_slow_end   = ONE_stock_int ( d_RUN_end ['slow_reservoir']   ) 
     618        SRF.RUN_mass_wat_stream_end = ONE_stock_int ( d_RUN_end ['stream_reservoir'] ) 
     619 
     620        SRF.RUN_mass_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  ) 
     621        SRF.RUN_mass_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   ) 
     622        SRF.RUN_mass_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   ) 
    733623 
    734624    if Routing == 'SECHIBA' : 
    735         RUN_mas_wat_fast_beg   = ONE_stock_int ( d_SRF_beg ['fastres']   ) 
    736         RUN_mas_wat_slow_beg   = ONE_stock_int ( d_SRF_beg ['slowres']   ) 
    737         RUN_mas_wat_stream_beg = ONE_stock_int ( d_SRF_beg ['streamres'] ) 
    738         RUN_mas_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  ) 
    739         RUN_mas_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   ) 
    740         RUN_mas_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   ) 
    741  
    742         RUN_mas_wat_fast_end   = ONE_stock_int ( d_SRF_end ['fastres']   ) 
    743         RUN_mas_wat_slow_end   = ONE_stock_int ( d_SRF_end ['slowres']   ) 
    744         RUN_mas_wat_stream_end = ONE_stock_int ( d_SRF_end ['streamres'] ) 
    745         RUN_mas_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  ) 
    746         RUN_mas_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   ) 
    747         RUN_mas_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   ) 
    748  
    749     RUN_mas_wat_beg = Sprec ( [RUN_mas_wat_fast_beg , RUN_mas_wat_slow_beg, RUN_mas_wat_stream_beg, 
    750                             RUN_mas_wat_flood_beg, RUN_mas_wat_lake_beg, RUN_mas_wat_pond_beg] ) 
    751  
    752     RUN_mas_wat_end = Sprec ( [RUN_mas_wat_fast_end  , RUN_mas_wat_slow_end , RUN_mas_wat_stream_end, 
    753                             RUN_mas_wat_flood_end , RUN_mas_wat_lake_end , RUN_mas_wat_pond_end] ) 
    754  
    755     dRUN_mas_wat_fast   = RUN_mas_wat_fast_end   - RUN_mas_wat_fast_beg 
    756     dRUN_mas_wat_slow   = RUN_mas_wat_slow_end   - RUN_mas_wat_slow_beg 
    757     dRUN_mas_wat_stream = RUN_mas_wat_stream_end - RUN_mas_wat_stream_beg 
    758     dRUN_mas_wat_flood  = RUN_mas_wat_flood_end  - RUN_mas_wat_flood_beg 
    759     dRUN_mas_wat_lake   = RUN_mas_wat_lake_end   - RUN_mas_wat_lake_beg 
    760     dRUN_mas_wat_pond   = RUN_mas_wat_pond_end   - RUN_mas_wat_pond_beg 
    761  
    762     dRUN_mas_wat        = RUN_mas_wat_end        - RUN_mas_wat_beg 
     625        SRF.RUN_mass_wat_fast_beg   = ONE_stock_int ( d_SRF_beg ['fastres']   ) 
     626        SRF.RUN_mass_wat_slow_beg   = ONE_stock_int ( d_SRF_beg ['slowres']   ) 
     627        SRF.RUN_mass_wat_stream_beg = ONE_stock_int ( d_SRF_beg ['streamres'] ) 
     628        SRF.RUN_mass_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  ) 
     629        SRF.RUN_mass_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   ) 
     630        SRF.RUN_mass_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   ) 
     631 
     632        SRF.RUN_mass_wat_fast_end   = ONE_stock_int ( d_SRF_end ['fastres']   ) 
     633        SRF.RUN_mass_wat_slow_end   = ONE_stock_int ( d_SRF_end ['slowres']   ) 
     634        SRF.RUN_mass_wat_stream_end = ONE_stock_int ( d_SRF_end ['streamres'] ) 
     635        SRF.RUN_mass_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  ) 
     636        SRF.RUN_mass_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   ) 
     637        SRF.RUN_mass_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   ) 
     638 
     639    SRF.RUN_mass_wat_beg = Sprec ( [SRF.RUN_mass_wat_fast_beg , SRF.RUN_mass_wat_slow_beg, SRF.RUN_mass_wat_stream_beg, 
     640                                   SRF.RUN_mass_wat_flood_beg, SRF.RUN_mass_wat_lake_beg, SRF.RUN_mass_wat_pond_beg] ) 
     641 
     642    SRF.RUN_mass_wat_end = Sprec ( [SRF.RUN_mass_wat_fast_end  , SRF.RUN_mass_wat_slow_end , SRF.RUN_mass_wat_stream_end, 
     643                                   SRF.RUN_mass_wat_flood_end , SRF.RUN_mass_wat_lake_end , SRF.RUN_mass_wat_pond_end] ) 
     644 
     645    SRF.RUN_diff_mass_wat_fast   = SRF.RUN_mass_wat_fast_end   - SRF.RUN_mass_wat_fast_beg 
     646    SRF.RUN_diff_mass_wat_slow   = SRF.RUN_mass_wat_slow_end   - SRF.RUN_mass_wat_slow_beg 
     647    SRF.RUN_diff_mass_wat_stream = SRF.RUN_mass_wat_stream_end - SRF.RUN_mass_wat_stream_beg 
     648    SRF.RUN_diff_mass_wat_flood  = SRF.RUN_mass_wat_flood_end  - SRF.RUN_mass_wat_flood_beg 
     649    SRF.RUN_diff_mass_wat_lake   = SRF.RUN_mass_wat_lake_end   - SRF.RUN_mass_wat_lake_beg 
     650    SRF.RUN_diff_mass_wat_pond   = SRF.RUN_mass_wat_pond_end   - SRF.RUN_mass_wat_pond_beg 
     651 
     652    SRF.RUN_diff_mass_wat        = SRF.RUN_mass_wat_end        - SRF.RUN_mass_wat_beg 
    763653 
    764654    echo ( f'\nRunoff reservoirs -- {Title} ') 
    765655    echo (  '------------------------------------------------------------------------------------' ) 
    766     echo ( f'RUN_mas_wat_fast_beg   = {RUN_mas_wat_fast_beg  :12.6e} kg | RUN_mas_wat_fast_end   = {RUN_mas_wat_fast_end  :12.6e} kg ' ) 
    767     echo ( f'RUN_mas_wat_slow_beg   = {RUN_mas_wat_slow_beg  :12.6e} kg | RUN_mas_wat_slow_end   = {RUN_mas_wat_slow_end  :12.6e} kg ' ) 
    768     echo ( f'RUN_mas_wat_stream_beg = {RUN_mas_wat_stream_beg:12.6e} kg | RUN_mas_wat_stream_end = {RUN_mas_wat_stream_end:12.6e} kg ' ) 
    769     echo ( f'RUN_mas_wat_flood_beg  = {RUN_mas_wat_flood_beg :12.6e} kg | RUN_mas_wat_flood_end  = {RUN_mas_wat_flood_end :12.6e} kg ' ) 
    770     echo ( f'RUN_mas_wat_lake_beg   = {RUN_mas_wat_lake_beg  :12.6e} kg | RUN_mas_wat_lake_end   = {RUN_mas_wat_lake_end  :12.6e} kg ' ) 
    771     echo ( f'RUN_mas_wat_pond_beg   = {RUN_mas_wat_pond_beg  :12.6e} kg | RUN_mas_wat_pond_end   = {RUN_mas_wat_pond_end  :12.6e} kg ' ) 
    772     echo ( f'RUN_mas_wat_beg        = {RUN_mas_wat_beg       :12.6e} kg | RUN_mas_wat_end        = {RUN_mas_wat_end       :12.6e} kg ' ) 
     656    echo ( f'SRF.RUN_mass_wat_fast_beg   = {SRF.RUN_mass_wat_fast_beg  :12.6e} kg | SRF.RUN_mass_wat_fast_end   = {SRF.RUN_mass_wat_fast_end  :12.6e} kg ' ) 
     657    echo ( f'SRF.RUN_mass_wat_slow_beg   = {SRF.RUN_mass_wat_slow_beg  :12.6e} kg | SRF.RUN_mass_wat_slow_end   = {SRF.RUN_mass_wat_slow_end  :12.6e} kg ' ) 
     658    echo ( f'SRF.RUN_mass_wat_stream_beg = {SRF.RUN_mass_wat_stream_beg:12.6e} kg | SRF.RUN_mass_wat_stream_end = {SRF.RUN_mass_wat_stream_end:12.6e} kg ' ) 
     659    echo ( f'SRF.RUN_mass_wat_flood_beg  = {SRF.RUN_mass_wat_flood_beg :12.6e} kg | SRF.RUN_mass_wat_flood_end  = {SRF.RUN_mass_wat_flood_end :12.6e} kg ' ) 
     660    echo ( f'SRF.RUN_mass_wat_lake_beg   = {SRF.RUN_mass_wat_lake_beg  :12.6e} kg | SRF.RUN_mass_wat_lake_end   = {SRF.RUN_mass_wat_lake_end  :12.6e} kg ' ) 
     661    echo ( f'SRF.RUN_mass_wat_pond_beg   = {SRF.RUN_mass_wat_pond_beg  :12.6e} kg | SRF.RUN_mass_wat_pond_end   = {SRF.RUN_mass_wat_pond_end  :12.6e} kg ' ) 
     662    echo ( f'SRF.RUN_mass_wat_beg        = {SRF.RUN_mass_wat_beg       :12.6e} kg | SRF.RUN_mass_wat_end        = {SRF.RUN_mass_wat_end       :12.6e} kg ' ) 
    773663 
    774664    echo ( '------------------------------------------------------------------------------------' ) 
    775665 
    776     prtFlux ( 'dMass (fast)   ', dRUN_mas_wat_fast  , 'e', True ) 
    777     prtFlux ( 'dMass (slow)   ', dRUN_mas_wat_slow  , 'e', True ) 
    778     prtFlux ( 'dMass (stream) ', dRUN_mas_wat_stream, 'e', True ) 
    779     prtFlux ( 'dMass (flood)  ', dRUN_mas_wat_flood , 'e', True ) 
    780     prtFlux ( 'dMass (lake)   ', dRUN_mas_wat_lake  , 'e', True ) 
    781     prtFlux ( 'dMass (pond)   ', dRUN_mas_wat_pond  , 'e', True ) 
    782     prtFlux ( 'dMass (all)    ', dRUN_mas_wat       , 'e', True ) 
     666    prtFlux ( 'dMass (fast)   ', SRF.RUN_diff_mass_wat_fast  , 'e', True ) 
     667    prtFlux ( 'dMass (slow)   ', SRF.RUN_diff_mass_wat_slow  , 'e', True ) 
     668    prtFlux ( 'dMass (stream) ', SRF.RUN_diff_mass_wat_stream, 'e', True ) 
     669    prtFlux ( 'dMass (flood)  ', SRF.RUN_diff_mass_wat_flood , 'e', True ) 
     670    prtFlux ( 'dMass (lake)   ', SRF.RUN_diff_mass_wat_lake  , 'e', True ) 
     671    prtFlux ( 'dMass (pond)   ', SRF.RUN_diff_mass_wat_pond  , 'e', True ) 
     672    prtFlux ( 'dMass (all)    ', SRF.RUN_diff_mass_wat       , 'e', True ) 
    783673 
    784674    echo ( f'\nWater content in routing  -- {Title} ' ) 
    785675    echo ( '------------------------------------------------------------------------------------' ) 
    786     echo ( f'RUN_mas_wat_beg = {RUN_mas_wat_end:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg' ) 
    787     prtFlux ( 'dMass (routing) ', dRUN_mas_wat , 'e', True   ) 
     676    echo ( f'SRF.RUN_mass_wat_beg = {SRF.RUN_mass_wat_end:12.6e} kg | SRF.RUN_mass_wat_end = {SRF.RUN_mass_wat_end:12.6e} kg' ) 
     677    prtFlux ( 'dMass (routing) ', SRF.RUN_diff_mass_wat , 'e', True   ) 
    788678 
    789679    echo ( '\n====================================================================================' ) 
    790     print ( 'Reading SRF restart') 
    791     SRF_tot_watveg_beg  = d_SRF_beg['tot_watveg_beg']  ; SRF_tot_watveg_beg  = SRF_tot_watveg_beg .where (SRF_tot_watveg_beg  < 1E15, 0.) 
    792     SRF_tot_watsoil_beg = d_SRF_beg['tot_watsoil_beg'] ; SRF_tot_watsoil_beg = SRF_tot_watsoil_beg.where (SRF_tot_watsoil_beg < 1E15, 0.) 
    793     SRF_snow_beg        = d_SRF_beg['snow_beg']        ; SRF_snow_beg        = SRF_snow_beg       .where (SRF_snow_beg        < 1E15, 0.) 
    794     SRF_lakeres_beg     = d_SRF_beg['lakeres']         ; SRF_lakeres_beg     = SRF_lakeres_beg    .where (SRF_lakeres_beg     < 1E15, 0.) 
    795  
    796     SRF_tot_watveg_end  = d_SRF_end['tot_watveg_beg']  ; SRF_tot_watveg_end  = SRF_tot_watveg_end .where (SRF_tot_watveg_end  < 1E15, 0.) 
    797     SRF_tot_watsoil_end = d_SRF_end['tot_watsoil_beg'] ; SRF_tot_watsoil_end = SRF_tot_watsoil_end.where (SRF_tot_watsoil_end < 1E15, 0.) 
    798     SRF_snow_end        = d_SRF_end['snow_beg']        ; SRF_snow_end        = SRF_snow_end       .where (SRF_snow_end        < 1E15, 0.) 
    799     SRF_lakeres_end     = d_SRF_end['lakeres']         ; SRF_lakeres_end     = SRF_lakeres_end    .where (SRF_lakeres_end     < 1E15, 0.) 
     680    print ( 'Reading SECHIBA restart') 
     681    SRF.tot_watveg_beg  = d_SRF_beg['tot_watveg_beg']  ; SRF.tot_watveg_beg  = SRF.tot_watveg_beg .where (SRF.tot_watveg_beg  < 1E15, 0.) 
     682    SRF.tot_watsoil_beg = d_SRF_beg['tot_watsoil_beg'] ; SRF.tot_watsoil_beg = SRF.tot_watsoil_beg.where (SRF.tot_watsoil_beg < 1E15, 0.) 
     683    SRF.snow_beg        = d_SRF_beg['snow_beg']        ; SRF.snow_beg        = SRF.snow_beg       .where (SRF.snow_beg        < 1E15, 0.) 
     684    SRF.lakeres_beg     = d_SRF_beg['lakeres']         ; SRF.lakeres_beg     = SRF.lakeres_beg    .where (SRF.lakeres_beg     < 1E15, 0.) 
     685 
     686    SRF.tot_watveg_end  = d_SRF_end['tot_watveg_beg']  ; SRF.tot_watveg_end  = SRF.tot_watveg_end .where (SRF.tot_watveg_end  < 1E15, 0.) 
     687    SRF.tot_watsoil_end = d_SRF_end['tot_watsoil_beg'] ; SRF.tot_watsoil_end = SRF.tot_watsoil_end.where (SRF.tot_watsoil_end < 1E15, 0.) 
     688    SRF.snow_end        = d_SRF_end['snow_beg']        ; SRF.snow_end        = SRF.snow_end       .where (SRF.snow_end        < 1E15, 0.) 
     689    SRF.lakeres_end     = d_SRF_end['lakeres']         ; SRF.lakeres_end     = SRF.lakeres_end    .where (SRF.lakeres_end     < 1E15, 0.) 
    800690 
    801691    if LMDZ : 
    802         SRF_tot_watveg_beg  = lmdz.geo2point (SRF_tot_watveg_beg , dim1d='cell') 
    803         SRF_tot_watsoil_beg = lmdz.geo2point (SRF_tot_watsoil_beg, dim1d='cell') 
    804         SRF_snow_beg        = lmdz.geo2point (SRF_snow_beg       , dim1d='cell') 
    805         SRF_lakeres_beg     = lmdz.geo2point (SRF_lakeres_beg    , dim1d='cell') 
    806         SRF_tot_watveg_end  = lmdz.geo2point (SRF_tot_watveg_end , dim1d='cell') 
    807         SRF_tot_watsoil_end = lmdz.geo2point (SRF_tot_watsoil_end, dim1d='cell') 
    808         SRF_snow_end        = lmdz.geo2point (SRF_snow_end       , dim1d='cell') 
    809         SRF_lakeres_end     = lmdz.geo2point (SRF_lakeres_end    , dim1d='cell') 
     692        SRF.tot_watveg_beg  = lmdz.geo2point (SRF.tot_watveg_beg , dim1d='cell') 
     693        SRF.tot_watsoil_beg = lmdz.geo2point (SRF.tot_watsoil_beg, dim1d='cell') 
     694        SRF.snow_beg        = lmdz.geo2point (SRF.snow_beg       , dim1d='cell') 
     695        SRF.lakeres_beg     = lmdz.geo2point (SRF.lakeres_beg    , dim1d='cell') 
     696        SRF.tot_watveg_end  = lmdz.geo2point (SRF.tot_watveg_end , dim1d='cell') 
     697        SRF.tot_watsoil_end = lmdz.geo2point (SRF.tot_watsoil_end, dim1d='cell') 
     698        SRF.snow_end        = lmdz.geo2point (SRF.snow_end       , dim1d='cell') 
     699        SRF.lakeres_end     = lmdz.geo2point (SRF.lakeres_end    , dim1d='cell') 
    810700 
    811701 
    812702    # Stock dSoilHum dInterce dSWE dStream dFastR dSlowR dLake dPond dFlood 
    813703 
    814     SRF_wat_beg = SRF_tot_watveg_beg + SRF_tot_watsoil_beg + SRF_snow_beg 
    815     SRF_wat_end = SRF_tot_watveg_end + SRF_tot_watsoil_end + SRF_snow_end 
     704    SRF.wat_beg = SRF.tot_watveg_beg + SRF.tot_watsoil_beg + SRF.snow_beg 
     705    SRF.wat_end = SRF.tot_watveg_end + SRF.tot_watsoil_end + SRF.snow_end 
    816706 
    817707    echo ( '\n====================================================================================' ) 
     
    819709 
    820710    print ( ' 1/8', end='' ) ; sys.stdout.flush () 
    821     SRF_mas_watveg_beg   = SRF_stock_int ( SRF_tot_watveg_beg  ) 
     711    SRF.mass_watveg_beg   = SRF_stock_int ( SRF.tot_watveg_beg  ) 
    822712    print ( ' 2/8', end='' ) ; sys.stdout.flush () 
    823     SRF_mas_watsoil_beg  = SRF_stock_int ( SRF_tot_watsoil_beg ) 
     713    SRF.mass_watsoil_beg  = SRF_stock_int ( SRF.tot_watsoil_beg ) 
    824714    print ( ' 3/8', end='' ) ; sys.stdout.flush () 
    825     SRF_mas_snow_beg     = SRF_stock_int ( SRF_snow_beg        ) 
     715    SRF.mass_snow_beg     = SRF_stock_int ( SRF.snow_beg        ) 
    826716    print ( ' 4/8', end='' ) ; sys.stdout.flush () 
    827     SRF_mas_lake_beg     = ONE_stock_int ( SRF_lakeres_beg     ) 
     717    SRF.mass_lake_beg     = ONE_stock_int ( SRF.lakeres_beg     ) 
    828718    print ( ' 5/8', end='' ) ; sys.stdout.flush () 
    829719 
    830     SRF_mas_watveg_end   = SRF_stock_int ( SRF_tot_watveg_end  ) 
     720    SRF.mass_watveg_end   = SRF_stock_int ( SRF.tot_watveg_end  ) 
    831721    print ( ' 6/8', end='' ) ; sys.stdout.flush () 
    832     SRF_mas_watsoil_end  = SRF_stock_int ( SRF_tot_watsoil_end ) 
     722    SRF.mass_watsoil_end  = SRF_stock_int ( SRF.tot_watsoil_end ) 
    833723    print ( ' 7/8', end='' ) ; sys.stdout.flush () 
    834     SRF_mas_snow_end     = SRF_stock_int ( SRF_snow_end        ) 
     724    SRF.mass_snow_end     = SRF_stock_int ( SRF.snow_end        ) 
    835725    print ( ' 8/8', end='' ) ; sys.stdout.flush () 
    836     SRF_mas_lake_end     = ONE_stock_int ( SRF_lakeres_end     ) 
     726    SRF.mass_lake_end     = ONE_stock_int ( SRF.lakeres_end     ) 
    837727 
    838728    print (' -- ') ; sys.stdout.flush () 
    839729 
    840     dSRF_mas_watveg   = Sprec ( [SRF_mas_watveg_end , -SRF_mas_watveg_beg] ) 
    841     dSRF_mas_watsoil  = Sprec ( [SRF_mas_watsoil_end, -SRF_mas_watsoil_beg] ) 
    842     dSRF_mas_snow     = Sprec ( [SRF_mas_snow_end   , -SRF_mas_snow_beg] ) 
    843     dSRF_mas_lake     = Sprec ( [SRF_mas_lake_end   , -SRF_mas_lake_beg] ) 
     730    SRF.diff_mass_watveg   = Sprec ( [SRF.mass_watveg_end , -SRF.mass_watveg_beg] ) 
     731    SRF.diff_mass_watsoil  = Sprec ( [SRF.mass_watsoil_end, -SRF.mass_watsoil_beg] ) 
     732    SRF.diff_mass_snow     = Sprec ( [SRF.mass_snow_end   , -SRF.mass_snow_beg] ) 
     733    SRF.diff_mass_lake     = Sprec ( [SRF.mass_lake_end   , -SRF.mass_lake_beg] ) 
    844734 
    845735    echo ( '------------------------------------------------------------------------------------' ) 
    846736    echo ( f'\nSurface reservoirs  -- {Title} ') 
    847     echo ( f'SRF_mas_watveg_beg   = {SRF_mas_watveg_beg :12.6e} kg | SRF_mas_watveg_end   = {SRF_mas_watveg_end :12.6e} kg ' ) 
    848     echo ( f'SRF_mas_watsoil_beg  = {SRF_mas_watsoil_beg:12.6e} kg | SRF_mas_watsoil_end  = {SRF_mas_watsoil_end:12.6e} kg ' ) 
    849     echo ( f'SRF_mas_snow_beg     = {SRF_mas_snow_beg   :12.6e} kg | SRF_mas_snow_end     = {SRF_mas_snow_end   :12.6e} kg ' ) 
    850     echo ( f'SRF_mas_lake_beg     = {SRF_mas_lake_beg   :12.6e} kg | SRF_mas_lake_end     = {SRF_mas_lake_end   :12.6e} kg ' ) 
    851  
    852     prtFlux ( 'dMass (watveg) ', dSRF_mas_watveg , 'e' , True ) 
    853     prtFlux ( 'dMass (watsoil)', dSRF_mas_watsoil, 'e' , True ) 
    854     prtFlux ( 'dMass (snow)   ', dSRF_mas_snow   , 'e' , True ) 
    855     prtFlux ( 'dMass (lake)   ', dSRF_mas_lake   , 'e' , True ) 
    856  
    857     SRF_mas_wat_beg = Sprec ( [SRF_mas_watveg_beg , SRF_mas_watsoil_beg, SRF_mas_snow_beg] ) 
    858     SRF_mas_wat_end = Sprec ( [SRF_mas_watveg_end , SRF_mas_watsoil_end, SRF_mas_snow_end] ) 
    859  
    860     dSRF_mas_wat    = Sprec ( [+SRF_mas_watveg_end , +SRF_mas_watsoil_end, +SRF_mas_snow_end, 
    861                                -SRF_mas_watveg_beg , -SRF_mas_watsoil_beg, -SRF_mas_snow_beg]  ) 
     737    echo ( f'SRF.mass_watveg_beg   = {SRF.mass_watveg_beg :12.6e} kg | SRF.mass_watveg_end   = {SRF.mass_watveg_end :12.6e} kg ' ) 
     738    echo ( f'SRF.mass_watsoil_beg  = {SRF.mass_watsoil_beg:12.6e} kg | SRF.mass_watsoil_end  = {SRF.mass_watsoil_end:12.6e} kg ' ) 
     739    echo ( f'SRF.mass_snow_beg     = {SRF.mass_snow_beg   :12.6e} kg | SRF.mass_snow_end     = {SRF.mass_snow_end   :12.6e} kg ' ) 
     740    echo ( f'SRF.mass_lake_beg     = {SRF.mass_lake_beg   :12.6e} kg | SRF.mass_lake_end     = {SRF.mass_lake_end   :12.6e} kg ' ) 
     741 
     742    prtFlux ( 'dMass (watveg) ', SRF.diff_mass_watveg , 'e' , True ) 
     743    prtFlux ( 'dMass (watsoil)', SRF.diff_mass_watsoil, 'e' , True ) 
     744    prtFlux ( 'dMass (snow)   ', SRF.diff_mass_snow   , 'e' , True ) 
     745    prtFlux ( 'dMass (lake)   ', SRF.diff_mass_lake   , 'e' , True ) 
     746 
     747    SRF.mass_wat_beg = Sprec ( [SRF.mass_watveg_beg , SRF.mass_watsoil_beg, SRF.mass_snow_beg] ) 
     748    SRF.mass_wat_end = Sprec ( [SRF.mass_watveg_end , SRF.mass_watsoil_end, SRF.mass_snow_end] ) 
     749 
     750    SRF.diff_mass_wat    = Sprec ( [+SRF.mass_watveg_end , +SRF.mass_watsoil_end, +SRF.mass_snow_end, 
     751                               -SRF.mass_watveg_beg , -SRF.mass_watsoil_beg, -SRF.mass_snow_beg]  ) 
    862752 
    863753    echo ( '------------------------------------------------------------------------------------' ) 
    864754    echo ( f'Water content in surface  -- {Title} ' ) 
    865     echo ( f'SRF_mas_wat_beg   = {SRF_mas_wat_beg:12.6e} kg | SRF_mas_wat_end  = {SRF_mas_wat_end:12.6e} kg ') 
    866     prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', True ) 
     755    echo ( f'SRF.mass_wat_beg   = {SRF.mass_wat_beg:12.6e} kg | SRF.mass_wat_end  = {SRF.mass_wat_end:12.6e} kg ') 
     756    prtFlux ( 'dMass (water srf)', SRF.diff_mass_wat, 'e', True ) 
    867757 
    868758    echo ( '------------------------------------------------------------------------------------' ) 
    869     echo ( 'Water content in  ATM + SRF + RUN + LAKE' ) 
    870     echo ( 'mas_wat_beg = {:12.6e} kg | mas_wat_end = {:12.6e} kg '. 
    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 ) ) 
    873     prtFlux ( 'dMass (water atm+srf+run+lake)', dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat + dSRF_mas_lake, 'e', True) 
     759    echo ( 'Water content in  ATM + SECHIBA + RUN + LAKE' ) 
     760    echo ( 'mass_wat_beg = {:12.6e} kg | mass_wat_end = {:12.6e} kg '. 
     761            format (DYN.mass_wat_beg + ATM.mass_sno_beg + SRF.RUN_mass_wat_beg + SRF.mass_wat_beg + SRF.mass_lake_beg , 
     762                    DYN.mass_wat_end + ATM.mass_sno_end + SRF.RUN_mass_wat_end + SRF.mass_wat_end + SRF.mass_lake_end ) ) 
     763    prtFlux ( 'dMass (water atm+srf+run+lake)', DYN.diff_mass_wat + ATM.diff_mass_sno + SRF.RUN_diff_mass_wat + SRF.diff_mass_wat + SRF.diff_mass_lake, 'e', True) 
    874764 
    875765echo ( '\n====================================================================================' ) 
     
    878768if ATM_HIS == 'latlon' : 
    879769    echo ( ' latlon case' ) 
    880     ATM_wbilo_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_oce']), dim1d='cell' ) 
    881     ATM_wbilo_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_sic']), dim1d='cell' ) 
    882     ATM_wbilo_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_ter']), dim1d='cell' ) 
    883     ATM_wbilo_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_lic']), dim1d='cell' ) 
    884     ATM_runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1d='cell' ) 
    885     ATM_fqcalving   = lmdz.geo2point ( rprec (d_ATM_his ['fqcalving']), dim1d='cell' ) 
    886     ATM_fqfonte     = lmdz.geo2point ( rprec (d_ATM_his ['fqfonte']  ), dim1d='cell' ) 
    887     ATM_precip      = lmdz.geo2point ( rprec (d_ATM_his ['precip']   ), dim1d='cell' ) 
    888     ATM_snowf       = lmdz.geo2point ( rprec (d_ATM_his ['snow']     ), dim1d='cell' ) 
    889     ATM_evap        = lmdz.geo2point ( rprec (d_ATM_his ['evap']     ), dim1d='cell' ) 
    890     ATM_wevap_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_ter']), dim1d='cell' ) 
    891     ATM_wevap_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_oce']), dim1d='cell' ) 
    892     ATM_wevap_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_lic']), dim1d='cell' ) 
    893     ATM_wevap_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_sic']), dim1d='cell' ) 
    894     ATM_wrain_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_ter']), dim1d='cell' ) 
    895     ATM_wrain_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_oce']), dim1d='cell' ) 
    896     ATM_wrain_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_lic']), dim1d='cell' ) 
    897     ATM_wrain_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_sic']), dim1d='cell' ) 
    898     ATM_wsnow_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_ter']), dim1d='cell' ) 
    899     ATM_wsnow_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_oce']), dim1d='cell' ) 
    900     ATM_wsnow_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_lic']), dim1d='cell' ) 
    901     ATM_wsnow_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_sic']), dim1d='cell' ) 
    902     ATM_runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1d='cell' ) 
     770    ATM.wbilo_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_oce']), dim1d='cell' ) 
     771    ATM.wbilo_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_sic']), dim1d='cell' ) 
     772    ATM.wbilo_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_ter']), dim1d='cell' ) 
     773    ATM.wbilo_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_lic']), dim1d='cell' ) 
     774    ATM.runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1d='cell' ) 
     775    ATM.fqcalving   = lmdz.geo2point ( rprec (d_ATM_his ['fqcalving']), dim1d='cell' ) 
     776    ATM.fqfonte     = lmdz.geo2point ( rprec (d_ATM_his ['fqfonte']  ), dim1d='cell' ) 
     777    ATM.precip      = lmdz.geo2point ( rprec (d_ATM_his ['precip']   ), dim1d='cell' ) 
     778    ATM.snowf       = lmdz.geo2point ( rprec (d_ATM_his ['snow']     ), dim1d='cell' ) 
     779    ATM.evap        = lmdz.geo2point ( rprec (d_ATM_his ['evap']     ), dim1d='cell' ) 
     780    ATM.wevap_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_ter']), dim1d='cell' ) 
     781    ATM.wevap_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_oce']), dim1d='cell' ) 
     782    ATM.wevap_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_lic']), dim1d='cell' ) 
     783    ATM.wevap_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_sic']), dim1d='cell' ) 
     784    ATM.wrain_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_ter']), dim1d='cell' ) 
     785    ATM.wrain_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_oce']), dim1d='cell' ) 
     786    ATM.wrain_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_lic']), dim1d='cell' ) 
     787    ATM.wrain_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_sic']), dim1d='cell' ) 
     788    ATM.wsnow_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_ter']), dim1d='cell' ) 
     789    ATM.wsnow_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_oce']), dim1d='cell' ) 
     790    ATM.wsnow_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_lic']), dim1d='cell' ) 
     791    ATM.wsnow_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_sic']), dim1d='cell' ) 
     792    ATM.runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1d='cell' ) 
    903793    echo ( 'End of LATLON case') 
    904794 
    905795if ATM_HIS == 'ico' : 
    906796    echo (' ico case') 
    907     ATM_wbilo_oce   = rprec (d_ATM_his ['wbilo_oce']) 
    908     ATM_wbilo_sic   = rprec (d_ATM_his ['wbilo_sic']) 
    909     ATM_wbilo_ter   = rprec (d_ATM_his ['wbilo_ter']) 
    910     ATM_wbilo_lic   = rprec (d_ATM_his ['wbilo_lic']) 
    911     ATM_runofflic   = rprec (d_ATM_his ['runofflic']) 
    912     ATM_fqcalving   = rprec (d_ATM_his ['fqcalving']) 
    913     ATM_fqfonte     = rprec (d_ATM_his ['fqfonte']  ) 
    914     ATM_precip      = rprec (d_ATM_his ['precip']   ) 
    915     ATM_snowf       = rprec (d_ATM_his ['snow']     ) 
    916     ATM_evap        = rprec (d_ATM_his ['evap']     ) 
    917     ATM_wevap_ter   = rprec (d_ATM_his ['wevap_ter']) 
    918     ATM_wevap_oce   = rprec (d_ATM_his ['wevap_oce']) 
    919     ATM_wevap_lic   = rprec (d_ATM_his ['wevap_lic']) 
    920     ATM_wevap_sic   = rprec (d_ATM_his ['wevap_sic']) 
    921     ATM_runofflic   = rprec (d_ATM_his ['runofflic']) 
    922     ATM_wevap_ter   = rprec (d_ATM_his ['wevap_ter']) 
    923     ATM_wevap_oce   = rprec (d_ATM_his ['wevap_oce']) 
    924     ATM_wevap_lic   = rprec (d_ATM_his ['wevap_lic']) 
    925     ATM_wevap_sic   = rprec (d_ATM_his ['wevap_sic']) 
    926     ATM_wrain_ter   = rprec (d_ATM_his ['wrain_ter']) 
    927     ATM_wrain_oce   = rprec (d_ATM_his ['wrain_oce']) 
    928     ATM_wrain_lic   = rprec (d_ATM_his ['wrain_lic']) 
    929     ATM_wrain_sic   = rprec (d_ATM_his ['wrain_sic']) 
    930     ATM_wsnow_ter   = rprec (d_ATM_his ['wsnow_ter']) 
    931     ATM_wsnow_oce   = rprec (d_ATM_his ['wsnow_oce']) 
    932     ATM_wsnow_lic   = rprec (d_ATM_his ['wsnow_lic']) 
    933     ATM_wsnow_sic   = rprec (d_ATM_his ['wsnow_sic']) 
     797    ATM.wbilo_oce   = rprec (d_ATM_his ['wbilo_oce']) 
     798    ATM.wbilo_sic   = rprec (d_ATM_his ['wbilo_sic']) 
     799    ATM.wbilo_ter   = rprec (d_ATM_his ['wbilo_ter']) 
     800    ATM.wbilo_lic   = rprec (d_ATM_his ['wbilo_lic']) 
     801    ATM.runofflic   = rprec (d_ATM_his ['runofflic']) 
     802    ATM.fqcalving   = rprec (d_ATM_his ['fqcalving']) 
     803    ATM.fqfonte     = rprec (d_ATM_his ['fqfonte']  ) 
     804    ATM.precip      = rprec (d_ATM_his ['precip']   ) 
     805    ATM.snowf       = rprec (d_ATM_his ['snow']     ) 
     806    ATM.evap        = rprec (d_ATM_his ['evap']     ) 
     807    ATM.wevap_ter   = rprec (d_ATM_his ['wevap_ter']) 
     808    ATM.wevap_oce   = rprec (d_ATM_his ['wevap_oce']) 
     809    ATM.wevap_lic   = rprec (d_ATM_his ['wevap_lic']) 
     810    ATM.wevap_sic   = rprec (d_ATM_his ['wevap_sic']) 
     811    ATM.runofflic   = rprec (d_ATM_his ['runofflic']) 
     812    ATM.wevap_ter   = rprec (d_ATM_his ['wevap_ter']) 
     813    ATM.wevap_oce   = rprec (d_ATM_his ['wevap_oce']) 
     814    ATM.wevap_lic   = rprec (d_ATM_his ['wevap_lic']) 
     815    ATM.wevap_sic   = rprec (d_ATM_his ['wevap_sic']) 
     816    ATM.wrain_ter   = rprec (d_ATM_his ['wrain_ter']) 
     817    ATM.wrain_oce   = rprec (d_ATM_his ['wrain_oce']) 
     818    ATM.wrain_lic   = rprec (d_ATM_his ['wrain_lic']) 
     819    ATM.wrain_sic   = rprec (d_ATM_his ['wrain_sic']) 
     820    ATM.wsnow_ter   = rprec (d_ATM_his ['wsnow_ter']) 
     821    ATM.wsnow_oce   = rprec (d_ATM_his ['wsnow_oce']) 
     822    ATM.wsnow_lic   = rprec (d_ATM_his ['wsnow_lic']) 
     823    ATM.wsnow_sic   = rprec (d_ATM_his ['wsnow_sic']) 
    934824    echo ( 'End of ico case ') 
    935825 
    936826echo ( 'ATM wprecip_oce' ) 
    937 ATM_wprecip_oce = ATM_wrain_oce + ATM_wsnow_oce 
    938 ATM_wprecip_ter = ATM_wrain_ter + ATM_wsnow_ter 
    939 ATM_wprecip_sic = ATM_wrain_sic + ATM_wsnow_sic 
    940 ATM_wprecip_lic = ATM_wrain_lic + ATM_wsnow_lic 
    941  
    942 ATM_wbilo       = ATM_wbilo_oce   + ATM_wbilo_sic   + ATM_wbilo_ter   + ATM_wbilo_lic 
    943 ATM_wevap       = ATM_wevap_oce   + ATM_wevap_sic   + ATM_wevap_ter   + ATM_wevap_lic 
    944 ATM_wprecip     = ATM_wprecip_oce + ATM_wprecip_sic + ATM_wprecip_ter + ATM_wprecip_lic 
    945 ATM_wsnow       = ATM_wsnow_oce   + ATM_wsnow_sic   + ATM_wsnow_ter   + ATM_wsnow_lic 
    946 ATM_wrain       = ATM_wrain_oce   + ATM_wrain_sic   + ATM_wrain_ter   + ATM_wrain_lic 
    947 ATM_wemp        = ATM_wevap - ATM_wprecip 
    948 ATM_emp         = ATM_evap  - ATM_precip 
    949  
    950 ATM_wprecip_sea = ATM_wprecip_oce + ATM_wprecip_sic 
    951 ATM_wsnow_sea   = ATM_wsnow_oce   + ATM_wsnow_sic 
    952 ATM_wrain_sea   = ATM_wrain_oce   + ATM_wrain_sic 
    953 ATM_wbilo_sea   = ATM_wbilo_oce   + ATM_wbilo_sic 
    954 ATM_wevap_sea   = ATM_wevap_sic   + ATM_wevap_oce 
    955  
    956 ATM_wemp_ter    = ATM_wevap_ter - ATM_wprecip_ter 
    957 ATM_wemp_oce    = ATM_wevap_oce - ATM_wprecip_oce 
    958 ATM_wemp_sic    = ATM_wevap_sic - ATM_wprecip_sic 
    959 ATM_wemp_lic    = ATM_wevap_lic - ATM_wprecip_lic 
    960 ATM_wemp_sea    = ATM_wevap_sic - ATM_wprecip_oce 
    961  
    962 if SRF : 
     827ATM.wprecip_oce = ATM.wrain_oce + ATM.wsnow_oce 
     828ATM.wprecip_ter = ATM.wrain_ter + ATM.wsnow_ter 
     829ATM.wprecip_sic = ATM.wrain_sic + ATM.wsnow_sic 
     830ATM.wprecip_lic = ATM.wrain_lic + ATM.wsnow_lic 
     831 
     832ATM.wbilo       = ATM.wbilo_oce   + ATM.wbilo_sic   + ATM.wbilo_ter   + ATM.wbilo_lic 
     833ATM.wevap       = ATM.wevap_oce   + ATM.wevap_sic   + ATM.wevap_ter   + ATM.wevap_lic 
     834ATM.wprecip     = ATM.wprecip_oce + ATM.wprecip_sic + ATM.wprecip_ter + ATM.wprecip_lic 
     835ATM.wsnow       = ATM.wsnow_oce   + ATM.wsnow_sic   + ATM.wsnow_ter   + ATM.wsnow_lic 
     836ATM.wrain       = ATM.wrain_oce   + ATM.wrain_sic   + ATM.wrain_ter   + ATM.wrain_lic 
     837ATM.wemp        = ATM.wevap - ATM.wprecip 
     838ATM.emp         = ATM.evap  - ATM.precip 
     839 
     840ATM.wprecip_sea = ATM.wprecip_oce + ATM.wprecip_sic 
     841ATM.wsnow_sea   = ATM.wsnow_oce   + ATM.wsnow_sic 
     842ATM.wrain_sea   = ATM.wrain_oce   + ATM.wrain_sic 
     843ATM.wbilo_sea   = ATM.wbilo_oce   + ATM.wbilo_sic 
     844ATM.wevap_sea   = ATM.wevap_sic   + ATM.wevap_oce 
     845 
     846ATM.wemp_ter    = ATM.wevap_ter - ATM.wprecip_ter 
     847ATM.wemp_oce    = ATM.wevap_oce - ATM.wprecip_oce 
     848ATM.wemp_sic    = ATM.wevap_sic - ATM.wprecip_sic 
     849ATM.wemp_lic    = ATM.wevap_lic - ATM.wprecip_lic 
     850ATM.wemp_sea    = ATM.wevap_sic - ATM.wprecip_oce 
     851 
     852if SECHIBA : 
    963853    if RUN_HIS == 'latlon' : 
    964854        echo ( 'RUN costalflow Grille LATLON' ) 
    965855        if TestInterp : 
    966856            echo ( 'RUN runoff TestInterp' ) 
    967             RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff_contfrac_interp']  )   , dim1d='cell' ) 
    968             RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage_contfrac_interp'])   , dim1d='cell' ) 
     857            SRF.RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff_contfrac_interp']  )   , dim1d='cell' ) 
     858            SRF.RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage_contfrac_interp'])   , dim1d='cell' ) 
    969859        else : 
    970860            echo ( 'RUN runoff' ) 
    971             RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff']         ), dim1d='cell' ) 
    972             RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage']       ), dim1d='cell' ) 
    973  
    974         RUN_coastalflow     = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow']    ), dim1d='cell' ) 
    975         RUN_riverflow       = lmdz.geo2point ( rprec (d_RUN_his ['riverflow']      ), dim1d='cell' ) 
    976         RUN_riversret       = lmdz.geo2point ( rprec (d_RUN_his ['riversret']      ), dim1d='cell' ) 
    977         RUN_coastalflow_cpl = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow_cpl']), dim1d='cell' ) 
    978         RUN_riverflow_cpl   = lmdz.geo2point ( rprec (d_RUN_his ['riverflow_cpl']  ), dim1d='cell' ) 
     861            SRF.RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff']         ), dim1d='cell' ) 
     862            SRF.RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage']       ), dim1d='cell' ) 
     863 
     864        SRF.RUN_coastalflow     = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow']    ), dim1d='cell' ) 
     865        SRF.RUN_riverflow       = lmdz.geo2point ( rprec (d_RUN_his ['riverflow']      ), dim1d='cell' ) 
     866        SRF.RUN_riversret       = lmdz.geo2point ( rprec (d_RUN_his ['riversret']      ), dim1d='cell' ) 
     867        SRF.RUN_coastalflow_cpl = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow_cpl']), dim1d='cell' ) 
     868        SRF.RUN_riverflow_cpl   = lmdz.geo2point ( rprec (d_RUN_his ['riverflow_cpl']  ), dim1d='cell' ) 
    979869 
    980870    if RUN_HIS == 'ico' : 
    981871        echo ( 'RUN costalflow Grille ICO' ) 
    982         RUN_coastalflow =  rprec (d_RUN_his ['coastalflow']) 
    983         RUN_riverflow   =  rprec (d_RUN_his ['riverflow']  ) 
    984         RUN_runoff      =  rprec (d_RUN_his ['runoff']     ) 
    985         RUN_drainage    =  rprec (d_RUN_his ['drainage']   ) 
    986         RUN_riversret   =  rprec (d_RUN_his ['riversret']  ) 
    987  
    988         RUN_coastalflow_cpl = rprec (d_RUN_his ['coastalflow_cpl']) 
    989         RUN_riverflow_cpl   = rprec (d_RUN_his ['riverflow_cpl']  ) 
     872        SRF.RUN_coastalflow =  rprec (d_RUN_his ['coastalflow']) 
     873        SRF.RUN_riverflow   =  rprec (d_RUN_his ['riverflow']  ) 
     874        SRF.RUN_runoff      =  rprec (d_RUN_his ['runoff']     ) 
     875        SRF.RUN_drainage    =  rprec (d_RUN_his ['drainage']   ) 
     876        SRF.RUN_riversret   =  rprec (d_RUN_his ['riversret']  ) 
     877 
     878        SRF.RUN_coastalflow_cpl = rprec (d_RUN_his ['coastalflow_cpl']) 
     879        SRF.RUN_riverflow_cpl   = rprec (d_RUN_his ['riverflow_cpl']  ) 
    990880 
    991881    Step = 0 
     
    993883    if SRF_HIS == 'latlon' : 
    994884        if TestInterp : 
    995             echo ( 'SRF rain TestInterp' ) 
    996             SRF_rain     = lmdz.geo2point ( rprec (d_SRF_his ['rain_contfrac_interp'] ), dim1d='cell') 
    997             SRF_evap     = lmdz.geo2point ( rprec (d_SRF_his ['evap_contfrac_interp'] ), dim1d='cell') 
    998             SRF_snowf    = lmdz.geo2point ( rprec (d_SRF_his ['snow_contfrac_interp'] ), dim1d='cell') 
    999             SRF_subli    = lmdz.geo2point ( rprec (d_SRF_his ['subli_contfrac_interp']), dim1d='cell') 
    1000             SRF_transpir = lmdz.geo2point ( rprec (d_SRF_his ['transpir_contfrac_interp']).sum(dim='veget'), dim1d='cell' ) 
    1001             #SRF_rain.attrs.update     ( d_SRF_his ['rain_contfrac_interp'].attrs ) 
    1002             #SRF_evap.attrs.update     ( d_SRF_his ['evap_contfrac_interp'].attrs ) 
    1003             #SRF_snowf.attrs.update    ( d_SRF_his ['snow_contfrac_interp'].attrs ) 
    1004             #SRF_subli.attrs.update    ( d_SRF_his ['subli_contfrac_interp'].attrs ) 
    1005             #SRF_transpir.attrs.update ( d_SRF_his ['transpir_contfrac_interp'].attrs ) 
     885            echo ( 'SECHIBA rain TestInterp' ) 
     886            SRF.rain     = lmdz.geo2point ( rprec (d_SRF_his ['rain_contfrac_interp'] ), dim1d='cell') 
     887            SRF.evap     = lmdz.geo2point ( rprec (d_SRF_his ['evap_contfrac_interp'] ), dim1d='cell') 
     888            SRF.snowf    = lmdz.geo2point ( rprec (d_SRF_his ['snow_contfrac_interp'] ), dim1d='cell') 
     889            SRF.subli    = lmdz.geo2point ( rprec (d_SRF_his ['subli_contfrac_interp']), dim1d='cell') 
     890            SRF.transpir = lmdz.geo2point ( rprec (d_SRF_his ['transpir_contfrac_interp']).sum(dim='veget'), dim1d='cell' ) 
     891            #SRF.rain.attrs.update     ( d_SRF_his ['rain_contfrac_interp'].attrs ) 
     892            #SRF.evap.attrs.update     ( d_SRF_his ['evap_contfrac_interp'].attrs ) 
     893            #SRF.snowf.attrs.update    ( d_SRF_his ['snow_contfrac_interp'].attrs ) 
     894            #SRF.subli.attrs.update    ( d_SRF_his ['subli_contfrac_interp'].attrs ) 
     895            #SRF.transpir.attrs.update ( d_SRF_his ['transpir_contfrac_interp'].attrs ) 
    1006896        else : 
    1007             echo ( 'SRF rain' ) 
    1008             SRF_rain     = lmdz.geo2point ( rprec (d_SRF_his ['rain'] ) , dim1d='cell') 
    1009             SRF_evap     = lmdz.geo2point ( rprec (d_SRF_his ['evap'] ) , dim1d='cell') 
    1010             SRF_snowf    = lmdz.geo2point ( rprec (d_SRF_his ['snowf']) , dim1d='cell') 
    1011             SRF_subli    = lmdz.geo2point ( rprec (d_SRF_his ['subli']) , dim1d='cell') 
    1012             SRF_transpir = lmdz.geo2point ( rprec (d_SRF_his ['transpir']).sum(dim='veget'), dim1d='cell' ) 
     897            echo ( 'SECHIBA rain' ) 
     898            SRF.rain     = lmdz.geo2point ( rprec (d_SRF_his ['rain'] ) , dim1d='cell') 
     899            SRF.evap     = lmdz.geo2point ( rprec (d_SRF_his ['evap'] ) , dim1d='cell') 
     900            SRF.snowf    = lmdz.geo2point ( rprec (d_SRF_his ['snowf']) , dim1d='cell') 
     901            SRF.subli    = lmdz.geo2point ( rprec (d_SRF_his ['subli']) , dim1d='cell') 
     902            SRF.transpir = lmdz.geo2point ( rprec (d_SRF_his ['transpir']).sum(dim='veget'), dim1d='cell' ) 
    1013903 
    1014904    if SRF_HIS == 'ico' : 
    1015         echo ( 'SRF rain') 
    1016         SRF_rain     = rprec (d_SRF_his ['rain'] ) 
    1017         SRF_evap     = rprec (d_SRF_his ['evap'] ) 
    1018         SRF_snowf    = rprec (d_SRF_his ['snowf']) 
    1019         SRF_subli    = rprec (d_SRF_his ['subli']) 
    1020         SRF_transpir = rprec (d_SRF_his ['transpir']).sum(dim='veget') 
    1021  
    1022     echo ( 'SRF emp' ) 
    1023     SRF_transpir.attrs['units'] = d_SRF_his ['transpir'].attrs['units'] 
    1024     SRF_emp      = SRF_evap - SRF_rain - SRF_snowf ; SRF_emp.attrs['units'] = SRF_rain.attrs['units'] 
     905        echo ( 'SECHIBA rain') 
     906        SRF.rain     = rprec (d_SRF_his ['rain'] ) 
     907        SRF.evap     = rprec (d_SRF_his ['evap'] ) 
     908        SRF.snowf    = rprec (d_SRF_his ['snowf']) 
     909        SRF.subli    = rprec (d_SRF_his ['subli']) 
     910        SRF.transpir = rprec (d_SRF_his ['transpir']).sum(dim='veget') 
     911 
     912    echo ( 'SECHIBA emp' ) 
     913    SRF.transpir.attrs['units'] = d_SRF_his ['transpir'].attrs['units'] 
     914    SRF.emp      = SRF.evap - SRF.rain - SRF.snowf ; SRF.emp.attrs['units'] = SRF.rain.attrs['units'] 
    1025915 
    1026916    ## Correcting units of SECHIBA variables 
    1027     def mmd2si ( pvar ) : 
     917    def mmd2SI ( Var ) : 
    1028918        '''Change unit from mm/d or m^3/s to kg/s if needed''' 
    1029         if 'units' in pvar.attrs : 
    1030             if pvar.attrs['units'] in ['m^3/s', 'm3/s', 'm3.s-1',] : 
    1031                 pvar.values = pvar.values  * ATM_RHO                    ;  pvar.attrs['units'] = 'kg/s' 
    1032             if pvar.attrs['units'] == 'mm/d' : 
    1033                 pvar.values = pvar.values  * ATM_RHO * (1e-3/lmdz.RDAY) ;  pvar.attrs['units'] = 'kg/s' 
    1034             if pvar.attrs['units'] in ['m^3', 'm3', ] : 
    1035                 pvar.values = pvar.values  * ATM_RHO                    ;  pvar.attrs['units'] = 'kg' 
    1036  
    1037     for var in [ 'runoff', 'drainage', 'riversret', 'coastalflow', 'riverflow', 'coastalflow_cpl', 'riverflow_cpl' ] : 
    1038         zvar = locals()['RUN_' + var] 
    1039         mmd2si (zvar) 
    1040  
    1041     for var in ['evap', 'snowf', 'subli', 'transpir', 'rain', 'emp' ] : 
    1042         zvar = locals()['SRF_' + var] 
    1043         mmd2si (zvar) 
     919        if 'units' in Var.attrs :  
     920            if Var.attrs['units'] in ['m^3/s', 'm3/s', 'm3.s-1',] : 
     921                Var.values = Var.values  * ATM_RHO                 ;  Var.attrs['units'] = 'kg/s' 
     922            if Var.attrs['units'] == 'mm/d' : 
     923                Var.values = Var.values  * ATM_RHO * (1e-3/86400.) ;  Var.attrs['units'] = 'kg/s' 
     924            if Var.attrs['units'] in ['m^3', 'm3', ] : 
     925                Var.values = Var.values  * ATM_RHO                 ;  Var.attrs['units'] = 'kg' 
     926        return Var 
     927     
     928    SRF.RUN_coastalflow     = mmd2SI ( SRF.RUN_coastalflow ) 
     929    SRF.RUN_coastalflow_cpl = mmd2SI ( SRF.RUN_coastalflow_cpl ) 
     930    SRF.RUN_drainage        = mmd2SI ( SRF.RUN_drainage ) 
     931    SRF.RUN_riverflow       = mmd2SI ( SRF.RUN_riverflow ) 
     932    SRF.RUN_riverflow_cpl   = mmd2SI ( SRF.RUN_riverflow_cpl ) 
     933    SRF.RUN_riversret       = mmd2SI ( SRF.RUN_riversret ) 
     934    SRF.RUN_runoff          = mmd2SI ( SRF.RUN_runoff ) 
     935             
     936    SRF.evap     = mmd2SI ( SRF.evap ) 
     937    SRF.snowf    = mmd2SI ( SRF.snowf ) 
     938    SRF.subli    = mmd2SI ( SRF.subli ) 
     939    SRF.transpir = mmd2SI ( SRF.transpir ) 
     940    SRF.rain     = mmd2SI ( SRF.rain ) 
     941    SRF.emp      = mmd2SI ( SRF.emp ) 
    1044942 
    1045943    echo ( 'RUN input' ) 
    1046     RUN_input  = RUN_runoff      + RUN_drainage 
    1047     RUN_output = RUN_coastalflow + RUN_riverflow 
     944    SRF.RUN_input  = SRF.RUN_runoff      + SRF.RUN_drainage 
     945    SRF.RUN_output = SRF.RUN_coastalflow + SRF.RUN_riverflow 
    1048946 
    1049947echo ( 'ATM flw_wbilo' ) 
    1050 ATM_flx_wbilo       = ATM_flux_int ( ATM_wbilo      ) 
    1051 ATM_flx_wevap       = ATM_flux_int ( ATM_wevap      ) 
    1052 ATM_flx_wprecip     = ATM_flux_int ( ATM_wprecip    ) 
    1053 ATM_flx_wsnow       = ATM_flux_int ( ATM_wsnow      ) 
    1054 ATM_flx_wrain       = ATM_flux_int ( ATM_wrain      ) 
    1055 ATM_flx_wemp        = ATM_flux_int ( ATM_wemp       ) 
    1056  
    1057 ATM_flx_wbilo_lic   = ATM_flux_int ( ATM_wbilo_lic  ) 
    1058 ATM_flx_wbilo_oce   = ATM_flux_int ( ATM_wbilo_oce  ) 
    1059 ATM_flx_wbilo_sea   = ATM_flux_int ( ATM_wbilo_sea  ) 
    1060 ATM_flx_wbilo_sic   = ATM_flux_int ( ATM_wbilo_sic  ) 
    1061 ATM_flx_wbilo_ter   = ATM_flux_int ( ATM_wbilo_ter  ) 
     948ATM.flx_wbilo       = ATM_flux_int ( ATM.wbilo      ) 
     949ATM.flx_wevap       = ATM_flux_int ( ATM.wevap      ) 
     950ATM.flx_wprecip     = ATM_flux_int ( ATM.wprecip    ) 
     951ATM.flx_wsnow       = ATM_flux_int ( ATM.wsnow      ) 
     952ATM.flx_wrain       = ATM_flux_int ( ATM.wrain      ) 
     953ATM.flx_wemp        = ATM_flux_int ( ATM.wemp       ) 
     954 
     955ATM.flx_wbilo_lic   = ATM_flux_int ( ATM.wbilo_lic  ) 
     956ATM.flx_wbilo_oce   = ATM_flux_int ( ATM.wbilo_oce  ) 
     957ATM.flx_wbilo_sea   = ATM_flux_int ( ATM.wbilo_sea  ) 
     958ATM.flx_wbilo_sic   = ATM_flux_int ( ATM.wbilo_sic  ) 
     959ATM.flx_wbilo_ter   = ATM_flux_int ( ATM.wbilo_ter  ) 
    1062960# Type d'integration a verifier 
    1063 ATM_flx_calving     = ATM_flux_int ( ATM_fqcalving  ) 
    1064 ATM_flx_fqfonte     = ATM_flux_int ( ATM_fqfonte    ) 
    1065  
    1066 LIC_flx_calving     = LIC_flux_int ( ATM_fqcalving  ) 
    1067 LIC_flx_fqfonte     = LIC_flux_int ( ATM_fqfonte    ) 
     961ATM.flx_calving     = ATM_flux_int ( ATM.fqcalving  ) 
     962ATM.flx_fqfonte     = ATM_flux_int ( ATM.fqfonte    ) 
     963 
     964ATM.LIC_flx_calving     = LIC_flux_int ( ATM.fqcalving  ) 
     965ATM.LIC_flx_fqfonte     = LIC_flux_int ( ATM.fqfonte    ) 
    1068966 
    1069967echo ( 'ATM flx precip' ) 
    1070 ATM_flx_precip      = ATM_flux_int ( ATM_precip     ) 
    1071 ATM_flx_snowf       = ATM_flux_int ( ATM_snowf      ) 
    1072 ATM_flx_evap        = ATM_flux_int ( ATM_evap       ) 
    1073 ATM_flx_runlic      = ATM_flux_int ( ATM_runofflic  ) 
    1074  
    1075 LIC_flx_precip      = LIC_flux_int ( ATM_precip     ) 
    1076 LIC_flx_snowf       = LIC_flux_int ( ATM_snowf      ) 
    1077 LIC_flx_evap        = LIC_flux_int ( ATM_evap       ) 
    1078 LIC_flx_runlic      = LIC_flux_int ( ATM_runofflic  ) 
     968ATM.flx_precip      = ATM_flux_int ( ATM.precip     ) 
     969ATM.flx_snowf       = ATM_flux_int ( ATM.snowf      ) 
     970ATM.flx_evap        = ATM_flux_int ( ATM.evap       ) 
     971ATM.flx_runlic      = ATM_flux_int ( ATM.runofflic  ) 
     972 
     973ATM.LIC_flx_precip      = LIC_flux_int ( ATM.precip     ) 
     974ATM.LIC_flx_snowf       = LIC_flux_int ( ATM.snowf      ) 
     975ATM.LIC_flx_evap        = LIC_flux_int ( ATM.evap       ) 
     976ATM.LIC_flx_runlic      = LIC_flux_int ( ATM.runofflic  ) 
    1079977 
    1080978echo ( 'ATM flx_wrain_ter' ) 
    1081 ATM_flx_wrain_ter    = ATM_flux_int ( ATM_wrain_ter ) 
    1082 ATM_flx_wrain_oce    = ATM_flux_int ( ATM_wrain_oce ) 
    1083 ATM_flx_wrain_lic    = ATM_flux_int ( ATM_wrain_lic ) 
    1084 ATM_flx_wrain_sic    = ATM_flux_int ( ATM_wrain_sic ) 
    1085 ATM_flx_wrain_sea    = ATM_flux_int ( ATM_wrain_sea ) 
    1086  
    1087 ATM_flx_wsnow_ter    = ATM_flux_int ( ATM_wsnow_ter ) 
    1088 ATM_flx_wsnow_oce    = ATM_flux_int ( ATM_wsnow_oce ) 
    1089 ATM_flx_wsnow_lic    = ATM_flux_int ( ATM_wsnow_lic ) 
    1090 ATM_flx_wsnow_sic    = ATM_flux_int ( ATM_wsnow_sic ) 
    1091 ATM_flx_wsnow_sea    = ATM_flux_int ( ATM_wsnow_sea ) 
     979ATM.flx_wrain_ter    = ATM_flux_int ( ATM.wrain_ter ) 
     980ATM.flx_wrain_oce    = ATM_flux_int ( ATM.wrain_oce ) 
     981ATM.flx_wrain_lic    = ATM_flux_int ( ATM.wrain_lic ) 
     982ATM.flx_wrain_sic    = ATM_flux_int ( ATM.wrain_sic ) 
     983ATM.flx_wrain_sea    = ATM_flux_int ( ATM.wrain_sea ) 
     984 
     985ATM.flx_wsnow_ter    = ATM_flux_int ( ATM.wsnow_ter ) 
     986ATM.flx_wsnow_oce    = ATM_flux_int ( ATM.wsnow_oce ) 
     987ATM.flx_wsnow_lic    = ATM_flux_int ( ATM.wsnow_lic ) 
     988ATM.flx_wsnow_sic    = ATM_flux_int ( ATM.wsnow_sic ) 
     989ATM.flx_wsnow_sea    = ATM_flux_int ( ATM.wsnow_sea ) 
    1092990 
    1093991echo ( 'ATM flx_evap_ter' ) 
    1094 ATM_flx_wevap_ter    = ATM_flux_int ( ATM_wevap_ter ) 
    1095 ATM_flx_wevap_oce    = ATM_flux_int ( ATM_wevap_oce ) 
    1096 ATM_flx_wevap_lic    = ATM_flux_int ( ATM_wevap_lic ) 
    1097 ATM_flx_wevap_sic    = ATM_flux_int ( ATM_wevap_sic ) 
    1098 ATM_flx_wevap_sea    = ATM_flux_int ( ATM_wevap_sea ) 
    1099 ATM_flx_wprecip_lic  = ATM_flux_int ( ATM_wprecip_lic ) 
    1100 ATM_flx_wprecip_oce  = ATM_flux_int ( ATM_wprecip_oce ) 
    1101 ATM_flx_wprecip_sic  = ATM_flux_int ( ATM_wprecip_sic ) 
    1102 ATM_flx_wprecip_ter  = ATM_flux_int ( ATM_wprecip_ter ) 
    1103 ATM_flx_wprecip_sea  = ATM_flux_int ( ATM_wprecip_sea ) 
    1104 ATM_flx_wemp_lic     = ATM_flux_int ( ATM_wemp_lic ) 
    1105 ATM_flx_wemp_oce     = ATM_flux_int ( ATM_wemp_oce ) 
    1106 ATM_flx_wemp_sic     = ATM_flux_int ( ATM_wemp_sic ) 
    1107 ATM_flx_wemp_ter     = ATM_flux_int ( ATM_wemp_ter ) 
    1108 ATM_flx_wemp_sea     = ATM_flux_int ( ATM_wemp_sea ) 
    1109  
    1110 ATM_flx_emp          = ATM_flux_int ( ATM_emp ) 
    1111  
    1112 if SRF : 
     992ATM.flx_wevap_ter    = ATM_flux_int ( ATM.wevap_ter ) 
     993ATM.flx_wevap_oce    = ATM_flux_int ( ATM.wevap_oce ) 
     994ATM.flx_wevap_lic    = ATM_flux_int ( ATM.wevap_lic ) 
     995ATM.flx_wevap_sic    = ATM_flux_int ( ATM.wevap_sic ) 
     996ATM.flx_wevap_sea    = ATM_flux_int ( ATM.wevap_sea ) 
     997ATM.flx_wprecip_lic  = ATM_flux_int ( ATM.wprecip_lic ) 
     998ATM.flx_wprecip_oce  = ATM_flux_int ( ATM.wprecip_oce ) 
     999ATM.flx_wprecip_sic  = ATM_flux_int ( ATM.wprecip_sic ) 
     1000ATM.flx_wprecip_ter  = ATM_flux_int ( ATM.wprecip_ter ) 
     1001ATM.flx_wprecip_sea  = ATM_flux_int ( ATM.wprecip_sea ) 
     1002ATM.flx_wemp_lic     = ATM_flux_int ( ATM.wemp_lic ) 
     1003ATM.flx_wemp_oce     = ATM_flux_int ( ATM.wemp_oce ) 
     1004ATM.flx_wemp_sic     = ATM_flux_int ( ATM.wemp_sic ) 
     1005ATM.flx_wemp_ter     = ATM_flux_int ( ATM.wemp_ter ) 
     1006ATM.flx_wemp_sea     = ATM_flux_int ( ATM.wemp_sea ) 
     1007 
     1008ATM.flx_emp          = ATM_flux_int ( ATM.emp ) 
     1009 
     1010if SECHIBA : 
    11131011    echo ( 'RUN flx_coastal' ) 
    1114     RUN_flx_coastal     = ONE_flux_int ( RUN_coastalflow) 
     1012    SRF.RUN_flx_coastal     = ONE_flux_int ( SRF.RUN_coastalflow) 
    11151013    echo ( 'RUN flx_river' ) 
    1116     RUN_flx_river       = ONE_flux_int ( RUN_riverflow  ) 
     1014    SRF.RUN_flx_river       = ONE_flux_int ( SRF.RUN_riverflow  ) 
    11171015    echo ( 'RUN flx_coastal_cpl' ) 
    1118     RUN_flx_coastal_cpl = ONE_flux_int ( RUN_coastalflow_cpl) 
     1016    SRF.RUN_flx_coastal_cpl = ONE_flux_int ( SRF.RUN_coastalflow_cpl) 
    11191017    echo ( 'RUN flx_river_cpl' ) 
    1120     RUN_flx_river_cpl   = ONE_flux_int ( RUN_riverflow_cpl  ) 
     1018    SRF.RUN_flx_river_cpl   = ONE_flux_int ( SRF.RUN_riverflow_cpl  ) 
    11211019    echo ( 'RUN flx_drainage' ) 
    1122     RUN_flx_drainage    = SRF_flux_int ( RUN_drainage   ) 
     1020    SRF.RUN_flx_drainage    = SRF_flux_int ( SRF.RUN_drainage   ) 
    11231021    echo ( 'RUN flx_riversset' ) 
    1124     RUN_flx_riversret   = SRF_flux_int ( RUN_riversret  ) 
     1022    SRF.RUN_flx_riversret   = SRF_flux_int ( SRF.RUN_riversret  ) 
    11251023    echo ( 'RUN flx_runoff' ) 
    1126     RUN_flx_runoff      = SRF_flux_int ( RUN_runoff     ) 
     1024    SRF.RUN_flx_runoff      = SRF_flux_int ( SRF.RUN_runoff     ) 
    11271025    echo ( 'RUN flx_input' ) 
    1128     RUN_flx_input       = SRF_flux_int ( RUN_input      ) 
     1026    SRF.RUN_flx_input       = SRF_flux_int ( SRF.RUN_input      ) 
    11291027    echo ( 'RUN flx_output' ) 
    1130     RUN_flx_output      = ONE_flux_int ( RUN_output     ) 
     1028    SRF.RUN_flx_output      = ONE_flux_int ( SRF.RUN_output     ) 
    11311029 
    11321030    echo ( 'RUN flx_bil' ) ; Step += 1 
    1133     #RUN_flx_bil    = RUN_flx_input   - RUN_flx_output 
    1134     #RUN_flx_rivcoa = RUN_flx_coastal + RUN_flx_river 
    1135  
    1136     RUN_flx_bil    = ONE_flux_int ( RUN_input       - RUN_output) 
    1137     RUN_flx_rivcoa = ONE_flux_int ( RUN_coastalflow + RUN_riverflow) 
    1138  
    1139 prtFlux ('wbilo_oce            ', ATM_flx_wbilo_oce     , 'f' ) 
    1140 prtFlux ('wbilo_sic            ', ATM_flx_wbilo_sic     , 'f' ) 
    1141 prtFlux ('wbilo_sic+oce        ', ATM_flx_wbilo_sea     , 'f' ) 
    1142 prtFlux ('wbilo_ter            ', ATM_flx_wbilo_ter     , 'f' ) 
    1143 prtFlux ('wbilo_lic            ', ATM_flx_wbilo_lic     , 'f' ) 
    1144 prtFlux ('Sum wbilo_*          ', ATM_flx_wbilo         , 'f', True) 
    1145 prtFlux ('E-P                  ', ATM_flx_emp           , 'f', True) 
    1146 prtFlux ('calving              ', ATM_flx_calving       , 'f' ) 
    1147 prtFlux ('fqfonte              ', ATM_flx_fqfonte       , 'f' ) 
    1148 prtFlux ('precip               ', ATM_flx_precip        , 'f' ) 
    1149 prtFlux ('snowf                ', ATM_flx_snowf         , 'f' ) 
    1150 prtFlux ('evap                 ', ATM_flx_evap          , 'f' ) 
    1151 prtFlux ('runoff lic           ', ATM_flx_runlic        , 'f' ) 
    1152  
    1153 prtFlux ('ATM_flx_wevap*       ', ATM_flx_wevap         , 'f' ) 
    1154 prtFlux ('ATM_flx_wrain*       ', ATM_flx_wrain         , 'f' ) 
    1155 prtFlux ('ATM_flx_wsnow*       ', ATM_flx_wsnow         , 'f' ) 
    1156 prtFlux ('ATM_flx_wprecip*     ', ATM_flx_wprecip       , 'f' ) 
    1157 prtFlux ('ATM_flx_wemp*        ', ATM_flx_wemp          , 'f', True ) 
     1031    #SRF.RUN_flx_bil    = SRF.RUN_flx_input   - SRF.RUN_flx_output 
     1032    #SRF.RUN_flx_rivcoa = SRF.RUN_flx_coastal + SRF.RUN_flx_river 
     1033 
     1034    SRF.RUN_flx_bil    = ONE_flux_int ( SRF.RUN_input       - SRF.RUN_output) 
     1035    SRF.RUN_flx_rivcoa = ONE_flux_int ( SRF.RUN_coastalflow + SRF.RUN_riverflow) 
     1036 
     1037prtFlux ('wbilo_oce            ', ATM.flx_wbilo_oce     , 'f' ) 
     1038prtFlux ('wbilo_sic            ', ATM.flx_wbilo_sic     , 'f' ) 
     1039prtFlux ('wbilo_sic+oce        ', ATM.flx_wbilo_sea     , 'f' ) 
     1040prtFlux ('wbilo_ter            ', ATM.flx_wbilo_ter     , 'f' ) 
     1041prtFlux ('wbilo_lic            ', ATM.flx_wbilo_lic     , 'f' ) 
     1042prtFlux ('Sum wbilo_*          ', ATM.flx_wbilo         , 'f', True) 
     1043prtFlux ('E-P                  ', ATM.flx_emp           , 'f', True) 
     1044prtFlux ('calving              ', ATM.flx_calving       , 'f' ) 
     1045prtFlux ('fqfonte              ', ATM.flx_fqfonte       , 'f' ) 
     1046prtFlux ('precip               ', ATM.flx_precip        , 'f' ) 
     1047prtFlux ('snowf                ', ATM.flx_snowf         , 'f' ) 
     1048prtFlux ('evap                 ', ATM.flx_evap          , 'f' ) 
     1049prtFlux ('runoff lic           ', ATM.flx_runlic        , 'f' ) 
     1050 
     1051prtFlux ('ATM.flx_wevap*       ', ATM.flx_wevap         , 'f' ) 
     1052prtFlux ('ATM.flx_wrain*       ', ATM.flx_wrain         , 'f' ) 
     1053prtFlux ('ATM.flx_wsnow*       ', ATM.flx_wsnow         , 'f' ) 
     1054prtFlux ('ATM.flx_wprecip*     ', ATM.flx_wprecip       , 'f' ) 
     1055prtFlux ('ATM.flx_wemp*        ', ATM.flx_wemp          , 'f', True ) 
    11581056 
    11591057echo ( 'Errors <field> vs. wbil_<field>' ) 
    1160 prtFlux ('ERROR evap           ', ATM_flx_wevap   - ATM_flx_evap  , 'e', True ) 
    1161 prtFlux ('ERROR precip         ', ATM_flx_wprecip - ATM_flx_precip, 'e', True ) 
    1162 prtFlux ('ERROR snow           ', ATM_flx_wsnow   - ATM_flx_snowf , 'e', True ) 
    1163 prtFlux ('ERROR emp            ', ATM_flx_wemp    - ATM_flx_emp   , 'e', True ) 
    1164  
    1165 if SRF : 
     1058prtFlux ('ERROR evap           ', ATM.flx_wevap   - ATM.flx_evap  , 'e', True ) 
     1059prtFlux ('ERROR precip         ', ATM.flx_wprecip - ATM.flx_precip, 'e', True ) 
     1060prtFlux ('ERROR snow           ', ATM.flx_wsnow   - ATM.flx_snowf , 'e', True ) 
     1061prtFlux ('ERROR emp            ', ATM.flx_wemp    - ATM.flx_emp   , 'e', True ) 
     1062 
     1063if SECHIBA : 
    11661064    echo ( '\n====================================================================================' ) 
    11671065    echo ( f'-- RUNOFF Fluxes  -- {Title} ' ) 
    1168     prtFlux ('coastalflow   ', RUN_flx_coastal    , 'f' ) 
    1169     prtFlux ('riverflow     ', RUN_flx_river      , 'f' ) 
    1170     prtFlux ('coastal_cpl   ', RUN_flx_coastal_cpl, 'f' ) 
    1171     prtFlux ('riverf_cpl    ', RUN_flx_river_cpl  , 'f' ) 
    1172     prtFlux ('river+coastal ', RUN_flx_rivcoa     , 'f' ) 
    1173     prtFlux ('drainage      ', RUN_flx_drainage   , 'f' ) 
    1174     prtFlux ('riversret     ', RUN_flx_riversret  , 'f' ) 
    1175     prtFlux ('runoff        ', RUN_flx_runoff     , 'f' ) 
    1176     prtFlux ('river in      ', RUN_flx_input      , 'f' ) 
    1177     prtFlux ('river out     ', RUN_flx_output     , 'f' ) 
    1178     prtFlux ('river bil     ', RUN_flx_bil        , 'f' ) 
    1179  
    1180 ATM_flx_budget = -ATM_flx_wbilo + ATM_flx_calving + ATM_flx_runlic #+# ATM_flx_fqfonte #+ RUN_flx_river 
     1066    prtFlux ('coastalflow   ', SRF.RUN_flx_coastal    , 'f' ) 
     1067    prtFlux ('riverflow     ', SRF.RUN_flx_river      , 'f' ) 
     1068    prtFlux ('coastal_cpl   ', SRF.RUN_flx_coastal_cpl, 'f' ) 
     1069    prtFlux ('riverf_cpl    ', SRF.RUN_flx_river_cpl  , 'f' ) 
     1070    prtFlux ('river+coastal ', SRF.RUN_flx_rivcoa     , 'f' ) 
     1071    prtFlux ('drainage      ', SRF.RUN_flx_drainage   , 'f' ) 
     1072    prtFlux ('riversret     ', SRF.RUN_flx_riversret  , 'f' ) 
     1073    prtFlux ('runoff        ', SRF.RUN_flx_runoff     , 'f' ) 
     1074    prtFlux ('river in      ', SRF.RUN_flx_input      , 'f' ) 
     1075    prtFlux ('river out     ', SRF.RUN_flx_output     , 'f' ) 
     1076    prtFlux ('river bil     ', SRF.RUN_flx_bil        , 'f' ) 
     1077 
     1078ATM.flx_budget = -ATM.flx_wbilo + ATM.flx_calving + ATM.flx_runlic #+# ATM.flx_fqfonte #+ SRF.RUN_flx_river 
    11811079 
    11821080 
    11831081echo ('') 
    1184 #echo ('  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 )) 
    1185  
    1186 #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 )) 
    1187  
    1188 ATM_flx_toSRF = -ATM_flx_wbilo_ter 
     1082#echo ('  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 )) 
     1083 
     1084#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 )) 
     1085 
     1086ATM.flx_toSRF = -ATM.flx_wbilo_ter 
    11891087 
    11901088echo (' ') 
    11911089echo ( '\n====================================================================================' ) 
    11921090echo ( f'--  Atmosphere  -- {Title} ' ) 
    1193 echo ( f'Mass begin = {DYN_mas_wat_beg:12.6e} kg | Mass end = {DYN_mas_wat_end:12.6e} kg' ) 
    1194 prtFlux ( 'dmass (atm)  = ', dDYN_mas_wat , 'e', True ) 
    1195 prtFlux ( 'Sum wbilo_*  = ', ATM_flx_wbilo, 'e', True ) 
    1196 prtFlux ( 'E-P          = ', ATM_flx_emp  , 'e', True ) 
     1091echo ( f'Mass begin = {DYN.mass_wat_beg:12.6e} kg | Mass end = {DYN.mass_wat_end:12.6e} kg' ) 
     1092prtFlux ( 'dmass (atm) ', DYN.diff_mass_wat , 'e', True ) 
     1093prtFlux ( 'Sum wbilo_* ', ATM.flx_wbilo, 'e', True ) 
     1094prtFlux ( 'E-P         ', ATM.flx_emp  , 'e', True ) 
    11971095echo ( ' ' ) 
    1198 prtFlux ( 'Water loss atm from wbil_*', ATM_flx_wbilo - dDYN_mas_wat, 'f', True ) 
    1199 echo ( 'Water loss atm = {:12.3e} (rel)  '.format ( (ATM_flx_wbilo - dDYN_mas_wat)/dDYN_mas_wat  ) ) 
     1096prtFlux ( 'Water loss atm from wbil_*', ATM.flx_wbilo - DYN.diff_mass_wat, 'f', True ) 
     1097echo ( 'Water loss atm = {:12.3e} (rel)  '.format ( (ATM.flx_wbilo - DYN.diff_mass_wat)/DYN.diff_mass_wat  ) ) 
    12001098 
    12011099echo (' ') 
    1202 prtFlux ( 'Water loss atm from E-P', ATM_flx_emp  - dDYN_mas_wat , 'f', True ) 
    1203 echo ( 'Water loss atm = {:12.3e} (rel)  '.format ( (ATM_flx_emp-dDYN_mas_wat)/dDYN_mas_wat  ) ) 
     1100prtFlux ( 'Water loss atm from E-P', ATM.flx_emp  - DYN.diff_mass_wat , 'f', True ) 
     1101echo ( 'Water loss atm = {:12.3e} (rel)  '.format ( (ATM.flx_emp-DYN.diff_mass_wat)/DYN.diff_mass_wat  ) ) 
    12041102echo (' ') 
    12051103 
    1206 ATM_error =  ATM_flx_emp  - dDYN_mas_wat 
     1104ATM.error =  ATM.flx_emp  - DYN.diff_mass_wat 
    12071105 
    12081106 
     
    12101108echo ( '\n====================================================================================' ) 
    12111109 
    1212 LIC_flx_budget1 = Sprec ( [-ATM_flx_wemp_lic  , -LIC_flx_calving , -LIC_flx_fqfonte] ) 
    1213 LIC_flx_budget2 = Sprec ( [-ATM_flx_wbilo_lic , -LIC_flx_calving , -LIC_flx_fqfonte] ) 
    1214 LIC_flx_budget3 = Sprec ( [-ATM_flx_wbilo_lic , -LIC_flx_runlic] ) 
    1215 LIC_flx_budget4 = Sprec ( [-ATM_flx_wemp_lic  , -LIC_flx_runlic] ) 
     1110ATM.LIC_flx_budget1 = Sprec ( [-ATM.flx_wemp_lic  , -ATM.LIC_flx_calving , -ATM.LIC_flx_fqfonte] ) 
     1111ATM.LIC_flx_budget2 = Sprec ( [-ATM.flx_wbilo_lic , -ATM.LIC_flx_calving , -ATM.LIC_flx_fqfonte] ) 
     1112ATM.LIC_flx_budget3 = Sprec ( [-ATM.flx_wbilo_lic , -ATM.LIC_flx_runlic] ) 
     1113ATM.LIC_flx_budget4 = Sprec ( [-ATM.flx_wemp_lic  , -ATM.LIC_flx_runlic] ) 
    12161114 
    12171115echo ( f'--  LIC  -- {Title} ' ) 
    1218 echo ( f'Mass total   begin = {LIC_mas_wat_beg    :12.6e} kg | Mass end = {LIC_mas_wat_end    :12.6e} kg' ) 
    1219 echo ( f'Mass snow    begin = {LIC_mas_sno_beg    :12.6e} kg | Mass end = {LIC_mas_sno_end    :12.6e} kg' ) 
    1220 echo ( f'Mass qs      begin = {LIC_mas_qs_beg     :12.6e} kg | Mass end = {LIC_mas_qs_end     :12.6e} kg' ) 
    1221 echo ( f'Mass runlic0 begin = {LIC_mas_runlic0_beg:12.6e} kg | Mass end = {LIC_mas_runlic0_end:12.6e} kg' ) 
    1222 prtFlux ( 'dmass (LIC sno)       ', dLIC_mas_sno          , 'f', True, width=45 ) 
    1223 prtFlux ( 'dmass (LIC qs)        ', dLIC_mas_qs           , 'e', True, width=45 ) 
    1224 prtFlux ( 'dmass (LIC wat)       ', dLIC_mas_wat          , 'f', True, width=45 ) 
    1225 prtFlux ( 'dmass (LIC runlic0)   ', dLIC_mas_runlic0      , 'e', True, width=45 ) 
    1226 prtFlux ( 'dmass (LIC total)     ', dLIC_mas_wat          , 'e', True, width=45 ) 
    1227 prtFlux ( 'LIC ATM_flx_wemp_lic  ', ATM_flx_wemp_lic      , 'f', True, width=45 ) 
    1228 prtFlux ( 'LIC LIC_flx_fqfonte   ', LIC_flx_fqfonte       , 'f', True, width=45 ) 
    1229 prtFlux ( 'LIC LIC_flx_calving   ', LIC_flx_calving       , 'f', True, width=45 ) 
    1230 prtFlux ( 'LIC LIC_flx_runofflic ', LIC_flx_runlic        , 'f', True, width=45 ) 
    1231 prtFlux ( 'LIC fqfonte + calving ', LIC_flx_calving+LIC_flx_fqfonte , 'f', True, width=45 ) 
    1232 prtFlux ( 'LIC fluxes 1 ( wemp_lic   - fqcalving - fqfonte)) ', LIC_flx_budget1              , 'f', True, width=45 ) 
    1233 prtFlux ( 'LIC fluxes 2 (-wbilo_lic - fqcalving - fqfonte)   ', LIC_flx_budget2              , 'f', True, width=45 ) 
    1234 prtFlux ( 'LIC fluxes 3 (-wbilo_lic - runofflic*frac_lic)    ', LIC_flx_budget3              , 'f', True, width=45 ) 
    1235 prtFlux ( 'LIC fluxes 3 ( wemp_lic  - runofflic*frac_lic)    ', LIC_flx_budget4              , 'f', True, width=45 ) 
    1236 prtFlux ( 'LIC error 1                                       ', LIC_flx_budget1-dLIC_mas_wat , 'e', True, width=45 ) 
    1237 prtFlux ( 'LIC error 2                                       ', LIC_flx_budget2-dLIC_mas_wat , 'e', True, width=45 ) 
    1238 prtFlux ( 'LIC error 3                                       ', LIC_flx_budget3-dLIC_mas_wat , 'e', True, width=45 ) 
    1239 echo ( 'LIC error (wevap - precip*frac_lic  - fqcalving - fqfonte)    = {:12.4e} (rel) '.format ( (LIC_flx_budget1-dLIC_mas_wat)/dLIC_mas_wat) ) 
    1240 echo ( 'LIC error (-wbilo_lic - fqcalving - fqfonte)                  = {:12.4e} (rel) '.format ( (LIC_flx_budget2-dLIC_mas_wat)/dLIC_mas_wat) ) 
    1241 echo ( 'LIC error (-wbilo_lic - runofflic*frac_lic)                   = {:12.4e} (rel) '.format ( (LIC_flx_budget3-dLIC_mas_wat)/dLIC_mas_wat) ) 
    1242  
    1243 if SRF : 
     1116echo ( f'Mass total   begin = {ATM.LIC_mass_wat_beg    :12.6e} kg | Mass end = {ATM.LIC_mass_wat_end    :12.6e} kg' ) 
     1117echo ( f'Mass snow    begin = {ATM.LIC_mass_sno_beg    :12.6e} kg | Mass end = {ATM.LIC_mass_sno_end    :12.6e} kg' ) 
     1118echo ( f'Mass qs      begin = {ATM.LIC_mass_qs_beg     :12.6e} kg | Mass end = {ATM.LIC_mass_qs_end     :12.6e} kg' ) 
     1119echo ( f'Mass runlic0 begin = {ATM.LIC_mass_runlic0_beg:12.6e} kg | Mass end = {ATM.LIC_mass_runlic0_end:12.6e} kg' ) 
     1120prtFlux ( 'dmass (LIC sno)       ', ATM.LIC_diff_mass_sno          , 'f', True, width=45 ) 
     1121prtFlux ( 'dmass (LIC qs)        ', ATM.LIC_diff_mass_qs           , 'e', True, width=45 ) 
     1122prtFlux ( 'dmass (LIC wat)       ', ATM.LIC_diff_mass_wat          , 'f', True, width=45 ) 
     1123prtFlux ( 'dmass (LIC runlic0)   ', ATM.LIC_diff_mass_runlic0      , 'e', True, width=45 ) 
     1124prtFlux ( 'dmass (LIC total)     ', ATM.LIC_diff_mass_wat          , 'e', True, width=45 ) 
     1125prtFlux ( 'LIC ATM.flx_wemp_lic  ', ATM.flx_wemp_lic          , 'f', True, width=45 ) 
     1126prtFlux ( 'LIC LIC_flx_fqfonte   ', ATM.LIC_flx_fqfonte       , 'f', True, width=45 ) 
     1127prtFlux ( 'LIC LIC_flx_calving   ', ATM.LIC_flx_calving       , 'f', True, width=45 ) 
     1128prtFlux ( 'LIC LIC_flx_runofflic ', ATM.LIC_flx_runlic        , 'f', True, width=45 ) 
     1129prtFlux ( 'LIC fqfonte + calving ', ATM.LIC_flx_calving+ATM.LIC_flx_fqfonte , 'f', True, width=45 ) 
     1130prtFlux ( 'LIC fluxes 1 ( wemp_lic   - fqcalving - fqfonte)) ', ATM.LIC_flx_budget1              , 'f', True, width=45 ) 
     1131prtFlux ( 'LIC fluxes 2 (-wbilo_lic - fqcalving - fqfonte)   ', ATM.LIC_flx_budget2              , 'f', True, width=45 ) 
     1132prtFlux ( 'LIC fluxes 3 (-wbilo_lic - runofflic*frac_lic)    ', ATM.LIC_flx_budget3              , 'f', True, width=45 ) 
     1133prtFlux ( 'LIC fluxes 3 ( wemp_lic  - runofflic*frac_lic)    ', ATM.LIC_flx_budget4              , 'f', True, width=45 ) 
     1134prtFlux ( 'LIC error 1                                       ', ATM.LIC_flx_budget1-ATM.LIC_diff_mass_wat , 'e', True, width=45 ) 
     1135prtFlux ( 'LIC error 2                                       ', ATM.LIC_flx_budget2-ATM.LIC_diff_mass_wat , 'e', True, width=45 ) 
     1136prtFlux ( 'LIC error 3                                       ', ATM.LIC_flx_budget3-ATM.LIC_diff_mass_wat , 'e', True, width=45 ) 
     1137echo ( 'LIC error (wevap - precip*frac_lic  - fqcalving - fqfonte)    = {:12.4e} (rel) '.format ( (ATM.LIC_flx_budget1-ATM.LIC_diff_mass_wat)/ATM.LIC_diff_mass_wat) ) 
     1138echo ( 'LIC error (-wbilo_lic - fqcalving - fqfonte)                  = {:12.4e} (rel) '.format ( (ATM.LIC_flx_budget2-ATM.LIC_diff_mass_wat)/ATM.LIC_diff_mass_wat) ) 
     1139echo ( 'LIC error (-wbilo_lic - runofflic*frac_lic)                   = {:12.4e} (rel) '.format ( (ATM.LIC_flx_budget3-ATM.LIC_diff_mass_wat)/ATM.LIC_diff_mass_wat) ) 
     1140 
     1141if SECHIBA : 
    12441142    echo ( '\n====================================================================================' ) 
    12451143    echo ( f'-- SECHIBA fluxes  -- {Title} ' ) 
    12461144 
    1247     SRF_flx_rain     = SRF_flux_int ( SRF_rain     ) 
    1248     SRF_flx_evap     = SRF_flux_int ( SRF_evap     ) 
    1249     SRF_flx_snowf    = SRF_flux_int ( SRF_snowf    ) 
    1250     SRF_flx_subli    = SRF_flux_int ( SRF_subli    ) 
    1251     SRF_flx_transpir = SRF_flux_int ( SRF_transpir ) 
    1252     SRF_flx_emp      = SRF_flux_int ( SRF_emp      ) 
    1253  
    1254     RUN_flx_torouting   = SRF_flux_int  ( RUN_runoff       + RUN_drainage) 
    1255     RUN_flx_fromrouting = ONE_flux_int  ( RUN_coastalflow + RUN_riverflow ) 
    1256  
    1257     SRF_flx_all =  SRF_flux_int  ( SRF_rain + SRF_snowf - SRF_evap - RUN_runoff - RUN_drainage ) 
    1258  
    1259     prtFlux ('rain         ', SRF_flx_rain       , 'f' ) 
    1260     prtFlux ('evap         ', SRF_flx_evap       , 'f' ) 
    1261     prtFlux ('snowf        ', SRF_flx_snowf      , 'f' ) 
    1262     prtFlux ('E-P          ', SRF_flx_emp        , 'f' ) 
    1263     prtFlux ('subli        ', SRF_flx_subli      , 'f' ) 
    1264     prtFlux ('transpir     ', SRF_flx_transpir   , 'f' ) 
    1265     prtFlux ('to routing   ', RUN_flx_torouting  , 'f' ) 
    1266     prtFlux ('budget       ', SRF_flx_all        , 'f', small=True ) 
     1145    SRF.flx_rain     = SRF_flux_int ( SRF.rain     ) 
     1146    SRF.flx_evap     = SRF_flux_int ( SRF.evap     ) 
     1147    SRF.flx_snowf    = SRF_flux_int ( SRF.snowf    ) 
     1148    SRF.flx_subli    = SRF_flux_int ( SRF.subli    ) 
     1149    SRF.flx_transpir = SRF_flux_int ( SRF.transpir ) 
     1150    SRF.flx_emp      = SRF_flux_int ( SRF.emp      ) 
     1151 
     1152    SRF.RUN_flx_torouting   = SRF_flux_int  ( SRF.RUN_runoff      + SRF.RUN_drainage) 
     1153    SRF.RUN_flx_fromrouting = ONE_flux_int  ( SRF.RUN_coastalflow + SRF.RUN_riverflow ) 
     1154 
     1155    SRF.flx_all =  SRF_flux_int  ( SRF.rain + SRF.snowf - SRF.evap - SRF.RUN_runoff - SRF.RUN_drainage ) 
     1156 
     1157    prtFlux ('rain         ', SRF.flx_rain       , 'f' ) 
     1158    prtFlux ('evap         ', SRF.flx_evap       , 'f' ) 
     1159    prtFlux ('snowf        ', SRF.flx_snowf      , 'f' ) 
     1160    prtFlux ('E-P          ', SRF.flx_emp        , 'f' ) 
     1161    prtFlux ('subli        ', SRF.flx_subli      , 'f' ) 
     1162    prtFlux ('transpir     ', SRF.flx_transpir   , 'f' ) 
     1163    prtFlux ('to routing   ', SRF.RUN_flx_torouting  , 'f' ) 
     1164    prtFlux ('budget       ', SRF.flx_all        , 'f', small=True ) 
    12671165 
    12681166    echo ( '\n------------------------------------------------------------------------------------' ) 
    12691167    echo ( 'Water content in surface ' ) 
    1270     echo ( f'SRF_mas_wat_beg  = {SRF_mas_wat_beg:12.6e} kg | SRF_mas_wat_end  = {SRF_mas_wat_end:12.6e} kg ' ) 
    1271     prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', small=True) 
    1272     prtFlux ( 'Error            ',  SRF_flx_all-dSRF_mas_wat, 'e', small=True ) 
    1273     echo ( 'dMass (water srf) = {:12.4e} (rel)   '.format ( (SRF_flx_all-dSRF_mas_wat)/dSRF_mas_wat) ) 
     1168    echo ( f'SRF.mass_wat_beg  = {SRF.mass_wat_beg:12.6e} kg | SRF.mass_wat_end  = {SRF.mass_wat_end:12.6e} kg ' ) 
     1169    prtFlux ( 'dMass (water srf)', SRF.diff_mass_wat, 'e', small=True) 
     1170    prtFlux ( 'Error            ', SRF.flx_all-SRF.diff_mass_wat, 'e', small=True ) 
     1171    echo ( 'dMass (water srf) = {:12.4e} (rel)   '.format ( (SRF.flx_all-SRF.diff_mass_wat)/SRF.diff_mass_wat) ) 
    12741172 
    12751173    echo ( '\n====================================================================================' ) 
    1276     echo ( f'-- Check ATM vs. SRF -- {Title} ' ) 
    1277     prtFlux ('E-P ATM       ', ATM_flx_wemp_ter   , 'f' ) 
    1278     prtFlux ('wbilo ter     ', ATM_flx_wbilo_ter  , 'f' ) 
    1279     prtFlux ('E-P SRF       ', SRF_flx_emp        , 'f' ) 
    1280     prtFlux ('SRF/ATM error ', ATM_flx_wbilo_ter - SRF_flx_emp, 'e', True) 
    1281     echo ( 'SRF/ATM error {:12.3e} (rel)  '.format ( (ATM_flx_wbilo_ter - SRF_flx_emp)/SRF_flx_emp ) ) 
     1174    echo ( f'-- Check ATM vs. SECHIBA -- {Title} ' ) 
     1175    prtFlux ('E-P ATM       ', ATM.flx_wemp_ter   , 'f' ) 
     1176    prtFlux ('wbilo ter     ', ATM.flx_wbilo_ter  , 'f' ) 
     1177    prtFlux ('E-P SECHIBA   ', SRF.flx_emp        , 'f' ) 
     1178    prtFlux ('SRF/ATM error ', ATM.flx_wbilo_ter - SRF.flx_emp, 'e', True) 
     1179    echo ( 'SRF/ATM error {:12.3e} (rel)  '.format ( (ATM.flx_wbilo_ter - SRF.flx_emp)/SRF.flx_emp ) ) 
    12821180 
    12831181    echo ('') 
    12841182    echo ( '\n====================================================================================' ) 
    12851183    echo ( f'-- RUNOFF fluxes  -- {Title} ' ) 
    1286     RUN_flx_all = RUN_flx_torouting - RUN_flx_river - RUN_flx_coastal 
    1287     prtFlux ('runoff    ', RUN_flx_runoff      , 'f' ) 
    1288     prtFlux ('drainage  ', RUN_flx_drainage    , 'f' ) 
    1289     prtFlux ('run+drain ', RUN_flx_torouting   , 'f' ) 
    1290     prtFlux ('river     ', RUN_flx_river       , 'f' ) 
    1291     prtFlux ('coastal   ', RUN_flx_coastal     , 'f' ) 
    1292     prtFlux ('riv+coa   ', RUN_flx_fromrouting , 'f' ) 
    1293     prtFlux ('budget    ', RUN_flx_all         , 'f' , small=True) 
     1184    SRF.RUN_flx_all = SRF.RUN_flx_torouting - SRF.RUN_flx_river - SRF.RUN_flx_coastal 
     1185    prtFlux ('runoff    ', SRF.RUN_flx_runoff      , 'f' ) 
     1186    prtFlux ('drainage  ', SRF.RUN_flx_drainage    , 'f' ) 
     1187    prtFlux ('run+drain ', SRF.RUN_flx_torouting   , 'f' ) 
     1188    prtFlux ('river     ', SRF.RUN_flx_river       , 'f' ) 
     1189    prtFlux ('coastal   ', SRF.RUN_flx_coastal     , 'f' ) 
     1190    prtFlux ('riv+coa   ', SRF.RUN_flx_fromrouting , 'f' ) 
     1191    prtFlux ('budget    ', SRF.RUN_flx_all         , 'f' , small=True) 
    12941192 
    12951193    echo ( '\n------------------------------------------------------------------------------------' ) 
    12961194    echo ( f'Water content in routing+lake  -- {Title} ' ) 
    1297     echo ( f'RUN_mas_wat_beg  = {RUN_mas_wat_beg:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg ' ) 
    1298     prtFlux ( 'dMass (routing)  ', dRUN_mas_wat+dSRF_mas_lake, 'f', small=True) 
    1299     prtFlux ( 'Routing error    ', RUN_flx_all+dSRF_mas_lake-dRUN_mas_wat, 'e', small=True ) 
    1300     echo ( 'Routing error : {:12.3e} (rel)'.format ( (RUN_flx_all-dSRF_mas_lake-dRUN_mas_wat)/(dRUN_mas_wat+dSRF_mas_lake) ) ) 
     1195    echo ( f'SRF.RUN_mass_wat_beg  = {SRF.RUN_mass_wat_beg:12.6e} kg | SRF.RUN_mass_wat_end = {SRF.RUN_mass_wat_end:12.6e} kg ' ) 
     1196    prtFlux ( 'dMass (routing)  ', SRF.RUN_diff_mass_wat+SRF.diff_mass_lake, 'f', small=True) 
     1197    prtFlux ( 'Routing error    ', SRF.RUN_flx_all+SRF.diff_mass_lake-SRF.RUN_diff_mass_wat, 'e', small=True ) 
     1198    echo ( 'Routing error : {:12.3e} (rel)'.format ( (SRF.RUN_flx_all-SRF.diff_mass_lake-SRF.RUN_diff_mass_wat)/(SRF.RUN_diff_mass_wat+SRF.diff_mass_lake) ) ) 
    13011199 
    13021200    echo ( '\n------------------------------------------------------------------------------------' ) 
    13031201    echo ( f'Water content in routing  -- {Title} ' ) 
    1304     echo ( f'RUN_mas_wat_beg  = {RUN_mas_wat_beg:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg ' ) 
    1305     prtFlux ( 'dMass (routing) ', dRUN_mas_wat, 'f', small=True  ) 
    1306     prtFlux ( 'Routing error   ', RUN_flx_all-dRUN_mas_wat, 'e', small=True) 
    1307     echo ( 'Routing error : {:12.3e} (rel)'.format ( (RUN_flx_all-dRUN_mas_wat)/dRUN_mas_wat ) ) 
     1202    echo ( f'SRF.RUN_mass_wat_beg  = {SRF.RUN_mass_wat_beg:12.6e} kg | SRF.RUN_mass_wat_end = {SRF.RUN_mass_wat_end:12.6e} kg ' ) 
     1203    prtFlux ( 'dMass (routing) ', SRF.RUN_diff_mass_wat, 'f', small=True  ) 
     1204    prtFlux ( 'Routing error   ', SRF.RUN_flx_all-SRF.RUN_diff_mass_wat, 'e', small=True) 
     1205    echo ( 'Routing error : {:12.3e} (rel)'.format ( (SRF.RUN_flx_all-SRF.RUN_diff_mass_wat)/SRF.RUN_diff_mass_wat ) ) 
    13081206 
    13091207echo ( ' ' ) 
Note: See TracChangeset for help on using the changeset viewer.