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

source: branches/UKMO/AMM15_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 7976

Last change on this file since 7976 was 7976, checked in by jgraham, 7 years ago

Changes to available diagnostics:
Kara MLD diagnostics:
NEMOGCM/CONFIG/SHARED/namelist_ref
NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
Diagnostics match those introduced in the following branch:
http://fcm3/projects/NEMO.xm/log/branches/UKMO/dev_r5107_mld_zint

Output rhop:
NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
Allow for output of rhop without key_diaar5.

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