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/r12377_ticket2386/src/OCE/DIA – NEMO

source: NEMO/branches/2020/r12377_ticket2386/src/OCE/DIA/diawri.F90 @ 13540

Last change on this file since 13540 was 13540, checked in by andmirek, 4 years ago

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