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