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

source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 4747

Last change on this file since 4747 was 4747, checked in by mathiot, 10 years ago

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