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_r12072_MERGE_OPTION2_2019/src/OCE/DIA – NEMO

source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/DIA/diawri.F90 @ 12251

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

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