Changeset 6764


Ignore:
Timestamp:
03/19/24 13:00:07 (6 weeks ago)
Author:
omamce
Message:

O.M. : water budget - version mars 2024

Location:
TOOLS/WATER_BUDGET
Files:
1 added
5 edited

Legend:

Unmodified
Added
Removed
  • TOOLS/WATER_BUDGET/ATM_waterbudget.py

    r6688 r6764  
    5050    FullIniFile = 'ATM_' + IniFile 
    5151 
    52 print ("Output file : ", FullIniFile ) 
     52print ("FullIniFile : ", FullIniFile ) 
    5353 
    5454## Experiment parameters : read in .ini file 
    5555## ----------------------------------------- 
    56 dpar = wu.ReadConfig ( IniFile ) 
     56dpar = wu.ReadConfig (IniFile) 
    5757 
    5858## Configure all needed parameter from existant parameters 
    5959## ------------------------------------------------------- 
    60 dpar = wu.SetDatesAndFiles ( dpar ) 
     60dpar = wu.SetDatesAndFiles (dpar) 
    6161 
    6262## Output file with water budget diagnostics 
     
    7777## Useful functions 
    7878## ---------------- 
    79 if repr(readPrec) == "<class 'numpy.float64'>" or readPrec == float : 
     79if repr (readPrec) == "<class 'numpy.float64'>" or readPrec == float : 
    8080    def rprec (ptab) : 
    8181        '''This version does nothing 
     
    121121## Open history files 
    122122## ------------------ 
    123 d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 
    124 if SECHIBA : 
    125     d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 
     123d_ATM_his = xr.open_dataset ( file_ATM_his, use_cftime=True, decode_times=True, decode_cf=True )#.squeeze() 
     124if SECHIBA : 
     125    d_SRF_his = xr.open_dataset ( file_SRF_his, use_cftime=True, decode_times=True, decode_cf=True )#.squeeze() 
    126126    if Routing == 'SECHIBA' : d_RUN_his = d_SRF_his 
    127     if Routing == 'SIMPLE'  : d_RUN_his = xr.open_dataset ( file_RUN_his, use_cftime=True, decode_times=True, decode_cf=True ).squeeze() 
     127    if Routing == 'SIMPLE'  : d_RUN_his = xr.open_dataset ( file_RUN_his, use_cftime=True, decode_times=True, decode_cf=True )#.squeeze() 
    128128 
    129129echo ( f'{file_ATM_his = }' ) 
     
    136136echo ( f'\nRun length : {(dtime/np.timedelta64(1, "D")).values:8.2f} days' ) 
    137137dtime_sec = (dtime/np.timedelta64(1, "s")).values.item() # Convert in seconds 
    138 dpar['Experiment']['dtime_sec'] = dtime_sec 
    139138 
    140139##-- Compute length of each period 
    141 dtime_per = (d_ATM_his.time_counter_bounds[:,-1] -  d_ATM_his.time_counter_bounds[:,0] ) 
     140dtime_per = (d_ATM_his.time_counter_bounds[...,-1] -  d_ATM_his.time_counter_bounds[...,0] ) 
    142141echo ( '\nPeriods lengths (days) : {:} days'.format ( (dtime_per/np.timedelta64(1, "D")).values ) ) 
    143142dtime_per_sec = (dtime_per/np.timedelta64(1, "s")).values # In seconds 
    144 dtime_per_sec = xr.DataArray (dtime_per_sec, dims=["time_counter", ], coords=[d_ATM_his.time_counter,] ) 
     143dtime_per_sec = xr.DataArray (dtime_per_sec, dims=("time_counter",), coords=(d_ATM_his.time_counter,) ) 
    145144dtime_per_sec.attrs['unit'] = 's' 
    146 dpar['Experiment']['dtime_per_sec'] = dtime_per_sec 
    147145 
    148146##-- Number of years (approximative) 
     
    393391 
    394392echo ( '\n====================================================================================' ) 
     393echo ( f'-- ATM checking GRAV  -- {Title} ' ) 
     394## Estimation de GRAV 
     395zGrav =  d_ATM_his.psol / d_ATM_his.mass.sum(dim='presnivs') 
     396echo ( f'Grav estimated {zGrav.min().values} {zGrav.max().values} versus specified {GRAV}' ) 
     397 
     398if GRAV > zGrav.max() or GRAV < zGrav.min() : 
     399    raise ValueError ( 'Error of GRAV value' ) 
     400 
     401echo ( '\n====================================================================================' ) 
    395402echo ( f'-- ATM changes in stores  -- {Title} ' ) 
    396403 
     
    401408ATM.Ahyb = d_ATM_his['Ahyb'].squeeze() 
    402409ATM.Bhyb = d_ATM_his['Bhyb'].squeeze() 
     410 
     411#DYN.Ahyb_beg = d_DYN_beg['ap'] 
     412#DYN.Bhyb_beg = d_DYN_beg['bp'] 
     413#DYN.Ahyb_end = d_DYN_end['ap'] 
     414#DYN.Bhyb_end = d_DYN_end['bp'] 
    403415 
    404416echo ( 'Surface pressure' ) 
     
    413425DYN.pres_beg = ATM.Ahyb + ATM.Bhyb * DYN.psol_beg 
    414426DYN.pres_end = ATM.Ahyb + ATM.Bhyb * DYN.psol_end 
     427 
     428#DYN.pres_beg = DYN.Ahyb_beg + DYN.Bhyb_beg *  DYN.psol_beg  
     429#DYN.pres_end = DYN.Ahyb_end + DYN.Bhyb_end *  DYN.psol_end 
    415430 
    416431echo ( 'Check computation of pressure levels' ) 
     
    427442ind[7] = (DYN.pres_end[-1]).max().item() 
    428443 
    429 if any ( ind != 0) : 
     444if any (ind != 0) : 
    430445    echo ( 'All values should be zero' ) 
    431446    echo ( f'(DYN.pres_beg[ 0]-DYN.psol_beg).min().item() = {ind[0]}' ) 
     
    457472 
    458473echo ( 'Atmosphere mass (kg)' ) 
    459 DYN.mass_air_beg = DYN_stock_int ( DYN.mass_beg_2D ) 
    460 DYN.mass_air_end = DYN_stock_int ( DYN.mass_end_2D ) 
     474DYN.mass_air_beg  = DYN_stock_int ( DYN.mass_beg_2D ) 
     475DYN.mass_air_end  = DYN_stock_int ( DYN.mass_end_2D ) 
     476 
     477DYN.mass_psol_beg = DYN_stock_int ( DYN.psol_beg ) / GRAV 
     478DYN.mass_psol_end = DYN_stock_int ( DYN.psol_end ) / GRAV 
    461479 
    462480echo ( 'Atmosphere mass change (kg)' ) 
    463 DYN.diff_mass_air = DYN.mass_air_end - DYN.mass_air_beg 
     481DYN.diff_mass_air  = DYN.mass_air_end  - DYN.mass_air_beg 
     482DYN.diff_mass_psol = DYN.mass_psol_end - DYN.mass_psol_beg 
     483 
     484echo ( f'\nChk psol vs. mass (dynamics)  -- {Title} ' ) 
     485echo ( '------------------------------------------------------------------------------------' ) 
     486echo ( 'DYN.mass_air_beg   = {:12.6e} kg | DYN.mass_psol_beg = {:12.6e} kg | diff: {:12.6e}'.format (DYN.mass_air_beg, DYN.mass_psol_beg, DYN.mass_air_beg-DYN.mass_psol_beg) ) 
     487echo ( 'DYN.mass_air_end   = {:12.6e} kg | DYN.mass_psol_end = {:12.6e} kg | diff: {:12.6e}'.format (DYN.mass_air_end, DYN.mass_psol_end, DYN.mass_air_end-DYN.mass_psol_end) ) 
     488echo ( 'DYN.diff_mass_air  = {:12.6e} kg'.format (DYN.diff_mass_air ) ) 
     489echo ( 'DYN.diff_mass_psol = {:12.6e} kg'.format (DYN.diff_mass_psol) ) 
     490echo ( 'diff               = {:12.6e} kg'.format (DYN.diff_mass_psol-DYN.diff_mass_air) ) 
    464491 
    465492echo ( f'\nChange of atmosphere mass (dynamics)  -- {Title} ' ) 
     
    469496echo ( f'dMass (atm) : {(DYN.diff_mass_air/DYN.mass_air_beg):9.4e} (-)' ) 
    470497 
     498echo ( f'\nChange of atmosphere mass from psol (dynamics)  -- {Title} ' ) 
     499echo ( '------------------------------------------------------------------------------------' ) 
     500echo ( 'DYN.mass_psol_beg = {:12.6e} kg | DYN.mass_psol_end = {:12.6e} kg'.format (DYN.mass_psol_beg, DYN.mass_psol_end) ) 
     501echo ( f'dMass (atm) : {DYN.diff_mass_psol:9.4e} kg' ) 
     502echo ( f'dMass (atm) : {(DYN.diff_mass_psol/DYN.mass_psol_beg):9.4e} (-)' ) 
     503 
    471504echo ( 'Vertical and horizontal integral, and sum of liquid, solid and vapor water phases' ) 
    472505if LMDZ : 
    473506    if 'H2Ov' in d_DYN_beg.variables : 
    474507        echo ('reading LATLON : H2Ov, H2Ol, H2Oi' ) 
    475         DYN.wat_beg = lmdz.geo3point ( (d_DYN_beg['H2Ov']   + d_DYN_beg['H2Ol']  + d_DYN_beg['H2Oi']).isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
    476         DYN.wat_end = lmdz.geo3point ( (d_DYN_end['H2Ov']   + d_DYN_end['H2Ol']  + d_DYN_end['H2Oi']).isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
     508        DYN.wat_beg_gas = lmdz.geo2point ( d_DYN_beg['H2Ov'].isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
     509        DYN.wat_beg_liq = lmdz.geo2point ( d_DYN_beg['H2Ol'].isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
     510        DYN.wat_beg_sol = lmdz.geo2point ( d_DYN_beg['H2Oi'].isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
     511        DYN.wat_end_gas = lmdz.geo2point ( d_DYN_end['H2Ov'].isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
     512        DYN.wat_end_liq = lmdz.geo2point ( d_DYN_end['H2Ol'].isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
     513        DYN.wat_end_sol = lmdz.geo2point ( d_DYN_end['H2Oi'].isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
    477514    if 'H2O_g' in d_DYN_beg.variables : 
    478515        echo ('reading LATLON : H2O_g, H2O_l, H2O_s' ) 
    479         DYN.wat_beg = lmdz.geo3point ( (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s']).isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
    480         DYN.wat_end = lmdz.geo3point ( (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s']).isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
     516        DYN.wat_beg_gas = lmdz.geo2point ( d_DYN_beg['H2O_g'].isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
     517        DYN.wat_beg_liq = lmdz.geo2point ( d_DYN_beg['H2O_l'].isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
     518        DYN.wat_beg_sol = lmdz.geo2point ( d_DYN_beg['H2O_s'].isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
     519        DYN.wat_end_gas = lmdz.geo2point ( d_DYN_end['H2O_g'].isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
     520        DYN.wat_end_liq = lmdz.geo2point ( d_DYN_end['H2O_l'].isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
     521        DYN.wat_end_sol = lmdz.geo2point ( d_DYN_end['H2O_i'].isel(rlonv=slice(0,-1) ), dim1d='cell' ) 
    481522if ICO : 
    482     if   'H2Ov_g' in d_DYN_beg.variables : 
     523    if 'H2Ov_g' in d_DYN_beg.variables : 
    483524        echo ('reading ICO : H2Ov_g, H2Ov_l, H2Ov_s' ) 
    484         DYN.wat_beg = d_DYN_beg['H2Ov_g'] + d_DYN_beg['H2Ov_l'] + d_DYN_beg['H2Ov_s'] 
    485         DYN.wat_end = d_DYN_end['H2Ov_g'] + d_DYN_end['H2Ov_l'] + d_DYN_end['H2Ov_s'] 
     525        DYN.wat_beg_gas = d_DYN_beg['H2Ov_g'] 
     526        DYN.wat_beg_liq = d_DYN_beg['H2Ov_l'] 
     527        DYN.wat_beg_sol = d_DYN_beg['H2Ov_s'] 
     528        DYN.wat_end_gas = d_DYN_end['H2Ov_g'] 
     529        DYN.wat_end_liq = d_DYN_end['H2Ov_l'] 
     530        DYN.wat_end_sol = d_DYN_end['H2Ov_s'] 
    486531    if 'H2O_g' in d_DYN_beg.variables : 
    487532        echo ('reading ICO : H2O_g, H2O_l, H2O_s' ) 
    488         DYN.wat_beg = d_DYN_beg['H2O_g']  + d_DYN_beg['H2O_l']  + d_DYN_beg['H2O_s'] 
    489         DYN.wat_end = d_DYN_end['H2O_g']  + d_DYN_end['H2O_l']  + d_DYN_end['H2O_s'] 
    490     if 'q' in d_DYN_beg.variables : 
    491         echo ('reading ICO : q' ) 
    492         DYN.wat_beg = d_DYN_beg['q'].isel(nq=0) + d_DYN_beg['q'].isel(nq=1) + d_DYN_beg['q'].isel(nq=2) 
    493         DYN.wat_end = d_DYN_end['q'].isel(nq=0) + d_DYN_end['q'].isel(nq=1) + d_DYN_end['q'].isel(nq=2) 
    494  
    495 if 'lev' in DYN.wat_beg.dims : 
    496     DYN.wat_beg = DYN.wat_beg.rename ( {'lev':'sigs'} ) 
    497     DYN.wat_end = DYN.wat_end.rename ( {'lev':'sigs'} ) 
    498  
     533        DYN.wat_beg_gas = d_DYN_beg['H2O_g'] 
     534        DYN.wat_beg_liq = d_DYN_beg['H2O_l'] 
     535        DYN.wat_beg_sol = d_DYN_beg['H2O_s'] 
     536        DYN.wat_end_gas = d_DYN_end['H2O_g'] 
     537        DYN.wat_end_liq = d_DYN_end['H2O_l'] 
     538        DYN.wat_end_sol = d_DYN_end['H2O_s'] 
     539    #if 'q' in d_DYN_beg.variables : 
     540    #    echo ('reading ICO : q' ) 
     541    #    DYN.wat_beg_gas = d_DYN_beg['q'].isel(nq=0) 
     542    #    DYN.wat_beg_liq = d_DYN_beg['q'].isel(nq=1) 
     543    #    DYN.wat_beg_sol = d_DYN_beg['q'].isel(nq=2) 
     544    #    DYN.wat_end_gas = d_DYN_end['q'].isel(nq=0) 
     545    #    DYN.wat_end_liq = d_DYN_end['q'].isel(nq=1) 
     546    #    DYN.wat_end_sol = d_DYN_end['q'].isel(nq=2) 
     547         
     548if 'lev' in DYN.wat_beg_gas.dims : 
     549    DYN.wat_beg_gas = DYN.wat_beg_gas.rename ( {'lev':'sigs'} ) 
     550    DYN.wat_beg_liq = DYN.wat_beg_liq.rename ( {'lev':'sigs'} ) 
     551    DYN.wat_beg_sol = DYN.wat_beg_sol.rename ( {'lev':'sigs'} ) 
     552    DYN.wat_end_gas = DYN.wat_end_gas.rename ( {'lev':'sigs'} ) 
     553    DYN.wat_end_liq = DYN.wat_end_liq.rename ( {'lev':'sigs'} ) 
     554    DYN.wat_end_sol = DYN.wat_end_sol.rename ( {'lev':'sigs'} ) 
     555 
     556DYN.wat_beg = DYN.wat_beg_gas +  DYN.wat_beg_liq +  DYN.wat_beg_sol 
     557DYN.wat_end = DYN.wat_end_gas +  DYN.wat_end_liq +  DYN.wat_end_sol 
     558     
    499559echo ( 'Compute water content : vertical and horizontal integral' ) 
    500  
    501 DYN.wat_beg_2D =  (DYN.mass_beg * DYN.wat_beg).sum (dim='sigs') 
    502 DYN.wat_end_2D =  (DYN.mass_end * DYN.wat_end).sum (dim='sigs') 
    503  
    504 DYN.mass_wat_beg = DYN_stock_int ( DYN.wat_beg_2D ) 
    505 DYN.mass_wat_end = DYN_stock_int ( DYN.wat_end_2D ) 
     560DYN.wat_beg_gas_2D =  (DYN.mass_beg * DYN.wat_beg_gas).sum (dim='sigs') 
     561DYN.wat_beg_liq_2D =  (DYN.mass_beg * DYN.wat_beg_liq).sum (dim='sigs') 
     562DYN.wat_beg_sol_2D =  (DYN.mass_beg * DYN.wat_beg_sol).sum (dim='sigs') 
     563DYN.wat_end_gas_2D =  (DYN.mass_end * DYN.wat_end_gas).sum (dim='sigs') 
     564DYN.wat_end_liq_2D =  (DYN.mass_end * DYN.wat_end_liq).sum (dim='sigs') 
     565DYN.wat_end_sol_2D =  (DYN.mass_end * DYN.wat_end_sol).sum (dim='sigs') 
     566DYN.wat_beg_2D     =  (DYN.mass_beg * DYN.wat_beg    ).sum (dim='sigs') 
     567DYN.wat_end_2D     =  (DYN.mass_end * DYN.wat_end    ).sum (dim='sigs') 
     568 
     569DYN.mass_wat_beg_gas = DYN_stock_int ( DYN.wat_beg_gas_2D ) 
     570DYN.mass_wat_beg_liq = DYN_stock_int ( DYN.wat_beg_liq_2D ) 
     571DYN.mass_wat_beg_sol = DYN_stock_int ( DYN.wat_beg_sol_2D ) 
     572DYN.mass_wat_end_gas = DYN_stock_int ( DYN.wat_end_gas_2D ) 
     573DYN.mass_wat_end_liq = DYN_stock_int ( DYN.wat_end_liq_2D ) 
     574DYN.mass_wat_end_sol = DYN_stock_int ( DYN.wat_end_sol_2D ) 
     575DYN.mass_wat_beg     = DYN_stock_int ( DYN.wat_beg_2D     ) 
     576DYN.mass_wat_end     = DYN_stock_int ( DYN.wat_end_2D     ) 
    506577 
    507578echo ( 'Variation of water content' ) 
    508 DYN.diff_mass_wat = DYN.mass_wat_end - DYN.mass_wat_beg 
    509  
     579DYN.diff_mass_wat_gas = DYN.mass_wat_end_gas - DYN.mass_wat_beg_gas 
     580DYN.diff_mass_wat_liq = DYN.mass_wat_end_liq - DYN.mass_wat_beg_liq 
     581DYN.diff_mass_wat_sol = DYN.mass_wat_end_sol - DYN.mass_wat_beg_sol 
     582DYN.diff_mass_wat     = DYN.mass_wat_end     - DYN.mass_wat_beg 
     583 
     584if 'prw_end' in d_ATM_his : 
     585    echo ( 'Check with pr variables' ) 
     586    if ATM_HIS == 'latlon' : 
     587        echo ( ' latlon case' ) 
     588        ATM.wat_prw     = lmdz.geo2point ( rprec (d_ATM_his ['prw_end']  [-1]), dim1d='cell' ) 
     589        ATM.wat_prlw    = lmdz.geo2point ( rprec (d_ATM_his ['prlw_end'] [-1]), dim1d='cell' ) 
     590        ATM.wat_prsw    = lmdz.geo2point ( rprec (d_ATM_his ['prsw_end'] [-1]), dim1d='cell' ) 
     591        ATM.wat_prbsw   = lmdz.geo2point ( rprec (d_ATM_his ['prbsw_end'][-1]), dim1d='cell' ) 
     592    if ATM_HIS == 'ico' : 
     593        echo ( ' ico case' ) 
     594        ATM.wat_prw     = rprec (d_ATM_his ['prw_end']  [-1]) 
     595        ATM.wat_prlw    = rprec (d_ATM_his ['prlw_end'] [-1]) 
     596        ATM.wat_prsw    = rprec (d_ATM_his ['prsw_end'] [-1]) 
     597        ATM.wat_prbsw   = rprec (d_ATM_his ['prbsw_end'][-1]) 
     598 
     599    ATM.mass_wat_prw   = DYN_stock_int ( ATM.wat_prw   ) 
     600    ATM.mass_wat_prlw  = DYN_stock_int ( ATM.wat_prlw  ) 
     601    ATM.mass_wat_prsw  = DYN_stock_int ( ATM.wat_prsw  ) 
     602    ATM.mass_wat_prbsw = DYN_stock_int ( ATM.wat_prbsw ) 
     603 
     604    ATM.mass_wat_pr = ATM.mass_wat_prw + ATM.mass_wat_prlw + ATM.mass_wat_prsw + ATM.mass_wat_prbsw 
     605         
     606else : 
     607    ATM.wat_pr    = 0 
     608    ATM.wat_prw   = 0 
     609    ATM.wat_prlw  = 0 
     610    ATM.wat_prsw  = 0 
     611    ATM.wat_prbsw = 0 
     612    ATM.mass_wat_prw   = 0 
     613    ATM.mass_wat_prlw  = 0 
     614    ATM.mass_wat_prsw  = 0 
     615    ATM.mass_wat_prbsw = 0 
     616    ATM.mass_wat_pr    = 0 
     617 
     618echo ( f'\nCheck atmosphere water content  -- {Title} ' ) 
     619echo ( '------------------------------------------------------------------------------------' ) 
     620echo ( 'DYN.mass_wat_end_gas = {:12.6e} kg | ATM.mass_wat_prw   = {:12.6e} kg'.format (DYN.mass_wat_end_gas, ATM.mass_wat_prw  ) ) 
     621echo ( 'DYN.mass_wat_end_liq = {:12.6e} kg | ATM.mass_wat_prlw  = {:12.6e} kg'.format (DYN.mass_wat_end_liq, ATM.mass_wat_prlw ) ) 
     622echo ( 'DYN.mass_wat_end_sol = {:12.6e} kg | ATM.mass_wat_prsw  = {:12.6e} kg'.format (DYN.mass_wat_end_sol, ATM.mass_wat_prsw ) ) 
     623echo ( 'DYN.mass_wat_end_sno = {:12.6e} kg | ATM.mass_wat_prbsw = {:12.6e} kg'.format (   0.               , ATM.mass_wat_prbsw) ) 
     624echo ( 'DYN.mass_wat_end     = {:12.6e} kg | ATM.mass_wat_pr    = {:12.6e} kg'.format (DYN.mass_wat_end    , ATM.mass_wat_pr   ) ) 
     625prtFlux ( 'Error water mass (dyn - atm)', DYN.mass_wat_end -  ATM.mass_wat_pr, 'e', True ) 
     626     
    510627echo ( f'\nChange of atmosphere water content (dynamics)  -- {Title} ' ) 
    511628echo ( '------------------------------------------------------------------------------------' ) 
  • TOOLS/WATER_BUDGET/OCE_waterbudget.py

    r6688 r6764  
    1 #!/usr/bin/env python3 
     1!/usr/bin/env python3 
    22### 
    33### Script to check water conservation in NEMO 
     
    5151 
    5252 
    53 print ("Output file : ", FullIniFile ) 
     53print ("Output parameter file : ", FullIniFile) 
    5454 
    5555## Experiment parameters 
     
    5959## Configure all needed parameter from existant parameters 
    6060## ------------------------------------------------------- 
    61 dpar = wu.SetDatesAndFiles ( dpar) 
     61dpar = wu.SetDatesAndFiles (dpar) 
    6262 
    6363## Output file with water budget diagnostics 
     
    8080 
    8181# Degrades precision 
    82 if repr(readPrec) == "<class 'numpy.float64'>" or readPrec == float : 
     82if repr (readPrec) == "<class 'numpy.float64'>" or readPrec == float : 
    8383    def rprec (ptab) : 
    8484        '''This version does nothing 
     
    271271 
    272272ErrorCount = 0 
    273  
    274273ErrorCount += extract_and_rebuild ( file_name=file_OCE_beg, tar_restart=tar_restart_beg, file_dir_comp=FileDirOCE ) 
    275274ErrorCount += extract_and_rebuild ( file_name=file_OCE_end, tar_restart=tar_restart_end, file_dir_comp=FileDirOCE ) 
  • TOOLS/WATER_BUDGET/WaterUtils.py

    r6688 r6764  
    3131import lmdz 
    3232 
    33 def ReadConfig ( ini_file, default_ini_file='defaults.ini' ) : 
     33def ReadConfig (ini_file, default_ini_file='defaults.ini') : 
    3434    ''' 
    3535    Reads experiment parameters 
     
    9292    return config2dict (params) 
    9393     
    94 def SetDatesAndFiles ( pdict ) : 
     94def SetDatesAndFiles (pdict) : 
    9595    ''' 
    9696    From readed experiment parameters, set all needed parameters 
     
    126126     
    127127    pdict['Files']['TmpDir'] = mm['TmpDir'] 
    128     pdict['libIGCM'].update (mm) 
     128    pdict['libIGCM'].update (mm.__dict__) 
    129129         
    130130    ## Complete configuration 
     
    156156            pdict['Config']['readPrec'] = np.float16 
    157157 
    158  
    159158    ## Defines begining and end of experiment 
    160159    ## -------------------------------------- 
     
    162161        pdict['Experiment']['DateBegin'] = f"{pdict['Experiment']['YearBegin']}0101" 
    163162    else : 
    164         YearBegin, MonthBegin, DayBegin = ldate.GetYearMonthDay ( pdict['Experiment']['DateBegin'] ) 
    165         DateBegin = FormatToGregorian (pdict['Experiment']['DateBegin']) 
    166         pdict['Experiment']['YearBegin'] = YearBegin 
     163        pdict['Experiment']['YearBegin'],pdict['Experiment']['MonthBegin'], pdict['Experiment']['DayBegin'] = ldate.GetYearMonthDay ( pdict['Experiment']['DateBegin'] ) 
     164        pdict['Experiment']['DateBegin'] = ldate.ConvertFormatToGregorian (pdict['Experiment']['DateBegin']) 
    167165         
    168166    if not pdict['Experiment']['DateEnd']  : 
    169167        pdict['Experiment']['DateEnd'] = f"{pdict['Experiment']['YearEnd']}1231" 
    170168    else : 
    171         YearEnd, MonthEnd, DayEnd = ldate.GetYearMonthDay ( pdict['Experiment']['DateEnd'] ) 
    172         DateEnd   = FormatToGregorian (pdict['Experiment']['DateEnd']) 
    173         pdict['Experiment']['DateEnd'] = DateEnd 
     169        pdict['Experiment']['YearEnd'], pdict['Experiment']['MonthEnd'], pdict['Experiment']['DayEnd'] = ldate.GetYearMonthDay ( pdict['Experiment']['DateEnd'] ) 
     170        pdict['Experiment']['DateEnd'] = ldate.ConvertFormatToGregorian (pdict['Experiment']['DateEnd']) 
    174171         
    175172    if not pdict['Experiment']['PackFrequency'] : 
    176         PackFrequencypdict['Experiment']['PackFrequency'] = f"{YearEnd - YearBegin + 1}YE" 
     173        pdict['Experiment']['PackFrequency'] = f"{pdict['Experiment']['YearEnd'] - pdict['Experiment']['YearBegin'] + 1}YE" 
    177174    
    178  
    179175    ## Set directory to extract files 
    180176    ## ------------------------------ 
    181177    if not pdict['Files']['FileDir'] : 
    182178         pdict['Files']['FileDir'] = os.path.join ( pdict['Files']['TmpDir'], f"WATER_{pdict['Experiment']['JobName']}" ) 
    183      
    184179    if not os.path.isdir ( pdict['Files']['FileDir'] ) : os.makedirs ( pdict['Files']['FileDir'] ) 
    185180 
     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         
    186189    FileOut = str2value (pdict['Files']['FileOut']) 
    187190    if not FileOut : 
    188         FileOut = f"ATM_waterbudget_{pdict['Experiment']['JobName']}_{pdict['Experiment']['YearBegin']}_{pdict['Experiment']['YearEnd']}" 
     191        FileOut = f"ATM_waterbudget_{pdict['Experiment']['JobName']}_{pdict['Experiment']['DateBegin']}_{pdict['Experiment']['DateEnd']}" 
    189192        if pdict['Experiment']['ICO'] : 
    190193            if pdict['Experiment']['ATM_HIS'] == 'latlon' : FileOut = f'{FileOut}_LATLON' 
     
    203206    echo ( f"FileDir : {pdict['Files']['FileDir']}"      , f_out=f_out ) 
    204207     
    205     echo ( f"\nDealing with {pdict['libIGCM']['L_EXP']}", f_out ) 
     208    echo ( f"\nDealing with {pdict['libIGCM']['L_EXP']=}", f_out ) 
    206209 
    207210    ## Creates model output directory names 
     
    280283    return pdict 
    281284 
    282 def SetRestartNames ( pdict, f_out) : 
     285def SetRestartNames (pdict, f_out) : 
    283286    ''' 
    284287    Defines dates for restart files 
     
    309312        pdict['Files']['TarRestartPeriod_end_DateBeg'] = ldate.DateAddPeriod ( pdict['Files']['TarRestartPeriod_end_DateEnd'], 
    310313                                                                                   f"-{pdict['Experiment']['PackFrequency']}" ) 
    311         pdict['Files']['TarRestartPeriod_end_DateBeg'] = ldate.AddOneDayToDate ( pdict['Files']['TarRestartPeriod_end_DateBeg'], 
    312                                                                                      pdict['Experiment']['CalendarType'] ) 
     314        pdict['Files']['TarRestartPeriod_end_DateBeg'] = ldate.AddOneDayToDate ( pdict['Files']['TarRestartPeriod_end_DateBeg'], pdict['Experiment']['CalendarType'] ) 
    313315         
    314316        pdict['Files']['TarRestartPeriod_end'] = f"{pdict['Files']['TarRestartPeriod_end_DateBeg']}_{pdict['Files']['TarRestartPeriod_end_DateEnd']}" 
     
    624626    The dictionary must have two levels (see configparser) 
    625627    ''' 
    626     zconf = configparser.ConfigParser ( interpolation=configparser.ExtendedInterpolation() ) 
     628    zconf = configparser.ConfigParser (interpolation=configparser.ExtendedInterpolation ()) 
    627629    zconf.optionxform = str # To keep capitals 
    628630    zconf.read_dict(dict2str(pdict)) 
  • TOOLS/WATER_BUDGET/libIGCM_sys.py

    r6688 r6764  
    88  environment variables and commands. 
    99  All those definitions depend on host particularities. 
    10   It manages a stack mechanism and test validity of operations. 
    1110 
    1211  This software is governed by the CeCILL  license under French law and 
     
    4746#     return unDefined 
    4847 
    49 def Mach ( long=False) : 
     48def Mach (long=False) : 
    5049    ''' 
    5150    Find the computer we are on 
     
    8483    return zmach 
    8584 
    86 def config ( JobName=None , TagName=None  , SpaceName=None, ExperimentName=None, 
    87              LongName=None, ModelName=None, ShortName=None, 
    88              Source=None  , Host=None     , ConfigCard=None, User=None, Group=None, 
    89              TGCC_User='p86mart', TGCC_Group='gen12006', IDRIS_User='rces009', IDRIS_Group='ces', 
    90              ARCHIVE=None, SCRATCHDIR=None, STORAGE=None, R_IN=None, R_OUT=None, 
    91              R_FIG=None, L_EXP=None, 
    92              R_SAVE=None , R_FIGR=None, R_BUF=None , R_BUFR=None, R_BUF_KSH=None, 
    93              REBUILD_DIR=None, POST_DIR=None, 
    94              ThreddsPrefix=None, R_GRAF=None, DB=None, 
    95              IGCM_OUT=None, IGCM_OUT_name=None, rebuild=None, TmpDir=None, 
    96              Out='dict') : 
    97     ''' 
    98     Defines the libIGCM directories - Returns a dictionnary or a namespace 
    99  
    100     Source : for Spip 
    101           Possibilities : 
    102           Local        : local (~/Data/IGCMG/...) 
    103           TGCC_sshfs   : TGCC disks mounted via sshfs 
    104           TGCC_thredds : thredds TGCC via IPSL 
    105                           ('https://thredds-su.ipsl.fr/thredds/dodsC/tgcc_thredds/store/...) 
    106  
    107     Exemple of use : 
    108     libIGCM = libIGCM_sys.config ( .... ) 
    109     globals().update ( libIGCM ) 
    110  
    111     ''' 
    112     #if not IGCM_OUT_name : 
    113     #    if Source == 'TGCC_thredds' : IGCM_OUT_name = '' 
    114     #    else                        : IGCM_OUT_name = 'IGCM_OUT' 
    115  
    116     if not Host : Host = Mach (long=False) 
    117  
    118     LocalUser = os.environ ['USER'] 
    119  
    120     # =========================================================================================== 
    121     # Reads config.card if available 
    122     if ConfigCard : 
    123         # Text existence of ConfigCard 
    124         if not os.path.exists (ConfigCard ) : 
    125             raise FileExistsError ( f"File not found : {ConfigCard = }" ) 
    126  
    127         ## Creates parser for reading .ini input file 
    128         MyReader = configparser.ConfigParser (interpolation=configparser.ExtendedInterpolation() ) 
    129         MyReader.optionxform = str # To keep capitals 
    130  
    131         MyReader.read (ConfigCard) 
    132  
    133         JobName        = MyReader['UserChoices']['JobName'] 
    134         #----- Short Name of Experiment 
    135         ExperimentName = MyReader['UserChoices']['ExperimentName'] 
    136         #----- DEVT TEST PROD 
    137         SpaceName      = MyReader['UserChoices']['SpaceName'] 
    138         LongName       = MyReader['UserChoices']['LongName'] 
    139         ModelName      = MyReader['UserChoices']['ModelName'] 
    140         TagName        = MyReader['UserChoices']['TagName'] 
    141  
    142     ### =========================================================================================== 
    143     ## Part specific to access by OpenDAP/Thredds server 
    144  
    145     # =========================================================================================== 
    146     if Source == 'TGCC_thredds' : 
    147         if not User  and TGCC_User  : User  = TGCC_User 
    148         if not Group and TGCC_Group : Group = TGCC_Group 
    149  
    150         if not ThreddsPrefix : 
    151             ThreddsPrefix = 'https://thredds-su.ipsl.fr/thredds/dodsC/tgcc_thredds' 
    152  
    153         if not ARCHIVE : ARCHIVE = f'{ThreddsPrefix}/store/{TGCC_User}' 
    154         if not R_FIG   : R_FIG   = f'{ThreddsPrefix}/work/{TGCC_User}' 
    155         if not R_IN    : R_IN    = f'{ThreddsPrefix}/work/igcmg/IGCM' 
    156         if not R_GRAF  : R_GRAF  = f'{ThreddsPrefix}/work/p86mart/GRAF/DATA' 
    157  
    158     # =========================================================================================== 
    159     if Source == 'IDRIS_thredds' : 
    160         if not User  and IDRIS_User  : User  = IDRIS_User 
    161         if not Group and IDRIS_Group : Group = IDRIS_Group 
    162  
    163         if not ThreddsPrefix : 
    164             ThreddsPrefix = 'https://thredds-su.ipsl.fr/thredds/catalog/idris_thredds' 
    165  
    166         if not ARCHIVE : ARCHIVE = f'{ThreddsPrefix}/store/{IDRIS_User}' 
    167         if not R_FIG   : R_FIG   = f'{ThreddsPrefix}/work/{IDRIS_User}' 
    168         if not R_IN    : R_IN    = f'{ThreddsPrefix}/work/igcmg/IGCM' 
    169         if not R_GRAF  : R_GRAF  = \ 
    170             'https://thredds-su.ipsl.fr/thredds/dodsC/tgcc_thredds/work/p86mart/GRAF/DATA' 
    171  
    172     ### =========================================================================================== 
    173     ## Machine dependant part 
    174  
    175     # =========================================================================================== 
    176     if Host == 'Spip' : 
    177         if not User : User = LocalUser 
    178         IGCM_OUT_name = 'IGCM_OUT' 
    179         if not ARCHIVE    : ARCHIVE     = os.path.join ( Path.home (), 'Data' ) 
    180         if not SCRATCHDIR : SCRATCHDIR  = os.path.join ( Path.home (), 'Scratch' ) 
    181         if not R_BUF      : R_BUF       = os.path.join ( Path.home (), 'Scratch' ) 
    182         if not R_FIG      : R_FIG       = os.path.join ( Path.home (), 'Data'    ) 
    183  
    184         if not STORAGE    : STORAGE     = ARCHIVE 
    185         if not R_IN       : R_IN        = os.path.join ( '/', 'Users', 'marti', 'Data', 'IGCM' ) 
    186         if not R_GRAF     : R_GRAF      = os.path.join ( '/', 'Users', 'marti', 'GRAF', 'DATA' ) 
    187         if not DB         : DB          = os.path.join ( '/', 'Users', 'marti', 'GRAF', 'DB'   ) 
    188         if not TmpDir     : TmpDir      = os.path.join ( '/', 'Users', LocalUser, 'Scratch' ) 
    189  
    190     # =========================================================================================== 
    191     if ( 'Irene' in Host ) or ( 'Rome' in Host ) : 
    192         
    193         LocalHome  = subprocess.getoutput ( 'ccc_home --ccchome' ) 
    194         LocalGroup = os.path.basename ( os.path.dirname (LocalHome)) 
    195         if not User or User == 'marti' : 
    196             if not TGCC_User : User = LocalUser 
    197             else             : User = TGCC_User 
    198  
    199         if not Group : 
    200             if TGCC_Group : Group = TGCC_Group 
    201             else          : Group = LocalGroup 
    202  
    203         IGCM_OUT_name = 'IGCM_OUT' 
    204  
    205         if not R_IN        : 
    206             R_IN       = os.path.join ( subprocess.getoutput ( 
    207                 'ccc_home --cccwork -d igcmg -u igcmg' ), 'IGCM') 
    208         if not ARCHIVE     : 
    209             ARCHIVE    = subprocess.getoutput ( 
    210                 f'ccc_home --cccstore   -u {User} -d {Group}' ) 
    211         if not STORAGE     : 
    212             STORAGE    = subprocess.getoutput ( 
    213                 f'ccc_home --cccwork    -u {User} -d {Group}' ) 
    214         if not SCRATCHDIR  : 
    215             SCRATCHDIR = subprocess.getoutput ( 
    216                 f'ccc_home --cccscratch -u {User} -d {Group}' ) 
    217         if not R_BUF       : 
    218             R_BUF      = subprocess.getoutput ( 
    219                 f'ccc_home --cccscratch -u {User} -d {Group}' ) 
    220         if not R_FIG       : 
    221             R_FIG      = subprocess.getoutput ( 
    222                 f'ccc_home --cccwork    -u {User} -d {Group}' ) 
    223         if not R_GRAF      : 
    224             R_GRAF     = os.path.join ( subprocess.getoutput ( 
    225                 'ccc_home --cccwork -d drf -u p86mart'), 'GRAF', 'DATA' ) 
    226         if not DB          : 
    227             DB         = os.path.join ( subprocess.getoutput ( 
    228                 'ccc_home --cccwork -d igcmg -u igcmg'), 'database' ) 
    229         if not rebuild : 
     85class config : 
     86    def __str__ (self) : 
     87        return str (self.__dict__) 
     88 
     89    def __repr__ (self) : 
     90        return str (self.__dict__) 
     91 
     92    def __name__ (self) : 
     93        return self.__class__.__name__ 
     94 
     95    def __getitem__ (self, key): 
     96        return getattr (self, key) 
     97 
     98    def __iter__ (self) : 
     99        return self 
     100 
     101    def __next__ (self) : 
     102        if self.index == 0: 
     103            raise StopIteration 
     104        self.index = self.index - 1 
     105        return self.data[self.index] 
     106 
     107    def __init__ (self, JobName=None, TagName=None, SpaceName=None, ExperimentName=None, 
     108                  LongName=None, ModelName=None, ShortName=None, 
     109                  Source=None, Host=None, ConfigCard=None, User=None, Group=None, 
     110                  TGCC_User='p86mart', TGCC_Group='gen12006', IDRIS_User='rces009', IDRIS_Group='ces', 
     111                  ARCHIVE=None, SCRATCHDIR=None, STORAGE=None, R_IN=None, R_OUT=None, 
     112                  R_FIG=None, L_EXP=None, 
     113                  R_SAVE=None, R_FIGR=None, R_BUF=None, R_BUFR=None, R_BUF_KSH=None, 
     114                  REBUILD_DIR=None, POST_DIR=None, 
     115                  ThreddsPrefix=None, R_GRAF=None, DB=None, 
     116                  IGCM_OUT=None, IGCM_OUT_name=None, rebuild=None, TmpDir=None) : 
     117        ''' 
     118        Defines the libIGCM directories - Returns a dictionary or a namespace 
     119         
     120        Source : for Spip 
     121        Possibilities : 
     122        Local        : local (~/Data/IGCMG/...) 
     123        TGCC_sshfs   : TGCC disks mounted via sshfs 
     124        TGCC_thredds : thredds TGCC via IPSL 
     125        ('https://thredds-su.ipsl.fr/thredds/dodsC/tgcc_thredds/store/...) 
     126         
     127        Exemple of use : 
     128        libIGCM = libIGCM_sys.config ( .... ) 
     129        globals().update ( libIGCM ) 
     130        ''' 
     131         
     132        if not Host : Host = Mach (long=False) 
     133             
     134        LocalUser = os.environ ['USER'] 
     135             
     136        # =========================================================================================== 
     137        # Reads config.card if available 
     138        if ConfigCard : 
     139            # Text existence of ConfigCard 
     140            if not os.path.exists (ConfigCard ) : 
     141                raise FileExistsError ( f"File not found : {ConfigCard = }" ) 
     142                 
     143            ## Creates parser for reading .ini input file 
     144            MyReader = configparser.ConfigParser (interpolation=configparser.ExtendedInterpolation() ) 
     145            MyReader.optionxform = str # To keep capitals 
     146             
     147            MyReader.read (ConfigCard) 
     148             
     149            JobName        = MyReader['UserChoices']['JobName'] 
     150            #----- Short Name of Experiment 
     151            ExperimentName = MyReader['UserChoices']['ExperimentName'] 
     152            #----- DEVT TEST PROD 
     153            SpaceName      = MyReader['UserChoices']['SpaceName'] 
     154            LongName       = MyReader['UserChoices']['LongName'] 
     155            ModelName      = MyReader['UserChoices']['ModelName'] 
     156            TagName        = MyReader['UserChoices']['TagName'] 
     157             
     158        ### =========================================================================================== 
     159        ## Part specific to access by OpenDAP/Thredds server 
     160         
     161        # =========================================================================================== 
     162        if Source == 'TGCC_thredds' : 
     163            if not User  and TGCC_User  : User  = TGCC_User 
     164            if not Group and TGCC_Group : Group = TGCC_Group 
     165                 
     166            if not ThreddsPrefix : 
     167                ThreddsPrefix = 'https://thredds-su.ipsl.fr/thredds/dodsC/tgcc_thredds' 
     168                 
     169            if not ARCHIVE : ARCHIVE = f'{ThreddsPrefix}/store/{TGCC_User}' 
     170            if not R_FIG   : R_FIG   = f'{ThreddsPrefix}/work/{TGCC_User}' 
     171            if not R_IN    : R_IN    = f'{ThreddsPrefix}/work/igcmg/IGCM' 
     172            if not R_GRAF  : R_GRAF  = f'{ThreddsPrefix}/work/p86mart/GRAF/DATA' 
     173                 
     174        # =========================================================================================== 
     175        if Source == 'IDRIS_thredds' : 
     176            if not User  and IDRIS_User  : User  = IDRIS_User 
     177            if not Group and IDRIS_Group : Group = IDRIS_Group 
     178                 
     179            if not ThreddsPrefix : 
     180                ThreddsPrefix = 'https://thredds-su.ipsl.fr/thredds/catalog/idris_thredds' 
     181                 
     182            if not ARCHIVE : ARCHIVE = f'{ThreddsPrefix}/store/{IDRIS_User}' 
     183            if not R_FIG   : R_FIG   = f'{ThreddsPrefix}/work/{IDRIS_User}' 
     184            if not R_IN    : R_IN    = f'{ThreddsPrefix}/work/igcmg/IGCM' 
     185            if not R_GRAF  : R_GRAF  = 'https://thredds-su.ipsl.fr/thredds/dodsC/tgcc_thredds/work/p86mart/GRAF/DATA' 
     186               
     187        ### =========================================================================================== 
     188        ## Machine dependant part 
     189         
     190        # =========================================================================================== 
     191        if Host == 'Spip' : 
     192            if not User : User = LocalUser 
     193            IGCM_OUT_name = 'IGCM_OUT' 
     194            if not ARCHIVE    : ARCHIVE     = os.path.join ( Path.home (), 'Data' ) 
     195            if not SCRATCHDIR : SCRATCHDIR  = os.path.join ( Path.home (), 'Scratch' ) 
     196            if not R_BUF      : R_BUF       = os.path.join ( Path.home (), 'Scratch' ) 
     197            if not R_FIG      : R_FIG       = os.path.join ( Path.home (), 'Data'    ) 
     198                 
     199            if not STORAGE    : STORAGE     = ARCHIVE 
     200            if not R_IN       : R_IN        = os.path.join ( '/', 'Users', 'marti', 'Data', 'IGCM' ) 
     201            if not R_GRAF     : R_GRAF      = os.path.join ( '/', 'Users', 'marti', 'GRAF', 'DATA' ) 
     202            if not DB         : DB          = os.path.join ( '/', 'Users', 'marti', 'GRAF', 'DB'   ) 
     203            if not TmpDir     : TmpDir      = os.path.join ( Path.home (), 'Scratch' ) 
     204                 
     205        # =========================================================================================== 
     206        if ( 'Irene' in Host ) or ( 'Rome' in Host ) : 
     207         
     208            LocalHome  = subprocess.getoutput ( 'ccc_home --ccchome' ) 
     209            LocalGroup = os.path.basename ( os.path.dirname (LocalHome)) 
     210            if not User or User == 'marti' : 
     211                if not TGCC_User : User = LocalUser 
     212                else             : User = TGCC_User 
     213                     
     214            if not Group : 
     215                if TGCC_Group : Group = TGCC_Group 
     216                else          : Group = LocalGroup 
     217                     
     218            IGCM_OUT_name = 'IGCM_OUT' 
     219             
     220            if not R_IN        : 
     221                R_IN       = os.path.join ( subprocess.getoutput ( 
     222                    'ccc_home --cccwork -d igcmg -u igcmg' ), 'IGCM') 
     223            if not ARCHIVE     : 
     224                ARCHIVE    = subprocess.getoutput ( 
     225                    f'ccc_home --cccstore   -u {User} -d {Group}' ) 
     226            if not STORAGE     : 
     227                STORAGE    = subprocess.getoutput ( 
     228                    f'ccc_home --cccwork    -u {User} -d {Group}' ) 
     229            if not SCRATCHDIR  : 
     230                SCRATCHDIR = subprocess.getoutput ( 
     231                    f'ccc_home --cccscratch -u {User} -d {Group}' ) 
     232            if not R_BUF       : 
     233                R_BUF      = subprocess.getoutput ( 
     234                    f'ccc_home --cccscratch -u {User} -d {Group}' ) 
     235            if not R_FIG       : 
     236                R_FIG      = subprocess.getoutput ( 
     237                    f'ccc_home --cccwork    -u {User} -d {Group}' ) 
     238            if not R_GRAF      : 
     239                R_GRAF     = os.path.join ( subprocess.getoutput ( 
     240                    'ccc_home --cccwork -d drf -u p86mart'), 'GRAF', 'DATA' ) 
     241            if not DB          : 
     242                DB         = os.path.join ( subprocess.getoutput ( 
     243                    'ccc_home --cccwork -d igcmg -u igcmg'), 'database' ) 
     244            if not rebuild : 
    230245                rebuild = os.path.join ( 
    231246                    subprocess.getoutput ( 'ccc_home --ccchome -d igcmg -u igcmg' ), 
    232247                    'Tools', 'irene', 'rebuild_nemo', 'bin', 'rebuild_nemo' ) 
    233248                 
    234         if not TmpDir : TmpDir = subprocess.getoutput ( f'ccc_home --cccscratch' ) 
    235  
    236     # =========================================================================================== 
    237     if Host == 'SpiritJ' : 
    238         if not User  : 
    239             if TGCC_User  : User = TGCC_User 
    240             else          : User = LocalUser 
    241         if not ARCHIVE    : 
    242             ARCHIVE    = os.path.join ( '/', 'thredds', 'tgcc', 'store', User ) 
    243         if not  STORAGE   : 
    244             STORAGE    = os.path.join ( '/', 'thredds', 'tgcc', 'work' , User ) 
    245         #if not SCRATCHDIR : SCRATCHDIR = os.path.join ( '/', 'thredds'  , 'tgcc', 'store', User ) 
    246         if not R_IN       : 
    247             R_IN       = os.path.join ( '/', 'projsu', 'igcmg', 'IGCM' ) 
    248         #if not R_GRAF     : R_GRAF     = os.path.join ('/', 'data', 'omamce', 'GRAF', 'DATA' ) 
    249         if not R_GRAF     : 
    250             R_GRAF     = os.path.join  ( '/', 'thredds', 'tgcc', 'work', 'p86mart', 'GRAF', 'DATA' ) 
    251         if not DB         : 
    252             DB         = os.path.join  ( '/', 'data', 'igcmg', 'database' ) 
    253         if not TmpDir     : 
    254             TmpDir = os.path.join ( '/', 'data', LocalUser ) 
    255  
    256     # =========================================================================================== 
    257     if Host == 'SpiritX' : 
    258         if not User  : 
    259             if TGCC_User  : User  = TGCC_User 
    260             else          : User = os.environ ['USER'] 
    261         if not ARCHIVE    : 
    262             ARCHIVE    = os.path.join ( '/', 'thredds', 'tgcc', 'store', User ) 
    263         if not  STORAGE   : 
    264             STORAGE    = os.path.join ( '/', 'thredds', 'tgcc', 'work' , User ) 
    265         #if not SCRATCHDIR : 
    266         #  SCRATCHDIR = os.path.join ( '/', 'thredds', 'tgcc', 'store', User ) 
    267         if not R_IN       : 
    268             R_IN       = os.path.join ( '/', 'projsu', 'igcmg', 'IGCM' ) 
    269         #if not R_GRAF     : 
    270         #    R_GRAF     = os.path.join ('/', 'data', 'omamce', 'GRAF', 'DATA' ) 
    271         if not R_GRAF     : 
    272             R_GRAF     = os.path.join  ( '/', 'thredds', 'tgcc', 'work' , 'p86mart', 'GRAF', 'DATA' ) 
    273         if not DB         : 
    274             DB         = os.path.join  ( '/', 'data', 'igcmg', 'database' ) 
    275  
    276     # =========================================================================================== 
    277     if Host == 'Jean-Zay' : 
    278         if not User  : User  = os.environ ['USER'] 
    279         LocalGroup = os.path.basename ( os.path.dirname ( Path.home () )) 
    280         if not Group : Group = LocalGroup 
    281  
    282         if not ARCHIVE    : 
    283             ARCHIVE    = os.path.join ( '/', 'gpfsstore'  , 'rech', Group, User ) 
    284         if not STORAGE    : 
    285             STORAGE    = os.path.join ( '/', 'gpfswork'   , 'rech', Group, User ) 
    286         if not SCRATCHDIR : 
    287             SCRATCHDIR = os.path.join ( '/', 'gpfsscratch', 'rech', Group, User ) 
    288         if not R_FIG      : 
    289             R_FIG      = os.path.join ( '/', 'cccwork'    , 'rech', Group, User ) 
    290         if not R_BUF      : 
    291             R_BUF      = os.path.join ( '/', 'gpfsscratch', 'rech', Group, User ) 
    292         if not R_IN       : 
    293             R_IN       = os.path.join ( '/', 'gpfswork'   , 'rech', 'psl', 'commun', 'IGCM' ) 
    294         if not R_GRAF     : 
    295             R_GRAF     = os.path.join ( '/', 'gpfswork'   , 'rech', Group, User, 'GRAF', 'DATA' ) 
    296         if not DB         : 
    297             DB         = os.path.join ( '/', 'gpfswork'   , 'rech', 'psl', 'commun', 'database' ) 
    298         if not rebuild : 
    299             rebuild = os.path.join ( '/', 'gpfswork', 'rech', 'psl', 'commun', 'Tools', 
    300                                          'rebuild', 'modipsl_IOIPSL_PLUS_v2_2_4', 'bin', 'rebuild' ) 
    301         if not TmpDir : TmpDir = os.path.join ( '/', 'gpfsscratch', 'rech', 
    302                             os.path.basename ( os.path.dirname ( Path.home () )), LocalUser ) 
    303  
    304     ### =========================================================================================== 
    305     ### The construction of the following variables is not machine dependant 
    306     ### =========================================================================================== 
    307     if SpaceName == 'TEST' : 
    308         if SCRATCHDIR and not R_OUT : R_OUT = SCRATCHDIR 
    309         if SCRATCHDIR and not R_FIG : R_FIG = SCRATCHDIR 
    310     else  : 
    311         if ARCHIVE    and not R_OUT : R_OUT = ARCHIVE 
    312         if STORAGE    and not R_FIG : R_FIG = STORAGE 
    313     if IGCM_OUT_name : 
    314         R_OUT = os.path.join ( R_OUT, IGCM_OUT_name ) 
    315         R_FIG = os.path.join ( R_FIG, IGCM_OUT_name ) 
    316  
    317     if SCRATCHDIR and not R_BUF : R_BUF  = os.path.join ( SCRATCHDIR, IGCM_OUT_name ) 
    318     if not IGCM_OUT : IGCM_OUT = R_OUT 
    319  
    320     if TagName and SpaceName and ExperimentName and JobName : 
    321         if not L_EXP : L_EXP = os.path.join ( TagName, SpaceName, ExperimentName, JobName ) 
    322  
    323         if R_OUT and not R_SAVE : R_SAVE      = os.path.join ( R_OUT  , L_EXP ) 
     249            if not TmpDir : TmpDir = subprocess.getoutput ( f'ccc_home --cccscratch' ) 
     250                 
     251        # =========================================================================================== 
     252        if Host == 'SpiritJ' : 
     253            if not User  : 
     254                if TGCC_User  : User = TGCC_User 
     255                else          : User = LocalUser 
     256            if not ARCHIVE    : 
     257                ARCHIVE    = os.path.join ( '/', 'thredds', 'tgcc', 'store', User ) 
     258            if not  STORAGE   : 
     259                STORAGE    = os.path.join ( '/', 'thredds', 'tgcc', 'work' , User ) 
     260            #if not SCRATCHDIR : SCRATCHDIR = os.path.join ( '/', 'thredds'  , 'tgcc', 'store', User ) 
     261            if not R_IN       : 
     262                R_IN       = os.path.join ( '/', 'projsu', 'igcmg', 'IGCM' ) 
     263            #if not R_GRAF     : R_GRAF     = os.path.join ('/', 'data', 'omamce', 'GRAF', 'DATA' ) 
     264            if not R_GRAF     : 
     265                R_GRAF     = os.path.join  ( '/', 'thredds', 'tgcc', 'work', 'p86mart', 'GRAF', 'DATA' ) 
     266            if not DB         : 
     267                DB         = os.path.join  ( '/', 'data', 'igcmg', 'database' ) 
     268            if not TmpDir     : 
     269                TmpDir = os.path.join ( '/', 'data', LocalUser ) 
     270                 
     271        # =========================================================================================== 
     272        if Host == 'SpiritX' : 
     273            if not User  : 
     274                if TGCC_User  : User  = TGCC_User 
     275                else          : User = os.environ ['USER'] 
     276            if not ARCHIVE    : 
     277                ARCHIVE    = os.path.join ( '/', 'thredds', 'tgcc', 'store', User ) 
     278            if not  STORAGE   : 
     279                STORAGE    = os.path.join ( '/', 'thredds', 'tgcc', 'work' , User ) 
     280            #if not SCRATCHDIR : 
     281            #  SCRATCHDIR = os.path.join ( '/', 'thredds', 'tgcc', 'store', User ) 
     282            if not R_IN       : 
     283                R_IN       = os.path.join ( '/', 'projsu', 'igcmg', 'IGCM' ) 
     284            if not R_GRAF     : 
     285                R_GRAF     = os.path.join  ( '/', 'thredds', 'tgcc', 'work' , 'p86mart', 'GRAF', 'DATA' ) 
     286            if not DB         : 
     287                DB         = os.path.join  ( '/', 'data', 'igcmg', 'database' ) 
     288                 
     289        # =========================================================================================== 
     290        if Host == 'Jean-Zay' : 
     291            if not User  : User  = os.environ ['USER'] 
     292            LocalGroup = os.path.basename ( os.path.dirname ( Path.home () )) 
     293            if not Group : Group = LocalGroup 
     294                 
     295            if not ARCHIVE    : 
     296                ARCHIVE    = os.path.join ( '/', 'gpfsstore'  , 'rech', Group, User ) 
     297            if not STORAGE    : 
     298                STORAGE    = os.path.join ( '/', 'gpfswork'   , 'rech', Group, User ) 
     299            if not SCRATCHDIR : 
     300                SCRATCHDIR = os.path.join ( '/', 'gpfsscratch', 'rech', Group, User ) 
     301            if not R_FIG      : 
     302                R_FIG      = os.path.join ( '/', 'cccwork'    , 'rech', Group, User ) 
     303            if not R_BUF      : 
     304                R_BUF      = os.path.join ( '/', 'gpfsscratch', 'rech', Group, User ) 
     305            if not R_IN       : 
     306                R_IN       = os.path.join ( '/', 'gpfswork'   , 'rech', 'psl', 'commun', 'IGCM' ) 
     307            if not R_GRAF     : 
     308                R_GRAF     = os.path.join ( '/', 'gpfswork'   , 'rech', Group, User, 'GRAF', 'DATA' ) 
     309            if not DB         : 
     310                DB         = os.path.join ( '/', 'gpfswork'   , 'rech', 'psl', 'commun', 'database' ) 
     311            if not rebuild : 
     312                rebuild = os.path.join ( '/', 'gpfswork', 'rech', 'psl', 'commun', 'Tools', 
     313                                            'rebuild', 'modipsl_IOIPSL_PLUS_v2_2_4', 'bin', 'rebuild' ) 
     314            if not TmpDir : TmpDir = os.path.join ( '/', 'gpfsscratch', 'rech', 
     315                                os.path.basename ( os.path.dirname ( Path.home () )), LocalUser ) 
     316                 
     317        ### =========================================================================================== 
     318        ### The construction of the following variables is not machine dependant 
     319        ### =========================================================================================== 
     320        if SpaceName == 'TEST' : 
     321            if SCRATCHDIR and not R_OUT : R_OUT = SCRATCHDIR 
     322            if SCRATCHDIR and not R_FIG : R_FIG = SCRATCHDIR 
     323        else  : 
     324            if ARCHIVE    and not R_OUT : R_OUT = ARCHIVE 
     325            if STORAGE    and not R_FIG : R_FIG = STORAGE 
    324326        if IGCM_OUT_name : 
    325             if STORAGE and not R_FIGR : R_FIGR      = os.path.join ( STORAGE, IGCM_OUT_name, L_EXP ) 
    326         else : 
    327             if STORAGE and not R_FIGR  : R_FIGR      = os.path.join ( STORAGE, L_EXP ) 
    328         if R_BUF   and not R_BUFR      : R_BUFR      = os.path.join ( R_BUF  , L_EXP ) 
    329         if R_BUFR  and not R_BUF_KSH   : R_BUF_KSH   = os.path.join ( R_BUFR , 'Out' ) 
    330         if R_BUF   and not REBUILD_DIR : REBUILD_DIR = os.path.join ( R_BUF  , L_EXP, 'REBUILD' ) 
    331         if R_BUF   and not POST_DIR    : POST_DIR    = os.path.join ( R_BUF  , L_EXP, 'Out' ) 
    332  
    333     ### =========================================================================================== 
    334     ## Builds the final dictionnary 
    335  
    336     libIGCM = { 'TagName'       : TagName  , 
    337                 'SpaceName'     : SpaceName , 
    338                 'ExperimentName': ExperimentName, 
    339                 'JobName'       : JobName  , 
    340                 'LongName'      : LongName  , 
    341                 'ModelName'     : ModelName    , 
    342                 'ShortName'     : ShortName, 
    343                 'ConfigCard'    : ConfigCard, 
    344                 'ARCHIVE'       : ARCHIVE  , 
    345                 'STORAGE'       : STORAGE   , 
    346                 'SCRATCHDIR'    : SCRATCHDIR    , 
    347                 'R_OUT'         : R_OUT    , 
    348                 'R_BUF'         : R_BUFR    , 
    349                 'R_GRAF'        : R_GRAF        , 
    350                 'DB'            : DB       , 
    351                 'IGCM_OUT'      : IGCM_OUT  , 
    352                 'R_SAVE'        : R_SAVE        , 
    353                 'R_FIGR'        : R_FIGR   , 
    354                 'R_BUFR'        : R_BUFR, 
    355                 'REBUILD_DIR'   : REBUILD_DIR, 
    356                 'POST_DIR'      : POST_DIR , 
    357                 'R_IN'          : R_IN      , 
    358                 'Mach'          : Host     , 
    359                 'Source'        : Source    , 
    360                 'User'          : User    , 
    361                 'Group'         : Group, 
    362                 'IGCM_OUT_name' : IGCM_OUT_name, 
    363                 'rebuild'       : rebuild, 
    364                 'TmpDir'        : TmpDir, 
    365                 'TGCC_User'     : TGCC_User, 
    366                 'TGCC_Group'    : TGCC_Group, 
    367                 'Thredds'       : ThreddsPrefix, 
    368                 'IDRIS_User'    : IDRIS_User, 
    369                 'IDRIS_Group'   : IDRIS_Group } 
    370  
    371     if Out in ['namespace', 'space', 'name', 'Space', 'Name', 'Names', 'NameSpace'] : 
    372         libIGCM = types.SimpleNamespace (**libIGCM) 
    373  
    374     return libIGCM 
     327            R_OUT = os.path.join ( R_OUT, IGCM_OUT_name ) 
     328            R_FIG = os.path.join ( R_FIG, IGCM_OUT_name ) 
     329             
     330        if SCRATCHDIR and not R_BUF : R_BUF  = os.path.join ( SCRATCHDIR, IGCM_OUT_name ) 
     331        if not IGCM_OUT : IGCM_OUT = R_OUT 
     332             
     333        if TagName and SpaceName and ExperimentName and JobName : 
     334            if not L_EXP : L_EXP = os.path.join ( TagName, SpaceName, ExperimentName, JobName ) 
     335                 
     336            if R_OUT and not R_SAVE : R_SAVE      = os.path.join ( R_OUT  , L_EXP ) 
     337            if IGCM_OUT_name : 
     338                if STORAGE and not R_FIGR : R_FIGR      = os.path.join ( STORAGE, IGCM_OUT_name, L_EXP ) 
     339            else : 
     340                if STORAGE and not R_FIGR  : R_FIGR      = os.path.join ( STORAGE, L_EXP ) 
     341            if R_BUF   and not R_BUFR      : R_BUFR      = os.path.join ( R_BUF  , L_EXP ) 
     342            if R_BUFR  and not R_BUF_KSH   : R_BUF_KSH   = os.path.join ( R_BUFR , 'Out' ) 
     343            if R_BUF   and not REBUILD_DIR : REBUILD_DIR = os.path.join ( R_BUF  , L_EXP, 'REBUILD' ) 
     344            if R_BUF   and not POST_DIR    : POST_DIR    = os.path.join ( R_BUF  , L_EXP, 'Out' ) 
     345                 
     346        ### =========================================================================================== 
     347        ## Builds the class components 
     348         
     349        self.TagName        = TagName 
     350        self.SpaceName      = SpaceName 
     351        self.ExperimentName = ExperimentName 
     352        self.JobName        = JobName 
     353        self.LongName       = LongName 
     354        self.ModelName      = ModelName 
     355        self.ShortName      = ShortName 
     356        self.ConfigCard     = ConfigCard 
     357        self.ARCHIVE        = ARCHIVE 
     358        self.STORAGE        = STORAGE 
     359        self.SCRATCHDIR     = SCRATCHDIR 
     360        self.R_OUT          = R_OUT 
     361        self.R_BUF          = R_BUFR 
     362        self.R_GRAF         = R_GRAF 
     363        self.DB             = DB 
     364        self.IGCM_OUT       = IGCM_OUT 
     365        self.R_SAVE         = R_SAVE 
     366        self.R_FIGR         = R_FIGR 
     367        self.R_BUFR         = R_BUFR 
     368        self.REBUILD_DIR    = REBUILD_DIR 
     369        self.POST_DIR       = POST_DIR 
     370        self.R_IN           = R_IN 
     371        self.Mach           = Host 
     372        self.Source         = Source 
     373        self.User           = User 
     374        self.Group          = Group 
     375        self.IGCM_OUT_name  = IGCM_OUT_name 
     376        self.rebuild        = rebuild 
     377        self.TmpDir         = TmpDir 
     378        self.TGCC_User      = TGCC_User 
     379        self.TGCC_Group     = TGCC_Group 
     380        self.Thredds        = ThreddsPrefix 
     381        self.IDRIS_User     = IDRIS_User 
     382        self.IDRIS_Group    = IDRIS_Group 
     383 
     384  
  • TOOLS/WATER_BUDGET/lmdz.py

    r6665 r6764  
    178178 
    179179    if verbose : 
    180         print ( f'{ou_dims   =}'     ) 
     180        print ( f'{ou_dims   =}' ) 
    181181        print ( f'{new_coords=}' ) 
    182182 
     
    261261    return pout 
    262262 
    263 def point2geo (p1d) : 
    264     '''From 1D (restart type) to 2D 
     263def point2geo (p1d, verbose=False, lon=False, lat=False, jpi=0, jpj=0, share_pole=False) : 
     264    ''' 
     265    From 1D (..., points_physiques) (restart type) to 2D (..., lat, lon) 
    265266    ''' 
    266267    math = __mmath__ (p1d) 
    267268 
    268     # Get the dimensions 
     269    # Get the horizontal dimension 
    269270    jpn = p1d.shape[-1] 
    270  
    271     if len (p1d.shape) > 1 : 
    272         jpt = p1d.shape[0] 
    273     else : 
    274         jpt = 0 
    275  
    276     if jpn ==  9026 : 
    277         jpi, jpj =  96, 96 
    278     if jpn == 17858 : 
    279         jpi, jpj = 144, 144 
    280     if jpn == 20306 : 
    281         jpi, jpj = 144, 143 
    282  
    283     if jpt > 0 : 
    284         p2d = np.zeros ( (jpt, jpj, jpi) ) 
    285         p2d [:, 1:jpj-1, :] = np.reshape (p1d [:,1:jpn-1], (jpt, jpj-2, jpi) ) 
    286         p2d [:, 0    , : ] = p1d[:,   0  ] 
    287         p2d [:, jpj-1, : ] = p1d[:, jpn-1] 
    288  
    289     else : 
    290         p2d = np.zeros ( (jpj, jpi) ) 
    291         p2d [1:jpj-1, :] = np.reshape (np.float64 (p1d [1:jpn-1]), (jpj-2, jpi) ) 
    292         p2d [0    , : ] = p1d[   0 ] 
    293         p2d [jpj-1, : ] = p1d[jpn-1] 
     271    # Get other dimension 
     272    form1 = list(p1d.shape [0:-1])  
     273 
     274    if jpi==0 and jpj>0 :  
     275        jpi = (jpn-2)//(jpj-2) 
     276    if jpi>0  and jpj == 0 : 
     277        jpj = (jpn-2)//jpi +2 
     278 
     279    def c_jpn (jpj, jpi) : return jpi*(jpj-2) + 2 
     280             
     281    if jpi==0 and jpi==0 : 
     282        if jpn == c_jpn( 95, 96) : 
     283            jpj, jpi =  95,  96 
     284        if jpn == c_jpn (96, 96) : 
     285            jpj, jpi =  96,  96 
     286        if jpn == c_jpn (144, 144) : 
     287            jpj, jpi = 144, 144 
     288        if jpn == c_jpn (143, 144) : 
     289            jpj, jpi = 143, 144 
     290         
     291    form2 = form1 + [jpj, jpi,] 
     292    form_shape = form1 + [jpj-2, jpi] 
     293    p2d = np.empty ( form2 ) 
     294     
     295    p2d [..., 1:jpj-1, :] = np.reshape (np.float64 (p1d [..., 1:jpn-1]), form_shape ) 
     296    if share_pole : 
     297        for ji in np.arange ( p2d.shape[-1] ) : 
     298            p2d [..., 0      , ji]  = p1d[...,    0 ] / float(jpi) 
     299            p2d [..., jpj-1  , ji]  = p1d[..., jpn-1] / float(jpi) 
     300    else :  
     301        for ji in np.arange ( p2d.shape[-1] ) : 
     302            p2d [..., 0      , ji]  = p1d[...,    0 ] 
     303            p2d [..., jpj-1  , ji]  = p1d[..., jpn-1] 
    294304 
    295305    if math == xr : 
    296306        p2d = xr.DataArray (p2d) 
    297307        p2d.attrs.update ( p1d.attrs ) 
    298         p2d = p2d.rename ( {p2d.dims[0]:p1d.dims[0], p2d.dims[-1]:'x', p2d.dims[-2]:'y'} ) 
     308        for idim in np.arange ( len(p1d.shape [0:-1]) ): 
     309            dim = p1d.dims[idim] 
     310            p2d = p2d.rename        ( {p2d.dims[idim]:p1d.dims[idim]}  ) 
     311            p2d = p2d.assign_coords ( {p2d.dims[idim]:p1d.coords[dim]} ) 
     312        p2d = p2d.rename ( {p2d.dims[-1]:'x', p2d.dims[-2]:'y'} ) 
     313        p2d['x'].attrs.update ( {'axis':'X'} ) 
     314        p2d['y'].attrs.update ( {'axis':'Y'} ) 
     315 
     316    if type(lon) == bool : 
     317        if lon : 
     318            lon = np.linspace ( -180, 180, jpi, endpoint=False) 
     319            #if math == xr : 
     320            #    lon = xr.DataArray ( lon, dims=('lon',), coords=(lon,) ) 
     321    if type(lat) == bool : 
     322        if lat : 
     323            lat = np.linspace ( 90, -90, jpj, endpoint=True)        
     324            #if math == xr : 
     325            #    lat = xr.DataArray ( lat, dims=('lat',), coords=(lat,) ) 
     326    if type(lon) != bool and type(lat) != bool : 
     327        if  __mmath__ (lat) == xr : 
     328            lat_name = lat.name 
     329        else : 
     330            lat_name = 'lat' 
     331        if __mmath__ (lon) == xr :  
     332            lon_name = lon.name 
     333        else : 
     334            lon_name = 'lon' 
     335        p2d = p2d.rename ( {p2d.dims[-1]:lon_name, p2d.dims[-2]:lat_name} ) 
     336        p2d = p2d.assign_coords ( {lon_name:lon, lat_name:lat} ) 
     337        p2d[lon_name].attrs.update ( { 'units':'degrees_east' , 'long_name':'Longitude', 'standard_name':'longitude', 'axis':'X' } ) 
     338        p2d[lat_name].attrs.update ( { 'units':'degrees_north', 'long_name':'Latitude' , 'standard_name':'latitude' , 'axis':'Y' }  ) 
    299339 
    300340    return p2d 
    301341 
    302342def geo2point ( p2d, cumul_poles=False, dim1d='points_physiques' ) : 
    303     '''From 2D (lat, lon) to 1D (points_phyiques) 
     343    ''' 
     344    From 2D (..., lat, lon) to 1D (..., points_phyiques) 
    304345    ''' 
    305346    math = __mmath__ (p2d) 
    306347    # 
    307     # Get the dimension 
    308  
     348    # Get the horizontal dimensions 
    309349    (jpj, jpi) = p2d.shape[-2:] 
    310  
    311     if len (p2d.shape) > 2 : 
    312         jpt = p2d.shape[0] 
     350    jpn = jpi*(jpj-2) + 2 
     351 
     352    # Get other dimensions 
     353    form1 = list(p2d.shape [0:-2])  
     354    form2 = form1 + [jpn,] 
     355         
     356    p1d = np.empty ( form2 ) 
     357    form_shape = form1 + [jpn-2,] 
     358     
     359    if math == xr : 
     360            p1d[..., 1:-1] = np.reshape ( np.float64 (p2d[..., 1:-1, :].values).ravel(), form_shape ) 
    313361    else : 
    314         jpt = -1 
    315  
    316     jpn = jpi*(jpj-2) + 2 
    317  
    318     if jpt == -1 : 
    319         p1d = np.zeros ( (jpn,) ) 
    320         p1d[1:-1] = np.float64(p2d[1:-1, :]).ravel() 
    321         p1d[ 0]   = p2d[ 0, 0] 
    322         p1d[-1]   = p2d[-1, 0] 
    323  
    324     else : 
    325         p1d = np.zeros ( (jpt, jpn) ) 
    326         if math == xr : 
    327             p1d[:, 1:-1] = np.reshape ( np.float64 (p2d[:, 1:-1, :].values).ravel(), (jpt, jpn-2) ) 
    328         else : 
    329             p1d[:, 1:-1] = np.reshape ( np.float64 (p2d[:, 1:-1, :]       ).ravel(), (jpt, jpn-2) ) 
    330         p1d[:,  0  ] = p2d[:,  0, 0] 
    331         p1d[:, -1  ] = p2d[:, -1, 0] 
     362            p1d[..., 1:-1] = np.reshape ( np.float64 (p2d[..., 1:-1, :]       ).ravel(), form_shape ) 
     363    p1d[...,  0  ] = p2d[...,  0, 0] 
     364    p1d[..., -1  ] = p2d[..., -1, 0] 
    332365 
    333366    if math == xr : 
    334367        p1d = xr.DataArray (p1d) 
    335368        p1d.attrs.update ( p2d.attrs ) 
    336         p1d = p1d.rename ( {p1d.dims[0]:p2d.dims[0], p1d.dims[-1]:dim1d} ) 
     369        if len(p2d.shape [0:-2]) > 0 :  
     370            for idim in np.arange ( len(p2d.shape [0:-2]) ): 
     371                dim=p2d.dims[idim] 
     372                p1d = p1d.rename        ( {p1d.dims[idim]:p2d.dims[idim]}   ) 
     373                p1d = p1d.assign_coords ( {p1d.dims[idim] :p2d.coords[dim].values} ) 
     374        p1d = p1d.rename ( {p1d.dims[-1]:dim1d} ) 
    337375 
    338376    if cumul_poles : 
     
    341379 
    342380    return p1d 
    343  
    344 def geo3point ( p3d, cumul_poles=False, dim1d='points_physiques' ) : 
    345     '''From 3D (lev, lat, lon) to 2D (lev, points_phyiques) 
    346     ''' 
    347     math = __mmath__ (p3d) 
    348     # 
    349     # Get the dimensions 
    350  
    351     (jpk, jpj, jpi) = p3d.shape[-3:] 
    352  
    353     if len (p3d.shape) > 3 : 
    354         jpt = p3d.shape[0] 
    355     else : 
    356         jpt = -1 
    357  
    358     jpn = jpi*(jpj-2) + 2 
    359  
    360     if jpt == -1 : 
    361         p2d = np.zeros ( (jpk, jpn,) ) 
    362         for jk in np.arange (jpk) : 
    363             p2d [jk, :] = geo2point ( p3d [jk,:,:], cumul_poles, dim1d ) 
    364     else : 
    365         p2d = np.zeros ( (jpt, jpk, jpn) ) 
    366         for jk in np.arange (jpk) : 
    367             p2d [:, jk, :] = geo2point ( p3d [:, jk,:,:], cumul_poles, dim1d  ) 
    368  
    369     if math == xr : 
    370         p2d       = xr.DataArray (p2d) 
    371         p2d.attrs.update ( p3d.attrs ) 
    372         p2d       = p2d.rename ( {p2d.dims[-1]:dim1d, p2d.dims[-2]:p3d.dims[-3]} ) 
    373  
    374     return p2d 
    375381 
    376382def geo2en (pxx, pyy, pzz, glam, gphi) : 
Note: See TracChangeset for help on using the changeset viewer.