New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diawri.F90 in branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

source: branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 6999

Last change on this file since 6999 was 6999, checked in by kingr, 8 years ago

Fix to typo in previous commmit.

File size: 56.8 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   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   dia_wri       : create the standart output files
23   !!   dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
24   !!----------------------------------------------------------------------
25   USE oce             ! ocean dynamics and tracers
26   USE dom_oce         ! ocean space and time domain
27   USE dynadv, ONLY: ln_dynadv_vec
28   USE zdf_oce         ! ocean vertical physics
29   USE ldftra_oce      ! ocean active tracers: lateral physics
30   USE ldfdyn_oce      ! ocean dynamics: lateral physics
31   USE traldf_iso_grif, ONLY : psix_eiv, psiy_eiv
32   USE sol_oce         ! solver variables
33   USE sbc_oce         ! Surface boundary condition: ocean fields
34   USE sbc_ice         ! Surface boundary condition: ice fields
35   USE icb_oce         ! Icebergs
36   USE icbdia          ! Iceberg budgets
37   USE sbcssr          ! restoring term toward SST/SSS climatology
38   USE phycst          ! physical constants
39   USE zdfmxl          ! mixed layer
40   USE dianam          ! build name of file (routine)
41   USE zdfddm          ! vertical  physics: double diffusion
42   USE diahth          ! thermocline diagnostics
43   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
44   USE in_out_manager  ! I/O manager
45   USE diadimg         ! dimg direct access file format output
46   USE diatmb          ! Top,middle,bottom output
47   USE dia25h          ! 25h Mean output
48   USE iom
49   USE ioipsl
50   USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities     
51
52#if defined key_lim2
53   USE limwri_2 
54#elif defined key_lim3
55   USE limwri 
56#endif
57   USE lib_mpp         ! MPP library
58   USE timing          ! preformance summary
59   USE 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) = 0.5_wp * 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#if defined key_dynspg_ts
211      CALL iom_put(  "ubar", un_adv(:,:)      )    ! barotropic i-current
212#else
213      CALL iom_put(  "ubar", un_b(:,:)        )    ! barotropic i-current
214#endif
215     
216      CALL iom_put( "voce", vn(:,:,:)         )    ! 3D j-current
217      CALL iom_put(  "ssv", vn(:,:,1)         )    ! surface j-current
218      IF ( iom_use("sbv") ) THEN
219         DO jj = 1, jpj
220            DO ji = 1, jpi
221               jkbot = mbkv(ji,jj)
222               z2d(ji,jj) = vn(ji,jj,jkbot)
223            END DO
224         END DO
225         CALL iom_put( "sbv", z2d )                ! bottom j-current
226      ENDIF
227#if defined key_dynspg_ts
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 * e12t(:,:)
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 ) / e1u(ji,jj) + ( zztmp - tsn(ji-1,jj  ,1,jp_tem) ) / e1u(ji-1,jj  )
253               zztmpy     = ( tsn(ji  ,jj+1,1,jp_tem) - zztmp ) / e2v(ji,jj) + ( zztmp - tsn(ji  ,jj-1,1,jp_tem) ) / 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      !!
415      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
416      !!
417      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
418      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
419      INTEGER  ::   inum = 11                                ! temporary logical unit
420      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
421      INTEGER  ::   ierr                                     ! error code return from allocation
422      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
423      INTEGER  ::   jn, ierror                               ! local integers
424      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt           ! local scalars
425      !!
426      REAL(wp), POINTER, DIMENSION(:,:)   :: zw2d       ! 2D workspace
427      REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d       ! 3D workspace
428      !!----------------------------------------------------------------------
429      !
430      IF( nn_timing == 1 )   CALL timing_start('dia_wri')
431      !
432      CALL wrk_alloc( jpi , jpj      , zw2d )
433      IF ( ln_traldf_gdia .OR. lk_vvl )  call wrk_alloc( jpi , jpj , jpk  , zw3d )
434      !
435      ! Output the initial state and forcings
436      IF( ninist == 1 ) THEN                       
437         CALL dia_wri_state( 'output.init', kt )
438         ninist = 0
439      ENDIF
440      !
441      ! 0. Initialisation
442      ! -----------------
443
444      ! local variable for debugging
445      ll_print = .FALSE.
446      ll_print = ll_print .AND. lwp
447
448      ! Define frequency of output and means
449      zdt = rdt
450      IF( nacc == 1 ) zdt = rdtmin
451      clop = "x"         ! no use of the mask value (require less cpu time, and otherwise the model crashes)
452#if defined key_diainstant
453      zsto = nwrite * zdt
454      clop = "inst("//TRIM(clop)//")"
455#else
456      zsto=zdt
457      clop = "ave("//TRIM(clop)//")"
458#endif
459      zout = nwrite * zdt
460      zmax = ( nitend - nit000 + 1 ) * zdt
461
462      ! Define indices of the horizontal output zoom and vertical limit storage
463      iimi = 1      ;      iima = jpi
464      ijmi = 1      ;      ijma = jpj
465      ipk = jpk
466
467      ! define time axis
468      it = kt
469      itmod = kt - nit000 + 1
470
471
472      ! 1. Define NETCDF files and fields at beginning of first time step
473      ! -----------------------------------------------------------------
474
475      IF( kt == nit000 ) THEN
476
477         ! Define the NETCDF files (one per grid)
478
479         ! Compute julian date from starting date of the run
480         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
481         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
482         IF(lwp)WRITE(numout,*)
483         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
484            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
485         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
486                                 ' limit storage in depth = ', ipk
487
488         ! WRITE root name in date.file for use by postpro
489         IF(lwp) THEN
490            CALL dia_nam( clhstnam, nwrite,' ' )
491            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
492            WRITE(inum,*) clhstnam
493            CLOSE(inum)
494         ENDIF
495
496         ! Define the T grid FILE ( nid_T )
497
498         CALL dia_nam( clhstnam, nwrite, 'grid_T' )
499         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
500         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
501            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
502            &          nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
503         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
504            &           "m", ipk, gdept_1d, nz_T, "down" )
505         !                                                            ! Index of ocean points
506         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
507         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
508         !
509         IF( ln_icebergs ) THEN
510            !
511            !! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
512            !! that routine is called from nemogcm, so do it here immediately before its needed
513            ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
514            IF( lk_mpp )   CALL mpp_sum( ierror )
515            IF( ierror /= 0 ) THEN
516               CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
517               RETURN
518            ENDIF
519            !
520            !! iceberg vertical coordinate is class number
521            CALL histvert( nid_T, "class", "Iceberg class",      &  ! Vertical grid: class
522               &           "number", nclasses, class_num, nb_T )
523            !
524            !! each class just needs the surface index pattern
525            ndim_bT = 3
526            DO jn = 1,nclasses
527               ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
528            ENDDO
529            !
530         ENDIF
531
532         ! Define the U grid FILE ( nid_U )
533
534         CALL dia_nam( clhstnam, nwrite, 'grid_U' )
535         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
536         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
537            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
538            &          nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
539         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
540            &           "m", ipk, gdept_1d, nz_U, "down" )
541         !                                                            ! Index of ocean points
542         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
543         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
544
545         ! Define the V grid FILE ( nid_V )
546
547         CALL dia_nam( clhstnam, nwrite, 'grid_V' )                   ! filename
548         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
549         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
550            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
551            &          nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
552         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
553            &          "m", ipk, gdept_1d, nz_V, "down" )
554         !                                                            ! Index of ocean points
555         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
556         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
557
558         ! Define the W grid FILE ( nid_W )
559
560         CALL dia_nam( clhstnam, nwrite, 'grid_W' )                   ! filename
561         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
562         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
563            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
564            &          nit000-1, zjulian, zdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
565         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
566            &          "m", ipk, gdepw_1d, nz_W, "down" )
567
568
569         ! Declare all the output fields as NETCDF variables
570
571         !                                                                                      !!! nid_T : 3D
572         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
573            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
574         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
575            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
576         IF(  lk_vvl  ) THEN
577            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n
578            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
579            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t_n
580            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
581            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t_n
582            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
583         ENDIF
584         !                                                                                      !!! nid_T : 2D
585         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
586            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
587         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
588            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
589         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
590            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
591         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
592            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
593         CALL histdef( nid_T, "sorunoff", "River runoffs"                      , "Kg/m2/s",   &  ! runoffs
594            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
595         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
596            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
597         IF(  .NOT. lk_vvl  ) THEN
598            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem)
599            &                                                                  , "KgC/m2/s",  &  ! sosst_cd
600            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
601            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * tsn(:,:,1,jp_sal)
602            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
603            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
604         ENDIF
605         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
606            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
607         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
608            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
609         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
610            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
611         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
612            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
613         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
614            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
615         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
616            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
617!
618         IF( ln_icebergs ) THEN
619            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
620               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
621            CALL histdef( nid_T, "calving_heat"        , "calving heat flux"                        , "XXXX"   , &
622               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
623            CALL histdef( nid_T, "berg_floating_melt"  , "Melt rate of icebergs + bits"             , "kg/m2/s", &
624               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
625            CALL histdef( nid_T, "berg_stored_ice"     , "Accumulated ice mass by class"            , "kg"     , &
626               &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
627            IF( ln_bergdia ) THEN
628               CALL histdef( nid_T, "berg_melt"           , "Melt rate of icebergs"                    , "kg/m2/s", &
629                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
630               CALL histdef( nid_T, "berg_buoy_melt"      , "Buoyancy component of iceberg melt rate"  , "kg/m2/s", &
631                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
632               CALL histdef( nid_T, "berg_eros_melt"      , "Erosion component of iceberg melt rate"   , "kg/m2/s", &
633                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
634               CALL histdef( nid_T, "berg_conv_melt"      , "Convective component of iceberg melt rate", "kg/m2/s", &
635                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
636               CALL histdef( nid_T, "berg_virtual_area"   , "Virtual coverage by icebergs"             , "m2"     , &
637                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
638               CALL histdef( nid_T, "bits_src"           , "Mass source of bergy bits"                , "kg/m2/s", &
639                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
640               CALL histdef( nid_T, "bits_melt"          , "Melt rate of bergy bits"                  , "kg/m2/s", &
641                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
642               CALL histdef( nid_T, "bits_mass"          , "Bergy bit density field"                  , "kg/m2"  , &
643                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
644               CALL histdef( nid_T, "berg_mass"           , "Iceberg density field"                    , "kg/m2"  , &
645                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
646               CALL histdef( nid_T, "berg_real_calving"   , "Calving into iceberg class"               , "kg/s"   , &
647                  &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
648            ENDIF
649         ENDIF
650
651         IF( .NOT. ln_cpl ) THEN
652            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
653               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
654            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
655               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
656            CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
657               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
658         ENDIF
659
660         IF( ln_cpl .AND. nn_ice <= 1 ) THEN
661            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
662               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
663            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
664               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
665            CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping"         , "Kg/m2/s",   &  ! erp * sn
666               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
667         ENDIF
668         
669         clmx ="l_max(only(x))"    ! max index on a period
670         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
671            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
672#if defined key_diahth
673         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
674            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
675         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
676            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
677         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
678            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
679         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "W"      ,   & ! htc3
680            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
681#endif
682
683         IF( ln_cpl .AND. nn_ice == 2 ) THEN
684            CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature"            , "K"      ,   &  ! tn_ice
685               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
686            CALL histdef( nid_T,"soicealb" , "Ice Albedo"                         , "[0,1]"  ,   &  ! alb_ice
687               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
688         ENDIF
689
690         CALL histend( nid_T, snc4chunks=snc4set )
691
692         !                                                                                      !!! nid_U : 3D
693         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
694            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
695         IF( ln_traldf_gdia ) THEN
696            CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
697                 &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
698         ELSE
699#if defined key_diaeiv
700            CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current"                  , "m/s"    ,   &  ! u_eiv
701            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
702#endif
703         END IF
704         !                                                                                      !!! nid_U : 2D
705         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
706            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
707
708         CALL histend( nid_U, snc4chunks=snc4set )
709
710         !                                                                                      !!! nid_V : 3D
711         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
712            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
713         IF( ln_traldf_gdia ) THEN
714            CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
715                 &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
716         ELSE 
717#if defined key_diaeiv
718            CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current"             , "m/s"    ,   &  ! v_eiv
719            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
720#endif
721         END IF
722         !                                                                                      !!! nid_V : 2D
723         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
724            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
725
726         CALL histend( nid_V, snc4chunks=snc4set )
727
728         !                                                                                      !!! nid_W : 3D
729         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn
730            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
731         IF( ln_traldf_gdia ) THEN
732            CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
733                 &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
734         ELSE
735#if defined key_diaeiv
736            CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity"              , "m/s"    ,   &  ! w_eiv
737                 &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
738#endif
739         END IF
740         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
741            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
742         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avmu
743            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
744
745         IF( lk_zdfddm ) THEN
746            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
747               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
748         ENDIF
749         !                                                                                      !!! nid_W : 2D
750#if defined key_traldf_c2d
751         CALL histdef( nid_W, "soleahtw", "lateral eddy diffusivity"           , "m2/s"   ,   &  ! ahtw
752            &          jpi, jpj, nh_W, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
753# if defined key_traldf_eiv 
754            CALL histdef( nid_W, "soleaeiw", "eddy induced vel. coeff. at w-point", "m2/s",   &  ! aeiw
755               &       jpi, jpj, nh_W, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
756# endif
757#endif
758
759         CALL histend( nid_W, snc4chunks=snc4set )
760
761         IF(lwp) WRITE(numout,*)
762         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
763         IF(ll_print) CALL FLUSH(numout )
764
765      ENDIF
766
767      ! 2. Start writing data
768      ! ---------------------
769
770      ! ndex(1) est utilise ssi l'avant dernier argument est different de
771      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
772      ! donne le nombre d'elements, et ndex la liste des indices a sortir
773
774      IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
775         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
776         WRITE(numout,*) '~~~~~~ '
777      ENDIF
778
779      IF( lk_vvl ) THEN
780         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
781         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
782         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
783         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * fse3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
784      ELSE
785         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature
786         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity
787         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature
788         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
789      ENDIF
790      IF( lk_vvl ) THEN
791         zw3d(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
792         CALL histwrite( nid_T, "vovvle3t", it, fse3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
793         CALL histwrite( nid_T, "vovvldep", it, fsdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
794         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
795      ENDIF
796      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
797      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
798      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
799      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
800                                                                                  ! (includes virtual salt flux beneath ice
801                                                                                  ! in linear free surface case)
802      IF( .NOT. lk_vvl ) THEN
803         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
804         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
805         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
806         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
807      ENDIF
808      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
809      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
810      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
811      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
812      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
813      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
814!
815      IF( ln_icebergs ) THEN
816         !
817         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT ) 
818         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
819         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT ) 
820         !
821         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
822         !
823         IF( ln_bergdia ) THEN
824            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   ) 
825            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   ) 
826            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   ) 
827            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   ) 
828            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   ) 
829            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   ) 
830            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   ) 
831            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   ) 
832            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   ) 
833            !
834            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
835         ENDIF
836      ENDIF
837
838      IF( .NOT. ln_cpl ) THEN
839         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
840         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
841         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
842         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
843      ENDIF
844      IF( ln_cpl .AND. nn_ice <= 1 ) THEN
845         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
846         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
847         IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
848         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
849      ENDIF
850!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
851!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
852
853#if defined key_diahth
854      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
855      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
856      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
857      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
858#endif
859
860      IF( ln_cpl .AND. nn_ice == 2 ) THEN
861         CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT )   ! surf. ice temperature
862         CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT )   ! ice albedo
863      ENDIF
864
865      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
866      IF( ln_traldf_gdia ) THEN
867         IF (.not. ALLOCATED(psix_eiv))THEN
868            ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr )
869            IF( lk_mpp   )   CALL mpp_sum ( ierr )
870            IF( ierr > 0 )   CALL ctl_stop('STOP', 'diawri: unable to allocate psi{x,y}_eiv')
871            psix_eiv(:,:,:) = 0.0_wp
872            psiy_eiv(:,:,:) = 0.0_wp
873         ENDIF
874         DO jk=1,jpkm1
875            zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk)  ! u_eiv = -dpsix/dz
876         END DO
877         zw3d(:,:,jpk) = 0._wp
878         CALL histwrite( nid_U, "vozoeivu", it, zw3d, ndim_U , ndex_U )           ! i-eiv current
879      ELSE
880#if defined key_diaeiv
881         CALL histwrite( nid_U, "vozoeivu", it, u_eiv, ndim_U , ndex_U )          ! i-eiv current
882#endif
883      ENDIF
884      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
885
886      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
887      IF( ln_traldf_gdia ) THEN
888         DO jk=1,jpk-1
889            zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk)  ! v_eiv = -dpsiy/dz
890         END DO
891         zw3d(:,:,jpk) = 0._wp
892         CALL histwrite( nid_V, "vomeeivv", it, zw3d, ndim_V , ndex_V )           ! j-eiv current
893      ELSE
894#if defined key_diaeiv
895         CALL histwrite( nid_V, "vomeeivv", it, v_eiv, ndim_V , ndex_V )          ! j-eiv current
896#endif
897      ENDIF
898      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
899
900      CALL histwrite( nid_W, "vovecrtz", it, wn             , ndim_T, ndex_T )    ! vert. current
901      IF( ln_traldf_gdia ) THEN
902         DO jk=1,jpk-1
903            DO jj = 2, jpjm1
904               DO ji = fs_2, fs_jpim1  ! vector opt.
905                  zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2v(ji,jj) + &
906                       &    (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1u(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx
907               END DO
908            END DO
909         END DO
910         zw3d(:,:,jpk) = 0._wp
911         CALL histwrite( nid_W, "voveeivw", it, zw3d          , ndim_T, ndex_T )    ! vert. eiv current
912      ELSE
913#   if defined key_diaeiv
914         CALL histwrite( nid_W, "voveeivw", it, w_eiv          , ndim_T, ndex_T )    ! vert. eiv current
915#   endif
916      ENDIF
917      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
918      CALL histwrite( nid_W, "votkeavm", it, avmu           , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
919      IF( lk_zdfddm ) THEN
920         CALL histwrite( nid_W, "voddmavs", it, fsavs(:,:,:), ndim_T, ndex_T )    ! S vert. eddy diff. coef.
921      ENDIF
922#if defined key_traldf_c2d
923      CALL histwrite( nid_W, "soleahtw", it, ahtw          , ndim_hT, ndex_hT )   ! lateral eddy diff. coef.
924# if defined key_traldf_eiv
925         CALL histwrite( nid_W, "soleaeiw", it, aeiw       , ndim_hT, ndex_hT )   ! EIV coefficient at w-point
926# endif
927#endif
928
929      ! 3. Close all files
930      ! ---------------------------------------
931      IF( kt == nitend ) THEN
932         CALL histclo( nid_T )
933         CALL histclo( nid_U )
934         CALL histclo( nid_V )
935         CALL histclo( nid_W )
936      ENDIF
937      !
938      CALL wrk_dealloc( jpi , jpj      , zw2d )
939      IF ( ln_traldf_gdia .OR. lk_vvl )  call wrk_dealloc( jpi , jpj , jpk  , zw3d )
940      !
941      IF( nn_timing == 1 )   CALL timing_stop('dia_wri')
942      !
943   END SUBROUTINE dia_wri
944# endif
945
946#endif
947
948   SUBROUTINE dia_wri_state( cdfile_name, kt )
949      !!---------------------------------------------------------------------
950      !!                 ***  ROUTINE dia_wri_state  ***
951      !!       
952      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
953      !!      the instantaneous ocean state and forcing fields.
954      !!        Used to find errors in the initial state or save the last
955      !!      ocean state in case of abnormal end of a simulation
956      !!
957      !! ** Method  :   NetCDF files using ioipsl
958      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
959      !!      File 'output.abort.nc' is created in case of abnormal job end
960      !!----------------------------------------------------------------------
961      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
962      INTEGER           , INTENT( in ) ::   kt               ! ocean time-step index
963      !!
964      CHARACTER (len=32) :: clname
965      CHARACTER (len=40) :: clop
966      INTEGER  ::   id_i , nz_i, nh_i       
967      INTEGER, DIMENSION(1) ::   idex             ! local workspace
968      REAL(wp) ::   zsto, zout, zmax, zjulian, zdt
969      !!----------------------------------------------------------------------
970      !
971!     IF( nn_timing == 1 )   CALL timing_start('dia_wri_state') ! not sure this works for routines not called in first timestep
972
973      ! 0. Initialisation
974      ! -----------------
975
976      ! Define name, frequency of output and means
977      clname = cdfile_name
978      IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
979      zdt  = rdt
980      zsto = rdt
981      clop = "inst(x)"           ! no use of the mask value (require less cpu time)
982      zout = rdt
983      zmax = ( nitend - nit000 + 1 ) * zdt
984
985      IF(lwp) WRITE(numout,*)
986      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
987      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
988      IF(lwp) WRITE(numout,*) '                and named :', clname, '.nc'
989
990
991      ! 1. Define NETCDF files and fields at beginning of first time step
992      ! -----------------------------------------------------------------
993
994      ! Compute julian date from starting date of the run
995      CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )         ! time axis
996      zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
997      CALL histbeg( clname, jpi, glamt, jpj, gphit,   &
998          1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
999      CALL histvert( id_i, "deptht", "Vertical T levels",   &    ! Vertical grid : gdept
1000          "m", jpk, gdept_1d, nz_i, "down")
1001
1002      ! Declare all the output fields as NetCDF variables
1003
1004      CALL histdef( id_i, "vosaline", "Salinity"              , "PSU"    ,   &   ! salinity
1005         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1006      CALL histdef( id_i, "votemper", "Temperature"           , "C"      ,   &   ! temperature
1007         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1008      CALL histdef( id_i, "sossheig", "Sea Surface Height"    , "m"      ,   &  ! ssh
1009         &          jpi, jpj, nh_i, 1  , 1, 1  , nz_i, 32, clop, zsto, zout )
1010      CALL histdef( id_i, "vozocrtx", "Zonal Current"         , "m/s"    ,   &   ! zonal current
1011         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1012      CALL histdef( id_i, "vomecrty", "Meridional Current"    , "m/s"    ,   &   ! meridonal current
1013         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
1014      CALL histdef( id_i, "vovecrtz", "Vertical Velocity"     , "m/s"    ,   &   ! vertical current
1015         &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout ) 
1016      CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S",   &   ! net freshwater
1017         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1018      CALL histdef( id_i, "sohefldo", "Net Downward Heat Flux", "W/m2"   ,   &   ! net heat flux
1019         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1020      CALL histdef( id_i, "soshfldo", "Shortwave Radiation"   , "W/m2"   ,   &   ! solar flux
1021         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1022      CALL histdef( id_i, "soicecov", "Ice fraction"          , "[0,1]"  ,   &   ! fr_i
1023         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1024      CALL histdef( id_i, "sozotaux", "Zonal Wind Stress"     , "N/m2"   ,   &   ! i-wind stress
1025         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1026      CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2"   ,   &   ! j-wind stress
1027         &          jpi, jpj, nh_i, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
1028      IF( lk_vvl ) THEN
1029         CALL histdef( id_i, "vovvldep", "T point depth"         , "m"      ,   &   ! t-point depth
1030            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1031         CALL histdef( id_i, "vovvle3t", "T point thickness"         , "m"      ,   &   ! t-point depth
1032            &          jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
1033      END IF
1034
1035#if defined key_lim2
1036      CALL lim_wri_state_2( kt, id_i, nh_i )
1037#elif defined key_lim3
1038      CALL lim_wri_state( kt, id_i, nh_i )
1039#else
1040      CALL histend( id_i, snc4chunks=snc4set )
1041#endif
1042
1043      ! 2. Start writing data
1044      ! ---------------------
1045      ! idex(1) est utilise ssi l'avant dernier argument est diffferent de
1046      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
1047      ! donne le nombre d'elements, et idex la liste des indices a sortir
1048      idex(1) = 1   ! init to avoid compil warning
1049
1050      ! Write all fields on T grid
1051      CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpi*jpj*jpk, idex )    ! now temperature
1052      CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpi*jpj*jpk, idex )    ! now salinity
1053      CALL histwrite( id_i, "sossheig", kt, sshn             , jpi*jpj    , idex )    ! sea surface height
1054      CALL histwrite( id_i, "vozocrtx", kt, un               , jpi*jpj*jpk, idex )    ! now i-velocity
1055      CALL histwrite( id_i, "vomecrty", kt, vn               , jpi*jpj*jpk, idex )    ! now j-velocity
1056      CALL histwrite( id_i, "vovecrtz", kt, wn               , jpi*jpj*jpk, idex )    ! now k-velocity
1057      CALL histwrite( id_i, "sowaflup", kt, (emp-rnf )       , jpi*jpj    , idex )    ! freshwater budget
1058      CALL histwrite( id_i, "sohefldo", kt, qsr + qns        , jpi*jpj    , idex )    ! total heat flux
1059      CALL histwrite( id_i, "soshfldo", kt, qsr              , jpi*jpj    , idex )    ! solar heat flux
1060      CALL histwrite( id_i, "soicecov", kt, fr_i             , jpi*jpj    , idex )    ! ice fraction
1061      CALL histwrite( id_i, "sozotaux", kt, utau             , jpi*jpj    , idex )    ! i-wind stress
1062      CALL histwrite( id_i, "sometauy", kt, vtau             , jpi*jpj    , idex )    ! j-wind stress
1063      IF( lk_vvl ) THEN
1064         CALL histwrite( id_i, "vovvldep", kt, fsdept_n(:,:,:), jpi*jpj*jpk, idex )!  T-cell depth       
1065         CALL histwrite( id_i, "vovvle3t", kt, fse3t_n (:,:,:), jpi*jpj*jpk, idex )!  T-cell thickness 
1066      END IF
1067
1068      ! 3. Close the file
1069      ! -----------------
1070      CALL histclo( id_i )
1071#if ! defined key_iomput && ! defined key_dimgout
1072      IF( ninist /= 1  ) THEN
1073         CALL histclo( nid_T )
1074         CALL histclo( nid_U )
1075         CALL histclo( nid_V )
1076         CALL histclo( nid_W )
1077      ENDIF
1078#endif
1079       
1080!     IF( nn_timing == 1 )   CALL timing_stop('dia_wri_state') ! not sure this works for routines not called in first timestep
1081      !
1082
1083   END SUBROUTINE dia_wri_state
1084   !!======================================================================
1085END MODULE diawri
Note: See TracBrowser for help on using the repository browser.