source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OCE_SRC/DIA/diawri.F90 @ 9570

Last change on this file since 9570 was 9570, checked in by nicolasmartin, 3 years ago

Global renaming for core routines (./NEMO)

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