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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diawri.F90 in branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 8093

Last change on this file since 8093 was 8093, checked in by gm, 7 years ago

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