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 @ 14179

Last change on this file since 14179 was 14179, checked in by techene, 4 years ago

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