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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/CANAL/MY_SRC/diawri.F90 @ 11771

Last change on this file since 11771 was 11771, checked in by acc, 5 years ago

Branch 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. Updates to MY_SRC files in CANAL test case (previously untested on this branch).
  • Property svn:keywords set to Id
File size: 54.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 nwrite 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( ln_timing )   CALL timing_start('dia_wri')
553      !
554      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==!
555         CALL dia_wri_state( Kmm, 'output.init' )
556         ninist = 0
557      ENDIF
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 = nwrite * rdt
569      clop = "inst("//TRIM(clop)//")"
570#else
571      zsto=rdt
572      clop = "ave("//TRIM(clop)//")"
573#endif
574      zout = nwrite * 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, nwrite,' ' )
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, nwrite, '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, nwrite, '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, nwrite, '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, nwrite, '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(:,:,:,Kmm)
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(:,:,:,Kmm)
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(:,:,:,Kmm)
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 * ts(:,:,1,jp_tem,Kmm)
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 * ts(:,:,1,jp_sal,Kmm)
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( .NOT. ln_cpl ) 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         IF( ln_cpl .AND. nn_ice <= 1 ) THEN
776            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
777               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
778            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
779               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
780            CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn
781               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
782         ENDIF
783         
784         clmx ="l_max(only(x))"    ! max index on a period
785!         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
786!            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
787#if defined key_diahth
788         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
789            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
790         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
791            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
792         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
793            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
794         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "J/m2"   ,   & ! htc3
795            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
796#endif
797
798         CALL histend( nid_T, snc4chunks=snc4set )
799
800         !                                                                                      !!! nid_U : 3D
801         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! uu(:,:,:,Kmm)
802            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
803         IF( ln_wave .AND. ln_sdw) THEN
804            CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current"         , "m/s"    ,   &  ! usd
805               &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
806         ENDIF
807         !                                                                                      !!! nid_U : 2D
808         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
809            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
810
811         CALL histend( nid_U, snc4chunks=snc4set )
812
813         !                                                                                      !!! nid_V : 3D
814         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vv(:,:,:,Kmm)
815            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
816         IF( ln_wave .AND. ln_sdw) THEN
817            CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current"    , "m/s"    ,   &  ! vsd
818               &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
819         ENDIF
820         !                                                                                      !!! nid_V : 2D
821         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
822            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
823
824         CALL histend( nid_V, snc4chunks=snc4set )
825
826         !                                                                                      !!! nid_W : 3D
827         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! ww
828            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
829         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
830            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
831         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avm
832            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
833
834         IF( ln_zdfddm ) THEN
835            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
836               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
837         ENDIF
838         
839         IF( ln_wave .AND. ln_sdw) THEN
840            CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current"   , "m/s"    ,   &  ! wsd
841               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
842         ENDIF
843         !                                                                                      !!! nid_W : 2D
844         CALL histend( nid_W, snc4chunks=snc4set )
845
846         IF(lwp) WRITE(numout,*)
847         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
848         IF(ll_print) CALL FLUSH(numout )
849
850      ENDIF
851
852      ! 2. Start writing data
853      ! ---------------------
854
855      ! ndex(1) est utilise ssi l'avant dernier argument est different de
856      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
857      ! donne le nombre d'elements, et ndex la liste des indices a sortir
858
859      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
860         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
861         WRITE(numout,*) '~~~~~~ '
862      ENDIF
863
864      IF( .NOT.ln_linssh ) THEN
865         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T  )   ! heat content
866         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * e3t(:,:,:,Kmm) , ndim_T , ndex_T  )   ! salt content
867         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface heat content
868         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * e3t(:,:,1,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity content
869      ELSE
870         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T  )   ! temperature
871         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) , ndim_T , ndex_T  )   ! salinity
872         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) , ndim_hT, ndex_hT )   ! sea surface temperature
873         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity
874      ENDIF
875      IF( .NOT.ln_linssh ) THEN
876         zw3d(:,:,:) = ( ( e3t(:,:,:,Kmm) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
877         CALL histwrite( nid_T, "vovvle3t", it, e3t (:,:,:,Kmm) , ndim_T , ndex_T  )   ! level thickness
878         CALL histwrite( nid_T, "vovvldep", it, gdept(:,:,:,Kmm) , ndim_T , ndex_T  )   ! t-point depth
879         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
880      ENDIF
881      CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm)          , ndim_hT, ndex_hT )   ! sea surface height
882      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
883      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
884      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
885                                                                                  ! (includes virtual salt flux beneath ice
886                                                                                  ! in linear free surface case)
887      IF( ln_linssh ) THEN
888         zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_tem,Kmm)
889         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
890         zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_sal,Kmm)
891         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
892      ENDIF
893      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
894      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
895      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
896      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
897      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
898      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
899!
900      IF( ln_icebergs ) THEN
901         !
902         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT ) 
903         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
904         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT ) 
905         !
906         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
907         !
908         IF( ln_bergdia ) THEN
909            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   ) 
910            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   ) 
911            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   ) 
912            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   ) 
913            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   ) 
914            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   ) 
915            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   ) 
916            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   ) 
917            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   ) 
918            !
919            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
920         ENDIF
921      ENDIF
922
923      IF( .NOT. ln_cpl ) THEN
924         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
925         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
926         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * ts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
927         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
928      ENDIF
929      IF( ln_cpl .AND. nn_ice <= 1 ) THEN
930         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
931         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
932         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * ts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
933         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
934      ENDIF
935!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
936!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
937
938#if defined key_diahth
939      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
940      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
941      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
942      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
943#endif
944
945      CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm)            , ndim_U , ndex_U )    ! i-current
946      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
947
948      CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm)            , ndim_V , ndex_V  )   ! j-current
949      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
950
951      CALL histwrite( nid_W, "vovecrtz", it, ww             , ndim_T, ndex_T )    ! vert. current
952      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
953      CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
954      IF( ln_zdfddm ) THEN
955         CALL histwrite( nid_W, "voddmavs", it, avs         , ndim_T, ndex_T )    ! S vert. eddy diff. coef.
956      ENDIF
957
958      IF( ln_wave .AND. ln_sdw ) THEN
959         CALL histwrite( nid_U, "sdzocrtx", it, usd         , ndim_U , ndex_U )    ! i-StokesDrift-current
960         CALL histwrite( nid_V, "sdmecrty", it, vsd         , ndim_V , ndex_V )    ! j-StokesDrift-current
961         CALL histwrite( nid_W, "sdvecrtz", it, wsd         , ndim_T , ndex_T )    ! StokesDrift vert. current
962      ENDIF
963
964      ! 3. Close all files
965      ! ---------------------------------------
966      IF( kt == nitend ) THEN
967         CALL histclo( nid_T )
968         CALL histclo( nid_U )
969         CALL histclo( nid_V )
970         CALL histclo( nid_W )
971      ENDIF
972      !
973      IF( ln_timing )   CALL timing_stop('dia_wri')
974      !
975   END SUBROUTINE dia_wri
976#endif
977
978   SUBROUTINE dia_wri_state( Kmm, cdfile_name )
979      !!---------------------------------------------------------------------
980      !!                 ***  ROUTINE dia_wri_state  ***
981      !!       
982      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
983      !!      the instantaneous ocean state and forcing fields.
984      !!        Used to find errors in the initial state or save the last
985      !!      ocean state in case of abnormal end of a simulation
986      !!
987      !! ** Method  :   NetCDF files using ioipsl
988      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
989      !!      File 'output.abort.nc' is created in case of abnormal job end
990      !!----------------------------------------------------------------------
991      INTEGER           , INTENT( in ) ::   Kmm              ! time level index
992      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
993      !!
994      INTEGER :: inum
995      !!----------------------------------------------------------------------
996      !
997      IF(lwp) WRITE(numout,*)
998      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
999      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
1000      IF(lwp) WRITE(numout,*) '                and named :', cdfile_name, '...nc'
1001
1002#if defined key_si3
1003     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
1004#else
1005     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
1006#endif
1007
1008      CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) )    ! now temperature
1009      CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) )    ! now salinity
1010      CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm)              )    ! sea surface height
1011      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)                )    ! now i-velocity
1012      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)                )    ! now j-velocity
1013      CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww                )    ! now k-velocity
1014      IF( ALLOCATED(ahtu) ) THEN
1015         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point
1016         CALL iom_rstput( 0, 0, inum,  'ahtv', ahtv              )    ! aht at v-point
1017      ENDIF
1018      IF( ALLOCATED(ahmt) ) THEN
1019         CALL iom_rstput( 0, 0, inum,  'ahmt', ahmt              )    ! ahmt at u-point
1020         CALL iom_rstput( 0, 0, inum,  'ahmf', ahmf              )    ! ahmf at v-point
1021      ENDIF
1022      CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf         )    ! freshwater budget
1023      CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns         )    ! total heat flux
1024      CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr               )    ! solar heat flux
1025      CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i              )    ! ice fraction
1026      CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
1027      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
1028      IF(  .NOT.ln_linssh  ) THEN             
1029         CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept(:,:,:,Kmm)        )    !  T-cell depth
1030         CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t(:,:,:,Kmm)          )    !  T-cell thickness 
1031      END IF
1032      IF( ln_wave .AND. ln_sdw ) THEN
1033         CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd            )    ! now StokesDrift i-velocity
1034         CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd            )    ! now StokesDrift j-velocity
1035         CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity
1036      ENDIF
1037 
1038#if defined key_si3
1039      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid
1040         CALL ice_wri_state( inum )
1041      ENDIF
1042#endif
1043      !
1044      CALL iom_close( inum )
1045      !
1046   END SUBROUTINE dia_wri_state
1047
1048   !!======================================================================
1049END MODULE diawri
Note: See TracBrowser for help on using the repository browser.