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/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/UKMO/AMM15_v3_6_STABLE_package_collate_coupling/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 10478

Last change on this file since 10478 was 10478, checked in by jcastill, 6 years ago

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