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

source: branches/2014/dev_r4650_UKMO13_CICE_changes_take2/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 4921

Last change on this file since 4921 was 4921, checked in by timgraham, 9 years ago

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