source: TOOLS/WATER_BUDGET/waterbudget.sh

Last change on this file was 6651, checked in by omamce, 7 months ago

O.M. : TOOLS/MOSAIX - Cosmetic changes

  • Property svn:executable set to *
File size: 29.2 KB
Line 
1#!/bin/bash
2#set -v
3
4#--script to check water conservation in the IPSL coupled model
5#--author: Olivier Boucher
6#--30/11/2016
7#--updated 10/05/2019
8#--updated July 2020 => Gurvan says dOCE_vol has to be computed with rho_oce/rho_liq
9
10#--only work for NEMO 1 degree ! mods needed for NEMO025
11
12#--direxp: directory of experiment up to experiment type
13#--exp: name of experiment
14#--yearbeg: first year of budget analysis (skip first year of experiment because of issues with some diagnostics)
15#--yearend: last year of budget analysis
16#--yearinc: time resolution (in year) of output files in OUTPUT
17#--for a good precision consider at least a period of 50 or 100 years
18
19direxp="/ccc/store/cont003/gencmip6/p86caub/IGCM_OUT/IPSLCM6/DEVT/piControl/"
20exp="VALID-CM622-LR.01"
21yearbeg=1890
22yearend=1899
23yearinc=10
24
25#--set maketimeseries to true to force rewriting temporary files
26#maketimeseries=false
27maketimeseries=true
28
29host=`echo $HOSTNAME |cut -c 1-3`
30
31if [ $host == 'ada' ]; then 
32  #--IDRIS
33  #--module load cdo
34  module unload cdo
35  module load cdo/1.6.5
36  module load netcdf
37  #--output directory where to store the temporary files needed for the analysis
38  dirout=$WORKDIR
39
40elif [ $host == 'ire' ]; then 
41  #--TGCC
42  #--assume cdo is loaded
43  #--output directory where to store the temporary files needed for the analysis
44  dirout=$CCCWORKDIR
45
46else
47  echo 'Machine not supported'
48  exit
49
50fi
51
52#--name of output file
53fileout='waterbudget_'${exp}'_'${yearbeg}'_'${yearend}'.out'
54echo Output will be redirected to $fileout
55rm -f ${fileout}
56
57echo 'Dealing with ' $direxp/$exp > ${fileout}
58echo >> ${fileout}
59
60#--here I assume there are 360 or 365 days per year on average
61#--a better script would diagnose the number of days from the netcdf files
62nbdayinyr=365
63nbmthinyr=12
64
65#--model output directories
66dir_atm_day=${direxp}${exp}"/ATM/Output/DA/"
67dir_atm_mth=${direxp}${exp}"/ATM/Output/MO/"
68dir_oce_day=${direxp}${exp}"/OCE/Output/DA/"
69dir_oce_mth=${direxp}${exp}"/OCE/Output/MO/"
70dir_ice_day=${direxp}${exp}"/ICE/Output/DA/"
71dir_ice_mth=${direxp}${exp}"/ICE/Output/MO/"
72dir_srf_day=${direxp}${exp}"/SRF/Output/DA/"
73dir_srf_mth=${direxp}${exp}"/SRF/Output/MO/"
74
75echo 'The analysis relies on files from the following model output directories' >> ${fileout}
76echo $dir_atm_day >> ${fileout}
77echo $dir_atm_mth >> ${fileout}
78echo $dir_oce_day >> ${fileout}
79echo $dir_oce_mth >> ${fileout}
80echo $dir_ice_day >> ${fileout}
81echo $dir_ice_mth >> ${fileout}
82echo $dir_srf_day >> ${fileout}
83echo $dir_srf_mth >> ${fileout}
84
85#--number of years for the budget analysis
86nbyear=$((${yearend}-${yearbeg}+1))
87
88#--water density (kg/m3) ORCHIDEE
89rho_liq=1000.0
90#--ice density (kg/m3) in LIM3
91rho_ice=917.0
92#--snow density (kg/m3) in LIM3
93rho_sno=330.0
94#--ocean water density (kg/m3) in NEMO
95rho_oce=1027.
96#--average ratio of salty water to fresh water
97#fresh2salt=1.035
98#--conversion from /day to /s
99day2sec=86400.
100#--approximated conversion factor from m sea to Sv
101ssh2sv=`python -c "print 4.*3.14159*6370000.*6370000.*0.7/1.e6"`
102
103#--temporary output files
104fileoceout=dq_oce_${exp}_extracted_${yearbeg}_${yearend}_short.nc
105fileatmout=dq_atm_${exp}_extracted_${yearbeg}_${yearend}_short.nc
106filesrfout=dq_srf_${exp}_extracted_${yearbeg}_${yearend}_short.nc
107
108#---PREPARE ATMOS FILE
109if [ ${maketimeseries} = true ] || [ ! -e ${dirout}/${fileatmout} ] ; then 
110  #--delete all yearly files that may be there from a previous execution of the script
111  rm -f ${dirout}/${fileatmout} ${dirout}/dq_atm_${exp}_extracted_????.nc
112  #--loop on years
113  for year1 in $(seq $yearbeg $yearinc $yearend) 
114  do
115    year2=$((year1+yearinc-1))
116    file_atm=$dir_atm_mth$exp'_'$year1'0101_'$year2'123?_1M_histmth.nc'
117    #--this diags may not be in the time series so I extract them
118    #--some are not needed though
119    echo cdo -L yearmonmean -selname,aire,pourc_lic,wevap_ter,wevap_oce,wevap_lic,wevap_sic,wsnow_lic,wsnow_sic,wsnow_oce,wsnow_ter,wrain_lic,wrain_sic,wrain_ter,wrain_oce,wbilo_ter,wbilo_oce,wbilo_lic,wbilo_sic,evap,precip,snow,msnow,runofflic ${file_atm} ${dirout}/dq_atm_${exp}_extracted_${year1}.nc
120    cdo -L yearmonmean -selname,aire,pourc_lic,wevap_ter,wevap_oce,wevap_lic,wevap_sic,wsnow_lic,wsnow_sic,wsnow_oce,wsnow_ter,wrain_lic,wrain_sic,wrain_ter,wrain_oce,wbilo_ter,wbilo_oce,wbilo_lic,wbilo_sic,evap,precip,snow,runofflic ${file_atm} ${dirout}/dq_atm_${exp}_extracted_${year1}.nc
121  done
122  cdo -L mergetime ${dirout}/dq_atm_${exp}_extracted_????.nc ${dirout}/${fileatmout}
123  rm -f ${dirout}/dq_atm_${exp}_extracted_????.nc
124fi
125
126#---PREPARE OCE FILE
127if [ ${maketimeseries} = true ] || [ ! -e ${dirout}/${fileoceout} ] ; then
128  #--delete all yearly files that may be there from a previous execution of the script
129  rm -f ${dirout}/${fileoceout} ${dirout}/dq_oce_${exp}_extracted_????.nc
130  #--loop on years
131  for year1 in $(seq $yearbeg $yearinc $yearend)
132  do
133    year2=$((year1+yearinc-1))
134    file_oce=$dir_oce_mth$exp'_'$year1'0101_'$year2'123?_1M_grid_T.nc'
135    #--this diags may not be in the time series so I extract them
136    cdo -L yearmonmean -selname,emp_oce,emp_ice,friver,calving,iceshelf,iceberg ${file_oce} ${dirout}/dq_oce_${exp}_extracted_${year1}.nc
137  done
138  cdo -L mergetime ${dirout}/dq_oce_${exp}_extracted_????.nc ${dirout}/${fileoceout}
139  rm -f ${dirout}/dq_oce_${exp}_extracted_????.nc
140fi
141
142#---PREPARE SURFACE FILE
143if [ ${maketimeseries} = true ] || [ ! -e ${dirout}/${filesrfout} ] ; then 
144  #--delete all yearly files that may be there from a previous execution of the script
145  rm -f ${dirout}/${filesrfout} ${dirout}/dq_srf_${exp}_extracted_????.nc
146  rm -f ${dirout}/dq0_srf_${exp}_extracted_????.nc ${dirout}/dq_contfrac_${exp}_extracted_????.nc
147  #--loop on years
148  for year1 in $(seq $yearbeg $yearinc $yearend)
149  do
150    year2=$((year1+yearinc-1))
151    file_srf=$dir_srf_mth$exp'_'$year1'0101_'$year2'123?_1M_sechiba_history.nc'
152    #--this diags may not be in the time series so I extract them
153    rm -f ${dirout}/contfrac.nc ${dirout}/maxveget.nc ${dirout}/summaxveget.nc
154    cdo -L selname,Contfrac ${file_srf} ${dirout}/contfrac.nc
155    cdo -L selname,maxvegetfrac ${file_srf} ${dirout}/maxveget.nc
156    cdo -L vertsum ${dirout}/maxveget.nc ${dirout}/summaxveget.nc
157    cdo -L yearmonmean ${dirout}/summaxveget.nc ${dirout}/summaxvegetmean.nc
158    cdo -L yearmonmean -selname,Areas,evap,evapnu,evapflo,inter,transpir,riverflow,coastalflow,rain,snowf,subli,snow,vegetfrac,DelSoilMoist_daily,DelIntercept_daily,DelSWE_daily,delstreamr_daily,delfastr_daily,delslowr_daily,delfloodr_daily,delpondr_daily,dellakevol_daily ${file_srf} ${dirout}/dq0_srf_${exp}_extracted_${year1}.nc
159    cdo -L merge ${dirout}/summaxvegetmean.nc ${dirout}/contfrac.nc ${dirout}/dq0_srf_${exp}_extracted_${year1}.nc ${dirout}/dq_srf_${exp}_extracted_${year1}.nc
160  done
161  cdo -L mergetime ${dirout}/dq_srf_${exp}_extracted_????.nc ${dirout}/${filesrfout}
162  rm -f ${dirout}/dq0_srf_${exp}_extracted_????.nc ${dirout}/dq_srf_${exp}_extracted_????.nc
163  rm -f ${dirout}/contfrac.nc ${dirout}/maxveget.nc ${dirout}/summaxveget.nc
164fi
165
166#---------------------------------------------------------------------------------
167echo >> ${fileout}
168echo 'NOTE all numbers are in Sv (1 Sv = 10^6 m3/s) of freshwater' >> ${fileout}
169echo 'NOTE 0.001 Sv = 8.8 mm sea level rise per century' >> ${fileout}
170echo >> ${fileout}
171echo 'We are dealing with years from' $yearbeg 'to' $yearend ', or' $nbyear 'years'>> ${fileout}
172if [ $nbyear -lt 30 ] ; then
173  echo 'WARNING: at least 30 years are needed for a good accuracy because of some approximations made' >> ${fileout}
174else
175  echo 'This is good, the longer the more acurate because of some approximations made' >> ${fileout}
176fi
177echo >> ${fileout}
178echo 'TER=land OCE=ocean SIC=sea-ice LIC=land ice' >> ${fileout}
179echo >> ${fileout}
180#
181#--using daily data when possible, otherwise monthly
182file_atm_daybeg=$dir_atm_day$exp'_'$yearbeg'0101_????123?_1D_histday.nc'
183file_atm_dayend=$dir_atm_day$exp'_????0101_'$yearend'123?_1D_histday.nc'
184file_atm_mthbeg=$dir_atm_mth$exp'_'$yearbeg'0101_????123?_1M_histmth.nc'
185file_atm_mthend=$dir_atm_mth$exp'_????0101_'$yearend'123?_1M_histmth.nc'
186#
187file_srf_daybeg=$dir_srf_day$exp'_'$yearbeg'0101_????123?_1D_sechiba_history.nc'
188file_srf_dayend=$dir_srf_day$exp'_????0101_'$yearend'123?_1D_sechiba_history.nc'
189file_srf_mthbeg=$dir_srf_mth$exp'_'$yearbeg'0101_????123?_1M_sechiba_history.nc'
190file_srf_mthend=$dir_srf_mth$exp'_????0101_'$yearend'123?_1M_sechiba_history.nc'
191#
192file_oce_daybeg=$dir_oce_day$exp'_'$yearbeg'0101_????123?_1D_scalar.nc'
193file_oce_dayend=$dir_oce_day$exp'_????0101_'$yearend'123?_1D_scalar.nc'
194file_oce_mthbeg=$dir_oce_mth$exp'_'$yearbeg'0101_????123?_1M_scalar.nc'
195file_oce_mthend=$dir_oce_mth$exp'_????0101_'$yearend'123?_1M_scalar.nc'
196#
197file_ice_daybeg=$dir_ice_day$exp'_'$yearbeg'0101_????123?_1D_icemod.nc'
198file_ice_dayend=$dir_ice_day$exp'_????0101_'$yearend'123?_1D_icemod.nc'
199file_ice_mthbeg=$dir_ice_mth$exp'_'$yearbeg'0101_????123?_1M_icemod.nc'
200file_ice_mthend=$dir_ice_mth$exp'_????0101_'$yearend'123?_1M_icemod.nc'
201#
202#---------------------------------------------------------------------------------
203echo '------------------------------------------------------------------------------------' >> ${fileout}
204echo '--NEMO change in stores (for the records)' >> ${fileout}
205#--note that the total number of days of the run should be diagnosed rather than assumed
206#--here the result is in Sv
207#
208#--change in ocean volume in freshwater equivalent
209if [ -e ${file_oce_daybeg} ]; then
210  maxtimestep=`ncdump -h ${file_oce_dayend} | grep currently | grep -o -E '[0-9]+'`
211  OCE_vol_daybeg=`cdo -L outputf,%12.10g,10 -seltimestep,1 -selvar,scvoltot ${file_oce_daybeg}`
212  OCE_vol_dayend=`cdo -L outputf,%12.10g,10 -seltimestep,${maxtimestep} -selvar,scvoltot ${file_oce_dayend}`
213  echo 'OCE vol beg =' ${OCE_vol_daybeg} >> ${fileout}
214  echo 'OCE vol end =' ${OCE_vol_dayend} >> ${fileout}
215  dOCE_vol=`python -c "print  ${rho_oce}/${rho_liq}*(${OCE_vol_dayend}-${OCE_vol_daybeg})/${nbyear}/${nbdayinyr}/${day2sec}/1.e6"`
216else
217  maxtimestep=`ncdump -h ${file_oce_mthend} | grep currently | grep -o -E '[0-9]+'`
218  OCE_vol_mthbeg=`cdo -L outputf,%12.10g,10 -seltimestep,1 -selvar,scvoltot ${file_oce_mthbeg}`
219  OCE_vol_mthend=`cdo -L outputf,%12.10g,10 -seltimestep,${maxtimestep} -selvar,scvoltot ${file_oce_mthend}`
220  echo 'OCE vol beg =' ${OCE_vol_mthbeg} >> ${fileout}
221  echo 'OCE vol end =' ${OCE_vol_mthend} >> ${fileout}
222  dOCE_vol=`python -c "print  ${rho_oce}/${rho_liq}*(${OCE_vol_mthend}-${OCE_vol_mthbeg})/${nbyear}/${nbdayinyr}/${day2sec}/1.e6"`
223fi
224echo 'These estimates of d OCE (converted to Sv) are approximate'  >> ${fileout}
225echo 'd OCE vol (Sv) =' ${dOCE_vol} >> ${fileout}
226#--
227#--computing sea-ice and snow-on-sea-ice volume
228#--it would be better to compute from daily files but for now only monthly files available
229#--here we get a volume of ice (m3 of ice) by multiplying depth with surface
230if [ -e ${file_ice_daybeg} ]; then
231  echo  >> ${fileout}
232else
233  maxtimestep=`ncdump -h ${file_ice_mthend} | grep currently | grep -o -E '[0-9]+'`
234  SIC_volice_mthbeg=`cdo -L outputf,%12.10g,10 -fldsum -mul -seltimestep,1 -selvar,sivolu ${file_ice_mthbeg} -gridarea ${file_ice_mthbeg}`
235  SIC_volice_mthend=`cdo -L outputf,%12.10g,10 -fldsum -mul -seltimestep,${maxtimestep} -selvar,sivolu ${file_ice_mthend} -gridarea ${file_ice_mthend}`
236  #--converting to freswater volume
237  #--note that the total number of days of the run should be diagnosed rather than assumed
238  dSIC_volice=`python -c "print ${rho_ice}/${rho_liq}*(${SIC_volice_mthend}-${SIC_volice_mthbeg})/${nbyear}/${nbdayinyr}/${day2sec}/1.e6"`
239  #--same for snow on sea-ice (m3 on snow)
240  SIC_volsno_mthbeg=`cdo -L outputf,%12.10g,10 -fldsum -mul -seltimestep,1 -selvar,snvolu ${file_ice_mthbeg} -gridarea ${file_ice_mthbeg}`
241  SIC_volsno_mthend=`cdo -L outputf,%12.10g,10 -fldsum -mul -seltimestep,${maxtimestep} -selvar,snvolu ${file_ice_mthend} -gridarea ${file_ice_mthend}`
242  #--converting to freswater volume
243  #--so we get change in sea-ice and snow in Sv
244  dSIC_volsno=`python -c "print ${rho_sno}/${rho_liq}*(${SIC_volsno_mthend}-${SIC_volsno_mthbeg})/${nbyear}/${nbdayinyr}/${day2sec}/1.e6"`
245fi
246echo 'd SIC ice (Sv)  =' ${dSIC_volice} >> ${fileout}
247echo 'd SIC snow (Sv) =' ${dSIC_volsno} >> ${fileout}
248echo '------------------------------------------------------------------------------------' >> ${fileout}
249echo '--NEMO fluxes (for the records)'  >> ${fileout}
250#
251file=$dirout'/'$fileoceout
252echo EMP_OCE cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -selindexbox,2,361,1,331 -mul -timavg -selname,emp_oce ${file} -gridarea ${file}
253EMP_OCE=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -selindexbox,2,361,1,331 -mul -timavg -selname,emp_oce ${file} -gridarea ${file}`
254EMP_SIC=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -selindexbox,2,361,1,331 -mul -timavg -selname,emp_ice ${file} -gridarea ${file}`
255ICEBERG=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -selindexbox,2,361,1,331 -mul -timavg -selname,iceberg ${file} -gridarea ${file}`
256ICESHELF=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -selindexbox,2,361,1,331 -mul -timavg -selname,iceshelf ${file} -gridarea ${file}`
257CALVING=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -selindexbox,2,361,1,331 -mul -timavg -selname,calving ${file} -gridarea ${file}`
258FRIVER=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -selindexbox,2,361,1,331 -mul -timavg -selname,friver ${file} -gridarea ${file}`
259echo 'EMP_OCE=' $EMP_OCE >> ${fileout}
260echo 'EMP_OCE + CALVING=' `python -c "print $EMP_OCE + $CALVING"` >> ${fileout}
261echo 'EMP_SIC=' $EMP_SIC >> ${fileout}
262echo 'EMP_OCE + EMP_SIC + CALVING =' `python -c "print $EMP_OCE + $EMP_SIC + $CALVING"` >> ${fileout}
263echo 'ICEBERG =' $ICEBERG >> ${fileout}
264echo 'ICESHELF=' $ICESHELF >> ${fileout}
265echo 'CALVING =' $CALVING >> ${fileout}
266echo 'FRIVER  =' $FRIVER >> ${fileout}
267echo '------------------------------------------------------------------------------------' >> ${fileout}
268CREATED_OCESIC=`python -c "print ${dOCE_vol}+${dSIC_volice}+${dSIC_volsno}+${EMP_OCE}+${EMP_SIC}-${FRIVER}-${ICEBERG}-${ICESHELF}"`
269echo 'WATER CREATED BY OCE+SIC is computed as' >> ${fileout}
270echo 'dOCE_vol + dSIC_volice + dSIC_volsno + EMP_OCE + EMP_SIC - FRIVER - ICEBERG - ICESHELF' >> ${fileout}
271echo 'and should be as small as possible, ideally ~ 0.001 Sv over a 50 to 100 year period'  >> ${fileout}
272echo 'WATER CREATED BY OCE+SIC=' ${CREATED_OCESIC} >> ${fileout}
273#
274echo '------------------------------------------------------------------------------------' >> ${fileout}
275echo '--LMDZ changes in stores (for the records)' >> ${fileout}
276echo 'should be less than 0.001 Sv for prw and even smaller for cloud water'  >> ${fileout}
277#--change in precipitable water from the atmosphere daily and monthly files
278#--compute sum weighted by gridcell area (kg/m2) then convert to Sv
279if [ -e ${file_atm_daybeg} ]; then
280  maxtimestep=`ncdump -h ${file_atm_dayend} | grep currently | grep -o -E '[0-9]+'`
281  ATM_prw_daybeg=`cdo -L outputf,%12.10g,10 -fldsum -mul -seltimestep,1 -selvar,prw ${file_atm_daybeg} -gridarea ${file_atm_daybeg}`
282  ATM_prw_dayend=`cdo -L outputf,%12.10g,10 -fldsum -mul -seltimestep,${maxtimestep} -selvar,prw ${file_atm_dayend} -gridarea ${file_atm_dayend}`
283  dATM_prw=`python -c "print (${ATM_prw_dayend}-${ATM_prw_daybeg})/${nbyear}/${nbdayinyr}/${day2sec}/1.e9"`
284  #--now same with cloud liquid water (small term)
285  ATM_prlw_daybeg=`cdo -L outputf,%12.10g,10 -fldsum -mul -seltimestep,1 -selvar,prlw ${file_atm_daybeg} -gridarea ${file_atm_daybeg}`
286  ATM_prlw_dayend=`cdo -L outputf,%12.10g,10 -fldsum -mul -seltimestep,${maxtimestep} -selvar,prlw ${file_atm_dayend} -gridarea ${file_atm_dayend}`
287  dATM_prlw=`python -c "print (${ATM_prlw_dayend}-${ATM_prlw_daybeg})/${nbyear}/${nbdayinyr}/${day2sec}/1.e9"`
288  #--now same with cloud ice water (small term)
289  ATM_prsw_daybeg=`cdo -L outputf,%12.10g,10 -fldsum -mul -seltimestep,1 -selvar,prsw ${file_atm_daybeg} -gridarea ${file_atm_daybeg}`
290  ATM_prsw_dayend=`cdo -L outputf,%12.10g,10 -fldsum -mul -seltimestep,${maxtimestep} -selvar,prsw ${file_atm_dayend} -gridarea ${file_atm_dayend}`
291  dATM_prsw=`python -c "print (${ATM_prsw_dayend}-${ATM_prsw_daybeg})/${nbyear}/${nbdayinyr}/${day2sec}/1.e9"`
292else
293  maxtimestep=`ncdump -h ${file_atm_mthend} | grep currently | grep -o -E '[0-9]+'`
294  ATM_prw_mthbeg=`cdo -L outputf,%12.10g,10 -fldsum -mul -seltimestep,1 -selvar,prw ${file_atm_mthbeg} -gridarea ${file_atm_mthbeg}`
295  ATM_prw_mthend=`cdo -L outputf,%12.10g,10 -fldsum -mul -seltimestep,${maxtimestep} -selvar,prw ${file_atm_mthend} -gridarea ${file_atm_mthend}`
296  dATM_prw=`python -c "print (${ATM_prw_mthend}-${ATM_prw_mthbeg})/${nbyear}/${nbdayinyr}/${day2sec}/1.e9"`
297  #
298  ATM_prlw_mthbeg=`cdo -L outputf,%12.10g,10 -fldsum -mul -seltimestep,1 -selvar,prlw ${file_atm_mthbeg} -gridarea ${file_atm_mthbeg}`
299  ATM_prlw_mthend=`cdo -L outputf,%12.10g,10 -fldsum -mul -seltimestep,${maxtimestep} -selvar,prlw ${file_atm_mthend} -gridarea ${file_atm_mthend}`
300  dATM_prlw=`python -c "print (${ATM_prlw_mthend}-${ATM_prlw_mthbeg})/${nbyear}/${nbdayinyr}/${day2sec}/1.e9"`
301  #
302  ATM_prsw_mthbeg=`cdo -L outputf,%12.10g,10 -fldsum -mul -seltimestep,1 -selvar,prsw ${file_atm_mthbeg} -gridarea ${file_atm_mthbeg}`
303  ATM_prsw_mthend=`cdo -L outputf,%12.10g,10 -fldsum -mul -seltimestep,${maxtimestep} -selvar,prsw ${file_atm_mthend} -gridarea ${file_atm_mthend}`
304  dATM_prsw=`python -c "print (${ATM_prsw_mthend}-${ATM_prsw_mthbeg})/${nbyear}/${nbdayinyr}/${day2sec}/1.e9"`
305  #
306  ATM_snowlic_mthbeg=`cdo -L outputf,%12.10g,10 -fldsum -mul -seltimestep,1 -expr,'msnow*pourc_lic/100.' ${file_atm_mthbeg} -gridarea ${file_atm_mthbeg}`
307  ATM_snowlic_mthend=`cdo -L outputf,%12.10g,10 -fldsum -mul -seltimestep,${maxtimestep} -expr,'msnow*pourc_lic/100.' ${file_atm_mthend} -gridarea ${file_atm_mthend}`
308  dATM_snowlic=`python -c "print (${ATM_snowlic_mthend}-${ATM_snowlic_mthbeg})/${nbyear}/${nbdayinyr}/${day2sec}/1.e9"`
309fi
310echo 'd ATM from prw  (Sv) =' ${dATM_prw} >> ${fileout}
311echo 'd ATM from prlw (Sv) =' ${dATM_prlw} >> ${fileout}
312echo 'd ATM from prsw (Sv) =' ${dATM_prsw} >> ${fileout}
313dATM_tot=`python -c "print ${dATM_prw}+${dATM_prlw}+${dATM_prsw}"`
314#
315echo '------------------------------------------------------------------------------------' >> ${fileout}
316echo '--LMDZ fluxes (for the records)' >> ${fileout}
317#--defining atmospheric file
318file=${dirout}/${fileatmout}
319#
320#--atmospheric fluxes (converted to Sv)
321PREC=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -expr,'toto=precip' ${file} -gridarea ${file}`
322RAIN_OCE=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -expr,'toto=wrain_oce' ${file} -gridarea ${file}`
323RAIN_SIC=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -expr,'toto=wrain_sic' ${file} -gridarea ${file}`
324RAIN_TER=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -expr,'toto=wrain_ter' ${file} -gridarea ${file}`
325RAIN_LIC=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -expr,'toto=wrain_lic' ${file} -gridarea ${file}`
326SNOW=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -expr,'toto=snow' ${file} -gridarea ${file}`
327SNOW_OCE=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -expr,'toto=wsnow_oce' ${file} -gridarea ${file}`
328SNOW_SIC=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -expr,'toto=wsnow_sic' ${file} -gridarea ${file}`
329SNOW_TER=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -expr,'toto=wsnow_ter' ${file} -gridarea ${file}`
330SNOW_LIC=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -expr,'toto=wsnow_lic' ${file} -gridarea ${file}`
331EVAP=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -expr,'toto=evap' ${file} -gridarea ${file}`
332EVAP_OCE=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -expr,'toto=wevap_oce' ${file} -gridarea ${file}`
333EVAP_SIC=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -expr,'toto=wevap_sic' ${file} -gridarea ${file}`
334EVAP_TER=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -expr,'toto=wevap_ter' ${file} -gridarea ${file}`
335EVAP_LIC=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -expr,'toto=wevap_lic' ${file} -gridarea ${file}`
336#--better to compute precip-evap then sum (rather than the other way around)
337PRECMINUSEVAP=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -expr,'toto=precip-evap' ${file} -gridarea ${file}`
338echo 'PRECIP - EVAP should be < 0.0001 Sv'  >> ${fileout}
339echo 'PRECIP - EVAP (Sv) =' ${PREC} '-' ${EVAP} '=' ${PRECMINUSEVAP} >> ${fileout}
340echo 'Next three lines are optional' >> ${fileout} 
341echo 'Numbers come in this order OCE + SIC + TER + LIC = TOTAL' >> ${fileout}
342echo 'RAIN by tile and total (Sv) =' ${RAIN_OCE} '+' ${RAIN_SIC} '+' ${RAIN_TER} '+' ${RAIN_LIC} '=' \
343                                    `python -c "print ${RAIN_OCE}+${RAIN_SIC}+${RAIN_TER}+${RAIN_LIC}"` >> ${fileout}
344echo 'SNOW by tile and total (Sv) =' ${SNOW_OCE} '+' ${SNOW_SIC} '+' ${SNOW_TER} '+' ${SNOW_LIC} '=' \
345                                    `python -c "print ${SNOW_OCE}+${SNOW_SIC}+${SNOW_TER}+${SNOW_LIC}"` >> ${fileout}
346echo 'EVAP by tile and total (Sv) =' ${EVAP_OCE} '+' ${EVAP_SIC} '+' ${EVAP_TER} '+' ${EVAP_LIC} '=' \
347                                    `python -c "print ${EVAP_OCE}+${EVAP_SIC}+${EVAP_TER}+${EVAP_LIC}"` >> ${fileout}
348#--water budget by surface type as seen from atmosphere
349echo 'Net water budget (evap-precip) by land surface type in LMDz' >> ${fileout}
350BILO_TER=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -selvar,wbilo_ter ${file} -gridarea ${file}`
351echo 'BILO TER (Sv)=' $BILO_TER >> ${fileout}
352BILO_OCE=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -selvar,wbilo_oce ${file} -gridarea ${file}`
353echo 'BILO OCE (Sv)=' $BILO_OCE >> ${fileout}
354BILO_SIC=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -selvar,wbilo_sic ${file} -gridarea ${file}`
355echo 'BILO SIC (Sv)=' $BILO_SIC >> ${fileout}
356BILO_LIC=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -selvar,wbilo_lic ${file} -gridarea ${file}`
357echo 'BILO LIC (Sv)=' $BILO_LIC >> ${fileout}
358BILO_tot=`python -c "print $BILO_TER+$BILO_OCE+$BILO_SIC+$BILO_LIC"` >> ${fileout}
359echo 'The BILO tot diag below should be similar (+/-5e-6) to PRECIP-EVAP above but with opposite sign' >> ${fileout}
360echo 'BILO tot (Sv)=' $BILO_tot >> ${fileout}
361echo 'The BILO OCE+SIC diag below should be the same (within <0.001 Sv) as EMP_OCE+EMP_ICE+CALVING above ' >> ${fileout}
362echo 'BILO OCE+SIC (Sv)=' `python -c "print $BILO_OCE+$BILO_SIC"` >> ${fileout}
363#--here we add PRLW and PRSW when they exist but these are small
364CREATED_ATM=`python -c "print ${PRECMINUSEVAP}+${dATM_prw}+${dATM_prlw}+${dATM_prsw}"`
365echo '------------------------------------------------------------------------------------' >> ${fileout}
366echo 'Water created in LMDZ should be < 0.001 Sv'  >> ${fileout}
367echo 'WATER_CREATED_IN_ATM (Sv)=' ${CREATED_ATM} >> ${fileout}
368echo '------------------------------------------------------------------------------------' >> ${fileout}
369#
370echo '--LIC change in store and fluxes' >> ${fileout}
371#--water going to the ocean
372RUNOFFLIC=`cdo -L outputf,%12.8g,8 -divc,1.e9 -fldsum -mul -timavg -expr,'toto=runofflic*pourc_lic/100.' ${file} -gridarea ${file}`
373echo 'Atmospheric calving can be different from that reaching the ocean because there is a buffer'  >> ${fileout}
374echo 'CALVING atm (Sv)=' $RUNOFFLIC >> ${fileout}
375dLIC_tot=`python -c "print -${BILO_LIC}-${RUNOFFLIC}"`
376echo 'd LIC (Sv) = ' ${dLIC_tot} >> ${fileout}
377#--by construction
378CREATED_LIC=0
379echo 'd LIC diagnosed from msnow (Sv) = ' ${dATM_snowlic} >> ${fileout}
380#
381#--------------------------------------------------------------------------------------------------------
382#--SECHIBA budgets
383file=${dirout}/${filesrfout}
384echo '------------------------------------------------------------------------------------' >> ${fileout}
385echo '--SECHIBA change in stores' >> ${fileout}
386#--changes in reservoirs - from tendencies
387dSoilHum_in_Sv=`cdo -L outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=(maxvegetfrac*DelSoilMoist_daily)*Contfrac' ${file} -gridarea ${file}`
388dInterce_in_Sv=`cdo -L outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=DelIntercept_daily*Contfrac' ${file} -gridarea ${file}`
389dSWE_in_Sv=`cdo -L outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=DelSWE_daily*Contfrac' ${file} -gridarea ${file}`
390dStream_in_Sv=`cdo -L outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=delstreamr_daily*Contfrac' ${file} -gridarea ${file}`
391dFastR_in_Sv=`cdo -L outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=delfastr_daily*Contfrac' ${file} -gridarea ${file}`
392dSlowR_in_Sv=`cdo -L outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=delslowr_daily*Contfrac' ${file} -gridarea ${file}`
393dFlood_in_Sv=`cdo -L outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=delfloodr_daily*Contfrac' ${file} -gridarea ${file}`
394dPond_in_Sv=`cdo -L outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=delpondr_daily*Contfrac' ${file} -gridarea ${file}`
395dLake_in_Sv=`cdo -L outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=dellakevol_daily*Contfrac' ${file} -gridarea ${file}`
396echo 'dSTOCK Soil Hum (Sv)  =' $dSoilHum_in_Sv >> ${fileout}
397echo 'dSTOCK Intercept (Sv) =' $dInterce_in_Sv >> ${fileout}
398echo 'dSTOCK Snow (Sv)      =' $dSWE_in_Sv >> ${fileout}
399echo 'dSTOCK Stream (Sv)    =' $dStream_in_Sv >> ${fileout}
400echo 'dSTOCK FastR (Sv)     =' $dFastR_in_Sv >> ${fileout}
401echo 'dSTOCK SlowR (Sv)     =' $dSlowR_in_Sv >> ${fileout}
402echo 'dSTOCK Lake (Sv)      =' $dLake_in_Sv >> ${fileout}
403echo 'dSTOCK Pond (Sv)      =' $dPond_in_Sv >> ${fileout}
404echo 'dSTOCK Flood (Sv)     =' $dFlood_in_Sv >> ${fileout}
405dSRF_tot=`python -c "print $dSoilHum_in_Sv+$dInterce_in_Sv+$dSWE_in_Sv+$dStream_in_Sv+$dFastR_in_Sv+$dSlowR_in_Sv+$dLake_in_Sv+$dPond_in_Sv+$dFlood_in_Sv"`
406echo 'dSTOCK Total (Sv)     =' ${dSRF_tot} >> ${fileout}
407#
408echo '------------------------------------------------------------------------------------' >> ${fileout}
409echo '--SECHIBA fluxes' >> ${fileout}
410EVAP_SRF=`cdo -L outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -selvar,evap ${file} -gridarea ${file}`
411echo 'EVAP SRF (Sv)=' ${EVAP_SRF} >> ${fileout}
412EVAPNU_SRF=`cdo -L outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -selvar,evapnu ${file} -gridarea ${file}`
413EVAPFL_SRF=`cdo -L outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -selvar,evapflo ${file} -gridarea ${file}`
414SUBLIM_SRF=`cdo -L outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -selvar,subli ${file} -gridarea ${file}`
415TRANSP_SRF=`cdo -L outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -vertsum -selvar,transpir ${file} -gridarea ${file}`
416INTERC_SRF=`cdo -L outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -vertsum -selvar,inter ${file} -gridarea ${file}`
417echo 'This diagnoses how different is EVAP from the sum of known terms in SECHIBA' >> ${fileout}
418echo 'The terms are evap bare soil, evap flooding plains, sublimation, transpiration, interception' >> ${fileout}
419echo 'The difference with EVAP SRF should be small <0.0001 Sv' >> ${fileout}
420echo 'EVAP TERMS SRF=' $EVAPNU_SRF '+' $EVAPFL_SRF '+' $SUBLIM_SRF '+' $TRANSP_SRF '+' $INTERC_SRF  >> ${fileout}
421echo 'EVAP TOTAL SRF=' `python -c "print ${EVAPNU_SRF}+${EVAPFL_SRF}+${SUBLIM_SRF}+${TRANSP_SRF}+${INTERC_SRF}"` >> ${fileout}
422#
423echo 'RAINFALL SRF=' `cdo -L outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=rain*Contfrac' ${file} -gridarea ${file}`  >> ${fileout}
424echo 'SNOWFALL SRF=' `cdo -L outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=snowf*Contfrac' ${file} -gridarea ${file}` >> ${fileout}
425IN_SRF=`cdo -L outputf,%12.8g,8 -divc,1.e9 -divc,${day2sec} -fldsum -mul -timavg -expr,'toto=(rain+snowf-evap)*Contfrac' ${file} -gridarea ${file}` 
426echo 'The SECH NET SRF (precip - evap over land) should be the same as BILO TER above but opposite sign' >> ${fileout}
427echo 'SECH NET SRF=' ${IN_SRF} >> ${fileout}
428#
429#--flow to ocean (Sv)
430FLOW=`cdo -L outputf,%12.8g,8 -divc,1.e6 -fldsum -timavg -expr,'toto=coastalflow+riverflow' ${file}`
431echo 'FLOW to the ocean should the same (<0.0001 Sv) as FRIVER above' >> ${fileout}
432echo 'FLOW to ocean (Sv)=' ${FLOW} >> ${fileout}
433echo '------------------------------------------------------------------------------------' >> ${fileout}
434#--budget of SRF alltogether
435CREATED_SRF=`python -c "print ${FLOW}-${IN_SRF}+${dSRF_tot}"`
436echo 'Created water by ORCHIDEE should be < 0.0001 Sv' >> ${fileout}
437echo 'WATER CREATED BY ORCHIDEE=' ${FLOW} '-' ${IN_SRF} '+' ${dSRF_tot} '=' ${CREATED_SRF} >> ${fileout}
438#
439echo '------------------------------------------------------------------------------------' >> ${fileout}
440echo 'SUMMARY' >> ${fileout}
441echo 'd ATM (Sv) =' ${dATM_tot}         >> ${fileout}
442echo 'd LIC (Sv) =' ${dLIC_tot}         >> ${fileout}
443echo 'd SRF (Sv) =' ${dSRF_tot}         >> ${fileout}
444echo 'd OCE (Sv) =' ${dOCE_vol}         >> ${fileout}
445echo 'd SIC ice  (Sv) =' ${dSIC_volice} >> ${fileout}
446echo 'd SIC snow (Sv) =' ${dSIC_volsno} >> ${fileout}
447CREATED_TOT1=`python -c "print ${dATM_tot}+${dLIC_tot}+${dSRF_tot}+${dOCE_vol}+${dSIC_volice}+${dSIC_volsno}"`
448CREATED_TOT2=`python -c "print ${CREATED_OCESIC}+${CREATED_SRF}+${CREATED_ATM}+${CREATED_LIC}"`
449echo 'Two different diagnostics of created water rate in the coupled model' >> ${fileout}
450echo 'should be < 0.003 Sv in the current state of affairs' >> ${fileout}
451echo 'WATER CREATED TOTAL (Sv) =' ${CREATED_TOT1} >> ${fileout}
452echo 'WATER CREATED TOTAL (Sv) =' ${CREATED_TOT2} >> ${fileout}
453echo '------------------------------------------------------------------------------------' >> ${fileout}
Note: See TracBrowser for help on using the repository browser.