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
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
MODULE stp2d
!!======================================================================
!! *** MODULE stp2d ***
!! Time-stepping : manager of the ocean, tracer and ice time stepping
!! using a 3rd order Rung Kuta with fixed or quasi-eulerian coordinate
!!======================================================================
!! History : 4.5 ! 2021-01 (S. Techene, G. Madec, N. Ducousso, F. Lemarie) Original code
!! NEMO
!!----------------------------------------------------------------------
#if defined key_qco || defined key_linssh
!!----------------------------------------------------------------------
!! 'key_qco' Quasi-Eulerian vertical coordinate
!! OR
!! 'key_linssh Fixed in time vertical coordinate
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! stp_2D : RK3 case
!!----------------------------------------------------------------------
USE step_oce ! time stepping used modules
USE domqco ! quasi-eulerian coordinate (dom_qco_r3c routine)
USE dynadv_cen2 ! centred flux form advection (dyn_adv_cen2 routine)
USE dynadv_ubs ! UBS flux form advection (dyn_adv_ubs routine)
USE dynkeg ! kinetic energy gradient (dyn_keg routine)
USE dynspg_ts ! 2D mode integration
USE sbc_ice , ONLY : snwice_mass, snwice_mass_b
USE sbcapr ! surface boundary condition: atmospheric pressure
USE sbcwave, ONLY : bhd_wave
#if defined key_agrif
USE agrif_oce_interp
USE agrif_oce_sponge
#endif
PRIVATE
PUBLIC stp_2D ! called by nemogcm.F90
REAL (wp) :: r1_2 = 0.5_wp
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: step.F90 12377 2020-02-12 14:39:06Z acc $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE stp_2D( kt, Kbb, Kmm, Kaa, Krhs )
!!----------------------------------------------------------------------
!! *** ROUTINE stp_2D ***
!!
!! ** Purpose : - Compute sea-surface height and barotropic velocity at Kaa
!! in single 1st RK3.
!!
!! ** Method : -1- Compute the 3D to 2D forcing
!! * Momentum (Ue,Ve)_rhs :
!! 3D to 2D dynamics, i.e. the vertical sum of :
!! - Hor. adv. : KEG + RVO in vector form
!! : ADV_h + MET in flux form
!! - LDF Lateral mixing
!! - HPG Hor. pressure gradient
!! External forcings
!! - baroclinic drag
!! - wind
!! - atmospheric pressure
!! - snow+ice load
!! - surface wave load
!! * ssh (sshe_rhs) :
!! Net column average freshwater flux
!!
!! -2- Solve the external mode Eqs. using sub-time step
!! by a call to dyn_spg_ts (will be renamed dyn_2D or stp_2D)
!!
!! ** action : ssh : N+1 sea surface height (Kaa=N+1)
!! (uu_b,vv_b) : N+1 barotropic velocity
!! (un_adv,vn_adv): barotropic transport from N to N+1
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: kt, Kbb, Kmm, Kaa, Krhs ! ocean time-step and time-level indices
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zg_2, zintp, zgrho0r, zld, zztmp ! local scalars
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zpice ! 2D workspace
!! ---------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('stp_2D')
!
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'stp_2D : barotropic field in single first '
IF(lwp) WRITE(numout,*) '~~~~~~'
ENDIF
!
IF( ln_linssh ) THEN !== Compute ww(:,:,1) ==! (needed for momentum advection)
!!gm only in Flux Form, Vector Form dzU_z=0 assumed to be zero
!!gm ww(k=1) = div_h(uu_b) ==> modif dans dynadv <<<=== TO BE DONE
ENDIF
ALLOCATE( sshe_rhs(jpi,jpj) , Ue_rhs(jpi,jpj) , Ve_rhs(jpi,jpj) , CdU_u(jpi,jpj) , CdU_v(jpi,jpj) )
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! RHS of barotropic momentum Eq.
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! !======================================!
! !== Dynamics 2D RHS from 3D trends ==! (HADV + LDF + HPG) (No Coriolis trend)
! !======================================!
uu(:,:,:,Krhs) = 0._wp ! set dynamics trends to zero
vv(:,:,:,Krhs) = 0._wp
SELECT CASE( n_dynadv ) !* compute horizontal advection only *!
CASE( np_VEC_c2 ) !- vector form ( HADV = KEG + VOR )
CALL dyn_keg( kt, nn_dynkeg, Kbb, uu, vv, Krhs ) ! grad_h(KE)
CALL dyn_vor( kt, Kbb, uu, vv, Krhs, np_RVO ) ! relative vorticity
CASE( np_FLX_c2 ) !- flux form ( HADV = ADV_h + MET )
CALL dyn_adv_cen2( kt , Kbb, uu, vv, Krhs, no_zad = 1 ) ! 2nd order centered scheme
CALL dyn_vor ( kt , Kbb, uu, vv, Krhs, np_MET ) ! metric term
CASE( np_FLX_ubs )
CALL dyn_adv_ubs ( kt , Kbb, Kbb, uu, vv, Krhs, no_zad = 1) ! 3rd order UBS scheme (UP3)
CALL dyn_vor ( kt , Kbb, uu, vv, Krhs, np_MET ) ! metric term
END SELECT
!
! !* lateral viscosity *!
CALL dyn_ldf( kt, Kbb, Kbb, uu, vv, Krhs )
#if defined key_agrif
IF(.NOT. Agrif_Root() ) THEN !* AGRIF: sponge *!
CALL Agrif_Sponge_dyn
ENDIF
#endif
!
! !* hydrostatic pressure gradient *!
CALL eos ( ts , Kbb, rhd ) ! in situ density anomaly at Kbb
CALL dyn_hpg( kt , Kbb , uu, vv, Krhs ) ! horizontal gradient of Hydrostatic pressure
!
! !* vertical averaging *!
Ue_rhs(:,:) = SUM( e3u_0(:,:,:) * uu(:,:,:,Krhs) * umask(:,:,:), DIM=3 ) * r1_hu_0(:,:)
Ve_rhs(:,:) = SUM( e3v_0(:,:,:) * vv(:,:,:,Krhs) * vmask(:,:,:), DIM=3 ) * r1_hv_0(:,:)
!
!
! !===========================!
! !== external 2D forcing ==!
! !===========================!
!
! !* baroclinic drag forcing *! (also provide the barotropic drag coeff.)
!
CALL dyn_drg_init( Kbb, Kbb, uu, vv, uu_b, vv_b, Ue_rhs, Ve_rhs, CdU_u, CdU_v )
!
! !* wind forcing *!
!!st ATTENTION stoke drift !!
IF( ln_bt_fw ) THEN
DO_2D( 0, 0, 0, 0 )
Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kbb)
Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kbb)
END_2D
ELSE
zztmp = r1_rho0 * r1_2
DO_2D( 0, 0, 0, 0 )
Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kbb)
Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kbb)
END_2D
ENDIF
!
! !* atmospheric pressure forcing *!
IF( ln_apr_dyn ) THEN
IF( ln_bt_fw ) THEN ! FORWARD integration: use kt+1/2 pressure (NOW+1/2)
DO_2D( 0, 0, 0, 0 )
Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
END_2D
ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW)
zztmp = grav * r1_2
DO_2D( 0, 0, 0, 0 )
Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &
& + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj)
Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + zztmp * ( 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
ENDIF
!
! !* snow+ice load *! (embedded sea ice)
IF( ln_ice_embd ) THEN
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 )
Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
END_2D
DEALLOCATE( zpice )
ENDIF
!
! !* surface wave load *! (Bernoulli head)
!
IF( ln_wave .AND. ln_bern_srfc ) THEN
DO_2D( 0, 0, 0, 0 )
Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) * r1_e1u(ji,jj) !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3]
Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) * r1_e1u(ji,jj)
END_2D
ENDIF
!
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! RHS of see surface height Eq.
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!
! !== Net water flux forcing ==! (applied to a water column)
!
IF (ln_bt_fw) THEN ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2)
sshe_rhs(:,:) = r1_rho0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) )
ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW)
zztmp = r1_rho0 * r1_2
sshe_rhs(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) &
& - rnf(:,:) - rnf_b(:,:) &
& + fwfisf_cav(:,:) + fwfisf_cav_b(:,:) &
& + fwfisf_par(:,:) + fwfisf_par_b(:,:) )
ENDIF
!
! !== Stokes drift divergence ==! (if exist)
!
IF( ln_sdw ) sshe_rhs(:,:) = sshe_rhs(:,:) + div_sd(:,:)
!
!
! !== ice sheet coupling ==!
!
IF( ln_isf .AND. ln_isfcpl ) THEN
IF( ln_rstart .AND. kt == nit000 ) sshe_rhs(:,:) = sshe_rhs(:,:) + risfcpl_ssh(:,:)
IF( ln_isfcpl_cons ) sshe_rhs(:,:) = sshe_rhs(:,:) + risfcpl_cons_ssh(:,:)
ENDIF
!
#if defined key_asminc
! !== Add the IAU weighted SSH increment ==!
!
IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) sshe_rhs(:,:) = sshe_rhs(:,:) - ssh_iau(:,:)
#endif
!
#if defined key_agrif
! !== AGRIF : fill boundary data arrays (on both )
IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt )
#endif
!
!
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! Compute ssh and (uu_b,vv_b) at N+1 (Kaa)
!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
! using a split-explicit time integration in forward mode
! ( ABM3-AM4 time-integration Shchepetkin et al. OM2005) with temporal diffusion (Demange et al. JCP2019) )
CALL dyn_spg_ts( kt, Kbb, Kbb, Krhs, uu, vv, ssh, uu_b, vv_b, Kaa ) ! time-splitting
DEALLOCATE( sshe_rhs , Ue_rhs , Ve_rhs , CdU_u , CdU_v )
!!gm this is useless I guess : RK3, done in each stages
!
! IF( ln_dynspg_ts ) THEN ! With split-explicit free surface, since now transports have been updated and ssh(:,:,Krhs)
! ! as well as vertical scale factors and vertical velocity need to be updated
! CALL div_hor ( kstp, Kbb, Kmm ) ! Horizontal divergence (2nd call in time-split case)
! IF(.NOT.lk_linssh) CALL dom_qco_r3c( ssh(:,:,Kaa), r3t(:,:,Kaa), r3u(:,:,Kaa), r3v(:,:,Kaa), r3f(:,:) ) ! update ssh/h_0 ratio at t,u,v,f pts
! ENDIF
!!gm
!
IF( ln_timing ) CALL timing_stop('stp_2D')
!
END SUBROUTINE stp_2D
#else
!!----------------------------------------------------------------------
!! default option EMPTY MODULE qco not activated
!!----------------------------------------------------------------------
#endif
!!======================================================================
END MODULE stp2d