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

Last change on this file since 1525 was 1525, checked in by mmaipsl, 13 years ago

Keep name of each history file.

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