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

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

trunk: avoid implicit loops in diawri
  • Property svn:keywords set to Id
File size: 65.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_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    ) :: z2d     ! 2D workspace
586      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d     ! 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      ! 1. Define NETCDF files and fields at beginning of first time step
628      ! -----------------------------------------------------------------
629
630      IF( kt == nit000 ) THEN
631
632         ! Define the NETCDF files (one per grid)
633
634         ! Compute julian date from starting date of the run
635         CALL ymds2ju( nyear, nmonth, nday, rn_Dt, zjulian )
636         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
637         IF(lwp)WRITE(numout,*)
638         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
639            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
640         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
641                                 ' limit storage in depth = ', ipk
642
643         ! WRITE root name in date.file for use by postpro
644         IF(lwp) THEN
645            CALL dia_nam( clhstnam, nn_write,' ' )
646            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
647            WRITE(inum,*) clhstnam
648            CLOSE(inum)
649         ENDIF
650
651         ! Define the T grid FILE ( nid_T )
652
653         CALL dia_nam( clhstnam, nn_write, 'grid_T' )
654         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
655         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
656            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
657            &          nit000-1, zjulian, rn_Dt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
658         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
659            &           "m", ipk, gdept_1d, nz_T, "down" )
660         !                                                            ! Index of ocean points
661         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
662         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
663         !
664         IF( ln_icebergs ) THEN
665            !
666            !! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
667            !! that routine is called from nemogcm, so do it here immediately before its needed
668            ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
669            CALL mpp_sum( 'diawri', ierror )
670            IF( ierror /= 0 ) THEN
671               CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
672               RETURN
673            ENDIF
674            !
675            !! iceberg vertical coordinate is class number
676            CALL histvert( nid_T, "class", "Iceberg class",      &  ! Vertical grid: class
677               &           "number", nclasses, class_num, nb_T )
678            !
679            !! each class just needs the surface index pattern
680            ndim_bT = 3
681            DO jn = 1,nclasses
682               ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
683            ENDDO
684            !
685         ENDIF
686
687         ! Define the U grid FILE ( nid_U )
688
689         CALL dia_nam( clhstnam, nn_write, 'grid_U' )
690         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
691         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
692            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
693            &          nit000-1, zjulian, rn_Dt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
694         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
695            &           "m", ipk, gdept_1d, nz_U, "down" )
696         !                                                            ! Index of ocean points
697         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
698         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
699
700         ! Define the V grid FILE ( nid_V )
701
702         CALL dia_nam( clhstnam, nn_write, 'grid_V' )                   ! filename
703         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
704         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
705            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
706            &          nit000-1, zjulian, rn_Dt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
707         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
708            &          "m", ipk, gdept_1d, nz_V, "down" )
709         !                                                            ! Index of ocean points
710         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
711         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
712
713         ! Define the W grid FILE ( nid_W )
714
715         CALL dia_nam( clhstnam, nn_write, 'grid_W' )                   ! filename
716         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
717         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
718            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
719            &          nit000-1, zjulian, rn_Dt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
720         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
721            &          "m", ipk, gdepw_1d, nz_W, "down" )
722
723         IF( ln_abl ) THEN 
724         ! Define the ABL grid FILE ( nid_A )
725            CALL dia_nam( clhstnam, nn_write, 'grid_ABL' )
726            IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
727            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
728               &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
729               &          nit000-1, zjulian, rn_Dt, nh_A, nid_A, domain_id=nidom, snc4chunks=snc4set )
730            CALL histvert( nid_A, "ght_abl", "Vertical T levels",      &  ! Vertical grid: gdept
731               &           "m", ipka, ght_abl(2:jpka), nz_A, "up" )
732            !                                                            ! Index of ocean points
733         ALLOCATE( zw3d_abl(jpi,jpj,ipka) ) 
734         zw3d_abl(:,:,:) = 1._wp 
735         CALL wheneq( jpi*jpj*ipka, zw3d_abl, 1, 1., ndex_A , ndim_A  )      ! volume
736            CALL wheneq( jpi*jpj     , zw3d_abl, 1, 1., ndex_hA, ndim_hA )      ! surface
737         DEALLOCATE(zw3d_abl)
738         ENDIF
739         !
740
741         ! Declare all the output fields as NETCDF variables
742
743         !                                                                                      !!! nid_T : 3D
744         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
745            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
746         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
747            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
748         IF(  .NOT.ln_linssh  ) THEN
749            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t n
750            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
751            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t n
752            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
753            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t n
754            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
755         ENDIF
756         !                                                                                      !!! nid_T : 2D
757         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
758            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
759         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
760            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
761         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
762            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
763         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
764            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
765         CALL histdef( nid_T, "sorunoff", "River runoffs"                      , "Kg/m2/s",   &  ! runoffs
766            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
767         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
768            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
769         IF(  ln_linssh  ) THEN
770            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * ts(:,:,1,jp_tem,Kmm)
771            &                                                                  , "KgC/m2/s",  &  ! sosst_cd
772            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
773            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * ts(:,:,1,jp_sal,Kmm)
774            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
775            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
776         ENDIF
777         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
778            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
779         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
780            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
781         IF( ALLOCATED(hmld) ) THEN   ! zdf_mxl not called by SWE
782            CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
783               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
784            CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
785               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
786         ENDIF
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         DO_3D( 0, 0, 0, 0, 1, jpk )
940            z3d(ji,jj,jk) = ts(ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm)
941         END_3D
942         CALL histwrite( nid_T, "votemper", it, z3d, ndim_T , ndex_T  )   ! heat content
943         DO_3D( 0, 0, 0, 0, 1, jpk )
944            z3d(ji,jj,jk) = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm)
945         END_3D
946         CALL histwrite( nid_T, "vosaline", it, z3d, ndim_T , ndex_T  )   ! salt content
947         DO_2D( 0, 0, 0, 0 )
948            z2d(ji,jj   ) = ts(ji,jj, 1,jp_tem,Kmm) * e3t(ji,jj, 1,Kmm)
949         END_2D
950         CALL histwrite( nid_T, "sosstsst", it, z2d, ndim_hT, ndex_hT )   ! sea surface heat content
951         DO_2D( 0, 0, 0, 0 )
952            z2d(ji,jj   ) = ts(ji,jj, 1,jp_sal,Kmm) * e3t(ji,jj, 1,Kmm)
953         END_2D
954         CALL histwrite( nid_T, "sosaline", it, z2d, ndim_hT, ndex_hT )   ! sea surface salinity content
955      ELSE
956         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T  )   ! temperature
957         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) , ndim_T , ndex_T  )   ! salinity
958         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) , ndim_hT, ndex_hT )   ! sea surface temperature
959         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity
960      ENDIF
961      IF( .NOT.ln_linssh ) THEN
962         DO_3D( 0, 0, 0, 0, 1, jpk )
963           z3d(ji,jj,jk) = e3t(ji,jj,jk,Kmm)     ! 3D workspace for qco substitution
964         END_3D
965         CALL histwrite( nid_T, "vovvle3t", it, z3d        , ndim_T , ndex_T  )   ! level thickness
966         DO_3D( 0, 0, 0, 0, 1, jpk )
967           z3d(ji,jj,jk) = gdept(ji,jj,jk,Kmm)   ! 3D workspace for qco substitution
968         END_3D
969         CALL histwrite( nid_T, "vovvldep", it, z3d        , ndim_T , ndex_T  )   ! t-point depth
970         DO_3D( 0, 0, 0, 0, 1, jpk )
971            z3d(ji,jj,jk) = ( ( e3t(ji,jj,jk,Kmm) - e3t_0(ji,jj,jk) ) / e3t_0(ji,jj,jk) * 100._wp * tmask(ji,jj,jk) ) ** 2
972         END_3D         
973         CALL histwrite( nid_T, "vovvldef", it, z3d        , ndim_T , ndex_T  )   ! level thickness deformation
974      ENDIF
975      CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm)  , ndim_hT, ndex_hT )   ! sea surface height
976      DO_2D( 0, 0, 0, 0 )
977         z2d(ji,jj) = emp(ji,jj) - rnf(ji,jj)
978      END_2D
979      CALL histwrite( nid_T, "sowaflup", it, z2d           , ndim_hT, ndex_hT )   ! upward water flux
980      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
981      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
982                                                                                  ! (includes virtual salt flux beneath ice
983                                                                                  ! in linear free surface case)
984      IF( ln_linssh ) THEN
985         DO_2D( 0, 0, 0, 0 )
986            z2d(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_tem,Kmm)
987         END_2D
988         CALL histwrite( nid_T, "sosst_cd", it, z2d, ndim_hT, ndex_hT )          ! c/d term on sst
989         DO_2D( 0, 0, 0, 0 )
990            z2d(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_sal,Kmm)
991         END_2D
992         CALL histwrite( nid_T, "sosss_cd", it, z2d, ndim_hT, ndex_hT )          ! c/d term on sss
993      ENDIF
994      DO_2D( 0, 0, 0, 0 )
995         z2d(ji,jj) = qsr(ji,jj) + qns(ji,jj)
996      END_2D
997      CALL histwrite( nid_T, "sohefldo", it, z2d           , ndim_hT, ndex_hT )   ! total heat flux
998      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
999      IF( ALLOCATED(hmld) ) THEN   ! zdf_mxl not called by SWE
1000         CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
1001         CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
1002      ENDIF
1003      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
1004      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
1005      !
1006      IF( ln_abl ) THEN
1007         ALLOCATE( zw3d_abl(jpi,jpj,jpka) )
1008         IF( ln_mskland )   THEN
1009            DO jk=1,jpka
1010               zw3d_abl(:,:,jk) = tmask(:,:,1)
1011            END DO       
1012         ELSE
1013            zw3d_abl(:,:,:) = 1._wp     
1014         ENDIF       
1015         CALL histwrite( nid_A,  "pblh"   , it, pblh(:,:)                  *zw3d_abl(:,:,1     ), ndim_hA, ndex_hA )   ! pblh
1016         CALL histwrite( nid_A,  "u_abl"  , it, u_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! u_abl
1017         CALL histwrite( nid_A,  "v_abl"  , it, v_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! v_abl
1018         CALL histwrite( nid_A,  "t_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! t_abl
1019         CALL histwrite( nid_A,  "q_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! q_abl       
1020         CALL histwrite( nid_A,  "tke_abl", it, tke_abl (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! tke_abl
1021         CALL histwrite( nid_A,  "avm_abl", it, avm_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avm_abl
1022         CALL histwrite( nid_A,  "avt_abl", it, avt_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avt_abl
1023#if defined key_si3
1024         CALL histwrite( nid_A,  "oce_frac"   , it, ato_i(:,:)                                  , ndim_hA, ndex_hA )   ! ato_i
1025#endif
1026         DEALLOCATE(zw3d_abl)
1027      ENDIF
1028      !
1029      IF( ln_icebergs ) THEN
1030         !
1031         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT ) 
1032         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
1033         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT ) 
1034         !
1035         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
1036         !
1037         IF( ln_bergdia ) THEN
1038            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   ) 
1039            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   ) 
1040            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   ) 
1041            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   ) 
1042            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   ) 
1043            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   ) 
1044            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   ) 
1045            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   ) 
1046            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   ) 
1047            !
1048            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
1049         ENDIF
1050      ENDIF
1051
1052      IF( ln_ssr ) THEN
1053         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
1054         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
1055         DO_2D( 0, 0, 0, 0 )
1056            z2d(ji,jj) = erp(ji,jj) * ts(ji,jj,1,jp_sal,Kmm) * tmask(ji,jj,1)
1057         END_2D
1058         CALL histwrite( nid_T, "sosafldp", it, z2d           , ndim_hT, ndex_hT )   ! salt flux damping
1059      ENDIF
1060!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
1061!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
1062
1063#if defined key_diahth
1064      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
1065      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
1066      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
1067      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
1068#endif
1069
1070      CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm) , ndim_U , ndex_U )    ! i-current
1071      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
1072
1073      CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm) , ndim_V , ndex_V  )   ! j-current
1074      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
1075
1076      IF( ln_zad_Aimp ) THEN
1077         DO_3D( 0, 0, 0, 0, 1, jpk )
1078           z3d(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
1079         END_3D         
1080         CALL histwrite( nid_W, "vovecrtz", it, z3d         , ndim_T, ndex_T )    ! vert. current
1081      ELSE
1082         CALL histwrite( nid_W, "vovecrtz", it, ww          , ndim_T, ndex_T )    ! vert. current
1083      ENDIF
1084      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
1085      CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
1086      IF( ln_zdfddm ) THEN
1087         CALL histwrite( nid_W, "voddmavs", it, avs         , ndim_T, ndex_T )    ! S vert. eddy diff. coef.
1088      ENDIF
1089
1090      IF( ln_wave .AND. ln_sdw ) THEN
1091         CALL histwrite( nid_U, "sdzocrtx", it, usd         , ndim_U , ndex_U )    ! i-StokesDrift-current
1092         CALL histwrite( nid_V, "sdmecrty", it, vsd         , ndim_V , ndex_V )    ! j-StokesDrift-current
1093         CALL histwrite( nid_W, "sdvecrtz", it, wsd         , ndim_T , ndex_T )    ! StokesDrift vert. current
1094      ENDIF
1095
1096      ! 3. Close all files
1097      ! ---------------------------------------
1098      IF( kt == nitend ) THEN
1099         CALL histclo( nid_T )
1100         CALL histclo( nid_U )
1101         CALL histclo( nid_V )
1102         CALL histclo( nid_W )
1103         IF(ln_abl) CALL histclo( nid_A )
1104      ENDIF
1105      !
1106      IF( ln_timing )   CALL timing_stop('dia_wri')
1107      !
1108   END SUBROUTINE dia_wri
1109#endif
1110
1111   SUBROUTINE dia_wri_state( Kmm, cdfile_name )
1112      !!---------------------------------------------------------------------
1113      !!                 ***  ROUTINE dia_wri_state  ***
1114      !!       
1115      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
1116      !!      the instantaneous ocean state and forcing fields.
1117      !!        Used to find errors in the initial state or save the last
1118      !!      ocean state in case of abnormal end of a simulation
1119      !!
1120      !! ** Method  :   NetCDF files using ioipsl
1121      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
1122      !!      File 'output.abort.nc' is created in case of abnormal job end
1123      !!----------------------------------------------------------------------
1124      INTEGER           , INTENT( in ) ::   Kmm              ! time level index
1125      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
1126      !!
1127      INTEGER ::   ji, jj, jk       ! dummy loop indices
1128      INTEGER ::   inum
1129      REAL(wp), DIMENSION(jpi,jpj)     :: z2d     
1130      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d     
1131      !!----------------------------------------------------------------------
1132      !
1133      IF(lwp) THEN
1134         WRITE(numout,*)
1135         WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
1136         WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
1137         WRITE(numout,*) '                and named :', cdfile_name, '...nc'
1138      ENDIF 
1139      !
1140      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
1141      !
1142      CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) )    ! now temperature
1143      CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) )    ! now salinity
1144      CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm)         )    ! sea surface height
1145      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)        )    ! now i-velocity
1146      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)        )    ! now j-velocity
1147      IF( ln_zad_Aimp ) THEN
1148         DO_3D( 0, 0, 0, 0, 1, jpk )
1149           z3d(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
1150         END_3D
1151         CALL iom_rstput( 0, 0, inum, 'vovecrtz', z3d            )    ! now k-velocity
1152      ELSE
1153         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww             )    ! now k-velocity
1154      ENDIF
1155      CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep            )
1156      CALL iom_rstput( 0, 0, inum, 'ht'     , ht(:,:)            )    ! now water column height
1157      !
1158      IF ( ln_isf ) THEN
1159         IF (ln_isfcav_mlt) THEN
1160            CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav          )
1161            CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav    )
1162            CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav    )
1163            CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,wp) )
1164            CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,wp) )
1165            CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,wp), ktype = jp_i1 )
1166         END IF
1167         IF (ln_isfpar_mlt) THEN
1168            CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,wp) )
1169            CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par          )
1170            CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par    )
1171            CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par    )
1172            CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,wp) )
1173            CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,wp) )
1174            CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,wp), ktype = jp_i1 )
1175         END IF
1176      END IF
1177      !
1178      IF( ALLOCATED(ahtu) ) THEN
1179         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point
1180         CALL iom_rstput( 0, 0, inum,  'ahtv', ahtv              )    ! aht at v-point
1181      ENDIF
1182      IF( ALLOCATED(ahmt) ) THEN
1183         CALL iom_rstput( 0, 0, inum,  'ahmt', ahmt              )    ! ahmt at u-point
1184         CALL iom_rstput( 0, 0, inum,  'ahmf', ahmf              )    ! ahmf at v-point
1185      ENDIF
1186      DO_2D( 0, 0, 0, 0 )
1187         z2d(ji,jj) = emp(ji,jj) - rnf(ji,jj)
1188      END_2D
1189      CALL iom_rstput( 0, 0, inum, 'sowaflup', z2d               )    ! freshwater budget
1190      DO_2D( 0, 0, 0, 0 )
1191         z2d(ji,jj) = qsr(ji,jj) + qns(ji,jj)
1192      END_2D
1193      CALL iom_rstput( 0, 0, inum, 'sohefldo', z2d               )    ! total heat flux
1194      CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr               )    ! solar heat flux
1195      CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i              )    ! ice fraction
1196      CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
1197      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
1198      IF(  .NOT.ln_linssh  ) THEN
1199         DO_3D( 0, 0, 0, 0, 1, jpk )
1200           z3d(ji,jj,jk) = gdept(ji,jj,jk,Kmm)   ! 3D workspace for qco substitution
1201         END_3D
1202         CALL iom_rstput( 0, 0, inum, 'vovvldep', z3d            )    !  T-cell depth
1203         DO_3D( 0, 0, 0, 0, 1, jpk )
1204           z3d(ji,jj,jk) = e3t(ji,jj,jk,Kmm)     ! 3D workspace for qco substitution
1205         END_3D
1206         CALL iom_rstput( 0, 0, inum, 'vovvle3t', z3d            )    !  T-cell thickness 
1207      END IF
1208      IF( ln_wave .AND. ln_sdw ) THEN
1209         CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd            )    ! now StokesDrift i-velocity
1210         CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd            )    ! now StokesDrift j-velocity
1211         CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity
1212      ENDIF
1213      IF ( ln_abl ) THEN
1214         CALL iom_rstput ( 0, 0, inum, "uz1_abl",   u_abl(:,:,2,nt_a  ) )   ! now first level i-wind
1215         CALL iom_rstput ( 0, 0, inum, "vz1_abl",   v_abl(:,:,2,nt_a  ) )   ! now first level j-wind
1216         CALL iom_rstput ( 0, 0, inum, "tz1_abl",  tq_abl(:,:,2,nt_a,1) )   ! now first level temperature
1217         CALL iom_rstput ( 0, 0, inum, "qz1_abl",  tq_abl(:,:,2,nt_a,2) )   ! now first level humidity
1218      ENDIF
1219      IF( ln_zdfosm ) THEN
1220         CALL iom_rstput( 0, 0, inum, 'hbl', hbl*tmask(:,:,1)  )      ! now boundary-layer depth
1221         CALL iom_rstput( 0, 0, inum, 'hml', hml*tmask(:,:,1)  )      ! now mixed-layer depth
1222         CALL iom_rstput( 0, 0, inum, 'avt_k', avt_k*wmask     )      ! w-level diffusion
1223         CALL iom_rstput( 0, 0, inum, 'avm_k', avm_k*wmask     )      ! now w-level viscosity
1224         CALL iom_rstput( 0, 0, inum, 'ghamt', ghamt*wmask     )      ! non-local t forcing
1225         CALL iom_rstput( 0, 0, inum, 'ghams', ghams*wmask     )      ! non-local s forcing
1226         CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu*umask     )      ! non-local u forcing
1227         CALL iom_rstput( 0, 0, inum, 'ghamv', ghamv*vmask     )      ! non-local v forcing
1228         IF( ln_osm_mle ) THEN
1229            CALL iom_rstput( 0, 0, inum, 'hmle', hmle*tmask(:,:,1)  ) ! now transition-layer depth
1230         END IF
1231      ENDIF
1232      !
1233      CALL iom_close( inum )
1234      !
1235#if defined key_si3
1236      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid
1237         CALL iom_open( TRIM(cdfile_name)//'_ice', inum, ldwrt = .TRUE., kdlev = jpl, cdcomp = 'ICE' )
1238         CALL ice_wri_state( inum )
1239         CALL iom_close( inum )
1240      ENDIF
1241      !
1242#endif
1243   END SUBROUTINE dia_wri_state
1244
1245   !!======================================================================
1246END MODULE diawri
Note: See TracBrowser for help on using the repository browser.