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

source: branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 6035

Last change on this file since 6035 was 6035, checked in by timgraham, 8 years ago

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