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

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

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