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/2019/dev_r11943_MERGE_2019/src/OCE/DIA – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DIA/diawri.F90 @ 12150

Last change on this file since 12150 was 12150, checked in by davestorkey, 4 years ago

2019/dev_r11943_MERGE_2019: Merge in UKMO_MERGE_2019.

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