source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 4378

Last change on this file since 4378 was 4378, checked in by rfurner, 7 years ago

mask values before outputing in vvl case, see ticket 1213

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