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

Last change on this file since 3790 was 3786, checked in by jbclement, 5 weeks ago

PEM:
Optimization (computing time is devided by 2.2) of the program "reshape_XIOS_output" to convert XIOS output onto the PCM grid.
JBC
File size: 16.7 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
[3579]13use version_info_mod, only: print_version_info
[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
41character(100)                     :: varName, arg
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! Handle command‐line argument "version"
[3579]75if (command_argument_count() > 0) then ! Get the number of command-line arguments
76    call get_command_argument(1,arg) ! Read the argument given to the program
77    select case (trim(adjustl(arg)))
78        case('version')
79            call print_version_info()
80            stop
81        case default
82            error stop 'The argument given to the program is unknown!'
83    end select
84endif
85
[3786]86! Main loop: two PCM years
[3076]87do numyear = 1,2
[3786]88    write(str,'(I1.1)') numyear
89    write(*,*) "> Reshaping variables from ""data2reshape_Y"//trim(str)//".nc""..."
[3181]90
[3786]91    ! Open input file (read‐only)
92    state = nf90_open("data2reshape_Y"//trim(str)//".nc",mode = nf90_nowrite,ncid = ncid_in)
[3076]93    if (state /= nf90_noerr) call handle_err(state)
[2850]94
[3786]95    ! If output exists, delete it
96    inquire(file = "data_PCM_Y"//trim(str)//".nc",exist = exists)
97    if (exists) then
98        call execute_command_line("rm data_PCM_Y"//trim(str)//".nc",cmdstat = cstat)
[3106]99        if (cstat > 0) then
100            error stop 'Command exection failed!'
101        else if (cstat < 0) then
102            error stop 'Command execution not supported!'
103        endif
104    endif
[3786]105
106    ! Create output file in define mode
107    state = nf90_create("data_PCM_Y"//trim(str)//".nc",cmode = or(nf90_noclobber,nf90_64bit_offset),ncid = ncid_out)
[3076]108    if (state /= nf90_noerr) call handle_err(state)
[2850]109
[3786]110    ! Inquire input for dims, vars, global atts, unlimited dim ID
111    state = nf90_inquire(ncid_in,ndims,nvars,nGlobalAtts,unlimDimID)
[3076]112    if (state /= nf90_noerr) call handle_err(state)
[2850]113
[3786]114    ! Allocate arrays for dim IDs, var IDs, names, lengths
115    allocate(dimids_in(ndims),varids_in(nvars),dimids_out(ndims),varids_out(nvars),dimNames(ndims),dimLens(ndims))
[2850]116
[3786]117    ! Get the dimension IDs and then query each for its name and length
118    state = nf90_inq_dimids(ncid_in,ndims,dimids_in,unlimDimID)
[3076]119    if (state /= nf90_noerr) call handle_err(state)
[2850]120
[3076]121    do i = 1,ndims
[3786]122        state = nf90_inquire_dimension(ncid_in,dimids_in(i),dimNames(i),dimLens(i))
[3076]123        if (state /= nf90_noerr) call handle_err(state)
[3786]124
125        select case (trim(dimNames(i)))
126            case ("lon","longitude")
127                idx_lon_in = i
128                len_lon_in = dimLens(i)
129            case ("lat","latitude")
130                idx_lat_in = i
131                len_lat_in = dimLens(i)
132            case ("time_counter","Time")
133                idx_time_in = i
134                len_time_in = dimLens(i)
135            case ("soil_layers","subsurface_layers")
136                idx_soil_in = i
137                len_soil_in = dimLens(i)
138            case default
139                ! nothing special
140        end select
141
142        ! Define the same dimension in the output, except lon becoming (len_lon_in + 1)
143        if (i == idx_lon_in) then
144            len_lon_out = len_lon_in + 1
145            state = nf90_def_dim(ncid_out,trim(dimNames(i)),len_lon_out,dimids_out(i))
146        else
147            state = nf90_def_dim(ncid_out,trim(dimNames(i)),dimLens(i),dimids_out(i))
[3076]148        endif
149        if (state /= nf90_noerr) call handle_err(state)
150    enddo
[2850]151
[3786]152    ! Ensure mandatory dims exist
153    if (idx_lon_in < 0 .or. idx_lat_in < 0) error stop "Input is missing mandatory 'lon' or 'lat' dimension."
154    if (idx_time_in < 0) len_time_in = 1
155    if (idx_soil_in < 0) len_soil_in = 1
156
157    ! Allocate only the "lon‐first" buffers (max‐sized) once
158    allocate(buf1D_in(len_lon_in),buf1D_out(len_lon_out))
159    allocate(buf2D_in(len_lon_in,len_lat_in),buf2D_out(len_lon_out, len_lat_in))
160    allocate(buf3D_in(len_lon_in,len_lat_in,len_time_in),buf3D_out(len_lon_out,len_lat_in,len_time_in))
161    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))
162
163    ! Get all variable IDs
164    state = nf90_inq_varids(ncid_in,nvars,varids_in)
165    if (state /= nf90_noerr) call handle_err(state)
166
167    ! Loop over each variable to define it in the output
[3076]168    do i = 1,nvars
[3786]169        ! Inquire name, xtype, ndims, natts
170        state = nf90_inquire_variable(ncid_in,varids_in(i),name = varName,xtype = xtypeVar,ndims = numDimsVar,natts = numAttsVar)
[3076]171        if (state /= nf90_noerr) call handle_err(state)
[3786]172        write(*,*) 'Treatment of '//varName
173
174        allocate(dimids_var_in(numDimsVar))
175        state = nf90_inquire_variable(ncid_in,varids_in(i),name = varName,xtype = xtypeVar,ndims = numDimsVar,dimids = dimids_var_in,natts = numAttsVar)
[3076]176        if (state /= nf90_noerr) call handle_err(state)
[3786]177
178        ! Detect if this variable first dimension is "lon"
179        if (numDimsVar >= 1 .and. dimids_var_in(1) == dimids_in(idx_lon_in)) then
180            uses_lon_first = .true.
181        else
182            uses_lon_first = .false.
[3076]183        endif
[3786]184
185        ! Build the output‐dimids list: replace the first dim with the output lon if needed
186        if (uses_lon_first) dimids_var_in(1) = dimids_out(idx_lon_in)
187        do j = 2,numDimsVar
188        ! Map each subsequent input dim to its output dim
189            do k = 1,ndims
190                if (dimids_var_in(j) == dimids_in(k)) then
191                    dimids_var_in(j) = dimids_out(k)
192                    exit
193                endif
194            enddo
195        enddo
196
197        ! Define this variable (same name, same xtype, but new dimids)
198        state = nf90_def_var(ncid_out,trim(varName),xtypeVar,dimids_var_in,varids_out(i))
199        if (state /= nf90_noerr) call handle_err(state)
200
201        deallocate(dimids_var_in)
[3076]202    enddo
[2850]203
[3786]204    ! Done defining all dims and vars exit define mode exactly once
205    state = nf90_enddef(ncid_out)
[3076]206    if (state /= nf90_noerr) call handle_err(state)
[3786]207
208    ! Loop over each variable to read from input and write to output
209    do i = 1,nvars
210        ! Re‐inquire metadata so we know dimids_var_in and numDimsVar
211        state = nf90_inquire_variable(ncid_in,varids_in(i),name = varName,xtype = xtypeVar,ndims = numDimsVar,natts = numAttsVar)
212        if (state /= nf90_noerr) call handle_err(state)
213
214        allocate(dimids_var_in(numDimsVar))
215        state = nf90_inquire_variable(ncid_in, varids_in(i),name = varName,xtype = xtypeVar,ndims = numDimsVar,dimids = dimids_var_in,natts = numAttsVar)
216        if (state /= nf90_noerr) call handle_err(state)
217
218        ! Detect again if first dim = lon
219        if (numDimsVar >= 1 .and. dimids_var_in(1) == dimids_in(idx_lon_in)) then
220            uses_lon_first = .true.
221        else
222            uses_lon_first = .false.
223        endif
224
225        select case (numDimsVar)
226            case (1)
227                if (uses_lon_first) then
228                    ! 1D lon sequence: read len_lon_in, extend to len_lon_out
229                    state = nf90_get_var(ncid_in,varids_in(i),buf1D_in)
230                    if (state /= nf90_noerr) call handle_err(state)
231
232                    buf1D_out(1:len_lon_in) = buf1D_in(1:len_lon_in)
233                    buf1D_out(len_lon_out)  = buf1D_in(1)  ! repeat first lon at end
234
235                    state = nf90_put_var(ncid_out,varids_out(i),buf1D_out)
236                    if (state /= nf90_noerr) call handle_err(state)
237
238                else
239                    ! Some other 1D (e.g. lat or time). Allocate exact 1D temp:
240                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(1),tmpDimName, thisLen)
241                    if (state /= nf90_noerr) call handle_err(state)
242
243                    allocate(tmp1D(thisLen))
244                    state = nf90_get_var(ncid_in,varids_in(i),tmp1D(1:thisLen))
245                    if (state /= nf90_noerr) call handle_err(state)
246
247                    state = nf90_put_var(ncid_out,varids_out(i),tmp1D(1:thisLen))
248                    if (state /= nf90_noerr) call handle_err(state)
249
250                    deallocate(tmp1D)
251                endif
252
253            case (2)
254                if (uses_lon_first) then
255                    ! 2D with first dim = lon (len_lon_in × lenDim2)
256                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(2),tmpDimName,lenDim2)
257                    if (state /= nf90_noerr) call handle_err(state)
258
259                    state = nf90_get_var(ncid_in,varids_in(i),buf2D_in(1:len_lon_in,1:lenDim2))
260                    if (state /= nf90_noerr) call handle_err(state)
261
262                    buf2D_out(1:len_lon_in,1:lenDim2) = buf2D_in(1:len_lon_in,1:lenDim2)
263                    buf2D_out(len_lon_out,1:lenDim2) = buf2D_in(1,1:lenDim2)
264
265                    state = nf90_put_var(ncid_out,varids_out(i),buf2D_out(1:len_lon_out,1:lenDim2))
266                    if (state /= nf90_noerr) call handle_err(state)
267
268                else
269                    ! Some other 2D (no lon‐extension). Allocate exact 2D temp:
270                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(1),tmpDimName,len1)
271                    if (state /= nf90_noerr) call handle_err(state)
272                    state = nf90_inquire_dimension(ncid_in, dimids_var_in(2),tmpDimName,len2)
273                    if (state /= nf90_noerr) call handle_err(state)
274
275                    allocate(tmp2D(len1,len2))
276                    state = nf90_get_var(ncid_in,varids_in(i),tmp2D(1:len1,1:len2))
277                    if (state /= nf90_noerr) call handle_err(state)
278
279                    state = nf90_put_var(ncid_out, varids_out(i), tmp2D(1:len1,1:len2))
280                    if (state /= nf90_noerr) call handle_err(state)
281
282                    deallocate(tmp2D)
283                endif
284
285            case (3)
286                if (uses_lon_first) then
287                    ! 3D with first dim = lon (len_lon_in × lenDim2 × lenDim3)
288                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(2),tmpDimName,lenDim2)
289                    if (state /= nf90_noerr) call handle_err(state)
290                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(3),tmpDimName,lenDim3)
291                    if (state /= nf90_noerr) call handle_err(state)
292
293                    state = nf90_get_var(ncid_in,varids_in(i),buf3D_in(1:len_lon_in,1:lenDim2,1:lenDim3))
294                    if (state /= nf90_noerr) call handle_err(state)
295
296                    buf3D_out(1:len_lon_in,1:lenDim2,1:lenDim3) = buf3D_in(1:len_lon_in,1:lenDim2,1:lenDim3)
297                    buf3D_out(len_lon_out,1:lenDim2,1:lenDim3) = buf3D_in(1,1:lenDim2,1:lenDim3)
298
299                    state = nf90_put_var(ncid_out,varids_out(i),buf3D_out(1:len_lon_out,1:lenDim2,1:lenDim3))
300                    if (state /= nf90_noerr) call handle_err(state)
301
302                else
303                    ! Some other 3D (no lon‐extension). Allocate exact 3D temp:
304                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(1),tmpDimName,len1)
305                    if (state /= nf90_noerr) call handle_err(state)
306                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(2),tmpDimName,len2)
307                    if (state /= nf90_noerr) call handle_err(state)
308                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(3),tmpDimName,len3)
309                    if (state /= nf90_noerr) call handle_err(state)
310
311                    allocate(tmp3D(len1,len2,len3))
312                    state = nf90_get_var(ncid_in,varids_in(i),tmp3D(1:len1,1:len2,1:len3))
313                    if (state /= nf90_noerr) call handle_err(state)
314
315                    state = nf90_put_var(ncid_out,varids_out(i),tmp3D(1:len1,1:len2,1:len3))
316                    if (state /= nf90_noerr) call handle_err(state)
317
318                    deallocate(tmp3D)
319                endif
320
321            case (4)
322                if (uses_lon_first) then ! 4D with first dim = lon (len_lon_in × lenDim2 × lenDim3 × lenDim4)
323                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(2),tmpDimName,lenDim2)
324                    if (state /= nf90_noerr) call handle_err(state)
325                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(3),tmpDimName,lenDim3)
326                    if (state /= nf90_noerr) call handle_err(state)
327                    state = nf90_inquire_dimension(ncid_in,dimids_var_in(4),tmpDimName,lenDim4)
328                    if (state /= nf90_noerr) call handle_err(state)
329
330                    state = nf90_get_var(ncid_in,varids_in(i),buf4D_in(1:len_lon_in,1:lenDim2,1:lenDim3,1:lenDim4))
331                    if (state /= nf90_noerr) call handle_err(state)
332
333                    buf4D_out(1:len_lon_in,1:lenDim2,1:lenDim3,1:lenDim4) = buf4D_in(1:len_lon_in, 1:lenDim2,1:lenDim3,1:lenDim4)
334                    buf4D_out(len_lon_out,1:lenDim2,1:lenDim3,1:lenDim4) = buf4D_in(1,1:lenDim2,1:lenDim3,1:lenDim4)
335
336                    state = nf90_put_var(ncid_out,varids_out(i),buf4D_out(1:len_lon_out,1:lenDim2,1:lenDim3,1:lenDim4))
337                    if (state /= nf90_noerr) call handle_err(state)
338
339                else ! Some other 4D (no lon‐extension). Allocate exact 4D temp:
340                    state = nf90_inquire_dimension(ncid_in, dimids_var_in(1),tmpDimName,len1)
341                    if (state /= nf90_noerr) call handle_err(state)
342                    state = nf90_inquire_dimension(ncid_in, dimids_var_in(2),tmpDimName,len2)
343                    if (state /= nf90_noerr) call handle_err(state)
344                    state = nf90_inquire_dimension(ncid_in, dimids_var_in(3),tmpDimName,len3)
345                    if (state /= nf90_noerr) call handle_err(state)
346                    state = nf90_inquire_dimension(ncid_in, dimids_var_in(4),tmpDimName,len4)
347                    if (state /= nf90_noerr) call handle_err(state)
348
349                    allocate(tmp4D(len1,len2,len3,len4))
350                    state = nf90_get_var(ncid_in,varids_in(i),tmp4D(1:len1,1:len2,1:len3,1:len4))
351                    if (state /= nf90_noerr) call handle_err(state)
352
353                    state = nf90_put_var(ncid_out,varids_out(i),tmp4D(1:len1,1:len2,1:len3,1:len4))
354                    if (state /= nf90_noerr) call handle_err(state)
355
356                    deallocate(tmp4D)
357                endif
358
359            case default
360                cycle ! Skip variables with 0 dims
361        end select
362
363        deallocate(dimids_var_in)
364    enddo
365
366    ! Close both NetCDF files
367    state = nf90_close(ncid_in)
[3076]368    if (state /= nf90_noerr) call handle_err(state)
[3786]369    state = nf90_close(ncid_out)
[3076]370    if (state /= nf90_noerr) call handle_err(state)
[2850]371
[3786]372    ! Deallocate everything
373    deallocate(dimids_in,dimids_out,varids_in,varids_out,dimNames,dimLens)
374    deallocate(buf1D_in,buf1D_out,buf2D_in,buf2D_out,buf3D_in,buf3D_out,buf4D_in,buf4D_out)
375
376    write(*,*) "> ""data2reshape_Y"//trim(str)//".nc"" processed!"
[2850]377enddo
378
[3076]379END PROGRAM reshape_XIOS_output
Note: See TracBrowser for help on using the repository browser.