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

Last change on this file since 3792 was 3786, checked in by jbclement, 6 months 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
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 version_info_mod, only: print_version_info
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, 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"
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
86! Main loop: two PCM years
87do numyear = 1,2
88    write(str,'(I1.1)') numyear
89    write(*,*) "> Reshaping variables from ""data2reshape_Y"//trim(str)//".nc""..."
90
91    ! Open input file (read‐only)
92    state = nf90_open("data2reshape_Y"//trim(str)//".nc",mode = nf90_nowrite,ncid = ncid_in)
93    if (state /= nf90_noerr) call handle_err(state)
94
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)
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
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)
108    if (state /= nf90_noerr) call handle_err(state)
109
110    ! Inquire input for dims, vars, global atts, unlimited dim ID
111    state = nf90_inquire(ncid_in,ndims,nvars,nGlobalAtts,unlimDimID)
112    if (state /= nf90_noerr) call handle_err(state)
113
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))
116
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)
119    if (state /= nf90_noerr) call handle_err(state)
120
121    do i = 1,ndims
122        state = nf90_inquire_dimension(ncid_in,dimids_in(i),dimNames(i),dimLens(i))
123        if (state /= nf90_noerr) call handle_err(state)
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))
148        endif
149        if (state /= nf90_noerr) call handle_err(state)
150    enddo
151
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
168    do i = 1,nvars
169        ! Inquire name, xtype, ndims, natts
170        state = nf90_inquire_variable(ncid_in,varids_in(i),name = varName,xtype = xtypeVar,ndims = numDimsVar,natts = numAttsVar)
171        if (state /= nf90_noerr) call handle_err(state)
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)
176        if (state /= nf90_noerr) call handle_err(state)
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.
183        endif
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)
202    enddo
203
204    ! Done defining all dims and vars exit define mode exactly once
205    state = nf90_enddef(ncid_out)
206    if (state /= nf90_noerr) call handle_err(state)
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)
368    if (state /= nf90_noerr) call handle_err(state)
369    state = nf90_close(ncid_out)
370    if (state /= nf90_noerr) call handle_err(state)
371
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!"
377enddo
378
379END PROGRAM reshape_XIOS_output
Note: See TracBrowser for help on using the repository browser.