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/2020/temporary_r4_trunk/tests/CANAL/MY_SRC – NEMO

source: NEMO/branches/2020/temporary_r4_trunk/tests/CANAL/MY_SRC/diawri.F90 @ 13469

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

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