source: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 8059

Last change on this file since 8059 was 8059, checked in by jgraham, 3 years ago

Merged branches required for AMM15 simulations, see ticket #1904.
Merged branches include:
branches/UKMO/CO6_KD490
branches/UKMO/CO6_Restartable_Tidal_Analysis
branches/UKMO/AMM15_v3_6_STABLE

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