source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 19 months ago

GMED 450 add flush after prints

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