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

Last change on this file since 953 was 953, checked in by bellier, 13 years ago

Suppress a little misstype in debug mode

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