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/2020/dev_r12527_Gurvan_ShallowWater/src/SWE – NEMO

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/diawri.F90 @ 12667

Last change on this file since 12667 was 12667, checked in by gm, 4 years ago

Shallow Water Eq. update

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