[14547] | 1 | MODULE 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 |
---|
[14782] | 29 | #if defined key_agrif |
---|
| 30 | USE agrif_oce_interp |
---|
[14800] | 31 | USE agrif_oce_sponge |
---|
[14782] | 32 | #endif |
---|
[14547] | 33 | |
---|
| 34 | PRIVATE |
---|
| 35 | |
---|
| 36 | PUBLIC stp_2D ! called by nemogcm.F90 |
---|
[14782] | 37 | REAL (wp) :: r1_2 = 0.5_wp |
---|
[14547] | 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 | !!---------------------------------------------------------------------- |
---|
| 47 | CONTAINS |
---|
| 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 ) |
---|
[14800] | 126 | #if defined key_agrif |
---|
| 127 | IF(.NOT. Agrif_Root() ) THEN !* AGRIF: sponge *! |
---|
| 128 | CALL Agrif_Sponge_dyn |
---|
| 129 | ENDIF |
---|
| 130 | #endif |
---|
[14547] | 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 | ! |
---|
[14949] | 235 | IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) sshe_rhs(:,:) = sshe_rhs(:,:) - ssh_iau(:,:) |
---|
[14547] | 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 | !!====================================================================== |
---|
| 277 | END MODULE stp2d |
---|