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_r12558_HPC-08_epico_Extra_Halo/tests/CANAL/MY_SRC – NEMO

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/tests/CANAL/MY_SRC/diawri.F90 @ 12939

Last change on this file since 12939 was 12939, checked in by smasson, 4 years ago

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