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 @ 4161

Last change on this file since 4161 was 4161, checked in by cetlod, 10 years ago

dev_LOCEAN_2013 : merge in the 3rd dev branch dev_r4028_CNRS_LIM3, see ticket #1169

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