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