source: NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/CANAL/MY_SRC/diawri.F90 @ 12443

Last change on this file since 12443 was 12443, checked in by davestorkey, 9 months ago

2020/KERNEL-03_Storkey_Coward_RK3_stage2: More variable renaming:
atfp → rn_atfp (use namelist parameter everywhere)
rdtbt → rDt_e
nn_baro → nn_e
rn_scal_load → rn_load
rau0 → rho0

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