source: TOOLS/WATER_BUDGET/WaterUtils.py @ 6709

Last change on this file since 6709 was 6688, checked in by omamce, 7 months ago

O.M. :

WATER_BUDGET : put function common to ATM and OCE in WAterUtils library

use of dictionary for all input parameters
  • Property svn:keywords set to Date Revision HeadURL Author
File size: 37.4 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)
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
159    ## Defines begining and end of experiment
160    ## --------------------------------------
161    if not pdict['Experiment']['DateBegin'] :
162        pdict['Experiment']['DateBegin'] = f"{pdict['Experiment']['YearBegin']}0101"
163    else :
164        YearBegin, MonthBegin, DayBegin = ldate.GetYearMonthDay ( pdict['Experiment']['DateBegin'] )
165        DateBegin = FormatToGregorian (pdict['Experiment']['DateBegin'])
166        pdict['Experiment']['YearBegin'] = YearBegin
167       
168    if not pdict['Experiment']['DateEnd']  :
169        pdict['Experiment']['DateEnd'] = f"{pdict['Experiment']['YearEnd']}1231"
170    else :
171        YearEnd, MonthEnd, DayEnd = ldate.GetYearMonthDay ( pdict['Experiment']['DateEnd'] )
172        DateEnd   = FormatToGregorian (pdict['Experiment']['DateEnd'])
173        pdict['Experiment']['DateEnd'] = DateEnd
174       
175    if not pdict['Experiment']['PackFrequency'] :
176        PackFrequencypdict['Experiment']['PackFrequency'] = f"{YearEnd - YearBegin + 1}YE"
177   
178
179    ## Set directory to extract files
180    ## ------------------------------
181    if not pdict['Files']['FileDir'] :
182         pdict['Files']['FileDir'] = os.path.join ( pdict['Files']['TmpDir'], f"WATER_{pdict['Experiment']['JobName']}" )
183   
184    if not os.path.isdir ( pdict['Files']['FileDir'] ) : os.makedirs ( pdict['Files']['FileDir'] )
185
186    FileOut = str2value (pdict['Files']['FileOut'])
187    if not FileOut :
188        FileOut = f"ATM_waterbudget_{pdict['Experiment']['JobName']}_{pdict['Experiment']['YearBegin']}_{pdict['Experiment']['YearEnd']}"
189        if pdict['Experiment']['ICO'] :
190            if pdict['Experiment']['ATM_HIS'] == 'latlon' : FileOut = f'{FileOut}_LATLON'
191            if pdict['Experiment']['ATM_HIS'] == 'ico'    : FileOut = f'{FileOut}_ICO'
192            if pdict['Config']['readPrec'] == np.float32  : FileOut = f'{FileOut}_float32'
193        FileOut = f'{FileOut}.out'
194       
195    pdict['Files']['FileOut'] = FileOut
196    f_out = open ( FileOut, mode = 'w', encoding="utf-8" )
197    pdict['Files']['f_out'] = f_out
198
199    echo (' ', f_out)
200    echo ( f"JobName : {pdict['Experiment']['JobName']}" , f_out=f_out )
201    echo ( f"Comment : {pdict['Experiment']['Comment']}" , f_out=f_out )
202    echo ( f"TmpDir  : {pdict['Files']['TmpDir']}"       , f_out=f_out )
203    echo ( f"FileDir : {pdict['Files']['FileDir']}"      , f_out=f_out )
204   
205    echo ( f"\nDealing with {pdict['libIGCM']['L_EXP']}", f_out )
206
207    ## Creates model output directory names
208    ## ------------------------------------
209    if not pdict['Files']['FreqDir'] : 
210        if pdict['Experiment']['Freq'] == "MO" :
211            pdict['Files']['FreqDir'] = os.path.join ('Output' , 'MO' )
212        if pdict['Experiment']['Freq'] == "SE" :
213            pdict['Files']['FreqDir'] = os.path.join ('Analyse', 'SE' )
214   
215    if not pdict['Files']['dir_ATM_his'] :
216        pdict['Files']['dir_ATM_his'] = os.path.join ( pdict['libIGCM']['R_SAVE'], "ATM", pdict['Files']['FreqDir'] )
217    if not pdict['Files']['dir_SRF_his'] :
218        pdict['Files']['dir_SRF_his'] = os.path.join ( pdict['libIGCM']['R_SAVE'], "SRF", pdict['Files']['FreqDir'] )
219       
220    if not pdict['Files']['dir_OCE_his'] :
221        pdict['Files']['dir_OCE_his'] = os.path.join ( pdict['libIGCM']['R_SAVE'], "OCE", pdict['Files']['FreqDir'] )
222    if not pdict['Files']['dir_ICE_his'] :
223        pdict['Files']['dir_ICE_his'] = os.path.join ( pdict['libIGCM']['R_SAVE'], "ICE", pdict['Files']['FreqDir'] )
224       
225    echo (  'The analysis relies on files from the following model output directories : ', f_out )
226    echo ( f"{pdict['Files']['dir_ATM_his'] = }" , f_out )
227    echo ( f"{pdict['Files']['dir_SRF_his'] = }" , f_out )
228    echo ( f"{pdict['Files']['dir_OCE_his'] = }" , f_out )
229    echo ( f"{pdict['Files']['dir_ICE_his'] = }" , f_out )
230   
231    ## Creates files names
232    ## --------------------
233    if not pdict['Experiment']['Period'] :
234        if pdict['Experiment']['Freq'] == 'MO' :
235            pdict['Experiment']['Period'] = f"{pdict['Experiment']['DateBegin']}_{pdict['Experiment']['DateEnd']}_1M"
236        if pdict['Experiment']['Freq'] == 'SE' :
237            pdict['Experiment']['Period'] = f"SE_{pdict['Files']['DateBegin']}_{pdict['Files']['DateEnd']}_1M"
238   
239    echo ( f"Period   : {pdict['Experiment']['Period']}", f_out )
240   
241    if not pdict['Files']['FileCommon']  :
242        pdict['Files']['FileCommon'] = f"{pdict['Experiment']['JobName']}_{pdict['Experiment']['Period']}"
243       
244    if not pdict['Files']['Title'] :
245        pdict['Files']['Title'] = f"{pdict['Experiment']['JobName']} : {pdict['Experiment']['Freq']} : {pdict['Experiment']['DateBegin']} - {pdict['Experiment']['DateEnd']}"
246       
247    echo ('\nOpen history files', f_out )
248    if not pdict['Files']['file_ATM_his'] :
249        if pdict['Experiment']['ATM_HIS'] == 'latlon' :
250            pdict['Files']['file_ATM_his'] = os.path.join ( pdict['Files']['dir_ATM_his'], f"{pdict['Files']['FileCommon']}_histmth.nc" )
251        if pdict['Experiment']['ATM_HIS'] == 'ico' :
252            pdict['Files']['file_ATM_his'] = os.path.join ( pdict['Files']['dir_ATM_his'], f"{pdict['Files']['FileCommon']}_histmth_ico.nc" )
253    if not pdict['Files']['file_SRF_his'] :
254        if pdict['Experiment']['ATM_HIS'] == 'latlon' :
255            pdict['Files']['file_SRF_his'] = os.path.join ( pdict['Files']['dir_SRF_his'], f"{pdict['Files']['FileCommon']}_sechiba_history.nc" )
256        if pdict['Experiment']['ATM_HIS'] == 'ico' :
257            pdict['Files']['file_SRF_his'] = os.path.join ( pdict['Files']['dir_SRF_his'], f"{pdict['Files']['FileCommon']}_sechiba_history_ico.nc" )
258       
259    if pdict['Experiment']['SECHIBA'] and pdict['Experiment']['Routing'] == 'SIMPLE' :
260        if not pdict['Files']['file_RUN_his'] :
261            if pdict['Experiment']['RUN_HIS'] == 'latlon' :
262                pdict['Files']['file_RUN_his'] = os.path.join ( pdict['Files']['dir_SRF_his'], f"{pdict['Files']['FileCommon']}_sechiba_history.nc" )
263            if pdict['Experiment']['SRF_HIS'] == 'ico' :
264                pdict['Files']['file_RUN_his'] = os.path.join ( pdict['Files']['dir_SRF_his'], f"{pdict['Files']['FileCommon']}_sechiba_history_ico.nc" )
265           
266    echo ( f"{pdict['Files']['file_ATM_his'] = }", f_out )
267    if pdict['Experiment']['SECHIBA'] :
268        echo ( f"{pdict['Files']['file_SRF_his'] = }", f_out )
269        if pdict['Experiment']['Routing'] == 'SIMPLE' : echo ( f"{pdict['Files']['file_RUN_his'] = }", f_out )
270
271    if not pdict['Files']['file_OCE_his'] :
272        pdict['Files']['file_OCE_his'] = os.path.join ( pdict['Files']['dir_OCE_his'], f"{pdict['Files']['FileCommon']}_grid_T.nc" )
273    if not pdict['Files']['file_OCE_sca'] :
274        pdict['Files']['file_OCE_sca'] = os.path.join ( pdict['Files']['dir_OCE_his'], f"{pdict['Files']['FileCommon']}_scalar.nc" )
275    if not pdict['Files']['file_OCE_srf'] :
276        pdict['Files']['file_OCE_srf'] = os.path.join ( pdict['Files']['dir_OCE_his'], f"{pdict['Files']['FileCommon']}_sbc.nc" )
277    if not pdict['Files']['file_ICE_his'] :
278        pdict['Files']['file_ICE_his'] = os.path.join ( pdict['Files']['dir_ICE_his'], f"{pdict['Files']['FileCommon']}_icemod.nc" )
279   
280    return pdict
281
282def SetRestartNames ( pdict, f_out) :
283    '''
284    Defines dates for restart files
285    Define names of tar files with restart
286
287    Returns a dictionary
288    '''
289   
290    if not pdict['Files']['TarRestartDate_beg'] :
291        pdict['Files']['TarRestartDate_beg'] = ldate.SubOneDayToDate ( pdict['Experiment']['DateBegin'], pdict['Experiment']['CalendarType'] )
292    if not pdict['Files']['TarRestartDate_end'] :
293        pdict['Files']['TarRestartDate_end'] = ldate.ConvertFormatToGregorian ( pdict['Experiment']['DateEnd'] )
294   
295    if not pdict['Files']['TarRestartPeriod_beg'] :
296       
297        pdict['Files']['TarRestartPeriod_beg_DateEnd'] = pdict['Files']['TarRestartDate_beg']
298        pdict['Files']['TarRestartPeriod_beg_DateBeg'] = ldate.DateAddPeriod ( pdict['Files']['TarRestartPeriod_beg_DateEnd'],
299                                                                                   f"-{pdict['Experiment']['PackFrequency']}")
300        pdict['Files']['TarRestartPeriod_beg_DateBeg'] = ldate.AddOneDayToDate ( pdict['Files']['TarRestartPeriod_beg_DateBeg'],
301                                                                                     pdict['Experiment']['CalendarType']  )
302       
303        pdict['Files']['TarRestartPeriod_beg'] = f"{pdict['Files']['TarRestartPeriod_beg_DateBeg']}_{pdict['Files']['TarRestartPeriod_beg_DateEnd']}"
304        echo ( f"Tar period for initial restart : {pdict['Files']['TarRestartPeriod_beg']}", f_out)
305       
306    if not pdict['Files']['TarRestartPeriod_end'] :
307
308        pdict['Files']['TarRestartPeriod_end_DateEnd'] = pdict['Files']['TarRestartDate_end']
309        pdict['Files']['TarRestartPeriod_end_DateBeg'] = ldate.DateAddPeriod ( pdict['Files']['TarRestartPeriod_end_DateEnd'],
310                                                                                   f"-{pdict['Experiment']['PackFrequency']}" )
311        pdict['Files']['TarRestartPeriod_end_DateBeg'] = ldate.AddOneDayToDate ( pdict['Files']['TarRestartPeriod_end_DateBeg'],
312                                                                                     pdict['Experiment']['CalendarType'] )
313       
314        pdict['Files']['TarRestartPeriod_end'] = f"{pdict['Files']['TarRestartPeriod_end_DateBeg']}_{pdict['Files']['TarRestartPeriod_end_DateEnd']}"
315        echo ( f"Tar period for final restart : {pdict['Files']['TarRestartPeriod_end']}", f_out )
316       
317    echo ( f"Restart dates - Start : {pdict['Files']['TarRestartPeriod_beg']}  /  End : {pdict['Files']['TarRestartPeriod_end']}", f_out)
318       
319    if not pdict['Files']['tar_restart_beg'] :
320        pdict['Files']['tar_restart_beg'] = os.path.join ( pdict['libIGCM']['R_SAVE'], 'RESTART', f"{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartPeriod_beg']}_restart.tar" )
321    if not  pdict['Files']['tar_restart_end'] :
322        pdict['Files']['tar_restart_end'] = os.path.join ( pdict['libIGCM']['R_SAVE'], 'RESTART', f"{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartPeriod_end']}_restart.tar" )
323       
324    echo ( f"{ pdict['Files']['tar_restart_beg'] = }", f_out )
325    echo ( f"{ pdict['Files']['tar_restart_end'] = }", f_out )
326
327    ##-- Names of tar files with restarts
328   
329    if not pdict['Files']['tar_restart_beg_ATM'] :
330        pdict['Files']['tar_restart_beg_ATM'] = pdict['Files']['tar_restart_beg']
331    if not pdict['Files']['tar_restart_beg_DYN'] :
332        pdict['Files']['tar_restart_beg_DYN'] = pdict['Files']['tar_restart_beg']
333    if not pdict['Files']['tar_restart_beg_RUN'] :
334        pdict['Files']['tar_restart_beg_RUN'] = pdict['Files']['tar_restart_beg']
335    if not pdict['Files']['tar_restart_beg_OCE'] :
336        pdict['Files']['tar_restart_beg_OCE'] = pdict['Files']['tar_restart_beg']
337    if not pdict['Files']['tar_restart_beg_ICE'] :
338        pdict['Files']['tar_restart_beg_ICE'] = pdict['Files']['tar_restart_beg']
339       
340    if not pdict['Files']['tar_restart_end_ATM'] :
341        pdict['Files']['tar_restart_end_ATM'] = pdict['Files']['tar_restart_end']
342    if not pdict['Files']['tar_restart_end_DYN'] :
343        pdict['Files']['tar_restart_end_DYN'] = pdict['Files']['tar_restart_end']
344    if not pdict['Files']['tar_restart_end_RUN'] :
345        pdict['Files']['tar_restart_end_RUN'] = pdict['Files']['tar_restart_end']
346    if not pdict['Files']['tar_restart_end_OCE'] :
347        pdict['Files']['tar_restart_end_OCE'] = pdict['Files']['tar_restart_end']
348    if not pdict['Files']['tar_restart_end_ICE'] :
349        pdict['Files']['tar_restart_end_ICE'] = pdict['Files']['tar_restart_end']
350       
351    if not pdict['Files']['file_DYN_beg'] :
352        if pdict['Experiment']['LMDZ'] :
353            pdict['Files']['file_DYN_beg'] = f"{pdict['Files']['FileDir']}/ATM_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_beg']}_restart.nc"
354        if pdict['Experiment']['ICO']  :
355            pdict['Files']['file_DYN_beg'] = f"{pdict['Files']['FileDir']}/ICO_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_beg']}_restart.nc"
356       
357    if not pdict['Files']['file_DYN_end'] :
358        if pdict['Experiment']['LMDZ'] :
359            pdict['Files']['file_DYN_end'] = f"{pdict['Files']['FileDir']}/ATM_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_end']}_restart.nc"
360        if pdict['Experiment']['ICO']  :
361            pdict['Files']['file_DYN_end'] = f"{pdict['Files']['FileDir']}/ICO_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_end']}_restart.nc"
362       
363    if pdict['Experiment']['SECHIBA'] :
364        if not pdict['Files']['file_SRF_beg'] :
365            pdict['Files']['file_SRF_beg'] = f"{pdict['Files']['FileDir']}/SRF_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_beg']}_sechiba_rest.nc"
366        if not pdict['Files']['file_SRF_end'] :
367            pdict['Files']['file_SRF_end'] = f"{pdict['Files']['FileDir']}/SRF_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_end']}_sechiba_rest.nc"
368       
369    if pdict['Experiment']['SECHIBA'] :
370        if not pdict['Files']['tar_restart_beg_SRF'] :
371            pdict['Files']['tar_restart_beg_SRF'] = pdict['Files']['tar_restart_beg']
372        if not pdict['Files']['tar_restart_end_SRF'] :
373            pdict['Files']['tar_restart_end_SRF'] = pdict['Files']['tar_restart_end']
374           
375    if not pdict['Files']['file_ATM_beg'] :
376        pdict['Files']['file_ATM_beg'] = f"{pdict['Files']['FileDir']}/ATM_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_beg']}_restartphy.nc"
377    if not pdict['Files']['file_ATM_end'] :
378        pdict['Files']['file_ATM_end'] = f"{pdict['Files']['FileDir']}/ATM_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_end']}_restartphy.nc"
379       
380    if pdict['Experiment']['SECHIBA'] :
381        echo ( f"{pdict['Files']['file_SRF_beg'] = }", f_out)
382        echo ( f"{pdict['Files']['file_SRF_end'] = }", f_out)
383       
384    if pdict['Experiment']['ICO'] :
385        if not pdict['Files']['file_DYN_aire'] :
386            pdict['Files']['file_DYN_aire'] = os.path.join ( pdict['libIGCM']['R_IN'], 'ATM', 'GRID',  f"{pdict['Experiment']['ResolAtm']}_grid.nc" )
387       
388    if pdict['Experiment']['SECHIBA'] and pdict['Experiment']['Routing'] == 'SIMPLE' :
389        if not pdict['Files']['file_RUN_beg'] :
390            pdict['Files']['file_RUN_beg'] = f"{ pdict['Files']['FileDir']}/SRF_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_beg']}_routing_restart.nc"
391        if not pdict['Files']['file_RUN_end'] :
392            pdict['Files']['file_RUN_end'] = f"{ pdict['Files']['FileDir']}/SRF_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_end']}_routing_restart.nc"
393
394    if not pdict['Files']['tar_restart_beg_OCE'] :
395        pdict['Files']['tar_restart_beg_OCE'] = pdict['Files']['tar_restart_beg']
396    if not pdict['Files']['tar_restart_beg_ICE'] :
397        pdict['Files']['tar_restart_beg_ICE'] = pdict['Files']['tar_restart_beg']
398       
399    if not pdict['Files']['tar_restart_end_OCE'] :
400        pdict['Files']['tar_restart_end_OCE'] = pdict['Files']['tar_restart_end']
401    if not pdict['Files']['tar_restart_end_ICE'] :
402        pdict['Files']['tar_restart_end_ICE'] = pdict['Files']['tar_restart_end']
403       
404    if not pdict['Files']['file_OCE_beg'] :
405        pdict['Files']['file_OCE_beg'] = f"{pdict['Files']['FileDir']}/OCE_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_beg']}_restart.nc"
406    if not pdict['Files']['file_OCE_end'] :
407        pdict['Files']['file_OCE_end'] = f"{pdict['Files']['FileDir']}/OCE_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_end']}_restart.nc"
408    if not pdict['Files']['file_ICE_beg'] :
409        pdict['Files']['file_ICE_beg'] = f"{pdict['Files']['FileDir']}/ICE_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_beg']}_restart_icemod.nc"
410    if not pdict['Files']['file_ICE_end'] :
411        pdict['Files']['file_ICE_end'] = f"{pdict['Files']['FileDir']}/ICE_{pdict['Experiment']['JobName']}_{pdict['Files']['TarRestartDate_end']}_restart_icemod.nc"
412           
413    return pdict
414
415def ComputeGridATM ( dpar, d_ATM_his ) :
416    readPrec = dpar['Config']['readPrec']
417    if repr(readPrec) in [ "<class 'numpy.float64'>" , float ] :
418        def rprec (ptab) :
419            return ptab
420    else                 :
421        def rprec (ptab) :
422            return ptab.astype(readPrec).astype(float)
423   
424    ATM = types.SimpleNamespace ()
425
426    # ATM grid with cell surfaces
427    if dpar['Experiment']['LMDZ'] :
428        #echo ('ATM grid with cell surfaces : LMDZ', f_out)
429        ATM.lat       = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1d='cell' )
430        ATM.lon       = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1d='cell' )
431        ATM.aire      = lmdz.geo2point ( rprec (d_ATM_his ['aire']     [0]), cumul_poles=True, dim1d='cell' )
432        ATM.fter      = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1d='cell' )
433        ATM.foce      = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1d='cell' )
434        ATM.fsic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1d='cell' )
435        ATM.flic      = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1d='cell' )
436       
437    if  dpar['Experiment']['ICO'] :
438        if  dpar['Experiment']['ATM_HIS'] == 'latlon' :
439            #echo ( 'ATM areas and fractions on LATLON grid' )
440            if 'lat_dom_out' in d_ATM_his.variables :
441                ATM.lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat_dom_out'])+0*rprec (d_ATM_his ['lon_dom_out']), dim1d='cell' )
442                ATM.lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat_dom_out'])+  rprec (d_ATM_his ['lon_dom_out']), dim1d='cell' )
443            else :
444                ATM.lat  = lmdz.geo2point (   rprec (d_ATM_his ['lat'])+0*rprec (d_ATM_his ['lon']), dim1d='cell' )
445                ATM.lon  = lmdz.geo2point ( 0*rprec (d_ATM_his ['lat'])+  rprec (d_ATM_his ['lon']), dim1d='cell' )
446            ATM.aire = lmdz.geo2point ( rprec (d_ATM_his ['aire'][0]).squeeze(), cumul_poles=True, dim1d='cell' )
447            ATM.fter = lmdz.geo2point ( rprec (d_ATM_his ['fract_ter'][0]), dim1d='cell' )
448            ATM.foce = lmdz.geo2point ( rprec (d_ATM_his ['fract_oce'][0]), dim1d='cell' )
449            ATM.fsic = lmdz.geo2point ( rprec (d_ATM_his ['fract_sic'][0]), dim1d='cell' )
450            ATM.flic = lmdz.geo2point ( rprec (d_ATM_his ['fract_lic'][0]), dim1d='cell' )
451
452           
453        if  dpar['Experiment']['ATM_HIS'] == 'ico' :
454            #echo ( 'ATM areas and fractions on ICO grid' )
455            ATM.aire =  rprec (d_ATM_his ['aire']     [0]).squeeze()
456            ATM.lat  =  rprec (d_ATM_his ['lat']         )
457            ATM.lon  =  rprec (d_ATM_his ['lon']         )
458            ATM.fter =  rprec (d_ATM_his ['fract_ter'][0])
459            ATM.foce =  rprec (d_ATM_his ['fract_oce'][0])
460            ATM.fsic =  rprec (d_ATM_his ['fract_sic'][0])
461            ATM.flic =  rprec (d_ATM_his ['fract_lic'][0])
462           
463       
464    ATM.fsea      = ATM.foce + ATM.fsic
465    ATM.flnd      = ATM.fter + ATM.flic
466    ATM.aire_fter = ATM.aire * ATM.fter
467    ATM.aire_flic = ATM.aire * ATM.flic
468    ATM.aire_fsic = ATM.aire * ATM.fsic
469    ATM.aire_foce = ATM.aire * ATM.foce
470    ATM.aire_flnd = ATM.aire * ATM.flnd
471    ATM.aire_fsea = ATM.aire * ATM.fsea
472
473    ATM.aire_sea     = ATM.aire * ATM.fsea
474    ATM.aire_tot     = P1sum ( ATM.aire      )
475    ATM.aire_sea_tot = P1sum ( ATM.aire_fsea )
476    ATM.aire_ter_tot = P1sum ( ATM.aire_fter )
477    ATM.aire_lic_tot = P1sum ( ATM.aire_flic )
478   
479    return dpar, ATM
480
481def ComputeGridSRF ( dpar, d_SRF_his ) :
482    readPrec = dpar['Config']['readPrec']
483    if repr(readPrec) in [ "<class 'numpy.float64'>" , float ] :
484        def rprec (ptab) :
485            return ptab
486    else                 :
487        def rprec (ptab) :
488            return ptab.astype(readPrec).astype(float)
489   
490    SRF = types.SimpleNamespace ()
491   
492    if dpar['Experiment']['SECHIBA'] : 
493        if dpar['Experiment']['LMDZ'] :
494            SRF.lat       = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1d='cell' )
495            SRF.lon       = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1d='cell' )
496            SRF.aire      = lmdz.geo2point ( rprec (d_SRF_his ['Areas']) * rprec (d_SRF_his ['Contfrac']), dim1d='cell', cumul_poles=True )
497            SRF.areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas'])  ,  dim1d='cell', cumul_poles=True )
498            SRF.contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac']), dim1d='cell' )
499   
500        if dpar['Experiment']['ICO'] :
501            if dpar['Experiment']['SRF_HIS'] == 'latlon' :
502                #echo ( 'SRF areas and fractions on LATLON grid' )
503                if 'lat_domain_landpoints_out' in d_SRF_his  :
504                    SRF.lat = lmdz.geo2point (   rprec (d_SRF_his ['lat_domain_landpoints_out'])+0*rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1d='cell' )
505                    SRF.lon = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_domain_landpoints_out'])+  rprec (d_SRF_his ['lon_domain_landpoints_out']), dim1d='cell' )
506                else :
507                    if 'lat_domain_landpoints_out' in d_SRF_his :
508                        SRF.lat = lmdz.geo2point (   rprec (d_SRF_his ['lat_dom_out'])+0*rprec (d_SRF_his ['lon_dom_out']), dim1d='cell' )
509                        SRF.lon = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat_dom_out'])+  rprec (d_SRF_his ['lon_dom_out']), dim1d='cell' )
510                    else :
511                        SRF.lat = lmdz.geo2point (   rprec (d_SRF_his ['lat'])+0*rprec (d_SRF_his ['lon']), dim1d='cell' )
512                        SRF.lon = lmdz.geo2point ( 0*rprec (d_SRF_his ['lat'])+  rprec (d_SRF_his ['lon']), dim1d='cell' )
513   
514                SRF.areas     = lmdz.geo2point ( rprec (d_SRF_his ['Areas']   )   , dim1d='cell', cumul_poles=True )
515                SRF.areafrac  = lmdz.geo2point ( rprec (d_SRF_his ['AreaFrac'])   , dim1d='cell', cumul_poles=True )
516                SRF.contfrac  = lmdz.geo2point ( rprec (d_SRF_his ['Contfrac'])   , dim1d='cell', cumul_poles=True )
517                SRF.aire      = SRF.areafrac
518   
519            if dpar['Experiment']['SRF_HIS'] == 'ico' :
520                #echo ( 'SRF areas and fractions on ICO grid' )
521                SRF.lat       =  rprec (d_SRF_his ['lat']     )
522                SRF.lon       =  rprec (d_SRF_his ['lon']     )
523                SRF.areas     =  rprec (d_SRF_his ['Areas']   )
524                SRF.contfrac  =  rprec (d_SRF_his ['Contfrac'])
525                SRF.aire      =  SRF.areas * SRF.contfrac
526
527        SRF.aire_tot = P1sum ( SRF.aire )
528
529    return dpar, SRF
530
531def ComputeGridDYN ( dpar, ATM, d_DYN_aire, d_ATM_beg ) : 
532    readPrec = dpar['Config']['readPrec']
533    if repr(readPrec) in [ "<class 'numpy.float64'>" , float ] :
534        def rprec (ptab) :
535            return ptab
536    else                 :
537        def rprec (ptab) :
538            return ptab.astype(readPrec).astype(float)
539       
540    DYN = types.SimpleNamespace ()
541    if dpar['Experiment']['ICO'] :
542        if dpar['Config']['SortIco'] :
543            # Creation d'une clef de tri pour le fichier aire
544            DYN.aire_keysort = np.lexsort ( (d_DYN_aire['lat'], d_DYN_aire['lon']) )
545        else :
546            DYN.aire_keysort = np.arange ( len ( d_DYN_aire['lat'] ) )
547           
548        DYN.lat = d_DYN_aire['lat']
549        DYN.lon = d_DYN_aire['lon']
550       
551        DYN.aire = d_DYN_aire['aire']
552        DYN.fsea = d_DYN_aire['fract_oce'] + d_DYN_aire['fract_sic']
553       
554        DYN.flnd = 1.0 - DYN.fsea
555        DYN.fter = d_ATM_beg['FTER']
556        DYN.flic = d_ATM_beg['FLIC']
557        DYN.foce = d_ATM_beg['FOCE']
558        DYN.aire_fter = DYN.aire * DYN.fter
559       
560    if dpar['Experiment']['ATM_HIS'] == 'ico' :
561        DYN.aire      = ATM.aire
562        DYN.foce      = ATM.foce
563        DYN.fsic      = ATM.fsic
564        DYN.flic      = ATM.flic
565        DYN.fter      = ATM.fter
566        DYN.fsea      = ATM.fsea
567        DYN.flnd      = ATM.flnd
568        DYN.aire_fter = ATM.aire_fter
569        DYN.aire_flic = ATM.aire_flic
570        DYN.aire_fsic = ATM.aire_fsic
571        DYN.aire_foce = ATM.aire_foce
572        DYN.aire_flnd = ATM.aire_flnd
573        DYN.aire_fsea = ATM.aire_fsea
574   
575    if dpar['Experiment']['LMDZ'] :
576        # Area on lon/lat grid
577        DYN.aire = ATM.aire
578        DYN.fsea = ATM.fsea
579        DYN.flnd = ATM.flnd
580        DYN.fter = rprec (d_ATM_beg['FTER'])
581        DYN.flic = rprec (d_ATM_beg['FLIC'])
582        DYN.aire_fter = DYN.aire * DYN.fter
583       
584    DYN.aire_tot =  P1sum ( DYN.aire ) 
585   
586    return dpar, DYN
587
588def ComputeGridOCE ( dpar, d_OCE_his, nperio=None ) :
589    OCE = types.SimpleNamespace ()
590
591    # Get mask and surfaces
592    sos = d_OCE_his ['sos'][0].squeeze()
593    OCE.OCE_msk = nemo.lbc_mask ( xr.where ( sos>0., 1., 0. ), cd_type = 'T', nperio=nperio, sval = 0. )
594   
595    so = sos = d_OCE_his ['sos'][0].squeeze()
596    OCE.msk3 = nemo.lbc_mask ( xr.where ( so>0., 1., 0. ), cd_type = 'T', nperio=nperio, sval = 0. )
597   
598    # lbc_mask removes the duplicate points (periodicity and north fold)
599    OCE.OCE_aire = nemo.lbc_mask ( d_OCE_his ['area'] * OCE.OCE_msk, cd_type = 'T', nperio=nperio, sval = 0.0 )
600    OCE.ICE_aire = OCE.OCE_aire
601
602    OCE.OCE_aire_tot = P1sum ( OCE.OCE_aire )
603    OCE.ICE_aire_tot = P1sum ( OCE.ICE_aire )
604
605    return dpar, OCE
606
607def config2dict ( pconf ) :
608    '''
609    Convert a config parser object into a dictionary
610    '''
611    pdict = {}
612    for section in pconf.keys () :
613        pdict[section] = {}
614        for key in pconf[section].keys() :
615            zz = str2value ( pconf[section][key] )
616            pdict[section].update ( {key:zz} )       
617    return pdict
618
619def dict2config ( pdict ):
620    '''
621    Convert a dictionary into a configparser object
622    All values are converted to strings
623
624    The dictionary must have two levels (see configparser)
625    '''
626    zconf = configparser.ConfigParser ( interpolation=configparser.ExtendedInterpolation() )
627    zconf.optionxform = str # To keep capitals
628    zconf.read_dict(dict2str(pdict))
629   
630    return zconf
631
632def updateConf ( pconf, pdict ):
633    '''
634    Update a config parser with a dictionary
635    All values are converted to strings
636
637    The dictionary must have two levels (see configparser)
638    '''
639    pconf.read_dict (dict1str(pdict))
640    #for section in pdict.keys() :
641    #    if section not in pconf.keys() : pconf[section] = {}
642    #    for key in pdict[section].keys() :
643    #        pconf[section][key] = str(pdict[section][key])
644    return pconf
645
646def dict2str (pdict) :
647    '''
648    Convert all values of a dictionary to strings
649    '''
650    zout = {}
651    for section in pdict.keys() :
652        zout[section] = {}
653        for key in pdict[section].keys() :
654            zout[section][key] = str(pdict[section][key])
655    return zout
656
657def echo (string, f_out, end='\n') :
658    '''Function to print to stdout *and* output file'''
659    print ( str(string), end=end  )
660    sys.stdout.flush ()
661    f_out.write ( str(string) + end )
662    f_out.flush ()
663
664def str2value ( pz ) :
665    '''
666    Tries to convert a string into integer, real, boolean or None
667    '''
668    zz = pz
669    zz = setBool (zz)
670    zz = setNum  (zz)
671    zz = setNone (zz)
672    return zz
673
674def setBool (chars) :
675    '''Convert specific char string in boolean if possible'''
676    z_set_bool = chars
677    if isinstance (chars, str) :
678        for key in configparser.ConfigParser.BOOLEAN_STATES.keys () :
679            if chars.lower() == key :
680                z_set_bool = configparser.ConfigParser.BOOLEAN_STATES[key]
681    return z_set_bool
682
683def setNum (chars) :
684    '''Convert specific char string in integer or real if possible'''
685    zset_num = chars
686    if isinstance (chars, str) :
687        realnum   = re.compile ("^[-+]?[0-9]*\.?[0-9]+(e[-+]?[0-9]+)?$")
688        is_number = realnum.match(chars.strip()) != None
689        is_int    = chars.strip().isdigit()
690        if is_number :
691            if is_int : zset_num = int   (chars)
692            else      : zset_num = float (chars)
693   
694    return zset_num
695
696def setNone (chars) :
697    '''Convert specific char string to None if possible'''
698    zset_none = chars
699    if isinstance (chars, str) :
700        if chars in ['None', 'NONE', 'none'] : zset_none = None
701    return zset_none
702
703def unDefined (char) :
704    '''True if a variable is not defined, or set to None'''
705    if char in globals () :
706        if globals()[char] :
707            zun_defined = False
708        else : zun_defined = True
709    else :     zun_defined = True
710    return zun_defined
711
712def Defined (char) :
713    '''True if a variable is defined and not equal to None'''
714    if char in globals () :
715        if globals()[char] is None :
716            zdefined = False
717        else            : zdefined = True
718    else                : zdefined = False
719    return zdefined
720
721def Ksum (tab) :
722    '''Kahan summation algorithm, also known as compensated summation
723
724    https://en.wikipedia.org/wiki/Kahan_summation_algorithm
725    '''
726    # Prepare the accumulator.
727    ksum = 0.0
728    # A running compensation for lost low-order bits.
729    comp = 0.0
730
731    for xx in np.where ( np.isnan(tab), 0., tab ) :
732        # comp is zero the first time around.
733        yy = xx - comp
734        # Alas, sum is big, y small, so low-order digits of y are lost.
735        tt = ksum + yy
736        # (tt - Ksum) cancels the high-order part of y; subtracting y recovers negative (low part of yy)
737        comp = (tt - ksum) - yy
738        # Algebraically, comp should always be zero. Beware overly-aggressive optimizing compilers !
739        ksum = tt
740        # Next time around, the lost low part will be added to y in a fresh attempt.
741
742    return ksum
743
744def Ssum (tab) :
745    '''Precision summation by sorting by absolute values'''
746    ssum = np.sum ( tab[np.argsort(np.abs(tab))] )
747    return ssum
748
749def KSsum (tab) :
750    '''Precision summation by sorting by absolute value, and applying Kahan summation'''
751    kssum = Ksum ( tab[np.argsort(np.abs(tab))] )
752    return kssum
753
754# Choosing algorithm for precise summation
755Psum = KSsum
756
757def P1sum (ptab) :
758    return Psum ( ptab.to_masked_array().ravel() )
759
760class Timer :
761    '''Measures execution time of a function'''
762    def __str__ (self):
763        return str (self.__class__)
764   
765    def __name__ (self):
766        return self.__class__.__name__
767
768    def __init__ (self, func, debug=False, timer=True) :
769        functools.update_wrapper (self, func)
770        self.func       = func
771        self.num_calls  = 0
772        self.cumul_time = 0.
773        self.debug      = debug
774        self.timer      = timer
775
776    def __call__ (self, *args, **kwargs) :
777        self.num_calls += 1
778        if self.debug :
779            print ( f'>-- Calling {self.__name__} --------------------------------------------------' )
780            args_repr   = [f"{repr (a)}" for a in args]
781            kwargs_repr = [f"{k}={v!r}" for k, v in kwargs.items ()]
782            signature   = ", ".join (args_repr + kwargs_repr)
783            print ( f"Signature : ({signature}")
784        start_time = time.perf_counter ()     
785        values     = self.func (*args, **kwargs)
786        end_time   = time.perf_counter ()
787        run_time   = end_time - start_time
788        self.cumul_time += run_time
789        if self.timer : 
790            print ( f"--> Calling {self.__name__!r} : {run_time:.3f} secs / cumul {self.cumul_time:.3f} # {self.num_calls:2d} calls")
791        if self.debug :
792            print ( '<------------------------------------------------------------' )
793        return values
794   
Note: See TracBrowser for help on using the repository browser.