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 @ 11334

Last change on this file since 11334 was 11334, checked in by flemarie, 5 years ago

Follow-on to the ABL implementation (see ticket #2131)

  • Modification of diawri.F90 to output ABL variables with IOIPSL (i.e. without key_iomput)
  • Modification of ice_var_agg in icevar.F90 to allow the use of tm_su in ABL and bulk
  • Error handling in case ln_abl = .true. and ABL sources were not compiled
  • ABL now works with SI3 considering an average over ice categories
  • Update reference namelist to avoid runtime warnings due to nam_tide
  • Sanity checks with ORCA2 for the following configurations :

1 - ABL src + ln_blk = .true.
2 - ABL src + ln_abl = .true.
3 - no ABL src + ln_blk = .true.

All configurations are MPP reproducible and configurations 1 and 3 results are identical

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