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 trunk/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 7 years ago

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