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 NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DIA – NEMO

source: NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DIA/diawri.F90 @ 10164

Last change on this file since 10164 was 10164, checked in by jchanut, 6 years ago

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