source: TOOLS/WATER_BUDGET/CPL_waterbudget.py @ 6647

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

O.M. :

New version of WATER_BUDGET

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