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

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

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