source: branches/2015/dev_r5094_UKMO_ISFCLEAN/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 5104

Last change on this file since 5104 was 5104, checked in by mathiot, 6 years ago

correction of minor bug with isf + add some extra test on ln_isfcav + restore ssu/ssv

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