- Timestamp:
- 2019-12-05T12:06:36+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/DYN/dynspg_ts.F90
r10860 r12065 63 63 USE diatmb ! Top,middle,bottom output 64 64 65 USE iom ! to remove 66 65 67 IMPLICIT NONE 66 68 PRIVATE … … 103 105 ! 104 106 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) ) 105 !106 107 IF( ln_dynvor_een .OR. ln_dynvor_eeT ) & 107 & ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 108 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) 108 & ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , ftsw(jpi,jpj) , ftse(jpi,jpj), STAT=ierr(2) ) 109 109 ! 110 110 ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) ) … … 148 148 LOGICAL :: ll_fw_start ! =T : forward integration 149 149 LOGICAL :: ll_init ! =T : special startup of 2d equations 150 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D 151 INTEGER :: ikbu, iktu, noffset ! local integers 152 INTEGER :: ikbv, iktv ! - - 153 REAL(wp) :: r1_2dt_b, z2dt_bf ! local scalars 154 REAL(wp) :: zx1, zx2, zu_spg, zhura, z1_hu ! - - 155 REAL(wp) :: zy1, zy2, zv_spg, zhvra, z1_hv ! - - 150 INTEGER :: noffset ! local integers : time offset for bdy update 151 REAL(wp) :: r1_2dt_b, z1_hu, z1_hv ! local scalars 156 152 REAL(wp) :: za0, za1, za2, za3 ! - - 157 REAL(wp) :: zmdi, zztmp , z1_ht ! - - 158 REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf 159 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc 160 REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv 161 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e, zhtp2_e 162 REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e 153 REAL(wp) :: zmdi, zztmp, zldg ! - - 154 REAL(wp) :: zhu_bck, zhv_bck, zhdiv ! - - 155 REAL(wp) :: zun_save, zvn_save ! - - 156 REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zu_frc, zu_spg, zssh_frc 157 REAL(wp), DIMENSION(jpi,jpj) :: zv_trd, zv_frc, zv_spg 158 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhtp2_e 159 REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zsshp2_e 163 160 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points 161 REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV ! fluxes 164 162 ! 165 163 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 ) 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 ) 721 426 ! Update tide potential at the beginning of current time substep 722 427 IF( ln_tide_pot .AND. ln_tide ) THEN … … 725 430 END IF 726 431 ! 727 ! Set extrapolation coefficients for predictor step: 432 ! !== extrapolation at mid-step ==! (jn+1/2) 433 ! 434 ! !* Set extrapolation coefficients for predictor step: 728 435 IF ((jn<3).AND.ll_init) THEN ! Forward 729 436 za1 = 1._wp … … 735 442 za3 = 0.281105_wp ! za3 = bet 736 443 ENDIF 737 738 ! Extrapolate barotropic velocities at step jit+0.5: 444 ! 445 ! !* Extrapolate barotropic velocities at mid-step (jn+1/2) 446 !-- m+1/2 m m-1 m-2 --! 447 !-- u = (3/2+beta) u -(1/2+2beta) u + beta u --! 448 !-------------------------------------------------------------------------! 739 449 ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 740 450 va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) … … 743 453 ! ! ------------------ 744 454 ! Extrapolate Sea Level at step jit+0.5: 455 !-- m+1/2 m m-1 m-2 --! 456 !-- ssh = (3/2+beta) ssh -(1/2+2beta) ssh + beta ssh --! 457 !--------------------------------------------------------------------------------! 745 458 zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 746 459 747 ! set wetting & drying mask at tracer points for this barotropic sub-step 748 IF ( ln_wd_dl ) THEN 749 ! 750 IF ( ln_wd_dl_rmp ) THEN 751 DO jj = 1, jpj 752 DO ji = 1, jpi ! vector opt. 753 IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 754 ! IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 755 ztwdmask(ji,jj) = 1._wp 756 ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 757 ztwdmask(ji,jj) = (tanh(50._wp*( ( zsshp2_e(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1)) ) 758 ELSE 759 ztwdmask(ji,jj) = 0._wp 760 END IF 761 END DO 762 END DO 763 ELSE 764 DO jj = 1, jpj 765 DO ji = 1, jpi ! vector opt. 766 IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 767 ztwdmask(ji,jj) = 1._wp 768 ELSE 769 ztwdmask(ji,jj) = 0._wp 770 ENDIF 771 END DO 772 END DO 773 ENDIF 774 ! 775 ENDIF 460 ! set wetting & drying mask at tracer points for this barotropic mid-step 461 IF( ln_wd_dl ) CALL wad_tmsk( zsshp2_e, ztwdmask ) 776 462 ! 777 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 778 DO ji = 2, fs_jpim1 ! Vector opt. 779 zwx(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 780 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 781 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 782 zwy(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 783 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 784 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 785 END DO 786 END DO 787 CALL lbc_lnk_multi( 'dynspg_ts', zwx, 'U', 1._wp, zwy, 'V', 1._wp ) 463 ! ! ocean t-depth at mid-step 464 zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:) 788 465 ! 789 zhup2_e(:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 790 zhvp2_e(:,:) = hv_0(:,:) + zwy(:,:) 791 zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:) 792 ELSE 793 zhup2_e(:,:) = hu_n(:,:) 794 zhvp2_e(:,:) = hv_n(:,:) 795 zhtp2_e(:,:) = ht_n(:,:) 796 ENDIF 797 ! !* after ssh 798 ! ! ----------- 799 ! 800 ! Enforce volume conservation at open boundaries: 466 ! ! ocean u- and v-depth at mid-step (separate DO-loops remove the need of a lbc_lnk) 467 DO jj = 1, jpj 468 DO ji = 1, jpim1 ! not jpi-column 469 zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) & 470 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 471 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 472 END DO 473 END DO 474 DO jj = 1, jpjm1 ! not jpj-row 475 DO ji = 1, jpi 476 zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) & 477 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 478 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 479 END DO 480 END DO 481 ! 482 ENDIF 483 ! 484 ! !== after SSH ==! (jn+1) 485 ! 486 ! ! update (ua_e,va_e) to enforce volume conservation at open boundaries 487 ! ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d 801 488 IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e ) 802 489 ! 803 zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5 804 zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 490 ! ! resulting flux at mid-step (not over the full domain) 491 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 492 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 805 493 ! 806 494 #if defined key_agrif … … 809 497 IF((nbondi == -1).OR.(nbondi == 2)) THEN 810 498 DO jj = 1, jpj 811 z wx(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj)812 z wy(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj)499 zhU(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj) 500 zhV(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj) 813 501 END DO 814 502 ENDIF 815 503 IF((nbondi == 1).OR.(nbondi == 2)) THEN 816 504 DO jj=1,jpj 817 z wx(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj)818 z wy(nlci-nbghostcells :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells :nlci-1,jj)505 zhU(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj) 506 zhV(nlci-nbghostcells :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells :nlci-1,jj) 819 507 END DO 820 508 ENDIF 821 509 IF((nbondj == -1).OR.(nbondj == 2)) THEN 822 510 DO ji=1,jpi 823 z wy(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1)824 z wx(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1)511 zhV(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1) 512 zhU(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1) 825 513 END DO 826 514 ENDIF 827 515 IF((nbondj == 1).OR.(nbondj == 2)) THEN 828 516 DO ji=1,jpi 829 z wy(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2)830 z wx(ji,nlcj-nbghostcells :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells :nlcj-1)517 zhV(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2) 518 zhU(ji,nlcj-nbghostcells :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells :nlcj-1) 831 519 END DO 832 520 ENDIF 833 521 ENDIF 834 522 #endif 835 IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 836 837 IF ( ln_wd_dl ) THEN 838 ! 839 ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells 840 ! 841 DO jj = 1, jpjm1 842 DO ji = 1, jpim1 843 IF ( zwx(ji,jj) > 0.0 ) THEN 844 zuwdmask(ji, jj) = ztwdmask(ji ,jj) 845 ELSE 846 zuwdmask(ji, jj) = ztwdmask(ji+1,jj) 847 END IF 848 zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj) 849 un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj) 850 851 IF ( zwy(ji,jj) > 0.0 ) THEN 852 zvwdmask(ji, jj) = ztwdmask(ji, jj ) 853 ELSE 854 zvwdmask(ji, jj) = ztwdmask(ji, jj+1) 855 END IF 856 zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj) 857 vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj) 858 END DO 859 END DO 523 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 524 525 IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where 526 ! ! the direction of the flow is from dry cells 527 CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) ! not jpi colomn for U, not jpj row for V 860 528 ! 861 529 ENDIF 862 863 ! Sum over sub-time-steps to compute advective velocities 864 za2 = wgtbtp2(jn) 865 un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 866 vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 867 868 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True) 530 ! 531 ! 532 ! Compute Sea Level at step jit+1 533 !-- m+1 m m+1/2 --! 534 !-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --! 535 !-------------------------------------------------------------------------! 536 DO jj = 2, jpjm1 ! INNER domain 537 DO ji = 2, jpim1 538 zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj) 539 ssha_e(ji,jj) = ( sshn_e(ji,jj) - rdtbt * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj) 540 END DO 541 END DO 542 ! 543 CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp ) 544 ! 545 ! ! Sum over sub-time-steps to compute advective velocities 546 za2 = wgtbtp2(jn) ! zhU, zhV hold fluxes extrapolated at jn+0.5 547 un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:) 548 vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:) 549 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True) 869 550 IF ( ln_wd_dl_bc ) THEN 870 zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:) 871 zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:) 872 END IF 873 874 ! Set next sea level: 875 DO jj = 2, jpjm1 876 DO ji = fs_2, fs_jpim1 ! vector opt. 877 zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) & 878 & + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj) 879 END DO 880 END DO 881 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 882 883 CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T', 1._wp ) 884 551 zuwdav2(1:jpim1,1:jpj ) = zuwdav2(1:jpim1,1:jpj ) + za2 * zuwdmask(1:jpim1,1:jpj ) ! not jpi-column 552 zvwdav2(1:jpi ,1:jpjm1) = zvwdav2(1:jpi ,1:jpjm1) + za2 * zvwdmask(1:jpi ,1:jpjm1) ! not jpj-row 553 END IF 554 ! 885 555 ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) 886 556 IF( ln_bdy ) CALL bdy_ssh( ssha_e ) … … 891 561 ! Sea Surface Height at u-,v-points (vvl case only) 892 562 IF( .NOT.ln_linssh ) THEN 893 DO jj = 2, jpjm1 563 DO jj = 2, jpjm1 ! INNER domain, will be extended to whole domain later 894 564 DO ji = 2, jpim1 ! NO Vector Opt. 895 565 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & … … 901 571 END DO 902 572 END DO 903 CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp )904 573 ENDIF 905 ! 906 ! Half-step back interpolation of SSH for surface pressure computation: 907 !---------------------------------------------------------------------- 908 IF ((jn==1).AND.ll_init) THEN 909 za0=1._wp ! Forward-backward 910 za1=0._wp 911 za2=0._wp 912 za3=0._wp 913 ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12 914 za0= 1.0833333333333_wp ! za0 = 1-gam-eps 915 za1=-0.1666666666666_wp ! za1 = gam 916 za2= 0.0833333333333_wp ! za2 = eps 917 za3= 0._wp 918 ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 919 IF (rn_bt_alpha==0._wp) THEN 920 za0=0.614_wp ! za0 = 1/2 + gam + 2*eps 921 za1=0.285_wp ! za1 = 1/2 - 2*gam - 3*eps 922 za2=0.088_wp ! za2 = gam 923 za3=0.013_wp ! za3 = eps 924 ELSE 925 zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha 926 zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha 927 za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon 928 za1 = 1._wp - za0 - zgamma - zepsilon 929 za2 = zgamma 930 za3 = zepsilon 931 ENDIF 932 ENDIF 933 ! 574 ! 575 ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2 576 !-- m+1/2 m+1 m m-1 m-2 --! 577 !-- ssh' = za0 * ssh + za1 * ssh + za2 * ssh + za3 * ssh --! 578 !------------------------------------------------------------------------------------------! 579 CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 ) ! coeficients of the interpolation 934 580 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 935 581 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 936 937 IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters 938 DO jj = 2, jpjm1 939 DO ji = 2, jpim1 940 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 941 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 942 & MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) & 943 & > rn_wdmin1 + rn_wdmin2 944 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji+1,jj)) > 1.E-12 ).AND.( & 945 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 946 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 947 948 IF(ll_tmp1) THEN 949 zcpx(ji,jj) = 1.0_wp 950 ELSE IF(ll_tmp2) THEN 951 ! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here 952 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 953 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) ) 954 ELSE 955 zcpx(ji,jj) = 0._wp 956 ENDIF 957 ! 958 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 959 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 960 & MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) & 961 & > rn_wdmin1 + rn_wdmin2 962 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji,jj+1)) > 1.E-12 ).AND.( & 963 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 964 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 965 966 IF(ll_tmp1) THEN 967 zcpy(ji,jj) = 1.0_wp 968 ELSEIF(ll_tmp2) THEN 969 ! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here 970 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 971 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) ) 972 ELSE 973 zcpy(ji,jj) = 0._wp 974 ENDIF 975 END DO 976 END DO 977 ENDIF 978 ! 979 ! Compute associated depths at U and V points: 980 IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form 981 ! 982 DO jj = 2, jpjm1 983 DO ji = 2, jpim1 984 zx1 = r1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) & 985 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) & 986 & + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) ) 987 zy1 = r1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) & 988 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) & 989 & + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) ) 990 zhust_e(ji,jj) = hu_0(ji,jj) + zx1 991 zhvst_e(ji,jj) = hv_0(ji,jj) + zy1 992 END DO 993 END DO 994 ! 582 ! 583 ! ! Surface pressure gradient 584 zldg = ( 1._wp - rn_scal_load ) * grav ! local factor 585 DO jj = 2, jpjm1 586 DO ji = 2, jpim1 587 zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 588 zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 589 END DO 590 END DO 591 IF( ln_wd_il ) THEN ! W/D : gravity filters applied on pressure gradient 592 CALL wad_spg( zsshp2_e, zcpx, zcpy ) ! Calculating W/D gravity filters 593 zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1) 594 zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1) 995 595 ENDIF 996 596 ! … … 998 598 ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated 999 599 ! at each time step. We however keep them constant here for optimization. 1000 ! Recall that zwx and zwy arrays hold fluxes at this stage: 1001 ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5 1002 ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 1003 ! 1004 SELECT CASE( nvor_scheme ) 1005 CASE( np_ENT ) ! energy conserving scheme (t-point) 1006 DO jj = 2, jpjm1 1007 DO ji = 2, jpim1 ! vector opt. 1008 1009 z1_hu = ssumask(ji,jj) / ( zhup2_e(ji,jj) + 1._wp - ssumask(ji,jj) ) 1010 z1_hv = ssvmask(ji,jj) / ( zhvp2_e(ji,jj) + 1._wp - ssvmask(ji,jj) ) 1011 1012 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu & 1013 & * ( 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) ) & 1014 & + e1e2t(ji ,jj)*zhtp2_e(ji ,jj)*ff_t(ji ,jj) * ( va_e(ji ,jj) + va_e(ji ,jj-1) ) ) 1015 ! 1016 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1017 & * ( 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) ) & 1018 & + e1e2t(ji,jj )*zhtp2_e(ji,jj )*ff_t(ji,jj ) * ( ua_e(ji,jj ) + ua_e(ji-1,jj ) ) ) 1019 END DO 1020 END DO 1021 ! 1022 CASE( np_ENE, np_MIX ) ! energy conserving scheme (f-point) 1023 DO jj = 2, jpjm1 1024 DO ji = fs_2, fs_jpim1 ! vector opt. 1025 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 1026 zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 1027 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 1028 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 1029 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 1030 zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 1031 END DO 1032 END DO 1033 ! 1034 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 1035 DO jj = 2, jpjm1 1036 DO ji = fs_2, fs_jpim1 ! vector opt. 1037 zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 1038 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 1039 zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 1040 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 1041 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 1042 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 1043 END DO 1044 END DO 1045 ! 1046 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1047 DO jj = 2, jpjm1 1048 DO ji = fs_2, fs_jpim1 ! vector opt. 1049 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 1050 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 1051 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 1052 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 1053 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 1054 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 1055 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 1056 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 1057 END DO 1058 END DO 1059 ! 1060 END SELECT 600 ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated) 601 CALL dyn_cor_2d( zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd ) 1061 602 ! 1062 603 ! Add tidal astronomical forcing if defined … … 1064 605 DO jj = 2, jpjm1 1065 606 DO ji = fs_2, fs_jpim1 ! vector opt. 1066 zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 1067 zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 1068 zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg 1069 zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg 607 zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 608 zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 1070 609 END DO 1071 610 END DO … … 1081 620 END DO 1082 621 END DO 1083 ENDIF 1084 ! 1085 ! Surface pressure trend: 1086 IF( ln_wd_il ) THEN 1087 DO jj = 2, jpjm1 1088 DO ji = 2, jpim1 1089 ! Add surface pressure gradient 1090 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 1091 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 1092 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj) 1093 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj) 1094 END DO 1095 END DO 1096 ELSE 1097 DO jj = 2, jpjm1 1098 DO ji = fs_2, fs_jpim1 ! vector opt. 1099 ! Add surface pressure gradient 1100 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 1101 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 1102 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg 1103 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg 1104 END DO 1105 END DO 1106 END IF 1107 622 ENDIF 1108 623 ! 1109 624 ! Set next velocities: 625 ! Compute barotropic speeds at step jit+1 (h : total height of the water colomn) 626 !-- VECTOR FORM 627 !-- m+1 m / m+1/2 \ --! 628 !-- u = u + delta_t' * \ (1-r)*g * grad_x( ssh') - f * k vect u + frc / --! 629 !-- --! 630 !-- FLUX FORM --! 631 !-- m+1 __1__ / m m / m+1/2 m+1/2 m+1/2 n \ \ --! 632 !-- u = m+1 | h * u + delta_t' * \ h * (1-r)*g * grad_x( ssh') - h * f * k vect u + h * frc / | --! 633 !-- h \ / --! 634 !------------------------------------------------------------------------------------------------------------------------! 1110 635 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 1111 636 DO jj = 2, jpjm1 1112 637 DO ji = fs_2, fs_jpim1 ! vector opt. 1113 638 ua_e(ji,jj) = ( un_e(ji,jj) & 1114 & + rdtbt * ( zwx(ji,jj) &639 & + rdtbt * ( zu_spg(ji,jj) & 1115 640 & + zu_trd(ji,jj) & 1116 641 & + zu_frc(ji,jj) ) & … … 1118 643 1119 644 va_e(ji,jj) = ( vn_e(ji,jj) & 1120 & + rdtbt * ( zwy(ji,jj) &645 & + rdtbt * ( zv_spg(ji,jj) & 1121 646 & + zv_trd(ji,jj) & 1122 647 & + zv_frc(ji,jj) ) & 1123 648 & ) * ssvmask(ji,jj) 1124 1125 649 END DO 1126 650 END DO … … 1128 652 ELSE !* Flux form 1129 653 DO jj = 2, jpjm1 1130 DO ji = fs_2, fs_jpim1 ! vector opt. 1131 1132 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 1133 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 1134 1135 zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj)) 1136 zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj)) 1137 1138 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & 1139 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 1140 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 1141 & + hu_n(ji,jj) * zu_frc(ji,jj) ) & 1142 & ) * zhura 1143 1144 va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) & 1145 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 1146 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & 1147 & + hv_n(ji,jj) * zv_frc(ji,jj) ) & 1148 & ) * zhvra 654 DO ji = 2, jpim1 655 ! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2 656 ! ! backward interpolated depth used in spg terms at jn+1/2 657 zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 658 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 659 zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 660 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 661 ! ! inverse depth at jn+1 662 z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) 663 z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) ) 664 ! 665 ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) & 666 & + rdtbt * ( zhu_bck * zu_spg (ji,jj) & ! 667 & + zhup2_e(ji,jj) * zu_trd (ji,jj) & ! 668 & + hu_n (ji,jj) * zu_frc (ji,jj) ) ) * z1_hu 669 ! 670 va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) & 671 & + rdtbt * ( zhv_bck * zv_spg (ji,jj) & ! 672 & + zhvp2_e(ji,jj) * zv_trd (ji,jj) & ! 673 & + hv_n (ji,jj) * zv_frc (ji,jj) ) ) * z1_hv 1149 674 END DO 1150 675 END DO … … 1159 684 END DO 1160 685 ENDIF 1161 1162 1163 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 1164 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 1165 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 1166 hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 1167 hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) 1168 ! 1169 ENDIF 1170 ! !* domain lateral boundary 1171 CALL lbc_lnk_multi( 'dynspg_ts', ua_e, 'U', -1._wp, va_e , 'V', -1._wp ) 686 687 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 688 hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1) 689 hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) ) 690 hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1) 691 hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) ) 692 CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp & 693 & , hu_e , 'U', -1._wp, hv_e , 'V', -1._wp & 694 & , hur_e, 'U', -1._wp, hvr_e, 'V', -1._wp ) 695 ELSE 696 CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp ) 697 ENDIF 698 ! 1172 699 ! 1173 700 ! ! open boundaries … … 1217 744 ! Set advection velocity correction: 1218 745 IF (ln_bt_fw) THEN 1219 zwx(:,:) = un_adv(:,:)1220 zwy(:,:) = vn_adv(:,:)1221 746 IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN 1222 un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) ) 1223 vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) ) 1224 ! 1225 ! Update corrective fluxes for next time step: 1226 un_bf(:,:) = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:)) 1227 vn_bf(:,:) = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:)) 747 DO jj = 1, jpj 748 DO ji = 1, jpi 749 zun_save = un_adv(ji,jj) 750 zvn_save = vn_adv(ji,jj) 751 ! ! apply the previously computed correction 752 un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) ) 753 vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) ) 754 ! ! Update corrective fluxes for next time step 755 un_bf(ji,jj) = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) ) 756 vn_bf(ji,jj) = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) ) 757 ! ! Save integrated transport for next computation 758 ub2_b(ji,jj) = zun_save 759 vb2_b(ji,jj) = zvn_save 760 END DO 761 END DO 1228 762 ELSE 1229 un_bf(:,:) = 0._wp 1230 vn_bf(:,:) = 0._wp 1231 END IF 1232 ! Save integrated transport for next computation 1233 ub2_b(:,:) = zwx(:,:) 1234 vb2_b(:,:) = zwy(:,:) 763 un_bf(:,:) = 0._wp ! corrective fluxes for next time step set to zero 764 vn_bf(:,:) = 0._wp 765 ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation 766 vb2_b(:,:) = vn_adv(:,:) 767 END IF 1235 768 ENDIF 1236 769 … … 1274 807 1275 808 IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN 809 ! need to set lbc here because not done prior time averaging 810 CALL lbc_lnk_multi( 'dynspg_ts', zuwdav2, 'U', 1._wp, zvwdav2, 'V', 1._wp) 1276 811 DO jk = 1, jpkm1 1277 812 un(:,:,jk) = ( un_adv(:,:)*r1_hu_n(:,:) & … … 1477 1012 REAL(wp) :: zxr2, zyr2, zcmax ! local scalar 1478 1013 REAL(wp), DIMENSION(jpi,jpj) :: zcu 1014 INTEGER :: inum 1479 1015 !!---------------------------------------------------------------------- 1480 1016 ! … … 1583 1119 END SUBROUTINE dyn_spg_ts_init 1584 1120 1121 1122 SUBROUTINE dyn_cor_2d_init 1123 !!--------------------------------------------------------------------- 1124 !! *** ROUTINE dyn_cor_2d_init *** 1125 !! 1126 !! ** Purpose : Set time splitting options 1127 !! Set arrays to remove/compute coriolis trend. 1128 !! Do it once during initialization if volume is fixed, else at each long time step. 1129 !! Note that these arrays are also used during barotropic loop. These are however frozen 1130 !! although they should be updated in the variable volume case. Not a big approximation. 1131 !! To remove this approximation, copy lines below inside barotropic loop 1132 !! and update depths at T-F points (ht and zhf resp.) at each barotropic time step 1133 !! 1134 !! Compute zwz = f / ( height of the water colomn ) 1135 !!---------------------------------------------------------------------- 1136 INTEGER :: ji ,jj, jk ! dummy loop indices 1137 REAL(wp) :: z1_ht 1138 REAL(wp), DIMENSION(jpi,jpj) :: zhf 1139 !!---------------------------------------------------------------------- 1140 ! 1141 SELECT CASE( nvor_scheme ) 1142 CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme) 1143 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point 1144 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 1145 DO jj = 1, jpjm1 1146 DO ji = 1, jpim1 1147 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 1148 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp 1149 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1150 END DO 1151 END DO 1152 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 1153 DO jj = 1, jpjm1 1154 DO ji = 1, jpim1 1155 zwz(ji,jj) = ( ht_n (ji ,jj+1) + ht_n (ji+1,jj+1) & 1156 & + ht_n (ji ,jj ) + ht_n (ji+1,jj ) ) & 1157 & / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) & 1158 & + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) ) 1159 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1160 END DO 1161 END DO 1162 END SELECT 1163 CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp ) 1164 ! 1165 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1166 DO jj = 2, jpj 1167 DO ji = 2, jpi 1168 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 1169 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 1170 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 1171 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 1172 END DO 1173 END DO 1174 ! 1175 CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme) 1176 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1177 DO jj = 2, jpj 1178 DO ji = 2, jpi 1179 z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) ) 1180 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht 1181 ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht 1182 ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht 1183 ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht 1184 END DO 1185 END DO 1186 ! 1187 CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT ! 1188 ! 1189 zwz(:,:) = 0._wp 1190 zhf(:,:) = 0._wp 1191 1192 !!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed 1193 !!gm A priori a better value should be something like : 1194 !!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1) 1195 !!gm divided by the sum of the corresponding mask 1196 !!gm 1197 !! 1198 IF( .NOT.ln_sco ) THEN 1199 1200 !!gm agree the JC comment : this should be done in a much clear way 1201 1202 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 1203 ! Set it to zero for the time being 1204 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level 1205 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth 1206 ! ENDIF 1207 ! zhf(:,:) = gdepw_0(:,:,jk+1) 1208 ! 1209 ELSE 1210 ! 1211 !zhf(:,:) = hbatf(:,:) 1212 DO jj = 1, jpjm1 1213 DO ji = 1, jpim1 1214 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) & 1215 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) & 1216 & / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) & 1217 & + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp ) 1218 END DO 1219 END DO 1220 ENDIF 1221 ! 1222 DO jj = 1, jpjm1 1223 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 1224 END DO 1225 ! 1226 DO jk = 1, jpkm1 1227 DO jj = 1, jpjm1 1228 zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 1229 END DO 1230 END DO 1231 CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp ) 1232 ! JC: TBC. hf should be greater than 0 1233 DO jj = 1, jpj 1234 DO ji = 1, jpi 1235 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) 1236 END DO 1237 END DO 1238 zwz(:,:) = ff_f(:,:) * zwz(:,:) 1239 END SELECT 1240 1241 END SUBROUTINE dyn_cor_2d_init 1242 1243 1244 1245 SUBROUTINE dyn_cor_2d( hu_n, hv_n, un_b, vn_b, zhU, zhV, zu_trd, zv_trd ) 1246 !!--------------------------------------------------------------------- 1247 !! *** ROUTINE dyn_cor_2d *** 1248 !! 1249 !! ** Purpose : Compute u and v coriolis trends 1250 !!---------------------------------------------------------------------- 1251 INTEGER :: ji ,jj ! dummy loop indices 1252 REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - - 1253 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: hu_n, hv_n, un_b, vn_b, zhU, zhV 1254 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd 1255 !!---------------------------------------------------------------------- 1256 SELECT CASE( nvor_scheme ) 1257 CASE( np_ENT ) ! enstrophy conserving scheme (f-point) 1258 DO jj = 2, jpjm1 1259 DO ji = 2, jpim1 1260 z1_hu = ssumask(ji,jj) / ( hu_n(ji,jj) + 1._wp - ssumask(ji,jj) ) 1261 z1_hv = ssvmask(ji,jj) / ( hv_n(ji,jj) + 1._wp - ssvmask(ji,jj) ) 1262 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu & 1263 & * ( 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) ) & 1264 & + e1e2t(ji ,jj)*ht_n(ji ,jj)*ff_t(ji ,jj) * ( vn_b(ji ,jj) + vn_b(ji ,jj-1) ) ) 1265 ! 1266 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1267 & * ( 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) ) & 1268 & + e1e2t(ji,jj )*ht_n(ji,jj )*ff_t(ji,jj ) * ( un_b(ji,jj ) + un_b(ji-1,jj ) ) ) 1269 END DO 1270 END DO 1271 ! 1272 CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX 1273 DO jj = 2, jpjm1 1274 DO ji = fs_2, fs_jpim1 ! vector opt. 1275 zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj) 1276 zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) 1277 zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj) 1278 zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj) 1279 ! energy conserving formulation for planetary vorticity term 1280 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 1281 zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 1282 END DO 1283 END DO 1284 ! 1285 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 1286 DO jj = 2, jpjm1 1287 DO ji = fs_2, fs_jpim1 ! vector opt. 1288 zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) & 1289 & + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) 1290 zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) & 1291 & + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj) 1292 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 1293 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 1294 END DO 1295 END DO 1296 ! 1297 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1298 DO jj = 2, jpjm1 1299 DO ji = fs_2, fs_jpim1 ! vector opt. 1300 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) & 1301 & + ftnw(ji+1,jj) * zhV(ji+1,jj ) & 1302 & + ftse(ji,jj ) * zhV(ji ,jj-1) & 1303 & + ftsw(ji+1,jj) * zhV(ji+1,jj-1) ) 1304 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) & 1305 & + ftse(ji,jj+1) * zhU(ji ,jj+1) & 1306 & + ftnw(ji,jj ) * zhU(ji-1,jj ) & 1307 & + ftne(ji,jj ) * zhU(ji ,jj ) ) 1308 END DO 1309 END DO 1310 ! 1311 END SELECT 1312 ! 1313 END SUBROUTINE dyn_cor_2D 1314 1315 1316 SUBROUTINE wad_tmsk( pssh, ptmsk ) 1317 !!---------------------------------------------------------------------- 1318 !! *** ROUTINE wad_lmt *** 1319 !! 1320 !! ** Purpose : set wetting & drying mask at tracer points 1321 !! for the current barotropic sub-step 1322 !! 1323 !! ** Method : ??? 1324 !! 1325 !! ** Action : ptmsk : wetting & drying t-mask 1326 !!---------------------------------------------------------------------- 1327 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh ! 1328 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: ptmsk ! 1329 ! 1330 INTEGER :: ji, jj ! dummy loop indices 1331 !!---------------------------------------------------------------------- 1332 ! 1333 IF( ln_wd_dl_rmp ) THEN 1334 DO jj = 1, jpj 1335 DO ji = 1, jpi 1336 IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 1337 ! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 1338 ptmsk(ji,jj) = 1._wp 1339 ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 1340 ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1) ) 1341 ELSE 1342 ptmsk(ji,jj) = 0._wp 1343 ENDIF 1344 END DO 1345 END DO 1346 ELSE 1347 DO jj = 1, jpj 1348 DO ji = 1, jpi 1349 IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp 1350 ELSE ; ptmsk(ji,jj) = 0._wp 1351 ENDIF 1352 END DO 1353 END DO 1354 ENDIF 1355 ! 1356 END SUBROUTINE wad_tmsk 1357 1358 1359 SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk ) 1360 !!---------------------------------------------------------------------- 1361 !! *** ROUTINE wad_lmt *** 1362 !! 1363 !! ** Purpose : set wetting & drying mask at tracer points 1364 !! for the current barotropic sub-step 1365 !! 1366 !! ** Method : ??? 1367 !! 1368 !! ** Action : ptmsk : wetting & drying t-mask 1369 !!---------------------------------------------------------------------- 1370 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pTmsk ! W & D t-mask 1371 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports 1372 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask 1373 ! 1374 INTEGER :: ji, jj ! dummy loop indices 1375 !!---------------------------------------------------------------------- 1376 ! 1377 DO jj = 1, jpj 1378 DO ji = 1, jpim1 ! not jpi-column 1379 IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj) 1380 ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj) 1381 ENDIF 1382 phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj) 1383 pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj) 1384 END DO 1385 END DO 1386 ! 1387 DO jj = 1, jpjm1 ! not jpj-row 1388 DO ji = 1, jpi 1389 IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj ) 1390 ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1) 1391 ENDIF 1392 phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 1393 pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj) 1394 END DO 1395 END DO 1396 ! 1397 END SUBROUTINE wad_Umsk 1398 1399 1400 SUBROUTINE wad_spg( sshn, zcpx, zcpy ) 1401 !!--------------------------------------------------------------------- 1402 !! *** ROUTINE wad_sp *** 1403 !! 1404 !! ** Purpose : 1405 !!---------------------------------------------------------------------- 1406 INTEGER :: ji ,jj ! dummy loop indices 1407 LOGICAL :: ll_tmp1, ll_tmp2 1408 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: sshn 1409 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy 1410 !!---------------------------------------------------------------------- 1411 DO jj = 2, jpjm1 1412 DO ji = 2, jpim1 1413 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1414 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 1415 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 1416 & > rn_wdmin1 + rn_wdmin2 1417 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 1418 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1419 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 1420 IF(ll_tmp1) THEN 1421 zcpx(ji,jj) = 1.0_wp 1422 ELSEIF(ll_tmp2) THEN 1423 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 1424 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 1425 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 1426 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 1427 ELSE 1428 zcpx(ji,jj) = 0._wp 1429 ENDIF 1430 ! 1431 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1432 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 1433 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 1434 & > rn_wdmin1 + rn_wdmin2 1435 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( & 1436 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1437 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1438 1439 IF(ll_tmp1) THEN 1440 zcpy(ji,jj) = 1.0_wp 1441 ELSE IF(ll_tmp2) THEN 1442 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 1443 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 1444 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 1445 zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) ) 1446 ELSE 1447 zcpy(ji,jj) = 0._wp 1448 ENDIF 1449 END DO 1450 END DO 1451 1452 END SUBROUTINE wad_spg 1453 1454 1455 1456 SUBROUTINE dyn_drg_init( pu_RHSi, pv_RHSi, pCdU_u, pCdU_v ) 1457 !!---------------------------------------------------------------------- 1458 !! *** ROUTINE dyn_drg_init *** 1459 !! 1460 !! ** Purpose : - add the baroclinic top/bottom drag contribution to 1461 !! the baroclinic part of the barotropic RHS 1462 !! - compute the barotropic drag coefficients 1463 !! 1464 !! ** Method : computation done over the INNER domain only 1465 !!---------------------------------------------------------------------- 1466 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS 1467 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients 1468 ! 1469 INTEGER :: ji, jj ! dummy loop indices 1470 INTEGER :: ikbu, ikbv, iktu, iktv 1471 REAL(wp) :: zztmp 1472 REAL(wp), DIMENSION(jpi,jpj) :: zu_i, zv_i 1473 !!---------------------------------------------------------------------- 1474 ! 1475 ! !== Set the barotropic drag coef. ==! 1476 ! 1477 IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities) 1478 1479 DO jj = 2, jpjm1 1480 DO ji = 2, jpim1 ! INNER domain 1481 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 1482 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) 1483 END DO 1484 END DO 1485 ELSE ! bottom friction only 1486 DO jj = 2, jpjm1 1487 DO ji = 2, jpim1 ! INNER domain 1488 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 1489 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 1490 END DO 1491 END DO 1492 ENDIF 1493 ! 1494 ! !== BOTTOM stress contribution from baroclinic velocities ==! 1495 ! 1496 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities 1497 1498 DO jj = 2, jpjm1 1499 DO ji = 2, jpim1 ! INNER domain 1500 ikbu = mbku(ji,jj) 1501 ikbv = mbkv(ji,jj) 1502 zu_i(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) 1503 zv_i(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj) 1504 END DO 1505 END DO 1506 ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities 1507 1508 DO jj = 2, jpjm1 1509 DO ji = 2, jpim1 ! INNER domain 1510 ikbu = mbku(ji,jj) 1511 ikbv = mbkv(ji,jj) 1512 zu_i(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) 1513 zv_i(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj) 1514 END DO 1515 END DO 1516 ENDIF 1517 ! 1518 IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please ! 1519 zztmp = -1._wp / rdtbt 1520 DO jj = 2, jpjm1 1521 DO ji = 2, jpim1 ! INNER domain 1522 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & 1523 & r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) 1524 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( & 1525 & r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) 1526 END DO 1527 END DO 1528 ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation) 1529 1530 DO jj = 2, jpjm1 1531 DO ji = 2, jpim1 ! INNER domain 1532 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) 1533 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) 1534 END DO 1535 END DO 1536 END IF 1537 ! 1538 ! !== TOP stress contribution from baroclinic velocities ==! (no W/D case) 1539 ! 1540 IF( ln_isfcav ) THEN 1541 ! 1542 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity 1543 1544 DO jj = 2, jpjm1 1545 DO ji = 2, jpim1 ! INNER domain 1546 iktu = miku(ji,jj) 1547 iktv = mikv(ji,jj) 1548 zu_i(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) 1549 zv_i(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 1550 END DO 1551 END DO 1552 ELSE ! CENTRED integration: use BEFORE top baroclinic velocity 1553 1554 DO jj = 2, jpjm1 1555 DO ji = 2, jpim1 ! INNER domain 1556 iktu = miku(ji,jj) 1557 iktv = mikv(ji,jj) 1558 zu_i(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) 1559 zv_i(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 1560 END DO 1561 END DO 1562 ENDIF 1563 ! 1564 ! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation) 1565 1566 DO jj = 2, jpjm1 1567 DO ji = 2, jpim1 ! INNER domain 1568 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) 1569 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) 1570 END DO 1571 END DO 1572 ! 1573 ENDIF 1574 ! 1575 END SUBROUTINE dyn_drg_init 1576 1577 SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in 1578 & za0, za1, za2, za3 ) ! ==> out 1579 !!---------------------------------------------------------------------- 1580 INTEGER ,INTENT(in ) :: jn ! index of sub time step 1581 LOGICAL ,INTENT(in ) :: ll_init ! 1582 REAL(wp),INTENT( out) :: za0, za1, za2, za3 ! Half-step back interpolation coefficient 1583 ! 1584 REAL(wp) :: zepsilon, zgamma ! - - 1585 !!---------------------------------------------------------------------- 1586 ! ! set Half-step back interpolation coefficient 1587 IF ( jn==1 .AND. ll_init ) THEN !* Forward-backward 1588 za0 = 1._wp 1589 za1 = 0._wp 1590 za2 = 0._wp 1591 za3 = 0._wp 1592 ELSEIF( jn==2 .AND. ll_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12 1593 za0 = 1.0833333333333_wp ! za0 = 1-gam-eps 1594 za1 =-0.1666666666666_wp ! za1 = gam 1595 za2 = 0.0833333333333_wp ! za2 = eps 1596 za3 = 0._wp 1597 ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 1598 IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion 1599 za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps 1600 za1 = 0.285_wp ! za1 = 1/2 - 2*gam - 3*eps 1601 za2 = 0.088_wp ! za2 = gam 1602 za3 = 0.013_wp ! za3 = eps 1603 ELSE ! no time diffusion 1604 zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha 1605 zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha 1606 za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon 1607 za1 = 1._wp - za0 - zgamma - zepsilon 1608 za2 = zgamma 1609 za3 = zepsilon 1610 ENDIF 1611 ENDIF 1612 END SUBROUTINE ts_bck_interp 1613 1614 1585 1615 !!====================================================================== 1586 1616 END MODULE dynspg_ts
Note: See TracChangeset
for help on using the changeset viewer.
