1 | program expandstartfi |
---|
2 | |
---|
3 | ! program which takes a GCM startfi.nc file and puts it on the |
---|
4 | ! corresponding lonxlat grid (so it can be displayed using Grads, Ferret, etc.) |
---|
5 | ! usage: |
---|
6 | ! expandstartfi [infile.nc] [outfile.nc] |
---|
7 | ! if outfile is not specified it is built as "infile_ex.nc" |
---|
8 | ! if infile is not specified, "startfi.nc" is used as default |
---|
9 | |
---|
10 | use netcdf |
---|
11 | implicit none |
---|
12 | |
---|
13 | ! Input and output files: |
---|
14 | character(len=256) :: infile="startfi.nc" ! default input file |
---|
15 | character(len=256) :: outfile="startfi_ex.nc" ! default output file |
---|
16 | |
---|
17 | ! NetCDF stuff |
---|
18 | integer :: status ! NetCDF return code |
---|
19 | integer :: inid,outid ! NetCDF file IDs |
---|
20 | integer :: varid ! to store the ID of a variable |
---|
21 | integer :: datashape(4) ! to store dimension IDs of a given dataset |
---|
22 | integer :: corner(4),edges(4) ! to read data with a time axis |
---|
23 | character(len=90) :: varname ! name of a variable |
---|
24 | character(len=90) :: varatt ! name of attribute of a variable |
---|
25 | character(len=90) :: varattcontent ! content of the attribute |
---|
26 | ! Input file dimension IDs: |
---|
27 | integer :: physical_points_dimid |
---|
28 | integer :: subsurface_layers_dimid |
---|
29 | integer :: nlayer_plus_1_dimid |
---|
30 | integer :: nslope_dimid |
---|
31 | integer :: number_of_advected_fields_dimid |
---|
32 | integer :: time_dimid |
---|
33 | integer :: nbindims ! number of dimensions in input file |
---|
34 | integer :: nbinvars ! number of variables in input file |
---|
35 | integer :: invarid ! to store ID of an input variable |
---|
36 | !Input file variables |
---|
37 | integer :: physical_points |
---|
38 | integer :: subsurface_layers |
---|
39 | integer :: nlayer_plus_1 |
---|
40 | integer :: nslope |
---|
41 | integer :: number_of_advected_fields |
---|
42 | integer :: timelen |
---|
43 | real,allocatable :: surf_field(:) ! to store a 1D field of physical_points elements |
---|
44 | real,allocatable :: subsurf_field(:,:) ! to store subsurface (2D field) |
---|
45 | real,allocatable :: subslope_field(:,:) |
---|
46 | real,allocatable :: nlayer_plus_1_field(:,:) |
---|
47 | real,allocatable :: subsurf_subslope_field(:,:,:) |
---|
48 | |
---|
49 | ! Output file dimensions: |
---|
50 | integer :: latlen |
---|
51 | integer :: lonlen |
---|
52 | ! Output file variables: |
---|
53 | real,allocatable :: longitude(:) |
---|
54 | real,allocatable :: latitude(:) |
---|
55 | real,allocatable :: depth(:) |
---|
56 | real,allocatable :: out_surf_field(:,:) |
---|
57 | real,allocatable :: out_subsurf_field(:,:,:) |
---|
58 | real,allocatable :: out_subslope_field(:,:,:) |
---|
59 | real,allocatable :: out_nlayer_plus_1_field(:,:,:) |
---|
60 | real,allocatable :: out_subsurf_subslope_field(:,:,:,:) |
---|
61 | ! Output file dimension IDs |
---|
62 | integer :: lon_dimid |
---|
63 | integer :: lat_dimid |
---|
64 | integer :: depth_dimid |
---|
65 | integer :: nslope_out_dimid |
---|
66 | integer :: nlayer_plus_1_out_dimid |
---|
67 | ! IDs of Output file variables |
---|
68 | integer :: lon_varid |
---|
69 | integer :: lat_varid |
---|
70 | integer :: depth_varid |
---|
71 | |
---|
72 | integer :: i,j,k,ig0,ivar |
---|
73 | integer :: nbdim,nbatt,shape(4) |
---|
74 | integer :: nbarg ! # of program arguments |
---|
75 | character(len=256) :: arg ! to store a program argument |
---|
76 | real :: pi ! 3.14159... |
---|
77 | |
---|
78 | pi=2.*asin(1.) |
---|
79 | |
---|
80 | ! 0. Preliminary stuff, check arguments (input and output files) |
---|
81 | ! get number of arguments this program was called with |
---|
82 | nbarg=command_argument_count() |
---|
83 | |
---|
84 | if (nbarg.ge.1) then |
---|
85 | call get_command_argument(1,arg) ! get argument # 1 |
---|
86 | infile=trim(arg) |
---|
87 | if (nbarg.eq.2) then |
---|
88 | call get_command_argument(2,arg) ! get argument # 2 |
---|
89 | outfile=trim(arg) |
---|
90 | else |
---|
91 | ! build outfile from infile (replace ".nc" with "_ex.nc" |
---|
92 | outfile=trim(infile(1:index(infile,".nc",back=.true.)-1))//"_ex.nc" |
---|
93 | endif |
---|
94 | if (nbarg.ge.3) then |
---|
95 | write(*,*) ' Warning: Too many arguments...' |
---|
96 | write(*,*) ' will only use the first 2 ' |
---|
97 | endif |
---|
98 | endif |
---|
99 | |
---|
100 | ! 1. open input file |
---|
101 | status=nf90_open(infile,NF90_NOWRITE,inid) |
---|
102 | if (status.ne.nf90_noerr) then |
---|
103 | write(*,*)"Failed to open datafile ",trim(infile) |
---|
104 | write(*,*)trim(nf90_strerror(status)) |
---|
105 | stop |
---|
106 | endif |
---|
107 | write(*,*) "Reading input file: ",trim(infile) |
---|
108 | |
---|
109 | ! 1.2 Identify/load dimensions in input file |
---|
110 | status=nf90_inq_dimid(inid,"physical_points",physical_points_dimid) |
---|
111 | if (status.ne.nf90_noerr) then |
---|
112 | write(*,*)"Failed to find physical_points dimension" |
---|
113 | write(*,*)trim(nf90_strerror(status)) |
---|
114 | stop |
---|
115 | endif |
---|
116 | status=nf90_inquire_dimension(inid,physical_points_dimid,len=physical_points) |
---|
117 | if (status.ne.nf90_noerr) then |
---|
118 | write(*,*)"Failed to find physical_points value" |
---|
119 | write(*,*)trim(nf90_strerror(status)) |
---|
120 | stop |
---|
121 | else |
---|
122 | write(*,*) " physical_points = ",physical_points |
---|
123 | endif |
---|
124 | |
---|
125 | status=nf90_inq_dimid(inid,"subsurface_layers",subsurface_layers_dimid) |
---|
126 | if (status.ne.nf90_noerr) then |
---|
127 | write(*,*)"Failed to find subsurface_layers dimension" |
---|
128 | write(*,*)trim(nf90_strerror(status)) |
---|
129 | stop |
---|
130 | endif |
---|
131 | status=nf90_inquire_dimension(inid,subsurface_layers_dimid,len=subsurface_layers) |
---|
132 | if (status.ne.nf90_noerr) then |
---|
133 | write(*,*)"Failed to find subsurface_layers value" |
---|
134 | write(*,*)trim(nf90_strerror(status)) |
---|
135 | stop |
---|
136 | else |
---|
137 | write(*,*) " subsurface_layers = ",subsurface_layers |
---|
138 | endif |
---|
139 | |
---|
140 | status=nf90_inq_dimid(inid,"nlayer_plus_1",nlayer_plus_1_dimid) |
---|
141 | if (status.ne.nf90_noerr) then |
---|
142 | write(*,*)"Failed to find nlayer_plus_1 dimension" |
---|
143 | write(*,*)trim(nf90_strerror(status)) |
---|
144 | stop |
---|
145 | endif |
---|
146 | status=nf90_inquire_dimension(inid,nlayer_plus_1_dimid,len=nlayer_plus_1) |
---|
147 | if (status.ne.nf90_noerr) then |
---|
148 | write(*,*)"Failed to find nlayer_plus_1 value" |
---|
149 | write(*,*)trim(nf90_strerror(status)) |
---|
150 | stop |
---|
151 | else |
---|
152 | write(*,*) " nlayer_plus_1 = ",nlayer_plus_1 |
---|
153 | endif |
---|
154 | |
---|
155 | status=nf90_inq_dimid(inid,"nslope",nslope_dimid) |
---|
156 | if (status.ne.nf90_noerr) then |
---|
157 | write(*,*)"Failed to find slope dimension, old startfi file" |
---|
158 | nslope=0 |
---|
159 | write(*,*)trim(nf90_strerror(status)) |
---|
160 | else |
---|
161 | status=nf90_inquire_dimension(inid,nslope_dimid,len=nslope) |
---|
162 | if (status.ne.nf90_noerr) then |
---|
163 | write(*,*)"Failed to find nslope value" |
---|
164 | write(*,*)trim(nf90_strerror(status)) |
---|
165 | stop |
---|
166 | else |
---|
167 | write(*,*) " nslope = ",nslope |
---|
168 | endif |
---|
169 | endif |
---|
170 | |
---|
171 | status=nf90_inq_dimid(inid,"Time",time_dimid) |
---|
172 | if (status.ne.nf90_noerr) then |
---|
173 | write(*,*)"Failed to find Time dimension" |
---|
174 | write(*,*)trim(nf90_strerror(status)) |
---|
175 | timelen = 0 |
---|
176 | else |
---|
177 | status=nf90_inquire_dimension(inid,time_dimid,len=timelen) |
---|
178 | if (status.ne.nf90_noerr) then |
---|
179 | write(*,*)"Failed to read Time dimension" |
---|
180 | write(*,*)trim(nf90_strerror(status)) |
---|
181 | stop |
---|
182 | else |
---|
183 | write(*,*) " time length = ",timelen |
---|
184 | endif |
---|
185 | endif |
---|
186 | |
---|
187 | ! 1.3 Allocate memory for input fields |
---|
188 | allocate(surf_field(physical_points)) |
---|
189 | allocate(subsurf_field(physical_points,subsurface_layers)) |
---|
190 | allocate(subslope_field(physical_points,nslope)) |
---|
191 | allocate(nlayer_plus_1_field(physical_points,nlayer_plus_1)) |
---|
192 | allocate(subsurf_subslope_field(physical_points,subsurface_layers,nslope)) |
---|
193 | |
---|
194 | ! 2.1. Create output file |
---|
195 | status=nf90_create(outfile,NF90_CLOBBER,outid) |
---|
196 | if (status.ne.nf90_noerr) then |
---|
197 | write(*,*) "Failed to create output file: ",trim(outfile) |
---|
198 | write(*,*)trim(nf90_strerror(status)) |
---|
199 | stop |
---|
200 | endif |
---|
201 | write(*,*) " " |
---|
202 | write(*,*) "Created output file: ",trim(outfile) |
---|
203 | |
---|
204 | ! 2.2. Build dimensions for output file |
---|
205 | |
---|
206 | ! 2.2.1 Compute lonlen and latlen from physical_points |
---|
207 | ! load "longitude(physical_points)" from input file |
---|
208 | status=nf90_inq_varid(inid,"longitude",varid) |
---|
209 | if (status.ne.nf90_noerr) then |
---|
210 | write(*,*) "Failed to find longitude ID" |
---|
211 | write(*,*)trim(nf90_strerror(status)) |
---|
212 | stop |
---|
213 | endif |
---|
214 | ! read longitude |
---|
215 | status=nf90_get_var(inid,varid,surf_field) |
---|
216 | if (status.ne.nf90_noerr) then |
---|
217 | write(*,*) "Failed to load longitude" |
---|
218 | write(*,*)trim(nf90_strerror(status)) |
---|
219 | stop |
---|
220 | endif |
---|
221 | |
---|
222 | ! count the number of points before longitude(i) wraps around |
---|
223 | i=3 |
---|
224 | lonlen=1 |
---|
225 | !write(*,*) "longitude(2)=",surf_field(2) |
---|
226 | do while (surf_field(i).ne.surf_field(2)) |
---|
227 | !write(*,*) "i:",i,"surf_field(i)=",surf_field(i) |
---|
228 | i=i+1 |
---|
229 | lonlen=lonlen+1 |
---|
230 | enddo |
---|
231 | ! and add 1 because we want a redundant lon=180 point |
---|
232 | lonlen=lonlen+1 |
---|
233 | write(*,*) " => lonlen=",lonlen |
---|
234 | |
---|
235 | allocate(longitude(lonlen)) |
---|
236 | ! fill longitude(1:lonlen) |
---|
237 | longitude(1:lonlen-1)=surf_field(2:lonlen) |
---|
238 | longitude(lonlen)=-longitude(1) ! redundant +Pi/2 |
---|
239 | ! convert to degrees |
---|
240 | longitude(:)=longitude(:)*(180./pi) |
---|
241 | |
---|
242 | ! knowing lonlen, compute latlen |
---|
243 | latlen=2+(physical_points-2)/(lonlen-1) |
---|
244 | write(*,*) " => latlen=",latlen |
---|
245 | |
---|
246 | allocate(latitude(latlen)) |
---|
247 | ! get field "latitude(physical_points)" from infile |
---|
248 | status=nf90_inq_varid(inid,"latitude",varid) |
---|
249 | if (status.ne.nf90_noerr) then |
---|
250 | write(*,*) "Failed to find latitude ID" |
---|
251 | write(*,*)trim(nf90_strerror(status)) |
---|
252 | stop |
---|
253 | endif |
---|
254 | ! read latitude |
---|
255 | status=nf90_get_var(inid,varid,surf_field) |
---|
256 | if (status.ne.nf90_noerr) then |
---|
257 | write(*,*) "Failed to load latitude" |
---|
258 | write(*,*)trim(nf90_strerror(status)) |
---|
259 | stop |
---|
260 | endif |
---|
261 | |
---|
262 | latitude(1)=surf_field(1) |
---|
263 | !write(*,*) "latitude(1)=",latitude(1) |
---|
264 | do i=2,latlen-1 |
---|
265 | latitude(i)=surf_field((i-1)*(lonlen-1)) |
---|
266 | ! write(*,*) "i:",i,"latitude(i)=",latitude(i) |
---|
267 | enddo |
---|
268 | latitude(latlen)=surf_field(physical_points) |
---|
269 | !write(*,*) "latitude(latlen)=",latitude(latlen) |
---|
270 | ! convert to degrees |
---|
271 | latitude(:)=latitude(:)*(180./pi) |
---|
272 | |
---|
273 | ! depth: |
---|
274 | allocate(depth(subsurface_layers)) |
---|
275 | ! load "soildepth(subsurface_layers)" from input file |
---|
276 | status=nf90_inq_varid(inid,"soildepth",varid) |
---|
277 | if (status.ne.nf90_noerr) then |
---|
278 | write(*,*) "Failed to find soildepth ID" |
---|
279 | write(*,*)trim(nf90_strerror(status)) |
---|
280 | stop |
---|
281 | endif |
---|
282 | ! read longitude |
---|
283 | status=nf90_get_var(inid,varid,depth) |
---|
284 | if (status.ne.nf90_noerr) then |
---|
285 | write(*,*) "Failed to load soildepth" |
---|
286 | write(*,*)trim(nf90_strerror(status)) |
---|
287 | stop |
---|
288 | endif |
---|
289 | |
---|
290 | ! 2.2.2 define dimensions to output file |
---|
291 | ! longitude: |
---|
292 | status=nf90_def_dim(outid,"longitude",lonlen,lon_dimid) |
---|
293 | if (status.ne.nf90_noerr) then |
---|
294 | write(*,*) "Failed creating longitude dimension" |
---|
295 | write(*,*)trim(nf90_strerror(status)) |
---|
296 | stop |
---|
297 | endif |
---|
298 | ! latitude: |
---|
299 | status=nf90_def_dim(outid,"latitude",latlen,lat_dimid) |
---|
300 | if (status.ne.nf90_noerr) then |
---|
301 | write(*,*) "Failed creating longitude dimension" |
---|
302 | write(*,*)trim(nf90_strerror(status)) |
---|
303 | stop |
---|
304 | endif |
---|
305 | ! depth: |
---|
306 | status=nf90_def_dim(outid,"depth",subsurface_layers,depth_dimid) |
---|
307 | if (status.ne.nf90_noerr) then |
---|
308 | write(*,*) "Failed creating depth dimension" |
---|
309 | write(*,*)trim(nf90_strerror(status)) |
---|
310 | stop |
---|
311 | endif |
---|
312 | ! nslope: |
---|
313 | status=nf90_def_dim(outid,"nslope",nslope,nslope_out_dimid) |
---|
314 | if (status.ne.nf90_noerr) then |
---|
315 | write(*,*) "Failed creating nslope dimension" |
---|
316 | write(*,*)trim(nf90_strerror(status)) |
---|
317 | stop |
---|
318 | endif |
---|
319 | ! nslope: |
---|
320 | status=nf90_def_dim(outid,"nlayer_plus_1",nlayer_plus_1,nlayer_plus_1_out_dimid) |
---|
321 | if (status.ne.nf90_noerr) then |
---|
322 | write(*,*) "Failed creating nslope dimension" |
---|
323 | write(*,*)trim(nf90_strerror(status)) |
---|
324 | stop |
---|
325 | endif |
---|
326 | |
---|
327 | !2.3. write variables to output file |
---|
328 | ! 2.3.1. define 1D variables |
---|
329 | ! longitude |
---|
330 | datashape(1)=lon_dimid |
---|
331 | status=nf90_def_var(outid,"longitude",NF90_REAL,lon_dimid,lon_varid) |
---|
332 | if (status.ne.nf90_noerr) then |
---|
333 | write(*,*) "Failed creating longitude variable" |
---|
334 | write(*,*)trim(nf90_strerror(status)) |
---|
335 | stop |
---|
336 | endif |
---|
337 | ! longitude attributes |
---|
338 | status=nf90_put_att(outid,lon_varid,"long_name","East longitude") |
---|
339 | status=nf90_put_att(outid,lon_varid,"units","degrees_east") |
---|
340 | |
---|
341 | !latitude |
---|
342 | datashape(2)=lat_dimid |
---|
343 | status=nf90_def_var(outid,"latitude",NF90_REAL,lat_dimid,lat_varid) |
---|
344 | if (status.ne.nf90_noerr) then |
---|
345 | write(*,*) "Failed creating latitude variable" |
---|
346 | write(*,*)trim(nf90_strerror(status)) |
---|
347 | stop |
---|
348 | endif |
---|
349 | ! latitude attributes |
---|
350 | status=nf90_put_att(outid,lat_varid,"long_name","North latitude") |
---|
351 | status=nf90_put_att(outid,lat_varid,"units","degrees_north") |
---|
352 | |
---|
353 | ! depth |
---|
354 | status=nf90_def_var(outid,"depth",NF90_REAL,depth_dimid,depth_varid) |
---|
355 | if (status.ne.nf90_noerr) then |
---|
356 | write(*,*) "Failed creating depth variable" |
---|
357 | write(*,*)trim(nf90_strerror(status)) |
---|
358 | stop |
---|
359 | endif |
---|
360 | !depth atributes |
---|
361 | status=nf90_put_att(outid,depth_varid,"long_name","Soil mid-layer depth") |
---|
362 | status=nf90_put_att(outid,depth_varid,"units","m") |
---|
363 | status=nf90_put_att(outid,depth_varid,"positive","down") |
---|
364 | |
---|
365 | ! 2.3.2 write 1D variable |
---|
366 | |
---|
367 | ! swich out of NetCDF define mode |
---|
368 | status=nf90_enddef(outid) |
---|
369 | if (status.ne.nf90_noerr) then |
---|
370 | write(*,*) "Failed to swich out of define mode" |
---|
371 | write(*,*)trim(nf90_strerror(status)) |
---|
372 | stop |
---|
373 | endif |
---|
374 | |
---|
375 | ! write longitude |
---|
376 | status=nf90_put_var(outid,lon_varid,longitude) |
---|
377 | if (status.ne.nf90_noerr) then |
---|
378 | write(*,*) "Failed writing longitude" |
---|
379 | write(*,*)trim(nf90_strerror(status)) |
---|
380 | stop |
---|
381 | endif |
---|
382 | |
---|
383 | ! write latitude |
---|
384 | status=nf90_put_var(outid,lat_varid,latitude) |
---|
385 | if (status.ne.nf90_noerr) then |
---|
386 | write(*,*) "Failed writing latitude" |
---|
387 | write(*,*)trim(nf90_strerror(status)) |
---|
388 | stop |
---|
389 | endif |
---|
390 | |
---|
391 | ! write depth |
---|
392 | status=nf90_put_var(outid,depth_varid,depth) |
---|
393 | if (status.ne.nf90_noerr) then |
---|
394 | write(*,*) "Failed writing depth" |
---|
395 | write(*,*)trim(nf90_strerror(status)) |
---|
396 | stop |
---|
397 | endif |
---|
398 | |
---|
399 | ! 2.3. 2D (surface) variables |
---|
400 | ! First find out how many variables there are in input file |
---|
401 | status=nf90_inquire(inid,nbindims,nbinvars) |
---|
402 | if (status.ne.nf90_noerr) then |
---|
403 | write(*,*) "Failed obtaining nbindims and nbinvars from input file" |
---|
404 | write(*,*)trim(nf90_strerror(status)) |
---|
405 | stop |
---|
406 | endif |
---|
407 | |
---|
408 | allocate(out_surf_field(lonlen,latlen)) |
---|
409 | allocate(out_subsurf_field(lonlen,latlen,subsurface_layers)) |
---|
410 | allocate(out_subslope_field(lonlen,latlen,nslope)) |
---|
411 | allocate(out_nlayer_plus_1_field(lonlen,latlen,nlayer_plus_1)) |
---|
412 | allocate(out_subsurf_subslope_field(lonlen,latlen,subsurface_layers,nslope)) |
---|
413 | |
---|
414 | do ivar=1,nbinvars ! loop on all input variables |
---|
415 | shape(:) = 0 |
---|
416 | ! find out what dimensions are linked to this variable |
---|
417 | status=nf90_inquire_variable(inid,ivar,name=varname,ndims=nbdim,& |
---|
418 | dimids=shape,natts=nbatt) |
---|
419 | |
---|
420 | if (((nbdim==1).and.(shape(1)==physical_points_dimid))& |
---|
421 | .or.((nbdim==2).and.(shape(1)==physical_points_dimid)& |
---|
422 | .and.(shape(2)==time_dimid))) then |
---|
423 | |
---|
424 | corner(1) = 1 |
---|
425 | corner(2) = timelen |
---|
426 | edges(1) = physical_points |
---|
427 | edges(2) = 1 |
---|
428 | |
---|
429 | ! skip "longitude" and "latitude" |
---|
430 | if (trim(varname)=="longitude") cycle |
---|
431 | if (trim(varname)=="latitude") cycle |
---|
432 | |
---|
433 | write(*,*) " processing: ",trim(varname) |
---|
434 | |
---|
435 | ! load input data: |
---|
436 | status=nf90_inq_varid(inid,varname,invarid) |
---|
437 | status=nf90_get_var(inid,invarid,surf_field,corner,edges) |
---|
438 | |
---|
439 | ! switch output file to to define mode |
---|
440 | status=nf90_redef(outid) |
---|
441 | if (status.ne.nf90_noerr) then |
---|
442 | write(*,*) "Failed to swich to define mode" |
---|
443 | write(*,*)trim(nf90_strerror(status)) |
---|
444 | stop |
---|
445 | endif |
---|
446 | !define the variable |
---|
447 | status=nf90_def_var(outid,trim(varname),NF90_REAL,& |
---|
448 | (/lon_dimid,lat_dimid/),varid) |
---|
449 | if (status.ne.nf90_noerr) then |
---|
450 | write(*,*) "Failed creating variable ",trim(varname) |
---|
451 | write(*,*)trim(nf90_strerror(status)) |
---|
452 | stop |
---|
453 | endif |
---|
454 | |
---|
455 | ! variable attributes |
---|
456 | do k=1,nbatt |
---|
457 | status=nf90_inq_attname(inid,invarid,k,varatt) |
---|
458 | if (status.ne.nf90_noerr) then |
---|
459 | write(*,*) "Failed getting attribute number",k," for ",trim(varname) |
---|
460 | write(*,*)trim(nf90_strerror(status)) |
---|
461 | stop |
---|
462 | endif |
---|
463 | status=nf90_get_att(inid,invarid,trim(varatt),varattcontent) |
---|
464 | if (status.ne.nf90_noerr) then |
---|
465 | write(*,*) "Failed loading attribute ",trim(varatt) |
---|
466 | write(*,*)trim(nf90_strerror(status)) |
---|
467 | stop |
---|
468 | endif |
---|
469 | |
---|
470 | ! write the attribute and its value to output |
---|
471 | status=nf90_put_att(outid,varid,trim(varatt),trim(varattcontent)) |
---|
472 | if (status.ne.nf90_noerr) then |
---|
473 | write(*,*) "Failed writing attribute ",trim(varatt) |
---|
474 | write(*,*)trim(nf90_strerror(status)) |
---|
475 | stop |
---|
476 | endif |
---|
477 | enddo ! do k=1,nbatt |
---|
478 | |
---|
479 | ! swich out of NetCDF define mode |
---|
480 | status=nf90_enddef(outid) |
---|
481 | if (status.ne.nf90_noerr) then |
---|
482 | write(*,*) "Failed to swich out of define mode" |
---|
483 | write(*,*)trim(nf90_strerror(status)) |
---|
484 | stop |
---|
485 | endif |
---|
486 | |
---|
487 | ! "convert" from surf_field(:) to out_surf_field(:,:) |
---|
488 | do i=1,lonlen |
---|
489 | out_surf_field(i,1)=surf_field(1) ! North Pole |
---|
490 | out_surf_field(i,latlen)=surf_field(physical_points) ! South Pole |
---|
491 | enddo |
---|
492 | do j=2,latlen-1 |
---|
493 | ig0=1+(j-2)*(lonlen-1) |
---|
494 | do i=1,lonlen-1 |
---|
495 | out_surf_field(i,j)=surf_field(ig0+i) |
---|
496 | enddo |
---|
497 | out_surf_field(lonlen,j)=out_surf_field(1,j) ! redundant lon=180 |
---|
498 | enddo |
---|
499 | |
---|
500 | ! write the variable |
---|
501 | status=nf90_put_var(outid,varid,out_surf_field) |
---|
502 | if (status.ne.nf90_noerr) then |
---|
503 | write(*,*) "Failed to write ",trim(varname) |
---|
504 | write(*,*)trim(nf90_strerror(status)) |
---|
505 | stop |
---|
506 | endif |
---|
507 | endif ! of if ((nbdim==1).and.(shape(1)==physical_points_dimid)) |
---|
508 | |
---|
509 | if ((nbdim==2).and.(shape(1)==physical_points_dimid)& |
---|
510 | .and.(shape(2)==subsurface_layers_dimid)) then |
---|
511 | |
---|
512 | corner(1) = 1 |
---|
513 | corner(2) = 1 |
---|
514 | corner(3) = timelen |
---|
515 | edges(1) = physical_points |
---|
516 | edges(2) = subsurface_layers |
---|
517 | edges(3) = 1 |
---|
518 | |
---|
519 | write(*,*) " processing: ",trim(varname) |
---|
520 | |
---|
521 | ! load input data: |
---|
522 | status=nf90_inq_varid(inid,varname,invarid) |
---|
523 | status=nf90_get_var(inid,invarid,subsurf_field,corner,edges) |
---|
524 | |
---|
525 | ! switch output file to to define mode |
---|
526 | status=nf90_redef(outid) |
---|
527 | if (status.ne.nf90_noerr) then |
---|
528 | write(*,*) "Failed to swich to define mode" |
---|
529 | write(*,*)trim(nf90_strerror(status)) |
---|
530 | stop |
---|
531 | endif |
---|
532 | !define the variable |
---|
533 | status=nf90_def_var(outid,trim(varname),NF90_REAL,& |
---|
534 | (/lon_dimid,lat_dimid,depth_dimid/),varid) |
---|
535 | if (status.ne.nf90_noerr) then |
---|
536 | write(*,*) "Failed creating variable ",trim(varname) |
---|
537 | write(*,*)trim(nf90_strerror(status)) |
---|
538 | stop |
---|
539 | endif |
---|
540 | |
---|
541 | ! variable attributes |
---|
542 | do k=1,nbatt |
---|
543 | status=nf90_inq_attname(inid,invarid,k,varatt) |
---|
544 | if (status.ne.nf90_noerr) then |
---|
545 | write(*,*) "Failed getting attribute number",k," for ",trim(varname) |
---|
546 | write(*,*)trim(nf90_strerror(status)) |
---|
547 | stop |
---|
548 | endif |
---|
549 | status=nf90_get_att(inid,invarid,trim(varatt),varattcontent) |
---|
550 | if (status.ne.nf90_noerr) then |
---|
551 | write(*,*) "Failed loading attribute ",trim(varatt) |
---|
552 | write(*,*)trim(nf90_strerror(status)) |
---|
553 | stop |
---|
554 | endif |
---|
555 | |
---|
556 | ! write the attribute and its value to output |
---|
557 | status=nf90_put_att(outid,varid,trim(varatt),trim(varattcontent)) |
---|
558 | if (status.ne.nf90_noerr) then |
---|
559 | write(*,*) "Failed writing attribute ",trim(varatt) |
---|
560 | write(*,*)trim(nf90_strerror(status)) |
---|
561 | stop |
---|
562 | endif |
---|
563 | enddo ! do k=1,nbatt |
---|
564 | |
---|
565 | ! swich out of NetCDF define mode |
---|
566 | status=nf90_enddef(outid) |
---|
567 | if (status.ne.nf90_noerr) then |
---|
568 | write(*,*) "Failed to swich out of define mode" |
---|
569 | write(*,*)trim(nf90_strerror(status)) |
---|
570 | stop |
---|
571 | endif |
---|
572 | |
---|
573 | ! "convert" from subsurf_field(:,:) to out_subsurf_field(:,:,:) |
---|
574 | do i=1,lonlen |
---|
575 | out_subsurf_field(i,1,:)=subsurf_field(1,:) ! North Pole |
---|
576 | out_subsurf_field(i,latlen,:)=subsurf_field(physical_points,:) ! South Pole |
---|
577 | enddo |
---|
578 | do j=2,latlen-1 |
---|
579 | ig0=1+(j-2)*(lonlen-1) |
---|
580 | do i=1,lonlen-1 |
---|
581 | out_subsurf_field(i,j,:)=subsurf_field(ig0+i,:) |
---|
582 | enddo |
---|
583 | out_subsurf_field(lonlen,j,:)=out_subsurf_field(1,j,:) ! redundant lon=180 |
---|
584 | enddo |
---|
585 | |
---|
586 | ! write the variable |
---|
587 | status=nf90_put_var(outid,varid,out_subsurf_field) |
---|
588 | if (status.ne.nf90_noerr) then |
---|
589 | write(*,*) "Failed to write ",trim(varname) |
---|
590 | write(*,*)trim(nf90_strerror(status)) |
---|
591 | stop |
---|
592 | endif |
---|
593 | endif ! of if ((nbdim==1).and.(shape(1)==physical_points_dimid)...) |
---|
594 | |
---|
595 | |
---|
596 | |
---|
597 | |
---|
598 | if ((nbdim==3).and.(shape(1)==physical_points_dimid)& |
---|
599 | .and.(shape(2)==nslope_dimid)) then |
---|
600 | |
---|
601 | corner(1) = 1 |
---|
602 | corner(2) = 1 |
---|
603 | corner(3) = timelen |
---|
604 | edges(1) = physical_points |
---|
605 | edges(2) = nslope |
---|
606 | edges(3) = 1 |
---|
607 | |
---|
608 | write(*,*) " processing: ",trim(varname) |
---|
609 | |
---|
610 | ! load input data: |
---|
611 | status=nf90_inq_varid(inid,varname,invarid) |
---|
612 | status=nf90_get_var(inid,invarid,subslope_field,corner,edges) |
---|
613 | |
---|
614 | ! switch output file to to define mode |
---|
615 | status=nf90_redef(outid) |
---|
616 | if (status.ne.nf90_noerr) then |
---|
617 | write(*,*) "Failed to swich to define mode" |
---|
618 | write(*,*)trim(nf90_strerror(status)) |
---|
619 | stop |
---|
620 | endif |
---|
621 | !define the variable |
---|
622 | status=nf90_def_var(outid,trim(varname),NF90_REAL,& |
---|
623 | (/lon_dimid,lat_dimid,nslope_out_dimid/),varid) |
---|
624 | if (status.ne.nf90_noerr) then |
---|
625 | write(*,*) "Failed creating variable ",trim(varname) |
---|
626 | write(*,*)trim(nf90_strerror(status)) |
---|
627 | stop |
---|
628 | endif |
---|
629 | |
---|
630 | ! variable attributes |
---|
631 | do k=1,nbatt |
---|
632 | status=nf90_inq_attname(inid,invarid,k,varatt) |
---|
633 | if (status.ne.nf90_noerr) then |
---|
634 | write(*,*) "Failed getting attribute number",k," for ",trim(varname) |
---|
635 | write(*,*)trim(nf90_strerror(status)) |
---|
636 | stop |
---|
637 | endif |
---|
638 | status=nf90_get_att(inid,invarid,trim(varatt),varattcontent) |
---|
639 | if (status.ne.nf90_noerr) then |
---|
640 | write(*,*) "Failed loading attribute ",trim(varatt) |
---|
641 | write(*,*)trim(nf90_strerror(status)) |
---|
642 | stop |
---|
643 | endif |
---|
644 | |
---|
645 | ! write the attribute and its value to output |
---|
646 | status=nf90_put_att(outid,varid,trim(varatt),trim(varattcontent)) |
---|
647 | if (status.ne.nf90_noerr) then |
---|
648 | write(*,*) "Failed writing attribute ",trim(varatt) |
---|
649 | write(*,*)trim(nf90_strerror(status)) |
---|
650 | stop |
---|
651 | endif |
---|
652 | enddo ! do k=1,nbatt |
---|
653 | |
---|
654 | ! swich out of NetCDF define mode |
---|
655 | status=nf90_enddef(outid) |
---|
656 | if (status.ne.nf90_noerr) then |
---|
657 | write(*,*) "Failed to swich out of define mode" |
---|
658 | write(*,*)trim(nf90_strerror(status)) |
---|
659 | stop |
---|
660 | endif |
---|
661 | |
---|
662 | ! "convert" from subsurf_field(:,:) to out_subslope_field(:,:,:) |
---|
663 | do i=1,lonlen |
---|
664 | out_subslope_field(i,1,:)=subslope_field(1,:) ! North Pole |
---|
665 | out_subslope_field(i,latlen,:)=subslope_field(physical_points,:) ! South Pole |
---|
666 | enddo |
---|
667 | do j=2,latlen-1 |
---|
668 | ig0=1+(j-2)*(lonlen-1) |
---|
669 | do i=1,lonlen-1 |
---|
670 | out_subslope_field(i,j,:)=subslope_field(ig0+i,:) |
---|
671 | enddo |
---|
672 | out_subslope_field(lonlen,j,:)=out_subslope_field(1,j,:) ! redundant lon=180 |
---|
673 | enddo |
---|
674 | |
---|
675 | ! write the variable |
---|
676 | status=nf90_put_var(outid,varid,out_subslope_field) |
---|
677 | if (status.ne.nf90_noerr) then |
---|
678 | write(*,*) "Failed to write ",trim(varname) |
---|
679 | write(*,*)trim(nf90_strerror(status)) |
---|
680 | stop |
---|
681 | endif |
---|
682 | endif ! of if ((nbdim==1).and.(shape(1)==physical_points_dimid)...) |
---|
683 | |
---|
684 | |
---|
685 | |
---|
686 | |
---|
687 | |
---|
688 | if ((nbdim==3).and.(shape(1)==physical_points_dimid)& |
---|
689 | .and.(shape(2)==nlayer_plus_1_dimid)) then |
---|
690 | |
---|
691 | corner(1) = 1 |
---|
692 | corner(2) = 1 |
---|
693 | corner(3) = timelen |
---|
694 | edges(1) = physical_points |
---|
695 | edges(2) = nlayer_plus_1 |
---|
696 | edges(3) = 1 |
---|
697 | |
---|
698 | write(*,*) " processing: ",trim(varname) |
---|
699 | |
---|
700 | ! load input data: |
---|
701 | status=nf90_inq_varid(inid,varname,invarid) |
---|
702 | status=nf90_get_var(inid,invarid,nlayer_plus_1_field,corner,edges) |
---|
703 | |
---|
704 | ! switch output file to to define mode |
---|
705 | status=nf90_redef(outid) |
---|
706 | if (status.ne.nf90_noerr) then |
---|
707 | write(*,*) "Failed to swich to define mode" |
---|
708 | write(*,*)trim(nf90_strerror(status)) |
---|
709 | stop |
---|
710 | endif |
---|
711 | !define the variable |
---|
712 | status=nf90_def_var(outid,trim(varname),NF90_REAL,& |
---|
713 | (/lon_dimid,lat_dimid,nlayer_plus_1_out_dimid/),varid) |
---|
714 | if (status.ne.nf90_noerr) then |
---|
715 | write(*,*) "Failed creating variable ",trim(varname) |
---|
716 | write(*,*)trim(nf90_strerror(status)) |
---|
717 | stop |
---|
718 | endif |
---|
719 | |
---|
720 | ! variable attributes |
---|
721 | do k=1,nbatt |
---|
722 | status=nf90_inq_attname(inid,invarid,k,varatt) |
---|
723 | if (status.ne.nf90_noerr) then |
---|
724 | write(*,*) "Failed getting attribute number",k," for ",trim(varname) |
---|
725 | write(*,*)trim(nf90_strerror(status)) |
---|
726 | stop |
---|
727 | endif |
---|
728 | status=nf90_get_att(inid,invarid,trim(varatt),varattcontent) |
---|
729 | if (status.ne.nf90_noerr) then |
---|
730 | write(*,*) "Failed loading attribute ",trim(varatt) |
---|
731 | write(*,*)trim(nf90_strerror(status)) |
---|
732 | stop |
---|
733 | endif |
---|
734 | |
---|
735 | ! write the attribute and its value to output |
---|
736 | status=nf90_put_att(outid,varid,trim(varatt),trim(varattcontent)) |
---|
737 | if (status.ne.nf90_noerr) then |
---|
738 | write(*,*) "Failed writing attribute ",trim(varatt) |
---|
739 | write(*,*)trim(nf90_strerror(status)) |
---|
740 | stop |
---|
741 | endif |
---|
742 | enddo ! do k=1,nbatt |
---|
743 | |
---|
744 | ! swich out of NetCDF define mode |
---|
745 | status=nf90_enddef(outid) |
---|
746 | if (status.ne.nf90_noerr) then |
---|
747 | write(*,*) "Failed to swich out of define mode" |
---|
748 | write(*,*)trim(nf90_strerror(status)) |
---|
749 | stop |
---|
750 | endif |
---|
751 | |
---|
752 | ! "convert" from subsurf_field(:,:) to out_nlayer_plus_1_field(:,:,:) |
---|
753 | do i=1,lonlen |
---|
754 | out_nlayer_plus_1_field(i,1,:)=nlayer_plus_1_field(1,:) ! North Pole |
---|
755 | out_nlayer_plus_1_field(i,latlen,:)=nlayer_plus_1_field(physical_points,:) ! South Pole |
---|
756 | enddo |
---|
757 | do j=2,latlen-1 |
---|
758 | ig0=1+(j-2)*(lonlen-1) |
---|
759 | do i=1,lonlen-1 |
---|
760 | out_nlayer_plus_1_field(i,j,:)=nlayer_plus_1_field(ig0+i,:) |
---|
761 | enddo |
---|
762 | out_nlayer_plus_1_field(lonlen,j,:)=out_nlayer_plus_1_field(1,j,:) ! redundant lon=180 |
---|
763 | enddo |
---|
764 | |
---|
765 | ! write the variable |
---|
766 | status=nf90_put_var(outid,varid,out_nlayer_plus_1_field) |
---|
767 | if (status.ne.nf90_noerr) then |
---|
768 | write(*,*) "Failed to write ",trim(varname) |
---|
769 | write(*,*)trim(nf90_strerror(status)) |
---|
770 | stop |
---|
771 | endif |
---|
772 | endif ! of if ((nbdim==1).and.(shape(1)==physical_points_dimid)...) |
---|
773 | |
---|
774 | |
---|
775 | |
---|
776 | |
---|
777 | if (nbdim==4) then |
---|
778 | |
---|
779 | corner(1) = 1 |
---|
780 | corner(2) = 1 |
---|
781 | corner(3) = 1 |
---|
782 | corner(4) = timelen |
---|
783 | edges(1) = physical_points |
---|
784 | edges(2) = subsurface_layers |
---|
785 | edges(3) = nslope |
---|
786 | edges(4) = 1 |
---|
787 | |
---|
788 | write(*,*) " processing: ",trim(varname) |
---|
789 | |
---|
790 | ! load input data: |
---|
791 | status=nf90_inq_varid(inid,varname,invarid) |
---|
792 | status=nf90_get_var(inid,invarid,subsurf_subslope_field,corner,edges) |
---|
793 | |
---|
794 | ! switch output file to to define mode |
---|
795 | status=nf90_redef(outid) |
---|
796 | if (status.ne.nf90_noerr) then |
---|
797 | write(*,*) "Failed to swich to define mode" |
---|
798 | write(*,*)trim(nf90_strerror(status)) |
---|
799 | stop |
---|
800 | endif |
---|
801 | !define the variable |
---|
802 | print *, "subsurface_layers_dimid", subsurface_layers_dimid |
---|
803 | print *, "nslope_out_dimid", nslope_out_dimid |
---|
804 | status=nf90_def_var(outid,trim(varname),NF90_REAL,& |
---|
805 | (/lon_dimid,lat_dimid,depth_dimid,nslope_out_dimid/),varid) |
---|
806 | if (status.ne.nf90_noerr) then |
---|
807 | write(*,*) "Failed creating variable ",trim(varname) |
---|
808 | write(*,*)trim(nf90_strerror(status)) |
---|
809 | stop |
---|
810 | endif |
---|
811 | |
---|
812 | ! variable attributes |
---|
813 | do k=1,nbatt |
---|
814 | status=nf90_inq_attname(inid,invarid,k,varatt) |
---|
815 | if (status.ne.nf90_noerr) then |
---|
816 | write(*,*) "Failed getting attribute number",k," for ",trim(varname) |
---|
817 | write(*,*)trim(nf90_strerror(status)) |
---|
818 | stop |
---|
819 | endif |
---|
820 | status=nf90_get_att(inid,invarid,trim(varatt),varattcontent) |
---|
821 | if (status.ne.nf90_noerr) then |
---|
822 | write(*,*) "Failed loading attribute ",trim(varatt) |
---|
823 | write(*,*)trim(nf90_strerror(status)) |
---|
824 | stop |
---|
825 | endif |
---|
826 | |
---|
827 | ! write the attribute and its value to output |
---|
828 | status=nf90_put_att(outid,varid,trim(varatt),trim(varattcontent)) |
---|
829 | if (status.ne.nf90_noerr) then |
---|
830 | write(*,*) "Failed writing attribute ",trim(varatt) |
---|
831 | write(*,*)trim(nf90_strerror(status)) |
---|
832 | stop |
---|
833 | endif |
---|
834 | enddo ! do k=1,nbatt |
---|
835 | |
---|
836 | ! swich out of NetCDF define mode |
---|
837 | status=nf90_enddef(outid) |
---|
838 | if (status.ne.nf90_noerr) then |
---|
839 | write(*,*) "Failed to swich out of define mode" |
---|
840 | write(*,*)trim(nf90_strerror(status)) |
---|
841 | stop |
---|
842 | endif |
---|
843 | |
---|
844 | ! "convert" from subsurf_field(:,:) to out_subsurf_subslope_field(:,:,:) |
---|
845 | do i=1,lonlen |
---|
846 | out_subsurf_subslope_field(i,1,:,:)=subsurf_subslope_field(1,:,:) ! North Pole |
---|
847 | out_subsurf_subslope_field(i,latlen,:,:)=subsurf_subslope_field(physical_points,:,:) ! South Pole |
---|
848 | enddo |
---|
849 | do j=2,latlen-1 |
---|
850 | ig0=1+(j-2)*(lonlen-1) |
---|
851 | do i=1,lonlen-1 |
---|
852 | out_subsurf_subslope_field(i,j,:,:)=subsurf_subslope_field(ig0+i,:,:) |
---|
853 | enddo |
---|
854 | out_subsurf_subslope_field(lonlen,j,:,:)=out_subsurf_subslope_field(1,j,:,:) ! redundant lon=180 |
---|
855 | enddo |
---|
856 | |
---|
857 | ! write the variable |
---|
858 | status=nf90_put_var(outid,varid,out_subsurf_subslope_field) |
---|
859 | if (status.ne.nf90_noerr) then |
---|
860 | write(*,*) "Failed to write ",trim(varname) |
---|
861 | write(*,*)trim(nf90_strerror(status)) |
---|
862 | stop |
---|
863 | endif |
---|
864 | endif ! of if ((nbdim==1).and.(shape(1)==physical_points_dimid)...) |
---|
865 | |
---|
866 | |
---|
867 | |
---|
868 | enddo ! of do i=1,nbinvars |
---|
869 | |
---|
870 | ! 3. Finish things and cleanup |
---|
871 | ! Close output file |
---|
872 | status=nf90_close(outid) |
---|
873 | if (status.ne.nf90_noerr) then |
---|
874 | write(*,*)"Failed to close file: ",trim(outfile) |
---|
875 | write(*,*)trim(nf90_strerror(status)) |
---|
876 | stop |
---|
877 | endif |
---|
878 | |
---|
879 | write(*,*)"Done writing file ",trim(outfile) |
---|
880 | |
---|
881 | end program expandstartfi |
---|