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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 4570

Last change on this file since 4570 was 4570, checked in by vichi, 10 years ago

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