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.
stp2d.F90 in NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE – NEMO

source: NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/stp2d.F90 @ 14949

Last change on this file since 14949 was 14949, checked in by acc, 3 years ago

#2605 : bug fixes to allow all SETTE tests to compile and run on this branch without key_RK3. This is a missed variable name change in stp2d.F90 and an uninecessarily duplicated routine in stprk3.F90. The duplicated routine (finalize_lbc) is unused but AGRIF compilations will pick up the incompatible first version from stpmlf.F90 (AGRIF does not respect scope). The copy in stprk3.F90 has been removed.
File size: 13.5 KB
Line 
1MODULE stp2d
2   !!======================================================================
3   !!                       ***  MODULE stp2d  ***
4   !! Time-stepping   : manager of the ocean, tracer and ice time stepping
5   !!                   using a 3rd order Rung Kuta  with fixed or quasi-eulerian coordinate
6   !!======================================================================
7   !! History :  4.5  !  2021-01  (S. Techene, G. Madec, N. Ducousso, F. Lemarie)  Original code
8   !!   NEMO     
9   !!----------------------------------------------------------------------
10#if defined key_qco   ||   defined key_linssh
11   !!----------------------------------------------------------------------
12   !!   'key_qco'                        Quasi-Eulerian vertical coordinate
13   !!                          OR
14   !!   'key_linssh                       Fixed in time vertical coordinate
15   !!----------------------------------------------------------------------
16
17   !!----------------------------------------------------------------------
18   !!    stp_2D   : RK3 case
19   !!----------------------------------------------------------------------
20   USE step_oce       ! time stepping used modules
21   USE domqco         ! quasi-eulerian coordinate      (dom_qco_r3c routine)
22   USE dynadv_cen2    ! centred flux form advection    (dyn_adv_cen2 routine)
23   USE dynadv_ubs     ! UBS flux form advection        (dyn_adv_ubs  routine)
24   USE dynkeg         ! kinetic energy gradient        (dyn_keg      routine)
25   USE dynspg_ts      ! 2D mode integration
26   USE sbc_ice , ONLY : snwice_mass, snwice_mass_b
27   USE sbcapr         ! surface boundary condition: atmospheric pressure
28   USE sbcwave,  ONLY : bhd_wave
29#if defined key_agrif
30   USE agrif_oce_interp
31   USE agrif_oce_sponge
32#endif
33
34   PRIVATE
35
36   PUBLIC   stp_2D   ! called by nemogcm.F90
37   REAL (wp) :: r1_2 = 0.5_wp 
38
39   !! * Substitutions
40#  include "do_loop_substitute.h90"
41#  include "domzgr_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
44   !! $Id: step.F90 12377 2020-02-12 14:39:06Z acc $
45   !! Software governed by the CeCILL license (see ./LICENSE)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE stp_2D( kt, Kbb, Kmm, Kaa, Krhs )
50      !!----------------------------------------------------------------------
51      !!                     ***  ROUTINE stp_2D  ***
52      !!
53      !! ** Purpose : - Compute sea-surface height and barotropic velocity at Kaa
54      !!                in single 1st RK3.
55      !!
56      !! ** Method  : -1- Compute the 3D to 2D forcing
57      !!                 * Momentum (Ue,Ve)_rhs :
58      !!                      3D to 2D dynamics, i.e. the vertical sum of :
59      !!                        - Hor. adv. : KEG   + RVO in vector form
60      !!                                    : ADV_h + MET in flux   form
61      !!                        - LDF Lateral mixing
62      !!                        - HPG Hor. pressure gradient
63      !!                      External forcings
64      !!                        - baroclinic drag
65      !!                        - wind
66      !!                        - atmospheric pressure
67      !!                        - snow+ice load
68      !!                        - surface wave load
69      !!                 * ssh (sshe_rhs) :
70      !!                      Net column average freshwater flux
71      !!
72      !!              -2- Solve the external mode Eqs. using sub-time step
73      !!                  by a call to dyn_spg_ts (will be renamed dyn_2D or stp_2D)
74      !!
75      !! ** action  :   ssh            : N+1 sea surface height (Kaa=N+1)
76      !!                (uu_b,vv_b)    : N+1 barotropic velocity
77      !!                (un_adv,vn_adv): barotropic transport from N to N+1
78      !!----------------------------------------------------------------------
79      INTEGER, INTENT(in) ::   kt, Kbb, Kmm, Kaa, Krhs    ! ocean time-step and time-level indices
80      !
81      INTEGER  ::   ji, jj, jk   ! dummy loop indices
82      REAL(wp) ::   zg_2, zintp, zgrho0r, zld, zztmp     ! local scalars
83      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zpice   ! 2D workspace
84      !! ---------------------------------------------------------------------
85      !
86      IF( ln_timing )   CALL timing_start('stp_2D')
87      !
88      IF( kt == nit000 ) THEN
89         IF(lwp) WRITE(numout,*)
90         IF(lwp) WRITE(numout,*) 'stp_2D : barotropic field in single first '
91         IF(lwp) WRITE(numout,*) '~~~~~~'
92      ENDIF
93      !
94      IF( ln_linssh ) THEN    !==  Compute ww(:,:,1)  ==!   (needed for momentum advection)
95!!gm  only in Flux Form, Vector Form  dzU_z=0 assumed to be zero
96!!gm  ww(k=1) = div_h(uu_b) ==> modif dans dynadv                        <<<=== TO BE DONE
97      ENDIF
98
99      ALLOCATE( sshe_rhs(jpi,jpj) , Ue_rhs(jpi,jpj) , Ve_rhs(jpi,jpj) , CdU_u(jpi,jpj) , CdU_v(jpi,jpj) )
100
101      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
102      !   RHS of barotropic momentum  Eq.
103      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
104
105      !                       !======================================!
106      !                       !==  Dynamics 2D RHS from 3D trends  ==!   (HADV + LDF + HPG) (No Coriolis trend)
107      !                       !======================================!
108
109      uu(:,:,:,Krhs) = 0._wp        ! set dynamics trends to zero
110      vv(:,:,:,Krhs) = 0._wp
111
112      SELECT CASE( n_dynadv )       !*  compute horizontal advection  only  *!
113      CASE( np_VEC_c2  )                                    !- vector form  ( HADV = KEG + VOR )
114         CALL dyn_keg( kt, nn_dynkeg, Kbb, uu, vv, Krhs )                  ! grad_h(KE)
115         CALL dyn_vor( kt,            Kbb, uu, vv, Krhs, np_RVO )          ! relative vorticity
116      CASE( np_FLX_c2  )                                    !-  flux form   ( HADV = ADV_h + MET )
117         CALL dyn_adv_cen2( kt      , Kbb, uu, vv, Krhs, no_zad = 1 )      ! 2nd order centered scheme
118         CALL dyn_vor     ( kt      , Kbb, uu, vv, Krhs, np_MET )          ! metric term
119      CASE( np_FLX_ubs ) 
120         CALL dyn_adv_ubs ( kt      , Kbb, Kbb, uu, vv, Krhs, no_zad = 1)  ! 3rd order UBS scheme (UP3)
121         CALL dyn_vor     ( kt      , Kbb, uu, vv, Krhs, np_MET )          ! metric term
122      END SELECT
123      !
124      !                             !*  lateral viscosity  *!
125      CALL dyn_ldf( kt,   Kbb, Kbb, uu, vv, Krhs )
126#if defined key_agrif
127      IF(.NOT. Agrif_Root() ) THEN  !*  AGRIF: sponge *!
128         CALL Agrif_Sponge_dyn
129      ENDIF
130#endif
131      !
132      !                             !*  hydrostatic pressure gradient  *! 
133      CALL eos    ( ts , Kbb, rhd )                          ! in situ density anomaly at Kbb
134      CALL dyn_hpg( kt  , Kbb     , uu, vv, Krhs )           ! horizontal gradient of Hydrostatic pressure
135      !
136      !                             !*  vertical averaging  *!
137      Ue_rhs(:,:) = SUM( e3u_0(:,:,:) * uu(:,:,:,Krhs) * umask(:,:,:), DIM=3 ) * r1_hu_0(:,:)
138      Ve_rhs(:,:) = SUM( e3v_0(:,:,:) * vv(:,:,:,Krhs) * vmask(:,:,:), DIM=3 ) * r1_hv_0(:,:)
139      !
140      !
141      !                       !===========================!
142      !                       !==  external 2D forcing  ==!
143      !                       !===========================!
144      !
145      !            !* baroclinic drag forcing *!   (also provide the barotropic drag coeff.)
146      !
147      CALL dyn_drg_init( Kbb, Kbb, uu, vv, uu_b, vv_b, Ue_rhs, Ve_rhs, CdU_u, CdU_v )
148      !
149      !                             !* wind forcing *!
150!!st ATTENTION stoke drift !!
151      IF( ln_bt_fw ) THEN
152         DO_2D( 0, 0, 0, 0 )
153            Ue_rhs(ji,jj) =  Ue_rhs(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kbb)
154            Ve_rhs(ji,jj) =  Ve_rhs(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kbb)
155         END_2D
156      ELSE
157         zztmp = r1_rho0 * r1_2
158         DO_2D( 0, 0, 0, 0 )
159            Ue_rhs(ji,jj) =  Ue_rhs(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kbb)
160            Ve_rhs(ji,jj) =  Ve_rhs(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kbb)
161         END_2D
162      ENDIF
163      !
164      !                             !* atmospheric pressure forcing *!
165      IF( ln_apr_dyn ) THEN
166         IF( ln_bt_fw ) THEN                          ! FORWARD integration: use kt+1/2 pressure (NOW+1/2)
167            DO_2D( 0, 0, 0, 0 )
168               Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + grav * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)
169               Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + grav * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)
170            END_2D
171         ELSE                                         ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW)
172            zztmp = grav * r1_2
173            DO_2D( 0, 0, 0, 0 )
174               Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + zztmp * (  ssh_ib (ji+1,jj  ) - ssh_ib (ji,jj)  &
175                    &                                   + ssh_ibb(ji+1,jj  ) - ssh_ibb(ji,jj)  ) * r1_e1u(ji,jj)
176               Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + zztmp * (  ssh_ib (ji  ,jj+1) - ssh_ib (ji,jj)  &
177                    &                                   + ssh_ibb(ji  ,jj+1) - ssh_ibb(ji,jj)  ) * r1_e2v(ji,jj)
178            END_2D
179         ENDIF
180      ENDIF
181      !
182      !                             !* snow+ice load *!   (embedded sea ice)
183      IF( ln_ice_embd ) THEN
184         ALLOCATE( zpice(jpi,jpj) )
185         zintp = REAL( MOD( kt-1, nn_fsbc ) ) / REAL( nn_fsbc )
186         zgrho0r     = - grav * r1_rho0
187         zpice(:,:) = (  zintp * snwice_mass(:,:) + ( 1.- zintp ) * snwice_mass_b(:,:)  ) * zgrho0r
188         DO_2D( 0, 0, 0, 0 )
189            Ue_rhs(ji,jj) = Ue_rhs(ji,jj) + ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
190            Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
191         END_2D
192         DEALLOCATE( zpice )
193      ENDIF
194      !
195      !                             !* surface wave load *!   (Bernoulli head)
196      !
197      IF( ln_wave .AND. ln_bern_srfc ) THEN
198         DO_2D( 0, 0, 0, 0 )
199            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]
200            Ve_rhs(ji,jj) = Ve_rhs(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) * r1_e1u(ji,jj)
201         END_2D
202      ENDIF
203      !
204      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
205      !   RHS of see surface height  Eq.
206      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
207      !
208      !     !==  Net water flux forcing  ==!  (applied to a water column)
209      !
210      IF (ln_bt_fw) THEN                          ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2)
211         sshe_rhs(:,:) = r1_rho0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) )
212      ELSE                                        ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW)
213         zztmp = r1_rho0 * r1_2
214         sshe_rhs(:,:) = zztmp * (   emp(:,:)        + emp_b(:,:)          &
215            &                      - rnf(:,:)        - rnf_b(:,:)          &
216            &                      + fwfisf_cav(:,:) + fwfisf_cav_b(:,:)   &
217            &                      + fwfisf_par(:,:) + fwfisf_par_b(:,:)   )
218      ENDIF
219      !
220      !     !==  Stokes drift divergence  ==!   (if exist)
221      !
222      IF( ln_sdw )    sshe_rhs(:,:) = sshe_rhs(:,:) + div_sd(:,:)
223      !
224      !
225      !     !==  ice sheet coupling  ==!
226      !
227      IF( ln_isf .AND. ln_isfcpl ) THEN
228         IF( ln_rstart .AND. kt == nit000 )   sshe_rhs(:,:) = sshe_rhs(:,:) + risfcpl_ssh(:,:)
229         IF( ln_isfcpl_cons               )   sshe_rhs(:,:) = sshe_rhs(:,:) + risfcpl_cons_ssh(:,:)
230      ENDIF
231      !
232#if defined key_asminc
233      !                 !==  Add the IAU weighted SSH increment  ==!
234      !
235      IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau )   sshe_rhs(:,:) = sshe_rhs(:,:) - ssh_iau(:,:)
236#endif
237      !
238#if defined key_agrif
239      !                 !==  AGRIF : fill boundary data arrays (on both )
240         IF( .NOT.Agrif_Root() )   CALL agrif_dta_ts( kt )
241#endif
242      !
243      !
244      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
245      !             Compute ssh and (uu_b,vv_b)  at N+1  (Kaa)
246      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
247
248      !    using a split-explicit time integration in forward mode
249      !    ( ABM3-AM4 time-integration Shchepetkin et al. OM2005) with temporal diffusion (Demange et al. JCP2019) )
250
251      CALL dyn_spg_ts( kt, Kbb, Kbb, Krhs, uu, vv, ssh, uu_b, vv_b, Kaa ) ! time-splitting
252                   
253
254      DEALLOCATE( sshe_rhs , Ue_rhs , Ve_rhs , CdU_u , CdU_v )
255
256!!gm  this is useless I guess : RK3,  done in each stages
257!
258!      IF( ln_dynspg_ts ) THEN      ! With split-explicit free surface, since now transports have been updated and ssh(:,:,Krhs)
259!                                   ! as well as vertical scale factors and vertical velocity need to be updated
260!                            CALL div_hor    ( kstp, Kbb, Kmm )                ! Horizontal divergence  (2nd call in time-split case)
261!         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
262!      ENDIF
263!!gm   
264      !
265      IF( ln_timing )   CALL timing_stop('stp_2D')
266      !
267   END SUBROUTINE stp_2D
268
269
270#else
271   !!----------------------------------------------------------------------
272   !!   default option             EMPTY MODULE           qco not activated
273   !!----------------------------------------------------------------------
274#endif
275   
276   !!======================================================================
277END MODULE stp2d
Note: See TracBrowser for help on using the repository browser.