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/releases/release-4.0-HEAD/src/OCE/DIA – NEMO

source: NEMO/releases/release-4.0-HEAD/src/OCE/DIA/diawri.F90 @ 12494

Last change on this file since 12494 was 12494, checked in by smasson, 4 years ago

release-4.0-HEAD: minor optimisation in diawri, see #2393
  • Property svn:keywords set to Id
File size: 49.1 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 dom_oce        ! ocean space and time domain
29   USE phycst         ! physical constants
30   USE dianam         ! build name of file (routine)
31   USE diahth         ! thermocline diagnostics
32   USE dynadv   , ONLY: ln_dynadv_vec
33   USE icb_oce        ! Icebergs
34   USE icbdia         ! Iceberg budgets
35   USE ldftra         ! lateral physics: eddy diffusivity coef.
36   USE ldfdyn         ! lateral physics: eddy viscosity   coef.
37   USE sbc_oce        ! Surface boundary condition: ocean fields
38   USE sbc_ice        ! Surface boundary condition: ice fields
39   USE sbcssr         ! restoring term toward SST/SSS climatology
40   USE sbcwave        ! wave parameters
41   USE wet_dry        ! wetting and drying
42   USE zdf_oce        ! ocean vertical physics
43   USE zdfdrg         ! ocean vertical physics: top/bottom friction
44   USE zdfmxl         ! mixed layer
45   !
46   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
47   USE in_out_manager ! I/O manager
48   USE dia25h         ! 25h Mean output
49   USE iom            !
50   USE ioipsl         !
51
52#if defined key_si3
53   USE ice 
54   USE icewri 
55#endif
56   USE lib_mpp         ! MPP library
57   USE timing          ! preformance summary
58   USE diurnal_bulk    ! diurnal warm layer
59   USE cool_skin       ! Cool skin
60
61   IMPLICIT NONE
62   PRIVATE
63
64   PUBLIC   dia_wri                 ! routines called by step.F90
65   PUBLIC   dia_wri_state
66   PUBLIC   dia_wri_alloc           ! Called by nemogcm module
67
68   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
69   INTEGER ::          nb_T              , ndim_bT   ! grid_T file
70   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
71   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
72   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file
73   INTEGER ::   ndex(1)                              ! ???
74   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
75   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
76   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
77
78   !! * Substitutions
79#  include "vectopt_loop_substitute.h90"
80   !!----------------------------------------------------------------------
81   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
82   !! $Id$
83   !! Software governed by the CeCILL license (see ./LICENSE)
84   !!----------------------------------------------------------------------
85CONTAINS
86
87#if defined key_iomput
88   !!----------------------------------------------------------------------
89   !!   'key_iomput'                                        use IOM library
90   !!----------------------------------------------------------------------
91   INTEGER FUNCTION dia_wri_alloc()
92      !
93      dia_wri_alloc = 0
94      !
95   END FUNCTION dia_wri_alloc
96
97   
98   SUBROUTINE dia_wri( kt )
99      !!---------------------------------------------------------------------
100      !!                  ***  ROUTINE dia_wri  ***
101      !!                   
102      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
103      !!      NETCDF format is used by default
104      !!
105      !! ** Method  :  use iom_put
106      !!----------------------------------------------------------------------
107      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
108      !!
109      INTEGER ::   ji, jj, jk       ! dummy loop indices
110      INTEGER ::   ikbot            ! local integer
111      REAL(wp)::   zztmp , zztmpx   ! local scalar
112      REAL(wp)::   zztmp2, zztmpy   !   -      -
113      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace
114      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace
115      !!----------------------------------------------------------------------
116      !
117      IF( ln_timing )   CALL timing_start('dia_wri')
118      !
119      ! Output the initial state and forcings
120      IF( ninist == 1 ) THEN                       
121         CALL dia_wri_state( 'output.init' )
122         ninist = 0
123      ENDIF
124
125      ! Output of initial vertical scale factor
126      CALL iom_put("e3t_0", e3t_0(:,:,:) )
127      CALL iom_put("e3u_0", e3u_0(:,:,:) )
128      CALL iom_put("e3v_0", e3v_0(:,:,:) )
129      !
130      CALL iom_put( "e3t" , e3t_n(:,:,:) )
131      CALL iom_put( "e3u" , e3u_n(:,:,:) )
132      CALL iom_put( "e3v" , e3v_n(:,:,:) )
133      CALL iom_put( "e3w" , e3w_n(:,:,:) )
134      IF( iom_use("e3tdef") )   &
135         CALL iom_put( "e3tdef"  , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
136
137      IF( ll_wd ) THEN
138         CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying)
139      ELSE
140         CALL iom_put( "ssh" , sshn )              ! sea surface height
141      ENDIF
142
143      IF( iom_use("wetdep") )   &                  ! wet depth
144         CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) )
145     
146      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
147      CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature
148      IF ( iom_use("sbt") ) THEN
149         DO jj = 1, jpj
150            DO ji = 1, jpi
151               ikbot = mbkt(ji,jj)
152               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem)
153            END DO
154         END DO
155         CALL iom_put( "sbt", z2d )                ! bottom temperature
156      ENDIF
157     
158      CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity
159      CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity
160      IF ( iom_use("sbs") ) THEN
161         DO jj = 1, jpj
162            DO ji = 1, jpi
163               ikbot = mbkt(ji,jj)
164               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal)
165            END DO
166         END DO
167         CALL iom_put( "sbs", z2d )                ! bottom salinity
168      ENDIF
169
170      IF ( iom_use("taubot") ) THEN                ! bottom stress
171         zztmp = rau0 * 0.25
172         z2d(:,:) = 0._wp
173         DO jj = 2, jpjm1
174            DO ji = fs_2, fs_jpim1   ! vector opt.
175               zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * un(ji  ,jj,mbku(ji  ,jj))  )**2   &
176                  &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj))  )**2   &
177                  &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vn(ji,jj  ,mbkv(ji,jj  ))  )**2   &
178                  &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1))  )**2
179               z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 
180               !
181            END DO
182         END DO
183         CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
184         CALL iom_put( "taubot", z2d )           
185      ENDIF
186         
187      CALL iom_put( "uoce", un(:,:,:) )            ! 3D i-current
188      CALL iom_put(  "ssu", un(:,:,1) )            ! surface i-current
189      IF ( iom_use("sbu") ) THEN
190         DO jj = 1, jpj
191            DO ji = 1, jpi
192               ikbot = mbku(ji,jj)
193               z2d(ji,jj) = un(ji,jj,ikbot)
194            END DO
195         END DO
196         CALL iom_put( "sbu", z2d )                ! bottom i-current
197      ENDIF
198     
199      CALL iom_put( "voce", vn(:,:,:) )            ! 3D j-current
200      CALL iom_put(  "ssv", vn(:,:,1) )            ! surface j-current
201      IF ( iom_use("sbv") ) THEN
202         DO jj = 1, jpj
203            DO ji = 1, jpi
204               ikbot = mbkv(ji,jj)
205               z2d(ji,jj) = vn(ji,jj,ikbot)
206            END DO
207         END DO
208         CALL iom_put( "sbv", z2d )                ! bottom j-current
209      ENDIF
210
211      IF( ln_zad_Aimp ) wn = wn + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output
212      !
213      CALL iom_put( "woce", wn )                   ! vertical velocity
214      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
215         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
216         z2d(:,:) = rau0 * e1e2t(:,:)
217         DO jk = 1, jpk
218            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
219         END DO
220         CALL iom_put( "w_masstr" , z3d ) 
221         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
222      ENDIF
223      !
224      IF( ln_zad_Aimp ) wn = wn - wi               ! Remove implicit part of vertical velocity that was added for diagnostic output
225
226      CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef.
227      CALL iom_put( "avs" , avs )                  ! S vert. eddy diff. coef.
228      CALL iom_put( "avm" , avm )                  ! T vert. eddy visc. coef.
229
230      IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
231      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
232
233      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
234         DO jj = 2, jpjm1                                    ! sst gradient
235            DO ji = fs_2, fs_jpim1   ! vector opt.
236               zztmp  = tsn(ji,jj,1,jp_tem)
237               zztmpx = ( tsn(ji+1,jj,1,jp_tem) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,1,jp_tem) ) * r1_e1u(ji-1,jj)
238               zztmpy = ( tsn(ji,jj+1,1,jp_tem) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - tsn(ji  ,jj-1,1,jp_tem) ) * r1_e2v(ji,jj-1)
239               z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
240                  &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
241            END DO
242         END DO
243         CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
244         CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient
245         z2d(:,:) = SQRT( z2d(:,:) )
246         CALL iom_put( "sstgrad" ,  z2d )          ! module of sst gradient
247      ENDIF
248         
249      ! heat and salt contents
250      IF( iom_use("heatc") ) THEN
251         z2d(:,:)  = 0._wp 
252         DO jk = 1, jpkm1
253            DO jj = 1, jpj
254               DO ji = 1, jpi
255                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
256               END DO
257            END DO
258         END DO
259         CALL iom_put( "heatc", rau0_rcp * z2d )   ! vertically integrated heat content (J/m2)
260      ENDIF
261
262      IF( iom_use("saltc") ) THEN
263         z2d(:,:)  = 0._wp 
264         DO jk = 1, jpkm1
265            DO jj = 1, jpj
266               DO ji = 1, jpi
267                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
268               END DO
269            END DO
270         END DO
271         CALL iom_put( "saltc", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
272      ENDIF
273      !
274      IF ( iom_use("eken") ) THEN
275         z3d(:,:,jpk) = 0._wp 
276         DO jk = 1, jpkm1
277            DO jj = 2, jpjm1
278               DO ji = fs_2, fs_jpim1   ! vector opt.
279                  zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
280                  z3d(ji,jj,jk) = zztmp * (  un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)   &
281                     &                     + un(ji  ,jj,jk)**2 * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)   &
282                     &                     + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)   &
283                     &                     + vn(ji,jj  ,jk)**2 * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk)   )
284               END DO
285            END DO
286         END DO
287         CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
288         CALL iom_put( "eken", z3d )                 ! kinetic energy
289      ENDIF
290      !
291      CALL iom_put( "hdiv", hdivn )                  ! Horizontal divergence
292      !
293      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
294         z3d(:,:,jpk) = 0.e0
295         z2d(:,:) = 0.e0
296         DO jk = 1, jpkm1
297            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)
298            z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
299         END DO
300         CALL iom_put( "u_masstr"     , z3d )         ! mass transport in i-direction
301         CALL iom_put( "u_masstr_vint", z2d )         ! mass transport in i-direction vertical sum
302      ENDIF
303     
304      IF( iom_use("u_heattr") ) THEN
305         z2d(:,:) = 0._wp 
306         DO jk = 1, jpkm1
307            DO jj = 2, jpjm1
308               DO ji = fs_2, fs_jpim1   ! vector opt.
309                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
310               END DO
311            END DO
312         END DO
313         CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
314         CALL iom_put( "u_heattr", 0.5*rcp * z2d )    ! heat transport in i-direction
315      ENDIF
316
317      IF( iom_use("u_salttr") ) THEN
318         z2d(:,:) = 0.e0 
319         DO jk = 1, jpkm1
320            DO jj = 2, jpjm1
321               DO ji = fs_2, fs_jpim1   ! vector opt.
322                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
323               END DO
324            END DO
325         END DO
326         CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
327         CALL iom_put( "u_salttr", 0.5 * z2d )        ! heat transport in i-direction
328      ENDIF
329
330     
331      IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
332         z3d(:,:,jpk) = 0.e0
333         DO jk = 1, jpkm1
334            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)
335         END DO
336         CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction
337      ENDIF
338     
339      IF( iom_use("v_heattr") ) THEN
340         z2d(:,:) = 0.e0 
341         DO jk = 1, jpkm1
342            DO jj = 2, jpjm1
343               DO ji = fs_2, fs_jpim1   ! vector opt.
344                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
345               END DO
346            END DO
347         END DO
348         CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
349         CALL iom_put( "v_heattr", 0.5*rcp * z2d )    !  heat transport in j-direction
350      ENDIF
351
352      IF( iom_use("v_salttr") ) THEN
353         z2d(:,:) = 0._wp 
354         DO jk = 1, jpkm1
355            DO jj = 2, jpjm1
356               DO ji = fs_2, fs_jpim1   ! vector opt.
357                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
358               END DO
359            END DO
360         END DO
361         CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
362         CALL iom_put( "v_salttr", 0.5 * z2d )        !  heat transport in j-direction
363      ENDIF
364
365      IF( iom_use("tosmint") ) THEN
366         z2d(:,:) = 0._wp
367         DO jk = 1, jpkm1
368            DO jj = 2, jpjm1
369               DO ji = fs_2, fs_jpim1   ! vector opt.
370                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem)
371               END DO
372            END DO
373         END DO
374         CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
375         CALL iom_put( "tosmint", rau0 * z2d )        ! Vertical integral of temperature
376      ENDIF
377      IF( iom_use("somint") ) THEN
378         z2d(:,:)=0._wp
379         DO jk = 1, jpkm1
380            DO jj = 2, jpjm1
381               DO ji = fs_2, fs_jpim1   ! vector opt.
382                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
383               END DO
384            END DO
385         END DO
386         CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
387         CALL iom_put( "somint", rau0 * z2d )         ! Vertical integral of salinity
388      ENDIF
389
390      CALL iom_put( "bn2", rn2 )                      ! Brunt-Vaisala buoyancy frequency (N^2)
391      !
392         
393      IF (ln_dia25h)   CALL dia_25h( kt )             ! 25h averaging
394
395      IF( ln_timing )   CALL timing_stop('dia_wri')
396      !
397   END SUBROUTINE dia_wri
398
399#else
400   !!----------------------------------------------------------------------
401   !!   Default option                                  use IOIPSL  library
402   !!----------------------------------------------------------------------
403
404   INTEGER FUNCTION dia_wri_alloc()
405      !!----------------------------------------------------------------------
406      INTEGER, DIMENSION(2) :: ierr
407      !!----------------------------------------------------------------------
408      IF( nn_write == -1 ) THEN
409         dia_wri_alloc = 0
410      ELSE   
411         ierr = 0
412         ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
413            &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
414            &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
415         !
416         dia_wri_alloc = MAXVAL(ierr)
417         CALL mpp_sum( 'diawri', dia_wri_alloc )
418         !
419      ENDIF
420      !
421   END FUNCTION dia_wri_alloc
422
423   
424   SUBROUTINE dia_wri( kt )
425      !!---------------------------------------------------------------------
426      !!                  ***  ROUTINE dia_wri  ***
427      !!                   
428      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
429      !!      NETCDF format is used by default
430      !!
431      !! ** Method  :   At the beginning of the first time step (nit000),
432      !!      define all the NETCDF files and fields
433      !!      At each time step call histdef to compute the mean if ncessary
434      !!      Each nn_write time step, output the instantaneous or mean fields
435      !!----------------------------------------------------------------------
436      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
437      !
438      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
439      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
440      INTEGER  ::   inum = 11                                ! temporary logical unit
441      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
442      INTEGER  ::   ierr                                     ! error code return from allocation
443      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
444      INTEGER  ::   jn, ierror                               ! local integers
445      REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
446      !
447      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace
448      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace
449      !!----------------------------------------------------------------------
450      !
451      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==!
452         CALL dia_wri_state( 'output.init' )
453         ninist = 0
454      ENDIF
455      !
456      IF( nn_write == -1 )   RETURN   ! we will never do any output
457      !
458      IF( ln_timing )   CALL timing_start('dia_wri')
459      !
460      ! 0. Initialisation
461      ! -----------------
462
463      ll_print = .FALSE.                  ! local variable for debugging
464      ll_print = ll_print .AND. lwp
465
466      ! Define frequency of output and means
467      clop = "x"         ! no use of the mask value (require less cpu time and otherwise the model crashes)
468#if defined key_diainstant
469      zsto = nn_write * rdt
470      clop = "inst("//TRIM(clop)//")"
471#else
472      zsto=rdt
473      clop = "ave("//TRIM(clop)//")"
474#endif
475      zout = nn_write * rdt
476      zmax = ( nitend - nit000 + 1 ) * rdt
477
478      ! Define indices of the horizontal output zoom and vertical limit storage
479      iimi = 1      ;      iima = jpi
480      ijmi = 1      ;      ijma = jpj
481      ipk = jpk
482
483      ! define time axis
484      it = kt
485      itmod = kt - nit000 + 1
486
487
488      ! 1. Define NETCDF files and fields at beginning of first time step
489      ! -----------------------------------------------------------------
490
491      IF( kt == nit000 ) THEN
492
493         ! Define the NETCDF files (one per grid)
494
495         ! Compute julian date from starting date of the run
496         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
497         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
498         IF(lwp)WRITE(numout,*)
499         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
500            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
501         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
502                                 ' limit storage in depth = ', ipk
503
504         ! WRITE root name in date.file for use by postpro
505         IF(lwp) THEN
506            CALL dia_nam( clhstnam, nn_write,' ' )
507            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
508            WRITE(inum,*) clhstnam
509            CLOSE(inum)
510         ENDIF
511
512         ! Define the T grid FILE ( nid_T )
513
514         CALL dia_nam( clhstnam, nn_write, 'grid_T' )
515         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
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, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
519         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
520            &           "m", ipk, gdept_1d, nz_T, "down" )
521         !                                                            ! Index of ocean points
522         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
523         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
524         !
525         IF( ln_icebergs ) THEN
526            !
527            !! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
528            !! that routine is called from nemogcm, so do it here immediately before its needed
529            ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
530            CALL mpp_sum( 'diawri', ierror )
531            IF( ierror /= 0 ) THEN
532               CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
533               RETURN
534            ENDIF
535            !
536            !! iceberg vertical coordinate is class number
537            CALL histvert( nid_T, "class", "Iceberg class",      &  ! Vertical grid: class
538               &           "number", nclasses, class_num, nb_T )
539            !
540            !! each class just needs the surface index pattern
541            ndim_bT = 3
542            DO jn = 1,nclasses
543               ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
544            ENDDO
545            !
546         ENDIF
547
548         ! Define the U grid FILE ( nid_U )
549
550         CALL dia_nam( clhstnam, nn_write, 'grid_U' )
551         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
552         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
553            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
554            &          nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
555         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
556            &           "m", ipk, gdept_1d, nz_U, "down" )
557         !                                                            ! Index of ocean points
558         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
559         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
560
561         ! Define the V grid FILE ( nid_V )
562
563         CALL dia_nam( clhstnam, nn_write, 'grid_V' )                   ! filename
564         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
565         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
566            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
567            &          nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
568         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
569            &          "m", ipk, gdept_1d, nz_V, "down" )
570         !                                                            ! Index of ocean points
571         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
572         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
573
574         ! Define the W grid FILE ( nid_W )
575
576         CALL dia_nam( clhstnam, nn_write, 'grid_W' )                   ! filename
577         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
578         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
579            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
580            &          nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
581         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
582            &          "m", ipk, gdepw_1d, nz_W, "down" )
583
584
585         ! Declare all the output fields as NETCDF variables
586
587         !                                                                                      !!! nid_T : 3D
588         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
589            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
590         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
591            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
592         IF(  .NOT.ln_linssh  ) THEN
593            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n
594            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
595            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t_n
596            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
597            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t_n
598            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
599         ENDIF
600         !                                                                                      !!! nid_T : 2D
601         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
602            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
603         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
604            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
605         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
606            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
607         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
608            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
609         CALL histdef( nid_T, "sorunoff", "River runoffs"                      , "Kg/m2/s",   &  ! runoffs
610            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
611         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
612            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
613         IF(  ln_linssh  ) THEN
614            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem)
615            &                                                                  , "KgC/m2/s",  &  ! sosst_cd
616            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
617            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * tsn(:,:,1,jp_sal)
618            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
619            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
620         ENDIF
621         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
622            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
623         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
624            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
625         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
626            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
627         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
628            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
629         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
630            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
631         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
632            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
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"    ,   &  ! un
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"    ,   &  ! vn
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"    ,   &  ! wn
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
756      IF( .NOT.ln_linssh ) THEN
757         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
758         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
759         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
760         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
761      ELSE
762         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature
763         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity
764         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature
765         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
766      ENDIF
767      IF( .NOT.ln_linssh ) THEN
768         zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
769         CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
770         CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
771         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
772      ENDIF
773      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
774      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
775      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
776      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
777                                                                                  ! (includes virtual salt flux beneath ice
778                                                                                  ! in linear free surface case)
779      IF( ln_linssh ) THEN
780         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
781         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
782         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
783         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
784      ENDIF
785      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
786      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
787      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
788      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
789      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
790      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
791!
792      IF( ln_icebergs ) THEN
793         !
794         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT ) 
795         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
796         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT ) 
797         !
798         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
799         !
800         IF( ln_bergdia ) THEN
801            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   ) 
802            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   ) 
803            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   ) 
804            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   ) 
805            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   ) 
806            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   ) 
807            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   ) 
808            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   ) 
809            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   ) 
810            !
811            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
812         ENDIF
813      ENDIF
814
815      IF( ln_ssr ) THEN
816         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
817         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
818         zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
819         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
820      ENDIF
821!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
822!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
823
824#if defined key_diahth
825      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
826      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
827      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
828      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
829#endif
830
831      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
832      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
833
834      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
835      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
836
837      IF( ln_zad_Aimp ) THEN
838         CALL histwrite( nid_W, "vovecrtz", it, wn + wi     , ndim_T, ndex_T )    ! vert. current
839      ELSE
840         CALL histwrite( nid_W, "vovecrtz", it, wn          , ndim_T, ndex_T )    ! vert. current
841      ENDIF
842      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
843      CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
844      IF( ln_zdfddm ) THEN
845         CALL histwrite( nid_W, "voddmavs", it, avs         , ndim_T, ndex_T )    ! S vert. eddy diff. coef.
846      ENDIF
847
848      IF( ln_wave .AND. ln_sdw ) THEN
849         CALL histwrite( nid_U, "sdzocrtx", it, usd         , ndim_U , ndex_U )    ! i-StokesDrift-current
850         CALL histwrite( nid_V, "sdmecrty", it, vsd         , ndim_V , ndex_V )    ! j-StokesDrift-current
851         CALL histwrite( nid_W, "sdvecrtz", it, wsd         , ndim_T , ndex_T )    ! StokesDrift vert. current
852      ENDIF
853
854      ! 3. Close all files
855      ! ---------------------------------------
856      IF( kt == nitend ) THEN
857         CALL histclo( nid_T )
858         CALL histclo( nid_U )
859         CALL histclo( nid_V )
860         CALL histclo( nid_W )
861      ENDIF
862      !
863      IF( ln_timing )   CALL timing_stop('dia_wri')
864      !
865   END SUBROUTINE dia_wri
866#endif
867
868   SUBROUTINE dia_wri_state( cdfile_name )
869      !!---------------------------------------------------------------------
870      !!                 ***  ROUTINE dia_wri_state  ***
871      !!       
872      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
873      !!      the instantaneous ocean state and forcing fields.
874      !!        Used to find errors in the initial state or save the last
875      !!      ocean state in case of abnormal end of a simulation
876      !!
877      !! ** Method  :   NetCDF files using ioipsl
878      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
879      !!      File 'output.abort.nc' is created in case of abnormal job end
880      !!----------------------------------------------------------------------
881      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
882      !!
883      INTEGER :: inum
884      !!----------------------------------------------------------------------
885      !
886      IF(lwp) WRITE(numout,*)
887      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
888      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
889      IF(lwp) WRITE(numout,*) '                and named :', cdfile_name, '...nc'
890
891#if defined key_si3
892     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
893#else
894     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
895#endif
896
897      CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) )    ! now temperature
898      CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) )    ! now salinity
899      CALL iom_rstput( 0, 0, inum, 'sossheig', sshn              )    ! sea surface height
900      CALL iom_rstput( 0, 0, inum, 'vozocrtx', un                )    ! now i-velocity
901      CALL iom_rstput( 0, 0, inum, 'vomecrty', vn                )    ! now j-velocity
902      IF( ln_zad_Aimp ) THEN
903         CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn + wi        )    ! now k-velocity
904      ELSE
905         CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn             )    ! now k-velocity
906      ENDIF
907      IF( ALLOCATED(ahtu) ) THEN
908         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point
909         CALL iom_rstput( 0, 0, inum,  'ahtv', ahtv              )    ! aht at v-point
910      ENDIF
911      IF( ALLOCATED(ahmt) ) THEN
912         CALL iom_rstput( 0, 0, inum,  'ahmt', ahmt              )    ! ahmt at u-point
913         CALL iom_rstput( 0, 0, inum,  'ahmf', ahmf              )    ! ahmf at v-point
914      ENDIF
915      CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf         )    ! freshwater budget
916      CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns         )    ! total heat flux
917      CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr               )    ! solar heat flux
918      CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i              )    ! ice fraction
919      CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
920      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
921      IF(  .NOT.ln_linssh  ) THEN             
922         CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept_n        )    !  T-cell depth
923         CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n          )    !  T-cell thickness 
924      END IF
925      IF( ln_wave .AND. ln_sdw ) THEN
926         CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd            )    ! now StokesDrift i-velocity
927         CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd            )    ! now StokesDrift j-velocity
928         CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity
929      ENDIF
930 
931#if defined key_si3
932      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid
933         CALL ice_wri_state( inum )
934      ENDIF
935#endif
936      !
937      CALL iom_close( inum )
938      !
939   END SUBROUTINE dia_wri_state
940
941   !!======================================================================
942END MODULE diawri
Note: See TracBrowser for help on using the repository browser.