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/2019/dev_ASINTER-01-05_merged/src/OCE/DIA – NEMO

source: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/DIA/diawri.F90 @ 12165

Last change on this file since 12165 was 11874, checked in by gsamson, 4 years ago

dev_r11265_ABL: output surface now ABL fields in output.init.nc file when nn_istate=1 (#2131)

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