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

source: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 3488

Last change on this file since 3488 was 3488, checked in by acc, 12 years ago

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