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

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

source: branches/UKMO/dev_r5518_GO6_package_OMP/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 9176

Last change on this file since 9176 was 9176, checked in by andmirek, 6 years ago

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