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/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DIA – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.2_ENHANCE-02_ISF_nemo/src/OCE/DIA/diawri.F90 @ 12708

Last change on this file since 12708 was 12708, checked in by mathiot, 4 years ago

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