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_r13296_HPC-07_mocavero_mpi3/src/SWE – NEMO

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/SWE/diawri.F90 @ 13630

Last change on this file since 13630 was 13630, checked in by mocavero, 4 years ago

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