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

Last change on this file since 14442 was 14442, checked in by davestorkey, 3 years ago

UKMO/NEMO_4.0.4_GO8_package : bug fix for previous commit.

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