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

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

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