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

source: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 5989

Last change on this file since 5989 was 5989, checked in by deazer, 8 years ago

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