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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diawri.F90 in NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/OCE/DIA – NEMO

source: NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/OCE/DIA/diawri.F90 @ 13942

Last change on this file since 13942 was 13942, checked in by jchanut, 4 years ago

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