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

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

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diawri.F90 @ 12622

Last change on this file since 12622 was 12622, checked in by techene, 4 years ago

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