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
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
!
! !* compute advection *!
!
CALL dyn_adv( kstp, Kbb, Kbb , uu, vv, Krhs) !- vector form KEG+ZAD
! !- flux form ADV
SELECT CASE( n_dynadv )
CASE( np_VEC_c2 ) !- vector form (add relative vorticity term)
CALL dyn_vor( kt, Kbb, uu, vv, Krhs, np_RVO )
CASE( np_FLX_c2 , np_FLX_ubs ) !- flux form (add metric term)
CALL dyn_vor ( kt , Kbb, uu, vv, Krhs, np_MET )
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(:,:)
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
! !===========================!
! !== 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)
Sibylle TECHENE
committed
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(:,:) &
Sibylle TECHENE
committed
& - 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
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
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! 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