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

Last change on this file since 4147 was 4147, checked in by jbclement, 8 days ago

PEM:

  • Simplification of subroutines to convert data between the physical and the dynamical/lon-lat grids + making them more robust.
  • Correction for air mass to give back to the PCM. The variable is extensive so poles must be treated specifically.
  • Making the PEM able to do 0 year.
  • Explicit information about the frost values computed by the PEM + enforcing positivity of yearly minima.

JBC
File size: 19.1 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)
387call vect2dyngrd(latitudes,var_dyn)
388call put_var_nc('latitude',var_dyn)
389call put_var_nc('ap',ap)
390call put_var_nc('bp',bp)
391call put_var_nc('soildepth',mlayer)
392call vect2dyngrd(cell_area,var_dyn)
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_in)
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_in' is mandatory.
424!       'dimids_in' 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!-----------------------------------------------------------------------
428
429! DEPENDENCIES
430! ------------
431use io_netcdf, only: open_nc, close_nc, def_var_nc, put_var_nc, diagevo_name
432use evolution, only: idt, n_yr_run
433use stoppage,  only: stop_clean
434use geometry,  only: vect2dyngrd, ngrid, nlon, nlat
435
436! DECLARATION
437! -----------
438implicit none
439
440! ARGUMENTS
441! ---------
442character(*),               intent(in)           :: var_name, title, units
443real(dp),    dimension(..), intent(in)           :: var ! Assumed-rank array
444integer(di), dimension(:),  intent(in), optional :: dimids_in
445
446! LOCAL VARIABLES
447! ---------------
448integer(di), dimension(:),       allocatable :: dimids
449integer(di)                                  :: idim_ngrid, nlon_loc, ndim, i, j
450real(dp),    dimension(:,:),     allocatable :: var_dyn
451real(dp),    dimension(:,:,:),   allocatable :: var1d_dyn
452real(dp),    dimension(:,:,:,:), allocatable :: var2d_dyn
453integer(di)                                  :: itime ! Current time index to record variables
454
455
456! CODE
457! ----
458! If no diagdef file, no output at all
459if (.not. is_diagdef) return
460
461! If variable not selected, return
462! Empty diagdef file means all the variables are ouputted
463if (.not. is_output_selected(var_name)) return
464
465! If output timing not met, return
466if (abs(mod(idt,output_rate)) > minieps) return
467
468! If it is the first writing of the current timestep
469if (writing_idt /= idt) then
470    writing_idt = idt
471    call open_nc(diagevo_name,'write',itime) ! Get the current time index
472    call put_var_nc('Time',n_yr_run,itime)
473end if
474
475! Write the variable
476call open_nc(diagevo_name,'write') ! Open diagevo_name only it is not done yet
477select rank (var) ! Assumed-rank array removes compile-time rank information which is required by the following subroutines using generic interface overload
478    rank(0) ! Scalar
479        call def_var_nc(var_name,title,units,(/dim_time/)) ! Define the variable only it is not done yet
480        call put_var_nc(var_name,var,itime)
481    rank default ! Arrays
482        if (.not. present(dimids_in)) call stop_clean(__FILE__,__LINE__,'The argument ''dimids_in'' must be present to output arrays!',1)
483        ndim = size(dimids_in)
484        if (any(dimids_in == dim_ngrid)) then
485            idim_ngrid = findloc(dimids_in,dim_ngrid,dim = 1)
486            if (idim_ngrid /= 1) call stop_clean(__FILE__,__LINE__,'The dimension ''ngrid'' must be the first one to output the array!',1)
487            allocate(dimids(ndim + 2))
488            dimids(1) = dim_nlon
489            dimids(2) = dim_nlat
490            dimids(3:ndim + 1) = dimids_in(2:ndim)
491            dimids(ndim + 2) = dim_time
492            if (ngrid == 1) then ! 1D
493                nlon_loc = 1
494            else ! 3D
495                nlon_loc = nlon + 1
496            end if
497            call def_var_nc(var_name,title,units,dimids) ! Define the variable only it is not done yet
498            select rank (var) ! Assumed-rank array removes compile-time rank information which is required by the following subroutines using generic interface overload
499                rank(1)
500                    allocate(var_dyn(nlon_loc,nlat))
501                    call vect2dyngrd(var,var_dyn)
502                    call put_var_nc(var_name,var_dyn,itime)
503                    deallocate(var_dyn)
504                rank(2)
505                    allocate(var1d_dyn(nlon_loc,nlat,size(var,2)))
506                    do i = 1,size(var,2)
507                        call vect2dyngrd(var(:,i),var1d_dyn(:,:,i))
508                    end do
509                    call put_var_nc(var_name,var1d_dyn,itime)
510                    deallocate(var1d_dyn)
511                rank(3)
512                    allocate(var2d_dyn(nlon_loc,nlat,size(var,2),size(var,3)))
513                    do i = 1,size(var,2)
514                        do j = 1,size(var,3)
515                            call vect2dyngrd(var(:,i,j),var2d_dyn(:,:,i,j))
516                        end do
517                    end do
518                    call put_var_nc(var_name,var2d_dyn,itime)
519                    deallocate(var2d_dyn)
520                rank default
521                    call stop_clean(__FILE__,__LINE__,'unsupported rank for the input array!',1)
522            end select
523        else ! No 'ngrid' dimension
524            allocate(dimids(ndim + 1))
525            dimids(1:ndim) = dimids_in(:)
526            dimids(ndim + 1) = dim_time
527            call def_var_nc(var_name,title,units,dimids) ! Define the variable only it is not done yet
528            select rank (var) ! Assumed-rank array removes compile-time rank information which is required by the following subroutines using generic interface overload
529                rank(1)
530                    call put_var_nc(var_name,var_dyn,itime)
531                rank(2)
532                    call put_var_nc(var_name,var1d_dyn,itime)
533                rank(3)
534                    call put_var_nc(var_name,var2d_dyn,itime)
535                rank(4)
536                    call put_var_nc(var_name,var2d_dyn,itime)
537                rank default
538                    call stop_clean(__FILE__,__LINE__,'unsupported rank for the input array!',1)
539            end select
540        end if
541        deallocate(dimids)
542end select
543call close_nc(diagevo_name)
544
545END SUBROUTINE write_diagevo
546!=================================================================
547
548END MODULE output
Note: See TracBrowser for help on using the repository browser.