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

source: NEMO/branches/UKMO/NEMO_4.0.4_GO6_mixing/src/OCE/DIA/diawri.F90 @ 15668

Last change on this file since 15668 was 15668, checked in by davestorkey, 2 years ago

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