New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diawri.F90 in NEMO/branches/2019/ENHANCE-02_ISF_nemo_TEST_MERGE/src/OCE/DIA – NEMO

source: NEMO/branches/2019/ENHANCE-02_ISF_nemo_TEST_MERGE/src/OCE/DIA/diawri.F90 @ 11970

Last change on this file since 11970 was 11970, checked in by davestorkey, 4 years ago

2019/ENHANCE-02_ISF_nemo_TEST_MERGE : copy changes from Pierre's branch.

  • Property svn:keywords set to Id
File size: 50.8 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
29   USE isfcpl
30   USE dom_oce        ! ocean space and time domain
31   USE phycst         ! physical constants
32   USE dianam         ! build name of file (routine)
33   USE diahth         ! thermocline diagnostics
34   USE dynadv   , ONLY: ln_dynadv_vec
35   USE icb_oce        ! Icebergs
36   USE icbdia         ! Iceberg budgets
37   USE ldftra         ! lateral physics: eddy diffusivity coef.
38   USE ldfdyn         ! lateral physics: eddy viscosity   coef.
39   USE sbc_oce        ! Surface boundary condition: ocean fields
40   USE sbc_ice        ! Surface boundary condition: ice fields
41   USE sbcssr         ! restoring term toward SST/SSS climatology
42   USE sbcwave        ! wave parameters
43   USE wet_dry        ! wetting and drying
44   USE zdf_oce        ! ocean vertical physics
45   USE zdfdrg         ! ocean vertical physics: top/bottom friction
46   USE zdfmxl         ! mixed layer
47   !
48   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
49   USE in_out_manager ! I/O manager
50   USE diatmb         ! Top,middle,bottom output
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 diurnal_bulk    ! diurnal warm layer
62   USE cool_skin       ! 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
71   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
72   INTEGER ::          nb_T              , ndim_bT   ! grid_T file
73   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
74   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
75   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file
76   INTEGER ::   ndex(1)                              ! ???
77   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
78   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
79   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
80
81   !! * Substitutions
82#  include "vectopt_loop_substitute.h90"
83   !!----------------------------------------------------------------------
84   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
85   !! $Id$
86   !! Software governed by the CeCILL license (see ./LICENSE)
87   !!----------------------------------------------------------------------
88CONTAINS
89
90#if defined key_iomput
91   !!----------------------------------------------------------------------
92   !!   'key_iomput'                                        use IOM library
93   !!----------------------------------------------------------------------
94   INTEGER FUNCTION dia_wri_alloc()
95      !
96      dia_wri_alloc = 0
97      !
98   END FUNCTION dia_wri_alloc
99
100   
101   SUBROUTINE dia_wri( kt )
102      !!---------------------------------------------------------------------
103      !!                  ***  ROUTINE dia_wri  ***
104      !!                   
105      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
106      !!      NETCDF format is used by default
107      !!
108      !! ** Method  :  use iom_put
109      !!----------------------------------------------------------------------
110      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
111      !!
112      INTEGER ::   ji, jj, jk       ! dummy loop indices
113      INTEGER ::   ikbot            ! local integer
114      REAL(wp)::   zztmp , zztmpx   ! local scalar
115      REAL(wp)::   zztmp2, zztmpy   !   -      -
116      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace
117      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace
118      !!----------------------------------------------------------------------
119      !
120      IF( ln_timing )   CALL timing_start('dia_wri')
121      !
122      ! Output the initial state and forcings
123      IF( ninist == 1 ) THEN                       
124         CALL dia_wri_state( 'output.init' )
125         ninist = 0
126      ENDIF
127
128      ! Output of initial vertical scale factor
129      CALL iom_put("e3t_0", e3t_0(:,:,:) )
130      CALL iom_put("e3u_0", e3u_0(:,:,:) )
131      CALL iom_put("e3v_0", e3v_0(:,:,:) )
132      !
133      CALL iom_put( "e3t" , e3t_n(:,:,:) )
134      CALL iom_put( "e3u" , e3u_n(:,:,:) )
135      CALL iom_put( "e3v" , e3v_n(:,:,:) )
136      CALL iom_put( "e3w" , e3w_n(:,:,:) )
137      IF( iom_use("e3tdef") )   &
138         CALL iom_put( "e3tdef"  , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
139
140      IF( ll_wd ) THEN
141         CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying)
142      ELSE
143         CALL iom_put( "ssh" , sshn )              ! sea surface height
144      ENDIF
145
146      IF( iom_use("wetdep") )   &                  ! wet depth
147         CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) )
148     
149      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
150      CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature
151      IF ( iom_use("sbt") ) THEN
152         DO jj = 1, jpj
153            DO ji = 1, jpi
154               ikbot = mbkt(ji,jj)
155               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem)
156            END DO
157         END DO
158         CALL iom_put( "sbt", z2d )                ! bottom temperature
159      ENDIF
160     
161      CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity
162      CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity
163      IF ( iom_use("sbs") ) THEN
164         DO jj = 1, jpj
165            DO ji = 1, jpi
166               ikbot = mbkt(ji,jj)
167               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal)
168            END DO
169         END DO
170         CALL iom_put( "sbs", z2d )                ! bottom salinity
171      ENDIF
172
173      IF ( iom_use("taubot") ) THEN                ! bottom stress
174         zztmp = rau0 * 0.25
175         z2d(:,:) = 0._wp
176         DO jj = 2, jpjm1
177            DO ji = fs_2, fs_jpim1   ! vector opt.
178               zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * un(ji  ,jj,mbku(ji  ,jj))  )**2   &
179                  &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj))  )**2   &
180                  &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vn(ji,jj  ,mbkv(ji,jj  ))  )**2   &
181                  &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1))  )**2
182               z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 
183               !
184            END DO
185         END DO
186         CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
187         CALL iom_put( "taubot", z2d )           
188      ENDIF
189         
190      CALL iom_put( "uoce", un(:,:,:) )            ! 3D i-current
191      CALL iom_put(  "ssu", un(:,:,1) )            ! surface i-current
192      IF ( iom_use("sbu") ) THEN
193         DO jj = 1, jpj
194            DO ji = 1, jpi
195               ikbot = mbku(ji,jj)
196               z2d(ji,jj) = un(ji,jj,ikbot)
197            END DO
198         END DO
199         CALL iom_put( "sbu", z2d )                ! bottom i-current
200      ENDIF
201     
202      CALL iom_put( "voce", vn(:,:,:) )            ! 3D j-current
203      CALL iom_put(  "ssv", vn(:,:,1) )            ! surface j-current
204      IF ( iom_use("sbv") ) THEN
205         DO jj = 1, jpj
206            DO ji = 1, jpi
207               ikbot = mbkv(ji,jj)
208               z2d(ji,jj) = vn(ji,jj,ikbot)
209            END DO
210         END DO
211         CALL iom_put( "sbv", z2d )                ! bottom j-current
212      ENDIF
213
214      IF( ln_zad_Aimp ) wn = wn + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output
215      !
216      CALL iom_put( "woce", wn )                   ! vertical velocity
217      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
218         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
219         z2d(:,:) = rau0 * e1e2t(:,:)
220         DO jk = 1, jpk
221            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
222         END DO
223         CALL iom_put( "w_masstr" , z3d ) 
224         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
225      ENDIF
226      !
227      IF( ln_zad_Aimp ) wn = wn - wi               ! Remove implicit part of vertical velocity that was added for diagnostic output
228
229      CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef.
230      CALL iom_put( "avs" , avs )                  ! S vert. eddy diff. coef.
231      CALL iom_put( "avm" , avm )                  ! T vert. eddy visc. coef.
232
233      IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
234      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
235
236      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
237         DO jj = 2, jpjm1                                    ! sst gradient
238            DO ji = fs_2, fs_jpim1   ! vector opt.
239               zztmp  = tsn(ji,jj,1,jp_tem)
240               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)
241               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)
242               z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
243                  &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
244            END DO
245         END DO
246         CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
247         CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient
248         z2d(:,:) = SQRT( z2d(:,:) )
249         CALL iom_put( "sstgrad" ,  z2d )          ! module of sst gradient
250      ENDIF
251         
252      ! heat and salt contents
253      IF( iom_use("heatc") ) THEN
254         z2d(:,:)  = 0._wp 
255         DO jk = 1, jpkm1
256            DO jj = 1, jpj
257               DO ji = 1, jpi
258                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
259               END DO
260            END DO
261         END DO
262         CALL iom_put( "heatc", rau0_rcp * z2d )   ! vertically integrated heat content (J/m2)
263      ENDIF
264
265      IF( iom_use("saltc") ) THEN
266         z2d(:,:)  = 0._wp 
267         DO jk = 1, jpkm1
268            DO jj = 1, jpj
269               DO ji = 1, jpi
270                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
271               END DO
272            END DO
273         END DO
274         CALL iom_put( "saltc", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
275      ENDIF
276      !
277      IF ( iom_use("eken") ) THEN
278         z3d(:,:,jpk) = 0._wp 
279         DO jk = 1, jpkm1
280            DO jj = 2, jpjm1
281               DO ji = fs_2, fs_jpim1   ! vector opt.
282                  zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
283                  z3d(ji,jj,jk) = zztmp * (  un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)   &
284                     &                     + un(ji  ,jj,jk)**2 * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)   &
285                     &                     + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)   &
286                     &                     + vn(ji,jj  ,jk)**2 * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk)   )
287               END DO
288            END DO
289         END DO
290         CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
291         CALL iom_put( "eken", z3d )                 ! kinetic energy
292      ENDIF
293      !
294      CALL iom_put( "hdiv", hdivn )                  ! Horizontal divergence
295      !
296      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
297         z3d(:,:,jpk) = 0.e0
298         z2d(:,:) = 0.e0
299         DO jk = 1, jpkm1
300            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)
301            z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
302         END DO
303         CALL iom_put( "u_masstr"     , z3d )         ! mass transport in i-direction
304         CALL iom_put( "u_masstr_vint", z2d )         ! mass transport in i-direction vertical sum
305      ENDIF
306     
307      IF( iom_use("u_heattr") ) THEN
308         z2d(:,:) = 0._wp 
309         DO jk = 1, jpkm1
310            DO jj = 2, jpjm1
311               DO ji = fs_2, fs_jpim1   ! vector opt.
312                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
313               END DO
314            END DO
315         END DO
316         CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
317         CALL iom_put( "u_heattr", 0.5*rcp * z2d )    ! heat transport in i-direction
318      ENDIF
319
320      IF( iom_use("u_salttr") ) THEN
321         z2d(:,:) = 0.e0 
322         DO jk = 1, jpkm1
323            DO jj = 2, jpjm1
324               DO ji = fs_2, fs_jpim1   ! vector opt.
325                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
326               END DO
327            END DO
328         END DO
329         CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
330         CALL iom_put( "u_salttr", 0.5 * z2d )        ! heat transport in i-direction
331      ENDIF
332
333     
334      IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
335         z3d(:,:,jpk) = 0.e0
336         DO jk = 1, jpkm1
337            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)
338         END DO
339         CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction
340      ENDIF
341     
342      IF( iom_use("v_heattr") ) THEN
343         z2d(:,:) = 0.e0 
344         DO jk = 1, jpkm1
345            DO jj = 2, jpjm1
346               DO ji = fs_2, fs_jpim1   ! vector opt.
347                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
348               END DO
349            END DO
350         END DO
351         CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
352         CALL iom_put( "v_heattr", 0.5*rcp * z2d )    !  heat transport in j-direction
353      ENDIF
354
355      IF( iom_use("v_salttr") ) THEN
356         z2d(:,:) = 0._wp 
357         DO jk = 1, jpkm1
358            DO jj = 2, jpjm1
359               DO ji = fs_2, fs_jpim1   ! vector opt.
360                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
361               END DO
362            END DO
363         END DO
364         CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
365         CALL iom_put( "v_salttr", 0.5 * z2d )        !  heat transport in j-direction
366      ENDIF
367
368      IF( iom_use("tosmint") ) THEN
369         z2d(:,:) = 0._wp
370         DO jk = 1, jpkm1
371            DO jj = 2, jpjm1
372               DO ji = fs_2, fs_jpim1   ! vector opt.
373                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem)
374               END DO
375            END DO
376         END DO
377         CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
378         CALL iom_put( "tosmint", rau0 * z2d )        ! Vertical integral of temperature
379      ENDIF
380      IF( iom_use("somint") ) THEN
381         z2d(:,:)=0._wp
382         DO jk = 1, jpkm1
383            DO jj = 2, jpjm1
384               DO ji = fs_2, fs_jpim1   ! vector opt.
385                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
386               END DO
387            END DO
388         END DO
389         CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
390         CALL iom_put( "somint", rau0 * z2d )         ! Vertical integral of salinity
391      ENDIF
392
393      CALL iom_put( "bn2", rn2 )                      ! Brunt-Vaisala buoyancy frequency (N^2)
394      !
395
396      IF (ln_diatmb)   CALL dia_tmb                   ! tmb values
397         
398      IF (ln_dia25h)   CALL dia_25h( kt )             ! 25h averaging
399
400      IF( ln_timing )   CALL timing_stop('dia_wri')
401      !
402   END SUBROUTINE dia_wri
403
404#else
405   !!----------------------------------------------------------------------
406   !!   Default option                                  use IOIPSL  library
407   !!----------------------------------------------------------------------
408
409   INTEGER FUNCTION dia_wri_alloc()
410      !!----------------------------------------------------------------------
411      INTEGER, DIMENSION(2) :: ierr
412      !!----------------------------------------------------------------------
413      ierr = 0
414      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
415         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
416         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
417         !
418      dia_wri_alloc = MAXVAL(ierr)
419      CALL mpp_sum( 'diawri', dia_wri_alloc )
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      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zisfdebug
883      !!
884      INTEGER :: inum, jk
885      !!----------------------------------------------------------------------
886      !
887      IF(lwp) WRITE(numout,*)
888      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
889      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
890      IF(lwp) WRITE(numout,*) '                and named :', cdfile_name, '...nc'
891
892#if defined key_si3
893     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
894#else
895     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
896#endif
897
898      CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) )    ! now temperature
899      CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) )    ! now salinity
900      CALL iom_rstput( 0, 0, inum, 'sossheig', sshn              )    ! sea surface height
901      CALL iom_rstput( 0, 0, inum, 'vozocrtx', un                )    ! now i-velocity
902      CALL iom_rstput( 0, 0, inum, 'vomecrty', vn                )    ! now j-velocity
903      IF( ln_zad_Aimp ) THEN
904         CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn + wi        )    ! now k-velocity
905      ELSE
906         CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn             )    ! now k-velocity
907      ENDIF
908      CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep            )    ! now k-velocity
909      CALL iom_rstput( 0, 0, inum, 'ht_n'   , ht_n               )    ! now k-velocity
910      IF ( ln_isf ) THEN
911         IF (ln_isfcav_mlt) THEN
912            CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav          )    ! now k-velocity
913            CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav    )    ! now k-velocity
914            CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav    )    ! now k-velocity
915            CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,8)    )    ! now k-velocity
916            CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,8)    )    ! now k-velocity
917         END IF
918         IF (ln_isfpar_mlt) THEN
919            CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,8)  )    ! now k-velocity
920            CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par          )    ! now k-velocity
921            CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par    )    ! now k-velocity
922            CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par    )    ! now k-velocity
923            CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,8)    )    ! now k-velocity
924            CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,8)    )    ! now k-velocity
925         END IF
926      END IF
927
928      IF ( ln_isf ) THEN
929         IF (ln_isfcav_mlt) CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,8), ktype = jp_i1 )
930         IF (ln_isfpar_mlt) CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,8), ktype = jp_i1 )
931      END IF
932     
933      IF( ALLOCATED(ahtu) ) THEN
934         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point
935         CALL iom_rstput( 0, 0, inum,  'ahtv', ahtv              )    ! aht at v-point
936      ENDIF
937      IF( ALLOCATED(ahmt) ) THEN
938         CALL iom_rstput( 0, 0, inum,  'ahmt', ahmt              )    ! ahmt at u-point
939         CALL iom_rstput( 0, 0, inum,  'ahmf', ahmf              )    ! ahmf at v-point
940      ENDIF
941      CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf         )    ! freshwater budget
942      CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns         )    ! total heat flux
943      CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr               )    ! solar heat flux
944      CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i              )    ! ice fraction
945      CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
946      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
947      IF(  .NOT.ln_linssh  ) THEN             
948         CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept_n        )    !  T-cell depth
949         CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n          )    !  T-cell thickness 
950      END IF
951      IF( ln_wave .AND. ln_sdw ) THEN
952         CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd            )    ! now StokesDrift i-velocity
953         CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd            )    ! now StokesDrift j-velocity
954         CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity
955      ENDIF
956 
957#if defined key_si3
958      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid
959         CALL ice_wri_state( inum )
960      ENDIF
961#endif
962      !
963      CALL iom_close( inum )
964      !
965   END SUBROUTINE dia_wri_state
966
967   !!======================================================================
968END MODULE diawri
Note: See TracBrowser for help on using the repository browser.