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

source: branches/UKMO/dev_r5518_MO_couple_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 6542

Last change on this file since 6542 was 6542, checked in by frrh, 7 years ago

Strip out svn keywords and properties.

File size: 56.3 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_heattr") .OR. iom_use("u_salttr") ) THEN
311         z3d(:,:,jpk) = 0.e0
312         DO jk = 1, jpkm1
313            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) * umask(:,:,jk)
314         END DO
315         CALL iom_put( "u_masstr", z3d )                  ! mass transport in i-direction
316      ENDIF
317     
318      IF( iom_use("u_heattr") ) THEN
319         z2d(:,:) = 0.e0 
320         DO jk = 1, jpkm1
321            DO jj = 2, jpjm1
322               DO ji = fs_2, fs_jpim1   ! vector opt.
323                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
324               END DO
325            END DO
326         END DO
327         CALL lbc_lnk( z2d, 'U', -1. )
328         CALL iom_put( "u_heattr", (0.5 * rcp) * z2d )    ! heat transport in i-direction
329      ENDIF
330
331      IF( iom_use("u_salttr") ) THEN
332         z2d(:,:) = 0.e0 
333         DO jk = 1, jpkm1
334            DO jj = 2, jpjm1
335               DO ji = fs_2, fs_jpim1   ! vector opt.
336                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
337               END DO
338            END DO
339         END DO
340         CALL lbc_lnk( z2d, 'U', -1. )
341         CALL iom_put( "u_salttr", 0.5 * z2d )            ! heat transport in i-direction
342      ENDIF
343
344     
345      IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
346         z3d(:,:,jpk) = 0.e0
347         DO jk = 1, jpkm1
348            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) * vmask(:,:,jk)
349         END DO
350         CALL iom_put( "v_masstr", z3d )                  ! mass transport in j-direction
351      ENDIF
352     
353      IF( iom_use("v_heattr") ) THEN
354         z2d(:,:) = 0.e0 
355         DO jk = 1, jpkm1
356            DO jj = 2, jpjm1
357               DO ji = fs_2, fs_jpim1   ! vector opt.
358                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
359               END DO
360            END DO
361         END DO
362         CALL lbc_lnk( z2d, 'V', -1. )
363         CALL iom_put( "v_heattr", (0.5 * rcp) * z2d )    !  heat transport in j-direction
364      ENDIF
365
366      IF( iom_use("v_salttr") ) THEN
367         z2d(:,:) = 0.e0 
368         DO jk = 1, jpkm1
369            DO jj = 2, jpjm1
370               DO ji = fs_2, fs_jpim1   ! vector opt.
371                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
372               END DO
373            END DO
374         END DO
375         CALL lbc_lnk( z2d, 'V', -1. )
376         CALL iom_put( "v_salttr", 0.5 * z2d )            !  heat transport in j-direction
377      ENDIF
378      !
379      CALL wrk_dealloc( jpi , jpj      , z2d )
380      CALL wrk_dealloc( jpi , jpj, jpk , z3d )
381      !
382      IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
383      !
384   END SUBROUTINE dia_wri
385
386#else
387   !!----------------------------------------------------------------------
388   !!   Default option                                  use IOIPSL  library
389   !!----------------------------------------------------------------------
390
391   SUBROUTINE dia_wri( kt )
392      !!---------------------------------------------------------------------
393      !!                  ***  ROUTINE dia_wri  ***
394      !!                   
395      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
396      !!      NETCDF format is used by default
397      !!
398      !! ** Method  :   At the beginning of the first time step (nit000),
399      !!      define all the NETCDF files and fields
400      !!      At each time step call histdef to compute the mean if ncessary
401      !!      Each nwrite time step, output the instantaneous or mean fields
402      !!----------------------------------------------------------------------
403      !!
404      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
405      !!
406      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
407      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
408      INTEGER  ::   inum = 11                                ! temporary logical unit
409      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
410      INTEGER  ::   ierr                                     ! error code return from allocation
411      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
412      INTEGER  ::   jn, ierror                               ! local integers
413      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt           ! local scalars
414      !!
415      REAL(wp), POINTER, DIMENSION(:,:)   :: zw2d       ! 2D workspace
416      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d       ! 3D workspace
417      !!----------------------------------------------------------------------
418      !
419      IF( nn_timing == 1 )   CALL timing_start('dia_wri')
420      !
421      CALL wrk_alloc( jpi , jpj      , zw2d )
422      IF ( ln_traldf_gdia .OR. lk_vvl )  call wrk_alloc( jpi , jpj , jpk  , zw3d )
423      !
424      ! Output the initial state and forcings
425      IF( ninist == 1 ) THEN                       
426         CALL dia_wri_state( 'output.init', kt )
427         ninist = 0
428      ENDIF
429      !
430      ! 0. Initialisation
431      ! -----------------
432
433      ! local variable for debugging
434      ll_print = .FALSE.
435      ll_print = ll_print .AND. lwp
436
437      ! Define frequency of output and means
438      zdt = rdt
439      IF( nacc == 1 ) zdt = rdtmin
440      IF( ln_mskland )   THEN   ;   clop = "only(x)"   ! put 1.e+20 on land (very expensive!!)
441      ELSE                      ;   clop = "x"         ! no use of the mask value (require less cpu time)
442      ENDIF
443#if defined key_diainstant
444      zsto = nwrite * zdt
445      clop = "inst("//TRIM(clop)//")"
446#else
447      zsto=zdt
448      clop = "ave("//TRIM(clop)//")"
449#endif
450      zout = nwrite * zdt
451      zmax = ( nitend - nit000 + 1 ) * zdt
452
453      ! Define indices of the horizontal output zoom and vertical limit storage
454      iimi = 1      ;      iima = jpi
455      ijmi = 1      ;      ijma = jpj
456      ipk = jpk
457
458      ! define time axis
459      it = kt
460      itmod = kt - nit000 + 1
461
462
463      ! 1. Define NETCDF files and fields at beginning of first time step
464      ! -----------------------------------------------------------------
465
466      IF( kt == nit000 ) THEN
467
468         ! Define the NETCDF files (one per grid)
469
470         ! Compute julian date from starting date of the run
471         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
472         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
473         IF(lwp)WRITE(numout,*)
474         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
475            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
476         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
477                                 ' limit storage in depth = ', ipk
478
479         ! WRITE root name in date.file for use by postpro
480         IF(lwp) THEN
481            CALL dia_nam( clhstnam, nwrite,' ' )
482            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
483            WRITE(inum,*) clhstnam
484            CLOSE(inum)
485         ENDIF
486
487         ! Define the T grid FILE ( nid_T )
488
489         CALL dia_nam( clhstnam, nwrite, 'grid_T' )
490         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
491         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
492            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
493            &          nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
494         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
495            &           "m", ipk, gdept_1d, nz_T, "down" )
496         !                                                            ! Index of ocean points
497         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
498         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
499         !
500         IF( ln_icebergs ) THEN
501            !
502            !! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
503            !! that routine is called from nemogcm, so do it here immediately before its needed
504            ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
505            IF( lk_mpp )   CALL mpp_sum( ierror )
506            IF( ierror /= 0 ) THEN
507               CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
508               RETURN
509            ENDIF
510            !
511            !! iceberg vertical coordinate is class number
512            CALL histvert( nid_T, "class", "Iceberg class",      &  ! Vertical grid: class
513               &           "number", nclasses, class_num, nb_T )
514            !
515            !! each class just needs the surface index pattern
516            ndim_bT = 3
517            DO jn = 1,nclasses
518               ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
519            ENDDO
520            !
521         ENDIF
522
523         ! Define the U grid FILE ( nid_U )
524
525         CALL dia_nam( clhstnam, nwrite, 'grid_U' )
526         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
527         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
528            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
529            &          nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
530         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
531            &           "m", ipk, gdept_1d, nz_U, "down" )
532         !                                                            ! Index of ocean points
533         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
534         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
535
536         ! Define the V grid FILE ( nid_V )
537
538         CALL dia_nam( clhstnam, nwrite, 'grid_V' )                   ! filename
539         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
540         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
541            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
542            &          nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
543         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
544            &          "m", ipk, gdept_1d, nz_V, "down" )
545         !                                                            ! Index of ocean points
546         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
547         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
548
549         ! Define the W grid FILE ( nid_W )
550
551         CALL dia_nam( clhstnam, nwrite, 'grid_W' )                   ! filename
552         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
553         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
554            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
555            &          nit000-1, zjulian, zdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
556         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
557            &          "m", ipk, gdepw_1d, nz_W, "down" )
558
559
560         ! Declare all the output fields as NETCDF variables
561
562         !                                                                                      !!! nid_T : 3D
563         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
564            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
565         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
566            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
567         IF(  lk_vvl  ) THEN
568            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n
569            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
570            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t_n
571            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
572            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t_n
573            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
574         ENDIF
575         !                                                                                      !!! nid_T : 2D
576         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
577            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
578         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
579            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
580         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
581            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
582         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
583            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
584         CALL histdef( nid_T, "sorunoff", "River runoffs"                      , "Kg/m2/s",   &  ! runoffs
585            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
586         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
587            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
588         IF(  .NOT. lk_vvl  ) THEN
589            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem)
590            &                                                                  , "KgC/m2/s",  &  ! sosst_cd
591            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
592            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * tsn(:,:,1,jp_sal)
593            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
594            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
595         ENDIF
596         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
597            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
598         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
599            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
600         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
601            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
602         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
603            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
604         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
605            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
606         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
607            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
608!
609         IF( ln_icebergs ) THEN
610            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
611               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
612            CALL histdef( nid_T, "calving_heat"        , "calving heat flux"                        , "XXXX"   , &
613               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
614            CALL histdef( nid_T, "berg_floating_melt"  , "Melt rate of icebergs + bits"             , "kg/m2/s", &
615               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
616            CALL histdef( nid_T, "berg_stored_ice"     , "Accumulated ice mass by class"            , "kg"     , &
617               &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
618            IF( ln_bergdia ) THEN
619               CALL histdef( nid_T, "berg_melt"           , "Melt rate of icebergs"                    , "kg/m2/s", &
620                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
621               CALL histdef( nid_T, "berg_buoy_melt"      , "Buoyancy component of iceberg melt rate"  , "kg/m2/s", &
622                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
623               CALL histdef( nid_T, "berg_eros_melt"      , "Erosion component of iceberg melt rate"   , "kg/m2/s", &
624                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
625               CALL histdef( nid_T, "berg_conv_melt"      , "Convective component of iceberg melt rate", "kg/m2/s", &
626                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
627               CALL histdef( nid_T, "berg_virtual_area"   , "Virtual coverage by icebergs"             , "m2"     , &
628                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
629               CALL histdef( nid_T, "bits_src"           , "Mass source of bergy bits"                , "kg/m2/s", &
630                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
631               CALL histdef( nid_T, "bits_melt"          , "Melt rate of bergy bits"                  , "kg/m2/s", &
632                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
633               CALL histdef( nid_T, "bits_mass"          , "Bergy bit density field"                  , "kg/m2"  , &
634                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
635               CALL histdef( nid_T, "berg_mass"           , "Iceberg density field"                    , "kg/m2"  , &
636                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
637               CALL histdef( nid_T, "berg_real_calving"   , "Calving into iceberg class"               , "kg/s"   , &
638                  &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
639            ENDIF
640         ENDIF
641
642         IF( .NOT. ln_cpl ) THEN
643            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
644               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
645            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
646               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
647            CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
648               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
649         ENDIF
650
651         IF( ln_cpl .AND. nn_ice <= 1 ) THEN
652            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
653               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
654            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
655               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
656            CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn
657               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
658         ENDIF
659         
660         clmx ="l_max(only(x))"    ! max index on a period
661         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
662            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
663#if defined key_diahth
664         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
665            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
666         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
667            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
668         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
669            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
670         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "W"      ,   & ! htc3
671            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
672#endif
673
674         IF( ln_cpl .AND. nn_ice == 2 ) THEN
675            CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature"            , "K"      ,   &  ! tn_ice
676               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
677            CALL histdef( nid_T,"soicealb" , "Ice Albedo"                         , "[0,1]"  ,   &  ! alb_ice
678               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
679         ENDIF
680
681         CALL histend( nid_T, snc4chunks=snc4set )
682
683         !                                                                                      !!! nid_U : 3D
684         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
685            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
686         IF( ln_traldf_gdia ) THEN
687            CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
688                 &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
689         ELSE
690#if defined key_diaeiv
691            CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
692            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
693#endif
694         END IF
695         !                                                                                      !!! nid_U : 2D
696         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
697            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
698
699         CALL histend( nid_U, snc4chunks=snc4set )
700
701         !                                                                                      !!! nid_V : 3D
702         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
703            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
704         IF( ln_traldf_gdia ) THEN
705            CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
706                 &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
707         ELSE 
708#if defined key_diaeiv
709            CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
710            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
711#endif
712         END IF
713         !                                                                                      !!! nid_V : 2D
714         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
715            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
716
717         CALL histend( nid_V, snc4chunks=snc4set )
718
719         !                                                                                      !!! nid_W : 3D
720         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn
721            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
722         IF( ln_traldf_gdia ) THEN
723            CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
724                 &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
725         ELSE
726#if defined key_diaeiv
727            CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
728                 &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
729#endif
730         END IF
731         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
732            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
733         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avmu
734            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
735
736         IF( lk_zdfddm ) THEN
737            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
738               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
739         ENDIF
740         !                                                                                      !!! nid_W : 2D
741#if defined key_traldf_c2d
742         CALL histdef( nid_W, "soleahtw", "lateral eddy diffusivity"           , "m2/s"   ,   &  ! ahtw
743            &          jpi, jpj, nh_W, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
744# if defined key_traldf_eiv 
745            CALL histdef( nid_W, "soleaeiw", "eddy induced vel. coeff. at w-point", "m2/s",   &  ! aeiw
746               &       jpi, jpj, nh_W, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
747# endif
748#endif
749
750         CALL histend( nid_W, snc4chunks=snc4set )
751
752         IF(lwp) WRITE(numout,*)
753         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
754         IF(ll_print) CALL FLUSH(numout )
755
756      ENDIF
757
758      ! 2. Start writing data
759      ! ---------------------
760
761      ! ndex(1) est utilise ssi l'avant dernier argument est different de
762      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
763      ! donne le nombre d'elements, et ndex la liste des indices a sortir
764
765      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
766         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
767         WRITE(numout,*) '~~~~~~ '
768      ENDIF
769
770      IF( lk_vvl ) THEN
771         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
772         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
773         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
774         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
775      ELSE
776         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature
777         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity
778         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature
779         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
780      ENDIF
781      IF( lk_vvl ) THEN
782         zw3d(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
783         CALL histwrite( nid_T, "vovvle3t", it, fse3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
784         CALL histwrite( nid_T, "vovvldep", it, fsdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
785         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
786      ENDIF
787      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
788      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
789      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
790      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
791                                                                                  ! (includes virtual salt flux beneath ice
792                                                                                  ! in linear free surface case)
793      IF( .NOT. lk_vvl ) THEN
794         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
795         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
796         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
797         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
798      ENDIF
799      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
800      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
801      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
802      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
803      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
804      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
805!
806      IF( ln_icebergs ) THEN
807         !
808         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT ) 
809         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
810         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT ) 
811         !
812         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
813         !
814         IF( ln_bergdia ) THEN
815            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   ) 
816            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   ) 
817            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   ) 
818            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   ) 
819            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   ) 
820            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   ) 
821            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   ) 
822            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   ) 
823            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   ) 
824            !
825            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
826         ENDIF
827      ENDIF
828
829      IF( .NOT. ln_cpl ) THEN
830         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
831         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
832         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
833         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
834      ENDIF
835      IF( ln_cpl .AND. nn_ice <= 1 ) THEN
836         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
837         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
838         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
839         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
840      ENDIF
841!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
842!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
843
844#if defined key_diahth
845      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
846      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
847      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
848      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
849#endif
850
851      IF( ln_cpl .AND. nn_ice == 2 ) THEN
852         CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature
853         CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo
854      ENDIF
855
856      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
857      IF( ln_traldf_gdia ) THEN
858         IF (.not. ALLOCATED(psix_eiv))THEN
859            ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr )
860            IF( lk_mpp   )   CALL mpp_sum ( ierr )
861            IF( ierr > 0 )   CALL ctl_stop('STOP', 'diawri: unable to allocate psi{x,y}_eiv')
862            psix_eiv(:,:,:) = 0.0_wp
863            psiy_eiv(:,:,:) = 0.0_wp
864         ENDIF
865         DO jk=1,jpkm1
866            zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk)  ! u_eiv = -dpsix/dz
867         END DO
868         zw3d(:,:,jpk) = 0._wp
869         CALL histwrite( nid_U, "vozoeivu", it, zw3d, ndim_U , ndex_U )           ! i-eiv current
870      ELSE
871#if defined key_diaeiv
872         CALL histwrite( nid_U, "vozoeivu", it, u_eiv, ndim_U , ndex_U )          ! i-eiv current
873#endif
874      ENDIF
875      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
876
877      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
878      IF( ln_traldf_gdia ) THEN
879         DO jk=1,jpk-1
880            zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk)  ! v_eiv = -dpsiy/dz
881         END DO
882         zw3d(:,:,jpk) = 0._wp
883         CALL histwrite( nid_V, "vomeeivv", it, zw3d, ndim_V , ndex_V )           ! j-eiv current
884      ELSE
885#if defined key_diaeiv
886         CALL histwrite( nid_V, "vomeeivv", it, v_eiv, ndim_V , ndex_V )          ! j-eiv current
887#endif
888      ENDIF
889      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
890
891      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
892      IF( ln_traldf_gdia ) THEN
893         DO jk=1,jpk-1
894            DO jj = 2, jpjm1
895               DO ji = fs_2, fs_jpim1  ! vector opt.
896                  zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2v(ji,jj) + &
897                       &    (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1u(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx
898               END DO
899            END DO
900         END DO
901         zw3d(:,:,jpk) = 0._wp
902         CALL histwrite( nid_W, "voveeivw", it, zw3d          , ndim_T, ndex_T )    ! vert. eiv current
903      ELSE
904#   if defined key_diaeiv
905         CALL histwrite( nid_W, "voveeivw", it, w_eiv          , ndim_T, ndex_T )    ! vert. eiv current
906#   endif
907      ENDIF
908      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
909      CALL histwrite( nid_W, "votkeavm", it, avmu           , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
910      IF( lk_zdfddm ) THEN
911         CALL histwrite( nid_W, "voddmavs", it, fsavs(:,:,:), ndim_T, ndex_T )    ! S vert. eddy diff. coef.
912      ENDIF
913#if defined key_traldf_c2d
914      CALL histwrite( nid_W, "soleahtw", it, ahtw          , ndim_hT, ndex_hT )   ! lateral eddy diff. coef.
915# if defined key_traldf_eiv
916         CALL histwrite( nid_W, "soleaeiw", it, aeiw       , ndim_hT, ndex_hT )   ! EIV coefficient at w-point
917# endif
918#endif
919
920      ! 3. Close all files
921      ! ---------------------------------------
922      IF( kt == nitend ) THEN
923         CALL histclo( nid_T )
924         CALL histclo( nid_U )
925         CALL histclo( nid_V )
926         CALL histclo( nid_W )
927      ENDIF
928      !
929      CALL wrk_dealloc( jpi , jpj      , zw2d )
930      IF ( ln_traldf_gdia .OR. lk_vvl )  call wrk_dealloc( jpi , jpj , jpk  , zw3d )
931      !
932      IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
933      !
934   END SUBROUTINE dia_wri
935# endif
936
937#endif
938
939   SUBROUTINE dia_wri_state( cdfile_name, kt )
940      !!---------------------------------------------------------------------
941      !!                 ***  ROUTINE dia_wri_state  ***
942      !!       
943      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
944      !!      the instantaneous ocean state and forcing fields.
945      !!        Used to find errors in the initial state or save the last
946      !!      ocean state in case of abnormal end of a simulation
947      !!
948      !! ** Method  :   NetCDF files using ioipsl
949      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
950      !!      File 'output.abort.nc' is created in case of abnormal job end
951      !!----------------------------------------------------------------------
952      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
953      INTEGER           , INTENT( in ) ::   kt               ! ocean time-step index
954      !!
955      CHARACTER (len=32) :: clname
956      CHARACTER (len=40) :: clop
957      INTEGER  ::   id_i , nz_i, nh_i       
958      INTEGER, DIMENSION(1) ::   idex             ! local workspace
959      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt
960      !!----------------------------------------------------------------------
961      !
962!     IF( nn_timing == 1 )   CALL timing_start('dia_wri_state') ! not sure this works for routines not called in first timestep
963
964      ! 0. Initialisation
965      ! -----------------
966
967      ! Define name, frequency of output and means
968      clname = cdfile_name
969      IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
970      zdt  = rdt
971      zsto = rdt
972      clop = "inst(x)"           ! no use of the mask value (require less cpu time)
973      zout = rdt
974      zmax = ( nitend - nit000 + 1 ) * zdt
975
976      IF(lwp) WRITE(numout,*)
977      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
978      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
979      IF(lwp) WRITE(numout,*) '                and named :', clname, '.nc'
980
981
982      ! 1. Define NETCDF files and fields at beginning of first time step
983      ! -----------------------------------------------------------------
984
985      ! Compute julian date from starting date of the run
986      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )         ! time axis
987      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
988      CALL histbeg( clname, jpi, glamt, jpj, gphit,   &
989          1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
990      CALL histvert( id_i, "deptht", "Vertical T levels",   &    ! Vertical grid : gdept
991          "m", jpk, gdept_1d, nz_i, "down")
992
993      ! Declare all the output fields as NetCDF variables
994
995      CALL histdef( id_i, "vosaline", "Salinity"              , "PSU"    ,   &   ! salinity
996         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
997      CALL histdef( id_i, "votemper", "Temperature"           , "C"      ,   &   ! temperature
998         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
999      CALL histdef( id_i, "sossheig", "Sea Surface Height"    , "m"      ,   &  ! ssh
1000         &          jpi, jpj, nh_i, 1  , 1, 1  , nz_i, 32, clop, zsto, zout )
1001      CALL histdef( id_i, "vozocrtx", "Zonal Current"         , "m/s"    ,   &   ! zonal current
1002         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1003      CALL histdef( id_i, "vomecrty", "Meridional Current"    , "m/s"    ,   &   ! meridonal current
1004         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
1005      CALL histdef( id_i, "vovecrtz", "Vertical Velocity"     , "m/s"    ,   &   ! vertical current
1006         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
1007      CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S",   &   ! net freshwater
1008         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1009      CALL histdef( id_i, "sohefldo", "Net Downward Heat Flux", "W/m2"   ,   &   ! net heat flux
1010         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1011      CALL histdef( id_i, "soshfldo", "Shortwave Radiation"   , "W/m2"   ,   &   ! solar flux
1012         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1013      CALL histdef( id_i, "soicecov", "Ice fraction"          , "[0,1]"  ,   &   ! fr_i
1014         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1015      CALL histdef( id_i, "sozotaux", "Zonal Wind Stress"     , "N/m2"   ,   &   ! i-wind stress
1016         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1017      CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2"   ,   &   ! j-wind stress
1018         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1019      IF( lk_vvl ) THEN
1020         CALL histdef( id_i, "vovvldep", "T point depth"         , "m"      ,   &   ! t-point depth
1021            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1022      END IF
1023
1024#if defined key_lim2
1025      CALL lim_wri_state_2( kt, id_i, nh_i )
1026#elif defined key_lim3
1027      CALL lim_wri_state( kt, id_i, nh_i )
1028#else
1029      CALL histend( id_i, snc4chunks=snc4set )
1030#endif
1031
1032      ! 2. Start writing data
1033      ! ---------------------
1034      ! idex(1) est utilise ssi l'avant dernier argument est diffferent de
1035      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
1036      ! donne le nombre d'elements, et idex la liste des indices a sortir
1037      idex(1) = 1   ! init to avoid compil warning
1038
1039      ! Write all fields on T grid
1040      CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature
1041      CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity
1042      CALL histwrite( id_i, "sossheig", kt, sshn             , jpi*jpj    , idex )    ! sea surface height
1043      CALL histwrite( id_i, "vozocrtx", kt, un               , jpi*jpj*jpk, idex )    ! now i-velocity
1044      CALL histwrite( id_i, "vomecrty", kt, vn               , jpi*jpj*jpk, idex )    ! now j-velocity
1045      CALL histwrite( id_i, "vovecrtz", kt, wn               , jpi*jpj*jpk, idex )    ! now k-velocity
1046      CALL histwrite( id_i, "sowaflup", kt, (emp-rnf )       , jpi*jpj    , idex )    ! freshwater budget
1047      CALL histwrite( id_i, "sohefldo", kt, qsr + qns        , jpi*jpj    , idex )    ! total heat flux
1048      CALL histwrite( id_i, "soshfldo", kt, qsr              , jpi*jpj    , idex )    ! solar heat flux
1049      CALL histwrite( id_i, "soicecov", kt, fr_i             , jpi*jpj    , idex )    ! ice fraction
1050      CALL histwrite( id_i, "sozotaux", kt, utau             , jpi*jpj    , idex )    ! i-wind stress
1051      CALL histwrite( id_i, "sometauy", kt, vtau             , jpi*jpj    , idex )    ! j-wind stress
1052
1053      ! 3. Close the file
1054      ! -----------------
1055      CALL histclo( id_i )
1056#if ! defined key_iomput && ! defined key_dimgout
1057      IF( ninist /= 1  ) THEN
1058         CALL histclo( nid_T )
1059         CALL histclo( nid_U )
1060         CALL histclo( nid_V )
1061         CALL histclo( nid_W )
1062      ENDIF
1063#endif
1064       
1065!     IF( nn_timing == 1 )   CALL timing_stop('dia_wri_state') ! not sure this works for routines not called in first timestep
1066      !
1067
1068   END SUBROUTINE dia_wri_state
1069   !!======================================================================
1070END MODULE diawri
Note: See TracBrowser for help on using the repository browser.