source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 8400

Last change on this file since 8400 was 8400, checked in by timgraham, 4 years ago

GMED ticket 335:

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