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

Last change on this file since 12731 was 12731, checked in by techene, 4 years ago

replace h. and gde. in case key_qco is activated - quick and dirty

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