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/cfgs/penzAM98/MY_SRC – NEMO

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/penzAM98/MY_SRC/diawri.F90 @ 13562

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

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