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 branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/NERC/dev_r5518_GO6_under_ice_relax/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

Last change on this file was 10047, checked in by jpalmier, 6 years ago

merge with GO6_package_branch 9385-10020 ; plus debug OMIP_DIC

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