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

source: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 3953

Last change on this file since 3953 was 3884, checked in by acc, 11 years ago

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