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