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 branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/UKMO/dev_r6501_GO6_package_trunk/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 6825

Last change on this file since 6825 was 6825, checked in by timgraham, 8 years ago

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