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

source: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 10958

Last change on this file since 10958 was 10958, checked in by rrenshaw, 5 years ago

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