- Timestamp:
- 2019-11-26T15:11:43+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/ENHANCE-02_ISF_nemo_TEST_MERGE/src/OCE/DYN/dynspg_ts.F90
r11541 r11967 65 65 USE diatmb ! Top,middle,bottom output 66 66 67 USE iom ! to remove 68 67 69 IMPLICIT NONE 68 70 PRIVATE … … 105 107 ! 106 108 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) ) 107 !108 109 IF( ln_dynvor_een .OR. ln_dynvor_eeT ) & 109 & ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 110 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) 110 & ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , ftsw(jpi,jpj) , ftse(jpi,jpj), STAT=ierr(2) ) 111 111 ! 112 112 ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) ) … … 150 150 LOGICAL :: ll_fw_start ! =T : forward integration 151 151 LOGICAL :: ll_init ! =T : special startup of 2d equations 152 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D 153 INTEGER :: ikbu, iktu, noffset ! local integers 154 INTEGER :: ikbv, iktv ! - - 155 REAL(wp) :: r1_2dt_b, z2dt_bf ! local scalars 156 REAL(wp) :: zx1, zx2, zu_spg, zhura, z1_hu ! - - 157 REAL(wp) :: zy1, zy2, zv_spg, zhvra, z1_hv ! - - 152 INTEGER :: noffset ! local integers : time offset for bdy update 153 REAL(wp) :: r1_2dt_b, z1_hu, z1_hv ! local scalars 158 154 REAL(wp) :: za0, za1, za2, za3 ! - - 159 REAL(wp) :: zmdi, zztmp , z1_ht ! - - 160 REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf 161 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc 162 REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv 163 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e, zhtp2_e 164 REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e 155 REAL(wp) :: zmdi, zztmp, zldg ! - - 156 REAL(wp) :: zhu_bck, zhv_bck, zhdiv ! - - 157 REAL(wp) :: zun_save, zvn_save ! - - 158 REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zu_frc, zu_spg, zssh_frc 159 REAL(wp), DIMENSION(jpi,jpj) :: zv_trd, zv_frc, zv_spg 160 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhtp2_e 161 REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zsshp2_e 165 162 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points 163 REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV ! fluxes 166 164 ! 167 165 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. … … 183 181 zwdramp = r_rn_wdmin1 ! simplest ramp 184 182 ! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 185 ! ! reciprocal of baroclinic time step 186 IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt 187 ELSE ; z2dt_bf = 2.0_wp * rdt 188 ENDIF 189 r1_2dt_b = 1.0_wp / z2dt_bf 183 ! ! inverse of baroclinic time step 184 IF( kt == nit000 .AND. neuler == 0 ) THEN ; r1_2dt_b = 1._wp / ( rdt ) 185 ELSE ; r1_2dt_b = 1._wp / ( 2._wp * rdt ) 186 ENDIF 190 187 ! 191 188 ll_init = ln_bt_av ! if no time averaging, then no specific restart … … 211 208 ll_fw_start =.FALSE. 212 209 ENDIF 213 ! 214 ! Set averaging weights and cycle length: 210 ! ! Set averaging weights and cycle length: 215 211 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 216 212 ! 217 ENDIF218 !219 IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities)220 DO jj = 2, jpjm1221 DO ji = fs_2, fs_jpim1 ! vector opt.222 zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )223 zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )224 END DO225 END DO226 ELSE ! bottom friction only227 DO jj = 2, jpjm1228 DO ji = fs_2, fs_jpim1 ! vector opt.229 zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )230 zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )231 END DO232 END DO233 ENDIF234 !235 ! Set arrays to remove/compute coriolis trend.236 ! Do it once at kt=nit000 if volume is fixed, else at each long time step.237 ! Note that these arrays are also used during barotropic loop. These are however frozen238 ! although they should be updated in the variable volume case. Not a big approximation.239 ! To remove this approximation, copy lines below inside barotropic loop240 ! and update depths at T-F points (ht and zhf resp.) at each barotropic time step241 !242 IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN243 !244 SELECT CASE( nvor_scheme )245 CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme)246 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point247 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)248 DO jj = 1, jpjm1249 DO ji = 1, jpim1250 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + &251 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp252 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)253 END DO254 END DO255 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)256 DO jj = 1, jpjm1257 DO ji = 1, jpim1258 zwz(ji,jj) = ( ht_n (ji ,jj+1) + ht_n (ji+1,jj+1) &259 & + ht_n (ji ,jj ) + ht_n (ji+1,jj ) ) &260 & / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) &261 & + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) )262 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)263 END DO264 END DO265 END SELECT266 CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )267 !268 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp269 DO jj = 2, jpj270 DO ji = 2, jpi271 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)272 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )273 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)274 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )275 END DO276 END DO277 !278 CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme)279 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp280 DO jj = 2, jpj281 DO ji = 2, jpi282 z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) )283 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht284 ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht285 ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht286 ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht287 END DO288 END DO289 !290 CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT !291 !292 zwz(:,:) = 0._wp293 zhf(:,:) = 0._wp294 295 !!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed296 !!gm A priori a better value should be something like :297 !!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)298 !!gm divided by the sum of the corresponding mask299 !!gm300 !!301 IF( .NOT.ln_sco ) THEN302 303 !!gm agree the JC comment : this should be done in a much clear way304 305 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case306 ! Set it to zero for the time being307 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level308 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth309 ! ENDIF310 ! zhf(:,:) = gdepw_0(:,:,jk+1)311 !312 ELSE313 !314 !zhf(:,:) = hbatf(:,:)315 DO jj = 1, jpjm1316 DO ji = 1, jpim1317 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) &318 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) &319 & / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) &320 & + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp )321 END DO322 END DO323 ENDIF324 !325 DO jj = 1, jpjm1326 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))327 END DO328 !329 DO jk = 1, jpkm1330 DO jj = 1, jpjm1331 zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)332 END DO333 END DO334 CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )335 ! JC: TBC. hf should be greater than 0336 DO jj = 1, jpj337 DO ji = 1, jpi338 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) ! zhf is actually hf here but it saves an array339 END DO340 END DO341 zwz(:,:) = ff_f(:,:) * zwz(:,:)342 END SELECT343 213 ENDIF 344 214 ! … … 349 219 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 350 220 ENDIF 221 ! 351 222 352 223 ! ----------------------------------------------------------------------------- … … 355 226 ! 356 227 ! 357 ! !* e3*d/dt(Ua) (Vertically integrated) 358 ! ! -------------------------------------------------- 359 zu_frc(:,:) = 0._wp 360 zv_frc(:,:) = 0._wp 361 ! 362 DO jk = 1, jpkm1 363 zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 364 zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 365 END DO 366 ! 367 zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:) 368 zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:) 369 ! 370 ! 371 ! !* baroclinic momentum trend (remove the vertical mean trend) 372 DO jk = 1, jpkm1 ! ----------------------------------------------------------- 373 DO jj = 2, jpjm1 374 DO ji = fs_2, fs_jpim1 ! vector opt. 375 ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk) 376 va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk) 377 END DO 378 END DO 228 ! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends) 229 ! ! --------------------------- ! 230 zu_frc(:,:) = SUM( e3u_n(:,:,:) * ua(:,:,:) * umask(:,:,:) , DIM=3 ) * r1_hu_n(:,:) 231 zv_frc(:,:) = SUM( e3v_n(:,:,:) * va(:,:,:) * vmask(:,:,:) , DIM=3 ) * r1_hv_n(:,:) 232 ! 233 ! 234 ! != Ua => baroclinic trend =! (remove its vertical mean) 235 DO jk = 1, jpkm1 ! ------------------------ ! 236 ua(:,:,jk) = ( ua(:,:,jk) - zu_frc(:,:) ) * umask(:,:,jk) 237 va(:,:,jk) = ( va(:,:,jk) - zv_frc(:,:) ) * vmask(:,:,jk) 379 238 END DO 380 239 … … 382 241 !!gm Is it correct to do so ? I think so... 383 242 384 ! 385 ! !* barotropic Coriolis trends (vorticity scheme dependent) 386 ! ! -------------------------------------------------------- 387 ! 388 zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes 389 zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) 390 ! 391 SELECT CASE( nvor_scheme ) 392 CASE( np_ENT ) ! enstrophy conserving scheme (f-point) 393 DO jj = 2, jpjm1 394 DO ji = 2, jpim1 ! vector opt. 395 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu_n(ji,jj) & 396 & * ( 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) ) & 397 & + e1e2t(ji ,jj)*ht_n(ji ,jj)*ff_t(ji ,jj) * ( vn_b(ji ,jj) + vn_b(ji ,jj-1) ) ) 398 ! 399 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * r1_hv_n(ji,jj) & 400 & * ( 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) ) & 401 & + e1e2t(ji,jj )*ht_n(ji,jj )*ff_t(ji,jj ) * ( un_b(ji,jj ) + un_b(ji-1,jj ) ) ) 402 END DO 403 END DO 404 ! 405 CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX 406 DO jj = 2, jpjm1 407 DO ji = fs_2, fs_jpim1 ! vector opt. 408 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 409 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 410 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 411 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 412 ! energy conserving formulation for planetary vorticity term 413 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 414 zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 415 END DO 416 END DO 417 ! 418 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 419 DO jj = 2, jpjm1 420 DO ji = fs_2, fs_jpim1 ! vector opt. 421 zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 422 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 423 zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 424 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 425 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 426 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 427 END DO 428 END DO 429 ! 430 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 431 DO jj = 2, jpjm1 432 DO ji = fs_2, fs_jpim1 ! vector opt. 433 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 434 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 435 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 436 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 437 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 438 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 439 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 440 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 441 END DO 442 END DO 443 ! 444 END SELECT 445 ! 446 ! !* Right-Hand-Side of the barotropic momentum equation 447 ! ! ---------------------------------------------------- 448 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 449 IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters 243 ! != remove 2D Coriolis and pressure gradient trends =! 244 ! ! ------------------------------------------------- ! 245 ! 246 IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2D_init ! Set zwz, the barotropic Coriolis force coefficient 247 ! ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes 248 ! 249 ! !* 2D Coriolis trends 250 zhU(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes 251 zhV(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) ! NB: FULL domain : put a value in last row and column 252 ! 253 CALL dyn_cor_2d( hu_n, hv_n, un_b, vn_b, zhU, zhV, & ! <<== in 254 & zu_trd, zv_trd ) ! ==>> out 255 ! 256 IF( .NOT.ln_linssh ) THEN !* surface pressure gradient (variable volume only) 257 ! 258 IF( ln_wd_il ) THEN ! W/D : limiter applied to spgspg 259 CALL wad_spg( sshn, zcpx, zcpy ) ! Calculating W/D gravity filters, zcpx and zcpy 450 260 DO jj = 2, jpjm1 451 DO ji = 2, jpim1 452 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 453 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 454 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 455 & > rn_wdmin1 + rn_wdmin2 456 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 457 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 458 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 459 IF(ll_tmp1) THEN 460 zcpx(ji,jj) = 1.0_wp 461 ELSEIF(ll_tmp2) THEN 462 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 463 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 464 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 465 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 466 ELSE 467 zcpx(ji,jj) = 0._wp 468 ENDIF 469 ! 470 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 471 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 472 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 473 & > rn_wdmin1 + rn_wdmin2 474 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( & 475 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 476 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 477 478 IF(ll_tmp1) THEN 479 zcpy(ji,jj) = 1.0_wp 480 ELSE IF(ll_tmp2) THEN 481 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 482 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 483 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 484 zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) ) 485 ELSE 486 zcpy(ji,jj) = 0._wp 487 ENDIF 488 END DO 489 END DO 490 ! 491 DO jj = 2, jpjm1 492 DO ji = 2, jpim1 261 DO ji = 2, jpim1 ! SPG with the application of W/D gravity filters 493 262 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 494 263 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth … … 497 266 END DO 498 267 END DO 499 ! 500 ELSE 501 ! 268 ELSE ! now suface pressure gradient 502 269 DO jj = 2, jpjm1 503 270 DO ji = fs_2, fs_jpim1 ! vector opt. … … 517 284 END DO 518 285 ! 519 ! ! Add bottom stress contribution from baroclinic velocities: 520 IF (ln_bt_fw) THEN 521 DO jj = 2, jpjm1 522 DO ji = fs_2, fs_jpim1 ! vector opt. 523 ikbu = mbku(ji,jj) 524 ikbv = mbkv(ji,jj) 525 zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities 526 zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj) 527 END DO 528 END DO 529 ELSE 530 DO jj = 2, jpjm1 531 DO ji = fs_2, fs_jpim1 ! vector opt. 532 ikbu = mbku(ji,jj) 533 ikbv = mbkv(ji,jj) 534 zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities 535 zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj) 536 END DO 537 END DO 538 ENDIF 539 ! 540 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 541 IF( ln_wd_il ) THEN 542 zztmp = -1._wp / rdtbt 543 DO jj = 2, jpjm1 544 DO ji = fs_2, fs_jpim1 ! vector opt. 545 zu_frc(ji,jj) = zu_frc(ji,jj) + & 546 & MAX(r1_hu_n(ji,jj) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ), zztmp ) * zwx(ji,jj) * wdrampu(ji,jj) 547 zv_frc(ji,jj) = zv_frc(ji,jj) + & 548 & MAX(r1_hv_n(ji,jj) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ), zztmp ) * zwy(ji,jj) * wdrampv(ji,jj) 549 END DO 550 END DO 551 ELSE 552 DO jj = 2, jpjm1 553 DO ji = fs_2, fs_jpim1 ! vector opt. 554 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zwx(ji,jj) 555 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zwy(ji,jj) 556 END DO 557 END DO 558 END IF 559 ! 560 IF( ln_isfcav ) THEN ! Add TOP stress contribution from baroclinic velocities: 561 IF( ln_bt_fw ) THEN 562 DO jj = 2, jpjm1 286 ! != Add bottom stress contribution from baroclinic velocities =! 287 ! ! ----------------------------------------------------------- ! 288 CALL dyn_drg_init( zu_frc, zv_frc, zCdU_u, zCdU_v ) ! also provide the barotropic drag coefficients 289 ! 290 ! != Add atmospheric pressure forcing =! 291 ! ! ---------------------------------- ! 292 IF( ln_apr_dyn ) THEN 293 IF( ln_bt_fw ) THEN ! FORWARD integration: use kt+1/2 pressure (NOW+1/2) 294 DO jj = 2, jpjm1 563 295 DO ji = fs_2, fs_jpim1 ! vector opt. 564 iktu = miku(ji,jj) 565 iktv = mikv(ji,jj) 566 zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities 567 zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 568 END DO 569 END DO 570 ELSE 571 DO jj = 2, jpjm1 296 zu_frc(ji,jj) = zu_frc(ji,jj) + grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 297 zv_frc(ji,jj) = zv_frc(ji,jj) + grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 298 END DO 299 END DO 300 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW) 301 zztmp = grav * r1_2 302 DO jj = 2, jpjm1 572 303 DO ji = fs_2, fs_jpim1 ! vector opt. 573 iktu = miku(ji,jj) 574 iktv = mikv(ji,jj) 575 zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities 576 zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 577 END DO 578 END DO 579 ENDIF 580 ! 581 ! Note that the "unclipped" top friction parameter is used even with explicit drag 582 DO jj = 2, jpjm1 583 DO ji = fs_2, fs_jpim1 ! vector opt. 584 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zwx(ji,jj) 585 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zwy(ji,jj) 586 END DO 587 END DO 588 ENDIF 589 ! 304 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 305 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 306 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 307 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 308 END DO 309 END DO 310 ENDIF 311 ENDIF 312 ! 313 ! != Add atmospheric pressure forcing =! 314 ! ! ---------------------------------- ! 590 315 IF( ln_bt_fw ) THEN ! Add wind forcing 591 316 DO jj = 2, jpjm1 … … 605 330 ENDIF 606 331 ! 607 IF( ln_apr_dyn ) THEN ! Add atm pressure forcing 608 IF( ln_bt_fw ) THEN 609 DO jj = 2, jpjm1 610 DO ji = fs_2, fs_jpim1 ! vector opt. 611 zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 612 zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 613 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 614 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg 615 END DO 616 END DO 617 ELSE 618 zztmp = grav * r1_2 619 DO jj = 2, jpjm1 620 DO ji = fs_2, fs_jpim1 ! vector opt. 621 zu_spg = zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 622 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 623 zv_spg = zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 624 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 625 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 626 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg 627 END DO 628 END DO 629 ENDIF 630 ENDIF 631 ! !* Right-Hand-Side of the barotropic ssh equation 632 ! ! ----------------------------------------------- 633 ! ! Surface net water flux, rivers and ice shelves 634 IF (ln_bt_fw) THEN 635 zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) 636 ELSE 332 ! !----------------! 333 ! !== sssh_frc ==! Right-Hand-Side of the barotropic ssh equation (over the FULL domain) 334 ! !----------------! 335 ! != Net water flux forcing applied to a water column =! 336 ! ! --------------------------------------------------- ! 337 IF (ln_bt_fw) THEN ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2) 338 zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) 339 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW) 637 340 zztmp = r1_rau0 * r1_2 638 341 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) & … … 641 344 & + fwfisf_par(:,:) + fwfisf_par_b(:,:) ) 642 345 ENDIF 643 ! 644 IF( ln_sdw ) THEN ! Stokes drift divergence added if necessary346 ! != Add Stokes drift divergence =! (if exist) 347 IF( ln_sdw ) THEN ! ----------------------------- ! 645 348 zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:) 646 349 ENDIF … … 662 365 ! 663 366 #if defined key_asminc 664 ! ! Include the IAU weighted SSH increment 367 ! != Add the IAU weighted SSH increment =! 368 ! ! ------------------------------------ ! 665 369 IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 666 370 zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:) 667 371 ENDIF 668 372 #endif 669 ! ! *Fill boundary data arrays for AGRIF373 ! != Fill boundary data arrays for AGRIF 670 374 ! ! ------------------------------------ 671 375 #if defined key_agrif … … 689 393 vb_e (:,:) = 0._wp 690 394 ENDIF 691 395 ! 396 IF( ln_linssh ) THEN ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0) 397 zhup2_e(:,:) = hu_n(:,:) 398 zhvp2_e(:,:) = hv_n(:,:) 399 zhtp2_e(:,:) = ht_n(:,:) 400 ENDIF 692 401 ! 693 402 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields … … 711 420 ENDIF 712 421 ! 713 !714 !715 422 ! Initialize sums: 716 423 ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form) … … 732 439 ! 733 440 l_full_nf_update = jn == icycle ! false: disable full North fold update (performances) for jn = 1 to icycle-1 734 ! ! ------------------ 735 ! !* Update the forcing (BDY and tides) 736 ! ! ------------------ 737 ! Update only tidal forcing at open boundaries 738 IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 ) 739 IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset ) 740 ! 741 ! Set extrapolation coefficients for predictor step: 441 ! 442 ! !== Update the forcing ==! (BDY and tides) 443 ! 444 IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, kt_offset= noffset+1 ) 445 IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, kt_offset= noffset ) 446 ! 447 ! !== extrapolation at mid-step ==! (jn+1/2) 448 ! 449 ! !* Set extrapolation coefficients for predictor step: 742 450 IF ((jn<3).AND.ll_init) THEN ! Forward 743 451 za1 = 1._wp … … 749 457 za3 = 0.281105_wp ! za3 = bet 750 458 ENDIF 751 752 ! Extrapolate barotropic velocities at step jit+0.5: 459 ! 460 ! !* Extrapolate barotropic velocities at mid-step (jn+1/2) 461 !-- m+1/2 m m-1 m-2 --! 462 !-- u = (3/2+beta) u -(1/2+2beta) u + beta u --! 463 !-------------------------------------------------------------------------! 753 464 ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 754 465 va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) … … 757 468 ! ! ------------------ 758 469 ! Extrapolate Sea Level at step jit+0.5: 470 !-- m+1/2 m m-1 m-2 --! 471 !-- ssh = (3/2+beta) ssh -(1/2+2beta) ssh + beta ssh --! 472 !--------------------------------------------------------------------------------! 759 473 zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 760 474 761 ! set wetting & drying mask at tracer points for this barotropic sub-step 762 IF ( ln_wd_dl ) THEN 763 ! 764 IF ( ln_wd_dl_rmp ) THEN 765 DO jj = 1, jpj 766 DO ji = 1, jpi ! vector opt. 767 IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 768 ! IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 769 ztwdmask(ji,jj) = 1._wp 770 ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 771 ztwdmask(ji,jj) = (tanh(50._wp*( ( zsshp2_e(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1)) ) 772 ELSE 773 ztwdmask(ji,jj) = 0._wp 774 END IF 775 END DO 776 END DO 777 ELSE 778 DO jj = 1, jpj 779 DO ji = 1, jpi ! vector opt. 780 IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 781 ztwdmask(ji,jj) = 1._wp 782 ELSE 783 ztwdmask(ji,jj) = 0._wp 784 ENDIF 785 END DO 786 END DO 787 ENDIF 788 ! 789 ENDIF 475 ! set wetting & drying mask at tracer points for this barotropic mid-step 476 IF( ln_wd_dl ) CALL wad_tmsk( zsshp2_e, ztwdmask ) 790 477 ! 791 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 792 DO ji = 2, fs_jpim1 ! Vector opt. 793 zwx(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 794 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 795 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 796 zwy(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 797 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 798 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 799 END DO 800 END DO 801 CALL lbc_lnk_multi( 'dynspg_ts', zwx, 'U', 1._wp, zwy, 'V', 1._wp ) 478 ! ! ocean t-depth at mid-step 479 zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:) 802 480 ! 803 zhup2_e(:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 804 zhvp2_e(:,:) = hv_0(:,:) + zwy(:,:) 805 zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:) 806 ELSE 807 zhup2_e(:,:) = hu_n(:,:) 808 zhvp2_e(:,:) = hv_n(:,:) 809 zhtp2_e(:,:) = ht_n(:,:) 810 ENDIF 811 ! !* after ssh 812 ! ! ----------- 813 ! 814 ! Enforce volume conservation at open boundaries: 481 ! ! ocean u- and v-depth at mid-step (separate DO-loops remove the need of a lbc_lnk) 482 DO jj = 1, jpj 483 DO ji = 1, jpim1 ! not jpi-column 484 zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) & 485 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 486 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 487 END DO 488 END DO 489 DO jj = 1, jpjm1 ! not jpj-row 490 DO ji = 1, jpi 491 zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) & 492 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 493 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 494 END DO 495 END DO 496 ! 497 ENDIF 498 ! 499 ! !== after SSH ==! (jn+1) 500 ! 501 ! ! update (ua_e,va_e) to enforce volume conservation at open boundaries 502 ! ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d 815 503 IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e ) 816 504 ! 817 zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5 818 zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 505 ! ! resulting flux at mid-step (not over the full domain) 506 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 507 zhV(1:jpi ,1:jpjm1) = e1v(1:jpi ,1:jpjm1) * va_e(1:jpi ,1:jpjm1) * zhvp2_e(1:jpi ,1:jpjm1) ! not jpj-row 819 508 ! 820 509 #if defined key_agrif … … 823 512 IF((nbondi == -1).OR.(nbondi == 2)) THEN 824 513 DO jj = 1, jpj 825 z wx(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj)826 z wy(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj)514 zhU(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj) 515 zhV(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj) 827 516 END DO 828 517 ENDIF 829 518 IF((nbondi == 1).OR.(nbondi == 2)) THEN 830 519 DO jj=1,jpj 831 z wx(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj)832 z wy(nlci-nbghostcells :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells :nlci-1,jj)520 zhU(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj) 521 zhV(nlci-nbghostcells :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells :nlci-1,jj) 833 522 END DO 834 523 ENDIF 835 524 IF((nbondj == -1).OR.(nbondj == 2)) THEN 836 525 DO ji=1,jpi 837 z wy(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1)838 z wx(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1)526 zhV(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1) 527 zhU(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1) 839 528 END DO 840 529 ENDIF 841 530 IF((nbondj == 1).OR.(nbondj == 2)) THEN 842 531 DO ji=1,jpi 843 z wy(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2)844 z wx(ji,nlcj-nbghostcells :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells :nlcj-1)532 zhV(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2) 533 zhU(ji,nlcj-nbghostcells :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells :nlcj-1) 845 534 END DO 846 535 ENDIF 847 536 ENDIF 848 537 #endif 849 IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 850 851 IF ( ln_wd_dl ) THEN 852 ! 853 ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells 854 ! 855 DO jj = 1, jpjm1 856 DO ji = 1, jpim1 857 IF ( zwx(ji,jj) > 0.0 ) THEN 858 zuwdmask(ji, jj) = ztwdmask(ji ,jj) 859 ELSE 860 zuwdmask(ji, jj) = ztwdmask(ji+1,jj) 861 END IF 862 zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj) 863 un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj) 864 865 IF ( zwy(ji,jj) > 0.0 ) THEN 866 zvwdmask(ji, jj) = ztwdmask(ji, jj ) 867 ELSE 868 zvwdmask(ji, jj) = ztwdmask(ji, jj+1) 869 END IF 870 zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj) 871 vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj) 872 END DO 873 END DO 538 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 539 540 IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where 541 ! ! the direction of the flow is from dry cells 542 CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) ! not jpi colomn for U, not jpj row for V 874 543 ! 875 544 ENDIF 876 877 ! Sum over sub-time-steps to compute advective velocities 878 za2 = wgtbtp2(jn) 879 un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 880 vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 881 882 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True) 545 ! 546 ! 547 ! Compute Sea Level at step jit+1 548 !-- m+1 m m+1/2 --! 549 !-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --! 550 !-------------------------------------------------------------------------! 551 DO jj = 2, jpjm1 ! INNER domain 552 DO ji = 2, jpim1 553 zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj) 554 ssha_e(ji,jj) = ( sshn_e(ji,jj) - rdtbt * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj) 555 END DO 556 END DO 557 ! 558 CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp ) 559 ! 560 ! ! Sum over sub-time-steps to compute advective velocities 561 za2 = wgtbtp2(jn) ! zhU, zhV hold fluxes extrapolated at jn+0.5 562 un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:) 563 vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:) 564 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True) 883 565 IF ( ln_wd_dl_bc ) THEN 884 zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:) 885 zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:) 886 END IF 887 888 ! Set next sea level: 889 DO jj = 2, jpjm1 890 DO ji = fs_2, fs_jpim1 ! vector opt. 891 zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) & 892 & + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj) 893 END DO 894 END DO 895 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 896 897 CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T', 1._wp ) 898 566 zuwdav2(1:jpim1,1:jpj ) = zuwdav2(1:jpim1,1:jpj ) + za2 * zuwdmask(1:jpim1,1:jpj ) ! not jpi-column 567 zvwdav2(1:jpi ,1:jpjm1) = zvwdav2(1:jpi ,1:jpjm1) + za2 * zvwdmask(1:jpi ,1:jpjm1) ! not jpj-row 568 END IF 569 ! 899 570 ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) 900 571 IF( ln_bdy ) CALL bdy_ssh( ssha_e ) … … 905 576 ! Sea Surface Height at u-,v-points (vvl case only) 906 577 IF( .NOT.ln_linssh ) THEN 907 DO jj = 2, jpjm1 578 DO jj = 2, jpjm1 ! INNER domain, will be extended to whole domain later 908 579 DO ji = 2, jpim1 ! NO Vector Opt. 909 580 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & … … 915 586 END DO 916 587 END DO 917 CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )918 588 ENDIF 919 ! 920 ! Half-step back interpolation of SSH for surface pressure computation: 921 !---------------------------------------------------------------------- 922 IF ((jn==1).AND.ll_init) THEN 923 za0=1._wp ! Forward-backward 924 za1=0._wp 925 za2=0._wp 926 za3=0._wp 927 ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12 928 za0= 1.0833333333333_wp ! za0 = 1-gam-eps 929 za1=-0.1666666666666_wp ! za1 = gam 930 za2= 0.0833333333333_wp ! za2 = eps 931 za3= 0._wp 932 ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 933 IF (rn_bt_alpha==0._wp) THEN 934 za0=0.614_wp ! za0 = 1/2 + gam + 2*eps 935 za1=0.285_wp ! za1 = 1/2 - 2*gam - 3*eps 936 za2=0.088_wp ! za2 = gam 937 za3=0.013_wp ! za3 = eps 938 ELSE 939 zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha 940 zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha 941 za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon 942 za1 = 1._wp - za0 - zgamma - zepsilon 943 za2 = zgamma 944 za3 = zepsilon 945 ENDIF 946 ENDIF 947 ! 589 ! 590 ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2 591 !-- m+1/2 m+1 m m-1 m-2 --! 592 !-- ssh' = za0 * ssh + za1 * ssh + za2 * ssh + za3 * ssh --! 593 !------------------------------------------------------------------------------------------! 594 CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 ) ! coeficients of the interpolation 948 595 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 949 596 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 950 951 IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters 952 DO jj = 2, jpjm1 953 DO ji = 2, jpim1 954 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 955 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 956 & MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) & 957 & > rn_wdmin1 + rn_wdmin2 958 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji+1,jj)) > 1.E-12 ).AND.( & 959 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 960 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 961 962 IF(ll_tmp1) THEN 963 zcpx(ji,jj) = 1.0_wp 964 ELSE IF(ll_tmp2) THEN 965 ! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here 966 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 967 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) ) 968 ELSE 969 zcpx(ji,jj) = 0._wp 970 ENDIF 971 ! 972 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 973 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 974 & MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) & 975 & > rn_wdmin1 + rn_wdmin2 976 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji,jj+1)) > 1.E-12 ).AND.( & 977 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 978 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 979 980 IF(ll_tmp1) THEN 981 zcpy(ji,jj) = 1.0_wp 982 ELSEIF(ll_tmp2) THEN 983 ! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here 984 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 985 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) ) 986 ELSE 987 zcpy(ji,jj) = 0._wp 988 ENDIF 989 END DO 990 END DO 991 ENDIF 992 ! 993 ! Compute associated depths at U and V points: 994 IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form 995 ! 996 DO jj = 2, jpjm1 997 DO ji = 2, jpim1 998 zx1 = r1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) & 999 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) & 1000 & + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) ) 1001 zy1 = r1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) & 1002 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) & 1003 & + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) ) 1004 zhust_e(ji,jj) = hu_0(ji,jj) + zx1 1005 zhvst_e(ji,jj) = hv_0(ji,jj) + zy1 1006 END DO 1007 END DO 1008 ! 597 ! 598 ! ! Surface pressure gradient 599 zldg = ( 1._wp - rn_scal_load ) * grav ! local factor 600 DO jj = 2, jpjm1 601 DO ji = 2, jpim1 602 zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 603 zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 604 END DO 605 END DO 606 IF( ln_wd_il ) THEN ! W/D : gravity filters applied on pressure gradient 607 CALL wad_spg( zsshp2_e, zcpx, zcpy ) ! Calculating W/D gravity filters 608 zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1) 609 zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1) 1009 610 ENDIF 1010 611 ! … … 1012 613 ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated 1013 614 ! at each time step. We however keep them constant here for optimization. 1014 ! Recall that zwx and zwy arrays hold fluxes at this stage: 1015 ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5 1016 ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 1017 ! 1018 SELECT CASE( nvor_scheme ) 1019 CASE( np_ENT ) ! energy conserving scheme (t-point) 1020 DO jj = 2, jpjm1 1021 DO ji = 2, jpim1 ! vector opt. 1022 1023 z1_hu = ssumask(ji,jj) / ( zhup2_e(ji,jj) + 1._wp - ssumask(ji,jj) ) 1024 z1_hv = ssvmask(ji,jj) / ( zhvp2_e(ji,jj) + 1._wp - ssvmask(ji,jj) ) 1025 1026 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu & 1027 & * ( e1e2t(ji+1,jj)*zhtp2_e(ji+1,jj)*ff_t(ji+1,jj) * ( va_e(ji+1,jj) + va_e(ji+1,jj-1) ) & 1028 & + e1e2t(ji ,jj)*zhtp2_e(ji ,jj)*ff_t(ji ,jj) * ( va_e(ji ,jj) + va_e(ji ,jj-1) ) ) 1029 ! 1030 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1031 & * ( e1e2t(ji,jj+1)*zhtp2_e(ji,jj+1)*ff_t(ji,jj+1) * ( ua_e(ji,jj+1) + ua_e(ji-1,jj+1) ) & 1032 & + e1e2t(ji,jj )*zhtp2_e(ji,jj )*ff_t(ji,jj ) * ( ua_e(ji,jj ) + ua_e(ji-1,jj ) ) ) 1033 END DO 1034 END DO 1035 ! 1036 CASE( np_ENE, np_MIX ) ! energy conserving scheme (f-point) 1037 DO jj = 2, jpjm1 1038 DO ji = fs_2, fs_jpim1 ! vector opt. 1039 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 1040 zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 1041 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 1042 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 1043 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 1044 zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 1045 END DO 1046 END DO 1047 ! 1048 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 1049 DO jj = 2, jpjm1 1050 DO ji = fs_2, fs_jpim1 ! vector opt. 1051 zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 1052 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 1053 zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 1054 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 1055 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 1056 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 1057 END DO 1058 END DO 1059 ! 1060 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1061 DO jj = 2, jpjm1 1062 DO ji = fs_2, fs_jpim1 ! vector opt. 1063 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 1064 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 1065 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 1066 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 1067 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 1068 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 1069 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 1070 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 1071 END DO 1072 END DO 1073 ! 1074 END SELECT 615 ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated) 616 CALL dyn_cor_2d( zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd ) 1075 617 ! 1076 618 ! Add tidal astronomical forcing if defined … … 1078 620 DO jj = 2, jpjm1 1079 621 DO ji = fs_2, fs_jpim1 ! vector opt. 1080 zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 1081 zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 1082 zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg 1083 zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg 622 zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 623 zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 1084 624 END DO 1085 625 END DO … … 1095 635 END DO 1096 636 END DO 1097 ENDIF 1098 ! 1099 ! Surface pressure trend: 1100 IF( ln_wd_il ) THEN 1101 DO jj = 2, jpjm1 1102 DO ji = 2, jpim1 1103 ! Add surface pressure gradient 1104 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 1105 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 1106 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj) 1107 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj) 1108 END DO 1109 END DO 1110 ELSE 1111 DO jj = 2, jpjm1 1112 DO ji = fs_2, fs_jpim1 ! vector opt. 1113 ! Add surface pressure gradient 1114 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 1115 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 1116 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg 1117 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg 1118 END DO 1119 END DO 1120 END IF 1121 637 ENDIF 1122 638 ! 1123 639 ! Set next velocities: 640 ! Compute barotropic speeds at step jit+1 (h : total height of the water colomn) 641 !-- VECTOR FORM 642 !-- m+1 m / m+1/2 \ --! 643 !-- u = u + delta_t' * \ (1-r)*g * grad_x( ssh') - f * k vect u + frc / --! 644 !-- --! 645 !-- FLUX FORM --! 646 !-- m+1 __1__ / m m / m+1/2 m+1/2 m+1/2 n \ \ --! 647 !-- u = m+1 | h * u + delta_t' * \ h * (1-r)*g * grad_x( ssh') - h * f * k vect u + h * frc / | --! 648 !-- h \ / --! 649 !------------------------------------------------------------------------------------------------------------------------! 1124 650 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 1125 651 DO jj = 2, jpjm1 1126 652 DO ji = fs_2, fs_jpim1 ! vector opt. 1127 653 ua_e(ji,jj) = ( un_e(ji,jj) & 1128 & + rdtbt * ( zwx(ji,jj) &654 & + rdtbt * ( zu_spg(ji,jj) & 1129 655 & + zu_trd(ji,jj) & 1130 656 & + zu_frc(ji,jj) ) & … … 1132 658 1133 659 va_e(ji,jj) = ( vn_e(ji,jj) & 1134 & + rdtbt * ( zwy(ji,jj) &660 & + rdtbt * ( zv_spg(ji,jj) & 1135 661 & + zv_trd(ji,jj) & 1136 662 & + zv_frc(ji,jj) ) & 1137 663 & ) * ssvmask(ji,jj) 1138 1139 664 END DO 1140 665 END DO … … 1142 667 ELSE !* Flux form 1143 668 DO jj = 2, jpjm1 1144 DO ji = fs_2, fs_jpim1 ! vector opt. 1145 1146 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 1147 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 1148 1149 zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj)) 1150 zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj)) 1151 1152 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & 1153 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 1154 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 1155 & + hu_n(ji,jj) * zu_frc(ji,jj) ) & 1156 & ) * zhura 1157 1158 va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) & 1159 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 1160 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & 1161 & + hv_n(ji,jj) * zv_frc(ji,jj) ) & 1162 & ) * zhvra 669 DO ji = 2, jpim1 670 ! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2 671 ! ! backward interpolated depth used in spg terms at jn+1/2 672 zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 673 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 674 zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 675 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 676 ! ! inverse depth at jn+1 677 z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) 678 z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) ) 679 ! 680 ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) & 681 & + rdtbt * ( zhu_bck * zu_spg (ji,jj) & ! 682 & + zhup2_e(ji,jj) * zu_trd (ji,jj) & ! 683 & + hu_n (ji,jj) * zu_frc (ji,jj) ) ) * z1_hu 684 ! 685 va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) & 686 & + rdtbt * ( zhv_bck * zv_spg (ji,jj) & ! 687 & + zhvp2_e(ji,jj) * zv_trd (ji,jj) & ! 688 & + hv_n (ji,jj) * zv_frc (ji,jj) ) ) * z1_hv 1163 689 END DO 1164 690 END DO … … 1173 699 END DO 1174 700 ENDIF 1175 1176 1177 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 1178 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 1179 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 1180 hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 1181 hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) 1182 ! 1183 ENDIF 1184 ! !* domain lateral boundary 1185 CALL lbc_lnk_multi( 'dynspg_ts', ua_e, 'U', -1._wp, va_e , 'V', -1._wp ) 701 702 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 703 hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1) 704 hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) ) 705 hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1) 706 hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) ) 707 CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp & 708 & , hu_e , 'U', -1._wp, hv_e , 'V', -1._wp & 709 & , hur_e, 'U', -1._wp, hvr_e, 'V', -1._wp ) 710 ELSE 711 CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp ) 712 ENDIF 713 ! 1186 714 ! 1187 715 ! ! open boundaries … … 1231 759 ! Set advection velocity correction: 1232 760 IF (ln_bt_fw) THEN 1233 zwx(:,:) = un_adv(:,:)1234 zwy(:,:) = vn_adv(:,:)1235 761 IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN 1236 un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) ) 1237 vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) ) 1238 ! 1239 ! Update corrective fluxes for next time step: 1240 un_bf(:,:) = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:)) 1241 vn_bf(:,:) = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:)) 762 DO jj = 1, jpj 763 DO ji = 1, jpi 764 zun_save = un_adv(ji,jj) 765 zvn_save = vn_adv(ji,jj) 766 ! ! apply the previously computed correction 767 un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) ) 768 vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) ) 769 ! ! Update corrective fluxes for next time step 770 un_bf(ji,jj) = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) ) 771 vn_bf(ji,jj) = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) ) 772 ! ! Save integrated transport for next computation 773 ub2_b(ji,jj) = zun_save 774 vb2_b(ji,jj) = zvn_save 775 END DO 776 END DO 1242 777 ELSE 1243 un_bf(:,:) = 0._wp 1244 vn_bf(:,:) = 0._wp 1245 END IF 1246 ! Save integrated transport for next computation 1247 ub2_b(:,:) = zwx(:,:) 1248 vb2_b(:,:) = zwy(:,:) 778 un_bf(:,:) = 0._wp ! corrective fluxes for next time step set to zero 779 vn_bf(:,:) = 0._wp 780 ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation 781 vb2_b(:,:) = vn_adv(:,:) 782 END IF 1249 783 ENDIF 1250 784 … … 1491 1025 REAL(wp) :: zxr2, zyr2, zcmax ! local scalar 1492 1026 REAL(wp), DIMENSION(jpi,jpj) :: zcu 1027 INTEGER :: inum 1493 1028 !!---------------------------------------------------------------------- 1494 1029 ! … … 1597 1132 END SUBROUTINE dyn_spg_ts_init 1598 1133 1134 1135 SUBROUTINE dyn_cor_2d_init 1136 !!--------------------------------------------------------------------- 1137 !! *** ROUTINE dyn_cor_2d_init *** 1138 !! 1139 !! ** Purpose : Set time splitting options 1140 !! Set arrays to remove/compute coriolis trend. 1141 !! Do it once during initialization if volume is fixed, else at each long time step. 1142 !! Note that these arrays are also used during barotropic loop. These are however frozen 1143 !! although they should be updated in the variable volume case. Not a big approximation. 1144 !! To remove this approximation, copy lines below inside barotropic loop 1145 !! and update depths at T-F points (ht and zhf resp.) at each barotropic time step 1146 !! 1147 !! Compute zwz = f / ( height of the water colomn ) 1148 !!---------------------------------------------------------------------- 1149 INTEGER :: ji ,jj, jk ! dummy loop indices 1150 REAL(wp) :: z1_ht 1151 REAL(wp), DIMENSION(jpi,jpj) :: zhf 1152 !!---------------------------------------------------------------------- 1153 ! 1154 SELECT CASE( nvor_scheme ) 1155 CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme) 1156 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point 1157 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 1158 DO jj = 1, jpjm1 1159 DO ji = 1, jpim1 1160 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 1161 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp 1162 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1163 END DO 1164 END DO 1165 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 1166 DO jj = 1, jpjm1 1167 DO ji = 1, jpim1 1168 zwz(ji,jj) = ( ht_n (ji ,jj+1) + ht_n (ji+1,jj+1) & 1169 & + ht_n (ji ,jj ) + ht_n (ji+1,jj ) ) & 1170 & / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) & 1171 & + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) ) 1172 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1173 END DO 1174 END DO 1175 END SELECT 1176 CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp ) 1177 ! 1178 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1179 DO jj = 2, jpj 1180 DO ji = 2, jpi 1181 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 1182 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 1183 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 1184 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 1185 END DO 1186 END DO 1187 ! 1188 CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme) 1189 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1190 DO jj = 2, jpj 1191 DO ji = 2, jpi 1192 z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) ) 1193 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht 1194 ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht 1195 ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht 1196 ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht 1197 END DO 1198 END DO 1199 ! 1200 CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT ! 1201 ! 1202 zwz(:,:) = 0._wp 1203 zhf(:,:) = 0._wp 1204 1205 !!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed 1206 !!gm A priori a better value should be something like : 1207 !!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1) 1208 !!gm divided by the sum of the corresponding mask 1209 !!gm 1210 !! 1211 IF( .NOT.ln_sco ) THEN 1212 1213 !!gm agree the JC comment : this should be done in a much clear way 1214 1215 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 1216 ! Set it to zero for the time being 1217 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level 1218 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth 1219 ! ENDIF 1220 ! zhf(:,:) = gdepw_0(:,:,jk+1) 1221 ! 1222 ELSE 1223 ! 1224 !zhf(:,:) = hbatf(:,:) 1225 DO jj = 1, jpjm1 1226 DO ji = 1, jpim1 1227 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) & 1228 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) & 1229 & / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) & 1230 & + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp ) 1231 END DO 1232 END DO 1233 ENDIF 1234 ! 1235 DO jj = 1, jpjm1 1236 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 1237 END DO 1238 ! 1239 DO jk = 1, jpkm1 1240 DO jj = 1, jpjm1 1241 zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 1242 END DO 1243 END DO 1244 CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp ) 1245 ! JC: TBC. hf should be greater than 0 1246 DO jj = 1, jpj 1247 DO ji = 1, jpi 1248 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) 1249 END DO 1250 END DO 1251 zwz(:,:) = ff_f(:,:) * zwz(:,:) 1252 END SELECT 1253 1254 END SUBROUTINE dyn_cor_2d_init 1255 1256 1257 1258 SUBROUTINE dyn_cor_2d( hu_n, hv_n, un_b, vn_b, zhU, zhV, zu_trd, zv_trd ) 1259 !!--------------------------------------------------------------------- 1260 !! *** ROUTINE dyn_cor_2d *** 1261 !! 1262 !! ** Purpose : Compute u and v coriolis trends 1263 !!---------------------------------------------------------------------- 1264 INTEGER :: ji ,jj ! dummy loop indices 1265 REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - - 1266 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: hu_n, hv_n, un_b, vn_b, zhU, zhV 1267 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd 1268 !!---------------------------------------------------------------------- 1269 SELECT CASE( nvor_scheme ) 1270 CASE( np_ENT ) ! enstrophy conserving scheme (f-point) 1271 DO jj = 2, jpjm1 1272 DO ji = 2, jpim1 1273 z1_hu = ssumask(ji,jj) / ( hu_n(ji,jj) + 1._wp - ssumask(ji,jj) ) 1274 z1_hv = ssvmask(ji,jj) / ( hv_n(ji,jj) + 1._wp - ssvmask(ji,jj) ) 1275 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu & 1276 & * ( 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) ) & 1277 & + e1e2t(ji ,jj)*ht_n(ji ,jj)*ff_t(ji ,jj) * ( vn_b(ji ,jj) + vn_b(ji ,jj-1) ) ) 1278 ! 1279 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1280 & * ( 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) ) & 1281 & + e1e2t(ji,jj )*ht_n(ji,jj )*ff_t(ji,jj ) * ( un_b(ji,jj ) + un_b(ji-1,jj ) ) ) 1282 END DO 1283 END DO 1284 ! 1285 CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX 1286 DO jj = 2, jpjm1 1287 DO ji = fs_2, fs_jpim1 ! vector opt. 1288 zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj) 1289 zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) 1290 zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj) 1291 zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj) 1292 ! energy conserving formulation for planetary vorticity term 1293 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 1294 zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 1295 END DO 1296 END DO 1297 ! 1298 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 1299 DO jj = 2, jpjm1 1300 DO ji = fs_2, fs_jpim1 ! vector opt. 1301 zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) & 1302 & + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) 1303 zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) & 1304 & + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj) 1305 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 1306 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 1307 END DO 1308 END DO 1309 ! 1310 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1311 DO jj = 2, jpjm1 1312 DO ji = fs_2, fs_jpim1 ! vector opt. 1313 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) & 1314 & + ftnw(ji+1,jj) * zhV(ji+1,jj ) & 1315 & + ftse(ji,jj ) * zhV(ji ,jj-1) & 1316 & + ftsw(ji+1,jj) * zhV(ji+1,jj-1) ) 1317 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) & 1318 & + ftse(ji,jj+1) * zhU(ji ,jj+1) & 1319 & + ftnw(ji,jj ) * zhU(ji-1,jj ) & 1320 & + ftne(ji,jj ) * zhU(ji ,jj ) ) 1321 END DO 1322 END DO 1323 ! 1324 END SELECT 1325 ! 1326 END SUBROUTINE dyn_cor_2D 1327 1328 1329 SUBROUTINE wad_tmsk( pssh, ptmsk ) 1330 !!---------------------------------------------------------------------- 1331 !! *** ROUTINE wad_lmt *** 1332 !! 1333 !! ** Purpose : set wetting & drying mask at tracer points 1334 !! for the current barotropic sub-step 1335 !! 1336 !! ** Method : ??? 1337 !! 1338 !! ** Action : ptmsk : wetting & drying t-mask 1339 !!---------------------------------------------------------------------- 1340 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh ! 1341 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: ptmsk ! 1342 ! 1343 INTEGER :: ji, jj ! dummy loop indices 1344 !!---------------------------------------------------------------------- 1345 ! 1346 IF( ln_wd_dl_rmp ) THEN 1347 DO jj = 1, jpj 1348 DO ji = 1, jpi 1349 IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 1350 ! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 1351 ptmsk(ji,jj) = 1._wp 1352 ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 1353 ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1) ) 1354 ELSE 1355 ptmsk(ji,jj) = 0._wp 1356 ENDIF 1357 END DO 1358 END DO 1359 ELSE 1360 DO jj = 1, jpj 1361 DO ji = 1, jpi 1362 IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp 1363 ELSE ; ptmsk(ji,jj) = 0._wp 1364 ENDIF 1365 END DO 1366 END DO 1367 ENDIF 1368 ! 1369 END SUBROUTINE wad_tmsk 1370 1371 1372 SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk ) 1373 !!---------------------------------------------------------------------- 1374 !! *** ROUTINE wad_lmt *** 1375 !! 1376 !! ** Purpose : set wetting & drying mask at tracer points 1377 !! for the current barotropic sub-step 1378 !! 1379 !! ** Method : ??? 1380 !! 1381 !! ** Action : ptmsk : wetting & drying t-mask 1382 !!---------------------------------------------------------------------- 1383 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pTmsk ! W & D t-mask 1384 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports 1385 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask 1386 ! 1387 INTEGER :: ji, jj ! dummy loop indices 1388 !!---------------------------------------------------------------------- 1389 ! 1390 DO jj = 1, jpj 1391 DO ji = 1, jpim1 ! not jpi-column 1392 IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj) 1393 ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj) 1394 ENDIF 1395 phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj) 1396 pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj) 1397 END DO 1398 END DO 1399 ! 1400 DO jj = 1, jpjm1 ! not jpj-row 1401 DO ji = 1, jpi 1402 IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj ) 1403 ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1) 1404 ENDIF 1405 phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 1406 pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj) 1407 END DO 1408 END DO 1409 ! 1410 END SUBROUTINE wad_Umsk 1411 1412 1413 SUBROUTINE wad_spg( sshn, zcpx, zcpy ) 1414 !!--------------------------------------------------------------------- 1415 !! *** ROUTINE wad_sp *** 1416 !! 1417 !! ** Purpose : 1418 !!---------------------------------------------------------------------- 1419 INTEGER :: ji ,jj ! dummy loop indices 1420 LOGICAL :: ll_tmp1, ll_tmp2 1421 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: sshn 1422 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy 1423 !!---------------------------------------------------------------------- 1424 DO jj = 2, jpjm1 1425 DO ji = 2, jpim1 1426 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1427 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 1428 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 1429 & > rn_wdmin1 + rn_wdmin2 1430 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 1431 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1432 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 1433 IF(ll_tmp1) THEN 1434 zcpx(ji,jj) = 1.0_wp 1435 ELSEIF(ll_tmp2) THEN 1436 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 1437 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 1438 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 1439 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 1440 ELSE 1441 zcpx(ji,jj) = 0._wp 1442 ENDIF 1443 ! 1444 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1445 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 1446 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 1447 & > rn_wdmin1 + rn_wdmin2 1448 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( & 1449 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1450 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1451 1452 IF(ll_tmp1) THEN 1453 zcpy(ji,jj) = 1.0_wp 1454 ELSE IF(ll_tmp2) THEN 1455 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 1456 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 1457 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 1458 zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) ) 1459 ELSE 1460 zcpy(ji,jj) = 0._wp 1461 ENDIF 1462 END DO 1463 END DO 1464 1465 END SUBROUTINE wad_spg 1466 1467 1468 1469 SUBROUTINE dyn_drg_init( pu_RHSi, pv_RHSi, pCdU_u, pCdU_v ) 1470 !!---------------------------------------------------------------------- 1471 !! *** ROUTINE dyn_drg_init *** 1472 !! 1473 !! ** Purpose : - add the baroclinic top/bottom drag contribution to 1474 !! the baroclinic part of the barotropic RHS 1475 !! - compute the barotropic drag coefficients 1476 !! 1477 !! ** Method : computation done over the INNER domain only 1478 !!---------------------------------------------------------------------- 1479 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS 1480 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients 1481 ! 1482 INTEGER :: ji, jj ! dummy loop indices 1483 INTEGER :: ikbu, ikbv, iktu, iktv 1484 REAL(wp) :: zztmp 1485 REAL(wp), DIMENSION(jpi,jpj) :: zu_i, zv_i 1486 !!---------------------------------------------------------------------- 1487 ! 1488 ! !== Set the barotropic drag coef. ==! 1489 ! 1490 IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities) 1491 1492 DO jj = 2, jpjm1 1493 DO ji = 2, jpim1 ! INNER domain 1494 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 1495 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) 1496 END DO 1497 END DO 1498 ELSE ! bottom friction only 1499 DO jj = 2, jpjm1 1500 DO ji = 2, jpim1 ! INNER domain 1501 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 1502 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 1503 END DO 1504 END DO 1505 ENDIF 1506 ! 1507 ! !== BOTTOM stress contribution from baroclinic velocities ==! 1508 ! 1509 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities 1510 1511 DO jj = 2, jpjm1 1512 DO ji = 2, jpim1 ! INNER domain 1513 ikbu = mbku(ji,jj) 1514 ikbv = mbkv(ji,jj) 1515 zu_i(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) 1516 zv_i(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj) 1517 END DO 1518 END DO 1519 ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities 1520 1521 DO jj = 2, jpjm1 1522 DO ji = 2, jpim1 ! INNER domain 1523 ikbu = mbku(ji,jj) 1524 ikbv = mbkv(ji,jj) 1525 zu_i(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) 1526 zv_i(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj) 1527 END DO 1528 END DO 1529 ENDIF 1530 ! 1531 IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please ! 1532 zztmp = -1._wp / rdtbt 1533 DO jj = 2, jpjm1 1534 DO ji = 2, jpim1 ! INNER domain 1535 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & 1536 & r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) 1537 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( & 1538 & r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) 1539 END DO 1540 END DO 1541 ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation) 1542 1543 DO jj = 2, jpjm1 1544 DO ji = 2, jpim1 ! INNER domain 1545 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) 1546 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) 1547 END DO 1548 END DO 1549 END IF 1550 ! 1551 ! !== TOP stress contribution from baroclinic velocities ==! (no W/D case) 1552 ! 1553 IF( ln_isfcav ) THEN 1554 ! 1555 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity 1556 1557 DO jj = 2, jpjm1 1558 DO ji = 2, jpim1 ! INNER domain 1559 iktu = miku(ji,jj) 1560 iktv = mikv(ji,jj) 1561 zu_i(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) 1562 zv_i(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 1563 END DO 1564 END DO 1565 ELSE ! CENTRED integration: use BEFORE top baroclinic velocity 1566 1567 DO jj = 2, jpjm1 1568 DO ji = 2, jpim1 ! INNER domain 1569 iktu = miku(ji,jj) 1570 iktv = mikv(ji,jj) 1571 zu_i(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) 1572 zv_i(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 1573 END DO 1574 END DO 1575 ENDIF 1576 ! 1577 ! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation) 1578 1579 DO jj = 2, jpjm1 1580 DO ji = 2, jpim1 ! INNER domain 1581 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) 1582 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) 1583 END DO 1584 END DO 1585 ! 1586 ENDIF 1587 ! 1588 END SUBROUTINE dyn_drg_init 1589 1590 SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in 1591 & za0, za1, za2, za3 ) ! ==> out 1592 !!---------------------------------------------------------------------- 1593 INTEGER ,INTENT(in ) :: jn ! index of sub time step 1594 LOGICAL ,INTENT(in ) :: ll_init ! 1595 REAL(wp),INTENT( out) :: za0, za1, za2, za3 ! Half-step back interpolation coefficient 1596 ! 1597 REAL(wp) :: zepsilon, zgamma ! - - 1598 !!---------------------------------------------------------------------- 1599 ! ! set Half-step back interpolation coefficient 1600 IF ( jn==1 .AND. ll_init ) THEN !* Forward-backward 1601 za0 = 1._wp 1602 za1 = 0._wp 1603 za2 = 0._wp 1604 za3 = 0._wp 1605 ELSEIF( jn==2 .AND. ll_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12 1606 za0 = 1.0833333333333_wp ! za0 = 1-gam-eps 1607 za1 =-0.1666666666666_wp ! za1 = gam 1608 za2 = 0.0833333333333_wp ! za2 = eps 1609 za3 = 0._wp 1610 ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 1611 IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion 1612 za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps 1613 za1 = 0.285_wp ! za1 = 1/2 - 2*gam - 3*eps 1614 za2 = 0.088_wp ! za2 = gam 1615 za3 = 0.013_wp ! za3 = eps 1616 ELSE ! no time diffusion 1617 zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha 1618 zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha 1619 za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon 1620 za1 = 1._wp - za0 - zgamma - zepsilon 1621 za2 = zgamma 1622 za3 = zepsilon 1623 ENDIF 1624 ENDIF 1625 END SUBROUTINE ts_bck_interp 1626 1627 1599 1628 !!====================================================================== 1600 1629 END MODULE dynspg_ts
Note: See TracChangeset
for help on using the changeset viewer.
