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

source: branches/2015/dev_r5056_CMCC4_simplification/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 5825

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