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

source: branches/2015/dev_CMCC_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 6051

Last change on this file since 6051 was 6051, checked in by lovato, 9 years ago

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