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_r11514_HPC-02_single-core-extrahalo/tests/CANAL/MY_SRC – NEMO

source: NEMO/branches/2019/dev_r11514_HPC-02_single-core-extrahalo/tests/CANAL/MY_SRC/diawri.F90 @ 11692

Last change on this file since 11692 was 11692, checked in by francesca, 4 years ago

Update branch to integrate the development starting from the current v4.01 ready trunk

  • Property svn:keywords set to Id
File size: 52.9 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 diatmb         ! Top,middle,bottom output
49   USE dia25h         ! 25h Mean output
50   USE iom            !
51   USE ioipsl         !
52
53#if defined key_si3
54   USE ice 
55   USE icewri 
56#endif
57   USE lib_mpp         ! MPP library
58   USE timing          ! preformance summary
59   USE diurnal_bulk    ! diurnal warm layer
60   USE cool_skin       ! Cool skin
61
62   IMPLICIT NONE
63   PRIVATE
64
65   PUBLIC   dia_wri                 ! routines called by step.F90
66   PUBLIC   dia_wri_state
67   PUBLIC   dia_wri_alloc           ! Called by nemogcm module
68
69   INTEGER ::   nid_T, nz_T, nh_T, ndim_T, ndim_hT   ! grid_T file
70   INTEGER ::          nb_T              , ndim_bT   ! grid_T file
71   INTEGER ::   nid_U, nz_U, nh_U, ndim_U, ndim_hU   ! grid_U file
72   INTEGER ::   nid_V, nz_V, nh_V, ndim_V, ndim_hV   ! grid_V file
73   INTEGER ::   nid_W, nz_W, nh_W                    ! grid_W file
74   INTEGER ::   ndex(1)                              ! ???
75   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
76   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
77   INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
78
79   !! * Substitutions
80#  include "vectopt_loop_substitute.h90"
81   !!----------------------------------------------------------------------
82   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
83   !! $Id$
84   !! Software governed by the CeCILL license (see ./LICENSE)
85   !!----------------------------------------------------------------------
86CONTAINS
87
88#if defined key_iomput
89   !!----------------------------------------------------------------------
90   !!   'key_iomput'                                        use IOM library
91   !!----------------------------------------------------------------------
92   INTEGER FUNCTION dia_wri_alloc()
93      !
94      dia_wri_alloc = 0
95      !
96   END FUNCTION dia_wri_alloc
97
98   
99   SUBROUTINE dia_wri( kt )
100      !!---------------------------------------------------------------------
101      !!                  ***  ROUTINE dia_wri  ***
102      !!                   
103      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
104      !!      NETCDF format is used by default
105      !!
106      !! ** Method  :  use iom_put
107      !!----------------------------------------------------------------------
108      INTEGER, INTENT( in ) ::   kt      ! ocean time-step 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      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   bu, bv   ! volume of u- and v-boxes
118      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   r1_bt    ! inverse of t-box volume
119      !!----------------------------------------------------------------------
120      !
121      IF( ln_timing )   CALL timing_start('dia_wri')
122      !
123      ! Output the initial state and forcings
124      IF( ninist == 1 ) THEN                       
125         CALL dia_wri_state( 'output.init' )
126         ninist = 0
127      ENDIF
128
129      ! Output of initial vertical scale factor
130      CALL iom_put("e3t_0", e3t_0(:,:,:) )
131      CALL iom_put("e3u_0", e3u_0(:,:,:) )
132      CALL iom_put("e3v_0", e3v_0(:,:,:) )
133      !
134      CALL iom_put( "e3t" , e3t_n(:,:,:) )
135      CALL iom_put( "e3u" , e3u_n(:,:,:) )
136      CALL iom_put( "e3v" , e3v_n(:,:,:) )
137      CALL iom_put( "e3w" , e3w_n(:,:,:) )
138      IF( iom_use("e3tdef") )   &
139         CALL iom_put( "e3tdef"  , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
140
141      IF( ll_wd ) THEN
142         CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying)
143      ELSE
144         CALL iom_put( "ssh" , sshn )              ! sea surface height
145      ENDIF
146
147      IF( iom_use("wetdep") )   &                  ! wet depth
148         CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) )
149     
150      CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature
151      CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature
152      IF ( iom_use("sbt") ) THEN
153         DO jj = 1, jpj
154            DO ji = 1, jpi
155               ikbot = mbkt(ji,jj)
156               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem)
157            END DO
158         END DO
159         CALL iom_put( "sbt", z2d )                ! bottom temperature
160      ENDIF
161     
162      CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity
163      CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity
164      IF ( iom_use("sbs") ) THEN
165         DO jj = 1, jpj
166            DO ji = 1, jpi
167               ikbot = mbkt(ji,jj)
168               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal)
169            END DO
170         END DO
171         CALL iom_put( "sbs", z2d )                ! bottom salinity
172      ENDIF
173
174      IF ( iom_use("taubot") ) THEN                ! bottom stress
175         zztmp = rau0 * 0.25
176         z2d(:,:) = 0._wp
177         DO jj = 2, jpjm1
178            DO ji = fs_2, fs_jpim1   ! vector opt.
179               zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * un(ji  ,jj,mbku(ji  ,jj))  )**2   &
180                  &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj))  )**2   &
181                  &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vn(ji,jj  ,mbkv(ji,jj  ))  )**2   &
182                  &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1))  )**2
183               z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 
184               !
185            END DO
186         END DO
187         CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
188         CALL iom_put( "taubot", z2d )           
189      ENDIF
190         
191      CALL iom_put( "uoce", un(:,:,:) )            ! 3D i-current
192      CALL iom_put(  "ssu", un(:,:,1) )            ! surface i-current
193      IF ( iom_use("sbu") ) THEN
194         DO jj = 1, jpj
195            DO ji = 1, jpi
196               ikbot = mbku(ji,jj)
197               z2d(ji,jj) = un(ji,jj,ikbot)
198            END DO
199         END DO
200         CALL iom_put( "sbu", z2d )                ! bottom i-current
201      ENDIF
202     
203      CALL iom_put( "voce", vn(:,:,:) )            ! 3D j-current
204      CALL iom_put(  "ssv", vn(:,:,1) )            ! surface j-current
205      IF ( iom_use("sbv") ) THEN
206         DO jj = 1, jpj
207            DO ji = 1, jpi
208               ikbot = mbkv(ji,jj)
209               z2d(ji,jj) = vn(ji,jj,ikbot)
210            END DO
211         END DO
212         CALL iom_put( "sbv", z2d )                ! bottom j-current
213      ENDIF
214
215      CALL iom_put( "woce", wn )                   ! vertical velocity
216      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
217         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
218         z2d(:,:) = rau0 * e1e2t(:,:)
219         DO jk = 1, jpk
220            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
221         END DO
222         CALL iom_put( "w_masstr" , z3d ) 
223         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
224      ENDIF
225
226      CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef.
227      CALL iom_put( "avs" , avs )                  ! S vert. eddy diff. coef.
228      CALL iom_put( "avm" , avm )                  ! T vert. eddy visc. coef.
229
230      IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
231      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
232
233      IF ( iom_use("salgrad") .OR. iom_use("salgrad2") ) THEN
234         z3d(:,:,jpk) = 0.
235         DO jk = 1, jpkm1
236            DO jj = 2, jpjm1                                    ! sal gradient
237               DO ji = fs_2, fs_jpim1   ! vector opt.
238                  zztmp  = tsn(ji,jj,jk,jp_sal)
239                  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)
240                  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)
241                  z3d(ji,jj,jk) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
242                     &                 * umask(ji,jj,jk) * umask(ji-1,jj,jk) * vmask(ji,jj,jk) * umask(ji,jj-1,jk)
243               END DO
244            END DO
245         END DO
246         CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
247         CALL iom_put( "salgrad2",  z3d )          ! square of module of sal gradient
248         z3d(:,:,:) = SQRT( z3d(:,:,:) )
249         CALL iom_put( "salgrad" ,  z3d )          ! module of sal gradient
250      ENDIF
251         
252      IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
253         DO jj = 2, jpjm1                                    ! sst gradient
254            DO ji = fs_2, fs_jpim1   ! vector opt.
255               zztmp  = tsn(ji,jj,1,jp_tem)
256               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)
257               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)
258               z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
259                  &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
260            END DO
261         END DO
262         CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
263         CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient
264         z2d(:,:) = SQRT( z2d(:,:) )
265         CALL iom_put( "sstgrad" ,  z2d )          ! module of sst gradient
266      ENDIF
267         
268      ! heat and salt contents
269      IF( iom_use("heatc") ) THEN
270         z2d(:,:)  = 0._wp 
271         DO jk = 1, jpkm1
272            DO jj = 1, jpj
273               DO ji = 1, jpi
274                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
275               END DO
276            END DO
277         END DO
278         CALL iom_put( "heatc", rau0_rcp * z2d )   ! vertically integrated heat content (J/m2)
279      ENDIF
280
281      IF( iom_use("saltc") ) THEN
282         z2d(:,:)  = 0._wp 
283         DO jk = 1, jpkm1
284            DO jj = 1, jpj
285               DO ji = 1, jpi
286                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
287               END DO
288            END DO
289         END DO
290         CALL iom_put( "saltc", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
291      ENDIF
292      !
293      IF( iom_use("salt2c") ) THEN
294         z2d(:,:)  = 0._wp 
295         DO jk = 1, jpkm1
296            DO jj = 1, jpj
297               DO ji = 1, jpi
298                  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)
299               END DO
300            END DO
301         END DO
302         CALL iom_put( "salt2c", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
303      ENDIF
304      !
305      IF ( iom_use("eken") ) THEN
306         z3d(:,:,jpk) = 0._wp 
307         DO jk = 1, jpkm1
308            DO jj = 2, jpj
309               DO ji = 2, jpi
310                  zztmpx = 0.5 * ( un(ji-1,jj  ,jk) + un(ji,jj,jk) )
311                  zztmpy = 0.5 * ( vn(ji  ,jj-1,jk) + vn(ji,jj,jk) )
312                  z3d(ji,jj,jk) = 0.5 * ( zztmpx*zztmpx + zztmpy*zztmpy )
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 * ( un(ji  ,jj,jk) * un(ji  ,jj,jk) * e1e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)  &
329                     &                      + un(ji-1,jj,jk) * un(ji-1,jj,jk) * e1e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)  &
330                     &                      + vn(ji,jj  ,jk) * vn(ji,jj  ,jk) * e1e2v(ji,jj  ) * e3v_n(ji,jj  ,jk)  &
331                     &                      + vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1e2v(ji,jj-1) * e3v_n(ji,jj-1,jk)  )  &
332                     &                    * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * 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_n(ji,jj,jk) * 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", hdivn )                  ! 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  ) * vn(ji+1,jj  ,jk) - e2v(ji,jj) * vn(ji,jj,jk)    &
361                     &              - e1u(ji  ,jj+1) * un(ji  ,jj+1,jk) + e1u(ji,jj) * un(ji,jj,jk)  ) * 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_n(ji,jj+1,jk)*tmask(ji,jj+1,jk) + e3t_n(ji+1,jj+1,jk)*tmask(ji+1,jj+1,jk)   &
381                     &    + e3t_n(ji,jj  ,jk)*tmask(ji,jj  ,jk) + e3t_n(ji+1,jj  ,jk)*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 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * 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) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
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) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
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 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * 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) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
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) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
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_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem)
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_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
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_diatmb)   CALL dia_tmb                   ! tmb values
496         
497      IF (ln_dia25h)   CALL dia_25h( kt )             ! 25h averaging
498
499      IF( ln_timing )   CALL timing_stop('dia_wri')
500      !
501   END SUBROUTINE dia_wri
502
503#else
504   !!----------------------------------------------------------------------
505   !!   Default option                                  use IOIPSL  library
506   !!----------------------------------------------------------------------
507
508   INTEGER FUNCTION dia_wri_alloc()
509      !!----------------------------------------------------------------------
510      INTEGER, DIMENSION(2) :: ierr
511      !!----------------------------------------------------------------------
512      ierr = 0
513      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
514         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
515         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
516         !
517      dia_wri_alloc = MAXVAL(ierr)
518      CALL mpp_sum( 'diawri', dia_wri_alloc )
519      !
520   END FUNCTION dia_wri_alloc
521
522   
523   SUBROUTINE dia_wri( kt )
524      !!---------------------------------------------------------------------
525      !!                  ***  ROUTINE dia_wri  ***
526      !!                   
527      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
528      !!      NETCDF format is used by default
529      !!
530      !! ** Method  :   At the beginning of the first time step (nit000),
531      !!      define all the NETCDF files and fields
532      !!      At each time step call histdef to compute the mean if ncessary
533      !!      Each nn_write time step, output the instantaneous or mean fields
534      !!----------------------------------------------------------------------
535      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
536      !
537      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
538      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
539      INTEGER  ::   inum = 11                                ! temporary logical unit
540      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
541      INTEGER  ::   ierr                                     ! error code return from allocation
542      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
543      INTEGER  ::   jn, ierror                               ! local integers
544      REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
545      !
546      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace
547      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace
548      !!----------------------------------------------------------------------
549      !
550      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==!
551         CALL dia_wri_state( 'output.init' )
552         ninist = 0
553      ENDIF
554      !
555      IF( nn_write == -1 )   RETURN   ! we will never do any output
556      !
557      IF( ln_timing )   CALL timing_start('dia_wri')
558      !
559      ! 0. Initialisation
560      ! -----------------
561
562      ll_print = .FALSE.                  ! local variable for debugging
563      ll_print = ll_print .AND. lwp
564
565      ! Define frequency of output and means
566      clop = "x"         ! no use of the mask value (require less cpu time and otherwise the model crashes)
567#if defined key_diainstant
568      zsto = nn_write * rdt
569      clop = "inst("//TRIM(clop)//")"
570#else
571      zsto=rdt
572      clop = "ave("//TRIM(clop)//")"
573#endif
574      zout = nn_write * rdt
575      zmax = ( nitend - nit000 + 1 ) * rdt
576
577      ! Define indices of the horizontal output zoom and vertical limit storage
578      iimi = 1      ;      iima = jpi
579      ijmi = 1      ;      ijma = jpj
580      ipk = jpk
581
582      ! define time axis
583      it = kt
584      itmod = kt - nit000 + 1
585
586
587      ! 1. Define NETCDF files and fields at beginning of first time step
588      ! -----------------------------------------------------------------
589
590      IF( kt == nit000 ) THEN
591
592         ! Define the NETCDF files (one per grid)
593
594         ! Compute julian date from starting date of the run
595         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
596         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
597         IF(lwp)WRITE(numout,*)
598         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
599            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
600         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
601                                 ' limit storage in depth = ', ipk
602
603         ! WRITE root name in date.file for use by postpro
604         IF(lwp) THEN
605            CALL dia_nam( clhstnam, nn_write,' ' )
606            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
607            WRITE(inum,*) clhstnam
608            CLOSE(inum)
609         ENDIF
610
611         ! Define the T grid FILE ( nid_T )
612
613         CALL dia_nam( clhstnam, nn_write, 'grid_T' )
614         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
615         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
616            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
617            &          nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
618         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
619            &           "m", ipk, gdept_1d, nz_T, "down" )
620         !                                                            ! Index of ocean points
621         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
622         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
623         !
624         IF( ln_icebergs ) THEN
625            !
626            !! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
627            !! that routine is called from nemogcm, so do it here immediately before its needed
628            ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
629            CALL mpp_sum( 'diawri', ierror )
630            IF( ierror /= 0 ) THEN
631               CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
632               RETURN
633            ENDIF
634            !
635            !! iceberg vertical coordinate is class number
636            CALL histvert( nid_T, "class", "Iceberg class",      &  ! Vertical grid: class
637               &           "number", nclasses, class_num, nb_T )
638            !
639            !! each class just needs the surface index pattern
640            ndim_bT = 3
641            DO jn = 1,nclasses
642               ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
643            ENDDO
644            !
645         ENDIF
646
647         ! Define the U grid FILE ( nid_U )
648
649         CALL dia_nam( clhstnam, nn_write, 'grid_U' )
650         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
651         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
652            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
653            &          nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
654         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
655            &           "m", ipk, gdept_1d, nz_U, "down" )
656         !                                                            ! Index of ocean points
657         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
658         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
659
660         ! Define the V grid FILE ( nid_V )
661
662         CALL dia_nam( clhstnam, nn_write, 'grid_V' )                   ! filename
663         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
664         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
665            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
666            &          nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
667         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
668            &          "m", ipk, gdept_1d, nz_V, "down" )
669         !                                                            ! Index of ocean points
670         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
671         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
672
673         ! Define the W grid FILE ( nid_W )
674
675         CALL dia_nam( clhstnam, nn_write, 'grid_W' )                   ! filename
676         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
677         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
678            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
679            &          nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
680         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
681            &          "m", ipk, gdepw_1d, nz_W, "down" )
682
683
684         ! Declare all the output fields as NETCDF variables
685
686         !                                                                                      !!! nid_T : 3D
687         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
688            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
689         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
690            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
691         IF(  .NOT.ln_linssh  ) THEN
692            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n
693            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
694            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t_n
695            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
696            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t_n
697            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
698         ENDIF
699         !                                                                                      !!! nid_T : 2D
700         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
701            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
702         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
703            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
704         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
705            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
706         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
707            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
708         CALL histdef( nid_T, "sorunoff", "River runoffs"                      , "Kg/m2/s",   &  ! runoffs
709            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
710         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
711            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
712         IF(  ln_linssh  ) THEN
713            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem)
714            &                                                                  , "KgC/m2/s",  &  ! sosst_cd
715            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
716            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * tsn(:,:,1,jp_sal)
717            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
718            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
719         ENDIF
720         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
721            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
722         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
723            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
724         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
725            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
726         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
727            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
728         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
729            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
730         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
731            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
732!
733         IF( ln_icebergs ) THEN
734            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
735               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
736            CALL histdef( nid_T, "calving_heat"        , "calving heat flux"                        , "XXXX"   , &
737               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
738            CALL histdef( nid_T, "berg_floating_melt"  , "Melt rate of icebergs + bits"             , "kg/m2/s", &
739               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
740            CALL histdef( nid_T, "berg_stored_ice"     , "Accumulated ice mass by class"            , "kg"     , &
741               &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
742            IF( ln_bergdia ) THEN
743               CALL histdef( nid_T, "berg_melt"           , "Melt rate of icebergs"                    , "kg/m2/s", &
744                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
745               CALL histdef( nid_T, "berg_buoy_melt"      , "Buoyancy component of iceberg melt rate"  , "kg/m2/s", &
746                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
747               CALL histdef( nid_T, "berg_eros_melt"      , "Erosion component of iceberg melt rate"   , "kg/m2/s", &
748                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
749               CALL histdef( nid_T, "berg_conv_melt"      , "Convective component of iceberg melt rate", "kg/m2/s", &
750                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
751               CALL histdef( nid_T, "berg_virtual_area"   , "Virtual coverage by icebergs"             , "m2"     , &
752                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
753               CALL histdef( nid_T, "bits_src"           , "Mass source of bergy bits"                , "kg/m2/s", &
754                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
755               CALL histdef( nid_T, "bits_melt"          , "Melt rate of bergy bits"                  , "kg/m2/s", &
756                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
757               CALL histdef( nid_T, "bits_mass"          , "Bergy bit density field"                  , "kg/m2"  , &
758                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
759               CALL histdef( nid_T, "berg_mass"           , "Iceberg density field"                    , "kg/m2"  , &
760                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
761               CALL histdef( nid_T, "berg_real_calving"   , "Calving into iceberg class"               , "kg/s"   , &
762                  &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
763            ENDIF
764         ENDIF
765
766         IF( ln_ssr ) THEN
767            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
768               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
769            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
770               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
771            CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
772               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
773         ENDIF
774       
775         clmx ="l_max(only(x))"    ! max index on a period
776!         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
777!            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
778#if defined key_diahth
779         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
780            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
781         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
782            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
783         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
784            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
785         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "J/m2"   ,   & ! htc3
786            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
787#endif
788
789         CALL histend( nid_T, snc4chunks=snc4set )
790
791         !                                                                                      !!! nid_U : 3D
792         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
793            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
794         IF( ln_wave .AND. ln_sdw) THEN
795            CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current"         , "m/s"    ,   &  ! usd
796               &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
797         ENDIF
798         !                                                                                      !!! nid_U : 2D
799         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
800            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
801
802         CALL histend( nid_U, snc4chunks=snc4set )
803
804         !                                                                                      !!! nid_V : 3D
805         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
806            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
807         IF( ln_wave .AND. ln_sdw) THEN
808            CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current"    , "m/s"    ,   &  ! vsd
809               &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
810         ENDIF
811         !                                                                                      !!! nid_V : 2D
812         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
813            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
814
815         CALL histend( nid_V, snc4chunks=snc4set )
816
817         !                                                                                      !!! nid_W : 3D
818         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn
819            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
820         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
821            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
822         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avm
823            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
824
825         IF( ln_zdfddm ) THEN
826            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
827               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
828         ENDIF
829         
830         IF( ln_wave .AND. ln_sdw) THEN
831            CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current"   , "m/s"    ,   &  ! wsd
832               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
833         ENDIF
834         !                                                                                      !!! nid_W : 2D
835         CALL histend( nid_W, snc4chunks=snc4set )
836
837         IF(lwp) WRITE(numout,*)
838         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
839         IF(ll_print) CALL FLUSH(numout )
840
841      ENDIF
842
843      ! 2. Start writing data
844      ! ---------------------
845
846      ! ndex(1) est utilise ssi l'avant dernier argument est different de
847      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
848      ! donne le nombre d'elements, et ndex la liste des indices a sortir
849
850      IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN
851         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
852         WRITE(numout,*) '~~~~~~ '
853      ENDIF
854
855      IF( .NOT.ln_linssh ) THEN
856         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
857         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
858         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
859         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
860      ELSE
861         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature
862         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity
863         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature
864         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
865      ENDIF
866      IF( .NOT.ln_linssh ) THEN
867         zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
868         CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
869         CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
870         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
871      ENDIF
872      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
873      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
874      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
875      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
876                                                                                  ! (includes virtual salt flux beneath ice
877                                                                                  ! in linear free surface case)
878      IF( ln_linssh ) THEN
879         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
880         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
881         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
882         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
883      ENDIF
884      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
885      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
886      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
887      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
888      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
889      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
890!
891      IF( ln_icebergs ) THEN
892         !
893         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT ) 
894         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
895         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT ) 
896         !
897         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
898         !
899         IF( ln_bergdia ) THEN
900            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   ) 
901            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   ) 
902            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   ) 
903            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   ) 
904            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   ) 
905            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   ) 
906            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   ) 
907            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   ) 
908            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   ) 
909            !
910            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
911         ENDIF
912      ENDIF
913
914      IF( ln_ssr ) THEN
915         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
916         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
917         zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
918         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
919      ENDIF
920!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
921!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
922
923#if defined key_diahth
924      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
925      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
926      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
927      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
928#endif
929
930      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
931      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
932
933      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
934      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
935
936      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
937      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
938      CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
939      IF( ln_zdfddm ) THEN
940         CALL histwrite( nid_W, "voddmavs", it, avs         , ndim_T, ndex_T )    ! S vert. eddy diff. coef.
941      ENDIF
942
943      IF( ln_wave .AND. ln_sdw ) THEN
944         CALL histwrite( nid_U, "sdzocrtx", it, usd         , ndim_U , ndex_U )    ! i-StokesDrift-current
945         CALL histwrite( nid_V, "sdmecrty", it, vsd         , ndim_V , ndex_V )    ! j-StokesDrift-current
946         CALL histwrite( nid_W, "sdvecrtz", it, wsd         , ndim_T , ndex_T )    ! StokesDrift vert. current
947      ENDIF
948
949      ! 3. Close all files
950      ! ---------------------------------------
951      IF( kt == nitend ) THEN
952         CALL histclo( nid_T )
953         CALL histclo( nid_U )
954         CALL histclo( nid_V )
955         CALL histclo( nid_W )
956      ENDIF
957      !
958      IF( ln_timing )   CALL timing_stop('dia_wri')
959      !
960   END SUBROUTINE dia_wri
961#endif
962
963   SUBROUTINE dia_wri_state( cdfile_name )
964      !!---------------------------------------------------------------------
965      !!                 ***  ROUTINE dia_wri_state  ***
966      !!       
967      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
968      !!      the instantaneous ocean state and forcing fields.
969      !!        Used to find errors in the initial state or save the last
970      !!      ocean state in case of abnormal end of a simulation
971      !!
972      !! ** Method  :   NetCDF files using ioipsl
973      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
974      !!      File 'output.abort.nc' is created in case of abnormal job end
975      !!----------------------------------------------------------------------
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', tsn(:,:,:,jp_tem) )    ! now temperature
993      CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) )    ! now salinity
994      CALL iom_rstput( 0, 0, inum, 'sossheig', sshn              )    ! sea surface height
995      CALL iom_rstput( 0, 0, inum, 'vozocrtx', un                )    ! now i-velocity
996      CALL iom_rstput( 0, 0, inum, 'vomecrty', vn                )    ! now j-velocity
997      CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn                )    ! 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_n        )    !  T-cell depth
1014         CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n          )    !  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.