source: trunk/LMDZ.MARS/util/expandstartfi.F90 @ 492

Last change on this file since 492 was 275, checked in by emillour, 13 years ago

Utilitaire pour GCM Mars: "expandstartfi" pour mettre le contenu d'un startfi.nc
sur une grille lonxlat pour faire des visus.
EM

File size: 16.7 KB
Line 
1program 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
10use netcdf
11implicit none
12
13! Input and output files:
14character(len=256) :: infile="startfi.nc"      ! default input file
15character(len=256) :: outfile="startfi_ex.nc"  ! default output file
16
17! NetCDF stuff
18integer :: status ! NetCDF return code
19integer :: inid,outid ! NetCDF file IDs
20integer :: varid ! to store the ID of a variable
21integer :: datashape(4) ! to store dimension IDs of a given dataset
22character(len=90) :: varname ! name of a variable
23character(len=90) :: varatt ! name of attribute of a variable
24character(len=90) :: varattcontent ! content of the attribute
25! Input file dimension IDs:
26integer :: physical_points_dimid
27integer :: subsurface_layers_dimid
28integer :: nlayer_plus_1_dimid
29integer :: number_of_advected_fields_dimid
30integer :: nbindims ! number of dimensions in input file
31integer :: nbinvars ! number of variables in input file
32integer :: invarid ! to store ID of an input variable
33!Input file variables
34integer :: physical_points
35integer :: subsurface_layers
36integer :: nlayer_plus_1
37integer :: number_of_advected_fields
38real,allocatable :: surf_field(:) ! to store a 1D field of physical_points elements
39real,allocatable :: subsurf_field(:,:) ! to store subsurface (2D field)
40
41! Output file dimensions:
42integer :: latlen
43integer :: lonlen
44! Output file variables:
45real,allocatable :: longitude(:)
46real,allocatable :: latitude(:)
47real,allocatable :: depth(:)
48real,allocatable :: out_surf_field(:,:)
49real,allocatable :: out_subsurf_field(:,:,:)
50! Output file dimension IDs
51integer :: lon_dimid
52integer :: lat_dimid
53integer :: depth_dimid
54! IDs of Output file variables
55integer :: lon_varid
56integer :: lat_varid
57integer :: depth_varid
58
59integer :: i,j,k,ig0,ivar
60integer :: nbdim,nbatt,shape(4)
61integer :: nbarg ! # of program arguments
62character(len=256) :: arg ! to store a program argument
63real :: pi ! 3.14159...
64
65pi=2.*asin(1.)
66
67! 0. Preliminary stuff, check arguments (input and output files)
68! get number of arguments this program was called with
69nbarg=command_argument_count()
70
71if (nbarg.ge.1) then
72  call get_command_argument(1,arg) ! get argument # 1
73  infile=trim(arg)
74  if (nbarg.eq.2) then
75    call get_command_argument(2,arg) ! get argument # 2
76    outfile=trim(arg)
77  else
78   ! build outfile from infile (replace ".nc" with "_ex.nc"
79   outfile=trim(infile(1:index(infile,".nc",back=.true.)-1))//"_ex.nc"
80  endif
81  if (nbarg.ge.3) then
82    write(*,*) ' Warning: Too many arguments...'
83    write(*,*) '         will only use the first 2 '
84  endif
85endif
86
87! 1. open input file
88status=nf90_open(infile,NF90_NOWRITE,inid)
89if (status.ne.nf90_noerr) then
90  write(*,*)"Failed to open datafile ",trim(infile)
91  write(*,*)trim(nf90_strerror(status))
92 stop
93endif
94write(*,*) "Reading input file: ",trim(infile)
95
96! 1.2 Identify/load dimensions in input file
97status=nf90_inq_dimid(inid,"physical_points",physical_points_dimid)
98if (status.ne.nf90_noerr) then
99  write(*,*)"Failed to find physical_points dimension"
100  write(*,*)trim(nf90_strerror(status))
101 stop
102endif
103status=nf90_inquire_dimension(inid,physical_points_dimid,len=physical_points)
104if (status.ne.nf90_noerr) then
105  write(*,*)"Failed to find physical_points value"
106  write(*,*)trim(nf90_strerror(status))
107 stop
108else
109  write(*,*) " physical_points = ",physical_points
110endif
111
112status=nf90_inq_dimid(inid,"subsurface_layers",subsurface_layers_dimid)
113if (status.ne.nf90_noerr) then
114  write(*,*)"Failed to find subsurface_layers dimension"
115  write(*,*)trim(nf90_strerror(status))
116 stop
117endif
118status=nf90_inquire_dimension(inid,subsurface_layers_dimid,len=subsurface_layers)
119if (status.ne.nf90_noerr) then
120  write(*,*)"Failed to find subsurface_layers value"
121  write(*,*)trim(nf90_strerror(status))
122 stop
123else
124  write(*,*) " subsurface_layers = ",subsurface_layers
125endif
126
127status=nf90_inq_dimid(inid,"nlayer_plus_1",nlayer_plus_1_dimid)
128if (status.ne.nf90_noerr) then
129  write(*,*)"Failed to find nlayer_plus_1 dimension"
130  write(*,*)trim(nf90_strerror(status))
131 stop
132endif
133status=nf90_inquire_dimension(inid,nlayer_plus_1_dimid,len=nlayer_plus_1)
134if (status.ne.nf90_noerr) then
135  write(*,*)"Failed to find nlayer_plus_1 value"
136  write(*,*)trim(nf90_strerror(status))
137 stop
138else
139  write(*,*) " nlayer_plus_1 = ",nlayer_plus_1
140endif
141
142status=nf90_inq_dimid(inid,"number_of_advected_fields",number_of_advected_fields_dimid)
143if (status.ne.nf90_noerr) then
144  write(*,*)"Failed to find number_of_advected_fields dimension"
145  write(*,*)trim(nf90_strerror(status))
146 stop
147endif
148status=nf90_inquire_dimension(inid,number_of_advected_fields_dimid,len=number_of_advected_fields)
149if (status.ne.nf90_noerr) then
150  write(*,*)"Failed to find number_of_advected_fields value"
151  write(*,*)trim(nf90_strerror(status))
152 stop
153else
154  write(*,*) " number_of_advected_fields = ",number_of_advected_fields
155endif
156
157! 1.3 Allocate memory for input fields
158allocate(surf_field(physical_points))
159allocate(subsurf_field(physical_points,subsurface_layers))
160
161! 2.1. Create output file
162status=nf90_create(outfile,NF90_CLOBBER,outid)
163if (status.ne.nf90_noerr) then
164  write(*,*) "Failed to create output file: ",trim(outfile)
165  write(*,*)trim(nf90_strerror(status))
166  stop
167endif
168write(*,*) " "
169write(*,*) "Created output file: ",trim(outfile)
170
171! 2.2. Build dimensions for output file
172
173! 2.2.1 Compute lonlen and latlen from  physical_points
174! load "longitude(physical_points)" from input file
175status=nf90_inq_varid(inid,"longitude",varid)
176if (status.ne.nf90_noerr) then
177  write(*,*) "Failed to find longitude ID"
178  write(*,*)trim(nf90_strerror(status))
179  stop
180endif
181! read longitude
182status=nf90_get_var(inid,varid,surf_field)
183if (status.ne.nf90_noerr) then
184  write(*,*) "Failed to load longitude"
185  write(*,*)trim(nf90_strerror(status))
186  stop
187endif
188
189! count the number of points before longitude(i) wraps around
190i=3
191lonlen=1
192!write(*,*) "longitude(2)=",surf_field(2)
193do while (surf_field(i).ne.surf_field(2))
194!write(*,*) "i:",i,"surf_field(i)=",surf_field(i)
195 i=i+1
196 lonlen=lonlen+1
197enddo
198! and add 1 because we want a redundant lon=180 point
199lonlen=lonlen+1
200write(*,*) " => lonlen=",lonlen
201
202allocate(longitude(lonlen))
203! fill longitude(1:lonlen)
204longitude(1:lonlen-1)=surf_field(2:lonlen)
205longitude(lonlen)=-longitude(1) ! redundant +Pi/2
206! convert to degrees
207longitude(:)=longitude(:)*(180./pi)
208
209! knowing lonlen, compute latlen
210latlen=2+(physical_points-2)/(lonlen-1)
211write(*,*) " => latlen=",latlen
212
213allocate(latitude(latlen))
214! get field "latitude(physical_points)" from infile
215status=nf90_inq_varid(inid,"latitude",varid)
216if (status.ne.nf90_noerr) then
217  write(*,*) "Failed to find latitude ID"
218  write(*,*)trim(nf90_strerror(status))
219  stop
220endif
221! read latitude
222status=nf90_get_var(inid,varid,surf_field)
223if (status.ne.nf90_noerr) then
224  write(*,*) "Failed to load latitude"
225  write(*,*)trim(nf90_strerror(status))
226  stop
227endif
228
229latitude(1)=surf_field(1)
230!write(*,*) "latitude(1)=",latitude(1)
231do i=2,latlen-1
232  latitude(i)=surf_field((i-1)*(lonlen-1))
233!  write(*,*) "i:",i,"latitude(i)=",latitude(i)
234enddo
235latitude(latlen)=surf_field(physical_points)
236!write(*,*) "latitude(latlen)=",latitude(latlen)
237! convert to degrees
238latitude(:)=latitude(:)*(180./pi)
239
240! depth:
241allocate(depth(subsurface_layers))
242! load "soildepth(subsurface_layers)" from input file
243status=nf90_inq_varid(inid,"soildepth",varid)
244if (status.ne.nf90_noerr) then
245  write(*,*) "Failed to find soildepth ID"
246  write(*,*)trim(nf90_strerror(status))
247  stop
248endif
249! read longitude
250status=nf90_get_var(inid,varid,depth)
251if (status.ne.nf90_noerr) then
252  write(*,*) "Failed to load soildepth"
253  write(*,*)trim(nf90_strerror(status))
254  stop
255endif
256
257! 2.2.2 define dimensions to output file
258! longitude:
259status=nf90_def_dim(outid,"longitude",lonlen,lon_dimid)
260if (status.ne.nf90_noerr) then
261  write(*,*) "Failed creating longitude dimension"
262  write(*,*)trim(nf90_strerror(status))
263  stop
264endif
265! latitude:
266status=nf90_def_dim(outid,"latitude",latlen,lat_dimid)
267if (status.ne.nf90_noerr) then
268  write(*,*) "Failed creating longitude dimension"
269  write(*,*)trim(nf90_strerror(status))
270  stop
271endif
272! depth:
273status=nf90_def_dim(outid,"depth",subsurface_layers,depth_dimid)
274if (status.ne.nf90_noerr) then
275  write(*,*) "Failed creating depth dimension"
276  write(*,*)trim(nf90_strerror(status))
277  stop
278endif
279
280!2.3. write variables to output file
281! 2.3.1. define 1D variables
282! longitude
283datashape(1)=lon_dimid
284status=nf90_def_var(outid,"longitude",NF90_REAL,lon_dimid,lon_varid)
285if (status.ne.nf90_noerr) then
286  write(*,*) "Failed creating longitude variable"
287  write(*,*)trim(nf90_strerror(status))
288  stop
289endif
290! longitude attributes
291status=nf90_put_att(outid,lon_varid,"long_name","East longitude")
292status=nf90_put_att(outid,lon_varid,"units","degrees_east")
293
294!latitude
295datashape(2)=lat_dimid
296status=nf90_def_var(outid,"latitude",NF90_REAL,lat_dimid,lat_varid)
297if (status.ne.nf90_noerr) then
298  write(*,*) "Failed creating latitude variable"
299  write(*,*)trim(nf90_strerror(status))
300  stop
301endif
302! latitude attributes
303status=nf90_put_att(outid,lat_varid,"long_name","North latitude")
304status=nf90_put_att(outid,lat_varid,"units","degrees_north")
305
306! depth
307status=nf90_def_var(outid,"depth",NF90_REAL,depth_dimid,depth_varid)
308if (status.ne.nf90_noerr) then
309  write(*,*) "Failed creating depth variable"
310  write(*,*)trim(nf90_strerror(status))
311  stop
312endif
313!depth atributes
314status=nf90_put_att(outid,depth_varid,"long_name","Soil mid-layer depth")
315status=nf90_put_att(outid,depth_varid,"units","m")
316status=nf90_put_att(outid,depth_varid,"positive","down")
317
318! 2.3.2 write 1D variable
319
320! swich out of NetCDF define mode
321status=nf90_enddef(outid)
322if (status.ne.nf90_noerr) then
323  write(*,*) "Failed to swich out of define mode"
324  write(*,*)trim(nf90_strerror(status))
325  stop
326endif
327
328! write longitude
329status=nf90_put_var(outid,lon_varid,longitude)
330if (status.ne.nf90_noerr) then
331  write(*,*) "Failed writing longitude"
332  write(*,*)trim(nf90_strerror(status))
333  stop
334endif
335
336! write latitude
337status=nf90_put_var(outid,lat_varid,latitude)
338if (status.ne.nf90_noerr) then
339  write(*,*) "Failed writing latitude"
340  write(*,*)trim(nf90_strerror(status))
341  stop
342endif
343
344! write depth
345status=nf90_put_var(outid,depth_varid,depth)
346if (status.ne.nf90_noerr) then
347  write(*,*) "Failed writing depth"
348  write(*,*)trim(nf90_strerror(status))
349  stop
350endif
351
352! 2.3. 2D (surface) variables
353! First find out how many variables there are in input file
354status=nf90_inquire(inid,nbindims,nbinvars)
355if (status.ne.nf90_noerr) then
356  write(*,*) "Failed obtaining nbindims and nbinvars from input file"
357  write(*,*)trim(nf90_strerror(status))
358  stop
359endif
360
361allocate(out_surf_field(lonlen,latlen))
362
363do ivar=1,nbinvars ! loop on all input variables
364  ! find out what dimensions are linked to this variable
365  status=nf90_inquire_variable(inid,ivar,name=varname,ndims=nbdim,&
366                               dimids=shape,natts=nbatt)
367  if ((nbdim==1).and.(shape(1)==physical_points_dimid)) then
368   
369    ! skip "longitude" and "latitude"
370    if (trim(varname)=="longitude") cycle
371    if (trim(varname)=="latitude") cycle
372   
373    write(*,*) " processing: ",trim(varname)
374
375    ! load input data:
376    status=nf90_inq_varid(inid,varname,invarid)
377    status=nf90_get_var(inid,invarid,surf_field)
378   
379    ! switch output file to to define mode
380    status=nf90_redef(outid)
381    if (status.ne.nf90_noerr) then
382      write(*,*) "Failed to swich to define mode"
383      write(*,*)trim(nf90_strerror(status))
384      stop
385    endif
386    !define the variable
387    status=nf90_def_var(outid,trim(varname),NF90_REAL,&
388                         (/lon_dimid,lat_dimid/),varid)
389    if (status.ne.nf90_noerr) then
390      write(*,*) "Failed creating variable ",trim(varname)
391      write(*,*)trim(nf90_strerror(status))
392      stop
393    endif
394
395    ! variable attributes
396    do k=1,nbatt
397      status=nf90_inq_attname(inid,invarid,k,varatt)
398      if (status.ne.nf90_noerr) then
399        write(*,*) "Failed getting attribute number",k," for ",trim(varname)
400        write(*,*)trim(nf90_strerror(status))
401      stop
402      endif
403      status=nf90_get_att(inid,invarid,trim(varatt),varattcontent)
404      if (status.ne.nf90_noerr) then
405        write(*,*) "Failed loading attribute ",trim(varatt)
406        write(*,*)trim(nf90_strerror(status))
407      stop
408      endif
409
410      ! write the attribute and its value to output
411      status=nf90_put_att(outid,varid,trim(varatt),trim(varattcontent))
412      if (status.ne.nf90_noerr) then
413        write(*,*) "Failed writing attribute ",trim(varatt)
414        write(*,*)trim(nf90_strerror(status))
415      stop
416      endif
417    enddo ! do k=1,nbatt
418
419    ! swich out of NetCDF define mode
420    status=nf90_enddef(outid)
421    if (status.ne.nf90_noerr) then
422      write(*,*) "Failed to swich out of define mode"
423      write(*,*)trim(nf90_strerror(status))
424      stop
425    endif
426   
427    ! "convert" from surf_field(:) to out_surf_field(:,:)
428    do i=1,lonlen
429      out_surf_field(i,1)=surf_field(1) ! North Pole
430      out_surf_field(i,latlen)=surf_field(1) ! South Pole
431    enddo
432    do j=2,latlen-1
433      ig0=1+(j-2)*(lonlen-1)
434      do i=1,lonlen-1
435        out_surf_field(i,j)=surf_field(ig0+i)
436      enddo
437      out_surf_field(lonlen,j)=out_surf_field(1,j) ! redundant lon=180
438    enddo
439   
440    ! write the variable
441    status=nf90_put_var(outid,varid,out_surf_field)
442    if (status.ne.nf90_noerr) then
443      write(*,*) "Failed to write ",trim(varname)
444      write(*,*)trim(nf90_strerror(status))
445      stop
446    endif
447  endif ! of if ((nbdim==1).and.(shape(1)==physical_points_dimid))
448enddo ! of do i=1,nbinvars
449
450! 2.4. 3D (subsurface) variables
451allocate(out_subsurf_field(lonlen,latlen,subsurface_layers))
452
453do ivar=1,nbinvars ! loop on all input variables
454  ! find out what dimensions are linked to this variable
455  status=nf90_inquire_variable(inid,ivar,name=varname,ndims=nbdim,&
456                               dimids=shape,natts=nbatt)
457  if ((nbdim==2).and.(shape(1)==physical_points_dimid) &
458                .and.(shape(2)==subsurface_layers_dimid)) then
459   
460    write(*,*) " processing: ",trim(varname)
461
462    ! load input data:
463    status=nf90_inq_varid(inid,varname,invarid)
464    status=nf90_get_var(inid,invarid,subsurf_field)
465   
466    ! switch output file to to define mode
467    status=nf90_redef(outid)
468    if (status.ne.nf90_noerr) then
469      write(*,*) "Failed to swich to define mode"
470      write(*,*)trim(nf90_strerror(status))
471      stop
472    endif
473    !define the variable
474    status=nf90_def_var(outid,trim(varname),NF90_REAL,&
475                         (/lon_dimid,lat_dimid,depth_dimid/),varid)
476    if (status.ne.nf90_noerr) then
477      write(*,*) "Failed creating variable ",trim(varname)
478      write(*,*)trim(nf90_strerror(status))
479      stop
480    endif
481
482    ! variable attributes
483    do k=1,nbatt
484      status=nf90_inq_attname(inid,invarid,k,varatt)
485      if (status.ne.nf90_noerr) then
486        write(*,*) "Failed getting attribute number",k," for ",trim(varname)
487        write(*,*)trim(nf90_strerror(status))
488      stop
489      endif
490      status=nf90_get_att(inid,invarid,trim(varatt),varattcontent)
491      if (status.ne.nf90_noerr) then
492        write(*,*) "Failed loading attribute ",trim(varatt)
493        write(*,*)trim(nf90_strerror(status))
494      stop
495      endif
496
497      ! write the attribute and its value to output
498      status=nf90_put_att(outid,varid,trim(varatt),trim(varattcontent))
499      if (status.ne.nf90_noerr) then
500        write(*,*) "Failed writing attribute ",trim(varatt)
501        write(*,*)trim(nf90_strerror(status))
502      stop
503      endif
504    enddo ! do k=1,nbatt
505
506    ! swich out of NetCDF define mode
507    status=nf90_enddef(outid)
508    if (status.ne.nf90_noerr) then
509      write(*,*) "Failed to swich out of define mode"
510      write(*,*)trim(nf90_strerror(status))
511      stop
512    endif
513   
514    ! "convert" from subsurf_field(:,:) to out_subsurf_field(:,:,:)
515    do i=1,lonlen
516      out_subsurf_field(i,1,:)=subsurf_field(1,:) ! North Pole
517      out_subsurf_field(i,latlen,:)=subsurf_field(1,:) ! South Pole
518    enddo
519    do j=2,latlen-1
520      ig0=1+(j-2)*(lonlen-1)
521      do i=1,lonlen-1
522        out_subsurf_field(i,j,:)=subsurf_field(ig0+i,:)
523      enddo
524      out_subsurf_field(lonlen,j,:)=out_subsurf_field(1,j,:) ! redundant lon=180
525    enddo
526   
527    ! write the variable
528    status=nf90_put_var(outid,varid,out_subsurf_field)
529    if (status.ne.nf90_noerr) then
530      write(*,*) "Failed to write ",trim(varname)
531      write(*,*)trim(nf90_strerror(status))
532      stop
533    endif
534  endif ! of if ((nbdim==1).and.(shape(1)==physical_points_dimid)...)
535enddo ! of do i=1,nbinvars
536
537
538! 3. Finish things and cleanup
539! Close output file
540status=nf90_close(outid)
541if (status.ne.nf90_noerr) then
542  write(*,*)"Failed to close file: ",trim(outfile)
543  write(*,*)trim(nf90_strerror(status))
544  stop
545endif
546
547end program expandstartfi
Note: See TracBrowser for help on using the repository browser.