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

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 6315

Last change on this file since 6315 was 6315, checked in by cetlod, 8 years ago

3.6 stable : updates for XIO2 and iom_put of reference vertical scale factor and some PISCES diagnostics needed for CMIP6, see ticket #1678

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