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 NEMO/branches/UKMO/NEMO_4.0.1_GO8_package/src/OCE/DIA – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.1_GO8_package/src/OCE/DIA/diawri.F90 @ 14335

Last change on this file since 14335 was 12574, checked in by cguiavarch, 4 years ago

Add Equation of State tag to salinity and temperature outputs and cell_methods attribute for mean_nemo (mean_nemo_wrapper not needed anymore)
File size: 50.7 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 phycst         ! physical constants
30   USE dianam         ! build name of file (routine)
31   USE diahth         ! thermocline diagnostics
32   USE dynadv   , ONLY: ln_dynadv_vec
33   USE icb_oce        ! Icebergs
34   USE icbdia         ! Iceberg budgets
35   USE ldftra         ! lateral physics: eddy diffusivity coef.
36   USE ldfdyn         ! lateral physics: eddy viscosity   coef.
37   USE sbc_oce        ! Surface boundary condition: ocean fields
38   USE sbc_ice        ! Surface boundary condition: ice fields
39   USE sbcssr         ! restoring term toward SST/SSS climatology
40   USE sbcwave        ! wave parameters
41   USE wet_dry        ! wetting and drying
42   USE zdf_oce        ! ocean vertical physics
43   USE zdfdrg         ! ocean vertical physics: top/bottom friction
44   USE zdfmxl         ! mixed layer
45   !
46   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
47   USE in_out_manager ! I/O manager
48   USE diatmb         ! Top,middle,bottom output
49   USE dia25h         ! 25h Mean output
50   USE iom            !
51   USE ioipsl         !
52   USE eosbn2
53#if defined key_si3
54   USE ice 
55   USE icewri 
56#endif
57   USE lib_mpp         ! MPP library
58   USE timing          ! preformance summary
59   USE diurnal_bulk    ! diurnal warm layer
60   USE cool_skin       ! Cool skin
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 "vectopt_loop_substitute.h90"
81   !!----------------------------------------------------------------------
82   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
83   !! $Id$
84   !! Software governed by the CeCILL license (see ./LICENSE)
85   !!----------------------------------------------------------------------
86CONTAINS
87
88#if defined key_iomput
89   !!----------------------------------------------------------------------
90   !!   'key_iomput'                                        use IOM library
91   !!----------------------------------------------------------------------
92   INTEGER FUNCTION dia_wri_alloc()
93      !
94      dia_wri_alloc = 0
95      !
96   END FUNCTION dia_wri_alloc
97
98   
99   SUBROUTINE dia_wri( kt )
100      !!---------------------------------------------------------------------
101      !!                  ***  ROUTINE dia_wri  ***
102      !!                   
103      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
104      !!      NETCDF format is used by default
105      !!
106      !! ** Method  :  use iom_put
107      !!----------------------------------------------------------------------
108      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
109      !!
110      INTEGER ::   ji, jj, jk       ! dummy loop indices
111      INTEGER ::   ikbot            ! local integer
112      REAL(wp)::   zztmp , zztmpx   ! local scalar
113      REAL(wp)::   zztmp2, zztmpy   !   -      -
114      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace
115      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace
116      CHARACTER(len=4),SAVE :: ttype , stype           ! temperature and salinity type
117      !!----------------------------------------------------------------------
118      !
119      IF( kt == nit000 ) THEN
120         IF( ln_TEOS10 ) THEN
121            IF ( iom_use("toce_pot") .OR. iom_use("soce_pra") .OR. iom_use("sst_pot") .OR. iom_use("sss_pra") &
122                  & .OR. iom_use("sbt_pot") .OR. iom_use("sbs_pra") .OR. iom_use("sstgrad_pot") .OR. iom_use("sstgrad2_pot") &
123                  & .OR. iom_use("tosmint_pot") .OR. iom_use("somint_pra"))  THEN
124               CALL ctl_stop( 'diawri: potential temperature and practical salinity not available with ln_TEOS10' )
125            ELSE
126               ttype='con' ; stype='abs'   ! teos-10 using conservative temperature and absolute salinity
127            ENDIF
128         ELSE IF( ln_EOS80  ) THEN
129            IF ( iom_use("toce_con") .OR. iom_use("soce_abs") .OR. iom_use("sst_con") .OR. iom_use("sss_abs") &
130                  & .OR. iom_use("sbt_con") .OR. iom_use("sbs_abs") .OR. iom_use("sstgrad_con") .OR. iom_use("sstgrad2_con") &
131                  & .OR. iom_use("tosmint_con") .OR. iom_use("somint_abs"))  THEN
132               CALL ctl_stop( 'diawri: conservative temperature and absolute salinity not available with ln_EOS80' )
133            ELSE
134               ttype='pot' ; stype='pra'   ! eos-80 using potential temperature and practical salinity
135            ENDIF
136         ELSE IF ( ln_SEOS) THEN
137            ttype='seos' ; stype='seos' ! seos using Simplified Equation of state
138         ENDIF
139      ENDIF
140
141      IF( ln_timing )   CALL timing_start('dia_wri')
142      !
143      ! Output the initial state and forcings
144      IF( ninist == 1 ) THEN                       
145         CALL dia_wri_state( 'output.init' )
146         ninist = 0
147      ENDIF
148
149      ! Output of initial vertical scale factor
150      CALL iom_put("e3t_0", e3t_0(:,:,:) )
151      CALL iom_put("e3u_0", e3u_0(:,:,:) )
152      CALL iom_put("e3v_0", e3v_0(:,:,:) )
153      !
154      CALL iom_put( "e3t" , e3t_n(:,:,:) )
155      CALL iom_put( "e3u" , e3u_n(:,:,:) )
156      CALL iom_put( "e3v" , e3v_n(:,:,:) )
157      CALL iom_put( "e3w" , e3w_n(:,:,:) )
158      IF( iom_use("e3tdef") )   &
159         CALL iom_put( "e3tdef"  , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
160
161      IF( ll_wd ) THEN
162         CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) )   ! sea surface height (brought back to the reference used for wetting and drying)
163      ELSE
164         CALL iom_put( "ssh" , sshn )              ! sea surface height
165      ENDIF
166
167      IF( iom_use("wetdep") )   &                  ! wet depth
168         CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) )
169     
170      CALL iom_put( "toce_"//ttype, tsn(:,:,:,jp_tem) )    ! 3D temperature
171      CALL iom_put(  "sst_"//ttype, tsn(:,:,1,jp_tem) )    ! surface temperature
172      IF ( iom_use("sbt_"//ttype) ) THEN
173         DO jj = 1, jpj
174            DO ji = 1, jpi
175               ikbot = mbkt(ji,jj)
176               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem)
177            END DO
178         END DO
179         CALL iom_put( "sbt_"//ttype, z2d )                ! bottom temperature
180      ENDIF
181     
182      CALL iom_put( "soce_"//stype, tsn(:,:,:,jp_sal) )    ! 3D salinity
183      CALL iom_put(  "sss_"//stype, tsn(:,:,1,jp_sal) )    ! surface salinity
184      IF ( iom_use("sbs_"//stype) ) THEN
185         DO jj = 1, jpj
186            DO ji = 1, jpi
187               ikbot = mbkt(ji,jj)
188               z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal)
189            END DO
190         END DO
191         CALL iom_put( "sbs_"//stype, z2d )                ! bottom salinity
192      ENDIF
193
194      IF ( iom_use("taubot") ) THEN                ! bottom stress
195         zztmp = rau0 * 0.25
196         z2d(:,:) = 0._wp
197         DO jj = 2, jpjm1
198            DO ji = fs_2, fs_jpim1   ! vector opt.
199               zztmp2 = (  ( rCdU_bot(ji+1,jj)+rCdU_bot(ji  ,jj) ) * un(ji  ,jj,mbku(ji  ,jj))  )**2   &
200                  &   + (  ( rCdU_bot(ji  ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj))  )**2   &
201                  &   + (  ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj  ) ) * vn(ji,jj  ,mbkv(ji,jj  ))  )**2   &
202                  &   + (  ( rCdU_bot(ji,jj  )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1))  )**2
203               z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1) 
204               !
205            END DO
206         END DO
207         CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
208         CALL iom_put( "taubot", z2d )           
209      ENDIF
210         
211      CALL iom_put( "uoce", un(:,:,:) )            ! 3D i-current
212      CALL iom_put(  "ssu", un(:,:,1) )            ! surface i-current
213      IF ( iom_use("sbu") ) THEN
214         DO jj = 1, jpj
215            DO ji = 1, jpi
216               ikbot = mbku(ji,jj)
217               z2d(ji,jj) = un(ji,jj,ikbot)
218            END DO
219         END DO
220         CALL iom_put( "sbu", z2d )                ! bottom 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         DO jj = 1, jpj
227            DO ji = 1, jpi
228               ikbot = mbkv(ji,jj)
229               z2d(ji,jj) = vn(ji,jj,ikbot)
230            END DO
231         END DO
232         CALL iom_put( "sbv", z2d )                ! bottom j-current
233      ENDIF
234
235      IF( ln_zad_Aimp ) wn = wn + wi               ! Recombine explicit and implicit parts of vertical velocity for diagnostic output
236      !
237      CALL iom_put( "woce", wn )                   ! vertical velocity
238      IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN   ! vertical mass transport & its square value
239         ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
240         z2d(:,:) = rau0 * e1e2t(:,:)
241         DO jk = 1, jpk
242            z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
243         END DO
244         CALL iom_put( "w_masstr" , z3d ) 
245         IF( iom_use('w_masstr2') )   CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
246      ENDIF
247      !
248      IF( ln_zad_Aimp ) wn = wn - wi               ! Remove implicit part of vertical velocity that was added for diagnostic output
249
250      CALL iom_put( "avt" , avt )                  ! T vert. eddy diff. coef.
251      CALL iom_put( "avs" , avs )                  ! S vert. eddy diff. coef.
252      CALL iom_put( "avm" , avm )                  ! T vert. eddy visc. coef.
253
254      IF( iom_use('logavt') )   CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
255      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
256
257      IF ( iom_use("sstgrad_"//ttype) .OR. iom_use("sstgrad2_"//ttype) ) THEN
258         DO jj = 2, jpjm1                                    ! sst gradient
259            DO ji = fs_2, fs_jpim1   ! vector opt.
260               zztmp  = tsn(ji,jj,1,jp_tem)
261               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)
262               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)
263               z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy )   &
264                  &              * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
265            END DO
266         END DO
267         CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
268         CALL iom_put( "sstgrad2_"//ttype,  z2d )          ! square of module of sst gradient
269         z2d(:,:) = SQRT( z2d(:,:) )
270         CALL iom_put( "sstgrad_"//ttype ,  z2d )          ! module of sst gradient
271      ENDIF
272         
273      ! heat and salt contents
274      IF( iom_use("heatc") ) THEN
275         z2d(:,:)  = 0._wp 
276         DO jk = 1, jpkm1
277            DO jj = 1, jpj
278               DO ji = 1, jpi
279                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
280               END DO
281            END DO
282         END DO
283         CALL iom_put( "heatc", rau0_rcp * z2d )   ! vertically integrated heat content (J/m2)
284      ENDIF
285
286      IF( iom_use("saltc") ) THEN
287         z2d(:,:)  = 0._wp 
288         DO jk = 1, jpkm1
289            DO jj = 1, jpj
290               DO ji = 1, jpi
291                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
292               END DO
293            END DO
294         END DO
295         CALL iom_put( "saltc", rau0 * z2d )          ! vertically integrated salt content (PSU*kg/m2)
296      ENDIF
297      !
298      IF ( iom_use("eken") ) THEN
299         z3d(:,:,jpk) = 0._wp 
300         DO jk = 1, jpkm1
301            DO jj = 2, jpjm1
302               DO ji = fs_2, fs_jpim1   ! vector opt.
303                  zztmp  = 0.25_wp * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
304                  z3d(ji,jj,jk) = zztmp * (  un(ji-1,jj,jk)**2 * e2u(ji-1,jj) * e3u_n(ji-1,jj,jk)   &
305                     &                     + un(ji  ,jj,jk)**2 * e2u(ji  ,jj) * e3u_n(ji  ,jj,jk)   &
306                     &                     + vn(ji,jj-1,jk)**2 * e1v(ji,jj-1) * e3v_n(ji,jj-1,jk)   &
307                     &                     + vn(ji,jj  ,jk)**2 * e1v(ji,jj  ) * e3v_n(ji,jj  ,jk)   )
308               END DO
309            END DO
310         END DO
311         CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
312         CALL iom_put( "eken", z3d )                 ! kinetic energy
313      ENDIF
314      !
315      CALL iom_put( "hdiv", hdivn )                  ! Horizontal divergence
316      !
317      IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
318         z3d(:,:,jpk) = 0.e0
319         z2d(:,:) = 0.e0
320         DO jk = 1, jpkm1
321            z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)
322            z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
323         END DO
324         CALL iom_put( "u_masstr"     , z3d )         ! mass transport in i-direction
325         CALL iom_put( "u_masstr_vint", z2d )         ! mass transport in i-direction vertical sum
326      ENDIF
327     
328      IF( iom_use("u_heattr") ) THEN
329         z2d(:,:) = 0._wp 
330         DO jk = 1, jpkm1
331            DO jj = 2, jpjm1
332               DO ji = fs_2, fs_jpim1   ! vector opt.
333                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
334               END DO
335            END DO
336         END DO
337         CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
338         CALL iom_put( "u_heattr", 0.5*rcp * z2d )    ! heat transport in i-direction
339      ENDIF
340
341      IF( iom_use("u_salttr") ) THEN
342         z2d(:,:) = 0.e0 
343         DO jk = 1, jpkm1
344            DO jj = 2, jpjm1
345               DO ji = fs_2, fs_jpim1   ! vector opt.
346                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
347               END DO
348            END DO
349         END DO
350         CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
351         CALL iom_put( "u_salttr", 0.5 * z2d )        ! heat transport in i-direction
352      ENDIF
353
354     
355      IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
356         z3d(:,:,jpk) = 0.e0
357         DO jk = 1, jpkm1
358            z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)
359         END DO
360         CALL iom_put( "v_masstr", z3d )              ! mass transport in j-direction
361      ENDIF
362     
363      IF( iom_use("v_heattr") ) THEN
364         z2d(:,:) = 0.e0 
365         DO jk = 1, jpkm1
366            DO jj = 2, jpjm1
367               DO ji = fs_2, fs_jpim1   ! vector opt.
368                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
369               END DO
370            END DO
371         END DO
372         CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
373         CALL iom_put( "v_heattr", 0.5*rcp * z2d )    !  heat transport in j-direction
374      ENDIF
375
376      IF( iom_use("v_salttr") ) THEN
377         z2d(:,:) = 0._wp 
378         DO jk = 1, jpkm1
379            DO jj = 2, jpjm1
380               DO ji = fs_2, fs_jpim1   ! vector opt.
381                  z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
382               END DO
383            END DO
384         END DO
385         CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
386         CALL iom_put( "v_salttr", 0.5 * z2d )        !  heat transport in j-direction
387      ENDIF
388
389      IF( iom_use("tosmint_"//ttype) ) THEN
390         z2d(:,:) = 0._wp
391         DO jk = 1, jpkm1
392            DO jj = 2, jpjm1
393               DO ji = fs_2, fs_jpim1   ! vector opt.
394                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) *  tsn(ji,jj,jk,jp_tem)
395               END DO
396            END DO
397         END DO
398         CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
399         CALL iom_put( "tosmint_"//ttype, rau0 * z2d )        ! Vertical integral of temperature
400      ENDIF
401      IF( iom_use("somint_"//stype) ) THEN
402         z2d(:,:)=0._wp
403         DO jk = 1, jpkm1
404            DO jj = 2, jpjm1
405               DO ji = fs_2, fs_jpim1   ! vector opt.
406                  z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
407               END DO
408            END DO
409         END DO
410         CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
411         CALL iom_put( "somint_"//stype, rau0 * z2d )         ! Vertical integral of salinity
412      ENDIF
413
414      CALL iom_put( "bn2", rn2 )                      ! Brunt-Vaisala buoyancy frequency (N^2)
415      !
416
417      IF (ln_diatmb)   CALL dia_tmb                   ! tmb values
418         
419      IF (ln_dia25h)   CALL dia_25h( kt )             ! 25h averaging
420
421      IF( ln_timing )   CALL timing_stop('dia_wri')
422      !
423   END SUBROUTINE dia_wri
424
425#else
426   !!----------------------------------------------------------------------
427   !!   Default option                                  use IOIPSL  library
428   !!----------------------------------------------------------------------
429
430   INTEGER FUNCTION dia_wri_alloc()
431      !!----------------------------------------------------------------------
432      INTEGER, DIMENSION(2) :: ierr
433      !!----------------------------------------------------------------------
434      ierr = 0
435      ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) ,     &
436         &      ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) ,     &
437         &      ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
438         !
439      dia_wri_alloc = MAXVAL(ierr)
440      CALL mpp_sum( 'diawri', dia_wri_alloc )
441      !
442   END FUNCTION dia_wri_alloc
443
444   
445   SUBROUTINE dia_wri( kt )
446      !!---------------------------------------------------------------------
447      !!                  ***  ROUTINE dia_wri  ***
448      !!                   
449      !! ** Purpose :   Standard output of opa: dynamics and tracer fields
450      !!      NETCDF format is used by default
451      !!
452      !! ** Method  :   At the beginning of the first time step (nit000),
453      !!      define all the NETCDF files and fields
454      !!      At each time step call histdef to compute the mean if ncessary
455      !!      Each nn_write time step, output the instantaneous or mean fields
456      !!----------------------------------------------------------------------
457      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
458      !
459      LOGICAL ::   ll_print = .FALSE.                        ! =T print and flush numout
460      CHARACTER (len=40) ::   clhstnam, clop, clmx           ! local names
461      INTEGER  ::   inum = 11                                ! temporary logical unit
462      INTEGER  ::   ji, jj, jk                               ! dummy loop indices
463      INTEGER  ::   ierr                                     ! error code return from allocation
464      INTEGER  ::   iimi, iima, ipk, it, itmod, ijmi, ijma   ! local integers
465      INTEGER  ::   jn, ierror                               ! local integers
466      REAL(wp) ::   zsto, zout, zmax, zjulian                ! local scalars
467      !
468      REAL(wp), DIMENSION(jpi,jpj)   :: zw2d       ! 2D workspace
469      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d       ! 3D workspace
470      !!----------------------------------------------------------------------
471      !
472      IF( ninist == 1 ) THEN     !==  Output the initial state and forcings  ==!
473         CALL dia_wri_state( 'output.init' )
474         ninist = 0
475      ENDIF
476      !
477      IF( nn_write == -1 )   RETURN   ! we will never do any output
478      !
479      IF( ln_timing )   CALL timing_start('dia_wri')
480      !
481      ! 0. Initialisation
482      ! -----------------
483
484      ll_print = .FALSE.                  ! local variable for debugging
485      ll_print = ll_print .AND. lwp
486
487      ! Define frequency of output and means
488      clop = "x"         ! no use of the mask value (require less cpu time and otherwise the model crashes)
489#if defined key_diainstant
490      zsto = nn_write * rdt
491      clop = "inst("//TRIM(clop)//")"
492#else
493      zsto=rdt
494      clop = "ave("//TRIM(clop)//")"
495#endif
496      zout = nn_write * rdt
497      zmax = ( nitend - nit000 + 1 ) * rdt
498
499      ! Define indices of the horizontal output zoom and vertical limit storage
500      iimi = 1      ;      iima = jpi
501      ijmi = 1      ;      ijma = jpj
502      ipk = jpk
503
504      ! define time axis
505      it = kt
506      itmod = kt - nit000 + 1
507
508
509      ! 1. Define NETCDF files and fields at beginning of first time step
510      ! -----------------------------------------------------------------
511
512      IF( kt == nit000 ) THEN
513
514         ! Define the NETCDF files (one per grid)
515
516         ! Compute julian date from starting date of the run
517         CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
518         zjulian = zjulian - adatrj   !   set calendar origin to the beginning of the experiment
519         IF(lwp)WRITE(numout,*)
520         IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear,   &
521            &                    ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
522         IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma,   &
523                                 ' limit storage in depth = ', ipk
524
525         ! WRITE root name in date.file for use by postpro
526         IF(lwp) THEN
527            CALL dia_nam( clhstnam, nn_write,' ' )
528            CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
529            WRITE(inum,*) clhstnam
530            CLOSE(inum)
531         ENDIF
532
533         ! Define the T grid FILE ( nid_T )
534
535         CALL dia_nam( clhstnam, nn_write, 'grid_T' )
536         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
537         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
538            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
539            &          nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
540         CALL histvert( nid_T, "deptht", "Vertical T levels",      &  ! Vertical grid: gdept
541            &           "m", ipk, gdept_1d, nz_T, "down" )
542         !                                                            ! Index of ocean points
543         CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T  )      ! volume
544         CALL wheneq( jpi*jpj    , tmask, 1, 1., ndex_hT, ndim_hT )      ! surface
545         !
546         IF( ln_icebergs ) THEN
547            !
548            !! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
549            !! that routine is called from nemogcm, so do it here immediately before its needed
550            ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
551            CALL mpp_sum( 'diawri', ierror )
552            IF( ierror /= 0 ) THEN
553               CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
554               RETURN
555            ENDIF
556            !
557            !! iceberg vertical coordinate is class number
558            CALL histvert( nid_T, "class", "Iceberg class",      &  ! Vertical grid: class
559               &           "number", nclasses, class_num, nb_T )
560            !
561            !! each class just needs the surface index pattern
562            ndim_bT = 3
563            DO jn = 1,nclasses
564               ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
565            ENDDO
566            !
567         ENDIF
568
569         ! Define the U grid FILE ( nid_U )
570
571         CALL dia_nam( clhstnam, nn_write, 'grid_U' )
572         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam    ! filename
573         CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu,           &  ! Horizontal grid: glamu and gphiu
574            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
575            &          nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
576         CALL histvert( nid_U, "depthu", "Vertical U levels",      &  ! Vertical grid: gdept
577            &           "m", ipk, gdept_1d, nz_U, "down" )
578         !                                                            ! Index of ocean points
579         CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U  )      ! volume
580         CALL wheneq( jpi*jpj    , umask, 1, 1., ndex_hU, ndim_hU )      ! surface
581
582         ! Define the V grid FILE ( nid_V )
583
584         CALL dia_nam( clhstnam, nn_write, 'grid_V' )                   ! filename
585         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
586         CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv,           &  ! Horizontal grid: glamv and gphiv
587            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
588            &          nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
589         CALL histvert( nid_V, "depthv", "Vertical V levels",      &  ! Vertical grid : gdept
590            &          "m", ipk, gdept_1d, nz_V, "down" )
591         !                                                            ! Index of ocean points
592         CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V  )      ! volume
593         CALL wheneq( jpi*jpj    , vmask, 1, 1., ndex_hV, ndim_hV )      ! surface
594
595         ! Define the W grid FILE ( nid_W )
596
597         CALL dia_nam( clhstnam, nn_write, 'grid_W' )                   ! filename
598         IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
599         CALL histbeg( clhstnam, jpi, glamt, jpj, gphit,           &  ! Horizontal grid: glamt and gphit
600            &          iimi, iima-iimi+1, ijmi, ijma-ijmi+1,       &
601            &          nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
602         CALL histvert( nid_W, "depthw", "Vertical W levels",      &  ! Vertical grid: gdepw
603            &          "m", ipk, gdepw_1d, nz_W, "down" )
604
605
606         ! Declare all the output fields as NETCDF variables
607
608         !                                                                                      !!! nid_T : 3D
609         CALL histdef( nid_T, "votemper", "Temperature"                        , "C"      ,   &  ! tn
610            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
611         CALL histdef( nid_T, "vosaline", "Salinity"                           , "PSU"    ,   &  ! sn
612            &          jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
613         IF(  .NOT.ln_linssh  ) THEN
614            CALL histdef( nid_T, "vovvle3t", "Level thickness"                    , "m"      ,&  ! e3t_n
615            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
616            CALL histdef( nid_T, "vovvldep", "T point depth"                      , "m"      ,&  ! e3t_n
617            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
618            CALL histdef( nid_T, "vovvldef", "Squared level deformation"          , "%^2"    ,&  ! e3t_n
619            &             jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
620         ENDIF
621         !                                                                                      !!! nid_T : 2D
622         CALL histdef( nid_T, "sosstsst", "Sea Surface temperature"            , "C"      ,   &  ! sst
623            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
624         CALL histdef( nid_T, "sosaline", "Sea Surface Salinity"               , "PSU"    ,   &  ! sss
625            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
626         CALL histdef( nid_T, "sossheig", "Sea Surface Height"                 , "m"      ,   &  ! ssh
627            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
628         CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! (emp-rnf)
629            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
630         CALL histdef( nid_T, "sorunoff", "River runoffs"                      , "Kg/m2/s",   &  ! runoffs
631            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
632         CALL histdef( nid_T, "sosfldow", "downward salt flux"                 , "PSU/m2/s",  &  ! sfx
633            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
634         IF(  ln_linssh  ) THEN
635            CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature"     &  ! emp * tsn(:,:,1,jp_tem)
636            &                                                                  , "KgC/m2/s",  &  ! sosst_cd
637            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
638            CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity"        &  ! emp * tsn(:,:,1,jp_sal)
639            &                                                                  , "KgPSU/m2/s",&  ! sosss_cd
640            &             jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
641         ENDIF
642         CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux"             , "W/m2"   ,   &  ! qns + qsr
643            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
644         CALL histdef( nid_T, "soshfldo", "Shortwave Radiation"                , "W/m2"   ,   &  ! qsr
645            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
646         CALL histdef( nid_T, "somixhgt", "Turbocline Depth"                   , "m"      ,   &  ! hmld
647            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
648         CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01"             , "m"      ,   &  ! hmlp
649            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
650         CALL histdef( nid_T, "soicecov", "Ice fraction"                       , "[0,1]"  ,   &  ! fr_i
651            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
652         CALL histdef( nid_T, "sowindsp", "wind speed at 10m"                  , "m/s"    ,   &  ! wndm
653            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
654!
655         IF( ln_icebergs ) THEN
656            CALL histdef( nid_T, "calving"             , "calving mass input"                       , "kg/s"   , &
657               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
658            CALL histdef( nid_T, "calving_heat"        , "calving heat flux"                        , "XXXX"   , &
659               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
660            CALL histdef( nid_T, "berg_floating_melt"  , "Melt rate of icebergs + bits"             , "kg/m2/s", &
661               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
662            CALL histdef( nid_T, "berg_stored_ice"     , "Accumulated ice mass by class"            , "kg"     , &
663               &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
664            IF( ln_bergdia ) THEN
665               CALL histdef( nid_T, "berg_melt"           , "Melt rate of icebergs"                    , "kg/m2/s", &
666                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
667               CALL histdef( nid_T, "berg_buoy_melt"      , "Buoyancy component of iceberg melt rate"  , "kg/m2/s", &
668                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
669               CALL histdef( nid_T, "berg_eros_melt"      , "Erosion component of iceberg melt rate"   , "kg/m2/s", &
670                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
671               CALL histdef( nid_T, "berg_conv_melt"      , "Convective component of iceberg melt rate", "kg/m2/s", &
672                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
673               CALL histdef( nid_T, "berg_virtual_area"   , "Virtual coverage by icebergs"             , "m2"     , &
674                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
675               CALL histdef( nid_T, "bits_src"           , "Mass source of bergy bits"                , "kg/m2/s", &
676                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
677               CALL histdef( nid_T, "bits_melt"          , "Melt rate of bergy bits"                  , "kg/m2/s", &
678                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
679               CALL histdef( nid_T, "bits_mass"          , "Bergy bit density field"                  , "kg/m2"  , &
680                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
681               CALL histdef( nid_T, "berg_mass"           , "Iceberg density field"                    , "kg/m2"  , &
682                  &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
683               CALL histdef( nid_T, "berg_real_calving"   , "Calving into iceberg class"               , "kg/s"   , &
684                  &          jpi, jpj, nh_T, nclasses  , 1, nclasses  , nb_T , 32, clop, zsto, zout )
685            ENDIF
686         ENDIF
687
688         IF( ln_ssr ) THEN
689            CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping"         , "W/m2"   ,   &  ! qrp
690               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
691            CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping"        , "Kg/m2/s",   &  ! erp
692               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
693            CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping"         , "Kg/m2/s",   &  ! erp * sn
694               &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
695         ENDIF
696       
697         clmx ="l_max(only(x))"    ! max index on a period
698!         CALL histdef( nid_T, "sobowlin", "Bowl Index"                         , "W-point",   &  ! bowl INDEX
699!            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clmx, zsto, zout )
700#if defined key_diahth
701         CALL histdef( nid_T, "sothedep", "Thermocline Depth"                  , "m"      ,   & ! hth
702            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
703         CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm"              , "m"      ,   & ! hd20
704            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
705         CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm"              , "m"      ,   & ! hd28
706            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
707         CALL histdef( nid_T, "sohtc300", "Heat content 300 m"                 , "J/m2"   ,   & ! htc3
708            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
709#endif
710
711         CALL histend( nid_T, snc4chunks=snc4set )
712
713         !                                                                                      !!! nid_U : 3D
714         CALL histdef( nid_U, "vozocrtx", "Zonal Current"                      , "m/s"    ,   &  ! un
715            &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
716         IF( ln_wave .AND. ln_sdw) THEN
717            CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current"         , "m/s"    ,   &  ! usd
718               &          jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
719         ENDIF
720         !                                                                                      !!! nid_U : 2D
721         CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis"           , "N/m2"   ,   &  ! utau
722            &          jpi, jpj, nh_U, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
723
724         CALL histend( nid_U, snc4chunks=snc4set )
725
726         !                                                                                      !!! nid_V : 3D
727         CALL histdef( nid_V, "vomecrty", "Meridional Current"                 , "m/s"    ,   &  ! vn
728            &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
729         IF( ln_wave .AND. ln_sdw) THEN
730            CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current"    , "m/s"    ,   &  ! vsd
731               &          jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
732         ENDIF
733         !                                                                                      !!! nid_V : 2D
734         CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis"           , "N/m2"   ,   &  ! vtau
735            &          jpi, jpj, nh_V, 1  , 1, 1  , - 99, 32, clop, zsto, zout )
736
737         CALL histend( nid_V, snc4chunks=snc4set )
738
739         !                                                                                      !!! nid_W : 3D
740         CALL histdef( nid_W, "vovecrtz", "Vertical Velocity"                  , "m/s"    ,   &  ! wn
741            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
742         CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity"          , "m2/s"   ,   &  ! avt
743            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
744         CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity"             , "m2/s"  ,   &  ! avm
745            &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
746
747         IF( ln_zdfddm ) THEN
748            CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity"    , "m2/s"   ,   &  ! avs
749               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
750         ENDIF
751         
752         IF( ln_wave .AND. ln_sdw) THEN
753            CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current"   , "m/s"    ,   &  ! wsd
754               &          jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
755         ENDIF
756         !                                                                                      !!! nid_W : 2D
757         CALL histend( nid_W, snc4chunks=snc4set )
758
759         IF(lwp) WRITE(numout,*)
760         IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
761         IF(ll_print) CALL FLUSH(numout )
762
763      ENDIF
764
765      ! 2. Start writing data
766      ! ---------------------
767
768      ! ndex(1) est utilise ssi l'avant dernier argument est different de
769      ! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
770      ! donne le nombre d'elements, et ndex la liste des indices a sortir
771
772      IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN
773         WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
774         WRITE(numout,*) '~~~~~~ '
775      ENDIF
776
777      IF( .NOT.ln_linssh ) THEN
778         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! heat content
779         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T  )   ! salt content
780         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface heat content
781         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT )   ! sea surface salinity content
782      ELSE
783         CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T  )   ! temperature
784         CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T  )   ! salinity
785         CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT )   ! sea surface temperature
786         CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT )   ! sea surface salinity
787      ENDIF
788      IF( .NOT.ln_linssh ) THEN
789         zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
790         CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T  )   ! level thickness
791         CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T  )   ! t-point depth
792         CALL histwrite( nid_T, "vovvldef", it, zw3d             , ndim_T , ndex_T  )   ! level thickness deformation
793      ENDIF
794      CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
795      CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf )   , ndim_hT, ndex_hT )   ! upward water flux
796      CALL histwrite( nid_T, "sorunoff", it, rnf           , ndim_hT, ndex_hT )   ! river runoffs
797      CALL histwrite( nid_T, "sosfldow", it, sfx           , ndim_hT, ndex_hT )   ! downward salt flux
798                                                                                  ! (includes virtual salt flux beneath ice
799                                                                                  ! in linear free surface case)
800      IF( ln_linssh ) THEN
801         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
802         CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sst
803         zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
804         CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT )          ! c/d term on sss
805      ENDIF
806      CALL histwrite( nid_T, "sohefldo", it, qns + qsr     , ndim_hT, ndex_hT )   ! total heat flux
807      CALL histwrite( nid_T, "soshfldo", it, qsr           , ndim_hT, ndex_hT )   ! solar heat flux
808      CALL histwrite( nid_T, "somixhgt", it, hmld          , ndim_hT, ndex_hT )   ! turbocline depth
809      CALL histwrite( nid_T, "somxl010", it, hmlp          , ndim_hT, ndex_hT )   ! mixed layer depth
810      CALL histwrite( nid_T, "soicecov", it, fr_i          , ndim_hT, ndex_hT )   ! ice fraction   
811      CALL histwrite( nid_T, "sowindsp", it, wndm          , ndim_hT, ndex_hT )   ! wind speed   
812!
813      IF( ln_icebergs ) THEN
814         !
815         CALL histwrite( nid_T, "calving"             , it, berg_grid%calving      , ndim_hT, ndex_hT ) 
816         CALL histwrite( nid_T, "calving_heat"        , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )         
817         CALL histwrite( nid_T, "berg_floating_melt"  , it, berg_grid%floating_melt, ndim_hT, ndex_hT ) 
818         !
819         CALL histwrite( nid_T, "berg_stored_ice"     , it, berg_grid%stored_ice   , ndim_bT, ndex_bT )
820         !
821         IF( ln_bergdia ) THEN
822            CALL histwrite( nid_T, "berg_melt"           , it, berg_melt        , ndim_hT, ndex_hT   ) 
823            CALL histwrite( nid_T, "berg_buoy_melt"      , it, buoy_melt        , ndim_hT, ndex_hT   ) 
824            CALL histwrite( nid_T, "berg_eros_melt"      , it, eros_melt        , ndim_hT, ndex_hT   ) 
825            CALL histwrite( nid_T, "berg_conv_melt"      , it, conv_melt        , ndim_hT, ndex_hT   ) 
826            CALL histwrite( nid_T, "berg_virtual_area"   , it, virtual_area     , ndim_hT, ndex_hT   ) 
827            CALL histwrite( nid_T, "bits_src"            , it, bits_src         , ndim_hT, ndex_hT   ) 
828            CALL histwrite( nid_T, "bits_melt"           , it, bits_melt        , ndim_hT, ndex_hT   ) 
829            CALL histwrite( nid_T, "bits_mass"           , it, bits_mass        , ndim_hT, ndex_hT   ) 
830            CALL histwrite( nid_T, "berg_mass"           , it, berg_mass        , ndim_hT, ndex_hT   ) 
831            !
832            CALL histwrite( nid_T, "berg_real_calving"   , it, real_calving     , ndim_bT, ndex_bT   )
833         ENDIF
834      ENDIF
835
836      IF( ln_ssr ) THEN
837         CALL histwrite( nid_T, "sohefldp", it, qrp           , ndim_hT, ndex_hT )   ! heat flux damping
838         CALL histwrite( nid_T, "sowafldp", it, erp           , ndim_hT, ndex_hT )   ! freshwater flux damping
839         zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
840         CALL histwrite( nid_T, "sosafldp", it, zw2d          , ndim_hT, ndex_hT )   ! salt flux damping
841      ENDIF
842!      zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
843!      CALL histwrite( nid_T, "sobowlin", it, zw2d          , ndim_hT, ndex_hT )   ! ???
844
845#if defined key_diahth
846      CALL histwrite( nid_T, "sothedep", it, hth           , ndim_hT, ndex_hT )   ! depth of the thermocline
847      CALL histwrite( nid_T, "so20chgt", it, hd20          , ndim_hT, ndex_hT )   ! depth of the 20 isotherm
848      CALL histwrite( nid_T, "so28chgt", it, hd28          , ndim_hT, ndex_hT )   ! depth of the 28 isotherm
849      CALL histwrite( nid_T, "sohtc300", it, htc3          , ndim_hT, ndex_hT )   ! first 300m heaat content
850#endif
851
852      CALL histwrite( nid_U, "vozocrtx", it, un            , ndim_U , ndex_U )    ! i-current
853      CALL histwrite( nid_U, "sozotaux", it, utau          , ndim_hU, ndex_hU )   ! i-wind stress
854
855      CALL histwrite( nid_V, "vomecrty", it, vn            , ndim_V , ndex_V  )   ! j-current
856      CALL histwrite( nid_V, "sometauy", it, vtau          , ndim_hV, ndex_hV )   ! j-wind stress
857
858      IF( ln_zad_Aimp ) THEN
859         CALL histwrite( nid_W, "vovecrtz", it, wn + wi     , ndim_T, ndex_T )    ! vert. current
860      ELSE
861         CALL histwrite( nid_W, "vovecrtz", it, wn          , ndim_T, ndex_T )    ! vert. current
862      ENDIF
863      CALL histwrite( nid_W, "votkeavt", it, avt            , ndim_T, ndex_T )    ! T vert. eddy diff. coef.
864      CALL histwrite( nid_W, "votkeavm", it, avm            , ndim_T, ndex_T )    ! T vert. eddy visc. coef.
865      IF( ln_zdfddm ) THEN
866         CALL histwrite( nid_W, "voddmavs", it, avs         , ndim_T, ndex_T )    ! S vert. eddy diff. coef.
867      ENDIF
868
869      IF( ln_wave .AND. ln_sdw ) THEN
870         CALL histwrite( nid_U, "sdzocrtx", it, usd         , ndim_U , ndex_U )    ! i-StokesDrift-current
871         CALL histwrite( nid_V, "sdmecrty", it, vsd         , ndim_V , ndex_V )    ! j-StokesDrift-current
872         CALL histwrite( nid_W, "sdvecrtz", it, wsd         , ndim_T , ndex_T )    ! StokesDrift vert. current
873      ENDIF
874
875      ! 3. Close all files
876      ! ---------------------------------------
877      IF( kt == nitend ) THEN
878         CALL histclo( nid_T )
879         CALL histclo( nid_U )
880         CALL histclo( nid_V )
881         CALL histclo( nid_W )
882      ENDIF
883      !
884      IF( ln_timing )   CALL timing_stop('dia_wri')
885      !
886   END SUBROUTINE dia_wri
887#endif
888
889   SUBROUTINE dia_wri_state( cdfile_name )
890      !!---------------------------------------------------------------------
891      !!                 ***  ROUTINE dia_wri_state  ***
892      !!       
893      !! ** Purpose :   create a NetCDF file named cdfile_name which contains
894      !!      the instantaneous ocean state and forcing fields.
895      !!        Used to find errors in the initial state or save the last
896      !!      ocean state in case of abnormal end of a simulation
897      !!
898      !! ** Method  :   NetCDF files using ioipsl
899      !!      File 'output.init.nc'  is created if ninist = 1 (namelist)
900      !!      File 'output.abort.nc' is created in case of abnormal job end
901      !!----------------------------------------------------------------------
902      CHARACTER (len=* ), INTENT( in ) ::   cdfile_name      ! name of the file created
903      !!
904      INTEGER :: inum
905      !!----------------------------------------------------------------------
906      !
907      IF(lwp) WRITE(numout,*)
908      IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
909      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   and forcing fields file created '
910      IF(lwp) WRITE(numout,*) '                and named :', cdfile_name, '...nc'
911
912#if defined key_si3
913     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
914#else
915     CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
916#endif
917
918      CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) )    ! now temperature
919      CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) )    ! now salinity
920      CALL iom_rstput( 0, 0, inum, 'sossheig', sshn              )    ! sea surface height
921      CALL iom_rstput( 0, 0, inum, 'vozocrtx', un                )    ! now i-velocity
922      CALL iom_rstput( 0, 0, inum, 'vomecrty', vn                )    ! now j-velocity
923      IF( ln_zad_Aimp ) THEN
924         CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn + wi        )    ! now k-velocity
925      ELSE
926         CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn             )    ! now k-velocity
927      ENDIF
928      IF( ALLOCATED(ahtu) ) THEN
929         CALL iom_rstput( 0, 0, inum,  'ahtu', ahtu              )    ! aht at u-point
930         CALL iom_rstput( 0, 0, inum,  'ahtv', ahtv              )    ! aht at v-point
931      ENDIF
932      IF( ALLOCATED(ahmt) ) THEN
933         CALL iom_rstput( 0, 0, inum,  'ahmt', ahmt              )    ! ahmt at u-point
934         CALL iom_rstput( 0, 0, inum,  'ahmf', ahmf              )    ! ahmf at v-point
935      ENDIF
936      CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf         )    ! freshwater budget
937      CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns         )    ! total heat flux
938      CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr               )    ! solar heat flux
939      CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i              )    ! ice fraction
940      CALL iom_rstput( 0, 0, inum, 'sozotaux', utau              )    ! i-wind stress
941      CALL iom_rstput( 0, 0, inum, 'sometauy', vtau              )    ! j-wind stress
942      IF(  .NOT.ln_linssh  ) THEN             
943         CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept_n        )    !  T-cell depth
944         CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n          )    !  T-cell thickness 
945      END IF
946      IF( ln_wave .AND. ln_sdw ) THEN
947         CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd            )    ! now StokesDrift i-velocity
948         CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd            )    ! now StokesDrift j-velocity
949         CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity
950      ENDIF
951#if defined key_si3
952      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid
953         CALL ice_wri_state( inum )
954      ENDIF
955#endif
956      !
957      CALL iom_close( inum )
958      !
959   END SUBROUTINE dia_wri_state
960
961   !!======================================================================
962END MODULE diawri
Note: See TracBrowser for help on using the repository browser.