New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diawri.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

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

Last change on this file since 5407 was 5407, checked in by smasson, 10 years ago

merge dev_r5218_CNRS17_coupling into the trunk

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