source: NEMO/trunk/src/OCE/DIA/diawri.F90 @ 11418

Last change on this file since 11418 was 11418, checked in by acc, 2 years ago

Minor correction to make sure full vertical velocity is used in diagnostic output from diawri.F90 when using adaptive-implicit vertical advection (ln_zad_Aimp=.true.)

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