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/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA – NEMO

source: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/DIA/diawri.F90 @ 12928

Last change on this file since 12928 was 12928, checked in by smueller, 4 years ago

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