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

source: branches/UKMO/dev_r5518_GO6_package_fix_rnf_MOCI_TEST_SUITE/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 8191

Last change on this file since 8191 was 8191, checked in by andmirek, 7 years ago

merge with XIOS restart read branch

File size: 58.8 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
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<<<<<<< .working
160      CALL iom_put("tpt_dep", fsdept_n(:,:,:) )
161      CALL iom_put("wpt_dep", fsdepw_n(:,:,:) )
162=======
163      CALL iom_put("tpt_dep", fsdept_n(:,:,:) )
164>>>>>>> .merge-right.r7923
165
166
167      CALL iom_put( "ssh" , sshn )                 ! sea surface height
168     
169      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
170      CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature
171      IF ( iom_use("sbt") ) THEN
172         DO jj = 1, jpj
173            DO ji = 1, jpi
174               jkbot = mbkt(ji,jj)
175               z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem)
176            END DO
177         END DO
178         CALL iom_put( "sbt", z2d )                ! bottom temperature
179      ENDIF
180     
181      CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity
182      CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity
183      IF ( iom_use("sbs") ) THEN
184         DO jj = 1, jpj
185            DO ji = 1, jpi
186               jkbot = mbkt(ji,jj)
187               z2d(ji,jj) = tsn(ji,jj,jkbot,jp_sal)
188            END DO
189         END DO
190         CALL iom_put( "sbs", z2d )                ! bottom salinity
191      ENDIF
192
193      IF ( iom_use("taubot") ) THEN                ! bottom stress
194         z2d(:,:) = 0._wp
195         DO jj = 2, jpjm1
196            DO ji = fs_2, fs_jpim1   ! vector opt.
197               zztmpx = (  bfrua(ji  ,jj) * un(ji  ,jj,mbku(ji  ,jj))  &
198                      &  + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj))  )     
199               zztmpy = (  bfrva(ji,  jj) * vn(ji,jj  ,mbkv(ji,jj  ))  &
200                      &  + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1))  ) 
201               z2d(ji,jj) = rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1) 
202               !
203            ENDDO
204         ENDDO
205         CALL lbc_lnk( z2d, 'T', 1. )
206         CALL iom_put( "taubot", z2d )           
207      ENDIF
208         
209      CALL iom_put( "uoce", un(:,:,:)         )    ! 3D i-current
210      CALL iom_put(  "ssu", un(:,:,1)         )    ! surface i-current
211      IF ( iom_use("sbu") ) THEN
212         DO jj = 1, jpj
213            DO ji = 1, jpi
214               jkbot = mbku(ji,jj)
215               z2d(ji,jj) = un(ji,jj,jkbot)
216            END DO
217         END DO
218         CALL iom_put( "sbu", z2d )                ! bottom i-current
219      ENDIF
220#if defined key_dynspg_ts
221      CALL iom_put(  "ubar", un_adv(:,:)      )    ! barotropic i-current
222#else
223      CALL iom_put(  "ubar", un_b(:,:)        )    ! barotropic i-current
224#endif
225     
226      CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current
227      CALL iom_put(  "ssv", vn(:,:,1)         )    ! surface j-current
228      IF ( iom_use("sbv") ) THEN
229         DO jj = 1, jpj
230            DO ji = 1, jpi
231               jkbot = mbkv(ji,jj)
232               z2d(ji,jj) = vn(ji,jj,jkbot)
233            END DO
234         END DO
235         CALL iom_put( "sbv", z2d )                ! bottom j-current
236      ENDIF
237#if defined key_dynspg_ts
238      CALL iom_put(  "vbar", vn_adv(:,:)      )    ! barotropic j-current
239#else
240      CALL iom_put(  "vbar", vn_b(:,:)        )    ! barotropic j-current
241#endif
242
243      CALL iom_put( "woce", wn )                   ! vertical velocity
244      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
245         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
246         z2d(:,:) = rau0 * e12t(:,:)
247         DO jk = 1, jpk
248            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
249         END DO
250         CALL 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   ||   defined key_esopa
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 iom_put( "u_masstr", z3d )                  ! mass transport in i-direction
341         CALL iom_put( "u_masstr_vint", z2d )             ! mass transport in i-direction vertical sum
342      ENDIF
343     
344      IF( iom_use("u_heattr") ) THEN
345         z2d(:,:) = 0.e0 
346         DO jk = 1, jpkm1
347            DO jj = 2, jpjm1
348               DO ji = fs_2, fs_jpim1   ! vector opt.
349                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
350               END DO
351            END DO
352         END DO
353         CALL lbc_lnk( z2d, 'U', -1. )
354         CALL iom_put( "u_heattr", (0.5 * rcp) * z2d )    ! heat transport in i-direction
355      ENDIF
356
357      IF( iom_use("u_salttr") ) THEN
358         z2d(:,:) = 0.e0 
359         DO jk = 1, jpkm1
360            DO jj = 2, jpjm1
361               DO ji = fs_2, fs_jpim1   ! vector opt.
362                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
363               END DO
364            END DO
365         END DO
366         CALL lbc_lnk( z2d, 'U', -1. )
367         CALL iom_put( "u_salttr", 0.5 * z2d )            ! heat transport in i-direction
368      ENDIF
369
370     
371      IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
372         z3d(:,:,jpk) = 0.e0
373         DO jk = 1, jpkm1
374            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) * vmask(:,:,jk)
375         END DO
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!     IF( nn_timing == 1 )   CALL timing_stop('dia_wri_state') ! not sure this works for routines not called in first timestep
1126      !
1127
1128   END SUBROUTINE dia_wri_state
1129   !!======================================================================
1130END MODULE diawri
Note: See TracBrowser for help on using the repository browser.