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

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

JB: some cleaning (-> fortran 90)

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