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 @ 13111

Last change on this file since 13111 was 13111, checked in by mattmartin, 4 years ago

Removed blank line after Jennie's review.

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