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/releases/r4.0/r4.0-HEAD/src/OCE/DIA – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/DIA/diawri.F90

Last change on this file was 15593, checked in by cetlod, 2 years ago

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