Changeset 6265


Ignore:
Timestamp:
11/02/22 13:45:56 (18 months ago)
Author:
omamce
Message:

O.M. : WATER_BUDGET :

Add nemo and lmdz libraries.
Small changes in routines : ATM water budget still not correct
OCE water budget matches fluxes for NEMO 4.2

Location:
TOOLS/WATER_BUDGET
Files:
2 added
3 edited

Legend:

Unmodified
Added
Removed

  • TOOLS/WATER_BUDGET/ATM_waterbudget.py

    r6264 r6265  
    8989    mpi4py.rc.initialize = False 
    9090         
    91     ## Access to the nemo and lmdz module 
    92     sys.path.append ( os.path.join (subprocess.getoutput ( f'ccc_home -u p86mart' ), 'Python', 'Library')) 
    93  
    9491    ## Creates output directory 
    9592    TmpDir = os.path.join ( subprocess.getoutput ( 'ccc_home --cccscratch' ), f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) 

  • TOOLS/WATER_BUDGET/OCE_waterbudget.py

    r6264 r6265  
    9797    mpi4py.rc.initialize = False 
    9898         
    99     ## Access to the nemo module 
    100     sys.path.append ( os.path.join (subprocess.getoutput ( f'ccc_home -u p86mart' ), 'Python', 'Library') ) 
    101  
    10299    ## Creates output directory 
    103100    TmpDir = os.path.join ( subprocess.getoutput ( 'ccc_home --cccscratch' ), f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) 

  • TOOLS/WATER_BUDGET/waterbudget.py

    r6264 r6265  
    2121 
    2222## Define Experiment 
    23 if True :  
     23if False :  
    2424    JobName="TEST-CM72-SIMPLE-ROUTING.12" ; TagName="IPSLCM7" ; SpaceName="DEVT" ; ExperimentName="piControl" ; User="p86caub" ; Project="gencmip6"  
    2525    Freq = 'MO' ; YearBegin = 1970 ; YearEnd = 1979 ; PackFrequency = 10 
    … …  
    3737    ORCA = 'eORCA1.2' ; ATM = 'LMD144142' ; Routing = 'ORCHIDEE' ; NEMO=3.6 ; OCE_relax = False ; OCE_icb = False 
    3838     
    39 if False : 
     39if True : 
    4040    JobName="CM65v420-LR-SKL-pi-05" ; TagName="IPSLCM6" ; SpaceName="DEVT" ; ExperimentName="piControl" ; User="p48ethe" ; Project="gencmip6"  
    41     Freq = 'MO' ;  YearBegin = 1890 ; YearEnd = 1899 ; PackFrequency = 10 
     41    #Freq = 'MO' ;  YearBegin = 1890 ; YearEnd = 1899 ; PackFrequency = 10 
     42    Freq = 'MO' ;  YearBegin = 2340 ; YearEnd = 2349 ; PackFrequency = 10 
    4243    ORCA = 'eORCA1.2'  ; ATM = 'LMD144142' ; Routing = 'ORCHIDEE' ; NEMO=4.2 ; OCE_relax = False ; OCE_icb = False ; Coupled = True 
    4344     
    … …  
    9293    mpi4py.rc.initialize = False 
    9394         
    94     ## Access to the nemo and lmdz module 
    95     sys.path.append ( os.path.join (subprocess.getoutput ( f'ccc_home -u p86mart' ), 'Python', 'Library')) 
    96  
    9795    ## Creates output directory 
    98     TmpDir = os.path.join ( subprocess.getoutput ( 'ccc_home --cccscratch' ), f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) 
     96    #TmpDir = os.path.join ( subprocess.getoutput ( 'ccc_home --cccscratch' ), f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) 
     97    TmpDir = os.path.join ( '/ccc/scratch/cont003/drf/p86mart', f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) 
    9998 
    10099if IDRIS : 
    … …  
    147146dir_ICE_his = os.path.join ( R_SAVE, "ICE", FreqDir ) 
    148147dir_SRF_his = os.path.join ( R_SAVE, "SRF", FreqDir ) 
    149  
    150 #-- Atmosphere grid file 
    151 file_DYN_grid = os.path.join ( R_IN,  
    152148 
    153149echo ( f'The analysis relies on files from the following model output directories : ' ) 
    … …  
    355351    ATM_aire = d_ATM_his ['aire'][0] 
    356352    ATM_fsea = d_ATM_his ['fract_oce'][0] + d_ATM_his ['fract_sic'][0] 
    357     ATM_fter = d_ATM_his ['fract_ter'][0] 
     353    ATM_flnd = d_ATM_his ['fract_ter'][0] + d_ATM_his ['fract_lic'][0] 
     354     
    358355if LMDZ : 
    359     ATM_aire = lmdz.geo2point ( d_ATM_his ['aire'][0] ) 
     356    ATM_aire  = lmdz.geo2point ( d_ATM_his ['aire'][0], cumulPoles=True ) 
     357    ATM_aire2 = lmdz.geo2point ( d_ATM_his ['aire'][0], cumulPoles=False ) 
    360358    ATM_fsea = lmdz.geo2point ( d_ATM_his ['fract_oce'][0] + d_ATM_his ['fract_sic'][0] ) 
    361     ATM_fter = lmdz.geo2point ( d_ATM_his ['fract_ter'][0]  ) 
    362     ATM_aire[0]  = np.sum ( d_ATM_his ['aire'][0, 0] ) 
    363     ATM_aire[-1] = np.sum ( d_ATM_his ['aire'][0,-1] ) 
    364  
    365 SRF_aire = ATM_aire * ATM_fter 
     359    ATM_flnd = lmdz.geo2point ( d_ATM_his ['fract_ter'][0]  ) 
     360    #ATM_aire[0]  = np.sum ( d_ATM_his ['aire'][0, 0] ) 
     361    #ATM_aire[-1] = np.sum ( d_ATM_his ['aire'][0,-1] ) 
     362 
     363SRF_aire = ATM_aire * ATM_flnd 
    366364 
    367365if ICO : 
     366    # Area on icosahedron grid 
    368367    file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ATM+'_grid.nc' ) 
    369368    d_DYN_aire = xr.open_dataset ( file_DYN_aire, decode_times=False).squeeze() 
    370     d_DYN_aire = d_DYN_aire.rename( {'cell':'cell_mesh'}) 
     369    d_DYN_aire = d_DYN_aire.rename ( {'cell':'cell_mesh'} ) 
    371370    DYN_aire   = d_DYN_aire['aire'] 
     371 
     372    DYN_fsea = d_DYN_aire ['fract_oce'] + d_DYN_aire ['fract_sic'] 
     373    DYN_flnd = 1.0 - DYN_fsea 
     374     
    372375if LMDZ : 
    373376    DYN_aire = ATM_aire 
     377    DYN_fsea = ATM_fsea 
     378    DYN_flnd = ATM_flnd 
    374379 
    375380#if LMDZ : 
    … …  
    434439     
    435440if NEMO == 4.0 or NEMO == 4.2 : 
    436     ICE_ice_beg = d_ICE_beg 'v_i']   ; ICE_ice_end = d_ICE_end ['v_i'] 
     441    ICE_ice_beg = d_ICE_beg ['v_i']   ; ICE_ice_end = d_ICE_end ['v_i'] 
    437442    ICE_sno_beg = d_ICE_beg ['v_s']  ; ICE_sno_end = d_ICE_end ['v_s'] 
    438443    ICE_pnd_beg = d_ICE_beg ['v_ip'] ; ICE_pnd_end = d_ICE_end ['v_ip'] 
    … …  
    502507 
    503508echo ( '\n------------------------------------------------------------------------------------' ) 
    504 echo ( '-- LMDZ changes in stores ' ) 
     509echo ( '-- ATM changes in stores ' ) 
    505510 
    506511#-- Change in precipitable water from the atmosphere daily and monthly files 
    … …  
    514519# Surface pressure 
    515520if ICO :  
    516     ATM_ps_beg = d_DYN_beg['ps'] 
    517     ATM_ps_end = d_DYN_end['ps'] 
     521    DYN_ps_beg = d_DYN_beg['ps'] 
     522    DYN_ps_end = d_DYN_end['ps'] 
    518523     
    519524if LMDZ :  
    520     ATM_ps_beg =  lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)) ) 
    521     ATM_ps_end =  lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)) ) 
     525    DYN_ps_beg =  lmdz.geo2point ( d_DYN_beg['ps'].isel(rlonv=slice(0,-1)) ) 
     526    DYN_ps_end =  lmdz.geo2point ( d_DYN_end['ps'].isel(rlonv=slice(0,-1)) ) 
    522527     
    523528# 3D Pressure 
    524 ATM_p_beg = ATM_Ahyb + ATM_Bhyb * ATM_ps_beg 
    525 ATM_p_end = ATM_Ahyb + ATM_Bhyb * ATM_ps_end 
     529DYN_p_beg = ATM_Ahyb + ATM_Bhyb * DYN_ps_beg 
     530DYN_p_end = ATM_Ahyb + ATM_Bhyb * DYN_ps_end 
    526531 
    527532# Layer thickness 
    528 ATM_sigma_beg = ATM_p_beg[0:-1]*0. 
    529 ATM_sigma_end = ATM_p_end[0:-1]*0.  
     533DYN_sigma_beg = DYN_p_beg[0:-1]*0. 
     534DYN_sigma_end = DYN_p_end[0:-1]*0.  
    530535 
    531536for k in np.arange (klevp1-1) : 
    532     ATM_sigma_beg[k,:] = (ATM_p_beg[k,:] - ATM_p_beg[k+1,:]) / g 
    533     ATM_sigma_end[k,:] = (ATM_p_end[k,:] - ATM_p_end[k+1,:]) / g 
    534  
    535 ATM_sigma_beg = ATM_sigma_beg.rename ( {'klevp1':'sigs'} ) 
    536 ATM_sigma_end = ATM_sigma_end.rename ( {'klevp1':'sigs'} ) 
     537    DYN_sigma_beg[k,:] = (DYN_p_beg[k,:] - DYN_p_beg[k+1,:]) / g 
     538    DYN_sigma_end[k,:] = (DYN_p_end[k,:] - DYN_p_end[k+1,:]) / g 
     539 
     540DYN_sigma_beg = DYN_sigma_beg.rename ( {'klevp1':'sigs'} ) 
     541DYN_sigma_end = DYN_sigma_end.rename ( {'klevp1':'sigs'} ) 
    537542 
    538543##-- Vertical and horizontal integral, and sum of liquid, solid and vapor water phases 
    539544if LMDZ : 
    540545    try :  
    541         ATM_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2Ov'] + d_DYN_beg['H2Ol'] + d_DYN_beg['H2Oi']).isel(rlonv=slice(0,-1) ) ) 
    542         ATM_wat_end = lmdz.geo3point ( (d_DYN_end['H2Ov'] + d_DYN_end['H2Ol'] + d_DYN_end['H2Oi']).isel(rlonv=slice(0,-1) ) ) 
     546        DYN_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2Ov'] + d_DYN_beg['H2Ol'] + d_DYN_beg['H2Oi']).isel(rlonv=slice(0,-1) ) ) 
     547        DYN_wat_end = lmdz.geo3point ( (d_DYN_end['H2Ov'] + d_DYN_end['H2Ol'] + d_DYN_end['H2Oi']).isel(rlonv=slice(0,-1) ) ) 
    543548    except : 
    544         ATM_wat_beg = lmdz.geo3point ( (d_DYN_beg['H2O_g'] + d_DYN_beg['H2O_l'] + d_DYN_beg['H2O_s']).isel(rlonv=slice(0,-1) ) ) 
    545         ATM_wat_end = lmdz.geo3point ( (d_DYN_end['H2O_g'] + d_DYN_end['H2O_l'] + d_DYN_end['H2O_s']).isel(rlonv=slice(0,-1) ) ) 
     549        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) ) ) 
     550        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) ) ) 
    546551if ICO : 
    547     ATM_wat_beg = ( d_DYN_beg['q'][0:4] ).sum (dim = ['nq',] ).rename ( {'lev':'sigs'} ) 
    548     ATM_wat_end = ( d_DYN_end['q'][0:4] ).sum (dim = ['nq',] ).rename ( {'lev':'sigs'} ) 
    549      
    550 ATM_mas_wat_beg = np.sum (ATM_sigma_beg * ATM_wat_beg * ATM_aire).values.item() 
    551 ATM_mas_wat_end = np.sum (ATM_sigma_end * ATM_wat_end * ATM_aire).values.item() 
     552    DYN_wat_beg = ( d_DYN_beg['q'][0:4] ).sum (dim = ['nq',] ).rename ( {'lev':'sigs'} ) 
     553    DYN_wat_end = ( d_DYN_end['q'][0:4] ).sum (dim = ['nq',] ).rename ( {'lev':'sigs'} ) 
     554     
     555ATM_mas_wat_beg = np.sum (DYN_sigma_beg * DYN_wat_beg * DYN_aire).values.item() 
     556ATM_mas_wat_end = np.sum (DYN_sigma_end * DYN_wat_end * DYN_aire).values.item() 
    552557 
    553558dATM_mas_wat = ATM_mas_wat_end - ATM_mas_wat_beg 
    554559 
    555 echo ( '\nVariation du contenu en eau atmosphere ' ) 
    556 echo ( 'ATM_mass_beg = {:12.6e} kg - ATM_mass_end = {:12.6e} kg'.format (ATM_mas_wat_beg, ATM_mas_wat_end) ) 
     560echo ( '\nVariation du contenu en eau atmosphere (dynamique) ' ) 
     561echo ( 'ATM_mas_beg = {:12.6e} kg - ATM_mas_end = {:12.6e} kg'.format (ATM_mas_wat_beg, ATM_mas_wat_end) ) 
    557562echo ( 'dMass(atm)   = {:12.3e} kg '.format (dATM_mas_wat) ) 
    558563echo ( 'dMass(atm)   = {:12.3e} Sv '.format (dATM_mas_wat/dtime_sec*1.E-9) ) 
    … …  
    561566LIC_sno_beg = d_ATM_beg['SNOW01']*d_ATM_beg['FTER']+d_ATM_beg['SNOW02']*d_ATM_beg['FLIC']+d_ATM_beg['SNOW03']*d_ATM_beg['FOCE']+d_ATM_beg['SNOW04']*d_ATM_beg['FSIC'] 
    562567LIC_sno_end = d_ATM_end['SNOW01']*d_ATM_end['FTER']+d_ATM_end['SNOW02']*d_ATM_end['FLIC']+d_ATM_end['SNOW03']*d_ATM_end['FOCE']+d_ATM_end['SNOW04']*d_ATM_end['FSIC'] 
     568 
    563569if ICO : 
    564570   LIC_sno_beg = LIC_sno_beg.rename ( {'points_physiques':'cell_mesh'} ) 
    565571   LIC_sno_end = LIC_sno_end.rename ( {'points_physiques':'cell_mesh'} ) 
    566572 
    567 LIC_mas_sno_beg = np.sum ( LIC_sno_beg * ATM_aire ).values.item() 
    568 LIC_mas_sno_end = np.sum ( LIC_sno_end * ATM_aire ).values.item() 
     573LIC_mas_sno_beg = np.sum ( LIC_sno_beg * DYN_aire ).values.item() 
     574LIC_mas_sno_end = np.sum ( LIC_sno_end * DYN_aire ).values.item() 
    569575 
    570576dLIC_mas_sno = LIC_mas_sno_end - LIC_mas_sno_beg 
    … …  
    610616snow_end        = d_SRF_end['snow_beg']        ; snow_end        = snow_end       .where (snow_end       < 1E10, 0.) 
    611617 
    612 SRF_mas_wat_beg = tot_watveg_beg + tot_watsoil_beg + snow_beg 
    613 SRF_mas_wat_end = tot_watveg_end + tot_watsoil_end + snow_end 
     618SRF_wat_beg = tot_watveg_beg + tot_watsoil_beg + snow_beg 
     619SRF_wat_end = tot_watveg_end + tot_watsoil_end + snow_end 
    614620 
    615621#SRF_mas_wat_beg = d_SRF_beg['tot_watveg_beg']+d_SRF_beg['tot_watsoil_beg'] + d_SRF_beg['snow_beg'] 
    … …  
    623629    tot_watsoil_end = tot_watsoil_end.rename ( {'y':'points_phyiques'} ) 
    624630    snow_end        = snow_end       .rename ( {'y':'points_phyiques'} ) 
    625     SRF_mas_wat_beg = SRF_mas_wat_beg.rename ( {'y':'points_phyiques'} ) 
    626     SRF_mas_wat_end = SRF_mas_wat_end.rename ( {'y':'points_phyiques'} ) 
     631    SRF_wat_beg     = SRF_wat_beg    .rename ( {'y':'points_phyiques'} ) 
     632    SRF_wat_end     = SRF_wat_end    .rename ( {'y':'points_phyiques'} ) 
    627633if ICO : 
    628634    tot_watveg_beg  = tot_watveg_beg .rename ( {'y':'cell_mesh'} ) 
    … …  
    632638    tot_watsoil_end = tot_watsoil_end.rename ( {'y':'cell_mesh'} ) 
    633639    snow_end        = snow_end       .rename ( {'y':'cell_mesh'} ) 
    634     SRF_mas_wat_end = SRF_mas_wat_end.rename ( {'y':'cell_mesh'} ) 
    635     SRF_mas_wat_end = SRF_mas_wat_end.rename ( {'y':'cell_mesh'} )   
     640    SRF_wat_beg     = SRF_wat_beg    .rename ( {'y':'cell_mesh'} ) 
     641    SRF_wat_end     = SRF_wat_end    .rename ( {'y':'cell_mesh'} )   
    636642 
    637643SRF_mas_watveg_beg  = np.sum ( tot_watveg_beg  * SRF_aire).values.item() 
    … …  
    649655echo ( 'SRF_mas_watveg_beg   = {:12.6e} kg - SRF_mas_watveg_end   = {:12.6e} kg '.format (SRF_mas_watveg_beg , SRF_mas_watveg_end) ) 
    650656echo ( 'SRF_mas_watsoil_beg  = {:12.6e} kg - SRF_mas_watsoil_end  = {:12.6e} kg '.format (SRF_mas_watsoil_beg, SRF_mas_watsoil_end) ) 
    651 echo ( 'SRF_mas_snow_beg     = {:12.6e} kg - SRF_mas_snow_end     = {:12.6e} kg '.format (SRF_mas_snow_beg    , SRF_mas_snow_end) ) 
    652  
    653 echo ( 'dMass(watveg)  = {:12.3e} kg '.format (dSRF_mas_watveg) ) 
    654 echo ( 'dMass(watsoil) = {:12.3e} kg '.format (dSRF_mas_watsoil) ) 
    655 echo ( 'dMass(sno)     = {:12.3e} kg '.format (dSRF_mas_snow) ) 
    656  
    657  
    658 SRF_mas_wat_beg = np.sum ( SRF_mas_wat_beg * SRF_aire).values.item() 
    659 SRF_mas_wat_end = np.sum ( SRF_mas_wat_end * SRF_aire).values.item() 
     657echo ( 'SRF_mas_snow_beg     = {:12.6e} kg - SRF_mas_snow_end     = {:12.6e} kg '.format (SRF_mas_snow_beg   , SRF_mas_snow_end) ) 
     658 
     659echo ( 'dMass(watveg)  = {:12.3e} kg {:12.2e} Sv {:12.2e} m '.format (dSRF_mas_watveg , dSRF_mas_watveg /dtime_sec*1E-9, dSRF_mas_watveg /OCE_aire_tot*1E-3) ) 
     660echo ( 'dMass(watsoil) = {:12.3e} kg {:12.2e} Sv {:12.2e} m '.format (dSRF_mas_watsoil, dSRF_mas_watsoil/dtime_sec*1E-9, dSRF_mas_watsoil/OCE_aire_tot*1E-3) ) 
     661echo ( 'dMass(sno)     = {:12.3e} kg {:12.2e} Sv {:12.2e} m '.format (dSRF_mas_snow   , dSRF_mas_snow   /dtime_sec*1E-9, dSRF_mas_snow   /OCE_aire_tot*1E-3  ) ) 
     662 
     663 
     664SRF_mas_wat_beg = np.sum ( SRF_wat_beg * SRF_aire).values.item() 
     665SRF_mas_wat_end = np.sum ( SRF_wat_end * SRF_aire).values.item() 
    660666 
    661667dSRF_mas_wat = SRF_mas_wat_end - SRF_mas_wat_beg 
Note: See TracChangeset for help on using the changeset viewer.