New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diawri.F90 in NEMO/trunk/src/OCE/DIA – NEMO

source: NEMO/trunk/src/OCE/DIA/diawri.F90 @ 14834

Last change on this file since 14834 was 14476, checked in by techene, 3 years ago

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