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/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/tests/DOME/MY_SRC – NEMO

source: NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/tests/DOME/MY_SRC/diawri.F90 @ 14552

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

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