source: NEMO/branches/UKMO/NEMO_4.0.2_mirror/src/OCE/DIA/diawri.F90 @ 12658

Last change on this file since 12658 was 12658, checked in by cguiavarch, 10 months ago

UKMO/NEMO_4.0.2_mirror : Remove SVN keywords.

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