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