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/UKMO/NEMO_4.0_FKOSM/src/OCE/DIA – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/DIA/diawri.F90 @ 12387

Last change on this file since 12387 was 12387, checked in by cguiavarch, 4 years ago

corrections to tramle & zdfosm to avoid division by 0 and restored OSMOSIS variables in output.abort

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