source: TOOLS/WATER_BUDGET/CPL_waterbudget.py @ 6652

Last change on this file since 6652 was 6651, checked in by omamce, 7 months ago

O.M. : TOOLS/MOSAIX - Cosmetic changes
  • Property svn:keywords set to Date Revision HeadURL Author Id
File size: 46.9 KB
Line 
1#!/usr/bin/env python3
2###
3### Script to check water conservation in the IPSL coupled model
4###
5##  Warning, to install, configure, run, use any of included software or
6##  to read the associated documentation you'll need at least one (1)
7##  brain in a reasonably working order. Lack of this implement will
8##  void any warranties (either express or implied).  Authors assumes
9##  no responsability for errors, omissions, data loss, or any other
10##  consequences caused directly or indirectly by the usage of his
11##  software by incorrectly or partially configured personal
12##
13##
14# SVN information
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
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), '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'{readPrec        = }' )
287
288echo ( f'{JobName         = }' )
289echo ( f'{ConfigCard      = }' )
290echo ( f'{libIGCM         = }' )     
291echo ( f'{User            = }' )       
292echo ( f'{Group           = }' )       
293echo ( f'{Freq            = }' )       
294echo ( f'{YearBegin       = }' )     
295echo ( f'{YearEnd         = }' )     
296echo ( f'{DateBegin       = }' )
297echo ( f'{DateEnd         = }' )
298echo ( f'{PackFrequency   = }' ) 
299echo ( f'{ATM             = }' )       
300echo ( f'{Routing         = }' )       
301echo ( f'{ORCA            = }' )     
302echo ( f'{NEMO            = }' )     
303echo ( f'{Coupled         = }' )     
304echo ( f'{ATM_HIS         = }' )     
305echo ( f'{SRF_HIS         = }' )     
306echo ( f'{RUN_HIS         = }' )
307
308## Set libIGCM directories
309## -----------------------
310if wu.unDefined ('R_OUT'      ) : R_OUT       = os.path.join ( ARCHIVE   , 'IGCM_OUT' )
311if wu.unDefined ('R_BUF'      ) : R_BUF       = os.path.join ( SCRATCHDIR, 'IGCM_OUT' )
312if wu.unDefined ('L_EXP'      ) : L_EXP       = os.path.join (TagName, SpaceName, ExperimentName, JobName)
313if wu.unDefined ('R_SAVE'     ) : R_SAVE      = os.path.join ( R_OUT, L_EXP )
314if wu.unDefined ('R_BUFR'     ) : R_BUFR      = os.path.join ( R_BUF, L_EXP )
315if wu.unDefined ('POST_DIR'   ) : POST_DIR    = os.path.join ( R_BUFR, 'Out' )
316if wu.unDefined ('REBUILD_DIR') : REBUILD_DIR = os.path.join ( R_BUFR, 'REBUILD' )
317if wu.unDefined ('R_BUF_KSH'  ) : R_BUF_KSH   = os.path.join ( R_BUFR, 'Out' )
318if wu.unDefined ('R_FIGR'     ) : R_FIGR      = os.path.join ( STORAGE, 'IGCM_OUT', L_EXP )
319
320config['libIGCM'].update ( { 'R_OUT':R_OUT, 'R_BUF':R_BUF, 'L_EXP':L_EXP, 'R_BUFR':R_BUFR, 'POST_DIR':POST_DIR,
321                      'REBUILD_DIR':REBUILD_DIR, 'R_BUF_KSH':R_BUF_KSH, 'R_FIGR':R_FIGR, 'rebuild':rebuild } )
322
323## Set directory to extract files
324## ------------------------------
325if wu.unDefined ( 'FileDir' ) : FileDir = os.path.join ( TmpDir, f'WATER_{JobName}' )
326config['Files']['FileDir'] = FileDir
327
328if not os.path.isdir ( FileDir ) : os.makedirs ( FileDir )
329
330##- Set directories to rebuild ocean and ice restart files
331if wu.unDefined ( 'FileDirOCE' ) : FileDirOCE = os.path.join ( FileDir, 'OCE' )
332if wu.unDefined ( 'FileDirICE' ) : FileDirICE = os.path.join ( FileDir, 'ICE' )
333if not os.path.exists ( FileDirOCE ) : os.mkdir ( FileDirOCE )
334if not os.path.exists ( FileDirICE ) : os.mkdir ( FileDirICE )
335
336echo (' ')
337echo ( f'JobName     : {JobName}'    )
338echo ( f'Comment     : {Comment}'    )
339echo ( f'TmpDir      : {TmpDir}'     )
340echo ( f'FileDir     : {FileDir}'    )
341echo ( f'FileDirOCE  : {FileDirOCE}' )
342echo ( f'FileDirICE  : {FileDirICE}' )
343
344echo ( f'\nDealing with {L_EXP}'  )
345
346echo (' ')
347echo ( f'JobName   : {JobName}'   )
348echo ( f'Comment   : {Comment}'   )
349echo ( f'TmpDir    : {TmpDir}'    )
350
351echo ( f'\nDealing with {L_EXP}'  )
352
353## Creates model output directory names
354## ------------------------------------
355if Freq == "MO" : FreqDir =  os.path.join ('Output' , 'MO' )
356if Freq == "SE" : FreqDir =  os.path.join ('Analyse', 'SE' )
357if dir_ATM_his == None :
358    dir_ATM_his = os.path.join ( R_SAVE, "ATM", FreqDir )
359    config['Files']['dir_ATM_his'] = dir_ATM_his
360if dir_SRF_his == None : 
361    dir_SRF_his = os.path.join ( R_SAVE, "SRF", FreqDir )
362    config['Files']['dir_SRF_his'] = dir_SRF_his
363if dir_OCE_his == None :
364    dir_OCE_his = os.path.join ( R_SAVE, "OCE", FreqDir )
365    config['Files']['dir_OCE_his'] = dir_OCE_his
366if dir_ICE_his == None : 
367    dir_ICE_his = os.path.join ( R_SAVE, "ICE", FreqDir )
368    config['Files']['dir_ICE_his'] = dir_ICE_his
369
370echo ( f'The analysis relies on files from the following model output directories : ' )
371echo ( f'{dir_ATM_his}' )
372echo ( f'{dir_OCE_his}' )
373echo ( f'{dir_ICE_his}' )
374echo ( f'{dir_SRF_his}' )
375
376##-- Creates files names
377if wu.unDefined ( 'Period' ) :
378    if Freq == 'MO' : Period = f'{DateBegin}_{DateEnd}_1M'
379    if Freq == 'SE' : Period = f'SE_{DateBegin}_{DateEnd}_1M'
380    config['Files']['Period'] = Period
381
382config['Files']['DateBegin'] = DateBegin
383config['Files']['DateBegin'] = DateEnd
384
385echo ( f'Period   : {Period}' )
386
387if wu.unDefined ( 'FileCommon' ) :
388    FileCommon = f'{JobName}_{Period}'
389    config['Files']['FileCommon'] = FileCommon
390
391if wu.unDefined ( 'Title' ) :
392    Title = f'{JobName} : {Freq} : {DateBegin} - {DateEnd}'
393    config['Files']['Title'] = Title
394echo ('\nOpen history files' )
395if wu.unDefined ( 'file_ATM_his' ) :
396    if ATM_HIS == 'latlon' :
397        file_ATM_his = os.path.join ( dir_ATM_his, f'{FileCommon}_histmth.nc' )
398    if ATM_HIS == 'ico' :
399        file_ATM_his = os.path.join ( dir_ATM_his, f'{FileCommon}_histmth_ico.nc' )
400    config['Files']['file_ATM_his'] = file_ATM_his
401if wu.unDefined ( 'file_SRF_his' ) :
402    if ATM_HIS == 'latlon' :
403        file_SRF_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' )
404    if ATM_HIS == 'ico' :
405        file_SRF_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history_ico.nc' )
406    config['Files']['file_SRF_his'] = file_SRF_his
407
408if Routing == 'SIMPLE' :
409    if file_RUN_his == None :
410        if ATM_HIS == 'latlon' :
411            file_RUN_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history.nc' )
412        if ATM_HIS == 'ico' :
413            file_RUN_his = os.path.join ( dir_SRF_his, f'{FileCommon}_sechiba_history_ico.nc' )
414        config['Files']['file_RUN_his'] = file_RUN_his
415
416echo ( f'{file_ATM_his = }' )
417echo ( f'{file_SRF_his = }' )
418if Routing == 'SIMPLE' : echo ( f'{file_RUN_his = }' )
419       
420d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
421d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
422if Routing == 'SECHIBA' : d_RUN_his = d_SRF_his
423if Routing == 'SIMPLE'  : d_RUN_his = xr.open_dataset ( file_RUN_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
424   
425if wu.unDefined ('file_OCE_his' ) :
426    file_OCE_his = os.path.join ( dir_OCE_his, f'{FileCommon}_grid_T.nc' )
427    file_OCE_his = file_OCE_his
428if wu.unDefined ('file_OCE_sca' ) :   
429    file_OCE_sca = os.path.join ( dir_OCE_his, f'{FileCommon}_scalar.nc' )
430    config['Files']['file_OCE_sca'] = file_OCE_sca
431if wu.unDefined ('file_OCE_srf' ) :   
432    file_OCE_srf = os.path.join ( dir_OCE_his, f'{FileCommon}_sbc.nc' )
433    config['Files']['file_OCE_srf'] = file_OCE_srf
434if wu.unDefined ( 'file_ICE_hi' ) : 
435    file_ICE_his = os.path.join ( dir_ICE_his, f'{FileCommon}_icemod.nc' )
436    config['Files']['file_ICE_his'] = file_ICE_his
437
438d_OCE_his = xr.open_dataset ( file_OCE_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
439d_OCE_sca = xr.open_dataset ( file_OCE_sca, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
440#d_OCE_srf = xr.open_dataset ( file_OCE_srf, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
441d_ICE_his = xr.open_dataset ( file_ICE_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze()
442if NEMO == '3.6' :d_ICE_his = d_ICE_his.rename ( {'y_grid_T':'y', 'x_grid_T':'x'} )
443
444echo ( f'{file_OCE_his = }' )
445echo ( f'{file_ICE_his = }' )
446echo ( f'{file_OCE_sca = }' )
447echo ( f'{file_OCE_srf = }' )
448
449## Compute run length
450## ------------------
451dtime = ( d_ATM_his.time_counter_bounds.max() - d_ATM_his.time_counter_bounds.min() )
452echo ('\nRun length : {:8.2f} days'.format ( (dtime/np.timedelta64(1, "D")).values ) )
453dtime_sec = (dtime/np.timedelta64(1, "s")).values.item() # Convert in seconds
454
455## Compute length of each period
456dtime_per = (d_ATM_his.time_counter_bounds[:,-1] -  d_ATM_his.time_counter_bounds[:,0] )
457echo ('\nPeriods lengths (days) : {:} days'.format ( (dtime_per/np.timedelta64(1, "D")).values ) )
458dtime_per_sec = (dtime_per/np.timedelta64(1, "s")).values # In seconds
459dtime_per_sec = xr.DataArray (dtime_per_sec, dims=["time_counter", ], coords=[d_ATM_his.time_counter,] )
460dtime_per_sec.attrs['unit'] = 's'
461
462# Number of years
463NbYear = dtime_sec / YearLength
464
465## Write the full configuration
466config_out = open (FullIniFile, 'w')
467config.write (config_out )
468config_out.close ()
469
470# ATM grid with cell surfaces
471if LMDZ :
472    echo ('ATM grid with cell surfaces : LMDZ')
473    ATM_lat       = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1D='cell' )
474    ATM_lon       = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1D='cell' )
475    ATM_aire      = lmdz.geo2point ( rprec (d_ATM_his ['aire'] )    [0], cumulPoles=True, dim1D='cell' )
476    ATM_fter      = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1D='cell' )
477    ATM_foce      = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1D='cell' )
478    ATM_fsic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1D='cell' )
479    ATM_flic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1D='cell' )
480    SRF_lat       = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1D='cell' )
481    SRF_lon       = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1D='cell' )
482    SRF_aire      = lmdz.geo2point ( rprec (d_SRF_his ['Areas']) * rprec (d_SRF_his ['Contfrac']), dim1D='cell', cumulPoles=True )
483    SRF_areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas'])  ,  dim1D='cell', cumulPoles=True )
484    SRF_contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac']), dim1D='cell' )
485if ICO :
486    if ATM_HIS == 'latlon' :
487        echo ( 'ATM areas and fractions on latlon grid' )
488        if 'lat_dom_out' in d_ATM_his.variables :
489            ATM_lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat_dom_out'])+0*rprec (d_ATM_his ['lon_dom_out']), dim1D='cell' )
490            ATM_lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat_dom_out'])+  rprec (d_ATM_his ['lon_dom_out']), dim1D='cell' )
491        else :
492            ATM_lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1D='cell' )
493            ATM_lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1D='cell' )
494        ATM_aire = lmdz.geo2point ( rprec (d_ATM_his ['aire'][0]).squeeze(), cumulPoles=True, dim1D='cell' )
495        ATM_fter = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1D='cell' )
496        ATM_foce = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1D='cell' )
497        ATM_fsic = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1D='cell' )
498        ATM_flic = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1D='cell' )
499
500    if ATM_HIS == 'ico' :
501        echo ( 'ATM areas and fractions on ICO grid' )
502        ATM_aire =  rprec (d_ATM_his ['aire']     [0]).squeeze()
503        ATM_lat  =  rprec (d_ATM_his ['lat']         )
504        ATM_lon  =  rprec (d_ATM_his ['lon']         )
505        ATM_fter =  rprec (d_ATM_his ['fract_ter'][0])
506        ATM_foce =  rprec (d_ATM_his ['fract_oce'][0])
507        ATM_fsic =  rprec (d_ATM_his ['fract_sic'][0])
508        ATM_flic =  rprec (d_ATM_his ['fract_lic'][0])
509
510    if SRF_HIS == 'latlon' :
511        echo ( 'SRF areas and fractions on latlon grid' )
512        if 'lat_domain_landpoints_out' in d_SRF_his  : 
513            SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat_domain_landpoints_out'])+0*rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1D='cell' )
514            SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_domain_landpoints_out'])+  rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1D='cell' )
515        else : 
516            if 'lat_domain_landpoints_out' in d_SRF_his :
517                SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat_dom_out'])+0*rprec (d_SRF_his ['lon_dom_out']), dim1D='cell' )
518                SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_dom_out'])+  rprec (d_SRF_his ['lon_dom_out']), dim1D='cell' )
519            else :
520                SRF_lat  = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1D='cell' )
521                SRF_lon  = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1D='cell' )
522               
523        SRF_areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas']   )   , dim1D='cell', cumulPoles=True )
524        SRF_areafrac  = lmdz.geo2point ( rprec (d_SRF_his ['AreaFrac'])   , dim1D='cell', cumulPoles=True )
525        SRF_contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac'])   , dim1D='cell', cumulPoles=True )
526        SRF_aire = SRF_areafrac
527 
528    if SRF_HIS == 'ico' :
529        echo ( 'SRF areas and fractions on latlon grid' )
530        SRF_lat       =  rprec (d_SRF_his ['lat']     )
531        SRF_lon       =  rprec (d_SRF_his ['lon']     )
532        SRF_areas     =  rprec (d_SRF_his ['Areas']   ) 
533        SRF_contfrac  =  rprec (d_SRF_his ['Contfrac'])
534        SRF_aire      =  SRF_areas * SRF_contfrac
535ATM_fsea      = ATM_foce + ATM_fsic
536ATM_flnd      = ATM_fter + ATM_flic
537ATM_aire_fter = ATM_aire * ATM_fter
538ATM_aire_flic = ATM_aire * ATM_flic
539ATM_aire_fsic = ATM_aire * ATM_fsic
540ATM_aire_foce = ATM_aire * ATM_foce
541ATM_aire_flnd = ATM_aire * ATM_flnd
542ATM_aire_fsea = ATM_aire * ATM_fsea
543
544#SRF_aire = SRF_aire.where ( SRF_aire < 1E15, 0.)
545
546# if ICO :
547#     if wu.unDefined ('file_DYN_aire') : file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ATM+'_grid.nc' )
548#     config['Files']['file_DYN_aire'] = file_DYN_aire
549
550# if ICO :
551#     # Area on icosahedron grid
552#     d_DYN_aire = xr.open_dataset ( file_DYN_aire, decode_times=False ).squeeze()
553           
554#     DYN_lat = d_DYN_aire['lat']
555#     DYN_lon = d_DYN_aire['lon']
556
557#     DYN_aire = d_DYN_aire['aire']
558#     DYN_fsea = d_DYN_aire['fract_oce'] + d_DYN_aire['fract_sic']
559#     DYN_flnd = 1.0 - DYN_fsea
560#     DYN_fter = d_ATM_beg['FTER']
561#     DYN_flic = d_ATM_beg['FLIC']
562#     DYN_aire_fter = DYN_aire * DYN_fter
563   
564# if LMDZ :
565#     # Area on lon/lat grid
566#     DYN_aire = ATM_aire
567#     DYN_fsea = ATM_fsea
568#     DYN_flnd = ATM_flnd
569#     DYN_fter = rprec (d_ATM_beg['FTER'])
570#     DYN_flic = rprec (d_ATM_beg['FLIC'])
571#     DYN_aire_fter = DYN_aire * DYN_fter
572
573# Functions computing integrals and sum
574# def ATM_stock_int (stock) :
575#     '''Integrate (* surface) stock on atmosphere grid'''
576#     ATM_stock_int  = wu.Psum ( (stock * DYN_aire).to_masked_array().ravel() )
577#     return ATM_stock_int
578
579def ATM_flux_int (flux) :
580    '''Integrate (* time * surface) flux on atmosphere grid'''
581    ATM_flux_int  = wu.Psum ( (flux * dtime_per_sec * ATM_aire).to_masked_array().ravel() )
582    return ATM_flux_int
583
584# def SRF_stock_int (stock) :
585#     '''Integrate (* surface) stock on land grid'''
586#     ATM_stock_int  = wu.Ksum (  ( (stock * DYN_aire_fter).to_masked_array().ravel()) )
587#     return ATM_stock_int
588
589def SRF_flux_int (flux) :
590    '''Integrate (* time * surface) flux on land grid'''
591    SRF_flux_int  = wu.Psum (  (flux * dtime_per_sec * SRF_aire).to_masked_array().ravel() )
592    return SRF_flux_int
593
594def LIC_flux_int (flux) :
595    '''Integrate (* time * surface) flux on land ice grid'''
596    LIC_flux_int  = wu.Psum ( (flux * dtime_per_sec * ATM_aire_flic).to_masked_array().ravel() )
597    return LIC_flux_int
598
599# def OCE_stock_int (stock) :
600#     '''Integrate stock on ocean grid'''
601#     OCE_stock_int = np.sum (  np.sort ( (stock * OCE_aire ).to_masked_array().ravel()) )
602#     return OCE_stock_int
603
604def ONE_stock_int (stock) :
605    '''Sum stock'''
606    ONE_stock_int =  wu.Psum ( stock.to_masked_array().ravel() )
607    return ONE_stock_int
608
609def OCE_flux_int (flux) :
610    '''Integrate flux on oce grid'''
611    OCE_flux_int = wu.Psum ( (flux * OCE_aire * dtime_per_sec).to_masked_array().ravel() )
612    return OCE_flux_int
613
614def ONE_flux_int (flux) :
615    '''Integrate flux on oce grid'''
616    OCE_flux_int = wu.Psum ( (flux * dtime_per_sec).to_masked_array().ravel() )
617    return OCE_flux_int
618
619# Get mask and surfaces
620sos = d_OCE_his ['sos'][0].squeeze()
621OCE_msk = nemo.lbc_mask ( xr.where ( sos>0, 1., 0.0 ), cd_type = 'T' )
622so = sos = d_OCE_his ['sos'][0].squeeze()
623OCE_msk3 = nemo.lbc_mask ( xr.where ( so>0., 1., 0. ), cd_type = 'T', sval = 0. )
624
625# lbc_mask removes the duplicate points (periodicity and north fold)
626OCE_aire = nemo.lbc_mask ( d_OCE_his ['area'] * OCE_msk, cd_type = 'T', sval = 0.0 )
627ICE_aire = OCE_aire
628
629ATM_aire_tot = ONE_stock_int (ATM_aire)
630SRF_aire_tot = ONE_stock_int (SRF_aire)
631OCE_aire_tot = ONE_stock_int (OCE_aire)
632ICE_aire_tot = ONE_stock_int (ICE_aire)
633
634ATM_aire_sea     = ATM_aire * ATM_fsea
635ATM_aire_sea_tot = ONE_stock_int ( ATM_aire_sea )
636
637echo ( '\n====================================================================================' )
638echo ( f'-- ATM Fluxes  -- {Title} ' )
639
640if ATM_HIS == 'latlon' :
641    echo ( ' latlon case' )
642    ATM_wbilo_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_oce']), dim1D='cell' )
643    ATM_wbilo_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_sic']), dim1D='cell' )
644    ATM_wbilo_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_ter']), dim1D='cell' )
645    ATM_wbilo_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wbilo_lic']), dim1D='cell' )
646    ATM_runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1D='cell' )
647    ATM_fqcalving   = lmdz.geo2point ( rprec (d_ATM_his ['fqcalving']), dim1D='cell' )
648    ATM_fqfonte     = lmdz.geo2point ( rprec (d_ATM_his ['fqfonte']  ), dim1D='cell' )
649    ATM_precip      = lmdz.geo2point ( rprec (d_ATM_his ['precip']   ), dim1D='cell' )
650    ATM_snowf       = lmdz.geo2point ( rprec (d_ATM_his ['snow']     ), dim1D='cell' )
651    ATM_evap        = lmdz.geo2point ( rprec (d_ATM_his ['evap']     ), dim1D='cell' )
652    ATM_wevap_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_ter']), dim1D='cell' )
653    ATM_wevap_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_oce']), dim1D='cell' )
654    ATM_wevap_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_lic']), dim1D='cell' )
655    ATM_wevap_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wevap_sic']), dim1D='cell' )
656    ATM_wrain_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_ter']), dim1D='cell' )
657    ATM_wrain_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_oce']), dim1D='cell' )
658    ATM_wrain_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_lic']), dim1D='cell' )
659    ATM_wrain_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wrain_sic']), dim1D='cell' )
660    ATM_wsnow_ter   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_ter']), dim1D='cell' )
661    ATM_wsnow_oce   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_oce']), dim1D='cell' )
662    ATM_wsnow_lic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_lic']), dim1D='cell' )
663    ATM_wsnow_sic   = lmdz.geo2point ( rprec (d_ATM_his ['wsnow_sic']), dim1D='cell' )
664    ATM_runofflic   = lmdz.geo2point ( rprec (d_ATM_his ['runofflic']), dim1D='cell' )
665    echo ( f'End of LATLON case')
666   
667if ATM_HIS == 'ico' :
668    echo (' ico case')
669    ATM_wbilo_oce   = rprec (d_ATM_his ['wbilo_oce'])
670    ATM_wbilo_sic   = rprec (d_ATM_his ['wbilo_sic'])
671    ATM_wbilo_ter   = rprec (d_ATM_his ['wbilo_ter'])
672    ATM_wbilo_lic   = rprec (d_ATM_his ['wbilo_lic'])
673    ATM_runofflic   = rprec (d_ATM_his ['runofflic'])
674    ATM_fqcalving   = rprec (d_ATM_his ['fqcalving'])
675    ATM_fqfonte     = rprec (d_ATM_his ['fqfonte']  )
676    ATM_precip      = rprec (d_ATM_his ['precip']   )
677    ATM_snowf       = rprec (d_ATM_his ['snow']     )
678    ATM_evap        = rprec (d_ATM_his ['evap']     )
679    ATM_wevap_ter   = rprec (d_ATM_his ['wevap_ter'])
680    ATM_wevap_oce   = rprec (d_ATM_his ['wevap_oce'])
681    ATM_wevap_lic   = rprec (d_ATM_his ['wevap_lic'])
682    ATM_wevap_sic   = rprec (d_ATM_his ['wevap_sic'])
683    ATM_runofflic   = rprec (d_ATM_his ['runofflic'])
684    ATM_wevap_ter   = rprec (d_ATM_his ['wevap_ter'])
685    ATM_wevap_oce   = rprec (d_ATM_his ['wevap_oce'])
686    ATM_wevap_lic   = rprec (d_ATM_his ['wevap_lic'])
687    ATM_wevap_sic   = rprec (d_ATM_his ['wevap_sic'])
688    ATM_wrain_ter   = rprec (d_ATM_his ['wrain_ter'])
689    ATM_wrain_oce   = rprec (d_ATM_his ['wrain_oce'])
690    ATM_wrain_lic   = rprec (d_ATM_his ['wrain_lic'])
691    ATM_wrain_sic   = rprec (d_ATM_his ['wrain_sic'])
692    ATM_wsnow_ter   = rprec (d_ATM_his ['wsnow_ter'])
693    ATM_wsnow_oce   = rprec (d_ATM_his ['wsnow_oce'])
694    ATM_wsnow_lic   = rprec (d_ATM_his ['wsnow_lic'])
695    ATM_wsnow_sic   = rprec (d_ATM_his ['wsnow_sic'])
696    echo ( f'End of ico case ')
697
698
699   
700echo ( 'ATM wprecip_oce' )
701ATM_wprecip_oce = ATM_wrain_oce + ATM_wsnow_oce
702ATM_wprecip_ter = ATM_wrain_ter + ATM_wsnow_ter
703ATM_wprecip_sic = ATM_wrain_sic + ATM_wsnow_sic
704ATM_wprecip_lic = ATM_wrain_lic + ATM_wsnow_lic
705
706ATM_wbilo       = ATM_wbilo_oce   + ATM_wbilo_sic   + ATM_wbilo_ter   + ATM_wbilo_lic
707ATM_wevap       = ATM_wevap_oce   + ATM_wevap_sic   + ATM_wevap_ter   + ATM_wevap_lic
708ATM_wprecip     = ATM_wprecip_oce + ATM_wprecip_sic + ATM_wprecip_ter + ATM_wprecip_lic
709ATM_wsnow       = ATM_wsnow_oce   + ATM_wsnow_sic   + ATM_wsnow_ter   + ATM_wsnow_lic
710ATM_wrain       = ATM_wrain_oce   + ATM_wrain_sic   + ATM_wrain_ter   + ATM_wrain_lic
711ATM_wemp        = ATM_wevap - ATM_wprecip
712ATM_emp         = ATM_evap  - ATM_precip
713
714ATM_wprecip_sea = ATM_wprecip_oce + ATM_wprecip_sic
715ATM_wsnow_sea   = ATM_wsnow_oce   + ATM_wsnow_sic
716ATM_wrain_sea   = ATM_wrain_oce   + ATM_wrain_sic
717ATM_wbilo_sea   = ATM_wbilo_oce   + ATM_wbilo_sic
718ATM_wevap_sea   = ATM_wevap_sic   + ATM_wevap_oce
719
720ATM_wemp_ter    = ATM_wevap_ter - ATM_wprecip_ter
721ATM_wemp_oce    = ATM_wevap_oce - ATM_wprecip_oce
722ATM_wemp_sic    = ATM_wevap_sic - ATM_wprecip_sic
723ATM_wemp_lic    = ATM_wevap_lic - ATM_wprecip_lic
724ATM_wemp_sea    = ATM_wevap_sic - ATM_wprecip_oce
725
726if RUN_HIS == 'latlon' :
727    echo ( f'RUN costalflow Grille LATLON' )
728    if TestInterp :
729        echo ( f'RUN runoff TestInterp' )
730        RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff_contfrac_interp']  )   , dim1D='cell' )
731        RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage_contfrac_interp'])   , dim1D='cell' )
732    else : 
733        echo ( f'RUN runoff' )
734        RUN_runoff      = lmdz.geo2point ( rprec (d_RUN_his ['runoff']         ), dim1D='cell' )
735        RUN_drainage    = lmdz.geo2point ( rprec (d_RUN_his ['drainage']       ), dim1D='cell' )
736       
737    RUN_coastalflow     = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow']    ), dim1D='cell' ) 
738    RUN_riverflow       = lmdz.geo2point ( rprec (d_RUN_his ['riverflow']      ), dim1D='cell' )
739    RUN_riversret       = lmdz.geo2point ( rprec (d_RUN_his ['riversret']      ), dim1D='cell' )
740    RUN_coastalflow_cpl = lmdz.geo2point ( rprec (d_RUN_his ['coastalflow_cpl']), dim1D='cell' ) 
741    RUN_riverflow_cpl   = lmdz.geo2point ( rprec (d_RUN_his ['riverflow_cpl']  ), dim1D='cell' )
742
743if RUN_HIS == 'ico' :
744    echo ( f'RUN costalflow Grille ICO' )
745    RUN_coastalflow =  rprec (d_RUN_his ['coastalflow'])
746    RUN_riverflow   =  rprec (d_RUN_his ['riverflow']  )
747    RUN_runoff      =  rprec (d_RUN_his ['runoff']     )
748    RUN_drainage    =  rprec (d_RUN_his ['drainage']   )
749    RUN_riversret   =  rprec (d_RUN_his ['riversret']  )
750   
751    RUN_coastalflow_cpl = rprec (d_RUN_his ['coastalflow_cpl'])
752    RUN_riverflow_cpl   = rprec (d_RUN_his ['riverflow_cpl']  )
753
754## Correcting units of SECHIBA variables
755def mmd2SI ( Var ) :
756    '''Change unit from mm/d or m^3/s to kg/s if needed'''
757    if 'units' in VarT.attrs : 
758        if VarT.attrs['units'] in ['m^3/s', 'm3/s', 'm3.s-1',] :
759            VarT.values = VarT.values  * ATM_rho                 ;  VarT.attrs['units'] = 'kg/s'
760        if VarT.attrs['units'] == 'mm/d' :
761            VarT.values = VarT.values  * ATM_rho * (1e-3/86400.) ;  VarT.attrs['units'] = 'kg/s'
762        if VarT.attrs['units'] in ['m^3', 'm3', ] :
763            VarT.values = VarT.values  * ATM_rho                 ;  VarT.attrs['units'] = 'kg'
764           
765for var in [ 'runoff', 'drainage', 'riversret', 'coastalflow', 'riverflow', 'coastalflow_cpl', 'riverflow_cpl' ] :
766    VarT = locals()['RUN_' + var]
767    mmd2SI (VarT)
768
769#for var in ['evap', 'snowf', 'subli', 'transpir', 'rain', 'emp' ] :
770#    VarT = locals()['SRF_' + var]
771#    mmd2SI (VarT)
772echo ( f'RUN input' )
773RUN_input  = RUN_runoff      + RUN_drainage
774RUN_output = RUN_coastalflow + RUN_riverflow
775
776echo ( f'ATM flw_wbilo' )
777ATM_flx_wbilo       = ATM_flux_int ( ATM_wbilo      )
778ATM_flx_wevap       = ATM_flux_int ( ATM_wevap      )
779ATM_flx_wprecip     = ATM_flux_int ( ATM_wprecip    )
780ATM_flx_wsnow       = ATM_flux_int ( ATM_wsnow      )
781ATM_flx_wrain       = ATM_flux_int ( ATM_wrain      )
782ATM_flx_wemp        = ATM_flux_int ( ATM_wemp       )
783
784ATM_flx_wbilo_lic   = ATM_flux_int ( ATM_wbilo_lic  )
785ATM_flx_wbilo_oce   = ATM_flux_int ( ATM_wbilo_oce  )
786ATM_flx_wbilo_sea   = ATM_flux_int ( ATM_wbilo_sea  )
787ATM_flx_wbilo_sic   = ATM_flux_int ( ATM_wbilo_sic  )
788ATM_flx_wbilo_ter   = ATM_flux_int ( ATM_wbilo_ter  )
789ATM_flx_calving     = ATM_flux_int ( ATM_fqcalving  )
790ATM_flx_fqfonte     = ATM_flux_int ( ATM_fqfonte    )
791
792LIC_flx_calving     = LIC_flux_int ( ATM_fqcalving  )
793LIC_flx_fqfonte     = LIC_flux_int ( ATM_fqfonte    )
794
795ATM_flx_precip      = ATM_flux_int ( ATM_precip     )
796ATM_flx_snowf       = ATM_flux_int ( ATM_snowf      )
797ATM_flx_evap        = ATM_flux_int ( ATM_evap       )
798ATM_flx_runlic      = ATM_flux_int ( ATM_runofflic  )
799
800LIC_flx_precip      = LIC_flux_int ( ATM_precip     )
801LIC_flx_snowf       = LIC_flux_int ( ATM_snowf      )
802LIC_flx_evap        = LIC_flux_int ( ATM_evap       )
803LIC_flx_runlic      = LIC_flux_int ( ATM_runofflic  )
804
805ATM_flx_wrain_ter   = ATM_flux_int ( ATM_wrain_ter )
806ATM_flx_wrain_oce   = ATM_flux_int ( ATM_wrain_oce )
807ATM_flx_wrain_lic   = ATM_flux_int ( ATM_wrain_lic )
808ATM_flx_wrain_sic   = ATM_flux_int ( ATM_wrain_sic )
809ATM_flx_wrain_sea   = ATM_flux_int ( ATM_wrain_sea )
810
811ATM_flx_wsnow_ter   = ATM_flux_int ( ATM_wsnow_ter )
812ATM_flx_wsnow_oce   = ATM_flux_int ( ATM_wsnow_oce )
813ATM_flx_wsnow_lic   = ATM_flux_int ( ATM_wsnow_lic )
814ATM_flx_wsnow_sic   = ATM_flux_int ( ATM_wsnow_sic )
815ATM_flx_wsnow_sea   = ATM_flux_int ( ATM_wsnow_sea )
816
817ATM_flx_wevap_ter   = ATM_flux_int ( ATM_wevap_ter )
818ATM_flx_wevap_oce   = ATM_flux_int ( ATM_wevap_oce )
819ATM_flx_wevap_lic   = ATM_flux_int ( ATM_wevap_lic )
820ATM_flx_wevap_sic   = ATM_flux_int ( ATM_wevap_sic )
821ATM_flx_wevap_sea   = ATM_flux_int ( ATM_wevap_sea )
822ATM_flx_wprecip_lic = ATM_flux_int ( ATM_wprecip_lic )
823ATM_flx_wprecip_oce = ATM_flux_int ( ATM_wprecip_oce )
824ATM_flx_wprecip_sic = ATM_flux_int ( ATM_wprecip_sic )
825ATM_flx_wprecip_ter = ATM_flux_int ( ATM_wprecip_ter )
826ATM_flx_wprecip_sea = ATM_flux_int ( ATM_wprecip_sea )
827ATM_flx_wemp_lic    = ATM_flux_int ( ATM_wemp_lic )
828ATM_flx_wemp_oce    = ATM_flux_int ( ATM_wemp_oce )
829ATM_flx_wemp_sic    = ATM_flux_int ( ATM_wemp_sic )
830ATM_flx_wemp_ter    = ATM_flux_int ( ATM_wemp_ter )
831ATM_flx_wemp_sea    = ATM_flux_int ( ATM_wemp_sea )
832
833ATM_flx_emp          = ATM_flux_int ( ATM_emp )
834
835RUN_flx_coastal     = ONE_flux_int ( RUN_coastalflow)
836RUN_flx_river       = ONE_flux_int ( RUN_riverflow  )
837RUN_flx_coastal_cpl = ONE_flux_int ( RUN_coastalflow_cpl)
838RUN_flx_river_cpl   = ONE_flux_int ( RUN_riverflow_cpl  )
839RUN_flx_drainage    = SRF_flux_int ( RUN_drainage   )
840RUN_flx_riversret   = SRF_flux_int ( RUN_riversret  )
841RUN_flx_runoff      = SRF_flux_int ( RUN_runoff     )
842RUN_flx_input       = SRF_flux_int ( RUN_input      )
843RUN_flx_output      = ONE_flux_int ( RUN_output     )
844
845RUN_flx_bil    = ONE_flux_int ( RUN_input       - RUN_output)
846RUN_flx_rivcoa = ONE_flux_int ( RUN_coastalflow + RUN_riverflow)
847
848prtFlux ('wbilo_oce            ', ATM_flx_wbilo_oce     , 'f' )         
849prtFlux ('wbilo_sic            ', ATM_flx_wbilo_sic     , 'f' )         
850prtFlux ('wbilo_sic+oce        ', ATM_flx_wbilo_sea     , 'f' )         
851prtFlux ('wbilo_ter            ', ATM_flx_wbilo_ter     , 'f' )         
852prtFlux ('wbilo_lic            ', ATM_flx_wbilo_lic     , 'f' )         
853prtFlux ('Sum wbilo_*          ', ATM_flx_wbilo         , 'f', True) 
854prtFlux ('E-P                  ', ATM_flx_emp           , 'f', True) 
855prtFlux ('calving              ', ATM_flx_calving       , 'f' ) 
856prtFlux ('fqfonte              ', ATM_flx_fqfonte       , 'f' )       
857prtFlux ('precip               ', ATM_flx_precip        , 'f' )       
858prtFlux ('snowf                ', ATM_flx_snowf         , 'f' )       
859prtFlux ('evap                 ', ATM_flx_evap          , 'f' )
860prtFlux ('runoff lic           ', ATM_flx_runlic        , 'f' )
861
862prtFlux ('ATM_flx_wevap*       ', ATM_flx_wevap         , 'f' )
863prtFlux ('ATM_flx_wrain*       ', ATM_flx_wrain         , 'f' )
864prtFlux ('ATM_flx_wsnow*       ', ATM_flx_wsnow         , 'f' )
865prtFlux ('ATM_flx_wprecip*     ', ATM_flx_wprecip       , 'f' )
866prtFlux ('ATM_flx_wemp*        ', ATM_flx_wemp          , 'f', True )
867
868prtFlux ('ERROR evap           ', ATM_flx_wevap   - ATM_flx_evap  , 'e', True )
869prtFlux ('ERROR precip         ', ATM_flx_wprecip - ATM_flx_precip, 'e', True )
870prtFlux ('ERROR snow           ', ATM_flx_wsnow   - ATM_flx_snowf , 'e', True )
871prtFlux ('ERROR emp            ', ATM_flx_wemp    - ATM_flx_emp   , 'e', True )
872
873
874echo ( '\n====================================================================================' )
875echo ( f'-- RUNOFF Fluxes  -- {Title} ' )
876prtFlux ('coastalflow   ', RUN_flx_coastal    , 'f' ) 
877prtFlux ('riverflow     ', RUN_flx_river      , 'f' )       
878prtFlux ('coastal_cpl   ', RUN_flx_coastal_cpl, 'f' ) 
879prtFlux ('riverf_cpl    ', RUN_flx_river_cpl  , 'f' )   
880prtFlux ('river+coastal ', RUN_flx_rivcoa     , 'f' )   
881prtFlux ('drainage      ', RUN_flx_drainage   , 'f' )   
882prtFlux ('riversret     ', RUN_flx_riversret  , 'f' )   
883prtFlux ('runoff        ', RUN_flx_runoff     , 'f' )   
884prtFlux ('river in      ', RUN_flx_input      , 'f' )   
885prtFlux ('river out     ', RUN_flx_output     , 'f' )   
886prtFlux ('river bil     ', RUN_flx_bil        , 'f' ) 
887   
888echo ( '\n====================================================================================' )
889echo ( f'-- OCE Fluxes  -- {Title} ' )
890
891# Read variable and computes integral over space and time
892OCE_empmr      = rprec (d_OCE_his['wfo']     )    ; OCE_mas_empmr    = OCE_flux_int ( OCE_empmr    )
893OCE_wfob       = rprec (d_OCE_his['wfob']    )    ; OCE_mas_wfob     = OCE_flux_int ( OCE_wfob     )
894OCE_emp_oce    = rprec (d_OCE_his['emp_oce'] )    ; OCE_mas_emp_oce  = OCE_flux_int ( OCE_emp_oce  )
895OCE_emp_ice    = rprec (d_OCE_his['emp_ice'] )    ; OCE_mas_emp_ice  = OCE_flux_int ( OCE_emp_ice  )
896OCE_iceshelf   = rprec (d_OCE_his['iceshelf'])    ; OCE_mas_iceshelf = OCE_flux_int ( OCE_iceshelf )
897OCE_calving    = rprec (d_OCE_his['calving'] )    ; OCE_mas_calving  = OCE_flux_int ( OCE_calving  )
898OCE_iceberg    = rprec (d_OCE_his['iceberg'] )    ; OCE_mas_iceberg  = OCE_flux_int ( OCE_iceberg  )
899OCE_friver     = rprec (d_OCE_his['friver']  )    ; OCE_mas_friver   = OCE_flux_int ( OCE_friver   )
900OCE_runoffs    = rprec (d_OCE_his['runoffs'] )    ; OCE_mas_runoffs  = OCE_flux_int ( OCE_runoffs  )
901if NEMO == 4.0 or NEMO == 4.2 :
902    OCE_wfxice     = rprec (d_OCE_his['vfxice']) ; OCE_mas_wfxice   = OCE_flux_int ( OCE_wfxice )
903    OCE_wfxsnw     = rprec (d_OCE_his['vfxsnw']) ; OCE_mas_wfxsnw   = OCE_flux_int ( OCE_wfxsnw )
904    OCE_wfxsub     = rprec (d_OCE_his['vfxsub']) ; OCE_mas_wfxsub   = OCE_flux_int ( OCE_wfxsub )
905if NEMO == 3.6 :
906    OCE_wfxice     = rprec (d_OCE_his['vfxice'])/86400.*ICE_rho_ice ; OCE_mas_wfxice   = OCE_flux_int ( OCE_wfxice )
907    OCE_wfxsnw     = rprec (d_OCE_his['vfxsnw'])/86400.*ICE_rho_sno ; OCE_mas_wfxsnw   = OCE_flux_int ( OCE_wfxsnw )
908    OCE_wfxsub     = rprec (d_OCE_his['vfxsub'])/86400.*ICE_rho_sno ; OCE_mas_wfxsub   = OCE_flux_int ( OCE_wfxsub )
909# Additional checks
910OCE_evap_oce   = rprec (d_OCE_his['evap_ao_cea']) ; OCE_mas_evap_oce   = OCE_flux_int ( OCE_evap_oce )
911ICE_evap_ice   = rprec (d_OCE_his['subl_ai_cea']) ; ICE_mas_evap_ice   = OCE_flux_int ( ICE_evap_ice )
912OCE_snow_oce   = rprec (d_OCE_his['snow_ao_cea']) ; OCE_mas_snow_oce   = OCE_flux_int ( OCE_snow_oce )
913OCE_snow_ice   = rprec (d_OCE_his['snow_ai_cea']) ; OCE_mas_snow_ice   = OCE_flux_int ( OCE_snow_ice )
914OCE_rain       = rprec (d_OCE_his['rain']       ) ; OCE_mas_rain       = OCE_flux_int ( OCE_rain     )
915ICE_wfxsub_err = rprec (d_ICE_his['vfxsub_err'] ) ; ICE_mas_wfxsub_err = OCE_flux_int ( ICE_wfxsub_err )
916if NEMO == 4.0 or NEMO == 4.2 :
917    ICE_wfxpnd     = rprec (d_ICE_his['vfxpnd']    ) ; ICE_mas_wfxpnd     = OCE_flux_int ( ICE_wfxpnd     )
918    ICE_wfxsnw_sub = rprec (d_ICE_his['vfxsnw_sub']) ; ICE_mas_wfxsnw_sub = OCE_flux_int ( ICE_wfxsnw_sub )
919    ICE_wfxsnw_pre = rprec (d_ICE_his['vfxsnw_pre']) ; ICE_mas_wfxsnw_pre = OCE_flux_int ( ICE_wfxsnw_pre )
920if NEMO == 3.6 :
921    ICE_wfxpnd = 0.0 ; ICE_mas_wfxpnd = 0.0
922    ICE_wfxsnw_sub = rprec (d_ICE_his['vfxsub'])/86400.*ICE_rho_sno  ; ICE_mas_wfxsnw_sub = OCE_flux_int ( ICE_wfxsnw_sub )
923    ICE_wfxsnw_pre = rprec (d_ICE_his['vfxspr'])/86400.*ICE_rho_sno  ; ICE_mas_wfxsnw_pre = OCE_flux_int ( ICE_wfxsnw_pre     )
924
925OCE_wfcorr    = rprec (d_OCE_his['wfcorr']) ; OCE_mas_wfcorr  = OCE_flux_int ( OCE_wfcorr )
926if OCE_relax  :
927    # ssr and fwb are included in emp=>empmr but not in emp_oce (outputed by sea-ice)
928    OCE_vflx_fwb = rprec (d_OCE_his['vflx_fwb']) ; OCE_mas_vflx_fwb   = OCE_flux_int ( OCE_vflx_fwb )
929    OCE_vflx_ssr = rprec (d_OCE_his['vflx_ssr']) ; OCE_mas_vflx_ssr   = OCE_flux_int ( OCE_vflx_ssr )
930else : 
931    OCE_fwb = 0.0 ; OCE_mas_fwb = 0.0
932    OCE_ssr = 0.0 ; OCE_mas_ssr = 0.0
933if OCE_icb :
934    OCE_berg_icb    = rprec (d_OCE_his['berg_floating_melt']) ; OCE_mas_berg_icb = OCE_flux_int ( OCE_berg_icb    )
935    OCE_calving_icb = rprec (d_OCE_his['calving_icb']       ) ; OCE_mas_calv_icb = OCE_flux_int ( OCE_calving_icb )
936else :
937    OCE_berg_icb = 0. ; OCE_mas_berg_icb = 0.
938    OCE_calv_icb = 0. ; OCE_mas_calv_icb = 0.
939
940OCE_mas_emp = OCE_mas_emp_oce - OCE_mas_wfxice  - OCE_mas_wfxsnw  - ICE_mas_wfxpnd - ICE_mas_wfxsub_err
941OCE_mas_all = OCE_mas_emp_oce + OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf
942
943prtFlux ('OCE+ICE budget ', OCE_mas_all       , 'e', True)
944prtFlux ('  EMPMR        ', OCE_mas_empmr     , 'e', True)
945prtFlux ('  WFOB         ', OCE_mas_wfob      , 'e', True)
946prtFlux ('  EMP_OCE      ', OCE_mas_emp_oce   , 'e', True)
947prtFlux ('  EMP_ICE      ', OCE_mas_emp_ice   , 'e', True)
948prtFlux ('  EMP          ', OCE_mas_emp       , 'e', True)
949prtFlux ('  ICEBERG      ', OCE_mas_iceberg   , 'e',     )
950prtFlux ('  ICESHELF     ', OCE_mas_iceshelf  , 'e', True)
951prtFlux ('  CALVING      ', OCE_mas_calving   , 'e', True)
952prtFlux ('  FRIVER       ', OCE_mas_friver    , 'e',     ) 
953prtFlux ('  RUNOFFS      ', OCE_mas_runoffs   , 'e', True)
954prtFlux ('  WFXICE       ', OCE_mas_wfxice    , 'e', True)
955prtFlux ('  WFXSNW       ', OCE_mas_wfxsnw    , 'e', True)
956prtFlux ('  WFXSUB       ', OCE_mas_wfxsub    , 'e', True)
957prtFlux ('  WFXPND       ', ICE_mas_wfxpnd    , 'e', True)
958prtFlux ('  WFXSNW_SUB   ', ICE_mas_wfxsnw_sub, 'e', True)
959prtFlux ('  WFXSNW_PRE   ', ICE_mas_wfxsnw_pre, 'e', True)
960prtFlux ('  WFXSUB_ERR   ', ICE_mas_wfxsub_err, 'e', True)
961prtFlux ('  EVAP_OCE     ', OCE_mas_evap_oce  , 'e'      )
962prtFlux ('  EVAP_ICE     ', ICE_mas_evap_ice  , 'e', True)
963prtFlux ('  SNOW_OCE     ', OCE_mas_snow_oce  , 'e', True)
964prtFlux ('  SNOW_ICE     ', OCE_mas_snow_ice  , 'e', True)
965prtFlux ('  RAIN         ', OCE_mas_rain      , 'e'      )
966prtFlux ('  FWB          ', OCE_mas_fwb       , 'e', True)
967prtFlux ('  SSR          ', OCE_mas_ssr       , 'e', True)
968prtFlux ('  WFCORR       ', OCE_mas_wfcorr    , 'e', True)
969prtFlux ('  BERG_ICB     ', OCE_mas_berg_icb  , 'e', True)
970prtFlux ('  CALV_ICB     ', OCE_mas_calv_icb  , 'e', True)
971
972
973echo (' ')
974
975prtFlux ( 'wbilo sea  ', ATM_flux_int (ATM_wbilo_sea), 'e',   )
976prtFlux ( 'costalflow ', ONE_flux_int (RUN_coastalflow), 'e',   )
977prtFlux ( 'riverflow  ', RUN_flx_river  , 'e',   )
978prtFlux ( 'costalflow ', RUN_flx_coastal, 'e',   )
979prtFlux ( 'runoff     ', RUN_flx_river+RUN_flx_coastal, 'e',   )
980
981ATM_to_OCE   =  ATM_flux_int (ATM_wbilo_sea) - RUN_flx_river - RUN_flx_coastal - ATM_flx_calving
982#OCE_from_ATM = -OCE_mas_emp_oce - OCE_mas_emp_ice + OCE_mas_runoffs + OCE_mas_iceberg + OCE_mas_calving + OCE_mas_iceshelf
983OCE_from_ATM = OCE_mas_all
984
985prtFlux ( 'ATM_to_OCE  ', ATM_to_OCE  , 'e', True )
986prtFlux ( 'OCE_from_ATM', OCE_from_ATM, 'e', True )
Note: See TracBrowser for help on using the repository browser.