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

source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 8316

Last change on this file since 8316 was 8316, checked in by clem, 7 years ago

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