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

source: NEMO/branches/UKMO/NEMO_4.0.2_GO8_package/src/OCE/DIA/diawri.F90 @ 12660

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

UKMO/NEMO_4.0.2_GO8_package: copy over changes from NEMO_4.0.1_GO8_package branch.

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