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

source: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 5038

Last change on this file since 5038 was 5038, checked in by jamesharle, 9 years ago

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