New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diawri.F90 in NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DIA – NEMO

source: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DIA/diawri.F90 @ 11521

Last change on this file since 11521 was 11521, checked in by mathiot, 5 years ago

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