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

source: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 11639

Last change on this file since 11639 was 11639, checked in by rrenshaw, 4 years ago

regional mean and transport diagnostics added for reanalysis

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