source: trunk/LMDZ.COMMON/libf/evolution/reshape_XIOS_output.F90 @ 3836

Last change on this file since 3836 was 3836, checked in by jbclement, 4 days ago

COMMON:
Rework related to the command-line parsing:

  • Replace ad-hoc argument parsing with a unified 'parse_args' subroutine, allowing easier extension to other programs and options across models;;
  • Use of '--version' (with ab optional output file) to print compilation/version details;
  • Addition of 'job_timelimit_mod' module to handle SLURM/PBS job time-limit via '--jobid' (currently only used in the PEM), allowing easier extension to other programs;
  • Replace manual SSO handling with 'parse_args' for the Mars start2archive;
  • Clean-up related legacy code in the programs supporting the version option.

JBC
File size: 16.3 KB
Line 
1PROGRAM reshape_XIOS_output
2
3!=======================================================================
4! Purpose: Read XIOS NetCDF files and convert them onto the PCM grid.
5!          XIOS longitudes run from -180 to +180 (exclusive). So we append
6!          the first longitude value again at the end in the output to
7!          complete the grid. Done for the two PCM years.
8!
9! Authors: RV & LL (original), JBC (optimized)
10!=======================================================================
11
12use netcdf
13use parse_args_mod, only: parse_args
14
15implicit none
16
17! Variables for NetCDF I/O and bookkeeping
18integer                            :: state
19integer                            :: ncid_in, ncid_out
20integer                            :: ndims, nvars, nGlobalAtts, unlimDimID
21integer, allocatable, dimension(:) :: dimids_in, varids_in
22integer, allocatable, dimension(:) :: dimids_out, varids_out
23
24! Store each input dimension name and length
25character(30), allocatable, dimension(:) :: dimNames
26integer,       allocatable, dimension(:) :: dimLens
27
28! Which input‐index corresponds to lon/lat/time/soil (–1 if not present)
29integer :: idx_lon_in  = -1
30integer :: idx_lat_in  = -1
31integer :: idx_time_in = -1
32integer :: idx_soil_in = -1
33
34! Lengths of key dims (input), plus output lon length
35integer :: len_lon_in, len_lat_in, len_time_in, len_soil_in
36integer :: len_lon_out
37
38! Loop and variable bookkeeping
39integer                            :: i, j, k
40integer                            :: numDimsVar, numAttsVar
41character(100)                     :: varName
42integer                            :: xtypeVar
43integer, allocatable, dimension(:) :: dimids_var_in
44
45! Buffers for reading/writing when first‐dim = lon (max‐sized)
46real, allocatable, dimension(:)       :: buf1D_in, buf1D_out
47real, allocatable, dimension(:,:)     :: buf2D_in, buf2D_out
48real, allocatable, dimension(:,:,:)   :: buf3D_in, buf3D_out
49real, allocatable, dimension(:,:,:,:) :: buf4D_in, buf4D_out
50
51! Temporaries for "non‐lon‐first" variables
52real, allocatable, dimension(:)       :: tmp1D
53real, allocatable, dimension(:,:)     :: tmp2D
54real, allocatable, dimension(:,:,:)   :: tmp3D
55real, allocatable, dimension(:,:,:,:) :: tmp4D
56
57! Temporaries for dimension inquiries
58integer       :: thisLen
59integer       :: len1, len2, len3, len4
60integer       :: lenDim2, lenDim3, lenDim4
61character(30) :: tmpDimName
62
63logical :: uses_lon_first
64
65! For looping over two "years"
66integer      :: numyear
67character(4) :: str
68
69! For deleting existing output
70integer :: cstat
71logical :: exists
72
73! CODE
74! Parse command-line options
75call parse_args()
76
77! Main loop: two PCM years
78do numyear = 1,2
79    write(str,'(I1.1)') numyear
80    write(*,*) "> Reshaping variables from ""data2reshape_Y"//trim(str)//".nc""..."
81
82    ! Open input file (read‐only)
83    state = nf90_open("data2reshape_Y"//trim(str)//".nc",mode = nf90_nowrite,ncid = ncid_in)
84    if (state /= nf90_noerr) call handle_err(state)
85
86    ! If output exists, delete it
87    inquire(file = "data_PCM_Y"//trim(str)//".nc",exist = exists)
88    if (exists) then
89        call execute_command_line("rm data_PCM_Y"//trim(str)//".nc",cmdstat = cstat)
90        if (cstat > 0) then
91            error stop 'Command exection failed!'
92        else if (cstat < 0) then
93            error stop 'Command execution not supported!'
94        endif
95    endif
96
97    ! Create output file in define mode
98    state = nf90_create("data_PCM_Y"//trim(str)//".nc",cmode = or(nf90_noclobber,nf90_64bit_offset),ncid = ncid_out)
99    if (state /= nf90_noerr) call handle_err(state)
100
101    ! Inquire input for dims, vars, global atts, unlimited dim ID
102    state = nf90_inquire(ncid_in,ndims,nvars,nGlobalAtts,unlimDimID)
103    if (state /= nf90_noerr) call handle_err(state)
104
105    ! Allocate arrays for dim IDs, var IDs, names, lengths
106    allocate(dimids_in(ndims),varids_in(nvars),dimids_out(ndims),varids_out(nvars),dimNames(ndims),dimLens(ndims))
107
108    ! Get the dimension IDs and then query each for its name and length
109    state = nf90_inq_dimids(ncid_in,ndims,dimids_in,unlimDimID)
110    if (state /= nf90_noerr) call handle_err(state)
111
112    do i = 1,ndims
113        state = nf90_inquire_dimension(ncid_in,dimids_in(i),dimNames(i),dimLens(i))
114        if (state /= nf90_noerr) call handle_err(state)
115
116        select case (trim(dimNames(i)))
117            case ("lon","longitude")
118                idx_lon_in = i
119                len_lon_in = dimLens(i)
120            case ("lat","latitude")
121                idx_lat_in = i
122                len_lat_in = dimLens(i)
123            case ("time_counter","Time")
124                idx_time_in = i
125                len_time_in = dimLens(i)
126            case ("soil_layers","subsurface_layers")
127                idx_soil_in = i
128                len_soil_in = dimLens(i)
129            case default
130                ! nothing special
131        end select
132
133        ! Define the same dimension in the output, except lon becoming (len_lon_in + 1)
134        if (i == idx_lon_in) then
135            len_lon_out = len_lon_in + 1
136            state = nf90_def_dim(ncid_out,trim(dimNames(i)),len_lon_out,dimids_out(i))
137        else
138            state = nf90_def_dim(ncid_out,trim(dimNames(i)),dimLens(i),dimids_out(i))
139        endif
140        if (state /= nf90_noerr) call handle_err(state)
141    enddo
142
143    ! Ensure mandatory dims exist
144    if (idx_lon_in < 0 .or. idx_lat_in < 0) error stop "Input is missing mandatory 'lon' or 'lat' dimension."
145    if (idx_time_in < 0) len_time_in = 1
146    if (idx_soil_in < 0) len_soil_in = 1
147
148    ! Allocate only the "lon‐first" buffers (max‐sized) once
149    allocate(buf1D_in(len_lon_in),buf1D_out(len_lon_out))
150    allocate(buf2D_in(len_lon_in,len_lat_in),buf2D_out(len_lon_out, len_lat_in))
151    allocate(buf3D_in(len_lon_in,len_lat_in,len_time_in),buf3D_out(len_lon_out,len_lat_in,len_time_in))
152    allocate(buf4D_in(len_lon_in,len_lat_in,len_soil_in,len_time_in),buf4D_out(len_lon_out,len_lat_in,len_soil_in,len_time_in))
153
154    ! Get all variable IDs
155    state = nf90_inq_varids(ncid_in,nvars,varids_in)
156    if (state /= nf90_noerr) call handle_err(state)
157
158    ! Loop over each variable to define it in the output
159    do i = 1,nvars
160        ! Inquire name, xtype, ndims, natts
161        state = nf90_inquire_variable(ncid_in,varids_in(i),name = varName,xtype = xtypeVar,ndims = numDimsVar,natts = numAttsVar)
162        if (state /= nf90_noerr) call handle_err(state)
163        write(*,*) 'Treatment of '//varName
164
165        allocate(dimids_var_in(numDimsVar))
166        state = nf90_inquire_variable(ncid_in,varids_in(i),name = varName,xtype = xtypeVar,ndims = numDimsVar,dimids = dimids_var_in,natts = numAttsVar)
167        if (state /= nf90_noerr) call handle_err(state)
168
169        ! Detect if this variable first dimension is "lon"
170        if (numDimsVar >= 1 .and. dimids_var_in(1) == dimids_in(idx_lon_in)) then
171            uses_lon_first = .true.
172        else
173            uses_lon_first = .false.
174        endif
175
176        ! Build the output‐dimids list: replace the first dim with the output lon if needed
177        if (uses_lon_first) dimids_var_in(1) = dimids_out(idx_lon_in)
178        do j = 2,numDimsVar
179        ! Map each subsequent input dim to its output dim
180            do k = 1,ndims
181                if (dimids_var_in(j) == dimids_in(k)) then
182                    dimids_var_in(j) = dimids_out(k)
183                    exit
184                endif
185            enddo
186        enddo
187
188        ! Define this variable (same name, same xtype, but new dimids)
189        state = nf90_def_var(ncid_out,trim(varName),xtypeVar,dimids_var_in,varids_out(i))
190        if (state /= nf90_noerr) call handle_err(state)
191
192        deallocate(dimids_var_in)
193    enddo
194
195    ! Done defining all dims and vars exit define mode exactly once
196    state = nf90_enddef(ncid_out)
197    if (state /= nf90_noerr) call handle_err(state)
198
199    ! Loop over each variable to read from input and write to output
200    do i = 1,nvars
201        ! Re‐inquire metadata so we know dimids_var_in and numDimsVar
202        state = nf90_inquire_variable(ncid_in,varids_in(i),name = varName,xtype = xtypeVar,ndims = numDimsVar,natts = numAttsVar)
203        if (state /= nf90_noerr) call handle_err(state)
204
205        allocate(dimids_var_in(numDimsVar))
206        state = nf90_inquire_variable(ncid_in, varids_in(i),name = varName,xtype = xtypeVar,ndims = numDimsVar,dimids = dimids_var_in,natts = numAttsVar)
207        if (state /= nf90_noerr) call handle_err(state)
208
209        ! Detect again if first dim = lon
210        if (numDimsVar >= 1 .and. dimids_var_in(1) == dimids_in(idx_lon_in)) then
211            uses_lon_first = .true.
212        else
213            uses_lon_first = .false.
214        endif
215
216        select case (numDimsVar)
217            case (1)
218                if (uses_lon_first) then
219                    ! 1D lon sequence: read len_lon_in, extend to len_lon_out
220                    state = nf90_get_var(ncid_in,varids_in(i),buf1D_in)
221                    if (state /= nf90_noerr) call handle_err(state)
222
223                    buf1D_out(1:len_lon_in) = buf1D_in(1:len_lon_in)
224                    buf1D_out(len_lon_out)  = buf1D_in(1)  ! repeat first lon at end
225
226                    state = nf90_put_var(ncid_out,varids_out(i),buf1D_out)
227                    if (state /= nf90_noerr) call handle_err(state)
228
229                else
230                    ! Some other 1D (e.g. lat or time). Allocate exact 1D temp:
231                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(1),tmpDimName, thisLen)
232                    if (state /= nf90_noerr) call handle_err(state)
233
234                    allocate(tmp1D(thisLen))
235                    state = nf90_get_var(ncid_in,varids_in(i),tmp1D(1:thisLen))
236                    if (state /= nf90_noerr) call handle_err(state)
237
238                    state = nf90_put_var(ncid_out,varids_out(i),tmp1D(1:thisLen))
239                    if (state /= nf90_noerr) call handle_err(state)
240
241                    deallocate(tmp1D)
242                endif
243
244            case (2)
245                if (uses_lon_first) then
246                    ! 2D with first dim = lon (len_lon_in × lenDim2)
247                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(2),tmpDimName,lenDim2)
248                    if (state /= nf90_noerr) call handle_err(state)
249
250                    state = nf90_get_var(ncid_in,varids_in(i),buf2D_in(1:len_lon_in,1:lenDim2))
251                    if (state /= nf90_noerr) call handle_err(state)
252
253                    buf2D_out(1:len_lon_in,1:lenDim2) = buf2D_in(1:len_lon_in,1:lenDim2)
254                    buf2D_out(len_lon_out,1:lenDim2) = buf2D_in(1,1:lenDim2)
255
256                    state = nf90_put_var(ncid_out,varids_out(i),buf2D_out(1:len_lon_out,1:lenDim2))
257                    if (state /= nf90_noerr) call handle_err(state)
258
259                else
260                    ! Some other 2D (no lon‐extension). Allocate exact 2D temp:
261                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(1),tmpDimName,len1)
262                    if (state /= nf90_noerr) call handle_err(state)
263                    state = nf90_inquire_dimension(ncid_in, dimids_var_in(2),tmpDimName,len2)
264                    if (state /= nf90_noerr) call handle_err(state)
265
266                    allocate(tmp2D(len1,len2))
267                    state = nf90_get_var(ncid_in,varids_in(i),tmp2D(1:len1,1:len2))
268                    if (state /= nf90_noerr) call handle_err(state)
269
270                    state = nf90_put_var(ncid_out, varids_out(i), tmp2D(1:len1,1:len2))
271                    if (state /= nf90_noerr) call handle_err(state)
272
273                    deallocate(tmp2D)
274                endif
275
276            case (3)
277                if (uses_lon_first) then
278                    ! 3D with first dim = lon (len_lon_in × lenDim2 × lenDim3)
279                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(2),tmpDimName,lenDim2)
280                    if (state /= nf90_noerr) call handle_err(state)
281                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(3),tmpDimName,lenDim3)
282                    if (state /= nf90_noerr) call handle_err(state)
283
284                    state = nf90_get_var(ncid_in,varids_in(i),buf3D_in(1:len_lon_in,1:lenDim2,1:lenDim3))
285                    if (state /= nf90_noerr) call handle_err(state)
286
287                    buf3D_out(1:len_lon_in,1:lenDim2,1:lenDim3) = buf3D_in(1:len_lon_in,1:lenDim2,1:lenDim3)
288                    buf3D_out(len_lon_out,1:lenDim2,1:lenDim3) = buf3D_in(1,1:lenDim2,1:lenDim3)
289
290                    state = nf90_put_var(ncid_out,varids_out(i),buf3D_out(1:len_lon_out,1:lenDim2,1:lenDim3))
291                    if (state /= nf90_noerr) call handle_err(state)
292
293                else
294                    ! Some other 3D (no lon‐extension). Allocate exact 3D temp:
295                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(1),tmpDimName,len1)
296                    if (state /= nf90_noerr) call handle_err(state)
297                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(2),tmpDimName,len2)
298                    if (state /= nf90_noerr) call handle_err(state)
299                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(3),tmpDimName,len3)
300                    if (state /= nf90_noerr) call handle_err(state)
301
302                    allocate(tmp3D(len1,len2,len3))
303                    state = nf90_get_var(ncid_in,varids_in(i),tmp3D(1:len1,1:len2,1:len3))
304                    if (state /= nf90_noerr) call handle_err(state)
305
306                    state = nf90_put_var(ncid_out,varids_out(i),tmp3D(1:len1,1:len2,1:len3))
307                    if (state /= nf90_noerr) call handle_err(state)
308
309                    deallocate(tmp3D)
310                endif
311
312            case (4)
313                if (uses_lon_first) then ! 4D with first dim = lon (len_lon_in × lenDim2 × lenDim3 × lenDim4)
314                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(2),tmpDimName,lenDim2)
315                    if (state /= nf90_noerr) call handle_err(state)
316                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(3),tmpDimName,lenDim3)
317                    if (state /= nf90_noerr) call handle_err(state)
318                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(4),tmpDimName,lenDim4)
319                    if (state /= nf90_noerr) call handle_err(state)
320
321                    state = nf90_get_var(ncid_in,varids_in(i),buf4D_in(1:len_lon_in,1:lenDim2,1:lenDim3,1:lenDim4))
322                    if (state /= nf90_noerr) call handle_err(state)
323
324                    buf4D_out(1:len_lon_in,1:lenDim2,1:lenDim3,1:lenDim4) = buf4D_in(1:len_lon_in, 1:lenDim2,1:lenDim3,1:lenDim4)
325                    buf4D_out(len_lon_out,1:lenDim2,1:lenDim3,1:lenDim4) = buf4D_in(1,1:lenDim2,1:lenDim3,1:lenDim4)
326
327                    state = nf90_put_var(ncid_out,varids_out(i),buf4D_out(1:len_lon_out,1:lenDim2,1:lenDim3,1:lenDim4))
328                    if (state /= nf90_noerr) call handle_err(state)
329
330                else ! Some other 4D (no lon‐extension). Allocate exact 4D temp:
331                    state = nf90_inquire_dimension(ncid_in, dimids_var_in(1),tmpDimName,len1)
332                    if (state /= nf90_noerr) call handle_err(state)
333                    state = nf90_inquire_dimension(ncid_in, dimids_var_in(2),tmpDimName,len2)
334                    if (state /= nf90_noerr) call handle_err(state)
335                    state = nf90_inquire_dimension(ncid_in, dimids_var_in(3),tmpDimName,len3)
336                    if (state /= nf90_noerr) call handle_err(state)
337                    state = nf90_inquire_dimension(ncid_in, dimids_var_in(4),tmpDimName,len4)
338                    if (state /= nf90_noerr) call handle_err(state)
339
340                    allocate(tmp4D(len1,len2,len3,len4))
341                    state = nf90_get_var(ncid_in,varids_in(i),tmp4D(1:len1,1:len2,1:len3,1:len4))
342                    if (state /= nf90_noerr) call handle_err(state)
343
344                    state = nf90_put_var(ncid_out,varids_out(i),tmp4D(1:len1,1:len2,1:len3,1:len4))
345                    if (state /= nf90_noerr) call handle_err(state)
346
347                    deallocate(tmp4D)
348                endif
349
350            case default
351                cycle ! Skip variables with 0 dims
352        end select
353
354        deallocate(dimids_var_in)
355    enddo
356
357    ! Close both NetCDF files
358    state = nf90_close(ncid_in)
359    if (state /= nf90_noerr) call handle_err(state)
360    state = nf90_close(ncid_out)
361    if (state /= nf90_noerr) call handle_err(state)
362
363    ! Deallocate everything
364    deallocate(dimids_in,dimids_out,varids_in,varids_out,dimNames,dimLens)
365    deallocate(buf1D_in,buf1D_out,buf2D_in,buf2D_out,buf3D_in,buf3D_out,buf4D_in,buf4D_out)
366
367    write(*,*) "> ""data2reshape_Y"//trim(str)//".nc"" processed!"
368enddo
369
370END PROGRAM reshape_XIOS_output
Note: See TracBrowser for help on using the repository browser.