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/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA – NEMO

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DIA/diawri.F90 @ 13193

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

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