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

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 6348

Last change on this file since 6348 was 6348, checked in by cetlod, 8 years ago

bugfix: move the output of scale factor before time swapping and output some variables needs for offline, see ticket #1682

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