source: trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 7753

Last change on this file since 7753 was 7753, checked in by mocavero, 4 years ago

Reverting trunk to remove OpenMP

  • Property svn:keywords set to Id
File size: 57.0 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 dynadv, ONLY: ln_dynadv_vec
30   USE zdf_oce         ! ocean vertical physics
31   USE ldftra          ! lateral physics: eddy diffusivity coef.
32   USE ldfdyn          ! lateral physics: eddy viscosity   coef.
33   USE sbc_oce         ! Surface boundary condition: ocean fields
34   USE sbc_ice         ! Surface boundary condition: ice fields
35   USE icb_oce         ! Icebergs
36   USE icbdia          ! Iceberg budgets
37   USE sbcssr          ! restoring term toward SST/SSS climatology
38   USE phycst          ! physical constants
39   USE zdfmxl          ! mixed layer
40   USE dianam          ! build name of file (routine)
41   USE zdfddm          ! vertical  physics: double diffusion
42   USE diahth          ! thermocline diagnostics
43   USE wet_dry         ! wetting and drying
44   USE sbcwave         ! wave parameters
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_lim2
54   USE limwri_2 
55#elif defined key_lim3
56   USE limwri 
57#endif
58   USE lib_mpp         ! MPP library
59   USE timing          ! preformance summary
60   USE diurnal_bulk    ! diurnal warm layer
61   USE cool_skin       ! Cool skin
62   USE wrk_nemo        ! working array
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 "zdfddm_substitute.h90"
83#  include "vectopt_loop_substitute.h90"
84   !!----------------------------------------------------------------------
85   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
86   !! $Id$
87   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
88   !!----------------------------------------------------------------------
89CONTAINS
90
91   INTEGER FUNCTION dia_wri_alloc()
92      !!----------------------------------------------------------------------
93      INTEGER, DIMENSION(2) :: ierr
94      !!----------------------------------------------------------------------
95      ierr = 0
96      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
97         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
98         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
99         !
100      dia_wri_alloc = MAXVAL(ierr)
101      IF( lk_mpp )   CALL mpp_sum( dia_wri_alloc )
102      !
103  END FUNCTION dia_wri_alloc
104
105   !!----------------------------------------------------------------------
106   !!   Default option                                   NetCDF output file
107   !!----------------------------------------------------------------------
108#if defined key_iomput
109   !!----------------------------------------------------------------------
110   !!   'key_iomput'                                        use IOM library
111   !!----------------------------------------------------------------------
112
113   SUBROUTINE dia_wri( kt )
114      !!---------------------------------------------------------------------
115      !!                  ***  ROUTINE dia_wri  ***
116      !!                   
117      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
118      !!      NETCDF format is used by default
119      !!
120      !! ** Method  :  use iom_put
121      !!----------------------------------------------------------------------
122      !!
123      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
124      !!
125      INTEGER                      ::   ji, jj, jk              ! dummy loop indices
126      INTEGER                      ::   jkbot                   !
127      REAL(wp)                     ::   zztmp, zztmpx, zztmpy   !
128      !!
129      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace
130      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace
131      !!----------------------------------------------------------------------
132      !
133      IF( nn_timing == 1 )   CALL timing_start('dia_wri')
134      !
135      CALL wrk_alloc( jpi , jpj      , z2d )
136      CALL wrk_alloc( jpi , jpj, jpk , z3d )
137      !
138      ! Output the initial state and forcings
139      IF( ninist == 1 ) THEN                       
140         CALL dia_wri_state( 'output.init', kt )
141         ninist = 0
142      ENDIF
143
144      ! Output of initial vertical scale factor
145      CALL iom_put("e3t_0", e3t_0(:,:,:) )
146      CALL iom_put("e3u_0", e3t_0(:,:,:) )
147      CALL iom_put("e3v_0", e3t_0(:,:,:) )
148      !
149      CALL iom_put( "e3t" , e3t_n(:,:,:) )
150      CALL iom_put( "e3u" , e3u_n(:,:,:) )
151      CALL iom_put( "e3v" , e3v_n(:,:,:) )
152      CALL iom_put( "e3w" , e3w_n(:,:,:) )
153      IF( iom_use("e3tdef") )   &
154         CALL iom_put( "e3tdef"  , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
155
156      CALL iom_put( "ssh" , sshn )                 ! sea surface height
157      IF( iom_use("wetdep") )   &                  ! wet depth
158         CALL iom_put( "wetdep" , ht_wd(:,:) + sshn(:,:) )
159     
160      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
161      CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature
162      IF ( iom_use("sbt") ) THEN
163         DO jj = 1, jpj
164            DO ji = 1, jpi
165               jkbot = mbkt(ji,jj)
166               z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem)
167            END DO
168         END DO
169         CALL iom_put( "sbt", z2d )                ! bottom temperature
170      ENDIF
171     
172      CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity
173      CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity
174      IF ( iom_use("sbs") ) THEN
175         DO jj = 1, jpj
176            DO ji = 1, jpi
177               jkbot = mbkt(ji,jj)
178               z2d(ji,jj) = tsn(ji,jj,jkbot,jp_sal)
179            END DO
180         END DO
181         CALL iom_put( "sbs", z2d )                ! bottom salinity
182      ENDIF
183
184      IF ( iom_use("taubot") ) THEN                ! bottom stress
185         z2d(:,:) = 0._wp
186         DO jj = 2, jpjm1
187            DO ji = fs_2, fs_jpim1   ! vector opt.
188               zztmpx = (  bfrua(ji  ,jj) * un(ji  ,jj,mbku(ji  ,jj))  &
189                      &  + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj))  )     
190               zztmpy = (  bfrva(ji,  jj) * vn(ji,jj  ,mbkv(ji,jj  ))  &
191                      &  + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1))  ) 
192               z2d(ji,jj) = rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1) 
193               !
194            ENDDO
195         ENDDO
196         CALL lbc_lnk( z2d, 'T', 1. )
197         CALL iom_put( "taubot", z2d )           
198      ENDIF
199         
200      CALL iom_put( "uoce", un(:,:,:)         )    ! 3D i-current
201      CALL iom_put(  "ssu", un(:,:,1)         )    ! surface i-current
202      IF ( iom_use("sbu") ) THEN
203         DO jj = 1, jpj
204            DO ji = 1, jpi
205               jkbot = mbku(ji,jj)
206               z2d(ji,jj) = un(ji,jj,jkbot)
207            END DO
208         END DO
209         CALL iom_put( "sbu", z2d )                ! bottom i-current
210      ENDIF
211     
212      CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current
213      CALL iom_put(  "ssv", vn(:,:,1)         )    ! surface j-current
214      IF ( iom_use("sbv") ) THEN
215         DO jj = 1, jpj
216            DO ji = 1, jpi
217               jkbot = mbkv(ji,jj)
218               z2d(ji,jj) = vn(ji,jj,jkbot)
219            END DO
220         END DO
221         CALL iom_put( "sbv", z2d )                ! bottom j-current
222      ENDIF
223
224      CALL iom_put( "woce", wn )                   ! vertical velocity
225      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
226         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
227         z2d(:,:) = rau0 * e1e2t(:,:)
228         DO jk = 1, jpk
229            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
230         END DO
231         CALL iom_put( "w_masstr" , z3d ) 
232         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
233      ENDIF
234
235      CALL iom_put( "avt" , avt                        )    ! T vert. eddy diff. coef.
236      CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef.
237      CALL iom_put( "avs" , fsavs(:,:,:)               )    ! S vert. eddy diff. coef. (useful only with key_zdfddm)
238
239      IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt  (:,:,:) ) ) )
240      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, fsavs(:,:,:) ) ) )
241
242      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
243         DO jj = 2, jpjm1                                    ! sst gradient
244            DO ji = fs_2, fs_jpim1   ! vector opt.
245               zztmp  = tsn(ji,jj,1,jp_tem)
246               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)
247               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)
248               z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
249                  &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
250            END DO
251         END DO
252         CALL lbc_lnk( z2d, 'T', 1. )
253         CALL iom_put( "sstgrad2",  z2d               )    ! square of module of sst gradient
254         z2d(:,:) = SQRT( z2d(:,:) )
255         CALL iom_put( "sstgrad" ,  z2d               )    ! module of sst gradient
256      ENDIF
257         
258      ! clem: heat and salt content
259      IF( iom_use("heatc") ) 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_tem) * tmask(ji,jj,jk)
265               END DO
266            END DO
267         END DO
268         CALL iom_put( "heatc", (rau0 * rcp) * z2d )    ! vertically integrated heat content (J/m2)
269      ENDIF
270
271      IF( iom_use("saltc") ) THEN
272         z2d(:,:)  = 0._wp 
273         DO jk = 1, jpkm1
274            DO jj = 1, jpj
275               DO ji = 1, jpi
276                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
277               END DO
278            END DO
279         END DO
280         CALL iom_put( "saltc", rau0 * z2d )   ! vertically integrated salt content (PSU*kg/m2)
281      ENDIF
282      !
283      IF ( iom_use("eken") ) THEN
284         rke(:,:,jk) = 0._wp                               !      kinetic energy
285         DO jk = 1, jpkm1
286            DO jj = 2, jpjm1
287               DO ji = fs_2, fs_jpim1   ! vector opt.
288                  zztmp   = 1._wp / ( e1e2t(ji,jj) * e3t_n(ji,jj,jk) )
289                  zztmpx  = 0.5 * (  un(ji-1,jj,jk) * un(ji-1,jj,jk) * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)    &
290                     &             + un(ji  ,jj,jk) * un(ji  ,jj,jk) * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk) )  &
291                     &          *  zztmp 
292                  !
293                  zztmpy  = 0.5 * (  vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)    &
294                     &             + vn(ji,jj  ,jk) * vn(ji,jj  ,jk) * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) )  &
295                     &          *  zztmp 
296                  !
297                  rke(ji,jj,jk) = 0.5_wp * ( zztmpx + zztmpy )
298                  !
299               ENDDO
300            ENDDO
301         ENDDO
302         CALL lbc_lnk( rke, 'T', 1. )
303         CALL iom_put( "eken", rke )           
304      ENDIF
305      !
306      CALL iom_put( "hdiv", hdivn )                  ! Horizontal divergence
307      !
308      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
309         z3d(:,:,jpk) = 0.e0
310         z2d(:,:) = 0.e0
311         DO jk = 1, jpkm1
312            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)
313            z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
314         END DO
315         CALL iom_put( "u_masstr", z3d )                  ! mass transport in i-direction
316         CALL iom_put( "u_masstr_vint", z2d )             ! mass transport in i-direction vertical sum
317      ENDIF
318     
319      IF( iom_use("u_heattr") ) THEN
320         z2d(:,:) = 0.e0 
321         DO jk = 1, jpkm1
322            DO jj = 2, jpjm1
323               DO ji = fs_2, fs_jpim1   ! vector opt.
324                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
325               END DO
326            END DO
327         END DO
328         CALL lbc_lnk( z2d, 'U', -1. )
329         CALL iom_put( "u_heattr", (0.5 * rcp) * z2d )    ! heat transport in i-direction
330      ENDIF
331
332      IF( iom_use("u_salttr") ) THEN
333         z2d(:,:) = 0.e0 
334         DO jk = 1, jpkm1
335            DO jj = 2, jpjm1
336               DO ji = fs_2, fs_jpim1   ! vector opt.
337                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
338               END DO
339            END DO
340         END DO
341         CALL lbc_lnk( z2d, 'U', -1. )
342         CALL iom_put( "u_salttr", 0.5 * z2d )            ! heat transport in i-direction
343      ENDIF
344
345     
346      IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
347         z3d(:,:,jpk) = 0.e0
348         DO jk = 1, jpkm1
349            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)
350         END DO
351         CALL iom_put( "v_masstr", z3d )                  ! mass transport in j-direction
352      ENDIF
353     
354      IF( iom_use("v_heattr") ) THEN
355         z2d(:,:) = 0.e0 
356         DO jk = 1, jpkm1
357            DO jj = 2, jpjm1
358               DO ji = fs_2, fs_jpim1   ! vector opt.
359                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
360               END DO
361            END DO
362         END DO
363         CALL lbc_lnk( z2d, 'V', -1. )
364         CALL iom_put( "v_heattr", (0.5 * rcp) * z2d )    !  heat transport in j-direction
365      ENDIF
366
367      IF( iom_use("v_salttr") ) THEN
368         z2d(:,:) = 0.e0 
369         DO jk = 1, jpkm1
370            DO jj = 2, jpjm1
371               DO ji = fs_2, fs_jpim1   ! vector opt.
372                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
373               END DO
374            END DO
375         END DO
376         CALL lbc_lnk( z2d, 'V', -1. )
377         CALL iom_put( "v_salttr", 0.5 * z2d )            !  heat transport in j-direction
378      ENDIF
379
380      ! Vertical integral of temperature
381      IF( iom_use("tosmint") ) THEN
382         z2d(:,:)=0._wp
383         DO jk = 1, jpkm1
384            DO jj = 2, jpjm1
385               DO ji = fs_2, fs_jpim1   ! vector opt.
386                  z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem)
387               END DO
388            END DO
389         END DO
390         CALL lbc_lnk( z2d, 'T', -1. )
391         CALL iom_put( "tosmint", z2d ) 
392      ENDIF
393
394      ! Vertical integral of salinity
395      IF( iom_use("somint") ) THEN
396         z2d(:,:)=0._wp
397         DO jk = 1, jpkm1
398            DO jj = 2, jpjm1
399               DO ji = fs_2, fs_jpim1   ! vector opt.
400                  z2d(ji,jj) = z2d(ji,jj) + rau0 * e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
401               END DO
402            END DO
403         END DO
404         CALL lbc_lnk( z2d, 'T', -1. )
405         CALL iom_put( "somint", z2d ) 
406      ENDIF
407
408      CALL iom_put( "bn2", rn2 )  !Brunt-Vaisala buoyancy frequency (N^2)
409      !
410      CALL wrk_dealloc( jpi , jpj      , z2d )
411      CALL wrk_dealloc( jpi , jpj, jpk , z3d )
412      !
413      ! If we want tmb values
414
415      IF (ln_diatmb) THEN
416         CALL dia_tmb 
417      ENDIF
418      IF (ln_dia25h) THEN
419         CALL dia_25h( kt )
420      ENDIF
421
422      IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
423      !
424   END SUBROUTINE dia_wri
425
426#else
427   !!----------------------------------------------------------------------
428   !!   Default option                                  use IOIPSL  library
429   !!----------------------------------------------------------------------
430
431   SUBROUTINE dia_wri( kt )
432      !!---------------------------------------------------------------------
433      !!                  ***  ROUTINE dia_wri  ***
434      !!                   
435      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
436      !!      NETCDF format is used by default
437      !!
438      !! ** Method  :   At the beginning of the first time step (nit000),
439      !!      define all the NETCDF files and fields
440      !!      At each time step call histdef to compute the mean if ncessary
441      !!      Each nwrite time step, output the instantaneous or mean fields
442      !!----------------------------------------------------------------------
443      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
444      !
445      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
446      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
447      INTEGER  ::   inum = 11                                ! temporary logical unit
448      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
449      INTEGER  ::   ierr                                     ! error code return from allocation
450      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
451      INTEGER  ::   jn, ierror                               ! local integers
452      REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
453      !
454      REAL(wp), POINTER, DIMENSION(:,:)   :: zw2d       ! 2D workspace
455      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d       ! 3D workspace
456      !!----------------------------------------------------------------------
457      !
458      IF( nn_timing == 1 )   CALL timing_start('dia_wri')
459      !
460                             CALL wrk_alloc( jpi,jpj      , zw2d )
461      IF( .NOT.ln_linssh )   CALL wrk_alloc( jpi,jpj,jpk  , zw3d )
462      !
463      ! Output the initial state and forcings
464      IF( ninist == 1 ) THEN                       
465         CALL dia_wri_state( 'output.init', kt )
466         ninist = 0
467      ENDIF
468      !
469      ! 0. Initialisation
470      ! -----------------
471
472      ! local variable for debugging
473      ll_print = .FALSE.
474      ll_print = ll_print .AND. lwp
475
476      ! Define frequency of output and means
477      clop = "x"         ! no use of the mask value (require less cpu time and otherwise the model crashes)
478#if defined key_diainstant
479      zsto = nwrite * rdt
480      clop = "inst("//TRIM(clop)//")"
481#else
482      zsto=rdt
483      clop = "ave("//TRIM(clop)//")"
484#endif
485      zout = nwrite * rdt
486      zmax = ( nitend - nit000 + 1 ) * rdt
487
488      ! Define indices of the horizontal output zoom and vertical limit storage
489      iimi = 1      ;      iima = jpi
490      ijmi = 1      ;      ijma = jpj
491      ipk = jpk
492
493      ! define time axis
494      it = kt
495      itmod = kt - nit000 + 1
496
497
498      ! 1. Define NETCDF files and fields at beginning of first time step
499      ! -----------------------------------------------------------------
500
501      IF( kt == nit000 ) THEN
502
503         ! Define the NETCDF files (one per grid)
504
505         ! Compute julian date from starting date of the run
506         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
507         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
508         IF(lwp)WRITE(numout,*)
509         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
510            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
511         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
512                                 ' limit storage in depth = ', ipk
513
514         ! WRITE root name in date.file for use by postpro
515         IF(lwp) THEN
516            CALL dia_nam( clhstnam, nwrite,' ' )
517            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
518            WRITE(inum,*) clhstnam
519            CLOSE(inum)
520         ENDIF
521
522         ! Define the T grid FILE ( nid_T )
523
524         CALL dia_nam( clhstnam, nwrite, 'grid_T' )
525         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
526         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
527            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
528            &          nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
529         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
530            &           "m", ipk, gdept_1d, nz_T, "down" )
531         !                                                            ! Index of ocean points
532         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
533         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
534         !
535         IF( ln_icebergs ) THEN
536            !
537            !! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
538            !! that routine is called from nemogcm, so do it here immediately before its needed
539            ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
540            IF( lk_mpp )   CALL mpp_sum( ierror )
541            IF( ierror /= 0 ) THEN
542               CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
543               RETURN
544            ENDIF
545            !
546            !! iceberg vertical coordinate is class number
547            CALL histvert( nid_T, "class", "Iceberg class",      &  ! Vertical grid: class
548               &           "number", nclasses, class_num, nb_T )
549            !
550            !! each class just needs the surface index pattern
551            ndim_bT = 3
552            DO jn = 1,nclasses
553               ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
554            ENDDO
555            !
556         ENDIF
557
558         ! Define the U grid FILE ( nid_U )
559
560         CALL dia_nam( clhstnam, nwrite, 'grid_U' )
561         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
562         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
563            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
564            &          nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
565         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
566            &           "m", ipk, gdept_1d, nz_U, "down" )
567         !                                                            ! Index of ocean points
568         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
569         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
570
571         ! Define the V grid FILE ( nid_V )
572
573         CALL dia_nam( clhstnam, nwrite, 'grid_V' )                   ! filename
574         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
575         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
576            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
577            &          nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
578         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
579            &          "m", ipk, gdept_1d, nz_V, "down" )
580         !                                                            ! Index of ocean points
581         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
582         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
583
584         ! Define the W grid FILE ( nid_W )
585
586         CALL dia_nam( clhstnam, nwrite, 'grid_W' )                   ! filename
587         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
588         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
589            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
590            &          nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
591         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
592            &          "m", ipk, gdepw_1d, nz_W, "down" )
593
594
595         ! Declare all the output fields as NETCDF variables
596
597         !                                                                                      !!! nid_T : 3D
598         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
599            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
600         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
601            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
602         IF(  .NOT.ln_linssh  ) THEN
603            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n
604            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
605            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t_n
606            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
607            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t_n
608            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
609         ENDIF
610         !                                                                                      !!! nid_T : 2D
611         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
612            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
613         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
614            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
615         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
616            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
617         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
618            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
619         CALL histdef( nid_T, "sorunoff", "River runoffs"                      , "Kg/m2/s",   &  ! runoffs
620            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
621         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
622            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
623         IF(  ln_linssh  ) THEN
624            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem)
625            &                                                                  , "KgC/m2/s",  &  ! sosst_cd
626            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
627            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * tsn(:,:,1,jp_sal)
628            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
629            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
630         ENDIF
631         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
632            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
633         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
634            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
635         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
636            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
637         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
638            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
639         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
640            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
641         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
642            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
643!
644         IF( ln_icebergs ) THEN
645            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
646               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
647            CALL histdef( nid_T, "calving_heat"        , "calving heat flux"                        , "XXXX"   , &
648               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
649            CALL histdef( nid_T, "berg_floating_melt"  , "Melt rate of icebergs + bits"             , "kg/m2/s", &
650               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
651            CALL histdef( nid_T, "berg_stored_ice"     , "Accumulated ice mass by class"            , "kg"     , &
652               &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
653            IF( ln_bergdia ) THEN
654               CALL histdef( nid_T, "berg_melt"           , "Melt rate of icebergs"                    , "kg/m2/s", &
655                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
656               CALL histdef( nid_T, "berg_buoy_melt"      , "Buoyancy component of iceberg melt rate"  , "kg/m2/s", &
657                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
658               CALL histdef( nid_T, "berg_eros_melt"      , "Erosion component of iceberg melt rate"   , "kg/m2/s", &
659                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
660               CALL histdef( nid_T, "berg_conv_melt"      , "Convective component of iceberg melt rate", "kg/m2/s", &
661                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
662               CALL histdef( nid_T, "berg_virtual_area"   , "Virtual coverage by icebergs"             , "m2"     , &
663                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
664               CALL histdef( nid_T, "bits_src"           , "Mass source of bergy bits"                , "kg/m2/s", &
665                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
666               CALL histdef( nid_T, "bits_melt"          , "Melt rate of bergy bits"                  , "kg/m2/s", &
667                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
668               CALL histdef( nid_T, "bits_mass"          , "Bergy bit density field"                  , "kg/m2"  , &
669                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
670               CALL histdef( nid_T, "berg_mass"           , "Iceberg density field"                    , "kg/m2"  , &
671                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
672               CALL histdef( nid_T, "berg_real_calving"   , "Calving into iceberg class"               , "kg/s"   , &
673                  &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
674            ENDIF
675         ENDIF
676
677         IF( .NOT. ln_cpl ) THEN
678            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
679               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
680            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
681               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
682            CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
683               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
684         ENDIF
685
686         IF( ln_cpl .AND. nn_ice <= 1 ) THEN
687            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
688               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
689            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
690               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
691            CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn
692               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
693         ENDIF
694         
695         clmx ="l_max(only(x))"    ! max index on a period
696!         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
697!            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
698#if defined key_diahth
699         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
700            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
701         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
702            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
703         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
704            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
705         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "J/m2"   ,   & ! htc3
706            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
707#endif
708
709         IF( ln_cpl .AND. nn_ice == 2 ) THEN
710            CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature"            , "K"      ,   &  ! tn_ice
711               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
712            CALL histdef( nid_T,"soicealb" , "Ice Albedo"                         , "[0,1]"  ,   &  ! alb_ice
713               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
714         ENDIF
715
716         CALL histend( nid_T, snc4chunks=snc4set )
717
718         !                                                                                      !!! nid_U : 3D
719         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
720            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
721         IF( ln_wave .AND. ln_sdw) THEN
722            CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current"         , "m/s"    ,   &  ! usd
723               &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
724         ENDIF
725         !                                                                                      !!! nid_U : 2D
726         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
727            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
728
729         CALL histend( nid_U, snc4chunks=snc4set )
730
731         !                                                                                      !!! nid_V : 3D
732         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
733            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
734         IF( ln_wave .AND. ln_sdw) THEN
735            CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current"    , "m/s"    ,   &  ! vsd
736               &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
737         ENDIF
738         !                                                                                      !!! nid_V : 2D
739         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
740            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
741
742         CALL histend( nid_V, snc4chunks=snc4set )
743
744         !                                                                                      !!! nid_W : 3D
745         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn
746            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
747         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
748            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
749         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avmu
750            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
751
752         IF( lk_zdfddm ) THEN
753            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
754               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
755         ENDIF
756         
757         IF( ln_wave .AND. ln_sdw) THEN
758            CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current"   , "m/s"    ,   &  ! wsd
759               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
760         ENDIF
761         !                                                                                      !!! nid_W : 2D
762         CALL histend( nid_W, snc4chunks=snc4set )
763
764         IF(lwp) WRITE(numout,*)
765         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
766         IF(ll_print) CALL FLUSH(numout )
767
768      ENDIF
769
770      ! 2. Start writing data
771      ! ---------------------
772
773      ! ndex(1) est utilise ssi l'avant dernier argument est different de
774      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
775      ! donne le nombre d'elements, et ndex la liste des indices a sortir
776
777      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
778         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
779         WRITE(numout,*) '~~~~~~ '
780      ENDIF
781
782      IF( .NOT.ln_linssh ) THEN
783         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
784         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
785         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
786         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
787      ELSE
788         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature
789         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity
790         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature
791         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
792      ENDIF
793      IF( .NOT.ln_linssh ) THEN
794         zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
795         CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
796         CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
797         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
798      ENDIF
799      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
800      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
801      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
802      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
803                                                                                  ! (includes virtual salt flux beneath ice
804                                                                                  ! in linear free surface case)
805      IF( ln_linssh ) THEN
806         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
807         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
808         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
809         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
810      ENDIF
811      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
812      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
813      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
814      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
815      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
816      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
817!
818      IF( ln_icebergs ) THEN
819         !
820         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT ) 
821         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
822         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT ) 
823         !
824         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
825         !
826         IF( ln_bergdia ) THEN
827            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   ) 
828            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   ) 
829            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   ) 
830            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   ) 
831            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   ) 
832            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   ) 
833            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   ) 
834            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   ) 
835            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   ) 
836            !
837            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
838         ENDIF
839      ENDIF
840
841      IF( .NOT. ln_cpl ) THEN
842         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
843         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
844         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
845         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
846      ENDIF
847      IF( ln_cpl .AND. nn_ice <= 1 ) THEN
848         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
849         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
850         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
851         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
852      ENDIF
853!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
854!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
855
856#if defined key_diahth
857      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
858      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
859      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
860      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
861#endif
862
863      IF( ln_cpl .AND. nn_ice == 2 ) THEN
864         CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature
865         CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo
866      ENDIF
867
868      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
869      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
870
871      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
872      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
873
874      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
875      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
876      CALL histwrite( nid_W, "votkeavm", it, avmu           , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
877      IF( lk_zdfddm ) THEN
878         CALL histwrite( nid_W, "voddmavs", it, fsavs(:,:,:), ndim_T, ndex_T )    ! S vert. eddy diff. coef.
879      ENDIF
880
881      IF( ln_wave .AND. ln_sdw ) THEN
882         CALL histwrite( nid_U, "sdzocrtx", it, usd           , ndim_U , ndex_U )    ! i-StokesDrift-current
883         CALL histwrite( nid_V, "sdmecrty", it, vsd           , ndim_V , ndex_V )    ! j-StokesDrift-current
884         CALL histwrite( nid_W, "sdvecrtz", it, wsd           , ndim_T , ndex_T )    ! StokesDrift vert. current
885      ENDIF
886
887      ! 3. Close all files
888      ! ---------------------------------------
889      IF( kt == nitend ) THEN
890         CALL histclo( nid_T )
891         CALL histclo( nid_U )
892         CALL histclo( nid_V )
893         CALL histclo( nid_W )
894      ENDIF
895      !
896                             CALL wrk_dealloc( jpi , jpj        , zw2d )
897      IF( .NOT.ln_linssh )   CALL wrk_dealloc( jpi , jpj , jpk  , zw3d )
898      !
899      IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
900      !
901   END SUBROUTINE dia_wri
902#endif
903
904   SUBROUTINE dia_wri_state( cdfile_name, kt )
905      !!---------------------------------------------------------------------
906      !!                 ***  ROUTINE dia_wri_state  ***
907      !!       
908      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
909      !!      the instantaneous ocean state and forcing fields.
910      !!        Used to find errors in the initial state or save the last
911      !!      ocean state in case of abnormal end of a simulation
912      !!
913      !! ** Method  :   NetCDF files using ioipsl
914      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
915      !!      File 'output.abort.nc' is created in case of abnormal job end
916      !!----------------------------------------------------------------------
917      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
918      INTEGER           , INTENT( in ) ::   kt               ! ocean time-step index
919      !!
920      CHARACTER (len=32) :: clname
921      CHARACTER (len=40) :: clop
922      INTEGER  ::   id_i , nz_i, nh_i       
923      INTEGER, DIMENSION(1) ::   idex             ! local workspace
924      REAL(wp) ::   zsto, zout, zmax, zjulian
925      !!----------------------------------------------------------------------
926      !
927      ! 0. Initialisation
928      ! -----------------
929
930      ! Define name, frequency of output and means
931      clname = cdfile_name
932      IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
933      zsto = rdt
934      clop = "inst(x)"           ! no use of the mask value (require less cpu time)
935      zout = rdt
936      zmax = ( nitend - nit000 + 1 ) * rdt
937
938      IF(lwp) WRITE(numout,*)
939      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
940      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
941      IF(lwp) WRITE(numout,*) '                and named :', clname, '.nc'
942
943
944      ! 1. Define NETCDF files and fields at beginning of first time step
945      ! -----------------------------------------------------------------
946
947      ! Compute julian date from starting date of the run
948      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )         ! time axis
949      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
950      CALL histbeg( clname, jpi, glamt, jpj, gphit,   &
951          1, jpi, 1, jpj, nit000-1, zjulian, rdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
952      CALL histvert( id_i, "deptht", "Vertical T levels",   &    ! Vertical grid : gdept
953          "m", jpk, gdept_1d, nz_i, "down")
954
955      ! Declare all the output fields as NetCDF variables
956
957      CALL histdef( id_i, "vosaline", "Salinity"              , "PSU"    ,   &   ! salinity
958         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
959      CALL histdef( id_i, "votemper", "Temperature"           , "C"      ,   &   ! temperature
960         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
961      CALL histdef( id_i, "sossheig", "Sea Surface Height"    , "m"      ,   &  ! ssh
962         &          jpi, jpj, nh_i, 1  , 1, 1  , nz_i, 32, clop, zsto, zout )
963      CALL histdef( id_i, "vozocrtx", "Zonal Current"         , "m/s"    ,   &   ! zonal current
964         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
965      CALL histdef( id_i, "vomecrty", "Meridional Current"    , "m/s"    ,   &   ! meridonal current
966         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
967      CALL histdef( id_i, "vovecrtz", "Vertical Velocity"     , "m/s"    ,   &   ! vertical current
968         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
969         !
970      IF( ALLOCATED(ahtu) ) THEN
971         CALL histdef( id_i, "ahtu"    , "u-eddy diffusivity"    , "m2/s"    ,   &   ! zonal current
972            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
973         CALL histdef( id_i, "ahtv"    , "v-eddy diffusivity"    , "m2/s"    ,   &   ! meridonal current
974            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
975      ENDIF
976      IF( ALLOCATED(ahmt) ) THEN
977         CALL histdef( id_i, "ahmt"    , "t-eddy viscosity"      , "m2/s"    ,   &   ! zonal current
978            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
979         CALL histdef( id_i, "ahmf"    , "f-eddy viscosity"      , "m2/s"    ,   &   ! meridonal current
980            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
981      ENDIF
982         !
983      CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S",   &   ! net freshwater
984         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
985      CALL histdef( id_i, "sohefldo", "Net Downward Heat Flux", "W/m2"   ,   &   ! net heat flux
986         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
987      CALL histdef( id_i, "soshfldo", "Shortwave Radiation"   , "W/m2"   ,   &   ! solar flux
988         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
989      CALL histdef( id_i, "soicecov", "Ice fraction"          , "[0,1]"  ,   &   ! fr_i
990         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
991      CALL histdef( id_i, "sozotaux", "Zonal Wind Stress"     , "N/m2"   ,   &   ! i-wind stress
992         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
993      CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2"   ,   &   ! j-wind stress
994         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
995      IF( .NOT.ln_linssh ) THEN
996         CALL histdef( id_i, "vovvldep", "T point depth"         , "m"      , &   ! t-point depth
997            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
998         CALL histdef( id_i, "vovvle3t", "T point thickness"     , "m"      , &   ! t-point depth
999            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1000      ENDIF
1001      !
1002      IF( ln_wave .AND. ln_sdw ) THEN
1003         CALL histdef( id_i, "sdzocrtx", "Stokes Drift Zonal"    , "m/s"    , &   ! StokesDrift zonal current
1004            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1005         CALL histdef( id_i, "sdmecrty", "Stokes Drift Merid"    , "m/s"    , &   ! StokesDrift meridonal current
1006            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1007         CALL histdef( id_i, "sdvecrtz", "Stokes Drift Vert"     , "m/s"    , &   ! StokesDrift vertical current
1008            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1009      ENDIF
1010
1011#if defined key_lim2
1012      CALL lim_wri_state_2( kt, id_i, nh_i )
1013#elif defined key_lim3
1014      CALL lim_wri_state( kt, id_i, nh_i )
1015#else
1016      CALL histend( id_i, snc4chunks=snc4set )
1017#endif
1018
1019      ! 2. Start writing data
1020      ! ---------------------
1021      ! idex(1) est utilise ssi l'avant dernier argument est diffferent de
1022      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
1023      ! donne le nombre d'elements, et idex la liste des indices a sortir
1024      idex(1) = 1   ! init to avoid compil warning
1025
1026      ! Write all fields on T grid
1027      CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature
1028      CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity
1029      CALL histwrite( id_i, "sossheig", kt, sshn             , jpi*jpj    , idex )    ! sea surface height
1030      CALL histwrite( id_i, "vozocrtx", kt, un               , jpi*jpj*jpk, idex )    ! now i-velocity
1031      CALL histwrite( id_i, "vomecrty", kt, vn               , jpi*jpj*jpk, idex )    ! now j-velocity
1032      CALL histwrite( id_i, "vovecrtz", kt, wn               , jpi*jpj*jpk, idex )    ! now k-velocity
1033      !
1034      IF( ALLOCATED(ahtu) ) THEN
1035         CALL histwrite( id_i, "ahtu"    , kt, ahtu             , jpi*jpj*jpk, idex )    ! aht at u-point
1036         CALL histwrite( id_i, "ahtv"    , kt, ahtv             , jpi*jpj*jpk, idex )    !  -  at v-point
1037      ENDIF
1038      IF( ALLOCATED(ahmt) ) THEN
1039         CALL histwrite( id_i, "ahmt"    , kt, ahmt             , jpi*jpj*jpk, idex )    ! ahm at t-point
1040         CALL histwrite( id_i, "ahmf"    , kt, ahmf             , jpi*jpj*jpk, idex )    !  -  at f-point
1041      ENDIF
1042      !
1043      CALL histwrite( id_i, "sowaflup", kt, emp - rnf        , jpi*jpj    , idex )    ! freshwater budget
1044      CALL histwrite( id_i, "sohefldo", kt, qsr + qns        , jpi*jpj    , idex )    ! total heat flux
1045      CALL histwrite( id_i, "soshfldo", kt, qsr              , jpi*jpj    , idex )    ! solar heat flux
1046      CALL histwrite( id_i, "soicecov", kt, fr_i             , jpi*jpj    , idex )    ! ice fraction
1047      CALL histwrite( id_i, "sozotaux", kt, utau             , jpi*jpj    , idex )    ! i-wind stress
1048      CALL histwrite( id_i, "sometauy", kt, vtau             , jpi*jpj    , idex )    ! j-wind stress
1049
1050      IF(  .NOT.ln_linssh  ) THEN             
1051         CALL histwrite( id_i, "vovvldep", kt, gdept_n(:,:,:), jpi*jpj*jpk, idex )!  T-cell depth
1052         CALL histwrite( id_i, "vovvle3t", kt, e3t_n (:,:,:) , jpi*jpj*jpk, idex )!  T-cell thickness 
1053      END IF
1054 
1055      IF( ln_wave .AND. ln_sdw ) THEN
1056         CALL histwrite( id_i, "sdzocrtx", kt, usd           , jpi*jpj*jpk, idex)     ! now StokesDrift i-velocity
1057         CALL histwrite( id_i, "sdmecrty", kt, vsd           , jpi*jpj*jpk, idex)     ! now StokesDrift j-velocity
1058         CALL histwrite( id_i, "sdvecrtz", kt, wsd           , jpi*jpj*jpk, idex)     ! now StokesDrift k-velocity
1059      ENDIF
1060
1061      ! 3. Close the file
1062      ! -----------------
1063      CALL histclo( id_i )
1064#if ! defined key_iomput
1065      IF( ninist /= 1  ) THEN
1066         CALL histclo( nid_T )
1067         CALL histclo( nid_U )
1068         CALL histclo( nid_V )
1069         CALL histclo( nid_W )
1070      ENDIF
1071#endif
1072      !
1073   END SUBROUTINE dia_wri_state
1074
1075   !!======================================================================
1076END MODULE diawri
Note: See TracBrowser for help on using the repository browser.