source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 3 years ago

All wrk_alloc removed

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