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

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

source: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/DIA/diawri.F90 @ 8367

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

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