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

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

source: branches/UKMO/dev_r5518_GO6_package_XIOS_read/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 8161

Last change on this file since 8161 was 8161, checked in by andmirek, 7 years ago

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