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

source: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 6351

Last change on this file since 6351 was 6351, checked in by cetlod, 8 years ago

trunk: move the output of scale factor before time swapping and output some variables needs for offline, see ticket #1682
  • Property svn:keywords set to Id
File size: 53.4 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 sbc_oce         ! Surface boundary condition: ocean fields
34   USE sbc_ice         ! Surface boundary condition: ice fields
35   USE icb_oce         ! Icebergs
36   USE icbdia          ! Iceberg budgets
37   USE sbcssr          ! restoring term toward SST/SSS climatology
38   USE phycst          ! physical constants
39   USE zdfmxl          ! mixed layer
40   USE dianam          ! build name of file (routine)
41   USE zdfddm          ! vertical  physics: double diffusion
42   USE diahth          ! thermocline diagnostics
43   !
44   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
45   USE in_out_manager  ! I/O manager
46   USE diatmb          ! Top,middle,bottom output
47   USE dia25h          ! 25h Mean output
48   USE iom
49   USE ioipsl
50
51#if defined key_lim2
52   USE limwri_2 
53#elif defined key_lim3
54   USE limwri 
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   USE wrk_nemo        ! working array
61
62   IMPLICIT NONE
63   PRIVATE
64
65   PUBLIC   dia_wri                 ! routines called by step.F90
66   PUBLIC   dia_wri_state
67   PUBLIC   dia_wri_alloc           ! Called by nemogcm module
68
69   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
70   INTEGER ::          nb_T              , ndim_bT   ! grid_T file
71   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
72   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
73   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file
74   INTEGER ::   ndex(1)                              ! ???
75   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
76   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
77   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
78
79   !! * Substitutions
80#  include "zdfddm_substitute.h90"
81#  include "vectopt_loop_substitute.h90"
82   !!----------------------------------------------------------------------
83   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
84   !! $Id$
85   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
86   !!----------------------------------------------------------------------
87CONTAINS
88
89   INTEGER FUNCTION dia_wri_alloc()
90      !!----------------------------------------------------------------------
91      INTEGER, DIMENSION(2) :: ierr
92      !!----------------------------------------------------------------------
93      ierr = 0
94      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
95         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
96         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
97         !
98      dia_wri_alloc = MAXVAL(ierr)
99      IF( lk_mpp )   CALL mpp_sum( dia_wri_alloc )
100      !
101  END FUNCTION dia_wri_alloc
102
103   !!----------------------------------------------------------------------
104   !!   Default option                                   NetCDF output file
105   !!----------------------------------------------------------------------
106#if defined key_iomput
107   !!----------------------------------------------------------------------
108   !!   'key_iomput'                                        use IOM library
109   !!----------------------------------------------------------------------
110
111   SUBROUTINE dia_wri( kt )
112      !!---------------------------------------------------------------------
113      !!                  ***  ROUTINE dia_wri  ***
114      !!                   
115      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
116      !!      NETCDF format is used by default
117      !!
118      !! ** Method  :  use iom_put
119      !!----------------------------------------------------------------------
120      !!
121      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
122      !!
123      INTEGER                      ::   ji, jj, jk              ! dummy loop indices
124      INTEGER                      ::   jkbot                   !
125      REAL(wp)                     ::   zztmp, zztmpx, zztmpy   !
126      !!
127      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace
128      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace
129      !!----------------------------------------------------------------------
130      !
131      IF( nn_timing == 1 )   CALL timing_start('dia_wri')
132      !
133      CALL wrk_alloc( jpi , jpj      , z2d )
134      CALL wrk_alloc( jpi , jpj, jpk , z3d )
135      !
136      ! Output the initial state and forcings
137      IF( ninist == 1 ) THEN                       
138         CALL dia_wri_state( 'output.init', kt )
139         ninist = 0
140      ENDIF
141
142      ! Output of initial vertical scale factor
143      CALL iom_put("e3t_0", e3t_0(:,:,:) )
144      CALL iom_put("e3u_0", e3t_0(:,:,:) )
145      CALL iom_put("e3v_0", e3t_0(:,:,:) )
146      !
147      CALL iom_put( "e3t" , fse3t_n(:,:,:) )
148      CALL iom_put( "e3u" , fse3u_n(:,:,:) )
149      CALL iom_put( "e3v" , fse3v_n(:,:,:) )
150      CALL iom_put( "e3w" , fse3w_n(:,:,:) )
151      IF( iom_use("e3tdef") )   &
152         CALL iom_put( "e3tdef"  , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
153
154      CALL iom_put( "ssh" , sshn )                 ! sea surface height
155     
156      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
157      CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature
158      IF ( iom_use("sbt") ) THEN
159         DO jj = 1, jpj
160            DO ji = 1, jpi
161               jkbot = mbkt(ji,jj)
162               z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem)
163            END DO
164         END DO
165         CALL iom_put( "sbt", z2d )                ! bottom temperature
166      ENDIF
167     
168      CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity
169      CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity
170      IF ( iom_use("sbs") ) THEN
171         DO jj = 1, jpj
172            DO ji = 1, jpi
173               jkbot = mbkt(ji,jj)
174               z2d(ji,jj) = tsn(ji,jj,jkbot,jp_sal)
175            END DO
176         END DO
177         CALL iom_put( "sbs", z2d )                ! bottom salinity
178      ENDIF
179
180      IF ( iom_use("taubot") ) THEN                ! bottom stress
181         z2d(:,:) = 0._wp
182         DO jj = 2, jpjm1
183            DO ji = fs_2, fs_jpim1   ! vector opt.
184               zztmpx = (  bfrua(ji  ,jj) * un(ji  ,jj,mbku(ji  ,jj))  &
185                      &  + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj))  )     
186               zztmpy = (  bfrva(ji,  jj) * vn(ji,jj  ,mbkv(ji,jj  ))  &
187                      &  + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1))  ) 
188               z2d(ji,jj) = rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1) 
189               !
190            ENDDO
191         ENDDO
192         CALL lbc_lnk( z2d, 'T', 1. )
193         CALL iom_put( "taubot", z2d )           
194      ENDIF
195         
196      CALL iom_put( "uoce", un(:,:,:)         )    ! 3D i-current
197      CALL iom_put(  "ssu", un(:,:,1)         )    ! surface i-current
198      IF ( iom_use("sbu") ) THEN
199         DO jj = 1, jpj
200            DO ji = 1, jpi
201               jkbot = mbku(ji,jj)
202               z2d(ji,jj) = un(ji,jj,jkbot)
203            END DO
204         END DO
205         CALL iom_put( "sbu", z2d )                ! bottom i-current
206      ENDIF
207     
208      CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current
209      CALL iom_put(  "ssv", vn(:,:,1)         )    ! surface j-current
210      IF ( iom_use("sbv") ) THEN
211         DO jj = 1, jpj
212            DO ji = 1, jpi
213               jkbot = mbkv(ji,jj)
214               z2d(ji,jj) = vn(ji,jj,jkbot)
215            END DO
216         END DO
217         CALL iom_put( "sbv", z2d )                ! bottom j-current
218      ENDIF
219
220      CALL iom_put( "woce", wn )                   ! vertical velocity
221      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
222         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
223         z2d(:,:) = rau0 * e1e2t(:,:)
224         DO jk = 1, jpk
225            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
226         END DO
227         CALL iom_put( "w_masstr" , z3d ) 
228         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
229      ENDIF
230
231      CALL iom_put( "avt" , avt                        )    ! T vert. eddy diff. coef.
232      CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef.
233      CALL iom_put( "avs" , fsavs(:,:,:)               )    ! S vert. eddy diff. coef. (useful only with key_zdfddm)
234
235      IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt  (:,:,:) ) ) )
236      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, fsavs(:,:,:) ) ) )
237
238      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
239         DO jj = 2, jpjm1                                    ! sst gradient
240            DO ji = fs_2, fs_jpim1   ! vector opt.
241               zztmp  = tsn(ji,jj,1,jp_tem)
242               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)
243               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)
244               z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
245                  &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
246            END DO
247         END DO
248         CALL lbc_lnk( z2d, 'T', 1. )
249         CALL iom_put( "sstgrad2",  z2d               )    ! square of module of sst gradient
250         z2d(:,:) = SQRT( z2d(:,:) )
251         CALL iom_put( "sstgrad" ,  z2d               )    ! module of sst gradient
252      ENDIF
253         
254      ! clem: heat and salt content
255      IF( iom_use("heatc") ) THEN
256         z2d(:,:)  = 0._wp 
257         DO jk = 1, jpkm1
258            DO jj = 1, jpj
259               DO ji = 1, jpi
260                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
261               END DO
262            END DO
263         END DO
264         CALL iom_put( "heatc", (rau0 * rcp) * z2d )    ! vertically integrated heat content (J/m2)
265      ENDIF
266
267      IF( iom_use("saltc") ) THEN
268         z2d(:,:)  = 0._wp 
269         DO jk = 1, jpkm1
270            DO jj = 1, jpj
271               DO ji = 1, jpi
272                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
273               END DO
274            END DO
275         END DO
276         CALL iom_put( "saltc", rau0 * z2d )   ! vertically integrated salt content (PSU*kg/m2)
277      ENDIF
278      !
279      IF ( iom_use("eken") ) THEN
280         rke(:,:,jk) = 0._wp                               !      kinetic energy
281         DO jk = 1, jpkm1
282            DO jj = 2, jpjm1
283               DO ji = fs_2, fs_jpim1   ! vector opt.
284                  zztmp   = 1._wp / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) )
285                  zztmpx  = 0.5 * (  un(ji-1,jj,jk) * un(ji-1,jj,jk) * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)    &
286                     &             + un(ji  ,jj,jk) * un(ji  ,jj,jk) * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) )  &
287                     &          *  zztmp 
288                  !
289                  zztmpy  = 0.5 * (  vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)    &
290                     &             + vn(ji,jj  ,jk) * vn(ji,jj  ,jk) * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) )  &
291                     &          *  zztmp 
292                  !
293                  rke(ji,jj,jk) = 0.5_wp * ( zztmpx + zztmpy )
294                  !
295               ENDDO
296            ENDDO
297         ENDDO
298         CALL lbc_lnk( rke, 'T', 1. )
299         CALL iom_put( "eken", rke )           
300      ENDIF
301      !
302      CALL iom_put( "hdiv", hdivn )                  ! Horizontal divergence
303      !
304      IF( iom_use("u_masstr") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
305         z3d(:,:,jpk) = 0.e0
306         DO jk = 1, jpkm1
307            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)
308         END DO
309         CALL iom_put( "u_masstr", z3d )                  ! mass transport in i-direction
310      ENDIF
311     
312      IF( iom_use("u_heattr") ) THEN
313         z2d(:,:) = 0.e0 
314         DO jk = 1, jpkm1
315            DO jj = 2, jpjm1
316               DO ji = fs_2, fs_jpim1   ! vector opt.
317                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
318               END DO
319            END DO
320         END DO
321         CALL lbc_lnk( z2d, 'U', -1. )
322         CALL iom_put( "u_heattr", (0.5 * rcp) * z2d )    ! heat transport in i-direction
323      ENDIF
324
325      IF( iom_use("u_salttr") ) THEN
326         z2d(:,:) = 0.e0 
327         DO jk = 1, jpkm1
328            DO jj = 2, jpjm1
329               DO ji = fs_2, fs_jpim1   ! vector opt.
330                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
331               END DO
332            END DO
333         END DO
334         CALL lbc_lnk( z2d, 'U', -1. )
335         CALL iom_put( "u_salttr", 0.5 * z2d )            ! heat transport in i-direction
336      ENDIF
337
338     
339      IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
340         z3d(:,:,jpk) = 0.e0
341         DO jk = 1, jpkm1
342            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)
343         END DO
344         CALL iom_put( "v_masstr", z3d )                  ! mass transport in j-direction
345      ENDIF
346     
347      IF( iom_use("v_heattr") ) THEN
348         z2d(:,:) = 0.e0 
349         DO jk = 1, jpkm1
350            DO jj = 2, jpjm1
351               DO ji = fs_2, fs_jpim1   ! vector opt.
352                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
353               END DO
354            END DO
355         END DO
356         CALL lbc_lnk( z2d, 'V', -1. )
357         CALL iom_put( "v_heattr", (0.5 * rcp) * z2d )    !  heat transport in j-direction
358      ENDIF
359
360      IF( iom_use("v_salttr") ) THEN
361         z2d(:,:) = 0.e0 
362         DO jk = 1, jpkm1
363            DO jj = 2, jpjm1
364               DO ji = fs_2, fs_jpim1   ! vector opt.
365                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
366               END DO
367            END DO
368         END DO
369         CALL lbc_lnk( z2d, 'V', -1. )
370         CALL iom_put( "v_salttr", 0.5 * z2d )            !  heat transport in j-direction
371      ENDIF
372      !
373      CALL wrk_dealloc( jpi , jpj      , z2d )
374      CALL wrk_dealloc( jpi , jpj, jpk , z3d )
375      !
376      ! If we want tmb values
377
378      IF (ln_diatmb) THEN
379         CALL dia_tmb 
380      ENDIF
381      IF (ln_dia25h) THEN
382         CALL dia_25h( kt )
383      ENDIF
384
385      IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
386      !
387   END SUBROUTINE dia_wri
388
389#else
390   !!----------------------------------------------------------------------
391   !!   Default option                                  use IOIPSL  library
392   !!----------------------------------------------------------------------
393
394   SUBROUTINE dia_wri( kt )
395      !!---------------------------------------------------------------------
396      !!                  ***  ROUTINE dia_wri  ***
397      !!                   
398      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
399      !!      NETCDF format is used by default
400      !!
401      !! ** Method  :   At the beginning of the first time step (nit000),
402      !!      define all the NETCDF files and fields
403      !!      At each time step call histdef to compute the mean if ncessary
404      !!      Each nwrite time step, output the instantaneous or mean fields
405      !!----------------------------------------------------------------------
406      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
407      !
408      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
409      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
410      INTEGER  ::   inum = 11                                ! temporary logical unit
411      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
412      INTEGER  ::   ierr                                     ! error code return from allocation
413      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
414      INTEGER  ::   jn, ierror                               ! local integers
415      REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
416      !
417      REAL(wp), POINTER, DIMENSION(:,:)   :: zw2d       ! 2D workspace
418      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d       ! 3D workspace
419      !!----------------------------------------------------------------------
420      !
421      IF( nn_timing == 1 )   CALL timing_start('dia_wri')
422      !
423                             CALL wrk_alloc( jpi,jpj      , zw2d )
424      IF( .NOT.ln_linssh )   CALL wrk_alloc( jpi,jpj,jpk  , zw3d )
425      !
426      ! Output the initial state and forcings
427      IF( ninist == 1 ) THEN                       
428         CALL dia_wri_state( 'output.init', kt )
429         ninist = 0
430      ENDIF
431      !
432      ! 0. Initialisation
433      ! -----------------
434
435      ! local variable for debugging
436      ll_print = .FALSE.
437      ll_print = ll_print .AND. lwp
438
439      ! Define frequency of output and means
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 * rdt
443      clop = "inst("//TRIM(clop)//")"
444#else
445      zsto=rdt
446      clop = "ave("//TRIM(clop)//")"
447#endif
448      zout = nwrite * rdt
449      zmax = ( nitend - nit000 + 1 ) * rdt
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, rdt, 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, rdt, 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, rdt, 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, rdt, 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   SUBROUTINE dia_wri_state( cdfile_name, kt )
849      !!---------------------------------------------------------------------
850      !!                 ***  ROUTINE dia_wri_state  ***
851      !!       
852      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
853      !!      the instantaneous ocean state and forcing fields.
854      !!        Used to find errors in the initial state or save the last
855      !!      ocean state in case of abnormal end of a simulation
856      !!
857      !! ** Method  :   NetCDF files using ioipsl
858      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
859      !!      File 'output.abort.nc' is created in case of abnormal job end
860      !!----------------------------------------------------------------------
861      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
862      INTEGER           , INTENT( in ) ::   kt               ! ocean time-step index
863      !!
864      CHARACTER (len=32) :: clname
865      CHARACTER (len=40) :: clop
866      INTEGER  ::   id_i , nz_i, nh_i       
867      INTEGER, DIMENSION(1) ::   idex             ! local workspace
868      REAL(wp) ::   zsto, zout, zmax, zjulian
869      !!----------------------------------------------------------------------
870      !
871      ! 0. Initialisation
872      ! -----------------
873
874      ! Define name, frequency of output and means
875      clname = cdfile_name
876      IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
877      zsto = rdt
878      clop = "inst(x)"           ! no use of the mask value (require less cpu time)
879      zout = rdt
880      zmax = ( nitend - nit000 + 1 ) * rdt
881
882      IF(lwp) WRITE(numout,*)
883      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
884      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
885      IF(lwp) WRITE(numout,*) '                and named :', clname, '.nc'
886
887
888      ! 1. Define NETCDF files and fields at beginning of first time step
889      ! -----------------------------------------------------------------
890
891      ! Compute julian date from starting date of the run
892      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )         ! time axis
893      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
894      CALL histbeg( clname, jpi, glamt, jpj, gphit,   &
895          1, jpi, 1, jpj, nit000-1, zjulian, rdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
896      CALL histvert( id_i, "deptht", "Vertical T levels",   &    ! Vertical grid : gdept
897          "m", jpk, gdept_1d, nz_i, "down")
898
899      ! Declare all the output fields as NetCDF variables
900
901      CALL histdef( id_i, "vosaline", "Salinity"              , "PSU"    ,   &   ! salinity
902         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
903      CALL histdef( id_i, "votemper", "Temperature"           , "C"      ,   &   ! temperature
904         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
905      CALL histdef( id_i, "sossheig", "Sea Surface Height"    , "m"      ,   &  ! ssh
906         &          jpi, jpj, nh_i, 1  , 1, 1  , nz_i, 32, clop, zsto, zout )
907      CALL histdef( id_i, "vozocrtx", "Zonal Current"         , "m/s"    ,   &   ! zonal current
908         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
909      CALL histdef( id_i, "vomecrty", "Meridional Current"    , "m/s"    ,   &   ! meridonal current
910         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
911      CALL histdef( id_i, "vovecrtz", "Vertical Velocity"     , "m/s"    ,   &   ! vertical current
912         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
913         !
914      CALL histdef( id_i, "ahtu"    , "u-eddy diffusivity"    , "m2/s"    ,   &   ! zonal current
915         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
916      CALL histdef( id_i, "ahtv"    , "v-eddy diffusivity"    , "m2/s"    ,   &   ! meridonal current
917         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
918      CALL histdef( id_i, "ahmt"    , "t-eddy viscosity"      , "m2/s"    ,   &   ! zonal current
919         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
920      CALL histdef( id_i, "ahmf"    , "f-eddy viscosity"      , "m2/s"    ,   &   ! meridonal current
921         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
922         !
923      CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S",   &   ! net freshwater
924         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
925      CALL histdef( id_i, "sohefldo", "Net Downward Heat Flux", "W/m2"   ,   &   ! net heat flux
926         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
927      CALL histdef( id_i, "soshfldo", "Shortwave Radiation"   , "W/m2"   ,   &   ! solar flux
928         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
929      CALL histdef( id_i, "soicecov", "Ice fraction"          , "[0,1]"  ,   &   ! fr_i
930         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
931      CALL histdef( id_i, "sozotaux", "Zonal Wind Stress"     , "N/m2"   ,   &   ! i-wind stress
932         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
933      CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2"   ,   &   ! j-wind stress
934         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
935      IF( .NOT.ln_linssh ) THEN
936         CALL histdef( id_i, "vovvldep", "T point depth"         , "m"      , &   ! t-point depth
937            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
938         CALL histdef( id_i, "vovvle3t", "T point thickness"     , "m"      , &   ! t-point depth
939            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
940      ENDIF
941
942#if defined key_lim2
943      CALL lim_wri_state_2( kt, id_i, nh_i )
944#elif defined key_lim3
945      CALL lim_wri_state( kt, id_i, nh_i )
946#else
947      CALL histend( id_i, snc4chunks=snc4set )
948#endif
949
950      ! 2. Start writing data
951      ! ---------------------
952      ! idex(1) est utilise ssi l'avant dernier argument est diffferent de
953      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
954      ! donne le nombre d'elements, et idex la liste des indices a sortir
955      idex(1) = 1   ! init to avoid compil warning
956
957      ! Write all fields on T grid
958      CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature
959      CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity
960      CALL histwrite( id_i, "sossheig", kt, sshn             , jpi*jpj    , idex )    ! sea surface height
961      CALL histwrite( id_i, "vozocrtx", kt, un               , jpi*jpj*jpk, idex )    ! now i-velocity
962      CALL histwrite( id_i, "vomecrty", kt, vn               , jpi*jpj*jpk, idex )    ! now j-velocity
963      CALL histwrite( id_i, "vovecrtz", kt, wn               , jpi*jpj*jpk, idex )    ! now k-velocity
964      !
965      CALL histwrite( id_i, "ahtu"    , kt, ahtu             , jpi*jpj*jpk, idex )    ! aht at u-point
966      CALL histwrite( id_i, "ahtv"    , kt, ahtv             , jpi*jpj*jpk, idex )    !  -  at v-point
967      CALL histwrite( id_i, "ahmt"    , kt, ahmt             , jpi*jpj*jpk, idex )    ! ahm at t-point
968      CALL histwrite( id_i, "ahmf"    , kt, ahmf             , jpi*jpj*jpk, idex )    !  -  at f-point
969      !
970      CALL histwrite( id_i, "sowaflup", kt, emp-rnf          , jpi*jpj    , idex )    ! freshwater budget
971      CALL histwrite( id_i, "sohefldo", kt, qsr + qns        , jpi*jpj    , idex )    ! total heat flux
972      CALL histwrite( id_i, "soshfldo", kt, qsr              , jpi*jpj    , idex )    ! solar heat flux
973      CALL histwrite( id_i, "soicecov", kt, fr_i             , jpi*jpj    , idex )    ! ice fraction
974      CALL histwrite( id_i, "sozotaux", kt, utau             , jpi*jpj    , idex )    ! i-wind stress
975      CALL histwrite( id_i, "sometauy", kt, vtau             , jpi*jpj    , idex )    ! j-wind stress
976
977      IF(  .NOT.ln_linssh  ) THEN             
978         CALL histwrite( id_i, "vovvldep", kt, gdept_n(:,:,:), jpi*jpj*jpk, idex )!  T-cell depth
979         CALL histwrite( id_i, "vovvle3t", kt, e3t_n (:,:,:) , jpi*jpj*jpk, idex )!  T-cell thickness 
980      END IF 
981      ! 3. Close the file
982      ! -----------------
983      CALL histclo( id_i )
984#if ! defined key_iomput
985      IF( ninist /= 1  ) THEN
986         CALL histclo( nid_T )
987         CALL histclo( nid_U )
988         CALL histclo( nid_V )
989         CALL histclo( nid_W )
990      ENDIF
991#endif
992      !
993   END SUBROUTINE dia_wri_state
994
995   !!======================================================================
996END MODULE diawri
Note: See TracBrowser for help on using the repository browser.