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

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 4147

Last change on this file since 4147 was 3704, checked in by smasson, 11 years ago

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