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

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

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

Last change on this file since 8561 was 8561, checked in by jgraham, 7 years ago

Updates for operational diagnostics:
25h mean diagnostics - bottom temperature (and insitu temp)
Operational foam diagnostics - diaopfoam and DIU routines added.

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