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/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DIA – NEMO

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DIA/diawri.F90 @ 14644

Last change on this file since 14644 was 14644, checked in by sparonuz, 3 years ago

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