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

source: branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 4902

Last change on this file since 4902 was 4902, checked in by cetlod, 9 years ago

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