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

Last change on this file since 3836 was 3836, checked in by jbclement, 3 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
RevLine 
[3076]1PROGRAM reshape_XIOS_output
[2855]2
3!=======================================================================
[3786]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.
[2855]8!
[3786]9! Authors: RV & LL (original), JBC (optimized)
[2855]10!=======================================================================
[3786]11
[3076]12use netcdf
[3836]13use parse_args_mod, only: parse_args
[2850]14
[3076]15implicit none
[2850]16
[3786]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
[2850]23
[3786]24! Store each input dimension name and length
25character(30), allocatable, dimension(:) :: dimNames
26integer,       allocatable, dimension(:) :: dimLens
[3579]27
[3786]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
[3836]41character(100)                     :: varName
[3786]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
[3836]74! Parse command-line options
75call parse_args()
[3579]76
[3786]77! Main loop: two PCM years
[3076]78do numyear = 1,2
[3786]79    write(str,'(I1.1)') numyear
80    write(*,*) "> Reshaping variables from ""data2reshape_Y"//trim(str)//".nc""..."
[3181]81
[3786]82    ! Open input file (read‐only)
83    state = nf90_open("data2reshape_Y"//trim(str)//".nc",mode = nf90_nowrite,ncid = ncid_in)
[3076]84    if (state /= nf90_noerr) call handle_err(state)
[2850]85
[3786]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)
[3106]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
[3786]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)
[3076]99    if (state /= nf90_noerr) call handle_err(state)
[2850]100
[3786]101    ! Inquire input for dims, vars, global atts, unlimited dim ID
102    state = nf90_inquire(ncid_in,ndims,nvars,nGlobalAtts,unlimDimID)
[3076]103    if (state /= nf90_noerr) call handle_err(state)
[2850]104
[3786]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))
[2850]107
[3786]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)
[3076]110    if (state /= nf90_noerr) call handle_err(state)
[2850]111
[3076]112    do i = 1,ndims
[3786]113        state = nf90_inquire_dimension(ncid_in,dimids_in(i),dimNames(i),dimLens(i))
[3076]114        if (state /= nf90_noerr) call handle_err(state)
[3786]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))
[3076]139        endif
140        if (state /= nf90_noerr) call handle_err(state)
141    enddo
[2850]142
[3786]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
[3076]159    do i = 1,nvars
[3786]160        ! Inquire name, xtype, ndims, natts
161        state = nf90_inquire_variable(ncid_in,varids_in(i),name = varName,xtype = xtypeVar,ndims = numDimsVar,natts = numAttsVar)
[3076]162        if (state /= nf90_noerr) call handle_err(state)
[3786]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)
[3076]167        if (state /= nf90_noerr) call handle_err(state)
[3786]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.
[3076]174        endif
[3786]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)
[3076]193    enddo
[2850]194
[3786]195    ! Done defining all dims and vars exit define mode exactly once
196    state = nf90_enddef(ncid_out)
[3076]197    if (state /= nf90_noerr) call handle_err(state)
[3786]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)
[3076]359    if (state /= nf90_noerr) call handle_err(state)
[3786]360    state = nf90_close(ncid_out)
[3076]361    if (state /= nf90_noerr) call handle_err(state)
[2850]362
[3786]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!"
[2850]368enddo
369
[3076]370END PROGRAM reshape_XIOS_output
Note: See TracBrowser for help on using the repository browser.