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_r13327_KERNEL-06_2_techene_e3/src/OCE/DIA – NEMO

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/DIA/diawri.F90 @ 13696

Last change on this file since 13696 was 13605, checked in by techene, 4 years ago

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