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

Last change on this file since 13565 was 13565, checked in by jchanut, 5 months ago

#2222, 1) Added parent bathymetry volume consistency check 2) Added velocity extrapolation in update 3) Corrected bdy issue #2519

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