Ignore:
Timestamp:
10/10/23 12:58:04 (7 months ago)
Author:
omamce
Message:

O.M. :

New version of WATER_BUDGET

  • Conservation in NEMO are OK
  • Conservation in Sechiba are OK, for both ICO and latlon grids
  • Problems in atmosphere, LIC surface and ocean/atmosphere coherence
File:
1 edited

Legend:

Unmodified
Added
Removed

  • TOOLS/WATER_BUDGET/OCE_waterbudget.py

    r6620 r6647  
    3939import nemo, lmdz 
    4040 
     41from 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 
     42 
    4143## Creates parser for reading .ini input file 
    4244## ------------------------------------------- 
    … …  
    4749## --------------------- 
    4850ATM=None ; ATM_HIS='latlon' ; SRF_HIS='latlon' ; RUN_HIS='latlon' ; ORCA=None ; NEMO=None ; OCE_relax=False 
    49 OCE_icb=False ; Coupled=False ; Routing=None 
     51OCE_icb=False ; Coupled=False ; Routing=None ; TestInterp=None 
    5052TarRestartPeriod_beg=None ; TarRestartPeriod_end=None ; Comment=None ; Period=None ; Title=None 
     53YearBegin=None ; YearEnd=None ; DateBegin=None ; DateEnd=None 
    5154 
    5255## 
    … …  
    6669tar_restart_end_RUN=None ; tar_restart_end_OCE=None ; tar_restart_end_ICE=None 
    6770ContinueOnError=False ; ErrorCount=0 ; SortIco = False 
     71 
     72## 
     73## Precision of history file reading 
     74## --------------------------------- 
     75# Default is float (full precision). Degrade the precision by using np.float32 
     76# Restart file are always read at the full precision 
     77readPrec=float 
    6878 
    6979## Read command line arguments 
    … …  
    137147 
    138148print ( f'\nConfig file readed : {IniFile} ' ) 
    139              
     149 
     150## 
     151## Reading prec 
     152if wu.unDefined ( 'readPrec' )  : 
     153    readPrec = np.float64 
     154else : 
     155    if readPrec in ["float", "float64", "r8", "double", "<class 'float'>"         ] : readPrec = float 
     156    if readPrec in [         "float32", "r4", "single", "<class 'numpy.float32'>" ] : readPrec = np.float32 
     157    if readPrec in [         "float16", "r2", "half"  , "<class 'numpy.float16'>" ] : readPrec = np.float16 
     158               
    140159## Some physical constants 
    141160## ======================= 
    … …  
    155174if not 'Files' in config.keys () : config['Files'] = {} 
    156175 
    157 config['Physics'] = { 'Ra':Ra, 'Grav':Grav, 'ICE_rho_ice':ICE_rho_ice, 'ICE_rho_sno':ICE_rho_sno, 
    158                       'OCE_rho_liq':OCE_rho_liq, 'ATM_rho':ATM_rho, 'SRF_rho':SRF_rho, 'RUN_rho':RUN_rho} 
    159  
    160 config['Config'] = { 'ContinueOnError':ContinueOnError, 'SortIco':SortIco } 
     176config['Physics'] = { 'Ra':str(Ra), 'Grav':str(Grav), 'ICE_rho_ice':str(ICE_rho_ice), 'ICE_rho_sno':str(ICE_rho_sno), 
     177                      'OCE_rho_liq':str(OCE_rho_liq), 'ATM_rho':str(ATM_rho), 'SRF_rho':str(SRF_rho), 'RUN_rho':str(RUN_rho)} 
     178 
     179config['Config'] = { 'ContinueOnError':ContinueOnError, 'SortIco':SortIco, 'TestInterp':str(TestInterp), 'readPrec':str(readPrec) } 
    161180 
    162181## -------------------------- 
    … …  
    178197    config['Experiment']['DateBegin'] = DateBegin 
    179198else : 
    180     YearBegin, MonthBegin, DayBegin = wu.DateSplit ( DateBegin ) 
     199    YearBegin, MonthBegin, DayBegin = wu.SplitDate ( DateBegin ) 
    181200    DateBegin = wu.FormatToGregorian ( DateBegin) 
    182     config['Experiment']['YearBegin'] = YearBegin 
     201    config['Experiment']['YearBegin'] = str(YearBegin) 
    183202 
    184203if wu.unDefined ( 'DateEnd' )   : 
    185204    DateEnd   = f'{YearEnd}1231' 
    186     config['Experiment']['DateEnd']   = DateEnd 
     205    config['Experiment']['DateEnd']   = str(DateEnd) 
    187206else : 
    188     YearEnd, MonthEnd, DayEnd = wu.DateSplit ( DateEnd ) 
     207    YearEnd, MonthEnd, DayEnd = wu.SplitDate ( DateEnd ) 
    189208    DateEnd   = wu.FormatToGregorian ( DateEnd) 
     209    config['Experiment']['DateEnd']   = str(DateEnd) 
    190210 
    191211if wu.unDefined ( 'PackFrequency' ) : 
    … …  
    204224## Output file with water budget diagnostics 
    205225## ----------------------------------------- 
    206 if wu.unDefined ( 'FileOut' ) : FileOut = f'ATM_waterbudget_{JobName}_{DateBegin}_{DateEnd}.out' 
    207 config['Files']['FileOut'] = FileOut 
     226if wu.unDefined ( 'FileOut' ) : 
     227    FileOut = f'OCE_waterbudget_{JobName}_{YearBegin}_{YearEnd}' 
     228    if readPrec == np.float32 : FileOut = f'{FileOut}_float32' 
     229    FileOut = f'{FileOut}.out' 
     230 
     231    config['Files']['FileOut'] = FileOut 
    208232 
    209233f_out = open ( FileOut, mode = 'w' ) 
    … …  
    211235## Useful functions 
    212236## ---------------- 
    213 def kg2Sv  (val, rho=OCE_rho_liq) : 
     237 
     238# Degrades precision  
     239if readPrec == float : 
     240    def rprec (tab) : return tab 
     241else :  
     242    def rprec (tab) : return tab.astype(readPrec).astype(float) 
     243         
     244def kg2Sv  (val, rho=ATM_rho) : 
    214245    '''From kg to Sverdrup''' 
    215246    return val/dtime_sec*1.0e-6/rho 
    216247 
    217 def kg2myear (val, rho=OCE_rho_liq) : 
     248def kg2myear (val, rho=ATM_rho) : 
    218249    '''From kg to m/year''' 
    219250    return val/OCE_aire_tot/rho/NbYear 
    220251 
    221 def var2prt (var, small=False, rho=OCE_rho_liq) : 
     252def var2prt (var, small=False, rho=ATM_rho) : 
    222253    if small :  return var , kg2Sv(var, rho=rho)*1000., kg2myear(var, rho=rho)*1000. 
    223254    else     :  return var , kg2Sv(var, rho=rho)      , kg2myear(var, rho=rho) 
    224255 
    225 def prtFlux (Desc, var, Form='F', small=False, rho=OCE_rho_liq, width=15) : 
     256def prtFlux (Desc, var, Form='F', small=False, rho=ATM_rho, width=15) : 
    226257    if small : 
    227258        if Form in ['f', 'F'] : ff=" {:14.6e} kg | {:12.4f} mSv | {:12.4f} mm/year " 
    … …  
    330361if NEMO == 3.6 :d_ICE_his = d_ICE_his.rename ( {'y_grid_T':'y', 'x_grid_T':'x'} ) 
    331362 
    332 echo ( file_OCE_his ) 
    333 echo ( file_OCE_sca ) 
    334 echo ( file_ICE_his ) 
     363echo ( f'{file_OCE_his = }' ) 
     364echo ( f'{file_OCE_sca = }' ) 
     365echo ( f'{file_ICE_his = }' ) 
    335366 
    336367## Compute run length 
    … …  
    666697 
    667698# Read variable and computes integral over space and time 
    668 OCE_empmr      = d_OCE_his['wfo']         ; OCE_mas_empmr    = OCE_flux_int ( OCE_empmr    ) 
    669 OCE_wfob       = d_OCE_his['wfob']        ; OCE_mas_wfob     = OCE_flux_int ( OCE_wfob     ) 
    670 OCE_emp_oce    = d_OCE_his['emp_oce']     ; OCE_mas_emp_oce  = OCE_flux_int ( OCE_emp_oce  ) 
    671 OCE_emp_ice    = d_OCE_his['emp_ice']     ; OCE_mas_emp_ice  = OCE_flux_int ( OCE_emp_ice  ) 
    672 OCE_iceshelf   = d_OCE_his['iceshelf']    ; OCE_mas_iceshelf = OCE_flux_int ( OCE_iceshelf ) 
    673 OCE_calving    = d_OCE_his['calving']     ; OCE_mas_calving  = OCE_flux_int ( OCE_calving  ) 
    674 OCE_iceberg    = d_OCE_his['iceberg']     ; OCE_mas_iceberg  = OCE_flux_int ( OCE_iceberg  ) 
    675 OCE_friver     = d_OCE_his['friver']      ; OCE_mas_friver   = OCE_flux_int ( OCE_friver   ) 
    676 OCE_runoffs    = d_OCE_his['runoffs']     ; OCE_mas_runoffs  = OCE_flux_int ( OCE_runoffs  ) 
     699OCE_empmr      = rprec (d_OCE_his['wfo']     )    ; OCE_mas_empmr    = OCE_flux_int ( OCE_empmr    ) 
     700OCE_wfob       = rprec (d_OCE_his['wfob']    )    ; OCE_mas_wfob     = OCE_flux_int ( OCE_wfob     ) 
     701OCE_emp_oce    = rprec (d_OCE_his['emp_oce'] )    ; OCE_mas_emp_oce  = OCE_flux_int ( OCE_emp_oce  ) 
     702OCE_emp_ice    = rprec (d_OCE_his['emp_ice'] )    ; OCE_mas_emp_ice  = OCE_flux_int ( OCE_emp_ice  ) 
     703OCE_iceshelf   = rprec (d_OCE_his['iceshelf'])    ; OCE_mas_iceshelf = OCE_flux_int ( OCE_iceshelf ) 
     704OCE_calving    = rprec (d_OCE_his['calving'] )    ; OCE_mas_calving  = OCE_flux_int ( OCE_calving  ) 
     705OCE_iceberg    = rprec (d_OCE_his['iceberg'] )    ; OCE_mas_iceberg  = OCE_flux_int ( OCE_iceberg  ) 
     706OCE_friver     = rprec (d_OCE_his['friver']  )    ; OCE_mas_friver   = OCE_flux_int ( OCE_friver   ) 
     707OCE_runoffs    = rprec (d_OCE_his['runoffs'] )    ; OCE_mas_runoffs  = OCE_flux_int ( OCE_runoffs  ) 
    677708if NEMO == 4.0 or NEMO == 4.2 : 
    678     OCE_wfxice     = d_OCE_his['vfxice'] ; OCE_mas_wfxice   = OCE_flux_int ( OCE_wfxice ) 
    679     OCE_wfxsnw     = d_OCE_his['vfxsnw'] ; OCE_mas_wfxsnw   = OCE_flux_int ( OCE_wfxsnw ) 
    680     OCE_wfxsub     = d_OCE_his['vfxsub'] ; OCE_mas_wfxsub   = OCE_flux_int ( OCE_wfxsub ) 
     709    OCE_wfxice     = rprec (d_OCE_his['vfxice']) ; OCE_mas_wfxice   = OCE_flux_int ( OCE_wfxice ) 
     710    OCE_wfxsnw     = rprec (d_OCE_his['vfxsnw']) ; OCE_mas_wfxsnw   = OCE_flux_int ( OCE_wfxsnw ) 
     711    OCE_wfxsub     = rprec (d_OCE_his['vfxsub']) ; OCE_mas_wfxsub   = OCE_flux_int ( OCE_wfxsub ) 
    681712if NEMO == 3.6 : 
    682     OCE_wfxice     = d_OCE_his['vfxice']/86400.*ICE_rho_ice ; OCE_mas_wfxice   = OCE_flux_int ( OCE_wfxice ) 
    683     OCE_wfxsnw     = d_OCE_his['vfxsnw']/86400.*ICE_rho_sno ; OCE_mas_wfxsnw   = OCE_flux_int ( OCE_wfxsnw ) 
    684     OCE_wfxsub     = d_OCE_his['vfxsub']/86400.*ICE_rho_sno ; OCE_mas_wfxsub   = OCE_flux_int ( OCE_wfxsub ) 
     713    OCE_wfxice     = rprec (d_OCE_his['vfxice'])/86400.*ICE_rho_ice ; OCE_mas_wfxice   = OCE_flux_int ( OCE_wfxice ) 
     714    OCE_wfxsnw     = rprec (d_OCE_his['vfxsnw'])/86400.*ICE_rho_sno ; OCE_mas_wfxsnw   = OCE_flux_int ( OCE_wfxsnw ) 
     715    OCE_wfxsub     = rprec (d_OCE_his['vfxsub'])/86400.*ICE_rho_sno ; OCE_mas_wfxsub   = OCE_flux_int ( OCE_wfxsub ) 
    685716# Additional checks 
    686 OCE_evap_oce   = d_OCE_his['evap_ao_cea'] ; OCE_mas_evap_oce   = OCE_flux_int ( OCE_evap_oce ) 
    687 ICE_evap_ice   = d_OCE_his['subl_ai_cea'] ; ICE_mas_evap_ice   = OCE_flux_int ( ICE_evap_ice ) 
    688 OCE_snow_oce   = d_OCE_his['snow_ao_cea'] ; OCE_mas_snow_oce   = OCE_flux_int ( OCE_snow_oce ) 
    689 OCE_snow_ice   = d_OCE_his['snow_ai_cea'] ; OCE_mas_snow_ice   = OCE_flux_int ( OCE_snow_ice ) 
    690 OCE_rain       = d_OCE_his['rain']        ; OCE_mas_rain       = OCE_flux_int ( OCE_rain     ) 
    691 ICE_wfxsub_err = d_ICE_his['vfxsub_err'] ; ICE_mas_wfxsub_err = OCE_flux_int ( ICE_wfxsub_err ) 
     717OCE_evap_oce   = rprec (d_OCE_his['evap_ao_cea']) ; OCE_mas_evap_oce   = OCE_flux_int ( OCE_evap_oce ) 
     718ICE_evap_ice   = rprec (d_OCE_his['subl_ai_cea']) ; ICE_mas_evap_ice   = OCE_flux_int ( ICE_evap_ice ) 
     719OCE_snow_oce   = rprec (d_OCE_his['snow_ao_cea']) ; OCE_mas_snow_oce   = OCE_flux_int ( OCE_snow_oce ) 
     720OCE_snow_ice   = rprec (d_OCE_his['snow_ai_cea']) ; OCE_mas_snow_ice   = OCE_flux_int ( OCE_snow_ice ) 
     721OCE_rain       = rprec (d_OCE_his['rain']       ) ; OCE_mas_rain       = OCE_flux_int ( OCE_rain     ) 
     722ICE_wfxsub_err = rprec (d_ICE_his['vfxsub_err'] ) ; ICE_mas_wfxsub_err = OCE_flux_int ( ICE_wfxsub_err ) 
    692723if NEMO == 4.0 or NEMO == 4.2 : 
    693     ICE_wfxpnd     = d_ICE_his['vfxpnd']    ; ICE_mas_wfxpnd     = OCE_flux_int ( ICE_wfxpnd     ) 
    694     ICE_wfxsnw_sub = d_ICE_his['vfxsnw_sub'] ; ICE_mas_wfxsnw_sub = OCE_flux_int ( ICE_wfxsnw_sub ) 
    695     ICE_wfxsnw_pre = d_ICE_his['vfxsnw_pre'] ; ICE_mas_wfxsnw_pre = OCE_flux_int ( ICE_wfxsnw_pre ) 
     724    ICE_wfxpnd     = rprec (d_ICE_his['vfxpnd']    ) ; ICE_mas_wfxpnd     = OCE_flux_int ( ICE_wfxpnd     ) 
     725    ICE_wfxsnw_sub = rprec (d_ICE_his['vfxsnw_sub']) ; ICE_mas_wfxsnw_sub = OCE_flux_int ( ICE_wfxsnw_sub ) 
     726    ICE_wfxsnw_pre = rprec (d_ICE_his['vfxsnw_pre']) ; ICE_mas_wfxsnw_pre = OCE_flux_int ( ICE_wfxsnw_pre ) 
    696727if NEMO == 3.6 : 
    697728    ICE_wfxpnd = 0.0 ; ICE_mas_wfxpnd = 0.0 
    698     ICE_wfxsnw_sub = d_ICE_his['vfxsub']/86400.*ICE_rho_sno  ; ICE_mas_wfxsnw_sub = OCE_flux_int ( ICE_wfxsnw_sub ) 
    699     ICE_wfxsnw_pre = d_ICE_his['vfxspr']/86400.*ICE_rho_sno  ; ICE_mas_wfxsnw_pre = OCE_flux_int ( ICE_wfxsnw_pre     ) 
    700  
    701 OCE_wfcorr    = d_OCE_his['wfcorr'] ; OCE_mas_wfcorr  = OCE_flux_int ( OCE_wfcorr ) 
     729    ICE_wfxsnw_sub = rprec (d_ICE_his['vfxsub'])/86400.*ICE_rho_sno  ; ICE_mas_wfxsnw_sub = OCE_flux_int ( ICE_wfxsnw_sub ) 
     730    ICE_wfxsnw_pre = rprec (d_ICE_his['vfxspr'])/86400.*ICE_rho_sno  ; ICE_mas_wfxsnw_pre = OCE_flux_int ( ICE_wfxsnw_pre     ) 
     731 
     732OCE_wfcorr    = rprec (d_OCE_his['wfcorr']) ; OCE_mas_wfcorr  = OCE_flux_int ( OCE_wfcorr ) 
    702733if OCE_relax  : 
    703734    # ssr and fwb are included in emp=>empmr but not in emp_oce (outputed by sea-ice) 
    704     OCE_vflx_fwb = d_OCE_his['vflx_fwb'] ; OCE_mas_vflx_fwb   = OCE_flux_int ( OCE_vflx_fwb ) 
    705     OCE_vflx_ssr = d_OCE_his['vflx_ssr'] ; OCE_mas_vflx_ssr   = OCE_flux_int ( OCE_vflx_ssr ) 
     735    OCE_vflx_fwb = rprec (d_OCE_his['vflx_fwb']) ; OCE_mas_vflx_fwb   = OCE_flux_int ( OCE_vflx_fwb ) 
     736    OCE_vflx_ssr = rprec (d_OCE_his['vflx_ssr']) ; OCE_mas_vflx_ssr   = OCE_flux_int ( OCE_vflx_ssr ) 
    706737else :  
    707738    OCE_fwb = 0.0 ; OCE_mas_fwb = 0.0 
    708739    OCE_ssr = 0.0 ; OCE_mas_ssr = 0.0 
    709740if OCE_icb : 
    710     OCE_berg_icb    = d_OCE_his['berg_floating_melt'] ; OCE_mas_berg_icb = OCE_flux_int ( OCE_berg_icb    ) 
    711     OCE_calving_icb = d_OCE_his['calving_icb']        ; OCE_mas_calv_icb = OCE_flux_int ( OCE_calving_icb ) 
     741    OCE_berg_icb    = rprec (d_OCE_his['berg_floating_melt']) ; OCE_mas_berg_icb = OCE_flux_int ( OCE_berg_icb    ) 
     742    OCE_calving_icb = rprec (d_OCE_his['calving_icb']       ) ; OCE_mas_calv_icb = OCE_flux_int ( OCE_calving_icb ) 
    712743else : 
    713744    OCE_berg_icb = 0. ; OCE_mas_berg_icb = 0. 
    … …  
    715746 
    716747OCE_mas_emp = OCE_mas_emp_oce - OCE_mas_wfxice - OCE_mas_wfxsnw - ICE_mas_wfxpnd - ICE_mas_wfxsub_err 
    717  
    718 OCE_mas_all = OCE_mas_emp_oce +OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf 
     748OCE_mas_all = OCE_mas_emp_oce + OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf 
    719749 
    720750echo ('\n-- Fields:' )  
    … …  
    748778prtFlux ('  CALV_ICB     ', OCE_mas_calv_icb  , 'e', True) 
    749779 
     780echo (' ') 
     781if Coupled : 
     782    prtFlux ( 'Bilan ocean :', -OCE_mas_emp_oce - OCE_mas_emp_ice + OCE_mas_runoffs + OCE_mas_iceshelf, 'e', True ) 
     783else : 
     784    prtFlux ( 'Bilan ocean :', -OCE_mas_emp_oce - OCE_mas_emp_ice + OCE_mas_runoffs + OCE_mas_iceberg + OCE_mas_calving + OCE_mas_iceshelf, 'e', True ) 
     785 
     786 
    750787echo     ('\n===>  Comparing ===>' )  
    751788echo     ('  WFX OCE                     = -empmr + iceshelf                                            = {:12.5e} (kg) '.format ( -OCE_mas_empmr + OCE_mas_iceshelf) ) 
    … …  
    770807echo ('                        = (dSEA_mas_wat + emp_oce + emp_ice - runoffs - iceshelf)*/abs(dSEA_mas_wat)      = {:12.1e} (-)  '.format ( (dSEA_mas_wat + OCE_mas_emp_oce + OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf)/abs(dSEA_mas_wat) ) ) 
    771808 
    772 prtFlux ('   Leak OCE             ', ( dOCE_mas_wat + OCE_mas_empmr - OCE_mas_iceshelf) , 'e' ) 
    773 prtFlux ('   Leak ICE+SNW+PND     ', ( dICE_mas_wat + OCE_mas_emp_ice + OCE_mas_wfxice + OCE_mas_wfxsnw + ICE_mas_wfxpnd + ICE_mas_wfxsub_err ) , 'e' ) 
    774 prtFlux ('   Leak OCE+ICE+SNW+PND ',  ( dSEA_mas_wat + OCE_mas_emp_oce +OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf ) , 'e' ) 
     809prtFlux ('   Leak OCE             ', ( dOCE_mas_wat + OCE_mas_empmr - OCE_mas_iceshelf) , 'e', True ) 
     810prtFlux ('   Leak ICE+SNW+PND     ', ( dICE_mas_wat + OCE_mas_emp_ice + OCE_mas_wfxice + OCE_mas_wfxsnw + ICE_mas_wfxpnd + ICE_mas_wfxsub_err ) , 'e', True ) 
     811prtFlux ('   Leak OCE+ICE+SNW+PND ',  ( dSEA_mas_wat + OCE_mas_emp_oce +OCE_mas_emp_ice - OCE_mas_runoffs - OCE_mas_iceshelf ) , 'e', True ) 
    775812 
    776813 
Note: See TracChangeset for help on using the changeset viewer.