import cdms2 import numpy as N from constantes import * from math import * import copy #------------------------------------------------------------------------------------------- # function write_nc def write_nc(listvar,listvar2,listvar3,val_lon,val_lat,year_start,year_end,name,name_path): timear = N.ma.arange(len(listvar[0]), dtype=N.float) timear = timear*1800 time = cdms2.createAxis(timear, id='tstep') time.units = "seconds since "+str(year_start)+"-01-01 00:00:00" time.origin = str(year_start)+"-01-01 00:00:00" time.calendar = "gregorian" zlevelar = N.ma.arange(1, dtype=N.float) zlevelar = zlevelar+2 zlevel = cdms2.createAxis(zlevelar, id='lev') zlevel.units = "m" zlevel.long_name = "Vertical levels" latar = N.ma.arange(1, dtype=N.float) latar = latar+val_lat lat = cdms2.createAxis(latar, id='lat') lat.units = "degrees_north" lat.long_name = "Latitude" lonar = N.ma.arange(1, dtype=N.float) lonar = latar+val_lon lon = cdms2.createAxis(lonar, id='lon') lon.units = "degrees_east" lon.long_name = "Longitude" g = cdms2.createUniformGrid(val_lat, 1, 0, val_lon, 1, 0) tair_ar = N.array(listvar[id_tair]) tair_ar.shape = (len(listvar[id_tair]),1,1,1) tair = cdms2.createVariable(tair_ar, id='Tair', axes = (time,zlevel,lat,lon)) tair.long_name = "Near surface air temperature" tair.units = "K" tair.missing_value = "1.e+20f" tair_qc_ar = N.array(listvar2[id_tair]) tair_qc_ar.shape = (len(listvar2[id_tair]),1,1,1) tair_qc = cdms2.createVariable(tair_qc_ar, id='tair_qc', axes = (time,zlevel,lat,lon)) tair_qc.long_name = "Quality check for Tair" tair_qc.units = "-" tair_qc.missing_value = "1.e+20f" psurf_ar = N.array(listvar[id_psurf]) psurf_ar.shape = (len(listvar[id_psurf]),1,1) psurf = cdms2.createVariable(psurf_ar, id='PSurf', axes = (time,lat,lon)) psurf.long_name = "Surface pressure" psurf.units = "Pa" psurf.missing_value = "1.e+20f" psurf_qc_ar = N.array(listvar2[id_psurf]) psurf_qc_ar.shape = (len(listvar2[id_psurf]),1,1) psurf_qc = cdms2.createVariable(psurf_qc_ar, id='psurf_qc', axes = (time,lat,lon)) psurf_qc.long_name = "quality check for Psurf" psurf_qc.units = "-" psurf_qc.missing_value = "1.e+20f" qair_ar = N.array(listvar[id_qair]) qair_ar.shape = (len(listvar[id_qair]),1,1,1) qair = cdms2.createVariable(qair_ar, id='Qair', axes = (time,zlevel,lat,lon)) qair.long_name = "Near Surface specific humidity" qair.units = "kg/kg" qair.missing_value = "1.e+20f" temp=listvar2[id_qair]*listvar2[id_tair]*listvar2[id_psurf] qair_qc_ar = N.array(temp) qair_qc_ar.shape = (len(temp),1,1,1) qair_qc = cdms2.createVariable(qair_qc_ar, id='qair_qc', axes = (time,zlevel,lat,lon)) qair_qc.long_name = "Quality check for qair" qair_qc.units = "-" qair_qc.missing_value = "1.e+20f" wind_n_ar = N.array(listvar[id_ws]) wind_n_ar.shape = (len(listvar[id_ws]),1,1,1) wind_n = cdms2.createVariable(wind_n_ar, id='Wind_N', axes = (time,zlevel,lat,lon)) wind_n.long_name = "Near surface northward wind component" wind_n.units = "m/s" wind_n.missing_value = "1.e+20f" wind_n_qc_ar= N.array(listvar2[id_ws]) wind_n_qc_ar.shape = (len(listvar2[id_ws]),1,1,1) wind_n_qc = cdms2.createVariable(wind_n_qc_ar, id='wind_n_qc', axes = (time,zlevel,lat,lon)) wind_n_qc.long_name = "Quality check for wind_n" wind_n_qc.units = "-" wind_n_qc.missing_value = "1.e+20f" wind_e_ar = N.zeros(len(listvar[id_ws])) wind_e_ar.shape = (len(listvar[id_ws]),1,1,1) wind_e = cdms2.createVariable(wind_e_ar, id='Wind_E', axes = (time,zlevel,lat,lon)) wind_e.long_name = "Near surface eastward wind component" wind_e.units = "m/s" wind_e.missing_value = "1.e+20f" wind_e_qc_ar = N.zeros(len(listvar2[id_ws])) wind_e_qc_ar.shape = (len(listvar2[id_ws]),1,1,1) wind_e_qc = cdms2.createVariable(wind_e_qc_ar, id='Wind_e_qc', axes = (time,zlevel,lat,lon)) wind_e_qc.long_name = "Quality check for wind_e" wind_e_qc.units = "-" wind_e_qc.missing_value = "1.e+20f" rainf_ar = N.array(listvar[id_rain]) rainf_ar.shape = (len(listvar[id_rain]),1,1) rainf = cdms2.createVariable(rainf_ar, id='Rainf', axes = (time,lat,lon)) rainf.long_name = "Rainfall rate" rainf.units = "kg/m^2s" rainf.missing_value = "1.e+20f" rainf_qc_ar = N.array(listvar2[id_precip]) rainf_qc_ar.shape = (len(listvar2[id_precip]),1,1) rainf_qc = cdms2.createVariable(rainf_qc_ar, id='Rainf_qc', axes = (time,lat,lon)) rainf_qc.long_name = "Quality check for Rainf" rainf_qc.units = "-" rainf_qc.missing_value = "1.e+20f" snowf_ar = N.array(listvar[id_snow]) snowf_ar.shape = (len(listvar[id_snow]),1,1) snowf = cdms2.createVariable(snowf_ar, id='Snowf', axes = (time,lat,lon)) snowf.long_name = "Snowfall rate" snowf.units = "kg/m^2s" snowf.missing_value = "1.e+20f" snowf_qc_ar = N.array(listvar2[id_precip]) snowf_qc_ar.shape = (len(listvar2[id_precip]),1,1) snowf_qc = cdms2.createVariable(snowf_qc_ar, id='Snowf_qc', axes = (time,lat,lon)) snowf_qc.long_name = "Quality check for Snowf" snowf_qc.units = "-" snowf_qc.missing_value = "1.e+20f" swdown_ar = N.array(listvar[id_swdown]) swdown_ar.shape = (len(listvar[id_swdown]),1,1) swdown = cdms2.createVariable(swdown_ar, id='SWdown', axes = (time,lat,lon)) swdown.long_name = "Surface incident shortwave radiation" swdown.units = "W/m^2" swdown.missing_value = "1.e+20f" swdown_qc_ar = N.array(listvar2[id_swdown]) swdown_qc_ar.shape = (len(listvar2[id_swdown]),1,1) swdown_qc = cdms2.createVariable(swdown_qc_ar, id='SWdown_qc', axes = (time,lat,lon)) swdown_qc.long_name = "Quality check for Swdown" swdown_qc.units = "-" swdown_qc.missing_value = "1.e+20f" lwdown_ar = N.array(listvar[id_lwdown]) lwdown_ar.shape = (len(listvar[id_lwdown]),1,1) lwdown = cdms2.createVariable(lwdown_ar, id='LWdown', axes = (time,lat,lon)) lwdown.long_name = "Surface incident longwave radiation" lwdown.units = "W/m^2" lwdown.missing_value = "1.e+20f" lwdown_qc_ar = N.array(listvar2[id_lwdown]) lwdown_qc_ar.shape = (len(listvar2[id_lwdown]),1,1) lwdown_qc = cdms2.createVariable(lwdown_qc_ar, id='LWdown_qc', axes = (time,lat,lon)) lwdown_qc.long_name = "Quality check for Lwdown" lwdown_qc.units = "-" lwdown_qc.missing_value = "1.e+20f" Ts1_f_ar = N.array(listvar3[id_Ts1_f]) Ts1_f_ar.shape = (len(listvar3[id_Ts1_f]),1,1) Ts1_f = cdms2.createVariable(Ts1_f_ar, id='TSOIL', axes = (time,lat,lon)) Ts1_f.long_name = "Soil Temperature" Ts1_f.units = "degC" Ts1_f.missing_value = -9999. Rh_ar = N.array(listvar3[id_Rh]) Rh_ar.shape = (len(listvar3[id_Rh]),1,1) Rh = cdms2.createVariable(Rh_ar, id='RH', axes = (time,lat,lon)) Rh.long_name = "Relative Humidity" Rh.units = "%" Rh.missing_value = -9999. NEE_f_ar = N.array(listvar3[id_NEE_f]) NEE_f_ar.shape = (len(listvar3[id_NEE_f]),1,1) NEE_f = cdms2.createVariable(NEE_f_ar, id='NEE', axes = (time,lat,lon)) NEE_f.long_name = "Net Ecosystem Exchange" NEE_f.units = "gC/m2/tstep" NEE_f.missing_value = -9999. H_f_ar = N.array(listvar3[id_H_f]) H_f_ar.shape = (len(listvar3[id_H_f]),1,1) H_f = cdms2.createVariable(H_f_ar, id='Fh', axes = (time,lat,lon)) H_f.long_name = "Sensible Heat Flux" H_f.units = "W/m2" H_f.missing_value = -9999. LE_f_ar = N.array(listvar3[id_LE_f]) LE_f_ar.shape = (len(listvar3[id_LE_f]),1,1) LE_f = cdms2.createVariable(LE_f_ar, id='Fle', axes = (time,lat,lon)) LE_f.long_name = "Latent Heat Flux" LE_f.units = "W/m2" LE_f.missing_value = -9999. Epot_f_ar = N.array(listvar3[id_Epot_f]) Epot_f_ar.shape = (len(listvar3[id_Epot_f]),1,1) Epot_f = cdms2.createVariable(Epot_f_ar, id='ET', axes = (time,lat,lon)) Epot_f.long_name = "Evapotranspiration" Epot_f.units = "mm/tstep" Epot_f.missing_value = -9999. SWC1_f_ar = N.array(listvar3[id_SWC1_f]) SWC1_f_ar.shape = (len(listvar3[id_SWC1_f]),1,1) SWC1_f = cdms2.createVariable(SWC1_f_ar, id='SWC', axes = (time,lat,lon)) SWC1_f.long_name = "Soil Water Content" SWC1_f.units = "VOL%" SWC1_f.missing_value = -9999. gsurf_f_ar = N.array(listvar3[id_gsurf_f]) gsurf_f_ar.shape = (len(listvar3[id_gsurf_f]),1,1) gsurf_f = cdms2.createVariable(gsurf_f_ar, id='GS', axes = (time,lat,lon)) gsurf_f.long_name = "Canopy conductance" gsurf_f.units = "mm/s" gsurf_f.missing_value = -9999. GPP_f_ar = N.array(listvar3[id_GPP_f]) GPP_f_ar.shape = (len(listvar3[id_GPP_f]),1,1) GPP_f = cdms2.createVariable(GPP_f_ar, id='GPP', axes = (time,lat,lon)) GPP_f.long_name = "Gross Primary Production" GPP_f.units = "gC/m2/tstep" GPP_f.missing_value = -9999. Reco_ar = N.array(listvar3[id_Reco]) Reco_ar.shape = (len(listvar3[id_Reco]),1,1) Reco = cdms2.createVariable(Reco_ar, id='Reco', axes = (time,lat,lon)) Reco.long_name = "Ecosystem Respiration" Reco.units = "gC/m2/tstep" Reco.missing_value = -9999. f = cdms2.open(name_path+name+'_'+str(year_start)+'-'+str(year_end)+'.nc', 'w') f.write(tair,axes=None) f.write(qair,axes=None) f.write(psurf,axes=None) f.write(wind_n,axes=None) f.write(wind_e,axes=None) f.write(rainf,axes=None) f.write(snowf,axes=None) f.write(swdown,axes=None) f.write(lwdown,axes=None) f.write(tair_qc,axes=None) f.write(qair_qc,axes=None) f.write(psurf_qc,axes=None) f.write(wind_n_qc,axes=None) f.write(wind_e_qc,axes=None) f.write(rainf_qc,axes=None) f.write(snowf_qc,axes=None) f.write(swdown_qc,axes=None) f.write(lwdown_qc,axes=None) f.write(Ts1_f,axes=None) f.write(Rh,axes=None) f.write(NEE_f,axes=None) f.write(H_f,axes=None) f.write(LE_f,axes=None) f.write(Epot_f,axes=None) f.write(SWC1_f,axes=None) f.write(gsurf_f,axes=None) f.write(GPP_f,axes=None) f.write(Reco,axes=None) f.close() # end function write_nc #------------------------------------------------------------------------------------------- def prepare_for_orchidee(weather_gapfill): weather_out=copy.deepcopy(weather_gapfill) # Conversion from VPD (Vapour Pressure Deficit, Pa) to Qair (Near Surface Specific Humidity, kg/kg) for i in range(len(weather_out[id_vpd])): # Magnus Tetens (Murray, 1967) http://cires.colorado.edu/~voemel/vp.html if(weather_out[id_tair][i]<273.15): eg=c74*exp(c70*(weather_out[id_tair][i]-273.15)/(weather_out[id_tair][i]-c71)) else: eg=c74*exp(c72*(weather_out[id_tair][i]-273.15)/(weather_out[id_tair][i]-c73)) weather_out[id_vpd][i]=(eg-weather_out[id_vpd][i])*c75/(weather_out[id_psurf][i]-(eg-weather_out[id_vpd][i])*c76) weather_out.append(N.zeros(len(weather_out[id_precip]))) for i in range(len(weather_out[id_precip])): if(weather_out[id_tair][i]<273.15): weather_out[id_snow][i]=weather_out[id_precip][i] weather_out[id_rain][i]=0. return weather_out #-------------------------------------------------------------------------------------------