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 @ 6005

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

Merged diurnal sst branch

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