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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 5460

Last change on this file since 5460 was 5460, checked in by jchanut, 9 years ago

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