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

source: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 3865

Last change on this file since 3865 was 3865, checked in by acc, 11 years ago

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