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/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 5433

Last change on this file since 5433 was 5260, checked in by deazer, 9 years ago

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