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 branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/UKMO/ROMS_WAD_7832/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 8870

Last change on this file since 8870 was 8870, checked in by deazer, 6 years ago

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