source: trunk/LMDZ.COMMON/libf/evolution/xios_data.F90

Last change on this file was 3991, checked in by jbclement, 4 weeks ago

PEM:
Apply documentation template everywhere: standardized headers format with short description, separators between functions/subroutines, normalized code sections, aligned dependencies/arguments/variables declaration.
JBC
File size: 15.5 KB
Line 
1MODULE xios_data
2!-----------------------------------------------------------------------
3! NAME
4!     xios_data
5!
6! DESCRIPTION
7!     Read XIOS output data and process it for PEM initialization.
8!
9! AUTHORS & DATE
10!     JB Clement, 2025
11!
12! NOTES
13!
14!-----------------------------------------------------------------------
15
16! DEPENDENCIES
17! ------------
18use netcdf, only: nf90_open, nf90_close, nf90_inquire_dimension, nf90_inq_dimid, nf90_noerr, nf90_nowrite, nf90_get_var, nf90_inq_varid
19
20! DECLARATION
21! -----------
22implicit none
23
24! MODULE VARIABLES
25! ----------------
26character(19), parameter :: file1_daily  = "Xoutdaily4pem_Y1.nc"   ! XIOS daily output file, year 1
27character(19), parameter :: file2_daily  = "Xoutdaily4pem_Y2.nc"   ! XIOS daily output file, year 2
28character(20), parameter :: file1_yearly = "Xoutyearly4pem_Y1.nc"  ! XIOS yearly output file, year 1
29character(20), parameter :: file2_yearly = "Xoutyearly4pem_Y2.nc"  ! XIOS yearly output file, year 2
30character(256)           :: msg                                    ! Message for reading errors
31integer                  :: fID, vID                               ! File and variable IDs for reading
32
33! INTERFACES
34! ----------
35interface get_var
36    module procedure get_var_1d, get_var_2d, get_var_3d, get_var_4d
37end interface get_var
38
39contains
40!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
41
42!=======================================================================
43SUBROUTINE load_xios_data(ngrid,nslope,nsoil_PCM,nsol,h2ofrost_PCM,co2frost_PCM,ps_avg,tsurf_avg,tsurf_avg_y1,tsoil_avg,tsoil_ts,watersurf_density_avg,d_h2oice,d_co2ice, &
44                         ps_ts,q_h2o_ts,q_co2_ts,watersoil_density_ts,min_h2oice,min_co2ice)
45!-----------------------------------------------------------------------
46! NAME
47!     load_xios_data
48!
49! DESCRIPTION
50!     Reads yearly and daily XIOS data, computes frost and ice tendencies.
51!
52! AUTHORS & DATE
53!     JB Clement, 2025
54!
55! NOTES
56!
57!-----------------------------------------------------------------------
58
59! DEPENDENCIES
60! ------------
61use grid_conversion, only: lonlat2vect
62use soil,            only: do_soil
63use tendencies,      only: compute_tend
64use metamorphism,    only: compute_frost
65
66! DECLARATION
67! -----------
68implicit none
69
70! ARGUMENTS
71! ---------
72integer,                                      intent(in)  :: ngrid, nslope, nsoil_PCM, nsol ! Grid dimensions
73real, dimension(ngrid,nslope),                intent(in)  :: h2ofrost_PCM, co2frost_PCM     ! PCM frost fields
74real, dimension(ngrid),                       intent(out) :: ps_avg                         ! Average surface pressure
75real, dimension(ngrid,nslope),                intent(out) :: tsurf_avg, tsurf_avg_y1        ! Surface temperature
76real, dimension(ngrid,nslope),                intent(out) :: watersurf_density_avg          ! Water density
77real, dimension(ngrid,nslope),                intent(out) :: d_h2oice, d_co2ice             ! Ice tendencies
78real, dimension(ngrid,nslope),                intent(out) :: min_h2oice, min_co2ice         ! Ice minima
79real, dimension(ngrid,nsoil_PCM,nslope),      intent(out) :: tsoil_avg                      ! Soil temperature
80real, dimension(ngrid,nsol),                  intent(out) :: ps_ts, q_h2o_ts, q_co2_ts      ! Time series
81real, dimension(ngrid,nsoil_PCM,nslope,nsol), intent(out) :: tsoil_ts, watersoil_density_ts ! Soil time series
82
83! LOCAL VARIABLES
84! ---------------
85integer                               :: islope, isoil, isol, nlon, nlat
86real, dimension(:,:),     allocatable :: var_read_2d
87real, dimension(:,:,:),   allocatable :: var_read_3d
88real, dimension(:,:,:,:), allocatable :: var_read_4d
89character(:),             allocatable :: num ! For reading slope variables
90real, dimension(ngrid,nslope,2)       :: min_h2operice, min_co2perice, min_h2ofrost, min_co2frost
91
92! CODE
93! ----
94! Initialization
95min_h2operice = 0.
96min_co2perice = 0.
97min_h2ofrost = 0.
98min_co2frost = 0.
99if (nslope == 1) then ! No slope
100    allocate(character(0) :: num)
101else ! Using slopes
102    allocate(character(8) :: num)
103endif
104
105!~~~~~~~~~~~~~~~~~~~~~~~~ Year 1 - Yearly data ~~~~~~~~~~~~~~~~~~~~~~~~~
106! Open the NetCDF file of XIOS outputs
107write(*,*) "Opening "//file1_yearly//"..."
108call error_msg(nf90_open(file1_yearly,nf90_nowrite,fID),'open',file1_yearly)
109
110! Get the dimensions
111call error_msg(nf90_inq_dimid(fID,'lon',vID),'inq',file1_yearly)
112call error_msg(nf90_inquire_dimension(fID,vID,len = nlon),'inq',file1_yearly)
113call error_msg(nf90_inq_dimid(fID,'lat',vID),'inq',file1_yearly)
114call error_msg(nf90_inquire_dimension(fID,vID,len = nlat),'inq',file1_yearly)
115
116! Allocate and read the variables
117allocate(var_read_2d(nlon,nlat),var_read_3d(nlon,nlat,nsoil_PCM))
118do islope = 1,nslope
119    if (nslope /= 1) then
120        num='  '
121        write(num,'(i2.2)') islope
122        num = '_slope'//num
123    endif
124    call get_var('co2ice'//num,var_read_2d)          ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2frost(:,islope,1))
125    call get_var('h2o_ice_s'//num,var_read_2d)       ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2ofrost(:,islope,1))
126    call get_var('watercap'//num,var_read_2d)        ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2operice(:,islope,1))
127    call get_var('perennial_co2ice'//num,var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2perice(:,islope,1))
128    call get_var('tsurf'//num,var_read_2d)           ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,tsurf_avg_y1(:,islope))
129enddo
130
131! Close the NetCDF file of XIOS outputs
132call error_msg(nf90_close(fID),"close",file1_yearly)
133write(*,*) '> Data from "'//file1_yearly//'" downloaded!'
134
135!~~~~~~~~~~~~~~~~~~~~~~~~ Year 2 - Yearly data ~~~~~~~~~~~~~~~~~~~~~~~~~
136! Open the NetCDF file of XIOS outputs
137write(*,*) "Opening "//file2_yearly//"..."
138call error_msg(nf90_open(file2_yearly,nf90_nowrite,fID),'open',file2_yearly)
139
140! Allocate and read the variables
141call get_var('ps',var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,ps_avg)
142do islope = 1,nslope
143    if (nslope /= 1) then
144        num='  '
145        write(num,'(i2.2)') islope
146        num = '_slope'//num
147    endif
148    call get_var('tsurf'//num,var_read_2d)           ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,tsurf_avg(:,islope))
149    call get_var('co2ice'//num,var_read_2d)          ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2frost(:,islope,2))
150    call get_var('h2o_ice_s'//num,var_read_2d)       ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2ofrost(:,islope,2))
151    call get_var('watercap'//num,var_read_2d)        ; call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_h2operice(:,islope,2))
152    call get_var('perennial_co2ice'//num,var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,min_co2perice(:,islope,2))
153    if (do_soil) then
154        call get_var('soiltemp'//num,var_read_3d)
155        do isoil = 1,nsoil_PCM
156            call lonlat2vect(nlon,nlat,ngrid,var_read_3d(:,:,isoil),tsoil_avg(:,isoil,islope))
157        enddo
158        call get_var('waterdensity_surface'//num,var_read_2d); call lonlat2vect(nlon,nlat,ngrid,var_read_2d,watersurf_density_avg(:,islope))
159    endif
160enddo
161! Close the NetCDF file of XIOS outputs
162call error_msg(nf90_close(fID),"close",file1_yearly)
163write(*,*) '> Data from "'//file1_yearly//'" downloaded!'
164
165!~~~~~~~~~~~~~~~~~~~~~~~~~ Year 2 - Daily data ~~~~~~~~~~~~~~~~~~~~~~~~~
166! Open the NetCDF file of XIOS outputs
167write(*,*) "Opening "//file2_daily//"..."
168call error_msg(nf90_open(file2_daily,nf90_nowrite,fID),'open',file2_daily)
169
170! Allocate and read the variables
171deallocate(var_read_2d,var_read_3d)
172allocate(var_read_3d(nlon,nlat,nsol),var_read_4d(nlon,nlat,nsoil_PCM,nsol))
173call get_var('ps',var_read_3d)
174do isol = 1,nsoil_PCM
175    call lonlat2vect(nlon,nlat,ngrid,var_read_3d(:,:,isol),ps_ts(:,isol))
176enddo
177call get_var('h2o_layer1',var_read_3d)
178do isol = 1,nsoil_PCM
179    call lonlat2vect(nlon,nlat,ngrid,var_read_3d(:,:,isol),q_h2o_ts(:,isol))
180enddo
181call get_var('co2_layer1',var_read_3d)
182do isol = 1,nsoil_PCM
183    call lonlat2vect(nlon,nlat,ngrid,var_read_3d(:,:,isol),q_co2_ts(:,isol))
184enddo
185if (do_soil) then
186    do islope = 1,nslope
187        if (nslope /= 1) then
188            num='  '
189            write(num,'(i2.2)') islope
190            num = '_slope'//num
191        endif
192        call get_var('soiltemp'//num,var_read_4d)
193        do isol = 1,nsol
194            do isoil = 1,nsoil_PCM
195                call lonlat2vect(nlon,nlat,ngrid,var_read_4d(:,:,isoil,isol),tsoil_ts(:,isoil,islope,isol))
196            enddo
197        enddo
198        call get_var('waterdensity_soil'//num,var_read_4d)
199        do isol = 1,nsol
200            do isoil = 1,nsoil_PCM
201                call lonlat2vect(nlon,nlat,ngrid,var_read_4d(:,:,isoil,isol),watersoil_density_ts(:,isoil,islope,isol))
202            enddo
203        enddo
204    enddo
205endif
206deallocate(var_read_3d,var_read_4d,num)
207
208! Close the NetCDF file of XIOS outputs
209call error_msg(nf90_close(fID),"close",file1_daily)
210write(*,*) '> Data from "'//file2_daily//'" downloaded!'
211
212!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
213! Compute frost from yearly minima
214call compute_frost(ngrid,nslope,h2ofrost_PCM,min_h2ofrost,co2frost_PCM,min_co2frost)
215
216!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
217! Compute ice tendencies from yearly minima
218write(*,*) '> Computing surface ice tendencies'
219call compute_tend(ngrid,nslope,min_h2operice + min_h2ofrost,d_h2oice)
220write(*,*) 'H2O ice tendencies (min/max):', minval(d_h2oice), maxval(d_h2oice)
221call compute_tend(ngrid,nslope,min_co2perice + min_co2frost,d_co2ice)
222write(*,*) 'CO2 ice tendencies (min/max):', minval(d_co2ice), maxval(d_co2ice)
223
224END SUBROUTINE load_xios_data
225!=======================================================================
226
227!=======================================================================
228SUBROUTINE get_timelen(filename,timelen)
229!-----------------------------------------------------------------------
230! NAME
231!     get_timelen
232!
233! DESCRIPTION
234!     Get time dimension length from a NetCDF file.
235!
236! AUTHORS & DATE
237!     JB Clement, 2025
238!
239! NOTES
240!
241!-----------------------------------------------------------------------
242
243! DEPENDENCIES
244! ------------
245use netcdf
246
247! DECLARATION
248! -----------
249implicit none
250
251! ARGUMENTS
252! ---------
253character(*), intent(in)  :: filename ! NetCDF filename
254integer,      intent(out) :: timelen  ! Length of time dimension
255
256! LOCAL VARIABLES
257! ---------------
258integer :: ncid  ! File ID
259integer :: dimid ! Dimension ID
260integer :: ierr  ! Return code
261
262! CODE
263! ----
264! Open the NetCDF file
265ierr = nf90_open(trim(filename),NF90_NOWRITE,ncid)
266if (ierr /= nf90_noerr) then
267    write(*,*) "Error opening file:", trim(nf90_strerror(ierr))
268    error stop
269endif
270
271! Get the dimension ID for 'time_counter'
272ierr = nf90_inq_dimid(ncid,"time_counter",dimid)
273if (ierr /= nf90_noerr) then
274    write(*,*) "Error getting dimid 'time_counter':", trim(nf90_strerror(ierr))
275    error stop
276endif
277
278! Get the size of the dimension 'time_counter'
279ierr = nf90_inquire_dimension(ncid,dimid,len = timelen)
280if (ierr /= nf90_noerr) then
281    write(*,*) "Error getting dimension length:", trim(nf90_strerror(ierr))
282    error stop
283endif
284
285! Close the file
286ierr = nf90_close(ncid)
287if (ierr /= nf90_noerr) then
288    write(*,*) "Error closing file:", trim(nf90_strerror(ierr))
289    error stop
290endif
291
292END SUBROUTINE get_timelen
293!=======================================================================
294
295!=======================================================================
296SUBROUTINE error_msg(ierr,typ,nam)
297!-----------------------------------------------------------------------
298! NAME
299!     error_msg
300!
301! DESCRIPTION
302!     Handle and report NetCDF errors.
303!
304! AUTHORS & DATE
305!     JB Clement, 2025
306!
307! NOTES
308!
309!-----------------------------------------------------------------------
310
311! DECLARATION
312! -----------
313implicit none
314
315! ARGUMENTS
316! ---------
317integer,      intent(in) :: ierr ! NetCDF error code
318character(*), intent(in) :: typ  ! Type of operation (inq, get, put, open, close)
319character(*), intent(in) :: nam  ! Field/file name
320
321! CODE
322! ----
323
324if (ierr == nf90_noerr) return
325select case(typ)
326    case('inq')  ; msg = "Dim/Field <"//trim(nam)//"> is missing"
327    case('get')  ; msg = "Reading failed for <"//trim(nam)//">"
328    case('put')  ; msg = "Writing failed for <"//trim(nam)//">"
329    case('open') ; msg = "File opening failed for <"//trim(nam)//">"
330    case('close'); msg = "File closing failed for <"//trim(nam)//">"
331    case default
332        write(*,*) 'There is no message for this error.'
333        error stop
334end select
335call abort_gcm(__FILE__,trim(msg),ierr)
336
337END SUBROUTINE error_msg
338!=======================================================================
339
340!=======================================================================
341SUBROUTINE get_var_1d(var,v)
342!-----------------------------------------------------------------------
343! NAME
344!     get_var_1d
345!
346! DESCRIPTION
347!     Read a 1D variable from open NetCDF file.
348!
349! AUTHORS & DATE
350!     JB Clement, 2025
351!
352! NOTES
353!
354!-----------------------------------------------------------------------
355
356! DECLARATION
357! -----------
358implicit none
359
360! ARGUMENTS
361! ---------
362character(*),       intent(in)  :: var ! Variable name
363real, dimension(:), intent(out) :: v   ! Output array
364
365! CODE
366! ----
367call error_msg(nf90_inq_varid(fID,var,vID),"inq",var)
368call error_msg(nf90_get_var(fID,vID,v),"get",var)
369
370END SUBROUTINE get_var_1d
371!=======================================================================
372
373!=======================================================================
374SUBROUTINE get_var_2d(var,v)
375!-----------------------------------------------------------------------
376! NAME
377!     get_var_2d
378!
379! DESCRIPTION
380!     Read a 2D variable from open NetCDF file.
381!
382! AUTHORS & DATE
383!     JB Clement, 2025
384!
385! NOTES
386!
387!-----------------------------------------------------------------------
388
389! DECLARATION
390! -----------
391implicit none
392
393! ARGUMENTS
394! ---------
395character(*),         intent(in)  :: var ! Variable name
396real, dimension(:,:), intent(out) :: v   ! Output array
397
398! CODE
399! ----
400
401call error_msg(nf90_inq_varid(fID,var,vID),"inq",var)
402call error_msg(nf90_get_var(fID,vID,v),"get",var)
403
404END SUBROUTINE get_var_2d
405!=======================================================================
406
407!=======================================================================
408SUBROUTINE get_var_3d(var,v)
409!-----------------------------------------------------------------------
410! NAME
411!     get_var_3d
412!
413! DESCRIPTION
414!     Read a 3D variable from open NetCDF file.
415!
416! AUTHORS & DATE
417!     JB Clement, 2025
418!
419! NOTES
420!
421!-----------------------------------------------------------------------
422
423! DECLARATION
424! -----------
425implicit none
426
427! ARGUMENTS
428! ---------
429character(*),           intent(in)  ::  var ! Variable name
430real, dimension(:,:,:), intent(out) ::  v   ! Output array
431
432! CODE
433! ----
434
435call error_msg(nf90_inq_varid(fID,var,vID),"inq",var)
436call error_msg(nf90_get_var(fID,vID,v),"get",var)
437
438END SUBROUTINE get_var_3d
439!=======================================================================
440
441!=======================================================================
442SUBROUTINE get_var_4d(var,v)
443!-----------------------------------------------------------------------
444! NAME
445!     get_var_4d
446!
447! DESCRIPTION
448!     Read a 4D variable from open NetCDF file.
449!
450! AUTHORS & DATE
451!     JB Clement, 2025
452!
453! NOTES
454!
455!-----------------------------------------------------------------------
456
457! DECLARATION
458! -----------
459implicit none
460
461! ARGUMENTS
462! ---------
463character(*),             intent(in)  ::  var ! Variable name
464real, dimension(:,:,:,:), intent(out) ::  v   ! Output array
465
466! CODE
467! ----
468call error_msg(nf90_inq_varid(fID,var,vID),"inq",var)
469call error_msg(nf90_get_var(fID,vID,v),"get",var)
470
471END SUBROUTINE get_var_4d
472!=======================================================================
473
474END MODULE xios_data
Note: See TracBrowser for help on using the repository browser.