source: TOOLS/WATER_BUDGET/CPL_waterbudget.py @ 6665

Last change on this file since 6665 was 6665, checked in by omamce, 9 months ago

O.M. : WATER BUDGET

Improved code with pylint analysis
  • Property svn:keywords set to Date Revision HeadURL Author Id
File size: 48.2 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
15#  $Author$
16#  $Date$
17#  $Revision$
18#  $Id$
19#  $HeadURL$
20
21# SVN Information
22SVN = {
23    'Author'  : "$Author$",
24    'Date'    : "$Date$",
25    'Revision': "$Revision$",
26    'Id'      : "$Id$",
27    'HeadURL' : "$HeadUrl:  $"
28    }
29###
30## Import system modules
31import sys, os, shutil, subprocess, platform
32import numpy as np
33import configparser, re
34from pathlib import Path
35
36###
37## Import system modules
38import sys, os, shutil#, subprocess, platform
39import configparser, re
40
41## Import needed scientific modules
42import numpy as np, xarray as xr
43
44## Import local modules
45import WaterUtils as wu
46import libIGCM_sys
47import nemo, lmdz
48
49from WaterUtils import RA, GRAV, ICE_RHO_ICE, ICE_RHO_SNO, OCE_RHO_LIQ, \
50                       ATM_RHO, SRF_RHO, RUN_RHO, ICE_RHO_PND, YEAR_LENGTH
51
52## Creates parser for reading .ini input file
53## -------------------------------------------
54config = configparser.ConfigParser ( interpolation=configparser.ExtendedInterpolation() )
55config.optionxform = str # To keep capitals
56
57## Experiment parameters
58## ---------------------
59ATM=None ; ATM_HIS='latlon' ; SRF_HIS='latlon' ; RUN_HIS='latlon' ; ORCA=None ; NEMO=None ; OCE_relax=False
60OCE_icb=False ; Coupled=False ; Routing=None ; TestInterp=None
61TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=None
62YearBegin=None ; YearEnd=None ; DateBegin=None ; DateEnd=None
63
64Timer=False ; Debug=False
65##
66ARCHIVE=None ; STORAGE=None ; SCRATCHDIR=None ; R_IN=None ; rebuild=None ; TmpDir=None
67FileDir=None ; FileOut=None ; R_OUT=None ; R_FIG=None ; R_FIGR=None ; R_BUFR=None ; R_SAVE=None
68R_BUF_KSH=None ; REBUILD_DIR=None ; POST_DIR=None ; L_EXP=None
69
70dir_ATM_his=None ; dir_SRF_his=None ; dir_OCE_his=None ; dir_ICE_his=None
71FileCommon=None ; file_ATM_his=None ; file_SRF_his=None ; file_RUN_his=None
72file_OCE_his=None ;  file_ICE_his=None ; file_OCE_sca=None ; file_OCE_srf=None
73tar_restart_beg=None ; tar_restart_end=None ; file_ATM_beg=None ; file_ATM_end=None ; file_DYN_beg=None
74file_DYN_end=None ; file_SRF_beg=None ; file_SRF_end=None
75file_RUN_beg=None ; file_RUN_end=None ; file_RUN_end=None ; file_OCE_beg=None
76file_ICE_beg=None ; file_OCE_beg=None
77file_OCE_end=None ; file_ICE_beg=None ; file_OCE_end=None ; file_ICE_end=None
78TarRestartDate_beg=None ; TarRestartDate_end=None
79file_DYN_aire=None
80tar_restart_beg_ATM=None ; tar_restart_beg_DYN=None ; tar_restart_beg_SRF=None
81tar_restart_beg_RUN=None ; tar_restart_beg_OCE=None ; tar_restart_beg_ICE=None
82tar_restart_end_ATM=None ; tar_restart_end_DYN=None ; tar_restart_end_SRF=None
83tar_restart_end_RUN=None ; tar_restart_end_OCE=None ; tar_restart_end_ICE=None
84ContinueOnError=False ; ErrorCount=0
85FileDirOCE=None ; FileDirATM=None ; FileDirICE=None ; FileDirSRF=None ; FileDirRUN=None
86
87##
88## Precision of history file reading
89## ---------------------------------
90# Default is float (full precision). Degrade the precision by using np.float32
91# Restart files are always read at the full precision
92readPrec=float
93
94## Read command line arguments
95## ---------------------------
96print ( "Name of Python script:", sys.argv[0] )
97IniFile = sys.argv[1]
98
99# Test existence of IniFile
100if not os.path.exists (IniFile ) :
101    raise FileExistsError ( f"File not found : {IniFile = }" )
102
103if 'full' in IniFile : FullIniFile = IniFile
104else                 : FullIniFile = 'full_' + IniFile
105   
106print ("Input file : ", IniFile )
107config.read (IniFile)
108FullIniFile = 'full_' + IniFile
109
110## Reading config.card if possible
111## -------------------------------
112ConfigCard = None
113
114if 'Experiment' in config.keys ()  : ## Read Experiment on Config file if possible
115    if 'ConfigCard' in config['Experiment'].keys () :
116        ConfigCard = config['Experiment']['ConfigCard']
117        print ( f'{ConfigCard=}' )
118
119if ConfigCard : ## Read config card if it exists
120    # Text existence of ConfigCard
121    if os.path.exists ( ConfigCard ) :
122        print ( f'Reading Config Card : {ConfigCard}' )
123        ## Creates parser for reading .ini input file
124        MyReader = configparser.ConfigParser (interpolation=configparser.ExtendedInterpolation() )
125        MyReader.optionxform = str # To keep capitals
126       
127        MyReader.read (ConfigCard)
128
129        for VarName in ['JobName', 'ExperimentName', 'SpaceName', 'LongName', 'ModelName', 'TagName'] :
130            if VarName in MyReader['UserChoices'].keys() :
131                locals()[VarName] = MyReader['UserChoices'][VarName]
132                exec  ( f'{VarName} = wu.setBool ({VarName})' )
133                exec  ( f'{VarName} = wu.setNum  ({VarName})' )
134                exec  ( f'{VarName} = wu.setNone ({VarName})' )
135                exec  ( f'wu.{VarName} = {VarName}' )
136                print ( f'    {VarName:21} set to : {locals()[VarName]:}' )
137               
138        for VarName in ['PackFrequency'] :
139            if VarName in MyReader['Post'].keys() :
140                locals()[VarName] = MyReader['Post'][VarName]
141                exec  ( f'{VarName} = wu.setBool ({VarName})' )
142                exec  ( f'{VarName} = wu.setNum  ({VarName})' )
143                exec  ( f'{VarName} = wu.setNone ({VarName})' )
144                exec  ( f'wu.{VarName} = {VarName}' )
145                print ( f'    {VarName:21} set to : {locals()[VarName]:}' )
146    else :
147        raise FileExistsError ( f"File not found : {ConfigCard = }" )
148
149   
150## Reading config file
151## -------------------
152# Each entry in the .ini file will create a Python variable with the same name
153for Section in config.keys () : 
154    print ( f'\nReading [{Section}]' )
155    for VarName in config[Section].keys() :
156        locals()[VarName] = config[Section][VarName]
157        exec  ( f'{VarName} = wu.setBool ({VarName})' )
158        exec  ( f'{VarName} = wu.setNum  ({VarName})' )
159        exec  ( f'{VarName} = wu.setNone ({VarName})' )
160        exec  ( f'wu.{VarName} = {VarName}' )
161        print ( f'    {VarName:21} set to : {locals()[VarName]}' )
162
163print ( f'\nConfig file readed : {IniFile} ' )
164
165##
166## Reading prec
167if not readPrec  :
168    readPrec = np.float64
169else :
170    if readPrec in ["float", "float64", "r8", "double", "<class 'float'>"         ] : readPrec = float
171    if readPrec in [         "float32", "r4", "single", "<class 'numpy.float32'>" ] : readPrec = np.float32
172    if readPrec in [         "float16", "r2", "half"  , "<class 'numpy.float16'>" ] : readPrec = np.float16
173   
174# Some physical constants
175## =======================
176if not RA           : RA          = wu.RA           #-- Earth Radius (m)
177if not GRAV         : GRAV        = wu.GRAV         #-- Gravity (m^2/s
178if not ICE_RHO_ICE  : ICE_RHO_ICE = wu.ICE_RHO_ICE  #-- Ice volumic mass (kg/m3) in LIM3
179if not ICE_RHO_SNO  : ICE_RHO_SNO = wu.ICE_RHO_SNO  #-- Snow volumic mass (kg/m3) in LIM3
180if not OCE_RHO_LIQ  : OCE_RHO_LIQ = wu.OCE_RHO_LIQ  #-- Ocean water volumic mass (kg/m3) in NEMO
181if not ATM_RHO      : ATM_RHO     = wu.ATM_RHO      #-- Water volumic mass in atmosphere (kg/m^3)
182if not SRF_RHO      : SRF_RHO     = wu.SRF_RHO      #-- Water volumic mass in surface reservoir (kg/m^3)
183if not RUN_RHO      : RUN_RHO     = wu.RUN_RHO      #-- Water volumic mass of rivers (kg/m^3)
184if not ICE_RHO_PND  : ICE_RHO_PND = wu.ICE_RHO_PND  #-- Water volumic mass in ice ponds (kg/m^3)
185if not YEAR_LENGTH   : YEAR_LENGTH  = wu.YEAR_LENGTH   #-- Year length (s)
186
187if not 'Files'   in config.keys () : config['Files']   = {}
188if not 'Physics' in config.keys () : config['Physics'] = {}
189if not 'Config'  in config.keys () : config['Physics'] = {}
190   
191
192config['Physics'].update ( { 'RA':str(RA), 'GRAV':str(GRAV), 'ICE_RHO_ICE':str(ICE_RHO_ICE), 'ICE_RHO_SNO':str(ICE_RHO_SNO),
193                      'OCE_RHO_LIQ':str(OCE_RHO_LIQ), 'ATM_RHO':str(ATM_RHO), 'SRF_RHO':str(SRF_RHO), 'RUN_RHO':str(RUN_RHO)} )
194
195config['Config'].update ( { 'ContinueOnError':str(ContinueOnError), 'TestInterp':str(TestInterp), 'readPrec':str(readPrec),
196                            'Debug':str(Debug), 'Timer':str(Timer) } )
197
198## --------------------------
199ICO  = ( 'ICO' in wu.ATM )
200LMDZ = ( 'LMD' in wu.ATM )
201
202mm = libIGCM_sys.config ( TagName=TagName, SpaceName=SpaceName, ExperimentName=ExperimentName, JobName=JobName, User=User, Group=Group,
203             ARCHIVE=ARCHIVE, SCRATCHDIR=SCRATCHDIR, STORAGE=STORAGE, R_IN=R_IN, R_OUT=R_OUT, R_FIG=R_FIG, rebuild=rebuild, TmpDir=TmpDir, 
204             R_SAVE=R_SAVE, R_FIGR=R_FIGR, R_BUFR=R_BUFR, R_BUF_KSH=R_BUF_KSH, REBUILD_DIR=REBUILD_DIR, POST_DIR=POST_DIR, L_EXP=L_EXP )
205globals().update(mm)
206
207config['Files']['TmpDir'] = TmpDir
208if not 'libIGCM'  in config.keys () : config['libIGCM'] = {}
209config['libIGCM'].update ( { 'ARCHIVE':str(ARCHIVE), 'STORAGE':str(STORAGE), 'TmpDir':str(TmpDir), 'R_IN':str(R_IN), 'rebuild':str(rebuild) } )
210
211## Debuging and timer
212Timer = wu.functools.partial (wu.Timer, debug=Debug, timer=Timer)
213
214## Defines begining and end of experiment
215## --------------------------------------
216if not DateBegin :
217    DateBegin = f'{YearBegin}0101'
218    config['Experiment']['DateBegin'] = str(DateBegin)
219else :
220    YearBegin, MonthBegin, DayBegin = wu.SplitDate ( DateBegin )
221    DateBegin = wu.FormatToGregorian (DateBegin)
222    config['Experiment']['YearBegin'] = str(YearBegin)
223
224if not DateEnd :
225    DateEnd   = f'{YearEnd}1231'
226    config['Experiment']['DateEnd'] = str(DateEnd)
227else :
228    YearEnd, MonthEnd, DayEnd = wu.SplitDate ( DateEnd )
229    DateEnd   = wu.FormatToGregorian (DateEnd)
230    config['Experiment']['DateEnd'] = str(DateEnd)
231
232if not PackFrequency : 
233    PackFrequency = YearEnd - YearBegin + 1
234    config['Experiment']['PackFrequency']   = f'{PackFrequency}'
235
236if type ( PackFrequency ) == str : 
237    if 'Y' in PackFrequency : PackFrequency = PackFrequency.replace ( 'Y', '') 
238    if 'M' in PackFrequency : PackFrequency = PackFrequency.replace ( 'M', '')
239    PackFrequency = int ( PackFrequency )
240   
241## Output file with water budget diagnostics
242## -----------------------------------------
243if not FileOut :
244    FileOut = f'CPL_waterbudget_{JobName}_{YearBegin}_{YearEnd}'
245    if ICO : 
246        if ATM_HIS == 'latlon' : FileOut = f'{FileOut}_LATLON'
247        if ATM_HIS == 'ico'    : FileOut = f'{FileOut}_ICO'
248    if readPrec == np.float32  : FileOut = f'{FileOut}_float32'
249    FileOut = f'{FileOut}.out'
250
251config['Files']['FileOut'] = FileOut
252
253f_out = open ( FileOut, mode = 'w' )
254   
255## Useful functions
256## ----------------
257if repr(readPrec) == "<class 'numpy.float64'>" or readPrec == float :
258    def rprec (ptab) :
259        '''This version does nothing
260
261        rprec may be use to reduce floating precision when reading history files
262        '''
263        return ptab
264else                 :
265    def rprec (ptab) :
266        '''Returns float with a different precision'''
267        return ptab.astype(readPrec).astype(float)
268   
269def kg2Sv    (val, rho=ATM_RHO) :
270    '''From kg to Sverdrup'''
271    return val/dtime_sec*1.0e-6/rho
272
273def kg2myear (val, rho=ATM_RHO) :
274    '''From kg to m/year'''
275    return val/ATM_aire_sea_tot/rho/NbYear
276
277def var2prt (var, small=False, rho=ATM_RHO) :
278    '''Formats value for printing'''
279    if small :  return var , kg2Sv(var, rho=rho)*1000., kg2myear(var, rho=rho)*1000
280    else     :  return var , kg2Sv(var, rho=rho)      , kg2myear(var, rho=rho)
281
282def prtFlux (Desc, var, Form='F', small=False, rho=ATM_RHO, width=15) :
283    '''Pretty print of formattd value'''
284    if small :
285        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} mSv | {:12.4f} mm/year "
286        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} mSv | {:12.4e} mm/year "
287    else :
288        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} Sv  | {:12.4f} m/year  "
289        if Form in ['e', 'E'] : ff=" {:14.6e} kg | {:12.4e} Sv  | {:12.4e} m/year  "
290    echo ( (' {:>{width}} = ' +ff).format (Desc, *var2prt(var, small=small, rho=rho), width=width ) )
291    return None
292
293def echo (string, end='\n') :
294    '''Function to print to stdout *and* output file'''
295    print ( str(string), end=end  )
296    sys.stdout.flush ()
297    f_out.write ( str(string) + end )
298    f_out.flush ()
299    return None
300
301echo ( f'{ContinueOnError = }' )
302echo ( f'{readPrec        = }' )
303
304echo ( f'{JobName         = }' )
305echo ( f'{ConfigCard      = }' )
306echo ( f'{libIGCM         = }' )     
307echo ( f'{User            = }' )       
308echo ( f'{Group           = }' )       
309echo ( f'{Freq            = }' )       
310echo ( f'{YearBegin       = }' )     
311echo ( f'{YearEnd         = }' )     
312echo ( f'{DateBegin       = }' )
313echo ( f'{DateEnd         = }' )
314echo ( f'{PackFrequency   = }' ) 
315echo ( f'{ATM             = }' )       
316echo ( f'{Routing         = }' )       
317echo ( f'{ORCA            = }' )     
318echo ( f'{NEMO            = }' )     
319echo ( f'{Coupled         = }' )     
320echo ( f'{ATM_HIS         = }' )     
321echo ( f'{SRF_HIS         = }' )     
322echo ( f'{RUN_HIS         = }' )
323
324## Set libIGCM directories
325## -----------------------
326if not R_OUT       : R_OUT       = os.path.join ( ARCHIVE   , 'IGCM_OUT' )
327if not R_BUF       : R_BUF       = os.path.join ( SCRATCHDIR, 'IGCM_OUT' )
328if not L_EXP       : L_EXP       = os.path.join (TagName, SpaceName, ExperimentName, JobName)
329if not R_SAVE      : R_SAVE      = os.path.join ( R_OUT, L_EXP )
330if not R_BUFR      : R_BUFR      = os.path.join ( R_BUF, L_EXP )
331if not POST_DIR    : POST_DIR    = os.path.join ( R_BUFR, 'Out' )
332if not REBUILD_DIR : REBUILD_DIR = os.path.join ( R_BUFR, 'REBUILD' )
333if not R_BUF_KSH   : R_BUF_KSH   = os.path.join ( R_BUFR, 'Out' )
334if not R_FIGR      : R_FIGR      = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP )
335
336config['libIGCM'].update ( { 'R_OUT':R_OUT, 'R_BUF':R_BUF, 'L_EXP':L_EXP, 'R_BUFR':R_BUFR, 'R_SAVE':R_SAVE, 'POST_DIR':POST_DIR,
337                             'REBUILD_DIR':REBUILD_DIR, 'R_BUF_KSH':R_BUF_KSH, 'R_FIGR':R_FIGR, 'rebuild':rebuild,
338                             'YEAR_LENGTH':str(YEAR_LENGTH)} )
339
340## Set directory to extract files
341## ------------------------------
342if not FileDir : FileDir = os.path.join ( TmpDir, f'WATER_{JobName}' )
343config['Files']['FileDir'] = FileDir
344
345if not os.path.isdir ( FileDir ) : os.makedirs ( FileDir )
346
347##- Set directories to rebuild ocean and ice restart files
348if not FileDirOCE : FileDirOCE = os.path.join ( FileDir, 'OCE' )
349if not FileDirICE : FileDirICE = os.path.join ( FileDir, 'ICE' )
350if not os.path.exists ( FileDirOCE ) : os.mkdir ( FileDirOCE )
351if not os.path.exists ( FileDirICE ) : os.mkdir ( FileDirICE )
352
353echo (' ')
354echo ( f'JobName     : {JobName}'    )
355echo ( f'Comment     : {Comment}'    )
356echo ( f'TmpDir      : {TmpDir}'     )
357echo ( f'FileDir     : {FileDir}'    )
358echo ( f'FileDirOCE  : {FileDirOCE}' )
359echo ( f'FileDirICE  : {FileDirICE}' )
360
361echo ( f'\nDealing with {L_EXP}'  )
362
363echo (' ')
364echo ( f'JobName   : {JobName}'   )
365echo ( f'Comment   : {Comment}'   )
366echo ( f'TmpDir    : {TmpDir}'    )
367
368echo ( f'\nDealing with {L_EXP}'  )
369
370## Creates model output directory names
371## ------------------------------------
372if Freq == "MO" : FreqDir =  os.path.join ('Output' , 'MO' )
373if Freq == "SE" : FreqDir =  os.path.join ('Analyse', 'SE' )
374if dir_ATM_his == None :
375    dir_ATM_his = os.path.join ( R_SAVE, "ATM", FreqDir )
376    config['Files']['dir_ATM_his'] = dir_ATM_his
377if SRF : 
378    if dir_SRF_his == None : 
379        dir_SRF_his = os.path.join ( R_SAVE, "SRF", FreqDir )
380        config['Files']['dir_SRF_his'] = dir_SRF_his
381if dir_OCE_his == None :
382    dir_OCE_his = os.path.join ( R_SAVE, "OCE", FreqDir )
383    config['Files']['dir_OCE_his'] = dir_OCE_his
384if dir_ICE_his == None : 
385    dir_ICE_his = os.path.join ( R_SAVE, "ICE", FreqDir )
386    config['Files']['dir_ICE_his'] = dir_ICE_his
387
388echo ( f'The analysis relies on files from the following model output directories : ' )
389echo ( f'{dir_ATM_his}' )
390echo ( f'{dir_OCE_his}' )
391echo ( f'{dir_ICE_his}' )
392if SRF : 
393    echo ( f'{dir_SRF_his}' )
394
395##-- Creates files names
396if not Period :
397    if Freq == 'MO' : Period = f'{DateBegin}_{DateEnd}_1M'
398    if Freq == 'SE' : Period = f'SE_{DateBegin}_{DateEnd}_1M'
399    config['Files']['Period'] = Period
400
401config['Files']['DateBegin'] = DateBegin
402config['Files']['DateBegin'] = DateEnd
403
404echo ( f'Period   : {Period}' )
405
406if not FileCommon :
407    FileCommon = f'{JobName}_{Period}'
408    config['Files']['FileCommon'] = FileCommon
409
410if not Title :
411    Title = f'{JobName} : {Freq} : {DateBegin} - {DateEnd}'
412    config['Files']['Title'] = Title
413echo ('\nOpen history files' )
414if not file_ATM_his :
415    if ATM_HIS == 'latlon' :
416        file_ATM_his = os.path.join ( dir_ATM_his, f'{FileCommon}_histmth.nc' )
417    if ATM_HIS == 'ico' :
418        file_ATM_his = os.path.join ( dir_ATM_his, f'{FileCommon}_histmth_ico.nc' )
419    config['Files']['file_ATM_his'] = file_ATM_his
420if SRF : 
421    if not file_SRF_his :
422        if ATM_HIS == 'latlon' :
423            file_SRF_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' )
424        if ATM_HIS == 'ico' :
425            file_SRF_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history_ico.nc' )
426        config['Files']['file_SRF_his'] = file_SRF_his
427
428    if Routing == 'SIMPLE' :
429        if file_RUN_his == None :
430            if ATM_HIS == 'latlon' :
431                file_RUN_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' )
432            if ATM_HIS == 'ico' :
433                file_RUN_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history_ico.nc' )
434            config['Files']['file_RUN_his'] = file_RUN_his
435
436echo ( f'{file_ATM_his = }' )
437if SRF : 
438    echo ( f'{file_SRF_his = }' )
439    if Routing == 'SIMPLE' : echo ( f'{file_RUN_his = }' )
440       
441d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
442if SRF : 
443    d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
444    if Routing == 'SECHIBA' : d_RUN_his = d_SRF_his
445    if Routing == 'SIMPLE'  : d_RUN_his = xr.open_dataset ( file_RUN_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
446   
447if not file_OCE_his :
448    file_OCE_his = os.path.join ( dir_OCE_his, f'{FileCommon}_grid_T.nc' )
449    file_OCE_his = file_OCE_his
450if not file_OCE_sca :   
451    file_OCE_sca = os.path.join ( dir_OCE_his, f'{FileCommon}_scalar.nc' )
452    config['Files']['file_OCE_sca'] = file_OCE_sca
453if not file_OCE_srf :   
454    file_OCE_srf = os.path.join ( dir_OCE_his, f'{FileCommon}_sbc.nc' )
455    config['Files']['file_OCE_srf'] = file_OCE_srf
456if not file_ICE_his : 
457    file_ICE_his = os.path.join ( dir_ICE_his, f'{FileCommon}_icemod.nc' )
458    config['Files']['file_ICE_his'] = file_ICE_his
459
460d_OCE_his = xr.open_dataset ( file_OCE_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
461d_OCE_sca = xr.open_dataset ( file_OCE_sca, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
462#d_OCE_srf = xr.open_dataset ( file_OCE_srf, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
463d_ICE_his = xr.open_dataset ( file_ICE_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
464if NEMO == '3.6' : d_ICE_his = d_ICE_his.rename ( {'y_grid_T':'y', 'x_grid_T':'x'} )
465
466echo ( f'{file_OCE_his = }' )
467echo ( f'{file_ICE_his = }' )
468echo ( f'{file_OCE_sca = }' )
469echo ( f'{file_OCE_srf = }' )
470
471## Compute run length
472## ------------------
473dtime = ( d_ATM_his.time_counter_bounds.max() - d_ATM_his.time_counter_bounds.min() )
474echo ('\nRun length : {:8.2f} days'.format ( (dtime/np.timedelta64(1, "D")).values ) )
475dtime_sec = (dtime/np.timedelta64(1, "s")).values.item() # Convert in seconds
476
477## Compute length of each period
478dtime_per = (d_ATM_his.time_counter_bounds[:,-1] -  d_ATM_his.time_counter_bounds[:,0] )
479echo ('\nPeriods lengths (days) : {:} days'.format ( (dtime_per/np.timedelta64(1, "D")).values ) )
480dtime_per_sec = (dtime_per/np.timedelta64(1, "s")).values # In seconds
481dtime_per_sec = xr.DataArray (dtime_per_sec, dims=["time_counter", ], coords=[d_ATM_his.time_counter,] )
482dtime_per_sec.attrs['unit'] = 's'
483
484# Number of years
485NbYear = dtime_sec / YEAR_LENGTH
486
487## Write the full configuration
488config_out = open (FullIniFile, 'w')
489config.write (config_out )
490config_out.close ()
491
492# ATM grid with cell surfaces
493if LMDZ :
494    echo ('ATM grid with cell surfaces : LMDZ')
495    ATM_lat       = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1D='cell' )
496    ATM_lon       = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1D='cell' )
497    ATM_aire      = lmdz.geo2point ( rprec (d_ATM_his ['aire'] )    [0], cumulPoles=True, dim1D='cell' )
498    ATM_fter      = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1D='cell' )
499    ATM_foce      = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1D='cell' )
500    ATM_fsic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1D='cell' )
501    ATM_flic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1D='cell' )
502    if SRF : 
503        SRF_lat       = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1D='cell' )
504        SRF_lon       = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1D='cell' )
505        SRF_aire      = lmdz.geo2point ( rprec (d_SRF_his ['Areas']) * rprec (d_SRF_his ['Contfrac']), dim1D='cell', cumulPoles=True )
506        SRF_areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas'])  ,  dim1D='cell', cumulPoles=True )
507        SRF_contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac']), dim1D='cell' )
508       
509if ICO :
510    if ATM_HIS == 'latlon' :
511        echo ( 'ATM areas and fractions on latlon grid' )
512        if 'lat_dom_out' in d_ATM_his.variables :
513            ATM_lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat_dom_out'])+0*rprec (d_ATM_his ['lon_dom_out']), dim1D='cell' )
514            ATM_lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat_dom_out'])+  rprec (d_ATM_his ['lon_dom_out']), dim1D='cell' )
515        else :
516            ATM_lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1D='cell' )
517            ATM_lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1D='cell' )
518        ATM_aire = lmdz.geo2point ( rprec (d_ATM_his ['aire'][0]).squeeze(), cumulPoles=True, dim1D='cell' )
519        ATM_fter = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1D='cell' )
520        ATM_foce = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1D='cell' )
521        ATM_fsic = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1D='cell' )
522        ATM_flic = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1D='cell' )
523
524    if ATM_HIS == 'ico' :
525        echo ( 'ATM areas and fractions on ICO grid' )
526        ATM_aire =  rprec (d_ATM_his ['aire']     [0]).squeeze()
527        ATM_lat  =  rprec (d_ATM_his ['lat']         )
528        ATM_lon  =  rprec (d_ATM_his ['lon']         )
529        ATM_fter =  rprec (d_ATM_his ['fract_ter'][0])
530        ATM_foce =  rprec (d_ATM_his ['fract_oce'][0])
531        ATM_fsic =  rprec (d_ATM_his ['fract_sic'][0])
532        ATM_flic =  rprec (d_ATM_his ['fract_lic'][0])
533
534    if SRF : 
535        if SRF_HIS == 'latlon' :
536            echo ( 'SRF areas and fractions on latlon grid' )
537            if 'lat_domain_landpoints_out' in d_SRF_his  : 
538                SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat_domain_landpoints_out'])+0*rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1D='cell' )
539                SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_domain_landpoints_out'])+  rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1D='cell' )
540            else : 
541                if 'lat_domain_landpoints_out' in d_SRF_his :
542                    SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat_dom_out'])+0*rprec (d_SRF_his ['lon_dom_out']), dim1D='cell' )
543                    SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_dom_out'])+  rprec (d_SRF_his ['lon_dom_out']), dim1D='cell' )
544                else :
545                    SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1D='cell' )
546                    SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1D='cell' )
547                   
548            SRF_areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas']   )   , dim1D='cell', cumulPoles=True )
549            SRF_areafrac  = lmdz.geo2point ( rprec (d_SRF_his ['AreaFrac'])   , dim1D='cell', cumulPoles=True )
550            SRF_contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac'])   , dim1D='cell', cumulPoles=True )
551            SRF_aire = SRF_areafrac
552 
553        if SRF_HIS == 'ico' :
554            echo ( 'SRF areas and fractions on latlon grid' )
555            SRF_lat       =  rprec (d_SRF_his ['lat']     )
556            SRF_lon       =  rprec (d_SRF_his ['lon']     )
557            SRF_areas     =  rprec (d_SRF_his ['Areas']   ) 
558            SRF_contfrac  =  rprec (d_SRF_his ['Contfrac'])
559            SRF_aire      =  SRF_areas * SRF_contfrac
560           
561ATM_fsea      = ATM_foce + ATM_fsic
562ATM_flnd      = ATM_fter + ATM_flic
563ATM_aire_fter = ATM_aire * ATM_fter
564ATM_aire_flic = ATM_aire * ATM_flic
565ATM_aire_fsic = ATM_aire * ATM_fsic
566ATM_aire_foce = ATM_aire * ATM_foce
567ATM_aire_flnd = ATM_aire * ATM_flnd
568ATM_aire_fsea = ATM_aire * ATM_fsea
569
570#SRF_aire = SRF_aire.where ( SRF_aire < 1E15, 0.)
571
572# if ICO :
573#     if wu.unDefined ('file_DYN_aire') : file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ATM+'_grid.nc' )
574#     config['Files']['file_DYN_aire'] = file_DYN_aire
575
576# if ICO :
577#     # Area on icosahedron grid
578#     d_DYN_aire = xr.open_dataset ( file_DYN_aire, decode_times=False ).squeeze()
579           
580#     DYN_lat = d_DYN_aire['lat']
581#     DYN_lon = d_DYN_aire['lon']
582
583#     DYN_aire = d_DYN_aire['aire']
584#     DYN_fsea = d_DYN_aire['fract_oce'] + d_DYN_aire['fract_sic']
585#     DYN_flnd = 1.0 - DYN_fsea
586#     DYN_fter = d_ATM_beg['FTER']
587#     DYN_flic = d_ATM_beg['FLIC']
588#     DYN_aire_fter = DYN_aire * DYN_fter
589   
590# if LMDZ :
591#     # Area on lon/lat grid
592#     DYN_aire = ATM_aire
593#     DYN_fsea = ATM_fsea
594#     DYN_flnd = ATM_flnd
595#     DYN_fter = rprec (d_ATM_beg['FTER'])
596#     DYN_flic = rprec (d_ATM_beg['FLIC'])
597#     DYN_aire_fter = DYN_aire * DYN_fter
598
599# Functions computing integrals and sum
600# def ATM_stock_int (stock) :
601#     '''Integrate (* surface) stock on atmosphere grid'''
602#     ATM_stock_int  = wu.Psum ( (stock * DYN_aire).to_masked_array().ravel() )
603#     return ATM_stock_int
604
605def ATM_flux_int (flux) :
606    '''Integrate (* time * surface) flux on atmosphere grid'''
607    ATM_flux_int  = wu.Psum ( (flux * dtime_per_sec * ATM_aire).to_masked_array().ravel() )
608    return ATM_flux_int
609
610# def SRF_stock_int (stock) :
611#     '''Integrate (* surface) stock on land grid'''
612#     ATM_stock_int  = wu.Ksum (  ( (stock * DYN_aire_fter).to_masked_array().ravel()) )
613#     return ATM_stock_int
614
615def SRF_flux_int (flux) :
616    '''Integrate (* time * surface) flux on land grid'''
617    SRF_flux_int  = wu.Psum (  (flux * dtime_per_sec * SRF_aire).to_masked_array().ravel() )
618    return SRF_flux_int
619
620# def LIC_flux_int (flux) :
621#     '''Integrate (* time * surface) flux on land ice grid'''
622#     LIC_flux_int  = wu.Psum ( (flux * dtime_per_sec * ATM_aire_flic).to_masked_array().ravel() )
623#     return LIC_flux_int
624
625# def OCE_stock_int (stock) :
626#     '''Integrate stock on ocean grid'''
627#     OCE_stock_int = np.sum (  np.sort ( (stock * OCE_aire ).to_masked_array().ravel()) )
628#     return OCE_stock_int
629
630def ONE_stock_int (stock) :
631    '''Sum stock'''
632    ONE_stock_int =  wu.Psum ( stock.to_masked_array().ravel() )
633    return ONE_stock_int
634
635def OCE_flux_int (flux) :
636    '''Integrate flux on oce grid'''
637    OCE_flux_int = wu.Psum ( (flux * OCE_aire * dtime_per_sec).to_masked_array().ravel() )
638    return OCE_flux_int
639
640def ONE_flux_int (flux) :
641    '''Integrate flux on oce grid'''
642    OCE_flux_int = wu.Psum ( (flux * dtime_per_sec).to_masked_array().ravel() )
643    return OCE_flux_int
644
645# Get mask and surfaces
646sos = d_OCE_his ['sos'][0].squeeze()
647OCE_msk = nemo.lbc_mask ( xr.where ( sos>0, 1., 0.0 ), cd_type = 'T' )
648so = sos = d_OCE_his ['sos'][0].squeeze()
649OCE_msk3 = nemo.lbc_mask ( xr.where ( so>0., 1., 0. ), cd_type = 'T', sval = 0. )
650
651# lbc_mask removes the duplicate points (periodicity and north fold)
652OCE_aire = nemo.lbc_mask ( d_OCE_his ['area'] * OCE_msk, cd_type = 'T', sval = 0.0 )
653ICE_aire = OCE_aire
654
655ATM_aire_tot = ONE_stock_int (ATM_aire)
656if SRF : 
657    SRF_aire_tot = ONE_stock_int (SRF_aire)
658OCE_aire_tot = ONE_stock_int (OCE_aire)
659ICE_aire_tot = ONE_stock_int (ICE_aire)
660
661ATM_aire_sea     = ATM_aire * ATM_fsea
662ATM_aire_sea_tot = ONE_stock_int ( ATM_aire_sea )
663
664echo ( '\n====================================================================================' )
665echo ( f'-- ATM Fluxes  -- {Title} ' )
666
667if ATM_HIS == 'latlon' :
668    echo ( ' latlon case' )
669    ATM_wbilo_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_oce']), dim1D='cell' )
670    ATM_wbilo_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_sic']), dim1D='cell' )
671    ATM_wbilo_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_ter']), dim1D='cell' )
672    ATM_wbilo_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_lic']), dim1D='cell' )
673    ATM_runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1D='cell' )
674    ATM_fqcalving   = lmdz.geo2point ( rprec (d_ATM_his ['fqcalving']), dim1D='cell' )
675    ATM_fqfonte     = lmdz.geo2point ( rprec (d_ATM_his ['fqfonte']  ), dim1D='cell' )
676    ATM_precip      = lmdz.geo2point ( rprec (d_ATM_his ['precip']   ), dim1D='cell' )
677    ATM_snowf       = lmdz.geo2point ( rprec (d_ATM_his ['snow']     ), dim1D='cell' )
678    ATM_evap        = lmdz.geo2point ( rprec (d_ATM_his ['evap']     ), dim1D='cell' )
679    ATM_wevap_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_ter']), dim1D='cell' )
680    ATM_wevap_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_oce']), dim1D='cell' )
681    ATM_wevap_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_lic']), dim1D='cell' )
682    ATM_wevap_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_sic']), dim1D='cell' )
683    ATM_wrain_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_ter']), dim1D='cell' )
684    ATM_wrain_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_oce']), dim1D='cell' )
685    ATM_wrain_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_lic']), dim1D='cell' )
686    ATM_wrain_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_sic']), dim1D='cell' )
687    ATM_wsnow_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_ter']), dim1D='cell' )
688    ATM_wsnow_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_oce']), dim1D='cell' )
689    ATM_wsnow_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_lic']), dim1D='cell' )
690    ATM_wsnow_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_sic']), dim1D='cell' )
691    ATM_runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1D='cell' )
692    echo ( f'End of LATLON case')
693   
694if ATM_HIS == 'ico' :
695    echo (' ico case')
696    ATM_wbilo_oce   = rprec (d_ATM_his ['wbilo_oce'])
697    ATM_wbilo_sic   = rprec (d_ATM_his ['wbilo_sic'])
698    ATM_wbilo_ter   = rprec (d_ATM_his ['wbilo_ter'])
699    ATM_wbilo_lic   = rprec (d_ATM_his ['wbilo_lic'])
700    ATM_runofflic   = rprec (d_ATM_his ['runofflic'])
701    ATM_fqcalving   = rprec (d_ATM_his ['fqcalving'])
702    ATM_fqfonte     = rprec (d_ATM_his ['fqfonte']  )
703    ATM_precip      = rprec (d_ATM_his ['precip']   )
704    ATM_snowf       = rprec (d_ATM_his ['snow']     )
705    ATM_evap        = rprec (d_ATM_his ['evap']     )
706    ATM_wevap_ter   = rprec (d_ATM_his ['wevap_ter'])
707    ATM_wevap_oce   = rprec (d_ATM_his ['wevap_oce'])
708    ATM_wevap_lic   = rprec (d_ATM_his ['wevap_lic'])
709    ATM_wevap_sic   = rprec (d_ATM_his ['wevap_sic'])
710    ATM_runofflic   = rprec (d_ATM_his ['runofflic'])
711    ATM_wevap_ter   = rprec (d_ATM_his ['wevap_ter'])
712    ATM_wevap_oce   = rprec (d_ATM_his ['wevap_oce'])
713    ATM_wevap_lic   = rprec (d_ATM_his ['wevap_lic'])
714    ATM_wevap_sic   = rprec (d_ATM_his ['wevap_sic'])
715    ATM_wrain_ter   = rprec (d_ATM_his ['wrain_ter'])
716    ATM_wrain_oce   = rprec (d_ATM_his ['wrain_oce'])
717    ATM_wrain_lic   = rprec (d_ATM_his ['wrain_lic'])
718    ATM_wrain_sic   = rprec (d_ATM_his ['wrain_sic'])
719    ATM_wsnow_ter   = rprec (d_ATM_his ['wsnow_ter'])
720    ATM_wsnow_oce   = rprec (d_ATM_his ['wsnow_oce'])
721    ATM_wsnow_lic   = rprec (d_ATM_his ['wsnow_lic'])
722    ATM_wsnow_sic   = rprec (d_ATM_his ['wsnow_sic'])
723    echo ( f'End of ico case ')
724
725
726   
727echo ( 'ATM wprecip_oce' )
728ATM_wprecip_oce = ATM_wrain_oce + ATM_wsnow_oce
729ATM_wprecip_ter = ATM_wrain_ter + ATM_wsnow_ter
730ATM_wprecip_sic = ATM_wrain_sic + ATM_wsnow_sic
731ATM_wprecip_lic = ATM_wrain_lic + ATM_wsnow_lic
732
733ATM_wbilo       = ATM_wbilo_oce   + ATM_wbilo_sic   + ATM_wbilo_ter   + ATM_wbilo_lic
734ATM_wevap       = ATM_wevap_oce   + ATM_wevap_sic   + ATM_wevap_ter   + ATM_wevap_lic
735ATM_wprecip     = ATM_wprecip_oce + ATM_wprecip_sic + ATM_wprecip_ter + ATM_wprecip_lic
736ATM_wsnow       = ATM_wsnow_oce   + ATM_wsnow_sic   + ATM_wsnow_ter   + ATM_wsnow_lic
737ATM_wrain       = ATM_wrain_oce   + ATM_wrain_sic   + ATM_wrain_ter   + ATM_wrain_lic
738ATM_wemp        = ATM_wevap - ATM_wprecip
739ATM_emp         = ATM_evap  - ATM_precip
740
741ATM_wprecip_sea = ATM_wprecip_oce + ATM_wprecip_sic
742ATM_wsnow_sea   = ATM_wsnow_oce   + ATM_wsnow_sic
743ATM_wrain_sea   = ATM_wrain_oce   + ATM_wrain_sic
744ATM_wbilo_sea   = ATM_wbilo_oce   + ATM_wbilo_sic
745ATM_wevap_sea   = ATM_wevap_sic   + ATM_wevap_oce
746
747ATM_wemp_ter    = ATM_wevap_ter - ATM_wprecip_ter
748ATM_wemp_oce    = ATM_wevap_oce - ATM_wprecip_oce
749ATM_wemp_sic    = ATM_wevap_sic - ATM_wprecip_sic
750ATM_wemp_lic    = ATM_wevap_lic - ATM_wprecip_lic
751ATM_wemp_sea    = ATM_wevap_sic - ATM_wprecip_oce
752
753if SRF : 
754    if RUN_HIS == 'latlon' :
755        echo ( f'RUN costalflow Grille LATLON' )
756        if TestInterp :
757            echo ( f'RUN runoff TestInterp' )
758            RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff_contfrac_interp']  )   , dim1D='cell' )
759            RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage_contfrac_interp'])   , dim1D='cell' )
760        else : 
761            echo ( f'RUN runoff' )
762            RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff']         ), dim1D='cell' )
763            RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage']       ), dim1D='cell' )
764           
765        RUN_coastalflow     = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow']    ), dim1D='cell' ) 
766        RUN_riverflow       = lmdz.geo2point ( rprec (d_RUN_his ['riverflow']      ), dim1D='cell' )
767        RUN_riversret       = lmdz.geo2point ( rprec (d_RUN_his ['riversret']      ), dim1D='cell' )
768        RUN_coastalflow_cpl = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow_cpl']), dim1D='cell' ) 
769        RUN_riverflow_cpl   = lmdz.geo2point ( rprec (d_RUN_his ['riverflow_cpl']  ), dim1D='cell' )
770       
771    if RUN_HIS == 'ico' :
772        echo ( f'RUN costalflow Grille ICO' )
773        RUN_coastalflow =  rprec (d_RUN_his ['coastalflow'])
774        RUN_riverflow   =  rprec (d_RUN_his ['riverflow']  )
775        RUN_runoff      =  rprec (d_RUN_his ['runoff']     )
776        RUN_drainage    =  rprec (d_RUN_his ['drainage']   )
777        RUN_riversret   =  rprec (d_RUN_his ['riversret']  )
778       
779        RUN_coastalflow_cpl = rprec (d_RUN_his ['coastalflow_cpl'])
780        RUN_riverflow_cpl   = rprec (d_RUN_his ['riverflow_cpl']  )
781       
782        ## Correcting units of SECHIBA variables
783    def mmd2SI ( Var ) :
784        '''Change unit from mm/d or m^3/s to kg/s if needed'''
785        if 'units' in VarT.attrs : 
786            if VarT.attrs['units'] in ['m^3/s', 'm3/s', 'm3.s-1',] :
787                VarT.values = VarT.values  * ATM_RHO                 ;  VarT.attrs['units'] = 'kg/s'
788            if VarT.attrs['units'] == 'mm/d' :
789                VarT.values = VarT.values  * ATM_RHO * (1e-3/86400.) ;  VarT.attrs['units'] = 'kg/s'
790            if VarT.attrs['units'] in ['m^3', 'm3', ] :
791                VarT.values = VarT.values  * ATM_RHO                 ;  VarT.attrs['units'] = 'kg'
792               
793    for var in [ 'runoff', 'drainage', 'riversret', 'coastalflow', 'riverflow', 'coastalflow_cpl', 'riverflow_cpl' ] :
794        VarT = locals()['RUN_' + var]
795        mmd2SI (VarT)
796       
797    #for var in ['evap', 'snowf', 'subli', 'transpir', 'rain', 'emp' ] :
798    #    VarT = locals()['SRF_' + var]
799    #    mmd2SI (VarT)
800    echo ( f'RUN input' )
801    RUN_input  = RUN_runoff      + RUN_drainage
802    RUN_output = RUN_coastalflow + RUN_riverflow
803   
804echo ( f'ATM flw_wbilo' )
805ATM_flx_wbilo       = ATM_flux_int ( ATM_wbilo      )
806ATM_flx_wevap       = ATM_flux_int ( ATM_wevap      )
807ATM_flx_wprecip     = ATM_flux_int ( ATM_wprecip    )
808ATM_flx_wsnow       = ATM_flux_int ( ATM_wsnow      )
809ATM_flx_wrain       = ATM_flux_int ( ATM_wrain      )
810ATM_flx_wemp        = ATM_flux_int ( ATM_wemp       )
811
812ATM_flx_wbilo_lic   = ATM_flux_int ( ATM_wbilo_lic  )
813ATM_flx_wbilo_oce   = ATM_flux_int ( ATM_wbilo_oce  )
814ATM_flx_wbilo_sea   = ATM_flux_int ( ATM_wbilo_sea  )
815ATM_flx_wbilo_sic   = ATM_flux_int ( ATM_wbilo_sic  )
816ATM_flx_wbilo_ter   = ATM_flux_int ( ATM_wbilo_ter  )
817ATM_flx_calving     = ATM_flux_int ( ATM_fqcalving  )
818ATM_flx_fqfonte     = ATM_flux_int ( ATM_fqfonte    )
819
820LIC_flx_calving     = ATM_flux_int ( ATM_fqcalving  )
821LIC_flx_fqfonte     = ATM_flux_int ( ATM_fqfonte    )
822
823ATM_flx_precip      = ATM_flux_int ( ATM_precip     )
824ATM_flx_snowf       = ATM_flux_int ( ATM_snowf      )
825ATM_flx_evap        = ATM_flux_int ( ATM_evap       )
826ATM_flx_runlic      = ATM_flux_int ( ATM_runofflic  )
827
828#LIC_flx_precip      = LIC_flux_int ( ATM_precip     )
829#LIC_flx_snowf       = LIC_flux_int ( ATM_snowf      )
830#LIC_flx_evap        = LIC_flux_int ( ATM_evap       )
831#LIC_flx_runlic      = LIC_flux_int ( ATM_runofflic  )
832
833ATM_flx_wrain_ter   = ATM_flux_int ( ATM_wrain_ter )
834ATM_flx_wrain_oce   = ATM_flux_int ( ATM_wrain_oce )
835ATM_flx_wrain_lic   = ATM_flux_int ( ATM_wrain_lic )
836ATM_flx_wrain_sic   = ATM_flux_int ( ATM_wrain_sic )
837ATM_flx_wrain_sea   = ATM_flux_int ( ATM_wrain_sea )
838
839ATM_flx_wsnow_ter   = ATM_flux_int ( ATM_wsnow_ter )
840ATM_flx_wsnow_oce   = ATM_flux_int ( ATM_wsnow_oce )
841ATM_flx_wsnow_lic   = ATM_flux_int ( ATM_wsnow_lic )
842ATM_flx_wsnow_sic   = ATM_flux_int ( ATM_wsnow_sic )
843ATM_flx_wsnow_sea   = ATM_flux_int ( ATM_wsnow_sea )
844
845ATM_flx_wevap_ter   = ATM_flux_int ( ATM_wevap_ter )
846ATM_flx_wevap_oce   = ATM_flux_int ( ATM_wevap_oce )
847ATM_flx_wevap_lic   = ATM_flux_int ( ATM_wevap_lic )
848ATM_flx_wevap_sic   = ATM_flux_int ( ATM_wevap_sic )
849ATM_flx_wevap_sea   = ATM_flux_int ( ATM_wevap_sea )
850ATM_flx_wprecip_lic = ATM_flux_int ( ATM_wprecip_lic )
851ATM_flx_wprecip_oce = ATM_flux_int ( ATM_wprecip_oce )
852ATM_flx_wprecip_sic = ATM_flux_int ( ATM_wprecip_sic )
853ATM_flx_wprecip_ter = ATM_flux_int ( ATM_wprecip_ter )
854ATM_flx_wprecip_sea = ATM_flux_int ( ATM_wprecip_sea )
855ATM_flx_wemp_lic    = ATM_flux_int ( ATM_wemp_lic )
856ATM_flx_wemp_oce    = ATM_flux_int ( ATM_wemp_oce )
857ATM_flx_wemp_sic    = ATM_flux_int ( ATM_wemp_sic )
858ATM_flx_wemp_ter    = ATM_flux_int ( ATM_wemp_ter )
859ATM_flx_wemp_sea    = ATM_flux_int ( ATM_wemp_sea )
860
861ATM_flx_emp          = ATM_flux_int ( ATM_emp )
862
863if SRF : 
864    RUN_flx_coastal     = ONE_flux_int ( RUN_coastalflow)
865    RUN_flx_river       = ONE_flux_int ( RUN_riverflow  )
866    RUN_flx_coastal_cpl = ONE_flux_int ( RUN_coastalflow_cpl)
867    RUN_flx_river_cpl   = ONE_flux_int ( RUN_riverflow_cpl  )
868    RUN_flx_drainage    = SRF_flux_int ( RUN_drainage   )
869    RUN_flx_riversret   = SRF_flux_int ( RUN_riversret  )
870    RUN_flx_runoff      = SRF_flux_int ( RUN_runoff     )
871    RUN_flx_input       = SRF_flux_int ( RUN_input      )
872    RUN_flx_output      = ONE_flux_int ( RUN_output     )
873   
874    RUN_flx_bil    = ONE_flux_int ( RUN_input       - RUN_output)
875    RUN_flx_rivcoa = ONE_flux_int ( RUN_coastalflow + RUN_riverflow)
876   
877prtFlux ('wbilo_oce            ', ATM_flx_wbilo_oce     , 'f' )         
878prtFlux ('wbilo_sic            ', ATM_flx_wbilo_sic     , 'f' )         
879prtFlux ('wbilo_sic+oce        ', ATM_flx_wbilo_sea     , 'f' )         
880prtFlux ('wbilo_ter            ', ATM_flx_wbilo_ter     , 'f' )         
881prtFlux ('wbilo_lic            ', ATM_flx_wbilo_lic     , 'f' )         
882prtFlux ('Sum wbilo_*          ', ATM_flx_wbilo         , 'f', True) 
883prtFlux ('E-P                  ', ATM_flx_emp           , 'f', True) 
884prtFlux ('calving              ', ATM_flx_calving       , 'f' ) 
885prtFlux ('fqfonte              ', ATM_flx_fqfonte       , 'f' )       
886prtFlux ('precip               ', ATM_flx_precip        , 'f' )       
887prtFlux ('snowf                ', ATM_flx_snowf         , 'f' )       
888prtFlux ('evap                 ', ATM_flx_evap          , 'f' )
889prtFlux ('runoff lic           ', ATM_flx_runlic        , 'f' )
890
891prtFlux ('ATM_flx_wevap*       ', ATM_flx_wevap         , 'f' )
892prtFlux ('ATM_flx_wrain*       ', ATM_flx_wrain         , 'f' )
893prtFlux ('ATM_flx_wsnow*       ', ATM_flx_wsnow         , 'f' )
894prtFlux ('ATM_flx_wprecip*     ', ATM_flx_wprecip       , 'f' )
895prtFlux ('ATM_flx_wemp*        ', ATM_flx_wemp          , 'f', True )
896
897prtFlux ('ERROR evap           ', ATM_flx_wevap   - ATM_flx_evap  , 'e', True )
898prtFlux ('ERROR precip         ', ATM_flx_wprecip - ATM_flx_precip, 'e', True )
899prtFlux ('ERROR snow           ', ATM_flx_wsnow   - ATM_flx_snowf , 'e', True )
900prtFlux ('ERROR emp            ', ATM_flx_wemp    - ATM_flx_emp   , 'e', True )
901
902if SRF : 
903    echo ( '\n====================================================================================' )
904    echo ( f'-- RUNOFF Fluxes  -- {Title} ' )
905    prtFlux ('coastalflow   ', RUN_flx_coastal    , 'f' ) 
906    prtFlux ('riverflow     ', RUN_flx_river      , 'f' )       
907    prtFlux ('coastal_cpl   ', RUN_flx_coastal_cpl, 'f' ) 
908    prtFlux ('riverf_cpl    ', RUN_flx_river_cpl  , 'f' )   
909    prtFlux ('river+coastal ', RUN_flx_rivcoa     , 'f' )   
910    prtFlux ('drainage      ', RUN_flx_drainage   , 'f' )   
911    prtFlux ('riversret     ', RUN_flx_riversret  , 'f' )   
912    prtFlux ('runoff        ', RUN_flx_runoff     , 'f' )   
913    prtFlux ('river in      ', RUN_flx_input      , 'f' )   
914    prtFlux ('river out     ', RUN_flx_output     , 'f' )   
915    prtFlux ('river bil     ', RUN_flx_bil        , 'f' ) 
916   
917echo ( '\n====================================================================================' )
918echo ( f'-- OCE Fluxes  -- {Title} ' )
919
920# Read variable and computes integral over space and time
921OCE_empmr      = rprec (d_OCE_his['wfo']     )    ; OCE_mas_empmr    = OCE_flux_int ( OCE_empmr    )
922OCE_wfob       = rprec (d_OCE_his['wfob']    )    ; OCE_mas_wfob     = OCE_flux_int ( OCE_wfob     )
923OCE_emp_oce    = rprec (d_OCE_his['emp_oce'] )    ; OCE_mas_emp_oce  = OCE_flux_int ( OCE_emp_oce  )
924OCE_emp_ice    = rprec (d_OCE_his['emp_ice'] )    ; OCE_mas_emp_ice  = OCE_flux_int ( OCE_emp_ice  )
925OCE_iceshelf   = rprec (d_OCE_his['iceshelf'])    ; OCE_mas_iceshelf = OCE_flux_int ( OCE_iceshelf )
926OCE_calving    = rprec (d_OCE_his['calving'] )    ; OCE_mas_calving  = OCE_flux_int ( OCE_calving  )
927OCE_iceberg    = rprec (d_OCE_his['iceberg'] )    ; OCE_mas_iceberg  = OCE_flux_int ( OCE_iceberg  )
928OCE_friver     = rprec (d_OCE_his['friver']  )    ; OCE_mas_friver   = OCE_flux_int ( OCE_friver   )
929OCE_runoffs    = rprec (d_OCE_his['runoffs'] )    ; OCE_mas_runoffs  = OCE_flux_int ( OCE_runoffs  )
930if NEMO == 4.0 or NEMO == 4.2 :
931    OCE_wfxice     = rprec (d_OCE_his['vfxice']) ; OCE_mas_wfxice   = OCE_flux_int ( OCE_wfxice )
932    OCE_wfxsnw     = rprec (d_OCE_his['vfxsnw']) ; OCE_mas_wfxsnw   = OCE_flux_int ( OCE_wfxsnw )
933    OCE_wfxsub     = rprec (d_OCE_his['vfxsub']) ; OCE_mas_wfxsub   = OCE_flux_int ( OCE_wfxsub )
934if NEMO == 3.6 :
935    OCE_wfxice     = rprec (d_OCE_his['vfxice'])/86400.*ICE_RHO_ICE ; OCE_mas_wfxice   = OCE_flux_int ( OCE_wfxice )
936    OCE_wfxsnw     = rprec (d_OCE_his['vfxsnw'])/86400.*ICE_RHO_SNO ; OCE_mas_wfxsnw   = OCE_flux_int ( OCE_wfxsnw )
937    OCE_wfxsub     = rprec (d_OCE_his['vfxsub'])/86400.*ICE_RHO_SNO ; OCE_mas_wfxsub   = OCE_flux_int ( OCE_wfxsub )
938# Additional checks
939OCE_evap_oce   = rprec (d_OCE_his['evap_ao_cea']) ; OCE_mas_evap_oce   = OCE_flux_int ( OCE_evap_oce )
940ICE_evap_ice   = rprec (d_OCE_his['subl_ai_cea']) ; ICE_mas_evap_ice   = OCE_flux_int ( ICE_evap_ice )
941OCE_snow_oce   = rprec (d_OCE_his['snow_ao_cea']) ; OCE_mas_snow_oce   = OCE_flux_int ( OCE_snow_oce )
942OCE_snow_ice   = rprec (d_OCE_his['snow_ai_cea']) ; OCE_mas_snow_ice   = OCE_flux_int ( OCE_snow_ice )
943OCE_rain       = rprec (d_OCE_his['rain']       ) ; OCE_mas_rain       = OCE_flux_int ( OCE_rain     )
944ICE_wfxsub_err = rprec (d_ICE_his['vfxsub_err'] ) ; ICE_mas_wfxsub_err = OCE_flux_int ( ICE_wfxsub_err )
945if NEMO == 4.0 or NEMO == 4.2 :
946    ICE_wfxpnd     = rprec (d_ICE_his['vfxpnd']    ) ; ICE_mas_wfxpnd     = OCE_flux_int ( ICE_wfxpnd     )
947    ICE_wfxsnw_sub = rprec (d_ICE_his['vfxsnw_sub']) ; ICE_mas_wfxsnw_sub = OCE_flux_int ( ICE_wfxsnw_sub )
948    ICE_wfxsnw_pre = rprec (d_ICE_his['vfxsnw_pre']) ; ICE_mas_wfxsnw_pre = OCE_flux_int ( ICE_wfxsnw_pre )
949if NEMO == 3.6 :
950    ICE_wfxpnd = 0.0 ; ICE_mas_wfxpnd = 0.0
951    ICE_wfxsnw_sub = rprec (d_ICE_his['vfxsub'])/86400.*ICE_RHO_SNO  ; ICE_mas_wfxsnw_sub = OCE_flux_int ( ICE_wfxsnw_sub )
952    ICE_wfxsnw_pre = rprec (d_ICE_his['vfxspr'])/86400.*ICE_RHO_SNO  ; ICE_mas_wfxsnw_pre = OCE_flux_int ( ICE_wfxsnw_pre     )
953
954OCE_wfcorr    = rprec (d_OCE_his['wfcorr']) ; OCE_mas_wfcorr  = OCE_flux_int ( OCE_wfcorr )
955if OCE_relax  :
956    # ssr and fwb are included in emp=>empmr but not in emp_oce (outputed by sea-ice)
957    OCE_vflx_fwb = rprec (d_OCE_his['vflx_fwb']) ; OCE_mas_vflx_fwb   = OCE_flux_int ( OCE_vflx_fwb )
958    OCE_vflx_ssr = rprec (d_OCE_his['vflx_ssr']) ; OCE_mas_vflx_ssr   = OCE_flux_int ( OCE_vflx_ssr )
959else : 
960    OCE_fwb = 0.0 ; OCE_mas_fwb = 0.0
961    OCE_ssr = 0.0 ; OCE_mas_ssr = 0.0
962if OCE_icb :
963    OCE_berg_icb    = rprec (d_OCE_his['berg_floating_melt']) ; OCE_mas_berg_icb = OCE_flux_int ( OCE_berg_icb    )
964    OCE_calving_icb = rprec (d_OCE_his['calving_icb']       ) ; OCE_mas_calv_icb = OCE_flux_int ( OCE_calving_icb )
965else :
966    OCE_berg_icb = 0. ; OCE_mas_berg_icb = 0.
967    OCE_calv_icb = 0. ; OCE_mas_calv_icb = 0.
968
969OCE_mas_emp = OCE_mas_emp_oce - OCE_mas_wfxice  - OCE_mas_wfxsnw  - ICE_mas_wfxpnd - ICE_mas_wfxsub_err
970OCE_mas_all = OCE_mas_emp_oce + OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf
971
972prtFlux ('OCE+ICE budget ', OCE_mas_all       , 'e', True)
973prtFlux ('  EMPMR        ', OCE_mas_empmr     , 'e', True)
974prtFlux ('  WFOB         ', OCE_mas_wfob      , 'e', True)
975prtFlux ('  EMP_OCE      ', OCE_mas_emp_oce   , 'e', True)
976prtFlux ('  EMP_ICE      ', OCE_mas_emp_ice   , 'e', True)
977prtFlux ('  EMP          ', OCE_mas_emp       , 'e', True)
978prtFlux ('  ICEBERG      ', OCE_mas_iceberg   , 'e',     )
979prtFlux ('  ICESHELF     ', OCE_mas_iceshelf  , 'e', True)
980prtFlux ('  CALVING      ', OCE_mas_calving   , 'e', True)
981prtFlux ('  FRIVER       ', OCE_mas_friver    , 'e',     ) 
982prtFlux ('  RUNOFFS      ', OCE_mas_runoffs   , 'e', True)
983prtFlux ('  WFXICE       ', OCE_mas_wfxice    , 'e', True)
984prtFlux ('  WFXSNW       ', OCE_mas_wfxsnw    , 'e', True)
985prtFlux ('  WFXSUB       ', OCE_mas_wfxsub    , 'e', True)
986prtFlux ('  WFXPND       ', ICE_mas_wfxpnd    , 'e', True)
987prtFlux ('  WFXSNW_SUB   ', ICE_mas_wfxsnw_sub, 'e', True)
988prtFlux ('  WFXSNW_PRE   ', ICE_mas_wfxsnw_pre, 'e', True)
989prtFlux ('  WFXSUB_ERR   ', ICE_mas_wfxsub_err, 'e', True)
990prtFlux ('  EVAP_OCE     ', OCE_mas_evap_oce  , 'e'      )
991prtFlux ('  EVAP_ICE     ', ICE_mas_evap_ice  , 'e', True)
992prtFlux ('  SNOW_OCE     ', OCE_mas_snow_oce  , 'e', True)
993prtFlux ('  SNOW_ICE     ', OCE_mas_snow_ice  , 'e', True)
994prtFlux ('  RAIN         ', OCE_mas_rain      , 'e'      )
995prtFlux ('  FWB          ', OCE_mas_fwb       , 'e', True)
996prtFlux ('  SSR          ', OCE_mas_ssr       , 'e', True)
997prtFlux ('  WFCORR       ', OCE_mas_wfcorr    , 'e', True)
998prtFlux ('  BERG_ICB     ', OCE_mas_berg_icb  , 'e', True)
999prtFlux ('  CALV_ICB     ', OCE_mas_calv_icb  , 'e', True)
1000
1001
1002echo (' ')
1003
1004prtFlux ( 'wbilo sea  ', ATM_flux_int (ATM_wbilo_sea), 'e',   )
1005if SRF : 
1006    prtFlux ( 'costalflow ', ONE_flux_int (RUN_coastalflow), 'e',   )
1007    prtFlux ( 'riverflow  ', RUN_flx_river  , 'e',   )
1008    prtFlux ( 'costalflow ', RUN_flx_coastal, 'e',   )
1009    prtFlux ( 'runoff     ', RUN_flx_river+RUN_flx_coastal, 'e',   )
1010
1011#ATM_to_OCE   =  ATM_flux_int (ATM_wbilo_sea) - RUN_flx_river - RUN_flx_coastal - ATM_flx_calving
1012ATM_to_OCE   =  ATM_flux_int (ATM_wbilo_sea) - ATM_flx_river - ATM_flx_coastal - ATM_flx_calving
1013#OCE_from_ATM = -OCE_mas_emp_oce - OCE_mas_emp_ice + OCE_mas_runoffs + OCE_mas_iceberg + OCE_mas_calving + OCE_mas_iceshelf
1014OCE_from_ATM = OCE_mas_all
1015
1016prtFlux ( 'ATM_to_OCE  ', ATM_to_OCE  , 'e', True )
1017prtFlux ( 'OCE_from_ATM', OCE_from_ATM, 'e', True )
1018
1019echo ( ' ' )
1020echo ( f'{Title = }' )
1021
1022echo ( 'SVN Information' )
1023for kk in SVN.keys():
1024    print ( SVN[kk] )
Note: See TracBrowser for help on using the repository browser.