source: TOOLS/WATER_BUDGET/WaterUtils.py @ 6823

Last change on this file since 6823 was 6764, checked in by omamce, 3 months ago

O.M. : water budget - version mars 2024
  • Property svn:keywords set to Date Revision HeadURL Author
File size: 38.0 KB
Line 
1#!/usr/bin/env python3
2'''
3  Script to check water conservation in the IPSL coupled model
4
5  Here, some common utilities to all scripts
6
7  Warning, to install, configure, run, use any of included software or
8  to read the associated documentation you'll need at least one (1)
9  brain in a reasonably working order. Lack of this implement will
10  void any warranties (either express or implied).  Authors assumes
11  no responsability for errors, omissions, data loss, or any other
12  consequences caused directly or indirectly by the usage of his
13  software by incorrectly or partially configured personal
14
15 SVN information
16 $Author$
17 $Date$
18 $Revision$
19 $Id:  $
20 $HeadURL$
21'''
22
23import os, sys, types
24import functools
25import time
26import configparser, re
27import numpy as np
28import xarray as xr
29import libIGCM_sys, libIGCM_date as ldate
30import nemo
31import lmdz
32
33def ReadConfig (ini_file, default_ini_file='defaults.ini') :
34    '''
35    Reads experiment parameters
36
37    Reads <default_ini_file> config file to set all defaults parameters
38    Reads <ini_file> config file to set all experiment parameters
39    Reads config.card files if specified in <ini_file>
40
41    Returns a dictionary with all parameters
42    '''
43    ## Read file with defaults values
44    ## ------------------------------
45    params = configparser.ConfigParser ( interpolation=configparser.ExtendedInterpolation() )
46    params.optionxform = str # To keep capitals
47    params.read (default_ini_file)
48    print ( f'\nConfig file readed : {default_ini_file} ' )
49
50    ## Read Experiment config file
51    ## ----------------------------
52    exp_params = configparser.ConfigParser ( interpolation=configparser.ExtendedInterpolation() )
53    exp_params.optionxform = str # To keep capitals
54    exp_params.read (ini_file)
55    print ( f'\nConfig file readed : {ini_file} ' )
56   
57    ## Reading config.card if possible
58    ## -------------------------------
59    if 'Experiment' in params.keys ()  : ## Read Experiment on parameter file if possible
60        if 'ConfigCard' in exp_params['Experiment'].keys () :
61            ConfigCard = exp_params['Experiment']['ConfigCard']
62            print ( f'{ConfigCard=}' )
63           
64    if ConfigCard : ## Read config card if it exists
65        # Text existence of ConfigCard
66        if os.path.exists ( ConfigCard ) :
67            print ( f'Reading Config Card : {ConfigCard}' )
68            ## Creates parser for reading .ini input file
69            config = configparser.ConfigParser (interpolation=configparser.ExtendedInterpolation() )
70            config.optionxform = str # To keep capitals
71            config.read (ConfigCard)
72           
73            for VarName in ['JobName', 'ExperimentName', 'SpaceName', 'LongName', 'ModelName', 'TagName',
74                            'ORCA_version', 'ExpType', 'PeriodLength', 'ResolAtm', 'ResolOce', 'CalendarType' ] :
75                if VarName in config['UserChoices'].keys() :
76                    exp_params['Experiment'][VarName] = config['UserChoices'][VarName]
77                   
78            if 'Post' in config.sections() :
79                if 'PackFrequency' in config['Post'] :
80                    PackFrequency = config['Post']['PackFrequency']
81                    exp_params['Experiment']['PackFrequency'] = PackFrequency
82                   
83        else :
84            raise FileExistsError ( f"File not found : {ConfigCard = }" )
85       
86    ## Reading config file
87    ## -------------------
88    for Section in exp_params.keys () :
89        params[Section].update ( exp_params[Section] )
90    print ( f'\nConfig file readed : {ini_file} ' )
91
92    return config2dict (params)
93   
94def SetDatesAndFiles (pdict) :
95    '''
96    From readed experiment parameters, set all needed parameters
97    '''
98       
99    ## Set libIGCM and machine dependant values
100    ## ----------------------------------------
101    pdict['Experiment']['ICO']  = 'ICO' in pdict['Experiment']['ResolAtm']
102    #pdict['Experiment']['LMDZ'] = 'LMD' in pdict['Experiment']['ResolAtm']
103    pdict['Experiment']['LMDZ'] = not pdict['Experiment']['ICO']
104
105    mm = libIGCM_sys.config (
106        JobName        = pdict['Experiment']['JobName'],
107        ExperimentName = pdict['Experiment']['ExperimentName'],
108        TagName        = pdict['Experiment']['TagName'],
109        SpaceName      = pdict['Experiment']['SpaceName'],
110        User           = pdict['Experiment']['User'],
111        Group          = pdict['Experiment']['Group'],
112        TmpDir         = pdict['Files']['TmpDir'],
113        ARCHIVE        = pdict['libIGCM']['ARCHIVE'],
114        SCRATCHDIR     = pdict['libIGCM']['SCRATCHDIR'],
115        STORAGE        = pdict['libIGCM']['STORAGE'],
116        R_IN           = pdict['libIGCM']['R_IN'],
117        R_OUT          = pdict['libIGCM']['R_OUT'],
118        R_FIG          = pdict['libIGCM']['R_FIG'],
119        R_SAVE         = pdict['libIGCM']['R_SAVE'],
120        R_FIGR         = pdict['libIGCM']['R_FIGR'],
121        R_BUFR         = pdict['libIGCM']['R_BUFR'],
122        R_BUF_KSH      = pdict['libIGCM']['R_BUF_KSH'],
123        POST_DIR       = pdict['libIGCM']['POST_DIR'],
124        L_EXP          = pdict['libIGCM']['L_EXP']
125        )
126   
127    pdict['Files']['TmpDir'] = mm['TmpDir']
128    pdict['libIGCM'].update (mm.__dict__)
129       
130    ## Complete configuration
131    ## ----------------------
132    if not pdict['Experiment']['NEMO'] and pdict['Experiment']['ORCA_version'] :
133        if '1.2' in pdict['Experiment']['ORCA_version'] :
134            pdict['Experiment']['NEMO'] = 3.6
135        if '4.0' in pdict['Experiment']['ORCA_version'] :
136            pdict['Experiment']['NEMO'] = 4.0
137        if '4.2' in pdict['Experiment']['ORCA_version'] :
138            pdict['Experiment']['NEMO'] = 4.2
139           
140    if pdict['Experiment']['ATM_HIS'] :
141        if not pdict['Experiment']['SRF_HIS'] :
142            pdict['Experiment']['SRF_HIS'] = pdict['Experiment']['ATM_HIS']
143        if not pdict['Experiment']['RUN_HIS'] :
144            pdict['Experiment']['RUN_HIS'] = pdict['Experiment']['ATM_HIS']
145           
146    ##
147    ## Reading prec
148    if not pdict['Config']['readPrec'] :
149        pdict['Config']['readPrec'] = np.float64
150    else :
151        if pdict['Config']['readPrec'] in ["float", "float64", "r8", "double", "<class 'float'>"         ] :
152            pdict['Config']['readPrec'] = float
153        if pdict['Config']['readPrec'] in [         "float32", "r4", "single", "<class 'numpy.float32'>" ] :
154            pdict['Config']['readPrec'] = np.float32
155        if pdict['Config']['readPrec'] in [         "float16", "r2", "half"  , "<class 'numpy.float16'>" ] :
156            pdict['Config']['readPrec'] = np.float16
157
158    ## Defines begining and end of experiment
159    ## --------------------------------------
160    if not pdict['Experiment']['DateBegin'] :
161        pdict['Experiment']['DateBegin'] = f"{pdict['Experiment']['YearBegin']}0101"
162    else :
163        pdict['Experiment']['YearBegin'],pdict['Experiment']['MonthBegin'], pdict['Experiment']['DayBegin'] = ldate.GetYearMonthDay ( pdict['Experiment']['DateBegin'] )
164        pdict['Experiment']['DateBegin'] = ldate.ConvertFormatToGregorian (pdict['Experiment']['DateBegin'])
165       
166    if not pdict['Experiment']['DateEnd']  :
167        pdict['Experiment']['DateEnd'] = f"{pdict['Experiment']['YearEnd']}1231"
168    else :
169        pdict['Experiment']['YearEnd'], pdict['Experiment']['MonthEnd'], pdict['Experiment']['DayEnd'] = ldate.GetYearMonthDay ( pdict['Experiment']['DateEnd'] )
170        pdict['Experiment']['DateEnd'] = ldate.ConvertFormatToGregorian (pdict['Experiment']['DateEnd'])
171       
172    if not pdict['Experiment']['PackFrequency'] :
173        pdict['Experiment']['PackFrequency'] = f"{pdict['Experiment']['YearEnd'] - pdict['Experiment']['YearBegin'] + 1}YE"
174   
175    ## Set directory to extract files
176    ## ------------------------------
177    if not pdict['Files']['FileDir'] :
178         pdict['Files']['FileDir'] = os.path.join ( pdict['Files']['TmpDir'], f"WATER_{pdict['Experiment']['JobName']}" )
179    if not os.path.isdir ( pdict['Files']['FileDir'] ) : os.makedirs ( pdict['Files']['FileDir'] )
180
181    if not pdict['Files']['FileDirOCE'] :
182         pdict['Files']['FileDirOCE'] = os.path.join ( pdict['Files']['FileDir'], "OCE" )
183    if not os.path.isdir ( pdict['Files']['FileDirOCE'] ) : os.makedirs ( pdict['Files']['FileDirOCE'] )
184
185    if not pdict['Files']['FileDirICE'] :
186         pdict['Files']['FileDirICE'] = os.path.join ( pdict['Files']['FileDir'], "ICE" )
187    if not os.path.isdir ( pdict['Files']['FileDirICE'] ) : os.makedirs ( pdict['Files']['FileDirICE'] )
188       
189    FileOut = str2value (pdict['Files']['FileOut'])
190    if not FileOut :
191        FileOut = f"ATM_waterbudget_{pdict['Experiment']['JobName']}_{pdict['Experiment']['DateBegin']}_{pdict['Experiment']['DateEnd']}"
192        if pdict['Experiment']['ICO'] :
193            if pdict['Experiment']['ATM_HIS'] == 'latlon' : FileOut = f'{FileOut}_LATLON'
194            if pdict['Experiment']['ATM_HIS'] == 'ico'    : FileOut = f'{FileOut}_ICO'
195            if pdict['Config']['readPrec'] == np.float32  : FileOut = f'{FileOut}_float32'
196        FileOut = f'{FileOut}.out'
197       
198    pdict['Files']['FileOut'] = FileOut
199    f_out = open ( FileOut, mode = 'w', encoding="utf-8" )
200    pdict['Files']['f_out'] = f_out
201
202    echo (' ', f_out)
203    echo ( f"JobName : {pdict['Experiment']['JobName']}" , f_out=f_out )
204    echo ( f"Comment : {pdict['Experiment']['Comment']}" , f_out=f_out )
205    echo ( f"TmpDir  : {pdict['Files']['TmpDir']}"       , f_out=f_out )
206    echo ( f"FileDir : {pdict['Files']['FileDir']}"      , f_out=f_out )
207   
208    echo ( f"\nDealing with {pdict['libIGCM']['L_EXP']=}", f_out )
209
210    ## Creates model output directory names
211    ## ------------------------------------
212    if not pdict['Files']['FreqDir'] : 
213        if pdict['Experiment']['Freq'] == "MO" :
214            pdict['Files']['FreqDir'] = os.path.join ('Output' , 'MO' )
215        if pdict['Experiment']['Freq'] == "SE" :
216            pdict['Files']['FreqDir'] = os.path.join ('Analyse', 'SE' )
217   
218    if not pdict['Files']['dir_ATM_his'] :
219        pdict['Files']['dir_ATM_his'] = os.path.join ( pdict['libIGCM']['R_SAVE'], "ATM", pdict['Files']['FreqDir'] )
220    if not pdict['Files']['dir_SRF_his'] :
221        pdict['Files']['dir_SRF_his'] = os.path.join ( pdict['libIGCM']['R_SAVE'], "SRF", pdict['Files']['FreqDir'] )
222       
223    if not pdict['Files']['dir_OCE_his'] :
224        pdict['Files']['dir_OCE_his'] = os.path.join ( pdict['libIGCM']['R_SAVE'], "OCE", pdict['Files']['FreqDir'] )
225    if not pdict['Files']['dir_ICE_his'] :
226        pdict['Files']['dir_ICE_his'] = os.path.join ( pdict['libIGCM']['R_SAVE'], "ICE", pdict['Files']['FreqDir'] )
227       
228    echo (  'The analysis relies on files from the following model output directories : ', f_out )
229    echo ( f"{pdict['Files']['dir_ATM_his'] = }" , f_out )
230    echo ( f"{pdict['Files']['dir_SRF_his'] = }" , f_out )
231    echo ( f"{pdict['Files']['dir_OCE_his'] = }" , f_out )
232    echo ( f"{pdict['Files']['dir_ICE_his'] = }" , f_out )
233   
234    ## Creates files names
235    ## --------------------
236    if not pdict['Experiment']['Period'] :
237        if pdict['Experiment']['Freq'] == 'MO' :
238            pdict['Experiment']['Period'] = f"{pdict['Experiment']['DateBegin']}_{pdict['Experiment']['DateEnd']}_1M"
239        if pdict['Experiment']['Freq'] == 'SE' :
240            pdict['Experiment']['Period'] = f"SE_{pdict['Files']['DateBegin']}_{pdict['Files']['DateEnd']}_1M"
241   
242    echo ( f"Period   : {pdict['Experiment']['Period']}", f_out )
243   
244    if not pdict['Files']['FileCommon']  :
245        pdict['Files']['FileCommon'] = f"{pdict['Experiment']['JobName']}_{pdict['Experiment']['Period']}"
246       
247    if not pdict['Files']['Title'] :
248        pdict['Files']['Title'] = f"{pdict['Experiment']['JobName']} : {pdict['Experiment']['Freq']} : {pdict['Experiment']['DateBegin']} - {pdict['Experiment']['DateEnd']}"
249       
250    echo ('\nOpen history files', f_out )
251    if not pdict['Files']['file_ATM_his'] :
252        if pdict['Experiment']['ATM_HIS'] == 'latlon' :
253            pdict['Files']['file_ATM_his'] = os.path.join ( pdict['Files']['dir_ATM_his'], f"{pdict['Files']['FileCommon']}_histmth.nc" )
254        if pdict['Experiment']['ATM_HIS'] == 'ico' :
255            pdict['Files']['file_ATM_his'] = os.path.join ( pdict['Files']['dir_ATM_his'], f"{pdict['Files']['FileCommon']}_histmth_ico.nc" )
256    if not pdict['Files']['file_SRF_his'] :
257        if pdict['Experiment']['ATM_HIS'] == 'latlon' :
258            pdict['Files']['file_SRF_his'] = os.path.join ( pdict['Files']['dir_SRF_his'], f"{pdict['Files']['FileCommon']}_sechiba_history.nc" )
259        if pdict['Experiment']['ATM_HIS'] == 'ico' :
260            pdict['Files']['file_SRF_his'] = os.path.join ( pdict['Files']['dir_SRF_his'], f"{pdict['Files']['FileCommon']}_sechiba_history_ico.nc" )
261       
262    if pdict['Experiment']['SECHIBA'] and pdict['Experiment']['Routing'] == 'SIMPLE' :
263        if not pdict['Files']['file_RUN_his'] :
264            if pdict['Experiment']['RUN_HIS'] == 'latlon' :
265                pdict['Files']['file_RUN_his'] = os.path.join ( pdict['Files']['dir_SRF_his'], f"{pdict['Files']['FileCommon']}_sechiba_history.nc" )
266            if pdict['Experiment']['SRF_HIS'] == 'ico' :
267                pdict['Files']['file_RUN_his'] = os.path.join ( pdict['Files']['dir_SRF_his'], f"{pdict['Files']['FileCommon']}_sechiba_history_ico.nc" )
268           
269    echo ( f"{pdict['Files']['file_ATM_his'] = }", f_out )
270    if pdict['Experiment']['SECHIBA'] :
271        echo ( f"{pdict['Files']['file_SRF_his'] = }", f_out )
272        if pdict['Experiment']['Routing'] == 'SIMPLE' : echo ( f"{pdict['Files']['file_RUN_his'] = }", f_out )
273
274    if not pdict['Files']['file_OCE_his'] :
275        pdict['Files']['file_OCE_his'] = os.path.join ( pdict['Files']['dir_OCE_his'], f"{pdict['Files']['FileCommon']}_grid_T.nc" )
276    if not pdict['Files']['file_OCE_sca'] :
277        pdict['Files']['file_OCE_sca'] = os.path.join ( pdict['Files']['dir_OCE_his'], f"{pdict['Files']['FileCommon']}_scalar.nc" )
278    if not pdict['Files']['file_OCE_srf'] :
279        pdict['Files']['file_OCE_srf'] = os.path.join ( pdict['Files']['dir_OCE_his'], f"{pdict['Files']['FileCommon']}_sbc.nc" )
280    if not pdict['Files']['file_ICE_his'] :
281        pdict['Files']['file_ICE_his'] = os.path.join ( pdict['Files']['dir_ICE_his'], f"{pdict['Files']['FileCommon']}_icemod.nc" )
282   
283    return pdict
284
285def SetRestartNames (pdict, f_out) :
286    '''
287    Defines dates for restart files
288    Define names of tar files with restart
289
290    Returns a dictionary
291    '''
292   
293    if not pdict['Files']['TarRestartDate_beg'] :
294        pdict['Files']['TarRestartDate_beg'] = ldate.SubOneDayToDate ( pdict['Experiment']['DateBegin'], pdict['Experiment']['CalendarType'] )
295    if not pdict['Files']['TarRestartDate_end'] :
296        pdict['Files']['TarRestartDate_end'] = ldate.ConvertFormatToGregorian ( pdict['Experiment']['DateEnd'] )
297   
298    if not pdict['Files']['TarRestartPeriod_beg'] :
299       
300        pdict['Files']['TarRestartPeriod_beg_DateEnd'] = pdict['Files']['TarRestartDate_beg']
301        pdict['Files']['TarRestartPeriod_beg_DateBeg'] = ldate.DateAddPeriod ( pdict['Files']['TarRestartPeriod_beg_DateEnd'],
302                                                                                   f"-{pdict['Experiment']['PackFrequency']}")
303        pdict['Files']['TarRestartPeriod_beg_DateBeg'] = ldate.AddOneDayToDate ( pdict['Files']['TarRestartPeriod_beg_DateBeg'],
304                                                                                     pdict['Experiment']['CalendarType']  )
305       
306        pdict['Files']['TarRestartPeriod_beg'] = f"{pdict['Files']['TarRestartPeriod_beg_DateBeg']}_{pdict['Files']['TarRestartPeriod_beg_DateEnd']}"
307        echo ( f"Tar period for initial restart : {pdict['Files']['TarRestartPeriod_beg']}", f_out)
308       
309    if not pdict['Files']['TarRestartPeriod_end'] :
310
311        pdict['Files']['TarRestartPeriod_end_DateEnd'] = pdict['Files']['TarRestartDate_end']
312        pdict['Files']['TarRestartPeriod_end_DateBeg'] = ldate.DateAddPeriod ( pdict['Files']['TarRestartPeriod_end_DateEnd'],
313                                                                                   f"-{pdict['Experiment']['PackFrequency']}" )
314        pdict['Files']['TarRestartPeriod_end_DateBeg'] = ldate.AddOneDayToDate ( pdict['Files']['TarRestartPeriod_end_DateBeg'], pdict['Experiment']['CalendarType'] )
315       
316        pdict['Files']['TarRestartPeriod_end'] = f"{pdict['Files']['TarRestartPeriod_end_DateBeg']}_{pdict['Files']['TarRestartPeriod_end_DateEnd']}"
317        echo ( f"Tar period for final restart : {pdict['Files']['TarRestartPeriod_end']}", f_out )
318       
319    echo ( f"Restart dates - Start : {pdict['Files']['TarRestartPeriod_beg']}  /  End : {pdict['Files']['TarRestartPeriod_end']}", f_out)
320       
321    if not pdict['Files']['tar_restart_beg'] :
322        pdict['Files']['tar_restart_beg'] = os.path.join ( pdict['libIGCM']['R_SAVE'], 'RESTART', f"{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartPeriod_beg']}_restart.tar" )
323    if not  pdict['Files']['tar_restart_end'] :
324        pdict['Files']['tar_restart_end'] = os.path.join ( pdict['libIGCM']['R_SAVE'], 'RESTART', f"{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartPeriod_end']}_restart.tar" )
325       
326    echo ( f"{ pdict['Files']['tar_restart_beg'] = }", f_out )
327    echo ( f"{ pdict['Files']['tar_restart_end'] = }", f_out )
328
329    ##-- Names of tar files with restarts
330   
331    if not pdict['Files']['tar_restart_beg_ATM'] :
332        pdict['Files']['tar_restart_beg_ATM'] = pdict['Files']['tar_restart_beg']
333    if not pdict['Files']['tar_restart_beg_DYN'] :
334        pdict['Files']['tar_restart_beg_DYN'] = pdict['Files']['tar_restart_beg']
335    if not pdict['Files']['tar_restart_beg_RUN'] :
336        pdict['Files']['tar_restart_beg_RUN'] = pdict['Files']['tar_restart_beg']
337    if not pdict['Files']['tar_restart_beg_OCE'] :
338        pdict['Files']['tar_restart_beg_OCE'] = pdict['Files']['tar_restart_beg']
339    if not pdict['Files']['tar_restart_beg_ICE'] :
340        pdict['Files']['tar_restart_beg_ICE'] = pdict['Files']['tar_restart_beg']
341       
342    if not pdict['Files']['tar_restart_end_ATM'] :
343        pdict['Files']['tar_restart_end_ATM'] = pdict['Files']['tar_restart_end']
344    if not pdict['Files']['tar_restart_end_DYN'] :
345        pdict['Files']['tar_restart_end_DYN'] = pdict['Files']['tar_restart_end']
346    if not pdict['Files']['tar_restart_end_RUN'] :
347        pdict['Files']['tar_restart_end_RUN'] = pdict['Files']['tar_restart_end']
348    if not pdict['Files']['tar_restart_end_OCE'] :
349        pdict['Files']['tar_restart_end_OCE'] = pdict['Files']['tar_restart_end']
350    if not pdict['Files']['tar_restart_end_ICE'] :
351        pdict['Files']['tar_restart_end_ICE'] = pdict['Files']['tar_restart_end']
352       
353    if not pdict['Files']['file_DYN_beg'] :
354        if pdict['Experiment']['LMDZ'] :
355            pdict['Files']['file_DYN_beg'] = f"{pdict['Files']['FileDir']}/ATM_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_beg']}_restart.nc"
356        if pdict['Experiment']['ICO']  :
357            pdict['Files']['file_DYN_beg'] = f"{pdict['Files']['FileDir']}/ICO_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_beg']}_restart.nc"
358       
359    if not pdict['Files']['file_DYN_end'] :
360        if pdict['Experiment']['LMDZ'] :
361            pdict['Files']['file_DYN_end'] = f"{pdict['Files']['FileDir']}/ATM_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_end']}_restart.nc"
362        if pdict['Experiment']['ICO']  :
363            pdict['Files']['file_DYN_end'] = f"{pdict['Files']['FileDir']}/ICO_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_end']}_restart.nc"
364       
365    if pdict['Experiment']['SECHIBA'] :
366        if not pdict['Files']['file_SRF_beg'] :
367            pdict['Files']['file_SRF_beg'] = f"{pdict['Files']['FileDir']}/SRF_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_beg']}_sechiba_rest.nc"
368        if not pdict['Files']['file_SRF_end'] :
369            pdict['Files']['file_SRF_end'] = f"{pdict['Files']['FileDir']}/SRF_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_end']}_sechiba_rest.nc"
370       
371    if pdict['Experiment']['SECHIBA'] :
372        if not pdict['Files']['tar_restart_beg_SRF'] :
373            pdict['Files']['tar_restart_beg_SRF'] = pdict['Files']['tar_restart_beg']
374        if not pdict['Files']['tar_restart_end_SRF'] :
375            pdict['Files']['tar_restart_end_SRF'] = pdict['Files']['tar_restart_end']
376           
377    if not pdict['Files']['file_ATM_beg'] :
378        pdict['Files']['file_ATM_beg'] = f"{pdict['Files']['FileDir']}/ATM_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_beg']}_restartphy.nc"
379    if not pdict['Files']['file_ATM_end'] :
380        pdict['Files']['file_ATM_end'] = f"{pdict['Files']['FileDir']}/ATM_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_end']}_restartphy.nc"
381       
382    if pdict['Experiment']['SECHIBA'] :
383        echo ( f"{pdict['Files']['file_SRF_beg'] = }", f_out)
384        echo ( f"{pdict['Files']['file_SRF_end'] = }", f_out)
385       
386    if pdict['Experiment']['ICO'] :
387        if not pdict['Files']['file_DYN_aire'] :
388            pdict['Files']['file_DYN_aire'] = os.path.join ( pdict['libIGCM']['R_IN'], 'ATM', 'GRID',  f"{pdict['Experiment']['ResolAtm']}_grid.nc" )
389       
390    if pdict['Experiment']['SECHIBA'] and pdict['Experiment']['Routing'] == 'SIMPLE' :
391        if not pdict['Files']['file_RUN_beg'] :
392            pdict['Files']['file_RUN_beg'] = f"{ pdict['Files']['FileDir']}/SRF_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_beg']}_routing_restart.nc"
393        if not pdict['Files']['file_RUN_end'] :
394            pdict['Files']['file_RUN_end'] = f"{ pdict['Files']['FileDir']}/SRF_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_end']}_routing_restart.nc"
395
396    if not pdict['Files']['tar_restart_beg_OCE'] :
397        pdict['Files']['tar_restart_beg_OCE'] = pdict['Files']['tar_restart_beg']
398    if not pdict['Files']['tar_restart_beg_ICE'] :
399        pdict['Files']['tar_restart_beg_ICE'] = pdict['Files']['tar_restart_beg']
400       
401    if not pdict['Files']['tar_restart_end_OCE'] :
402        pdict['Files']['tar_restart_end_OCE'] = pdict['Files']['tar_restart_end']
403    if not pdict['Files']['tar_restart_end_ICE'] :
404        pdict['Files']['tar_restart_end_ICE'] = pdict['Files']['tar_restart_end']
405       
406    if not pdict['Files']['file_OCE_beg'] :
407        pdict['Files']['file_OCE_beg'] = f"{pdict['Files']['FileDir']}/OCE_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_beg']}_restart.nc"
408    if not pdict['Files']['file_OCE_end'] :
409        pdict['Files']['file_OCE_end'] = f"{pdict['Files']['FileDir']}/OCE_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_end']}_restart.nc"
410    if not pdict['Files']['file_ICE_beg'] :
411        pdict['Files']['file_ICE_beg'] = f"{pdict['Files']['FileDir']}/ICE_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_beg']}_restart_icemod.nc"
412    if not pdict['Files']['file_ICE_end'] :
413        pdict['Files']['file_ICE_end'] = f"{pdict['Files']['FileDir']}/ICE_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_end']}_restart_icemod.nc"
414           
415    return pdict
416
417def ComputeGridATM ( dpar, d_ATM_his ) :
418    readPrec = dpar['Config']['readPrec']
419    if repr(readPrec) in [ "<class 'numpy.float64'>" , float ] :
420        def rprec (ptab) :
421            return ptab
422    else                 :
423        def rprec (ptab) :
424            return ptab.astype(readPrec).astype(float)
425   
426    ATM = types.SimpleNamespace ()
427
428    # ATM grid with cell surfaces
429    if dpar['Experiment']['LMDZ'] :
430        #echo ('ATM grid with cell surfaces : LMDZ', f_out)
431        ATM.lat       = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1d='cell' )
432        ATM.lon       = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1d='cell' )
433        ATM.aire      = lmdz.geo2point ( rprec (d_ATM_his ['aire']     [0]), cumul_poles=True, dim1d='cell' )
434        ATM.fter      = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1d='cell' )
435        ATM.foce      = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1d='cell' )
436        ATM.fsic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1d='cell' )
437        ATM.flic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1d='cell' )
438       
439    if  dpar['Experiment']['ICO'] :
440        if  dpar['Experiment']['ATM_HIS'] == 'latlon' :
441            #echo ( 'ATM areas and fractions on LATLON grid' )
442            if 'lat_dom_out' in d_ATM_his.variables :
443                ATM.lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat_dom_out'])+0*rprec (d_ATM_his ['lon_dom_out']), dim1d='cell' )
444                ATM.lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat_dom_out'])+  rprec (d_ATM_his ['lon_dom_out']), dim1d='cell' )
445            else :
446                ATM.lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1d='cell' )
447                ATM.lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1d='cell' )
448            ATM.aire = lmdz.geo2point ( rprec (d_ATM_his ['aire'][0]).squeeze(), cumul_poles=True, dim1d='cell' )
449            ATM.fter = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1d='cell' )
450            ATM.foce = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1d='cell' )
451            ATM.fsic = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1d='cell' )
452            ATM.flic = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1d='cell' )
453
454           
455        if  dpar['Experiment']['ATM_HIS'] == 'ico' :
456            #echo ( 'ATM areas and fractions on ICO grid' )
457            ATM.aire =  rprec (d_ATM_his ['aire']     [0]).squeeze()
458            ATM.lat  =  rprec (d_ATM_his ['lat']         )
459            ATM.lon  =  rprec (d_ATM_his ['lon']         )
460            ATM.fter =  rprec (d_ATM_his ['fract_ter'][0])
461            ATM.foce =  rprec (d_ATM_his ['fract_oce'][0])
462            ATM.fsic =  rprec (d_ATM_his ['fract_sic'][0])
463            ATM.flic =  rprec (d_ATM_his ['fract_lic'][0])
464           
465       
466    ATM.fsea      = ATM.foce + ATM.fsic
467    ATM.flnd      = ATM.fter + ATM.flic
468    ATM.aire_fter = ATM.aire * ATM.fter
469    ATM.aire_flic = ATM.aire * ATM.flic
470    ATM.aire_fsic = ATM.aire * ATM.fsic
471    ATM.aire_foce = ATM.aire * ATM.foce
472    ATM.aire_flnd = ATM.aire * ATM.flnd
473    ATM.aire_fsea = ATM.aire * ATM.fsea
474
475    ATM.aire_sea     = ATM.aire * ATM.fsea
476    ATM.aire_tot     = P1sum ( ATM.aire      )
477    ATM.aire_sea_tot = P1sum ( ATM.aire_fsea )
478    ATM.aire_ter_tot = P1sum ( ATM.aire_fter )
479    ATM.aire_lic_tot = P1sum ( ATM.aire_flic )
480   
481    return dpar, ATM
482
483def ComputeGridSRF ( dpar, d_SRF_his ) :
484    readPrec = dpar['Config']['readPrec']
485    if repr(readPrec) in [ "<class 'numpy.float64'>" , float ] :
486        def rprec (ptab) :
487            return ptab
488    else                 :
489        def rprec (ptab) :
490            return ptab.astype(readPrec).astype(float)
491   
492    SRF = types.SimpleNamespace ()
493   
494    if dpar['Experiment']['SECHIBA'] : 
495        if dpar['Experiment']['LMDZ'] :
496            SRF.lat       = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1d='cell' )
497            SRF.lon       = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1d='cell' )
498            SRF.aire      = lmdz.geo2point ( rprec (d_SRF_his ['Areas']) * rprec (d_SRF_his ['Contfrac']), dim1d='cell', cumul_poles=True )
499            SRF.areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas'])  ,  dim1d='cell', cumul_poles=True )
500            SRF.contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac']), dim1d='cell' )
501   
502        if dpar['Experiment']['ICO'] :
503            if dpar['Experiment']['SRF_HIS'] == 'latlon' :
504                #echo ( 'SRF areas and fractions on LATLON grid' )
505                if 'lat_domain_landpoints_out' in d_SRF_his  :
506                    SRF.lat = lmdz.geo2point (   rprec (d_SRF_his ['lat_domain_landpoints_out'])+0*rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1d='cell' )
507                    SRF.lon = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_domain_landpoints_out'])+  rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1d='cell' )
508                else :
509                    if 'lat_domain_landpoints_out' in d_SRF_his :
510                        SRF.lat = lmdz.geo2point (   rprec (d_SRF_his ['lat_dom_out'])+0*rprec (d_SRF_his ['lon_dom_out']), dim1d='cell' )
511                        SRF.lon = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_dom_out'])+  rprec (d_SRF_his ['lon_dom_out']), dim1d='cell' )
512                    else :
513                        SRF.lat = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1d='cell' )
514                        SRF.lon = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1d='cell' )
515   
516                SRF.areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas']   )   , dim1d='cell', cumul_poles=True )
517                SRF.areafrac  = lmdz.geo2point ( rprec (d_SRF_his ['AreaFrac'])   , dim1d='cell', cumul_poles=True )
518                SRF.contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac'])   , dim1d='cell', cumul_poles=True )
519                SRF.aire      = SRF.areafrac
520   
521            if dpar['Experiment']['SRF_HIS'] == 'ico' :
522                #echo ( 'SRF areas and fractions on ICO grid' )
523                SRF.lat       =  rprec (d_SRF_his ['lat']     )
524                SRF.lon       =  rprec (d_SRF_his ['lon']     )
525                SRF.areas     =  rprec (d_SRF_his ['Areas']   )
526                SRF.contfrac  =  rprec (d_SRF_his ['Contfrac'])
527                SRF.aire      =  SRF.areas * SRF.contfrac
528
529        SRF.aire_tot = P1sum ( SRF.aire )
530
531    return dpar, SRF
532
533def ComputeGridDYN ( dpar, ATM, d_DYN_aire, d_ATM_beg ) : 
534    readPrec = dpar['Config']['readPrec']
535    if repr(readPrec) in [ "<class 'numpy.float64'>" , float ] :
536        def rprec (ptab) :
537            return ptab
538    else                 :
539        def rprec (ptab) :
540            return ptab.astype(readPrec).astype(float)
541       
542    DYN = types.SimpleNamespace ()
543    if dpar['Experiment']['ICO'] :
544        if dpar['Config']['SortIco'] :
545            # Creation d'une clef de tri pour le fichier aire
546            DYN.aire_keysort = np.lexsort ( (d_DYN_aire['lat'], d_DYN_aire['lon']) )
547        else :
548            DYN.aire_keysort = np.arange ( len ( d_DYN_aire['lat'] ) )
549           
550        DYN.lat = d_DYN_aire['lat']
551        DYN.lon = d_DYN_aire['lon']
552       
553        DYN.aire = d_DYN_aire['aire']
554        DYN.fsea = d_DYN_aire['fract_oce'] + d_DYN_aire['fract_sic']
555       
556        DYN.flnd = 1.0 - DYN.fsea
557        DYN.fter = d_ATM_beg['FTER']
558        DYN.flic = d_ATM_beg['FLIC']
559        DYN.foce = d_ATM_beg['FOCE']
560        DYN.aire_fter = DYN.aire * DYN.fter
561       
562    if dpar['Experiment']['ATM_HIS'] == 'ico' :
563        DYN.aire      = ATM.aire
564        DYN.foce      = ATM.foce
565        DYN.fsic      = ATM.fsic
566        DYN.flic      = ATM.flic
567        DYN.fter      = ATM.fter
568        DYN.fsea      = ATM.fsea
569        DYN.flnd      = ATM.flnd
570        DYN.aire_fter = ATM.aire_fter
571        DYN.aire_flic = ATM.aire_flic
572        DYN.aire_fsic = ATM.aire_fsic
573        DYN.aire_foce = ATM.aire_foce
574        DYN.aire_flnd = ATM.aire_flnd
575        DYN.aire_fsea = ATM.aire_fsea
576   
577    if dpar['Experiment']['LMDZ'] :
578        # Area on lon/lat grid
579        DYN.aire = ATM.aire
580        DYN.fsea = ATM.fsea
581        DYN.flnd = ATM.flnd
582        DYN.fter = rprec (d_ATM_beg['FTER'])
583        DYN.flic = rprec (d_ATM_beg['FLIC'])
584        DYN.aire_fter = DYN.aire * DYN.fter
585       
586    DYN.aire_tot =  P1sum ( DYN.aire ) 
587   
588    return dpar, DYN
589
590def ComputeGridOCE ( dpar, d_OCE_his, nperio=None ) :
591    OCE = types.SimpleNamespace ()
592
593    # Get mask and surfaces
594    sos = d_OCE_his ['sos'][0].squeeze()
595    OCE.OCE_msk = nemo.lbc_mask ( xr.where ( sos>0., 1., 0. ), cd_type = 'T', nperio=nperio, sval = 0. )
596   
597    so = sos = d_OCE_his ['sos'][0].squeeze()
598    OCE.msk3 = nemo.lbc_mask ( xr.where ( so>0., 1., 0. ), cd_type = 'T', nperio=nperio, sval = 0. )
599   
600    # lbc_mask removes the duplicate points (periodicity and north fold)
601    OCE.OCE_aire = nemo.lbc_mask ( d_OCE_his ['area'] * OCE.OCE_msk, cd_type = 'T', nperio=nperio, sval = 0.0 )
602    OCE.ICE_aire = OCE.OCE_aire
603
604    OCE.OCE_aire_tot = P1sum ( OCE.OCE_aire )
605    OCE.ICE_aire_tot = P1sum ( OCE.ICE_aire )
606
607    return dpar, OCE
608
609def config2dict ( pconf ) :
610    '''
611    Convert a config parser object into a dictionary
612    '''
613    pdict = {}
614    for section in pconf.keys () :
615        pdict[section] = {}
616        for key in pconf[section].keys() :
617            zz = str2value ( pconf[section][key] )
618            pdict[section].update ( {key:zz} )       
619    return pdict
620
621def dict2config ( pdict ):
622    '''
623    Convert a dictionary into a configparser object
624    All values are converted to strings
625
626    The dictionary must have two levels (see configparser)
627    '''
628    zconf = configparser.ConfigParser (interpolation=configparser.ExtendedInterpolation ())
629    zconf.optionxform = str # To keep capitals
630    zconf.read_dict(dict2str(pdict))
631   
632    return zconf
633
634def updateConf ( pconf, pdict ):
635    '''
636    Update a config parser with a dictionary
637    All values are converted to strings
638
639    The dictionary must have two levels (see configparser)
640    '''
641    pconf.read_dict (dict1str(pdict))
642    #for section in pdict.keys() :
643    #    if section not in pconf.keys() : pconf[section] = {}
644    #    for key in pdict[section].keys() :
645    #        pconf[section][key] = str(pdict[section][key])
646    return pconf
647
648def dict2str (pdict) :
649    '''
650    Convert all values of a dictionary to strings
651    '''
652    zout = {}
653    for section in pdict.keys() :
654        zout[section] = {}
655        for key in pdict[section].keys() :
656            zout[section][key] = str(pdict[section][key])
657    return zout
658
659def echo (string, f_out, end='\n') :
660    '''Function to print to stdout *and* output file'''
661    print ( str(string), end=end  )
662    sys.stdout.flush ()
663    f_out.write ( str(string) + end )
664    f_out.flush ()
665
666def str2value ( pz ) :
667    '''
668    Tries to convert a string into integer, real, boolean or None
669    '''
670    zz = pz
671    zz = setBool (zz)
672    zz = setNum  (zz)
673    zz = setNone (zz)
674    return zz
675
676def setBool (chars) :
677    '''Convert specific char string in boolean if possible'''
678    z_set_bool = chars
679    if isinstance (chars, str) :
680        for key in configparser.ConfigParser.BOOLEAN_STATES.keys () :
681            if chars.lower() == key :
682                z_set_bool = configparser.ConfigParser.BOOLEAN_STATES[key]
683    return z_set_bool
684
685def setNum (chars) :
686    '''Convert specific char string in integer or real if possible'''
687    zset_num = chars
688    if isinstance (chars, str) :
689        realnum   = re.compile ("^[-+]?[0-9]*\.?[0-9]+(e[-+]?[0-9]+)?$")
690        is_number = realnum.match(chars.strip()) != None
691        is_int    = chars.strip().isdigit()
692        if is_number :
693            if is_int : zset_num = int   (chars)
694            else      : zset_num = float (chars)
695   
696    return zset_num
697
698def setNone (chars) :
699    '''Convert specific char string to None if possible'''
700    zset_none = chars
701    if isinstance (chars, str) :
702        if chars in ['None', 'NONE', 'none'] : zset_none = None
703    return zset_none
704
705def unDefined (char) :
706    '''True if a variable is not defined, or set to None'''
707    if char in globals () :
708        if globals()[char] :
709            zun_defined = False
710        else : zun_defined = True
711    else :     zun_defined = True
712    return zun_defined
713
714def Defined (char) :
715    '''True if a variable is defined and not equal to None'''
716    if char in globals () :
717        if globals()[char] is None :
718            zdefined = False
719        else            : zdefined = True
720    else                : zdefined = False
721    return zdefined
722
723def Ksum (tab) :
724    '''Kahan summation algorithm, also known as compensated summation
725
726    https://en.wikipedia.org/wiki/Kahan_summation_algorithm
727    '''
728    # Prepare the accumulator.
729    ksum = 0.0
730    # A running compensation for lost low-order bits.
731    comp = 0.0
732
733    for xx in np.where ( np.isnan(tab), 0., tab ) :
734        # comp is zero the first time around.
735        yy = xx - comp
736        # Alas, sum is big, y small, so low-order digits of y are lost.
737        tt = ksum + yy
738        # (tt - Ksum) cancels the high-order part of y; subtracting y recovers negative (low part of yy)
739        comp = (tt - ksum) - yy
740        # Algebraically, comp should always be zero. Beware overly-aggressive optimizing compilers !
741        ksum = tt
742        # Next time around, the lost low part will be added to y in a fresh attempt.
743
744    return ksum
745
746def Ssum (tab) :
747    '''Precision summation by sorting by absolute values'''
748    ssum = np.sum ( tab[np.argsort(np.abs(tab))] )
749    return ssum
750
751def KSsum (tab) :
752    '''Precision summation by sorting by absolute value, and applying Kahan summation'''
753    kssum = Ksum ( tab[np.argsort(np.abs(tab))] )
754    return kssum
755
756# Choosing algorithm for precise summation
757Psum = KSsum
758
759def P1sum (ptab) :
760    return Psum ( ptab.to_masked_array().ravel() )
761
762class Timer :
763    '''Measures execution time of a function'''
764    def __str__ (self):
765        return str (self.__class__)
766   
767    def __name__ (self):
768        return self.__class__.__name__
769
770    def __init__ (self, func, debug=False, timer=True) :
771        functools.update_wrapper (self, func)
772        self.func       = func
773        self.num_calls  = 0
774        self.cumul_time = 0.
775        self.debug      = debug
776        self.timer      = timer
777
778    def __call__ (self, *args, **kwargs) :
779        self.num_calls += 1
780        if self.debug :
781            print ( f'>-- Calling {self.__name__} --------------------------------------------------' )
782            args_repr   = [f"{repr (a)}" for a in args]
783            kwargs_repr = [f"{k}={v!r}" for k, v in kwargs.items ()]
784            signature   = ", ".join (args_repr + kwargs_repr)
785            print ( f"Signature : ({signature}")
786        start_time = time.perf_counter ()     
787        values     = self.func (*args, **kwargs)
788        end_time   = time.perf_counter ()
789        run_time   = end_time - start_time
790        self.cumul_time += run_time
791        if self.timer : 
792            print ( f"--> Calling {self.__name__!r} : {run_time:.3f} secs / cumul {self.cumul_time:.3f} # {self.num_calls:2d} calls")
793        if self.debug :
794            print ( '<------------------------------------------------------------' )
795        return values
796   
Note: See TracBrowser for help on using the repository browser.