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

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 7037

Last change on this file since 7037 was 7037, checked in by mocavero, 8 years ago

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