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

source: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 7158

Last change on this file since 7158 was 7158, checked in by clem, 7 years ago

debug branch
  • 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      clop = "x"         ! no use of the mask value (require less cpu time, and otherwise the model crashes)
450#if defined key_diainstant
451      zsto = nwrite * zdt
452      clop = "inst("//TRIM(clop)//")"
453#else
454      zsto=zdt
455      clop = "ave("//TRIM(clop)//")"
456#endif
457      zout = nwrite * zdt
458      zmax = ( nitend - nit000 + 1 ) * zdt
459
460      ! Define indices of the horizontal output zoom and vertical limit storage
461      iimi = 1      ;      iima = jpi
462      ijmi = 1      ;      ijma = jpj
463      ipk = jpk
464
465      ! define time axis
466      it = kt
467      itmod = kt - nit000 + 1
468
469
470      ! 1. Define NETCDF files and fields at beginning of first time step
471      ! -----------------------------------------------------------------
472
473      IF( kt == nit000 ) THEN
474
475         ! Define the NETCDF files (one per grid)
476
477         ! Compute julian date from starting date of the run
478         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
479         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
480         IF(lwp)WRITE(numout,*)
481         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
482            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
483         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
484                                 ' limit storage in depth = ', ipk
485
486         ! WRITE root name in date.file for use by postpro
487         IF(lwp) THEN
488            CALL dia_nam( clhstnam, nwrite,' ' )
489            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
490            WRITE(inum,*) clhstnam
491            CLOSE(inum)
492         ENDIF
493
494         ! Define the T grid FILE ( nid_T )
495
496         CALL dia_nam( clhstnam, nwrite, 'grid_T' )
497         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
498         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
499            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
500            &          nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
501         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
502            &           "m", ipk, gdept_1d, nz_T, "down" )
503         !                                                            ! Index of ocean points
504         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
505         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
506         !
507         IF( ln_icebergs ) THEN
508            !
509            !! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
510            !! that routine is called from nemogcm, so do it here immediately before its needed
511            ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
512            IF( lk_mpp )   CALL mpp_sum( ierror )
513            IF( ierror /= 0 ) THEN
514               CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
515               RETURN
516            ENDIF
517            !
518            !! iceberg vertical coordinate is class number
519            CALL histvert( nid_T, "class", "Iceberg class",      &  ! Vertical grid: class
520               &           "number", nclasses, class_num, nb_T )
521            !
522            !! each class just needs the surface index pattern
523            ndim_bT = 3
524            DO jn = 1,nclasses
525               ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
526            ENDDO
527            !
528         ENDIF
529
530         ! Define the U grid FILE ( nid_U )
531
532         CALL dia_nam( clhstnam, nwrite, 'grid_U' )
533         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
534         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
535            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
536            &          nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
537         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
538            &           "m", ipk, gdept_1d, nz_U, "down" )
539         !                                                            ! Index of ocean points
540         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
541         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
542
543         ! Define the V grid FILE ( nid_V )
544
545         CALL dia_nam( clhstnam, nwrite, 'grid_V' )                   ! filename
546         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
547         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
548            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
549            &          nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
550         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
551            &          "m", ipk, gdept_1d, nz_V, "down" )
552         !                                                            ! Index of ocean points
553         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
554         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
555
556         ! Define the W grid FILE ( nid_W )
557
558         CALL dia_nam( clhstnam, nwrite, 'grid_W' )                   ! filename
559         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
560         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
561            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
562            &          nit000-1, zjulian, zdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
563         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
564            &          "m", ipk, gdepw_1d, nz_W, "down" )
565
566
567         ! Declare all the output fields as NETCDF variables
568
569         !                                                                                      !!! nid_T : 3D
570         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
571            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
572         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
573            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
574         IF(  lk_vvl  ) THEN
575            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n
576            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
577            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t_n
578            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
579            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t_n
580            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
581         ENDIF
582         !                                                                                      !!! nid_T : 2D
583         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
584            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
585         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
586            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
587         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
588            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
589         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
590            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
591         CALL histdef( nid_T, "sorunoff", "River runoffs"                      , "Kg/m2/s",   &  ! runoffs
592            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
593         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
594            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
595         IF(  .NOT. lk_vvl  ) THEN
596            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem)
597            &                                                                  , "KgC/m2/s",  &  ! sosst_cd
598            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
599            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * tsn(:,:,1,jp_sal)
600            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
601            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
602         ENDIF
603         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
604            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
605         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
606            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
607         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
608            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
609         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
610            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
611         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
612            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
613         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
614            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
615!
616         IF( ln_icebergs ) THEN
617            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
618               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
619            CALL histdef( nid_T, "calving_heat"        , "calving heat flux"                        , "XXXX"   , &
620               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
621            CALL histdef( nid_T, "berg_floating_melt"  , "Melt rate of icebergs + bits"             , "kg/m2/s", &
622               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
623            CALL histdef( nid_T, "berg_stored_ice"     , "Accumulated ice mass by class"            , "kg"     , &
624               &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
625            IF( ln_bergdia ) THEN
626               CALL histdef( nid_T, "berg_melt"           , "Melt rate of icebergs"                    , "kg/m2/s", &
627                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
628               CALL histdef( nid_T, "berg_buoy_melt"      , "Buoyancy component of iceberg melt rate"  , "kg/m2/s", &
629                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
630               CALL histdef( nid_T, "berg_eros_melt"      , "Erosion component of iceberg melt rate"   , "kg/m2/s", &
631                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
632               CALL histdef( nid_T, "berg_conv_melt"      , "Convective component of iceberg melt rate", "kg/m2/s", &
633                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
634               CALL histdef( nid_T, "berg_virtual_area"   , "Virtual coverage by icebergs"             , "m2"     , &
635                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
636               CALL histdef( nid_T, "bits_src"           , "Mass source of bergy bits"                , "kg/m2/s", &
637                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
638               CALL histdef( nid_T, "bits_melt"          , "Melt rate of bergy bits"                  , "kg/m2/s", &
639                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
640               CALL histdef( nid_T, "bits_mass"          , "Bergy bit density field"                  , "kg/m2"  , &
641                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
642               CALL histdef( nid_T, "berg_mass"           , "Iceberg density field"                    , "kg/m2"  , &
643                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
644               CALL histdef( nid_T, "berg_real_calving"   , "Calving into iceberg class"               , "kg/s"   , &
645                  &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
646            ENDIF
647         ENDIF
648
649         IF( .NOT. ln_cpl ) THEN
650            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
651               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
652            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
653               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
654            CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
655               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
656         ENDIF
657
658         IF( ln_cpl .AND. nn_ice <= 1 ) THEN
659            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
660               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
661            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
662               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
663            CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn
664               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
665         ENDIF
666         
667         clmx ="l_max(only(x))"    ! max index on a period
668         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
669            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
670#if defined key_diahth
671         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
672            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
673         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
674            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
675         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
676            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
677         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "W"      ,   & ! htc3
678            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
679#endif
680
681         IF( ln_cpl .AND. nn_ice == 2 ) THEN
682            CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature"            , "K"      ,   &  ! tn_ice
683               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
684            CALL histdef( nid_T,"soicealb" , "Ice Albedo"                         , "[0,1]"  ,   &  ! alb_ice
685               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
686         ENDIF
687
688         CALL histend( nid_T, snc4chunks=snc4set )
689
690         !                                                                                      !!! nid_U : 3D
691         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
692            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
693         IF( ln_traldf_gdia ) THEN
694            CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
695                 &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
696         ELSE
697#if defined key_diaeiv
698            CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
699            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
700#endif
701         END IF
702         !                                                                                      !!! nid_U : 2D
703         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
704            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
705
706         CALL histend( nid_U, snc4chunks=snc4set )
707
708         !                                                                                      !!! nid_V : 3D
709         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
710            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
711         IF( ln_traldf_gdia ) THEN
712            CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
713                 &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
714         ELSE 
715#if defined key_diaeiv
716            CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
717            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
718#endif
719         END IF
720         !                                                                                      !!! nid_V : 2D
721         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
722            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
723
724         CALL histend( nid_V, snc4chunks=snc4set )
725
726         !                                                                                      !!! nid_W : 3D
727         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn
728            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
729         IF( ln_traldf_gdia ) THEN
730            CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
731                 &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
732         ELSE
733#if defined key_diaeiv
734            CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
735                 &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
736#endif
737         END IF
738         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
739            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
740         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avmu
741            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
742
743         IF( lk_zdfddm ) THEN
744            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
745               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
746         ENDIF
747         !                                                                                      !!! nid_W : 2D
748#if defined key_traldf_c2d
749         CALL histdef( nid_W, "soleahtw", "lateral eddy diffusivity"           , "m2/s"   ,   &  ! ahtw
750            &          jpi, jpj, nh_W, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
751# if defined key_traldf_eiv 
752            CALL histdef( nid_W, "soleaeiw", "eddy induced vel. coeff. at w-point", "m2/s",   &  ! aeiw
753               &       jpi, jpj, nh_W, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
754# endif
755#endif
756
757         CALL histend( nid_W, snc4chunks=snc4set )
758
759         IF(lwp) WRITE(numout,*)
760         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
761         IF(ll_print) CALL FLUSH(numout )
762
763      ENDIF
764
765      ! 2. Start writing data
766      ! ---------------------
767
768      ! ndex(1) est utilise ssi l'avant dernier argument est different de
769      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
770      ! donne le nombre d'elements, et ndex la liste des indices a sortir
771
772      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
773         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
774         WRITE(numout,*) '~~~~~~ '
775      ENDIF
776
777      IF( lk_vvl ) THEN
778         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
779         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
780         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
781         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
782      ELSE
783         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature
784         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity
785         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature
786         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
787      ENDIF
788      IF( lk_vvl ) THEN
789         zw3d(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
790         CALL histwrite( nid_T, "vovvle3t", it, fse3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
791         CALL histwrite( nid_T, "vovvldep", it, fsdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
792         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
793      ENDIF
794      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
795      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
796      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
797      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
798                                                                                  ! (includes virtual salt flux beneath ice
799                                                                                  ! in linear free surface case)
800      IF( .NOT. lk_vvl ) THEN
801         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
802         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
803         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
804         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
805      ENDIF
806      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
807      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
808      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
809      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
810      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
811      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
812!
813      IF( ln_icebergs ) THEN
814         !
815         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT ) 
816         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
817         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT ) 
818         !
819         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
820         !
821         IF( ln_bergdia ) THEN
822            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   ) 
823            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   ) 
824            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   ) 
825            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   ) 
826            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   ) 
827            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   ) 
828            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   ) 
829            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   ) 
830            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   ) 
831            !
832            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
833         ENDIF
834      ENDIF
835
836      IF( .NOT. ln_cpl ) THEN
837         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
838         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
839         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
840         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
841      ENDIF
842      IF( ln_cpl .AND. nn_ice <= 1 ) THEN
843         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
844         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
845         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
846         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
847      ENDIF
848!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
849!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
850
851#if defined key_diahth
852      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
853      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
854      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
855      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
856#endif
857
858      IF( ln_cpl .AND. nn_ice == 2 ) THEN
859         CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature
860         CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo
861      ENDIF
862
863      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
864      IF( ln_traldf_gdia ) THEN
865         IF (.not. ALLOCATED(psix_eiv))THEN
866            ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr )
867            IF( lk_mpp   )   CALL mpp_sum ( ierr )
868            IF( ierr > 0 )   CALL ctl_stop('STOP', 'diawri: unable to allocate psi{x,y}_eiv')
869            psix_eiv(:,:,:) = 0.0_wp
870            psiy_eiv(:,:,:) = 0.0_wp
871         ENDIF
872         DO jk=1,jpkm1
873            zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk)  ! u_eiv = -dpsix/dz
874         END DO
875         zw3d(:,:,jpk) = 0._wp
876         CALL histwrite( nid_U, "vozoeivu", it, zw3d, ndim_U , ndex_U )           ! i-eiv current
877      ELSE
878#if defined key_diaeiv
879         CALL histwrite( nid_U, "vozoeivu", it, u_eiv, ndim_U , ndex_U )          ! i-eiv current
880#endif
881      ENDIF
882      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
883
884      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
885      IF( ln_traldf_gdia ) THEN
886         DO jk=1,jpk-1
887            zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk)  ! v_eiv = -dpsiy/dz
888         END DO
889         zw3d(:,:,jpk) = 0._wp
890         CALL histwrite( nid_V, "vomeeivv", it, zw3d, ndim_V , ndex_V )           ! j-eiv current
891      ELSE
892#if defined key_diaeiv
893         CALL histwrite( nid_V, "vomeeivv", it, v_eiv, ndim_V , ndex_V )          ! j-eiv current
894#endif
895      ENDIF
896      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
897
898      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
899      IF( ln_traldf_gdia ) THEN
900         DO jk=1,jpk-1
901            DO jj = 2, jpjm1
902               DO ji = fs_2, fs_jpim1  ! vector opt.
903                  zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2v(ji,jj) + &
904                       &    (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1u(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx
905               END DO
906            END DO
907         END DO
908         zw3d(:,:,jpk) = 0._wp
909         CALL histwrite( nid_W, "voveeivw", it, zw3d          , ndim_T, ndex_T )    ! vert. eiv current
910      ELSE
911#   if defined key_diaeiv
912         CALL histwrite( nid_W, "voveeivw", it, w_eiv          , ndim_T, ndex_T )    ! vert. eiv current
913#   endif
914      ENDIF
915      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
916      CALL histwrite( nid_W, "votkeavm", it, avmu           , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
917      IF( lk_zdfddm ) THEN
918         CALL histwrite( nid_W, "voddmavs", it, fsavs(:,:,:), ndim_T, ndex_T )    ! S vert. eddy diff. coef.
919      ENDIF
920#if defined key_traldf_c2d
921      CALL histwrite( nid_W, "soleahtw", it, ahtw          , ndim_hT, ndex_hT )   ! lateral eddy diff. coef.
922# if defined key_traldf_eiv
923         CALL histwrite( nid_W, "soleaeiw", it, aeiw       , ndim_hT, ndex_hT )   ! EIV coefficient at w-point
924# endif
925#endif
926
927      ! 3. Close all files
928      ! ---------------------------------------
929      IF( kt == nitend ) THEN
930         CALL histclo( nid_T )
931         CALL histclo( nid_U )
932         CALL histclo( nid_V )
933         CALL histclo( nid_W )
934      ENDIF
935      !
936      CALL wrk_dealloc( jpi , jpj      , zw2d )
937      IF ( ln_traldf_gdia .OR. lk_vvl )  call wrk_dealloc( jpi , jpj , jpk  , zw3d )
938      !
939      IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
940      !
941   END SUBROUTINE dia_wri
942# endif
943
944#endif
945
946   SUBROUTINE dia_wri_state( cdfile_name, kt )
947      !!---------------------------------------------------------------------
948      !!                 ***  ROUTINE dia_wri_state  ***
949      !!       
950      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
951      !!      the instantaneous ocean state and forcing fields.
952      !!        Used to find errors in the initial state or save the last
953      !!      ocean state in case of abnormal end of a simulation
954      !!
955      !! ** Method  :   NetCDF files using ioipsl
956      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
957      !!      File 'output.abort.nc' is created in case of abnormal job end
958      !!----------------------------------------------------------------------
959      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
960      INTEGER           , INTENT( in ) ::   kt               ! ocean time-step index
961      !!
962      CHARACTER (len=32) :: clname
963      CHARACTER (len=40) :: clop
964      INTEGER  ::   id_i , nz_i, nh_i       
965      INTEGER, DIMENSION(1) ::   idex             ! local workspace
966      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt
967      !!----------------------------------------------------------------------
968      !
969!     IF( nn_timing == 1 )   CALL timing_start('dia_wri_state') ! not sure this works for routines not called in first timestep
970
971      ! 0. Initialisation
972      ! -----------------
973
974      ! Define name, frequency of output and means
975      clname = cdfile_name
976      IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
977      zdt  = rdt
978      zsto = rdt
979      clop = "inst(x)"           ! no use of the mask value (require less cpu time)
980      zout = rdt
981      zmax = ( nitend - nit000 + 1 ) * zdt
982
983      IF(lwp) WRITE(numout,*)
984      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
985      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
986      IF(lwp) WRITE(numout,*) '                and named :', clname, '.nc'
987
988
989      ! 1. Define NETCDF files and fields at beginning of first time step
990      ! -----------------------------------------------------------------
991
992      ! Compute julian date from starting date of the run
993      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )         ! time axis
994      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
995      CALL histbeg( clname, jpi, glamt, jpj, gphit,   &
996          1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
997      CALL histvert( id_i, "deptht", "Vertical T levels",   &    ! Vertical grid : gdept
998          "m", jpk, gdept_1d, nz_i, "down")
999
1000      ! Declare all the output fields as NetCDF variables
1001
1002      CALL histdef( id_i, "vosaline", "Salinity"              , "PSU"    ,   &   ! salinity
1003         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1004      CALL histdef( id_i, "votemper", "Temperature"           , "C"      ,   &   ! temperature
1005         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1006      CALL histdef( id_i, "sossheig", "Sea Surface Height"    , "m"      ,   &  ! ssh
1007         &          jpi, jpj, nh_i, 1  , 1, 1  , nz_i, 32, clop, zsto, zout )
1008      CALL histdef( id_i, "vozocrtx", "Zonal Current"         , "m/s"    ,   &   ! zonal current
1009         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1010      CALL histdef( id_i, "vomecrty", "Meridional Current"    , "m/s"    ,   &   ! meridonal current
1011         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
1012      CALL histdef( id_i, "vovecrtz", "Vertical Velocity"     , "m/s"    ,   &   ! vertical current
1013         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
1014      CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S",   &   ! net freshwater
1015         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1016      CALL histdef( id_i, "sohefldo", "Net Downward Heat Flux", "W/m2"   ,   &   ! net heat flux
1017         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1018      CALL histdef( id_i, "soshfldo", "Shortwave Radiation"   , "W/m2"   ,   &   ! solar flux
1019         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1020      CALL histdef( id_i, "soicecov", "Ice fraction"          , "[0,1]"  ,   &   ! fr_i
1021         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1022      CALL histdef( id_i, "sozotaux", "Zonal Wind Stress"     , "N/m2"   ,   &   ! i-wind stress
1023         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1024      CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2"   ,   &   ! j-wind stress
1025         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1026      IF( lk_vvl ) THEN
1027         CALL histdef( id_i, "vovvldep", "T point depth"         , "m"      ,   &   ! t-point depth
1028            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1029         CALL histdef( id_i, "vovvle3t", "T point thickness"         , "m"      ,   &   ! t-point depth
1030            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1031      END IF
1032
1033#if defined key_lim2
1034      CALL lim_wri_state_2( kt, id_i, nh_i )
1035#elif defined key_lim3
1036      CALL lim_wri_state( kt, id_i, nh_i )
1037#else
1038      CALL histend( id_i, snc4chunks=snc4set )
1039#endif
1040
1041      ! 2. Start writing data
1042      ! ---------------------
1043      ! idex(1) est utilise ssi l'avant dernier argument est diffferent de
1044      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
1045      ! donne le nombre d'elements, et idex la liste des indices a sortir
1046      idex(1) = 1   ! init to avoid compil warning
1047
1048      ! Write all fields on T grid
1049      CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature
1050      CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity
1051      CALL histwrite( id_i, "sossheig", kt, sshn             , jpi*jpj    , idex )    ! sea surface height
1052      CALL histwrite( id_i, "vozocrtx", kt, un               , jpi*jpj*jpk, idex )    ! now i-velocity
1053      CALL histwrite( id_i, "vomecrty", kt, vn               , jpi*jpj*jpk, idex )    ! now j-velocity
1054      CALL histwrite( id_i, "vovecrtz", kt, wn               , jpi*jpj*jpk, idex )    ! now k-velocity
1055      CALL histwrite( id_i, "sowaflup", kt, (emp-rnf )       , jpi*jpj    , idex )    ! freshwater budget
1056      CALL histwrite( id_i, "sohefldo", kt, qsr + qns        , jpi*jpj    , idex )    ! total heat flux
1057      CALL histwrite( id_i, "soshfldo", kt, qsr              , jpi*jpj    , idex )    ! solar heat flux
1058      CALL histwrite( id_i, "soicecov", kt, fr_i             , jpi*jpj    , idex )    ! ice fraction
1059      CALL histwrite( id_i, "sozotaux", kt, utau             , jpi*jpj    , idex )    ! i-wind stress
1060      CALL histwrite( id_i, "sometauy", kt, vtau             , jpi*jpj    , idex )    ! j-wind stress
1061      IF( lk_vvl ) THEN
1062         CALL histwrite( id_i, "vovvldep", kt, fsdept_n(:,:,:), jpi*jpj*jpk, idex )!  T-cell depth       
1063         CALL histwrite( id_i, "vovvle3t", kt, fse3t_n (:,:,:), jpi*jpj*jpk, idex )!  T-cell thickness 
1064      END IF
1065
1066      ! 3. Close the file
1067      ! -----------------
1068      CALL histclo( id_i )
1069#if ! defined key_iomput && ! defined key_dimgout
1070      IF( ninist /= 1  ) THEN
1071         CALL histclo( nid_T )
1072         CALL histclo( nid_U )
1073         CALL histclo( nid_V )
1074         CALL histclo( nid_W )
1075      ENDIF
1076#endif
1077       
1078!     IF( nn_timing == 1 )   CALL timing_stop('dia_wri_state') ! not sure this works for routines not called in first timestep
1079      !
1080
1081   END SUBROUTINE dia_wri_state
1082   !!======================================================================
1083END MODULE diawri
Note: See TracBrowser for help on using the repository browser.