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

Last change on this file since 416 was 416, checked in by bellier, 14 years ago

Changed an initialization to suppress a bug in a test (from JG)

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