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 @ 7508

Last change on this file since 7508 was 7508, checked in by mocavero, 7 years ago

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