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 branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 9019

Last change on this file since 9019 was 9019, checked in by timgraham, 6 years ago

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