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 @ 7931

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

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