source: trunk/LMDZ.COMMON/libf/evolution/output.F90 @ 4152

Last change on this file since 4152 was 4152, checked in by jbclement, 7 days ago

PEM:

  • Fix outputs "diagevo.nc" for 3D data.
  • Fix sign in computing exchanges due to adsorption/ice table and in balancing H2O flux from/into atmosphere.
  • Few cleanings.

JBC
File size: 19.5 KB
Line 
1MODULE output
2!-----------------------------------------------------------------------
3! NAME
4!     output
5!
6! DESCRIPTION
7!     Tools to write PEM diagnostic outputs.
8!
9! AUTHORS & DATE
10!     JB Clement, 12/2025
11!
12! NOTES
13!
14!-----------------------------------------------------------------------
15
16! DEPENDENCIES
17! ------------
18use numerics, only: dp, di, k4, minieps
19
20! DECLARATION
21! -----------
22implicit none
23
24! PARAMETERS
25! ----------
26character(12),                             parameter, private :: diagdef_name = 'diagevo.def'
27integer(di),                               parameter          :: dim_ngrid = -1_di    ! Dimension ID to identify 'ngrid' for "diagevo.nc" logic
28logical(k4),                               protected, private :: is_diagdef = .false. ! Flag to know if "diagevo.def" exists
29logical(k4),                               protected, private :: is_diagevo = .false. ! Flag to know if "diagevo.nc" exists
30character(256), dimension(:), allocatable, protected, private :: outputs_name         ! Names of output variables defined in "diagevo.def"
31integer(di),                               protected, private :: n_outputs            ! Number of output variables defined in "diagevo.def"
32integer(di),                               protected, private :: output_rate          ! Output rate
33integer(di),                               protected, private :: writing_idt = -1_di  ! Current writing number of timestep
34integer(di),                               protected          :: dim_nlon, dim_nlat, dim_nalt, dim_interlayer, dim_nsoil, dim_time ! Dimensions ID to write in "diagevo.nc"
35
36contains
37!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
38
39!=======================================================================
40SUBROUTINE ini_output()
41!-----------------------------------------------------------------------
42! NAME
43!     ini_output
44!
45! DESCRIPTION
46!     Initialize the parameters of module 'output'.
47!
48! AUTHORS & DATE
49!     JB Clement, 12/2025
50!
51! NOTES
52!     At first call.
53!-----------------------------------------------------------------------
54
55! DECLARATION
56! -----------
57implicit none
58
59! CODE
60! ----
61n_outputs = 0
62allocate(outputs_name(0))
63call read_diagdef()
64
65END SUBROUTINE ini_output
66!=======================================================================
67
68!=======================================================================
69SUBROUTINE end_output()
70!-----------------------------------------------------------------------
71! NAME
72!     end_output
73!
74! DESCRIPTION
75!     Deallocate output arrays.
76!
77! AUTHORS & DATE
78!     JB Clement, 01/2026
79!
80! NOTES
81!
82!-----------------------------------------------------------------------
83
84! DECLARATION
85! -----------
86implicit none
87
88! CODE
89! ----
90if (allocated(outputs_name)) deallocate(outputs_name)
91
92END SUBROUTINE end_output
93!=======================================================================
94
95!=======================================================================
96SUBROUTINE set_output_config(output_rate_in)
97!-----------------------------------------------------------------------
98! NAME
99!     set_output_config
100!
101! DESCRIPTION
102!     Setter for 'output' configuration parameters.
103!
104! AUTHORS & DATE
105!     JB Clement, 02/2026
106!
107! NOTES
108!
109!-----------------------------------------------------------------------
110
111! DEPENDENCIES
112! ------------
113use utility,  only: int2str
114use display,  only: print_msg, LVL_NFO
115use stoppage, only: stop_clean
116
117! DECLARATION
118! -----------
119implicit none
120
121! ARGUMENTS
122! ---------
123integer(di), intent(in) :: output_rate_in
124
125! CODE
126! ----
127output_rate = output_rate_in
128call print_msg('output_rate = '//int2str(output_rate),LVL_NFO)
129if (output_rate < 1) call stop_clean(__FILE__,__LINE__,'''output_rate'' must be positive!',1)
130
131END SUBROUTINE set_output_config
132!=======================================================================
133
134!=======================================================================
135SUBROUTINE append_name(var_name)
136!-----------------------------------------------------------------------
137! NAME
138!     append_name
139!
140! DESCRIPTION
141!     Append the variable name read in "diagevo.def" to 'outputs_name'.
142!
143! AUTHORS & DATE
144!     JB Clement, 01/2026
145!
146! NOTES
147!
148!-----------------------------------------------------------------------
149
150! DECLARATION
151! -----------
152implicit none
153
154! ARGUMENTS
155! ---------
156character(*), intent(in) :: var_name
157
158! LOCAL VARIABLES
159! ---------------
160character(256), dimension(:), allocatable :: tmp
161
162! CODE
163! ----
164allocate(tmp(n_outputs + 1))
165if (n_outputs > 0) tmp(1:n_outputs) = outputs_name
166tmp(n_outputs + 1) = var_name
167
168call move_alloc(tmp,outputs_name)
169n_outputs = n_outputs + 1
170
171END SUBROUTINE append_name
172!=======================================================================
173
174!=======================================================================
175SUBROUTINE read_diagdef()
176!-----------------------------------------------------------------------
177! NAME
178!     read_diagdef
179!
180! DESCRIPTION
181!     Read the variables names defined in "diagevo.def" to ouput.
182!
183! AUTHORS & DATE
184!     JB Clement, 01/2026
185!
186! NOTES
187!     > If no "diagevo.def", then no output.
188!     > If "diagevo.def" is empty, then all the possible coded outputs.
189!     > If "diagevo.def" holds a list of variable names, then these are
190!     outputted.
191!-----------------------------------------------------------------------
192
193! DEPENDENCIES
194! ------------
195use stoppage, only: stop_clean
196use utility,  only: is_id_1st_char, is_id_char, int2str
197use display,  only: print_msg, LVL_NFO, LVL_WRN
198
199! DECLARATION
200! -----------
201implicit none
202
203! LOCAL VARIABLES
204! ---------------
205integer(di)    :: i, iline, ierr, funit
206character(256) :: line
207
208! CODE
209! ----
210inquire(file = diagdef_name,exist = is_diagdef)
211if (is_diagdef) then
212    call print_msg('> Reading "'//diagdef_name//'"',LVL_NFO)
213    open(newunit = funit,file = diagdef_name,status = 'old',form = 'formatted',action = 'read',iostat = ierr)
214    if (ierr /= 0) call stop_clean(__FILE__,__LINE__,'error opening file "'//diagdef_name//'"!',ierr)
215    iline = 1
216    do
217        ! Read the line
218        read(funit,'(a)',iostat = ierr) line
219        if (ierr /= 0) exit ! Empty file or end of file
220        iline = iline + 1
221
222        ! Examine the line
223        line = trim(line)
224        if (len(line) == 0) cycle ! Ignore empty lines
225        if (.not. is_id_1st_char(line(1:1))) call print_msg('Invalid output name at line '//int2str(iline)//' in '//diagdef_name//'"!',LVL_WRN)
226        do i = 2,len(line)
227            if (.not. is_id_char(line(i:i))) call print_msg('Invalid output name at line '//int2str(iline)//' in '//diagdef_name//'"!',LVL_WRN)
228        end do
229
230        ! Fill in the table 'outputs_name'
231        call append_name(line)
232    end do
233    close(funit)
234    if (n_outputs == 0) call print_msg('There is a "'//diagdef_name//'" which is is empty. So all the possible coded variables are outputted.',LVL_WRN)
235    ! Create the output file
236    call create_diagevo()
237else
238    call print_msg('There is no "'//diagdef_name//'". So no variable is outputted.',LVL_WRN)
239end if
240
241END SUBROUTINE read_diagdef
242!=======================================================================
243
244!=======================================================================
245FUNCTION is_output_selected(var_name) RESULT(output_var)
246!-----------------------------------------------------------------------
247! NAME
248!     is_output_selected
249!
250! DESCRIPTION
251!     Read the variables names to output defined in "diagevo.def".
252!
253! AUTHORS & DATE
254!     JB Clement, 01/2026
255!
256! NOTES
257!     > If no "diagevo.def", then no output.
258!     > If "diagevo.def" is empty, then all the possible coded outputs.
259!     > If "diagevo.def" holds a list of variable names, then these are
260!     outputted.
261!-----------------------------------------------------------------------
262
263! DECLARATION
264! -----------
265implicit none
266
267! ARGUMENTS
268! ---------
269character(*), intent(in) :: var_name
270logical(k4)              :: output_var
271
272! LOCAL VARIABLES
273! ---------------
274integer(di) :: i
275
276! CODE
277! ----
278output_var = .false.
279
280! If diagdef exists but is empty, output everything
281if (n_outputs == 0) then
282    output_var = .true.
283    return
284end if
285
286do i = 1,n_outputs
287    if (trim(outputs_name(i)) == trim(var_name)) then
288        output_var = .true.
289        return
290    end if
291end do
292
293END FUNCTION is_output_selected
294!=======================================================================
295
296!=======================================================================
297SUBROUTINE create_diagevo()
298!-----------------------------------------------------------------------
299! NAME
300!     create_diagevo
301!
302! DESCRIPTION
303!     Create a NetCDF file "diagevo.nc" containing the diagnostic
304!     variables.
305!
306! AUTHORS & DATE
307!     JB Clement, 01/2026
308!
309! NOTES
310!
311!-----------------------------------------------------------------------
312
313! DEPENDENCIES
314! ------------
315use netcdf,     only: nf90_double, nf90_clobber, nf90_unlimited, nf90_global, &
316                      nf90_create, nf90_enddef, nf90_def_dim, nf90_def_var, nf90_put_att
317use io_netcdf,  only: check_nc, open_nc, close_nc, put_var_nc, diagevo_name, ncid_diagevo
318use geometry,   only: dim_init, ngrid, nlon, nlat, nlayer, nsoil, cell_area, longitudes, latitudes, vect2dyngrd
319use stoppage,   only: stop_clean
320use soil,       only: mlayer
321use atmosphere, only: ap, bp
322use display,    only: print_msg, LVL_NFO
323
324! DECLARATION
325! -----------
326implicit none
327
328! LOCAL VARIABLES
329! ---------------
330integer(di)                           :: varid ! Variable ID
331integer(di)                           :: nlon_loc
332real(dp), dimension(:,:), allocatable :: var_dyn
333
334! CODE
335! ----
336! Check if dimensions are well initialized
337if (.not. dim_init) call stop_clean(__FILE__,__LINE__,'dimensions are not initilized in the PEM!',1)
338if (ngrid == 1) then ! 1D
339    nlon_loc = 1
340else ! 3D
341    nlon_loc = nlon + 1
342end if
343
344! Create file
345call print_msg('> Creating "'//diagevo_name//'"',LVL_NFO)
346call check_nc(nf90_create(diagevo_name,nf90_clobber,ncid_diagevo),'creating '//diagevo_name) ! Enter define mode
347
348! Define dimensions
349call check_nc(nf90_def_dim(ncid_diagevo,'Time',NF90_UNLIMITED,dim_time),'defining dimension Time')
350call check_nc(nf90_def_dim(ncid_diagevo,'longitude',nlon_loc,dim_nlon),'defining dimension longitude')
351call check_nc(nf90_def_dim(ncid_diagevo,'latitude',nlat,dim_nlat),'defining dimension latitude')
352call check_nc(nf90_def_dim(ncid_diagevo,'altitude',nlayer,dim_nalt),'defining dimension altitude')
353call check_nc(nf90_def_dim(ncid_diagevo,'interlayer',nlayer + 1,dim_interlayer),'defining dimension interlayer')
354call check_nc(nf90_def_dim(ncid_diagevo,'subsurface_layers',nsoil,dim_nsoil),'defining dimension subsurface_layers')
355
356! Define variables
357call check_nc(nf90_def_var(ncid_diagevo,'Time',nf90_double,(/dim_time/),varid),'defining variable Time')
358call check_nc(nf90_put_att(ncid_diagevo,varid,'title','Year of run'),'putting title attribute for Time')
359call check_nc(nf90_put_att(ncid_diagevo,varid,'units','Planetary year'),'putting units attribute for Time')
360
361call check_nc(nf90_def_var(ncid_diagevo,'longitude',nf90_double,(/dim_nlon/),varid),'defining variable longitude')
362call check_nc(nf90_put_att(ncid_diagevo,varid,'title','Longitudes of the grid'),'putting title attribute for longitude')
363
364call check_nc(nf90_def_var(ncid_diagevo,'latitude',nf90_double,(/dim_nlat/),varid),'defining variable latitude')
365call check_nc(nf90_put_att(ncid_diagevo,varid,'title','Latitudes of the grid'),'putting title attribute for latitude')
366
367call check_nc(nf90_def_var(ncid_diagevo,'ap',nf90_double,(/dim_interlayer/),varid),'defining variable ap')
368call check_nc(nf90_put_att(ncid_diagevo,varid,'title','Hybrid pressure at interlayers'),'putting title attribute for ap')
369
370call check_nc(nf90_def_var(ncid_diagevo,'bp',nf90_double,(/dim_interlayer/),varid),'defining variable bp')
371call check_nc(nf90_put_att(ncid_diagevo,varid,'title','Hybrid sigma at interlayers'),'putting title attribute for bp')
372
373call check_nc(nf90_def_var(ncid_diagevo,'soildepth',nf90_double,(/dim_nsoil/),varid),'defining variable soildepth')
374call check_nc(nf90_put_att(ncid_diagevo,varid,'title','Depths of soil layers'),'putting title attribute for soildepth')
375
376call check_nc(nf90_def_var(ncid_diagevo,'cell_area',nf90_double,(/dim_nlon,dim_nlat/),varid),'defining variable cell_area')
377call check_nc(nf90_put_att(ncid_diagevo,varid,'title','Cell area'),'putting title attribute for cell_area')
378
379! Global attributes
380call check_nc(nf90_put_att(ncid_diagevo,nf90_global,'title','PEM diagnostic file'),'putting global attribute')
381
382! Leave define mode and putting variables defining dimensions
383call check_nc(nf90_enddef(ncid_diagevo),'leaving define mode')
384allocate(var_dyn(nlon_loc,nlat))
385call vect2dyngrd(longitudes,var_dyn)
386call put_var_nc('longitude',var_dyn(:,1))
387call vect2dyngrd(latitudes,var_dyn)
388call put_var_nc('latitude',var_dyn(1,:))
389call put_var_nc('ap',ap)
390call put_var_nc('bp',bp)
391call put_var_nc('soildepth',mlayer)
392call vect2dyngrd(cell_area,var_dyn,extensive = .true.)
393call put_var_nc('cell_area',var_dyn)
394deallocate(var_dyn)
395call close_nc(diagevo_name)
396
397! File creation done
398is_diagevo = .true.
399
400END SUBROUTINE create_diagevo
401!=======================================================================
402
403!=======================================================================
404SUBROUTINE write_diagevo(var_name,title,units,var,dimids,extensive)
405!-----------------------------------------------------------------------
406! NAME
407!     write_diagevo
408!
409! DESCRIPTION
410!     Write selected diagnostic variables to NetCDF file "diagevo.nc".
411!
412! AUTHORS & DATE
413!     JB Clement, 01/2026
414!
415! NOTES
416!     > Output rate controlled by 'output_rate' from "run.def" (PEM).
417!     > Variable selection controlled by the file "diagevo.def".
418!         - If no "diagevo.def", then no output.
419!         - If "diagevo.def" is empty, then all the possible coded
420!           outputs are done.
421!         - If "diagevo.def" holds a list of variable names, then these
422!           are outputted.
423!     > If 'var' is an array, then the argument 'dimids' is mandatory.
424!       'dimids' contains the dimension IDs of the variable in order.
425!       If the variable holds a dimension 'ngrid', then it must appear
426!       in the first position.
427!     > For extensive variables, the values at the poles must be
428!       averaged over the number of longitudes to be computed/outputted
429!       on the dynamical grid.
430!-----------------------------------------------------------------------
431
432! DEPENDENCIES
433! ------------
434use io_netcdf, only: open_nc, close_nc, def_var_nc, put_var_nc, diagevo_name
435use evolution, only: idt, n_yr_run
436use stoppage,  only: stop_clean
437use geometry,  only: vect2dyngrd, ngrid, nlon, nlat
438
439! DECLARATION
440! -----------
441implicit none
442
443! ARGUMENTS
444! ---------
445character(*),               intent(in)           :: var_name, title, units
446real(dp),    dimension(..), intent(in)           :: var ! Assumed-rank array
447integer(di), dimension(:),  intent(in), optional :: dimids
448logical(k4),                intent(in), optional :: extensive
449
450! LOCAL VARIABLES
451! ---------------
452integer(di), dimension(:),       allocatable :: dimids_nc
453integer(di)                                  :: idim_ngrid, nlon_loc, ndim, i, j
454real(dp),    dimension(:,:),     allocatable :: var_dyn
455real(dp),    dimension(:,:,:),   allocatable :: var1d_dyn
456real(dp),    dimension(:,:,:,:), allocatable :: var2d_dyn
457integer(di)                                  :: itime ! Current time index to record variables
458logical(k4)                                  :: extensive_l
459
460
461! CODE
462! ----
463! If no diagdef file, no output at all
464if (.not. is_diagdef) return
465
466! If variable not selected, return
467! Empty diagdef file means all the variables are ouputted
468if (.not. is_output_selected(var_name)) return
469
470! If output timing not met, return
471if (abs(mod(idt,output_rate)) > minieps) return
472
473! If it is the first writing of the current timestep
474if (writing_idt /= idt) then
475    writing_idt = idt
476    call open_nc(diagevo_name,'write',itime) ! Get the current time index
477    call put_var_nc('Time',n_yr_run,itime)
478end if
479
480! Write the variable
481extensive_l = .false.
482if (present(extensive)) extensive_l = extensive
483call open_nc(diagevo_name,'write') ! Open diagevo_name only it is not done yet
484select rank (var) ! Assumed-rank array removes compile-time rank information which is required by the following subroutines using generic interface overload
485    rank(0) ! Scalar
486        call def_var_nc(var_name,title,units,(/dim_time/)) ! Define the variable only it is not done yet
487        call put_var_nc(var_name,var,itime)
488    rank default ! Arrays
489        if (.not. present(dimids)) call stop_clean(__FILE__,__LINE__,'The argument ''dimids'' must be present to output arrays!',1)
490        ndim = size(dimids)
491        if (any(dimids == dim_ngrid)) then
492            idim_ngrid = findloc(dimids,dim_ngrid,dim = 1)
493            if (idim_ngrid /= 1) call stop_clean(__FILE__,__LINE__,'The dimension ''ngrid'' must be the first one to output the array!',1)
494            allocate(dimids_nc(ndim + 2))
495            dimids_nc(1) = dim_nlon
496            dimids_nc(2) = dim_nlat
497            dimids_nc(3:ndim + 1) = dimids(2:ndim)
498            dimids_nc(ndim + 2) = dim_time
499            if (ngrid == 1) then ! 1D
500                nlon_loc = 1
501            else ! 3D
502                nlon_loc = nlon + 1
503            end if
504            call def_var_nc(var_name,title,units,dimids_nc) ! Define the variable only it is not done yet
505            select rank (var) ! Assumed-rank array removes compile-time rank information which is required by the following subroutines using generic interface overload
506                rank(1)
507                    allocate(var_dyn(nlon_loc,nlat))
508                    call vect2dyngrd(var,var_dyn,extensive_l)
509                    call put_var_nc(var_name,var_dyn,itime)
510                    deallocate(var_dyn)
511                rank(2)
512                    allocate(var1d_dyn(nlon_loc,nlat,size(var,2)))
513                    do i = 1,size(var,2)
514                        call vect2dyngrd(var(:,i),var1d_dyn(:,:,i),extensive_l)
515                    end do
516                    call put_var_nc(var_name,var1d_dyn,itime)
517                    deallocate(var1d_dyn)
518                rank(3)
519                    allocate(var2d_dyn(nlon_loc,nlat,size(var,2),size(var,3)))
520                    do i = 1,size(var,2)
521                        do j = 1,size(var,3)
522                            call vect2dyngrd(var(:,i,j),var2d_dyn(:,:,i,j),extensive_l)
523                        end do
524                    end do
525                    call put_var_nc(var_name,var2d_dyn,itime)
526                    deallocate(var2d_dyn)
527                rank default
528                    call stop_clean(__FILE__,__LINE__,'unsupported rank for the input array!',1)
529            end select
530        else ! No 'ngrid' dimension
531            allocate(dimids_nc(ndim + 1))
532            dimids_nc(1:ndim) = dimids(:)
533            dimids_nc(ndim + 1) = dim_time
534            call def_var_nc(var_name,title,units,dimids_nc) ! Define the variable only it is not done yet
535            select rank (var) ! Assumed-rank array removes compile-time rank information which is required by the following subroutines using generic interface overload
536                rank(1)
537                    call put_var_nc(var_name,var_dyn,itime)
538                rank(2)
539                    call put_var_nc(var_name,var1d_dyn,itime)
540                rank(3)
541                    call put_var_nc(var_name,var2d_dyn,itime)
542                rank(4)
543                    call put_var_nc(var_name,var2d_dyn,itime)
544                rank default
545                    call stop_clean(__FILE__,__LINE__,'unsupported rank for the input array!',1)
546            end select
547        end if
548        deallocate(dimids_nc)
549end select
550call close_nc(diagevo_name)
551
552END SUBROUTINE write_diagevo
553!=================================================================
554
555END MODULE output
Note: See TracBrowser for help on using the repository browser.