New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diawri.F90 in branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/UKMO/r6232_INGV1_WAVE-coupling/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 7863

Last change on this file since 7863 was 7600, checked in by jcastill, 7 years ago

Changes as in branches/2016/dev_merge_2016@7562 (branches/2016/dev_INGV_UKMO_2016@7558)
File size: 58.8 KB
Line 
1MODULE diawri
2   !!======================================================================
3   !!                     ***  MODULE  diawri  ***
4   !! Ocean diagnostics :  write ocean output files
5   !!=====================================================================
6   !! History :  OPA  ! 1991-03  (M.-A. Foujols)  Original code
7   !!            4.0  ! 1991-11  (G. Madec)
8   !!                 ! 1992-06  (M. Imbard)  correction restart file
9   !!                 ! 1992-07  (M. Imbard)  split into diawri and rstwri
10   !!                 ! 1993-03  (M. Imbard)  suppress writibm
11   !!                 ! 1998-01  (C. Levy)  NETCDF format using ioipsl INTERFACE
12   !!                 ! 1999-02  (E. Guilyardi)  name of netCDF files + variables
13   !!            8.2  ! 2000-06  (M. Imbard)  Original code (diabort.F)
14   !!   NEMO     1.0  ! 2002-06  (A.Bozec, E. Durand)  Original code (diainit.F)
15   !!             -   ! 2002-09  (G. Madec)  F90: Free form and module
16   !!             -   ! 2002-12  (G. Madec)  merge of diabort and diainit, F90
17   !!                 ! 2005-11  (V. Garnier) Surface pressure gradient organization
18   !!            3.2  ! 2008-11  (B. Lemaire) creation from old diawri
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   dia_wri       : create the standart output files
23   !!   dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
24   !!----------------------------------------------------------------------
25   USE oce             ! ocean dynamics and tracers
26   USE dom_oce         ! ocean space and time domain
27   USE dynadv, ONLY: ln_dynadv_vec
28   USE zdf_oce         ! ocean vertical physics
29   USE ldftra_oce      ! ocean active tracers: lateral physics
30   USE ldfdyn_oce      ! ocean dynamics: lateral physics
31   USE traldf_iso_grif, ONLY : psix_eiv, psiy_eiv
32   USE sol_oce         ! solver variables
33   USE sbc_oce         ! Surface boundary condition: ocean fields
34   USE sbc_ice         ! Surface boundary condition: ice fields
35   USE icb_oce         ! Icebergs
36   USE icbdia          ! Iceberg budgets
37   USE sbcssr          ! restoring term toward SST/SSS climatology
38   USE phycst          ! physical constants
39   USE zdfmxl          ! mixed layer
40   USE dianam          ! build name of file (routine)
41   USE zdfddm          ! vertical  physics: double diffusion
42   USE diahth          ! thermocline diagnostics
43   USE sbcwave         ! wave parameters
44   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
45   USE in_out_manager  ! I/O manager
46   USE diadimg         ! dimg direct access file format output
47   USE iom
48   USE ioipsl
49   USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities     
50
51#if defined key_lim2
52   USE limwri_2 
53#elif defined key_lim3
54   USE limwri 
55#endif
56   USE lib_mpp         ! MPP library
57   USE timing          ! preformance summary
58   USE wrk_nemo        ! working array
59
60   IMPLICIT NONE
61   PRIVATE
62
63   PUBLIC   dia_wri                 ! routines called by step.F90
64   PUBLIC   dia_wri_state
65   PUBLIC   dia_wri_alloc           ! Called by nemogcm module
66
67   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
68   INTEGER ::          nb_T              , ndim_bT   ! grid_T file
69   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
70   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
71   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file
72   INTEGER ::   ndex(1)                              ! ???
73   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
74   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
75   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
76
77   !! * Substitutions
78#  include "zdfddm_substitute.h90"
79#  include "domzgr_substitute.h90"
80#  include "vectopt_loop_substitute.h90"
81   !!----------------------------------------------------------------------
82   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
83   !! $Id$
84   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
85   !!----------------------------------------------------------------------
86CONTAINS
87
88   INTEGER FUNCTION dia_wri_alloc()
89      !!----------------------------------------------------------------------
90      INTEGER, DIMENSION(2) :: ierr
91      !!----------------------------------------------------------------------
92      ierr = 0
93      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
94         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
95         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
96         !
97      dia_wri_alloc = MAXVAL(ierr)
98      IF( lk_mpp )   CALL mpp_sum( dia_wri_alloc )
99      !
100  END FUNCTION dia_wri_alloc
101
102#if defined key_dimgout
103   !!----------------------------------------------------------------------
104   !!   'key_dimgout'                                      DIMG output file
105   !!----------------------------------------------------------------------
106#   include "diawri_dimg.h90"
107
108#else
109   !!----------------------------------------------------------------------
110   !!   Default option                                   NetCDF output file
111   !!----------------------------------------------------------------------
112# if defined key_iomput
113   !!----------------------------------------------------------------------
114   !!   'key_iomput'                                        use IOM library
115   !!----------------------------------------------------------------------
116
117   SUBROUTINE dia_wri( kt )
118      !!---------------------------------------------------------------------
119      !!                  ***  ROUTINE dia_wri  ***
120      !!                   
121      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
122      !!      NETCDF format is used by default
123      !!
124      !! ** Method  :  use iom_put
125      !!----------------------------------------------------------------------
126      !!
127      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
128      !!
129      INTEGER                      ::   ji, jj, jk              ! dummy loop indices
130      INTEGER                      ::   jkbot                   !
131      REAL(wp)                     ::   zztmp, zztmpx, zztmpy   !
132      !!
133      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace
134      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace
135      !!----------------------------------------------------------------------
136      !
137      IF( nn_timing == 1 )   CALL timing_start('dia_wri')
138      !
139      CALL wrk_alloc( jpi , jpj      , z2d )
140      CALL wrk_alloc( jpi , jpj, jpk , z3d )
141      !
142      ! Output the initial state and forcings
143      IF( ninist == 1 ) THEN                       
144         CALL dia_wri_state( 'output.init', kt )
145         ninist = 0
146      ENDIF
147
148      IF( .NOT.lk_vvl ) THEN
149         CALL iom_put( "e3t" , fse3t_n(:,:,:) )
150         CALL iom_put( "e3u" , fse3u_n(:,:,:) )
151         CALL iom_put( "e3v" , fse3v_n(:,:,:) )
152         CALL iom_put( "e3w" , fse3w_n(:,:,:) )
153      ENDIF
154
155      CALL iom_put( "ssh" , sshn )                 ! sea surface height
156      if( iom_use('ssh2') )   CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of sea surface height
157     
158      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
159      CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature
160      IF ( iom_use("sbt") ) THEN
161         DO jj = 1, jpj
162            DO ji = 1, jpi
163               jkbot = mbkt(ji,jj)
164               z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem)
165            END DO
166         END DO
167         CALL iom_put( "sbt", z2d )                ! bottom temperature
168      ENDIF
169     
170      CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity
171      CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity
172      IF ( iom_use("sbs") ) THEN
173         DO jj = 1, jpj
174            DO ji = 1, jpi
175               jkbot = mbkt(ji,jj)
176               z2d(ji,jj) = tsn(ji,jj,jkbot,jp_sal)
177            END DO
178         END DO
179         CALL iom_put( "sbs", z2d )                ! bottom salinity
180      ENDIF
181
182      IF ( iom_use("taubot") ) THEN                ! bottom stress
183         z2d(:,:) = 0._wp
184         DO jj = 2, jpjm1
185            DO ji = fs_2, fs_jpim1   ! vector opt.
186               zztmpx = (  bfrua(ji  ,jj) * un(ji  ,jj,mbku(ji  ,jj))  &
187                      &  + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj))  )     
188               zztmpy = (  bfrva(ji,  jj) * vn(ji,jj  ,mbkv(ji,jj  ))  &
189                      &  + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1))  ) 
190               z2d(ji,jj) = rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1) 
191               !
192            ENDDO
193         ENDDO
194         CALL lbc_lnk( z2d, 'T', 1. )
195         CALL iom_put( "taubot", z2d )           
196      ENDIF
197         
198      CALL iom_put( "uoce", un(:,:,:)         )    ! 3D i-current
199      CALL iom_put(  "ssu", un(:,:,1)         )    ! surface i-current
200      IF ( iom_use("sbu") ) THEN
201         DO jj = 1, jpj
202            DO ji = 1, jpi
203               jkbot = mbku(ji,jj)
204               z2d(ji,jj) = un(ji,jj,jkbot)
205            END DO
206         END DO
207         CALL iom_put( "sbu", z2d )                ! bottom i-current
208      ENDIF
209#if defined key_dynspg_ts
210      CALL iom_put(  "ubar", un_adv(:,:)      )    ! barotropic i-current
211#else
212      CALL iom_put(  "ubar", un_b(:,:)        )    ! barotropic i-current
213#endif
214     
215      CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current
216      CALL iom_put(  "ssv", vn(:,:,1)         )    ! surface j-current
217      IF ( iom_use("sbv") ) THEN
218         DO jj = 1, jpj
219            DO ji = 1, jpi
220               jkbot = mbkv(ji,jj)
221               z2d(ji,jj) = vn(ji,jj,jkbot)
222            END DO
223         END DO
224         CALL iom_put( "sbv", z2d )                ! bottom j-current
225      ENDIF
226#if defined key_dynspg_ts
227      CALL iom_put(  "vbar", vn_adv(:,:)      )    ! barotropic j-current
228#else
229      CALL iom_put(  "vbar", vn_b(:,:)        )    ! barotropic j-current
230#endif
231
232      CALL iom_put( "woce", wn )                   ! vertical velocity
233      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
234         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
235         z2d(:,:) = rau0 * e12t(:,:)
236         DO jk = 1, jpk
237            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
238         END DO
239         CALL iom_put( "w_masstr" , z3d ) 
240         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
241      ENDIF
242
243      CALL iom_put( "avt" , avt                        )    ! T vert. eddy diff. coef.
244      CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef.
245      CALL iom_put( "avs" , fsavs(:,:,:)               )    ! S vert. eddy diff. coef. (useful only with key_zdfddm)
246
247      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
248         DO jj = 2, jpjm1                                    ! sst gradient
249            DO ji = fs_2, fs_jpim1   ! vector opt.
250               zztmp      = tsn(ji,jj,1,jp_tem)
251               zztmpx     = ( tsn(ji+1,jj  ,1,jp_tem) - zztmp ) / e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,1,jp_tem) ) / e1u(ji-1,jj  )
252               zztmpy     = ( tsn(ji  ,jj+1,1,jp_tem) - zztmp ) / e2v(ji,jj) + ( zztmp - tsn(ji  ,jj-1,1,jp_tem) ) / e2v(ji  ,jj-1)
253               z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
254                  &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
255            END DO
256         END DO
257         CALL lbc_lnk( z2d, 'T', 1. )
258         CALL iom_put( "sstgrad2",  z2d               )    ! square of module of sst gradient
259         z2d(:,:) = SQRT( z2d(:,:) )
260         CALL iom_put( "sstgrad" ,  z2d               )    ! module of sst gradient
261      ENDIF
262         
263      ! clem: heat and salt content
264      IF( iom_use("heatc") ) THEN
265         z2d(:,:)  = 0._wp 
266         DO jk = 1, jpkm1
267            DO jj = 1, jpj
268               DO ji = 1, jpi
269                  z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
270               END DO
271            END DO
272         END DO
273         CALL iom_put( "heatc", (rau0 * rcp) * z2d )    ! vertically integrated heat content (J/m2)
274      ENDIF
275
276      IF( iom_use("saltc") ) THEN
277         z2d(:,:)  = 0._wp 
278         DO jk = 1, jpkm1
279            DO jj = 1, jpj
280               DO ji = 1, jpi
281                  z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
282               END DO
283            END DO
284         END DO
285         CALL iom_put( "saltc", rau0 * z2d )   ! vertically integrated salt content (PSU*kg/m2)
286      ENDIF
287      !
288      IF ( iom_use("eken") ) THEN
289         rke(:,:,jk) = 0._wp                               !      kinetic energy
290         DO jk = 1, jpkm1
291            DO jj = 2, jpjm1
292               DO ji = fs_2, fs_jpim1   ! vector opt.
293                  zztmp   = 1._wp / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
294                  zztmpx  = 0.5 * (  un(ji-1,jj,jk) * un(ji-1,jj,jk) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk)    &
295                     &             + un(ji  ,jj,jk) * un(ji  ,jj,jk) * e2u(ji  ,jj) * fse3u(ji  ,jj,jk) )  &
296                     &          *  zztmp 
297                  !
298                  zztmpy  = 0.5 * (  vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1v(ji,jj-1) * fse3v(ji,jj-1,jk)    &
299                     &             + vn(ji,jj  ,jk) * vn(ji,jj  ,jk) * e1v(ji,jj  ) * fse3v(ji,jj  ,jk) )  &
300                     &          *  zztmp 
301                  !
302                  rke(ji,jj,jk) = 0.5_wp * ( zztmpx + zztmpy )
303                  !
304               ENDDO
305            ENDDO
306         ENDDO
307         CALL lbc_lnk( rke, 'T', 1. )
308         CALL iom_put( "eken", rke )           
309      ENDIF
310         
311      IF( iom_use("u_masstr") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
312         z3d(:,:,jpk) = 0.e0
313         DO jk = 1, jpkm1
314            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) * umask(:,:,jk)
315         END DO
316         CALL iom_put( "u_masstr", z3d )                  ! mass transport in i-direction
317      ENDIF
318     
319      IF( iom_use("u_heattr") ) THEN
320         z2d(:,:) = 0.e0 
321         DO jk = 1, jpkm1
322            DO jj = 2, jpjm1
323               DO ji = fs_2, fs_jpim1   ! vector opt.
324                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
325               END DO
326            END DO
327         END DO
328         CALL lbc_lnk( z2d, 'U', -1. )
329         CALL iom_put( "u_heattr", (0.5 * rcp) * z2d )    ! heat transport in i-direction
330      ENDIF
331
332      IF( iom_use("u_salttr") ) THEN
333         z2d(:,:) = 0.e0 
334         DO jk = 1, jpkm1
335            DO jj = 2, jpjm1
336               DO ji = fs_2, fs_jpim1   ! vector opt.
337                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
338               END DO
339            END DO
340         END DO
341         CALL lbc_lnk( z2d, 'U', -1. )
342         CALL iom_put( "u_salttr", 0.5 * z2d )            ! heat transport in i-direction
343      ENDIF
344
345     
346      IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
347         z3d(:,:,jpk) = 0.e0
348         DO jk = 1, jpkm1
349            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) * vmask(:,:,jk)
350         END DO
351         CALL iom_put( "v_masstr", z3d )                  ! mass transport in j-direction
352      ENDIF
353     
354      IF( iom_use("v_heattr") ) THEN
355         z2d(:,:) = 0.e0 
356         DO jk = 1, jpkm1
357            DO jj = 2, jpjm1
358               DO ji = fs_2, fs_jpim1   ! vector opt.
359                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
360               END DO
361            END DO
362         END DO
363         CALL lbc_lnk( z2d, 'V', -1. )
364         CALL iom_put( "v_heattr", (0.5 * rcp) * z2d )    !  heat transport in j-direction
365      ENDIF
366
367      IF( iom_use("v_salttr") ) THEN
368         z2d(:,:) = 0.e0 
369         DO jk = 1, jpkm1
370            DO jj = 2, jpjm1
371               DO ji = fs_2, fs_jpim1   ! vector opt.
372                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
373               END DO
374            END DO
375         END DO
376         CALL lbc_lnk( z2d, 'V', -1. )
377         CALL iom_put( "v_salttr", 0.5 * z2d )            !  heat transport in j-direction
378      ENDIF
379      !
380      CALL wrk_dealloc( jpi , jpj      , z2d )
381      CALL wrk_dealloc( jpi , jpj, jpk , z3d )
382      !
383      IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
384      !
385   END SUBROUTINE dia_wri
386
387#else
388   !!----------------------------------------------------------------------
389   !!   Default option                                  use IOIPSL  library
390   !!----------------------------------------------------------------------
391
392   SUBROUTINE dia_wri( kt )
393      !!---------------------------------------------------------------------
394      !!                  ***  ROUTINE dia_wri  ***
395      !!                   
396      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
397      !!      NETCDF format is used by default
398      !!
399      !! ** Method  :   At the beginning of the first time step (nit000),
400      !!      define all the NETCDF files and fields
401      !!      At each time step call histdef to compute the mean if ncessary
402      !!      Each nwrite time step, output the instantaneous or mean fields
403      !!----------------------------------------------------------------------
404      !!
405      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
406      !!
407      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
408      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
409      INTEGER  ::   inum = 11                                ! temporary logical unit
410      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
411      INTEGER  ::   ierr                                     ! error code return from allocation
412      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
413      INTEGER  ::   jn, ierror                               ! local integers
414      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt           ! local scalars
415      !!
416      REAL(wp), POINTER, DIMENSION(:,:)   :: zw2d       ! 2D workspace
417      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d       ! 3D workspace
418      !!----------------------------------------------------------------------
419      !
420      IF( nn_timing == 1 )   CALL timing_start('dia_wri')
421      !
422      CALL wrk_alloc( jpi , jpj      , zw2d )
423      IF ( ln_traldf_gdia .OR. lk_vvl )  call wrk_alloc( jpi , jpj , jpk  , zw3d )
424      !
425      ! Output the initial state and forcings
426      IF( ninist == 1 ) THEN                       
427         CALL dia_wri_state( 'output.init', kt )
428         ninist = 0
429      ENDIF
430      !
431      ! 0. Initialisation
432      ! -----------------
433
434      ! local variable for debugging
435      ll_print = .FALSE.
436      ll_print = ll_print .AND. lwp
437
438      ! Define frequency of output and means
439      zdt = rdt
440      IF( nacc == 1 ) zdt = rdtmin
441      clop = "x"         ! no use of the mask value (require less cpu time, and otherwise the model crashes)
442#if defined key_diainstant
443      zsto = nwrite * zdt
444      clop = "inst("//TRIM(clop)//")"
445#else
446      zsto=zdt
447      clop = "ave("//TRIM(clop)//")"
448#endif
449      zout = nwrite * zdt
450      zmax = ( nitend - nit000 + 1 ) * zdt
451
452      ! Define indices of the horizontal output zoom and vertical limit storage
453      iimi = 1      ;      iima = jpi
454      ijmi = 1      ;      ijma = jpj
455      ipk = jpk
456
457      ! define time axis
458      it = kt
459      itmod = kt - nit000 + 1
460
461
462      ! 1. Define NETCDF files and fields at beginning of first time step
463      ! -----------------------------------------------------------------
464
465      IF( kt == nit000 ) THEN
466
467         ! Define the NETCDF files (one per grid)
468
469         ! Compute julian date from starting date of the run
470         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
471         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
472         IF(lwp)WRITE(numout,*)
473         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
474            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
475         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
476                                 ' limit storage in depth = ', ipk
477
478         ! WRITE root name in date.file for use by postpro
479         IF(lwp) THEN
480            CALL dia_nam( clhstnam, nwrite,' ' )
481            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
482            WRITE(inum,*) clhstnam
483            CLOSE(inum)
484         ENDIF
485
486         ! Define the T grid FILE ( nid_T )
487
488         CALL dia_nam( clhstnam, nwrite, 'grid_T' )
489         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
490         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
491            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
492            &          nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
493         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
494            &           "m", ipk, gdept_1d, nz_T, "down" )
495         !                                                            ! Index of ocean points
496         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
497         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
498         !
499         IF( ln_icebergs ) THEN
500            !
501            !! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
502            !! that routine is called from nemogcm, so do it here immediately before its needed
503            ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
504            IF( lk_mpp )   CALL mpp_sum( ierror )
505            IF( ierror /= 0 ) THEN
506               CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
507               RETURN
508            ENDIF
509            !
510            !! iceberg vertical coordinate is class number
511            CALL histvert( nid_T, "class", "Iceberg class",      &  ! Vertical grid: class
512               &           "number", nclasses, class_num, nb_T )
513            !
514            !! each class just needs the surface index pattern
515            ndim_bT = 3
516            DO jn = 1,nclasses
517               ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
518            ENDDO
519            !
520         ENDIF
521
522         ! Define the U grid FILE ( nid_U )
523
524         CALL dia_nam( clhstnam, nwrite, 'grid_U' )
525         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
526         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
527            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
528            &          nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
529         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
530            &           "m", ipk, gdept_1d, nz_U, "down" )
531         !                                                            ! Index of ocean points
532         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
533         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
534
535         ! Define the V grid FILE ( nid_V )
536
537         CALL dia_nam( clhstnam, nwrite, 'grid_V' )                   ! filename
538         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
539         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
540            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
541            &          nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
542         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
543            &          "m", ipk, gdept_1d, nz_V, "down" )
544         !                                                            ! Index of ocean points
545         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
546         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
547
548         ! Define the W grid FILE ( nid_W )
549
550         CALL dia_nam( clhstnam, nwrite, 'grid_W' )                   ! filename
551         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
552         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
553            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
554            &          nit000-1, zjulian, zdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
555         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
556            &          "m", ipk, gdepw_1d, nz_W, "down" )
557
558
559         ! Declare all the output fields as NETCDF variables
560
561         !                                                                                      !!! nid_T : 3D
562         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
563            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
564         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
565            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
566         IF(  lk_vvl  ) THEN
567            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n
568            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
569            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t_n
570            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
571            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t_n
572            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
573         ENDIF
574         !                                                                                      !!! nid_T : 2D
575         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
576            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
577         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
578            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
579         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
580            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
581         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
582            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
583         CALL histdef( nid_T, "sorunoff", "River runoffs"                      , "Kg/m2/s",   &  ! runoffs
584            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
585         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
586            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
587         IF(  .NOT. lk_vvl  ) THEN
588            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem)
589            &                                                                  , "KgC/m2/s",  &  ! sosst_cd
590            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
591            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * tsn(:,:,1,jp_sal)
592            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
593            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
594         ENDIF
595         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
596            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
597         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
598            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
599         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
600            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
601         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
602            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
603         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
604            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
605         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
606            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
607!
608         IF( ln_icebergs ) THEN
609            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
610               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
611            CALL histdef( nid_T, "calving_heat"        , "calving heat flux"                        , "XXXX"   , &
612               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
613            CALL histdef( nid_T, "berg_floating_melt"  , "Melt rate of icebergs + bits"             , "kg/m2/s", &
614               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
615            CALL histdef( nid_T, "berg_stored_ice"     , "Accumulated ice mass by class"            , "kg"     , &
616               &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
617            IF( ln_bergdia ) THEN
618               CALL histdef( nid_T, "berg_melt"           , "Melt rate of icebergs"                    , "kg/m2/s", &
619                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
620               CALL histdef( nid_T, "berg_buoy_melt"      , "Buoyancy component of iceberg melt rate"  , "kg/m2/s", &
621                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
622               CALL histdef( nid_T, "berg_eros_melt"      , "Erosion component of iceberg melt rate"   , "kg/m2/s", &
623                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
624               CALL histdef( nid_T, "berg_conv_melt"      , "Convective component of iceberg melt rate", "kg/m2/s", &
625                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
626               CALL histdef( nid_T, "berg_virtual_area"   , "Virtual coverage by icebergs"             , "m2"     , &
627                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
628               CALL histdef( nid_T, "bits_src"           , "Mass source of bergy bits"                , "kg/m2/s", &
629                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
630               CALL histdef( nid_T, "bits_melt"          , "Melt rate of bergy bits"                  , "kg/m2/s", &
631                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
632               CALL histdef( nid_T, "bits_mass"          , "Bergy bit density field"                  , "kg/m2"  , &
633                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
634               CALL histdef( nid_T, "berg_mass"           , "Iceberg density field"                    , "kg/m2"  , &
635                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
636               CALL histdef( nid_T, "berg_real_calving"   , "Calving into iceberg class"               , "kg/s"   , &
637                  &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
638            ENDIF
639         ENDIF
640
641         IF( .NOT. ln_cpl ) THEN
642            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
643               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
644            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
645               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
646            CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
647               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
648         ENDIF
649
650         IF( ln_cpl .AND. nn_ice <= 1 ) THEN
651            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
652               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
653            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
654               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
655            CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn
656               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
657         ENDIF
658         
659         clmx ="l_max(only(x))"    ! max index on a period
660         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
661            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
662#if defined key_diahth
663         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
664            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
665         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
666            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
667         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
668            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
669         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "W"      ,   & ! htc3
670            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
671#endif
672
673         IF( ln_cpl .AND. nn_ice == 2 ) THEN
674            CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature"            , "K"      ,   &  ! tn_ice
675               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
676            CALL histdef( nid_T,"soicealb" , "Ice Albedo"                         , "[0,1]"  ,   &  ! alb_ice
677               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
678         ENDIF
679
680         CALL histend( nid_T, snc4chunks=snc4set )
681
682         !                                                                                      !!! nid_U : 3D
683         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
684            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
685         IF( ln_wave .AND. ln_sdw) THEN
686            CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current"         , "m/s"    ,   &  ! usd
687               &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout ) 
688         ENDIF
689         IF( ln_traldf_gdia ) THEN
690            CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
691                 &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
692         ELSE
693#if defined key_diaeiv
694            CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
695            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
696#endif
697         END IF
698         !                                                                                      !!! nid_U : 2D
699         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
700            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
701
702         CALL histend( nid_U, snc4chunks=snc4set )
703
704         !                                                                                      !!! nid_V : 3D
705         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
706            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
707         IF( ln_wave .AND. ln_sdw) THEN
708            CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current"    , "m/s"    ,   &  ! vsd
709               &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout ) 
710         ENDIF
711         IF( ln_traldf_gdia ) THEN
712            CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
713                 &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
714         ELSE 
715#if defined key_diaeiv
716            CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
717            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
718#endif
719         END IF
720         !                                                                                      !!! nid_V : 2D
721         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
722            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
723
724         CALL histend( nid_V, snc4chunks=snc4set )
725
726         !                                                                                      !!! nid_W : 3D
727         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn
728            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
729         IF( ln_traldf_gdia ) THEN
730            CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
731                 &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
732         ELSE
733#if defined key_diaeiv
734            CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
735                 &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
736#endif
737         END IF
738         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
739            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
740         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avmu
741            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
742
743         IF( lk_zdfddm ) THEN
744            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
745               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
746         ENDIF
747         
748         IF( ln_wave .AND. ln_sdw) THEN
749            CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current"   , "m/s"    ,   &  ! wsd
750               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout ) 
751         ENDIF 
752         !                                                                                      !!! nid_W : 2D
753#if defined key_traldf_c2d
754         CALL histdef( nid_W, "soleahtw", "lateral eddy diffusivity"           , "m2/s"   ,   &  ! ahtw
755            &          jpi, jpj, nh_W, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
756# if defined key_traldf_eiv 
757            CALL histdef( nid_W, "soleaeiw", "eddy induced vel. coeff. at w-point", "m2/s",   &  ! aeiw
758               &       jpi, jpj, nh_W, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
759# endif
760#endif
761
762         CALL histend( nid_W, snc4chunks=snc4set )
763
764         IF(lwp) WRITE(numout,*)
765         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
766         IF(ll_print) CALL FLUSH(numout )
767
768      ENDIF
769
770      ! 2. Start writing data
771      ! ---------------------
772
773      ! ndex(1) est utilise ssi l'avant dernier argument est different de
774      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
775      ! donne le nombre d'elements, et ndex la liste des indices a sortir
776
777      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
778         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
779         WRITE(numout,*) '~~~~~~ '
780      ENDIF
781
782      IF( lk_vvl ) THEN
783         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
784         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
785         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
786         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
787      ELSE
788         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature
789         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity
790         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature
791         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
792      ENDIF
793      IF( lk_vvl ) THEN
794         zw3d(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
795         CALL histwrite( nid_T, "vovvle3t", it, fse3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
796         CALL histwrite( nid_T, "vovvldep", it, fsdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
797         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
798      ENDIF
799      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
800      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
801      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
802      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
803                                                                                  ! (includes virtual salt flux beneath ice
804                                                                                  ! in linear free surface case)
805      IF( .NOT. lk_vvl ) THEN
806         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
807         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
808         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
809         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
810      ENDIF
811      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
812      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
813      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
814      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
815      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
816      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
817!
818      IF( ln_icebergs ) THEN
819         !
820         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT ) 
821         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
822         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT ) 
823         !
824         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
825         !
826         IF( ln_bergdia ) THEN
827            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   ) 
828            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   ) 
829            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   ) 
830            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   ) 
831            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   ) 
832            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   ) 
833            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   ) 
834            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   ) 
835            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   ) 
836            !
837            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
838         ENDIF
839      ENDIF
840
841      IF( .NOT. ln_cpl ) THEN
842         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
843         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
844         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
845         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
846      ENDIF
847      IF( ln_cpl .AND. nn_ice <= 1 ) THEN
848         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
849         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
850         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
851         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
852      ENDIF
853!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
854!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
855
856#if defined key_diahth
857      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
858      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
859      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
860      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
861#endif
862
863      IF( ln_cpl .AND. nn_ice == 2 ) THEN
864         CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature
865         CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo
866      ENDIF
867
868      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
869      IF( ln_traldf_gdia ) THEN
870         IF (.not. ALLOCATED(psix_eiv))THEN
871            ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr )
872            IF( lk_mpp   )   CALL mpp_sum ( ierr )
873            IF( ierr > 0 )   CALL ctl_stop('STOP', 'diawri: unable to allocate psi{x,y}_eiv')
874            psix_eiv(:,:,:) = 0.0_wp
875            psiy_eiv(:,:,:) = 0.0_wp
876         ENDIF
877         DO jk=1,jpkm1
878            zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk)  ! u_eiv = -dpsix/dz
879         END DO
880         zw3d(:,:,jpk) = 0._wp
881         CALL histwrite( nid_U, "vozoeivu", it, zw3d, ndim_U , ndex_U )           ! i-eiv current
882      ELSE
883#if defined key_diaeiv
884         CALL histwrite( nid_U, "vozoeivu", it, u_eiv, ndim_U , ndex_U )          ! i-eiv current
885#endif
886      ENDIF
887      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
888
889      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
890      IF( ln_traldf_gdia ) THEN
891         DO jk=1,jpk-1
892            zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk)  ! v_eiv = -dpsiy/dz
893         END DO
894         zw3d(:,:,jpk) = 0._wp
895         CALL histwrite( nid_V, "vomeeivv", it, zw3d, ndim_V , ndex_V )           ! j-eiv current
896      ELSE
897#if defined key_diaeiv
898         CALL histwrite( nid_V, "vomeeivv", it, v_eiv, ndim_V , ndex_V )          ! j-eiv current
899#endif
900      ENDIF
901      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
902
903      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
904      IF( ln_traldf_gdia ) THEN
905         DO jk=1,jpk-1
906            DO jj = 2, jpjm1
907               DO ji = fs_2, fs_jpim1  ! vector opt.
908                  zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2v(ji,jj) + &
909                       &    (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1u(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx
910               END DO
911            END DO
912         END DO
913         zw3d(:,:,jpk) = 0._wp
914         CALL histwrite( nid_W, "voveeivw", it, zw3d          , ndim_T, ndex_T )    ! vert. eiv current
915      ELSE
916#   if defined key_diaeiv
917         CALL histwrite( nid_W, "voveeivw", it, w_eiv          , ndim_T, ndex_T )    ! vert. eiv current
918#   endif
919      ENDIF
920      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
921      CALL histwrite( nid_W, "votkeavm", it, avmu           , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
922      IF( lk_zdfddm ) THEN
923         CALL histwrite( nid_W, "voddmavs", it, fsavs(:,:,:), ndim_T, ndex_T )    ! S vert. eddy diff. coef.
924      ENDIF
925#if defined key_traldf_c2d
926      CALL histwrite( nid_W, "soleahtw", it, ahtw          , ndim_hT, ndex_hT )   ! lateral eddy diff. coef.
927# if defined key_traldf_eiv
928         CALL histwrite( nid_W, "soleaeiw", it, aeiw       , ndim_hT, ndex_hT )   ! EIV coefficient at w-point
929# endif
930#endif
931
932      IF( ln_wave .AND. ln_sdw ) THEN
933         CALL histwrite( nid_U, "sdzocrtx", it, usd           , ndim_U , ndex_U )    ! i-StokesDrift-current
934         CALL histwrite( nid_V, "sdmecrty", it, vsd           , ndim_V , ndex_V )    ! j-StokesDrift-current
935         CALL histwrite( nid_W, "sdvecrtz", it, wsd           , ndim_T , ndex_T )    ! StokesDrift vert. current
936      ENDIF 
937     
938      ! 3. Close all files
939      ! ---------------------------------------
940      IF( kt == nitend ) THEN
941         CALL histclo( nid_T )
942         CALL histclo( nid_U )
943         CALL histclo( nid_V )
944         CALL histclo( nid_W )
945      ENDIF
946      !
947      CALL wrk_dealloc( jpi , jpj      , zw2d )
948      IF ( ln_traldf_gdia .OR. lk_vvl )  call wrk_dealloc( jpi , jpj , jpk  , zw3d )
949      !
950      IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
951      !
952   END SUBROUTINE dia_wri
953# endif
954
955#endif
956
957   SUBROUTINE dia_wri_state( cdfile_name, kt )
958      !!---------------------------------------------------------------------
959      !!                 ***  ROUTINE dia_wri_state  ***
960      !!       
961      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
962      !!      the instantaneous ocean state and forcing fields.
963      !!        Used to find errors in the initial state or save the last
964      !!      ocean state in case of abnormal end of a simulation
965      !!
966      !! ** Method  :   NetCDF files using ioipsl
967      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
968      !!      File 'output.abort.nc' is created in case of abnormal job end
969      !!----------------------------------------------------------------------
970      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
971      INTEGER           , INTENT( in ) ::   kt               ! ocean time-step index
972      !!
973      CHARACTER (len=32) :: clname
974      CHARACTER (len=40) :: clop
975      INTEGER  ::   id_i , nz_i, nh_i       
976      INTEGER, DIMENSION(1) ::   idex             ! local workspace
977      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt
978      !!----------------------------------------------------------------------
979      !
980!     IF( nn_timing == 1 )   CALL timing_start('dia_wri_state') ! not sure this works for routines not called in first timestep
981
982      ! 0. Initialisation
983      ! -----------------
984
985      ! Define name, frequency of output and means
986      clname = cdfile_name
987      IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
988      zdt  = rdt
989      zsto = rdt
990      clop = "inst(x)"           ! no use of the mask value (require less cpu time)
991      zout = rdt
992      zmax = ( nitend - nit000 + 1 ) * zdt
993
994      IF(lwp) WRITE(numout,*)
995      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
996      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
997      IF(lwp) WRITE(numout,*) '                and named :', clname, '.nc'
998
999
1000      ! 1. Define NETCDF files and fields at beginning of first time step
1001      ! -----------------------------------------------------------------
1002
1003      ! Compute julian date from starting date of the run
1004      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )         ! time axis
1005      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
1006      CALL histbeg( clname, jpi, glamt, jpj, gphit,   &
1007          1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
1008      CALL histvert( id_i, "deptht", "Vertical T levels",   &    ! Vertical grid : gdept
1009          "m", jpk, gdept_1d, nz_i, "down")
1010
1011      ! Declare all the output fields as NetCDF variables
1012
1013      CALL histdef( id_i, "vosaline", "Salinity"              , "PSU"    ,   &   ! salinity
1014         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1015      CALL histdef( id_i, "votemper", "Temperature"           , "C"      ,   &   ! temperature
1016         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1017      CALL histdef( id_i, "sossheig", "Sea Surface Height"    , "m"      ,   &  ! ssh
1018         &          jpi, jpj, nh_i, 1  , 1, 1  , nz_i, 32, clop, zsto, zout )
1019      CALL histdef( id_i, "vozocrtx", "Zonal Current"         , "m/s"    ,   &   ! zonal current
1020         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1021      CALL histdef( id_i, "vomecrty", "Meridional Current"    , "m/s"    ,   &   ! meridonal current
1022         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
1023      CALL histdef( id_i, "vovecrtz", "Vertical Velocity"     , "m/s"    ,   &   ! vertical current
1024         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
1025      CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S",   &   ! net freshwater
1026         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1027      CALL histdef( id_i, "sohefldo", "Net Downward Heat Flux", "W/m2"   ,   &   ! net heat flux
1028         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1029      CALL histdef( id_i, "soshfldo", "Shortwave Radiation"   , "W/m2"   ,   &   ! solar flux
1030         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1031      CALL histdef( id_i, "soicecov", "Ice fraction"          , "[0,1]"  ,   &   ! fr_i
1032         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1033      CALL histdef( id_i, "sozotaux", "Zonal Wind Stress"     , "N/m2"   ,   &   ! i-wind stress
1034         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1035      CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2"   ,   &   ! j-wind stress
1036         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1037      IF( lk_vvl ) THEN
1038         CALL histdef( id_i, "vovvldep", "T point depth"         , "m"      ,   &   ! t-point depth
1039            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1040         CALL histdef( id_i, "vovvle3t", "T point thickness"         , "m"      ,   &   ! t-point depth
1041            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1042      END IF
1043      !
1044      IF( ln_wave .AND. ln_sdw ) THEN
1045         CALL histdef( id_i, "sdzocrtx", "Stokes Drift Zonal", "m/s"    , &   ! StokesDrift zonal current
1046            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
1047         CALL histdef( id_i, "sdmecrty", "Stokes Drift Merid", "m/s"    , &   ! StokesDrift meridonal current
1048            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
1049         CALL histdef( id_i, "sdvecrtz", "Stokes Drift Vert", "m/s"    , &   ! StokesDrift vertical current
1050            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
1051      ENDIF 
1052
1053#if defined key_lim2
1054      CALL lim_wri_state_2( kt, id_i, nh_i )
1055#elif defined key_lim3
1056      CALL lim_wri_state( kt, id_i, nh_i )
1057#else
1058      CALL histend( id_i, snc4chunks=snc4set )
1059#endif
1060
1061      ! 2. Start writing data
1062      ! ---------------------
1063      ! idex(1) est utilise ssi l'avant dernier argument est diffferent de
1064      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
1065      ! donne le nombre d'elements, et idex la liste des indices a sortir
1066      idex(1) = 1   ! init to avoid compil warning
1067
1068      ! Write all fields on T grid
1069      CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature
1070      CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity
1071      CALL histwrite( id_i, "sossheig", kt, sshn             , jpi*jpj    , idex )    ! sea surface height
1072      CALL histwrite( id_i, "vozocrtx", kt, un               , jpi*jpj*jpk, idex )    ! now i-velocity
1073      CALL histwrite( id_i, "vomecrty", kt, vn               , jpi*jpj*jpk, idex )    ! now j-velocity
1074      CALL histwrite( id_i, "vovecrtz", kt, wn               , jpi*jpj*jpk, idex )    ! now k-velocity
1075      CALL histwrite( id_i, "sowaflup", kt, (emp-rnf )       , jpi*jpj    , idex )    ! freshwater budget
1076      CALL histwrite( id_i, "sohefldo", kt, qsr + qns        , jpi*jpj    , idex )    ! total heat flux
1077      CALL histwrite( id_i, "soshfldo", kt, qsr              , jpi*jpj    , idex )    ! solar heat flux
1078      CALL histwrite( id_i, "soicecov", kt, fr_i             , jpi*jpj    , idex )    ! ice fraction
1079      CALL histwrite( id_i, "sozotaux", kt, utau             , jpi*jpj    , idex )    ! i-wind stress
1080      CALL histwrite( id_i, "sometauy", kt, vtau             , jpi*jpj    , idex )    ! j-wind stress
1081      IF( lk_vvl ) THEN
1082         CALL histwrite( id_i, "vovvldep", kt, fsdept_n(:,:,:), jpi*jpj*jpk, idex )!  T-cell depth       
1083         CALL histwrite( id_i, "vovvle3t", kt, fse3t_n (:,:,:), jpi*jpj*jpk, idex )!  T-cell thickness 
1084      END IF
1085     
1086      IF( ln_wave .AND. ln_sdw ) THEN
1087         CALL histwrite( id_i, "sdzocrtx", kt, usd           , jpi*jpj*jpk, idex)     ! now StokesDrift i-velocity
1088         CALL histwrite( id_i, "sdmecrty", kt, vsd           , jpi*jpj*jpk, idex)     ! now StokesDrift j-velocity
1089         CALL histwrite( id_i, "sdvecrtz", kt, wsd           , jpi*jpj*jpk, idex)     ! now StokesDrift k-velocity
1090      ENDIF 
1091
1092      ! 3. Close the file
1093      ! -----------------
1094      CALL histclo( id_i )
1095#if ! defined key_iomput && ! defined key_dimgout
1096      IF( ninist /= 1  ) THEN
1097         CALL histclo( nid_T )
1098         CALL histclo( nid_U )
1099         CALL histclo( nid_V )
1100         CALL histclo( nid_W )
1101      ENDIF
1102#endif
1103       
1104!     IF( nn_timing == 1 )   CALL timing_stop('dia_wri_state') ! not sure this works for routines not called in first timestep
1105      !
1106
1107   END SUBROUTINE dia_wri_state
1108   !!======================================================================
1109END MODULE diawri
Note: See TracBrowser for help on using the repository browser.