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

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

branch KERNEL-06 : merge with trunk@12698 #2385 - in duplcated files : changes to comply to the new trunk variables and some loop bug fixes
  • Property svn:keywords set to Id
File size: 56.3 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       ! 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      END DO
482
483
484      ! 1. Define NETCDF files and fields at beginning of first time step
485      ! -----------------------------------------------------------------
486
487      IF( kt == nit000 ) THEN
488
489         ! Define the NETCDF files (one per grid)
490
491         ! Compute julian date from starting date of the run
492         CALL ymds2ju( nyear, nmonth, nday, rn_Dt, zjulian )
493         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
494         IF(lwp)WRITE(numout,*)
495         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
496            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
497         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
498                                 ' limit storage in depth = ', ipk
499
500         ! WRITE root name in date.file for use by postpro
501         IF(lwp) THEN
502            CALL dia_nam( clhstnam, nn_write,' ' )
503            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
504            WRITE(inum,*) clhstnam
505            CLOSE(inum)
506         ENDIF
507
508         ! Define the T grid FILE ( nid_T )
509
510         CALL dia_nam( clhstnam, nn_write, 'grid_T' )
511         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
512         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
513            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
514            &          nit000-1, zjulian, rn_Dt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
515         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
516            &           "m", ipk, gdept_1d, nz_T, "down" )
517         !                                                            ! Index of ocean points
518         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
519         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
520         !
521         IF( ln_icebergs ) THEN
522            !
523            !! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
524            !! that routine is called from nemogcm, so do it here immediately before its needed
525            ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
526            CALL mpp_sum( 'diawri', ierror )
527            IF( ierror /= 0 ) THEN
528               CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
529               RETURN
530            ENDIF
531            !
532            !! iceberg vertical coordinate is class number
533            CALL histvert( nid_T, "class", "Iceberg class",      &  ! Vertical grid: class
534               &           "number", nclasses, class_num, nb_T )
535            !
536            !! each class just needs the surface index pattern
537            ndim_bT = 3
538            DO jn = 1,nclasses
539               ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
540            ENDDO
541            !
542         ENDIF
543
544         ! Define the U grid FILE ( nid_U )
545
546         CALL dia_nam( clhstnam, nn_write, 'grid_U' )
547         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
548         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
549            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
550            &          nit000-1, zjulian, rn_Dt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
551         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
552            &           "m", ipk, gdept_1d, nz_U, "down" )
553         !                                                            ! Index of ocean points
554         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
555         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
556
557         ! Define the V grid FILE ( nid_V )
558
559         CALL dia_nam( clhstnam, nn_write, 'grid_V' )                   ! filename
560         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
561         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
562            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
563            &          nit000-1, zjulian, rn_Dt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
564         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
565            &          "m", ipk, gdept_1d, nz_V, "down" )
566         !                                                            ! Index of ocean points
567         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
568         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
569
570         ! Define the W grid FILE ( nid_W )
571
572         CALL dia_nam( clhstnam, nn_write, 'grid_W' )                   ! filename
573         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
574         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
575            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
576            &          nit000-1, zjulian, rn_Dt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
577         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
578            &          "m", ipk, gdepw_1d, nz_W, "down" )
579
580         IF( ln_abl ) THEN 
581         ! Define the ABL grid FILE ( nid_A )
582            CALL dia_nam( clhstnam, nn_write, 'grid_ABL' )
583            IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
584            CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
585               &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
586               &          nit000-1, zjulian, rn_Dt, nh_A, nid_A, domain_id=nidom, snc4chunks=snc4set )
587            CALL histvert( nid_A, "ght_abl", "Vertical T levels",      &  ! Vertical grid: gdept
588               &           "m", ipka, ght_abl(2:jpka), nz_A, "up" )
589            !                                                            ! Index of ocean points
590         ALLOCATE( zw3d_abl(jpi,jpj,ipka) ) 
591         zw3d_abl(:,:,:) = 1._wp 
592         CALL wheneq( jpi*jpj*ipka, zw3d_abl, 1, 1., ndex_A , ndim_A  )      ! volume
593            CALL wheneq( jpi*jpj     , zw3d_abl, 1, 1., ndex_hA, ndim_hA )      ! surface
594         DEALLOCATE(zw3d_abl)
595         ENDIF
596         !
597
598         ! Declare all the output fields as NETCDF variables
599
600         !                                                                                      !!! nid_T : 3D
601         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
602            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
603         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
604            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
605         IF(  .NOT.ln_linssh  ) THEN
606            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t n
607            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
608            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t n
609            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
610            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t n
611            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
612         ENDIF
613         !                                                                                      !!! nid_T : 2D
614         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
615            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
616         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
617            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
618         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
619            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
620         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
621            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
622         CALL histdef( nid_T, "sorunoff", "River runoffs"                      , "Kg/m2/s",   &  ! runoffs
623            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
624         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
625            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
626         IF(  ln_linssh  ) THEN
627            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * ts(:,:,1,jp_tem,Kmm)
628            &                                                                  , "KgC/m2/s",  &  ! sosst_cd
629            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
630            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * ts(:,:,1,jp_sal,Kmm)
631            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
632            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
633         ENDIF
634         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
635            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
636         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
637            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
638         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
639            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
640         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
641            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
642         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
643            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
644         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
645            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
646         !
647         IF( ln_abl ) THEN
648            CALL histdef( nid_A, "t_abl", "Potential Temperature"     , "K"        ,       &  ! t_abl
649               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
650            CALL histdef( nid_A, "q_abl", "Humidity"                  , "kg/kg"    ,       &  ! q_abl
651               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
652            CALL histdef( nid_A, "u_abl", "Atmospheric U-wind   "     , "m/s"        ,     &  ! u_abl
653               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
654            CALL histdef( nid_A, "v_abl", "Atmospheric V-wind   "     , "m/s"    ,         &  ! v_abl
655               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
656            CALL histdef( nid_A, "tke_abl", "Atmospheric TKE   "     , "m2/s2"    ,        &  ! tke_abl
657               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
658            CALL histdef( nid_A, "avm_abl", "Atmospheric turbulent viscosity", "m2/s"   ,  &  ! avm_abl
659               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
660            CALL histdef( nid_A, "avt_abl", "Atmospheric turbulent diffusivity", "m2/s2",  &  ! avt_abl
661               &          jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout ) 
662            CALL histdef( nid_A, "pblh", "Atmospheric boundary layer height "  , "m",      &  ! pblh
663               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout )                 
664#if defined key_si3
665            CALL histdef( nid_A, "oce_frac", "Fraction of open ocean"  , " ",      &  ! ato_i
666               &          jpi, jpj, nh_A,  1  , 1, 1   , -99 , 32, clop, zsto, zout )
667#endif
668            CALL histend( nid_A, snc4chunks=snc4set )
669         ENDIF
670         !
671         IF( ln_icebergs ) THEN
672            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
673               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
674            CALL histdef( nid_T, "calving_heat"        , "calving heat flux"                        , "XXXX"   , &
675               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
676            CALL histdef( nid_T, "berg_floating_melt"  , "Melt rate of icebergs + bits"             , "kg/m2/s", &
677               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
678            CALL histdef( nid_T, "berg_stored_ice"     , "Accumulated ice mass by class"            , "kg"     , &
679               &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
680            IF( ln_bergdia ) THEN
681               CALL histdef( nid_T, "berg_melt"           , "Melt rate of icebergs"                    , "kg/m2/s", &
682                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
683               CALL histdef( nid_T, "berg_buoy_melt"      , "Buoyancy component of iceberg melt rate"  , "kg/m2/s", &
684                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
685               CALL histdef( nid_T, "berg_eros_melt"      , "Erosion component of iceberg melt rate"   , "kg/m2/s", &
686                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
687               CALL histdef( nid_T, "berg_conv_melt"      , "Convective component of iceberg melt rate", "kg/m2/s", &
688                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
689               CALL histdef( nid_T, "berg_virtual_area"   , "Virtual coverage by icebergs"             , "m2"     , &
690                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
691               CALL histdef( nid_T, "bits_src"           , "Mass source of bergy bits"                , "kg/m2/s", &
692                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
693               CALL histdef( nid_T, "bits_melt"          , "Melt rate of bergy bits"                  , "kg/m2/s", &
694                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
695               CALL histdef( nid_T, "bits_mass"          , "Bergy bit density field"                  , "kg/m2"  , &
696                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
697               CALL histdef( nid_T, "berg_mass"           , "Iceberg density field"                    , "kg/m2"  , &
698                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
699               CALL histdef( nid_T, "berg_real_calving"   , "Calving into iceberg class"               , "kg/s"   , &
700                  &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
701            ENDIF
702         ENDIF
703
704         IF( ln_ssr ) THEN
705            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
706               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
707            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
708               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
709            CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
710               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
711         ENDIF
712       
713         clmx ="l_max(only(x))"    ! max index on a period
714!         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
715!            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
716#if defined key_diahth
717         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
718            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
719         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
720            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
721         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
722            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
723         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "J/m2"   ,   & ! htc3
724            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
725#endif
726
727         CALL histend( nid_T, snc4chunks=snc4set )
728
729         !                                                                                      !!! nid_U : 3D
730         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! uu(:,:,:,Kmm)
731            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
732         IF( ln_wave .AND. ln_sdw) THEN
733            CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current"         , "m/s"    ,   &  ! usd
734               &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
735         ENDIF
736         !                                                                                      !!! nid_U : 2D
737         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
738            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
739
740         CALL histend( nid_U, snc4chunks=snc4set )
741
742         !                                                                                      !!! nid_V : 3D
743         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vv(:,:,:,Kmm)
744            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
745         IF( ln_wave .AND. ln_sdw) THEN
746            CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current"    , "m/s"    ,   &  ! vsd
747               &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
748         ENDIF
749         !                                                                                      !!! nid_V : 2D
750         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
751            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
752
753         CALL histend( nid_V, snc4chunks=snc4set )
754
755         !                                                                                      !!! nid_W : 3D
756         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! ww
757            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
758         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
759            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
760         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avm
761            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
762
763         IF( ln_zdfddm ) THEN
764            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
765               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
766         ENDIF
767         
768         IF( ln_wave .AND. ln_sdw) THEN
769            CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current"   , "m/s"    ,   &  ! wsd
770               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
771         ENDIF
772         !                                                                                      !!! nid_W : 2D
773         CALL histend( nid_W, snc4chunks=snc4set )
774
775         IF(lwp) WRITE(numout,*)
776         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
777         IF(ll_print) CALL FLUSH(numout )
778
779      ENDIF
780
781      ! 2. Start writing data
782      ! ---------------------
783
784      ! ndex(1) est utilise ssi l'avant dernier argument est different de
785      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
786      ! donne le nombre d'elements, et ndex la liste des indices a sortir
787
788      IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN
789         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
790         WRITE(numout,*) '~~~~~~ '
791      ENDIF
792
793      IF( .NOT.ln_linssh ) THEN
794         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T  )   ! heat content
795         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T  )   ! salt content
796         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
797         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
798      ELSE
799         CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T  )   ! temperature
800         CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) , ndim_T , ndex_T  )   ! salinity
801         CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) , ndim_hT, ndex_hT )   ! sea surface temperature
802         CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) , ndim_hT, ndex_hT )   ! sea surface salinity
803      ENDIF
804      IF( .NOT.ln_linssh ) THEN
805         zw3d(:,:,:) = ( ( ze3t(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
806         CALL histwrite( nid_T, "vovvle3t", it, ze3t (:,:,:)     , ndim_T , ndex_T  )   ! level thickness
807         CALL histwrite( nid_T, "vovvldep", it, gdept(:,:,:,Kmm) , ndim_T , ndex_T  )   ! t-point depth
808         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
809      ENDIF
810      CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm)          , ndim_hT, ndex_hT )   ! sea surface height
811      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
812      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
813      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
814                                                                                  ! (includes virtual salt flux beneath ice
815                                                                                  ! in linear free surface case)
816      IF( ln_linssh ) THEN
817         zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_tem,Kmm)
818         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
819         zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_sal,Kmm)
820         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
821      ENDIF
822      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
823      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
824      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
825      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
826      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
827      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
828      !
829      IF( ln_abl ) THEN
830         ALLOCATE( zw3d_abl(jpi,jpj,jpka) )
831         IF( ln_mskland )   THEN
832            DO jk=1,jpka
833               zw3d_abl(:,:,jk) = tmask(:,:,1)
834            END DO       
835         ELSE
836            zw3d_abl(:,:,:) = 1._wp     
837         ENDIF       
838         CALL histwrite( nid_A,  "pblh"   , it, pblh(:,:)                  *zw3d_abl(:,:,1     ), ndim_hA, ndex_hA )   ! pblh
839         CALL histwrite( nid_A,  "u_abl"  , it, u_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! u_abl
840         CALL histwrite( nid_A,  "v_abl"  , it, v_abl   (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! v_abl
841         CALL histwrite( nid_A,  "t_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! t_abl
842         CALL histwrite( nid_A,  "q_abl"  , it, tq_abl  (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! q_abl       
843         CALL histwrite( nid_A,  "tke_abl", it, tke_abl (:,:,2:jpka,nt_n  )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! tke_abl
844         CALL histwrite( nid_A,  "avm_abl", it, avm_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avm_abl
845         CALL histwrite( nid_A,  "avt_abl", it, avt_abl (:,:,2:jpka       )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A  )   ! avt_abl
846#if defined key_si3
847         CALL histwrite( nid_A,  "oce_frac"   , it, ato_i(:,:)                                  , ndim_hA, ndex_hA )   ! ato_i
848#endif
849         DEALLOCATE(zw3d_abl)
850      ENDIF
851      !
852      IF( ln_icebergs ) THEN
853         !
854         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT ) 
855         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
856         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT ) 
857         !
858         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
859         !
860         IF( ln_bergdia ) THEN
861            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   ) 
862            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   ) 
863            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   ) 
864            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   ) 
865            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   ) 
866            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   ) 
867            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   ) 
868            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   ) 
869            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   ) 
870            !
871            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
872         ENDIF
873      ENDIF
874
875      IF( ln_ssr ) THEN
876         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
877         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
878         zw2d(:,:) = erp(:,:) * ts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
879         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
880      ENDIF
881!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
882!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
883
884#if defined key_diahth
885      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
886      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
887      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
888      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
889#endif
890
891      CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm)            , ndim_U , ndex_U )    ! i-current
892      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
893
894      CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm)            , ndim_V , ndex_V  )   ! j-current
895      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
896
897      IF( ln_zad_Aimp ) THEN
898         CALL histwrite( nid_W, "vovecrtz", it, ww + wi     , ndim_T, ndex_T )    ! vert. current
899      ELSE
900         CALL histwrite( nid_W, "vovecrtz", it, ww          , ndim_T, ndex_T )    ! vert. current
901      ENDIF
902      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
903      CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
904      IF( ln_zdfddm ) THEN
905         CALL histwrite( nid_W, "voddmavs", it, avs         , ndim_T, ndex_T )    ! S vert. eddy diff. coef.
906      ENDIF
907
908      IF( ln_wave .AND. ln_sdw ) THEN
909         CALL histwrite( nid_U, "sdzocrtx", it, usd         , ndim_U , ndex_U )    ! i-StokesDrift-current
910         CALL histwrite( nid_V, "sdmecrty", it, vsd         , ndim_V , ndex_V )    ! j-StokesDrift-current
911         CALL histwrite( nid_W, "sdvecrtz", it, wsd         , ndim_T , ndex_T )    ! StokesDrift vert. current
912      ENDIF
913
914      ! 3. Close all files
915      ! ---------------------------------------
916      IF( kt == nitend ) THEN
917         CALL histclo( nid_T )
918         CALL histclo( nid_U )
919         CALL histclo( nid_V )
920         CALL histclo( nid_W )
921         IF(ln_abl) CALL histclo( nid_A )
922      ENDIF
923      !
924      IF( ln_timing )   CALL timing_stop('dia_wri')
925      !
926   END SUBROUTINE dia_wri
927#endif
928
929   SUBROUTINE dia_wri_state( Kmm, cdfile_name )
930      !!---------------------------------------------------------------------
931      !!                 ***  ROUTINE dia_wri_state  ***
932      !!       
933      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
934      !!      the instantaneous ocean state and forcing fields.
935      !!        Used to find errors in the initial state or save the last
936      !!      ocean state in case of abnormal end of a simulation
937      !!
938      !! ** Method  :   NetCDF files using ioipsl
939      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
940      !!      File 'output.abort.nc' is created in case of abnormal job end
941      !!----------------------------------------------------------------------
942      INTEGER           , INTENT( in ) ::   Kmm              ! time level index
943      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
944      !!
945      INTEGER :: inum, jk
946      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t       ! 3D workspace
947      !!----------------------------------------------------------------------
948      !
949      IF(lwp) THEN
950         WRITE(numout,*)
951         WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
952         WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
953         WRITE(numout,*) '                and named :', cdfile_name, '...nc'
954      ENDIF 
955      !
956      DO jk = 1, jpk
957         ze3t(:,:,jk) =  e3t(:,:,jk,Kmm)
958      END DO
959      !
960      CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
961      !
962      CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) )    ! now temperature
963      CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) )    ! now salinity
964      CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm)              )    ! sea surface height
965      CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm)                )    ! now i-velocity
966      CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm)                )    ! now j-velocity
967      IF( ln_zad_Aimp ) THEN
968         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww + wi        )    ! now k-velocity
969      ELSE
970         CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww             )    ! now k-velocity
971      ENDIF
972      CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep            )    ! now k-velocity
973      CALL iom_rstput( 0, 0, inum, 'ht'     , ht                 )    ! now water column height
974      !
975      IF ( ln_isf ) THEN
976         IF (ln_isfcav_mlt) THEN
977            CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav          )    ! now k-velocity
978            CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav    )    ! now k-velocity
979            CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav    )    ! now k-velocity
980            CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,wp) )    ! now k-velocity
981            CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,wp) )    ! now k-velocity
982            CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,wp), ktype = jp_i1 )
983         END IF
984         IF (ln_isfpar_mlt) THEN
985            CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,wp) )    ! now k-velocity
986            CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par          )    ! now k-velocity
987            CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par    )    ! now k-velocity
988            CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par    )    ! now k-velocity
989            CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,wp) )    ! now k-velocity
990            CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,wp) )    ! now k-velocity
991            CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,wp), ktype = jp_i1 )
992         END IF
993      END IF
994      !
995      IF( ALLOCATED(ahtu) ) THEN
996         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point
997         CALL iom_rstput( 0, 0, inum,  'ahtv', ahtv              )    ! aht at v-point
998      ENDIF
999      IF( ALLOCATED(ahmt) ) THEN
1000         CALL iom_rstput( 0, 0, inum,  'ahmt', ahmt              )    ! ahmt at u-point
1001         CALL iom_rstput( 0, 0, inum,  'ahmf', ahmf              )    ! ahmf at v-point
1002      ENDIF
1003      CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf         )    ! freshwater budget
1004      CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns         )    ! total heat flux
1005      CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr               )    ! solar heat flux
1006      CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i              )    ! ice fraction
1007      CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
1008      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
1009      IF(  .NOT.ln_linssh  ) THEN             
1010         CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept(:,:,:,Kmm)        )    !  T-cell depth
1011         CALL iom_rstput( 0, 0, inum, 'vovvle3t', ze3t(:,:,:)          )    !  T-cell thickness 
1012      END IF
1013      IF( ln_wave .AND. ln_sdw ) THEN
1014         CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd            )    ! now StokesDrift i-velocity
1015         CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd            )    ! now StokesDrift j-velocity
1016         CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity
1017      ENDIF
1018      IF ( ln_abl ) THEN
1019         CALL iom_rstput ( 0, 0, inum, "uz1_abl",   u_abl(:,:,2,nt_a  ) )   ! now first level i-wind
1020         CALL iom_rstput ( 0, 0, inum, "vz1_abl",   v_abl(:,:,2,nt_a  ) )   ! now first level j-wind
1021         CALL iom_rstput ( 0, 0, inum, "tz1_abl",  tq_abl(:,:,2,nt_a,1) )   ! now first level temperature
1022         CALL iom_rstput ( 0, 0, inum, "qz1_abl",  tq_abl(:,:,2,nt_a,2) )   ! now first level humidity
1023      ENDIF
1024      !
1025      CALL iom_close( inum )
1026      !
1027#if defined key_si3
1028      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid
1029         CALL iom_open( TRIM(cdfile_name)//'_ice', inum, ldwrt = .TRUE., kdlev = jpl, cdcomp = 'ICE' )
1030         CALL ice_wri_state( inum )
1031         CALL iom_close( inum )
1032      ENDIF
1033#endif
1034
1035   END SUBROUTINE dia_wri_state
1036
1037   !!======================================================================
1038END MODULE diawri
Note: See TracBrowser for help on using the repository browser.