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

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14_ssh_dt/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 13074

Last change on this file since 13074 was 13074, checked in by mattmartin, 5 years ago

Include SSH tendency diagnostic.

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