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/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 5866

Last change on this file since 5866 was 5866, checked in by gm, 8 years ago

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