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

Last change on this file since 3419 was 3419, checked in by acc, 9 years ago

Branch: dev_r3385_NOCS04_HAMF; #665. Stage 4 of 2012 development: Changes to get LIM3 working with embedded sea-ice. Working with a reduced (halved) timestep but exhibiting the same stability problems as LIM2_EVP with standard ORCA2 settings.

  • 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!!$#if defined key_lim3 || defined key_lim2
375!!$         ! sowaflup = sowaflep + sorunoff + sowafldp + a term associated to
376!!$         !    internal damping to Levitus that can be diagnosed from others
377!!$         ! sowaflcd = sowaflep + sorunoff + sowafldp + iowaflup
378!!$         CALL histdef( nid_T, "iowaflup", "Ice=>ocean net freshwater"          , "kg/m2/s",   &  ! fsalt
379!!$            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
380!!$         CALL histdef( nid_T, "sowaflep", "atmos=>ocean net freshwater"        , "kg/m2/s",   &  ! fmass
381!!$            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
382!!$#endif
383         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
384            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
385!!$         CALL histdef( nid_T, "sorunoff", "Runoffs"                            , "Kg/m2/s",   &  ! runoffs
386!!$            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
387         CALL histdef( nid_T, "sowaflcd", "concentration/dilution water flux"  , "kg/m2/s",   &  ! (emps-rnf)
388            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
389         CALL histdef( nid_T, "sosalflx", "Surface Salt Flux"                  , "Kg/m2/s",   &  ! (emps-rnf) * sn
390            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
391         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
392            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
393         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
394            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
395         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
396            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
397         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
398            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
399         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
400            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
401         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
402            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
403#if ! defined key_coupled
404         CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
405            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
406         CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
407            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
408         CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
409            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
410#endif
411
412
413
414#if ( defined key_coupled && ! defined key_lim3 && ! defined key_lim2 )
415         CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
416            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
417         CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
418            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
419         CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn
420            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
421#endif
422         clmx ="l_max(only(x))"    ! max index on a period
423         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
424            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
425#if defined key_diahth
426         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
427            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
428         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
429            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
430         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
431            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
432         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "W"      ,   & ! htc3
433            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
434#endif
435
436#if defined key_coupled 
437# if defined key_lim3
438         Must be adapted to LIM3
439# else
440         CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature"            , "K"      ,   &  ! tn_ice
441            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
442         CALL histdef( nid_T,"soicealb" , "Ice Albedo"                         , "[0,1]"  ,   &  ! alb_ice
443            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
444# endif 
445#endif
446
447         CALL histend( nid_T, snc4chunks=snc4set )
448
449         !                                                                                      !!! nid_U : 3D
450         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
451            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
452         IF( ln_traldf_gdia ) THEN
453            CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
454                 &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
455         ELSE
456#if defined key_diaeiv
457            CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
458            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
459#endif
460         END IF
461         !                                                                                      !!! nid_U : 2D
462         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
463            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
464
465         CALL histend( nid_U, snc4chunks=snc4set )
466
467         !                                                                                      !!! nid_V : 3D
468         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
469            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
470         IF( ln_traldf_gdia ) THEN
471            CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
472                 &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
473         ELSE 
474#if defined key_diaeiv
475            CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
476            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
477#endif
478         END IF
479         !                                                                                      !!! nid_V : 2D
480         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
481            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
482
483         CALL histend( nid_V, snc4chunks=snc4set )
484
485         !                                                                                      !!! nid_W : 3D
486         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn
487            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
488         IF( ln_traldf_gdia ) THEN
489            CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
490                 &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
491         ELSE
492#if defined key_diaeiv
493            CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
494                 &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
495#endif
496         END IF
497         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
498            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
499         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avmu
500            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
501
502         IF( lk_zdfddm ) THEN
503            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
504               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
505         ENDIF
506         !                                                                                      !!! nid_W : 2D
507#if defined key_traldf_c2d
508         CALL histdef( nid_W, "soleahtw", "lateral eddy diffusivity"           , "m2/s"   ,   &  ! ahtw
509            &          jpi, jpj, nh_W, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
510# if defined key_traldf_eiv 
511            CALL histdef( nid_W, "soleaeiw", "eddy induced vel. coeff. at w-point", "m2/s",   &  ! aeiw
512               &       jpi, jpj, nh_W, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
513# endif
514#endif
515
516         CALL histend( nid_W, snc4chunks=snc4set )
517
518         IF(lwp) WRITE(numout,*)
519         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
520         IF(ll_print) CALL FLUSH(numout )
521
522      ENDIF
523
524      ! 2. Start writing data
525      ! ---------------------
526
527      ! ndex(1) est utilise ssi l'avant dernier argument est diffferent de
528      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
529      ! donne le nombre d'elements, et ndex la liste des indices a sortir
530
531      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
532         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
533         WRITE(numout,*) '~~~~~~ '
534      ENDIF
535
536      ! Write fields on T grid
537      CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem), ndim_T , ndex_T  )   ! temperature
538      CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal), ndim_T , ndex_T  )   ! salinity
539      CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem), ndim_hT, ndex_hT )   ! sea surface temperature
540      CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal), ndim_hT, ndex_hT )   ! sea surface salinity
541      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
542!!$#if  defined key_lim3 || defined key_lim2
543!!$      CALL histwrite( nid_T, "iowaflup", it, fsalt(:,:)    , ndim_hT, ndex_hT )   ! ice=>ocean water flux
544!!$      CALL histwrite( nid_T, "sowaflep", it, fmass(:,:)    , ndim_hT, ndex_hT )   ! atmos=>ocean water flux
545!!$#endif
546      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
547!!$      CALL histwrite( nid_T, "sorunoff", it, runoff        , ndim_hT, ndex_hT )   ! runoff
548      CALL histwrite( nid_T, "sowaflcd", it, ( sfx -rnf )  , ndim_hT, ndex_hT )   ! c/d water flux
549      zw2d(:,:) = ( sfx (:,:) - rnf(:,:) ) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
550      CALL histwrite( nid_T, "sosalflx", it, zw2d          , ndim_hT, ndex_hT )   ! c/d salt flux
551      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
552      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
553      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
554      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
555      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
556      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
557#if ! defined key_coupled
558      CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
559      CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
560      IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
561      CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
562#endif
563#if ( defined key_coupled && ! defined key_lim3 && ! defined key_lim2 )
564      CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
565      CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
566         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
567      CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
568#endif
569      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
570      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
571
572#if defined key_diahth
573      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
574      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
575      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
576      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
577#endif
578
579#if defined key_coupled 
580# if defined key_lim3
581      Must be adapted for LIM3
582      CALL histwrite( nid_T, "soicetem", it, tn_ice        , ndim_hT, ndex_hT )   ! surf. ice temperature
583      CALL histwrite( nid_T, "soicealb", it, alb_ice       , ndim_hT, ndex_hT )   ! ice albedo
584# else
585      CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature
586      CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo
587# endif
588#endif
589         ! Write fields on U grid
590      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
591      IF( ln_traldf_gdia ) THEN
592         IF (.not. ALLOCATED(psix_eiv))THEN
593            ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr )
594            IF( lk_mpp   )   CALL mpp_sum ( ierr )
595            IF( ierr > 0 )   CALL ctl_stop('STOP', 'diawri: unable to allocate psi{x,y}_eiv')
596            psix_eiv(:,:,:) = 0.0_wp
597            psiy_eiv(:,:,:) = 0.0_wp
598         ENDIF
599         DO jk=1,jpkm1
600            zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk)  ! u_eiv = -dpsix/dz
601         END DO
602         zw3d(:,:,jpk) = 0._wp
603         CALL histwrite( nid_U, "vozoeivu", it, zw3d, ndim_U , ndex_U )           ! i-eiv current
604      ELSE
605#if defined key_diaeiv
606         CALL histwrite( nid_U, "vozoeivu", it, u_eiv, ndim_U , ndex_U )          ! i-eiv current
607#endif
608      ENDIF
609      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
610
611         ! Write fields on V grid
612      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
613      IF( ln_traldf_gdia ) THEN
614         DO jk=1,jpk-1
615            zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk)  ! v_eiv = -dpsiy/dz
616         END DO
617         zw3d(:,:,jpk) = 0._wp
618         CALL histwrite( nid_V, "vomeeivv", it, zw3d, ndim_V , ndex_V )           ! j-eiv current
619      ELSE
620#if defined key_diaeiv
621         CALL histwrite( nid_V, "vomeeivv", it, v_eiv, ndim_V , ndex_V )          ! j-eiv current
622#endif
623      ENDIF
624      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
625
626         ! Write fields on W grid
627      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
628      IF( ln_traldf_gdia ) THEN
629         DO jk=1,jpk-1
630            DO jj = 2, jpjm1
631               DO ji = fs_2, fs_jpim1  ! vector opt.
632                  zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2v(ji,jj) + &
633                       &    (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1u(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx
634               END DO
635            END DO
636         END DO
637         zw3d(:,:,jpk) = 0._wp
638         CALL histwrite( nid_W, "voveeivw", it, zw3d          , ndim_T, ndex_T )    ! vert. eiv current
639      ELSE
640#   if defined key_diaeiv
641         CALL histwrite( nid_W, "voveeivw", it, w_eiv          , ndim_T, ndex_T )    ! vert. eiv current
642#   endif
643      ENDIF
644      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
645      CALL histwrite( nid_W, "votkeavm", it, avmu           , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
646      IF( lk_zdfddm ) THEN
647         CALL histwrite( nid_W, "voddmavs", it, fsavs(:,:,:), ndim_T, ndex_T )    ! S vert. eddy diff. coef.
648      ENDIF
649#if defined key_traldf_c2d
650      CALL histwrite( nid_W, "soleahtw", it, ahtw          , ndim_hT, ndex_hT )   ! lateral eddy diff. coef.
651# if defined key_traldf_eiv
652         CALL histwrite( nid_W, "soleaeiw", it, aeiw       , ndim_hT, ndex_hT )   ! EIV coefficient at w-point
653# endif
654#endif
655
656      ! 3. Close all files
657      ! ---------------------------------------
658      IF( kt == nitend ) THEN
659         CALL histclo( nid_T )
660         CALL histclo( nid_U )
661         CALL histclo( nid_V )
662         CALL histclo( nid_W )
663      ENDIF
664      !
665      CALL wrk_dealloc( jpi , jpj      , zw2d )
666      IF ( ln_traldf_gdia )  call wrk_dealloc( jpi , jpj , jpk  , zw3d )
667      !
668      IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
669      !
670   END SUBROUTINE dia_wri
671# endif
672
673#endif
674
675   SUBROUTINE dia_wri_state( cdfile_name, kt )
676      !!---------------------------------------------------------------------
677      !!                 ***  ROUTINE dia_wri_state  ***
678      !!       
679      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
680      !!      the instantaneous ocean state and forcing fields.
681      !!        Used to find errors in the initial state or save the last
682      !!      ocean state in case of abnormal end of a simulation
683      !!
684      !! ** Method  :   NetCDF files using ioipsl
685      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
686      !!      File 'output.abort.nc' is created in case of abnormal job end
687      !!----------------------------------------------------------------------
688      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
689      INTEGER           , INTENT( in ) ::   kt               ! ocean time-step index
690      !!
691      CHARACTER (len=32) :: clname
692      CHARACTER (len=40) :: clop
693      INTEGER  ::   id_i , nz_i, nh_i       
694      INTEGER, DIMENSION(1) ::   idex             ! local workspace
695      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt
696      !!----------------------------------------------------------------------
697      !
698!     IF( nn_timing == 1 )   CALL timing_start('dia_wri_state') ! not sure this works for routines not called in first timestep
699
700      ! 0. Initialisation
701      ! -----------------
702
703      ! Define name, frequency of output and means
704      clname = cdfile_name
705      IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
706      zdt  = rdt
707      zsto = rdt
708      clop = "inst(x)"           ! no use of the mask value (require less cpu time)
709      zout = rdt
710      zmax = ( nitend - nit000 + 1 ) * zdt
711
712      IF(lwp) WRITE(numout,*)
713      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
714      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
715      IF(lwp) WRITE(numout,*) '                and named :', clname, '.nc'
716
717
718      ! 1. Define NETCDF files and fields at beginning of first time step
719      ! -----------------------------------------------------------------
720
721      ! Compute julian date from starting date of the run
722      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )         ! time axis
723      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
724      CALL histbeg( clname, jpi, glamt, jpj, gphit,   &
725          1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
726      CALL histvert( id_i, "deptht", "Vertical T levels",   &    ! Vertical grid : gdept
727          "m", jpk, gdept_0, nz_i, "down")
728
729      ! Declare all the output fields as NetCDF variables
730
731      CALL histdef( id_i, "vosaline", "Salinity"              , "PSU"    ,   &   ! salinity
732         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
733      CALL histdef( id_i, "votemper", "Temperature"           , "C"      ,   &   ! temperature
734         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
735      CALL histdef( id_i, "sossheig", "Sea Surface Height"    , "m"      ,   &  ! ssh
736         &          jpi, jpj, nh_i, 1  , 1, 1  , nz_i, 32, clop, zsto, zout )
737      CALL histdef( id_i, "vozocrtx", "Zonal Current"         , "m/s"    ,   &   ! zonal current
738         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
739      CALL histdef( id_i, "vomecrty", "Meridional Current"    , "m/s"    ,   &   ! meridonal current
740         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
741      CALL histdef( id_i, "vovecrtz", "Vertical Velocity"     , "m/s"    ,   &   ! vertical current
742         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
743      CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S",   &   ! net freshwater
744         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
745      CALL histdef( id_i, "sohefldo", "Net Downward Heat Flux", "W/m2"   ,   &   ! net heat flux
746         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
747      CALL histdef( id_i, "soshfldo", "Shortwave Radiation"   , "W/m2"   ,   &   ! solar flux
748         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
749      CALL histdef( id_i, "soicecov", "Ice fraction"          , "[0,1]"  ,   &   ! fr_i
750         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
751      CALL histdef( id_i, "sozotaux", "Zonal Wind Stress"     , "N/m2"   ,   &   ! i-wind stress
752         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
753      CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2"   ,   &   ! j-wind stress
754         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
755
756#if defined key_lim2
757      CALL lim_wri_state_2( kt, id_i, nh_i )
758#else
759      CALL histend( id_i, snc4chunks=snc4set )
760#endif
761
762      ! 2. Start writing data
763      ! ---------------------
764      ! idex(1) est utilise ssi l'avant dernier argument est diffferent de
765      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
766      ! donne le nombre d'elements, et idex la liste des indices a sortir
767      idex(1) = 1   ! init to avoid compil warning
768
769      ! Write all fields on T grid
770      CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature
771      CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity
772      CALL histwrite( id_i, "sossheig", kt, sshn             , jpi*jpj    , idex )    ! sea surface height
773      CALL histwrite( id_i, "vozocrtx", kt, un               , jpi*jpj*jpk, idex )    ! now i-velocity
774      CALL histwrite( id_i, "vomecrty", kt, vn               , jpi*jpj*jpk, idex )    ! now j-velocity
775      CALL histwrite( id_i, "vovecrtz", kt, wn               , jpi*jpj*jpk, idex )    ! now k-velocity
776      CALL histwrite( id_i, "sowaflup", kt, (emp-rnf )       , jpi*jpj    , idex )    ! freshwater budget
777      CALL histwrite( id_i, "sohefldo", kt, qsr + qns        , jpi*jpj    , idex )    ! total heat flux
778      CALL histwrite( id_i, "soshfldo", kt, qsr              , jpi*jpj    , idex )    ! solar heat flux
779      CALL histwrite( id_i, "soicecov", kt, fr_i             , jpi*jpj    , idex )    ! ice fraction
780      CALL histwrite( id_i, "sozotaux", kt, utau             , jpi*jpj    , idex )    ! i-wind stress
781      CALL histwrite( id_i, "sometauy", kt, vtau             , jpi*jpj    , idex )    ! j-wind stress
782
783      ! 3. Close the file
784      ! -----------------
785      CALL histclo( id_i )
786#if ! defined key_iomput && ! defined key_dimgout
787      IF( ninist /= 1  ) THEN
788         CALL histclo( nid_T )
789         CALL histclo( nid_U )
790         CALL histclo( nid_V )
791         CALL histclo( nid_W )
792      ENDIF
793#endif
794       
795!     IF( nn_timing == 1 )   CALL timing_stop('dia_wri_state') ! not sure this works for routines not called in first timestep
796      !
797
798   END SUBROUTINE dia_wri_state
799   !!======================================================================
800END MODULE diawri
Note: See TracBrowser for help on using the repository browser.