source: IOIPSL/trunk/src/histcom.f90 @ 714

Last change on this file since 714 was 714, checked in by bellier, 15 years ago

Declare the time axis in double precision in NETCDF files

  • Property svn:keywords set to Id
File size: 79.4 KB
Line 
1MODULE histcom
2!-
3!$Id$
4!-
5! This software is governed by the CeCILL license
6! See IOIPSL/IOIPSL_License_CeCILL.txt
7!-
8  USE netcdf
9!-
10  USE stringop, ONLY : nocomma,cmpblank,findpos,find_str,strlowercase
11  USE mathelp,  ONLY : mathop,moycum,buildop
12  USE fliocom,  ONLY : flio_dom_file,flio_dom_att
13  USE calendar
14  USE errioipsl, ONLY : ipslerr,ipsldbg
15!-
16  IMPLICIT NONE
17!-
18  PRIVATE
19  PUBLIC :: histbeg,histdef,histhori,histvert,histend, &
20 &          histwrite,histclo,histsync,ioconf_modname
21!---------------------------------------------------------------------
22!- Some confusing vocabulary in this code !
23!- =========================================
24!-
25!- A REGULAR grid is a grid which is i,j indexes
26!- and thus it is stored in a 2D matrix.
27!- This is opposed to a IRREGULAR grid which is only in a vector
28!- and where we do not know which neighbors we have.
29!- As a consequence we need the bounds for each grid-cell.
30!-
31!- A RECTILINEAR grid is a special case of a regular grid
32!- in which all longitudes for i constant are equal
33!- and all latitudes for j constant.
34!- In other words we do not need the full 2D matrix
35!- to describe the grid, just two vectors.
36!---------------------------------------------------------------------
37  INTERFACE histwrite
38!---------------------------------------------------------------------
39!- The "histwrite" routines will give the data to the I/O system.
40!- It will trigger the operations to be performed,
41!- and the writting to the file if needed
42!-
43!- We test for the work to be done at this time here so that at a
44!- later stage we can call different operation and write subroutine
45!- for the REAL and INTEGER interfaces
46!-
47!- INPUT
48!- pfileid  : The ID of the file on which this variable is to be,
49!-            written. The variable should have been defined in
50!-            this file before.
51!- pvarname : The short name of the variable
52!- pitau    : Current timestep
53!- pdata    : The variable, I mean the real data !
54!- nbindex  : The number of indexes provided. If it is equal to
55!-            the size of the full field as provided in histdef
56!-            then nothing is done.
57!- nindex   : The indices used to expand the variable (pdata)
58!-            onto the full field.
59!---------------------------------------------------------------------
60!- histwrite - we have to prepare different type of fields :
61!-             real and integer, 1,2 or 3D
62    MODULE PROCEDURE histwrite_r1d,histwrite_r2d,histwrite_r3d
63  END INTERFACE
64!-
65  INTERFACE histbeg
66    MODULE PROCEDURE histbeg_totreg,histbeg_regular,histbeg_irregular
67  END INTERFACE
68!-
69  INTERFACE histhori
70    MODULE PROCEDURE histhori_regular,histhori_irregular
71  END INTERFACE
72!-
73! Fixed parameter
74!-
75  INTEGER,PARAMETER :: nb_files_max=20,nb_var_max=400, &
76 &                     nb_hax_max=5,nb_zax_max=10,nbopp_max=10
77  REAL,PARAMETER :: missing_val=nf90_fill_real
78!-
79  INTEGER :: bufftmp_max(nb_files_max) = 1
80!-
81! Time variables
82!-
83  INTEGER,SAVE :: itau0(nb_files_max)=0
84  REAL,DIMENSION(nb_files_max),SAVE ::date0,deltat
85!-
86! Counter of elements
87!-
88  INTEGER,SAVE :: nb_files=0
89  INTEGER,DIMENSION(nb_files_max),SAVE :: nb_var=0,nb_tax=0
90!-
91! DOMAIN IDs for files
92!-
93  INTEGER,DIMENSION(nb_files_max),SAVE :: dom_id_svg=-1
94!-
95! NETCDF IDs for files and axes
96!-
97  INTEGER,DIMENSION(nb_files_max),SAVE :: ncdf_ids,xid,yid,tid
98!-
99! General definitions in the NETCDF file
100!-
101  INTEGER,DIMENSION(nb_files_max,2),SAVE :: &
102 &   full_size=0,slab_ori,slab_sz
103!-
104! The horizontal axes
105!-
106  INTEGER,SAVE :: nb_hax(nb_files_max)=0
107  CHARACTER(LEN=25),SAVE :: hax_name(nb_files_max,nb_hax_max,2)
108!-
109! The vertical axes
110!-
111  INTEGER,SAVE :: nb_zax(nb_files_max)=0
112  INTEGER,DIMENSION(nb_files_max,nb_zax_max),SAVE :: &
113  &  zax_size,zax_ids
114  CHARACTER(LEN=20),SAVE :: zax_name(nb_files_max,nb_zax_max)
115!-
116! Informations on each variable
117!-
118  INTEGER,DIMENSION(nb_files_max,nb_var_max),SAVE :: &
119 &  nbopp
120  CHARACTER(LEN=20),DIMENSION(nb_files_max,nb_var_max),SAVE :: &
121 &  name,unit_name
122  CHARACTER(LEN=80),DIMENSION(nb_files_max,nb_var_max),SAVE :: &
123 &  title,fullop
124  CHARACTER(LEN=7),SAVE :: topp(nb_files_max,nb_var_max)
125  CHARACTER(LEN=7),SAVE :: sopps(nb_files_max,nb_var_max,nbopp_max)
126  REAL,SAVE :: scal(nb_files_max,nb_var_max,nbopp_max)
127!- Sizes of the associated grid and zommed area
128  INTEGER,DIMENSION(nb_files_max,nb_var_max,3),SAVE :: &
129 &   scsize,zorig,zsize
130!- Sizes for the data as it goes through the various math operations
131  INTEGER,SAVE :: datasz_in(nb_files_max,nb_var_max,3) = -1
132  INTEGER,SAVE :: datasz_max(nb_files_max,nb_var_max) = -1
133!-
134  INTEGER,DIMENSION(nb_files_max,nb_var_max),SAVE :: &
135 &  var_haxid,var_zaxid,var_axid,ncvar_ids
136!-
137  REAL,DIMENSION(nb_files_max,nb_var_max,2),SAVE :: hist_minmax
138  LOGICAL,DIMENSION(nb_files_max,nb_var_max),SAVE :: &
139 &  hist_calc_rng=.FALSE.,hist_wrt_rng=.FALSE.
140!-
141  REAL,DIMENSION(nb_files_max,nb_var_max),SAVE :: &
142 &  freq_opp,freq_wrt
143  INTEGER,DIMENSION(nb_files_max,nb_var_max),SAVE :: &
144 &  last_opp,last_wrt,last_opp_chk,last_wrt_chk,nb_opp,nb_wrt,point
145!-
146! Book keeping for the buffers
147!-
148  INTEGER,SAVE :: buff_pos=0
149  REAL,ALLOCATABLE,DIMENSION(:),SAVE :: buffer
150  LOGICAL,DIMENSION(nb_files_max),SAVE :: zoom=.FALSE.,regular=.TRUE.
151!-
152! Book keeping of the axes
153!-
154  INTEGER,DIMENSION(nb_files_max,nb_var_max),SAVE :: &
155 &  tdimid,tax_last
156  CHARACTER(LEN=40),DIMENSION(nb_files_max,nb_var_max),SAVE :: &
157 &  tax_name
158!-
159! A list of functions which require special action
160! (Needs to be updated when functions are added
161!  but they are well located here)
162!-
163  CHARACTER(LEN=120),SAVE :: &
164 &  indchfun = 'scatter, fill, gather, coll', &
165 &  fuchnbout = 'scatter, fill'
166!- Some configurable variables with locks
167  CHARACTER(LEN=80),SAVE :: model_name='An IPSL model'
168  LOGICAL,SAVE :: lock_modname=.FALSE.
169!-
170!===
171CONTAINS
172!===
173!-
174SUBROUTINE histbeg_totreg                 &
175 & (pfilename,pim,plon,pjm,plat,          &
176 &  par_orix,par_szx,par_oriy,par_szy,    &
177 &  pitau0,pdate0,pdeltat,phoriid,pfileid,domain_id)
178!---------------------------------------------------------------------
179!- This is just an interface for histbeg_regular in case when
180!- the user provides plon and plat as vectors.
181!- Obviously this can only be used for very regular grids.
182!-
183!- INPUT
184!-
185!- pfilename : Name of the netcdf file to be created
186!- pim       : Size of arrays in longitude direction
187!- plon      : Coordinates of points in longitude
188!- pjm       : Size of arrays in latitude direction
189!- plat      : Coordinates of points in latitude
190!-
191!- The next 4 arguments allow to define a horizontal zoom
192!- for this file. It is assumed that all variables to come
193!- have the same index space. This can not be assumed for
194!- the z axis and thus we define the zoom in histdef.
195!-
196!- par_orix  : Origin of the slab of data within the X axis (pim)
197!- par_szx   : Size of the slab of data in X
198!- par_oriy  : Origin of the slab of data within the Y axis (pjm)
199!- par_szy   : Size of the slab of data in Y
200!-
201!- pitau0    : time step at which the history tape starts
202!- pdate0    : The Julian date at which the itau was equal to 0
203!- pdeltat   : Time step in seconds. Time step of the counter itau
204!-             used in histwrite for instance
205!-
206!- OUTPUT
207!-
208!- phoriid   : ID of the horizontal grid
209!- pfileid   : ID of the netcdf file
210!-
211!- Optional INPUT arguments
212!-
213!- domain_id       : Domain identifier
214!-
215!- TO DO
216!-
217!- This package should be written in f90
218!- and use the following features :
219!-    - structures for the meta-data of the files and variables
220!-    - memory allocation as needed
221!-    - Pointers
222!-
223!- VERSION
224!-
225!---------------------------------------------------------------------
226  IMPLICIT NONE
227!-
228  CHARACTER(LEN=*) :: pfilename
229  INTEGER,INTENT(IN) :: pim,pjm
230  REAL,DIMENSION(pim),INTENT(IN) :: plon
231  REAL,DIMENSION(pjm),INTENT(IN) :: plat
232  INTEGER,INTENT(IN):: par_orix,par_szx,par_oriy,par_szy
233  INTEGER,INTENT(IN) :: pitau0
234  REAL,INTENT(IN) :: pdate0,pdeltat
235  INTEGER,INTENT(OUT) :: pfileid,phoriid
236  INTEGER,INTENT(IN),OPTIONAL :: domain_id
237!-
238  REAL,ALLOCATABLE,DIMENSION(:,:) :: lon_tmp,lat_tmp
239  LOGICAL :: l_dbg
240!---------------------------------------------------------------------
241  CALL ipsldbg (old_status=l_dbg)
242!-
243  IF (l_dbg) WRITE(*,*) "histbeg_totreg"
244!-
245  ALLOCATE(lon_tmp(pim,pjm),lat_tmp(pim,pjm))
246!-
247  lon_tmp(:,:) = SPREAD(plon(:),2,pjm)
248  lat_tmp(:,:) = SPREAD(plat(:),1,pim)
249!-
250  CALL histbeg_regular &
251 &  (pfilename,pim,lon_tmp,pjm,lat_tmp, &
252 &   par_orix,par_szx,par_oriy,par_szy, &
253 &   pitau0,pdate0,pdeltat,phoriid,pfileid, &
254 &   .TRUE.,domain_id)
255!-
256  DEALLOCATE(lon_tmp,lat_tmp)
257!----------------------------
258END SUBROUTINE histbeg_totreg
259!===
260SUBROUTINE histbeg_regular &
261 & (pfilename,pim,plon,pjm,plat,           &
262 &  par_orix,par_szx,par_oriy,par_szy,     &
263 &  pitau0,pdate0,pdeltat,phoriid,pfileid, &
264 &  opt_rectilinear,domain_id)
265!---------------------------------------------------------------------
266!- This subroutine initializes a netcdf file and returns the ID.
267!- It will set up the geographical space on which the data will be
268!- stored and offers the possibility of seting a zoom.
269!- It also gets the global parameters into the I/O subsystem.
270!-
271!- INPUT
272!-
273!- pfilename : Name of the netcdf file to be created
274!- pim       : Size of arrays in longitude direction
275!- plon      : Coordinates of points in longitude
276!- pjm       : Size of arrays in latitude direction
277!- plat      : Coordinates of points in latitude
278!-
279!- The next 4 arguments allow to define a horizontal zoom
280!- for this file. It is assumed that all variables to come
281!- have the same index space. This can not be assumed for
282!- the z axis and thus we define the zoom in histdef.
283!-
284!- par_orix  : Origin of the slab of data within the X axis (pim)
285!- par_szx   : Size of the slab of data in X
286!- par_oriy  : Origin of the slab of data within the Y axis (pjm)
287!- par_szy   : Size of the slab of data in Y
288!-
289!- pitau0    : time step at which the history tape starts
290!- pdate0    : The Julian date at which the itau was equal to 0
291!- pdeltat   : Time step in seconds. Time step of the counter itau
292!-             used in histwrite for instance
293!-
294!- OUTPUT
295!-
296!- phoriid   : ID of the horizontal grid
297!- pfileid   : ID of the netcdf file
298!-
299!- Optional INPUT arguments
300!-
301!- opt_rectilinear : If true we know the grid is rectilinear
302!- domain_id       : Domain identifier
303!-
304!- TO DO
305!-
306!- This package should be written in F90 and use the following
307!- feature :
308!-   - structures for the meta-data of the files and variables
309!-   - memory allocation as needed
310!-   - Pointers
311!-
312!- VERSION
313!-
314!---------------------------------------------------------------------
315  IMPLICIT NONE
316!-
317  CHARACTER(LEN=*) :: pfilename
318  INTEGER,INTENT(IN) :: pim,pjm
319  REAL,DIMENSION(pim,pjm),INTENT(IN) :: plon,plat
320  INTEGER,INTENT(IN):: par_orix,par_szx,par_oriy,par_szy
321  INTEGER,INTENT(IN) :: pitau0
322  REAL,INTENT(IN) :: pdate0,pdeltat
323  INTEGER,INTENT(OUT) :: pfileid,phoriid
324  LOGICAL,INTENT(IN),OPTIONAL :: opt_rectilinear
325  INTEGER,INTENT(IN),OPTIONAL :: domain_id
326!-
327  INTEGER :: ncid,iret
328  CHARACTER(LEN=120) :: file
329  CHARACTER(LEN=30) :: timenow
330  LOGICAL :: rectilinear
331  LOGICAL :: l_dbg
332!---------------------------------------------------------------------
333  CALL ipsldbg (old_status=l_dbg)
334!-
335  nb_files = nb_files+1
336  pfileid  = nb_files
337!-
338! 1.0 Transfering into the common for future use
339!-
340  IF (l_dbg) WRITE(*,*) "histbeg_regular 1.0"
341!-
342  itau0(pfileid) = pitau0
343  date0(pfileid) = pdate0
344  deltat(pfileid) = pdeltat
345!-
346  IF (PRESENT(opt_rectilinear)) THEN
347    rectilinear = opt_rectilinear
348  ELSE
349    rectilinear = .FALSE.
350  ENDIF
351!-
352! 2.0 Initializes all variables for this file
353!-
354  IF (l_dbg) WRITE(*,*) "histbeg_regular 2.0"
355!-
356  IF (nb_files > nb_files_max) THEN
357    CALL ipslerr (3,"histbeg", &
358   &  'Table of files too small. You should increase nb_files_max', &
359   &  'in histcom.f90 in order to accomodate all these files',' ')
360  ENDIF
361!-
362  nb_var(pfileid) = 0
363  nb_tax(pfileid) = 0
364  nb_hax(pfileid) = 0
365  nb_zax(pfileid) = 0
366!-
367  slab_ori(pfileid,1:2) = (/ par_orix,par_oriy /)
368  slab_sz(pfileid,1:2)  = (/ par_szx, par_szy  /)
369!-
370! 3.0 Opening netcdf file and defining dimensions
371!-
372  IF (l_dbg) WRITE(*,*) "histbeg_regular 3.0"
373!-
374! Add DOMAIN number and ".nc" suffix in file name if needed
375!-
376  file = pfilename
377  CALL flio_dom_file (file,domain_id)
378!-
379  iret = NF90_CREATE (file,NF90_CLOBBER,ncid)
380!-
381  IF (rectilinear) THEN
382    iret = NF90_DEF_DIM (ncid,'lon',par_szx,xid(nb_files))
383    iret = NF90_DEF_DIM (ncid,'lat',par_szy,yid(nb_files))
384  ELSE
385    iret = NF90_DEF_DIM (ncid,'x',par_szx,xid(nb_files))
386    iret = NF90_DEF_DIM (ncid,'y',par_szy,yid(nb_files))
387  ENDIF
388!-
389! 4.0 Declaring the geographical coordinates and other attributes
390!-
391  IF (l_dbg) WRITE(*,*) "histbeg_regular 4.0"
392!-
393! 4.3 Global attributes
394!-
395  iret = NF90_PUT_ATT (ncid,NF90_GLOBAL,'Conventions','CF-1.1')
396  iret = NF90_PUT_ATT (ncid,NF90_GLOBAL,'file_name',TRIM(file))
397  iret = NF90_PUT_ATT (ncid,NF90_GLOBAL,'production',TRIM(model_name))
398  lock_modname = .TRUE.
399  CALL ioget_timestamp (timenow)
400  iret = NF90_PUT_ATT (ncid,NF90_GLOBAL,'TimeStamp',TRIM(timenow))
401!-
402! 5.0 Saving some important information on this file in the common
403!-
404  IF (l_dbg) WRITE(*,*) "histbeg_regular 5.0"
405!-
406  IF (PRESENT(domain_id)) THEN
407    dom_id_svg(pfileid) = domain_id
408  ENDIF
409  ncdf_ids(pfileid) = ncid
410  full_size(pfileid,1:2) = (/ pim,pjm /)
411!-
412! 6.0 storing the geographical coordinates
413!-
414  zoom(pfileid) = (pim /= par_szx).OR.(pjm /= par_szy)
415  regular(pfileid)=.TRUE.
416!-
417  CALL histhori_regular (pfileid,pim,plon,pjm,plat, &
418 &  ' ' ,'Default grid',phoriid,rectilinear)
419!-----------------------------
420END SUBROUTINE histbeg_regular
421!===
422SUBROUTINE histbeg_irregular &
423 &  (pfilename,pim,plon,plon_bounds,plat,plat_bounds,   &
424 &   pitau0,pdate0,pdeltat,phoriid,pfileid,domain_id)
425!---------------------------------------------------------------------
426!- This subroutine initializes a netcdf file and returns the ID.
427!- This version is for totaly irregular grids. In this case all
428!- all the data comes in as vectors and for the grid we have
429!- the coordinates of the 4 corners.
430!- It also gets the global parameters into the I/O subsystem.
431!-
432!- INPUT
433!-
434!- pfilename   : Name of the netcdf file to be created
435!- pim         : Size of arrays in longitude direction
436!- plon        : Coordinates of points in longitude
437!- plon_bounds : The 2 corners of the grid in longitude
438!- plat        : Coordinates of points in latitude
439!- plat_bounds : The 2 corners of the grid in latitude
440!-
441!- pitau0    : time step at which the history tape starts
442!- pdate0    : The Julian date at which the itau was equal to 0
443!- pdeltat   : Time step in seconds. Time step of the counter itau
444!-             used in histwrite for instance
445!-
446!- OUTPUT
447!-
448!- phoriid   : ID of the horizontal grid
449!- pfileid   : ID of the netcdf file
450!-
451!- Optional INPUT arguments
452!-
453!- domain_id       : Domain identifier
454!-
455!- TO DO
456!-
457!- This package should be written in F90 and use the following
458!- feature :
459!-   - structures for the meta-data of the files and variables
460!-   - memory allocation as needed
461!-   - Pointers
462!-
463!- VERSION
464!-
465!---------------------------------------------------------------------
466  IMPLICIT NONE
467!-
468  CHARACTER(LEN=*) :: pfilename
469  INTEGER,INTENT(IN) :: pim
470  REAL,DIMENSION(pim),INTENT(IN) :: plon,plat
471  REAL,DIMENSION(:,:),INTENT(IN) :: plon_bounds,plat_bounds
472  INTEGER,INTENT(IN) :: pitau0
473  REAL,INTENT(IN) :: pdate0,pdeltat
474  INTEGER,INTENT(OUT) :: pfileid,phoriid
475  INTEGER,INTENT(IN),OPTIONAL :: domain_id
476!-
477  INTEGER :: ncid,iret
478  CHARACTER(LEN=120) :: file
479  CHARACTER(LEN=30) :: timenow
480  LOGICAL :: l_dbg
481!---------------------------------------------------------------------
482  CALL ipsldbg (old_status=l_dbg)
483!-
484  nb_files = nb_files+1
485  pfileid  = nb_files
486!-
487! 1.0 Transfering into the common for future use
488!-
489  IF (l_dbg) WRITE(*,*) "histbeg_irregular 1.0"
490!-
491  itau0(pfileid) = pitau0
492  date0(pfileid) = pdate0
493  deltat(pfileid) = pdeltat
494!-
495! 2.0 Initializes all variables for this file
496!-
497  IF (l_dbg) WRITE(*,*) "histbeg_irregular 2.0"
498!-
499  IF (nb_files > nb_files_max) THEN
500    CALL ipslerr (3,"histbeg", &
501   &  'Table of files too small. You should increase nb_files_max', &
502   &  'in histcom.f90 in order to accomodate all these files',' ')
503  ENDIF
504!-
505  nb_var(pfileid) = 0
506  nb_tax(pfileid) = 0
507  nb_hax(pfileid) = 0
508  nb_zax(pfileid) = 0
509!-
510  slab_ori(pfileid,1:2) = (/ 1,1 /)
511  slab_sz(pfileid,1:2)  = (/ pim,1 /)
512!-
513! 3.0 Opening netcdf file and defining dimensions
514!-
515  IF (l_dbg) WRITE(*,*) "histbeg_irregular 3.0"
516!-
517! Add DOMAIN number and ".nc" suffix in file name if needed
518!-
519  file  = pfilename
520  CALL flio_dom_file (file,domain_id)
521!-
522  iret = NF90_CREATE (file,NF90_CLOBBER,ncid)
523!-
524  iret = NF90_DEF_DIM (ncid,'x',pim,xid(nb_files))
525  yid(nb_files) = 0
526!-
527! 4.0 Declaring the geographical coordinates and other attributes
528!-
529  IF (l_dbg) WRITE(*,*) "histbeg_irregular 4.0"
530!-
531! 4.3 Global attributes
532!-
533  iret = NF90_PUT_ATT (ncid,NF90_GLOBAL,'Conventions','CF-1.1')
534  iret = NF90_PUT_ATT (ncid,NF90_GLOBAL,'file_name',TRIM(file))
535  iret = NF90_PUT_ATT (ncid,NF90_GLOBAL,'production',TRIM(model_name))
536  lock_modname = .TRUE.
537  CALL ioget_timestamp (timenow)
538  iret = NF90_PUT_ATT (ncid,NF90_GLOBAL,'TimeStamp',TRIM(timenow))
539!-
540! 5.0 Saving some important information on this file in the common
541!-
542  IF (l_dbg) WRITE(*,*) "histbeg_irregular 5.0"
543!-
544  IF (PRESENT(domain_id)) THEN
545    dom_id_svg(pfileid) = domain_id
546  ENDIF
547  ncdf_ids(pfileid) = ncid
548  full_size(pfileid,1:2) = (/ pim,1 /)
549!-
550! 6.0 storing the geographical coordinates
551!-
552  zoom(pfileid)=.FALSE.
553  regular(pfileid)=.FALSE.
554!-
555  CALL histhori_irregular &
556 &  (pfileid,pim,plon,plon_bounds,plat,plat_bounds, &
557 &   ' ' ,'Default grid',phoriid)
558!-------------------------------
559END SUBROUTINE histbeg_irregular
560!===
561SUBROUTINE histhori_regular &
562 &  (pfileid,pim,plon,pjm,plat,phname,phtitle,phid,opt_rectilinear)
563!---------------------------------------------------------------------
564!- This subroutine is made to declare a new horizontale grid.
565!- It has to have the same number of points as
566!- the original and thus in this routine we will only
567!- add two variable (longitude and latitude).
568!- Any variable in the file can thus point to this pair
569!- through an attribute. This routine is very usefull
570!- to allow staggered grids.
571!-
572!- INPUT
573!-
574!- pfileid : The id of the file to which the grid should be added
575!- pim     : Size in the longitude direction
576!- plon    : The longitudes
577!- pjm     : Size in the latitude direction
578!- plat    : The latitudes
579!- phname  : The name of grid
580!- phtitle : The title of the grid
581!-
582!- OUTPUT
583!-
584!- phid    : Id of the created grid
585!-
586!- OPTIONAL
587!-
588!- opt_rectilinear : If true we know the grid is rectilinear.
589!-
590!---------------------------------------------------------------------
591  IMPLICIT NONE
592!-
593  INTEGER,INTENT(IN) :: pfileid,pim,pjm
594  REAL,INTENT(IN),DIMENSION(pim,pjm) :: plon,plat
595  CHARACTER(LEN=*),INTENT(IN) :: phname,phtitle
596  INTEGER,INTENT(OUT) :: phid
597  LOGICAL,INTENT(IN),OPTIONAL :: opt_rectilinear
598!-
599  CHARACTER(LEN=25) :: lon_name,lat_name
600  CHARACTER(LEN=80) :: tmp_title,tmp_name
601  INTEGER :: ndim
602  INTEGER,DIMENSION(2) :: dims
603  INTEGER :: nlonid,nlatid
604  INTEGER :: orix,oriy,par_szx,par_szy
605  INTEGER :: iret,ncid
606  LOGICAL :: rectilinear
607  LOGICAL :: l_dbg
608!---------------------------------------------------------------------
609  CALL ipsldbg (old_status=l_dbg)
610!-
611! 1.0 Check that all fits in the buffers
612!-
613  IF (    (pim /= full_size(pfileid,1)) &
614    & .OR.(pjm /= full_size(pfileid,2)) ) THEN
615    CALL ipslerr (3,"histhori", &
616   &  'The new horizontal grid does not have the same size', &
617   &  'as the one provided to histbeg. This is not yet ', &
618   &  'possible in the hist package.')
619  ENDIF
620!-
621  IF (PRESENT(opt_rectilinear)) THEN
622    rectilinear = opt_rectilinear
623  ELSE
624    rectilinear = .FALSE.
625  ENDIF
626!-
627! 1.1 Create all the variables needed
628!-
629  IF (l_dbg) WRITE(*,*) "histhori_regular 1.0"
630!-
631  ncid = ncdf_ids(pfileid)
632!-
633  ndim = 2
634  dims(1:2) = (/ xid(pfileid),yid(pfileid) /)
635!-
636  tmp_name = phname
637  IF (rectilinear) THEN
638    IF (nb_hax(pfileid) == 0) THEN
639      lon_name = 'lon'
640      lat_name = 'lat'
641    ELSE
642      lon_name = 'lon_'//TRIM(tmp_name)
643      lat_name = 'lat_'//TRIM(tmp_name)
644    ENDIF
645  ELSE
646    IF (nb_hax(pfileid) == 0) THEN
647      lon_name = 'nav_lon'
648      lat_name = 'nav_lat'
649    ELSE
650      lon_name = 'nav_lon_'//TRIM(tmp_name)
651      lat_name = 'nav_lat_'//TRIM(tmp_name)
652    ENDIF
653  ENDIF
654!-
655! 1.2 Save the informations
656!-
657  phid =  nb_hax(pfileid)+1
658  nb_hax(pfileid) = phid
659!-
660  hax_name(pfileid,phid,1:2) = (/ lon_name,lat_name /)
661  tmp_title = phtitle
662!-
663! 2.0 Longitude
664!-
665  IF (l_dbg) WRITE(*,*) "histhori_regular 2.0"
666!-
667  IF (rectilinear) THEN
668    ndim = 1
669    dims(1:1) = (/ xid(pfileid) /)
670  ENDIF
671  iret = NF90_DEF_VAR (ncid,lon_name,NF90_FLOAT,dims(1:ndim),nlonid)
672  IF (rectilinear) THEN
673    iret = NF90_PUT_ATT (ncid,nlonid,'axis',"X")
674  ENDIF
675  iret = NF90_PUT_ATT (ncid,nlonid,'standard_name',"longitude")
676  iret = NF90_PUT_ATT (ncid,nlonid,'units',"degrees_east")
677  iret = NF90_PUT_ATT (ncid,nlonid,'valid_min', &
678 &                     REAL(MINVAL(plon),KIND=4))
679  iret = NF90_PUT_ATT (ncid,nlonid,'valid_max', &
680 &                     REAL(MAXVAL(plon),KIND=4))
681  iret = NF90_PUT_ATT (ncid,nlonid,'long_name',"Longitude")
682  iret = NF90_PUT_ATT (ncid,nlonid,'nav_model',TRIM(tmp_title))
683!-
684! 3.0 Latitude
685!-
686  IF (l_dbg) WRITE(*,*) "histhori_regular 3.0"
687!-
688  IF (rectilinear) THEN
689    ndim = 1
690    dims(1:1) = (/ yid(pfileid) /)
691  ENDIF
692  iret = NF90_DEF_VAR (ncid,lat_name,NF90_FLOAT,dims(1:ndim),nlatid)
693  IF (rectilinear) THEN
694    iret = NF90_PUT_ATT (ncid,nlatid,'axis',"Y")
695  ENDIF
696  iret = NF90_PUT_ATT (ncid,nlatid,'standard_name',"latitude")
697  iret = NF90_PUT_ATT (ncid,nlatid,'units',"degrees_north")
698  iret = NF90_PUT_ATT (ncid,nlatid,'valid_min', &
699 &                     REAL(MINVAL(plat),KIND=4))
700  iret = NF90_PUT_ATT (ncid,nlatid,'valid_max', &
701 &                     REAL(MAXVAL(plat),KIND=4))
702  iret = NF90_PUT_ATT (ncid,nlatid,'long_name',"Latitude")
703  iret = NF90_PUT_ATT (ncid,nlatid,'nav_model',TRIM(tmp_title))
704!-
705  iret = NF90_ENDDEF (ncid)
706!-
707! 4.0 storing the geographical coordinates
708!-
709  IF (l_dbg) WRITE(*,*) "histhori_regular 4.0"
710!-
711  orix = slab_ori(pfileid,1)
712  oriy = slab_ori(pfileid,2)
713  par_szx = slab_sz(pfileid,1)
714  par_szy = slab_sz(pfileid,2)
715!-
716! Transfer the longitude
717!-
718  IF (rectilinear) THEN
719    iret = NF90_PUT_VAR (ncid,nlonid,plon(orix:orix+par_szx-1,1))
720  ELSE
721    iret = NF90_PUT_VAR (ncid,nlonid, &
722 &           plon(orix:orix+par_szx-1,oriy:oriy+par_szy-1))
723  ENDIF
724!-
725! Transfer the latitude
726!-
727  IF (rectilinear) THEN
728    iret = NF90_PUT_VAR (ncid,nlatid,plat(1,oriy:oriy+par_szy-1))
729  ELSE
730    iret = NF90_PUT_VAR (ncid,nlatid, &
731 &           plat(orix:orix+par_szx-1,oriy:oriy+par_szy-1))
732  ENDIF
733!-
734  iret = NF90_REDEF (ncid)
735!------------------------------
736END SUBROUTINE histhori_regular
737!===
738SUBROUTINE histhori_irregular &
739 &  (pfileid,pim,plon,plon_bounds,plat,plat_bounds, &
740 &   phname,phtitle,phid)
741!---------------------------------------------------------------------
742!- This subroutine is made to declare a new horizontale grid.
743!- It has to have the same number of points as
744!- the original and thus in this routine we will only
745!- add two variable (longitude and latitude).
746!- Any variable in the file can thus point to this pair
747!- through an attribute. This routine is very usefull
748!- to allow staggered grids.
749!-
750!- INPUT
751!-
752!- pfileid     : The id of the file to which the grid should be added
753!- pim         : Size in the longitude direction
754!- plon        : The longitudes
755!- plon_bounds : The boundaries of the grid in longitude
756!- plat        : The latitudes
757!- plat_bounds : Boundaries of the grid in latitude
758!- phname      : The name of grid
759!- phtitle     : The title of the grid
760!-
761!- OUTPUT
762!-
763!- phid    : Id of the created grid
764!---------------------------------------------------------------------
765  IMPLICIT NONE
766!-
767  INTEGER,INTENT(IN) :: pfileid,pim
768  REAL,DIMENSION(pim),INTENT(IN) :: plon,plat
769  REAL,DIMENSION(:,:),INTENT(IN) :: plon_bounds,plat_bounds
770  CHARACTER(LEN=*),INTENT(IN) :: phname,phtitle
771  INTEGER,INTENT(OUT) :: phid
772!-
773  CHARACTER(LEN=25) :: lon_name,lat_name
774  CHARACTER(LEN=30) :: lonbound_name,latbound_name
775  CHARACTER(LEN=80) :: tmp_title,tmp_name,longname
776  INTEGER :: ndim,dims(2)
777  INTEGER :: ndimb,dimsb(2)
778  INTEGER :: nbbounds
779  INTEGER :: nlonid,nlatid,nlonidb,nlatidb
780  INTEGER :: iret,ncid,twoid
781  LOGICAL :: transp = .FALSE.
782  REAL,ALLOCATABLE,DIMENSION(:,:) :: bounds_trans
783  LOGICAL :: l_dbg
784!---------------------------------------------------------------------
785  CALL ipsldbg (old_status=l_dbg)
786!-
787! 1.0 Check that all fits in the buffers
788!-
789  IF (    (pim /= full_size(pfileid,1)) &
790    & .OR.(full_size(pfileid,2) /= 1) ) THEN
791    CALL ipslerr (3,"histhori", &
792   &  'The new horizontal grid does not have the same size', &
793   &  'as the one provided to histbeg. This is not yet ', &
794   &  'possible in the hist package.')
795  ENDIF
796!-
797! 1.1 Create all the variables needed
798!-
799  IF (l_dbg) WRITE(*,*) 'histhori_irregular 1.0'
800!-
801  ncid = ncdf_ids(pfileid)
802!-
803  IF     (SIZE(plon_bounds,DIM=1) == pim) THEN
804    nbbounds = SIZE(plon_bounds,DIM=2)
805    transp = .TRUE.
806  ELSEIF (SIZE(plon_bounds,DIM=2) == pim) THEN
807    nbbounds = SIZE(plon_bounds,DIM=1)
808    transp = .FALSE.
809  ELSE
810    CALL ipslerr (3,"histhori", &
811   &  'The boundary variable does not have any axis corresponding', &
812   &  'to the size of the longitude or latitude variable', &
813   &  '.')
814  ENDIF
815!-
816  ALLOCATE(bounds_trans(nbbounds,pim))
817!-
818  iret = NF90_DEF_DIM (ncid,'nbnd',nbbounds,twoid)
819  ndim = 1
820  dims(1) = xid(pfileid)
821  ndimb = 2
822  dimsb(1:2) = (/ twoid,xid(pfileid) /)
823!-
824  tmp_name = phname
825  IF (nb_hax(pfileid) == 0) THEN
826    lon_name = 'nav_lon'
827    lat_name = 'nav_lat'
828  ELSE
829    lon_name = 'nav_lon_'//TRIM(tmp_name)
830    lat_name = 'nav_lat_'//TRIM(tmp_name)
831  ENDIF
832  lonbound_name = TRIM(lon_name)//'_bounds'
833  latbound_name = TRIM(lat_name)//'_bounds'
834!-
835! 1.2 Save the informations
836!-
837  phid =  nb_hax(pfileid)+1
838  nb_hax(pfileid) = phid
839!-
840  hax_name(pfileid,phid,1:2) = (/ lon_name,lat_name /)
841  tmp_title = phtitle
842!-
843! 2.0 Longitude
844!-
845  IF (l_dbg) WRITE(*,*) "histhori_irregular 2.0"
846!-
847  iret = NF90_DEF_VAR (ncid,lon_name,NF90_FLOAT,dims(1:ndim),nlonid)
848  iret = NF90_PUT_ATT (ncid,nlonid,'standard_name',"longitude")
849  iret = NF90_PUT_ATT (ncid,nlonid,'units',"degrees_east")
850  iret = NF90_PUT_ATT (ncid,nlonid,'valid_min', &
851 &                     REAL(MINVAL(plon),KIND=4))
852  iret = NF90_PUT_ATT (ncid,nlonid,'valid_max', &
853 &                     REAL(MAXVAL(plon),KIND=4))
854  iret = NF90_PUT_ATT (ncid,nlonid,'long_name',"Longitude")
855  iret = NF90_PUT_ATT (ncid,nlonid,'nav_model',TRIM(tmp_title))
856!-
857! 2.1 Longitude bounds
858!-
859  iret = NF90_PUT_ATT (ncid,nlonid,'bounds',TRIM(lonbound_name))
860  iret = NF90_DEF_VAR (ncid,lonbound_name,NF90_FLOAT, &
861 &                     dimsb(1:ndimb),nlonidb)
862  longname = 'Boundaries for coordinate variable '//TRIM(lon_name)
863  iret = NF90_PUT_ATT (ncid,nlonidb,'long_name',TRIM(longname))
864!-
865! 3.0 Latitude
866!-
867  IF (l_dbg) WRITE(*,*) "histhori_irregular 3.0"
868!-
869  iret = NF90_DEF_VAR (ncid,lat_name,NF90_FLOAT,dims(1:ndim),nlatid)
870  iret = NF90_PUT_ATT (ncid,nlatid,'standard_name',"latitude")
871  iret = NF90_PUT_ATT (ncid,nlatid,'units',"degrees_north")
872  iret = NF90_PUT_ATT (ncid,nlatid,'valid_min', &
873 &                     REAL(MINVAL(plat),KIND=4))
874  iret = NF90_PUT_ATT (ncid,nlatid,'valid_max', &
875 &                     REAL(MAXVAL(plat),KIND=4))
876  iret = NF90_PUT_ATT (ncid,nlatid,'long_name',"Latitude")
877  iret = NF90_PUT_ATT (ncid,nlatid,'nav_model',TRIM(tmp_title))
878!-
879! 3.1 Latitude bounds
880!-
881  iret = NF90_PUT_ATT (ncid,nlatid,'bounds',TRIM(latbound_name))
882  iret = NF90_DEF_VAR (ncid,latbound_name,NF90_FLOAT, &
883 &                     dimsb(1:ndimb),nlatidb)
884  longname = 'Boundaries for coordinate variable '//TRIM(lat_name)
885  iret = NF90_PUT_ATT (ncid,nlatidb,'long_name',TRIM(longname))
886!-
887  iret = NF90_ENDDEF (ncid)
888!-
889! 4.0 storing the geographical coordinates
890!-
891  IF (l_dbg) WRITE(*,*) "histhori_irregular 4.0"
892!-
893! 4.1 Write the longitude
894!-
895  iret = NF90_PUT_VAR (ncid,nlonid,plon(1:pim))
896!-
897! 4.2 Write the longitude bounds
898!-
899  IF (transp) THEN
900    bounds_trans = TRANSPOSE(plon_bounds)
901  ELSE
902    bounds_trans = plon_bounds
903  ENDIF
904  iret = NF90_PUT_VAR (ncid,nlonidb,bounds_trans(1:nbbounds,1:pim))
905!-
906! 4.3 Write the latitude
907!-
908  iret = NF90_PUT_VAR (ncid,nlatid,plat(1:pim))
909!-
910! 4.4 Write the latitude bounds
911!-
912  IF (transp) THEN
913    bounds_trans = TRANSPOSE(plat_bounds)
914  ELSE
915    bounds_trans = plat_bounds
916  ENDIF
917  iret = NF90_PUT_VAR (ncid,nlatidb,bounds_trans(1:nbbounds,1:pim))
918!-
919  DEALLOCATE(bounds_trans)
920!-
921  iret = NF90_REDEF (ncid)
922!--------------------------------
923END SUBROUTINE histhori_irregular
924!===
925SUBROUTINE histvert (pfileid,pzaxname,pzaxtitle,pzaxunit, &
926 &                   pzsize,pzvalues,pzaxid,pdirect)
927!---------------------------------------------------------------------
928!- This subroutine defines a vertical axis and returns it s id.
929!- It gives the user the possibility to the user to define many
930!- different vertical axes. For each variable defined with histdef a
931!- vertical axis can be specified with by it s ID.
932!-
933!- INPUT
934!-
935!- pfileid  : ID of the file the variable should be archived in
936!- pzaxname : Name of the vertical axis
937!- pzaxtitle: title of the vertical axis
938!- pzaxunit : Units of the vertical axis (no units if blank string)
939!- pzsize   : size of the vertical axis
940!- pzvalues : Coordinate values of the vetical axis
941!-
942!- pdirect  : is an optional argument which allows to specify the
943!-            the positive direction of the axis : up or down.
944!- OUTPUT
945!-
946!- pzaxid   : Returns the ID of the axis.
947!-            Note that this is not the netCDF ID !
948!-
949!- VERSION
950!-
951!---------------------------------------------------------------------
952  IMPLICIT NONE
953!-
954  INTEGER,INTENT(IN) :: pfileid,pzsize
955  CHARACTER(LEN=*),INTENT(IN) :: pzaxname,pzaxunit,pzaxtitle
956  REAL,INTENT(IN) :: pzvalues(pzsize)
957  INTEGER,INTENT(OUT) :: pzaxid
958  CHARACTER(LEN=*),INTENT(IN),OPTIONAL :: pdirect
959!-
960  INTEGER :: pos,iv,zdimid,zaxid_tmp
961  CHARACTER(LEN=70) :: str71
962  CHARACTER(LEN=80) :: str80
963  CHARACTER(LEN=20) :: direction
964  INTEGER :: iret,leng,ncid
965  LOGICAL :: l_dbg
966!---------------------------------------------------------------------
967  CALL ipsldbg (old_status=l_dbg)
968!-
969! 1.0 Verifications :
970!    Do we have enough space for an extra axis ?
971!    Is the name already in use ?
972!-
973  IF (l_dbg) WRITE(*,*) "histvert : 1.0 Verifications", &
974 &                      pzaxname,'---',pzaxunit,'---',pzaxtitle
975!-
976! - Direction of axis. Can we get if from the user.
977!   If not we put unknown.
978!-
979  IF (PRESENT(pdirect)) THEN
980    direction = TRIM(pdirect)
981    CALL strlowercase (direction)
982  ELSE
983    direction = 'unknown'
984  ENDIF
985!-
986! Check the consistency of the attribute
987!-
988  IF (     (direction /= 'unknown') &
989 &    .AND.(direction /= 'up')      &
990 &    .AND.(direction /= 'down')   ) THEN
991    direction = 'unknown'
992    str80 = 'The specified axis was : '//TRIM(direction)
993    CALL ipslerr (2,"histvert",&
994   & "The specified direction for the vertical axis is not possible.",&
995   & "it is replaced by : unknown",str80)
996  ENDIF
997!-
998  IF (nb_zax(pfileid)+1 > nb_zax_max) THEN
999    CALL ipslerr (3,"histvert", &
1000   &  'Table of vertical axes too small. You should increase ',&
1001   &  'nb_zax_max in histcom.f90 in order to accomodate all ', &
1002   &  'these variables ')
1003  ENDIF
1004!-
1005  iv = nb_zax(pfileid)
1006  IF (iv > 1) THEN
1007    CALL find_str (zax_name(pfileid,1:iv-1),pzaxname,pos)
1008  ELSE
1009    pos = 0
1010  ENDIF
1011!-
1012  IF (pos > 0) THEN
1013    WRITE(str71,'("Check variable ",A," in file",I3)') &
1014 &    TRIM(pzaxname),pfileid
1015    CALL ipslerr (3,"histvert", &
1016 &    "Vertical axis already exists",TRIM(str71), &
1017 &    "Can also be a wrong file ID in another declaration")
1018  ENDIF
1019!-
1020  iv = nb_zax(pfileid)+1
1021!-
1022! 2.0 Add the information to the file
1023!-
1024  IF (l_dbg) &
1025 &  WRITE(*,*) "histvert : 2.0 Add the information to the file"
1026!-
1027  ncid = ncdf_ids(pfileid)
1028!-
1029  leng = MIN(LEN_TRIM(pzaxname),20)
1030  iret = NF90_DEF_DIM (ncid,pzaxname(1:leng),pzsize,zaxid_tmp)
1031  iret = NF90_DEF_VAR (ncid,pzaxname(1:leng),NF90_FLOAT, &
1032 &                     zaxid_tmp,zdimid)
1033  iret = NF90_PUT_ATT (ncid,zdimid,'axis',"Z")
1034  iret = NF90_PUT_ATT (ncid,zdimid,'standard_name',"model_level_number")
1035  leng = MIN(LEN_TRIM(pzaxunit),20)
1036  IF (leng > 0) THEN
1037    iret = NF90_PUT_ATT (ncid,zdimid,'units',pzaxunit(1:leng))
1038  ENDIF
1039  iret = NF90_PUT_ATT (ncid,zdimid,'positive',TRIM(direction))
1040  iret = NF90_PUT_ATT (ncid,zdimid,'valid_min', &
1041 &                     REAL(MINVAL(pzvalues(1:pzsize)),KIND=4))
1042  iret = NF90_PUT_ATT (ncid,zdimid,'valid_max', &
1043 &                     REAL(MAXVAL(pzvalues(1:pzsize)),KIND=4))
1044  leng = MIN(LEN_TRIM(pzaxname),20)
1045  iret = NF90_PUT_ATT (ncid,zdimid,'title',pzaxname(1:leng))
1046  leng = MIN(LEN_TRIM(pzaxtitle),80)
1047  iret = NF90_PUT_ATT (ncid,zdimid,'long_name',pzaxtitle(1:leng))
1048!-
1049  iret = NF90_ENDDEF (ncid)
1050!-
1051  iret = NF90_PUT_VAR (ncid,zdimid,pzvalues(1:pzsize))
1052!-
1053  iret = NF90_REDEF (ncid)
1054!-
1055!- 3.0 add the information to the common
1056!-
1057  IF (l_dbg) &
1058  &  WRITE(*,*) "histvert : 3.0 add the information to the common"
1059!-
1060  nb_zax(pfileid) = iv
1061  zax_size(pfileid,iv) = pzsize
1062  zax_name(pfileid,iv) = pzaxname
1063  zax_ids(pfileid,iv) = zaxid_tmp
1064  pzaxid =  iv
1065!----------------------
1066END SUBROUTINE histvert
1067!===
1068SUBROUTINE histdef (pfileid,pvarname,ptitle,punit, &
1069 &                  pxsize,pysize,phoriid,pzsize,  &
1070 &                  par_oriz,par_szz,pzid,         &
1071 &                  pnbbyt,popp,pfreq_opp,pfreq_wrt,var_range)
1072!---------------------------------------------------------------------
1073!- With this subroutine each variable to be archived on the history
1074!- tape should be declared.
1075!-
1076!- It gives the user the choise of operation
1077!- to be performed on the variables, the frequency of this operation
1078!- and finaly the frequency of the archiving.
1079!-
1080!- INPUT
1081!-
1082!- pfileid  : ID of the file the variable should be archived in
1083!- pvarname : Name of the variable, short and easy to remember
1084!- ptitle   : Full name of the variable
1085!- punit    : Units of the variable (no units if blank string)
1086!-
1087!- The next 3 arguments give the size of that data
1088!- that will be passed to histwrite. The zoom will be
1089!- done there with the horizontal information obtained
1090!- in histbeg and the vertical information to follow.
1091!-
1092!- pxsize   : Size in X direction (size of the data that will be
1093!-            given to histwrite)
1094!- pysize   : Size in Y direction
1095!- phoriid  : ID of the horizontal axis
1096!-
1097!- The next two arguments give the vertical zoom to use.
1098!-
1099!- pzsize   : Size in Z direction (If 1 then no axis is declared
1100!-            for this variable and pzid is not used)
1101!- par_oriz : Off set of the zoom
1102!- par_szz  : Size of the zoom
1103!-
1104!- pzid     : ID of the vertical axis to use. It has to have
1105!-            the size of the zoom.
1106!- pnbbyt   : Number of bytes on which to store in netCDF (Not opp.)
1107!- popp     : Operation to be performed. The following options
1108!-            exist today :
1109!- inst : keeps instantaneous values for writting
1110!- ave  : Computes the average from call between writes
1111!- pfreq_opp: Frequency of this operation (in seconds)
1112!- pfreq_wrt: Frequency at which the variable should be
1113!-            written (in seconds)
1114!- var_range: Range of the variable.
1115!-            If the minimum is greater than the maximum,
1116!-            the values will be calculated.
1117!-
1118!- VERSION
1119!---------------------------------------------------------------------
1120  IMPLICIT NONE
1121!-
1122  INTEGER,INTENT(IN) :: pfileid,pxsize,pysize,pzsize,pzid
1123  INTEGER,INTENT(IN) :: par_oriz,par_szz,pnbbyt,phoriid
1124  CHARACTER(LEN=*),INTENT(IN) :: pvarname,punit,popp,ptitle
1125  REAL,INTENT(IN) :: pfreq_opp,pfreq_wrt
1126  REAL,DIMENSION(2),OPTIONAL,INTENT(IN) :: var_range
1127!-
1128  INTEGER :: iv,i
1129  CHARACTER(LEN=70) :: str70,str71,str72
1130  CHARACTER(LEN=20) :: tmp_name
1131  CHARACTER(LEN=40) :: str40
1132  CHARACTER(LEN=10) :: str10
1133  CHARACTER(LEN=80) :: tmp_str80
1134  CHARACTER(LEN=7) :: tmp_topp,tmp_sopp(nbopp_max)
1135  CHARACTER(LEN=120) :: ex_topps
1136  REAL :: tmp_scal(nbopp_max),un_an,un_jour,test_fopp,test_fwrt
1137  INTEGER :: pos,buff_sz
1138  LOGICAL :: l_dbg
1139!---------------------------------------------------------------------
1140  CALL ipsldbg (old_status=l_dbg)
1141!-
1142  ex_topps = 'ave, inst, t_min, t_max, t_sum, once, never, l_max, l_min'
1143!-
1144  nb_var(pfileid) = nb_var(pfileid)+1
1145  iv = nb_var(pfileid)
1146!-
1147  IF (iv > nb_var_max) THEN
1148    CALL ipslerr (3,"histdef", &
1149   &  'Table of variables too small. You should increase nb_var_max',&
1150   &  'in histcom.f90 in order to accomodate all these variables', &
1151   &  ' ')
1152  ENDIF
1153!-
1154! 1.0 Transfer informations on the variable to the common
1155!     and verify that it does not already exist
1156!-
1157  IF (l_dbg) WRITE(*,*) "histdef : 1.0"
1158!-
1159  IF (iv > 1) THEN
1160    CALL find_str (name(pfileid,1:iv-1),pvarname,pos)
1161  ELSE
1162    pos = 0
1163  ENDIF
1164!-
1165  IF (pos > 0) THEN
1166    str70 = "Variable already exists"
1167    WRITE(str71,'("Check variable  ",a," in file",I3)') &
1168 &    TRIM(pvarname),pfileid
1169    str72 = "Can also be a wrong file ID in another declaration"
1170    CALL ipslerr (3,"histdef",str70,str71,str72)
1171  ENDIF
1172!-
1173  name(pfileid,iv) = pvarname
1174  title(pfileid,iv) = ptitle
1175  unit_name(pfileid,iv) = punit
1176  tmp_name =  name(pfileid,iv)
1177!-
1178! 1.1 decode the operations
1179!-
1180  fullop(pfileid,iv) = popp
1181  tmp_str80 = popp
1182  CALL buildop &
1183 &  (tmp_str80,ex_topps,tmp_topp,nbopp_max,missing_val, &
1184 &   tmp_sopp,tmp_scal,nbopp(pfileid,iv))
1185!-
1186  topp(pfileid,iv) = tmp_topp
1187  DO i=1,nbopp(pfileid,iv)
1188    sopps(pfileid,iv,i) = tmp_sopp(i)
1189    scal(pfileid,iv,i) = tmp_scal(i)
1190  ENDDO
1191!-
1192! 1.2 If we have an even number of operations
1193!     then we need to add identity
1194!-
1195  IF (2*INT(nbopp(pfileid,iv)/2.0) == nbopp(pfileid,iv)) THEN
1196    nbopp(pfileid,iv) = nbopp(pfileid,iv)+1
1197    sopps(pfileid,iv,nbopp(pfileid,iv)) = 'ident'
1198    scal(pfileid,iv,nbopp(pfileid,iv)) = missing_val
1199  ENDIF
1200!-
1201! 2.0 Put the size of the variable in the common and check
1202!-
1203  IF (l_dbg) &
1204 &  WRITE(*,*) "histdef : 2.0",pfileid,iv,nbopp(pfileid,iv), &
1205 &    sopps(pfileid,iv,1:nbopp(pfileid,iv)), &
1206 &    scal(pfileid,iv,1:nbopp(pfileid,iv))
1207!-
1208  scsize(pfileid,iv,1:3) = (/ pxsize,pysize,pzsize /)
1209!-
1210  zorig(pfileid,iv,1:3) = &
1211 &  (/ slab_ori(pfileid,1),slab_ori(pfileid,2),par_oriz /)
1212!-
1213  zsize(pfileid,iv,1:3) = &
1214 &  (/ slab_sz(pfileid,1),slab_sz(pfileid,2),par_szz /)
1215!-
1216! Is the size of the full array the same as that of the coordinates  ?
1217!-
1218  IF (    (pxsize > full_size(pfileid,1)) &
1219 &    .OR.(pysize > full_size(pfileid,2)) ) THEN
1220!-
1221    str70 = "The size of the variable is different "// &
1222 &          "from the one of the coordinates"
1223    WRITE(str71,'("Size of coordinates :",2I4)') &
1224 &   full_size(pfileid,1),full_size(pfileid,2)
1225    WRITE(str72,'("Size declared for variable ",a," :",2I4)') &
1226 &   TRIM(tmp_name),pxsize,pysize
1227    CALL ipslerr (3,"histdef",str70,str71,str72)
1228  ENDIF
1229!-
1230! Is the size of the zoom smaler than the coordinates ?
1231!-
1232  IF (    (full_size(pfileid,1) < slab_sz(pfileid,1)) &
1233 &    .OR.(full_size(pfileid,2) < slab_sz(pfileid,2)) ) THEN
1234    str70 = &
1235 &   "Size of variable should be greater or equal to those of the zoom"
1236    WRITE(str71,'("Size of XY zoom :",2I4)') &
1237 &   slab_sz(pfileid,1),slab_sz(pfileid,1)
1238    WRITE(str72,'("Size declared for variable ",a," :",2I4)') &
1239 &   TRIM(tmp_name),pxsize,pysize
1240    CALL ipslerr (3,"histdef",str70,str71,str72)
1241  ENDIF
1242!-
1243! 2.1 We store the horizontal grid information with minimal
1244!     and a fall back onto the default grid
1245!-
1246  IF ( (phoriid > 0).AND.(phoriid <= nb_hax(pfileid)) ) THEN
1247    var_haxid(pfileid,iv) = phoriid
1248  ELSE
1249    var_haxid(pfileid,iv) = 1
1250    CALL ipslerr (2,"histdef", &
1251   &  'We use the default grid for variable as an invalide',&
1252   &  'ID was provided for variable : ',TRIM(pvarname))
1253  ENDIF
1254!-
1255! 2.2 Check the vertical coordinates if needed
1256!-
1257  IF (par_szz > 1) THEN
1258!-
1259!-- Does the vertical coordinate exist ?
1260!-
1261    IF (pzid > nb_zax(pfileid)) THEN
1262      WRITE(str70, &
1263 &    '("The vertical coordinate chosen for variable ",a)') &
1264 &     TRIM(tmp_name)
1265      str71 = " Does not exist."
1266      CALL ipslerr (3,"histdef",str70,str71," ")
1267    ENDIF
1268!-
1269!-- Is the vertical size of the variable equal to that of the axis ?
1270!-
1271    IF (par_szz /= zax_size(pfileid,pzid)) THEN
1272      str70 = "The size of the zoom does not correspond "// &
1273 &            "to the size of the chosen vertical axis"
1274      WRITE(str71,'("Size of zoom in z :",I4)') par_szz
1275      WRITE(str72,'("Size declared for axis ",A," :",I4)') &
1276 &     TRIM(zax_name(pfileid,pzid)),zax_size(pfileid,pzid)
1277      CALL ipslerr (3,"histdef",str70,str71,str72)
1278    ENDIF
1279!-
1280!-- Is the zoom smaler that the total size of the variable ?
1281!-
1282    IF (pzsize < par_szz) THEN
1283      str70 = "The vertical size of variable "// &
1284 &            "is smaller than that of the zoom."
1285      WRITE(str71,'("Declared vertical size of data :",I5)') pzsize
1286      WRITE(str72,'("Size of zoom for variable ",a," = ",I5)') &
1287 &     TRIM(tmp_name),par_szz
1288      CALL ipslerr (3,"histdef",str70,str71,str72)
1289    ENDIF
1290    var_zaxid(pfileid,iv) = pzid
1291  ELSE
1292    var_zaxid(pfileid,iv) = -99
1293  ENDIF
1294!-
1295! 3.0 Determine the position of the variable in the buffer
1296!     If it is instantaneous output then we do not use the buffer
1297!-
1298  IF (l_dbg) WRITE(*,*) "histdef : 3.0"
1299!-
1300! 3.1 We get the size of the arrays histwrite will get and check
1301!     that they fit into the tmp_buffer
1302!-
1303  buff_sz = zsize(pfileid,iv,1)*zsize(pfileid,iv,2)*zsize(pfileid,iv,3)
1304!-
1305! 3.2 move the pointer of the buffer array for operation
1306!     which need bufferisation
1307!-
1308  IF (     (TRIM(tmp_topp) /= "inst") &
1309 &    .AND.(TRIM(tmp_topp) /= "once") &
1310 &    .AND.(TRIM(tmp_topp) /= "never") )THEN
1311    point(pfileid,iv) = buff_pos+1
1312    buff_pos = buff_pos+buff_sz
1313    IF (l_dbg) THEN
1314      WRITE(*,*) "histdef : 3.2 bufpos for iv = ",iv, &
1315 &               " pfileid = ",pfileid," is = ",point(pfileid,iv)
1316    ENDIF
1317  ENDIF
1318!-
1319! 4.0 Transfer the frequency of the operations and check
1320!     for validity. We have to pay attention to negative values
1321!     of the frequency which indicate monthly time-steps.
1322!     The strategy is to bring it back to seconds for the tests
1323!-
1324  IF (l_dbg) WRITE(*,*) "histdef : 4.0"
1325!-
1326  freq_opp(pfileid,iv) = pfreq_opp
1327  freq_wrt(pfileid,iv) = pfreq_wrt
1328!-
1329  CALL ioget_calendar(un_an,un_jour)
1330  IF (pfreq_opp < 0) THEN
1331    CALL ioget_calendar(un_an)
1332    test_fopp = pfreq_opp*(-1.)*un_an/12.*un_jour
1333  ELSE
1334    test_fopp = pfreq_opp
1335  ENDIF
1336  IF (pfreq_wrt < 0) THEN
1337    CALL ioget_calendar(un_an)
1338    test_fwrt = pfreq_wrt*(-1.)*un_an/12.*un_jour
1339  ELSE
1340    test_fwrt = pfreq_wrt
1341  ENDIF
1342!-
1343! 4.1 Frequency of operations and output should be larger than deltat !
1344!-
1345  IF (test_fopp < deltat(pfileid)) THEN
1346    str70 = 'Frequency of operations should be larger than deltat'
1347    WRITE(str71,'("It is not the case for variable ",a," :",F10.4)') &
1348 &   TRIM(tmp_name),pfreq_opp
1349    str72 = "PATCH : frequency set to deltat"
1350!-
1351    CALL ipslerr (2,"histdef",str70,str71,str72)
1352!-
1353    freq_opp(pfileid,iv) = deltat(pfileid)
1354  ENDIF
1355!-
1356  IF (test_fwrt < deltat(pfileid)) THEN
1357    str70 = 'Frequency of output should be larger than deltat'
1358    WRITE(str71,'("It is not the case for variable ",a," :",F10.4)') &
1359 &   TRIM(tmp_name),pfreq_wrt
1360    str72 = "PATCH : frequency set to deltat"
1361!-
1362    CALL ipslerr (2,"histdef",str70,str71,str72)
1363!-
1364    freq_wrt(pfileid,iv) = deltat(pfileid)
1365  ENDIF
1366!-
1367! 4.2 First the existence of the operation is tested and then
1368!     its compaticility with the choice of frequencies
1369!-
1370  IF (TRIM(tmp_topp) == "inst") THEN
1371    IF (test_fopp /= test_fwrt) THEN
1372      str70 = 'For instantaneous output the frequency '// &
1373 &            'of operations and output'
1374      WRITE(str71, &
1375 &     '("should be the same, this was not case for variable ",a)') &
1376 &      TRIM(tmp_name)
1377      str72 = "PATCH : The smalest frequency of both is used"
1378      CALL ipslerr (2,"histdef",str70,str71,str72)
1379      IF (test_fopp < test_fwrt) THEN
1380        freq_opp(pfileid,iv) = pfreq_opp
1381        freq_wrt(pfileid,iv) = pfreq_opp
1382      ELSE
1383        freq_opp(pfileid,iv) = pfreq_wrt
1384        freq_wrt(pfileid,iv) = pfreq_wrt
1385      ENDIF
1386    ENDIF
1387  ELSE IF (INDEX(ex_topps,TRIM(tmp_topp)) > 0) THEN
1388    IF (test_fopp > test_fwrt) THEN
1389      str70 = 'For averages the frequency of operations '// &
1390&             'should be smaller or equal'
1391      WRITE(str71, &
1392 &     '("to that of output. It is not the case for variable ",a)') &
1393 &     TRIM(tmp_name)
1394      str72 = 'PATCH : The output frequency is used for both'
1395      CALL ipslerr (2,"histdef",str70,str71,str72)
1396      freq_opp(pfileid,iv) = pfreq_wrt
1397    ENDIF
1398  ELSE
1399    WRITE (str70,'("Operation on variable ",a," is unknown")') &
1400&    TRIM(tmp_name)
1401    WRITE (str71,'("operation requested is :",a)') tmp_topp
1402    WRITE (str72,'("File ID :",I3)') pfileid
1403    CALL ipslerr (3,"histdef",str70,str71,str72)
1404  ENDIF
1405!-
1406! 5.0 Initialize other variables of the common
1407!-
1408  IF (l_dbg) WRITE(*,*) "histdef : 5.0"
1409!-
1410  hist_wrt_rng(pfileid,iv) = (PRESENT(var_range))
1411  IF (hist_wrt_rng(pfileid,iv)) THEN
1412    hist_calc_rng(pfileid,iv) = (var_range(1) > var_range(2))
1413    IF (hist_calc_rng(pfileid,iv)) THEN
1414      hist_minmax(pfileid,iv,1:2) = &
1415 &      (/ ABS(missing_val),-ABS(missing_val) /)
1416    ELSE
1417      hist_minmax(pfileid,iv,1:2) = var_range(1:2)
1418    ENDIF
1419  ENDIF
1420!-
1421! - freq_opp(pfileid,iv)/2./deltat(pfileid)
1422  last_opp(pfileid,iv) = itau0(pfileid)
1423! - freq_wrt(pfileid,iv)/2./deltat(pfileid)
1424  last_wrt(pfileid,iv) = itau0(pfileid)
1425! - freq_opp(pfileid,iv)/2./deltat(pfileid)
1426  last_opp_chk(pfileid,iv) = itau0(pfileid)
1427! - freq_wrt(pfileid,iv)/2./deltat(pfileid)
1428  last_wrt_chk(pfileid,iv) = itau0(pfileid)
1429  nb_opp(pfileid,iv) = 0
1430  nb_wrt(pfileid,iv) = 0
1431!-
1432! 6.0 Get the time axis for this variable
1433!-
1434  IF (l_dbg) WRITE(*,*) "histdef : 6.0"
1435!-
1436  IF (freq_wrt(pfileid,iv) > 0) THEN
1437    WRITE(str10,'(I8.8)') INT(freq_wrt(pfileid,iv))
1438    str40 = TRIM(tmp_topp)//"_"//TRIM(str10)
1439  ELSE
1440    WRITE(str10,'(I2.2,"month")') ABS(INT(freq_wrt(pfileid,iv)))
1441    str40 = TRIM(tmp_topp)//"_"//TRIM(str10)
1442  ENDIF
1443  CALL find_str (tax_name(pfileid,1:nb_tax(pfileid)),str40,pos)
1444!-
1445! No time axis for once, l_max, l_min or never operation
1446!-
1447  IF (     (TRIM(tmp_topp) /= 'once')  &
1448 &    .AND.(TRIM(tmp_topp) /= 'never') &
1449 &    .AND.(TRIM(tmp_topp) /= 'l_max') &
1450 &    .AND.(TRIM(tmp_topp) /= 'l_min') ) THEN
1451    IF (pos < 0) THEN
1452      nb_tax(pfileid) = nb_tax(pfileid)+1
1453      tax_name(pfileid,nb_tax(pfileid)) = str40
1454      tax_last(pfileid,nb_tax(pfileid)) = 0
1455      var_axid(pfileid,iv) = nb_tax(pfileid)
1456    ELSE
1457      var_axid(pfileid,iv) = pos
1458    ENDIF
1459  ELSE
1460    IF (l_dbg)   WRITE(*,*) "histdef : 7.0 ",TRIM(tmp_topp),'----'
1461    var_axid(pfileid,iv) = -99
1462  ENDIF
1463!-
1464! 7.0 prepare frequence of writing and operation
1465!     for never or once operation
1466!-
1467  IF (    (TRIM(tmp_topp) == 'once')  &
1468 &    .OR.(TRIM(tmp_topp) == 'never') ) THEN
1469    freq_opp(pfileid,iv) = 0.
1470    freq_wrt(pfileid,iv) = 0.
1471  ENDIF
1472!---------------------
1473END SUBROUTINE histdef
1474!===
1475SUBROUTINE histend (pfileid)
1476!---------------------------------------------------------------------
1477!- This subroutine end the decalaration of variables and sets the
1478!- time axes in the netcdf file and puts it into the write mode.
1479!-
1480!- INPUT
1481!-
1482!- pfileid : ID of the file to be worked on
1483!-
1484!- VERSION
1485!-
1486!---------------------------------------------------------------------
1487  IMPLICIT NONE
1488!-
1489  INTEGER,INTENT(IN) :: pfileid
1490!-
1491  INTEGER :: ncid,ncvarid,iret,ndim,iv,itx,ziv,itax,dim_cnt
1492  INTEGER,DIMENSION(4) :: dims
1493  INTEGER :: year,month,day,hours,minutes
1494  REAL :: sec
1495  REAL :: rtime0
1496  CHARACTER(LEN=30) :: str30
1497  CHARACTER(LEN=120) :: assoc
1498  CHARACTER(LEN=70) :: str70
1499  CHARACTER(LEN=3),DIMENSION(12) :: cal =   &
1500 &  (/ 'JAN','FEB','MAR','APR','MAY','JUN', &
1501 &     'JUL','AUG','SEP','OCT','NOV','DEC' /)
1502  CHARACTER(LEN=7) :: tmp_opp
1503  LOGICAL :: l_dbg
1504!---------------------------------------------------------------------
1505  CALL ipsldbg (old_status=l_dbg)
1506!-
1507  ncid = ncdf_ids(pfileid)
1508!-
1509! 1.0 Create the time axes
1510!-
1511  IF (l_dbg) WRITE(*,*) "histend : 1.0"
1512!---
1513  iret = NF90_DEF_DIM (ncid,'time_counter',NF90_UNLIMITED,tid(pfileid))
1514!-
1515! 1.1 Define all the time axes needed for this file
1516!-
1517  DO itx=1,nb_tax(pfileid)
1518    dims(1) = tid(pfileid)
1519    IF (nb_tax(pfileid) > 1) THEN
1520      str30 = "t_"//tax_name(pfileid,itx)
1521    ELSE
1522      str30 = "time_counter"
1523    ENDIF
1524    iret = NF90_DEF_VAR (ncid,str30,NF90_DOUBLE, &
1525 &                       dims(1),tdimid(pfileid,itx))
1526    IF (nb_tax(pfileid) <= 1) THEN
1527      iret = NF90_PUT_ATT (ncid,tdimid(pfileid,itx),'axis',"T")
1528    ENDIF
1529    iret = NF90_PUT_ATT (ncid,tdimid(pfileid,itx),'standard_name',"time")
1530!---
1531!   To transform the current itau into a real date and take it
1532!   as the origin of the file requires the time counter to change.
1533!   Thus it is an operation the user has to ask for.
1534!   This function should thus only be re-instated
1535!   if there is a ioconf routine to control it.
1536!---
1537!-- rtime0 = itau2date(itau0(pfileid),date0(pfileid),deltat(pfileid))
1538    rtime0 = date0(pfileid)
1539!-
1540    CALL ju2ymds(rtime0,year,month,day,sec)
1541!---
1542!   Catch any error induced by a change in calendar !
1543!---
1544    IF (year < 0) THEN
1545      year = 2000+year
1546    ENDIF
1547!-
1548    hours = INT(sec/(60.*60.))
1549    minutes = INT((sec-hours*60.*60.)/60.)
1550    sec = sec-(hours*60.*60.+minutes*60.)
1551!-
1552    WRITE (UNIT=str70, &
1553 &   FMT='(A,I4.4,"-",I2.2,"-",I2.2," ",I2.2,":",I2.2,":",I2.2)') &
1554 &    'seconds since ',year,month,day,hours,minutes,INT(sec)
1555    iret = NF90_PUT_ATT (ncid,tdimid(pfileid,itx),'units',TRIM(str70))
1556!-
1557    CALL ioget_calendar (str30)
1558    iret = NF90_PUT_ATT (ncid,tdimid(pfileid,itx), &
1559 &                       'calendar',TRIM(str30))
1560!-
1561    iret = NF90_PUT_ATT (ncid,tdimid(pfileid,itx),'title','Time')
1562!-
1563    iret = NF90_PUT_ATT (ncid,tdimid(pfileid,itx), &
1564 &                       'long_name','Time axis')
1565!-
1566    WRITE (UNIT=str70, &
1567 &   FMT='(" ",I4.4,"-",A3,"-",I2.2," ",I2.2,":",I2.2,":",I2.2)') &
1568 &    year,cal(month),day,hours,minutes,INT(sec)
1569    iret = NF90_PUT_ATT (ncid,tdimid(pfileid,itx), &
1570 &                       'time_origin',TRIM(str70))
1571  ENDDO
1572!-
1573! 2.0 declare the variables
1574!-
1575  IF (l_dbg) WRITE(*,*) "histend : 2.0"
1576!-
1577  DO iv=1,nb_var(pfileid)
1578!---
1579    itax = var_axid(pfileid,iv)
1580!---
1581    IF (regular(pfileid) ) THEN
1582      dims(1:2) = (/ xid(pfileid),yid(pfileid) /)
1583      dim_cnt = 2
1584    ELSE
1585      dims(1) = xid(pfileid)
1586      dim_cnt = 1
1587    ENDIF
1588!---
1589    tmp_opp = topp(pfileid,iv)
1590    ziv = var_zaxid(pfileid,iv)
1591!---
1592!   2.1 dimension of field
1593!---
1594    IF ((TRIM(tmp_opp) /= 'never')) THEN
1595      IF (     (TRIM(tmp_opp) /= 'once')  &
1596     &    .AND.(TRIM(tmp_opp) /= 'l_max') &
1597     &    .AND.(TRIM(tmp_opp) /= 'l_min') ) THEN
1598        IF (ziv == -99) THEN
1599          ndim = dim_cnt+1
1600          dims(dim_cnt+1:dim_cnt+2) = (/ tid(pfileid),0 /)
1601        ELSE
1602          ndim = dim_cnt+2
1603          dims(dim_cnt+1:dim_cnt+2) = (/zax_ids(pfileid,ziv),tid(pfileid)/)
1604        ENDIF
1605      ELSE
1606        IF (ziv == -99) THEN
1607          ndim = dim_cnt
1608          dims(dim_cnt+1:dim_cnt+2) = (/ 0,0 /)
1609        ELSE
1610          ndim = dim_cnt+1
1611          dims(dim_cnt+1:dim_cnt+2) = (/ zax_ids(pfileid,ziv),0 /)
1612        ENDIF
1613      ENDIF
1614!-
1615      iret = NF90_DEF_VAR (ncid,TRIM(name(pfileid,iv)),NF90_FLOAT, &
1616 &                         dims(1:ABS(ndim)),ncvarid)
1617!-
1618      ncvar_ids(pfileid,iv) = ncvarid
1619!-
1620      IF (LEN_TRIM(unit_name(pfileid,iv)) > 0) THEN
1621        iret = NF90_PUT_ATT (ncid,ncvarid,'units', &
1622 &                           TRIM(unit_name(pfileid,iv)))
1623      ENDIF
1624      iret = NF90_PUT_ATT (ncid,ncvarid,'standard_name', &
1625 &                         TRIM(title(pfileid,iv)))
1626!-
1627      iret = NF90_PUT_ATT (ncid,ncvarid,'_Fillvalue', &
1628 &                         REAL(missing_val,KIND=4))
1629      IF (hist_wrt_rng(pfileid,iv)) THEN
1630        iret = NF90_PUT_ATT (ncid,ncvarid,'valid_min', &
1631 &                           REAL(hist_minmax(pfileid,iv,1),KIND=4))
1632        iret = NF90_PUT_ATT (ncid,ncvarid,'valid_max', &
1633 &                           REAL(hist_minmax(pfileid,iv,2),KIND=4))
1634      ENDIF
1635      iret = NF90_PUT_ATT (ncid,ncvarid,'long_name', &
1636 &                         TRIM(title(pfileid,iv)))
1637      iret = NF90_PUT_ATT (ncid,ncvarid,'online_operation', &
1638 &                         TRIM(fullop(pfileid,iv)))
1639!-
1640      SELECT CASE(ndim)
1641      CASE(-3,2:4)
1642      CASE DEFAULT
1643        CALL ipslerr (3,"histend", &
1644       &  'less than 2 or more than 4 dimensions are not', &
1645       &  'allowed at this stage',' ')
1646      END SELECT
1647!-
1648      assoc=TRIM(hax_name(pfileid,var_haxid(pfileid,iv),2))  &
1649 &   //' '//TRIM(hax_name(pfileid,var_haxid(pfileid,iv),1))
1650!-
1651      ziv = var_zaxid(pfileid,iv)
1652      IF (ziv > 0) THEN
1653        str30 = zax_name(pfileid,ziv)
1654        assoc = TRIM(str30)//' '//TRIM(assoc)
1655      ENDIF
1656!-
1657      IF (itax > 0) THEN
1658        IF (nb_tax(pfileid) > 1) THEN
1659          str30 = "t_"//tax_name(pfileid,itax)
1660        ELSE
1661          str30 = "time_counter"
1662        ENDIF
1663        assoc = TRIM(str30)//' '//TRIM(assoc)
1664!-
1665        IF (l_dbg) THEN
1666          WRITE(*,*) "histend : 2.0.n, freq_opp, freq_wrt", &
1667 &                   freq_opp(pfileid,iv),freq_wrt(pfileid,iv)
1668        ENDIF
1669!-
1670        iret = NF90_PUT_ATT (ncid,ncvarid,'interval_operation', &
1671 &                           REAL(freq_opp(pfileid,iv),KIND=4))
1672        iret = NF90_PUT_ATT (ncid,ncvarid,'interval_write', &
1673 &                           REAL(freq_wrt(pfileid,iv),KIND=4))
1674      ENDIF
1675      iret = NF90_PUT_ATT (ncid,ncvarid,'coordinates',TRIM(assoc))
1676    ENDIF
1677  ENDDO
1678!-
1679! 2.2 Add DOMAIN attributes if needed
1680!-
1681  IF (dom_id_svg(pfileid) >= 0) THEN
1682    CALL flio_dom_att (ncid,dom_id_svg(pfileid))
1683  ENDIF
1684!-
1685! 3.0 Put the netcdf file into write mode
1686!-
1687  IF (l_dbg) WRITE(*,*) "histend : 3.0"
1688!-
1689  iret = NF90_ENDDEF (ncid)
1690!-
1691! 4.0 Give some informations to the user
1692!-
1693  IF (l_dbg) WRITE(*,*) "histend : 4.0"
1694!-
1695  WRITE(str70,'("All variables have been initialized on file :",I3)') pfileid
1696  CALL ipslerr (1,'histend',str70,'',' ')
1697!---------------------
1698END SUBROUTINE histend
1699!===
1700SUBROUTINE histwrite_r1d (pfileid,pvarname,pitau,pdata,nbindex,nindex)
1701!---------------------------------------------------------------------
1702  IMPLICIT NONE
1703!-
1704  INTEGER,INTENT(IN) :: pfileid,pitau,nbindex
1705  INTEGER,DIMENSION(nbindex),INTENT(IN) :: nindex
1706  REAL,DIMENSION(:),INTENT(IN) :: pdata
1707  CHARACTER(LEN=*),INTENT(IN) :: pvarname
1708!---------------------------------------------------------------------
1709  CALL histw_rnd (pfileid,pvarname,pitau,nbindex,nindex,pdata_1d=pdata)
1710!---------------------------
1711END SUBROUTINE histwrite_r1d
1712!===
1713SUBROUTINE histwrite_r2d (pfileid,pvarname,pitau,pdata,nbindex,nindex)
1714!---------------------------------------------------------------------
1715  IMPLICIT NONE
1716!-
1717  INTEGER,INTENT(IN) :: pfileid,pitau,nbindex
1718  INTEGER,DIMENSION(nbindex),INTENT(IN) :: nindex
1719  REAL,DIMENSION(:,:),INTENT(IN) :: pdata
1720  CHARACTER(LEN=*),INTENT(IN) :: pvarname
1721!---------------------------------------------------------------------
1722  CALL histw_rnd (pfileid,pvarname,pitau,nbindex,nindex,pdata_2d=pdata)
1723!---------------------------
1724END SUBROUTINE histwrite_r2d
1725!===
1726SUBROUTINE histwrite_r3d (pfileid,pvarname,pitau,pdata,nbindex,nindex)
1727!---------------------------------------------------------------------
1728  IMPLICIT NONE
1729!-
1730  INTEGER,INTENT(IN) :: pfileid,pitau,nbindex
1731  INTEGER,DIMENSION(nbindex),INTENT(IN) :: nindex
1732  REAL,DIMENSION(:,:,:),INTENT(IN) :: pdata
1733  CHARACTER(LEN=*),INTENT(IN) :: pvarname
1734!---------------------------------------------------------------------
1735  CALL histw_rnd (pfileid,pvarname,pitau,nbindex,nindex,pdata_3d=pdata)
1736!---------------------------
1737END SUBROUTINE histwrite_r3d
1738!===
1739SUBROUTINE histw_rnd (pfileid,pvarname,pitau,nbindex,nindex, &
1740  &                   pdata_1d,pdata_2d,pdata_3d)
1741!---------------------------------------------------------------------
1742  IMPLICIT NONE
1743!-
1744  INTEGER,INTENT(IN) :: pfileid,pitau,nbindex
1745  INTEGER,DIMENSION(nbindex),INTENT(IN) :: nindex
1746  CHARACTER(LEN=*),INTENT(IN) :: pvarname
1747  REAL,DIMENSION(:),INTENT(IN),OPTIONAL     :: pdata_1d
1748  REAL,DIMENSION(:,:),INTENT(IN),OPTIONAL   :: pdata_2d
1749  REAL,DIMENSION(:,:,:),INTENT(IN),OPTIONAL :: pdata_3d
1750!-
1751  LOGICAL :: do_oper,do_write,largebuf,l1d,l2d,l3d
1752  INTEGER :: varid,io,nbpt_out
1753  INTEGER              :: nbpt_in1
1754  INTEGER,DIMENSION(2) :: nbpt_in2
1755  INTEGER,DIMENSION(3) :: nbpt_in3
1756  REAL,ALLOCATABLE,DIMENSION(:),SAVE :: buff_tmp
1757  INTEGER,SAVE :: buff_tmp_sz
1758  CHARACTER(LEN=7) :: tmp_opp
1759  CHARACTER(LEN=13) :: c_nam
1760  LOGICAL :: l_dbg
1761!---------------------------------------------------------------------
1762  CALL ipsldbg (old_status=l_dbg)
1763!-
1764  l1d=PRESENT(pdata_1d); l2d=PRESENT(pdata_2d); l3d=PRESENT(pdata_3d);
1765  IF      (l1d) THEN
1766    c_nam = 'histwrite_r1d'
1767  ELSE IF (l2d) THEN
1768    c_nam = 'histwrite_r2d'
1769  ELSE IF (l3d) THEN
1770    c_nam = 'histwrite_r3d'
1771  ENDIF
1772!-
1773! 1.0 Try to catch errors like specifying the wrong file ID.
1774!     Thanks Marine for showing us what errors users can make !
1775!-
1776  IF ( (pfileid < 1).OR.(pfileid > nb_files) ) THEN
1777    CALL ipslerr (3,"histwrite", &
1778 &    'Illegal file ID in the histwrite of variable',pvarname,' ')
1779  ENDIF
1780!-
1781! 1.1 Find the id of the variable to be written and the real time
1782!-
1783  CALL histvar_seq (pfileid,pvarname,varid)
1784!-
1785! 2.0 do nothing for never operation
1786!-
1787  tmp_opp = topp(pfileid,varid)
1788!-
1789  IF (TRIM(tmp_opp) == "never") THEN
1790    last_opp_chk(pfileid,varid) = -99
1791    last_wrt_chk(pfileid,varid) = -99
1792  ENDIF
1793!-
1794! 3.0 We check if we need to do an operation
1795!-
1796  IF (last_opp_chk(pfileid,varid) == pitau) THEN
1797    CALL ipslerr (3,"histwrite", &
1798 &    'This variable has already been analysed at the present', &
1799 &    'time step',TRIM(pvarname))
1800  ENDIF
1801!-
1802  CALL isittime &
1803 &  (pitau,date0(pfileid),deltat(pfileid),freq_opp(pfileid,varid), &
1804 &   last_opp(pfileid,varid),last_opp_chk(pfileid,varid),do_oper)
1805!-
1806! 4.0 We check if we need to write the data
1807!-
1808  IF (last_wrt_chk(pfileid,varid) == pitau) THEN
1809    CALL ipslerr (3,"histwrite", &
1810 &    'This variable has already been written for the present', &
1811 &    'time step',TRIM(pvarname))
1812  ENDIF
1813!-
1814  CALL isittime &
1815 &  (pitau,date0(pfileid),deltat(pfileid),freq_wrt(pfileid,varid), &
1816 &   last_wrt(pfileid,varid),last_wrt_chk(pfileid,varid),do_write)
1817!-
1818! 5.0 histwrite called
1819!-
1820  IF (do_oper.OR.do_write) THEN
1821!-
1822!-- 5.1 Get the sizes of the data we will handle
1823!-
1824    IF (datasz_in(pfileid,varid,1) <= 0) THEN
1825!---- There is the risk here that the user has over-sized the array.
1826!---- But how can we catch this ?
1827!---- In the worst case we will do impossible operations
1828!---- on part of the data !
1829      datasz_in(pfileid,varid,1:3) = -1
1830      IF      (l1d) THEN
1831        datasz_in(pfileid,varid,1) = SIZE(pdata_1d)
1832      ELSE IF (l2d) THEN
1833        datasz_in(pfileid,varid,1) = SIZE(pdata_2d,DIM=1)
1834        datasz_in(pfileid,varid,2) = SIZE(pdata_2d,DIM=2)
1835      ELSE IF (l3d) THEN
1836        datasz_in(pfileid,varid,1) = SIZE(pdata_3d,DIM=1)
1837        datasz_in(pfileid,varid,2) = SIZE(pdata_3d,DIM=2)
1838        datasz_in(pfileid,varid,3) = SIZE(pdata_3d,DIM=3)
1839      ENDIF
1840    ENDIF
1841!-
1842!-- 5.2 The maximum size of the data will give the size of the buffer
1843!-
1844    IF (datasz_max(pfileid,varid) <= 0) THEN
1845      largebuf = .FALSE.
1846      DO io=1,nbopp(pfileid,varid)
1847        IF (INDEX(fuchnbout,sopps(pfileid,varid,io)) > 0) THEN
1848          largebuf = .TRUE.
1849        ENDIF
1850      ENDDO
1851      IF (largebuf) THEN
1852        datasz_max(pfileid,varid) = &
1853 &        scsize(pfileid,varid,1) &
1854 &       *scsize(pfileid,varid,2) &
1855 &       *scsize(pfileid,varid,3)
1856      ELSE
1857        IF      (l1d) THEN
1858          datasz_max(pfileid,varid) = &
1859 &          datasz_in(pfileid,varid,1)
1860        ELSE IF (l2d) THEN
1861          datasz_max(pfileid,varid) = &
1862 &          datasz_in(pfileid,varid,1) &
1863 &         *datasz_in(pfileid,varid,2)
1864        ELSE IF (l3d) THEN
1865          datasz_max(pfileid,varid) = &
1866 &          datasz_in(pfileid,varid,1) &
1867 &         *datasz_in(pfileid,varid,2) &
1868 &         *datasz_in(pfileid,varid,3)
1869        ENDIF
1870      ENDIF
1871    ENDIF
1872!-
1873    IF (.NOT.ALLOCATED(buff_tmp)) THEN
1874      IF (l_dbg) THEN
1875        WRITE(*,*) &
1876 &       c_nam//" : allocate buff_tmp for buff_sz = ", &
1877 &       datasz_max(pfileid,varid)
1878      ENDIF
1879      ALLOCATE(buff_tmp(datasz_max(pfileid,varid)))
1880      buff_tmp_sz = datasz_max(pfileid,varid)
1881    ELSE IF (datasz_max(pfileid,varid) > buff_tmp_sz) THEN
1882      IF (l_dbg) THEN
1883        WRITE(*,*) &
1884 &       c_nam//" : re-allocate buff_tmp for buff_sz = ", &
1885 &       datasz_max(pfileid,varid)
1886      ENDIF
1887      DEALLOCATE(buff_tmp)
1888      ALLOCATE(buff_tmp(datasz_max(pfileid,varid)))
1889      buff_tmp_sz = datasz_max(pfileid,varid)
1890    ENDIF
1891!-
1892!-- We have to do the first operation anyway.
1893!-- Thus we do it here and change the ranke
1894!-- of the data at the same time. This should speed up things.
1895!-
1896    nbpt_out = datasz_max(pfileid,varid)
1897    IF      (l1d) THEN
1898      nbpt_in1 = datasz_in(pfileid,varid,1)
1899      CALL mathop (sopps(pfileid,varid,1),nbpt_in1,pdata_1d, &
1900 &                 missing_val,nbindex,nindex, &
1901 &                 scal(pfileid,varid,1),nbpt_out,buff_tmp)
1902    ELSE IF (l2d) THEN
1903      nbpt_in2(1:2) = datasz_in(pfileid,varid,1:2)
1904      CALL mathop (sopps(pfileid,varid,1),nbpt_in2,pdata_2d, &
1905 &                 missing_val,nbindex,nindex, &
1906 &                 scal(pfileid,varid,1),nbpt_out,buff_tmp)
1907    ELSE IF (l3d) THEN
1908      nbpt_in3(1:3) = datasz_in(pfileid,varid,1:3)
1909      CALL mathop (sopps(pfileid,varid,1),nbpt_in3,pdata_3d, &
1910 &                 missing_val,nbindex,nindex, &
1911 &                 scal(pfileid,varid,1),nbpt_out,buff_tmp)
1912    ENDIF
1913    CALL histwrite_real (pfileid,varid,pitau,nbpt_out, &
1914 &            buff_tmp,nbindex,nindex,do_oper,do_write)
1915  ENDIF
1916!-
1917! 6.0 Manage time steps
1918!-
1919  IF ((TRIM(tmp_opp) /= "once").AND.(TRIM(tmp_opp) /= "never")) THEN
1920    last_opp_chk(pfileid,varid) = pitau
1921    last_wrt_chk(pfileid,varid) = pitau
1922  ELSE
1923    last_opp_chk(pfileid,varid) = -99
1924    last_wrt_chk(pfileid,varid) = -99
1925  ENDIF
1926!-----------------------
1927END SUBROUTINE histw_rnd
1928!===
1929SUBROUTINE histwrite_real &
1930 & (pfileid,varid,pitau,nbdpt,buff_tmp,nbindex,nindex,do_oper,do_write)
1931!---------------------------------------------------------------------
1932!- This subroutine is internal and does the calculations and writing
1933!- if needed. At a later stage it should be split into an operation
1934!- and writing subroutines.
1935!---------------------------------------------------------------------
1936  IMPLICIT NONE
1937!-
1938  INTEGER,INTENT(IN) :: pfileid,pitau,varid, &
1939 &                      nbindex,nindex(nbindex),nbdpt
1940  REAL,DIMENSION(:)  :: buff_tmp
1941  LOGICAL,INTENT(IN) :: do_oper,do_write
1942!-
1943  INTEGER :: tsz,ncid,ncvarid,i,iret,ipt,itax,io,nbin,nbout
1944  INTEGER :: nx,ny,nz,ky,kz,kt,kc
1945  INTEGER,DIMENSION(4) :: corner,edges
1946  INTEGER :: itime
1947!-
1948  REAL :: rtime
1949  CHARACTER(LEN=7) :: tmp_opp
1950  REAL,ALLOCATABLE,DIMENSION(:),SAVE :: buff_tmp2,buffer_used
1951  INTEGER,SAVE :: buff_tmp2_sz,buffer_sz
1952  LOGICAL :: l_dbg
1953!---------------------------------------------------------------------
1954  CALL ipsldbg (old_status=l_dbg)
1955!-
1956  IF (l_dbg) THEN
1957    WRITE(*,*) "histwrite 0.0 : VAR : ",name(pfileid,varid)
1958    WRITE(*,*) "histwrite 0.0 : nbindex,nindex :", &
1959 &  nbindex,nindex(1:MIN(3,nbindex)),'...',nindex(MAX(1,nbindex-3):nbindex)
1960  ENDIF
1961!-
1962! The sizes which can be encoutered
1963!-
1964  tsz =  zsize(pfileid,varid,1) &
1965 &      *zsize(pfileid,varid,2) &
1966 &      *zsize(pfileid,varid,3)
1967!-
1968! 1.0 We allocate the memory needed to store the data between write
1969!     and the temporary space needed for operations.
1970!     We have to keep precedent buffer if needed
1971!-
1972  IF (.NOT. ALLOCATED(buffer)) THEN
1973    IF (l_dbg) WRITE(*,*) "histwrite_real 1.0 allocate buffer ",buff_pos
1974    ALLOCATE(buffer(buff_pos))
1975    buffer_sz = buff_pos
1976    buffer(:)=0.0
1977  ELSE IF (buffer_sz < buff_pos) THEN
1978    IF (l_dbg) THEN
1979      WRITE(*,*) "histwrite_real 1.0.1 re-allocate buffer for ", &
1980 &      buff_pos," instead of ",SIZE(buffer)
1981    ENDIF
1982    IF (SUM(buffer)/=0.0) THEN
1983      IF (l_dbg) WRITE (*,*) 'histwrite : buffer has been used. ', &
1984 &      'We have to save it before re-allocating.'
1985      ALLOCATE(buffer_used(buffer_sz))
1986      buffer_used(:)=buffer(:)
1987      DEALLOCATE(buffer)
1988      ALLOCATE(buffer(buff_pos))
1989      buffer_sz = buff_pos
1990      buffer(:)=0.0
1991      buffer(:SIZE(buffer_used))=buffer_used
1992      DEALLOCATE(buffer_used)
1993    ELSE
1994      IF (l_dbg) WRITE (*,*) 'histwrite : buffer has not been used. ', &
1995 &      'We have just to re-allocate it.'
1996      DEALLOCATE(buffer)
1997      ALLOCATE(buffer(buff_pos))
1998      buffer_sz = buff_pos
1999      buffer(:)=0.0
2000    ENDIF
2001  ENDIF
2002!-
2003! The buffers are only deallocated when more space is needed. This
2004! reduces the umber of allocates but increases memory needs.
2005!-
2006  IF (.NOT.ALLOCATED(buff_tmp2)) THEN
2007    IF (l_dbg) THEN
2008      WRITE(*,*) "histwrite_real 1.1 allocate buff_tmp2 ",SIZE(buff_tmp)
2009    ENDIF
2010    ALLOCATE(buff_tmp2(datasz_max(pfileid,varid)))
2011    buff_tmp2_sz = datasz_max(pfileid,varid)
2012  ELSE IF (datasz_max(pfileid,varid) > buff_tmp2_sz) THEN
2013    IF (l_dbg) THEN
2014      WRITE(*,*) "histwrite_real 1.2 re-allocate buff_tmp2 : ", &
2015     & SIZE(buff_tmp)," instead of ",SIZE(buff_tmp2)
2016    ENDIF
2017    DEALLOCATE(buff_tmp2)
2018    ALLOCATE(buff_tmp2(datasz_max(pfileid,varid)))
2019    buff_tmp2_sz = datasz_max(pfileid,varid)
2020  ENDIF
2021!-
2022  rtime = pitau * deltat(pfileid)
2023  tmp_opp = topp(pfileid,varid)
2024!-
2025! 3.0 Do the operations or transfer the slab of data into buff_tmp
2026!-
2027  IF (l_dbg) WRITE(*,*) "histwrite: 3.0",pfileid
2028!-
2029! 3.1 DO the Operations only if needed
2030!-
2031  IF (do_oper) THEN
2032    i = pfileid
2033    nbout = nbdpt
2034!-
2035!-- 3.4 We continue the sequence of operations
2036!--     we started in the interface routine
2037!-
2038    DO io = 2,nbopp(i,varid),2
2039      nbin = nbout
2040      nbout = datasz_max(i,varid)
2041      CALL mathop(sopps(i,varid,io),nbin,buff_tmp,missing_val, &
2042 &      nbindex,nindex,scal(i,varid,io),nbout,buff_tmp2)
2043      IF (l_dbg) THEN
2044        WRITE(*,*) &
2045 &       "histwrite: 3.4a nbout : ",nbin,nbout,sopps(i,varid,io)
2046      ENDIF
2047!-
2048      nbin = nbout
2049      nbout = datasz_max(i,varid)
2050      CALL mathop(sopps(i,varid,io+1),nbin,buff_tmp2,missing_val, &
2051 &      nbindex,nindex,scal(i,varid,io+1),nbout,buff_tmp)
2052      IF (l_dbg) THEN
2053        WRITE(*,*) &
2054 &       "histwrite: 3.4b nbout : ",nbin,nbout,sopps(i,varid,io+1)
2055      ENDIF
2056    ENDDO
2057!-
2058!   3.5 Zoom into the data
2059!-
2060    IF (l_dbg) THEN
2061      WRITE(*,*) &
2062 &     "histwrite: 3.5 size(buff_tmp) : ",SIZE(buff_tmp)
2063      WRITE(*,*) &
2064 &     "histwrite: 3.5 slab in X :",zorig(i,varid,1),zsize(i,varid,1)
2065      WRITE(*,*) &
2066 &     "histwrite: 3.5 slab in Y :",zorig(i,varid,2),zsize(i,varid,2)
2067      WRITE(*,*) &
2068 &     "histwrite: 3.5 slab in Z :",zorig(i,varid,3),zsize(i,varid,3)
2069      WRITE(*,*) &
2070 &     "histwrite: 3.5 slab of input:", &
2071 &     scsize(i,varid,1),scsize(i,varid,2),scsize(i,varid,3)
2072    ENDIF
2073!---
2074!-- We have to consider blocks of contiguous data
2075!---
2076    nx=MAX(zsize(i,varid,1),1)
2077    ny=MAX(zsize(i,varid,2),1)
2078    nz=MAX(zsize(i,varid,3),1)
2079    IF     (     (zorig(i,varid,1) == 1) &
2080   &        .AND.(zsize(i,varid,1) == scsize(i,varid,1)) &
2081   &        .AND.(zorig(i,varid,2) == 1) &
2082   &        .AND.(zsize(i,varid,2) == scsize(i,varid,2))) THEN
2083      kt = (zorig(i,varid,3)-1)*nx*ny
2084      buff_tmp2(1:nx*ny*nz) = buff_tmp(kt+1:kt+nx*ny*nz)
2085    ELSEIF (     (zorig(i,varid,1) == 1) &
2086   &        .AND.(zsize(i,varid,1) == scsize(i,varid,1))) THEN
2087      kc = -nx*ny
2088      DO kz=zorig(i,varid,3),zorig(i,varid,3)+nz-1
2089        kc = kc+nx*ny
2090        kt = ((kz-1)*scsize(i,varid,2)+zorig(i,varid,2)-1)*nx
2091        buff_tmp2(kc+1:kc+nx*ny) = buff_tmp(kt+1:kt+nx*ny)
2092      ENDDO
2093    ELSE
2094      kc = -nx
2095      DO kz=zorig(i,varid,3),zorig(i,varid,3)+nz-1
2096        DO ky=zorig(i,varid,2),zorig(i,varid,2)+ny-1
2097          kc = kc+nx
2098          kt = ((kz-1)*scsize(i,varid,2)+ky-1)*scsize(i,varid,1) &
2099   &          +zorig(i,varid,1)-1
2100          buff_tmp2(kc+1:kc+nx) = buff_tmp(kt+1:kt+nx)
2101        ENDDO
2102      ENDDO
2103    ENDIF
2104!-
2105!-- 4.0 Get the min and max of the field (buff_tmp)
2106!-
2107    IF (l_dbg) WRITE(*,*) "histwrite: 4.0 buff_tmp",pfileid,varid, &
2108 &    TRIM(tmp_opp),' ---- ',LEN_TRIM(tmp_opp),nbindex
2109!-
2110    IF (hist_calc_rng(pfileid,varid)) THEN
2111      hist_minmax(pfileid,varid,1) = &
2112 &      MIN(hist_minmax(pfileid,varid,1), &
2113 &      MINVAL(buff_tmp2(1:tsz),MASK=buff_tmp2(1:tsz) /= missing_val))
2114      hist_minmax(pfileid,varid,2) = &
2115 &      MAX(hist_minmax(pfileid,varid,2), &
2116 &      MAXVAL(buff_tmp2(1:tsz),MASK=buff_tmp2(1:tsz) /= missing_val))
2117    ENDIF
2118!-
2119!-- 5.0 Do the operations if needed. In the case of instantaneous
2120!--     output we do not transfer to the buffer.
2121!-
2122    IF (l_dbg) WRITE(*,*) "histwrite: 5.0",pfileid,"tsz :",tsz
2123!-
2124    ipt = point(pfileid,varid)
2125!-
2126!   WRITE(*,*) 'OPE ipt, buffer :',pvarname,ipt,varid
2127!-
2128    IF (     (TRIM(tmp_opp) /= "inst") &
2129   &    .AND.(TRIM(tmp_opp) /= "once") ) THEN
2130      CALL moycum(tmp_opp,tsz,buffer(ipt:), &
2131     &       buff_tmp2,nb_opp(pfileid,varid))
2132    ENDIF
2133!-
2134    last_opp(pfileid,varid) = pitau
2135    nb_opp(pfileid,varid) = nb_opp(pfileid,varid)+1
2136!-
2137  ENDIF
2138!-
2139! 6.0 Write to file if needed
2140!-
2141  IF (l_dbg) WRITE(*,*) "histwrite: 6.0",pfileid
2142!-
2143  IF (do_write) THEN
2144!-
2145    ncvarid = ncvar_ids(pfileid,varid)
2146    ncid = ncdf_ids(pfileid)
2147!-
2148!-- 6.1 Do the operations that are needed before writting
2149!-
2150    IF (l_dbg) WRITE(*,*) "histwrite: 6.1",pfileid
2151!-
2152    IF (     (TRIM(tmp_opp) /= "inst") &
2153   &    .AND.(TRIM(tmp_opp) /= "once") ) THEN
2154      rtime = (rtime+last_wrt(pfileid,varid)*deltat(pfileid))/2.0
2155    ENDIF
2156!-
2157!-- 6.2 Add a value to the time axis of this variable if needed
2158!-
2159    IF (     (TRIM(tmp_opp) /= "l_max") &
2160   &    .AND.(TRIM(tmp_opp) /= "l_min") &
2161   &    .AND.(TRIM(tmp_opp) /= "once") ) THEN
2162!-
2163      IF (l_dbg) WRITE(*,*) "histwrite: 6.2",pfileid
2164!-
2165      itax = var_axid(pfileid,varid)
2166      itime = nb_wrt(pfileid,varid)+1
2167!-
2168      IF (tax_last(pfileid,itax) < itime) THEN
2169        iret = NF90_PUT_VAR (ncid,tdimid(pfileid,itax),(/ rtime /), &
2170 &                            start=(/ itime /),count=(/ 1 /))
2171        tax_last(pfileid,itax) = itime
2172      ENDIF
2173    ELSE
2174      itime=1
2175    ENDIF
2176!-
2177!-- 6.3 Write the data. Only in the case of instantaneous output
2178!       we do not write the buffer.
2179!-
2180    IF (l_dbg) THEN
2181      WRITE(*,*) "histwrite: 6.3",pfileid,ncid,ncvarid,varid,itime
2182    ENDIF
2183!-
2184    IF (scsize(pfileid,varid,3) == 1) THEN
2185      IF (regular(pfileid)) THEN
2186        corner(1:4) = (/ 1,1,itime,0 /)
2187        edges(1:4) = (/ zsize(pfileid,varid,1), &
2188 &                      zsize(pfileid,varid,2),1,0 /)
2189      ELSE
2190        corner(1:4) = (/ 1,itime,0,0 /)
2191        edges(1:4) = (/ zsize(pfileid,varid,1),1,0,0 /)
2192      ENDIF
2193    ELSE
2194      IF (regular(pfileid)) THEN
2195        corner(1:4) = (/ 1,1,1,itime /)
2196        edges(1:4) = (/ zsize(pfileid,varid,1), &
2197 &                      zsize(pfileid,varid,2), &
2198 &                      zsize(pfileid,varid,3),1 /)
2199      ELSE
2200        corner(1:4) = (/ 1,1,itime,0 /)
2201        edges(1:4) = (/ zsize(pfileid,varid,1), &
2202 &                      zsize(pfileid,varid,3),1,0 /)
2203      ENDIF
2204    ENDIF
2205!-
2206    ipt = point(pfileid,varid)
2207!-
2208    IF (     (TRIM(tmp_opp) /= "inst") &
2209 &      .AND.(TRIM(tmp_opp) /= "once") ) THEN
2210      iret = NF90_PUT_VAR (ncid,ncvarid,buffer(ipt:), &
2211 &                       start=corner(1:4),count=edges(1:4))
2212    ELSE
2213      iret = NF90_PUT_VAR (ncid,ncvarid,buff_tmp2, &
2214 &                       start=corner(1:4),count=edges(1:4))
2215    ENDIF
2216!-
2217    last_wrt(pfileid,varid) = pitau
2218    nb_wrt(pfileid,varid) = nb_wrt(pfileid,varid)+1
2219    nb_opp(pfileid,varid) = 0
2220!---
2221!   After the write the file can be synchronized so that no data is
2222!   lost in case of a crash. This feature gives up on the benefits of
2223!   buffering and should only be used in debuging mode. A flag is
2224!   needed here to switch to this mode.
2225!---
2226!   iret = NF90_SYNC (ncid)
2227!-
2228  ENDIF
2229!----------------------------
2230END SUBROUTINE histwrite_real
2231!===
2232SUBROUTINE histvar_seq (pfid,pvarname,pvid)
2233!---------------------------------------------------------------------
2234!- This subroutine optimized the search for the variable in the table.
2235!- In a first phase it will learn the succession of the variables
2236!- called and then it will use the table to guess what comes next.
2237!- It is the best solution to avoid lengthy searches through array
2238!- vectors.
2239!-
2240!- ARGUMENTS :
2241!-
2242!- pfid  : id of the file on which we work
2243!- pvarname : The name of the variable we are looking for
2244!- pvid     : The var id we found
2245!---------------------------------------------------------------------
2246  IMPLICIT NONE
2247!-
2248  INTEGER,INTENT(in)  :: pfid
2249  CHARACTER(LEN=*),INTENT(IN) :: pvarname
2250  INTEGER,INTENT(out) :: pvid
2251!-
2252  LOGICAL,SAVE :: learning(nb_files_max)=.TRUE.
2253  INTEGER,SAVE :: overlap(nb_files_max) = -1
2254  INTEGER,SAVE :: varseq(nb_files_max,nb_var_max*3)
2255  INTEGER,SAVE :: varseq_len(nb_files_max) = 0
2256  INTEGER,SAVE :: varseq_pos(nb_files_max)
2257  INTEGER,SAVE :: varseq_err(nb_files_max) = 0
2258  INTEGER      :: ib,sp,nn,pos
2259  CHARACTER(LEN=70) :: str70
2260  LOGICAL :: l_dbg
2261!---------------------------------------------------------------------
2262  CALL ipsldbg (old_status=l_dbg)
2263!-
2264  IF (l_dbg) THEN
2265    WRITE(*,*) 'histvar_seq, start of the subroutine :',learning(pfid)
2266  ENDIF
2267!-
2268  IF (learning(pfid)) THEN
2269!-
2270!-- 1.0 We compute the length over which we are going
2271!--     to check the overlap
2272!-
2273    IF (overlap(pfid) <= 0) THEN
2274      IF (nb_var(pfid) > 6) THEN
2275        overlap(pfid) = nb_var(pfid)/3*2
2276      ELSE
2277        overlap(pfid) = nb_var(pfid)
2278      ENDIF
2279    ENDIF
2280!-
2281!-- 1.1 Find the position of this string
2282!-
2283    CALL find_str (name(pfid,1:nb_var(pfid)),pvarname,pos)
2284    IF (pos > 0) THEN
2285      pvid = pos
2286    ELSE
2287      CALL ipslerr (3,"histvar_seq", &
2288 &      'The name of the variable you gave has not been declared', &
2289 &      'You should use subroutine histdef for declaring variable', &
2290 &      TRIM(pvarname))
2291    ENDIF
2292!-
2293!-- 1.2 If we have not given up we store the position
2294!--     in the sequence of calls
2295!-
2296    IF (varseq_err(pfid) >= 0) THEN
2297      sp = varseq_len(pfid)+1
2298      IF (sp <= nb_var_max*3) THEN
2299        varseq(pfid,sp) = pvid
2300        varseq_len(pfid) = sp
2301      ELSE
2302        CALL ipslerr (2,"histvar_seq",&
2303 &       'The learning process has failed and we give up. '// &
2304 &       'Either you sequence is',&
2305 &       'too complex or I am too dumb. '// &
2306 &       'This will only affect the efficiency',&
2307 &       'of your code. Thus if you wish to save time'// &
2308 &       ' contact the IOIPSL team. ')
2309        WRITE(*,*) 'The sequence we have found up to now :'
2310        WRITE(*,*) varseq(pfid,1:sp-1)
2311        varseq_err(pfid) = -1
2312      ENDIF
2313!-
2314!---- 1.3 Check if we have found the right overlap
2315!-
2316      IF (varseq_len(pfid) .GE. overlap(pfid)*2) THEN
2317!-
2318!------ We skip a few variables if needed as they could come
2319!------ from the initialisation of the model.
2320!-
2321        DO ib = 0,sp-overlap(pfid)*2
2322          IF ( learning(pfid) .AND.&
2323            & SUM(ABS(varseq(pfid,ib+1:ib+overlap(pfid)) -&
2324            & varseq(pfid,sp-overlap(pfid)+1:sp))) == 0 ) THEN
2325            learning(pfid) = .FALSE.
2326            varseq_len(pfid) = sp-overlap(pfid)-ib
2327            varseq_pos(pfid) = overlap(pfid)+ib
2328            varseq(pfid,1:varseq_len(pfid)) = &
2329 &            varseq(pfid,ib+1:ib+varseq_len(pfid))
2330          ENDIF
2331        ENDDO
2332      ENDIF
2333    ENDIF
2334  ELSE
2335!-
2336!-- 2.0 Now we know how the calls to histwrite are sequenced
2337!--     and we can get a guess at the var ID
2338!-
2339    nn = varseq_pos(pfid)+1
2340    IF (nn > varseq_len(pfid)) nn = 1
2341!-
2342    pvid = varseq(pfid,nn)
2343!-
2344    IF (TRIM(name(pfid,pvid)) /= TRIM(pvarname)) THEN
2345      CALL find_str (name(pfid,1:nb_var(pfid)),pvarname,pos)
2346      IF (pos > 0) THEN
2347        pvid = pos
2348      ELSE
2349        CALL ipslerr (3,"histvar_seq", &
2350 &  'The name of the variable you gave has not been declared',&
2351 &  'You should use subroutine histdef for declaring variable', &
2352 &  TRIM(pvarname))
2353      ENDIF
2354      varseq_err(pfid) = varseq_err(pfid)+1
2355    ELSE
2356!-
2357!---- We only keep the new position if we have found the variable
2358!---- this way. This way an out of sequence call to histwrite does
2359!---- not defeat the process.
2360!-
2361      varseq_pos(pfid) = nn
2362    ENDIF
2363!-
2364    IF (varseq_err(pfid) .GE. 10) THEN
2365      WRITE(str70,'("for file ",I3)') pfid
2366      CALL ipslerr (2,"histvar_seq", &
2367 &  'There were 10 errors in the learned sequence of variables',&
2368 &  str70,'This looks like a bug, please report it.')
2369         varseq_err(pfid) = 0
2370    ENDIF
2371  ENDIF
2372!-
2373  IF (l_dbg) THEN
2374    WRITE(*,*) &
2375 &   'histvar_seq, end of the subroutine :',TRIM(pvarname),pvid
2376  ENDIF
2377!-------------------------
2378END SUBROUTINE histvar_seq
2379!===
2380SUBROUTINE histsync (file)
2381!---------------------------------------------------------------------
2382!- This subroutine will synchronise all
2383!- (or one if defined) opened files.
2384!-
2385!- VERSION
2386!-
2387!---------------------------------------------------------------------
2388  IMPLICIT NONE
2389!-
2390! file  : optional argument for fileid
2391  INTEGER,INTENT(in),OPTIONAL :: file
2392!-
2393  INTEGER :: ifile,ncid,iret
2394!-
2395  LOGICAL :: file_exists
2396  LOGICAL :: l_dbg
2397!---------------------------------------------------------------------
2398  CALL ipsldbg (old_status=l_dbg)
2399!-
2400  IF (l_dbg) WRITE(*,*) 'Entering loop on files : ',nb_files
2401!-
2402! 1.The loop on files to synchronise
2403!-
2404  DO ifile = 1,nb_files
2405!-
2406    IF (PRESENT(file)) THEN
2407      file_exists = (ifile == file)
2408    ELSE
2409      file_exists = .TRUE.
2410    ENDIF
2411!-
2412    IF (file_exists) THEN
2413      IF (l_dbg) THEN
2414        WRITE(*,*) 'Synchronising specified file number :',file
2415      ENDIF
2416      ncid = ncdf_ids(ifile)
2417      iret = NF90_SYNC (ncid)
2418    ENDIF
2419!-
2420  ENDDO
2421!----------------------
2422END SUBROUTINE histsync
2423!===
2424SUBROUTINE histclo (fid)
2425!---------------------------------------------------------------------
2426!- This subroutine will close all (or one if defined) opened files
2427!-
2428!- VERSION
2429!-
2430!---------------------------------------------------------------------
2431  IMPLICIT NONE
2432!-
2433! fid  : optional argument for fileid
2434  INTEGER,INTENT(in),OPTIONAL :: fid
2435!-
2436  INTEGER :: ifile,ncid,iret,iv
2437  INTEGER :: start_loop,end_loop
2438  CHARACTER(LEN=70) :: str70
2439  LOGICAL :: l_dbg
2440!---------------------------------------------------------------------
2441  CALL ipsldbg (old_status=l_dbg)
2442!-
2443  IF (l_dbg) WRITE(*,*) 'Entering loop on files :',nb_files
2444!-
2445  IF (PRESENT(fid)) THEN
2446    start_loop = fid
2447    end_loop = fid
2448  ELSE
2449    start_loop = 1
2450    end_loop = nb_files
2451  ENDIF
2452!-
2453  DO ifile=start_loop,end_loop
2454    IF (l_dbg) WRITE(*,*) 'Closing specified file number :',ifile
2455    ncid = ncdf_ids(ifile)
2456    iret = NF90_REDEF (ncid)
2457!---
2458!-- 1. Loop on the number of variables to add some final information
2459!---
2460    IF (l_dbg) WRITE(*,*) 'Entering loop on vars : ',nb_var(ifile)
2461    DO iv=1,nb_var(ifile)
2462      IF (hist_wrt_rng(ifile,iv)) THEN
2463        IF (l_dbg) THEN
2464          WRITE(*,*) 'min value for file :',ifile,' var n. :',iv, &
2465         &           ' is : ',hist_minmax(ifile,iv,1)
2466          WRITE(*,*) 'max value for file :',ifile,' var n. :',iv, &
2467         &           ' is : ',hist_minmax(ifile,iv,2)
2468        ENDIF
2469        IF (hist_calc_rng(ifile,iv)) THEN
2470!-------- Put the min and max values on the file
2471          iret = NF90_PUT_ATT (ncid,ncvar_ids(ifile,iv),'valid_min', &
2472 &                             REAL(hist_minmax(ifile,iv,1),KIND=4))
2473          iret = NF90_PUT_ATT (ncid,ncvar_ids(ifile,iv),'valid_max', &
2474 &                             REAL(hist_minmax(ifile,iv,2),KIND=4))
2475        ENDIF
2476      ENDIF
2477    ENDDO
2478!---
2479!-- 2. Close the file
2480!---
2481    IF (l_dbg) WRITE(*,*) 'close file :',ncid
2482    iret = NF90_CLOSE (ncid)
2483    IF (iret /= NF90_NOERR) THEN
2484      WRITE(str70,'("This file has been already closed :",I3)') ifile
2485      CALL ipslerr (2,'histclo',str70,'','')
2486    ENDIF
2487  ENDDO
2488!---------------------
2489END SUBROUTINE histclo
2490!===
2491SUBROUTINE ioconf_modname (str)
2492!---------------------------------------------------------------------
2493!- This subroutine allows to configure the name
2494!- of the model written into the file
2495!---------------------------------------------------------------------
2496  IMPLICIT NONE
2497!-
2498  CHARACTER(LEN=*),INTENT(IN) :: str
2499!---------------------------------------------------------------------
2500  IF (.NOT.lock_modname) THEN
2501    model_name = str(1:MIN(LEN_TRIM(str),80))
2502    lock_modname = .TRUE.
2503  ELSE
2504    CALL ipslerr (2,"ioconf_modname", &
2505   &  'The model name can only be changed once and only', &
2506   &  'before it is used. It is now set to :',model_name)
2507  ENDIF
2508!----------------------------
2509END SUBROUTINE ioconf_modname
2510!-
2511!===
2512!-
2513!-----------------
2514END MODULE histcom
Note: See TracBrowser for help on using the repository browser.