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

source: branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 7557

Last change on this file since 7557 was 7557, checked in by emanuelaclementi, 7 years ago

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