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

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

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

Last change on this file since 5098 was 5098, checked in by mathiot, 9 years ago

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