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.3_GM_rossby_radius_ramp/src/OCE/DIA – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.3_GM_rossby_radius_ramp/src/OCE/DIA/diawri.F90

Last change on this file was 15770, checked in by davestorkey, 2 years ago

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