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_r5518_GC3p0_package_v2/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/UKMO/dev_r5518_GC3p0_package_v2/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 6679

Last change on this file since 6679 was 6679, checked in by malcolmroberts, 8 years ago

Merged in changes from v3_6_extra_CMIP6_diagnostics up to revision 6674
File size: 56.5 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   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   dia_wri       : create the standart output files
23   !!   dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
24   !!----------------------------------------------------------------------
25   USE oce             ! ocean dynamics and tracers
26   USE dom_oce         ! ocean space and time domain
27   USE dynadv, ONLY: ln_dynadv_vec
28   USE zdf_oce         ! ocean vertical physics
29   USE ldftra_oce      ! ocean active tracers: lateral physics
30   USE ldfdyn_oce      ! ocean dynamics: lateral physics
31   USE traldf_iso_grif, ONLY : psix_eiv, psiy_eiv
32   USE sol_oce         ! solver variables
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 lbclnk          ! ocean lateral boundary conditions (or mpp link)
44   USE in_out_manager  ! I/O manager
45   USE diadimg         ! dimg direct access file format output
46   USE iom
47   USE ioipsl
48   USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities     
49
50#if defined key_lim2
51   USE limwri_2 
52#elif defined key_lim3
53   USE limwri 
54#endif
55   USE lib_mpp         ! MPP library
56   USE timing          ! preformance summary
57   USE wrk_nemo        ! working array
58
59   IMPLICIT NONE
60   PRIVATE
61
62   PUBLIC   dia_wri                 ! routines called by step.F90
63   PUBLIC   dia_wri_state
64   PUBLIC   dia_wri_alloc           ! Called by nemogcm module
65
66   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
67   INTEGER ::          nb_T              , ndim_bT   ! grid_T file
68   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
69   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
70   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file
71   INTEGER ::   ndex(1)                              ! ???
72   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
73   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
74   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
75
76   !! * Substitutions
77#  include "zdfddm_substitute.h90"
78#  include "domzgr_substitute.h90"
79#  include "vectopt_loop_substitute.h90"
80   !!----------------------------------------------------------------------
81   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
82   !! $Id$
83   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
84   !!----------------------------------------------------------------------
85CONTAINS
86
87   INTEGER FUNCTION dia_wri_alloc()
88      !!----------------------------------------------------------------------
89      INTEGER, DIMENSION(2) :: ierr
90      !!----------------------------------------------------------------------
91      ierr = 0
92      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
93         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
94         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
95         !
96      dia_wri_alloc = MAXVAL(ierr)
97      IF( lk_mpp )   CALL mpp_sum( dia_wri_alloc )
98      !
99  END FUNCTION dia_wri_alloc
100
101#if defined key_dimgout
102   !!----------------------------------------------------------------------
103   !!   'key_dimgout'                                      DIMG output file
104   !!----------------------------------------------------------------------
105#   include "diawri_dimg.h90"
106
107#else
108   !!----------------------------------------------------------------------
109   !!   Default option                                   NetCDF output file
110   !!----------------------------------------------------------------------
111# if defined key_iomput
112   !!----------------------------------------------------------------------
113   !!   'key_iomput'                                        use IOM library
114   !!----------------------------------------------------------------------
115
116   SUBROUTINE dia_wri( kt )
117      !!---------------------------------------------------------------------
118      !!                  ***  ROUTINE dia_wri  ***
119      !!                   
120      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
121      !!      NETCDF format is used by default
122      !!
123      !! ** Method  :  use iom_put
124      !!----------------------------------------------------------------------
125      !!
126      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
127      !!
128      INTEGER                      ::   ji, jj, jk              ! dummy loop indices
129      INTEGER                      ::   jkbot                   !
130      REAL(wp)                     ::   zztmp, zztmpx, zztmpy   !
131      !!
132      REAL(wp), POINTER, DIMENSION(:,:)   :: z2d      ! 2D workspace
133      REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d      ! 3D workspace
134      !!----------------------------------------------------------------------
135      !
136      IF( nn_timing == 1 )   CALL timing_start('dia_wri')
137      !
138      CALL wrk_alloc( jpi , jpj      , z2d )
139      CALL wrk_alloc( jpi , jpj, jpk , z3d )
140      !
141      ! Output the initial state and forcings
142      IF( ninist == 1 ) THEN                       
143         CALL dia_wri_state( 'output.init', kt )
144         ninist = 0
145      ENDIF
146
147      IF( .NOT.lk_vvl ) THEN
148         CALL iom_put( "e3t" , fse3t_n(:,:,:) )
149         CALL iom_put( "e3u" , fse3u_n(:,:,:) )
150         CALL iom_put( "e3v" , fse3v_n(:,:,:) )
151         CALL iom_put( "e3w" , fse3w_n(:,:,:) )
152      ENDIF
153
154      CALL iom_put( "ssh" , sshn )                 ! sea surface height
155      if( iom_use('ssh2') )   CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) )   ! square of 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#if defined key_dynspg_ts
209      CALL iom_put(  "ubar", un_adv(:,:)      )    ! barotropic i-current
210#else
211      CALL iom_put(  "ubar", un_b(:,:)        )    ! barotropic i-current
212#endif
213     
214      CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current
215      CALL iom_put(  "ssv", vn(:,:,1)         )    ! surface j-current
216      IF ( iom_use("sbv") ) THEN
217         DO jj = 1, jpj
218            DO ji = 1, jpi
219               jkbot = mbkv(ji,jj)
220               z2d(ji,jj) = vn(ji,jj,jkbot)
221            END DO
222         END DO
223         CALL iom_put( "sbv", z2d )                ! bottom j-current
224      ENDIF
225#if defined key_dynspg_ts
226      CALL iom_put(  "vbar", vn_adv(:,:)      )    ! barotropic j-current
227#else
228      CALL iom_put(  "vbar", vn_b(:,:)        )    ! barotropic j-current
229#endif
230
231      CALL iom_put( "woce", wn )                   ! vertical velocity
232      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
233         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
234         z2d(:,:) = rau0 * e12t(:,:)
235         DO jk = 1, jpk
236            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
237         END DO
238         CALL iom_put( "w_masstr" , z3d ) 
239         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
240      ENDIF
241
242      CALL iom_put( "avt" , avt                        )    ! T vert. eddy diff. coef.
243      CALL iom_put( "avm" , avmu                       )    ! T vert. eddy visc. coef.
244      CALL iom_put( "avs" , fsavs(:,:,:)               )    ! S vert. eddy diff. coef. (useful only with key_zdfddm)
245
246      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
247         DO jj = 2, jpjm1                                    ! sst gradient
248            DO ji = fs_2, fs_jpim1   ! vector opt.
249               zztmp      = tsn(ji,jj,1,jp_tem)
250               zztmpx     = ( tsn(ji+1,jj  ,1,jp_tem) - zztmp ) / e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,1,jp_tem) ) / e1u(ji-1,jj  )
251               zztmpy     = ( tsn(ji  ,jj+1,1,jp_tem) - zztmp ) / e2v(ji,jj) + ( zztmp - tsn(ji  ,jj-1,1,jp_tem) ) / e2v(ji  ,jj-1)
252               z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
253                  &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
254            END DO
255         END DO
256         CALL lbc_lnk( z2d, 'T', 1. )
257         CALL iom_put( "sstgrad2",  z2d               )    ! square of module of sst gradient
258         z2d(:,:) = SQRT( z2d(:,:) )
259         CALL iom_put( "sstgrad" ,  z2d               )    ! module of sst gradient
260      ENDIF
261         
262      ! clem: heat and salt content
263      IF( iom_use("heatc") ) THEN
264         z2d(:,:)  = 0._wp 
265         DO jk = 1, jpkm1
266            DO jj = 1, jpj
267               DO ji = 1, jpi
268                  z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
269               END DO
270            END DO
271         END DO
272         CALL iom_put( "heatc", (rau0 * rcp) * z2d )    ! vertically integrated heat content (J/m2)
273      ENDIF
274
275      IF( iom_use("saltc") ) THEN
276         z2d(:,:)  = 0._wp 
277         DO jk = 1, jpkm1
278            DO jj = 1, jpj
279               DO ji = 1, jpi
280                  z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
281               END DO
282            END DO
283         END DO
284         CALL iom_put( "saltc", rau0 * z2d )   ! vertically integrated salt content (PSU*kg/m2)
285      ENDIF
286      !
287      IF ( iom_use("eken") ) THEN
288         rke(:,:,jk) = 0._wp                               !      kinetic energy
289         DO jk = 1, jpkm1
290            DO jj = 2, jpjm1
291               DO ji = fs_2, fs_jpim1   ! vector opt.
292                  zztmp   = 1._wp / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
293                  zztmpx  = 0.5 * (  un(ji-1,jj,jk) * un(ji-1,jj,jk) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk)    &
294                     &             + un(ji  ,jj,jk) * un(ji  ,jj,jk) * e2u(ji  ,jj) * fse3u(ji  ,jj,jk) )  &
295                     &          *  zztmp 
296                  !
297                  zztmpy  = 0.5 * (  vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1v(ji,jj-1) * fse3v(ji,jj-1,jk)    &
298                     &             + vn(ji,jj  ,jk) * vn(ji,jj  ,jk) * e1v(ji,jj  ) * fse3v(ji,jj  ,jk) )  &
299                     &          *  zztmp 
300                  !
301                  rke(ji,jj,jk) = 0.5_wp * ( zztmpx + zztmpy )
302                  !
303               ENDDO
304            ENDDO
305         ENDDO
306         CALL lbc_lnk( rke, 'T', 1. )
307         CALL iom_put( "eken", rke )           
308      ENDIF
309         
310      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
311         z3d(:,:,jpk) = 0.e0
312         z2d(:,:) = 0.e0
313         DO jk = 1, jpkm1
314            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) * umask(:,:,jk)
315            z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
316         END DO
317         CALL iom_put( "u_masstr", z3d )                  ! mass transport in i-direction
318         CALL iom_put( "u_masstr_vint", z2d )             ! mass transport in i-direction vertical sum
319      ENDIF
320     
321      IF( iom_use("u_heattr") ) THEN
322         z2d(:,:) = 0.e0 
323         DO jk = 1, jpkm1
324            DO jj = 2, jpjm1
325               DO ji = fs_2, fs_jpim1   ! vector opt.
326                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
327               END DO
328            END DO
329         END DO
330         CALL lbc_lnk( z2d, 'U', -1. )
331         CALL iom_put( "u_heattr", (0.5 * rcp) * z2d )    ! heat transport in i-direction
332      ENDIF
333
334      IF( iom_use("u_salttr") ) THEN
335         z2d(:,:) = 0.e0 
336         DO jk = 1, jpkm1
337            DO jj = 2, jpjm1
338               DO ji = fs_2, fs_jpim1   ! vector opt.
339                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
340               END DO
341            END DO
342         END DO
343         CALL lbc_lnk( z2d, 'U', -1. )
344         CALL iom_put( "u_salttr", 0.5 * z2d )            ! heat transport in i-direction
345      ENDIF
346
347     
348      IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
349         z3d(:,:,jpk) = 0.e0
350         DO jk = 1, jpkm1
351            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) * vmask(:,:,jk)
352         END DO
353         CALL iom_put( "v_masstr", z3d )                  ! mass transport in j-direction
354      ENDIF
355     
356      IF( iom_use("v_heattr") ) THEN
357         z2d(:,:) = 0.e0 
358         DO jk = 1, jpkm1
359            DO jj = 2, jpjm1
360               DO ji = fs_2, fs_jpim1   ! vector opt.
361                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
362               END DO
363            END DO
364         END DO
365         CALL lbc_lnk( z2d, 'V', -1. )
366         CALL iom_put( "v_heattr", (0.5 * rcp) * z2d )    !  heat transport in j-direction
367      ENDIF
368
369      IF( iom_use("v_salttr") ) THEN
370         z2d(:,:) = 0.e0 
371         DO jk = 1, jpkm1
372            DO jj = 2, jpjm1
373               DO ji = fs_2, fs_jpim1   ! vector opt.
374                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
375               END DO
376            END DO
377         END DO
378         CALL lbc_lnk( z2d, 'V', -1. )
379         CALL iom_put( "v_salttr", 0.5 * z2d )            !  heat transport in j-direction
380      ENDIF
381
382      CALL iom_put( "bn2", rn2 )  !Brunt-Vaisala buoyancy frequency (N^2)
383      !
384      CALL wrk_dealloc( jpi , jpj      , z2d )
385      CALL wrk_dealloc( jpi , jpj, jpk , z3d )
386      !
387      IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
388      !
389   END SUBROUTINE dia_wri
390
391#else
392   !!----------------------------------------------------------------------
393   !!   Default option                                  use IOIPSL  library
394   !!----------------------------------------------------------------------
395
396   SUBROUTINE dia_wri( kt )
397      !!---------------------------------------------------------------------
398      !!                  ***  ROUTINE dia_wri  ***
399      !!                   
400      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
401      !!      NETCDF format is used by default
402      !!
403      !! ** Method  :   At the beginning of the first time step (nit000),
404      !!      define all the NETCDF files and fields
405      !!      At each time step call histdef to compute the mean if ncessary
406      !!      Each nwrite time step, output the instantaneous or mean fields
407      !!----------------------------------------------------------------------
408      !!
409      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
410      !!
411      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
412      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
413      INTEGER  ::   inum = 11                                ! temporary logical unit
414      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
415      INTEGER  ::   ierr                                     ! error code return from allocation
416      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
417      INTEGER  ::   jn, ierror                               ! local integers
418      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt           ! local scalars
419      !!
420      REAL(wp), POINTER, DIMENSION(:,:)   :: zw2d       ! 2D workspace
421      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d       ! 3D workspace
422      !!----------------------------------------------------------------------
423      !
424      IF( nn_timing == 1 )   CALL timing_start('dia_wri')
425      !
426      CALL wrk_alloc( jpi , jpj      , zw2d )
427      IF ( ln_traldf_gdia .OR. lk_vvl )  call wrk_alloc( jpi , jpj , jpk  , zw3d )
428      !
429      ! Output the initial state and forcings
430      IF( ninist == 1 ) THEN                       
431         CALL dia_wri_state( 'output.init', kt )
432         ninist = 0
433      ENDIF
434      !
435      ! 0. Initialisation
436      ! -----------------
437
438      ! local variable for debugging
439      ll_print = .FALSE.
440      ll_print = ll_print .AND. lwp
441
442      ! Define frequency of output and means
443      zdt = rdt
444      IF( nacc == 1 ) zdt = rdtmin
445      IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!)
446      ELSE                      ;   clop = "x"         ! no use of the mask value (require less cpu time)
447      ENDIF
448#if defined key_diainstant
449      zsto = nwrite * zdt
450      clop = "inst("//TRIM(clop)//")"
451#else
452      zsto=zdt
453      clop = "ave("//TRIM(clop)//")"
454#endif
455      zout = nwrite * zdt
456      zmax = ( nitend - nit000 + 1 ) * zdt
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, zdt, 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, zdt, 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, zdt, 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, zdt, 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(  lk_vvl  ) 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(  .NOT. lk_vvl  ) 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         IF( ln_traldf_gdia ) THEN
692            CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
693                 &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
694         ELSE
695#if defined key_diaeiv
696            CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
697            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
698#endif
699         END IF
700         !                                                                                      !!! nid_U : 2D
701         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
702            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
703
704         CALL histend( nid_U, snc4chunks=snc4set )
705
706         !                                                                                      !!! nid_V : 3D
707         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
708            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
709         IF( ln_traldf_gdia ) THEN
710            CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
711                 &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
712         ELSE 
713#if defined key_diaeiv
714            CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
715            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
716#endif
717         END IF
718         !                                                                                      !!! nid_V : 2D
719         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
720            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
721
722         CALL histend( nid_V, snc4chunks=snc4set )
723
724         !                                                                                      !!! nid_W : 3D
725         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn
726            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
727         IF( ln_traldf_gdia ) THEN
728            CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
729                 &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
730         ELSE
731#if defined key_diaeiv
732            CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
733                 &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
734#endif
735         END IF
736         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
737            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
738         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avmu
739            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
740
741         IF( lk_zdfddm ) THEN
742            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
743               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
744         ENDIF
745         !                                                                                      !!! nid_W : 2D
746#if defined key_traldf_c2d
747         CALL histdef( nid_W, "soleahtw", "lateral eddy diffusivity"           , "m2/s"   ,   &  ! ahtw
748            &          jpi, jpj, nh_W, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
749# if defined key_traldf_eiv 
750            CALL histdef( nid_W, "soleaeiw", "eddy induced vel. coeff. at w-point", "m2/s",   &  ! aeiw
751               &       jpi, jpj, nh_W, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
752# endif
753#endif
754
755         CALL histend( nid_W, snc4chunks=snc4set )
756
757         IF(lwp) WRITE(numout,*)
758         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
759         IF(ll_print) CALL FLUSH(numout )
760
761      ENDIF
762
763      ! 2. Start writing data
764      ! ---------------------
765
766      ! ndex(1) est utilise ssi l'avant dernier argument est different de
767      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
768      ! donne le nombre d'elements, et ndex la liste des indices a sortir
769
770      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
771         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
772         WRITE(numout,*) '~~~~~~ '
773      ENDIF
774
775      IF( lk_vvl ) THEN
776         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
777         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
778         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
779         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
780      ELSE
781         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature
782         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity
783         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature
784         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
785      ENDIF
786      IF( lk_vvl ) THEN
787         zw3d(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
788         CALL histwrite( nid_T, "vovvle3t", it, fse3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
789         CALL histwrite( nid_T, "vovvldep", it, fsdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
790         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
791      ENDIF
792      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
793      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
794      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
795      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
796                                                                                  ! (includes virtual salt flux beneath ice
797                                                                                  ! in linear free surface case)
798      IF( .NOT. lk_vvl ) THEN
799         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
800         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
801         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
802         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
803      ENDIF
804      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
805      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
806      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
807      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
808      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
809      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
810!
811      IF( ln_icebergs ) THEN
812         !
813         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT ) 
814         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
815         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT ) 
816         !
817         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
818         !
819         IF( ln_bergdia ) THEN
820            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   ) 
821            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   ) 
822            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   ) 
823            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   ) 
824            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   ) 
825            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   ) 
826            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   ) 
827            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   ) 
828            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   ) 
829            !
830            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
831         ENDIF
832      ENDIF
833
834      IF( .NOT. ln_cpl ) THEN
835         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
836         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
837         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
838         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
839      ENDIF
840      IF( ln_cpl .AND. nn_ice <= 1 ) THEN
841         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
842         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
843         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
844         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
845      ENDIF
846!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
847!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
848
849#if defined key_diahth
850      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
851      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
852      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
853      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
854#endif
855
856      IF( ln_cpl .AND. nn_ice == 2 ) THEN
857         CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature
858         CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo
859      ENDIF
860
861      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
862      IF( ln_traldf_gdia ) THEN
863         IF (.not. ALLOCATED(psix_eiv))THEN
864            ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr )
865            IF( lk_mpp   )   CALL mpp_sum ( ierr )
866            IF( ierr > 0 )   CALL ctl_stop('STOP', 'diawri: unable to allocate psi{x,y}_eiv')
867            psix_eiv(:,:,:) = 0.0_wp
868            psiy_eiv(:,:,:) = 0.0_wp
869         ENDIF
870         DO jk=1,jpkm1
871            zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk)  ! u_eiv = -dpsix/dz
872         END DO
873         zw3d(:,:,jpk) = 0._wp
874         CALL histwrite( nid_U, "vozoeivu", it, zw3d, ndim_U , ndex_U )           ! i-eiv current
875      ELSE
876#if defined key_diaeiv
877         CALL histwrite( nid_U, "vozoeivu", it, u_eiv, ndim_U , ndex_U )          ! i-eiv current
878#endif
879      ENDIF
880      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
881
882      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
883      IF( ln_traldf_gdia ) THEN
884         DO jk=1,jpk-1
885            zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk)  ! v_eiv = -dpsiy/dz
886         END DO
887         zw3d(:,:,jpk) = 0._wp
888         CALL histwrite( nid_V, "vomeeivv", it, zw3d, ndim_V , ndex_V )           ! j-eiv current
889      ELSE
890#if defined key_diaeiv
891         CALL histwrite( nid_V, "vomeeivv", it, v_eiv, ndim_V , ndex_V )          ! j-eiv current
892#endif
893      ENDIF
894      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
895
896      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
897      IF( ln_traldf_gdia ) THEN
898         DO jk=1,jpk-1
899            DO jj = 2, jpjm1
900               DO ji = fs_2, fs_jpim1  ! vector opt.
901                  zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2v(ji,jj) + &
902                       &    (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1u(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx
903               END DO
904            END DO
905         END DO
906         zw3d(:,:,jpk) = 0._wp
907         CALL histwrite( nid_W, "voveeivw", it, zw3d          , ndim_T, ndex_T )    ! vert. eiv current
908      ELSE
909#   if defined key_diaeiv
910         CALL histwrite( nid_W, "voveeivw", it, w_eiv          , ndim_T, ndex_T )    ! vert. eiv current
911#   endif
912      ENDIF
913      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
914      CALL histwrite( nid_W, "votkeavm", it, avmu           , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
915      IF( lk_zdfddm ) THEN
916         CALL histwrite( nid_W, "voddmavs", it, fsavs(:,:,:), ndim_T, ndex_T )    ! S vert. eddy diff. coef.
917      ENDIF
918#if defined key_traldf_c2d
919      CALL histwrite( nid_W, "soleahtw", it, ahtw          , ndim_hT, ndex_hT )   ! lateral eddy diff. coef.
920# if defined key_traldf_eiv
921         CALL histwrite( nid_W, "soleaeiw", it, aeiw       , ndim_hT, ndex_hT )   ! EIV coefficient at w-point
922# endif
923#endif
924
925      ! 3. Close all files
926      ! ---------------------------------------
927      IF( kt == nitend ) THEN
928         CALL histclo( nid_T )
929         CALL histclo( nid_U )
930         CALL histclo( nid_V )
931         CALL histclo( nid_W )
932      ENDIF
933      !
934      CALL wrk_dealloc( jpi , jpj      , zw2d )
935      IF ( ln_traldf_gdia .OR. lk_vvl )  call wrk_dealloc( jpi , jpj , jpk  , zw3d )
936      !
937      IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
938      !
939   END SUBROUTINE dia_wri
940# endif
941
942#endif
943
944   SUBROUTINE dia_wri_state( cdfile_name, kt )
945      !!---------------------------------------------------------------------
946      !!                 ***  ROUTINE dia_wri_state  ***
947      !!       
948      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
949      !!      the instantaneous ocean state and forcing fields.
950      !!        Used to find errors in the initial state or save the last
951      !!      ocean state in case of abnormal end of a simulation
952      !!
953      !! ** Method  :   NetCDF files using ioipsl
954      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
955      !!      File 'output.abort.nc' is created in case of abnormal job end
956      !!----------------------------------------------------------------------
957      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
958      INTEGER           , INTENT( in ) ::   kt               ! ocean time-step index
959      !!
960      CHARACTER (len=32) :: clname
961      CHARACTER (len=40) :: clop
962      INTEGER  ::   id_i , nz_i, nh_i       
963      INTEGER, DIMENSION(1) ::   idex             ! local workspace
964      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt
965      !!----------------------------------------------------------------------
966      !
967!     IF( nn_timing == 1 )   CALL timing_start('dia_wri_state') ! not sure this works for routines not called in first timestep
968
969      ! 0. Initialisation
970      ! -----------------
971
972      ! Define name, frequency of output and means
973      clname = cdfile_name
974      IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
975      zdt  = rdt
976      zsto = rdt
977      clop = "inst(x)"           ! no use of the mask value (require less cpu time)
978      zout = rdt
979      zmax = ( nitend - nit000 + 1 ) * zdt
980
981      IF(lwp) WRITE(numout,*)
982      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
983      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
984      IF(lwp) WRITE(numout,*) '                and named :', clname, '.nc'
985
986
987      ! 1. Define NETCDF files and fields at beginning of first time step
988      ! -----------------------------------------------------------------
989
990      ! Compute julian date from starting date of the run
991      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )         ! time axis
992      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
993      CALL histbeg( clname, jpi, glamt, jpj, gphit,   &
994          1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
995      CALL histvert( id_i, "deptht", "Vertical T levels",   &    ! Vertical grid : gdept
996          "m", jpk, gdept_1d, nz_i, "down")
997
998      ! Declare all the output fields as NetCDF variables
999
1000      CALL histdef( id_i, "vosaline", "Salinity"              , "PSU"    ,   &   ! salinity
1001         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1002      CALL histdef( id_i, "votemper", "Temperature"           , "C"      ,   &   ! temperature
1003         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1004      CALL histdef( id_i, "sossheig", "Sea Surface Height"    , "m"      ,   &  ! ssh
1005         &          jpi, jpj, nh_i, 1  , 1, 1  , nz_i, 32, clop, zsto, zout )
1006      CALL histdef( id_i, "vozocrtx", "Zonal Current"         , "m/s"    ,   &   ! zonal current
1007         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1008      CALL histdef( id_i, "vomecrty", "Meridional Current"    , "m/s"    ,   &   ! meridonal current
1009         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
1010      CALL histdef( id_i, "vovecrtz", "Vertical Velocity"     , "m/s"    ,   &   ! vertical current
1011         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
1012      CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S",   &   ! net freshwater
1013         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1014      CALL histdef( id_i, "sohefldo", "Net Downward Heat Flux", "W/m2"   ,   &   ! net heat flux
1015         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1016      CALL histdef( id_i, "soshfldo", "Shortwave Radiation"   , "W/m2"   ,   &   ! solar flux
1017         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1018      CALL histdef( id_i, "soicecov", "Ice fraction"          , "[0,1]"  ,   &   ! fr_i
1019         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1020      CALL histdef( id_i, "sozotaux", "Zonal Wind Stress"     , "N/m2"   ,   &   ! i-wind stress
1021         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1022      CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2"   ,   &   ! j-wind stress
1023         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1024      IF( lk_vvl ) THEN
1025         CALL histdef( id_i, "vovvldep", "T point depth"         , "m"      ,   &   ! t-point depth
1026            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1027      END IF
1028
1029#if defined key_lim2
1030      CALL lim_wri_state_2( kt, id_i, nh_i )
1031#elif defined key_lim3
1032      CALL lim_wri_state( kt, id_i, nh_i )
1033#else
1034      CALL histend( id_i, snc4chunks=snc4set )
1035#endif
1036
1037      ! 2. Start writing data
1038      ! ---------------------
1039      ! idex(1) est utilise ssi l'avant dernier argument est diffferent de
1040      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
1041      ! donne le nombre d'elements, et idex la liste des indices a sortir
1042      idex(1) = 1   ! init to avoid compil warning
1043
1044      ! Write all fields on T grid
1045      CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature
1046      CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity
1047      CALL histwrite( id_i, "sossheig", kt, sshn             , jpi*jpj    , idex )    ! sea surface height
1048      CALL histwrite( id_i, "vozocrtx", kt, un               , jpi*jpj*jpk, idex )    ! now i-velocity
1049      CALL histwrite( id_i, "vomecrty", kt, vn               , jpi*jpj*jpk, idex )    ! now j-velocity
1050      CALL histwrite( id_i, "vovecrtz", kt, wn               , jpi*jpj*jpk, idex )    ! now k-velocity
1051      CALL histwrite( id_i, "sowaflup", kt, (emp-rnf )       , jpi*jpj    , idex )    ! freshwater budget
1052      CALL histwrite( id_i, "sohefldo", kt, qsr + qns        , jpi*jpj    , idex )    ! total heat flux
1053      CALL histwrite( id_i, "soshfldo", kt, qsr              , jpi*jpj    , idex )    ! solar heat flux
1054      CALL histwrite( id_i, "soicecov", kt, fr_i             , jpi*jpj    , idex )    ! ice fraction
1055      CALL histwrite( id_i, "sozotaux", kt, utau             , jpi*jpj    , idex )    ! i-wind stress
1056      CALL histwrite( id_i, "sometauy", kt, vtau             , jpi*jpj    , idex )    ! j-wind stress
1057
1058      ! 3. Close the file
1059      ! -----------------
1060      CALL histclo( id_i )
1061#if ! defined key_iomput && ! defined key_dimgout
1062      IF( ninist /= 1  ) THEN
1063         CALL histclo( nid_T )
1064         CALL histclo( nid_U )
1065         CALL histclo( nid_V )
1066         CALL histclo( nid_W )
1067      ENDIF
1068#endif
1069       
1070!     IF( nn_timing == 1 )   CALL timing_stop('dia_wri_state') ! not sure this works for routines not called in first timestep
1071      !
1072
1073   END SUBROUTINE dia_wri_state
1074   !!======================================================================
1075END MODULE diawri
Note: See TracBrowser for help on using the repository browser.