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/UKMO/NEMO_4.0_mirror_text_diagnostics/src/OCE/DIA – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_mirror_text_diagnostics/src/OCE/DIA/diawri.F90

Last change on this file was 10986, checked in by andmirek, 5 years ago

GMED 462 add flush
File size: 49.9 KB
Line 
1MODULE diawri
2   !!======================================================================
3   !!                     ***  MODULE  diawri  ***
4   !! Ocean diagnostics :  write ocean output files
5   !!=====================================================================
6   !! History :  OPA  ! 1991-03  (M.-A. Foujols)  Original code
7   !!            4.0  ! 1991-11  (G. Madec)
8   !!                 ! 1992-06  (M. Imbard)  correction restart file
9   !!                 ! 1992-07  (M. Imbard)  split into diawri and rstwri
10   !!                 ! 1993-03  (M. Imbard)  suppress writibm
11   !!                 ! 1998-01  (C. Levy)  NETCDF format using ioipsl INTERFACE
12   !!                 ! 1999-02  (E. Guilyardi)  name of netCDF files + variables
13   !!            8.2  ! 2000-06  (M. Imbard)  Original code (diabort.F)
14   !!   NEMO     1.0  ! 2002-06  (A.Bozec, E. Durand)  Original code (diainit.F)
15   !!             -   ! 2002-09  (G. Madec)  F90: Free form and module
16   !!             -   ! 2002-12  (G. Madec)  merge of diabort and diainit, F90
17   !!                 ! 2005-11  (V. Garnier) Surface pressure gradient organization
18   !!            3.2  ! 2008-11  (B. Lemaire) creation from old diawri
19   !!            3.7  ! 2014-01  (G. Madec) remove eddy induced velocity from no-IOM output
20   !!                 !                     change name of output variables in dia_wri_state
21   !!----------------------------------------------------------------------
22
23   !!----------------------------------------------------------------------
24   !!   dia_wri       : create the standart output files
25   !!   dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
26   !!----------------------------------------------------------------------
27   USE oce            ! ocean dynamics and tracers
28   USE 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 diatmb         ! Top,middle,bottom output
49   USE dia25h         ! 25h Mean output
50   USE iom            !
51   USE ioipsl         !
52
53#if defined key_si3
54   USE ice 
55   USE icewri 
56#endif
57   USE lib_mpp         ! MPP library
58   USE timing          ! preformance summary
59   USE diurnal_bulk    ! diurnal warm layer
60   USE cool_skin       ! Cool skin
61
62   IMPLICIT NONE
63   PRIVATE
64
65   PUBLIC   dia_wri                 ! routines called by step.F90
66   PUBLIC   dia_wri_state
67   PUBLIC   dia_wri_alloc           ! Called by nemogcm module
68
69   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
70   INTEGER ::          nb_T              , ndim_bT   ! grid_T file
71   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
72   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
73   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file
74   INTEGER ::   ndex(1)                              ! ???
75   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
76   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
77   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
78
79   !! * Substitutions
80#  include "vectopt_loop_substitute.h90"
81   !!----------------------------------------------------------------------
82   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
83   !! $Id$
84   !! Software governed by the CeCILL license (see ./LICENSE)
85   !!----------------------------------------------------------------------
86CONTAINS
87
88#if defined key_iomput
89   !!----------------------------------------------------------------------
90   !!   'key_iomput'                                        use IOM library
91   !!----------------------------------------------------------------------
92   INTEGER FUNCTION dia_wri_alloc()
93      !
94      dia_wri_alloc = 0
95      !
96   END FUNCTION dia_wri_alloc
97
98   
99   SUBROUTINE dia_wri( kt )
100      !!---------------------------------------------------------------------
101      !!                  ***  ROUTINE dia_wri  ***
102      !!                   
103      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
104      !!      NETCDF format is used by default
105      !!
106      !! ** Method  :  use iom_put
107      !!----------------------------------------------------------------------
108      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
109      !!
110      INTEGER ::   ji, jj, jk       ! dummy loop indices
111      INTEGER ::   ikbot            ! local integer
112      REAL(wp)::   zztmp , zztmpx   ! local scalar
113      REAL(wp)::   zztmp2, zztmpy   !   -      -
114      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace
115      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace
116      !!----------------------------------------------------------------------
117      !
118      IF( ln_timing )   CALL timing_start('dia_wri')
119      !
120      ! Output the initial state and forcings
121      IF( ninist == 1 ) THEN                       
122         CALL dia_wri_state( 'output.init' )
123         ninist = 0
124      ENDIF
125
126      ! Output of initial vertical scale factor
127      CALL iom_put("e3t_0", e3t_0(:,:,:) )
128      CALL iom_put("e3u_0", e3u_0(:,:,:) )
129      CALL iom_put("e3v_0", e3v_0(:,:,:) )
130      !
131      CALL iom_put( "e3t" , e3t_n(:,:,:) )
132      CALL iom_put( "e3u" , e3u_n(:,:,:) )
133      CALL iom_put( "e3v" , e3v_n(:,:,:) )
134      CALL iom_put( "e3w" , e3w_n(:,:,:) )
135      IF( iom_use("e3tdef") )   &
136         CALL iom_put( "e3tdef"  , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
137
138      IF( ll_wd ) THEN
139         CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying)
140      ELSE
141         CALL iom_put( "ssh" , sshn )              ! sea surface height
142      ENDIF
143
144      IF( iom_use("wetdep") )   &                  ! wet depth
145         CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) )
146     
147      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
148      CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature
149      IF ( iom_use("sbt") ) THEN
150         DO jj = 1, jpj
151            DO ji = 1, jpi
152               ikbot = mbkt(ji,jj)
153               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem)
154            END DO
155         END DO
156         CALL iom_put( "sbt", z2d )                ! bottom temperature
157      ENDIF
158     
159      CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity
160      CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity
161      IF ( iom_use("sbs") ) THEN
162         DO jj = 1, jpj
163            DO ji = 1, jpi
164               ikbot = mbkt(ji,jj)
165               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal)
166            END DO
167         END DO
168         CALL iom_put( "sbs", z2d )                ! bottom salinity
169      ENDIF
170
171      IF ( iom_use("taubot") ) THEN                ! bottom stress
172         zztmp = rau0 * 0.25
173         z2d(:,:) = 0._wp
174         DO jj = 2, jpjm1
175            DO ji = fs_2, fs_jpim1   ! vector opt.
176               zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * un(ji  ,jj,mbku(ji  ,jj))  )**2   &
177                  &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj))  )**2   &
178                  &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vn(ji,jj  ,mbkv(ji,jj  ))  )**2   &
179                  &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1))  )**2
180               z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 
181               !
182            END DO
183         END DO
184         CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
185         CALL iom_put( "taubot", z2d )           
186      ENDIF
187         
188      CALL iom_put( "uoce", un(:,:,:) )            ! 3D i-current
189      CALL iom_put(  "ssu", un(:,:,1) )            ! surface i-current
190      IF ( iom_use("sbu") ) THEN
191         DO jj = 1, jpj
192            DO ji = 1, jpi
193               ikbot = mbku(ji,jj)
194               z2d(ji,jj) = un(ji,jj,ikbot)
195            END DO
196         END DO
197         CALL iom_put( "sbu", z2d )                ! bottom i-current
198      ENDIF
199     
200      CALL iom_put( "voce", vn(:,:,:) )            ! 3D j-current
201      CALL iom_put(  "ssv", vn(:,:,1) )            ! surface j-current
202      IF ( iom_use("sbv") ) THEN
203         DO jj = 1, jpj
204            DO ji = 1, jpi
205               ikbot = mbkv(ji,jj)
206               z2d(ji,jj) = vn(ji,jj,ikbot)
207            END DO
208         END DO
209         CALL iom_put( "sbv", z2d )                ! bottom j-current
210      ENDIF
211
212      CALL iom_put( "woce", wn )                   ! vertical velocity
213      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
214         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
215         z2d(:,:) = rau0 * e1e2t(:,:)
216         DO jk = 1, jpk
217            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
218         END DO
219         CALL iom_put( "w_masstr" , z3d ) 
220         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
221      ENDIF
222
223      CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef.
224      CALL iom_put( "avs" , avs )                  ! S vert. eddy diff. coef.
225      CALL iom_put( "avm" , avm )                  ! T vert. eddy visc. coef.
226
227      IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
228      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
229
230      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
231         DO jj = 2, jpjm1                                    ! sst gradient
232            DO ji = fs_2, fs_jpim1   ! vector opt.
233               zztmp  = tsn(ji,jj,1,jp_tem)
234               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)
235               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)
236               z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
237                  &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
238            END DO
239         END DO
240         CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
241         CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient
242         z2d(:,:) = SQRT( z2d(:,:) )
243         CALL iom_put( "sstgrad" ,  z2d )          ! module of sst gradient
244      ENDIF
245         
246      ! heat and salt contents
247      IF( iom_use("heatc") ) THEN
248         z2d(:,:)  = 0._wp 
249         DO jk = 1, jpkm1
250            DO jj = 1, jpj
251               DO ji = 1, jpi
252                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
253               END DO
254            END DO
255         END DO
256         CALL iom_put( "heatc", rau0_rcp * z2d )   ! vertically integrated heat content (J/m2)
257      ENDIF
258
259      IF( iom_use("saltc") ) THEN
260         z2d(:,:)  = 0._wp 
261         DO jk = 1, jpkm1
262            DO jj = 1, jpj
263               DO ji = 1, jpi
264                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
265               END DO
266            END DO
267         END DO
268         CALL iom_put( "saltc", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
269      ENDIF
270      !
271      IF ( iom_use("eken") ) THEN
272         z3d(:,:,jpk) = 0._wp 
273         DO jk = 1, jpkm1
274            DO jj = 2, jpjm1
275               DO ji = fs_2, fs_jpim1   ! vector opt.
276                  zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
277                  z3d(ji,jj,jk) = zztmp * (  un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)   &
278                     &                     + un(ji  ,jj,jk)**2 * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)   &
279                     &                     + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)   &
280                     &                     + vn(ji,jj  ,jk)**2 * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk)   )
281               END DO
282            END DO
283         END DO
284         CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
285         CALL iom_put( "eken", z3d )                 ! kinetic energy
286      ENDIF
287      !
288      CALL iom_put( "hdiv", hdivn )                  ! Horizontal divergence
289      !
290      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
291         z3d(:,:,jpk) = 0.e0
292         z2d(:,:) = 0.e0
293         DO jk = 1, jpkm1
294            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)
295            z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
296         END DO
297         CALL iom_put( "u_masstr"     , z3d )         ! mass transport in i-direction
298         CALL iom_put( "u_masstr_vint", z2d )         ! mass transport in i-direction vertical sum
299      ENDIF
300     
301      IF( iom_use("u_heattr") ) THEN
302         z2d(:,:) = 0._wp 
303         DO jk = 1, jpkm1
304            DO jj = 2, jpjm1
305               DO ji = fs_2, fs_jpim1   ! vector opt.
306                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
307               END DO
308            END DO
309         END DO
310         CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
311         CALL iom_put( "u_heattr", 0.5*rcp * z2d )    ! heat transport in i-direction
312      ENDIF
313
314      IF( iom_use("u_salttr") ) THEN
315         z2d(:,:) = 0.e0 
316         DO jk = 1, jpkm1
317            DO jj = 2, jpjm1
318               DO ji = fs_2, fs_jpim1   ! vector opt.
319                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
320               END DO
321            END DO
322         END DO
323         CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
324         CALL iom_put( "u_salttr", 0.5 * z2d )        ! heat transport in i-direction
325      ENDIF
326
327     
328      IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
329         z3d(:,:,jpk) = 0.e0
330         DO jk = 1, jpkm1
331            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)
332         END DO
333         CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction
334      ENDIF
335     
336      IF( iom_use("v_heattr") ) THEN
337         z2d(:,:) = 0.e0 
338         DO jk = 1, jpkm1
339            DO jj = 2, jpjm1
340               DO ji = fs_2, fs_jpim1   ! vector opt.
341                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
342               END DO
343            END DO
344         END DO
345         CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
346         CALL iom_put( "v_heattr", 0.5*rcp * z2d )    !  heat transport in j-direction
347      ENDIF
348
349      IF( iom_use("v_salttr") ) THEN
350         z2d(:,:) = 0._wp 
351         DO jk = 1, jpkm1
352            DO jj = 2, jpjm1
353               DO ji = fs_2, fs_jpim1   ! vector opt.
354                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
355               END DO
356            END DO
357         END DO
358         CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
359         CALL iom_put( "v_salttr", 0.5 * z2d )        !  heat transport in j-direction
360      ENDIF
361
362      IF( iom_use("tosmint") ) THEN
363         z2d(:,:) = 0._wp
364         DO jk = 1, jpkm1
365            DO jj = 2, jpjm1
366               DO ji = fs_2, fs_jpim1   ! vector opt.
367                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem)
368               END DO
369            END DO
370         END DO
371         CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
372         CALL iom_put( "tosmint", rau0 * z2d )        ! Vertical integral of temperature
373      ENDIF
374      IF( iom_use("somint") ) THEN
375         z2d(:,:)=0._wp
376         DO jk = 1, jpkm1
377            DO jj = 2, jpjm1
378               DO ji = fs_2, fs_jpim1   ! vector opt.
379                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
380               END DO
381            END DO
382         END DO
383         CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
384         CALL iom_put( "somint", rau0 * z2d )         ! Vertical integral of salinity
385      ENDIF
386
387      CALL iom_put( "bn2", rn2 )                      ! Brunt-Vaisala buoyancy frequency (N^2)
388      !
389
390      IF (ln_diatmb)   CALL dia_tmb                   ! tmb values
391         
392      IF (ln_dia25h)   CALL dia_25h( kt )             ! 25h averaging
393
394      IF( ln_timing )   CALL timing_stop('dia_wri')
395      !
396   END SUBROUTINE dia_wri
397
398#else
399   !!----------------------------------------------------------------------
400   !!   Default option                                  use IOIPSL  library
401   !!----------------------------------------------------------------------
402
403   INTEGER FUNCTION dia_wri_alloc()
404      !!----------------------------------------------------------------------
405      INTEGER, DIMENSION(2) :: ierr
406      !!----------------------------------------------------------------------
407      ierr = 0
408      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
409         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
410         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
411         !
412      dia_wri_alloc = MAXVAL(ierr)
413      CALL mpp_sum( 'diawri', dia_wri_alloc )
414      !
415   END FUNCTION dia_wri_alloc
416
417   
418   SUBROUTINE dia_wri( kt )
419      !!---------------------------------------------------------------------
420      !!                  ***  ROUTINE dia_wri  ***
421      !!                   
422      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
423      !!      NETCDF format is used by default
424      !!
425      !! ** Method  :   At the beginning of the first time step (nit000),
426      !!      define all the NETCDF files and fields
427      !!      At each time step call histdef to compute the mean if ncessary
428      !!      Each nwrite time step, output the instantaneous or mean fields
429      !!----------------------------------------------------------------------
430      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
431      !
432      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
433      INTEGER  ::   inum = 11                                ! temporary logical unit
434      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
435      INTEGER  ::   ierr                                     ! error code return from allocation
436      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
437      INTEGER  ::   jn, ierror                               ! local integers
438      REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
439      !
440      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace
441      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace
442      !!----------------------------------------------------------------------
443      !
444      IF( ln_timing )   CALL timing_start('dia_wri')
445      !
446      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==!
447         CALL dia_wri_state( 'output.init' )
448         ninist = 0
449      ENDIF
450      !
451      ! 0. Initialisation
452      ! -----------------
453
454      ll_print = .FALSE.                  ! local variable for debugging
455      ll_print = ll_print .AND. lwp
456
457      ! Define frequency of output and means
458      clop = "x"         ! no use of the mask value (require less cpu time and otherwise the model crashes)
459#if defined key_diainstant
460      zsto = nwrite * rdt
461      clop = "inst("//TRIM(clop)//")"
462#else
463      zsto=rdt
464      clop = "ave("//TRIM(clop)//")"
465#endif
466      zout = nwrite * rdt
467      zmax = ( nitend - nit000 + 1 ) * rdt
468
469      ! Define indices of the horizontal output zoom and vertical limit storage
470      iimi = 1      ;      iima = jpi
471      ijmi = 1      ;      ijma = jpj
472      ipk = jpk
473
474      ! define time axis
475      it = kt
476      itmod = kt - nit000 + 1
477
478
479      ! 1. Define NETCDF files and fields at beginning of first time step
480      ! -----------------------------------------------------------------
481
482      IF( kt == nit000 ) THEN
483
484         ! Define the NETCDF files (one per grid)
485
486         ! Compute julian date from starting date of the run
487         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
488         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
489         IF(lwp) THEN
490            WRITE(numout,*)
491            WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
492            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
493            WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
494                                 ' limit storage in depth = ', ipk
495            IF(lflush) CALL FLUSH(numout)
496            ! WRITE root name in date.file for use by postpro
497            CALL dia_nam( clhstnam, nwrite,' ' )
498            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
499            WRITE(inum,*) clhstnam
500            CLOSE(inum)
501         ENDIF
502
503         ! Define the T grid FILE ( nid_T )
504
505         CALL dia_nam( clhstnam, nwrite, 'grid_T' )
506         IF(lwp) THEN
507            WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
508            IF(lflush) CALL FLUSH(numout)
509         ENDIF
510         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
511            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
512            &          nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
513         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
514            &           "m", ipk, gdept_1d, nz_T, "down" )
515         !                                                            ! Index of ocean points
516         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
517         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
518         !
519         IF( ln_icebergs ) THEN
520            !
521            !! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
522            !! that routine is called from nemogcm, so do it here immediately before its needed
523            ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
524            CALL mpp_sum( 'diawri', ierror )
525            IF( ierror /= 0 ) THEN
526               CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
527               RETURN
528            ENDIF
529            !
530            !! iceberg vertical coordinate is class number
531            CALL histvert( nid_T, "class", "Iceberg class",      &  ! Vertical grid: class
532               &           "number", nclasses, class_num, nb_T )
533            !
534            !! each class just needs the surface index pattern
535            ndim_bT = 3
536            DO jn = 1,nclasses
537               ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
538            ENDDO
539            !
540         ENDIF
541
542         ! Define the U grid FILE ( nid_U )
543
544         CALL dia_nam( clhstnam, nwrite, 'grid_U' )
545         IF(lwp) THEN
546            WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
547            IF(lflush) CALL FLUSH(numout)
548         ENDIF
549         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
550            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
551            &          nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
552         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
553            &           "m", ipk, gdept_1d, nz_U, "down" )
554         !                                                            ! Index of ocean points
555         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
556         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
557
558         ! Define the V grid FILE ( nid_V )
559
560         CALL dia_nam( clhstnam, nwrite, 'grid_V' )                   ! filename
561         IF(lwp) THEN
562            WRITE(numout,*) " Name of NETCDF file ", clhstnam
563            IF(lflush) CALL FLUSH(numout)
564         ENDIF
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, nwrite, 'grid_W' )                   ! filename
577         IF(lwp) THEN
578            WRITE(numout,*) " Name of NETCDF file ", clhstnam
579            IF(lflush) CALL FLUSH(numout)
580         ENDIF
581         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
582            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
583            &          nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
584         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
585            &          "m", ipk, gdepw_1d, nz_W, "down" )
586
587
588         ! Declare all the output fields as NETCDF variables
589
590         !                                                                                      !!! nid_T : 3D
591         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
592            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
593         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
594            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
595         IF(  .NOT.ln_linssh  ) THEN
596            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n
597            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
598            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t_n
599            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
600            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t_n
601            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
602         ENDIF
603         !                                                                                      !!! nid_T : 2D
604         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
605            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
606         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
607            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
608         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
609            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
610         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
611            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
612         CALL histdef( nid_T, "sorunoff", "River runoffs"                      , "Kg/m2/s",   &  ! runoffs
613            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
614         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
615            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
616         IF(  ln_linssh  ) THEN
617            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem)
618            &                                                                  , "KgC/m2/s",  &  ! sosst_cd
619            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
620            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * tsn(:,:,1,jp_sal)
621            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
622            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
623         ENDIF
624         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
625            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
626         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
627            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
628         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
629            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
630         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
631            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
632         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
633            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
634         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
635            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
636!
637         IF( ln_icebergs ) THEN
638            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
639               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
640            CALL histdef( nid_T, "calving_heat"        , "calving heat flux"                        , "XXXX"   , &
641               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
642            CALL histdef( nid_T, "berg_floating_melt"  , "Melt rate of icebergs + bits"             , "kg/m2/s", &
643               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
644            CALL histdef( nid_T, "berg_stored_ice"     , "Accumulated ice mass by class"            , "kg"     , &
645               &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
646            IF( ln_bergdia ) THEN
647               CALL histdef( nid_T, "berg_melt"           , "Melt rate of icebergs"                    , "kg/m2/s", &
648                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
649               CALL histdef( nid_T, "berg_buoy_melt"      , "Buoyancy component of iceberg melt rate"  , "kg/m2/s", &
650                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
651               CALL histdef( nid_T, "berg_eros_melt"      , "Erosion component of iceberg melt rate"   , "kg/m2/s", &
652                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
653               CALL histdef( nid_T, "berg_conv_melt"      , "Convective component of iceberg melt rate", "kg/m2/s", &
654                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
655               CALL histdef( nid_T, "berg_virtual_area"   , "Virtual coverage by icebergs"             , "m2"     , &
656                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
657               CALL histdef( nid_T, "bits_src"           , "Mass source of bergy bits"                , "kg/m2/s", &
658                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
659               CALL histdef( nid_T, "bits_melt"          , "Melt rate of bergy bits"                  , "kg/m2/s", &
660                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
661               CALL histdef( nid_T, "bits_mass"          , "Bergy bit density field"                  , "kg/m2"  , &
662                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
663               CALL histdef( nid_T, "berg_mass"           , "Iceberg density field"                    , "kg/m2"  , &
664                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
665               CALL histdef( nid_T, "berg_real_calving"   , "Calving into iceberg class"               , "kg/s"   , &
666                  &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
667            ENDIF
668         ENDIF
669
670         IF( .NOT. ln_cpl ) THEN
671            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
672               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
673            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
674               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
675            CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
676               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
677         ENDIF
678
679         IF( ln_cpl .AND. nn_ice <= 1 ) THEN
680            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
681               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
682            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
683               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
684            CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn
685               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
686         ENDIF
687         
688         clmx ="l_max(only(x))"    ! max index on a period
689!         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
690!            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
691#if defined key_diahth
692         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
693            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
694         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
695            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
696         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
697            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
698         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "J/m2"   ,   & ! htc3
699            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
700#endif
701
702         CALL histend( nid_T, snc4chunks=snc4set )
703
704         !                                                                                      !!! nid_U : 3D
705         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
706            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
707         IF( ln_wave .AND. ln_sdw) THEN
708            CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current"         , "m/s"    ,   &  ! usd
709               &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
710         ENDIF
711         !                                                                                      !!! nid_U : 2D
712         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
713            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
714
715         CALL histend( nid_U, snc4chunks=snc4set )
716
717         !                                                                                      !!! nid_V : 3D
718         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
719            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
720         IF( ln_wave .AND. ln_sdw) THEN
721            CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current"    , "m/s"    ,   &  ! vsd
722               &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
723         ENDIF
724         !                                                                                      !!! nid_V : 2D
725         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
726            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
727
728         CALL histend( nid_V, snc4chunks=snc4set )
729
730         !                                                                                      !!! nid_W : 3D
731         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn
732            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
733         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
734            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
735         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avm
736            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
737
738         IF( ln_zdfddm ) THEN
739            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
740               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
741         ENDIF
742         
743         IF( ln_wave .AND. ln_sdw) THEN
744            CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current"   , "m/s"    ,   &  ! wsd
745               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
746         ENDIF
747         !                                                                                      !!! nid_W : 2D
748         CALL histend( nid_W, snc4chunks=snc4set )
749
750         IF(lwp) THEN
751            WRITE(numout,*)
752            WRITE(numout,*) 'End of NetCDF Initialization'
753            IF(lflush) CALL FLUSH(numout)
754         ENDIF
755
756      ENDIF
757
758      ! 2. Start writing data
759      ! ---------------------
760
761      ! ndex(1) est utilise ssi l'avant dernier argument est different de
762      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
763      ! donne le nombre d'elements, et ndex la liste des indices a sortir
764
765      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
766         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
767         WRITE(numout,*) '~~~~~~ '
768         IF(lflush) CALL FLUSH(numout)
769      ENDIF
770
771      IF( .NOT.ln_linssh ) THEN
772         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
773         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
774         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
775         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
776      ELSE
777         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature
778         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity
779         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature
780         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
781      ENDIF
782      IF( .NOT.ln_linssh ) THEN
783         zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
784         CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
785         CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
786         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
787      ENDIF
788      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
789      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
790      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
791      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
792                                                                                  ! (includes virtual salt flux beneath ice
793                                                                                  ! in linear free surface case)
794      IF( ln_linssh ) THEN
795         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
796         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
797         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
798         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
799      ENDIF
800      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
801      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
802      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
803      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
804      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
805      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
806!
807      IF( ln_icebergs ) THEN
808         !
809         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT ) 
810         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
811         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT ) 
812         !
813         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
814         !
815         IF( ln_bergdia ) THEN
816            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   ) 
817            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   ) 
818            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   ) 
819            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   ) 
820            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   ) 
821            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   ) 
822            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   ) 
823            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   ) 
824            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   ) 
825            !
826            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
827         ENDIF
828      ENDIF
829
830      IF( .NOT. ln_cpl ) THEN
831         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
832         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
833         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
834         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
835      ENDIF
836      IF( ln_cpl .AND. nn_ice <= 1 ) THEN
837         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
838         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
839         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
840         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
841      ENDIF
842!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
843!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
844
845#if defined key_diahth
846      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
847      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
848      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
849      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
850#endif
851
852      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
853      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
854
855      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
856      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
857
858      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
859      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
860      CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
861      IF( ln_zdfddm ) THEN
862         CALL histwrite( nid_W, "voddmavs", it, avs         , ndim_T, ndex_T )    ! S vert. eddy diff. coef.
863      ENDIF
864
865      IF( ln_wave .AND. ln_sdw ) THEN
866         CALL histwrite( nid_U, "sdzocrtx", it, usd         , ndim_U , ndex_U )    ! i-StokesDrift-current
867         CALL histwrite( nid_V, "sdmecrty", it, vsd         , ndim_V , ndex_V )    ! j-StokesDrift-current
868         CALL histwrite( nid_W, "sdvecrtz", it, wsd         , ndim_T , ndex_T )    ! StokesDrift vert. current
869      ENDIF
870
871      ! 3. Close all files
872      ! ---------------------------------------
873      IF( kt == nitend ) THEN
874         CALL histclo( nid_T )
875         CALL histclo( nid_U )
876         CALL histclo( nid_V )
877         CALL histclo( nid_W )
878      ENDIF
879      !
880      IF( ln_timing )   CALL timing_stop('dia_wri')
881      !
882   END SUBROUTINE dia_wri
883#endif
884
885   SUBROUTINE dia_wri_state( cdfile_name )
886      !!---------------------------------------------------------------------
887      !!                 ***  ROUTINE dia_wri_state  ***
888      !!       
889      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
890      !!      the instantaneous ocean state and forcing fields.
891      !!        Used to find errors in the initial state or save the last
892      !!      ocean state in case of abnormal end of a simulation
893      !!
894      !! ** Method  :   NetCDF files using ioipsl
895      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
896      !!      File 'output.abort.nc' is created in case of abnormal job end
897      !!----------------------------------------------------------------------
898      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
899      !!
900      INTEGER :: inum
901      !!----------------------------------------------------------------------
902      !
903      IF(lwp) THEN
904         WRITE(numout,*)
905         WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
906         WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
907         WRITE(numout,*) '                and named :', cdfile_name, '...nc'
908         IF(lflush) CALL FLUSH(numout)
909      ENDIF
910
911#if defined key_si3
912     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
913#else
914     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
915#endif
916
917      CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) )    ! now temperature
918      CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) )    ! now salinity
919      CALL iom_rstput( 0, 0, inum, 'sossheig', sshn              )    ! sea surface height
920      CALL iom_rstput( 0, 0, inum, 'vozocrtx', un                )    ! now i-velocity
921      CALL iom_rstput( 0, 0, inum, 'vomecrty', vn                )    ! now j-velocity
922      CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn                )    ! now k-velocity
923      IF( ALLOCATED(ahtu) ) THEN
924         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point
925         CALL iom_rstput( 0, 0, inum,  'ahtv', ahtv              )    ! aht at v-point
926      ENDIF
927      IF( ALLOCATED(ahmt) ) THEN
928         CALL iom_rstput( 0, 0, inum,  'ahmt', ahmt              )    ! ahmt at u-point
929         CALL iom_rstput( 0, 0, inum,  'ahmf', ahmf              )    ! ahmf at v-point
930      ENDIF
931      CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf         )    ! freshwater budget
932      CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns         )    ! total heat flux
933      CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr               )    ! solar heat flux
934      CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i              )    ! ice fraction
935      CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
936      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
937      IF(  .NOT.ln_linssh  ) THEN             
938         CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept_n        )    !  T-cell depth
939         CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n          )    !  T-cell thickness 
940      END IF
941      IF( ln_wave .AND. ln_sdw ) THEN
942         CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd            )    ! now StokesDrift i-velocity
943         CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd            )    ! now StokesDrift j-velocity
944         CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity
945      ENDIF
946 
947#if defined key_si3
948      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid
949         CALL ice_wri_state( inum )
950      ENDIF
951#endif
952      !
953      CALL iom_close( inum )
954      !
955   END SUBROUTINE dia_wri_state
956
957   !!======================================================================
958END MODULE diawri
Note: See TracBrowser for help on using the repository browser.