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_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA – NEMO

source: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/DIA/diawri.F90 @ 11363

Last change on this file since 11363 was 11363, checked in by gsamson, 4 years ago

dev_r11265_ABL : see #2131

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