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

source: branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 6808

Last change on this file since 6808 was 6808, checked in by jamesharle, 8 years ago

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