New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diawri.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/tests/CANAL/MY_SRC – NEMO

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

Last change on this file since 11949 was 11949, checked in by acc, 4 years ago

Merge in changes from 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. This just creates a fresh copy of this branch to use as the merge base. See ticket #2341

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