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

Last change on this file since 752 was 752, checked in by bellier, 13 years ago

Added the possibility to choose the kind of NETCDF variables (R4/R8) by calling histdef

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