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

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

Last change on this file since 5836 was 5836, checked in by cetlod, 9 years ago

merge the simplification branch onto the trunk, see ticket #1612
  • Property svn:keywords set to Id
File size: 51.9 KB
Line 
1MODULE diawri
2   !!======================================================================
3   !!                     ***  MODULE  diawri  ***
4   !! Ocean diagnostics :  write ocean output files
5   !!=====================================================================
6   !! History :  OPA  ! 1991-03  (M.-A. Foujols)  Original code
7   !!            4.0  ! 1991-11  (G. Madec)
8   !!                 ! 1992-06  (M. Imbard)  correction restart file
9   !!                 ! 1992-07  (M. Imbard)  split into diawri and rstwri
10   !!                 ! 1993-03  (M. Imbard)  suppress writibm
11   !!                 ! 1998-01  (C. Levy)  NETCDF format using ioipsl INTERFACE
12   !!                 ! 1999-02  (E. Guilyardi)  name of netCDF files + variables
13   !!            8.2  ! 2000-06  (M. Imbard)  Original code (diabort.F)
14   !!   NEMO     1.0  ! 2002-06  (A.Bozec, E. Durand)  Original code (diainit.F)
15   !!             -   ! 2002-09  (G. Madec)  F90: Free form and module
16   !!             -   ! 2002-12  (G. Madec)  merge of diabort and diainit, F90
17   !!                 ! 2005-11  (V. Garnier) Surface pressure gradient organization
18   !!            3.2  ! 2008-11  (B. Lemaire) creation from old diawri
19   !!            3.7  ! 2014-01  (G. Madec) remove eddy induced velocity from no-IOM output
20   !!                 !                     change name of output variables in dia_wri_state
21   !!----------------------------------------------------------------------
22
23   !!----------------------------------------------------------------------
24   !!   dia_wri       : create the standart output files
25   !!   dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
26   !!----------------------------------------------------------------------
27   USE oce             ! ocean dynamics and tracers
28   USE dom_oce         ! ocean space and time domain
29   USE dynadv, ONLY: ln_dynadv_vec
30   USE zdf_oce         ! ocean vertical physics
31   USE ldftra          ! lateral physics: eddy diffusivity coef.
32   USE 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 * e1e2t(:,:)
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 ) * r1_e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,1,jp_tem) ) * r1_e1u(ji-1,jj)
251               zztmpy = ( tsn(ji,jj+1,1,jp_tem) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - tsn(ji  ,jj-1,1,jp_tem) ) * r1_e2v(ji,jj-1)
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      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
404      !
405      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
406      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
407      INTEGER  ::   inum = 11                                ! temporary logical unit
408      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
409      INTEGER  ::   ierr                                     ! error code return from allocation
410      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
411      INTEGER  ::   jn, ierror                               ! local integers
412      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt           ! local scalars
413      !
414      REAL(wp), POINTER, DIMENSION(:,:)   :: zw2d       ! 2D workspace
415      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d       ! 3D workspace
416      !!----------------------------------------------------------------------
417      !
418      IF( nn_timing == 1 )   CALL timing_start('dia_wri')
419      !
420                     CALL wrk_alloc( jpi,jpj      , zw2d )
421      IF( lk_vvl )   CALL wrk_alloc( jpi,jpj,jpk  , zw3d )
422      !
423      ! Output the initial state and forcings
424      IF( ninist == 1 ) THEN                       
425         CALL dia_wri_state( 'output.init', kt )
426         ninist = 0
427      ENDIF
428      !
429      ! 0. Initialisation
430      ! -----------------
431
432      ! local variable for debugging
433      ll_print = .FALSE.
434      ll_print = ll_print .AND. lwp
435
436      ! Define frequency of output and means
437      zdt = rdt
438      IF( nacc == 1 ) zdt = rdtmin
439      clop = "x"         ! no use of the mask value (require less cpu time and otherwise the model crashes)
440#if defined key_diainstant
441      zsto = nwrite * zdt
442      clop = "inst("//TRIM(clop)//")"
443#else
444      zsto=zdt
445      clop = "ave("//TRIM(clop)//")"
446#endif
447      zout = nwrite * zdt
448      zmax = ( nitend - nit000 + 1 ) * zdt
449
450      ! Define indices of the horizontal output zoom and vertical limit storage
451      iimi = 1      ;      iima = jpi
452      ijmi = 1      ;      ijma = jpj
453      ipk = jpk
454
455      ! define time axis
456      it = kt
457      itmod = kt - nit000 + 1
458
459
460      ! 1. Define NETCDF files and fields at beginning of first time step
461      ! -----------------------------------------------------------------
462
463      IF( kt == nit000 ) THEN
464
465         ! Define the NETCDF files (one per grid)
466
467         ! Compute julian date from starting date of the run
468         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
469         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
470         IF(lwp)WRITE(numout,*)
471         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
472            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
473         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
474                                 ' limit storage in depth = ', ipk
475
476         ! WRITE root name in date.file for use by postpro
477         IF(lwp) THEN
478            CALL dia_nam( clhstnam, nwrite,' ' )
479            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
480            WRITE(inum,*) clhstnam
481            CLOSE(inum)
482         ENDIF
483
484         ! Define the T grid FILE ( nid_T )
485
486         CALL dia_nam( clhstnam, nwrite, 'grid_T' )
487         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
488         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
489            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
490            &          nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
491         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
492            &           "m", ipk, gdept_1d, nz_T, "down" )
493         !                                                            ! Index of ocean points
494         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
495         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
496         !
497         IF( ln_icebergs ) THEN
498            !
499            !! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
500            !! that routine is called from nemogcm, so do it here immediately before its needed
501            ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
502            IF( lk_mpp )   CALL mpp_sum( ierror )
503            IF( ierror /= 0 ) THEN
504               CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
505               RETURN
506            ENDIF
507            !
508            !! iceberg vertical coordinate is class number
509            CALL histvert( nid_T, "class", "Iceberg class",      &  ! Vertical grid: class
510               &           "number", nclasses, class_num, nb_T )
511            !
512            !! each class just needs the surface index pattern
513            ndim_bT = 3
514            DO jn = 1,nclasses
515               ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
516            ENDDO
517            !
518         ENDIF
519
520         ! Define the U grid FILE ( nid_U )
521
522         CALL dia_nam( clhstnam, nwrite, 'grid_U' )
523         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
524         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
525            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
526            &          nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
527         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
528            &           "m", ipk, gdept_1d, nz_U, "down" )
529         !                                                            ! Index of ocean points
530         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
531         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
532
533         ! Define the V grid FILE ( nid_V )
534
535         CALL dia_nam( clhstnam, nwrite, 'grid_V' )                   ! filename
536         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
537         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
538            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
539            &          nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
540         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
541            &          "m", ipk, gdept_1d, nz_V, "down" )
542         !                                                            ! Index of ocean points
543         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
544         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
545
546         ! Define the W grid FILE ( nid_W )
547
548         CALL dia_nam( clhstnam, nwrite, 'grid_W' )                   ! filename
549         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
550         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
551            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
552            &          nit000-1, zjulian, zdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
553         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
554            &          "m", ipk, gdepw_1d, nz_W, "down" )
555
556
557         ! Declare all the output fields as NETCDF variables
558
559         !                                                                                      !!! nid_T : 3D
560         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
561            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
562         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
563            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
564         IF(  lk_vvl  ) THEN
565            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n
566            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
567            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t_n
568            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
569            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t_n
570            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
571         ENDIF
572         !                                                                                      !!! nid_T : 2D
573         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
574            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
575         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
576            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
577         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
578            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
579         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
580            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
581         CALL histdef( nid_T, "sorunoff", "River runoffs"                      , "Kg/m2/s",   &  ! runoffs
582            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
583         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
584            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
585         IF(  .NOT. lk_vvl  ) THEN
586            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem)
587            &                                                                  , "KgC/m2/s",  &  ! sosst_cd
588            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
589            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * tsn(:,:,1,jp_sal)
590            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
591            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
592         ENDIF
593         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
594            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
595         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
596            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
597         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
598            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
599         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
600            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
601         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
602            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
603         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
604            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
605!
606         IF( ln_icebergs ) THEN
607            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
608               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
609            CALL histdef( nid_T, "calving_heat"        , "calving heat flux"                        , "XXXX"   , &
610               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
611            CALL histdef( nid_T, "berg_floating_melt"  , "Melt rate of icebergs + bits"             , "kg/m2/s", &
612               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
613            CALL histdef( nid_T, "berg_stored_ice"     , "Accumulated ice mass by class"            , "kg"     , &
614               &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
615            IF( ln_bergdia ) THEN
616               CALL histdef( nid_T, "berg_melt"           , "Melt rate of icebergs"                    , "kg/m2/s", &
617                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
618               CALL histdef( nid_T, "berg_buoy_melt"      , "Buoyancy component of iceberg melt rate"  , "kg/m2/s", &
619                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
620               CALL histdef( nid_T, "berg_eros_melt"      , "Erosion component of iceberg melt rate"   , "kg/m2/s", &
621                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
622               CALL histdef( nid_T, "berg_conv_melt"      , "Convective component of iceberg melt rate", "kg/m2/s", &
623                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
624               CALL histdef( nid_T, "berg_virtual_area"   , "Virtual coverage by icebergs"             , "m2"     , &
625                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
626               CALL histdef( nid_T, "bits_src"           , "Mass source of bergy bits"                , "kg/m2/s", &
627                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
628               CALL histdef( nid_T, "bits_melt"          , "Melt rate of bergy bits"                  , "kg/m2/s", &
629                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
630               CALL histdef( nid_T, "bits_mass"          , "Bergy bit density field"                  , "kg/m2"  , &
631                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
632               CALL histdef( nid_T, "berg_mass"           , "Iceberg density field"                    , "kg/m2"  , &
633                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
634               CALL histdef( nid_T, "berg_real_calving"   , "Calving into iceberg class"               , "kg/s"   , &
635                  &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
636            ENDIF
637         ENDIF
638
639         IF( .NOT. ln_cpl ) THEN
640            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
641               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
642            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
643               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
644            CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
645               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
646         ENDIF
647
648         IF( ln_cpl .AND. nn_ice <= 1 ) THEN
649            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
650               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
651            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
652               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
653            CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn
654               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
655         ENDIF
656         
657         clmx ="l_max(only(x))"    ! max index on a period
658!         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
659!            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
660#if defined key_diahth
661         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
662            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
663         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
664            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
665         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
666            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
667         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "W"      ,   & ! htc3
668            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
669#endif
670
671         IF( ln_cpl .AND. nn_ice == 2 ) THEN
672            CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature"            , "K"      ,   &  ! tn_ice
673               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
674            CALL histdef( nid_T,"soicealb" , "Ice Albedo"                         , "[0,1]"  ,   &  ! alb_ice
675               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
676         ENDIF
677
678         CALL histend( nid_T, snc4chunks=snc4set )
679
680         !                                                                                      !!! nid_U : 3D
681         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
682            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
683         !                                                                                      !!! nid_U : 2D
684         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
685            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
686
687         CALL histend( nid_U, snc4chunks=snc4set )
688
689         !                                                                                      !!! nid_V : 3D
690         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
691            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
692         !                                                                                      !!! nid_V : 2D
693         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
694            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
695
696         CALL histend( nid_V, snc4chunks=snc4set )
697
698         !                                                                                      !!! nid_W : 3D
699         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn
700            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
701         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
702            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
703         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avmu
704            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
705
706         IF( lk_zdfddm ) THEN
707            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
708               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
709         ENDIF
710         !                                                                                      !!! nid_W : 2D
711         CALL histend( nid_W, snc4chunks=snc4set )
712
713         IF(lwp) WRITE(numout,*)
714         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
715         IF(ll_print) CALL FLUSH(numout )
716
717      ENDIF
718
719      ! 2. Start writing data
720      ! ---------------------
721
722      ! ndex(1) est utilise ssi l'avant dernier argument est different de
723      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
724      ! donne le nombre d'elements, et ndex la liste des indices a sortir
725
726      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
727         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
728         WRITE(numout,*) '~~~~~~ '
729      ENDIF
730
731      IF( lk_vvl ) THEN
732         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
733         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
734         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
735         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
736      ELSE
737         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature
738         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity
739         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature
740         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
741      ENDIF
742      IF( lk_vvl ) THEN
743         zw3d(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
744         CALL histwrite( nid_T, "vovvle3t", it, fse3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
745         CALL histwrite( nid_T, "vovvldep", it, fsdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
746         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
747      ENDIF
748      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
749      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
750      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
751      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
752                                                                                  ! (includes virtual salt flux beneath ice
753                                                                                  ! in linear free surface case)
754      IF( .NOT. lk_vvl ) THEN
755         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
756         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
757         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
758         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
759      ENDIF
760      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
761      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
762      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
763      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
764      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
765      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
766!
767      IF( ln_icebergs ) THEN
768         !
769         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT ) 
770         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
771         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT ) 
772         !
773         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
774         !
775         IF( ln_bergdia ) THEN
776            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   ) 
777            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   ) 
778            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   ) 
779            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   ) 
780            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   ) 
781            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   ) 
782            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   ) 
783            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   ) 
784            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   ) 
785            !
786            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
787         ENDIF
788      ENDIF
789
790      IF( .NOT. ln_cpl ) THEN
791         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
792         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
793         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
794         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
795      ENDIF
796      IF( ln_cpl .AND. nn_ice <= 1 ) THEN
797         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
798         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
799         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
800         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
801      ENDIF
802!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
803!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
804
805#if defined key_diahth
806      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
807      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
808      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
809      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
810#endif
811
812      IF( ln_cpl .AND. nn_ice == 2 ) THEN
813         CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature
814         CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo
815      ENDIF
816
817      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
818      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
819
820      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
821      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
822
823      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
824      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
825      CALL histwrite( nid_W, "votkeavm", it, avmu           , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
826      IF( lk_zdfddm ) THEN
827         CALL histwrite( nid_W, "voddmavs", it, fsavs(:,:,:), ndim_T, ndex_T )    ! S vert. eddy diff. coef.
828      ENDIF
829
830      ! 3. Close all files
831      ! ---------------------------------------
832      IF( kt == nitend ) THEN
833         CALL histclo( nid_T )
834         CALL histclo( nid_U )
835         CALL histclo( nid_V )
836         CALL histclo( nid_W )
837      ENDIF
838      !
839                     CALL wrk_dealloc( jpi , jpj        , zw2d )
840      IF( lk_vvl )   CALL wrk_dealloc( jpi , jpj , jpk  , zw3d )
841      !
842      IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
843      !
844   END SUBROUTINE dia_wri
845# endif
846
847#endif
848
849   SUBROUTINE dia_wri_state( cdfile_name, kt )
850      !!---------------------------------------------------------------------
851      !!                 ***  ROUTINE dia_wri_state  ***
852      !!       
853      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
854      !!      the instantaneous ocean state and forcing fields.
855      !!        Used to find errors in the initial state or save the last
856      !!      ocean state in case of abnormal end of a simulation
857      !!
858      !! ** Method  :   NetCDF files using ioipsl
859      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
860      !!      File 'output.abort.nc' is created in case of abnormal job end
861      !!----------------------------------------------------------------------
862      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
863      INTEGER           , INTENT( in ) ::   kt               ! ocean time-step index
864      !!
865      CHARACTER (len=32) :: clname
866      CHARACTER (len=40) :: clop
867      INTEGER  ::   id_i , nz_i, nh_i       
868      INTEGER, DIMENSION(1) ::   idex             ! local workspace
869      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt
870      !!----------------------------------------------------------------------
871      !
872      ! 0. Initialisation
873      ! -----------------
874
875      ! Define name, frequency of output and means
876      clname = cdfile_name
877      IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
878      zdt  = rdt
879      zsto = rdt
880      clop = "inst(x)"           ! no use of the mask value (require less cpu time)
881      zout = rdt
882      zmax = ( nitend - nit000 + 1 ) * zdt
883
884      IF(lwp) WRITE(numout,*)
885      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
886      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
887      IF(lwp) WRITE(numout,*) '                and named :', clname, '.nc'
888
889
890      ! 1. Define NETCDF files and fields at beginning of first time step
891      ! -----------------------------------------------------------------
892
893      ! Compute julian date from starting date of the run
894      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )         ! time axis
895      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
896      CALL histbeg( clname, jpi, glamt, jpj, gphit,   &
897          1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
898      CALL histvert( id_i, "deptht", "Vertical T levels",   &    ! Vertical grid : gdept
899          "m", jpk, gdept_1d, nz_i, "down")
900
901      ! Declare all the output fields as NetCDF variables
902
903      CALL histdef( id_i, "vosaline", "Salinity"              , "PSU"    ,   &   ! salinity
904         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
905      CALL histdef( id_i, "votemper", "Temperature"           , "C"      ,   &   ! temperature
906         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
907      CALL histdef( id_i, "sossheig", "Sea Surface Height"    , "m"      ,   &  ! ssh
908         &          jpi, jpj, nh_i, 1  , 1, 1  , nz_i, 32, clop, zsto, zout )
909      CALL histdef( id_i, "vozocrtx", "Zonal Current"         , "m/s"    ,   &   ! zonal current
910         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
911      CALL histdef( id_i, "vomecrty", "Meridional Current"    , "m/s"    ,   &   ! meridonal current
912         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
913      CALL histdef( id_i, "vovecrtz", "Vertical Velocity"     , "m/s"    ,   &   ! vertical current
914         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
915      CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S",   &   ! net freshwater
916         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
917      CALL histdef( id_i, "sohefldo", "Net Downward Heat Flux", "W/m2"   ,   &   ! net heat flux
918         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
919      CALL histdef( id_i, "soshfldo", "Shortwave Radiation"   , "W/m2"   ,   &   ! solar flux
920         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
921      CALL histdef( id_i, "soicecov", "Ice fraction"          , "[0,1]"  ,   &   ! fr_i
922         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
923      CALL histdef( id_i, "sozotaux", "Zonal Wind Stress"     , "N/m2"   ,   &   ! i-wind stress
924         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
925      CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2"   ,   &   ! j-wind stress
926         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
927      IF( lk_vvl ) THEN
928         CALL histdef( id_i, "vovvldep", "T point depth"         , "m"      ,   &   ! t-point depth
929            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
930      ENDIF
931
932#if defined key_lim2
933      CALL lim_wri_state_2( kt, id_i, nh_i )
934#elif defined key_lim3
935      CALL lim_wri_state( kt, id_i, nh_i )
936#else
937      CALL histend( id_i, snc4chunks=snc4set )
938#endif
939
940      ! 2. Start writing data
941      ! ---------------------
942      ! idex(1) est utilise ssi l'avant dernier argument est diffferent de
943      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
944      ! donne le nombre d'elements, et idex la liste des indices a sortir
945      idex(1) = 1   ! init to avoid compil warning
946
947      ! Write all fields on T grid
948      CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature
949      CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity
950      CALL histwrite( id_i, "sossheig", kt, sshn             , jpi*jpj    , idex )    ! sea surface height
951      CALL histwrite( id_i, "vozocrtx", kt, un               , jpi*jpj*jpk, idex )    ! now i-velocity
952      CALL histwrite( id_i, "vomecrty", kt, vn               , jpi*jpj*jpk, idex )    ! now j-velocity
953      CALL histwrite( id_i, "vovecrtz", kt, wn               , jpi*jpj*jpk, idex )    ! now k-velocity
954      CALL histwrite( id_i, "sowaflup", kt, emp-rnf          , jpi*jpj    , idex )    ! freshwater budget
955      CALL histwrite( id_i, "sohefldo", kt, qsr + qns        , jpi*jpj    , idex )    ! total heat flux
956      CALL histwrite( id_i, "soshfldo", kt, qsr              , jpi*jpj    , idex )    ! solar heat flux
957      CALL histwrite( id_i, "soicecov", kt, fr_i             , jpi*jpj    , idex )    ! ice fraction
958      CALL histwrite( id_i, "sozotaux", kt, utau             , jpi*jpj    , idex )    ! i-wind stress
959      CALL histwrite( id_i, "sometauy", kt, vtau             , jpi*jpj    , idex )    ! j-wind stress
960
961      ! 3. Close the file
962      ! -----------------
963      CALL histclo( id_i )
964#if ! defined key_iomput && ! defined key_dimgout
965      IF( ninist /= 1  ) THEN
966         CALL histclo( nid_T )
967         CALL histclo( nid_U )
968         CALL histclo( nid_V )
969         CALL histclo( nid_W )
970      ENDIF
971#endif
972      !
973   END SUBROUTINE dia_wri_state
974   !!======================================================================
975END MODULE diawri
Note: See TracBrowser for help on using the repository browser.