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

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

New handling of the time_buffer : one buffer for each variable

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