Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
MODULE diawri
!!======================================================================
!! *** MODULE diawri ***
!! Ocean diagnostics : write ocean output files
!!=====================================================================
!! History : OPA ! 1991-03 (M.-A. Foujols) Original code
!! 4.0 ! 1991-11 (G. Madec)
!! ! 1992-06 (M. Imbard) correction restart file
!! ! 1992-07 (M. Imbard) split into diawri and rstwri
!! ! 1993-03 (M. Imbard) suppress writibm
!! ! 1998-01 (C. Levy) NETCDF format using ioipsl INTERFACE
!! ! 1999-02 (E. Guilyardi) name of netCDF files + variables
!! 8.2 ! 2000-06 (M. Imbard) Original code (diabort.F)
!! NEMO 1.0 ! 2002-06 (A.Bozec, E. Durand) Original code (diainit.F)
!! - ! 2002-09 (G. Madec) F90: Free form and module
!! - ! 2002-12 (G. Madec) merge of diabort and diainit, F90
!! ! 2005-11 (V. Garnier) Surface pressure gradient organization
!! 3.2 ! 2008-11 (B. Lemaire) creation from old diawri
!! 3.7 ! 2014-01 (G. Madec) remove eddy induced velocity from no-IOM output
!! ! change name of output variables in dia_wri_state
!! 4.0 ! 2020-10 (A. Nasser, S. Techene) add diagnostic for SWE
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! dia_wri : create the standart output files
!! dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers
USE isf_oce
USE isfcpl
USE abl ! abl variables in case ln_abl = .true.
USE dom_oce ! ocean space and time domain
USE phycst ! physical constants
USE dianam ! build name of file (routine)
USE diahth ! thermocline diagnostics
USE dynadv , ONLY: ln_dynadv_vec
USE icb_oce ! Icebergs
USE icbdia ! Iceberg budgets
USE ldftra ! lateral physics: eddy diffusivity coef.
USE ldfdyn ! lateral physics: eddy viscosity coef.
USE sbc_oce ! Surface boundary condition: ocean fields
USE sbc_ice ! Surface boundary condition: ice fields
USE sbcssr ! restoring term toward SST/SSS climatology
USE sbcwave ! wave parameters
USE wet_dry ! wetting and drying
USE zdf_oce ! ocean vertical physics
USE zdfdrg ! ocean vertical physics: top/bottom friction
USE zdfmxl ! mixed layer
USE zdfosm ! mixed layer
!
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE in_out_manager ! I/O manager
USE dia25h ! 25h Mean output
USE iom !
USE ioipsl !
#if defined key_si3
USE ice
USE icewri
#endif
USE lib_mpp ! MPP library
USE timing ! preformance summary
USE diu_bulk ! diurnal warm layer
USE diu_coolskin ! Cool skin
IMPLICIT NONE
PRIVATE
PUBLIC dia_wri ! routines called by step.F90
PUBLIC dia_wri_state
PUBLIC dia_wri_alloc ! Called by nemogcm module
#if ! defined key_xios
PUBLIC dia_wri_alloc_abl ! Called by sbcabl module (if ln_abl = .true.)
#endif
INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file
INTEGER :: nb_T , ndim_bT ! grid_T file
INTEGER :: nid_U, nz_U, nh_U, ndim_U, ndim_hU ! grid_U file
INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file
INTEGER :: nid_W, nz_W, nh_W ! grid_W file
INTEGER :: nid_A, nz_A, nh_A, ndim_A, ndim_hA ! grid_ABL file
INTEGER :: ndex(1) ! ???
INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hA, ndex_A ! ABL
INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: diawri.F90 15141 2021-07-23 14:20:12Z smasson $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
#if defined key_xios
!!----------------------------------------------------------------------
!! 'key_xios' use IOM library
!!----------------------------------------------------------------------
INTEGER FUNCTION dia_wri_alloc()
!
dia_wri_alloc = 0
!
END FUNCTION dia_wri_alloc
SUBROUTINE dia_wri( kt, Kmm )
!!---------------------------------------------------------------------
!! *** ROUTINE dia_wri ***
!!
!! ** Purpose : Standard output of opa: dynamics and tracer fields
!! NETCDF format is used by default
!!
!! ** Method : use iom_put
!!----------------------------------------------------------------------
INTEGER, INTENT( in ) :: kt ! ocean time-step index
INTEGER, INTENT( in ) :: Kmm ! ocean time level index
!!
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ikbot ! local integer
REAL(wp):: zztmp , zztmpx ! local scalar
REAL(wp):: zztmp2, zztmpy ! - -
REAL(wp):: ze3
Clement Rousset
committed
REAL(wp), DIMENSION(A2D(0)) :: z2d0 ! 2D workspace
REAL(wp), DIMENSION(A2D(0) ,jpk) :: z3d0 ! 3D workspace
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z3d ! 3D workspace
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dia_wri')
!
! Output the initial state and forcings
IF( ninist == 1 ) THEN
CALL dia_wri_state( Kmm, 'output.init' )
ninist = 0
ENDIF
! initialize arrays
Clement Rousset
committed
z2d0(:,:) = 0._wp
z3d0(:,:,:) = 0._wp
z3d (:,:,:) = 0._wp
! Output of initial vertical scale factor
CALL iom_put("e3t_0", e3t_0(:,:,:) )
CALL iom_put("e3u_0", e3u_0(:,:,:) )
CALL iom_put("e3v_0", e3v_0(:,:,:) )
CALL iom_put("e3f_0", e3f_0(:,:,:) )
!
IF ( iom_use("tpt_dep") ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
Clement Rousset
committed
z3d0(ji,jj,jk) = gdept(ji,jj,jk,Kmm)
Clement Rousset
committed
CALL iom_put( "tpt_dep", z3d0 )
ENDIF
! --- vertical scale factors --- !
IF ( iom_use("e3t") .OR. iom_use("e3tdef") ) THEN ! time-varying e3t
Clement Rousset
committed
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
z3d(ji,jj,jk) = e3t(ji,jj,jk,Kmm)
END_3D
CALL iom_put( "e3t", z3d )
IF ( iom_use("e3tdef") ) THEN
Clement Rousset
committed
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
z3d(ji,jj,jk) = ( ( z3d(ji,jj,jk) - e3t_0(ji,jj,jk) ) / e3t_0(ji,jj,jk) * 100._wp * tmask(ji,jj,jk) ) ** 2
END_3D
CALL iom_put( "e3tdef", z3d )
ENDIF
ENDIF
IF ( iom_use("e3u") ) THEN ! time-varying e3u
Clement Rousset
committed
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
z3d(ji,jj,jk) = e3u(ji,jj,jk,Kmm)
END_3D
CALL iom_put( "e3u" , z3d )
ENDIF
IF ( iom_use("e3v") ) THEN ! time-varying e3v
Clement Rousset
committed
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
z3d(ji,jj,jk) = e3v(ji,jj,jk,Kmm)
END_3D
CALL iom_put( "e3v" , z3d )
ENDIF
IF ( iom_use("e3w") ) THEN ! time-varying e3w
Clement Rousset
committed
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
z3d(ji,jj,jk) = e3w(ji,jj,jk,Kmm)
END_3D
CALL iom_put( "e3w" , z3d )
ENDIF
IF ( iom_use("e3f") ) THEN ! time-varying e3f caution here at Kaa
Clement Rousset
committed
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
z3d(ji,jj,jk) = e3f(ji,jj,jk)
END_3D
CALL iom_put( "e3f" , z3d )
ENDIF
IF ( iom_use("ssh") ) THEN
IF( ll_wd ) THEN ! sea surface height (brought back to the reference used for wetting and drying)
CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*ssmask(:,:) )
ELSE
CALL iom_put( "ssh" , ssh(:,:,Kmm) ) ! sea surface height
ENDIF
ENDIF
IF( iom_use("wetdep") ) CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) ) ! wet depth
#if defined key_qco
IF( iom_use("ht") ) CALL iom_put( "ht" , ht(:,:) ) ! water column at t-point
IF( iom_use("hu") ) CALL iom_put( "hu" , hu(:,:,Kmm) ) ! water column at u-point
IF( iom_use("hv") ) CALL iom_put( "hv" , hv(:,:,Kmm) ) ! water column at v-point
IF( iom_use("hf") ) CALL iom_put( "hf" , hf_0(:,:)*( 1._wp + r3f(:,:) ) ) ! water column at f-point (caution here at Naa)
#endif
! --- tracers T&S --- !
CALL iom_put( "toce", ts(:,:,:,jp_tem,Kmm) ) ! 3D temperature
CALL iom_put( "sst", ts(:,:,1,jp_tem,Kmm) ) ! surface temperature
IF ( iom_use("sbt") ) THEN
DO_2D( 0, 0, 0, 0 )
ikbot = mbkt(ji,jj)
Clement Rousset
committed
z2d0(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm)
Clement Rousset
committed
CALL iom_put( "sbt", z2d0 ) ! bottom temperature
ENDIF
CALL iom_put( "soce", ts(:,:,:,jp_sal,Kmm) ) ! 3D salinity
CALL iom_put( "sss", ts(:,:,1,jp_sal,Kmm) ) ! surface salinity
IF ( iom_use("sbs") ) THEN
DO_2D( 0, 0, 0, 0 )
ikbot = mbkt(ji,jj)
Clement Rousset
committed
z2d0(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm)
Clement Rousset
committed
CALL iom_put( "sbs", z2d0 ) ! bottom salinity
ENDIF
IF( .NOT.lk_SWE ) CALL iom_put( "rhop", rhop(:,:,:) ) ! 3D potential density (sigma0)
! --- momentum --- !
IF ( iom_use("taubot") ) THEN ! bottom stress
zztmp = rho0 * 0.25_wp
Clement Rousset
committed
z2d0(:,:) = 0._wp
DO_2D( 0, 0, 0, 0 )
zztmp2 = ( ( rCdU_bot(ji+1,jj)+rCdU_bot(ji ,jj) ) * uu(ji ,jj,mbku(ji ,jj),Kmm) )**2 &
& + ( ( rCdU_bot(ji ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm) )**2 &
& + ( ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj ) ) * vv(ji,jj ,mbkv(ji,jj ),Kmm) )**2 &
& + ( ( rCdU_bot(ji,jj )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm) )**2
Clement Rousset
committed
z2d0(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)
Clement Rousset
committed
CALL iom_put( "taubot", z2d0 )
ENDIF
CALL iom_put( "uoce", uu(:,:,:,Kmm) ) ! 3D i-current
CALL iom_put( "ssu", uu(:,:,1,Kmm) ) ! surface i-current
IF ( iom_use("sbu") ) THEN
DO_2D( 0, 0, 0, 0 )
ikbot = mbku(ji,jj)
Clement Rousset
committed
z2d0(ji,jj) = uu(ji,jj,ikbot,Kmm)
Clement Rousset
committed
CALL iom_put( "sbu", z2d0 ) ! bottom i-current
ENDIF
CALL iom_put( "voce", vv(:,:,:,Kmm) ) ! 3D j-current
CALL iom_put( "ssv", vv(:,:,1,Kmm) ) ! surface j-current
IF ( iom_use("sbv") ) THEN
DO_2D( 0, 0, 0, 0 )
ikbot = mbkv(ji,jj)
Clement Rousset
committed
z2d0(ji,jj) = vv(ji,jj,ikbot,Kmm)
Clement Rousset
committed
CALL iom_put( "sbv", z2d0 ) ! bottom j-current
ENDIF
! ! vertical velocity
IF( ln_zad_Aimp ) THEN
IF( iom_use('woce') ) THEN
Clement Rousset
committed
CALL iom_put( "woce", ww+wi ) ! explicit plus implicit parts
ENDIF
ELSE
CALL iom_put( "woce", ww )
ENDIF
IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value
! ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
IF( ln_zad_Aimp ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
Clement Rousset
committed
z3d0(ji,jj,jk) = rho0 * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wi(ji,jj,jk) )
Clement Rousset
committed
z3d0(ji,jj,jk) = rho0 * e1e2t(ji,jj) * ww(ji,jj,jk)
Clement Rousset
committed
CALL iom_put( "w_masstr" , z3d0 )
IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d0 * z3d0 )
ENDIF
CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef.
CALL iom_put( "avs" , avs ) ! S vert. eddy diff. coef.
CALL iom_put( "avm" , avm ) ! T vert. eddy visc. coef.
IF( iom_use('logavt') ) CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
IF( iom_use('logavs') ) CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
IF ( iom_use("sssgrad") .OR. iom_use("sssgrad2") ) THEN
DO_2D( 0, 0, 0, 0 ) ! sss gradient
zztmp = ts(ji,jj,1,jp_sal,Kmm)
zztmpx = (ts(ji+1,jj,1,jp_sal,Kmm) - zztmp) * r1_e1u(ji,jj) + (zztmp - ts(ji-1,jj ,1,jp_sal,Kmm)) * r1_e1u(ji-1,jj)
zztmpy = (ts(ji,jj+1,1,jp_sal,Kmm) - zztmp) * r1_e2v(ji,jj) + (zztmp - ts(ji ,jj-1,1,jp_sal,Kmm)) * r1_e2v(ji,jj-1)
Clement Rousset
committed
z2d0(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * vmask(ji,jj-1,1)
END_2D
Clement Rousset
committed
CALL iom_put( "sssgrad2", z2d0 ) ! square of module of sss gradient
IF ( iom_use("sssgrad") ) THEN
DO_2D( 0, 0, 0, 0 )
Clement Rousset
committed
z2d0(ji,jj) = SQRT( z2d0(ji,jj) )
Clement Rousset
committed
CALL iom_put( "sssgrad", z2d0 ) ! module of sss gradient
ENDIF
ENDIF
IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
DO_2D( 0, 0, 0, 0 ) ! sst gradient
zztmp = ts(ji,jj,1,jp_tem,Kmm)
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)
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)
Clement Rousset
committed
z2d0(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * vmask(ji,jj-1,1)
END_2D
Clement Rousset
committed
CALL iom_put( "sstgrad2", z2d0 ) ! square of module of sst gradient
IF ( iom_use("sstgrad") ) THEN
DO_2D( 0, 0, 0, 0 )
Clement Rousset
committed
z2d0(ji,jj) = SQRT( z2d0(ji,jj) )
Clement Rousset
committed
CALL iom_put( "sstgrad", z2d0 ) ! module of sst gradient
ENDIF
ENDIF
! heat and salt contents
IF( iom_use("heatc") ) THEN
Clement Rousset
committed
z2d0(:,:) = 0._wp
Clement Rousset
committed
z2d0(ji,jj) = z2d0(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk)
Clement Rousset
committed
CALL iom_put( "heatc", rho0_rcp * z2d0 ) ! vertically integrated heat content (J/m2)
Clement Rousset
committed
z2d0(:,:) = 0._wp
Clement Rousset
committed
z2d0(ji,jj) = z2d0(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
Clement Rousset
committed
CALL iom_put( "saltc", rho0 * z2d0 ) ! vertically integrated salt content (PSU*kg/m2)
Clement Rousset
committed
z2d0(:,:) = 0._wp
Clement Rousset
committed
z2d0(ji,jj) = z2d0(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
Clement Rousset
committed
CALL iom_put( "salt2c", rho0 * z2d0 ) ! vertically integrated square of salt content (PSU2*kg/m2)
ENDIF
!
IF ( iom_use("ke") .OR. iom_use("ke_int") ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
zztmpx = uu(ji-1,jj ,jk,Kmm) + uu(ji,jj,jk,Kmm)
zztmpy = vv(ji ,jj-1,jk,Kmm) + vv(ji,jj,jk,Kmm)
Clement Rousset
committed
z3d0(ji,jj,jk) = 0.25_wp * ( zztmpx*zztmpx + zztmpy*zztmpy )
Clement Rousset
committed
CALL iom_put( "ke", z3d0 ) ! kinetic energy
Clement Rousset
committed
z2d0(:,:) = 0._wp
Clement Rousset
committed
z2d0(ji,jj) = z2d0(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d0(ji,jj,jk) * e1e2t(ji,jj) * tmask(ji,jj,jk)
Clement Rousset
committed
CALL iom_put( "ke_int", z2d0 ) ! vertically integrated kinetic energy
ENDIF
!
IF ( iom_use("sKE") ) THEN ! surface kinetic energy at T point
Clement Rousset
committed
z2d0(:,:) = 0._wp
Clement Rousset
committed
z2d0(ji,jj) = 0.25_wp * ( uu(ji ,jj,1,Kmm) * uu(ji ,jj,1,Kmm) * e1e2u(ji ,jj) * e3u(ji ,jj,1,Kmm) &
& + uu(ji-1,jj,1,Kmm) * uu(ji-1,jj,1,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,1,Kmm) &
& + vv(ji,jj ,1,Kmm) * vv(ji,jj ,1,Kmm) * e1e2v(ji,jj ) * e3v(ji,jj ,1,Kmm) &
& + vv(ji,jj-1,1,Kmm) * vv(ji,jj-1,1,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,1,Kmm) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,1,Kmm) * ssmask(ji,jj)
END_2D
Clement Rousset
committed
IF ( iom_use("sKE" ) ) CALL iom_put( "sKE" , z2d0 )
ENDIF
!
IF ( iom_use("ssKEf") ) THEN ! surface kinetic energy at F point
Clement Rousset
committed
z2d0(:,:) = 0._wp ! CAUTION : only valid in SWE, not with bathymetry
Clement Rousset
committed
z2d0(ji,jj) = 0.25_wp * ( uu(ji,jj ,1,Kmm) * uu(ji,jj ,1,Kmm) * e1e2u(ji,jj ) * e3u(ji,jj ,1,Kmm) &
& + uu(ji,jj+1,1,Kmm) * uu(ji,jj+1,1,Kmm) * e1e2u(ji,jj+1) * e3u(ji,jj+1,1,Kmm) &
& + vv(ji ,jj,1,Kmm) * vv(ji,jj ,1,Kmm) * e1e2v(ji ,jj) * e3v(ji ,jj,1,Kmm) &
& + vv(ji+1,jj,1,Kmm) * vv(ji+1,jj,1,Kmm) * e1e2v(ji+1,jj) * e3v(ji+1,jj,1,Kmm) ) &
& * r1_e1e2f(ji,jj) / e3f(ji,jj,1) * ssfmask(ji,jj)
END_2D
Clement Rousset
committed
CALL iom_put( "ssKEf", z2d0 )
ENDIF
!
CALL iom_put( "hdiv", hdiv ) ! Horizontal divergence
!
IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
Clement Rousset
committed
z3d0(ji,jj,jk) = rho0 * uu(ji,jj,jk,Kmm) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
Clement Rousset
committed
CALL iom_put( "u_masstr" , z3d0 ) ! mass transport in i-direction
Clement Rousset
committed
z2d0(:,:) = 0._wp
Clement Rousset
committed
z2d0(ji,jj) = z2d0(ji,jj) + z3d0(ji,jj,jk)
Clement Rousset
committed
CALL iom_put( "u_masstr_vint", z2d0 ) ! mass transport in i-direction vertical sum
Clement Rousset
committed
z2d0(:,:) = 0._wp
zztmp = 0.5_wp * rcp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
Clement Rousset
committed
z2d0(ji,jj) = z2d0(ji,jj) + zztmp * z3d0(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) )
Clement Rousset
committed
CALL iom_put( "u_heattr", z2d0 ) ! heat transport in i-direction
Clement Rousset
committed
z2d0(:,:) = 0._wp
Clement Rousset
committed
z2d0(ji,jj) = z2d0(ji,jj) + 0.5 * z3d0(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) )
Clement Rousset
committed
CALL iom_put( "u_salttr", z2d0 ) ! heat transport in i-direction
ENDIF
ENDIF
IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
Clement Rousset
committed
z3d0(ji,jj,jk) = rho0 * vv(ji,jj,jk,Kmm) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
Clement Rousset
committed
CALL iom_put( "v_masstr", z3d0 ) ! mass transport in j-direction
Clement Rousset
committed
z2d0(:,:) = 0._wp
zztmp = 0.5_wp * rcp
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
Clement Rousset
committed
z2d0(ji,jj) = z2d0(ji,jj) + zztmp * z3d0(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) )
Clement Rousset
committed
CALL iom_put( "v_heattr", z2d0 ) ! heat transport in j-direction
Clement Rousset
committed
z2d0(:,:) = 0._wp
Clement Rousset
committed
z2d0(ji,jj) = z2d0(ji,jj) + 0.5 * z3d0(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) )
Clement Rousset
committed
CALL iom_put( "v_salttr", z2d0 ) ! heat transport in j-direction
ENDIF
ENDIF
IF( iom_use("tosmint") ) THEN
Clement Rousset
committed
z2d0(:,:) = 0._wp
Clement Rousset
committed
z2d0(ji,jj) = z2d0(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)
Clement Rousset
committed
CALL iom_put( "tosmint", z2d0 ) ! Vertical integral of temperature
Clement Rousset
committed
z2d0(:,:) = 0._wp
Clement Rousset
committed
z2d0(ji,jj) = z2d0(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
Clement Rousset
committed
CALL iom_put( "somint", z2d0 ) ! Vertical integral of salinity
ENDIF
CALL iom_put( "bn2", rn2 ) ! Brunt-Vaisala buoyancy frequency (N^2)
IF (ln_dia25h) CALL dia_25h( kt, Kmm ) ! 25h averaging
! Output of surface vorticity terms
!
Clement Rousset
committed
IF ( iom_use("ssplavor") .OR. iom_use("ssrelvor") .OR. iom_use("ssEns") .OR. &
& iom_use("ssrelpotvor") .OR. iom_use("ssabspotvor") ) THEN
!
DO_2D( 0, 0, 0, 0 )
Clement Rousset
committed
z2d0(ji,jj) = ff_f(ji,jj)
END_2D
CALL iom_put( "ssplavor", z2d0 ) ! planetary vorticity ( f )
DO_2D( 0, 0, 0, 0 )
z2d0(ji,jj) = ( e2v(ji+1,jj ) * vv(ji+1,jj ,1,Kmm) - e2v(ji,jj) * vv(ji,jj,1,Kmm) &
& - e1u(ji ,jj+1) * uu(ji ,jj+1,1,Kmm) + e1u(ji,jj) * uu(ji,jj,1,Kmm) ) * r1_e1e2f(ji,jj)
END_2D
Clement Rousset
committed
CALL iom_put( "ssrelvor", z2d0 ) ! relative vorticity ( zeta )
!
IF ( iom_use("ssEns") .OR. iom_use("ssrelpotvor") .OR. iom_use("ssabspotvor") ) THEN
DO_2D( 0, 0, 0, 0 )
ze3 = ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1) &
Clement Rousset
committed
& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3
ELSE ; ze3 = 0._wp
ENDIF
Clement Rousset
committed
z2d0(ji,jj) = ze3 * z2d0(ji,jj)
Clement Rousset
committed
CALL iom_put( "ssrelpotvor", z2d0 ) ! relative potential vorticity (zeta/h)
!
IF ( iom_use("ssEns") .OR. iom_use("ssabspotvor") ) THEN
DO_2D( 0, 0, 0, 0 )
ze3 = ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1) &
Clement Rousset
committed
& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3
ELSE ; ze3 = 0._wp
ENDIF
Clement Rousset
committed
z2d0(ji,jj) = ze3 * ff_f(ji,jj) + z2d0(ji,jj)
Clement Rousset
committed
CALL iom_put( "ssabspotvor", z2d0 ) ! absolute potential vorticity ( q )
!
IF ( iom_use("ssEns") ) THEN
DO_2D( 0, 0, 0, 0 )
Clement Rousset
committed
z2d0(ji,jj) = 0.5_wp * z2d0(ji,jj) * z2d0(ji,jj)
Clement Rousset
committed
CALL iom_put( "ssEns", z2d0 ) ! potential enstrophy ( 1/2*q2 )
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
ENDIF
ENDIF
ENDIF
ENDIF
IF( ln_timing ) CALL timing_stop('dia_wri')
!
END SUBROUTINE dia_wri
#else
!!----------------------------------------------------------------------
!! Default option use IOIPSL library
!!----------------------------------------------------------------------
INTEGER FUNCTION dia_wri_alloc()
!!----------------------------------------------------------------------
INTEGER, DIMENSION(2) :: ierr
!!----------------------------------------------------------------------
IF( nn_write == -1 ) THEN
dia_wri_alloc = 0
ELSE
ierr = 0
ALLOCATE( ndex_hT(jpi*jpj) , ndex_T(jpi*jpj*jpk) , &
& ndex_hU(jpi*jpj) , ndex_U(jpi*jpj*jpk) , &
& ndex_hV(jpi*jpj) , ndex_V(jpi*jpj*jpk) , STAT=ierr(1) )
!
dia_wri_alloc = MAXVAL(ierr)
CALL mpp_sum( 'diawri', dia_wri_alloc )
!
ENDIF
!
END FUNCTION dia_wri_alloc
INTEGER FUNCTION dia_wri_alloc_abl()
!!----------------------------------------------------------------------
ALLOCATE( ndex_hA(jpi*jpj), ndex_A (jpi*jpj*jpkam1), STAT=dia_wri_alloc_abl)
CALL mpp_sum( 'diawri', dia_wri_alloc_abl )
!
END FUNCTION dia_wri_alloc_abl
SUBROUTINE dia_wri( kt, Kmm )
!!---------------------------------------------------------------------
!! *** ROUTINE dia_wri ***
!!
!! ** Purpose : Standard output of opa: dynamics and tracer fields
!! NETCDF format is used by default
!!
!! ** Method : At the beginning of the first time step (nit000),
!! define all the NETCDF files and fields
!! At each time step call histdef to compute the mean if ncessary
!! Each nn_write time step, output the instantaneous or mean fields
!!----------------------------------------------------------------------
INTEGER, INTENT( in ) :: kt ! ocean time-step index
INTEGER, INTENT( in ) :: Kmm ! ocean time level index
!
LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
CHARACTER (len=40) :: clhstnam, clop, clmx ! local names
INTEGER :: inum = 11 ! temporary logical unit
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ierr ! error code return from allocation
INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers
INTEGER :: ipka ! ABL
INTEGER :: jn, ierror ! local integers
REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars
!
Clement Rousset
committed
REAL(wp), DIMENSION(jpi,jpj ) :: z2d0 ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d0 ! 3D workspace
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl ! ABL 3D workspace
!!----------------------------------------------------------------------
!
IF( ninist == 1 ) THEN !== Output the initial state and forcings ==!
CALL dia_wri_state( Kmm, 'output.init' )
ninist = 0
ENDIF
!
IF( nn_write == -1 ) RETURN ! we will never do any output
!
IF( ln_timing ) CALL timing_start('dia_wri')
!
! 0. Initialisation
! -----------------
ll_print = .FALSE. ! local variable for debugging
ll_print = ll_print .AND. lwp
! Define frequency of output and means
clop = "x" ! no use of the mask value (require less cpu time and otherwise the model crashes)
#if defined key_diainstant
zsto = nn_write * rn_Dt
clop = "inst("//TRIM(clop)//")"
#else
zsto=rn_Dt
clop = "ave("//TRIM(clop)//")"
#endif
zout = nn_write * rn_Dt
zmax = ( nitend - nit000 + 1 ) * rn_Dt
! Define indices of the horizontal output zoom and vertical limit storage
iimi = Nis0 ; iima = Nie0
ijmi = Njs0 ; ijma = Nje0
ipk = jpk
IF(ln_abl) ipka = jpkam1
! define time axis
it = kt
itmod = kt - nit000 + 1
! 1. Define NETCDF files and fields at beginning of first time step
! -----------------------------------------------------------------
IF( kt == nit000 ) THEN
! Define the NETCDF files (one per grid)
! Compute julian date from starting date of the run
CALL ymds2ju( nyear, nmonth, nday, rn_Dt, zjulian )
zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
IF(lwp)WRITE(numout,*)
IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear, &
& ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma, &
' limit storage in depth = ', ipk
! WRITE root name in date.file for use by postpro
IF(lwp) THEN
CALL dia_nam( clhstnam, nn_write,' ' )
CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
WRITE(inum,*) clhstnam
CLOSE(inum)
ENDIF
! Define the T grid FILE ( nid_T )
CALL dia_nam( clhstnam, nn_write, 'grid_T' )
IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
& nit000-1, zjulian, rn_Dt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
& "m", ipk, gdept_1d, nz_T, "down" )
! ! Index of ocean points
CALL wheneq( jpi*jpj*ipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume
CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
!
IF( ln_icebergs ) THEN
!
!! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
!! that routine is called from nemogcm, so do it here immediately before its needed
ALLOCATE( ndex_bT(jpi*jpj*nclasses), STAT=ierror )
CALL mpp_sum( 'diawri', ierror )
IF( ierror /= 0 ) THEN
CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
RETURN
ENDIF
!
!! iceberg vertical coordinate is class number
CALL histvert( nid_T, "class", "Iceberg class", & ! Vertical grid: class
& "number", nclasses, class_num, nb_T )
!
!! each class just needs the surface index pattern
ndim_bT = 3
DO jn = 1,nclasses
ndex_bT((jn-1)*jpi*jpj+1:jn*jpi*jpj) = ndex_hT(1:jpi*jpj)
ENDDO
!
ENDIF
! Define the U grid FILE ( nid_U )
CALL dia_nam( clhstnam, nn_write, 'grid_U' )
IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
& nit000-1, zjulian, rn_Dt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept
& "m", ipk, gdept_1d, nz_U, "down" )
! ! Index of ocean points
CALL wheneq( jpi*jpj*ipk, umask, 1, 1., ndex_U , ndim_U ) ! volume
CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
! Define the V grid FILE ( nid_V )
CALL dia_nam( clhstnam, nn_write, 'grid_V' ) ! filename
IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
& nit000-1, zjulian, rn_Dt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept
& "m", ipk, gdept_1d, nz_V, "down" )
! ! Index of ocean points
CALL wheneq( jpi*jpj*ipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume
CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
! Define the W grid FILE ( nid_W )
CALL dia_nam( clhstnam, nn_write, 'grid_W' ) ! filename
IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
& nit000-1, zjulian, rn_Dt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw
& "m", ipk, gdepw_1d, nz_W, "down" )
IF( ln_abl ) THEN
! Define the ABL grid FILE ( nid_A )
CALL dia_nam( clhstnam, nn_write, 'grid_ABL' )
IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
& nit000-1, zjulian, rn_Dt, nh_A, nid_A, domain_id=nidom, snc4chunks=snc4set )
CALL histvert( nid_A, "ght_abl", "Vertical T levels", & ! Vertical grid: gdept
& "m", ipka, ght_abl(2:jpka), nz_A, "up" )
! ! Index of ocean points
ALLOCATE( zw3d_abl(jpi,jpj,ipka) )
zw3d_abl(:,:,:) = 1._wp
CALL wheneq( jpi*jpj*ipka, zw3d_abl, 1, 1., ndex_A , ndim_A ) ! volume
CALL wheneq( jpi*jpj , zw3d_abl, 1, 1., ndex_hA, ndim_hA ) ! surface
DEALLOCATE(zw3d_abl)
ENDIF
!
! Declare all the output fields as NETCDF variables
! !!! nid_T : 3D
CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn
& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn
& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
IF( .NOT.ln_linssh ) THEN
CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t n
& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t n
& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t n
& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
ENDIF
! !!! nid_T : 2D
CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "sosaline", "Sea Surface Salinity" , "PSU" , & ! sss
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! ssh
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux" , "Kg/m2/s", & ! (emp-rnf)
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "sorunoff", "River runoffs" , "Kg/m2/s", & ! runoffs
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
IF( ln_linssh ) THEN
CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * ts(:,:,1,jp_tem,Kmm)
& , "KgC/m2/s", & ! sosst_cd
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * ts(:,:,1,jp_sal,Kmm)
& , "KgPSU/m2/s",& ! sosss_cd
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
ENDIF
CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! qsr
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
IF( ALLOCATED(hmld) ) THEN ! zdf_mxl not called by SWE
CALL histdef( nid_T, "somixhgt", "Turbocline Depth" , "m" , & ! hmld
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01" , "m" , & ! hmlp
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
ENDIF
CALL histdef( nid_T, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
!
IF( ln_abl ) THEN
CALL histdef( nid_A, "t_abl", "Potential Temperature" , "K" , & ! t_abl
& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
CALL histdef( nid_A, "q_abl", "Humidity" , "kg/kg" , & ! q_abl
& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
CALL histdef( nid_A, "u_abl", "Atmospheric U-wind " , "m/s" , & ! u_abl
& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
CALL histdef( nid_A, "v_abl", "Atmospheric V-wind " , "m/s" , & ! v_abl
& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
CALL histdef( nid_A, "tke_abl", "Atmospheric TKE " , "m2/s2" , & ! tke_abl
& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
CALL histdef( nid_A, "avm_abl", "Atmospheric turbulent viscosity", "m2/s" , & ! avm_abl
& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
CALL histdef( nid_A, "avt_abl", "Atmospheric turbulent diffusivity", "m2/s2", & ! avt_abl
& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
CALL histdef( nid_A, "pblh", "Atmospheric boundary layer height " , "m", & ! pblh
& jpi, jpj, nh_A, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
#if defined key_si3
CALL histdef( nid_A, "oce_frac", "Fraction of open ocean" , " ", & ! ato_i
& jpi, jpj, nh_A, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
#endif
CALL histend( nid_A, snc4chunks=snc4set )
ENDIF
!
IF( ln_icebergs ) THEN
CALL histdef( nid_T, "calving" , "calving mass input" , "kg/s" , &
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "calving_heat" , "calving heat flux" , "XXXX" , &
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "berg_floating_melt" , "Melt rate of icebergs + bits" , "kg/m2/s", &
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "berg_stored_ice" , "Accumulated ice mass by class" , "kg" , &
& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
IF( ln_bergdia ) THEN
CALL histdef( nid_T, "berg_melt" , "Melt rate of icebergs" , "kg/m2/s", &
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "berg_buoy_melt" , "Buoyancy component of iceberg melt rate" , "kg/m2/s", &
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "berg_eros_melt" , "Erosion component of iceberg melt rate" , "kg/m2/s", &
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "berg_conv_melt" , "Convective component of iceberg melt rate", "kg/m2/s", &
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "berg_virtual_area" , "Virtual coverage by icebergs" , "m2" , &
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "bits_src" , "Mass source of bergy bits" , "kg/m2/s", &
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "bits_melt" , "Melt rate of bergy bits" , "kg/m2/s", &
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "bits_mass" , "Bergy bit density field" , "kg/m2" , &
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "berg_mass" , "Iceberg density field" , "kg/m2" , &
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "berg_real_calving" , "Calving into iceberg class" , "kg/s" , &
& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
ENDIF
ENDIF
IF( ln_ssr ) THEN
CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping" , "Kg/m2/s", & ! erp * sn
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
ENDIF
clmx ="l_max(only(x))" ! max index on a period
! CALL histdef( nid_T, "sobowlin", "Bowl Index" , "W-point", & ! bowl INDEX
! & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clmx, zsto, zout )
#if defined key_diahth
CALL histdef( nid_T, "sothedep", "Thermocline Depth" , "m" , & ! hth
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm" , "m" , & ! hd20
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm" , "m" , & ! hd28
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "J/m2" , & ! htc3
& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
#endif
CALL histdef( nid_T, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
& jpi, jpj, nh_T, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
CALL histdef( nid_T, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
& jpi, jpj, nh_T, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
!
CALL histend( nid_T, snc4chunks=snc4set )
! !!! nid_U : 3D
CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! uu(:,:,:,Kmm)
& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
IF( ln_wave .AND. ln_sdw) THEN
CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current" , "m/s" , & ! usd
& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
ENDIF
CALL histend( nid_U, snc4chunks=snc4set )
! !!! nid_V : 3D
CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vv(:,:,:,Kmm)
& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
IF( ln_wave .AND. ln_sdw) THEN
CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current" , "m/s" , & ! vsd
& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
ENDIF
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
CALL histend( nid_V, snc4chunks=snc4set )
! !!! nid_W : 3D
CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! ww
& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt
& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avm
& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
IF( ln_zdfddm ) THEN
CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity" , "m2/s" , & ! avs
& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
ENDIF
IF( ln_wave .AND. ln_sdw) THEN
CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current" , "m/s" , & ! wsd
& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
ENDIF
! !!! nid_W : 2D
CALL histend( nid_W, snc4chunks=snc4set )
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
IF(ll_print) CALL FLUSH(numout )
ENDIF
! 2. Start writing data
! ---------------------
! ndex(1) est utilise ssi l'avant dernier argument est different de
! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
! donne le nombre d'elements, et ndex la liste des indices a sortir
IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN
WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
WRITE(numout,*) '~~~~~~ '
ENDIF
IF( .NOT.ln_linssh ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
Clement Rousset
committed
z3d0(ji,jj,jk) = ts(ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm)
Clement Rousset
committed
CALL histwrite( nid_T, "votemper", it, z3d0, ndim_T , ndex_T ) ! heat content
Clement Rousset
committed
z3d0(ji,jj,jk) = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm)
Clement Rousset
committed
CALL histwrite( nid_T, "vosaline", it, z3d0, ndim_T , ndex_T ) ! salt content
Clement Rousset
committed
z2d0(ji,jj ) = ts(ji,jj, 1,jp_tem,Kmm) * e3t(ji,jj, 1,Kmm)
Clement Rousset
committed
CALL histwrite( nid_T, "sosstsst", it, z2d0, ndim_hT, ndex_hT ) ! sea surface heat content
Clement Rousset
committed
z2d0(ji,jj ) = ts(ji,jj, 1,jp_sal,Kmm) * e3t(ji,jj, 1,Kmm)
Clement Rousset
committed
CALL histwrite( nid_T, "sosaline", it, z2d0, ndim_hT, ndex_hT ) ! sea surface salinity content
ELSE
CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T ) ! temperature
CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) , ndim_T , ndex_T ) ! salinity
CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) , ndim_hT, ndex_hT ) ! sea surface temperature
CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) , ndim_hT, ndex_hT ) ! sea surface salinity
ENDIF
IF( .NOT.ln_linssh ) THEN
DO_3D( 0, 0, 0, 0, 1, jpk )
Clement Rousset
committed
z3d0(ji,jj,jk) = e3t(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
Clement Rousset
committed
CALL histwrite( nid_T, "vovvle3t", it, z3d0 , ndim_T , ndex_T ) ! level thickness
Clement Rousset
committed
z3d0(ji,jj,jk) = gdept(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
Clement Rousset
committed
CALL histwrite( nid_T, "vovvldep", it, z3d0 , ndim_T , ndex_T ) ! t-point depth
Clement Rousset
committed
z3d0(ji,jj,jk) = ( ( e3t(ji,jj,jk,Kmm) - e3t_0(ji,jj,jk) ) / e3t_0(ji,jj,jk) * 100._wp * tmask(ji,jj,jk) ) ** 2
Clement Rousset
committed
CALL histwrite( nid_T, "vovvldef", it, z3d0 , ndim_T , ndex_T ) ! level thickness deformation
ENDIF
CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm) , ndim_hT, ndex_hT ) ! sea surface height
DO_2D( 0, 0, 0, 0 )
Clement Rousset
committed
z2d0(ji,jj) = emp(ji,jj) - rnf(ji,jj)
Clement Rousset
committed
CALL histwrite( nid_T, "sowaflup", it, z2d0 , ndim_hT, ndex_hT ) ! upward water flux
CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
! (includes virtual salt flux beneath ice
! in linear free surface case)
IF( ln_linssh ) THEN
DO_2D( 0, 0, 0, 0 )
Clement Rousset
committed
z2d0(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_tem,Kmm)
Clement Rousset
committed
CALL histwrite( nid_T, "sosst_cd", it, z2d0, ndim_hT, ndex_hT ) ! c/d term on sst
Clement Rousset
committed
z2d0(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_sal,Kmm)
Clement Rousset
committed
CALL histwrite( nid_T, "sosss_cd", it, z2d0, ndim_hT, ndex_hT ) ! c/d term on sss
Clement Rousset
committed
z2d0(ji,jj) = qsr(ji,jj) + qns(ji,jj)
Clement Rousset
committed
CALL histwrite( nid_T, "sohefldo", it, z2d0 , ndim_hT, ndex_hT ) ! total heat flux
CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
IF( ALLOCATED(hmld) ) THEN ! zdf_mxl not called by SWE
CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
CALL histwrite( nid_T, "somxl010", it, hmlp , ndim_hT, ndex_hT ) ! mixed layer depth