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

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

source: branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 11134

Last change on this file since 11134 was 11134, checked in by jcastill, 5 years ago

Full set of changes as in the original branch

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