Ignore:
Timestamp:
10/10/23 12:58:04 (7 months ago)
Author:
omamce
Message:

O.M. :

New version of WATER_BUDGET

  • Conservation in NEMO are OK
  • Conservation in Sechiba are OK, for both ICO and latlon grids
  • Problems in atmosphere, LIC surface and ocean/atmosphere coherence
File:
1 edited

Legend:

Unmodified
Added
Removed

  • TOOLS/WATER_BUDGET/ATM_waterbudget.py

    r6620 r6647  
    3434    raise Exception ( "Minimum Python version is 3.8" ) 
    3535 
    36 ## Import local module 
     36## Import local modules 
    3737import WaterUtils as wu 
    3838import libIGCM_sys 
    3939import nemo, lmdz 
     40 
     41from WaterUtils import VarInt, Rho, Ra, Grav, ICE_rho_ice, ICE_rho_sno, OCE_rho_liq, ATM_rho, SRF_rho, RUN_rho, ICE_rho_pnd, YearLength 
    4042 
    4143## Creates parser for reading .ini input file 
    … …  
    4749## --------------------- 
    4850ATM=None ; ATM_HIS='latlon' ; SRF_HIS='latlon' ; RUN_HIS='latlon' ; ORCA=None ; NEMO=None ; OCE_relax=False 
    49 OCE_icb=False ; Coupled=False ; Routing=None 
     51OCE_icb=False ; Coupled=False ; Routing=None ; TestInterp=None 
    5052TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=None 
     53YearBegin=None ; YearEnd=None ; DateBegin=None ; DateEnd=None 
    5154 
    5255## 
    … …  
    6669tar_restart_end_RUN=None ; tar_restart_end_OCE=None ; tar_restart_end_ICE=None 
    6770ContinueOnError=False ; ErrorCount=0 ; SortIco = False 
     71 
     72## 
     73## Precision of history file reading 
     74## --------------------------------- 
     75# Default is float (full precision). Degrade the precision by using np.float32 
     76# Restart file are always read at the full precision 
     77readPrec=float 
    6878 
    6979## Read command line arguments 
    … …  
    138148print ( f'\nConfig file readed : {IniFile} ' ) 
    139149 
     150## 
     151## Reading prec 
     152if wu.unDefined ( 'readPrec' )  : 
     153    readPrec = np.float64 
     154else : 
     155    if readPrec in ["float", "float64", "r8", "double", "<class 'float'>"         ] : readPrec = float 
     156    if readPrec in [         "float32", "r4", "single", "<class 'numpy.float32'>" ] : readPrec = np.float32 
     157    if readPrec in [         "float16", "r2", "half"  , "<class 'numpy.float16'>" ] : readPrec = np.float16 
     158     
    140159## Some physical constants 
    141160## ======================= 
    … …  
    155174if not 'Files' in config.keys () : config['Files'] = {} 
    156175 
    157 config['Physics'] = { 'Ra':Ra, 'Grav':Grav, 'ICE_rho_ice':ICE_rho_ice, 'ICE_rho_sno':ICE_rho_sno, 
    158                       'OCE_rho_liq':OCE_rho_liq, 'ATM_rho':ATM_rho, 'SRF_rho':SRF_rho, 'RUN_rho':RUN_rho} 
    159  
    160 config['Config'] = { 'ContinueOnError':ContinueOnError, 'SortIco':SortIco } 
     176config['Physics'] = { 'Ra':str(Ra), 'Grav':str(Grav), 'ICE_rho_ice':str(ICE_rho_ice), 'ICE_rho_sno':str(ICE_rho_sno), 
     177                      'OCE_rho_liq':str(OCE_rho_liq), 'ATM_rho':str(ATM_rho), 'SRF_rho':str(SRF_rho), 'RUN_rho':str(RUN_rho)} 
     178 
     179config['Config'] = { 'ContinueOnError':str(ContinueOnError), 'SortIco':str(SortIco), 'TestInterp':str(TestInterp), 'readPrec':str(readPrec) } 
    161180 
    162181## -------------------------- 
    … …  
    176195if wu.unDefined ( 'DateBegin' ) : 
    177196    DateBegin = f'{YearBegin}0101' 
    178     config['Experiment']['DateBegin'] = DateBegin 
     197    config['Experiment']['DateBegin'] = str(DateBegin) 
    179198else : 
    180     YearBegin, MonthBegin, DayBegin = wu.DateSplit ( DateBegin ) 
    181     DateBegin = wu.FormatToGregorian ( DateBegin) 
    182     config['Experiment']['YearBegin'] = YearBegin 
     199    YearBegin, MonthBegin, DayBegin = wu.SplitDate ( DateBegin ) 
     200    DateBegin = wu.FormatToGregorian (DateBegin) 
     201    config['Experiment']['YearBegin'] = str(YearBegin) 
    183202 
    184203if wu.unDefined ( 'DateEnd' )   : 
    185204    DateEnd   = f'{YearEnd}1231' 
    186     config['Experiment']['DateEnd']   = DateEnd 
     205    config['Experiment']['DateEnd'] = str(DateEnd) 
    187206else : 
    188     YearEnd, MonthEnd, DayEnd = wu.DateSplit ( DateEnd ) 
    189     DateEnd   = wu.FormatToGregorian ( DateEnd) 
     207    YearEnd, MonthEnd, DayEnd = wu.SplitDate ( DateEnd ) 
     208    DateEnd   = wu.FormatToGregorian (DateEnd) 
     209    config['Experiment']['DateEnd'] = str(DateEnd) 
    190210 
    191211if wu.unDefined ( 'PackFrequency' ) : 
    … …  
    200220## Output file with water budget diagnostics 
    201221## ----------------------------------------- 
    202 if wu.unDefined ( 'FileOut' ) : FileOut = f'ATM_waterbudget_{JobName}_{DateBegin}_{DateEnd}.out' 
     222if wu.unDefined ( 'FileOut' ) : 
     223    FileOut = f'ATM_waterbudget_{JobName}_{YearBegin}_{YearEnd}' 
     224    if ICO :  
     225        if ATM_HIS == 'latlon' : FileOut = f'{FileOut}_LATLON' 
     226        if ATM_HIS == 'ico'    : FileOut = f'{FileOut}_ICO' 
     227    if readPrec == np.float32  : FileOut = f'{FileOut}_float32' 
     228    FileOut = f'{FileOut}.out' 
     229 
    203230config['Files']['FileOut'] = FileOut 
    204231 
    … …  
    207234## Useful functions 
    208235## ---------------- 
     236if readPrec == float : 
     237    def rprec (tab) : return tab 
     238else :  
     239    def rprec (tab) : return tab.astype(readPrec).astype(float) 
     240     
    209241def kg2Sv    (val, rho=ATM_rho) : 
    210242    '''From kg to Sverdrup''' 
    … …  
    216248 
    217249def var2prt (var, small=False, rho=ATM_rho) : 
    218     if small :  return var , kg2Sv(var, rho=rho)*1000., kg2myear(var, rho=rho)*100# 
     250    if small :  return var , kg2Sv(var, rho=rho)*1000., kg2myear(var, rho=rho)*1000 
    219251    else     :  return var , kg2Sv(var, rho=rho)      , kg2myear(var, rho=rho) 
    220252 
    … …  
    236268    f_out.flush () 
    237269    return None 
    238      
     270 
     271echo ( f'{ContinueOnError = }' ) 
     272echo ( f'{SortIco         = }' ) 
     273echo ( f'{readPrec        = }' ) 
     274 
     275echo ( f'{JobName         = }' ) 
     276echo ( f'{ConfigCard      = }' ) 
     277echo ( f'{libIGCM         = }' )      
     278echo ( f'{User            = }' )        
     279echo ( f'{Group           = }' )        
     280echo ( f'{Freq            = }' )        
     281echo ( f'{YearBegin       = }' )      
     282echo ( f'{YearEnd         = }' )      
     283echo ( f'{DateBegin       = }' ) 
     284echo ( f'{DateEnd         = }' ) 
     285echo ( f'{PackFrequency   = }' )  
     286echo ( f'{ATM             = }' )        
     287echo ( f'{Routing         = }' )        
     288echo ( f'{ORCA            = }' )       
     289echo ( f'{NEMO            = }' )       
     290echo ( f'{Coupled         = }' )       
     291echo ( f'{ATM_HIS         = }' )       
     292echo ( f'{SRF_HIS         = }' )       
     293echo ( f'{RUN_HIS         = }' )       
     294 
    239295## Set libIGCM directories 
    240296## ----------------------- 
    … …  
    266322 
    267323echo (' ') 
    268 echo ( f'JobName    : {JobName}'   ) 
    269 echo ( f'Comment    : {Comment}'   ) 
    270 echo ( f'TmpDir     : {TmpDir}'    ) 
     324echo ( f'JobName     : {JobName}'    ) 
     325echo ( f'Comment     : {Comment}'    ) 
     326echo ( f'TmpDir      : {TmpDir}'     ) 
    271327echo ( f'FileDir     : {FileDir}'    ) 
    272328echo ( f'FileDirOCE  : {FileDirOCE}' ) 
    … …  
    598654# ATM grid with cell surfaces 
    599655if LMDZ : 
    600     ATM_lat  = lmdz.geo2point (   d_ATM_his ['lat']+0*d_ATM_his ['lon'], dim1D='cell_latlon' ) 
    601     ATM_lon  = lmdz.geo2point ( 0*d_ATM_his ['lat']+  d_ATM_his ['lon'], dim1D='cell_latlon' ) 
    602     ATM_aire = lmdz.geo2point ( d_ATM_his ['aire']     [0], cumulPoles=True, dim1D='cell_latlon' ) 
    603     ATM_fter = lmdz.geo2point ( d_ATM_his ['fract_ter'][0], dim1D='cell_latlon' ) 
    604     ATM_foce = lmdz.geo2point ( d_ATM_his ['fract_oce'][0], dim1D='cell_latlon' ) 
    605     ATM_fsic = lmdz.geo2point ( d_ATM_his ['fract_sic'][0], dim1D='cell_latlon' ) 
    606     ATM_flic = lmdz.geo2point ( d_ATM_his ['fract_lic'][0], dim1D='cell_latlon' ) 
    607     SRF_lat  = lmdz.geo2point (   d_SRF_his ['lat']+0*d_SRF_his ['lon'], dim1D='cell_latlon' ) 
    608     SRF_lon  = lmdz.geo2point ( 0*d_SRF_his ['lat']+  d_SRF_his ['lon'], dim1D='cell_latlon' ) 
    609     SRF_aire      = lmdz.geo2point ( d_SRF_his ['Areas'] * d_SRF_his ['Contfrac'], dim1D='cell_latlonl', cumulPoles=True ) 
    610     SRF_areas     = lmdz.geo2point ( d_SRF_his ['Areas'],  dim1D='cell_latlonl', cumulPoles=True ) 
    611     SRF_contfrac  = lmdz.geo2point ( d_SRF_his ['Contfrac'], dim1D='cell_latlonl' ) 
    612  
    613     #SRF_res_aire = SRF_aire 
    614      
     656    echo ('ATM grid with cell surfaces : LMDZ') 
     657    ATM_lat       = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1D='cell' ) 
     658    ATM_lon       = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1D='cell' ) 
     659    ATM_aire      = lmdz.geo2point ( rprec (d_ATM_his ['aire'] )    [0], cumulPoles=True, dim1D='cell' ) 
     660    ATM_fter      = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1D='cell' ) 
     661    ATM_foce      = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1D='cell' ) 
     662    ATM_fsic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1D='cell' ) 
     663    ATM_flic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1D='cell' ) 
     664    SRF_lat       = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1D='cell' ) 
     665    SRF_lon       = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1D='cell' ) 
     666    SRF_aire      = lmdz.geo2point ( rprec (d_SRF_his ['Areas']) * rprec (d_SRF_his ['Contfrac']), dim1D='cell', cumulPoles=True ) 
     667    SRF_areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas'])  ,  dim1D='cell', cumulPoles=True ) 
     668    SRF_contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac']), dim1D='cell' ) 
     669 
    615670if ICO : 
    616671    if ATM_HIS == 'latlon' : 
    617         jpja, jpia = d_ATM_his['aire'][0].shape 
    618         if wu.unDefined ( 'file_ATM_aire' ) :  
    619             file_ATM_aire = os.path.join ( R_IN, 'ATM', 'GRID', f'aire_{ATM}_to_{jpia}x{jpja}.nc' ) 
    620             config['Files']['file_ATM_aire'] = file_ATM_aire 
    621         echo ( f'Aire sur grille reguliere : {file_ATM_aire = }' ) 
    622         d_ATM_aire = xr.open_dataset ( file_ATM_aire, decode_times=False ).squeeze() 
    623         try : 
    624             ATM_lat  = lmdz.geo2point (   d_ATM_his ['lat']+0*d_ATM_his ['lon'], dim1D='cell_latlon' ) 
    625             ATM_lon  = lmdz.geo2point ( 0*d_ATM_his ['lat']+  d_ATM_his ['lon'], dim1D='cell_latlon' ) 
    626         except : 
    627             ATM_lat  = lmdz.geo2point (   d_ATM_his ['lat_dom_out']+0*d_ATM_his ['lon_dom_out'], dim1D='cell_latlon' ) 
    628             ATM_lon  = lmdz.geo2point ( 0*d_ATM_his ['lat_dom_out']+  d_ATM_his ['lon_dom_out'], dim1D='cell_latlon' ) 
    629         ATM_aire = lmdz.geo2point ( d_ATM_aire ['aire'].squeeze(), cumulPoles=True, dim1D='cell_latlon' ) 
    630         ATM_fter = lmdz.geo2point ( d_ATM_his ['fract_ter'][0], dim1D='cell_latlon' ) 
    631         ATM_foce = lmdz.geo2point ( d_ATM_his ['fract_oce'][0], dim1D='cell_latlon' ) 
    632         ATM_fsic = lmdz.geo2point ( d_ATM_his ['fract_sic'][0], dim1D='cell_latlon' ) 
    633         ATM_flic = lmdz.geo2point ( d_ATM_his ['fract_lic'][0], dim1D='cell_latlon' ) 
     672        echo ( 'ATM areas and fractions on latlon grid' ) 
     673        if 'lat_dom_out' in d_ATM_his.variables : 
     674            ATM_lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat_dom_out'])+0*rprec (d_ATM_his ['lon_dom_out']), dim1D='cell' ) 
     675            ATM_lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat_dom_out'])+  rprec (d_ATM_his ['lon_dom_out']), dim1D='cell' ) 
     676        else : 
     677            ATM_lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1D='cell' ) 
     678            ATM_lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1D='cell' ) 
     679        ATM_aire = lmdz.geo2point ( rprec (d_ATM_his ['aire'][0]).squeeze(), cumulPoles=True, dim1D='cell' ) 
     680        ATM_fter = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1D='cell' ) 
     681        ATM_foce = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1D='cell' ) 
     682        ATM_fsic = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1D='cell' ) 
     683        ATM_flic = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1D='cell' ) 
     684 
     685    if ATM_HIS == 'ico' : 
     686        echo ( 'ATM areas and fractions on ICO grid' ) 
     687        ATM_aire =  rprec (d_ATM_his ['aire']     [0])[ATM_his_keysort].squeeze() 
     688        ATM_lat  =  rprec (d_ATM_his ['lat']         )[ATM_his_keysort] 
     689        ATM_lon  =  rprec (d_ATM_his ['lon']         )[ATM_his_keysort] 
     690        ATM_fter =  rprec (d_ATM_his ['fract_ter'][0])[ATM_his_keysort] 
     691        ATM_foce =  rprec (d_ATM_his ['fract_oce'][0])[ATM_his_keysort] 
     692        ATM_fsic =  rprec (d_ATM_his ['fract_sic'][0])[ATM_his_keysort] 
     693        ATM_flic =  rprec (d_ATM_his ['fract_lic'][0])[ATM_his_keysort] 
    634694 
    635695    if SRF_HIS == 'latlon' : 
    636         try : 
    637             SRF_lat  = lmdz.geo2point (   d_SRF_his ['lat']+0*d_SRF_his ['lon'], dim1D='cell_latlon' ) 
    638             SRF_lon  = lmdz.geo2point ( 0*d_SRF_his ['lat']+  d_SRF_his ['lon'], dim1D='cell_latlon' ) 
    639         except : 
    640             try :  
    641                 SRF_lat  = lmdz.geo2point (   d_SRF_his ['lat_dom_out']+0*d_SRF_his ['lon_dom_out'], dim1D='cell_latlon' ) 
    642                 SRF_lon  = lmdz.geo2point ( 0*d_SRF_his ['lat_dom_out']+  d_SRF_his ['lon_dom_out'], dim1D='cell_latlon' ) 
    643             except :  
    644                 SRF_lat  = lmdz.geo2point (   d_SRF_his ['lat_domain_landpoints_out']+0*d_SRF_his ['lon_domain_landpoints_out'], dim1D='cell_latlon' ) 
    645                 SRF_lon  = lmdz.geo2point ( 0*d_SRF_his ['lat_domain_landpoints_out']+  d_SRF_his ['lon_domain_landpoints_out'], dim1D='cell_latlon' ) 
    646         SRF_areas    = d_SRF_his ['Areas'].values 
    647         SRF_areafrac = d_SRF_his ['AreaFrac'].values 
    648         SRF_areas    = xr.DataArray ( SRF_areas   , coords=d_SRF_his ['Contfrac'].coords, dims=d_SRF_his ['Contfrac'].dims) 
    649         SRF_areafrac = xr.DataArray ( SRF_areafrac, coords=d_SRF_his ['Contfrac'].coords, dims=d_SRF_his ['Contfrac'].dims) 
    650  
    651         SRF_areas     = lmdz.geo2point ( SRF_areas   , dim1D='cell_latlon', cumulPoles=True ) 
    652         SRF_areafrac  = lmdz.geo2point ( SRF_areafrac, dim1D='cell_latlon', cumulPoles=True ) 
    653         SRF_contfrac  = SRF_areafrac / SRF_areas 
    654         SRF_aire      = SRF_areafrac 
    655  
    656     if ATM_HIS == 'ico' : 
    657         echo ( f'Grille icosaedre' ) 
    658         ATM_aire =  d_ATM_his ['aire']     [0][ATM_his_keysort].squeeze() 
    659         ATM_lat  =  d_ATM_his ['lat']         [ATM_his_keysort] 
    660         ATM_lon  =  d_ATM_his ['lon']         [ATM_his_keysort] 
    661         ATM_fter =  d_ATM_his ['fract_ter'][0][ATM_his_keysort] 
    662         ATM_foce =  d_ATM_his ['fract_oce'][0][ATM_his_keysort] 
    663         ATM_fsic =  d_ATM_his ['fract_sic'][0][ATM_his_keysort] 
    664         ATM_flic =  d_ATM_his ['fract_lic'][0][ATM_his_keysort] 
    665          
    666     if SRF_HIS == 'ico' :  
    667         SRF_lat       =  d_SRF_his ['lat']      [SRF_his_keysort] 
    668         SRF_lon       =  d_SRF_his ['lon']      [SRF_his_keysort] 
    669         SRF_areas     =  d_SRF_his ['Areas']    [SRF_his_keysort]  
    670         #SRF_areafrac =  d_SRF_his ['AreaFrac'] [SRF_his_keysort]  
    671         SRF_contfrac  =  d_SRF_his ['Contfrac'][SRF_his_keysort] 
     696        echo ( 'SRF areas and fractions on latlon grid' ) 
     697        if 'lat_domain_landpoints_out' in d_SRF_his  :  
     698            SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat_domain_landpoints_out'])+0*rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1D='cell' ) 
     699            SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_domain_landpoints_out'])+  rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1D='cell' ) 
     700        else :  
     701            if 'lat_domain_landpoints_out' in d_SRF_his : 
     702                SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat_dom_out'])+0*rprec (d_SRF_his ['lon_dom_out']), dim1D='cell' ) 
     703                SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_dom_out'])+  rprec (d_SRF_his ['lon_dom_out']), dim1D='cell' ) 
     704            else : 
     705                SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1D='cell' ) 
     706                SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1D='cell' ) 
     707                 
     708        SRF_areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas']   )   , dim1D='cell', cumulPoles=True ) 
     709        SRF_areafrac  = lmdz.geo2point ( rprec (d_SRF_his ['AreaFrac'])   , dim1D='cell', cumulPoles=True ) 
     710        SRF_contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac'])   , dim1D='cell', cumulPoles=True ) 
     711        SRF_aire = SRF_areafrac 
     712  
     713    if SRF_HIS == 'ico' : 
     714        echo ( 'SRF areas and fractions on latlon grid' ) 
     715        SRF_lat       =  rprec (d_SRF_his ['lat']     ) [SRF_his_keysort] 
     716        SRF_lon       =  rprec (d_SRF_his ['lon']     ) [SRF_his_keysort] 
     717        SRF_areas     =  rprec (d_SRF_his ['Areas']   ) [SRF_his_keysort]  
     718        SRF_contfrac  =  rprec (d_SRF_his ['Contfrac']) [SRF_his_keysort] 
    672719        SRF_aire      =  SRF_areas * SRF_contfrac 
    673720 
    … …  
    681728ATM_aire_fsea = ATM_aire * ATM_fsea 
    682729 
    683 SRF_aire = SRF_aire.where ( np.abs (SRF_aire) < 1E15, 0. ) 
     730#SRF_aire = SRF_aire.where ( np.abs (SRF_aire) < 1E15, 0. ) 
    684731 
    685732## Write the full configuration 
    … …  
    713760    DYN_fsea = ATM_fsea 
    714761    DYN_flnd = ATM_flnd 
    715     DYN_fter = d_ATM_beg['FTER'] 
    716     DYN_flic = d_ATM_beg['FLIC'] 
     762    DYN_fter = rprec (d_ATM_beg['FTER']) 
     763    DYN_flic = rprec (d_ATM_beg['FLIC']) 
    717764    DYN_aire_fter = DYN_aire * DYN_fter 
    718765 
    … …  
    728775    return ATM_flux_int 
    729776 
     777def LIC_flux_int (flux) : 
     778    '''Integrate (* time * surface) flux on land ice grid''' 
     779    LIC_flux_int  = wu.Psum ( (flux * dtime_per_sec * ATM_aire_flic).to_masked_array().ravel() ) 
     780    return LIC_flux_int 
     781 
    730782def SRF_stock_int (stock) : 
    731783    '''Integrate (* surface) stock on land grid''' 
    … …  
    749801    return ONE_flux_int 
    750802 
     803def Sprec ( tlist ) : 
     804    '''Accurate sum of list of scalar elements''' 
     805    return wu.Psum ( np.array ( tlist) ) 
     806 
    751807ATM_aire_sea     = ATM_aire * ATM_fsea 
     808 
    752809ATM_aire_tot     = ONE_stock_int (ATM_aire) 
    753 ATM_aire_sea_tot = ONE_stock_int (ATM_aire*ATM_fsea) 
    754  
    755 DYN_aire_tot     =  ONE_stock_int ( DYN_aire ) 
    756 SRF_aire_tot     =  ONE_stock_int ( SRF_aire ) 
     810ATM_aire_sea_tot = ONE_stock_int (ATM_aire_fsea) 
     811ATM_aire_ter_tot = ONE_stock_int (ATM_aire_fter) 
     812ATM_aire_lic_tot = ONE_stock_int (ATM_aire_flic) 
     813 
     814DYN_aire_tot     = ONE_stock_int ( DYN_aire ) 
     815SRF_aire_tot     = ONE_stock_int ( SRF_aire ) 
    757816 
    758817echo ('') 
    759 echo ( f'ATM DYN : Area of atmosphere : {DYN_aire_tot:18.9e}' ) 
    760 echo ( f'ATM HIS : Area of atmosphere : {ATM_aire_tot:18.9e}' ) 
    761 echo ( f'ATM SRF : Area of atmosphere : {SRF_aire_tot:18.9e}' ) 
     818echo ( f'ATM DYN : Area of atmosphere        : {DYN_aire_tot:18.9e}' ) 
     819echo ( f'ATM HIS : Area of atmosphere        : {ATM_aire_tot:18.9e}' ) 
     820echo ( f'ATM HIS : Area of ter in atmosphere : {ATM_aire_ter_tot:18.9e}' ) 
     821echo ( f'ATM HIS : Area of lic in atmosphere : {ATM_aire_lic_tot:18.9e}' ) 
     822echo ( f'ATM SRF : Area of atmosphere        : {SRF_aire_tot:18.9e}' ) 
    762823echo ('') 
    763 echo ( 'ATM DYN : Area of atmosphere/(4pi R^2) : {:18.9f}'.format(DYN_aire_tot/(Ra*Ra*4*np.pi) ) ) 
    764 echo ( 'ATM HIS : Area of atmosphere/(4pi R^2) : {:18.9f}'.format(ATM_aire_tot/(Ra*Ra*4*np.pi) ) ) 
    765 echo ( 'ATM SRF : Area of atmosphere/(4pi R^2) : {:18.9f}'.format(SRF_aire_tot/(Ra*Ra*4*np.pi) ) ) 
     824echo ( 'ATM DYN : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(DYN_aire_tot/(Ra*Ra*4*np.pi) ) ) 
     825echo ( 'ATM HIS : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(ATM_aire_tot/(Ra*Ra*4*np.pi) ) ) 
     826echo ( 'ATM HIS : Area of ter in atmosphere/(4pi R^2) : {:18.9f}'.format(ATM_aire_ter_tot/(Ra*Ra*4*np.pi) ) ) 
     827echo ( 'ATM SRF : Area of atmosphere/(4pi R^2)        : {:18.9f}'.format(SRF_aire_tot/(Ra*Ra*4*np.pi) ) ) 
    766828echo ('') 
    767829echo ( f'ATM SRF : Area of atmosphere (no contfrac): {ONE_stock_int (SRF_areas):18.9e}' ) 
    … …  
    786848    DYN_psol_end = d_DYN_end['ps'][DYN_end_keysort] 
    787849if LMDZ :  
    788     DYN_psol_beg = lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)), dim1D='cell_latlon' ) 
    789     DYN_psol_end = lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)), dim1D='cell_latlon' ) 
     850    DYN_psol_beg = lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)), dim1D='cell' ) 
     851    DYN_psol_end = lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)), dim1D='cell' ) 
    790852     
    791853echo ( '3D Pressure at the interface layers (not scalar points)' ) 
    … …  
    793855DYN_pres_end = ATM_Ahyb + ATM_Bhyb * DYN_psol_end 
    794856 
     857echo ( 'Check computation of pressure levels' ) 
     858 
     859ind = np.empty (8) 
     860 
     861ind[0] = (DYN_pres_beg[ 0]-DYN_psol_beg).min().item() 
     862ind[1] = (DYN_pres_beg[ 0]-DYN_psol_beg).max().item() 
     863ind[2] = (DYN_pres_beg[-1]).min().item() 
     864ind[3] = (DYN_pres_beg[-1]).max().item() 
     865ind[4] = (DYN_pres_end[ 0]-DYN_psol_end).min().item() 
     866ind[5] = (DYN_pres_end[ 0]-DYN_psol_end).max().item() 
     867ind[6] = (DYN_pres_end[-1]).min().item() 
     868ind[7] = (DYN_pres_end[-1]).max().item() 
     869 
     870if any ( ind != 0) : 
     871    echo ( 'All values should be zero' ) 
     872    echo ( f'(DYN_pres_beg[ 0]-DYN_psol_beg).min().item() = {ind[0]}' ) 
     873    echo ( f'(DYN_pres_beg[ 0]-DYN_psol_beg).max().item() = {ind[1]}' ) 
     874    echo ( f'(DYN_pres_beg[-1]).min().item()              = {ind[2]}' ) 
     875    echo ( f'(DYN_pres_beg[-1]).max().item()              = {ind[3]}' ) 
     876    echo ( f'(DYN_pres_end[ 0]-DYN_psol_end).min().item() = {ind[4]}' ) 
     877    echo ( f'(DYN_pres_end[ 0]-DYN_psol_end).max().item() = {ind[5]}' ) 
     878    echo ( f'(DYN_pres_end[-1]).min().item()              = {ind[6]}' ) 
     879    echo ( f'(DYN_pres_end[-1]).max().item()              = {ind[7]}' ) 
     880    raise Exception 
     881     
    795882klevp1 = ATM_Bhyb.shape[-1] 
    796883cell   = DYN_psol_beg.shape[-1] 
    … …  
    798885 
    799886echo ( 'Layer thickness (pressure)' ) 
    800 DYN_dpres_beg = xr.DataArray ( np.empty( (klev, cell )), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) ) 
    801 DYN_dpres_end = xr.DataArray ( np.empty( (klev, cell )), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) ) 
    802      
    803 for k in np.arange (klevp1-1) : 
    804     DYN_dpres_beg[k,:] = DYN_pres_beg[k,:] - DYN_pres_beg[k+1,:] 
    805     DYN_dpres_end[k,:] = DYN_pres_end[k,:] - DYN_pres_end[k+1,:] 
     887DYN_mass_beg = xr.DataArray ( np.empty( (klev, cell)), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) ) 
     888DYN_mass_end = xr.DataArray ( np.empty( (klev, cell)), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) ) 
     889     
     890for k in np.arange (klev) : 
     891    DYN_mass_beg[k,:] = ( DYN_pres_beg[k,:] - DYN_pres_beg[k+1,:] ) / Grav 
     892    DYN_mass_end[k,:] = ( DYN_pres_end[k,:] - DYN_pres_end[k+1,:] ) / Grav 
     893 
     894DYN_mass_beg_2D = DYN_mass_beg.sum (dim='sigs') 
     895DYN_mass_end_2D = DYN_mass_end.sum (dim='sigs') 
     896 
     897DYN_mas_air_beg = DYN_stock_int ( DYN_mass_beg_2D ) 
     898DYN_mas_air_end = DYN_stock_int ( DYN_mass_end_2D ) 
    806899 
    807900echo ( 'Vertical and horizontal integral, and sum of liquid, solid and vapor water phases' ) 
    808901if LMDZ : 
    809     try :  
    810         DYN_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2Ov']  + d_DYN_beg['H2Ol']  + d_DYN_beg['H2Oi'] ).isel(rlonv=slice(0,-1) ), dim1D='cell' ) 
    811         DYN_wat_end = lmdz.geo3point ( (d_DYN_end['H2Ov']  + d_DYN_end['H2Ol']  + d_DYN_end['H2Oi'] ).isel(rlonv=slice(0,-1) ), dim1D='cell' ) 
    812     except : 
     902    if 'H2Ov' in d_DYN_beg.variables : 
     903        echo ('reading LATLON : H2Ov, H2Ol, H2Oi' ) 
     904        DYN_wat_beg = lmdz.geo3point ( d_DYN_beg['H2Ov'] + d_DYN_beg['H2Ol'] + d_DYN_beg['H2Oi'].isel(rlonv=slice(0,-1) ), dim1D='cell' ) 
     905        DYN_wat_end = lmdz.geo3point ( d_DYN_end['H2Ov'] + d_DYN_end['H2Ol'] + d_DYN_end['H2Oi'].isel(rlonv=slice(0,-1) ), dim1D='cell' ) 
     906    if 'H2Ov_g' in d_DYN_beg.variables : 
     907        echo ('reading LATLON : H2O_g, H2O_l, H2O_s' ) 
    813908        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' ) 
    814909        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' ) 
    815910if ICO : 
    816     try : 
     911    if 'H2Ov_g' in d_DYN_beg.variables : 
     912        echo ('reading ICO : H2O_g, H2O_l, H2O_s' ) 
    817913        DYN_wat_beg = (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s'])[..., DYN_beg_keysort] 
    818914        DYN_wat_end = (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s'])[..., DYN_beg_keysort] 
    819          
    820     except :  
     915    elif 'H2O_g' in d_DYN_beg.variables : 
     916        echo ('reading ICO : H2O_g, H2O_l, H2O_s' ) 
     917        DYN_wat_beg = (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s'])[..., DYN_beg_keysort] 
     918        DYN_wat_end = (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s'])[..., DYN_beg_keysort] 
     919    elif 'q' in d_DYN_beg.variables : 
     920        echo ('reading ICO : q' ) 
    821921        DYN_wat_beg = (d_DYN_beg['q'].isel(nq=0) + d_DYN_beg['q'].isel(nq=1) + d_DYN_beg['q'].isel(nq=2) )[..., DYN_beg_keysort] 
    822922        DYN_wat_end = (d_DYN_end['q'].isel(nq=0) + d_DYN_end['q'].isel(nq=1) + d_DYN_end['q'].isel(nq=2) )[..., DYN_beg_keysort] 
    823          
     923 
     924if 'lev' in DYN_wat_beg.dims :  
    824925    DYN_wat_beg = DYN_wat_beg.rename ( {'lev':'sigs'} ) 
    825     DYN_wat_beg = DYN_wat_end.rename ( {'lev':'sigs'} ) 
     926    DYN_wat_end = DYN_wat_end.rename ( {'lev':'sigs'} ) 
    826927     
    827928echo ( 'Compute water content : vertical and horizontal integral' ) 
    828 DYN_mas_wat_beg = DYN_stock_int ( (DYN_dpres_beg * DYN_wat_beg).sum (dim='sigs') ) / Grav 
    829 DYN_mas_wat_end = DYN_stock_int ( (DYN_dpres_end * DYN_wat_end).sum (dim='sigs') ) / Grav 
     929 
     930DYN_wat_beg_2D =  (DYN_mass_beg * DYN_wat_beg).sum (dim='sigs') 
     931DYN_wat_end_2D =  (DYN_mass_end * DYN_wat_end).sum (dim='sigs') 
     932 
     933DYN_mas_wat_beg = DYN_stock_int ( DYN_wat_beg_2D ) 
     934DYN_mas_wat_end = DYN_stock_int ( DYN_wat_end_2D ) 
    830935 
    831936echo ( 'Variation of water content' ) 
    832937dDYN_mas_wat = DYN_mas_wat_end - DYN_mas_wat_beg 
    833938 
    834 # \([a-z,A-Z,_]*\)/dtime_sec\*1e-9  kg2Sv(\1) 
    835 # \([a-z,A-Z,_]*\)\/ATM_aire_sea_tot\/ATM_rho\/NbYear  kg2myear(\1) 
    836  
    837939echo ( f'\nChange of atmosphere water content (dynamics)  -- {Title} ' ) 
    838940echo ( '------------------------------------------------------------------------------------' ) 
    839 echo ( 'DYN_mas_beg = {:12.6e} kg | DYN_mas_end = {:12.6e} kg'.format (DYN_mas_wat_beg, DYN_mas_wat_end) ) 
     941echo ( 'DYN_mas_air_beg = {:12.6e} kg | DYN_mas_air_end = {:12.6e} kg'.format (DYN_mas_air_beg, DYN_mas_air_end) ) 
     942echo ( 'DYN_mas_wat_beg = {:12.6e} kg | DYN_mas_wat_end = {:12.6e} kg'.format (DYN_mas_wat_beg, DYN_mas_wat_end) ) 
    840943prtFlux ( 'dMass (atm)  ', dDYN_mas_wat, 'e', True ) 
    841944 
    … …  
    850953              d_ATM_end['QS03']  *d_ATM_end['FOCE'] + d_ATM_end['QS04']  *d_ATM_end['FSIC'] 
    851954 
    852 ATM_qsol_beg = d_ATM_beg['QSOL'] 
    853 ATM_qsol_end = d_ATM_end['QSOL'] 
     955ATM_qsol_beg     = d_ATM_beg['QSOL'] 
     956ATM_qsol_end     = d_ATM_end['QSOL'] 
    854957 
    855958LIC_sno_beg      = d_ATM_beg['SNOW02']*d_ATM_beg['FLIC'] 
    … …  
    859962LIC_runlic0_end  = d_ATM_end['RUNOFFLIC0'] 
    860963 
    861 ATM_qs01_beg  = d_ATM_beg['QS01'] * d_ATM_beg['FTER'] 
    862 ATM_qs01_end  = d_ATM_end['QS01'] * d_ATM_end['FTER'] 
    863 ATM_qs02_beg  = d_ATM_beg['QS02'] * d_ATM_beg['FLIC'] 
    864 ATM_qs02_end  = d_ATM_end['QS02'] * d_ATM_end['FLIC'] 
    865 ATM_qs03_beg  = d_ATM_beg['QS03'] * d_ATM_beg['FOCE'] 
    866 ATM_qs03_end  = d_ATM_end['QS03'] * d_ATM_end['FOCE'] 
    867 ATM_qs04_beg  = d_ATM_beg['QS04'] * d_ATM_beg['FSIC'] 
    868 ATM_qs04_end  = d_ATM_end['QS04'] * d_ATM_end['FSIC']   
     964ATM_qs01_beg     = d_ATM_beg['QS01'] * d_ATM_beg['FTER'] 
     965ATM_qs02_beg     = d_ATM_beg['QS02'] * d_ATM_beg['FLIC'] 
     966ATM_qs03_beg     = d_ATM_beg['QS03'] * d_ATM_beg['FOCE'] 
     967ATM_qs04_beg     = d_ATM_beg['QS04'] * d_ATM_beg['FSIC'] 
     968 
     969ATM_qs01_end     = d_ATM_end['QS01'] * d_ATM_end['FTER'] 
     970ATM_qs02_end     = d_ATM_end['QS02'] * d_ATM_end['FLIC'] 
     971ATM_qs03_end     = d_ATM_end['QS03'] * d_ATM_end['FOCE'] 
     972ATM_qs04_end     = d_ATM_end['QS04'] * d_ATM_end['FSIC']   
    869973 
    870974if ICO : 
    871      ATM_sno_beg       = ATM_sno_beg    [ATM_beg_keysort] 
    872      ATM_sno_end       = ATM_sno_end    [ATM_end_keysort] 
    873      ATM_qs_beg        = ATM_qs_beg     [ATM_beg_keysort] 
    874      ATM_qs_end        = ATM_qs_end     [ATM_end_keysort] 
    875      ATM_qsol_beg      = ATM_qsol_beg   [ATM_beg_keysort] 
    876      ATM_qsol_end      = ATM_qsol_end   [ATM_end_keysort] 
    877      ATM_qs01_beg      = ATM_qs01_beg   [ATM_beg_keysort] 
    878      ATM_qs01_end      = ATM_qs01_end   [ATM_end_keysort] 
    879      ATM_qs02_beg      = ATM_qs02_beg   [ATM_beg_keysort] 
    880      ATM_qs02_end      = ATM_qs02_end   [ATM_end_keysort] 
    881      ATM_qs03_beg      = ATM_qs03_beg   [ATM_beg_keysort] 
    882      ATM_qs03_end      = ATM_qs03_end   [ATM_end_keysort] 
    883      ATM_qs04_beg      = ATM_qs04_beg   [ATM_beg_keysort] 
    884      ATM_qs04_end      = ATM_qs04_end   [ATM_end_keysort] 
    885      LIC_sno_beg       = LIC_sno_beg    [ATM_beg_keysort] 
    886      LIC_sno_end       = LIC_sno_end    [ATM_end_keysort] 
    887      LIC_runlic0_beg   = LIC_runlic0_beg[ATM_beg_keysort] 
    888      LIC_runlic0_end   = LIC_runlic0_end[ATM_end_keysort] 
     975     ATM_sno_beg     = ATM_sno_beg    [ATM_beg_keysort] 
     976     ATM_sno_end     = ATM_sno_end    [ATM_end_keysort] 
     977     ATM_qs_beg      = ATM_qs_beg     [ATM_beg_keysort] 
     978     ATM_qs_end      = ATM_qs_end     [ATM_end_keysort] 
     979     ATM_qsol_beg    = ATM_qsol_beg   [ATM_beg_keysort] 
     980     ATM_qs01_beg    = ATM_qs01_beg   [ATM_beg_keysort] 
     981     ATM_qs02_beg    = ATM_qs02_beg   [ATM_beg_keysort] 
     982     ATM_qs03_beg    = ATM_qs03_beg   [ATM_beg_keysort] 
     983     ATM_qs04_beg    = ATM_qs04_beg   [ATM_beg_keysort] 
     984     ATM_qsol_end    = ATM_qsol_end   [ATM_end_keysort] 
     985     ATM_qs01_end    = ATM_qs01_end   [ATM_end_keysort] 
     986     ATM_qs02_end    = ATM_qs02_end   [ATM_end_keysort] 
     987     ATM_qs03_end    = ATM_qs03_end   [ATM_end_keysort] 
     988     ATM_qs04_end    = ATM_qs04_end   [ATM_end_keysort] 
     989     LIC_sno_beg     = LIC_sno_beg    [ATM_beg_keysort] 
     990     LIC_sno_end     = LIC_sno_end    [ATM_end_keysort] 
     991     LIC_runlic0_beg = LIC_runlic0_beg[ATM_beg_keysort] 
     992     LIC_runlic0_end = LIC_runlic0_end[ATM_end_keysort] 
    889993    
    890994LIC_qs_beg = ATM_qs02_beg 
    … …  
    8971001ATM_mas_qs_end      = DYN_stock_int ( ATM_qs_end   ) 
    8981002ATM_mas_qsol_beg    = DYN_stock_int ( ATM_qsol_beg ) 
     1003ATM_mas_qs01_beg    = DYN_stock_int ( ATM_qs01_beg ) 
     1004ATM_mas_qs02_beg    = DYN_stock_int ( ATM_qs02_beg ) 
     1005ATM_mas_qs03_beg    = DYN_stock_int ( ATM_qs03_beg ) 
     1006ATM_mas_qs04_beg    = DYN_stock_int ( ATM_qs04_beg ) 
    8991007ATM_mas_qsol_end    = DYN_stock_int ( ATM_qsol_end ) 
    900 ATM_mas_qs01_beg    = DYN_stock_int ( ATM_qs01_beg ) 
    9011008ATM_mas_qs01_end    = DYN_stock_int ( ATM_qs01_end ) 
    902 ATM_mas_qs02_beg    = DYN_stock_int ( ATM_qs02_beg ) 
    9031009ATM_mas_qs02_end    = DYN_stock_int ( ATM_qs02_end ) 
    904 ATM_mas_qs03_beg    = DYN_stock_int ( ATM_qs03_beg ) 
    9051010ATM_mas_qs03_end    = DYN_stock_int ( ATM_qs03_end ) 
    906 ATM_mas_qs04_beg    = DYN_stock_int ( ATM_qs04_beg ) 
    9071011ATM_mas_qs04_end    = DYN_stock_int ( ATM_qs04_end ) 
    9081012 
    … …  
    9291033dLIC_mas_sno     = LIC_mas_sno_end     - LIC_mas_sno_beg 
    9301034dLIC_mas_runlic0 = LIC_mas_runlic0_end - LIC_mas_runlic0_beg 
    931 dLIC_mas_wat     = dLIC_mas_qs + dLIC_mas_sno 
     1035 
     1036dLIC_mas_wat     = dLIC_mas_qs + dLIC_mas_sno # + dLIC_mas_runlic0 
    9321037 
    9331038echo ( f'\nChange of atmosphere snow content (Land ice) -- {Title} ' ) 
    … …  
    9801085    RUN_mas_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   ) 
    9811086 
    982 RUN_mas_wat_beg = RUN_mas_wat_fast_beg  + RUN_mas_wat_slow_beg + RUN_mas_wat_stream_beg + \ 
    983                   RUN_mas_wat_flood_beg + RUN_mas_wat_lake_beg + RUN_mas_wat_pond_beg 
    984 RUN_mas_wat_end = RUN_mas_wat_fast_end  + RUN_mas_wat_slow_end + RUN_mas_wat_stream_end + \ 
    985                   RUN_mas_wat_flood_end + RUN_mas_wat_lake_end + RUN_mas_wat_pond_end 
     1087RUN_mas_wat_beg = Sprec ( [RUN_mas_wat_fast_beg , RUN_mas_wat_slow_beg, RUN_mas_wat_stream_beg, 
     1088                           RUN_mas_wat_flood_beg, RUN_mas_wat_lake_beg, RUN_mas_wat_pond_beg] ) 
     1089                   
     1090RUN_mas_wat_end = Sprec ( [RUN_mas_wat_fast_end  , RUN_mas_wat_slow_end , RUN_mas_wat_stream_end, 
     1091                           RUN_mas_wat_flood_end , RUN_mas_wat_lake_end , RUN_mas_wat_pond_end] ) 
    9861092 
    9871093dRUN_mas_wat_fast   = RUN_mas_wat_fast_end   - RUN_mas_wat_fast_beg 
    … …  
    10321138 
    10331139if LMDZ : 
    1034     SRF_tot_watveg_beg  = lmdz.geo2point (SRF_tot_watveg_beg , dim1D='cell_latlon') 
    1035     SRF_tot_watsoil_beg = lmdz.geo2point (SRF_tot_watsoil_beg, dim1D='cell_latlon') 
    1036     SRF_snow_beg        = lmdz.geo2point (SRF_snow_beg       , dim1D='cell_latlon') 
    1037     SRF_lakeres_beg     = lmdz.geo2point (SRF_lakeres_beg    , dim1D='cell_latlon') 
    1038     SRF_tot_watveg_end  = lmdz.geo2point (SRF_tot_watveg_end , dim1D='cell_latlon') 
    1039     SRF_tot_watsoil_end = lmdz.geo2point (SRF_tot_watsoil_end, dim1D='cell_latlon') 
    1040     SRF_snow_end        = lmdz.geo2point (SRF_snow_end       , dim1D='cell_latlon') 
    1041     SRF_lakeres_end     = lmdz.geo2point (SRF_lakeres_end    , dim1D='cell_latlon') 
     1140    SRF_tot_watveg_beg  = lmdz.geo2point (SRF_tot_watveg_beg , dim1D='cell') 
     1141    SRF_tot_watsoil_beg = lmdz.geo2point (SRF_tot_watsoil_beg, dim1D='cell') 
     1142    SRF_snow_beg        = lmdz.geo2point (SRF_snow_beg       , dim1D='cell') 
     1143    SRF_lakeres_beg     = lmdz.geo2point (SRF_lakeres_beg    , dim1D='cell') 
     1144    SRF_tot_watveg_end  = lmdz.geo2point (SRF_tot_watveg_end , dim1D='cell') 
     1145    SRF_tot_watsoil_end = lmdz.geo2point (SRF_tot_watsoil_end, dim1D='cell') 
     1146    SRF_snow_end        = lmdz.geo2point (SRF_snow_end       , dim1D='cell') 
     1147    SRF_lakeres_end     = lmdz.geo2point (SRF_lakeres_end    , dim1D='cell') 
    10421148   
    10431149if ICO : 
    … …  
    10531159# Stock dSoilHum dInterce dSWE dStream dFastR dSlowR dLake dPond dFlood 
    10541160 
    1055 SRF_wat_beg = SRF_tot_watveg_beg + SRF_tot_watsoil_beg + SRF_snow_beg  
     1161SRF_wat_beg = SRF_tot_watveg_beg + SRF_tot_watsoil_beg + SRF_snow_beg 
    10561162SRF_wat_end = SRF_tot_watveg_end + SRF_tot_watsoil_end + SRF_snow_end 
    10571163 
    … …  
    10791185print (' -- ') ; sys.stdout.flush () 
    10801186 
    1081 dSRF_mas_watveg   = SRF_mas_watveg_end  - SRF_mas_watveg_beg 
    1082 dSRF_mas_watsoil  = SRF_mas_watsoil_end - SRF_mas_watsoil_beg 
    1083 dSRF_mas_snow     = SRF_mas_snow_end    - SRF_mas_snow_beg 
    1084 dSRF_mas_lake     = SRF_mas_lake_end    - SRF_mas_lake_beg 
     1187dSRF_mas_watveg   = Sprec ( [SRF_mas_watveg_end , -SRF_mas_watveg_beg] ) 
     1188dSRF_mas_watsoil  = Sprec ( [SRF_mas_watsoil_end, -SRF_mas_watsoil_beg] ) 
     1189dSRF_mas_snow     = Sprec ( [SRF_mas_snow_end   , -SRF_mas_snow_beg] ) 
     1190dSRF_mas_lake     = Sprec ( [SRF_mas_lake_end   , -SRF_mas_lake_beg] ) 
    10851191 
    10861192echo ( '------------------------------------------------------------------------------------' ) 
    … …  
    10961202prtFlux ( 'dMass (lake)   ', dSRF_mas_lake   , 'e' , True ) 
    10971203 
    1098 SRF_mas_wat_beg = SRF_mas_watveg_beg + SRF_mas_watsoil_beg + SRF_mas_snow_beg  
    1099 SRF_mas_wat_end = SRF_mas_watveg_end + SRF_mas_watsoil_end + SRF_mas_snow_end  
    1100 dSRF_mas_wat    = SRF_mas_wat_end - SRF_mas_wat_beg 
     1204SRF_mas_wat_beg = Sprec ( [SRF_mas_watveg_beg , SRF_mas_watsoil_beg, SRF_mas_snow_beg] ) 
     1205SRF_mas_wat_end = Sprec ( [SRF_mas_watveg_end , SRF_mas_watsoil_end, SRF_mas_snow_end] ) 
     1206 
     1207dSRF_mas_wat    = Sprec ( [+SRF_mas_watveg_end , +SRF_mas_watsoil_end, +SRF_mas_snow_end,  
     1208                           -SRF_mas_watveg_beg , -SRF_mas_watsoil_beg, -SRF_mas_snow_beg]  ) 
    11011209 
    11021210echo ( '------------------------------------------------------------------------------------' ) 
    … …  
    11151223echo ( f'-- ATM Fluxes  -- {Title} ' ) 
    11161224 
    1117 Step=1 
    1118 print ( 'Step {Step}', end='' ) ; Step += 1 
    11191225if ATM_HIS == 'latlon' : 
    1120     ATM_wbilo_oce   = lmdz.geo2point ( d_ATM_his ['wbilo_oce'], dim1D='cell_latlon' ) 
    1121     ATM_wbilo_sic   = lmdz.geo2point ( d_ATM_his ['wbilo_sic'], dim1D='cell_latlon' ) 
    1122     ATM_wbilo_ter   = lmdz.geo2point ( d_ATM_his ['wbilo_ter'], dim1D='cell_latlon' ) 
    1123     ATM_wbilo_lic   = lmdz.geo2point ( d_ATM_his ['wbilo_lic'], dim1D='cell_latlon' ) 
    1124     ATM_runofflic   = lmdz.geo2point ( d_ATM_his ['runofflic'], dim1D='cell_latlon' ) 
    1125     ATM_fqcalving   = lmdz.geo2point ( d_ATM_his ['fqcalving'], dim1D='cell_latlon' ) 
    1126     ATM_fqfonte     = lmdz.geo2point ( d_ATM_his ['fqfonte']  , dim1D='cell_latlon' ) 
    1127     ATM_precip      = lmdz.geo2point ( d_ATM_his ['precip']   , dim1D='cell_latlon' ) 
    1128     ATM_snowf       = lmdz.geo2point ( d_ATM_his ['snow']     , dim1D='cell_latlon' ) 
    1129     ATM_evap        = lmdz.geo2point ( d_ATM_his ['evap']     , dim1D='cell_latlon' ) 
    1130     ATM_wevap_ter   = lmdz.geo2point ( d_ATM_his ['wevap_ter'], dim1D='cell_latlon' ) 
    1131     ATM_wevap_oce   = lmdz.geo2point ( d_ATM_his ['wevap_oce'], dim1D='cell_latlon' ) 
    1132     ATM_wevap_lic   = lmdz.geo2point ( d_ATM_his ['wevap_lic'], dim1D='cell_latlon' ) 
    1133     ATM_wevap_sic   = lmdz.geo2point ( d_ATM_his ['wevap_sic'], dim1D='cell_latlon' ) 
    1134  
    1135     ATM_wrain_ter   = lmdz.geo2point ( d_ATM_his ['wrain_ter'], dim1D='cell_latlon' ) 
    1136     ATM_wrain_oce   = lmdz.geo2point ( d_ATM_his ['wrain_oce'], dim1D='cell_latlon' ) 
    1137     ATM_wrain_lic   = lmdz.geo2point ( d_ATM_his ['wrain_lic'], dim1D='cell_latlon' ) 
    1138     ATM_wrain_sic   = lmdz.geo2point ( d_ATM_his ['wrain_sic'], dim1D='cell_latlon' ) 
    1139  
    1140     ATM_wsnow_ter   = lmdz.geo2point ( d_ATM_his ['wsnow_ter'], dim1D='cell_latlon' ) 
    1141     ATM_wsnow_oce   = lmdz.geo2point ( d_ATM_his ['wsnow_oce'], dim1D='cell_latlon' ) 
    1142     ATM_wsnow_lic   = lmdz.geo2point ( d_ATM_his ['wsnow_lic'], dim1D='cell_latlon' ) 
    1143     ATM_wsnow_sic   = lmdz.geo2point ( d_ATM_his ['wsnow_sic'], dim1D='cell_latlon' ) 
    1144  
    1145     ATM_wevap_ter   = lmdz.geo2point ( d_ATM_his ['wevap_ter'], dim1D='cell_latlon' ) 
    1146     ATM_wevap_oce   = lmdz.geo2point ( d_ATM_his ['wevap_oce'], dim1D='cell_latlon' ) 
    1147     ATM_wevap_lic   = lmdz.geo2point ( d_ATM_his ['wevap_lic'], dim1D='cell_latlon' ) 
    1148     ATM_wevap_sic   = lmdz.geo2point ( d_ATM_his ['wevap_sic'], dim1D='cell_latlon' ) 
    1149      
    1150     ATM_runofflic   = lmdz.geo2point ( d_ATM_his ['runofflic'], dim1D='cell_latlon' ) 
    1151      
    1152 print ( 'Step {Step}', end='' ) ; Step += 1 
     1226    echo ( ' latlon case' ) 
     1227    ATM_wbilo_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_oce']), dim1D='cell' ) 
     1228    ATM_wbilo_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_sic']), dim1D='cell' ) 
     1229    ATM_wbilo_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_ter']), dim1D='cell' ) 
     1230    ATM_wbilo_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_lic']), dim1D='cell' ) 
     1231    ATM_runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1D='cell' ) 
     1232    ATM_fqcalving   = lmdz.geo2point ( rprec (d_ATM_his ['fqcalving']), dim1D='cell' ) 
     1233    ATM_fqfonte     = lmdz.geo2point ( rprec (d_ATM_his ['fqfonte']  ), dim1D='cell' ) 
     1234    ATM_precip      = lmdz.geo2point ( rprec (d_ATM_his ['precip']   ), dim1D='cell' ) 
     1235    ATM_snowf       = lmdz.geo2point ( rprec (d_ATM_his ['snow']     ), dim1D='cell' ) 
     1236    ATM_evap        = lmdz.geo2point ( rprec (d_ATM_his ['evap']     ), dim1D='cell' ) 
     1237    ATM_wevap_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_ter']), dim1D='cell' ) 
     1238    ATM_wevap_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_oce']), dim1D='cell' ) 
     1239    ATM_wevap_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_lic']), dim1D='cell' ) 
     1240    ATM_wevap_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_sic']), dim1D='cell' ) 
     1241    ATM_wrain_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_ter']), dim1D='cell' ) 
     1242    ATM_wrain_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_oce']), dim1D='cell' ) 
     1243    ATM_wrain_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_lic']), dim1D='cell' ) 
     1244    ATM_wrain_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_sic']), dim1D='cell' ) 
     1245    ATM_wsnow_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_ter']), dim1D='cell' ) 
     1246    ATM_wsnow_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_oce']), dim1D='cell' ) 
     1247    ATM_wsnow_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_lic']), dim1D='cell' ) 
     1248    ATM_wsnow_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_sic']), dim1D='cell' ) 
     1249    ATM_runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1D='cell' ) 
     1250    echo ( f'End of LATLON case') 
     1251     
    11531252if ATM_HIS == 'ico' : 
    1154     ATM_wbilo_oce   = d_ATM_his ['wbilo_oce'][..., ATM_his_keysort] 
    1155     ATM_wbilo_sic   = d_ATM_his ['wbilo_sic'][..., ATM_his_keysort] 
    1156     ATM_wbilo_ter   = d_ATM_his ['wbilo_ter'][..., ATM_his_keysort] 
    1157     ATM_wbilo_lic   = d_ATM_his ['wbilo_lic'][..., ATM_his_keysort] 
    1158     ATM_runofflic   = d_ATM_his ['runofflic'][..., ATM_his_keysort] 
    1159     ATM_fqcalving   = d_ATM_his ['fqcalving'][..., ATM_his_keysort] 
    1160     ATM_fqfonte     = d_ATM_his ['fqfonte']  [..., ATM_his_keysort] 
    1161     ATM_precip      = d_ATM_his ['precip']   [..., ATM_his_keysort] 
    1162     ATM_snowf       = d_ATM_his ['snow']     [..., ATM_his_keysort] 
    1163     ATM_evap        = d_ATM_his ['evap']     [..., ATM_his_keysort] 
    1164     ATM_wevap_ter   = d_ATM_his ['wevap_ter'][..., ATM_his_keysort] 
    1165     ATM_wevap_oce   = d_ATM_his ['wevap_oce'][..., ATM_his_keysort] 
    1166     ATM_wevap_lic   = d_ATM_his ['wevap_lic'][..., ATM_his_keysort] 
    1167     ATM_wevap_sic   = d_ATM_his ['wevap_sic'][..., ATM_his_keysort] 
    1168     ATM_runofflic   = d_ATM_his ['runofflic'][..., ATM_his_keysort] 
    1169  
    1170     ATM_wevap_ter   = d_ATM_his ['wevap_ter'][..., ATM_his_keysort] 
    1171     ATM_wevap_oce   = d_ATM_his ['wevap_oce'][..., ATM_his_keysort] 
    1172     ATM_wevap_lic   = d_ATM_his ['wevap_lic'][..., ATM_his_keysort] 
    1173     ATM_wevap_sic   = d_ATM_his ['wevap_sic'][..., ATM_his_keysort] 
    1174  
    1175     ATM_wrain_ter   = d_ATM_his ['wrain_ter'][..., ATM_his_keysort] 
    1176     ATM_wrain_oce   = d_ATM_his ['wrain_oce'][..., ATM_his_keysort] 
    1177     ATM_wrain_lic   = d_ATM_his ['wrain_lic'][..., ATM_his_keysort] 
    1178     ATM_wrain_sic   = d_ATM_his ['wrain_sic'][..., ATM_his_keysort] 
    1179  
    1180     ATM_wsnow_ter   = d_ATM_his ['wsnow_ter'][..., ATM_his_keysort] 
    1181     ATM_wsnow_oce   = d_ATM_his ['wsnow_oce'][..., ATM_his_keysort] 
    1182     ATM_wsnow_lic   = d_ATM_his ['wsnow_lic'][..., ATM_his_keysort] 
    1183     ATM_wsnow_sic   = d_ATM_his ['wsnow_sic'][..., ATM_his_keysort] 
    1184  
    1185 print ( 'Step {Step}', end='' ) ; Step += 1 
    1186  
     1253    echo (' ico case') 
     1254    ATM_wbilo_oce   = rprec (d_ATM_his ['wbilo_oce'])[..., ATM_his_keysort] 
     1255    ATM_wbilo_sic   = rprec (d_ATM_his ['wbilo_sic'])[..., ATM_his_keysort] 
     1256    ATM_wbilo_ter   = rprec (d_ATM_his ['wbilo_ter'])[..., ATM_his_keysort] 
     1257    ATM_wbilo_lic   = rprec (d_ATM_his ['wbilo_lic'])[..., ATM_his_keysort] 
     1258    ATM_runofflic   = rprec (d_ATM_his ['runofflic'])[..., ATM_his_keysort] 
     1259    ATM_fqcalving   = rprec (d_ATM_his ['fqcalving'])[..., ATM_his_keysort] 
     1260    ATM_fqfonte     = rprec (d_ATM_his ['fqfonte']  )[..., ATM_his_keysort] 
     1261    ATM_precip      = rprec (d_ATM_his ['precip']   )[..., ATM_his_keysort] 
     1262    ATM_snowf       = rprec (d_ATM_his ['snow']     )[..., ATM_his_keysort] 
     1263    ATM_evap        = rprec (d_ATM_his ['evap']     )[..., ATM_his_keysort] 
     1264    ATM_wevap_ter   = rprec (d_ATM_his ['wevap_ter'])[..., ATM_his_keysort] 
     1265    ATM_wevap_oce   = rprec (d_ATM_his ['wevap_oce'])[..., ATM_his_keysort] 
     1266    ATM_wevap_lic   = rprec (d_ATM_his ['wevap_lic'])[..., ATM_his_keysort] 
     1267    ATM_wevap_sic   = rprec (d_ATM_his ['wevap_sic'])[..., ATM_his_keysort] 
     1268    ATM_runofflic   = rprec (d_ATM_his ['runofflic'])[..., ATM_his_keysort] 
     1269    ATM_wevap_ter   = rprec (d_ATM_his ['wevap_ter'])[..., ATM_his_keysort] 
     1270    ATM_wevap_oce   = rprec (d_ATM_his ['wevap_oce'])[..., ATM_his_keysort] 
     1271    ATM_wevap_lic   = rprec (d_ATM_his ['wevap_lic'])[..., ATM_his_keysort] 
     1272    ATM_wevap_sic   = rprec (d_ATM_his ['wevap_sic'])[..., ATM_his_keysort] 
     1273    ATM_wrain_ter   = rprec (d_ATM_his ['wrain_ter'])[..., ATM_his_keysort] 
     1274    ATM_wrain_oce   = rprec (d_ATM_his ['wrain_oce'])[..., ATM_his_keysort] 
     1275    ATM_wrain_lic   = rprec (d_ATM_his ['wrain_lic'])[..., ATM_his_keysort] 
     1276    ATM_wrain_sic   = rprec (d_ATM_his ['wrain_sic'])[..., ATM_his_keysort] 
     1277    ATM_wsnow_ter   = rprec (d_ATM_his ['wsnow_ter'])[..., ATM_his_keysort] 
     1278    ATM_wsnow_oce   = rprec (d_ATM_his ['wsnow_oce'])[..., ATM_his_keysort] 
     1279    ATM_wsnow_lic   = rprec (d_ATM_his ['wsnow_lic'])[..., ATM_his_keysort] 
     1280    ATM_wsnow_sic   = rprec (d_ATM_his ['wsnow_sic'])[..., ATM_his_keysort] 
     1281    echo ( f'End of ico case ') 
     1282 
     1283echo ( 'ATM wprecip_oce' ) 
    11871284ATM_wprecip_oce = ATM_wrain_oce + ATM_wsnow_oce 
    11881285ATM_wprecip_ter = ATM_wrain_ter + ATM_wsnow_ter 
    … …  
    11901287ATM_wprecip_lic = ATM_wrain_lic + ATM_wsnow_lic 
    11911288 
    1192 ATM_wbilo       = ATM_wbilo_oce + ATM_wbilo_sic + ATM_wbilo_ter + ATM_wbilo_lic 
    1193 ATM_emp         = ATM_evap - ATM_precip 
     1289ATM_wbilo       = ATM_wbilo_oce   + ATM_wbilo_sic   + ATM_wbilo_ter   + ATM_wbilo_lic 
     1290ATM_wevap       = ATM_wevap_oce   + ATM_wevap_sic   + ATM_wevap_ter   + ATM_wevap_lic 
     1291ATM_wprecip     = ATM_wprecip_oce + ATM_wprecip_sic + ATM_wprecip_ter + ATM_wprecip_lic 
     1292ATM_wsnow       = ATM_wsnow_oce   + ATM_wsnow_sic   + ATM_wsnow_ter   + ATM_wsnow_lic 
     1293ATM_wrain       = ATM_wrain_oce   + ATM_wrain_sic   + ATM_wrain_ter   + ATM_wrain_lic 
     1294ATM_wemp        = ATM_wevap - ATM_wprecip 
     1295ATM_emp         = ATM_evap  - ATM_precip 
     1296 
    11941297ATM_wprecip_sea = ATM_wprecip_oce + ATM_wprecip_sic 
    11951298ATM_wsnow_sea   = ATM_wsnow_oce   + ATM_wsnow_sic 
    … …  
    11981301ATM_wevap_sea   = ATM_wevap_sic   + ATM_wevap_oce 
    11991302 
    1200 #ATM_wemp_ter = ATM_wevap_ter - ATM_precip * ATM_fter 
    1201 #ATM_wemp_oce = ATM_wevap_oce - ATM_precip * ATM_foce 
    1202 #ATM_wemp_sic = ATM_wevap_sic - ATM_precip * ATM_fsic 
    1203 #ATM_wemp_lic = ATM_wevap_lic - ATM_precip * ATM_flic 
    1204 #ATM_wemp_sea = ATM_wevap_sic + ATM_wevap_oce - ATM_precip * ATM_fsea 
    1205  
    1206 ATM_wemp_ter = ATM_wevap_ter - ATM_wprecip_ter  
    1207 ATM_wemp_oce = ATM_wevap_oce - ATM_wprecip_oce 
    1208 ATM_wemp_sic = ATM_wevap_sic - ATM_wprecip_sic 
    1209 ATM_wemp_lic = ATM_wevap_lic - ATM_wprecip_lic 
    1210 ATM_wemp_sea = ATM_wevap_sic - ATM_wprecip_oce 
    1211  
    1212 print ( 'Step {Step}', end='' ) ; Step += 1 
    1213 if RUN_HIS == 'latlon' :  
    1214     RUN_coastalflow = lmdz.geo2point ( d_RUN_his ['coastalflow'], dim1D='cell_latlon' )  
    1215     RUN_riverflow   = lmdz.geo2point ( d_RUN_his ['riverflow']  , dim1D='cell_latlon' ) 
    1216     RUN_runoff      = lmdz.geo2point ( d_RUN_his ['runoff']     , dim1D='cell_latlon' ) 
    1217     RUN_drainage    = lmdz.geo2point ( d_RUN_his ['drainage']   , dim1D='cell_latlon' ) 
    1218     RUN_riversret   = lmdz.geo2point ( d_RUN_his ['riversret']  , dim1D='cell_latlon' ) 
    1219      
    1220     RUN_coastalflow_cpl = lmdz.geo2point ( d_RUN_his ['coastalflow_cpl'], dim1D='cell_latlon' )  
    1221     RUN_riverflow_cpl   = lmdz.geo2point ( d_RUN_his ['riverflow_cpl']  , dim1D='cell_latlon' ) 
    1222  
    1223 print ( 'Step {Step}', end='' ) ; Step += 1 
     1303ATM_wemp_ter    = ATM_wevap_ter - ATM_wprecip_ter 
     1304ATM_wemp_oce    = ATM_wevap_oce - ATM_wprecip_oce 
     1305ATM_wemp_sic    = ATM_wevap_sic - ATM_wprecip_sic 
     1306ATM_wemp_lic    = ATM_wevap_lic - ATM_wprecip_lic 
     1307ATM_wemp_sea    = ATM_wevap_sic - ATM_wprecip_oce 
     1308 
     1309if RUN_HIS == 'latlon' : 
     1310    echo ( f'RUN costalflow Grille LATLON' ) 
     1311    if TestInterp : 
     1312         echo ( f'RUN runoff TestInterp' ) 
     1313         RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff_contfrac_interp']  )   , dim1D='cell' ) 
     1314         RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage_contfrac_interp'])   , dim1D='cell' ) 
     1315    else :  
     1316        echo ( f'RUN runoff' ) 
     1317        RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff']         ), dim1D='cell' ) 
     1318        RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage']       ), dim1D='cell' ) 
     1319         
     1320    RUN_coastalflow     = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow']    ), dim1D='cell' )  
     1321    RUN_riverflow       = lmdz.geo2point ( rprec (d_RUN_his ['riverflow']      ), dim1D='cell' ) 
     1322    RUN_riversret       = lmdz.geo2point ( rprec (d_RUN_his ['riversret']      ), dim1D='cell' ) 
     1323    RUN_coastalflow_cpl = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow_cpl']), dim1D='cell' )  
     1324    RUN_riverflow_cpl   = lmdz.geo2point ( rprec (d_RUN_his ['riverflow_cpl']  ), dim1D='cell' ) 
     1325 
    12241326if RUN_HIS == 'ico' : 
    1225     RUN_coastalflow =  d_RUN_his ['coastalflow'][..., RUN_his_keysort] 
    1226     RUN_riverflow   =  d_RUN_his ['riverflow']  [..., RUN_his_keysort] 
    1227     RUN_runoff      =  d_RUN_his ['runoff']     [..., RUN_his_keysort] 
    1228     RUN_drainage    =  d_RUN_his ['drainage']   [..., RUN_his_keysort] 
    1229     RUN_riversret   =  d_RUN_his ['riversret']  [..., RUN_his_keysort] 
    1230      
    1231     RUN_coastalflow_cpl = d_RUN_his ['coastalflow_cpl'][..., RUN_his_keysort] 
    1232     RUN_riverflow_cpl   = d_RUN_his ['riverflow_cpl']  [..., RUN_his_keysort] 
    1233  
    1234 print ( 'Step {Step}', end='' ) ; Step += 1 
     1327    echo ( f'RUN costalflow Grille ICO' ) 
     1328    RUN_coastalflow =  rprec (d_RUN_his ['coastalflow'])[..., RUN_his_keysort] 
     1329    RUN_riverflow   =  rprec (d_RUN_his ['riverflow']  )[..., RUN_his_keysort] 
     1330    RUN_runoff      =  rprec (d_RUN_his ['runoff']     )[..., RUN_his_keysort] 
     1331    RUN_drainage    =  rprec (d_RUN_his ['drainage']   )[..., RUN_his_keysort] 
     1332    RUN_riversret   =  rprec (d_RUN_his ['riversret']  )[..., RUN_his_keysort] 
     1333     
     1334    RUN_coastalflow_cpl = rprec (d_RUN_his ['coastalflow_cpl'])[..., RUN_his_keysort] 
     1335    RUN_riverflow_cpl   = rprec (d_RUN_his ['riverflow_cpl']  )[..., RUN_his_keysort] 
     1336 
     1337Step = 0 
     1338 
    12351339if SRF_HIS == 'latlon' : 
    1236     SRF_rain     = lmdz.geo2point ( d_SRF_his ['rain']  , dim1D='cell_latlon') 
    1237     SRF_evap     = lmdz.geo2point ( d_SRF_his ['evap']  , dim1D='cell_latlon') 
    1238     SRF_snowf    = lmdz.geo2point ( d_SRF_his ['snowf'] , dim1D='cell_latlon') 
    1239     SRF_subli    = lmdz.geo2point ( d_SRF_his ['subli'] , dim1D='cell_latlon') 
    1240     SRF_transpir = lmdz.geo2point ( np.sum (d_SRF_his ['transpir'], axis=1), dim1D='cell_latlon' ) 
    1241  
    1242 print ( 'Step {Step}', end='' ) ; Step += 1 
     1340    if TestInterp : 
     1341         echo ( f'SRF rain TestInterp' ) 
     1342         SRF_rain     = lmdz.geo2point ( rprec (d_SRF_his ['rain_contfrac_interp'] ), dim1D='cell') 
     1343         SRF_evap     = lmdz.geo2point ( rprec (d_SRF_his ['evap_contfrac_interp'] ), dim1D='cell') 
     1344         SRF_snowf    = lmdz.geo2point ( rprec (d_SRF_his ['snow_contfrac_interp'] ), dim1D='cell') 
     1345         SRF_subli    = lmdz.geo2point ( rprec (d_SRF_his ['subli_contfrac_interp']), dim1D='cell') 
     1346         SRF_transpir = lmdz.geo2point ( rprec (d_SRF_his ['transpir_contfrac_interp']).sum(dim='veget'), dim1D='cell' ) 
     1347         SRF_rain.attrs     = d_SRF_his ['rain_contfrac_interp'].attrs 
     1348         SRF_evap.attrs     = d_SRF_his ['evap_contfrac_interp'].attrs 
     1349         SRF_snowf.attrs    = d_SRF_his ['snow_contfrac_interp'].attrs 
     1350         SRF_subli.attrs    = d_SRF_his ['subli_contfrac_interp'].attrs 
     1351         SRF_transpir.attrs = d_SRF_his ['transpir_contfrac_interp'].attrs 
     1352    else :  
     1353        echo ( f'SRF rain' ) 
     1354        SRF_rain     = lmdz.geo2point ( rprec (d_SRF_his ['rain'] ) , dim1D='cell') 
     1355        SRF_evap     = lmdz.geo2point ( rprec (d_SRF_his ['evap'] ) , dim1D='cell') 
     1356        SRF_snowf    = lmdz.geo2point ( rprec (d_SRF_his ['snowf']) , dim1D='cell') 
     1357        SRF_subli    = lmdz.geo2point ( rprec (d_SRF_his ['subli']) , dim1D='cell') 
     1358        SRF_transpir = lmdz.geo2point ( rprec (d_SRF_his ['transpir']).sum(dim='veget'), dim1D='cell' ) 
     1359 
    12431360if SRF_HIS == 'ico' : 
    1244     SRF_rain     = d_SRF_his ['rain'] [..., SRF_his_keysort] 
    1245     SRF_evap     = d_SRF_his ['evap'] [..., SRF_his_keysort] 
    1246     SRF_snowf    = d_SRF_his ['snowf'][..., SRF_his_keysort] 
    1247     SRF_subli    = d_SRF_his ['subli'][..., SRF_his_keysort] 
    1248     SRF_transpir = np.sum (d_SRF_his ['transpir'], axis=1)[..., SRF_his_keysort] 
    1249  
    1250 print ( 'Step {Step}', end='' ) ; Step += 1 
     1361    echo ( f'SRF rain') 
     1362    SRF_rain     = rprec (d_SRF_his ['rain'] )[..., SRF_his_keysort] 
     1363    SRF_evap     = rprec (d_SRF_his ['evap'] )[..., SRF_his_keysort] 
     1364    SRF_snowf    = rprec (d_SRF_his ['snowf'])[..., SRF_his_keysort] 
     1365    SRF_subli    = rprec (d_SRF_his ['subli'])[..., SRF_his_keysort] 
     1366    SRF_transpir = rprec (d_SRF_his ['transpir']).sum(dim='veget')[..., SRF_his_keysort] 
     1367 
     1368echo ( f'SRF emp' ) 
    12511369SRF_transpir.attrs['units'] = d_SRF_his ['transpir'].attrs['units'] 
    12521370SRF_emp      = SRF_evap - SRF_rain - SRF_snowf ; SRF_emp.attrs['units'] = SRF_rain.attrs['units'] 
    … …  
    12711389    mmd2SI (VarT) 
    12721390 
    1273 print ( 'Step {Step}', end='' ) ; Step += 1 
     1391echo ( f'RUN input' ) 
    12741392RUN_input  = RUN_runoff      + RUN_drainage 
    12751393RUN_output = RUN_coastalflow + RUN_riverflow 
    12761394 
    1277 print ( 'Step {Step}', end='' ) ; Step += 1 
     1395echo ( f'ATM flw_wbilo' ) 
    12781396ATM_flx_wbilo       = ATM_flux_int ( ATM_wbilo      ) 
     1397ATM_flx_wevap       = ATM_flux_int ( ATM_wevap      ) 
     1398ATM_flx_wprecip     = ATM_flux_int ( ATM_wprecip    ) 
     1399ATM_flx_wsnow       = ATM_flux_int ( ATM_wsnow      ) 
     1400ATM_flx_wrain       = ATM_flux_int ( ATM_wrain      ) 
     1401ATM_flx_wemp        = ATM_flux_int ( ATM_wemp       ) 
     1402 
    12791403ATM_flx_wbilo_lic   = ATM_flux_int ( ATM_wbilo_lic  ) 
    12801404ATM_flx_wbilo_oce   = ATM_flux_int ( ATM_wbilo_oce  ) 
    … …  
    12851409ATM_flx_fqfonte     = ATM_flux_int ( ATM_fqfonte    ) 
    12861410 
    1287 print ( 'Step {Step}', end='' ) ; Step += 1 
     1411LIC_flx_calving     = LIC_flux_int ( ATM_fqcalving  ) 
     1412LIC_flx_fqfonte     = LIC_flux_int ( ATM_fqfonte    ) 
     1413 
     1414echo ( f'ATM flx precip' ) 
    12881415ATM_flx_precip      = ATM_flux_int ( ATM_precip     ) 
    12891416ATM_flx_snowf       = ATM_flux_int ( ATM_snowf      ) 
    … …  
    12911418ATM_flx_runlic      = ATM_flux_int ( ATM_runofflic  ) 
    12921419 
    1293 print ( 'Step {Step}', end='' ) ; Step += 1 
     1420LIC_flx_precip      = LIC_flux_int ( ATM_precip     ) 
     1421LIC_flx_snowf       = LIC_flux_int ( ATM_snowf      ) 
     1422LIC_flx_evap        = LIC_flux_int ( ATM_evap       ) 
     1423LIC_flx_runlic      = LIC_flux_int ( ATM_runofflic  ) 
     1424 
     1425echo ( f'ATM flx_wrain_ter' ) 
    12941426ATM_flx_wrain_ter    = ATM_flux_int ( ATM_wrain_ter ) 
    12951427ATM_flx_wrain_oce    = ATM_flux_int ( ATM_wrain_oce ) 
    … …  
    13031435ATM_flx_wsnow_sic    = ATM_flux_int ( ATM_wsnow_sic ) 
    13041436ATM_flx_wsnow_sea    = ATM_flux_int ( ATM_wsnow_sea ) 
    1305 print ( 'Step {Step}', end='' ) ; Step += 1 
    1306 ATM_flx_wrain_ter    = ATM_flux_int ( ATM_wrain_ter ) 
    1307 ATM_flx_wrain_oce    = ATM_flux_int ( ATM_wrain_oce ) 
    1308 ATM_flx_wrain_lic    = ATM_flux_int ( ATM_wrain_lic ) 
    1309 ATM_flx_wrain_sic    = ATM_flux_int ( ATM_wrain_sic ) 
    1310 ATM_flx_wrain_sea    = ATM_flux_int ( ATM_wrain_sea ) 
     1437 
     1438echo ( f'ATM flx_evap_ter' )  
    13111439ATM_flx_wevap_ter    = ATM_flux_int ( ATM_wevap_ter ) 
    13121440ATM_flx_wevap_oce    = ATM_flux_int ( ATM_wevap_oce ) 
    … …  
    13251453ATM_flx_wemp_sea     = ATM_flux_int ( ATM_wemp_sea ) 
    13261454 
    1327 ATM_flx_emp          = ATM_flux_int ( ATM_emp) 
    1328  
    1329 print ( 'Step {Step}', end='' ) ; Step += 1 
     1455ATM_flx_emp          = ATM_flux_int ( ATM_emp ) 
     1456 
     1457echo ( f'RUN flx_coastal' ) 
    13301458RUN_flx_coastal     = ONE_flux_int ( RUN_coastalflow) 
     1459echo ( f'RUN flx_river' ) 
    13311460RUN_flx_river       = ONE_flux_int ( RUN_riverflow  ) 
     1461echo ( f'RUN flx_coastal_cpl' ) 
    13321462RUN_flx_coastal_cpl = ONE_flux_int ( RUN_coastalflow_cpl) 
     1463echo ( f'RUN flx_river_cpl' ) 
    13331464RUN_flx_river_cpl   = ONE_flux_int ( RUN_riverflow_cpl  ) 
     1465echo ( f'RUN flx_drainage' ) 
    13341466RUN_flx_drainage    = SRF_flux_int ( RUN_drainage   ) 
     1467echo ( f'RUN flx_riversset' ) 
    13351468RUN_flx_riversret   = SRF_flux_int ( RUN_riversret  ) 
     1469echo ( f'RUN flx_runoff' ) 
    13361470RUN_flx_runoff      = SRF_flux_int ( RUN_runoff     ) 
     1471echo ( f'RUN flx_input' ) 
    13371472RUN_flx_input       = SRF_flux_int ( RUN_input      ) 
     1473echo ( f'RUN flx_output' ) 
    13381474RUN_flx_output      = ONE_flux_int ( RUN_output     ) 
    13391475 
    1340 print ( 'Step {Step}', end='' ) ; Step += 1 
    1341 RUN_flx_bil    = RUN_flx_input   - RUN_flx_output 
    1342 RUN_flx_rivcoa = RUN_flx_coastal + RUN_flx_river 
    1343  
    1344 prtFlux ('wbilo_oce     ', ATM_flx_wbilo_oce  , 'f' )           
    1345 prtFlux ('wbilo_sic     ', ATM_flx_wbilo_sic  , 'f' )           
    1346 prtFlux ('wbilo_sic+oce ', ATM_flx_wbilo_sea  , 'f' )           
    1347 prtFlux ('wbilo_ter     ', ATM_flx_wbilo_ter  , 'f' )           
    1348 prtFlux ('wbilo_lic     ', ATM_flx_wbilo_lic  , 'f' )           
    1349 prtFlux ('Sum wbilo_*   ', ATM_flx_wbilo      , 'f', True)   
    1350 prtFlux ('E-P           ', ATM_flx_emp        , 'f', True)   
    1351 prtFlux ('calving       ', ATM_flx_calving    , 'f' )  
    1352 prtFlux ('fqfonte       ', ATM_flx_fqfonte    , 'f' )        
    1353 prtFlux ('precip        ', ATM_flx_precip     , 'f' )        
    1354 prtFlux ('snowf         ', ATM_flx_snowf      , 'f' )         
    1355 prtFlux ('evap          ', ATM_flx_evap       , 'f' ) 
    1356 prtFlux ('runoff lic    ', ATM_flx_runlic     , 'f' ) 
    1357  
    1358 prtFlux ('ATM_flx_wemp_lic     ', ATM_flx_wemp_lic      , 'f' ) 
    1359 prtFlux ('ATM_flx_wemp_oce     ', ATM_flx_wemp_oce      , 'f' ) 
    1360 prtFlux ('ATM_flx_wemp_sic     ', ATM_flx_wemp_sic      , 'f' ) 
    1361 prtFlux ('ATM_flx_wemp_ter     ', ATM_flx_wemp_ter      , 'f' ) 
    1362  
    1363 prtFlux ('ATM_flx_wprecip_lic  ', ATM_flx_wprecip_lic   , 'f' ) 
    1364 prtFlux ('ATM_flx_wprecip_oce  ', ATM_flx_wprecip_oce   , 'f' ) 
    1365 prtFlux ('ATM_flx_wprecip_sic  ', ATM_flx_wprecip_sic   , 'f' ) 
    1366 prtFlux ('ATM_flx_wprecip_ter  ', ATM_flx_wprecip_ter   , 'f' ) 
    1367  
    1368 prtFlux ('ATM_flx_wrain_lic    ', ATM_flx_wrain_lic     , 'f' ) 
    1369 prtFlux ('ATM_flx_wrain_oce    ', ATM_flx_wrain_oce     , 'f' ) 
    1370 prtFlux ('ATM_flx_wrain_sic    ', ATM_flx_wrain_sic     , 'f' ) 
    1371 prtFlux ('ATM_flx_wrain_ter    ', ATM_flx_wrain_ter     , 'f' ) 
    1372  
    1373 prtFlux ('ATM_flx_wsnow_lic    ', ATM_flx_wsnow_lic     , 'f' ) 
    1374 prtFlux ('ATM_flx_wsnow_oce    ', ATM_flx_wsnow_oce     , 'f' ) 
    1375 prtFlux ('ATM_flx_wsnow_sic    ', ATM_flx_wsnow_sic     , 'f' ) 
    1376 prtFlux ('ATM_flx_wsnow_ter    ', ATM_flx_wsnow_ter     , 'f' ) 
    1377  
    1378 prtFlux ('ATM_flx_wevap_lic    ', ATM_flx_wevap_lic     , 'f' ) 
    1379 prtFlux ('ATM_flx_wevap_oce    ', ATM_flx_wevap_oce     , 'f' ) 
    1380 prtFlux ('ATM_flx_wevap_sic    ', ATM_flx_wevap_sic     , 'f' ) 
    1381 prtFlux ('ATM_flx_wevap_ter    ', ATM_flx_wevap_ter     , 'f' ) 
    1382  
    1383 prtFlux ('ATM_flx_wevap_lic    ', ATM_flx_wevap_lic     , 'f' ) 
    1384 prtFlux ('ATM_flx_wevap_oce    ', ATM_flx_wevap_oce     , 'f' ) 
    1385 prtFlux ('ATM_flx_wevap_sic    ', ATM_flx_wevap_sic     , 'f' ) 
    1386 prtFlux ('ATM_flx_wevap_ter    ', ATM_flx_wevap_ter     , 'f' ) 
     1476echo ( f'RUN flx_bil' ) ; Step += 1 
     1477#RUN_flx_bil    = RUN_flx_input   - RUN_flx_output 
     1478#RUN_flx_rivcoa = RUN_flx_coastal + RUN_flx_river 
     1479 
     1480RUN_flx_bil    = ONE_flux_int ( RUN_input       - RUN_output) 
     1481RUN_flx_rivcoa = ONE_flux_int ( RUN_coastalflow + RUN_riverflow) 
     1482 
     1483prtFlux ('wbilo_oce            ', ATM_flx_wbilo_oce     , 'f' )           
     1484prtFlux ('wbilo_sic            ', ATM_flx_wbilo_sic     , 'f' )           
     1485prtFlux ('wbilo_sic+oce        ', ATM_flx_wbilo_sea     , 'f' )           
     1486prtFlux ('wbilo_ter            ', ATM_flx_wbilo_ter     , 'f' )           
     1487prtFlux ('wbilo_lic            ', ATM_flx_wbilo_lic     , 'f' )           
     1488prtFlux ('Sum wbilo_*          ', ATM_flx_wbilo         , 'f', True)   
     1489prtFlux ('E-P                  ', ATM_flx_emp           , 'f', True)   
     1490prtFlux ('calving              ', ATM_flx_calving       , 'f' )  
     1491prtFlux ('fqfonte              ', ATM_flx_fqfonte       , 'f' )        
     1492prtFlux ('precip               ', ATM_flx_precip        , 'f' )        
     1493prtFlux ('snowf                ', ATM_flx_snowf         , 'f' )         
     1494prtFlux ('evap                 ', ATM_flx_evap          , 'f' ) 
     1495prtFlux ('runoff lic           ', ATM_flx_runlic        , 'f' ) 
     1496 
     1497prtFlux ('ATM_flx_wevap*       ', ATM_flx_wevap         , 'f' ) 
     1498prtFlux ('ATM_flx_wrain*       ', ATM_flx_wrain         , 'f' ) 
     1499prtFlux ('ATM_flx_wsnow*       ', ATM_flx_wsnow         , 'f' ) 
     1500prtFlux ('ATM_flx_wprecip*     ', ATM_flx_wprecip       , 'f' ) 
     1501prtFlux ('ATM_flx_wemp*        ', ATM_flx_wemp          , 'f', True ) 
     1502 
     1503prtFlux ('ERROR evap           ', ATM_flx_wevap   - ATM_flx_evap  , 'e', True ) 
     1504prtFlux ('ERROR precip         ', ATM_flx_wprecip - ATM_flx_precip, 'e', True ) 
     1505prtFlux ('ERROR snow           ', ATM_flx_wsnow   - ATM_flx_snowf , 'e', True ) 
     1506prtFlux ('ERROR emp            ', ATM_flx_wemp    - ATM_flx_emp   , 'e', True ) 
     1507 
    13871508 
    13881509echo ( '\n====================================================================================' ) 
    … …  
    14011522 
    14021523ATM_flx_budget = -ATM_flx_wbilo + ATM_flx_calving + ATM_flx_runlic #+ ATM_flx_fqfonte + RUN_flx_river 
     1524 
     1525 
    14031526echo ('') 
    14041527#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 )) 
    … …  
    14271550echo ( '\n====================================================================================' ) 
    14281551 
    1429 ATM_flx_runofflic_lic = ATM_flux_int ( ATM_runofflic * ATM_flic ) 
    1430  
    1431 LIC_flx_budget1 = -ATM_flx_wemp_lic  - ATM_flx_calving - ATM_flx_fqfonte 
    1432 LIC_flx_budget2 = -ATM_flx_wbilo_lic - ATM_flx_calving - ATM_flx_fqfonte 
    1433 LIC_flx_budget3 = -ATM_flx_wbilo_lic - ATM_flx_runofflic_lic 
     1552LIC_flx_budget1 = Sprec ( [-ATM_flx_wemp_lic  , -LIC_flx_calving , -LIC_flx_fqfonte] ) 
     1553LIC_flx_budget2 = Sprec ( [-ATM_flx_wbilo_lic , -LIC_flx_calving , -LIC_flx_fqfonte] ) 
     1554LIC_flx_budget3 = Sprec ( [-ATM_flx_wbilo_lic , -LIC_flx_runlic] ) 
     1555LIC_flx_budget4 = Sprec ( [-ATM_flx_wemp_lic  , -LIC_flx_runlic] ) 
    14341556 
    14351557echo ( f'--  LIC  -- {Title} ' ) 
    … …  
    14381560echo ( f'Mass qs      begin = {LIC_mas_qs_beg     :12.6e} kg | Mass end = {LIC_mas_qs_end     :12.6e} kg' ) 
    14391561echo ( f'Mass runlic0 begin = {LIC_mas_runlic0_beg:12.6e} kg | Mass end = {LIC_mas_runlic0_end:12.6e} kg' ) 
    1440 prtFlux ( 'dmass (LIC sno)       ', dLIC_mas_sno          , 'f', True, width=65 ) 
    1441 prtFlux ( 'dmass (LIC qs)        ', dLIC_mas_qs           , 'e', True, width=65 ) 
    1442 prtFlux ( 'dmass (LIC wat)       ', dLIC_mas_wat          , 'f', True, width=65 ) 
    1443 prtFlux ( 'dmass (LIC runlic0)   ', dLIC_mas_runlic0      , 'e', True, width=65 ) 
    1444 prtFlux ( 'dmass (LIC total)     ', dLIC_mas_wat          , 'e', True, width=65 ) 
    1445 prtFlux ( 'LIC ATM_flx_wemp_lic  ', ATM_flx_wemp_lic      , 'f', True, width=65 ) 
    1446 prtFlux ( 'LIC ATM_flx_fqfonte   ', ATM_flx_fqfonte       , 'f', True, width=65 ) 
    1447 prtFlux ( 'LIC ATM_flx_calving   ', ATM_flx_calving       , 'f', True, width=65 ) 
    1448 prtFlux ( 'LIC ATM_flx_runofflic ', ATM_flx_runofflic_lic , 'f', True, width=65 ) 
    1449 prtFlux ( 'LIC fqfonte + calving ', ATM_flx_calving+ATM_flx_fqfonte , 'f', True, width=65 ) 
    1450 prtFlux ( 'LIC fluxes 1 (wevap - precip*frac_lic  - fqcalving - fqfonte) ', LIC_flx_budget1              , 'f', True, width=65 ) 
    1451 prtFlux ( 'LIC fluxes 2 (-wbilo_lic - fqcalving - fqfonte)               ', LIC_flx_budget2              , 'f', True, width=65 ) 
    1452 prtFlux ( 'LIC fluxes 3 (-wbilo_lic - runofflic*frac_lic)                ', LIC_flx_budget3              , 'f', True, width=65 ) 
    1453 prtFlux ( 'LIC error 1                                                   ', LIC_flx_budget1-dLIC_mas_wat , 'e', True, width=65 ) 
    1454 prtFlux ( 'LIC error 2                                                   ', LIC_flx_budget2-dLIC_mas_wat , 'e', True, width=65 ) 
    1455 prtFlux ( 'LIC error 3                                                   ', LIC_flx_budget3-dLIC_mas_wat , 'e', True, width=65 ) 
     1562prtFlux ( 'dmass (LIC sno)       ', dLIC_mas_sno          , 'f', True, width=45 ) 
     1563prtFlux ( 'dmass (LIC qs)        ', dLIC_mas_qs           , 'e', True, width=45 ) 
     1564prtFlux ( 'dmass (LIC wat)       ', dLIC_mas_wat          , 'f', True, width=45 ) 
     1565prtFlux ( 'dmass (LIC runlic0)   ', dLIC_mas_runlic0      , 'e', True, width=45 ) 
     1566prtFlux ( 'dmass (LIC total)     ', dLIC_mas_wat          , 'e', True, width=45 ) 
     1567prtFlux ( 'LIC ATM_flx_wemp_lic  ', ATM_flx_wemp_lic      , 'f', True, width=45 ) 
     1568prtFlux ( 'LIC LIC_flx_fqfonte   ', LIC_flx_fqfonte       , 'f', True, width=45 ) 
     1569prtFlux ( 'LIC LIC_flx_calving   ', LIC_flx_calving       , 'f', True, width=45 ) 
     1570prtFlux ( 'LIC LIC_flx_runofflic ', LIC_flx_runlic        , 'f', True, width=45 ) 
     1571prtFlux ( 'LIC fqfonte + calving ', LIC_flx_calving+LIC_flx_fqfonte , 'f', True, width=45 ) 
     1572prtFlux ( 'LIC fluxes 1 ( wemp_lic   - fqcalving - fqfonte)) ', LIC_flx_budget1              , 'f', True, width=45 ) 
     1573prtFlux ( 'LIC fluxes 2 (-wbilo_lic - fqcalving - fqfonte)   ', LIC_flx_budget2              , 'f', True, width=45 ) 
     1574prtFlux ( 'LIC fluxes 3 (-wbilo_lic - runofflic*frac_lic)    ', LIC_flx_budget3              , 'f', True, width=45 ) 
     1575prtFlux ( 'LIC fluxes 3 ( wemp_lic  - runofflic*frac_lic)    ', LIC_flx_budget4              , 'f', True, width=45 ) 
     1576prtFlux ( 'LIC error 1                                       ', LIC_flx_budget1-dLIC_mas_wat , 'e', True, width=45 ) 
     1577prtFlux ( 'LIC error 2                                       ', LIC_flx_budget2-dLIC_mas_wat , 'e', True, width=45 ) 
     1578prtFlux ( 'LIC error 3                                       ', LIC_flx_budget3-dLIC_mas_wat , 'e', True, width=45 ) 
    14561579echo ( 'LIC error (wevap - precip*frac_lic  - fqcalving - fqfonte)    = {:12.4e} (rel) '.format ( (LIC_flx_budget1-dLIC_mas_wat)/dLIC_mas_wat) ) 
    14571580echo ( 'LIC error (-wbilo_lic - fqcalving - fqfonte)                  = {:12.4e} (rel) '.format ( (LIC_flx_budget2-dLIC_mas_wat)/dLIC_mas_wat) ) 
    … …  
    14681591SRF_flx_emp      = SRF_flux_int ( SRF_emp      ) 
    14691592 
    1470 RUN_flx_torouting  =  RUN_flx_runoff  + RUN_flx_drainage 
    1471 RUN_flx_fromrouting = RUN_flx_coastal + RUN_flx_river 
    1472  
    1473 SRF_flx_all = SRF_flx_rain + SRF_flx_snowf - SRF_flx_evap - RUN_flx_runoff - RUN_flx_drainage 
     1593RUN_flx_torouting   = SRF_flux_int  ( RUN_runoff       + RUN_drainage) 
     1594RUN_flx_fromrouting = ONE_flux_int  ( RUN_coastalflow + RUN_riverflow ) 
     1595 
     1596SRF_flx_all =  SRF_flux_int  ( SRF_rain + SRF_snowf - SRF_evap - RUN_runoff - RUN_drainage ) 
    14741597 
    14751598prtFlux ('rain         ', SRF_flx_rain       , 'f' ) 
Note: See TracChangeset for help on using the changeset viewer.