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

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

source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 4924

Last change on this file since 4924 was 4924, checked in by mathiot, 9 years ago

UKM02_ice_shelves merged and SETTE tested with revision 4879 of trunk

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