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

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

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