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_r12527_Gurvan_ShallowWater/src/SWE – NEMO

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/diawri.F90 @ 13151

Last change on this file since 13151 was 13151, checked in by gm, 4 years ago

result from merge with qco r12983

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