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 @ 5461

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

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