Changeset 12377 for NEMO/trunk/src/OCE/DYN/dynspg_ts.F90
- Timestamp:
- 2020-02-12T15:39:06+01:00 (4 years ago)
- Location:
- NEMO/trunk
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/DYN/dynspg_ts.F90 (modified) (44 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEAD ext/AGRIF5 ^/vendors/AGRIF/dev_r11615_ENHANCE-04_namelists_as_internalfiles_agrif@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL
-
- Property svn:externals
-
NEMO/trunk/src/OCE/DYN/dynspg_ts.F90
r12206 r12377 1 1 MODULE dynspg_ts 2 2 3 !! Includes ROMS wd scheme with diagnostic outputs ; un and uaupdates are commented out !3 !! Includes ROMS wd scheme with diagnostic outputs ; puu(:,:,:,Kmm) and puu(:,:,:,Krhs) updates are commented out ! 4 4 5 5 !!====================================================================== … … 31 31 USE dom_oce ! ocean space and time domain 32 32 USE sbc_oce ! surface boundary condition: ocean 33 USE isf_oce ! ice shelf variable (fwfisf) 33 34 USE zdf_oce ! vertical physics: variables 34 35 USE zdfdrg ! vertical physics: top/bottom drag coef. 35 USE sbcisf ! ice shelf variable (fwfisf)36 36 USE sbcapr ! surface boundary condition: atmospheric pressure 37 37 USE dynadv , ONLY: ln_dynadv_vec … … 44 44 USE bdytides ! open boundary condition data 45 45 USE bdydyn2d ! open boundary conditions on barotropic variables 46 USE sbctide ! tides 47 USE updtide ! tide potential 46 USE tide_mod ! 48 47 USE sbcwave ! surface wave 49 48 #if defined key_agrif … … 87 86 88 87 !! * Substitutions 89 # include " vectopt_loop_substitute.h90"88 # include "do_loop_substitute.h90" 90 89 !!---------------------------------------------------------------------- 91 90 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 117 116 118 117 119 SUBROUTINE dyn_spg_ts( kt )118 SUBROUTINE dyn_spg_ts( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa ) 120 119 !!---------------------------------------------------------------------- 121 120 !! … … 132 131 !! 133 132 !! ** Action : 134 !! -Update the filtered free surface at step "n+1" : ssha135 !! -Update filtered barotropic velocities at step "n+1" : ua_b, va_b133 !! -Update the filtered free surface at step "n+1" : pssh(:,:,Kaa) 134 !! -Update filtered barotropic velocities at step "n+1" : puu_b(:,:,:,Kaa), vv_b(:,:,:,Kaa) 136 135 !! -Compute barotropic advective fluxes at step "n" : un_adv, vn_adv 137 136 !! These are used to advect tracers and are compliant with discrete 138 137 !! continuity equation taken at the baroclinic time steps. This 139 138 !! ensures tracers conservation. 140 !! - ( ua, va) momentum trend updated with barotropic component.139 !! - (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) momentum trend updated with barotropic component. 141 140 !! 142 141 !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005. 143 142 !!--------------------------------------------------------------------- 144 INTEGER, INTENT(in) :: kt ! ocean time-step index 143 INTEGER , INTENT( in ) :: kt ! ocean time-step index 144 INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices 145 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 146 REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(inout) :: pssh, puu_b, pvv_b ! SSH and barotropic velocities at main time levels 145 147 ! 146 148 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 168 170 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points 169 171 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2 ! averages over the sub-steps of zuwdmask and zvwdmask 172 REAL(wp) :: zt0substep ! Time of day at the beginning of the time substep 170 173 !!---------------------------------------------------------------------- 171 174 ! … … 223 226 ! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends) 224 227 ! ! --------------------------- ! 225 zu_frc(:,:) = SUM( e3u _n(:,:,:) * ua(:,:,:) * umask(:,:,:) , DIM=3 ) * r1_hu_n(:,:)226 zv_frc(:,:) = SUM( e3v _n(:,:,:) * va(:,:,:) * vmask(:,:,:) , DIM=3 ) * r1_hv_n(:,:)227 ! 228 ! 229 ! != U a=> baroclinic trend =! (remove its vertical mean)230 DO jk = 1, jpkm1 ! ------------------------ !231 u a(:,:,jk) = ( ua(:,:,jk) - zu_frc(:,:) ) * umask(:,:,jk)232 v a(:,:,jk) = ( va(:,:,jk) - zv_frc(:,:) ) * vmask(:,:,jk)228 zu_frc(:,:) = SUM( e3u(:,:,:,Kmm) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm) 229 zv_frc(:,:) = SUM( e3v(:,:,:,Kmm) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm) 230 ! 231 ! 232 ! != U(Krhs) => baroclinic trend =! (remove its vertical mean) 233 DO jk = 1, jpkm1 ! ----------------------------- ! 234 uu(:,:,jk,Krhs) = ( uu(:,:,jk,Krhs) - zu_frc(:,:) ) * umask(:,:,jk) 235 vv(:,:,jk,Krhs) = ( vv(:,:,jk,Krhs) - zv_frc(:,:) ) * vmask(:,:,jk) 233 236 END DO 234 237 … … 239 242 ! ! ------------------------------------------------- ! 240 243 ! 241 IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2D_init ! Set zwz, the barotropic Coriolis force coefficient244 IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2D_init( Kmm ) ! Set zwz, the barotropic Coriolis force coefficient 242 245 ! ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes 243 246 ! 244 247 ! !* 2D Coriolis trends 245 zhU(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes246 zhV(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) ! NB: FULL domain : put a value in last row and column247 ! 248 CALL dyn_cor_2d( hu _n, hv_n, un_b, vn_b, zhU, zhV, & ! <<== in249 & zu_trd, zv_trd ) ! ==>> out248 zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:) ! now fluxes 249 zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) ! NB: FULL domain : put a value in last row and column 250 ! 251 CALL dyn_cor_2d( hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in 252 & zu_trd, zv_trd ) ! ==>> out 250 253 ! 251 254 IF( .NOT.ln_linssh ) THEN !* surface pressure gradient (variable volume only) 252 255 ! 253 256 IF( ln_wd_il ) THEN ! W/D : limiter applied to spgspg 254 CALL wad_spg( sshn, zcpx, zcpy ) ! Calculating W/D gravity filters, zcpx and zcpy 255 DO jj = 2, jpjm1 256 DO ji = 2, jpim1 ! SPG with the application of W/D gravity filters 257 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 258 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth 259 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 260 & * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth 261 END DO 262 END DO 257 CALL wad_spg( pssh(:,:,Kmm), zcpx, zcpy ) ! Calculating W/D gravity filters, zcpx and zcpy 258 DO_2D_00_00 259 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj ,Kmm) - pssh(ji ,jj ,Kmm) ) & 260 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth 261 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( pssh(ji ,jj+1,Kmm) - pssh(ji ,jj ,Kmm) ) & 262 & * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth 263 END_2D 263 264 ELSE ! now suface pressure gradient 264 DO jj = 2, jpjm1 265 DO ji = fs_2, fs_jpim1 ! vector opt. 266 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj) 267 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj) 268 END DO 269 END DO 270 ENDIF 271 ! 272 ENDIF 273 ! 274 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 275 DO ji = fs_2, fs_jpim1 276 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 277 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 278 END DO 279 END DO 265 DO_2D_00_00 266 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj ,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e1u(ji,jj) 267 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( pssh(ji ,jj+1,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e2v(ji,jj) 268 END_2D 269 ENDIF 270 ! 271 ENDIF 272 ! 273 DO_2D_00_00 274 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 275 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 276 END_2D 280 277 ! 281 278 ! != Add bottom stress contribution from baroclinic velocities =! 282 279 ! ! ----------------------------------------------------------- ! 283 CALL dyn_drg_init( zu_frc, zv_frc, zCdU_u, zCdU_v ) ! also provide the barotropic drag coefficients 284 ! 280 CALL dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, zu_frc, zv_frc, zCdU_u, zCdU_v ) ! also provide the barotropic drag coefficients 285 281 ! != Add atmospheric pressure forcing =! 286 282 ! ! ---------------------------------- ! 287 283 IF( ln_apr_dyn ) THEN 288 284 IF( ln_bt_fw ) THEN ! FORWARD integration: use kt+1/2 pressure (NOW+1/2) 289 DO jj = 2, jpjm1 290 DO ji = fs_2, fs_jpim1 ! vector opt. 291 zu_frc(ji,jj) = zu_frc(ji,jj) + grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 292 zv_frc(ji,jj) = zv_frc(ji,jj) + grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 293 END DO 294 END DO 285 DO_2D_00_00 286 zu_frc(ji,jj) = zu_frc(ji,jj) + grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 287 zv_frc(ji,jj) = zv_frc(ji,jj) + grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 288 END_2D 295 289 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW) 296 290 zztmp = grav * r1_2 297 DO jj = 2, jpjm1 298 DO ji = fs_2, fs_jpim1 ! vector opt. 299 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 300 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 301 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 302 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 303 END DO 304 END DO 305 ENDIF 291 DO_2D_00_00 292 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 293 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 294 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 295 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 296 END_2D 297 ENDIF 306 298 ENDIF 307 299 ! … … 309 301 ! ! ---------------------------------- ! 310 302 IF( ln_bt_fw ) THEN ! Add wind forcing 311 DO jj = 2, jpjm1 312 DO ji = fs_2, fs_jpim1 ! vector opt. 313 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu_n(ji,jj) 314 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv_n(ji,jj) 315 END DO 316 END DO 303 DO_2D_00_00 304 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu(ji,jj,Kmm) 305 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm) 306 END_2D 317 307 ELSE 318 308 zztmp = r1_rau0 * r1_2 319 DO jj = 2, jpjm1 320 DO ji = fs_2, fs_jpim1 ! vector opt. 321 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj) 322 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj) 323 END DO 324 END DO 309 DO_2D_00_00 310 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm) 311 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm) 312 END_2D 325 313 ENDIF 326 314 ! … … 331 319 ! ! --------------------------------------------------- ! 332 320 IF (ln_bt_fw) THEN ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2) 333 zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:))321 zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) 334 322 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW) 335 323 zztmp = r1_rau0 * r1_2 336 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) + fwfisf(:,:) + fwfisf_b(:,:) ) 324 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) & 325 & - rnf(:,:) - rnf_b(:,:) & 326 & + fwfisf_cav(:,:) + fwfisf_cav_b(:,:) & 327 & + fwfisf_par(:,:) + fwfisf_par_b(:,:) ) 337 328 ENDIF 338 329 ! != Add Stokes drift divergence =! (if exist) … … 340 331 zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:) 341 332 ENDIF 333 ! 334 ! ! ice sheet coupling 335 IF ( ln_isf .AND. ln_isfcpl ) THEN 336 ! 337 ! ice sheet coupling 338 IF( ln_rstart .AND. kt == nit000 ) THEN 339 zssh_frc(:,:) = zssh_frc(:,:) + risfcpl_ssh(:,:) 340 END IF 341 ! 342 ! conservation option 343 IF( ln_isfcpl_cons ) THEN 344 zssh_frc(:,:) = zssh_frc(:,:) + risfcpl_cons_ssh(:,:) 345 END IF 346 ! 347 END IF 342 348 ! 343 349 #if defined key_asminc … … 372 378 ! 373 379 IF( ln_linssh ) THEN ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0) 374 zhup2_e(:,:) = hu _n(:,:)375 zhvp2_e(:,:) = hv _n(:,:)376 zhtp2_e(:,:) = ht _n(:,:)380 zhup2_e(:,:) = hu(:,:,Kmm) 381 zhvp2_e(:,:) = hv(:,:,Kmm) 382 zhtp2_e(:,:) = ht(:,:) 377 383 ENDIF 378 384 ! 379 385 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 380 sshn_e(:,:) = sshn(:,:)381 un_e (:,:) = un_b(:,:)382 vn_e (:,:) = vn_b(:,:)383 ! 384 hu_e (:,:) = hu _n(:,:)385 hv_e (:,:) = hv _n(:,:)386 hur_e (:,:) = r1_hu _n(:,:)387 hvr_e (:,:) = r1_hv _n(:,:)386 sshn_e(:,:) = pssh(:,:,Kmm) 387 un_e (:,:) = puu_b(:,:,Kmm) 388 vn_e (:,:) = pvv_b(:,:,Kmm) 389 ! 390 hu_e (:,:) = hu(:,:,Kmm) 391 hv_e (:,:) = hv(:,:,Kmm) 392 hur_e (:,:) = r1_hu(:,:,Kmm) 393 hvr_e (:,:) = r1_hv(:,:,Kmm) 388 394 ELSE ! CENTRED integration: start from BEFORE fields 389 sshn_e(:,:) = sshb(:,:)390 un_e (:,:) = ub_b(:,:)391 vn_e (:,:) = vb_b(:,:)392 ! 393 hu_e (:,:) = hu _b(:,:)394 hv_e (:,:) = hv _b(:,:)395 hur_e (:,:) = r1_hu _b(:,:)396 hvr_e (:,:) = r1_hv _b(:,:)395 sshn_e(:,:) = pssh(:,:,Kbb) 396 un_e (:,:) = puu_b(:,:,Kbb) 397 vn_e (:,:) = pvv_b(:,:,Kbb) 398 ! 399 hu_e (:,:) = hu(:,:,Kbb) 400 hv_e (:,:) = hv(:,:,Kbb) 401 hur_e (:,:) = r1_hu(:,:,Kbb) 402 hvr_e (:,:) = r1_hv(:,:,Kbb) 397 403 ENDIF 398 404 ! 399 405 ! Initialize sums: 400 ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form)401 va_b (:,:) = 0._wp402 ssha (:,:) = 0._wp ! Sum for after averaged sea level406 puu_b (:,:,Kaa) = 0._wp ! After barotropic velocities (or transport if flux form) 407 pvv_b (:,:,Kaa) = 0._wp 408 pssh (:,:,Kaa) = 0._wp ! Sum for after averaged sea level 403 409 un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop 404 410 vn_adv(:,:) = 0._wp … … 419 425 ! !== Update the forcing ==! (BDY and tides) 420 426 ! 421 IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, kt_offset= noffset+1 ) 422 IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, kt_offset= noffset ) 427 IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, pt_offset= REAL(noffset+1,wp) ) 428 ! Update tide potential at the beginning of current time substep 429 IF( ln_tide_pot .AND. ln_tide ) THEN 430 zt0substep = REAL(nsec_day, wp) - 0.5_wp*rdt + (jn + noffset - 1) * rdt / REAL(nn_baro, wp) 431 CALL upd_tide(zt0substep, Kmm) 432 END IF 423 433 ! 424 434 ! !== extrapolation at mid-step ==! (jn+1/2) … … 457 467 ! 458 468 ! ! ocean u- and v-depth at mid-step (separate DO-loops remove the need of a lbc_lnk) 459 DO jj = 1, jpj 460 DO ji = 1, jpim1 ! not jpi-column 461 zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) & 462 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 463 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 464 END DO 465 END DO 466 DO jj = 1, jpjm1 ! not jpj-row 467 DO ji = 1, jpi 468 zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) & 469 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 470 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 471 END DO 472 END DO 469 DO_2D_11_10 470 zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) & 471 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 472 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 473 END_2D 474 DO_2D_10_11 475 zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) & 476 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 477 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 478 END_2D 473 479 ! 474 480 ENDIF … … 479 485 ! ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d 480 486 IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e ) 481 ! 487 ! 482 488 ! ! resulting flux at mid-step (not over the full domain) 483 489 zhU(1:jpim1,1:jpj ) = e2u(1:jpim1,1:jpj ) * ua_e(1:jpim1,1:jpj ) * zhup2_e(1:jpim1,1:jpj ) ! not jpi-column … … 486 492 #if defined key_agrif 487 493 ! Set fluxes during predictor step to ensure volume conservation 488 IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN 489 IF((nbondi == -1).OR.(nbondi == 2)) THEN 490 DO jj = 1, jpj 491 zhU(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj) 492 zhV(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj) 493 END DO 494 ENDIF 495 IF((nbondi == 1).OR.(nbondi == 2)) THEN 496 DO jj=1,jpj 497 zhU(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj) 498 zhV(nlci-nbghostcells :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells :nlci-1,jj) 499 END DO 500 ENDIF 501 IF((nbondj == -1).OR.(nbondj == 2)) THEN 502 DO ji=1,jpi 503 zhV(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1) 504 zhU(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1) 505 END DO 506 ENDIF 507 IF((nbondj == 1).OR.(nbondj == 2)) THEN 508 DO ji=1,jpi 509 zhV(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2) 510 zhU(ji,nlcj-nbghostcells :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells :nlcj-1) 511 END DO 512 ENDIF 513 ENDIF 494 IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) CALL agrif_dyn_ts_flux( jn, zhU, zhV ) 514 495 #endif 515 496 IF( ln_wd_il ) CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rdtbt) !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV … … 526 507 !-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --! 527 508 !-------------------------------------------------------------------------! 528 DO jj = 2, jpjm1 ! INNER domain 529 DO ji = 2, jpim1 530 zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj) 531 ssha_e(ji,jj) = ( sshn_e(ji,jj) - rdtbt * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj) 532 END DO 533 END DO 509 DO_2D_00_00 510 zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj) 511 ssha_e(ji,jj) = ( sshn_e(ji,jj) - rdtbt * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj) 512 END_2D 534 513 ! 535 514 CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp ) … … 553 532 ! Sea Surface Height at u-,v-points (vvl case only) 554 533 IF( .NOT.ln_linssh ) THEN 555 DO jj = 2, jpjm1 ! INNER domain, will be extended to whole domain later 556 DO ji = 2, jpim1 ! NO Vector Opt. 557 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 558 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 559 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) 560 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 561 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 562 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) 563 END DO 564 END DO 534 DO_2D_00_00 535 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 536 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 537 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) 538 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 539 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 540 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) 541 END_2D 565 542 ENDIF 566 543 ! … … 575 552 ! ! Surface pressure gradient 576 553 zldg = ( 1._wp - rn_scal_load ) * grav ! local factor 577 DO jj = 2, jpjm1 578 DO ji = 2, jpim1 579 zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 580 zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 581 END DO 582 END DO 554 DO_2D_00_00 555 zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 556 zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 557 END_2D 583 558 IF( ln_wd_il ) THEN ! W/D : gravity filters applied on pressure gradient 584 559 CALL wad_spg( zsshp2_e, zcpx, zcpy ) ! Calculating W/D gravity filters … … 595 570 ! Add tidal astronomical forcing if defined 596 571 IF ( ln_tide .AND. ln_tide_pot ) THEN 597 DO jj = 2, jpjm1 598 DO ji = fs_2, fs_jpim1 ! vector opt. 599 zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 600 zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 601 END DO 602 END DO 572 DO_2D_00_00 573 zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 574 zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 575 END_2D 603 576 ENDIF 604 577 ! … … 606 579 !jth do implicitly instead 607 580 IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs 608 DO jj = 2, jpjm1 609 DO ji = fs_2, fs_jpim1 ! vector opt. 610 zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj) 611 zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj) 612 END DO 613 END DO 581 DO_2D_00_00 582 zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj) 583 zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj) 584 END_2D 614 585 ENDIF 615 586 ! … … 626 597 !------------------------------------------------------------------------------------------------------------------------! 627 598 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 628 DO jj = 2, jpjm1 629 DO ji = fs_2, fs_jpim1 ! vector opt. 630 ua_e(ji,jj) = ( un_e(ji,jj) & 631 & + rdtbt * ( zu_spg(ji,jj) & 632 & + zu_trd(ji,jj) & 633 & + zu_frc(ji,jj) ) & 634 & ) * ssumask(ji,jj) 635 636 va_e(ji,jj) = ( vn_e(ji,jj) & 637 & + rdtbt * ( zv_spg(ji,jj) & 638 & + zv_trd(ji,jj) & 639 & + zv_frc(ji,jj) ) & 640 & ) * ssvmask(ji,jj) 641 END DO 642 END DO 599 DO_2D_00_00 600 ua_e(ji,jj) = ( un_e(ji,jj) & 601 & + rdtbt * ( zu_spg(ji,jj) & 602 & + zu_trd(ji,jj) & 603 & + zu_frc(ji,jj) ) & 604 & ) * ssumask(ji,jj) 605 606 va_e(ji,jj) = ( vn_e(ji,jj) & 607 & + rdtbt * ( zv_spg(ji,jj) & 608 & + zv_trd(ji,jj) & 609 & + zv_frc(ji,jj) ) & 610 & ) * ssvmask(ji,jj) 611 END_2D 643 612 ! 644 613 ELSE !* Flux form 645 DO jj = 2, jpjm1 646 DO ji = 2, jpim1 647 ! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2 648 ! ! backward interpolated depth used in spg terms at jn+1/2 649 zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 650 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 651 zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 652 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 653 ! ! inverse depth at jn+1 654 z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) 655 z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) ) 656 ! 657 ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) & 658 & + rdtbt * ( zhu_bck * zu_spg (ji,jj) & ! 659 & + zhup2_e(ji,jj) * zu_trd (ji,jj) & ! 660 & + hu_n (ji,jj) * zu_frc (ji,jj) ) ) * z1_hu 661 ! 662 va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) & 663 & + rdtbt * ( zhv_bck * zv_spg (ji,jj) & ! 664 & + zhvp2_e(ji,jj) * zv_trd (ji,jj) & ! 665 & + hv_n (ji,jj) * zv_frc (ji,jj) ) ) * z1_hv 666 END DO 667 END DO 614 DO_2D_00_00 615 ! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2 616 ! ! backward interpolated depth used in spg terms at jn+1/2 617 zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 618 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 619 zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 620 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 621 ! ! inverse depth at jn+1 622 z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) 623 z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) ) 624 ! 625 ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) & 626 & + rdtbt * ( zhu_bck * zu_spg (ji,jj) & ! 627 & + zhup2_e(ji,jj) * zu_trd (ji,jj) & ! 628 & + hu(ji,jj,Kmm) * zu_frc (ji,jj) ) ) * z1_hu 629 ! 630 va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) & 631 & + rdtbt * ( zhv_bck * zv_spg (ji,jj) & ! 632 & + zhvp2_e(ji,jj) * zv_trd (ji,jj) & ! 633 & + hv(ji,jj,Kmm) * zv_frc (ji,jj) ) ) * z1_hv 634 END_2D 668 635 ENDIF 669 636 !jth implicit bottom friction: 670 637 IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs 671 DO jj = 2, jpjm1 672 DO ji = fs_2, fs_jpim1 ! vector opt. 673 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj)) 674 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj)) 675 END DO 676 END DO 638 DO_2D_00_00 639 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj)) 640 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj)) 641 END_2D 677 642 ENDIF 678 643 … … 713 678 za1 = wgtbtp1(jn) 714 679 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities 715 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)716 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)680 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) 681 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) 717 682 ELSE ! Sum transports 718 683 IF ( .NOT.ln_wd_dl ) THEN 719 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)720 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)684 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) 685 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) 721 686 ELSE 722 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)723 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)687 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:) 688 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:) 724 689 END IF 725 690 ENDIF 726 691 ! ! Sum sea level 727 ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)692 pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:) 728 693 729 694 ! ! ==================== ! … … 737 702 IF (ln_bt_fw) THEN 738 703 IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN 739 DO jj = 1, jpj 740 DO ji = 1, jpi 741 zun_save = un_adv(ji,jj) 742 zvn_save = vn_adv(ji,jj) 743 ! ! apply the previously computed correction 744 un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) ) 745 vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) ) 746 ! ! Update corrective fluxes for next time step 747 un_bf(ji,jj) = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) ) 748 vn_bf(ji,jj) = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) ) 749 ! ! Save integrated transport for next computation 750 ub2_b(ji,jj) = zun_save 751 vb2_b(ji,jj) = zvn_save 752 END DO 753 END DO 704 DO_2D_11_11 705 zun_save = un_adv(ji,jj) 706 zvn_save = vn_adv(ji,jj) 707 ! ! apply the previously computed correction 708 un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) ) 709 vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) ) 710 ! ! Update corrective fluxes for next time step 711 un_bf(ji,jj) = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) ) 712 vn_bf(ji,jj) = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) ) 713 ! ! Save integrated transport for next computation 714 ub2_b(ji,jj) = zun_save 715 vb2_b(ji,jj) = zvn_save 716 END_2D 754 717 ELSE 755 718 un_bf(:,:) = 0._wp ! corrective fluxes for next time step set to zero … … 765 728 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 766 729 DO jk=1,jpkm1 767 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b768 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b730 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_2dt_b 731 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_2dt_b 769 732 END DO 770 733 ELSE 771 ! At this stage, ssha has been corrected: compute new depths at velocity points 772 DO jj = 1, jpjm1 773 DO ji = 1, jpim1 ! NO Vector Opt. 774 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 775 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 776 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 777 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 778 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 779 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 780 END DO 781 END DO 734 ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points 735 DO_2D_10_10 736 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 737 & * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) & 738 & + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) 739 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 740 & * ( e1e2t(ji,jj ) * pssh(ji,jj ,Kaa) & 741 & + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) 742 END_2D 782 743 CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 783 744 ! 784 745 DO jk=1,jpkm1 785 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b786 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b746 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_2dt_b 747 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_2dt_b 787 748 END DO 788 749 ! Save barotropic velocities not transport: 789 ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )790 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )750 puu_b(:,:,Kaa) = puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 751 pvv_b(:,:,Kaa) = pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 791 752 ENDIF 792 753 … … 794 755 ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 795 756 DO jk = 1, jpkm1 796 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)*r1_hu_n(:,:) - un_b(:,:) ) * umask(:,:,jk)797 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)*r1_hv_n(:,:) - vn_b(:,:) ) * vmask(:,:,jk)757 puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)*r1_hu(:,:,Kmm) - puu_b(:,:,Kmm) ) * umask(:,:,jk) 758 pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) + vn_adv(:,:)*r1_hv(:,:,Kmm) - pvv_b(:,:,Kmm) ) * vmask(:,:,jk) 798 759 END DO 799 760 800 761 IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN 801 ! need to set lbc here because not done prior time averaging802 CALL lbc_lnk_multi( 'dynspg_ts', zuwdav2, 'U', 1._wp, zvwdav2, 'V', 1._wp)803 762 DO jk = 1, jpkm1 804 un(:,:,jk) = ( un_adv(:,:)*r1_hu_n(:,:) &805 & + zuwdav2(:,:)*( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:)) ) * umask(:,:,jk)806 vn(:,:,jk) = ( vn_adv(:,:)*r1_hv_n(:,:) &807 & + zvwdav2(:,:)*( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:)) ) * vmask(:,:,jk)763 puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) & 764 & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk) 765 pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) & 766 & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk) 808 767 END DO 809 768 END IF 810 769 811 770 812 CALL iom_put( "ubar", un_adv(:,:)*r1_hu _n(:,:) ) ! barotropic i-current813 CALL iom_put( "vbar", vn_adv(:,:)*r1_hv _n(:,:) ) ! barotropic i-current771 CALL iom_put( "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) ) ! barotropic i-current 772 CALL iom_put( "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) ) ! barotropic i-current 814 773 ! 815 774 #if defined key_agrif … … 834 793 IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 ) 835 794 ! 836 CALL iom_put( "baro_u" , un_b) ! Barotropic U Velocity837 CALL iom_put( "baro_v" , vn_b) ! Barotropic V Velocity795 CALL iom_put( "baro_u" , puu_b(:,:,Kmm) ) ! Barotropic U Velocity 796 CALL iom_put( "baro_v" , pvv_b(:,:,Kmm) ) ! Barotropic V Velocity 838 797 ! 839 798 END SUBROUTINE dyn_spg_ts … … 1002 961 REAL(wp) :: zxr2, zyr2, zcmax ! local scalar 1003 962 REAL(wp), DIMENSION(jpi,jpj) :: zcu 1004 INTEGER :: inum1005 963 !!---------------------------------------------------------------------- 1006 964 ! 1007 965 ! Max courant number for ext. grav. waves 1008 966 ! 1009 DO jj = 1, jpj 1010 DO ji =1, jpi 1011 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1012 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1013 zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) ) 1014 END DO 1015 END DO 967 DO_2D_11_11 968 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 969 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 970 zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) ) 971 END_2D 1016 972 ! 1017 973 zcmax = MAXVAL( zcu(:,:) ) … … 1110 1066 1111 1067 1112 SUBROUTINE dyn_cor_2 d_init1068 SUBROUTINE dyn_cor_2D_init( Kmm ) 1113 1069 !!--------------------------------------------------------------------- 1114 !! *** ROUTINE dyn_cor_2 d_init ***1070 !! *** ROUTINE dyn_cor_2D_init *** 1115 1071 !! 1116 1072 !! ** Purpose : Set time splitting options … … 1124 1080 !! Compute zwz = f / ( height of the water colomn ) 1125 1081 !!---------------------------------------------------------------------- 1082 INTEGER, INTENT(in) :: Kmm ! Time index 1126 1083 INTEGER :: ji ,jj, jk ! dummy loop indices 1127 1084 REAL(wp) :: z1_ht … … 1133 1090 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point 1134 1091 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 1135 DO jj = 1, jpjm1 1136 DO ji = 1, jpim1 1137 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 1138 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp 1139 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1140 END DO 1141 END DO 1092 DO_2D_10_10 1093 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 1094 & ht(ji ,jj ) + ht(ji+1,jj ) ) * 0.25_wp 1095 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1096 END_2D 1142 1097 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 1143 DO jj = 1, jpjm1 1144 DO ji = 1, jpim1 1145 zwz(ji,jj) = ( ht_n (ji ,jj+1) + ht_n (ji+1,jj+1) & 1146 & + ht_n (ji ,jj ) + ht_n (ji+1,jj ) ) & 1147 & / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) & 1148 & + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) ) 1149 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1150 END DO 1151 END DO 1098 DO_2D_10_10 1099 zwz(ji,jj) = ( ht (ji ,jj+1) + ht (ji+1,jj+1) & 1100 & + ht (ji ,jj ) + ht (ji+1,jj ) ) & 1101 & / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) & 1102 & + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) ) 1103 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1104 END_2D 1152 1105 END SELECT 1153 1106 CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp ) 1154 1107 ! 1155 1108 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1156 DO jj = 2, jpj 1157 DO ji = 2, jpi 1158 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 1159 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 1160 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 1161 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 1162 END DO 1163 END DO 1109 DO_2D_01_01 1110 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 1111 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 1112 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 1113 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 1114 END_2D 1164 1115 ! 1165 1116 CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme) 1166 1117 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1167 DO jj = 2, jpj 1168 DO ji = 2, jpi 1169 z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) ) 1170 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht 1171 ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht 1172 ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht 1173 ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht 1174 END DO 1175 END DO 1118 DO_2D_01_01 1119 z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) 1120 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht 1121 ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht 1122 ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht 1123 ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht 1124 END_2D 1176 1125 ! 1177 1126 CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT ! … … 1200 1149 ! 1201 1150 !zhf(:,:) = hbatf(:,:) 1202 DO jj = 1, jpjm1 1203 DO ji = 1, jpim1 1204 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) & 1205 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) & 1206 & / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) & 1207 & + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp ) 1208 END DO 1209 END DO 1151 DO_2D_10_10 1152 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) & 1153 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) & 1154 & / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) & 1155 & + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp ) 1156 END_2D 1210 1157 ENDIF 1211 1158 ! … … 1216 1163 DO jk = 1, jpkm1 1217 1164 DO jj = 1, jpjm1 1218 zhf(:,jj) = zhf(:,jj) + e3f _n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)1165 zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 1219 1166 END DO 1220 1167 END DO 1221 1168 CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp ) 1222 1169 ! JC: TBC. hf should be greater than 0 1223 DO jj = 1, jpj 1224 DO ji = 1, jpi 1225 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) 1226 END DO 1227 END DO 1170 DO_2D_11_11 1171 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) 1172 END_2D 1228 1173 zwz(:,:) = ff_f(:,:) * zwz(:,:) 1229 1174 END SELECT … … 1233 1178 1234 1179 1235 SUBROUTINE dyn_cor_2d( hu_n, hv_n, un_b, vn_b, zhU, zhV, zu_trd, zv_trd )1180 SUBROUTINE dyn_cor_2d( phu, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd ) 1236 1181 !!--------------------------------------------------------------------- 1237 1182 !! *** ROUTINE dyn_cor_2d *** … … 1241 1186 INTEGER :: ji ,jj ! dummy loop indices 1242 1187 REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - - 1243 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: hu_n, hv_n, un_b, vn_b, zhU, zhV1188 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: phu, phv, punb, pvnb, zhU, zhV 1244 1189 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd 1245 1190 !!---------------------------------------------------------------------- 1246 1191 SELECT CASE( nvor_scheme ) 1247 1192 CASE( np_ENT ) ! enstrophy conserving scheme (f-point) 1248 DO jj = 2, jpjm1 1249 DO ji = 2, jpim1 1250 z1_hu = ssumask(ji,jj) / ( hu_n(ji,jj) + 1._wp - ssumask(ji,jj) ) 1251 z1_hv = ssvmask(ji,jj) / ( hv_n(ji,jj) + 1._wp - ssvmask(ji,jj) ) 1252 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu & 1253 & * ( e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) ) & 1254 & + e1e2t(ji ,jj)*ht_n(ji ,jj)*ff_t(ji ,jj) * ( vn_b(ji ,jj) + vn_b(ji ,jj-1) ) ) 1255 ! 1256 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1257 & * ( e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) ) & 1258 & + e1e2t(ji,jj )*ht_n(ji,jj )*ff_t(ji,jj ) * ( un_b(ji,jj ) + un_b(ji-1,jj ) ) ) 1259 END DO 1260 END DO 1193 DO_2D_00_00 1194 z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) ) 1195 z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) ) 1196 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu & 1197 & * ( e1e2t(ji+1,jj)*ht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) ) & 1198 & + e1e2t(ji ,jj)*ht(ji ,jj)*ff_t(ji ,jj) * ( pvnb(ji ,jj) + pvnb(ji ,jj-1) ) ) 1199 ! 1200 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1201 & * ( e1e2t(ji,jj+1)*ht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) & 1202 & + e1e2t(ji,jj )*ht(ji,jj )*ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) ) 1203 END_2D 1261 1204 ! 1262 1205 CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX 1263 DO jj = 2, jpjm1 1264 DO ji = fs_2, fs_jpim1 ! vector opt. 1265 zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj) 1266 zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) 1267 zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj) 1268 zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj) 1269 ! energy conserving formulation for planetary vorticity term 1270 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 1271 zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 1272 END DO 1273 END DO 1206 DO_2D_00_00 1207 zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj) 1208 zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) 1209 zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj) 1210 zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj) 1211 ! energy conserving formulation for planetary vorticity term 1212 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 1213 zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 1214 END_2D 1274 1215 ! 1275 1216 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 1276 DO jj = 2, jpjm1 1277 DO ji = fs_2, fs_jpim1 ! vector opt. 1278 zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) & 1279 & + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) 1280 zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) & 1281 & + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj) 1282 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 1283 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 1284 END DO 1285 END DO 1217 DO_2D_00_00 1218 zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) & 1219 & + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) 1220 zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) & 1221 & + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj) 1222 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 1223 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 1224 END_2D 1286 1225 ! 1287 1226 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1288 DO jj = 2, jpjm1 1289 DO ji = fs_2, fs_jpim1 ! vector opt. 1290 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) & 1291 & + ftnw(ji+1,jj) * zhV(ji+1,jj ) & 1292 & + ftse(ji,jj ) * zhV(ji ,jj-1) & 1293 & + ftsw(ji+1,jj) * zhV(ji+1,jj-1) ) 1294 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) & 1295 & + ftse(ji,jj+1) * zhU(ji ,jj+1) & 1296 & + ftnw(ji,jj ) * zhU(ji-1,jj ) & 1297 & + ftne(ji,jj ) * zhU(ji ,jj ) ) 1298 END DO 1299 END DO 1227 DO_2D_00_00 1228 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) & 1229 & + ftnw(ji+1,jj) * zhV(ji+1,jj ) & 1230 & + ftse(ji,jj ) * zhV(ji ,jj-1) & 1231 & + ftsw(ji+1,jj) * zhV(ji+1,jj-1) ) 1232 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) & 1233 & + ftse(ji,jj+1) * zhU(ji ,jj+1) & 1234 & + ftnw(ji,jj ) * zhU(ji-1,jj ) & 1235 & + ftne(ji,jj ) * zhU(ji ,jj ) ) 1236 END_2D 1300 1237 ! 1301 1238 END SELECT … … 1322 1259 ! 1323 1260 IF( ln_wd_dl_rmp ) THEN 1324 DO jj = 1, jpj 1325 DO ji = 1, jpi 1326 IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 1327 ! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 1328 ptmsk(ji,jj) = 1._wp 1329 ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 1330 ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1) ) 1331 ELSE 1332 ptmsk(ji,jj) = 0._wp 1333 ENDIF 1334 END DO 1335 END DO 1261 DO_2D_11_11 1262 IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 1263 ! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 1264 ptmsk(ji,jj) = 1._wp 1265 ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 1266 ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1) ) 1267 ELSE 1268 ptmsk(ji,jj) = 0._wp 1269 ENDIF 1270 END_2D 1336 1271 ELSE 1337 DO jj = 1, jpj 1338 DO ji = 1, jpi 1339 IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp 1340 ELSE ; ptmsk(ji,jj) = 0._wp 1341 ENDIF 1342 END DO 1343 END DO 1272 DO_2D_11_11 1273 IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp 1274 ELSE ; ptmsk(ji,jj) = 0._wp 1275 ENDIF 1276 END_2D 1344 1277 ENDIF 1345 1278 ! … … 1365 1298 !!---------------------------------------------------------------------- 1366 1299 ! 1367 DO jj = 1, jpj 1368 DO ji = 1, jpim1 ! not jpi-column 1369 IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj) 1370 ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj) 1371 ENDIF 1372 phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj) 1373 pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj) 1374 END DO 1375 END DO 1376 ! 1377 DO jj = 1, jpjm1 ! not jpj-row 1378 DO ji = 1, jpi 1379 IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj ) 1380 ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1) 1381 ENDIF 1382 phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 1383 pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj) 1384 END DO 1385 END DO 1300 DO_2D_11_10 1301 IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj) 1302 ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj) 1303 ENDIF 1304 phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj) 1305 pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj) 1306 END_2D 1307 ! 1308 DO_2D_10_11 1309 IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj ) 1310 ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1) 1311 ENDIF 1312 phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 1313 pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj) 1314 END_2D 1386 1315 ! 1387 1316 END SUBROUTINE wad_Umsk 1388 1317 1389 1318 1390 SUBROUTINE wad_spg( sshn, zcpx, zcpy )1319 SUBROUTINE wad_spg( pshn, zcpx, zcpy ) 1391 1320 !!--------------------------------------------------------------------- 1392 1321 !! *** ROUTINE wad_sp *** … … 1396 1325 INTEGER :: ji ,jj ! dummy loop indices 1397 1326 LOGICAL :: ll_tmp1, ll_tmp2 1398 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: sshn1327 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pshn 1399 1328 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy 1400 1329 !!---------------------------------------------------------------------- 1401 DO jj = 2, jpjm1 1402 DO ji = 2, jpim1 1403 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1404 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 1405 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 1406 & > rn_wdmin1 + rn_wdmin2 1407 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 1408 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1409 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 1410 IF(ll_tmp1) THEN 1411 zcpx(ji,jj) = 1.0_wp 1412 ELSEIF(ll_tmp2) THEN 1413 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 1414 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 1415 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 1416 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 1417 ELSE 1418 zcpx(ji,jj) = 0._wp 1419 ENDIF 1420 ! 1421 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1422 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 1423 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 1424 & > rn_wdmin1 + rn_wdmin2 1425 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( & 1426 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1427 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1428 1429 IF(ll_tmp1) THEN 1430 zcpy(ji,jj) = 1.0_wp 1431 ELSE IF(ll_tmp2) THEN 1432 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 1433 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 1434 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 1435 zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) ) 1436 ELSE 1437 zcpy(ji,jj) = 0._wp 1438 ENDIF 1439 END DO 1440 END DO 1330 DO_2D_00_00 1331 ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji+1,jj) ) > & 1332 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 1333 & MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji+1,jj) + ht_0(ji+1,jj) ) & 1334 & > rn_wdmin1 + rn_wdmin2 1335 ll_tmp2 = ( ABS( pshn(ji+1,jj) - pshn(ji ,jj)) > 1.E-12 ).AND.( & 1336 & MAX( pshn(ji,jj) , pshn(ji+1,jj) ) > & 1337 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 1338 IF(ll_tmp1) THEN 1339 zcpx(ji,jj) = 1.0_wp 1340 ELSEIF(ll_tmp2) THEN 1341 ! no worries about pshn(ji+1,jj) - pshn(ji ,jj) = 0, it won't happen ! here 1342 zcpx(ji,jj) = ABS( (pshn(ji+1,jj) + ht_0(ji+1,jj) - pshn(ji,jj) - ht_0(ji,jj)) & 1343 & / (pshn(ji+1,jj) - pshn(ji ,jj)) ) 1344 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 1345 ELSE 1346 zcpx(ji,jj) = 0._wp 1347 ENDIF 1348 ! 1349 ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji,jj+1) ) > & 1350 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 1351 & MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji,jj+1) + ht_0(ji,jj+1) ) & 1352 & > rn_wdmin1 + rn_wdmin2 1353 ll_tmp2 = ( ABS( pshn(ji,jj) - pshn(ji,jj+1)) > 1.E-12 ).AND.( & 1354 & MAX( pshn(ji,jj) , pshn(ji,jj+1) ) > & 1355 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1356 1357 IF(ll_tmp1) THEN 1358 zcpy(ji,jj) = 1.0_wp 1359 ELSE IF(ll_tmp2) THEN 1360 ! no worries about pshn(ji,jj+1) - pshn(ji,jj ) = 0, it won't happen ! here 1361 zcpy(ji,jj) = ABS( (pshn(ji,jj+1) + ht_0(ji,jj+1) - pshn(ji,jj) - ht_0(ji,jj)) & 1362 & / (pshn(ji,jj+1) - pshn(ji,jj )) ) 1363 zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) ) 1364 ELSE 1365 zcpy(ji,jj) = 0._wp 1366 ENDIF 1367 END_2D 1441 1368 1442 1369 END SUBROUTINE wad_spg … … 1444 1371 1445 1372 1446 SUBROUTINE dyn_drg_init( pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )1373 SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v ) 1447 1374 !!---------------------------------------------------------------------- 1448 1375 !! *** ROUTINE dyn_drg_init *** … … 1454 1381 !! ** Method : computation done over the INNER domain only 1455 1382 !!---------------------------------------------------------------------- 1456 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS 1457 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients 1383 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 1384 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in ) :: puu, pvv ! ocean velocities and RHS of momentum equation 1385 REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(in ) :: puu_b, pvv_b ! barotropic velocities at main time levels 1386 REAL(wp), DIMENSION(jpi,jpj) , INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS 1387 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients 1458 1388 ! 1459 1389 INTEGER :: ji, jj ! dummy loop indices … … 1467 1397 IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities) 1468 1398 1469 DO jj = 2, jpjm1 1470 DO ji = 2, jpim1 ! INNER domain 1471 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 1472 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) 1473 END DO 1474 END DO 1399 DO_2D_00_00 1400 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 1401 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) 1402 END_2D 1475 1403 ELSE ! bottom friction only 1476 DO jj = 2, jpjm1 1477 DO ji = 2, jpim1 ! INNER domain 1478 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 1479 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 1480 END DO 1481 END DO 1404 DO_2D_00_00 1405 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 1406 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 1407 END_2D 1482 1408 ENDIF 1483 1409 ! … … 1486 1412 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities 1487 1413 1488 DO jj = 2, jpjm1 1489 DO ji = 2, jpim1 ! INNER domain 1490 ikbu = mbku(ji,jj) 1491 ikbv = mbkv(ji,jj) 1492 zu_i(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) 1493 zv_i(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj) 1494 END DO 1495 END DO 1414 DO_2D_00_00 1415 ikbu = mbku(ji,jj) 1416 ikbv = mbkv(ji,jj) 1417 zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm) 1418 zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm) 1419 END_2D 1496 1420 ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities 1497 1421 1498 DO jj = 2, jpjm1 1499 DO ji = 2, jpim1 ! INNER domain 1500 ikbu = mbku(ji,jj) 1501 ikbv = mbkv(ji,jj) 1502 zu_i(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) 1503 zv_i(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj) 1504 END DO 1505 END DO 1422 DO_2D_00_00 1423 ikbu = mbku(ji,jj) 1424 ikbv = mbkv(ji,jj) 1425 zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb) 1426 zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb) 1427 END_2D 1506 1428 ENDIF 1507 1429 ! 1508 1430 IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please ! 1509 1431 zztmp = -1._wp / rdtbt 1510 DO jj = 2, jpjm1 1511 DO ji = 2, jpim1 ! INNER domain 1512 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & 1513 & r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) 1514 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( & 1515 & r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) 1516 END DO 1517 END DO 1432 DO_2D_00_00 1433 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & 1434 & r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) 1435 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( & 1436 & r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) 1437 END_2D 1518 1438 ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation) 1519 1439 1520 DO jj = 2, jpjm1 1521 DO ji = 2, jpim1 ! INNER domain 1522 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj) 1523 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj) 1524 END DO 1525 END DO 1440 DO_2D_00_00 1441 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj) 1442 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj) 1443 END_2D 1526 1444 END IF 1527 1445 ! … … 1532 1450 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity 1533 1451 1534 DO jj = 2, jpjm1 1535 DO ji = 2, jpim1 ! INNER domain 1536 iktu = miku(ji,jj) 1537 iktv = mikv(ji,jj) 1538 zu_i(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) 1539 zv_i(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 1540 END DO 1541 END DO 1452 DO_2D_00_00 1453 iktu = miku(ji,jj) 1454 iktv = mikv(ji,jj) 1455 zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm) 1456 zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm) 1457 END_2D 1542 1458 ELSE ! CENTRED integration: use BEFORE top baroclinic velocity 1543 1459 1544 DO jj = 2, jpjm1 1545 DO ji = 2, jpim1 ! INNER domain 1546 iktu = miku(ji,jj) 1547 iktv = mikv(ji,jj) 1548 zu_i(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) 1549 zv_i(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 1550 END DO 1551 END DO 1460 DO_2D_00_00 1461 iktu = miku(ji,jj) 1462 iktv = mikv(ji,jj) 1463 zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb) 1464 zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb) 1465 END_2D 1552 1466 ENDIF 1553 1467 ! 1554 1468 ! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation) 1555 1469 1556 DO jj = 2, jpjm1 1557 DO ji = 2, jpim1 ! INNER domain 1558 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu_n(ji,jj) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj) 1559 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv_n(ji,jj) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj) 1560 END DO 1561 END DO 1470 DO_2D_00_00 1471 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj) 1472 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj) 1473 END_2D 1562 1474 ! 1563 1475 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
