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

source: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 4792

Last change on this file since 4792 was 4792, checked in by jamesharle, 10 years ago

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