source: branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 11883

Last change on this file since 11883 was 11883, checked in by timgraham, 11 months ago

Mergerd in changes from Chris' working copy (except for symlink changes)

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