source: TOOLS/WATER_BUDGET/ATM_waterbudget.py @ 6620

Last change on this file since 6620 was 6620, checked in by omamce, 10 months ago

WATER_BUDET (by OM) : OK for SECHIBA on LMD grid. Still need improvment for ICO.
  • Property svn:executable set to *
  • Property svn:keywords set to Date Revision HeadURL Author
File size: 76.9 KB
Line 
1#!/usr/bin/env python3
2###
3### Script to check water conservation in the IPSL coupled model
4###
5##  Warning, to install, configure, run, use any of included software or
6##  to read the associated documentation you'll need at least one (1)
7##  brain in a reasonably working order. Lack of this implement will
8##  void any warranties (either express or implied).  Authors assumes
9##  no responsability for errors, omissions, data loss, or any other
10##  consequences caused directly or indirectly by the usage of his
11##  software by incorrectly or partially configured personal
12##
13##
14# SVN Information
15SVN = {
16    'Author'  : "$Author$",
17    'Date'    : "$Date$",
18    'Revision': "$Revision$",
19    'Id'      : "$Id: ATM_waterbudget.py 6508 2023-06-13 10:58:38Z omamce $",
20    'HeadURL' : "$HeadUrl: svn+ssh://omamce@forge.ipsl.jussieu.fr/ipsl/forge/projets/igcmg/svn/TOOLS/WATER_BUDGET/ATM_waterbudget.py $"
21    }
22
23###
24## Import system modules
25import sys, os, shutil#, subprocess, platform
26import configparser, re
27
28## Import needed scientific modules
29import numpy as np, xarray as xr
30
31# Check python version
32if sys.version_info < (3, 8, 0) :
33    print ( f'Python version : {platform.python_version()}' ) 
34    raise Exception ( "Minimum Python version is 3.8" )
35
36## Import local module
37import WaterUtils as wu
38import libIGCM_sys
39import nemo, lmdz
40
41## Creates parser for reading .ini input file
42## -------------------------------------------
43config = configparser.ConfigParser ( interpolation=configparser.ExtendedInterpolation() )
44config.optionxform = str # To keep capitals
45
46## Experiment parameters
47## ---------------------
48ATM=None ; ATM_HIS='latlon' ; SRF_HIS='latlon' ; RUN_HIS='latlon' ; ORCA=None ; NEMO=None ; OCE_relax=False
49OCE_icb=False ; Coupled=False ; Routing=None
50TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=None
51
52##
53ARCHIVE=None ; STORAGE=None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None ; TmpDir=None
54FileDir=None ; FileOut=None
55dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None
56FileCommon=None ; file_ATM_his=None ; file_SRF_his=None ; file_RUN_his=None
57file_OCE_his=None ;  file_ICE_his=None ; file_OCE_sca=None
58tar_restart_beg=None ; tar_restart_end=None ; file_ATM_beg=None ; file_ATM_end=None ; file_DYN_beg=None
59file_DYN_end=None ; file_SRF_beg=None ; file_SRF_end=None
60file_RUN_beg=None ; file_RUN_end=None ; file_RUN_end=None ; file_OCE_beg=None
61file_ICE_beg=None ; file_OCE_beg=None
62file_OCE_end=None ; file_ICE_beg=None ; file_OCE_end=None ; file_ICE_end=None
63tar_restart_beg_ATM=None ; tar_restart_beg_DYN=None ; tar_restart_beg_SRF=None
64tar_restart_beg_RUN=None ; tar_restart_beg_OCE=None ; tar_restart_beg_ICE=None
65tar_restart_end_ATM=None ; tar_restart_end_DYN=None ; tar_restart_end_SRF=None
66tar_restart_end_RUN=None ; tar_restart_end_OCE=None ; tar_restart_end_ICE=None
67ContinueOnError=False ; ErrorCount=0 ; SortIco = False
68
69## Read command line arguments
70## ---------------------------
71print ( "Name of Python script:", sys.argv[0] )
72IniFile = sys.argv[1]
73
74# Test existence of IniFile
75if not os.path.exists (IniFile ) :
76    raise FileExistsError ( f"File not found : {IniFile = }" )
77
78if 'full' in IniFile : FullIniFile = IniFile
79else                 : FullIniFile = 'full_' + IniFile
80   
81print ("Input file : ", IniFile )
82config.read (IniFile)
83FullIniFile = 'full_' + IniFile
84
85## Reading config.card if possible
86## -------------------------------
87ConfigCard = None
88
89if 'Experiment' in config.keys ()  : ## Read Experiment on Config file if possible
90    if 'ConfigCard' in config['Experiment'].keys () :
91        ConfigCard = config['Experiment']['ConfigCard']
92        print ( f'{ConfigCard=}' )
93
94if ConfigCard : ## Read config card if it exists
95    # Text existence of ConfigCard
96    if os.path.exists ( ConfigCard ) :
97        print ( f'Reading Config Card : {ConfigCard}' )
98        ## Creates parser for reading .ini input file
99        MyReader = configparser.ConfigParser (interpolation=configparser.ExtendedInterpolation() )
100        MyReader.optionxform = str # To keep capitals
101       
102        MyReader.read (ConfigCard)
103
104        for VarName in ['JobName', 'ExperimentName', 'SpaceName', 'LongName', 'ModelName', 'TagName'] :
105            if VarName in MyReader['UserChoices'].keys() :
106                locals()[VarName] = MyReader['UserChoices'][VarName]
107                exec  ( f'{VarName} = wu.setBool ({VarName})' )
108                exec  ( f'{VarName} = wu.setNum  ({VarName})' )
109                exec  ( f'{VarName} = wu.setNone ({VarName})' )
110                exec  ( f'wu.{VarName} = {VarName}' )
111                print ( f'    {VarName:21} set to : {locals()[VarName]:}' )
112               
113        for VarName in ['PackFrequency'] :
114            if VarName in MyReader['Post'].keys() :
115                locals()[VarName] = MyReader['Post'][VarName]
116                exec  ( f'{VarName} = wu.setBool ({VarName})' )
117                exec  ( f'{VarName} = wu.setNum  ({VarName})' )
118                exec  ( f'{VarName} = wu.setNone ({VarName})' )
119                exec  ( f'wu.{VarName} = {VarName}' )
120                print ( f'    {VarName:21} set to : {locals()[VarName]:}' )
121    else :
122        raise FileExistsError ( f"File not found : {ConfigCard = }" )
123               
124## Reading config file
125## -------------------
126for Section in ['Config', 'Experiment', 'libIGCM', 'Files', 'Physics' ] :
127   if Section in config.keys () : 
128        print ( f'\nReading [{Section}]' )
129        for VarName in config[Section].keys() :
130            locals()[VarName] = config[Section][VarName]
131            exec  ( f'{VarName} = wu.setBool ({VarName})' )
132            exec  ( f'{VarName} = wu.setNum  ({VarName})' )
133            exec  ( f'{VarName} = wu.setNone ({VarName})' )
134            exec  ( f'wu.{VarName} = {VarName}' )
135            print ( f'    {VarName:21} set to : {locals()[VarName]}' )
136            #exec ( f'del {VarName}' )
137
138print ( f'\nConfig file readed : {IniFile} ' )
139
140## Some physical constants
141## =======================
142if wu.unDefined ( 'Ra' )           : Ra          = wu.Ra           #-- Earth Radius (m)
143if wu.unDefined ( 'Grav' )         : Grav        = wu.Grav         #-- Gravity (m^2/s
144if wu.unDefined ( 'ICE_rho_ice' )  : ICE_rho_ice = wu.ICE_rho_ice  #-- Ice volumic mass (kg/m3) in LIM3
145if wu.unDefined ( 'ICE_rho_sno')   : ICE_rho_sno = wu.ICE_rho_sno  #-- Snow volumic mass (kg/m3) in LIM3
146if wu.unDefined ( 'OCE_rho_liq' )  : OCE_rho_liq = wu.OCE_rho_liq  #-- Ocean water volumic mass (kg/m3) in NEMO
147if wu.unDefined ( 'ATM_rho' )      : ATM_rho     = wu.ATM_rho      #-- Water volumic mass in atmosphere (kg/m^3)
148if wu.unDefined ( 'SRF_rho' )      : SRF_rho     = wu.SRF_rho      #-- Water volumic mass in surface reservoir (kg/m^3)
149if wu.unDefined ( 'RUN_rho' )      : RUN_rho     = wu.RUN_rho      #-- Water volumic mass of rivers (kg/m^3)
150if wu.unDefined ( 'ICE_rho_pnd' )  : ICE_rho_pnd = wu.ICE_rho_pnd  #-- Water volumic mass in ice ponds (kg/m^3)
151if wu.unDefined ( 'YearLength' )   : YearLength  = wu.YearLength   #-- Year length (s)
152   
153## Set libIGCM and machine dependant values
154## ----------------------------------------
155if not 'Files' in config.keys () : config['Files'] = {}
156
157config['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
160config['Config'] = { 'ContinueOnError':ContinueOnError, 'SortIco':SortIco }
161
162## --------------------------
163ICO  = ( 'ICO' in wu.ATM )
164LMDZ = ( 'LMD' in wu.ATM )
165
166mm = libIGCM_sys.config ( TagName=TagName, SpaceName=SpaceName, ExperimentName=ExperimentName, JobName=JobName, User=User, Group=Group,
167             ARCHIVE=None, SCRATCHDIR=None, STORAGE=None, R_IN=None, R_OUT=None, R_FIG=None, rebuild=None, TmpDir=None, 
168             R_SAVE=None, R_FIGR=None, R_BUFR=None, R_BUF_KSH=None, REBUILD_DIR=None, POST_DIR=None )
169globals().update(mm)
170
171config['Files']['TmpDir'] = TmpDir
172config['libIGCM'] = { 'ARCHIVE':ARCHIVE, 'STORAGE':STORAGE, 'TmpDir':TmpDir, 'R_IN':R_IN, 'rebuild':rebuild }
173
174## Defines begining and end of experiment
175## --------------------------------------
176if wu.unDefined ( 'DateBegin' ) :
177    DateBegin = f'{YearBegin}0101'
178    config['Experiment']['DateBegin'] = DateBegin
179else :
180    YearBegin, MonthBegin, DayBegin = wu.DateSplit ( DateBegin )
181    DateBegin = wu.FormatToGregorian ( DateBegin)
182    config['Experiment']['YearBegin'] = YearBegin
183
184if wu.unDefined ( 'DateEnd' )   :
185    DateEnd   = f'{YearEnd}1231'
186    config['Experiment']['DateEnd']   = DateEnd
187else :
188    YearEnd, MonthEnd, DayEnd = wu.DateSplit ( DateEnd )
189    DateEnd   = wu.FormatToGregorian ( DateEnd)
190
191if wu.unDefined ( 'PackFrequency' ) :
192    PackFrequency = YearEnd - YearBegin + 1
193    config['Experiment']['PackFrequency']   = f'{PackFrequency}'
194
195if type ( PackFrequency ) == str : 
196    if 'Y' in PackFrequency : PackFrequency = PackFrequency.replace ( 'Y', '') 
197    if 'M' in PackFrequency : PackFrequency = PackFrequency.replace ( 'M', '')
198    PackFrequency = int ( PackFrequency )
199   
200## Output file with water budget diagnostics
201## -----------------------------------------
202if wu.unDefined ( 'FileOut' ) : FileOut = f'ATM_waterbudget_{JobName}_{DateBegin}_{DateEnd}.out'
203config['Files']['FileOut'] = FileOut
204
205f_out = open ( FileOut, mode = 'w' )
206
207## Useful functions
208## ----------------
209def kg2Sv    (val, rho=ATM_rho) :
210    '''From kg to Sverdrup'''
211    return val/dtime_sec*1.0e-6/rho
212
213def kg2myear (val, rho=ATM_rho) :
214    '''From kg to m/year'''
215    return val/ATM_aire_sea_tot/rho/NbYear
216
217def var2prt (var, small=False, rho=ATM_rho) :
218    if small :  return var , kg2Sv(var, rho=rho)*1000., kg2myear(var, rho=rho)*100#
219    else     :  return var , kg2Sv(var, rho=rho)      , kg2myear(var, rho=rho)
220
221def prtFlux (Desc, var, Form='F', small=False, rho=ATM_rho, width=15) :
222    if small :
223        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} mSv | {:12.4f} mm/year "
224        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} mSv | {:12.4e} mm/year "
225    else :
226        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} Sv  | {:12.4f} m/year  "
227        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} Sv  | {:12.4e} m/year  "
228    echo ( (' {:>{width}} = ' +ff).format (Desc, *var2prt(var, small=small, rho=rho), width=width ) )
229    return None
230
231def echo (string, end='\n') :
232    '''Function to print to stdout *and* output file'''
233    print ( str(string), end=end  )
234    sys.stdout.flush ()
235    f_out.write ( str(string) + end )
236    f_out.flush ()
237    return None
238   
239## Set libIGCM directories
240## -----------------------
241if wu.unDefined ('R_OUT'      ) : R_OUT       = os.path.join ( ARCHIVE   , 'IGCM_OUT' )
242if wu.unDefined ('R_BUF'      ) : R_BUF       = os.path.join ( SCRATCHDIR, 'IGCM_OUT' )
243if wu.unDefined ('L_EXP'      ) : L_EXP       = os.path.join (TagName, SpaceName, ExperimentName, JobName)
244if wu.unDefined ('R_SAVE'     ) : R_SAVE      = os.path.join ( R_OUT, L_EXP )
245if wu.unDefined ('R_BUFR'     ) : R_BUFR      = os.path.join ( R_BUF, L_EXP )
246if wu.unDefined ('POST_DIR'   ) : POST_DIR    = os.path.join ( R_BUFR, 'Out' )
247if wu.unDefined ('REBUILD_DIR') : REBUILD_DIR = os.path.join ( R_BUFR, 'REBUILD' )
248if wu.unDefined ('R_BUF_KSH'  ) : R_BUF_KSH   = os.path.join ( R_BUFR, 'Out' )
249if wu.unDefined ('R_FIGR'     ) : R_FIGR      = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP )
250
251config['libIGCM'].update ( { 'R_OUT':R_OUT, 'R_BUF':R_BUF, 'L_EXP':L_EXP, 'R_BUFR':R_BUFR, 'POST_DIR':POST_DIR,
252                      'REBUILD_DIR':REBUILD_DIR, 'R_BUF_KSH':R_BUF_KSH, 'R_FIGR':R_FIGR, 'rebuild':rebuild } )
253
254## Set directory to extract files
255## ------------------------------
256if wu.unDefined ( 'FileDir' ) : FileDir = os.path.join ( TmpDir, f'WATER_{JobName}' )
257config['Files']['FileDir'] = FileDir
258
259if not os.path.isdir ( FileDir ) : os.makedirs ( FileDir )
260
261##- Set directories to rebuild ocean and ice restart files
262if wu.unDefined ( 'FileDirOCE' ) : FileDirOCE = os.path.join ( FileDir, 'OCE' )
263if wu.unDefined ( 'FileDirICE' ) : FileDirICE = os.path.join ( FileDir, 'ICE' )
264if not os.path.exists ( FileDirOCE ) : os.mkdir ( FileDirOCE )
265if not os.path.exists ( FileDirICE ) : os.mkdir ( FileDirICE )
266
267echo (' ')
268echo ( f'JobName    : {JobName}'   )
269echo ( f'Comment    : {Comment}'   )
270echo ( f'TmpDir     : {TmpDir}'    )
271echo ( f'FileDir     : {FileDir}'    )
272echo ( f'FileDirOCE  : {FileDirOCE}' )
273echo ( f'FileDirICE  : {FileDirICE}' )
274
275echo ( f'\nDealing with {L_EXP}'  )
276
277## Creates model output directory names
278## ------------------------------------
279if Freq == "MO" : FreqDir = os.path.join ('Output' , 'MO' )
280if Freq == "SE" : FreqDir = os.path.join ('Analyse', 'SE' )
281if wu.unDefined ('dir_ATM_his' ) :
282    dir_ATM_his = os.path.join ( R_SAVE, "ATM", FreqDir )
283    config['Files']['dir_ATM_his'] = dir_ATM_his
284if wu.unDefined ( 'dir_SRF_his' ) : 
285    dir_SRF_his = os.path.join ( R_SAVE, "SRF", FreqDir )
286    config['Files']['dir_SRF_his'] = dir_SRF_his
287   
288echo ( f'The analysis relies on files from the following model output directories : ' )
289echo ( f'{dir_ATM_his = }' )
290echo ( f'{dir_SRF_his = }' )
291
292##-- Creates files names
293if wu.unDefined ( 'Period' ) :
294    if Freq == 'MO' : Period = f'{DateBegin}_{DateEnd}_1M'
295    if Freq == 'SE' : Period = f'SE_{DateBegin}_{DateEnd}_1M'
296    config['Files']['Period'] = Period
297
298config['Files']['DateBegin'] = DateBegin
299config['Files']['DateBegin'] = DateEnd
300       
301echo ( f'Period   : {Period}' )
302   
303if wu.unDefined ( 'FileCommon' ) :
304    FileCommon = f'{JobName}_{Period}'
305    config['Files']['FileCommon'] = FileCommon
306
307if wu.unDefined ( 'Title' ) :
308    Title = f'{JobName} : {Freq} : {DateBegin} - {DateEnd}'
309    config['Files']['Title'] = Title
310     
311echo ('\nOpen history files' )
312if wu.unDefined ( 'file_ATM_his' ) :
313    if ATM_HIS == 'latlon' :
314        file_ATM_his = os.path.join ( dir_ATM_his, f'{FileCommon}_histmth.nc' )
315    if ATM_HIS == 'ico' :
316        file_ATM_his = os.path.join ( dir_ATM_his, f'{FileCommon}_histmth_ico.nc' )
317    config['Files']['file_ATM_his'] = file_ATM_his
318if wu.unDefined ( 'file_SRF_his' ) :
319    if ATM_HIS == 'latlon' :
320        file_SRF_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' )
321    if ATM_HIS == 'ico' :
322        file_SRF_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history_ico.nc' )
323    config['Files']['file_SRF_his'] = file_SRF_his
324
325if Routing == 'SIMPLE' :
326    if file_RUN_his == None :
327        if ATM_HIS == 'latlon' :
328            file_RUN_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' )
329        if ATM_HIS == 'ico' :
330            file_RUN_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history_ico.nc' )
331        config['Files']['file_RUN_his'] = file_RUN_his
332
333d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
334d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
335if Routing == 'SECHIBA' : d_RUN_his = d_SRF_his
336if Routing == 'SIMPLE'  : d_RUN_his = xr.open_dataset ( file_RUN_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
337
338echo ( f'{file_ATM_his = }' )
339echo ( f'{file_SRF_his = }' )
340if Routing == 'SIMPLE' : echo ( f'{file_RUN_his = }' )
341
342## Compute run length
343## ------------------
344dtime = ( d_ATM_his.time_counter_bounds.max() - d_ATM_his.time_counter_bounds.min() )
345echo ('\nRun length : {:8.2f} days'.format ( (dtime/np.timedelta64(1, "D")).values ) )
346dtime_sec = (dtime/np.timedelta64(1, "s")).values.item() # Convert in seconds
347
348##-- Compute length of each period
349dtime_per = (d_ATM_his.time_counter_bounds[:,-1] -  d_ATM_his.time_counter_bounds[:,0] )
350echo ('\nPeriods lengths (days) : {:} days'.format ( (dtime_per/np.timedelta64(1, "D")).values ) )
351dtime_per_sec = (dtime_per/np.timedelta64(1, "s")).values # In seconds
352dtime_per_sec = xr.DataArray (dtime_per_sec, dims=["time_counter", ], coords=[d_ATM_his.time_counter,] )
353dtime_per_sec.attrs['unit'] = 's'
354
355##-- Number of years (approximative)
356NbYear = dtime_sec / YearLength
357
358##-- Extract restart files from tar
359
360if wu.unDefined ('TarRestartDate_beg' ) : TarRestartDate_beg = wu.DateMinusOneDay ( DateBegin )
361if wu.unDefined ('TarRestartDate_end' ) : TarRestartDate_end = wu.FormatToGregorian ( DateEnd )
362
363if wu.unDefined ( 'TarRestartPeriod_beg' ) :
364
365    TarRestartPeriod_beg_DateEnd = TarRestartDate_beg
366    TarRestartPeriod_beg_DateBeg = wu.DateAddYear ( TarRestartPeriod_beg_DateEnd, -PackFrequency )
367    TarRestartPeriod_beg_DateBeg = wu.DatePlusOneDay ( TarRestartPeriod_beg_DateBeg )
368   
369    TarRestartPeriod_beg = f'{TarRestartPeriod_beg_DateBeg}_{TarRestartPeriod_beg_DateEnd}'
370    echo (f'Tar period for initial restart : {TarRestartPeriod_beg}')
371    config['Files']['TarRestartPeriod_beg'] = TarRestartPeriod_beg
372
373if wu.unDefined ( 'TarRestartPeriod_end' ) :
374
375    TarRestartPeriod_end_DateEnd = TarRestartDate_end
376    TarRestartPeriod_end_DateBeg = wu.DateAddYear ( TarRestartPeriod_end_DateEnd, -PackFrequency )
377    TarRestartPeriod_end_DateBeg = wu.DatePlusOneDay ( TarRestartPeriod_end_DateBeg )
378     
379    TarRestartPeriod_end = f'{TarRestartPeriod_end_DateBeg}_{TarRestartPeriod_end_DateEnd}'
380    echo (f'Tar period for final restart : {TarRestartPeriod_end}')
381    config['Files']['TarRestartPeriod_end'] = TarRestartPeriod_end
382
383echo (f'Restart dates - Start : {TarRestartPeriod_beg}  /  End : {TarRestartPeriod_end}')
384
385if wu.unDefined ( 'tar_restart_beg' ) :
386    tar_restart_beg = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_beg}_restart.tar' )
387    config['Files']['tar_restart_beg'] = tar_restart_beg
388if wu.unDefined ( 'tar_restart_end' ) :
389    tar_restart_end = os.path.join ( R_SAVE, 'RESTART', f'{JobName}_{TarRestartPeriod_end}_restart.tar' )
390    config['Files']['tar_restart_end'] = tar_restart_end
391   
392echo ( f'{tar_restart_beg = }' )
393echo ( f'{tar_restart_end = }' )
394
395##-- Names of tar files with restarts
396if wu.unDefined ( 'SRF_HIS' ) : SRF_HIS = ATM_HIS
397
398if wu.unDefined ( 'tar_restart_beg_ATM' ) : tar_restart_beg_ATM = tar_restart_beg
399if wu.unDefined ( 'tar_restart_beg_DYN' ) : tar_restart_beg_DYN = tar_restart_beg
400if wu.unDefined ( 'tar_restart_beg_SRF' ) : tar_restart_beg_SRF = tar_restart_beg
401if wu.unDefined ( 'tar_restart_beg_RUN' ) : tar_restart_beg_RUN = tar_restart_beg
402if wu.unDefined ( 'tar_restart_beg_OCE' ) : tar_restart_beg_OCE = tar_restart_beg
403if wu.unDefined ( 'tar_restart_beg_ICE' ) : tar_restart_beg_ICE = tar_restart_beg
404
405if wu.unDefined ( 'tar_restart_end_ATM' ) : tar_restart_end_ATM = tar_restart_end
406if wu.unDefined ( 'tar_restart_end_DYN' ) : tar_restart_end_DYN = tar_restart_end
407if wu.unDefined ( 'tar_restart_end_SRF' ) : tar_restart_end_SRF = tar_restart_end
408if wu.unDefined ( 'tar_restart_end_RUN' ) : tar_restart_end_RUN = tar_restart_end
409if wu.unDefined ( 'tar_restart_end_OCE' ) : tar_restart_end_OCE = tar_restart_end
410if wu.unDefined ( 'tar_restart_end_ICE' ) : tar_restart_end_ICE = tar_restart_end
411
412if wu.unDefined ( 'file_ATM_beg' ) : 
413    file_ATM_beg = f'{FileDir}/ATM_{JobName}_{TarRestartDate_beg}_restartphy.nc'
414    config['Files']['file_ATM_beg'] = file_ATM_beg
415if wu.unDefined ( 'file_ATM_end' ) :
416    file_ATM_end = f'{FileDir}/ATM_{JobName}_{TarRestartDate_end}_restartphy.nc'
417    config['Files']['file_ATM_end'] = file_ATM_end
418
419liste_beg = [file_ATM_beg, ]
420liste_end = [file_ATM_end, ]
421
422if wu.unDefined ( 'file_DYN_beg' ) :
423    if LMDZ : file_DYN_beg = f'{FileDir}/ATM_{JobName}_{TarRestartDate_beg}_restart.nc'
424    if ICO  : file_DYN_beg = f'{FileDir}/ICO_{JobName}_{TarRestartDate_beg}_restart.nc'
425    liste_beg.append (file_DYN_beg)
426    config['Files']['file_DYN_beg'] = file_DYN_beg
427   
428if wu.unDefined ( 'file_DYN_end' ) : 
429    if LMDZ : file_DYN_end = f'{FileDir}/ATM_{JobName}_{TarRestartDate_end}_restart.nc'
430    if ICO  : file_DYN_end = f'{FileDir}/ICO_{JobName}_{TarRestartDate_end}_restart.nc'
431    liste_end.append ( file_DYN_end )
432    config['Files']['file_DYN_end'] = file_DYN_end
433
434if wu.unDefined ( 'file_SRF_beg' ) :
435    file_SRF_beg = f'{FileDir}/SRF_{JobName}_{TarRestartDate_beg}_sechiba_rest.nc'
436    config['Files']['file_SRF_beg'] = file_SRF_beg
437if wu.unDefined ( 'file_SRF_end' ) :
438    file_SRF_end = f'{FileDir}/SRF_{JobName}_{TarRestartDate_end}_sechiba_rest.nc'
439    config['Files']['file_SRF_end'] = file_SRF_end
440   
441liste_beg.append ( file_SRF_beg )
442liste_end.append ( file_SRF_end )
443
444echo ( f'{file_ATM_beg = }')
445echo ( f'{file_ATM_end = }')
446echo ( f'{file_DYN_beg = }')
447echo ( f'{file_DYN_end = }')
448echo ( f'{file_SRF_beg = }')
449echo ( f'{file_SRF_end = }')
450
451if ICO :
452    if wu.unDefined ('file_DYN_aire') : file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ATM+'_grid.nc' )
453    config['Files']['file_DYN_aire'] = file_DYN_aire
454
455if Routing == 'SIMPLE' :
456    if wu.unDefined ( 'file_RUN_beg' ) : 
457        file_RUN_beg = f'{FileDir}/SRF_{JobName}_{TarRestartDate_beg}_routing_restart.nc'
458        config['Files']['file_RUN_beg'] = file_RUN_beg
459    if wu.unDefined ( 'file_RUN_end' ) : 
460        file_RUN_end = f'{FileDir}/SRF_{JobName}_{TarRestartDate_end}_routing_restart.nc'
461        config['Files']['file_RUN_end'] = file_RUN_end
462       
463    liste_beg.append ( file_RUN_beg )
464    liste_end.append ( file_RUN_end )
465    echo ( f'{file_RUN_beg = }' )
466    echo ( f'{file_RUN_end = }' )
467
468echo ('\nExtract restart files from tar : ATM, ICO and SRF')
469
470def extract_and_rebuild ( file_name=file_ATM_beg, tar_restart=tar_restart_end, FileDirComp=FileDir, ErrorCount=ErrorCount ) :
471    echo ( f'----------')
472    echo ( f'file to extract : {file_name = }' )
473    if os.path.exists ( os.path.join (FileDir, file_name) ) :
474        echo ( f'file found : {file_name = }' )
475    else :
476        echo ( f'file not found : {file_name = }' )
477        base_resFile = os.path.basename (file_name)
478        if os.path.exists ( tar_restart ) :
479            command =  f'cd {FileDir} ; tar xf {tar_restart} {base_resFile}'
480            echo ( f'{command = }' )
481            try : os.system ( command )
482            except :
483                if ContinueOnError :
484                    ErrorCount += 1
485                    echo ( f'****** Command failed : {command}' )
486                    echo ( f'****** Trying to continue' )
487                    echo ( ' ')
488                else :
489                    raise Exception ( f'**** command failed : {command} - Stopping' )
490            else :
491                echo ( f'tar done : {base_resFile}' )
492        else :
493            echo ( f'****** Tar restart file {tar_restart = } not found ' )
494            if ContinueOnError :
495                ErrorCount += 1
496            else : 
497                raise Exception ( f'****** tar file not found {tar_restart = } - Stopping' )
498    return ErrorCount
499 
500ErrorCount += extract_and_rebuild ( file_name=file_ATM_beg, tar_restart=tar_restart_beg_ATM, FileDirComp=FileDir )
501ErrorCount += extract_and_rebuild ( file_name=file_DYN_beg, tar_restart=tar_restart_beg_DYN, FileDirComp=FileDir )
502ErrorCount += extract_and_rebuild ( file_name=file_SRF_beg, tar_restart=tar_restart_beg_SRF, FileDirComp=FileDir )
503
504ErrorCount += extract_and_rebuild ( file_name=file_ATM_end, tar_restart=tar_restart_end_ATM, FileDirComp=FileDir )
505ErrorCount += extract_and_rebuild ( file_name=file_DYN_end, tar_restart=tar_restart_end_DYN, FileDirComp=FileDir )
506ErrorCount += extract_and_rebuild ( file_name=file_SRF_end, tar_restart=tar_restart_end_SRF, FileDirComp=FileDir )
507
508if Routing == 'SIMPLE' :
509    ErrorCount += extract_and_rebuild ( file_name=file_RUN_beg, tar_restart=tar_restart_beg_RUN, FileDirComp=FileDir )
510    ErrorCount += extract_and_rebuild ( file_name=file_RUN_end, tar_restart=tar_restart_end_RUN, FileDirComp=FileDir )
511
512##-- Exit in case of error in the opening file phase
513if ErrorCount > 0 :
514    echo ( '  ' )
515    raise Exception ( f'**** Some files missing - Stopping - {ErrorCount = }' )
516
517##
518echo ('\nOpening ATM SRF and ICO restart files')
519d_ATM_beg = xr.open_dataset ( os.path.join (FileDir, file_ATM_beg), decode_times=False, decode_cf=True ).squeeze()
520d_ATM_end = xr.open_dataset ( os.path.join (FileDir, file_ATM_end), decode_times=False, decode_cf=True ).squeeze()
521d_SRF_beg = xr.open_dataset ( os.path.join (FileDir, file_SRF_beg), decode_times=False, decode_cf=True ).squeeze()
522d_SRF_end = xr.open_dataset ( os.path.join (FileDir, file_SRF_end), decode_times=False, decode_cf=True ).squeeze()
523d_DYN_beg = xr.open_dataset ( os.path.join (FileDir, file_DYN_beg), decode_times=False, decode_cf=True ).squeeze()
524d_DYN_end = xr.open_dataset ( os.path.join (FileDir, file_DYN_end), decode_times=False, decode_cf=True ).squeeze()
525
526for var in d_SRF_beg.variables :
527    d_SRF_beg[var] = d_SRF_beg[var].where ( d_SRF_beg[var]<1.e20, 0.)
528    d_SRF_end[var] = d_SRF_end[var].where ( d_SRF_end[var]<1.e20, 0.)
529
530if Routing == 'SIMPLE' :
531    d_RUN_beg = xr.open_dataset ( os.path.join (FileDir, file_RUN_beg), decode_times=False, decode_cf=True ).squeeze()
532    d_RUN_end = xr.open_dataset ( os.path.join (FileDir, file_RUN_end), decode_times=False, decode_cf=True ).squeeze()
533
534def to_cell ( dd, newname='cell' ) :
535    '''Set space dimension to 'cell' '''
536    for oldname in [ 'cell_mesh', 'y', 'points_physiques' ] :
537        try    : dd = dd.rename ( {oldname : newname} )
538        except : pass
539    return dd
540
541d_ATM_beg = to_cell ( d_ATM_beg )
542d_ATM_end = to_cell ( d_ATM_end )
543d_SRF_beg = to_cell ( d_SRF_beg )
544d_SRF_end = to_cell ( d_SRF_end )
545d_DYN_beg = to_cell ( d_DYN_beg )
546d_DYN_end = to_cell ( d_DYN_end )
547
548if Routing == 'SIMPLE' :
549    d_RUN_beg = to_cell ( d_RUN_beg )
550    d_RUN_end = to_cell ( d_RUN_end )
551
552d_ATM_his = to_cell ( d_ATM_his )
553d_SRF_his = to_cell ( d_SRF_his )
554   
555echo ( f'{file_ATM_beg = }' )
556echo ( f'{file_ATM_end = }' )
557echo ( f'{file_DYN_beg = }' )
558echo ( f'{file_DYN_end = }' )
559echo ( f'{file_SRF_beg = }' )
560echo ( f'{file_SRF_end = }' )
561if Routing == 'SIMPLE' :
562    echo ( f'{file_RUN_beg = }' )
563    echo ( f'{file_RUN_end = }' )
564
565## Write the full configuration
566config_out = open (FullIniFile, 'w')
567config.write ( config_out )
568config_out.close ()
569
570if ICO :
571    if SortIco : 
572        # Creation d'une clef de tri pour les restarts (pour chaque restart !)
573        DYN_beg_keysort = np.lexsort ( (d_DYN_beg['lat_mesh'], d_DYN_beg['lon_mesh'] ) )
574        ATM_beg_keysort = np.lexsort ( (d_ATM_beg['latitude'], d_ATM_beg['longitude']) )
575        SRF_beg_keysort = np.lexsort ( (d_SRF_beg['nav_lat' ], d_SRF_beg['nav_lon']  ) )
576        DYN_end_keysort = np.lexsort ( (d_DYN_end['lat_mesh'], d_DYN_end['lon_mesh'] ) )
577        ATM_end_keysort = np.lexsort ( (d_ATM_end['latitude'], d_ATM_end['longitude']) )
578        SRF_end_keysort = np.lexsort ( (d_SRF_end['nav_lat' ], d_SRF_end['nav_lon']  ) )
579       
580        if ATM_HIS == 'ico' :
581            ATM_his_keysort = np.lexsort ( (d_ATM_his['lat'], d_ATM_his['lon'] ) )
582        if SRF_HIS == 'ico' :
583            SRF_his_keysort = np.lexsort ( (d_SRF_his['lat'], d_SRF_his['lon'] ) )
584            RUN_his_keysort = SRF_his_keysort
585    else : 
586        DYN_beg_keysort = np.arange ( len ( d_DYN_beg['lat_mesh'] ) )
587        ATM_beg_keysort = DYN_beg_keysort
588        SRF_beg_keysort = DYN_beg_keysort
589
590        DYN_end_keysort = DYN_beg_keysort
591        ATM_end_keysort = DYN_beg_keysort
592        SRF_end_keysort = DYN_beg_keysort
593   
594        ATM_his_keysort = DYN_beg_keysort
595        SRF_his_keysort = DYN_beg_keysort
596        RUN_his_keysort = SRF_his_keysort
597       
598# ATM grid with cell surfaces
599if 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   
615if ICO :
616    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' )
634
635    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]
672        SRF_aire      =  SRF_areas * SRF_contfrac
673
674ATM_fsea      = ATM_foce + ATM_fsic
675ATM_flnd      = ATM_fter + ATM_flic
676ATM_aire_fter = ATM_aire * ATM_fter
677ATM_aire_flic = ATM_aire * ATM_flic
678ATM_aire_fsic = ATM_aire * ATM_fsic
679ATM_aire_foce = ATM_aire * ATM_foce
680ATM_aire_flnd = ATM_aire * ATM_flnd
681ATM_aire_fsea = ATM_aire * ATM_fsea
682
683SRF_aire = SRF_aire.where ( np.abs (SRF_aire) < 1E15, 0. )
684
685## Write the full configuration
686config_out = open (FullIniFile, 'w')
687config.write ( config_out )
688config_out.close ()
689           
690if ICO :
691    # Area on icosahedron grid
692    d_DYN_aire = xr.open_dataset ( file_DYN_aire, decode_times=False ).squeeze()
693
694    if SortIco : 
695        # Creation d'une clef de tri pour le fichier aire
696        DYN_aire_keysort = np.lexsort ( (d_DYN_aire['lat'], d_DYN_aire['lon']) )
697    else : 
698        DYN_aire_keysort = np.arange ( len ( d_DYN_aire['lat'] ) )
699           
700    DYN_lat = d_DYN_aire['lat'][DYN_aire_keysort]
701    DYN_lon = d_DYN_aire['lon'][DYN_aire_keysort]
702
703    DYN_aire = d_DYN_aire['aire'][DYN_aire_keysort]
704    DYN_fsea = d_DYN_aire['fract_oce'][DYN_aire_keysort] + d_DYN_aire['fract_sic'][DYN_aire_keysort]
705    DYN_flnd = 1.0 - DYN_fsea
706    DYN_fter = d_ATM_beg['FTER'][ATM_beg_keysort]
707    DYN_flic = d_ATM_beg['FLIC'][ATM_beg_keysort]
708    DYN_aire_fter = DYN_aire * DYN_fter
709   
710if LMDZ :
711    # Area on lon/lat grid
712    DYN_aire = ATM_aire
713    DYN_fsea = ATM_fsea
714    DYN_flnd = ATM_flnd
715    DYN_fter = d_ATM_beg['FTER']
716    DYN_flic = d_ATM_beg['FLIC']
717    DYN_aire_fter = DYN_aire * DYN_fter
718
719# Functions computing integrals and sum
720def DYN_stock_int (stock) :
721    '''Integrate (* surface) stock on atmosphere grid'''
722    DYN_stock_int  = wu.Psum ( (stock * DYN_aire).to_masked_array().ravel() ) 
723    return DYN_stock_int
724
725def ATM_flux_int (flux) :
726    '''Integrate (* time * surface) flux on atmosphere grid'''
727    ATM_flux_int  = wu.Psum ( (flux * dtime_per_sec * ATM_aire).to_masked_array().ravel() )
728    return ATM_flux_int
729
730def SRF_stock_int (stock) :
731    '''Integrate (* surface) stock on land grid'''
732    #SRF_stock_int  = wu.Ksum (  ( (stock * DYN_aire_fter).to_masked_array().ravel()) ) # Marche avec ICO et lon/lat na
733    SRF_stock_int  = wu.Ksum (  ( (stock * DYN_aire_fter).to_masked_array().ravel()) )
734    return SRF_stock_int
735
736def SRF_flux_int (flux) :
737    '''Integrate (* time * surface) flux on land grid'''
738    SRF_flux_int  = wu.Psum (  (flux * dtime_per_sec * SRF_aire).to_masked_array().ravel() )
739    return SRF_flux_int
740
741def ONE_stock_int (stock) :
742    '''Sum stock '''
743    ONE_stock_int  = wu.Psum ( stock.to_masked_array().ravel() )
744    return ONE_stock_int
745
746def ONE_flux_int (flux) :
747    '''Integrate (* time) flux on area=1 grid '''
748    ONE_flux_int  = wu.Psum ( (flux * dtime_per_sec ).to_masked_array().ravel() )
749    return ONE_flux_int
750
751ATM_aire_sea     = ATM_aire * ATM_fsea
752ATM_aire_tot     = ONE_stock_int (ATM_aire)
753ATM_aire_sea_tot = ONE_stock_int (ATM_aire*ATM_fsea)
754
755DYN_aire_tot     =  ONE_stock_int ( DYN_aire )
756SRF_aire_tot     =  ONE_stock_int ( SRF_aire )
757
758echo ('')
759echo ( f'ATM DYN : Area of atmosphere : {DYN_aire_tot:18.9e}' )
760echo ( f'ATM HIS : Area of atmosphere : {ATM_aire_tot:18.9e}' )
761echo ( f'ATM SRF : Area of atmosphere : {SRF_aire_tot:18.9e}' )
762echo ('')
763echo ( 'ATM DYN : Area of atmosphere/(4pi R^2) : {:18.9f}'.format(DYN_aire_tot/(Ra*Ra*4*np.pi) ) )
764echo ( 'ATM HIS : Area of atmosphere/(4pi R^2) : {:18.9f}'.format(ATM_aire_tot/(Ra*Ra*4*np.pi) ) )
765echo ( 'ATM SRF : Area of atmosphere/(4pi R^2) : {:18.9f}'.format(SRF_aire_tot/(Ra*Ra*4*np.pi) ) )
766echo ('')
767echo ( f'ATM SRF : Area of atmosphere (no contfrac): {ONE_stock_int (SRF_areas):18.9e}' )
768
769
770if (  np.abs (ATM_aire_tot/(Ra*Ra*4*np.pi) - 1.0) > 0.01 ) :
771    raise Exception ('Error of atmosphere surface interpolated on lon/lat grid')
772
773echo ( '\n====================================================================================' )
774echo ( f'-- ATM changes in stores  -- {Title} ' )
775
776#-- Change in precipitable water from the atmosphere daily and monthly files
777#-- Compute sum weighted by gridcell area (kg/m2) then convert to Sv
778
779echo ( 'ATM vertical grid' )
780ATM_Ahyb = d_ATM_his['Ahyb'].squeeze()
781ATM_Bhyb = d_ATM_his['Bhyb'].squeeze()
782
783echo ( 'Surface pressure' )
784if ICO :
785    DYN_psol_beg = d_DYN_beg['ps'][DYN_beg_keysort]
786    DYN_psol_end = d_DYN_end['ps'][DYN_end_keysort]
787if 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' )
790   
791echo ( '3D Pressure at the interface layers (not scalar points)' )
792DYN_pres_beg = ATM_Ahyb + ATM_Bhyb * DYN_psol_beg
793DYN_pres_end = ATM_Ahyb + ATM_Bhyb * DYN_psol_end
794
795klevp1 = ATM_Bhyb.shape[-1]
796cell   = DYN_psol_beg.shape[-1]
797klev   = klevp1 - 1
798
799echo ( 'Layer thickness (pressure)' )
800DYN_dpres_beg = xr.DataArray ( np.empty( (klev, cell )), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) )
801DYN_dpres_end = xr.DataArray ( np.empty( (klev, cell )), dims = ('sigs', 'cell'), coords=(np.arange(klev), np.arange(cell) ) )
802   
803for 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,:]
806
807echo ( 'Vertical and horizontal integral, and sum of liquid, solid and vapor water phases' )
808if 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 :
813        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' )
814        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' )
815if ICO :
816    try :
817        DYN_wat_beg = (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s'])[..., DYN_beg_keysort]
818        DYN_wat_end = (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s'])[..., DYN_beg_keysort]
819       
820    except : 
821        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]
822        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       
824    DYN_wat_beg = DYN_wat_beg.rename ( {'lev':'sigs'} )
825    DYN_wat_beg = DYN_wat_end.rename ( {'lev':'sigs'} )
826   
827echo ( 'Compute water content : vertical and horizontal integral' )
828DYN_mas_wat_beg = DYN_stock_int ( (DYN_dpres_beg * DYN_wat_beg).sum (dim='sigs') ) / Grav
829DYN_mas_wat_end = DYN_stock_int ( (DYN_dpres_end * DYN_wat_end).sum (dim='sigs') ) / Grav
830
831echo ( 'Variation of water content' )
832dDYN_mas_wat = DYN_mas_wat_end - DYN_mas_wat_beg
833
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
837echo ( f'\nChange of atmosphere water content (dynamics)  -- {Title} ' )
838echo ( '------------------------------------------------------------------------------------' )
839echo ( 'DYN_mas_beg = {:12.6e} kg | DYN_mas_end = {:12.6e} kg'.format (DYN_mas_wat_beg, DYN_mas_wat_end) )
840prtFlux ( 'dMass (atm)  ', dDYN_mas_wat, 'e', True )
841
842ATM_sno_beg = d_ATM_beg['SNOW01']*d_ATM_beg['FTER'] + d_ATM_beg['SNOW02']*d_ATM_beg['FLIC'] + \
843              d_ATM_beg['SNOW03']*d_ATM_beg['FOCE'] + d_ATM_beg['SNOW04']*d_ATM_beg['FSIC']
844ATM_sno_end = d_ATM_end['SNOW01']*d_ATM_end['FTER'] + d_ATM_end['SNOW02']*d_ATM_end['FLIC'] + \
845              d_ATM_end['SNOW03']*d_ATM_end['FOCE'] + d_ATM_end['SNOW04']*d_ATM_end['FSIC']
846
847ATM_qs_beg  = d_ATM_beg['QS01']  *d_ATM_beg['FTER'] + d_ATM_beg['QS02']  *d_ATM_beg['FLIC'] + \
848              d_ATM_beg['QS03']  *d_ATM_beg['FOCE'] + d_ATM_beg['QS04']  *d_ATM_beg['FSIC']
849ATM_qs_end  = d_ATM_end['QS01']  *d_ATM_end['FTER'] + d_ATM_end['QS02']  *d_ATM_end['FLIC'] + \
850              d_ATM_end['QS03']  *d_ATM_end['FOCE'] + d_ATM_end['QS04']  *d_ATM_end['FSIC']
851
852ATM_qsol_beg = d_ATM_beg['QSOL']
853ATM_qsol_end = d_ATM_end['QSOL']
854
855LIC_sno_beg      = d_ATM_beg['SNOW02']*d_ATM_beg['FLIC']
856LIC_sno_end      = d_ATM_end['SNOW02']*d_ATM_end['FLIC']
857
858LIC_runlic0_beg  = d_ATM_beg['RUNOFFLIC0']
859LIC_runlic0_end  = d_ATM_end['RUNOFFLIC0']
860
861ATM_qs01_beg  = d_ATM_beg['QS01'] * d_ATM_beg['FTER']
862ATM_qs01_end  = d_ATM_end['QS01'] * d_ATM_end['FTER']
863ATM_qs02_beg  = d_ATM_beg['QS02'] * d_ATM_beg['FLIC']
864ATM_qs02_end  = d_ATM_end['QS02'] * d_ATM_end['FLIC']
865ATM_qs03_beg  = d_ATM_beg['QS03'] * d_ATM_beg['FOCE']
866ATM_qs03_end  = d_ATM_end['QS03'] * d_ATM_end['FOCE']
867ATM_qs04_beg  = d_ATM_beg['QS04'] * d_ATM_beg['FSIC']
868ATM_qs04_end  = d_ATM_end['QS04'] * d_ATM_end['FSIC'] 
869
870if 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]
889   
890LIC_qs_beg = ATM_qs02_beg
891LIC_qs_end = ATM_qs02_end
892
893ATM_mas_sno_beg     = DYN_stock_int ( ATM_sno_beg  )
894ATM_mas_sno_end     = DYN_stock_int ( ATM_sno_end  )
895
896ATM_mas_qs_beg      = DYN_stock_int ( ATM_qs_beg   )
897ATM_mas_qs_end      = DYN_stock_int ( ATM_qs_end   )
898ATM_mas_qsol_beg    = DYN_stock_int ( ATM_qsol_beg )
899ATM_mas_qsol_end    = DYN_stock_int ( ATM_qsol_end )
900ATM_mas_qs01_beg    = DYN_stock_int ( ATM_qs01_beg )
901ATM_mas_qs01_end    = DYN_stock_int ( ATM_qs01_end )
902ATM_mas_qs02_beg    = DYN_stock_int ( ATM_qs02_beg )
903ATM_mas_qs02_end    = DYN_stock_int ( ATM_qs02_end )
904ATM_mas_qs03_beg    = DYN_stock_int ( ATM_qs03_beg )
905ATM_mas_qs03_end    = DYN_stock_int ( ATM_qs03_end )
906ATM_mas_qs04_beg    = DYN_stock_int ( ATM_qs04_beg )
907ATM_mas_qs04_end    = DYN_stock_int ( ATM_qs04_end )
908
909LIC_mas_sno_beg     = DYN_stock_int ( LIC_sno_beg     )
910LIC_mas_sno_end     = DYN_stock_int ( LIC_sno_end     )
911LIC_mas_runlic0_beg = DYN_stock_int ( LIC_runlic0_beg )
912LIC_mas_runlic0_end = DYN_stock_int ( LIC_runlic0_end )
913
914LIC_mas_qs_beg  = ATM_mas_qs02_beg
915LIC_mas_qs_end  = ATM_mas_qs02_end
916LIC_mas_wat_beg = LIC_mas_qs_beg + LIC_mas_sno_beg
917LIC_mas_wat_end = LIC_mas_qs_end + LIC_mas_sno_end
918
919dATM_mas_sno  = ATM_mas_sno_end  - ATM_mas_sno_beg
920dATM_mas_qs   = ATM_mas_qs_end   - ATM_mas_qs_beg
921dATM_mas_qsol = ATM_mas_qsol_end - ATM_mas_qsol_beg
922
923dATM_mas_qs01 = ATM_mas_qs01_end - ATM_mas_qs01_beg
924dATM_mas_qs02 = ATM_mas_qs02_end - ATM_mas_qs02_beg
925dATM_mas_qs03 = ATM_mas_qs03_end - ATM_mas_qs03_beg
926dATM_mas_qs04 = ATM_mas_qs04_end - ATM_mas_qs04_beg
927
928dLIC_mas_qs      = LIC_mas_qs_end      - LIC_mas_qs_beg
929dLIC_mas_sno     = LIC_mas_sno_end     - LIC_mas_sno_beg
930dLIC_mas_runlic0 = LIC_mas_runlic0_end - LIC_mas_runlic0_beg
931dLIC_mas_wat     = dLIC_mas_qs + dLIC_mas_sno
932
933echo ( f'\nChange of atmosphere snow content (Land ice) -- {Title} ' )
934echo ( '------------------------------------------------------------------------------------' )
935echo ( 'ATM_mas_sno_beg  = {:12.6e} kg | ATM_mas_sno_end  = {:12.6e} kg'.format (ATM_mas_sno_beg, ATM_mas_sno_end) )
936prtFlux ( 'dMass (neige atm) ', dATM_mas_sno  , 'e', True  )
937
938echo ( f'\nChange of soil humidity content -- {Title} ' )
939echo ( '------------------------------------------------------------------------------------' )
940echo ( 'ATM_mas_qs_beg  = {:12.6e} kg | ATM_mas_qs_end  = {:12.6e} kg'.format (ATM_mas_qs_beg, ATM_mas_qs_end) )
941prtFlux ( 'dMass (neige atm) ', dATM_mas_qs, 'e', True )
942
943echo ( f'\nChange of atmosphere water+snow content -- {Title}  ' )
944echo ( '------------------------------------------------------------------------------------' )
945prtFlux ( 'dMass (eau + neige atm) ', dDYN_mas_wat + dATM_mas_sno , 'e', True)
946
947echo ( '\n====================================================================================' )
948echo ( f'-- SRF changes  -- {Title} ' )
949
950if Routing == 'SIMPLE' :
951    RUN_mas_wat_fast_beg   = ONE_stock_int ( d_RUN_beg ['fast_reservoir']   )
952    RUN_mas_wat_slow_beg   = ONE_stock_int ( d_RUN_beg ['slow_reservoir']   )
953    RUN_mas_wat_stream_beg = ONE_stock_int ( d_RUN_beg ['stream_reservoir'] )
954   
955    RUN_mas_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  )
956    RUN_mas_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   )
957    RUN_mas_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   )
958   
959    RUN_mas_wat_fast_end   = ONE_stock_int ( d_RUN_end ['fast_reservoir']   )
960    RUN_mas_wat_slow_end   = ONE_stock_int ( d_RUN_end ['slow_reservoir']   )
961    RUN_mas_wat_stream_end = ONE_stock_int ( d_RUN_end ['stream_reservoir'] )
962   
963    RUN_mas_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  )
964    RUN_mas_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   )
965    RUN_mas_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   )
966
967if Routing == 'SECHIBA' :
968    RUN_mas_wat_fast_beg   = ONE_stock_int ( d_SRF_beg ['fastres']   )
969    RUN_mas_wat_slow_beg   = ONE_stock_int ( d_SRF_beg ['slowres']   )
970    RUN_mas_wat_stream_beg = ONE_stock_int ( d_SRF_beg ['streamres'] )
971    RUN_mas_wat_flood_beg  = ONE_stock_int ( d_SRF_beg ['floodres']  )
972    RUN_mas_wat_lake_beg   = ONE_stock_int ( d_SRF_beg ['lakeres']   )
973    RUN_mas_wat_pond_beg   = ONE_stock_int ( d_SRF_beg ['pondres']   )
974   
975    RUN_mas_wat_fast_end   = ONE_stock_int ( d_SRF_end ['fastres']   )
976    RUN_mas_wat_slow_end   = ONE_stock_int ( d_SRF_end ['slowres']   )
977    RUN_mas_wat_stream_end = ONE_stock_int ( d_SRF_end ['streamres'] )
978    RUN_mas_wat_flood_end  = ONE_stock_int ( d_SRF_end ['floodres']  )
979    RUN_mas_wat_lake_end   = ONE_stock_int ( d_SRF_end ['lakeres']   )
980    RUN_mas_wat_pond_end   = ONE_stock_int ( d_SRF_end ['pondres']   )
981
982RUN_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
984RUN_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
986
987dRUN_mas_wat_fast   = RUN_mas_wat_fast_end   - RUN_mas_wat_fast_beg
988dRUN_mas_wat_slow   = RUN_mas_wat_slow_end   - RUN_mas_wat_slow_beg
989dRUN_mas_wat_stream = RUN_mas_wat_stream_end - RUN_mas_wat_stream_beg
990dRUN_mas_wat_flood  = RUN_mas_wat_flood_end  - RUN_mas_wat_flood_beg
991dRUN_mas_wat_lake   = RUN_mas_wat_lake_end   - RUN_mas_wat_lake_beg
992dRUN_mas_wat_pond   = RUN_mas_wat_pond_end   - RUN_mas_wat_pond_beg
993
994dRUN_mas_wat        = RUN_mas_wat_end        - RUN_mas_wat_beg
995
996echo ( f'\nRunoff reservoirs -- {Title} ')
997echo ( f'------------------------------------------------------------------------------------' )
998echo ( 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 ' )
999echo ( 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 ' )
1000echo ( 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 ' )
1001echo ( 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 ' )
1002echo ( 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 ' )
1003echo ( 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 ' )
1004echo ( f'RUN_mas_wat_beg        = {RUN_mas_wat_beg       :12.6e} kg | RUN_mas_wat_end        = {RUN_mas_wat_end       :12.6e} kg ' )
1005
1006echo ( '------------------------------------------------------------------------------------' )
1007
1008prtFlux ( 'dMass (fast)   ', dRUN_mas_wat_fast  , 'e', True )
1009prtFlux ( 'dMass (slow)   ', dRUN_mas_wat_slow  , 'e', True )
1010prtFlux ( 'dMass (stream) ', dRUN_mas_wat_stream, 'e', True )
1011prtFlux ( 'dMass (flood)  ', dRUN_mas_wat_flood , 'e', True )
1012prtFlux ( 'dMass (lake)   ', dRUN_mas_wat_lake  , 'e', True )
1013prtFlux ( 'dMass (pond)   ', dRUN_mas_wat_pond  , 'e', True )
1014prtFlux ( 'dMass (all)    ', dRUN_mas_wat       , 'e', True )
1015
1016echo ( f'\nWater content in routing  -- {Title} ' )
1017echo ( '------------------------------------------------------------------------------------' )
1018echo ( f'RUN_mas_wat_beg = {RUN_mas_wat_end:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg' )
1019prtFlux ( 'dMass (routing) ', dRUN_mas_wat , 'e', True   )
1020
1021echo ( '\n====================================================================================' )
1022print (f'Reading SRF restart')
1023SRF_tot_watveg_beg  = d_SRF_beg['tot_watveg_beg']  ; SRF_tot_watveg_beg  = SRF_tot_watveg_beg .where (SRF_tot_watveg_beg  < 1E15, 0.)
1024SRF_tot_watsoil_beg = d_SRF_beg['tot_watsoil_beg'] ; SRF_tot_watsoil_beg = SRF_tot_watsoil_beg.where (SRF_tot_watsoil_beg < 1E15, 0.)
1025SRF_snow_beg        = d_SRF_beg['snow_beg']        ; SRF_snow_beg        = SRF_snow_beg       .where (SRF_snow_beg        < 1E15, 0.)
1026SRF_lakeres_beg     = d_SRF_beg['lakeres']         ; SRF_lakeres_beg     = SRF_lakeres_beg    .where (SRF_lakeres_beg     < 1E15, 0.)
1027
1028SRF_tot_watveg_end  = d_SRF_end['tot_watveg_beg']  ; SRF_tot_watveg_end  = SRF_tot_watveg_end .where (SRF_tot_watveg_end  < 1E15, 0.)
1029SRF_tot_watsoil_end = d_SRF_end['tot_watsoil_beg'] ; SRF_tot_watsoil_end = SRF_tot_watsoil_end.where (SRF_tot_watsoil_end < 1E15, 0.)
1030SRF_snow_end        = d_SRF_end['snow_beg']        ; SRF_snow_end        = SRF_snow_end       .where (SRF_snow_end        < 1E15, 0.)
1031SRF_lakeres_end     = d_SRF_end['lakeres']         ; SRF_lakeres_end     = SRF_lakeres_end    .where (SRF_lakeres_end     < 1E15, 0.)
1032
1033if 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')
1042 
1043if ICO :
1044    SRF_tot_watveg_beg  = SRF_tot_watveg_beg  [SRF_beg_keysort]
1045    SRF_tot_watsoil_beg = SRF_tot_watsoil_beg [SRF_beg_keysort]
1046    SRF_snow_beg        = SRF_snow_beg        [SRF_beg_keysort]
1047    SRF_lakeres_beg     = SRF_lakeres_beg     [SRF_beg_keysort]
1048    SRF_tot_watveg_end  = SRF_tot_watveg_end  [SRF_end_keysort]
1049    SRF_tot_watsoil_end = SRF_tot_watsoil_end [SRF_end_keysort]
1050    SRF_snow_end        = SRF_snow_end        [SRF_end_keysort]
1051    SRF_lakeres_end     = SRF_lakeres_end     [SRF_end_keysort]
1052
1053# Stock dSoilHum dInterce dSWE dStream dFastR dSlowR dLake dPond dFlood
1054
1055SRF_wat_beg = SRF_tot_watveg_beg + SRF_tot_watsoil_beg + SRF_snow_beg
1056SRF_wat_end = SRF_tot_watveg_end + SRF_tot_watsoil_end + SRF_snow_end
1057
1058echo ( '\n====================================================================================' )
1059print ('Computing integrals')
1060
1061print ( ' 1/8', end='' ) ; sys.stdout.flush ()
1062SRF_mas_watveg_beg   = SRF_stock_int ( SRF_tot_watveg_beg  )
1063print ( ' 2/8', end='' ) ; sys.stdout.flush ()
1064SRF_mas_watsoil_beg  = SRF_stock_int ( SRF_tot_watsoil_beg )
1065print ( ' 3/8', end='' ) ; sys.stdout.flush ()
1066SRF_mas_snow_beg     = SRF_stock_int ( SRF_snow_beg        )
1067print ( ' 4/8', end='' ) ; sys.stdout.flush ()
1068SRF_mas_lake_beg     = ONE_stock_int ( SRF_lakeres_beg     )
1069print ( ' 5/8', end='' ) ; sys.stdout.flush ()
1070
1071SRF_mas_watveg_end   = SRF_stock_int ( SRF_tot_watveg_end  )
1072print ( ' 6/8', end='' ) ; sys.stdout.flush ()
1073SRF_mas_watsoil_end  = SRF_stock_int ( SRF_tot_watsoil_end )
1074print ( ' 7/8', end='' ) ; sys.stdout.flush ()
1075SRF_mas_snow_end     = SRF_stock_int ( SRF_snow_end        )
1076print ( ' 8/8', end='' ) ; sys.stdout.flush ()
1077SRF_mas_lake_end     = ONE_stock_int ( SRF_lakeres_end     )
1078
1079print (' -- ') ; sys.stdout.flush ()
1080
1081dSRF_mas_watveg   = SRF_mas_watveg_end  - SRF_mas_watveg_beg
1082dSRF_mas_watsoil  = SRF_mas_watsoil_end - SRF_mas_watsoil_beg
1083dSRF_mas_snow     = SRF_mas_snow_end    - SRF_mas_snow_beg
1084dSRF_mas_lake     = SRF_mas_lake_end    - SRF_mas_lake_beg
1085
1086echo ( '------------------------------------------------------------------------------------' )
1087echo ( f'\nSurface reservoirs  -- {Title} ')
1088echo ( f'SRF_mas_watveg_beg   = {SRF_mas_watveg_beg :12.6e} kg | SRF_mas_watveg_end   = {SRF_mas_watveg_end :12.6e} kg ' )
1089echo ( f'SRF_mas_watsoil_beg  = {SRF_mas_watsoil_beg:12.6e} kg | SRF_mas_watsoil_end  = {SRF_mas_watsoil_end:12.6e} kg ' )
1090echo ( f'SRF_mas_snow_beg     = {SRF_mas_snow_beg   :12.6e} kg | SRF_mas_snow_end     = {SRF_mas_snow_end   :12.6e} kg ' )
1091echo ( f'SRF_mas_lake_beg     = {SRF_mas_lake_beg   :12.6e} kg | SRF_mas_lake_end     = {SRF_mas_lake_end   :12.6e} kg ' )
1092
1093prtFlux ( 'dMass (watveg) ', dSRF_mas_watveg , 'e' , True )
1094prtFlux ( 'dMass (watsoil)', dSRF_mas_watsoil, 'e' , True )
1095prtFlux ( 'dMass (snow)   ', dSRF_mas_snow   , 'e' , True )
1096prtFlux ( 'dMass (lake)   ', dSRF_mas_lake   , 'e' , True )
1097
1098SRF_mas_wat_beg = SRF_mas_watveg_beg + SRF_mas_watsoil_beg + SRF_mas_snow_beg
1099SRF_mas_wat_end = SRF_mas_watveg_end + SRF_mas_watsoil_end + SRF_mas_snow_end
1100dSRF_mas_wat    = SRF_mas_wat_end - SRF_mas_wat_beg
1101
1102echo ( '------------------------------------------------------------------------------------' )
1103echo ( f'Water content in surface  -- {Title} ' )
1104echo ( f'SRF_mas_wat_beg   = {SRF_mas_wat_beg:12.6e} kg | SRF_mas_wat_end  = {SRF_mas_wat_end:12.6e} kg ')
1105prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', True )
1106
1107echo ( '------------------------------------------------------------------------------------' )
1108echo ( 'Water content in  ATM + SRF + RUN + LAKE' )
1109echo ( 'mas_wat_beg = {:12.6e} kg | mas_wat_end = {:12.6e} kg '.
1110           format (DYN_mas_wat_beg + ATM_mas_sno_beg + RUN_mas_wat_beg + SRF_mas_wat_beg + SRF_mas_lake_beg ,
1111                   DYN_mas_wat_end + ATM_mas_sno_end + RUN_mas_wat_end + SRF_mas_wat_end + SRF_mas_lake_end ) )
1112prtFlux ( 'dMass (water atm+srf+run+lake)', dDYN_mas_wat + dATM_mas_sno + dRUN_mas_wat + dSRF_mas_wat + dSRF_mas_lake, 'e', True)
1113
1114echo ( '\n====================================================================================' )
1115echo ( f'-- ATM Fluxes  -- {Title} ' )
1116
1117Step=1
1118print ( 'Step {Step}', end='' ) ; Step += 1
1119if 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   
1152print ( 'Step {Step}', end='' ) ; Step += 1
1153if 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
1185print ( 'Step {Step}', end='' ) ; Step += 1
1186
1187ATM_wprecip_oce = ATM_wrain_oce + ATM_wsnow_oce
1188ATM_wprecip_ter = ATM_wrain_ter + ATM_wsnow_ter
1189ATM_wprecip_sic = ATM_wrain_sic + ATM_wsnow_sic
1190ATM_wprecip_lic = ATM_wrain_lic + ATM_wsnow_lic
1191
1192ATM_wbilo       = ATM_wbilo_oce + ATM_wbilo_sic + ATM_wbilo_ter + ATM_wbilo_lic
1193ATM_emp         = ATM_evap - ATM_precip
1194ATM_wprecip_sea = ATM_wprecip_oce + ATM_wprecip_sic
1195ATM_wsnow_sea   = ATM_wsnow_oce   + ATM_wsnow_sic
1196ATM_wrain_sea   = ATM_wrain_oce   + ATM_wrain_sic
1197ATM_wbilo_sea   = ATM_wbilo_oce   + ATM_wbilo_sic
1198ATM_wevap_sea   = ATM_wevap_sic   + ATM_wevap_oce
1199
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
1206ATM_wemp_ter = ATM_wevap_ter - ATM_wprecip_ter
1207ATM_wemp_oce = ATM_wevap_oce - ATM_wprecip_oce
1208ATM_wemp_sic = ATM_wevap_sic - ATM_wprecip_sic
1209ATM_wemp_lic = ATM_wevap_lic - ATM_wprecip_lic
1210ATM_wemp_sea = ATM_wevap_sic - ATM_wprecip_oce
1211
1212print ( 'Step {Step}', end='' ) ; Step += 1
1213if 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
1223print ( 'Step {Step}', end='' ) ; Step += 1
1224if 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
1234print ( 'Step {Step}', end='' ) ; Step += 1
1235if 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
1242print ( 'Step {Step}', end='' ) ; Step += 1
1243if 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
1250print ( 'Step {Step}', end='' ) ; Step += 1
1251SRF_transpir.attrs['units'] = d_SRF_his ['transpir'].attrs['units']
1252SRF_emp      = SRF_evap - SRF_rain - SRF_snowf ; SRF_emp.attrs['units'] = SRF_rain.attrs['units']
1253   
1254## Correcting units of SECHIBA variables
1255def mmd2SI ( Var ) :
1256    '''Change unit from mm/d or m^3/s to kg/s if needed'''
1257    if 'units' in VarT.attrs : 
1258        if VarT.attrs['units'] in ['m^3/s', 'm3/s', 'm3.s-1',] :
1259            VarT.values = VarT.values  * ATM_rho                 ;  VarT.attrs['units'] = 'kg/s'
1260        if VarT.attrs['units'] == 'mm/d' :
1261            VarT.values = VarT.values  * ATM_rho * (1e-3/86400.) ;  VarT.attrs['units'] = 'kg/s'
1262        if VarT.attrs['units'] in ['m^3', 'm3', ] :
1263            VarT.values = VarT.values  * ATM_rho                 ;  VarT.attrs['units'] = 'kg'
1264           
1265for var in [ 'runoff', 'drainage', 'riversret', 'coastalflow', 'riverflow', 'coastalflow_cpl', 'riverflow_cpl' ] :
1266    VarT = locals()['RUN_' + var]
1267    mmd2SI (VarT)
1268
1269for var in ['evap', 'snowf', 'subli', 'transpir', 'rain', 'emp' ] :
1270    VarT = locals()['SRF_' + var]
1271    mmd2SI (VarT)
1272
1273print ( 'Step {Step}', end='' ) ; Step += 1
1274RUN_input  = RUN_runoff      + RUN_drainage
1275RUN_output = RUN_coastalflow + RUN_riverflow
1276
1277print ( 'Step {Step}', end='' ) ; Step += 1
1278ATM_flx_wbilo       = ATM_flux_int ( ATM_wbilo      )
1279ATM_flx_wbilo_lic   = ATM_flux_int ( ATM_wbilo_lic  )
1280ATM_flx_wbilo_oce   = ATM_flux_int ( ATM_wbilo_oce  )
1281ATM_flx_wbilo_sea   = ATM_flux_int ( ATM_wbilo_sea  )
1282ATM_flx_wbilo_sic   = ATM_flux_int ( ATM_wbilo_sic  )
1283ATM_flx_wbilo_ter   = ATM_flux_int ( ATM_wbilo_ter  )
1284ATM_flx_calving     = ATM_flux_int ( ATM_fqcalving  )
1285ATM_flx_fqfonte     = ATM_flux_int ( ATM_fqfonte    )
1286
1287print ( 'Step {Step}', end='' ) ; Step += 1
1288ATM_flx_precip      = ATM_flux_int ( ATM_precip     )
1289ATM_flx_snowf       = ATM_flux_int ( ATM_snowf      )
1290ATM_flx_evap        = ATM_flux_int ( ATM_evap       )
1291ATM_flx_runlic      = ATM_flux_int ( ATM_runofflic  )
1292
1293print ( 'Step {Step}', end='' ) ; Step += 1
1294ATM_flx_wrain_ter    = ATM_flux_int ( ATM_wrain_ter )
1295ATM_flx_wrain_oce    = ATM_flux_int ( ATM_wrain_oce )
1296ATM_flx_wrain_lic    = ATM_flux_int ( ATM_wrain_lic )
1297ATM_flx_wrain_sic    = ATM_flux_int ( ATM_wrain_sic )
1298ATM_flx_wrain_sea    = ATM_flux_int ( ATM_wrain_sea )
1299
1300ATM_flx_wsnow_ter    = ATM_flux_int ( ATM_wsnow_ter )
1301ATM_flx_wsnow_oce    = ATM_flux_int ( ATM_wsnow_oce )
1302ATM_flx_wsnow_lic    = ATM_flux_int ( ATM_wsnow_lic )
1303ATM_flx_wsnow_sic    = ATM_flux_int ( ATM_wsnow_sic )
1304ATM_flx_wsnow_sea    = ATM_flux_int ( ATM_wsnow_sea )
1305print ( 'Step {Step}', end='' ) ; Step += 1
1306ATM_flx_wrain_ter    = ATM_flux_int ( ATM_wrain_ter )
1307ATM_flx_wrain_oce    = ATM_flux_int ( ATM_wrain_oce )
1308ATM_flx_wrain_lic    = ATM_flux_int ( ATM_wrain_lic )
1309ATM_flx_wrain_sic    = ATM_flux_int ( ATM_wrain_sic )
1310ATM_flx_wrain_sea    = ATM_flux_int ( ATM_wrain_sea )
1311ATM_flx_wevap_ter    = ATM_flux_int ( ATM_wevap_ter )
1312ATM_flx_wevap_oce    = ATM_flux_int ( ATM_wevap_oce )
1313ATM_flx_wevap_lic    = ATM_flux_int ( ATM_wevap_lic )
1314ATM_flx_wevap_sic    = ATM_flux_int ( ATM_wevap_sic )
1315ATM_flx_wevap_sea    = ATM_flux_int ( ATM_wevap_sea )
1316ATM_flx_wprecip_lic  = ATM_flux_int ( ATM_wprecip_lic )
1317ATM_flx_wprecip_oce  = ATM_flux_int ( ATM_wprecip_oce )
1318ATM_flx_wprecip_sic  = ATM_flux_int ( ATM_wprecip_sic )
1319ATM_flx_wprecip_ter  = ATM_flux_int ( ATM_wprecip_ter )
1320ATM_flx_wprecip_sea  = ATM_flux_int ( ATM_wprecip_sea )
1321ATM_flx_wemp_lic     = ATM_flux_int ( ATM_wemp_lic )
1322ATM_flx_wemp_oce     = ATM_flux_int ( ATM_wemp_oce )
1323ATM_flx_wemp_sic     = ATM_flux_int ( ATM_wemp_sic )
1324ATM_flx_wemp_ter     = ATM_flux_int ( ATM_wemp_ter )
1325ATM_flx_wemp_sea     = ATM_flux_int ( ATM_wemp_sea )
1326
1327ATM_flx_emp          = ATM_flux_int ( ATM_emp)
1328
1329print ( 'Step {Step}', end='' ) ; Step += 1
1330RUN_flx_coastal     = ONE_flux_int ( RUN_coastalflow)
1331RUN_flx_river       = ONE_flux_int ( RUN_riverflow  )
1332RUN_flx_coastal_cpl = ONE_flux_int ( RUN_coastalflow_cpl)
1333RUN_flx_river_cpl   = ONE_flux_int ( RUN_riverflow_cpl  )
1334RUN_flx_drainage    = SRF_flux_int ( RUN_drainage   )
1335RUN_flx_riversret   = SRF_flux_int ( RUN_riversret  )
1336RUN_flx_runoff      = SRF_flux_int ( RUN_runoff     )
1337RUN_flx_input       = SRF_flux_int ( RUN_input      )
1338RUN_flx_output      = ONE_flux_int ( RUN_output     )
1339
1340print ( 'Step {Step}', end='' ) ; Step += 1
1341RUN_flx_bil    = RUN_flx_input   - RUN_flx_output
1342RUN_flx_rivcoa = RUN_flx_coastal + RUN_flx_river
1343
1344prtFlux ('wbilo_oce     ', ATM_flx_wbilo_oce  , 'f' )         
1345prtFlux ('wbilo_sic     ', ATM_flx_wbilo_sic  , 'f' )         
1346prtFlux ('wbilo_sic+oce ', ATM_flx_wbilo_sea  , 'f' )         
1347prtFlux ('wbilo_ter     ', ATM_flx_wbilo_ter  , 'f' )         
1348prtFlux ('wbilo_lic     ', ATM_flx_wbilo_lic  , 'f' )         
1349prtFlux ('Sum wbilo_*   ', ATM_flx_wbilo      , 'f', True) 
1350prtFlux ('E-P           ', ATM_flx_emp        , 'f', True) 
1351prtFlux ('calving       ', ATM_flx_calving    , 'f' ) 
1352prtFlux ('fqfonte       ', ATM_flx_fqfonte    , 'f' )       
1353prtFlux ('precip        ', ATM_flx_precip     , 'f' )       
1354prtFlux ('snowf         ', ATM_flx_snowf      , 'f' )       
1355prtFlux ('evap          ', ATM_flx_evap       , 'f' )
1356prtFlux ('runoff lic    ', ATM_flx_runlic     , 'f' )
1357
1358prtFlux ('ATM_flx_wemp_lic     ', ATM_flx_wemp_lic      , 'f' )
1359prtFlux ('ATM_flx_wemp_oce     ', ATM_flx_wemp_oce      , 'f' )
1360prtFlux ('ATM_flx_wemp_sic     ', ATM_flx_wemp_sic      , 'f' )
1361prtFlux ('ATM_flx_wemp_ter     ', ATM_flx_wemp_ter      , 'f' )
1362
1363prtFlux ('ATM_flx_wprecip_lic  ', ATM_flx_wprecip_lic   , 'f' )
1364prtFlux ('ATM_flx_wprecip_oce  ', ATM_flx_wprecip_oce   , 'f' )
1365prtFlux ('ATM_flx_wprecip_sic  ', ATM_flx_wprecip_sic   , 'f' )
1366prtFlux ('ATM_flx_wprecip_ter  ', ATM_flx_wprecip_ter   , 'f' )
1367
1368prtFlux ('ATM_flx_wrain_lic    ', ATM_flx_wrain_lic     , 'f' )
1369prtFlux ('ATM_flx_wrain_oce    ', ATM_flx_wrain_oce     , 'f' )
1370prtFlux ('ATM_flx_wrain_sic    ', ATM_flx_wrain_sic     , 'f' )
1371prtFlux ('ATM_flx_wrain_ter    ', ATM_flx_wrain_ter     , 'f' )
1372
1373prtFlux ('ATM_flx_wsnow_lic    ', ATM_flx_wsnow_lic     , 'f' )
1374prtFlux ('ATM_flx_wsnow_oce    ', ATM_flx_wsnow_oce     , 'f' )
1375prtFlux ('ATM_flx_wsnow_sic    ', ATM_flx_wsnow_sic     , 'f' )
1376prtFlux ('ATM_flx_wsnow_ter    ', ATM_flx_wsnow_ter     , 'f' )
1377
1378prtFlux ('ATM_flx_wevap_lic    ', ATM_flx_wevap_lic     , 'f' )
1379prtFlux ('ATM_flx_wevap_oce    ', ATM_flx_wevap_oce     , 'f' )
1380prtFlux ('ATM_flx_wevap_sic    ', ATM_flx_wevap_sic     , 'f' )
1381prtFlux ('ATM_flx_wevap_ter    ', ATM_flx_wevap_ter     , 'f' )
1382
1383prtFlux ('ATM_flx_wevap_lic    ', ATM_flx_wevap_lic     , 'f' )
1384prtFlux ('ATM_flx_wevap_oce    ', ATM_flx_wevap_oce     , 'f' )
1385prtFlux ('ATM_flx_wevap_sic    ', ATM_flx_wevap_sic     , 'f' )
1386prtFlux ('ATM_flx_wevap_ter    ', ATM_flx_wevap_ter     , 'f' )
1387
1388echo ( '\n====================================================================================' )
1389echo ( f'-- RUNOFF Fluxes  -- {Title} ' )
1390prtFlux ('coastalflow   ', RUN_flx_coastal    , 'f' ) 
1391prtFlux ('riverflow     ', RUN_flx_river      , 'f' )       
1392prtFlux ('coastal_cpl   ', RUN_flx_coastal_cpl, 'f' ) 
1393prtFlux ('riverf_cpl    ', RUN_flx_river_cpl  , 'f' )   
1394prtFlux ('river+coastal ', RUN_flx_rivcoa     , 'f' )   
1395prtFlux ('drainage      ', RUN_flx_drainage   , 'f' )   
1396prtFlux ('riversret     ', RUN_flx_riversret  , 'f' )   
1397prtFlux ('runoff        ', RUN_flx_runoff     , 'f' )   
1398prtFlux ('river in      ', RUN_flx_input      , 'f' )   
1399prtFlux ('river out     ', RUN_flx_output     , 'f' )   
1400prtFlux ('river bil     ', RUN_flx_bil        , 'f' )   
1401
1402ATM_flx_budget = -ATM_flx_wbilo + ATM_flx_calving + ATM_flx_runlic #+ ATM_flx_fqfonte + RUN_flx_river
1403echo ('')
1404#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 ))
1405
1406#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 ))
1407
1408ATM_flx_toSRF = -ATM_flx_wbilo_ter
1409
1410echo (' ')
1411echo ( '\n====================================================================================' )
1412echo ( f'--  Atmosphere  -- {Title} ' )
1413echo ( f'Mass begin = {DYN_mas_wat_beg:12.6e} kg | Mass end = {DYN_mas_wat_end:12.6e} kg' )
1414prtFlux ( 'dmass (atm)  = ', dDYN_mas_wat , 'e', True )
1415prtFlux ( 'Sum wbilo_*  = ', ATM_flx_wbilo, 'e', True )
1416prtFlux ( 'E-P          = ', ATM_flx_emp  , 'e', True )
1417echo ( ' ' )
1418prtFlux ( 'Water loss atm', ATM_flx_wbilo - dDYN_mas_wat, 'f', True )
1419echo ( 'Water loss atm = {:12.3e} (rel)  '.format ( (ATM_flx_wbilo - dDYN_mas_wat)/dDYN_mas_wat  ) )
1420
1421echo (' ')
1422prtFlux ( 'Water loss atm', ATM_flx_emp  - dDYN_mas_wat , 'f', True )
1423echo ( 'Water loss atm = {:12.3e} (rel)  '.format ( (ATM_flx_emp-dDYN_mas_wat)/dDYN_mas_wat  ) )
1424echo (' ')
1425
1426echo (' ')
1427echo ( '\n====================================================================================' )
1428
1429ATM_flx_runofflic_lic = ATM_flux_int ( ATM_runofflic * ATM_flic )
1430
1431LIC_flx_budget1 = -ATM_flx_wemp_lic  - ATM_flx_calving - ATM_flx_fqfonte
1432LIC_flx_budget2 = -ATM_flx_wbilo_lic - ATM_flx_calving - ATM_flx_fqfonte
1433LIC_flx_budget3 = -ATM_flx_wbilo_lic - ATM_flx_runofflic_lic
1434
1435echo ( f'--  LIC  -- {Title} ' )
1436echo ( f'Mass total   begin = {LIC_mas_wat_beg    :12.6e} kg | Mass end = {LIC_mas_wat_end    :12.6e} kg' )
1437echo ( f'Mass snow    begin = {LIC_mas_sno_beg    :12.6e} kg | Mass end = {LIC_mas_sno_end    :12.6e} kg' )
1438echo ( f'Mass qs      begin = {LIC_mas_qs_beg     :12.6e} kg | Mass end = {LIC_mas_qs_end     :12.6e} kg' )
1439echo ( f'Mass runlic0 begin = {LIC_mas_runlic0_beg:12.6e} kg | Mass end = {LIC_mas_runlic0_end:12.6e} kg' )
1440prtFlux ( 'dmass (LIC sno)       ', dLIC_mas_sno          , 'f', True, width=65 )
1441prtFlux ( 'dmass (LIC qs)        ', dLIC_mas_qs           , 'e', True, width=65 )
1442prtFlux ( 'dmass (LIC wat)       ', dLIC_mas_wat          , 'f', True, width=65 )
1443prtFlux ( 'dmass (LIC runlic0)   ', dLIC_mas_runlic0      , 'e', True, width=65 )
1444prtFlux ( 'dmass (LIC total)     ', dLIC_mas_wat          , 'e', True, width=65 )
1445prtFlux ( 'LIC ATM_flx_wemp_lic  ', ATM_flx_wemp_lic      , 'f', True, width=65 )
1446prtFlux ( 'LIC ATM_flx_fqfonte   ', ATM_flx_fqfonte       , 'f', True, width=65 )
1447prtFlux ( 'LIC ATM_flx_calving   ', ATM_flx_calving       , 'f', True, width=65 )
1448prtFlux ( 'LIC ATM_flx_runofflic ', ATM_flx_runofflic_lic , 'f', True, width=65 )
1449prtFlux ( 'LIC fqfonte + calving ', ATM_flx_calving+ATM_flx_fqfonte , 'f', True, width=65 )
1450prtFlux ( 'LIC fluxes 1 (wevap - precip*frac_lic  - fqcalving - fqfonte) ', LIC_flx_budget1              , 'f', True, width=65 )
1451prtFlux ( 'LIC fluxes 2 (-wbilo_lic - fqcalving - fqfonte)               ', LIC_flx_budget2              , 'f', True, width=65 )
1452prtFlux ( 'LIC fluxes 3 (-wbilo_lic - runofflic*frac_lic)                ', LIC_flx_budget3              , 'f', True, width=65 )
1453prtFlux ( 'LIC error 1                                                   ', LIC_flx_budget1-dLIC_mas_wat , 'e', True, width=65 )
1454prtFlux ( 'LIC error 2                                                   ', LIC_flx_budget2-dLIC_mas_wat , 'e', True, width=65 )
1455prtFlux ( 'LIC error 3                                                   ', LIC_flx_budget3-dLIC_mas_wat , 'e', True, width=65 )
1456echo ( 'LIC error (wevap - precip*frac_lic  - fqcalving - fqfonte)    = {:12.4e} (rel) '.format ( (LIC_flx_budget1-dLIC_mas_wat)/dLIC_mas_wat) )
1457echo ( 'LIC error (-wbilo_lic - fqcalving - fqfonte)                  = {:12.4e} (rel) '.format ( (LIC_flx_budget2-dLIC_mas_wat)/dLIC_mas_wat) )
1458echo ( 'LIC error (-wbilo_lic - runofflic*frac_lic)                   = {:12.4e} (rel) '.format ( (LIC_flx_budget3-dLIC_mas_wat)/dLIC_mas_wat) )
1459
1460echo ( '\n====================================================================================' )
1461echo ( f'-- SECHIBA fluxes  -- {Title} ' )
1462
1463SRF_flx_rain     = SRF_flux_int ( SRF_rain     )
1464SRF_flx_evap     = SRF_flux_int ( SRF_evap     )
1465SRF_flx_snowf    = SRF_flux_int ( SRF_snowf    )
1466SRF_flx_subli    = SRF_flux_int ( SRF_subli    )
1467SRF_flx_transpir = SRF_flux_int ( SRF_transpir )
1468SRF_flx_emp      = SRF_flux_int ( SRF_emp      )
1469
1470RUN_flx_torouting  =  RUN_flx_runoff  + RUN_flx_drainage
1471RUN_flx_fromrouting = RUN_flx_coastal + RUN_flx_river
1472
1473SRF_flx_all = SRF_flx_rain + SRF_flx_snowf - SRF_flx_evap - RUN_flx_runoff - RUN_flx_drainage
1474
1475prtFlux ('rain         ', SRF_flx_rain       , 'f' )
1476prtFlux ('evap         ', SRF_flx_evap       , 'f' )
1477prtFlux ('snowf        ', SRF_flx_snowf      , 'f' )
1478prtFlux ('E-P          ', SRF_flx_emp        , 'f' )
1479prtFlux ('subli        ', SRF_flx_subli      , 'f' )
1480prtFlux ('transpir     ', SRF_flx_transpir   , 'f' )
1481prtFlux ('to routing   ', RUN_flx_torouting  , 'f' )
1482prtFlux ('budget       ', SRF_flx_all        , 'f', small=True )
1483
1484echo ( '\n------------------------------------------------------------------------------------' )
1485echo ( 'Water content in surface ' )
1486echo ( f'SRF_mas_wat_beg  = {SRF_mas_wat_beg:12.6e} kg | SRF_mas_wat_end  = {SRF_mas_wat_end:12.6e} kg ' )
1487prtFlux ( 'dMass (water srf)', dSRF_mas_wat, 'e', small=True)
1488prtFlux ( 'Error            ',  SRF_flx_all-dSRF_mas_wat, 'e', small=True )
1489echo ( 'dMass (water srf) = {:12.4e} (rel)   '.format ( (SRF_flx_all-dSRF_mas_wat)/dSRF_mas_wat) )
1490
1491echo ( '\n====================================================================================' )
1492echo ( f'-- Check ATM vs. SRF -- {Title} ' )
1493prtFlux ('E-P ATM       ', ATM_flx_wemp_ter   , 'f' )
1494prtFlux ('wbilo ter     ', ATM_flx_wbilo_ter  , 'f' )
1495prtFlux ('E-P SRF       ', SRF_flx_emp        , 'f' )
1496prtFlux ('SRF/ATM error ', ATM_flx_wbilo_ter - SRF_flx_emp, 'e', True)
1497echo ( 'SRF/ATM error {:12.3e} (rel)  '.format ( (ATM_flx_wbilo_ter - SRF_flx_emp)/SRF_flx_emp ) )
1498
1499echo ('')
1500echo ( '\n====================================================================================' )
1501echo ( f'-- RUNOFF fluxes  -- {Title} ' )
1502RUN_flx_all = RUN_flx_torouting - RUN_flx_river - RUN_flx_coastal
1503prtFlux ('runoff    ', RUN_flx_runoff      , 'f' ) 
1504prtFlux ('drainage  ', RUN_flx_drainage    , 'f' ) 
1505prtFlux ('run+drain ', RUN_flx_torouting   , 'f' ) 
1506prtFlux ('river     ', RUN_flx_river       , 'f' ) 
1507prtFlux ('coastal   ', RUN_flx_coastal     , 'f' ) 
1508prtFlux ('riv+coa   ', RUN_flx_fromrouting , 'f' ) 
1509prtFlux ('budget    ', RUN_flx_all         , 'f' , small=True)
1510
1511echo ( '\n------------------------------------------------------------------------------------' )
1512echo ( f'Water content in routing+lake  -- {Title} ' )
1513echo ( f'RUN_mas_wat_beg  = {RUN_mas_wat_beg:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg ' )
1514prtFlux ( 'dMass (routing)  ', dRUN_mas_wat+dSRF_mas_lake, 'f', small=True)
1515prtFlux ( 'Routing error    ', RUN_flx_all+dSRF_mas_lake-dRUN_mas_wat, 'e', small=True )
1516echo ( 'Routing error : {:12.3e} (rel)'.format ( (RUN_flx_all-dSRF_mas_lake-dRUN_mas_wat)/(dRUN_mas_wat+dSRF_mas_lake) ) )
1517
1518echo ( '\n------------------------------------------------------------------------------------' )
1519echo ( f'Water content in routing  -- {Title} ' )
1520echo ( f'RUN_mas_wat_beg  = {RUN_mas_wat_beg:12.6e} kg | RUN_mas_wat_end = {RUN_mas_wat_end:12.6e} kg ' )
1521prtFlux ( 'dMass (routing) ', dRUN_mas_wat, 'f', small=True  )
1522prtFlux ( 'Routing error   ', RUN_flx_all-dRUN_mas_wat, 'e', small=True)
1523echo ( 'Routing error : {:12.3e} (rel)'.format ( (RUN_flx_all-dRUN_mas_wat)/dRUN_mas_wat ) )
1524
1525echo ( ' ' )
1526echo ( f'{Title = }' )
1527
1528echo ( 'SVN Information' )
1529for clef in SVN.keys () : echo ( SVN[clef] )
Note: See TracBrowser for help on using the repository browser.