source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diawri.F90 @ 9939

Last change on this file since 9939 was 9939, checked in by gm, 2 years ago

#1911 (ENHANCE-04): RK3 branche phased with MLF@9937 branche

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