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/dev_r8814_surge_modelling_Nemo4/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/UKMO/dev_r8814_surge_modelling_Nemo4/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 9243

Last change on this file since 9243 was 9243, checked in by clne, 6 years ago

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