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

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

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