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

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

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