= Conservation de l'évaporation : où est l'erreur ? = 

== Exemple avec le MR025 - script cdo == 

Expérience : /ccc/cont003/home/gencmip6/oboucher/IPSLCM6.2/modipsl/config/IPSLCM6/CM61-MR025-pd-04

Sorties : /ccc/scratch/cont003/gencmip6/oboucher/IGCM_OUT/IPSLCM6/TEST/pdControl/CM61-MR025-pd-04/

Détail : 1 mois de simu 

Script cdo : on pondère par tmaskutil et l'aire des mailles pour l'océan. On pondère par l'aire des mailles et pourc_oce+pourc_sic pour l'atmosphère (moyenne mensuelle mais la somme des deux est un invariant). 

{{{ 

exp="CM61-MR025-pd-04"
yr1=1850
yr2=1850
dir=/ccc/scratch/cont003/gencmip6/oboucher/IGCM_OUT/IPSLCM6/TEST/pdControl/${exp}
fileatm=${dir}/ATM/Output/MO/${exp}_${yr1}0101_${yr2}0131_1M_histmth.nc
fileevap1=${dir}/CPL/Output/MO/${exp}_${yr1}0101_${yr2}0131_COTOTEVA.nc
fileevap2=${dir}/CPL/Output/MO/${exp}_${yr1}0101_${yr2}0131_OTotEvap.nc
filerain1=${dir}/CPL/Output/MO/${exp}_${yr1}0101_${yr2}0131_COTOTRAI.nc
filerain2=${dir}/CPL/Output/MO/${exp}_${yr1}0101_${yr2}0131_OTotRain.nc
filesnow1=${dir}/CPL/Output/MO/${exp}_${yr1}0101_${yr2}0131_COTOTSNO.nc
filesnow2=${dir}/CPL/Output/MO/${exp}_${yr1}0101_${yr2}0131_OTotSnow.nc
fileoce=${dir}/OCE/Output/MO/${exp}_${yr1}0101_${yr2}0131_1M_grid_T.nc
maskoce=/ccc/work/cont003/igcmg/igcmg/IGCM/OCE/NEMO/eORCA025.4/GRIDS/eORCA025_mesh_mask.nc

sumevap1=`cdo outputf,%12.10g,10 -divc,1.e9 -fldsum -mul -seltimestep,1/10 -selname,COTOTEVA ${fileevap1} -seltimestep,1 -expr,'toto=aire/100.*(pourc_oce+pourc_sic)' ${fileatm}`
sumevap2=`cdo outputf,%12.10g,10 -divc,1.e9 -fldsum -mul -mul -seltimestep,1/10 -selname,OTotEvap ${fileevap2} -selvar,tmaskutil ${maskoce} -gridarea ${fileoce}`
sumrain1=`cdo outputf,%12.10g,10 -divc,1.e9 -fldsum -mul -seltimestep,1/10 -selname,COTOTRAI ${filerain1} -seltimestep,1 -expr,'toto=aire/100.*(pourc_oce+pourc_sic)' ${fileatm}`
sumrain2=`cdo outputf,%12.10g,10 -divc,1.e9 -fldsum -mul -mul -seltimestep,1/10 -selname,OTotRain ${filerain2} -selvar,tmaskutil ${maskoce} -gridarea ${fileoce}`
sumsnow1=`cdo outputf,%12.10g,10 -divc,1.e9 -fldsum -mul -seltimestep,1/10 -selname,COTOTSNO ${filesnow1} -seltimestep,1 -expr,'toto=aire/100.*(pourc_oce+pourc_sic)' ${fileatm}`
sumsnow2=`cdo outputf,%12.10g,10 -divc,1.e9 -fldsum -mul -mul -seltimestep,1/10 -selname,OTotSnow ${filesnow2} -selvar,tmaskutil ${maskoce} -gridarea ${fileoce}`

echo 'evap in=' $sumevap1
echo 'evap out=' $sumevap2
echo 'rain in=' $sumrain1
echo 'rain out=' $sumrain2
echo 'snow in=' $sumsnow1
echo 'snow out=' $sumsnow2

 }}}

Résultats : 
{{{ 
evap in= 16.58377264 16.5842436 16.73348733 16.70182978 16.66627324 16.63719809 16.59929387 16.55878592 16.52697635 16.4860888
evap out= 16.55351427 16.55350019 16.70269787 16.67088149 16.63496352 16.60587634 16.56772969 16.52697626 16.49513478 16.45415434
rain in= 12.76418216 11.79569265 11.90206466 12.29730812 12.41094174 12.47357032 12.53198387 12.55976865 12.5649098 12.58444113
rain out= 12.74270543 11.77684503 11.88402995 12.27952764 12.395365 12.45955016 12.51887749 12.5477753 12.55303025 12.57239723
snow in= 0.6307305996 0.613739933 0.6188579544 0.6443643258 0.6604737915 0.6693078788 0.6729603416 0.6815849762 0.684608096 0.6879445765
snow out= 0.6216787337 0.6042128991 0.608490469 0.6332710421 0.6475660592 0.6557678674 0.6592571085 0.6672780528 0.6703978744 0.6741682585
}}}

Les écarts sont très importants : 4e décimale et 0.2% pour rain et evap, 2e décimale et 2% pour snow.  
Bien sûr, le biais est moindre pour P+S-E car il y a compensation d'erreur. 

== Exemple avec le VLR - script python == 

Expérience : /ccc/work/cont003/ra5424/p25khod/IPSLCM6.5_work_15DEC/modipsl/config/IPSLCM6/CM65-VLR-pd-tn-01-test

Fichier poids: /ccc/work/cont003/gencmip6/lebasn/INPUT/VLR/IPSLCM6_9695_ORCA2.3/grids_ORCA2.3xLMD9695_MOSAIX_v1.nc

Sorties : /ccc/store/cont003/gencmip6/p25khod/IGCM_OUT/IPSLCM6/DEVT/pdControl/CM65-VLR-pd-tn-01-test

Détails : 12 mois de simu, 18500101_18501231, cpl_old_calving=no, bug fixé dans LMDZ pour le calcul du flux cumulé par bande de latitude 

Problème : le flux de évap qui repart du coupleur diffère un peu trop du flux de évap qui arrive, alors que les flux de pluie et neige sont très bien conservés 

Fichier d'aires NEMO : ceux du fichier de sortie grid_T de NEMO 

Fichier masque NEMO : /ccc/cont003/home/lmd/oboucher/BILAN_EAU/maskutil_T.nc

Diagnostic 1 - script python 

{{{
import xarray as xr
import numpy as np

exp="CM65-VLR-pd-tn-01-test"
dir="/ccc/store/cont003/gencmip6/p25khod/IGCM_OUT/IPSLCM6/DEVT/pdControl/"+exp+"/"

file_atm=dir+"ATM/Output/MO/"+exp+"_18500101_18501231_1M_histmth.nc"

file_calving_1=dir+"CPL/Output/MO/"+exp+"_18500101_18501231_COCALVIN.nc"
file_calving_2=dir+"CPL/Output/MO/"+exp+"_18500101_18501231_OCalving.nc"

file_evap_1=dir+"CPL/Output/MO/"+exp+"_18500101_18501231_COTOTEVA.nc"
file_evap_2=dir+"CPL/Output/MO/"+exp+"_18500101_18501231_OTotEvap.nc"

file_rain_1=dir+"CPL/Output/MO/"+exp+"_18500101_18501231_COTOTRAI.nc"
file_rain_2=dir+"CPL/Output/MO/"+exp+"_18500101_18501231_OTotRain.nc"

file_snow_1=dir+"CPL/Output/MO/"+exp+"_18500101_18501231_COTOTSNO.nc"
file_snow_2=dir+"CPL/Output/MO/"+exp+"_18500101_18501231_OTotSnow.nc"

fileoce=dir+"OCE/Output/MO/"+exp+"_18500101_18501231_1M_grid_T.nc"
maskoce="/ccc/cont003/home/lmd/oboucher/BILAN_EAU/maskutil_T.nc"

#--ouverture des xarray
xr_atm=xr.open_dataset(file_atm)
xr_calving_1=xr.open_dataset(file_calving_1)
xr_calving_2=xr.open_dataset(file_calving_2)
xr_evap_1=xr.open_dataset(file_evap_1)
xr_evap_2=xr.open_dataset(file_evap_2)
xr_rain_1=xr.open_dataset(file_rain_1)
xr_rain_2=xr.open_dataset(file_rain_2)
xr_snow_1=xr.open_dataset(file_snow_1)
xr_snow_2=xr.open_dataset(file_snow_2)
xr_ocean=xr.open_dataset(fileoce)
xr_mask=xr.open_dataset(maskoce)

#--extraction des array numpy - conversion en longdouble pour les sommes
atm_aire0=np.longdouble(xr_atm['aire'].values)
atm_oce=np.longdouble(xr_atm['pourc_oce'].values)
atm_sic=np.longdouble(xr_atm['pourc_sic'].values)

wbilo_oce=np.longdouble(xr_atm['wbilo_oce'].values)
wbilo_sic=np.longdouble(xr_atm['wbilo_sic'].values)

atm_aire=atm_aire0*(atm_oce+atm_sic)/100.

calving_1=np.longdouble(xr_calving_1['COCALVIN'].values)
calving_2=np.longdouble(xr_calving_2['OCalving'].values)

evap_1=np.longdouble(xr_evap_1['COTOTEVA'].values)
evap_2=np.longdouble(xr_evap_2['OTotEvap'].values)

rain_1=np.longdouble(xr_rain_1['COTOTRAI'].values)
rain_2=np.longdouble(xr_rain_2['OTotRain'].values)

snow_1=np.longdouble(xr_snow_1['COTOTSNO'].values)
snow_2=np.longdouble(xr_snow_2['OTotSnow'].values)

oce_area=np.longdouble(xr_ocean['area'].values)

mask=np.longdouble(xr_mask['maskutil_T'][:,:].values)

#--totaux sur la grille - conversion en Sv
tot_calving_1=np.sum(calving_1,axis=(1,2))/1.e9
tot_calving_2=np.sum(calving_2*oce_area*mask,axis=(1,2))/1.e9
print('calving in  Sv =',tot_calving_1[0:10])
print('calving out Sv =',tot_calving_2[0:10])

}}}

Résultats pour le calving : 
{{{

('calving in  Sv =', array([ 0.0,  0.025110938,  0.07871894,  0.060956317,  0.066686016,
        0.032709929,  0.049833877,  0.033137503,  0.032376687,  0.059238586], dtype=float128))
('calving out Sv =', array([ 0.0,  0.025110938,  0.07871894,  0.060956317,  0.066686016,
        0.032709929,  0.049833878,  0.033137503,  0.032376687,  0.059238586], dtype=float128))

}}}

Tout va bien, c'est la même chose en input et en output 

{{{
tot_evap_1=np.sum(evap_1*atm_aire[0,:,:],axis=(1,2))/1.e9
tot_evap_2=np.sum(evap_2*oce_area*mask,axis=(1,2))/1.e9
tot_rain_1=np.sum(rain_1*atm_aire[0,:,:],axis=(1,2))/1.e9
tot_rain_2=np.sum(rain_2*oce_area*mask,axis=(1,2))/1.e9
tot_snow_1=np.sum(snow_1*atm_aire[0,:,:],axis=(1,2))/1.e9
tot_snow_2=np.sum(snow_2*oce_area*mask,axis=(1,2))/1.e9

print('oceanic area from grid_T file')
print('rain in  Sv =',tot_rain_1[0:10])
print('rain out Sv =',tot_rain_2[0:10])
print('rain diff mSv=',(tot_rain_1[1:10]-tot_rain_2[1:10])*1000.)
print('snow in  Sv =',tot_snow_1[0:10])
print('snow out Sv =',tot_snow_2[0:10])
print('snow diff mSv=',(tot_snow_1[1:10]-tot_snow_2[1:10])*1000.)
print('evap in  Sv =',tot_evap_1[0:10])
print('evap out Sv =',tot_evap_2[0:10])
print('evap diff mSv=',(tot_evap_1[1:10]-tot_evap_2[1:10])*1000.)

}}}

Résultats : 
{{{
('rain in  Sv =', array([ 0.0,  3.4585411,  8.8186534,  11.318086,  11.570535,  11.836867,
        11.357751,  11.295551,  11.242278,  11.49349], dtype=float128))
('rain out Sv =', array([ 0.0,  3.4585411,  8.8186534,  11.317966,  11.570457,  11.836863,
        11.35775,  11.295479,  11.242256,  11.493486], dtype=float128))
('rain diff mSv=', array([-3.9859023e-08,  4.4476907e-05,  0.11972505,  0.078397034,
        0.0045189078,  0.0016220226,  0.071562863,  0.022124081,
        0.0034941963], dtype=float128))
('snow in  Sv =', array([ 0.0,  0.43405517,  0.48368972,  0.41251261,  0.36736733,
        0.39909007,  0.46146092,  0.3871201,  0.33174664,  0.45059974], dtype=float128))
('snow out Sv =', array([ 0.0,  0.43405517,  0.48368972,  0.41251261,  0.36736733,
        0.39909007,  0.46146092,  0.38712011,  0.33174664,  0.45059974], dtype=float128))
('snow diff mSv=', array([-5.991792e-07, -2.3985875e-07,  4.5716792e-07,  3.592597e-07,
        1.0149842e-08, -1.356738e-08, -3.040086e-07,  3.1490164e-07,
       -3.3566433e-07], dtype=float128))
('evap in  Sv =', array([ 0.0,  12.638699,  12.227281,  13.689618,  14.567473,  14.940661,
        14.904491,  14.177904,  13.691065,  13.653616], dtype=float128))
('evap out Sv =', array([ 0.0,  12.638585,  12.226845,  13.689247,  14.566987,  14.940099,
        14.903804,  14.177131,  13.690155,  13.652478], dtype=float128))
('evap diff mSv=', array([ 0.11373978,  0.43656301,  0.37129745,  0.48537692,  0.5611987,
        0.68672948,  0.77286054,  0.91008484,  1.1375948], dtype=float128))

}}}

Snow est très bien conservé, Rain pas trop mal (<0.1 mSv) mais ce n'est pas le cas de Evap avec des différences jusqu'à 1 mSv.