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 NEMO/branches/2019/dev_r11943_MERGE_2019/tests/CANAL/MY_SRC – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/tests/CANAL/MY_SRC/diawri.F90 @ 12252

Last change on this file since 12252 was 12252, checked in by smasson, 4 years ago

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