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

Last change on this file since 4 was 4, checked in by rblod, 19 years ago

First import of IOIPSL sources

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