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 NEMO/trunk/src/OCE/DIA – NEMO

source: NEMO/trunk/src/OCE/DIA/diawri.F90 @ 15048

Last change on this file since 15048 was 15017, checked in by smasson, 3 years ago

trunk: introduce iom_put without halos, #2694
  • Property svn:keywords set to Id
File size: 63.9 KB
Line 
1MODULE diawri
2   !!======================================================================
3   !!                     ***  MODULE  diawri  ***
4   !! Ocean diagnostics :  write ocean output files
5   !!=====================================================================
6   !! History :  OPA  ! 1991-03  (M.-A. Foujols)  Original code
7   !!            4.0  ! 1991-11  (G. Madec)
8   !!                 ! 1992-06  (M. Imbard)  correction restart file
9   !!                 ! 1992-07  (M. Imbard)  split into diawri and rstwri
10   !!                 ! 1993-03  (M. Imbard)  suppress writibm
11   !!                 ! 1998-01  (C. Levy)  NETCDF format using ioipsl INTERFACE
12   !!                 ! 1999-02  (E. Guilyardi)  name of netCDF files + variables
13   !!            8.2  ! 2000-06  (M. Imbard)  Original code (diabort.F)
14   !!   NEMO     1.0  ! 2002-06  (A.Bozec, E. Durand)  Original code (diainit.F)
15   !!             -   ! 2002-09  (G. Madec)  F90: Free form and module
16   !!             -   ! 2002-12  (G. Madec)  merge of diabort and diainit, F90
17   !!                 ! 2005-11  (V. Garnier) Surface pressure gradient organization
18   !!            3.2  ! 2008-11  (B. Lemaire) creation from old diawri
19   !!            3.7  ! 2014-01  (G. Madec) remove eddy induced velocity from no-IOM output
20   !!                 !                     change name of output variables in dia_wri_state
21   !!            4.0  ! 2020-10  (A. Nasser, S. Techene) add diagnostic for SWE
22   !!----------------------------------------------------------------------
23
24   !!----------------------------------------------------------------------
25   !!   dia_wri       : create the standart output files
26   !!   dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
27   !!----------------------------------------------------------------------
28   USE oce            ! ocean dynamics and tracers
29   USE isf_oce
30   USE isfcpl
31   USE abl            ! abl variables in case ln_abl = .true.
32   USE dom_oce        ! ocean space and time domain
33   USE phycst         ! physical constants
34   USE dianam         ! build name of file (routine)
35   USE diahth         ! thermocline diagnostics
36   USE dynadv   , ONLY: ln_dynadv_vec
37   USE icb_oce        ! Icebergs
38   USE icbdia         ! Iceberg budgets
39   USE ldftra         ! lateral physics: eddy diffusivity coef.
40   USE ldfdyn         ! lateral physics: eddy viscosity   coef.
41   USE sbc_oce        ! Surface boundary condition: ocean fields
42   USE sbc_ice        ! Surface boundary condition: ice fields
43   USE sbcssr         ! restoring term toward SST/SSS climatology
44   USE sbcwave        ! wave parameters
45   USE wet_dry        ! wetting and drying
46   USE zdf_oce        ! ocean vertical physics
47   USE zdfdrg         ! ocean vertical physics: top/bottom friction
48   USE zdfmxl         ! mixed layer
49   USE zdfosm         ! mixed layer
50   !
51   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
52   USE in_out_manager ! I/O manager
53   USE dia25h         ! 25h Mean output
54   USE iom            !
55   USE ioipsl         !
56
57#if defined key_si3
58   USE ice 
59   USE icewri 
60#endif
61   USE lib_mpp         ! MPP library
62   USE timing          ! preformance summary
63   USE diu_bulk        ! diurnal warm layer
64   USE diu_coolskin    ! Cool skin
65
66   IMPLICIT NONE
67   PRIVATE
68
69   PUBLIC   dia_wri                 ! routines called by step.F90
70   PUBLIC   dia_wri_state
71   PUBLIC   dia_wri_alloc           ! Called by nemogcm module
72#if ! defined key_xios   
73   PUBLIC   dia_wri_alloc_abl       ! Called by sbcabl  module (if ln_abl = .true.)
74#endif
75   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
76   INTEGER ::          nb_T              , ndim_bT   ! grid_T file
77   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
78   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
79   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file
80   INTEGER ::   nid_A, nz_A, nh_A, ndim_A, ndim_hA   ! grid_ABL file   
81   INTEGER ::   ndex(1)                              ! ???
82   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
83   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hA, ndex_A ! ABL
84   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
85   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
86
87   !! * Substitutions
88#  include "do_loop_substitute.h90"
89#  include "domzgr_substitute.h90"
90   !!----------------------------------------------------------------------
91   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
92   !! $Id$
93   !! Software governed by the CeCILL license (see ./LICENSE)
94   !!----------------------------------------------------------------------
95CONTAINS
96
97#if defined key_xios
98   !!----------------------------------------------------------------------
99   !!   'key_xios'                                        use IOM library
100   !!----------------------------------------------------------------------
101   INTEGER FUNCTION dia_wri_alloc()
102      !
103      dia_wri_alloc = 0
104      !
105   END FUNCTION dia_wri_alloc
106
107   
108   SUBROUTINE dia_wri( kt, Kmm )
109      !!---------------------------------------------------------------------
110      !!                  ***  ROUTINE dia_wri  ***
111      !!                   
112      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
113      !!      NETCDF format is used by default
114      !!
115      !! ** Method  :  use iom_put
116      !!----------------------------------------------------------------------
117      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
118      INTEGER, INTENT( in ) ::   Kmm     ! ocean time level index
119      !!
120      INTEGER ::   ji, jj, jk       ! dummy loop indices
121      INTEGER ::   ikbot            ! local integer
122      REAL(wp)::   zztmp , zztmpx   ! local scalar
123      REAL(wp)::   zztmp2, zztmpy   !   -      -
124      REAL(wp)::   ze3
125      REAL(wp), DIMENSION(A2D(     0))     ::   z2d   ! 2D workspace
126      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   z3d   ! 3D workspace
127      !!----------------------------------------------------------------------
128      !
129      IF( ln_timing )   CALL timing_start('dia_wri')
130      !
131      ! Output the initial state and forcings
132      IF( ninist == 1 ) THEN                       
133         CALL dia_wri_state( Kmm, 'output.init' )
134         ninist = 0
135      ENDIF
136
137      ! Output of initial vertical scale factor
138      CALL iom_put("e3t_0", e3t_0(:,:,:) )
139      CALL iom_put("e3u_0", e3u_0(:,:,:) )
140      CALL iom_put("e3v_0", e3v_0(:,:,:) )
141      CALL iom_put("e3f_0", e3f_0(:,:,:) )
142      !
143      IF ( iom_use("tpt_dep") ) THEN
144         DO jk = 1, jpk
145            z3d(:,:,jk) = gdept(:,:,jk,Kmm)
146         END DO
147         CALL iom_put( "tpt_dep", z3d )
148      ENDIF
149
150      IF ( iom_use("e3t") .OR. iom_use("e3tdef") ) THEN  ! time-varying e3t
151         DO jk = 1, jpk
152            z3d(:,:,jk) =  e3t(:,:,jk,Kmm)
153         END DO
154         CALL iom_put( "e3t", z3d )
155         IF ( iom_use("e3tdef") ) THEN
156            z3d(:,:,:) = ( ( z3d(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100._wp * tmask(:,:,:) ) ** 2 
157            CALL iom_put( "e3tdef", z3d ) 
158         ENDIF
159      ENDIF
160      IF ( iom_use("e3u") ) THEN                         ! time-varying e3u
161         DO jk = 1, jpk
162            z3d(:,:,jk) =  e3u(:,:,jk,Kmm)
163         END DO
164         CALL iom_put( "e3u" , z3d )
165      ENDIF
166      IF ( iom_use("e3v") ) THEN                         ! time-varying e3v
167         DO jk = 1, jpk
168            z3d(:,:,jk) =  e3v(:,:,jk,Kmm)
169         END DO
170         CALL iom_put( "e3v" , z3d )
171      ENDIF
172      IF ( iom_use("e3w") ) THEN                         ! time-varying e3w
173         DO jk = 1, jpk
174            z3d(:,:,jk) =  e3w(:,:,jk,Kmm)
175         END DO
176         CALL iom_put( "e3w" , z3d )
177      ENDIF
178      IF ( iom_use("e3f") ) THEN                         ! time-varying e3f caution here at Kaa
179          DO jk = 1, jpk
180            z3d(:,:,jk) =  e3f(:,:,jk)
181         END DO
182         CALL iom_put( "e3f" , z3d )
183      ENDIF
184
185      IF ( iom_use("ssh") ) THEN
186         IF( ll_wd ) THEN                                ! sea surface height (brought back to the reference used for wetting and drying)
187            CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*ssmask(:,:) )
188         ELSE
189            CALL iom_put( "ssh" ,  ssh(:,:,Kmm) )        ! sea surface height
190         ENDIF
191      ENDIF
192
193      IF( iom_use("wetdep") )    CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) )   ! wet depth
194         
195#if defined key_qco
196      IF( iom_use("ht") )   CALL iom_put( "ht" , ht(:,:)     )   ! water column at t-point
197      IF( iom_use("hu") )   CALL iom_put( "hu" , hu(:,:,Kmm) )   ! water column at u-point
198      IF( iom_use("hv") )   CALL iom_put( "hv" , hv(:,:,Kmm) )   ! water column at v-point
199      IF( iom_use("hf") )   CALL iom_put( "hf" , hf_0(:,:)*( 1._wp + r3f(:,:) ) )   ! water column at f-point (caution here at Naa)
200#endif
201     
202      CALL iom_put( "toce", ts(:,:,:,jp_tem,Kmm) )    ! 3D temperature
203      CALL iom_put(  "sst", ts(:,:,1,jp_tem,Kmm) )    ! surface temperature
204      IF ( iom_use("sbt") ) THEN
205         DO_2D( 0, 0, 0, 0 )
206            ikbot = mbkt(ji,jj)
207            z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm)
208         END_2D
209         CALL iom_put( "sbt", z2d )                ! bottom temperature
210      ENDIF
211     
212      CALL iom_put( "soce", ts(:,:,:,jp_sal,Kmm) )    ! 3D salinity
213      CALL iom_put(  "sss", ts(:,:,1,jp_sal,Kmm) )    ! surface salinity
214      IF ( iom_use("sbs") ) THEN
215         DO_2D( 0, 0, 0, 0 )
216            ikbot = mbkt(ji,jj)
217            z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm)
218         END_2D
219         CALL iom_put( "sbs", z2d )                ! bottom salinity
220      ENDIF
221
222      IF( .NOT.lk_SWE )   CALL iom_put( "rhop", rhop(:,:,:) )          ! 3D potential density (sigma0)
223
224      IF ( iom_use("taubot") ) THEN                ! bottom stress
225         zztmp = rho0 * 0.25_wp
226         z2d(:,:) = 0._wp
227         DO_2D( 0, 0, 0, 0 )
228            zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * uu(ji  ,jj,mbku(ji  ,jj),Kmm)  )**2   &
229               &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm)  )**2   &
230               &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vv(ji,jj  ,mbkv(ji,jj  ),Kmm)  )**2   &
231               &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm)  )**2
232            z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 
233            !
234         END_2D
235         CALL iom_put( "taubot", z2d )           
236      ENDIF
237         
238      CALL iom_put( "uoce", uu(:,:,:,Kmm) )            ! 3D i-current
239      CALL iom_put(  "ssu", uu(:,:,1,Kmm) )            ! surface i-current
240      IF ( iom_use("sbu") ) THEN
241         DO_2D( 0, 0, 0, 0 )
242            ikbot = mbku(ji,jj)
243            z2d(ji,jj) = uu(ji,jj,ikbot,Kmm)
244         END_2D
245         CALL iom_put( "sbu", z2d )                ! bottom i-current
246      ENDIF
247     
248      CALL iom_put( "voce", vv(:,:,:,Kmm) )            ! 3D j-current
249      CALL iom_put(  "ssv", vv(:,:,1,Kmm) )            ! surface j-current
250      IF ( iom_use("sbv") ) THEN
251         DO_2D( 0, 0, 0, 0 )
252            ikbot = mbkv(ji,jj)
253            z2d(ji,jj) = vv(ji,jj,ikbot,Kmm)
254         END_2D
255         CALL iom_put( "sbv", z2d )                ! bottom j-current
256      ENDIF
257
258      !                                            ! vertical velocity
259      IF( ln_zad_Aimp ) THEN   ;   IF( iom_use('woce') )   CALL iom_put( "woce", ww + wi )   ! explicit plus implicit parts
260      ELSE                     ;                           CALL iom_put( "woce", ww )
261      ENDIF
262
263      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
264         !                     ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
265         DO jk = 1, jpk
266            IF( ln_zad_Aimp ) THEN
267               z3d(:,:,jk) = rho0 * e1e2t(:,:) * ( ww(:,:,jk) + wi(:,:,jk) )
268            ELSE
269               z3d(:,:,jk) = rho0 * e1e2t(:,:) * ww(:,:,jk)
270            ENDIF
271         END DO
272         CALL iom_put( "w_masstr" , z3d ) 
273         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d * z3d )
274      ENDIF
275
276      CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef.
277      CALL iom_put( "avs" , avs )                  ! S vert. eddy diff. coef.
278      CALL iom_put( "avm" , avm )                  ! T vert. eddy visc. coef.
279
280      IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
281      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
282
283      IF ( iom_use("sssgrad") .OR. iom_use("sssgrad2") ) THEN
284         DO_2D( 0, 0, 0, 0 )                       ! sss gradient
285            zztmp  = ts(ji,jj,1,jp_sal,Kmm)
286            zztmpx = (ts(ji+1,jj,1,jp_sal,Kmm) - zztmp) * r1_e1u(ji,jj) + (zztmp - ts(ji-1,jj  ,1,jp_sal,Kmm)) * r1_e1u(ji-1,jj)
287            zztmpy = (ts(ji,jj+1,1,jp_sal,Kmm) - zztmp) * r1_e2v(ji,jj) + (zztmp - ts(ji  ,jj-1,1,jp_sal,Kmm)) * r1_e2v(ji,jj-1)
288            z2d(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
289               &                 * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * vmask(ji,jj-1,1)
290         END_2D
291         CALL iom_put( "sssgrad2",  z2d )          ! square of module of sss gradient
292         IF ( iom_use("sssgrad") ) THEN
293            DO_2D( 0, 0, 0, 0 )
294               z2d(ji,jj) = SQRT( z2d(ji,jj) )
295            END_2D
296            CALL iom_put( "sssgrad",  z2d )        ! module of sss gradient
297         ENDIF
298      ENDIF
299         
300      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
301         DO_2D( 0, 0, 0, 0 )                       ! sst gradient
302            zztmp  = ts(ji,jj,1,jp_tem,Kmm)
303            zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj  ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj)
304            zztmpy = ( ts(ji,jj+1,1,jp_tem,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji  ,jj-1,1,jp_tem,Kmm) ) * r1_e2v(ji,jj-1)
305            z2d(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
306               &                 * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * vmask(ji,jj-1,1)
307         END_2D
308         CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient
309         IF ( iom_use("sstgrad") ) THEN
310            DO_2D( 0, 0, 0, 0 )
311               z2d(ji,jj) = SQRT( z2d(ji,jj) )
312            END_2D
313            CALL iom_put( "sstgrad",  z2d )        ! module of sst gradient
314         ENDIF
315      ENDIF
316         
317      ! heat and salt contents
318      IF( iom_use("heatc") ) THEN
319         z2d(:,:)  = 0._wp 
320         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
321            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk)
322         END_3D
323         CALL iom_put( "heatc", rho0_rcp * z2d )   ! vertically integrated heat content (J/m2)
324      ENDIF
325
326      IF( iom_use("saltc") ) THEN
327         z2d(:,:)  = 0._wp 
328         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
329            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
330         END_3D
331         CALL iom_put( "saltc", rho0 * z2d )       ! vertically integrated salt content (PSU*kg/m2)
332      ENDIF
333      !
334      IF( iom_use("salt2c") ) THEN
335         z2d(:,:)  = 0._wp 
336         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
337            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
338         END_3D
339         CALL iom_put( "salt2c", rho0 * z2d )      ! vertically integrated salt content (PSU*kg/m2)
340      ENDIF
341      !
342      IF ( iom_use("ke") .OR. iom_use("ke_int") ) THEN
343         z3d(:,:,jpk) = 0._wp 
344         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
345            zztmpx = uu(ji-1,jj  ,jk,Kmm) + uu(ji,jj,jk,Kmm)
346            zztmpy = vv(ji  ,jj-1,jk,Kmm) + vv(ji,jj,jk,Kmm)
347            z3d(ji,jj,jk) = 0.25_wp * ( zztmpx*zztmpx + zztmpy*zztmpy )
348         END_3D
349         CALL iom_put( "ke", z3d )                 ! kinetic energy
350
351         z2d(:,:)  = 0._wp 
352         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
353            z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * e1e2t(ji,jj) * tmask(ji,jj,jk)
354         END_3D
355         CALL iom_put( "ke_int", z2d )             ! vertically integrated kinetic energy
356      ENDIF
357      !
358      IF ( iom_use("sKE") ) THEN                   ! surface kinetic energy at T point
359         z2d(:,:) = 0._wp
360         DO_2D( 0, 0, 0, 0 )
361            z2d(ji,jj) = 0.25_wp * ( uu(ji  ,jj,1,Kmm) * uu(ji  ,jj,1,Kmm) * e1e2u(ji  ,jj) * e3u(ji  ,jj,1,Kmm)  &
362               &                   + uu(ji-1,jj,1,Kmm) * uu(ji-1,jj,1,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,1,Kmm)  &
363               &                   + vv(ji,jj  ,1,Kmm) * vv(ji,jj  ,1,Kmm) * e1e2v(ji,jj  ) * e3v(ji,jj  ,1,Kmm)  & 
364               &                   + vv(ji,jj-1,1,Kmm) * vv(ji,jj-1,1,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,1,Kmm)  )  &
365               &                 * r1_e1e2t(ji,jj) / e3t(ji,jj,1,Kmm) * ssmask(ji,jj)
366         END_2D
367         IF ( iom_use("sKE" ) )  CALL iom_put( "sKE" , z2d )   
368      ENDIF
369      !   
370      IF ( iom_use("ssKEf") ) THEN                 ! surface kinetic energy at F point
371         z2d(:,:) = 0._wp                          ! CAUTION : only valid in SWE, not with bathymetry
372         DO_2D( 0, 0, 0, 0 )
373            z2d(ji,jj) = 0.25_wp * ( uu(ji,jj  ,1,Kmm) * uu(ji,jj  ,1,Kmm) * e1e2u(ji,jj  ) * e3u(ji,jj  ,1,Kmm)  &
374               &                   + uu(ji,jj+1,1,Kmm) * uu(ji,jj+1,1,Kmm) * e1e2u(ji,jj+1) * e3u(ji,jj+1,1,Kmm)  &
375               &                   + vv(ji  ,jj,1,Kmm) * vv(ji,jj  ,1,Kmm) * e1e2v(ji  ,jj) * e3v(ji  ,jj,1,Kmm)  & 
376               &                   + vv(ji+1,jj,1,Kmm) * vv(ji+1,jj,1,Kmm) * e1e2v(ji+1,jj) * e3v(ji+1,jj,1,Kmm)  )  &
377               &                 * r1_e1e2f(ji,jj) / e3f(ji,jj,1) * ssfmask(ji,jj)
378         END_2D
379         CALL iom_put( "ssKEf", z2d )                     
380      ENDIF
381      !
382      CALL iom_put( "hdiv", hdiv )                 ! Horizontal divergence
383      !
384      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
385         
386         z3d(:,:,jpk) = 0._wp
387         DO jk = 1, jpkm1
388            z3d(:,:,jk) = rho0 * uu(:,:,jk,Kmm) * e2u(:,:) * e3u(:,:,jk,Kmm) * umask(:,:,jk)
389         END DO
390         CALL iom_put( "u_masstr"     , z3d )      ! mass transport in i-direction
391         
392         IF( iom_use("u_masstr_vint") ) THEN
393            z2d(:,:) = 0._wp 
394            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
395               z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
396            END_3D
397            CALL iom_put( "u_masstr_vint", z2d )   ! mass transport in i-direction vertical sum
398         ENDIF
399         IF( iom_use("u_heattr") ) THEN
400            z2d(:,:) = 0._wp 
401            zztmp = 0.5_wp * rcp
402            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
403               z2d(ji,jj) = z2d(ji,jj) + zztmp * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) )
404            END_3D
405            CALL iom_put( "u_heattr", z2d )        ! heat transport in i-direction
406         ENDIF
407         IF( iom_use("u_salttr") ) THEN
408            z2d(:,:) = 0._wp 
409            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
410               z2d(ji,jj) = z2d(ji,jj) +   0.5 * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) )
411            END_3D
412            CALL iom_put( "u_salttr", z2d )        ! heat transport in i-direction
413         ENDIF
414         
415      ENDIF
416     
417      IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
418         
419         z3d(:,:,jpk) = 0._wp
420         DO jk = 1, jpkm1
421            z3d(:,:,jk) = rho0 * vv(:,:,jk,Kmm) * e1v(:,:) * e3v(:,:,jk,Kmm) * vmask(:,:,jk)
422         END DO
423         CALL iom_put( "v_masstr", z3d )           ! mass transport in j-direction
424         
425         IF( iom_use("v_heattr") ) THEN
426            z2d(:,:) = 0._wp
427            zztmp = 0.5_wp * rcp
428            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
429               z2d(ji,jj) = z2d(ji,jj) + zztmp * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) )
430            END_3D
431            CALL iom_put( "v_heattr", z2d )        !  heat transport in j-direction
432         ENDIF
433         IF( iom_use("v_salttr") ) THEN
434            z2d(:,:) = 0._wp 
435            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
436               z2d(ji,jj) = z2d(ji,jj) +   0.5 * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) )
437            END_3D
438            CALL iom_put( "v_salttr", z2d )        !  heat transport in j-direction
439         ENDIF
440
441      ENDIF
442
443      IF( iom_use("tosmint") ) THEN
444         z2d(:,:) = 0._wp
445         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
446            z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)
447         END_3D
448         CALL iom_put( "tosmint", z2d )            ! Vertical integral of temperature
449      ENDIF
450      IF( iom_use("somint") ) THEN
451         z2d(:,:) = 0._wp
452         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
453            z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
454         END_3D
455         CALL iom_put( "somint", z2d )             ! Vertical integral of salinity
456      ENDIF
457
458      CALL iom_put( "bn2", rn2 )                   ! Brunt-Vaisala buoyancy frequency (N^2)
459     
460      IF (ln_dia25h)   CALL dia_25h( kt, Kmm )     ! 25h averaging
461     
462      ! Output of surface vorticity terms
463      !
464      CALL iom_put( "ssplavor", ff_f )             ! planetary vorticity ( f )
465      !
466      IF ( iom_use("ssrelvor")    .OR. iom_use("ssEns")    .OR.   &
467         & iom_use("ssrelpotvor") .OR. iom_use("ssabspotvor") ) THEN
468         !
469         z2d(:,:) = 0._wp 
470         DO_2D( 0, 0, 0, 0 )
471            z2d(ji,jj) = (   e2v(ji+1,jj  ) * vv(ji+1,jj  ,1,Kmm) - e2v(ji,jj) * vv(ji,jj,1,Kmm)    &
472            &              - e1u(ji  ,jj+1) * uu(ji  ,jj+1,1,Kmm) + e1u(ji,jj) * uu(ji,jj,1,Kmm)  ) * r1_e1e2f(ji,jj)
473         END_2D
474         CALL iom_put( "ssrelvor", z2d )           ! relative vorticity ( zeta )
475         !
476         IF ( iom_use("ssEns") .OR. iom_use("ssrelpotvor") .OR. iom_use("ssabspotvor") ) THEN
477            DO_2D( 0, 0, 0, 0 ) 
478               ze3 = (  e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1)    &
479                  &    + e3t(ji,jj  ,1,Kmm) * e1e2t(ji,jj  ) + e3t(ji+1,jj  ,1,Kmm) * e1e2t(ji+1,jj  )  ) * r1_e1e2f(ji,jj)
480               IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3
481               ELSE                      ;   ze3 = 0._wp
482               ENDIF
483               z2d(ji,jj) = ze3 * z2d(ji,jj) 
484            END_2D
485            CALL iom_put( "ssrelpotvor", z2d )     ! relative potential vorticity (zeta/h)
486            !
487            IF ( iom_use("ssEns") .OR. iom_use("ssabspotvor") ) THEN
488               DO_2D( 0, 0, 0, 0 )
489                  ze3 = (  e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1)    &
490                     &    + e3t(ji,jj  ,1,Kmm) * e1e2t(ji,jj  ) + e3t(ji+1,jj  ,1,Kmm) * e1e2t(ji+1,jj  )  ) * r1_e1e2f(ji,jj)
491                  IF( ze3 /= 0._wp ) THEN   ;   ze3 = 4._wp / ze3
492                  ELSE                      ;   ze3 = 0._wp
493                  ENDIF
494                  z2d(ji,jj) = ze3 * ff_f(ji,jj) + z2d(ji,jj) 
495               END_2D
496               CALL iom_put( "ssabspotvor", z2d )  ! absolute potential vorticity ( q )
497               !
498               IF ( iom_use("ssEns") ) THEN
499                  DO_2D( 0, 0, 0, 0 ) 
500                     z2d(ji,jj) = 0.5_wp * z2d(ji,jj) * z2d(ji,jj) 
501                  END_2D
502                  CALL iom_put( "ssEns", z2d )     ! potential enstrophy ( 1/2*q2 )
503               ENDIF
504            ENDIF
505         ENDIF
506      ENDIF
507
508      IF( ln_timing )   CALL timing_stop('dia_wri')
509      !
510   END SUBROUTINE dia_wri
511
512#else
513   !!----------------------------------------------------------------------
514   !!   Default option                                  use IOIPSL  library
515   !!----------------------------------------------------------------------
516
517   INTEGER FUNCTION dia_wri_alloc()
518      !!----------------------------------------------------------------------
519      INTEGER, DIMENSION(2) :: ierr
520      !!----------------------------------------------------------------------
521      IF( nn_write == -1 ) THEN
522         dia_wri_alloc = 0
523      ELSE   
524         ierr = 0
525         ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
526            &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
527            &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
528         !
529         dia_wri_alloc = MAXVAL(ierr)
530         CALL mpp_sum( 'diawri', dia_wri_alloc )
531         !
532      ENDIF
533      !
534   END FUNCTION dia_wri_alloc
535 
536   INTEGER FUNCTION dia_wri_alloc_abl()
537      !!----------------------------------------------------------------------
538     ALLOCATE(   ndex_hA(jpi*jpj), ndex_A (jpi*jpj*jpkam1), STAT=dia_wri_alloc_abl)
539      CALL mpp_sum( 'diawri', dia_wri_alloc_abl )
540      !
541   END FUNCTION dia_wri_alloc_abl
542
543   
544   SUBROUTINE dia_wri( kt, Kmm )
545      !!---------------------------------------------------------------------
546      !!                  ***  ROUTINE dia_wri  ***
547      !!                   
548      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
549      !!      NETCDF format is used by default
550      !!
551      !! ** Method  :   At the beginning of the first time step (nit000),
552      !!      define all the NETCDF files and fields
553      !!      At each time step call histdef to compute the mean if ncessary
554      !!      Each nn_write time step, output the instantaneous or mean fields
555      !!----------------------------------------------------------------------
556      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
557      INTEGER, INTENT( in ) ::   Kmm  ! ocean time level index
558      !
559      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
560      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
561      INTEGER  ::   inum = 11                                ! temporary logical unit
562      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
563      INTEGER  ::   ierr                                     ! error code return from allocation
564      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
565      INTEGER  ::   ipka                                     ! ABL
566      INTEGER  ::   jn, ierror                               ! local integers
567      REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
568      !
569      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace
570      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d, ze3t, zgdept       ! 3D workspace
571      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl   ! ABL 3D workspace
572      !!----------------------------------------------------------------------
573      !
574      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==!
575         CALL dia_wri_state( Kmm, 'output.init' )
576         ninist = 0
577      ENDIF
578      !
579      IF( nn_write == -1 )   RETURN   ! we will never do any output
580      !
581      IF( ln_timing )   CALL timing_start('dia_wri')
582      !
583      ! 0. Initialisation
584      ! -----------------
585
586      ll_print = .FALSE.                  ! local variable for debugging
587      ll_print = ll_print .AND. lwp
588
589      ! Define frequency of output and means
590      clop = "x"         ! no use of the mask value (require less cpu time and otherwise the model crashes)
591#if defined key_diainstant
592      zsto = nn_write * rn_Dt
593      clop = "inst("//TRIM(clop)//")"
594#else
595      zsto=rn_Dt
596      clop = "ave("//TRIM(clop)//")"
597#endif
598      zout = nn_write * rn_Dt
599      zmax = ( nitend - nit000 + 1 ) * rn_Dt
600
601      ! Define indices of the horizontal output zoom and vertical limit storage
602      iimi = Nis0   ;   iima = Nie0
603      ijmi = Njs0   ;   ijma = Nje0
604      ipk = jpk
605      IF(ln_abl) ipka = jpkam1
606
607      ! define time axis
608      it = kt
609      itmod = kt - nit000 + 1
610
611      ! store e3t for subsitute
612      DO jk = 1, jpk
613         ze3t  (:,:,jk) =  e3t  (:,:,jk,Kmm)
614         zgdept(:,:,jk) =  gdept(:,:,jk,Kmm)
615      END DO
616
617
618      ! 1. Define NETCDF files and fields at beginning of first time step
619      ! -----------------------------------------------------------------
620
621      IF( kt == nit000 ) THEN
622
623         ! Define the NETCDF files (one per grid)
624
625         ! Compute julian date from starting date of the run
626         CALL ymds2ju( nyear, nmonth, nday, rn_Dt, zjulian )
627         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
628         IF(lwp)WRITE(numout,*)
629         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
630            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
631         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
632                                 ' limit storage in depth = ', ipk
633
634         ! WRITE root name in date.file for use by postpro
635         IF(lwp) THEN
636            CALL dia_nam( clhstnam, nn_write,' ' )
637            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
638            WRITE(inum,*) clhstnam
639            CLOSE(inum)
640         ENDIF
641
642         ! Define the T grid FILE ( nid_T )
643
644         CALL dia_nam( clhstnam, nn_write, 'grid_T' )
645         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
646         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
647            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
648            &          nit000-1, zjulian, rn_Dt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
649         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
650            &           "m", ipk, gdept_1d, nz_T, "down" )
651         !                                                            ! Index of ocean points
652         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
653         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
654         !
655         IF( ln_icebergs ) THEN
656            !
657            !! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
658            !! that routine is called from nemogcm, so do it here immediately before its needed
659            ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
660            CALL mpp_sum( 'diawri', ierror )
661            IF( ierror /= 0 ) THEN
662               CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
663               RETURN
664            ENDIF
665            !
666            !! iceberg vertical coordinate is class number
667            CALL histvert( nid_T, "class", "Iceberg class",      &  ! Vertical grid: class
668               &           "number", nclasses, class_num, nb_T )
669            !
670            !! each class just needs the surface index pattern
671            ndim_bT = 3
672            DO jn = 1,nclasses
673               ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
674            ENDDO
675            !
676         ENDIF
677
678         ! Define the U grid FILE ( nid_U )
679
680         CALL dia_nam( clhstnam, nn_write, 'grid_U' )
681         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
682         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
683            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
684            &          nit000-1, zjulian, rn_Dt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
685         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
686            &           "m", ipk, gdept_1d, nz_U, "down" )
687         !                                                            ! Index of ocean points
688         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
689         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
690
691         ! Define the V grid FILE ( nid_V )
692
693         CALL dia_nam( clhstnam, nn_write, 'grid_V' )                   ! filename
694         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
695         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
696            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
697            &          nit000-1, zjulian, rn_Dt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
698         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
699            &          "m", ipk, gdept_1d, nz_V, "down" )
700         !                                                            ! Index of ocean points
701         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
702         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
703
704         ! Define the W grid FILE ( nid_W )
705
706         CALL dia_nam( clhstnam, nn_write, 'grid_W' )                   ! filename
707         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
708         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
709            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
710            &          nit000-1, zjulian, rn_Dt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
711         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
712            &          "m", ipk, gdepw_1d, nz_W, "down" )
713
714         IF( ln_abl ) THEN 
715         ! Define the ABL grid FILE ( nid_A )
716            CALL dia_nam( clhstnam, nn_write, 'grid_ABL' )
717            IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
718            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
719               &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
720               &          nit000-1, zjulian, rn_Dt, nh_A, nid_A, domain_id=nidom, snc4chunks=snc4set )
721            CALL histvert( nid_A, "ght_abl", "Vertical T levels",      &  ! Vertical grid: gdept
722               &           "m", ipka, ght_abl(2:jpka), nz_A, "up" )
723            !                                                            ! Index of ocean points
724         ALLOCATE( zw3d_abl(jpi,jpj,ipka) ) 
725         zw3d_abl(:,:,:) = 1._wp 
726         CALL wheneq( jpi*jpj*ipka, zw3d_abl, 1, 1., ndex_A , ndim_A  )      ! volume
727            CALL wheneq( jpi*jpj     , zw3d_abl, 1, 1., ndex_hA, ndim_hA )      ! surface
728         DEALLOCATE(zw3d_abl)
729         ENDIF
730         !
731
732         ! Declare all the output fields as NETCDF variables
733
734         !                                                                                      !!! nid_T : 3D
735         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
736            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
737         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
738            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
739         IF(  .NOT.ln_linssh  ) THEN
740            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t n
741            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
742            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t n
743            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
744            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t n
745            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
746         ENDIF
747         !                                                                                      !!! nid_T : 2D
748         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
749            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
750         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
751            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
752         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
753            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
754         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
755            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
756         CALL histdef( nid_T, "sorunoff", "River runoffs"                      , "Kg/m2/s",   &  ! runoffs
757            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
758         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
759            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
760         IF(  ln_linssh  ) THEN
761            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * ts(:,:,1,jp_tem,Kmm)
762            &                                                                  , "KgC/m2/s",  &  ! sosst_cd
763            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
764            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * ts(:,:,1,jp_sal,Kmm)
765            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
766            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
767         ENDIF
768         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
769            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
770         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
771            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
772         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
773            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
774         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
775            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
776         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
777            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
778         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
779            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
780         !
781         IF( ln_abl ) THEN
782            CALL histdef( nid_A, "t_abl", "Potential Temperature"     , "K"        ,       &  ! t_abl
783               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
784            CALL histdef( nid_A, "q_abl", "Humidity"                  , "kg/kg"    ,       &  ! q_abl
785               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
786            CALL histdef( nid_A, "u_abl", "Atmospheric U-wind   "     , "m/s"        ,     &  ! u_abl
787               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
788            CALL histdef( nid_A, "v_abl", "Atmospheric V-wind   "     , "m/s"    ,         &  ! v_abl
789               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
790            CALL histdef( nid_A, "tke_abl", "Atmospheric TKE   "     , "m2/s2"    ,        &  ! tke_abl
791               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
792            CALL histdef( nid_A, "avm_abl", "Atmospheric turbulent viscosity", "m2/s"   ,  &  ! avm_abl
793               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
794            CALL histdef( nid_A, "avt_abl", "Atmospheric turbulent diffusivity", "m2/s2",  &  ! avt_abl
795               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
796            CALL histdef( nid_A, "pblh", "Atmospheric boundary layer height "  , "m",      &  ! pblh
797               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout )                 
798#if defined key_si3
799            CALL histdef( nid_A, "oce_frac", "Fraction of open ocean"  , " ",      &  ! ato_i
800               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout )
801#endif
802            CALL histend( nid_A, snc4chunks=snc4set )
803         ENDIF
804         !
805         IF( ln_icebergs ) THEN
806            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
807               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
808            CALL histdef( nid_T, "calving_heat"        , "calving heat flux"                        , "XXXX"   , &
809               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
810            CALL histdef( nid_T, "berg_floating_melt"  , "Melt rate of icebergs + bits"             , "kg/m2/s", &
811               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
812            CALL histdef( nid_T, "berg_stored_ice"     , "Accumulated ice mass by class"            , "kg"     , &
813               &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
814            IF( ln_bergdia ) THEN
815               CALL histdef( nid_T, "berg_melt"           , "Melt rate of icebergs"                    , "kg/m2/s", &
816                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
817               CALL histdef( nid_T, "berg_buoy_melt"      , "Buoyancy component of iceberg melt rate"  , "kg/m2/s", &
818                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
819               CALL histdef( nid_T, "berg_eros_melt"      , "Erosion component of iceberg melt rate"   , "kg/m2/s", &
820                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
821               CALL histdef( nid_T, "berg_conv_melt"      , "Convective component of iceberg melt rate", "kg/m2/s", &
822                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
823               CALL histdef( nid_T, "berg_virtual_area"   , "Virtual coverage by icebergs"             , "m2"     , &
824                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
825               CALL histdef( nid_T, "bits_src"           , "Mass source of bergy bits"                , "kg/m2/s", &
826                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
827               CALL histdef( nid_T, "bits_melt"          , "Melt rate of bergy bits"                  , "kg/m2/s", &
828                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
829               CALL histdef( nid_T, "bits_mass"          , "Bergy bit density field"                  , "kg/m2"  , &
830                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
831               CALL histdef( nid_T, "berg_mass"           , "Iceberg density field"                    , "kg/m2"  , &
832                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
833               CALL histdef( nid_T, "berg_real_calving"   , "Calving into iceberg class"               , "kg/s"   , &
834                  &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
835            ENDIF
836         ENDIF
837
838         IF( ln_ssr ) THEN
839            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
840               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
841            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
842               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
843            CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
844               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
845         ENDIF
846       
847         clmx ="l_max(only(x))"    ! max index on a period
848!         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
849!            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
850#if defined key_diahth
851         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
852            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
853         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
854            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
855         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
856            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
857         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "J/m2"   ,   & ! htc3
858            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
859#endif
860
861         CALL histend( nid_T, snc4chunks=snc4set )
862
863         !                                                                                      !!! nid_U : 3D
864         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! uu(:,:,:,Kmm)
865            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
866         IF( ln_wave .AND. ln_sdw) THEN
867            CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current"         , "m/s"    ,   &  ! usd
868               &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
869         ENDIF
870         !                                                                                      !!! nid_U : 2D
871         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
872            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
873
874         CALL histend( nid_U, snc4chunks=snc4set )
875
876         !                                                                                      !!! nid_V : 3D
877         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vv(:,:,:,Kmm)
878            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
879         IF( ln_wave .AND. ln_sdw) THEN
880            CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current"    , "m/s"    ,   &  ! vsd
881               &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
882         ENDIF
883         !                                                                                      !!! nid_V : 2D
884         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
885            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
886
887         CALL histend( nid_V, snc4chunks=snc4set )
888
889         !                                                                                      !!! nid_W : 3D
890         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! ww
891            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
892         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
893            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
894         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avm
895            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
896
897         IF( ln_zdfddm ) THEN
898            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
899               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
900         ENDIF
901         
902         IF( ln_wave .AND. ln_sdw) THEN
903            CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current"   , "m/s"    ,   &  ! wsd
904               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
905         ENDIF
906         !                                                                                      !!! nid_W : 2D
907         CALL histend( nid_W, snc4chunks=snc4set )
908
909         IF(lwp) WRITE(numout,*)
910         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
911         IF(ll_print) CALL FLUSH(numout )
912
913      ENDIF
914
915      ! 2. Start writing data
916      ! ---------------------
917
918      ! ndex(1) est utilise ssi l'avant dernier argument est different de
919      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
920      ! donne le nombre d'elements, et ndex la liste des indices a sortir
921
922      IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN
923         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
924         WRITE(numout,*) '~~~~~~ '
925      ENDIF
926
927      IF( .NOT.ln_linssh ) THEN
928         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T  )   ! heat content
929         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T  )   ! salt content
930         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
931         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
932      ELSE
933         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T  )   ! temperature
934         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) , ndim_T , ndex_T  )   ! salinity
935         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) , ndim_hT, ndex_hT )   ! sea surface temperature
936         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity
937      ENDIF
938      IF( .NOT.ln_linssh ) THEN
939         zw3d(:,:,:) = ( ( ze3t(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
940         CALL histwrite( nid_T, "vovvle3t", it, ze3t (:,:,:)     , ndim_T , ndex_T  )   ! level thickness
941         CALL histwrite( nid_T, "vovvldep", it, zgdept , ndim_T , ndex_T  )   ! t-point depth
942         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
943      ENDIF
944      CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm)          , ndim_hT, ndex_hT )   ! sea surface height
945      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
946      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
947      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
948                                                                                  ! (includes virtual salt flux beneath ice
949                                                                                  ! in linear free surface case)
950      IF( ln_linssh ) THEN
951         zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_tem,Kmm)
952         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
953         zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_sal,Kmm)
954         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
955      ENDIF
956      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
957      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
958      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
959      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
960      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
961      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
962      !
963      IF( ln_abl ) THEN
964         ALLOCATE( zw3d_abl(jpi,jpj,jpka) )
965         IF( ln_mskland )   THEN
966            DO jk=1,jpka
967               zw3d_abl(:,:,jk) = tmask(:,:,1)
968            END DO       
969         ELSE
970            zw3d_abl(:,:,:) = 1._wp     
971         ENDIF       
972         CALL histwrite( nid_A,  "pblh"   , it, pblh(:,:)                  *zw3d_abl(:,:,1     ), ndim_hA, ndex_hA )   ! pblh
973         CALL histwrite( nid_A,  "u_abl"  , it, u_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! u_abl
974         CALL histwrite( nid_A,  "v_abl"  , it, v_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! v_abl
975         CALL histwrite( nid_A,  "t_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! t_abl
976         CALL histwrite( nid_A,  "q_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! q_abl       
977         CALL histwrite( nid_A,  "tke_abl", it, tke_abl (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! tke_abl
978         CALL histwrite( nid_A,  "avm_abl", it, avm_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avm_abl
979         CALL histwrite( nid_A,  "avt_abl", it, avt_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avt_abl
980#if defined key_si3
981         CALL histwrite( nid_A,  "oce_frac"   , it, ato_i(:,:)                                  , ndim_hA, ndex_hA )   ! ato_i
982#endif
983         DEALLOCATE(zw3d_abl)
984      ENDIF
985      !
986      IF( ln_icebergs ) THEN
987         !
988         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT ) 
989         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
990         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT ) 
991         !
992         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
993         !
994         IF( ln_bergdia ) THEN
995            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   ) 
996            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   ) 
997            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   ) 
998            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   ) 
999            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   ) 
1000            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   ) 
1001            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   ) 
1002            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   ) 
1003            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   ) 
1004            !
1005            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
1006         ENDIF
1007      ENDIF
1008
1009      IF( ln_ssr ) THEN
1010         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
1011         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
1012         zw2d(:,:) = erp(:,:) * ts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
1013         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
1014      ENDIF
1015!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
1016!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
1017
1018#if defined key_diahth
1019      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
1020      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
1021      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
1022      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
1023#endif
1024
1025      CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm)            , ndim_U , ndex_U )    ! i-current
1026      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
1027
1028      CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm)            , ndim_V , ndex_V  )   ! j-current
1029      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
1030
1031      IF( ln_zad_Aimp ) THEN
1032         CALL histwrite( nid_W, "vovecrtz", it, ww + wi     , ndim_T, ndex_T )    ! vert. current
1033      ELSE
1034         CALL histwrite( nid_W, "vovecrtz", it, ww          , ndim_T, ndex_T )    ! vert. current
1035      ENDIF
1036      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
1037      CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
1038      IF( ln_zdfddm ) THEN
1039         CALL histwrite( nid_W, "voddmavs", it, avs         , ndim_T, ndex_T )    ! S vert. eddy diff. coef.
1040      ENDIF
1041
1042      IF( ln_wave .AND. ln_sdw ) THEN
1043         CALL histwrite( nid_U, "sdzocrtx", it, usd         , ndim_U , ndex_U )    ! i-StokesDrift-current
1044         CALL histwrite( nid_V, "sdmecrty", it, vsd         , ndim_V , ndex_V )    ! j-StokesDrift-current
1045         CALL histwrite( nid_W, "sdvecrtz", it, wsd         , ndim_T , ndex_T )    ! StokesDrift vert. current
1046      ENDIF
1047
1048      ! 3. Close all files
1049      ! ---------------------------------------
1050      IF( kt == nitend ) THEN
1051         CALL histclo( nid_T )
1052         CALL histclo( nid_U )
1053         CALL histclo( nid_V )
1054         CALL histclo( nid_W )
1055         IF(ln_abl) CALL histclo( nid_A )
1056      ENDIF
1057      !
1058      IF( ln_timing )   CALL timing_stop('dia_wri')
1059      !
1060   END SUBROUTINE dia_wri
1061#endif
1062
1063   SUBROUTINE dia_wri_state( Kmm, cdfile_name )
1064      !!---------------------------------------------------------------------
1065      !!                 ***  ROUTINE dia_wri_state  ***
1066      !!       
1067      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
1068      !!      the instantaneous ocean state and forcing fields.
1069      !!        Used to find errors in the initial state or save the last
1070      !!      ocean state in case of abnormal end of a simulation
1071      !!
1072      !! ** Method  :   NetCDF files using ioipsl
1073      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
1074      !!      File 'output.abort.nc' is created in case of abnormal job end
1075      !!----------------------------------------------------------------------
1076      INTEGER           , INTENT( in ) ::   Kmm              ! time level index
1077      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
1078      !!
1079      INTEGER :: inum, jk
1080      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, zgdept       ! 3D workspace for qco substitution
1081      !!----------------------------------------------------------------------
1082      !
1083      IF(lwp) THEN
1084         WRITE(numout,*)
1085         WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
1086         WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
1087         WRITE(numout,*) '                and named :', cdfile_name, '...nc'
1088      ENDIF 
1089      !
1090      DO jk = 1, jpk
1091         ze3t(:,:,jk) =  e3t(:,:,jk,Kmm)
1092         zgdept(:,:,jk) =  gdept(:,:,jk,Kmm)
1093      END DO
1094      !
1095      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
1096      !
1097      CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) )    ! now temperature
1098      CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) )    ! now salinity
1099      CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm)              )    ! sea surface height
1100      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)                )    ! now i-velocity
1101      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)                )    ! now j-velocity
1102      IF( ln_zad_Aimp ) THEN
1103         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww + wi        )    ! now k-velocity
1104      ELSE
1105         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww             )    ! now k-velocity
1106      ENDIF
1107      CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep            )    ! now k-velocity
1108      CALL iom_rstput( 0, 0, inum, 'ht'     , ht(:,:)            )    ! now water column height
1109      !
1110      IF ( ln_isf ) THEN
1111         IF (ln_isfcav_mlt) THEN
1112            CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav          )    ! now k-velocity
1113            CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav    )    ! now k-velocity
1114            CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav    )    ! now k-velocity
1115            CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,wp) )    ! now k-velocity
1116            CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,wp) )    ! now k-velocity
1117            CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,wp), ktype = jp_i1 )
1118         END IF
1119         IF (ln_isfpar_mlt) THEN
1120            CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,wp) )    ! now k-velocity
1121            CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par          )    ! now k-velocity
1122            CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par    )    ! now k-velocity
1123            CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par    )    ! now k-velocity
1124            CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,wp) )    ! now k-velocity
1125            CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,wp) )    ! now k-velocity
1126            CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,wp), ktype = jp_i1 )
1127         END IF
1128      END IF
1129      !
1130      IF( ALLOCATED(ahtu) ) THEN
1131         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point
1132         CALL iom_rstput( 0, 0, inum,  'ahtv', ahtv              )    ! aht at v-point
1133      ENDIF
1134      IF( ALLOCATED(ahmt) ) THEN
1135         CALL iom_rstput( 0, 0, inum,  'ahmt', ahmt              )    ! ahmt at u-point
1136         CALL iom_rstput( 0, 0, inum,  'ahmf', ahmf              )    ! ahmf at v-point
1137      ENDIF
1138      CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf         )    ! freshwater budget
1139      CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns         )    ! total heat flux
1140      CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr               )    ! solar heat flux
1141      CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i              )    ! ice fraction
1142      CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
1143      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
1144      IF(  .NOT.ln_linssh  ) THEN             
1145         CALL iom_rstput( 0, 0, inum, 'vovvldep', zgdept        )    !  T-cell depth
1146         CALL iom_rstput( 0, 0, inum, 'vovvle3t', ze3t          )    !  T-cell thickness 
1147      END IF
1148      IF( ln_wave .AND. ln_sdw ) THEN
1149         CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd            )    ! now StokesDrift i-velocity
1150         CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd            )    ! now StokesDrift j-velocity
1151         CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity
1152      ENDIF
1153      IF ( ln_abl ) THEN
1154         CALL iom_rstput ( 0, 0, inum, "uz1_abl",   u_abl(:,:,2,nt_a  ) )   ! now first level i-wind
1155         CALL iom_rstput ( 0, 0, inum, "vz1_abl",   v_abl(:,:,2,nt_a  ) )   ! now first level j-wind
1156         CALL iom_rstput ( 0, 0, inum, "tz1_abl",  tq_abl(:,:,2,nt_a,1) )   ! now first level temperature
1157         CALL iom_rstput ( 0, 0, inum, "qz1_abl",  tq_abl(:,:,2,nt_a,2) )   ! now first level humidity
1158      ENDIF
1159      IF( ln_zdfosm ) THEN
1160         CALL iom_rstput( 0, 0, inum, 'hbl', hbl*tmask(:,:,1)  )      ! now boundary-layer depth
1161         CALL iom_rstput( 0, 0, inum, 'hml', hml*tmask(:,:,1)  )      ! now mixed-layer depth
1162         CALL iom_rstput( 0, 0, inum, 'avt_k', avt_k*wmask     )      ! w-level diffusion
1163         CALL iom_rstput( 0, 0, inum, 'avm_k', avm_k*wmask     )      ! now w-level viscosity
1164         CALL iom_rstput( 0, 0, inum, 'ghamt', ghamt*wmask     )      ! non-local t forcing
1165         CALL iom_rstput( 0, 0, inum, 'ghams', ghams*wmask     )      ! non-local s forcing
1166         CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu*umask     )      ! non-local u forcing
1167         CALL iom_rstput( 0, 0, inum, 'ghamv', ghamv*vmask     )      ! non-local v forcing
1168         IF( ln_osm_mle ) THEN
1169            CALL iom_rstput( 0, 0, inum, 'hmle', hmle*tmask(:,:,1)  ) ! now transition-layer depth
1170         END IF
1171      ENDIF
1172      !
1173      CALL iom_close( inum )
1174      !
1175#if defined key_si3
1176      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid
1177         CALL iom_open( TRIM(cdfile_name)//'_ice', inum, ldwrt = .TRUE., kdlev = jpl, cdcomp = 'ICE' )
1178         CALL ice_wri_state( inum )
1179         CALL iom_close( inum )
1180      ENDIF
1181      !
1182#endif
1183   END SUBROUTINE dia_wri_state
1184
1185   !!======================================================================
1186END MODULE diawri
Note: See TracBrowser for help on using the repository browser.