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
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
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
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
MODULE dynspg
!!======================================================================
!! *** MODULE dynspg ***
!! Ocean dynamics: surface pressure gradient control
!!======================================================================
!! History : 1.0 ! 2005-12 (C. Talandier, G. Madec, V. Garnier) Original code
!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
!! 4.2 ! 2020-12 (G. Madec, E. Clementi) add Bernoulli Head for
!! wave coupling
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! dyn_spg : update the dynamics trend with surface pressure gradient
!! dyn_spg_init: initialization, namelist read, and parameters control
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers variables
USE dom_oce ! ocean space and time domain variables
USE c1d ! 1D vertical configuration
USE phycst ! physical constants
USE sbc_oce ! surface boundary condition: ocean
USE sbc_ice , ONLY : snwice_mass, snwice_mass_b
USE sbcapr ! surface boundary condition: atmospheric pressure
USE sbcwave, ONLY : bhd_wave
USE dynspg_exp ! surface pressure gradient (dyn_spg_exp routine)
USE dynspg_ts ! surface pressure gradient (dyn_spg_ts routine)
USE tide_mod !
USE trd_oce ! trends: ocean variables
USE trddyn ! trend manager: dynamics
!
USE prtctl ! Print control (prt_ctl routine)
USE in_out_manager ! I/O manager
USE lib_mpp ! MPP library
USE timing ! Timing
IMPLICIT NONE
PRIVATE
PUBLIC dyn_spg ! routine called by step module
PUBLIC dyn_spg_init ! routine called by opa module
INTEGER :: nspg = 0 ! type of surface pressure gradient scheme defined from lk_dynspg_...
! ! Parameter to control the surface pressure gradient scheme
INTEGER, PARAMETER :: np_TS = 1 ! split-explicit time stepping (Time-Splitting)
INTEGER, PARAMETER :: np_EXP = 0 ! explicit time stepping
INTEGER, PARAMETER :: np_NO =-1 ! no surface pressure gradient, no scheme
!
REAL(wp) :: zt0step ! Time of day at the beginning of the time step
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynspg.F90 14225 2020-12-19 14:58:39Z smasson $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_spg( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa, k_only_ADV )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_spg ***
!!
!! ** Purpose : compute surface pressure gradient including the
!! atmospheric pressure forcing (ln_apr_dyn=T).
!!
!! ** Method : Two schemes:
!! - explicit : the spg is evaluated at now
!! - split-explicit : a time splitting technique is used
!!
!! ln_apr_dyn=T : the atmospheric pressure forcing is applied
!! as the gradient of the inverse barometer ssh:
!! apgu = - 1/rho0 di[apr] = 0.5*grav di[ssh_ib+ssh_ibb]
!! apgv = - 1/rho0 dj[apr] = 0.5*grav dj[ssh_ib+ssh_ibb]
!! Note that as all external forcing a time averaging over a two rn_Dt
!! period is used to prevent the divergence of odd and even time step.
!!----------------------------------------------------------------------
INTEGER , INTENT( in ) :: kt ! ocean time-step index
INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(inout) :: pssh, puu_b, pvv_b ! SSH and barotropic velocities at main time levels
INTEGER , OPTIONAL , INTENT( in ) :: k_only_ADV ! only Advection in the RHS
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: z2dt, zg_2, zintp, zgrho0r, zld ! local scalars
REAL(wp) , DIMENSION(jpi,jpj) :: zpgu, zpgv ! 2D workspace
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zpice
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dyn_spg')
!
IF( l_trddyn ) THEN ! temporary save of ta and sa trends
ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
ztrdu(:,:,:) = puu(:,:,:,Krhs)
ztrdv(:,:,:) = pvv(:,:,:,Krhs)
ENDIF
!
IF( ln_apr_dyn & ! atmos. pressure
.OR. ( .NOT.ln_dynspg_ts .AND. (ln_tide_pot .AND. ln_tide) ) & ! tide potential (no time slitting)
.OR. ln_ice_embd & ! embedded sea-ice
.OR. ( ln_wave .and. ln_bern_srfc ) ) THEN ! depth-independent Bernoulli head
!
DO_2D( 0, 0, 0, 0 )
zpgu(ji,jj) = 0._wp
zpgv(ji,jj) = 0._wp
END_2D
!
IF( ln_apr_dyn .AND. .NOT.ln_dynspg_ts ) THEN !== Atmospheric pressure gradient (added later in time-split case) ==!
zg_2 = grav * 0.5
DO_2D( 0, 0, 0, 0 ) ! gradient of Patm using inverse barometer ssh
zpgu(ji,jj) = zpgu(ji,jj) + zg_2 * ( ssh_ib (ji+1,jj) - ssh_ib (ji,jj) &
& + ssh_ibb(ji+1,jj) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj)
zpgv(ji,jj) = zpgv(ji,jj) + zg_2 * ( ssh_ib (ji,jj+1) - ssh_ib (ji,jj) &
& + ssh_ibb(ji,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj)
END_2D
ENDIF
!
! !== tide potential forcing term ==!
IF( .NOT.ln_dynspg_ts .AND. ( ln_tide_pot .AND. ln_tide ) ) THEN ! N.B. added directly at sub-time-step in ts-case
!
! Update tide potential at the beginning of current time step
zt0step = REAL(nsec_day, wp)-0.5_wp*rn_Dt
CALL upd_tide(zt0step, Kmm)
!
DO_2D( 0, 0, 0, 0 ) ! add tide potential forcing
zpgu(ji,jj) = zpgu(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj)
zpgv(ji,jj) = zpgv(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj)
END_2D
!
IF (ln_scal_load) THEN
zld = rn_scal_load * grav
DO_2D( 0, 0, 0, 0 ) ! add scalar approximation for load potential
zpgu(ji,jj) = zpgu(ji,jj) + zld * ( pssh(ji+1,jj,Kmm) - pssh(ji,jj,Kmm) ) * r1_e1u(ji,jj)
zpgv(ji,jj) = zpgv(ji,jj) + zld * ( pssh(ji,jj+1,Kmm) - pssh(ji,jj,Kmm) ) * r1_e2v(ji,jj)
END_2D
ENDIF
ENDIF
!
IF( ln_ice_embd ) THEN !== embedded sea ice: Pressure gradient due to snow-ice mass ==!
ALLOCATE( zpice(jpi,jpj) )
zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc )
zgrho0r = - grav * r1_rho0
zpice(:,:) = ( zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:) ) * zgrho0r
DO_2D( 0, 0, 0, 0 )
zpgu(ji,jj) = zpgu(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
zpgv(ji,jj) = zpgv(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
END_2D
DEALLOCATE( zpice )
ENDIF
!
IF( ln_wave .and. ln_bern_srfc ) THEN !== Add J terms: depth-independent Bernoulli head
DO_2D( 0, 0, 0, 0 )
zpgu(ji,jj) = zpgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) / e1u(ji,jj) !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3]
zpgv(ji,jj) = zpgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) / e2v(ji,jj)
END_2D
ENDIF
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !== Add all terms to the general trend
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zpgu(ji,jj)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zpgv(ji,jj)
END_3D
!
!!gm add here a call to dyn_trd for ice pressure gradient, the surf pressure trends ????
!
ENDIF
!
SELECT CASE ( nspg ) !== surface pressure gradient computed and add to the general trend ==!
CASE ( np_EXP ) ; CALL dyn_spg_exp( kt, Kmm, puu, pvv, Krhs ) ! explicit
CASE ( np_TS ) ; CALL dyn_spg_ts ( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa, k_only_ADV ) ! time-splitting
END SELECT
!
IF( l_trddyn ) THEN ! save the surface pressure gradient trends for further diagnostics
ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
CALL trd_dyn( ztrdu, ztrdv, jpdyn_spg, kt, Kmm )
DEALLOCATE( ztrdu , ztrdv )
ENDIF
! ! print mean trends (used for debugging)
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' spg - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
IF( ln_timing ) CALL timing_stop('dyn_spg')
!
END SUBROUTINE dyn_spg
SUBROUTINE dyn_spg_init
!!---------------------------------------------------------------------
!! *** ROUTINE dyn_spg_init ***
!!
!! ** Purpose : Control the consistency between namelist options for
!! surface pressure gradient schemes
!!----------------------------------------------------------------------
INTEGER :: ioptio, ios ! local integers
!
NAMELIST/namdyn_spg/ ln_dynspg_exp , ln_dynspg_ts, &
& ln_bt_fw, ln_bt_av , ln_bt_auto , &
& nn_e , rn_bt_cmax, nn_bt_flt, rn_bt_alpha
!!----------------------------------------------------------------------
!
IF(lwp) THEN
WRITE(numout,*)
WRITE(numout,*) 'dyn_spg_init : choice of the surface pressure gradient scheme'
WRITE(numout,*) '~~~~~~~~~~~~'
ENDIF
!
READ ( numnam_ref, namdyn_spg, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_spg in reference namelist' )
!
READ ( numnam_cfg, namdyn_spg, IOSTAT = ios, ERR = 902 )
902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_spg in configuration namelist' )
IF(lwm) WRITE ( numond, namdyn_spg )
!
IF(lwp) THEN ! Namelist print
WRITE(numout,*) ' Namelist : namdyn_spg '
WRITE(numout,*) ' Explicit free surface ln_dynspg_exp = ', ln_dynspg_exp
WRITE(numout,*) ' Free surface with time splitting ln_dynspg_ts = ', ln_dynspg_ts
ENDIF
! ! Control of surface pressure gradient scheme options
nspg = np_NO ; ioptio = 0
IF( ln_dynspg_exp ) THEN ; nspg = np_EXP ; ioptio = ioptio + 1 ; ENDIF
IF( ln_dynspg_ts ) THEN ; nspg = np_TS ; ioptio = ioptio + 1 ; ENDIF
!
IF( ioptio > 1 ) CALL ctl_stop( 'Choose only one surface pressure gradient scheme' )
IF( ioptio == 0 ) CALL ctl_warn( 'NO surface pressure gradient trend in momentum Eqs.' )
IF( ln_dynspg_exp .AND. ln_isfcav ) &
& CALL ctl_stop( ' dynspg_exp not tested with ice shelf cavity ' )
!
IF(lwp) THEN
WRITE(numout,*)
IF( nspg == np_EXP ) WRITE(numout,*) ' ==>>> explicit free surface'
IF( nspg == np_TS ) WRITE(numout,*) ' ==>>> free surface with time splitting scheme'
IF( nspg == np_NO ) WRITE(numout,*) ' ==>>> No surface surface pressure gradient trend in momentum Eqs.'
ENDIF
!
IF( nspg == np_TS ) THEN ! split-explicit scheme initialisation
CALL dyn_spg_ts_init ! do it first: set nn_e used to allocate some arrays later on
ENDIF
!
END SUBROUTINE dyn_spg_init
!!======================================================================
END MODULE dynspg