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

source: branches/UKMO/2015_CO6_CO5_shelfdiagnostic/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 5666

Last change on this file since 5666 was 5666, checked in by deazer, 9 years ago

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