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

Last change on this file since 2975 was 2975, checked in by romain.vande, 18 months ago

Mars util :

Adapt expandstartfi to subslope dimension + small correction regarding depreciated dimension "number_of_advected_field"

RV

File size: 26.9 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
22integer :: corner(4),edges(4) ! to read data with a time axis
23character(len=90) :: varname ! name of a variable
24character(len=90) :: varatt ! name of attribute of a variable
25character(len=90) :: varattcontent ! content of the attribute
26! Input file dimension IDs:
27integer :: physical_points_dimid
28integer :: subsurface_layers_dimid
29integer :: nlayer_plus_1_dimid
30integer :: nslope_dimid
31integer :: number_of_advected_fields_dimid
32integer :: time_dimid
33integer :: nbindims ! number of dimensions in input file
34integer :: nbinvars ! number of variables in input file
35integer :: invarid ! to store ID of an input variable
36!Input file variables
37integer :: physical_points
38integer :: subsurface_layers
39integer :: nlayer_plus_1
40integer :: nslope
41integer :: number_of_advected_fields
42integer :: timelen
43real,allocatable :: surf_field(:) ! to store a 1D field of physical_points elements
44real,allocatable :: subsurf_field(:,:) ! to store subsurface (2D field)
45real,allocatable :: subslope_field(:,:)
46real,allocatable :: nlayer_plus_1_field(:,:)
47real,allocatable :: subsurf_subslope_field(:,:,:)
48
49! Output file dimensions:
50integer :: latlen
51integer :: lonlen
52! Output file variables:
53real,allocatable :: longitude(:)
54real,allocatable :: latitude(:)
55real,allocatable :: depth(:)
56real,allocatable :: out_surf_field(:,:)
57real,allocatable :: out_subsurf_field(:,:,:)
58real,allocatable :: out_subslope_field(:,:,:)
59real,allocatable :: out_nlayer_plus_1_field(:,:,:)
60real,allocatable :: out_subsurf_subslope_field(:,:,:,:)
61! Output file dimension IDs
62integer :: lon_dimid
63integer :: lat_dimid
64integer :: depth_dimid
65integer :: nslope_out_dimid
66integer :: nlayer_plus_1_out_dimid
67! IDs of Output file variables
68integer :: lon_varid
69integer :: lat_varid
70integer :: depth_varid
71
72integer :: i,j,k,ig0,ivar
73integer :: nbdim,nbatt,shape(4)
74integer :: nbarg ! # of program arguments
75character(len=256) :: arg ! to store a program argument
76real :: pi ! 3.14159...
77
78pi=2.*asin(1.)
79
80! 0. Preliminary stuff, check arguments (input and output files)
81! get number of arguments this program was called with
82nbarg=command_argument_count()
83
84if (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
98endif
99
100! 1. open input file
101status=nf90_open(infile,NF90_NOWRITE,inid)
102if (status.ne.nf90_noerr) then
103  write(*,*)"Failed to open datafile ",trim(infile)
104  write(*,*)trim(nf90_strerror(status))
105 stop
106endif
107write(*,*) "Reading input file: ",trim(infile)
108
109! 1.2 Identify/load dimensions in input file
110status=nf90_inq_dimid(inid,"physical_points",physical_points_dimid)
111if (status.ne.nf90_noerr) then
112  write(*,*)"Failed to find physical_points dimension"
113  write(*,*)trim(nf90_strerror(status))
114 stop
115endif
116status=nf90_inquire_dimension(inid,physical_points_dimid,len=physical_points)
117if (status.ne.nf90_noerr) then
118  write(*,*)"Failed to find physical_points value"
119  write(*,*)trim(nf90_strerror(status))
120 stop
121else
122  write(*,*) " physical_points = ",physical_points
123endif
124
125status=nf90_inq_dimid(inid,"subsurface_layers",subsurface_layers_dimid)
126if (status.ne.nf90_noerr) then
127  write(*,*)"Failed to find subsurface_layers dimension"
128  write(*,*)trim(nf90_strerror(status))
129 stop
130endif
131status=nf90_inquire_dimension(inid,subsurface_layers_dimid,len=subsurface_layers)
132if (status.ne.nf90_noerr) then
133  write(*,*)"Failed to find subsurface_layers value"
134  write(*,*)trim(nf90_strerror(status))
135 stop
136else
137  write(*,*) " subsurface_layers = ",subsurface_layers
138endif
139
140status=nf90_inq_dimid(inid,"nlayer_plus_1",nlayer_plus_1_dimid)
141if (status.ne.nf90_noerr) then
142  write(*,*)"Failed to find nlayer_plus_1 dimension"
143  write(*,*)trim(nf90_strerror(status))
144 stop
145endif
146status=nf90_inquire_dimension(inid,nlayer_plus_1_dimid,len=nlayer_plus_1)
147if (status.ne.nf90_noerr) then
148  write(*,*)"Failed to find nlayer_plus_1 value"
149  write(*,*)trim(nf90_strerror(status))
150 stop
151else
152  write(*,*) " nlayer_plus_1 = ",nlayer_plus_1
153endif
154
155status=nf90_inq_dimid(inid,"nslope",nslope_dimid)
156if (status.ne.nf90_noerr) then
157  write(*,*)"Failed to find slope dimension, old startfi file"
158  nslope=0
159  write(*,*)trim(nf90_strerror(status))
160else
161status=nf90_inquire_dimension(inid,nslope_dimid,len=nslope)
162if (status.ne.nf90_noerr) then
163  write(*,*)"Failed to find nslope value"
164  write(*,*)trim(nf90_strerror(status))
165 stop
166else
167  write(*,*) " nslope = ",nslope
168endif
169endif
170
171status=nf90_inq_dimid(inid,"Time",time_dimid)
172if (status.ne.nf90_noerr) then
173  write(*,*)"Failed to find Time dimension"
174  write(*,*)trim(nf90_strerror(status))
175  timelen = 0
176else
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
185endif
186
187! 1.3 Allocate memory for input fields
188allocate(surf_field(physical_points))
189allocate(subsurf_field(physical_points,subsurface_layers))
190allocate(subslope_field(physical_points,nslope))
191allocate(nlayer_plus_1_field(physical_points,nlayer_plus_1))
192allocate(subsurf_subslope_field(physical_points,subsurface_layers,nslope))
193
194! 2.1. Create output file
195status=nf90_create(outfile,NF90_CLOBBER,outid)
196if (status.ne.nf90_noerr) then
197  write(*,*) "Failed to create output file: ",trim(outfile)
198  write(*,*)trim(nf90_strerror(status))
199  stop
200endif
201write(*,*) " "
202write(*,*) "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
208status=nf90_inq_varid(inid,"longitude",varid)
209if (status.ne.nf90_noerr) then
210  write(*,*) "Failed to find longitude ID"
211  write(*,*)trim(nf90_strerror(status))
212  stop
213endif
214! read longitude
215status=nf90_get_var(inid,varid,surf_field)
216if (status.ne.nf90_noerr) then
217  write(*,*) "Failed to load longitude"
218  write(*,*)trim(nf90_strerror(status))
219  stop
220endif
221
222! count the number of points before longitude(i) wraps around
223i=3
224lonlen=1
225!write(*,*) "longitude(2)=",surf_field(2)
226do 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
230enddo
231! and add 1 because we want a redundant lon=180 point
232lonlen=lonlen+1
233write(*,*) " => lonlen=",lonlen
234
235allocate(longitude(lonlen))
236! fill longitude(1:lonlen)
237longitude(1:lonlen-1)=surf_field(2:lonlen)
238longitude(lonlen)=-longitude(1) ! redundant +Pi/2
239! convert to degrees
240longitude(:)=longitude(:)*(180./pi)
241
242! knowing lonlen, compute latlen
243latlen=2+(physical_points-2)/(lonlen-1)
244write(*,*) " => latlen=",latlen
245
246allocate(latitude(latlen))
247! get field "latitude(physical_points)" from infile
248status=nf90_inq_varid(inid,"latitude",varid)
249if (status.ne.nf90_noerr) then
250  write(*,*) "Failed to find latitude ID"
251  write(*,*)trim(nf90_strerror(status))
252  stop
253endif
254! read latitude
255status=nf90_get_var(inid,varid,surf_field)
256if (status.ne.nf90_noerr) then
257  write(*,*) "Failed to load latitude"
258  write(*,*)trim(nf90_strerror(status))
259  stop
260endif
261
262latitude(1)=surf_field(1)
263!write(*,*) "latitude(1)=",latitude(1)
264do i=2,latlen-1
265  latitude(i)=surf_field((i-1)*(lonlen-1))
266!  write(*,*) "i:",i,"latitude(i)=",latitude(i)
267enddo
268latitude(latlen)=surf_field(physical_points)
269!write(*,*) "latitude(latlen)=",latitude(latlen)
270! convert to degrees
271latitude(:)=latitude(:)*(180./pi)
272
273! depth:
274allocate(depth(subsurface_layers))
275! load "soildepth(subsurface_layers)" from input file
276status=nf90_inq_varid(inid,"soildepth",varid)
277if (status.ne.nf90_noerr) then
278  write(*,*) "Failed to find soildepth ID"
279  write(*,*)trim(nf90_strerror(status))
280  stop
281endif
282! read longitude
283status=nf90_get_var(inid,varid,depth)
284if (status.ne.nf90_noerr) then
285  write(*,*) "Failed to load soildepth"
286  write(*,*)trim(nf90_strerror(status))
287  stop
288endif
289
290! 2.2.2 define dimensions to output file
291! longitude:
292status=nf90_def_dim(outid,"longitude",lonlen,lon_dimid)
293if (status.ne.nf90_noerr) then
294  write(*,*) "Failed creating longitude dimension"
295  write(*,*)trim(nf90_strerror(status))
296  stop
297endif
298! latitude:
299status=nf90_def_dim(outid,"latitude",latlen,lat_dimid)
300if (status.ne.nf90_noerr) then
301  write(*,*) "Failed creating longitude dimension"
302  write(*,*)trim(nf90_strerror(status))
303  stop
304endif
305! depth:
306status=nf90_def_dim(outid,"depth",subsurface_layers,depth_dimid)
307if (status.ne.nf90_noerr) then
308  write(*,*) "Failed creating depth dimension"
309  write(*,*)trim(nf90_strerror(status))
310  stop
311endif
312! nslope:
313status=nf90_def_dim(outid,"nslope",nslope,nslope_out_dimid)
314if (status.ne.nf90_noerr) then
315  write(*,*) "Failed creating nslope dimension"
316  write(*,*)trim(nf90_strerror(status))
317  stop
318endif
319! nslope:
320status=nf90_def_dim(outid,"nlayer_plus_1",nlayer_plus_1,nlayer_plus_1_out_dimid)
321if (status.ne.nf90_noerr) then
322  write(*,*) "Failed creating nslope dimension"
323  write(*,*)trim(nf90_strerror(status))
324  stop
325endif
326
327!2.3. write variables to output file
328! 2.3.1. define 1D variables
329! longitude
330datashape(1)=lon_dimid
331status=nf90_def_var(outid,"longitude",NF90_REAL,lon_dimid,lon_varid)
332if (status.ne.nf90_noerr) then
333  write(*,*) "Failed creating longitude variable"
334  write(*,*)trim(nf90_strerror(status))
335  stop
336endif
337! longitude attributes
338status=nf90_put_att(outid,lon_varid,"long_name","East longitude")
339status=nf90_put_att(outid,lon_varid,"units","degrees_east")
340
341!latitude
342datashape(2)=lat_dimid
343status=nf90_def_var(outid,"latitude",NF90_REAL,lat_dimid,lat_varid)
344if (status.ne.nf90_noerr) then
345  write(*,*) "Failed creating latitude variable"
346  write(*,*)trim(nf90_strerror(status))
347  stop
348endif
349! latitude attributes
350status=nf90_put_att(outid,lat_varid,"long_name","North latitude")
351status=nf90_put_att(outid,lat_varid,"units","degrees_north")
352
353! depth
354status=nf90_def_var(outid,"depth",NF90_REAL,depth_dimid,depth_varid)
355if (status.ne.nf90_noerr) then
356  write(*,*) "Failed creating depth variable"
357  write(*,*)trim(nf90_strerror(status))
358  stop
359endif
360!depth atributes
361status=nf90_put_att(outid,depth_varid,"long_name","Soil mid-layer depth")
362status=nf90_put_att(outid,depth_varid,"units","m")
363status=nf90_put_att(outid,depth_varid,"positive","down")
364
365! 2.3.2 write 1D variable
366
367! swich out of NetCDF define mode
368status=nf90_enddef(outid)
369if (status.ne.nf90_noerr) then
370  write(*,*) "Failed to swich out of define mode"
371  write(*,*)trim(nf90_strerror(status))
372  stop
373endif
374
375! write longitude
376status=nf90_put_var(outid,lon_varid,longitude)
377if (status.ne.nf90_noerr) then
378  write(*,*) "Failed writing longitude"
379  write(*,*)trim(nf90_strerror(status))
380  stop
381endif
382
383! write latitude
384status=nf90_put_var(outid,lat_varid,latitude)
385if (status.ne.nf90_noerr) then
386  write(*,*) "Failed writing latitude"
387  write(*,*)trim(nf90_strerror(status))
388  stop
389endif
390
391! write depth
392status=nf90_put_var(outid,depth_varid,depth)
393if (status.ne.nf90_noerr) then
394  write(*,*) "Failed writing depth"
395  write(*,*)trim(nf90_strerror(status))
396  stop
397endif
398
399! 2.3. 2D (surface) variables
400! First find out how many variables there are in input file
401status=nf90_inquire(inid,nbindims,nbinvars)
402if (status.ne.nf90_noerr) then
403  write(*,*) "Failed obtaining nbindims and nbinvars from input file"
404  write(*,*)trim(nf90_strerror(status))
405  stop
406endif
407
408allocate(out_surf_field(lonlen,latlen))
409allocate(out_subsurf_field(lonlen,latlen,subsurface_layers))
410allocate(out_subslope_field(lonlen,latlen,nslope))
411allocate(out_nlayer_plus_1_field(lonlen,latlen,nlayer_plus_1))
412allocate(out_subsurf_subslope_field(lonlen,latlen,subsurface_layers,nslope))
413
414do 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
868enddo ! of do i=1,nbinvars
869
870! 3. Finish things and cleanup
871! Close output file
872status=nf90_close(outid)
873if (status.ne.nf90_noerr) then
874  write(*,*)"Failed to close file: ",trim(outfile)
875  write(*,*)trim(nf90_strerror(status))
876  stop
877endif
878
879  write(*,*)"Done writing file ",trim(outfile)
880
881end program expandstartfi
Note: See TracBrowser for help on using the repository browser.