source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/tests/CANAL/MY_SRC/diawri.F90 @ 12756

Last change on this file since 12756 was 12756, checked in by clem, 9 months ago

fix xml again

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