source: input_nc_orchidee.py @ 123

Last change on this file since 123 was 123, checked in by nicolas.vuichard, 13 years ago

import .py files
File size: 11.0 KB
RevLine 
[123]1import cdms2
2import numpy as N
3from constantes import *
4from math import *
5import copy
6
7#-------------------------------------------------------------------------------------------
8# function write_nc
9def write_nc(listvar,listvar2,listvar3,val_lon,val_lat,year_start,year_end,name,name_path):
10   timear = N.ma.arange(len(listvar[0]), dtype=N.float)
11   timear = timear*1800
12   time = cdms2.createAxis(timear, id='tstep')
13   time.units = "seconds since "+str(year_start)+"-01-01 00:00:00"
14   time.origin = str(year_start)+"-01-01 00:00:00"
15   time.calendar = "gregorian"
16
17   zlevelar = N.ma.arange(1, dtype=N.float)
18   zlevelar = zlevelar+2
19   zlevel = cdms2.createAxis(zlevelar, id='lev')
20   zlevel.units = "m"
21   zlevel.long_name = "Vertical levels"
22   
23   latar = N.ma.arange(1, dtype=N.float)
24   latar = latar+val_lat
25   lat = cdms2.createAxis(latar, id='lat')
26   lat.units = "degrees_north"
27   lat.long_name = "Latitude"
28
29   
30   lonar = N.ma.arange(1, dtype=N.float)
31   lonar = latar+val_lon
32   lon = cdms2.createAxis(lonar, id='lon')
33   lon.units = "degrees_east"
34   lon.long_name = "Longitude"
35   
36   g = cdms2.createUniformGrid(val_lat, 1, 0, val_lon, 1, 0)
37
38   tair_ar = N.array(listvar[id_tair])
39   tair_ar.shape = (len(listvar[id_tair]),1,1,1)
40   tair = cdms2.createVariable(tair_ar, id='Tair', axes = (time,zlevel,lat,lon))
41   tair.long_name = "Near surface air temperature" 
42   tair.units = "K" 
43   tair.missing_value = "1.e+20f"
44
45   tair_qc_ar = N.array(listvar2[id_tair])
46   tair_qc_ar.shape = (len(listvar2[id_tair]),1,1,1)
47   tair_qc = cdms2.createVariable(tair_qc_ar, id='tair_qc', axes = (time,zlevel,lat,lon))
48   tair_qc.long_name = "Quality check for Tair" 
49   tair_qc.units = "-" 
50   tair_qc.missing_value = "1.e+20f"
51
52   psurf_ar = N.array(listvar[id_psurf])
53   psurf_ar.shape = (len(listvar[id_psurf]),1,1)
54   psurf = cdms2.createVariable(psurf_ar, id='PSurf', axes = (time,lat,lon))
55   psurf.long_name = "Surface pressure" 
56   psurf.units = "Pa" 
57   psurf.missing_value = "1.e+20f"
58
59   psurf_qc_ar = N.array(listvar2[id_psurf])
60   psurf_qc_ar.shape = (len(listvar2[id_psurf]),1,1)
61   psurf_qc = cdms2.createVariable(psurf_qc_ar, id='psurf_qc', axes = (time,lat,lon))
62   psurf_qc.long_name = "quality check for Psurf" 
63   psurf_qc.units = "-" 
64   psurf_qc.missing_value = "1.e+20f"
65
66   qair_ar = N.array(listvar[id_qair])
67   qair_ar.shape = (len(listvar[id_qair]),1,1,1)
68   qair = cdms2.createVariable(qair_ar, id='Qair', axes = (time,zlevel,lat,lon))
69   qair.long_name = "Near Surface specific humidity" 
70   qair.units = "kg/kg" 
71   qair.missing_value = "1.e+20f"
72
73   temp=listvar2[id_qair]*listvar2[id_tair]*listvar2[id_psurf]
74   qair_qc_ar = N.array(temp)
75   qair_qc_ar.shape = (len(temp),1,1,1)
76   qair_qc = cdms2.createVariable(qair_qc_ar, id='qair_qc', axes = (time,zlevel,lat,lon))
77   qair_qc.long_name = "Quality check for qair" 
78   qair_qc.units = "-" 
79   qair_qc.missing_value = "1.e+20f"
80
81   wind_n_ar = N.array(listvar[id_ws])
82   wind_n_ar.shape = (len(listvar[id_ws]),1,1,1)
83   wind_n = cdms2.createVariable(wind_n_ar, id='Wind_N', axes = (time,zlevel,lat,lon))
84   wind_n.long_name = "Near surface northward wind component" 
85   wind_n.units = "m/s" 
86   wind_n.missing_value = "1.e+20f"
87
88   wind_n_qc_ar= N.array(listvar2[id_ws])
89   wind_n_qc_ar.shape = (len(listvar2[id_ws]),1,1,1)
90   wind_n_qc = cdms2.createVariable(wind_n_qc_ar, id='wind_n_qc', axes = (time,zlevel,lat,lon))
91   wind_n_qc.long_name = "Quality check for wind_n" 
92   wind_n_qc.units = "-" 
93   wind_n_qc.missing_value = "1.e+20f"
94
95   wind_e_ar = N.zeros(len(listvar[id_ws]))
96   wind_e_ar.shape = (len(listvar[id_ws]),1,1,1)
97   wind_e = cdms2.createVariable(wind_e_ar, id='Wind_E', axes = (time,zlevel,lat,lon))
98   wind_e.long_name = "Near surface eastward wind component" 
99   wind_e.units = "m/s" 
100   wind_e.missing_value = "1.e+20f"
101
102   wind_e_qc_ar = N.zeros(len(listvar2[id_ws]))
103   wind_e_qc_ar.shape = (len(listvar2[id_ws]),1,1,1)
104   wind_e_qc = cdms2.createVariable(wind_e_qc_ar, id='Wind_e_qc', axes = (time,zlevel,lat,lon))
105   wind_e_qc.long_name = "Quality check for wind_e" 
106   wind_e_qc.units = "-" 
107   wind_e_qc.missing_value = "1.e+20f"
108
109   rainf_ar = N.array(listvar[id_rain])
110   rainf_ar.shape = (len(listvar[id_rain]),1,1)
111   rainf = cdms2.createVariable(rainf_ar, id='Rainf', axes = (time,lat,lon))
112   rainf.long_name = "Rainfall rate" 
113   rainf.units = "kg/m^2s" 
114   rainf.missing_value = "1.e+20f"
115
116   rainf_qc_ar = N.array(listvar2[id_precip])
117   rainf_qc_ar.shape = (len(listvar2[id_precip]),1,1)
118   rainf_qc = cdms2.createVariable(rainf_qc_ar, id='Rainf_qc', axes = (time,lat,lon))
119   rainf_qc.long_name = "Quality check for Rainf" 
120   rainf_qc.units = "-" 
121   rainf_qc.missing_value = "1.e+20f"
122
123   snowf_ar = N.array(listvar[id_snow])
124   snowf_ar.shape = (len(listvar[id_snow]),1,1)
125   snowf = cdms2.createVariable(snowf_ar, id='Snowf', axes = (time,lat,lon))
126   snowf.long_name = "Snowfall rate" 
127   snowf.units = "kg/m^2s" 
128   snowf.missing_value = "1.e+20f"
129
130   snowf_qc_ar = N.array(listvar2[id_precip])
131   snowf_qc_ar.shape = (len(listvar2[id_precip]),1,1)
132   snowf_qc = cdms2.createVariable(snowf_qc_ar, id='Snowf_qc', axes = (time,lat,lon))
133   snowf_qc.long_name = "Quality check for Snowf" 
134   snowf_qc.units = "-" 
135   snowf_qc.missing_value = "1.e+20f"
136
137   swdown_ar = N.array(listvar[id_swdown])
138   swdown_ar.shape = (len(listvar[id_swdown]),1,1)
139   swdown = cdms2.createVariable(swdown_ar, id='SWdown', axes = (time,lat,lon))
140   swdown.long_name = "Surface incident shortwave radiation" 
141   swdown.units = "W/m^2" 
142   swdown.missing_value = "1.e+20f"
143
144   swdown_qc_ar = N.array(listvar2[id_swdown])
145   swdown_qc_ar.shape = (len(listvar2[id_swdown]),1,1)
146   swdown_qc = cdms2.createVariable(swdown_qc_ar, id='SWdown_qc', axes = (time,lat,lon))
147   swdown_qc.long_name = "Quality check for Swdown" 
148   swdown_qc.units = "-" 
149   swdown_qc.missing_value = "1.e+20f"
150
151   lwdown_ar = N.array(listvar[id_lwdown])
152   lwdown_ar.shape = (len(listvar[id_lwdown]),1,1)
153   lwdown = cdms2.createVariable(lwdown_ar, id='LWdown', axes = (time,lat,lon))
154   lwdown.long_name = "Surface incident longwave radiation" 
155   lwdown.units = "W/m^2" 
156   lwdown.missing_value = "1.e+20f"
157
158   lwdown_qc_ar = N.array(listvar2[id_lwdown])
159   lwdown_qc_ar.shape = (len(listvar2[id_lwdown]),1,1)
160   lwdown_qc = cdms2.createVariable(lwdown_qc_ar, id='LWdown_qc', axes = (time,lat,lon))
161   lwdown_qc.long_name = "Quality check for Lwdown" 
162   lwdown_qc.units = "-" 
163   lwdown_qc.missing_value = "1.e+20f"
164
165   Ts1_f_ar = N.array(listvar3[id_Ts1_f])
166   Ts1_f_ar.shape = (len(listvar3[id_Ts1_f]),1,1)
167   Ts1_f = cdms2.createVariable(Ts1_f_ar, id='TSOIL', axes = (time,lat,lon))
168   Ts1_f.long_name = "Soil Temperature" 
169   Ts1_f.units = "degC" 
170   Ts1_f.missing_value = -9999.
171
172   Rh_ar = N.array(listvar3[id_Rh])
173   Rh_ar.shape = (len(listvar3[id_Rh]),1,1)
174   Rh = cdms2.createVariable(Rh_ar, id='RH', axes = (time,lat,lon))
175   Rh.long_name = "Relative Humidity" 
176   Rh.units = "%" 
177   Rh.missing_value = -9999.
178
179   NEE_f_ar = N.array(listvar3[id_NEE_f])
180   NEE_f_ar.shape = (len(listvar3[id_NEE_f]),1,1)
181   NEE_f = cdms2.createVariable(NEE_f_ar, id='NEE', axes = (time,lat,lon))
182   NEE_f.long_name = "Net Ecosystem Exchange" 
183   NEE_f.units = "gC/m2/tstep" 
184   NEE_f.missing_value = -9999.
185
186   H_f_ar = N.array(listvar3[id_H_f])
187   H_f_ar.shape = (len(listvar3[id_H_f]),1,1)
188   H_f = cdms2.createVariable(H_f_ar, id='Fh', axes = (time,lat,lon))
189   H_f.long_name = "Sensible Heat Flux" 
190   H_f.units = "W/m2" 
191   H_f.missing_value = -9999.
192
193   LE_f_ar = N.array(listvar3[id_LE_f])
194   LE_f_ar.shape = (len(listvar3[id_LE_f]),1,1)
195   LE_f = cdms2.createVariable(LE_f_ar, id='Fle', axes = (time,lat,lon))
196   LE_f.long_name = "Latent Heat Flux" 
197   LE_f.units = "W/m2" 
198   LE_f.missing_value = -9999.
199
200   Epot_f_ar = N.array(listvar3[id_Epot_f])
201   Epot_f_ar.shape = (len(listvar3[id_Epot_f]),1,1)
202   Epot_f = cdms2.createVariable(Epot_f_ar, id='ET', axes = (time,lat,lon))
203   Epot_f.long_name = "Evapotranspiration" 
204   Epot_f.units = "mm/tstep" 
205   Epot_f.missing_value = -9999.
206
207   SWC1_f_ar = N.array(listvar3[id_SWC1_f])
208   SWC1_f_ar.shape = (len(listvar3[id_SWC1_f]),1,1)
209   SWC1_f = cdms2.createVariable(SWC1_f_ar, id='SWC', axes = (time,lat,lon))
210   SWC1_f.long_name = "Soil Water Content" 
211   SWC1_f.units = "VOL%" 
212   SWC1_f.missing_value = -9999.
213
214   gsurf_f_ar = N.array(listvar3[id_gsurf_f])
215   gsurf_f_ar.shape = (len(listvar3[id_gsurf_f]),1,1)
216   gsurf_f = cdms2.createVariable(gsurf_f_ar, id='GS', axes = (time,lat,lon))
217   gsurf_f.long_name = "Canopy conductance" 
218   gsurf_f.units = "mm/s" 
219   gsurf_f.missing_value = -9999.
220
221   GPP_f_ar = N.array(listvar3[id_GPP_f])
222   GPP_f_ar.shape = (len(listvar3[id_GPP_f]),1,1)
223   GPP_f = cdms2.createVariable(GPP_f_ar, id='GPP', axes = (time,lat,lon))
224   GPP_f.long_name = "Gross Primary Production" 
225   GPP_f.units = "gC/m2/tstep" 
226   GPP_f.missing_value = -9999.
227
228   Reco_ar = N.array(listvar3[id_Reco])
229   Reco_ar.shape = (len(listvar3[id_Reco]),1,1)
230   Reco = cdms2.createVariable(Reco_ar, id='Reco', axes = (time,lat,lon))
231   Reco.long_name = "Ecosystem Respiration" 
232   Reco.units = "gC/m2/tstep" 
233   Reco.missing_value = -9999.
234   
235   f = cdms2.open(name_path+name+'_'+str(year_start)+'-'+str(year_end)+'.nc', 'w')
236   f.write(tair,axes=None)
237   f.write(qair,axes=None)
238   f.write(psurf,axes=None)
239   f.write(wind_n,axes=None)
240   f.write(wind_e,axes=None)
241   f.write(rainf,axes=None)
242   f.write(snowf,axes=None)
243   f.write(swdown,axes=None)
244   f.write(lwdown,axes=None)
245
246   f.write(tair_qc,axes=None)
247   f.write(qair_qc,axes=None)
248   f.write(psurf_qc,axes=None)
249   f.write(wind_n_qc,axes=None)
250   f.write(wind_e_qc,axes=None)
251   f.write(rainf_qc,axes=None)
252   f.write(snowf_qc,axes=None)
253   f.write(swdown_qc,axes=None)
254   f.write(lwdown_qc,axes=None)
255
256   f.write(Ts1_f,axes=None)
257   f.write(Rh,axes=None)
258   f.write(NEE_f,axes=None)
259   f.write(H_f,axes=None)
260   f.write(LE_f,axes=None)
261   f.write(Epot_f,axes=None)
262   f.write(SWC1_f,axes=None)
263   f.write(gsurf_f,axes=None)
264   f.write(GPP_f,axes=None)
265   f.write(Reco,axes=None)
266   
267   f.close()
268# end function write_nc
269#-------------------------------------------------------------------------------------------
270
271
272
273def prepare_for_orchidee(weather_gapfill):
274   weather_out=copy.deepcopy(weather_gapfill)
275   # Conversion from VPD (Vapour Pressure Deficit, Pa) to Qair (Near Surface Specific Humidity, kg/kg)
276   for i in range(len(weather_out[id_vpd])):
277      # Magnus Tetens (Murray, 1967) http://cires.colorado.edu/~voemel/vp.html     
278      if(weather_out[id_tair][i]<273.15):
279         eg=c74*exp(c70*(weather_out[id_tair][i]-273.15)/(weather_out[id_tair][i]-c71))
280      else:
281         eg=c74*exp(c72*(weather_out[id_tair][i]-273.15)/(weather_out[id_tair][i]-c73))
282         
283      weather_out[id_vpd][i]=(eg-weather_out[id_vpd][i])*c75/(weather_out[id_psurf][i]-(eg-weather_out[id_vpd][i])*c76)
284
285   weather_out.append(N.zeros(len(weather_out[id_precip])))
286   
287   for i in range(len(weather_out[id_precip])):
288      if(weather_out[id_tair][i]<273.15):
289         weather_out[id_snow][i]=weather_out[id_precip][i]
290         weather_out[id_rain][i]=0.
291
292   return weather_out
293
294
295
296#-------------------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.