New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diawri.F90 in branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 4616

Last change on this file since 4616 was 4616, checked in by gm, 10 years ago

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