Changeset 6265 for TOOLS/WATER_BUDGET
- Timestamp:
- 11/02/22 13:45:56 (20 months ago)
- Location:
- TOOLS/WATER_BUDGET
- Files:
-
- 2 added
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
TOOLS/WATER_BUDGET/ATM_waterbudget.py
r6264 r6265 89 89 mpi4py.rc.initialize = False 90 90 91 ## Access to the nemo and lmdz module92 sys.path.append ( os.path.join (subprocess.getoutput ( f'ccc_home -u p86mart' ), 'Python', 'Library'))93 94 91 ## Creates output directory 95 92 TmpDir = os.path.join ( subprocess.getoutput ( 'ccc_home --cccscratch' ), f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) -
TOOLS/WATER_BUDGET/OCE_waterbudget.py
r6264 r6265 97 97 mpi4py.rc.initialize = False 98 98 99 ## Access to the nemo module100 sys.path.append ( os.path.join (subprocess.getoutput ( f'ccc_home -u p86mart' ), 'Python', 'Library') )101 102 99 ## Creates output directory 103 100 TmpDir = os.path.join ( subprocess.getoutput ( 'ccc_home --cccscratch' ), f'WATER_{JobName}_{YearBegin}_{YearEnd}' ) -
TOOLS/WATER_BUDGET/waterbudget.py
r6264 r6265 21 21 22 22 ## Define Experiment 23 if True :23 if False : 24 24 JobName="TEST-CM72-SIMPLE-ROUTING.12" ; TagName="IPSLCM7" ; SpaceName="DEVT" ; ExperimentName="piControl" ; User="p86caub" ; Project="gencmip6" 25 25 Freq = 'MO' ; YearBegin = 1970 ; YearEnd = 1979 ; PackFrequency = 10 … … 37 37 ORCA = 'eORCA1.2' ; ATM = 'LMD144142' ; Routing = 'ORCHIDEE' ; NEMO=3.6 ; OCE_relax = False ; OCE_icb = False 38 38 39 if False :39 if True : 40 40 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 42 43 ORCA = 'eORCA1.2' ; ATM = 'LMD144142' ; Routing = 'ORCHIDEE' ; NEMO=4.2 ; OCE_relax = False ; OCE_icb = False ; Coupled = True 43 44 … … 92 93 mpi4py.rc.initialize = False 93 94 94 ## Access to the nemo and lmdz module95 sys.path.append ( os.path.join (subprocess.getoutput ( f'ccc_home -u p86mart' ), 'Python', 'Library'))96 97 95 ## 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}' ) 99 98 100 99 if IDRIS : … … 147 146 dir_ICE_his = os.path.join ( R_SAVE, "ICE", FreqDir ) 148 147 dir_SRF_his = os.path.join ( R_SAVE, "SRF", FreqDir ) 149 150 #-- Atmosphere grid file151 file_DYN_grid = os.path.join ( R_IN,152 148 153 149 echo ( f'The analysis relies on files from the following model output directories : ' ) … … 355 351 ATM_aire = d_ATM_his ['aire'][0] 356 352 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 358 355 if 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 ) 360 358 ATM_fsea = lmdz.geo2point ( d_ATM_his ['fract_oce'][0] + d_ATM_his ['fract_sic'][0] ) 361 ATM_f ter= 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_f ter359 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 363 SRF_aire = ATM_aire * ATM_flnd 366 364 367 365 if ICO : 366 # Area on icosahedron grid 368 367 file_DYN_aire = os.path.join ( R_IN, 'ATM', 'GRID', ATM+'_grid.nc' ) 369 368 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'} ) 371 370 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 372 375 if LMDZ : 373 376 DYN_aire = ATM_aire 377 DYN_fsea = ATM_fsea 378 DYN_flnd = ATM_flnd 374 379 375 380 #if LMDZ : … … 434 439 435 440 if 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'] 437 442 ICE_sno_beg = d_ICE_beg ['v_s'] ; ICE_sno_end = d_ICE_end ['v_s'] 438 443 ICE_pnd_beg = d_ICE_beg ['v_ip'] ; ICE_pnd_end = d_ICE_end ['v_ip'] … … 502 507 503 508 echo ( '\n------------------------------------------------------------------------------------' ) 504 echo ( '-- LMDZchanges in stores ' )509 echo ( '-- ATM changes in stores ' ) 505 510 506 511 #-- Change in precipitable water from the atmosphere daily and monthly files … … 514 519 # Surface pressure 515 520 if 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'] 518 523 519 524 if 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)) ) 522 527 523 528 # 3D Pressure 524 ATM_p_beg = ATM_Ahyb + ATM_Bhyb * ATM_ps_beg525 ATM_p_end = ATM_Ahyb + ATM_Bhyb * ATM_ps_end529 DYN_p_beg = ATM_Ahyb + ATM_Bhyb * DYN_ps_beg 530 DYN_p_end = ATM_Ahyb + ATM_Bhyb * DYN_ps_end 526 531 527 532 # Layer thickness 528 ATM_sigma_beg = ATM_p_beg[0:-1]*0.529 ATM_sigma_end = ATM_p_end[0:-1]*0.533 DYN_sigma_beg = DYN_p_beg[0:-1]*0. 534 DYN_sigma_end = DYN_p_end[0:-1]*0. 530 535 531 536 for k in np.arange (klevp1-1) : 532 ATM_sigma_beg[k,:] = (ATM_p_beg[k,:] - ATM_p_beg[k+1,:]) / g533 ATM_sigma_end[k,:] = (ATM_p_end[k,:] - ATM_p_end[k+1,:]) / g534 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 540 DYN_sigma_beg = DYN_sigma_beg.rename ( {'klevp1':'sigs'} ) 541 DYN_sigma_end = DYN_sigma_end.rename ( {'klevp1':'sigs'} ) 537 542 538 543 ##-- Vertical and horizontal integral, and sum of liquid, solid and vapor water phases 539 544 if LMDZ : 540 545 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) ) ) 543 548 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) ) ) 546 551 if 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 555 ATM_mas_wat_beg = np.sum (DYN_sigma_beg * DYN_wat_beg * DYN_aire).values.item() 556 ATM_mas_wat_end = np.sum (DYN_sigma_end * DYN_wat_end * DYN_aire).values.item() 552 557 553 558 dATM_mas_wat = ATM_mas_wat_end - ATM_mas_wat_beg 554 559 555 echo ( '\nVariation du contenu en eau atmosphere ' )556 echo ( 'ATM_mas s_beg = {:12.6e} kg - ATM_mass_end = {:12.6e} kg'.format (ATM_mas_wat_beg, ATM_mas_wat_end) )560 echo ( '\nVariation du contenu en eau atmosphere (dynamique) ' ) 561 echo ( 'ATM_mas_beg = {:12.6e} kg - ATM_mas_end = {:12.6e} kg'.format (ATM_mas_wat_beg, ATM_mas_wat_end) ) 557 562 echo ( 'dMass(atm) = {:12.3e} kg '.format (dATM_mas_wat) ) 558 563 echo ( 'dMass(atm) = {:12.3e} Sv '.format (dATM_mas_wat/dtime_sec*1.E-9) ) … … 561 566 LIC_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'] 562 567 LIC_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 563 569 if ICO : 564 570 LIC_sno_beg = LIC_sno_beg.rename ( {'points_physiques':'cell_mesh'} ) 565 571 LIC_sno_end = LIC_sno_end.rename ( {'points_physiques':'cell_mesh'} ) 566 572 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()573 LIC_mas_sno_beg = np.sum ( LIC_sno_beg * DYN_aire ).values.item() 574 LIC_mas_sno_end = np.sum ( LIC_sno_end * DYN_aire ).values.item() 569 575 570 576 dLIC_mas_sno = LIC_mas_sno_end - LIC_mas_sno_beg … … 610 616 snow_end = d_SRF_end['snow_beg'] ; snow_end = snow_end .where (snow_end < 1E10, 0.) 611 617 612 SRF_ mas_wat_beg = tot_watveg_beg + tot_watsoil_beg + snow_beg613 SRF_ mas_wat_end = tot_watveg_end + tot_watsoil_end + snow_end618 SRF_wat_beg = tot_watveg_beg + tot_watsoil_beg + snow_beg 619 SRF_wat_end = tot_watveg_end + tot_watsoil_end + snow_end 614 620 615 621 #SRF_mas_wat_beg = d_SRF_beg['tot_watveg_beg']+d_SRF_beg['tot_watsoil_beg'] + d_SRF_beg['snow_beg'] … … 623 629 tot_watsoil_end = tot_watsoil_end.rename ( {'y':'points_phyiques'} ) 624 630 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'} ) 627 633 if ICO : 628 634 tot_watveg_beg = tot_watveg_beg .rename ( {'y':'cell_mesh'} ) … … 632 638 tot_watsoil_end = tot_watsoil_end.rename ( {'y':'cell_mesh'} ) 633 639 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'} ) 636 642 637 643 SRF_mas_watveg_beg = np.sum ( tot_watveg_beg * SRF_aire).values.item() … … 649 655 echo ( 'SRF_mas_watveg_beg = {:12.6e} kg - SRF_mas_watveg_end = {:12.6e} kg '.format (SRF_mas_watveg_beg , SRF_mas_watveg_end) ) 650 656 echo ( '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 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()657 echo ( 'SRF_mas_snow_beg = {:12.6e} kg - SRF_mas_snow_end = {:12.6e} kg '.format (SRF_mas_snow_beg , SRF_mas_snow_end) ) 658 659 echo ( '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) ) 660 echo ( '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) ) 661 echo ( '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 664 SRF_mas_wat_beg = np.sum ( SRF_wat_beg * SRF_aire).values.item() 665 SRF_mas_wat_end = np.sum ( SRF_wat_end * SRF_aire).values.item() 660 666 661 667 dSRF_mas_wat = SRF_mas_wat_end - SRF_mas_wat_beg
Note: See TracChangeset
for help on using the changeset viewer.