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.4_GO8_package/src/OCE/DIA – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/OCE/DIA/diawri.F90 @ 14078

Last change on this file since 14078 was 14078, checked in by cguiavarch, 3 years ago

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