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

Last change on this file since 428 was 427, checked in by bellier, 16 years ago

minor modifications

  • Property svn:keywords set to Id
File size: 78.2 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 histwrite for instance
205!-
206!- OUTPUT
207!-
208!- phoriid   : ID of the horizontal grid
209!- pfileid   : ID of the netcdf file
210!-
211!- Optional INPUT arguments
212!-
213!- domain_id       : Domain identifier
214!-
215!- TO DO
216!-
217!- This package should be written in f90
218!- and use the following features :
219!-    - structures for the meta-data of the files and variables
220!-    - memory allocation as needed
221!-    - Pointers
222!-
223!- VERSION
224!-
225!---------------------------------------------------------------------
226  IMPLICIT NONE
227!-
228  CHARACTER(LEN=*) :: pfilename
229  INTEGER,INTENT(IN) :: pim,pjm
230  REAL,DIMENSION(pim),INTENT(IN) :: plon
231  REAL,DIMENSION(pjm),INTENT(IN) :: plat
232  INTEGER,INTENT(IN):: par_orix,par_szx,par_oriy,par_szy
233  INTEGER,INTENT(IN) :: pitau0
234  REAL,INTENT(IN) :: pdate0,pdeltat
235  INTEGER,INTENT(OUT) :: pfileid,phoriid
236  INTEGER,INTENT(IN),OPTIONAL :: domain_id
237!-
238  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 histwrite 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 histwrite 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 if 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      iret = NF90_PUT_ATT (ncid,ncvarid,'standard_name', &
1619 &                         TRIM(title(pfileid,iv)))
1620!-
1621      iret = NF90_PUT_ATT (ncid,ncvarid,'_Fillvalue', &
1622 &                         REAL(missing_val,KIND=4))
1623      IF (hist_wrt_rng(pfileid,iv)) THEN
1624        iret = NF90_PUT_ATT (ncid,ncvarid,'valid_min', &
1625 &                           REAL(hist_minmax(pfileid,iv,1),KIND=4))
1626        iret = NF90_PUT_ATT (ncid,ncvarid,'valid_max', &
1627 &                           REAL(hist_minmax(pfileid,iv,2),KIND=4))
1628      ENDIF
1629      iret = NF90_PUT_ATT (ncid,ncvarid,'long_name', &
1630 &                         TRIM(title(pfileid,iv)))
1631      iret = NF90_PUT_ATT (ncid,ncvarid,'online_operation', &
1632 &                         TRIM(fullop(pfileid,iv)))
1633!-
1634      SELECT CASE(ndim)
1635      CASE(-3,2:4)
1636      CASE DEFAULT
1637        CALL ipslerr (3,"histend", &
1638       &  'less than 2 or more than 4 dimensions are not', &
1639       &  'allowed at this stage',' ')
1640      END SELECT
1641!-
1642      assoc=TRIM(hax_name(pfileid,var_haxid(pfileid,iv),2))  &
1643 &   //' '//TRIM(hax_name(pfileid,var_haxid(pfileid,iv),1))
1644!-
1645      ziv = var_zaxid(pfileid,iv)
1646      IF (ziv > 0) THEN
1647        str30 = zax_name(pfileid,ziv)
1648        assoc = TRIM(str30)//' '//TRIM(assoc)
1649      ENDIF
1650!-
1651      IF (itax > 0) THEN
1652        IF (nb_tax(pfileid) > 1) THEN
1653          str30 = "t_"//tax_name(pfileid,itax)
1654        ELSE
1655          str30 = "time_counter"
1656        ENDIF
1657        assoc = TRIM(str30)//' '//TRIM(assoc)
1658!-
1659        IF (check) THEN
1660          WRITE(*,*) "histend : 2.0.n, freq_opp, freq_wrt", &
1661 &                   freq_opp(pfileid,iv),freq_wrt(pfileid,iv)
1662        ENDIF
1663!-
1664        iret = NF90_PUT_ATT (ncid,ncvarid,'interval_operation', &
1665 &                           REAL(freq_opp(pfileid,iv),KIND=4))
1666        iret = NF90_PUT_ATT (ncid,ncvarid,'interval_write', &
1667 &                           REAL(freq_wrt(pfileid,iv),KIND=4))
1668      ENDIF
1669      iret = NF90_PUT_ATT (ncid,ncvarid,'coordinates',TRIM(assoc))
1670    ENDIF
1671  ENDDO
1672!-
1673! 2.2 Add DOMAIN attributes if needed
1674!-
1675  IF (dom_id_svg(pfileid) >= 0) THEN
1676    CALL flio_dom_att (ncid,dom_id_svg(pfileid))
1677  ENDIF
1678!-
1679! 3.0 Put the netcdf file into write mode
1680!-
1681  IF (check) WRITE(*,*) "histend : 3.0"
1682!-
1683  iret = NF90_ENDDEF (ncid)
1684!-
1685! 4.0 Give some informations to the user
1686!-
1687  IF (check) WRITE(*,*) "histend : 4.0"
1688!-
1689  WRITE(str70,'("All variables have been initialized on file :",I3)') pfileid
1690  CALL ipslerr (1,'histend',str70,'',' ')
1691!---------------------
1692END SUBROUTINE histend
1693!===
1694SUBROUTINE histwrite_r1d (pfileid,pvarname,pitau,pdata,nbindex,nindex)
1695!---------------------------------------------------------------------
1696  IMPLICIT NONE
1697!-
1698  INTEGER,INTENT(IN) :: pfileid,pitau,nbindex
1699  INTEGER,DIMENSION(nbindex),INTENT(IN) :: nindex
1700  REAL,DIMENSION(:),INTENT(IN) :: pdata
1701  CHARACTER(LEN=*),INTENT(IN) :: pvarname
1702!---------------------------------------------------------------------
1703  CALL histw_rnd (pfileid,pvarname,pitau,nbindex,nindex,pdata_1d=pdata)
1704!---------------------------
1705END SUBROUTINE histwrite_r1d
1706!===
1707SUBROUTINE histwrite_r2d (pfileid,pvarname,pitau,pdata,nbindex,nindex)
1708!---------------------------------------------------------------------
1709  IMPLICIT NONE
1710!-
1711  INTEGER,INTENT(IN) :: pfileid,pitau,nbindex
1712  INTEGER,DIMENSION(nbindex),INTENT(IN) :: nindex
1713  REAL,DIMENSION(:,:),INTENT(IN) :: pdata
1714  CHARACTER(LEN=*),INTENT(IN) :: pvarname
1715!---------------------------------------------------------------------
1716  CALL histw_rnd (pfileid,pvarname,pitau,nbindex,nindex,pdata_2d=pdata)
1717!---------------------------
1718END SUBROUTINE histwrite_r2d
1719!===
1720SUBROUTINE histwrite_r3d (pfileid,pvarname,pitau,pdata,nbindex,nindex)
1721!---------------------------------------------------------------------
1722  IMPLICIT NONE
1723!-
1724  INTEGER,INTENT(IN) :: pfileid,pitau,nbindex
1725  INTEGER,DIMENSION(nbindex),INTENT(IN) :: nindex
1726  REAL,DIMENSION(:,:,:),INTENT(IN) :: pdata
1727  CHARACTER(LEN=*),INTENT(IN) :: pvarname
1728!---------------------------------------------------------------------
1729  CALL histw_rnd (pfileid,pvarname,pitau,nbindex,nindex,pdata_3d=pdata)
1730!---------------------------
1731END SUBROUTINE histwrite_r3d
1732!===
1733SUBROUTINE histw_rnd (pfileid,pvarname,pitau,nbindex,nindex, &
1734  &                   pdata_1d,pdata_2d,pdata_3d)
1735!---------------------------------------------------------------------
1736  IMPLICIT NONE
1737!-
1738  INTEGER,INTENT(IN) :: pfileid,pitau,nbindex
1739  INTEGER,DIMENSION(nbindex),INTENT(IN) :: nindex
1740  CHARACTER(LEN=*),INTENT(IN) :: pvarname
1741  REAL,DIMENSION(:),INTENT(IN),OPTIONAL     :: pdata_1d
1742  REAL,DIMENSION(:,:),INTENT(IN),OPTIONAL   :: pdata_2d
1743  REAL,DIMENSION(:,:,:),INTENT(IN),OPTIONAL :: pdata_3d
1744!-
1745  LOGICAL :: do_oper,do_write,largebuf,l1d,l2d,l3d
1746  INTEGER :: varid,io,nbpt_out
1747  INTEGER              :: nbpt_in1
1748  INTEGER,DIMENSION(2) :: nbpt_in2
1749  INTEGER,DIMENSION(3) :: nbpt_in3
1750  REAL,ALLOCATABLE,SAVE :: buff_tmp(:)
1751  INTEGER,SAVE :: buff_tmp_sz
1752  CHARACTER(LEN=7) :: tmp_opp
1753  CHARACTER(LEN=13) :: c_nam
1754!-
1755  LOGICAL :: check = .FALSE.
1756!---------------------------------------------------------------------
1757  l1d=PRESENT(pdata_1d); l2d=PRESENT(pdata_2d); l3d=PRESENT(pdata_3d);
1758  IF      (l1d) THEN
1759    c_nam = 'histwrite_r1d'
1760  ELSE IF (l2d) THEN
1761    c_nam = 'histwrite_r2d'
1762  ELSE IF (l3d) THEN
1763    c_nam = 'histwrite_r3d'
1764  ENDIF
1765!-
1766! 1.0 Try to catch errors like specifying the wrong file ID.
1767!     Thanks Marine for showing us what errors users can make !
1768!-
1769  IF ( (pfileid < 1).OR.(pfileid > nb_files) ) THEN
1770    CALL ipslerr (3,"histwrite", &
1771 &    'Illegal file ID in the histwrite of variable',pvarname,' ')
1772  ENDIF
1773!-
1774! 1.1 Find the id of the variable to be written and the real time
1775!-
1776  CALL histvar_seq (pfileid,pvarname,varid)
1777!-
1778! 2.0 do nothing for never operation
1779!-
1780  tmp_opp = topp(pfileid,varid)
1781!-
1782  IF (TRIM(tmp_opp) == "never") THEN
1783    last_opp_chk(pfileid,varid) = -99
1784    last_wrt_chk(pfileid,varid) = -99
1785  ENDIF
1786!-
1787! 3.0 We check if we need to do an operation
1788!-
1789  IF (last_opp_chk(pfileid,varid) == pitau) THEN
1790    CALL ipslerr (3,"histwrite", &
1791 &    'This variable as already been analysed at the present', &
1792 &    'time step',' ')
1793  ENDIF
1794!-
1795  CALL isittime &
1796 &  (pitau,date0(pfileid),deltat(pfileid),freq_opp(pfileid,varid), &
1797 &   last_opp(pfileid,varid),last_opp_chk(pfileid,varid),do_oper)
1798!-
1799! 4.0 We check if we need to write the data
1800!-
1801  IF (last_wrt_chk(pfileid,varid) == pitau) THEN
1802    CALL ipslerr (3,"histwrite", &
1803 &    'This variable as already been written for the present', &
1804 &    'time step',' ')
1805  ENDIF
1806!-
1807  CALL isittime &
1808 &  (pitau,date0(pfileid),deltat(pfileid),freq_wrt(pfileid,varid), &
1809 &   last_wrt(pfileid,varid),last_wrt_chk(pfileid,varid),do_write)
1810!-
1811! 5.0 histwrite called
1812!-
1813  IF (do_oper.OR.do_write) THEN
1814!-
1815!-- 5.1 Get the sizes of the data we will handle
1816!-
1817    IF (datasz_in(pfileid,varid,1) <= 0) THEN
1818!---- There is the risk here that the user has over-sized the array.
1819!---- But how can we catch this ?
1820!---- In the worst case we will do impossible operations
1821!---- on part of the data !
1822      datasz_in(pfileid,varid,1:3) = -1
1823      IF      (l1d) THEN
1824        datasz_in(pfileid,varid,1) = SIZE(pdata_1d)
1825      ELSE IF (l2d) THEN
1826        datasz_in(pfileid,varid,1) = SIZE(pdata_2d,DIM=1)
1827        datasz_in(pfileid,varid,2) = SIZE(pdata_2d,DIM=2)
1828      ELSE IF (l3d) THEN
1829        datasz_in(pfileid,varid,1) = SIZE(pdata_3d,DIM=1)
1830        datasz_in(pfileid,varid,2) = SIZE(pdata_3d,DIM=2)
1831        datasz_in(pfileid,varid,3) = SIZE(pdata_3d,DIM=3)
1832      ENDIF
1833    ENDIF
1834!-
1835!-- 5.2 The maximum size of the data will give the size of the buffer
1836!-
1837    IF (datasz_max(pfileid,varid) <= 0) THEN
1838      largebuf = .FALSE.
1839      DO io=1,nbopp(pfileid,varid)
1840        IF (INDEX(fuchnbout,sopps(pfileid,varid,io)) > 0) THEN
1841          largebuf = .TRUE.
1842        ENDIF
1843      ENDDO
1844      IF (largebuf) THEN
1845        datasz_max(pfileid,varid) = &
1846 &        scsize(pfileid,varid,1) &
1847 &       *scsize(pfileid,varid,2) &
1848 &       *scsize(pfileid,varid,3)
1849      ELSE
1850        IF      (l1d) THEN
1851          datasz_max(pfileid,varid) = &
1852 &          datasz_in(pfileid,varid,1)
1853        ELSE IF (l2d) THEN
1854          datasz_max(pfileid,varid) = &
1855 &          datasz_in(pfileid,varid,1) &
1856 &         *datasz_in(pfileid,varid,2)
1857        ELSE IF (l3d) THEN
1858          datasz_max(pfileid,varid) = &
1859 &          datasz_in(pfileid,varid,1) &
1860 &         *datasz_in(pfileid,varid,2) &
1861 &         *datasz_in(pfileid,varid,3)
1862        ENDIF
1863      ENDIF
1864    ENDIF
1865!-
1866    IF (.NOT.ALLOCATED(buff_tmp)) THEN
1867      IF (check) THEN
1868        WRITE(*,*) &
1869 &       c_nam//" : allocate buff_tmp for buff_sz = ", &
1870 &       datasz_max(pfileid,varid)
1871      ENDIF
1872      ALLOCATE (buff_tmp(datasz_max(pfileid,varid)))
1873      buff_tmp_sz = datasz_max(pfileid,varid)
1874    ELSE IF (datasz_max(pfileid,varid) > buff_tmp_sz) THEN
1875      IF (check) THEN
1876        WRITE(*,*) &
1877 &       c_nam//" : re-allocate buff_tmp for buff_sz = ", &
1878 &       datasz_max(pfileid,varid)
1879      ENDIF
1880      DEALLOCATE (buff_tmp)
1881      ALLOCATE (buff_tmp(datasz_max(pfileid,varid)))
1882      buff_tmp_sz = datasz_max(pfileid,varid)
1883    ENDIF
1884!-
1885!-- We have to do the first operation anyway.
1886!-- Thus we do it here and change the ranke
1887!-- of the data at the same time. This should speed up things.
1888!-
1889    nbpt_out = datasz_max(pfileid,varid)
1890    IF      (l1d) THEN
1891      nbpt_in1 = datasz_in(pfileid,varid,1)
1892      CALL mathop (sopps(pfileid,varid,1),nbpt_in1,pdata_1d, &
1893 &                 missing_val,nbindex,nindex, &
1894 &                 scal(pfileid,varid,1),nbpt_out,buff_tmp)
1895    ELSE IF (l2d) THEN
1896      nbpt_in2(1:2) = datasz_in(pfileid,varid,1:2)
1897      CALL mathop (sopps(pfileid,varid,1),nbpt_in2,pdata_2d, &
1898 &                 missing_val,nbindex,nindex, &
1899 &                 scal(pfileid,varid,1),nbpt_out,buff_tmp)
1900    ELSE IF (l3d) THEN
1901      nbpt_in3(1:3) = datasz_in(pfileid,varid,1:3)
1902      CALL mathop (sopps(pfileid,varid,1),nbpt_in3,pdata_3d, &
1903 &                 missing_val,nbindex,nindex, &
1904 &                 scal(pfileid,varid,1),nbpt_out,buff_tmp)
1905    ENDIF
1906    CALL histwrite_real (pfileid,varid,pitau,nbpt_out, &
1907 &            buff_tmp,nbindex,nindex,do_oper,do_write)
1908  ENDIF
1909!-
1910! 6.0 Manage time steps
1911!-
1912  IF ((TRIM(tmp_opp) /= "once").AND.(TRIM(tmp_opp) /= "never")) THEN
1913    last_opp_chk(pfileid,varid) = pitau
1914    last_wrt_chk(pfileid,varid) = pitau
1915  ELSE
1916    last_opp_chk(pfileid,varid) = -99
1917    last_wrt_chk(pfileid,varid) = -99
1918  ENDIF
1919!-----------------------
1920END SUBROUTINE histw_rnd
1921!===
1922SUBROUTINE histwrite_real &
1923 & (pfileid,varid,pitau,nbdpt,buff_tmp,nbindex,nindex,do_oper,do_write)
1924!---------------------------------------------------------------------
1925!- This subroutine is internal and does the calculations and writing
1926!- if needed. At a later stage it should be split into an operation
1927!- and writing subroutines.
1928!---------------------------------------------------------------------
1929  IMPLICIT NONE
1930!-
1931  INTEGER,INTENT(IN) :: pfileid,pitau,varid, &
1932 &                      nbindex,nindex(nbindex),nbdpt
1933  REAL,DIMENSION(:)  :: buff_tmp
1934  LOGICAL,INTENT(IN) :: do_oper,do_write
1935!-
1936  INTEGER :: tsz,ncid,ncvarid,i,iret,ipt,itax,io,nbin,nbout
1937  INTEGER,DIMENSION(4) :: corner,edges
1938  INTEGER :: itime
1939!-
1940  REAL :: rtime
1941  CHARACTER(LEN=7) :: tmp_opp
1942  REAL,ALLOCATABLE,SAVE :: buff_tmp2(:)
1943  INTEGER,SAVE          :: buff_tmp2_sz
1944  REAL,ALLOCATABLE,SAVE :: buffer_used(:)
1945  INTEGER,SAVE          :: buffer_sz
1946!-
1947  LOGICAL :: check = .FALSE.
1948!---------------------------------------------------------------------
1949  IF (check) THEN
1950    WRITE(*,*) "histwrite 0.0 :  VAR : ",name(pfileid,varid)
1951    WRITE(*,*) "histwrite 0.0 : nbindex,nindex :", &
1952    &  nbindex,nindex(1:MIN(3,nbindex)),'...',nindex(MAX(1,nbindex-3):nbindex)
1953  ENDIF
1954!-
1955! The sizes which can be encoutered
1956!-
1957  tsz = zsize(pfileid,varid,1)*zsize(pfileid,varid,2)*zsize(pfileid,varid,3)
1958!-
1959! 1.0 We allocate the memory needed to store the data between write
1960!     and the temporary space needed for operations.
1961!     We have to keep precedent buffer if needed
1962!-
1963  IF (.NOT. ALLOCATED(buffer)) THEN
1964    IF (check) WRITE(*,*) "histwrite_real 1.0 allocate buffer ",buff_pos
1965    ALLOCATE(buffer(buff_pos))
1966    buffer_sz = buff_pos
1967    buffer(:)=0.0
1968  ELSE IF (buffer_sz < buff_pos) THEN
1969    IF (check) WRITE(*,*) "histwrite_real 1.0.1 re-allocate buffer for ",buff_pos," instead of ",SIZE(buffer)
1970    IF (SUM(buffer)/=0.0) THEN
1971      IF (check) WRITE (*,*) ' histwrite : buffer has been used. We have to save it before re-allocating it '
1972      ALLOCATE (buffer_used(buffer_sz))
1973      buffer_used(:)=buffer(:)
1974      DEALLOCATE (buffer)
1975      ALLOCATE (buffer(buff_pos))
1976      buffer_sz = buff_pos
1977      buffer(:)=0.0
1978      buffer(:SIZE(buffer_used))=buffer_used
1979      DEALLOCATE (buffer_used)
1980    ELSE
1981      IF (check) WRITE (*,*) ' histwrite : buffer has not been used. We have just to re-allocating it '
1982      DEALLOCATE (buffer)
1983      ALLOCATE (buffer(buff_pos))
1984      buffer_sz = buff_pos
1985      buffer(:)=0.0
1986    ENDIF
1987  ENDIF
1988!-
1989! The buffers are only deallocated when more space is needed. This
1990! reduces the umber of allocates but increases memory needs.
1991!-
1992  IF (.NOT.ALLOCATED(buff_tmp2)) THEN
1993    IF (check) THEN
1994      WRITE(*,*) "histwrite_real 1.1 allocate buff_tmp2 ",SIZE(buff_tmp)
1995    ENDIF
1996    ALLOCATE (buff_tmp2(datasz_max(pfileid,varid)))
1997    buff_tmp2_sz = datasz_max(pfileid,varid)
1998  ELSE IF (datasz_max(pfileid,varid) > buff_tmp2_sz) THEN
1999    IF (check) THEN
2000      WRITE(*,*) "histwrite_real 1.2 re-allocate buff_tmp2 : ", &
2001     & SIZE(buff_tmp)," instead of ",SIZE(buff_tmp2)
2002    ENDIF
2003    DEALLOCATE (buff_tmp2)
2004    ALLOCATE (buff_tmp2(datasz_max(pfileid,varid)))
2005    buff_tmp2_sz = datasz_max(pfileid,varid)
2006  ENDIF
2007!-
2008  rtime = pitau * deltat(pfileid)
2009  tmp_opp = topp(pfileid,varid)
2010!-
2011! 3.0 Do the operations or transfer the slab of data into buff_tmp
2012!-
2013  IF (check) WRITE(*,*) "histwrite: 3.0",pfileid
2014!-
2015! 3.1 DO the Operations only if needed
2016!-
2017  IF (do_oper) THEN
2018    i = pfileid
2019    nbout = nbdpt
2020!-
2021!-- 3.4 We continue the sequence of operations
2022!--     we started in the interface routine
2023!-
2024    DO io = 2,nbopp(i,varid),2
2025      nbin = nbout
2026      nbout = datasz_max(i,varid)
2027      CALL mathop(sopps(i,varid,io),nbin,buff_tmp,missing_val, &
2028 &      nbindex,nindex,scal(i,varid,io),nbout,buff_tmp2)
2029      IF (check) THEN
2030        WRITE(*,*) &
2031 &       "histwrite: 3.4a nbout : ",nbin,nbout,sopps(i,varid,io)
2032      ENDIF
2033!-
2034      nbin = nbout
2035      nbout = datasz_max(i,varid)
2036      CALL mathop(sopps(i,varid,io+1),nbin,buff_tmp2,missing_val, &
2037 &      nbindex,nindex,scal(i,varid,io+1),nbout,buff_tmp)
2038      IF (check) THEN
2039        WRITE(*,*) &
2040 &       "histwrite: 3.4b nbout : ",nbin,nbout,sopps(i,varid,io+1)
2041      ENDIF
2042    ENDDO
2043!-
2044!   3.5 Zoom into the data
2045!-
2046    IF (check) THEN
2047      WRITE(*,*) &
2048 &     "histwrite: 3.5 size(buff_tmp) : ",SIZE(buff_tmp)
2049      WRITE(*,*) &
2050 &     "histwrite: 3.5 slab in X :",zorig(i,varid,1),zsize(i,varid,1)
2051      WRITE(*,*) &
2052 &     "histwrite: 3.5 slab in Y :",zorig(i,varid,2),zsize(i,varid,2)
2053      WRITE(*,*) &
2054 &     "histwrite: 3.5 slab in Z :",zorig(i,varid,3),zsize(i,varid,3)
2055      WRITE(*,*) &
2056 &     "histwrite: 3.5 slab of input:", &
2057 &     scsize(i,varid,1),scsize(i,varid,2),scsize(i,varid,3)
2058    ENDIF
2059    CALL trans_buff &
2060 &      (zorig(i,varid,1),zsize(i,varid,1), &
2061 &       zorig(i,varid,2),zsize(i,varid,2), &
2062 &       zorig(i,varid,3),zsize(i,varid,3), &
2063 &       scsize(i,varid,1),scsize(i,varid,2),scsize(i,varid,3), &
2064 &       buff_tmp,buff_tmp2_sz,buff_tmp2)
2065!-
2066!-- 4.0 Get the min and max of the field (buff_tmp)
2067!-
2068    IF (check) WRITE(*,*) "histwrite: 4.0 buff_tmp",pfileid,varid, &
2069 &    TRIM(tmp_opp),'----',LEN_TRIM(tmp_opp),nbindex
2070!-
2071    IF (hist_calc_rng(pfileid,varid)) THEN
2072      hist_minmax(pfileid,varid,1) = &
2073 &      MIN(hist_minmax(pfileid,varid,1), &
2074 &      MINVAL(buff_tmp2(1:tsz),MASK=buff_tmp2(1:tsz) /= missing_val))
2075      hist_minmax(pfileid,varid,2) = &
2076 &      MAX(hist_minmax(pfileid,varid,2), &
2077 &      MAXVAL(buff_tmp2(1:tsz),MASK=buff_tmp2(1:tsz) /= missing_val))
2078    ENDIF
2079!-
2080!-- 5.0 Do the operations if needed. In the case of instantaneous
2081!--     output we do not transfer to the buffer.
2082!-
2083    IF (check) WRITE(*,*) "histwrite: 5.0",pfileid,"tsz :",tsz
2084!-
2085    ipt = point(pfileid,varid)
2086!-
2087!   WRITE(*,*) 'OPE ipt, buffer :',pvarname,ipt,varid
2088!-
2089    IF (     (TRIM(tmp_opp) /= "inst") &
2090   &    .AND.(TRIM(tmp_opp) /= "once") ) THEN
2091      CALL moycum(tmp_opp,tsz,buffer(ipt:), &
2092     &       buff_tmp2,nb_opp(pfileid,varid))
2093    ENDIF
2094!-
2095    last_opp(pfileid,varid) = pitau
2096    nb_opp(pfileid,varid) = nb_opp(pfileid,varid)+1
2097!-
2098   ENDIF
2099!-
2100! 6.0 Write to file if needed
2101!-
2102  IF (check) WRITE(*,*) "histwrite: 6.0",pfileid
2103!-
2104  IF (do_write) THEN
2105!-
2106    ncvarid = ncvar_ids(pfileid,varid)
2107    ncid = ncdf_ids(pfileid)
2108!-
2109!-- 6.1 Do the operations that are needed before writting
2110!-
2111    IF (check) WRITE(*,*) "histwrite: 6.1",pfileid
2112!-
2113    IF (     (TRIM(tmp_opp) /= "inst") &
2114   &    .AND.(TRIM(tmp_opp) /= "once") ) THEN
2115      rtime = (rtime+last_wrt(pfileid,varid)*deltat(pfileid))/2.0
2116    ENDIF
2117!-
2118!-- 6.2 Add a value to the time axis of this variable if needed
2119!-
2120    IF (     (TRIM(tmp_opp) /= "l_max") &
2121   &    .AND.(TRIM(tmp_opp) /= "l_min") &
2122   &    .AND.(TRIM(tmp_opp) /= "once") ) THEN
2123!-
2124      IF (check) WRITE(*,*) "histwrite: 6.2",pfileid
2125!-
2126      itax = var_axid(pfileid,varid)
2127      itime = nb_wrt(pfileid,varid)+1
2128!-
2129      IF (tax_last(pfileid,itax) < itime) THEN
2130        iret = NF90_PUT_VAR (ncid,tdimid(pfileid,itax),(/ rtime /), &
2131 &                            start=(/ itime /),count=(/ 1 /))
2132        tax_last(pfileid,itax) = itime
2133      ENDIF
2134    ELSE
2135      itime=1
2136    ENDIF
2137!-
2138!-- 6.3 Write the data. Only in the case of instantaneous output
2139!       we do not write the buffer.
2140!-
2141    IF (check) THEN
2142      WRITE(*,*) "histwrite: 6.3",pfileid,ncid,ncvarid,varid,itime
2143    ENDIF
2144!-
2145    IF (scsize(pfileid,varid,3) == 1) THEN
2146      IF (regular(pfileid)) THEN
2147        corner(1:4) = (/ 1,1,itime,0 /)
2148        edges(1:4) = (/ zsize(pfileid,varid,1), &
2149 &                      zsize(pfileid,varid,2),1,0 /)
2150      ELSE
2151        corner(1:4) = (/ 1,itime,0,0 /)
2152        edges(1:4) = (/ zsize(pfileid,varid,1),1,0,0 /)
2153      ENDIF
2154    ELSE
2155      IF (regular(pfileid)) THEN
2156        corner(1:4) = (/ 1,1,1,itime /)
2157        edges(1:4) = (/ zsize(pfileid,varid,1), &
2158 &                      zsize(pfileid,varid,2), &
2159 &                      zsize(pfileid,varid,3),1 /)
2160      ELSE
2161        corner(1:4) = (/ 1,1,itime,0 /)
2162        edges(1:4) = (/ zsize(pfileid,varid,1), &
2163 &                      zsize(pfileid,varid,3),1,0 /)
2164      ENDIF
2165    ENDIF
2166!-
2167    ipt = point(pfileid,varid)
2168!-
2169    IF (     (TRIM(tmp_opp) /= "inst") &
2170 &      .AND.(TRIM(tmp_opp) /= "once") ) THEN
2171      iret = NF90_PUT_VAR (ncid,ncvarid,buffer(ipt:), &
2172 &                       start=corner(1:4),count=edges(1:4))
2173    ELSE
2174      iret = NF90_PUT_VAR (ncid,ncvarid,buff_tmp2, &
2175 &                       start=corner(1:4),count=edges(1:4))
2176    ENDIF
2177!-
2178    last_wrt(pfileid,varid) = pitau
2179    nb_wrt(pfileid,varid) = nb_wrt(pfileid,varid)+1
2180    nb_opp(pfileid,varid) = 0
2181!---
2182!   After the write the file can be synchronized so that no data is
2183!   lost in case of a crash. This feature gives up on the benefits of
2184!   buffering and should only be used in debuging mode. A flag is
2185!   needed here to switch to this mode.
2186!---
2187!   iret = NF90_SYNC (ncid)
2188!-
2189  ENDIF
2190!----------------------------
2191END SUBROUTINE histwrite_real
2192!===
2193SUBROUTINE histvar_seq (pfid,pvarname,pvid)
2194!---------------------------------------------------------------------
2195!- This subroutine optimized the search for the variable in the table.
2196!- In a first phase it will learn the succession of the variables
2197!- called and then it will use the table to guess what comes next.
2198!- It is the best solution to avoid lengthy searches through array
2199!- vectors.
2200!-
2201!- ARGUMENTS :
2202!-
2203!- pfid  : id of the file on which we work
2204!- pvarname : The name of the variable we are looking for
2205!- pvid     : The var id we found
2206!---------------------------------------------------------------------
2207  IMPLICIT NONE
2208!-
2209  INTEGER,INTENT(in)  :: pfid
2210  CHARACTER(LEN=*),INTENT(IN) :: pvarname
2211  INTEGER,INTENT(out) :: pvid
2212!-
2213  LOGICAL,SAVE :: learning(nb_files_max)=.TRUE.
2214  INTEGER,SAVE :: overlap(nb_files_max) = -1
2215  INTEGER,SAVE :: varseq(nb_files_max,nb_var_max*3)
2216  INTEGER,SAVE :: varseq_len(nb_files_max) = 0
2217  INTEGER,SAVE :: varseq_pos(nb_files_max)
2218  INTEGER,SAVE :: varseq_err(nb_files_max) = 0
2219  INTEGER      :: ib,sp,nx,pos
2220  CHARACTER(LEN=70) :: str70
2221!-
2222  LOGICAL :: check = .FALSE.
2223!---------------------------------------------------------------------
2224  IF (check) THEN
2225    WRITE(*,*) 'histvar_seq, start of the subroutine :',learning(pfid)
2226  ENDIF
2227!-
2228  IF (learning(pfid)) THEN
2229!-
2230!-- 1.0 We compute the length over which we are going
2231!--     to check the overlap
2232!-
2233    IF (overlap(pfid) <= 0) THEN
2234      IF (nb_var(pfid) > 6) THEN
2235        overlap(pfid) = nb_var(pfid)/3*2
2236      ELSE
2237        overlap(pfid) = nb_var(pfid)
2238      ENDIF
2239    ENDIF
2240!-
2241!-- 1.1 Find the position of this string
2242!-
2243    CALL find_str (name(pfid,1:nb_var(pfid)),pvarname,pos)
2244    IF (pos > 0) THEN
2245      pvid = pos
2246    ELSE
2247      CALL ipslerr (3,"histvar_seq", &
2248 &      'The name of the variable you gave has not been declared', &
2249 &      'You should use subroutine histdef for declaring variable', &
2250 &      TRIM(pvarname))
2251    ENDIF
2252!-
2253!-- 1.2 If we have not given up we store the position
2254!--     in the sequence of calls
2255!-
2256    IF (varseq_err(pfid) >= 0) THEN
2257      sp = varseq_len(pfid)+1
2258      IF (sp <= nb_var_max*3) THEN
2259        varseq(pfid,sp) = pvid
2260        varseq_len(pfid) = sp
2261      ELSE
2262        CALL ipslerr (2,"histvar_seq",&
2263 &       'The learning process has failed and we give up. '// &
2264 &       'Either you sequence is',&
2265 &       'too complex or I am too dumb. '// &
2266 &       'This will only affect the efficiency',&
2267 &       'of your code. Thus if you wish to save time'// &
2268 &       ' contact the IOIPSL team. ')
2269        WRITE(*,*) 'The sequence we have found up to now :'
2270        WRITE(*,*) varseq(pfid,1:sp-1)
2271        varseq_err(pfid) = -1
2272      ENDIF
2273!-
2274!---- 1.3 Check if we have found the right overlap
2275!-
2276      IF (varseq_len(pfid) .GE. overlap(pfid)*2) THEN
2277!-
2278!------ We skip a few variables if needed as they could come
2279!------ from the initialisation of the model.
2280!-
2281        DO ib = 0,sp-overlap(pfid)*2
2282          IF ( learning(pfid) .AND.&
2283            & SUM(ABS(varseq(pfid,ib+1:ib+overlap(pfid)) -&
2284            & varseq(pfid,sp-overlap(pfid)+1:sp))) == 0 ) THEN
2285            learning(pfid) = .FALSE.
2286            varseq_len(pfid) = sp-overlap(pfid)-ib
2287            varseq_pos(pfid) = overlap(pfid)+ib
2288            varseq(pfid,1:varseq_len(pfid)) = &
2289 &            varseq(pfid,ib+1:ib+varseq_len(pfid))
2290          ENDIF
2291        ENDDO
2292      ENDIF
2293    ENDIF
2294  ELSE
2295!-
2296!-- 2.0 Now we know how the calls to histwrite are sequenced
2297!--     and we can get a guess at the var ID
2298!-
2299    nx = varseq_pos(pfid)+1
2300    IF (nx > varseq_len(pfid)) nx = 1
2301!-
2302    pvid = varseq(pfid,nx)
2303!-
2304    IF (TRIM(name(pfid,pvid)) /= TRIM(pvarname)) THEN
2305      CALL find_str (name(pfid,1:nb_var(pfid)),pvarname,pos)
2306      IF (pos > 0) THEN
2307        pvid = pos
2308      ELSE
2309        CALL ipslerr (3,"histvar_seq", &
2310 &  'The name of the variable you gave has not been declared',&
2311 &  'You should use subroutine histdef for declaring variable', &
2312 &  TRIM(pvarname))
2313      ENDIF
2314      varseq_err(pfid) = varseq_err(pfid)+1
2315    ELSE
2316!-
2317!---- We only keep the new position if we have found the variable
2318!---- this way. This way an out of sequence call to histwrite does
2319!---- not defeat the process.
2320!-
2321      varseq_pos(pfid) = nx
2322    ENDIF
2323!-
2324    IF (varseq_err(pfid) .GE. 10) THEN
2325      WRITE(str70,'("for file ",I3)') pfid
2326      CALL ipslerr (2,"histvar_seq", &
2327 &  'There were 10 errors in the learned sequence of variables',&
2328 &  str70,'This looks like a bug, please report it.')
2329         varseq_err(pfid) = 0
2330    ENDIF
2331  ENDIF
2332!-
2333  IF (check) THEN
2334    WRITE(*,*) &
2335 &   'histvar_seq, end of the subroutine :',TRIM(pvarname),pvid
2336  ENDIF
2337!-------------------------
2338END SUBROUTINE histvar_seq
2339!===
2340SUBROUTINE histsync (file)
2341!---------------------------------------------------------------------
2342!- This subroutine will synchronise all
2343!- (or one if defined) opened files.
2344!-
2345!- VERSION
2346!-
2347!---------------------------------------------------------------------
2348  IMPLICIT NONE
2349!-
2350! file  : optional argument for fileid
2351  INTEGER,INTENT(in),OPTIONAL :: file
2352!-
2353  INTEGER :: ifile,ncid,iret
2354!-
2355  LOGICAL :: file_exists
2356  LOGICAL :: check = .FALSE.
2357!---------------------------------------------------------------------
2358  IF (check) WRITE(*,*) 'Entering loop on files :',nb_files
2359!-
2360! 1.The loop on files to synchronise
2361!-
2362  DO ifile = 1,nb_files
2363!-
2364    IF (PRESENT(file)) THEN
2365      file_exists = (ifile == file)
2366    ELSE
2367      file_exists = .TRUE.
2368    ENDIF
2369!-
2370    IF (file_exists) THEN
2371      IF (check) THEN
2372        WRITE(*,*) 'Synchronising specified file number :',file
2373      ENDIF
2374      ncid = ncdf_ids(ifile)
2375      iret = NF90_SYNC (ncid)
2376    ENDIF
2377!-
2378  ENDDO
2379!----------------------
2380END SUBROUTINE histsync
2381!===
2382SUBROUTINE histclo (fid)
2383!---------------------------------------------------------------------
2384!- This subroutine will close all (or one if defined) opened files
2385!-
2386!- VERSION
2387!-
2388!---------------------------------------------------------------------
2389  IMPLICIT NONE
2390!-
2391! fid  : optional argument for fileid
2392  INTEGER,INTENT(in),OPTIONAL :: fid
2393!-
2394  INTEGER :: ifile,ncid,iret,iv
2395  INTEGER :: start_loop,end_loop
2396  CHARACTER(LEN=70) :: str70
2397!-
2398  LOGICAL :: check=.FALSE.
2399!---------------------------------------------------------------------
2400  IF (check) WRITE(*,*) 'Entering loop on files :',nb_files
2401!-
2402  IF (PRESENT(fid)) THEN
2403    start_loop = fid
2404    end_loop = fid
2405  ELSE
2406    start_loop = 1
2407    end_loop = nb_files
2408  ENDIF
2409!-
2410  DO ifile=start_loop,end_loop
2411    IF (check) WRITE(*,*) 'Closing specified file number :',ifile
2412    ncid = ncdf_ids(ifile)
2413    iret = NF90_REDEF (ncid)
2414!---
2415!-- 1. Loop on the number of variables to add some final information
2416!---
2417    IF (check) WRITE(*,*) 'Entering loop on vars : ',nb_var(ifile)
2418    DO iv=1,nb_var(ifile)
2419      IF (hist_wrt_rng(ifile,iv)) THEN
2420        IF (check) THEN
2421          WRITE(*,*) 'min value for file :',ifile,' var n. :',iv, &
2422         &           ' is : ',hist_minmax(ifile,iv,1)
2423          WRITE(*,*) 'max value for file :',ifile,' var n. :',iv, &
2424         &           ' is : ',hist_minmax(ifile,iv,2)
2425        ENDIF
2426        IF (hist_calc_rng(ifile,iv)) THEN
2427!-------- Put the min and max values on the file
2428          iret = NF90_PUT_ATT (ncid,ncvar_ids(ifile,iv),'valid_min', &
2429 &                             REAL(hist_minmax(ifile,iv,1),KIND=4))
2430          iret = NF90_PUT_ATT (ncid,ncvar_ids(ifile,iv),'valid_max', &
2431 &                             REAL(hist_minmax(ifile,iv,2),KIND=4))
2432        ENDIF
2433      ENDIF
2434    ENDDO
2435!---
2436!-- 2. Close the file
2437!---
2438    IF (check) WRITE(*,*) 'close file :',ncid
2439    iret = NF90_CLOSE (ncid)
2440    IF (iret /= NF90_NOERR) THEN
2441      WRITE(str70,'("This file has been already closed :",I3)') ifile
2442      CALL ipslerr (2,'histclo',str70,'','')
2443    ENDIF
2444  ENDDO
2445!---------------------
2446END SUBROUTINE histclo
2447!===
2448SUBROUTINE ioconf_modname (str)
2449!---------------------------------------------------------------------
2450!- This subroutine allows to configure the name
2451!- of the model written into the file
2452!---------------------------------------------------------------------
2453  IMPLICIT NONE
2454!-
2455  CHARACTER(LEN=*),INTENT(IN) :: str
2456!---------------------------------------------------------------------
2457  IF (.NOT.lock_modname) THEN
2458    model_name = str(1:MIN(LEN_TRIM(str),80))
2459    lock_modname = .TRUE.
2460  ELSE
2461    CALL ipslerr (2,"ioconf_modname", &
2462   &  'The model name can only be changed once and only', &
2463   &  'before it is used. It is now set to :',model_name)
2464  ENDIF
2465!----------------------------
2466END SUBROUTINE ioconf_modname
2467!-
2468!===
2469!-
2470!-----------------
2471END MODULE histcom
Note: See TracBrowser for help on using the repository browser.