source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/CANAL/MY_SRC/diawri.F90 @ 10410

Last change on this file since 10410 was 10410, checked in by smasson, 21 months ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: mssing commit in test following [10380], see #2133

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