New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diawri.F90 in NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DIA – NEMO

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

Last change on this file since 12252 was 12252, checked in by smasson, 4 years ago

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