New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diawri.F90 in branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

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

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

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