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

source: branches/2016/dev_r7233_CMIP6_diags_trunk_version/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 7240

Last change on this file since 7240 was 7240, checked in by timgraham, 7 years ago

Added vertically integrated zonal mass transport, temperature and salinity
Added 3D ice shelf budget terms

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