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

source: branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 3625

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

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