Changeset 15403
- Timestamp:
- 2021-10-19T15:48:28+02:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90
r13681 r15403 79 79 !! 80 80 !! f1(u) = g1(v) 81 !! f2(v) = g2( v)81 !! f2(v) = g2(u) 82 82 !! 83 83 !! The right-hand side (RHS) is explicit … … 143 143 REAL(wp) :: zkt ! isotropic tensile strength for landfast ice 144 144 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV ! ice/snow mass and volume 145 REAL(wp) :: zd eltat, zds2, zdt, zdt2, zdiv, zdiv2! temporary scalars145 REAL(wp) :: zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 146 146 REAL(wp) :: zp_deltastar_f ! 147 147 REAL(wp) :: zu_cV, zv_cU ! … … 156 156 REAL(wp), DIMENSION(jpi,jpj) :: zmassU_t, zmassV_t ! Mass per unit area divided by time step 157 157 ! 158 REAL(wp), DIMENSION(jpi,jpj) :: zdelta star_t !Delta* at T-points158 REAL(wp), DIMENSION(jpi,jpj) :: zdeltat, zdeltastar_t ! Delta & Delta* at T-points 159 159 REAL(wp), DIMENSION(jpi,jpj) :: zten_i ! Tension 160 160 REAL(wp), DIMENSION(jpi,jpj) :: zp_deltastar_t ! P/delta* at T points … … 207 207 !!---------------------------------------------------------------------------------------------------------------------- 208 208 ! DEBUG put all forcing terms to zero 209 ! air-ice drag 210 utau_ice(:,:) = 0._wp 211 vtau_ice(:,:) = 0._wp 212 ! coriolis 213 ff_t(:,:) = 0._wp 214 ! ice-ocean drag 215 rn_cio = 0._wp 216 ! ssh 217 ! done line 330 !!! dont forget to act there 209 ! air-ice drag 210 !utau_ice(:,:) = 0._wp 211 !vtau_ice(:,:) = 0._wp 212 213 ! coriolis 214 !ff_t(:,:) = 0._wp 215 216 ! ice-ocean drag 217 !rn_cio = 0._wp 218 219 ! sea surface height 220 ! embedded sea ice: compute representative ice top surface 221 ! non-embedded sea ice: use ocean surface for slope calculation 222 zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) ! DEBUG 223 !zsshdyn(:,:) = 0._wp 224 218 225 ! END DEBUG 226 !!---------------------------------------------------------------------------------------------------------------------- 219 227 220 228 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_vp: VP sea-ice rheology (LSR solver)' … … 265 273 ENDIF 266 274 267 zs1_rhsu(:,:) = 0._wp; zs2_rhsu(:,:) = 0._wp; zs1_rhsv(:,:) = 0._wp; zs2_rhsv(:,:) = 0._wp 275 zs1_rhsu(:,:) = 0._wp; zs2_rhsu(:,:) = 0._wp 276 zs1_rhsv(:,:) = 0._wp; zs2_rhsv(:,:) = 0._wp 268 277 zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp; 269 278 zAV(:,:) = 0._wp; zBV(:,:) = 0._wp; zCV(:,:) = 0._wp; zDV(:,:) = 0._wp; zEV(:,:) = 0._wp; 270 zrhsu (:,:) = 0._wp; zrhsv(:,:) = 0._wp271 zf_rhsu(:,:) = 0._wp ; zf_rhsv(:,:) = 0._wp279 zrhsu (:,:) = 0._wp ; zrhsv (:,:) = 0._wp 280 zf_rhsu(:,:) = 0._wp ; zf_rhsv(:,:) = 0._wp 272 281 273 282 !------------------------------------------------------------------------------! … … 279 288 CALL ice_strength ! strength at T points 280 289 281 !--------------------------- ---282 ! -- F-mask 283 !--------------------------- ---290 !--------------------------- 291 ! -- F-mask (code from EVP) 292 !--------------------------- 284 293 ! MartinV: 285 294 ! In EVP routine, zfmask is applied on shear at F-points, in order to enforce the lateral boundary condition (no-slip, ..., free-slip) … … 289 298 DO jj = 1, jpj - 1 290 299 DO ji = 1, jpi - 1 291 zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)300 zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 292 301 END DO 293 302 END DO … … 326 335 !---------------------------------------------------------------------------------------------------------- 327 336 ! Compute all terms & factors independent of velocities, or only depending on velocities at previous time step 328 329 ! sea surface height330 ! embedded sea ice: compute representative ice top surface331 ! non-embedded sea ice: use ocean surface for slope calculation332 zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)333 zsshdyn(:,:) = 0._wp ! DEBUG CAREFUL !!!334 337 335 338 zmt(:,:) = rhos * vt_s(:,:) + rhoi * vt_i(:,:) ! Snow and ice mass at T-point … … 347 350 zm2 = zmt(ji+1,jj) 348 351 zm3 = zmt(ji,jj+1) 349 350 352 zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 351 353 zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) … … 382 384 383 385 ! MV TEST DEBUG 384 385 386 387 388 389 390 386 ! IF ( ( zmt(ji,jj) <= zmmin .OR. zmt(ji+1,jj) <= zmmin ) .AND. & 387 ! & ( at_i(ji,jj) <= zamin .OR. at_i(ji+1,jj) <= zamin ) ) THEN ; zmsk01x(ji,jj) = 0._wp 388 ! ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF 389 390 ! IF ( ( zmt(ji,jj) <= zmmin .OR. zmt(ji,jj+1) <= zmmin ) .AND. & 391 ! & ( at_i(ji,jj) <= zamin .OR. at_i(ji,jj+1) <= zamin ) ) THEN ; zmsk01y(ji,jj) = 0._wp 392 ! ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF 391 393 ! END MV TEST DEBUG 392 394 393 395 END DO 394 396 END DO 397 398 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zspgU , 'U', -1., zspgV , 'V', -1. ) 399 CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_ai, 'U', -1., ztauy_ai, 'V', -1. ) 395 400 396 401 CALL iom_put( 'zmsk00x' , zmsk00x ) ! MV DEBUG … … 431 436 ! In the outer loop, one needs to update all RHS terms 432 437 ! with explicit velocity dependencies (viscosities, coriolis, ocean stress) 433 ! as a function of uc 434 ! 438 ! as a function of "current" velocities (uc, vc) 435 439 436 440 !------------------------------------------ … … 451 455 452 456 CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! MV TEST could be un-necessary according to Gurvan 453 CALL iom_put( 'zds' , zds ) ! MV DEBUG 457 458 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zds' , zds ) ! MV DEBUG 454 459 455 460 IF( lwp ) WRITE(numout,*) ' outer loop 1a i_out : ', i_out … … 464 469 465 470 ! shear**2 at T points (doc eq. A16) 466 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) &467 & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) &468 & ) * 0.25_wp * r1_e1e2t(ji,jj)471 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & 472 & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & 473 & ) * 0.25_wp * r1_e1e2t(ji,jj) 469 474 470 475 ! divergence at T points … … 475 480 476 481 ! tension at T points 477 zdt = ( ( zu_c(ji,jj) * r1_e2u(ji,jj) - zu_c(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) &478 & - ( zv_c(ji,jj) * r1_e1v(ji,jj) - zv_c(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) &479 & ) * r1_e1e2t(ji,jj)480 zdt2 = zdt * zdt482 zdt = ( ( zu_c(ji,jj) * r1_e2u(ji,jj) - zu_c(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj) & 483 & - ( zv_c(ji,jj) * r1_e1v(ji,jj) - zv_c(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 484 & ) * r1_e1e2t(ji,jj) 485 zdt2 = zdt * zdt 481 486 482 487 ! delta at T points 483 zdeltat = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )488 zdeltat(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 484 489 485 490 ! delta* at T points (following Lemieux and Dupont, GMD 2020) 486 zdeltastar_t(ji,jj) = zdeltat+ rn_creepl491 zdeltastar_t(ji,jj) = zdeltat(ji,jj) + rn_creepl 487 492 488 493 ! P/delta at T-points … … 490 495 491 496 ! Temporary zzt and zet factors at T-points 492 zzt(ji,jj) = zp_deltastar_t(ji,jj) * r1_e1e2t(ji,jj)493 zet(ji,jj) = zzt(ji,jj) * z1_ecc2497 zzt(ji,jj) = zp_deltastar_t(ji,jj) * r1_e1e2t(ji,jj) 498 zet(ji,jj) = zzt(ji,jj) * z1_ecc2 494 499 495 500 END DO 496 501 END DO 497 498 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. ) 499 500 CALL iom_put( 'zzt' , zzt ) ! MV DEBUG 501 CALL iom_put( 'zet' , zet ) ! MV DEBUG 502 CALL iom_put( 'zp_deltastar_t', zp_deltastar_t ) ! MV DEBUG 502 503 ! zzt(:,:) = 0. ! test 504 ! zet(:,:) = 0. ! test 505 506 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. , zdeltat, 'T', 1.) 507 508 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zzt' , zzt ) ! MV DEBUG 509 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zet' , zet ) ! MV DEBUG 510 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zp_deltastar_t', zp_deltastar_t ) ! MV DEBUG 511 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zdeltat' , zdeltat ) ! MV DEBUG 503 512 504 513 IF( lwp ) WRITE(numout,*) ' outer loop 1b i_out : ', i_out … … 515 524 END DO 516 525 END DO 526 527 ! zef(:,:) = 0. ! test 517 528 518 529 CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. ) 519 CALL iom_put( 'zef' , zef ) ! MV DEBUG530 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zef' , zef ) ! MV DEBUG 520 531 IF( lwp ) WRITE(numout,*) ' outer loop 1c i_out : ', i_out 521 532 … … 551 562 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * zu_c(ji,jj ) + e2u(ji-1,jj ) * zu_c(ji-1,jj ) ) & 552 563 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * zu_c(ji,jj+1) + e2u(ji-1,jj+1) * zu_c(ji-1,jj+1) ) ) 553 564 554 565 END DO 555 566 END DO … … 560 571 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCorU, 'U', -1., zCorV, 'V', -1. ) 561 572 562 CALL iom_put( 'zCwU' , zCwU ) ! MV DEBUG563 CALL iom_put( 'zCwV' , zCwV ) ! MV DEBUG564 CALL iom_put( 'zCorU' , zCorU ) ! MV DEBUG565 CALL iom_put( 'zCorV' , zCorV ) ! MV DEBUG573 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zCwU' , zCwU ) ! MV DEBUG 574 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zCwV' , zCwV ) ! MV DEBUG 575 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zCorU' , zCorU ) ! MV DEBUG 576 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zCorV' , zCorV ) ! MV DEBUG 566 577 567 578 IF( lwp ) WRITE(numout,*) ' outer loop 1f i_out : ', i_out … … 575 586 576 587 ! --- Stress contributions at T-points 577 DO jj = 2, jpj ! loop to jpi,jpj to avoid making a communication for zs1 ,zs2,zs12588 DO jj = 2, jpj ! loop to jpi,jpj to avoid making a communication for zs1 & zs2 578 589 DO ji = 2, jpi ! 579 590 580 591 ! sig1 contribution to RHS of U-equation at T-points 581 zs1_rhsu(ji,jj) = zzt(ji,jj) * ( e1v(ji,jj) * zv_c(ji,jj) - e1v(ji,jj-1) * zv_c(ji,jj-1) - 1.0_wp ) 592 zs1_rhsu(ji,jj) = zzt(ji,jj) * ( e1v(ji,jj) * zv_c(ji,jj) - e1v(ji,jj-1) * zv_c(ji,jj-1) ) & 593 & - zp_deltastar_t(ji,jj) * zdeltat(ji,jj) 582 594 583 595 ! sig2 contribution to RHS of U-equation at T-points … … 585 597 586 598 ! sig1 contribution to RHS of V-equation at T-points 587 zs1_rhsv(ji,jj) = zzt(ji,jj) * ( e2u(ji,jj) * zu_c(ji,jj) - e2u(ji-1,jj) * zu_c(ji-1,jj) - 1.0_wp ) 599 zs1_rhsv(ji,jj) = zzt(ji,jj) * ( e2u(ji,jj) * zu_c(ji,jj) - e2u(ji-1,jj) * zu_c(ji-1,jj) ) & 600 & - zp_deltastar_t(ji,jj) * zdeltat(ji,jj) 588 601 589 602 ! sig2 contribution to RHS of V-equation at T-points … … 593 606 END DO 594 607 595 CALL iom_put( 'zs1_rhsu' , zs1_rhsu ) ! MV DEBUG 596 CALL iom_put( 'zs2_rhsu' , zs2_rhsu ) ! MV DEBUG 597 CALL iom_put( 'zs1_rhsv' , zs1_rhsv ) ! MV DEBUG 598 CALL iom_put( 'zs2_rhsv' , zs2_rhsv ) ! MV DEBUG 608 ! TEST 609 ! zs2_rhsv(:,:) = 0._wp 610 611 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs1_rhsu' , zs1_rhsu ) ! MV DEBUG 612 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs2_rhsu' , zs2_rhsu ) ! MV DEBUG 613 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs1_rhsv' , zs1_rhsv ) ! MV DEBUG 614 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs2_rhsv' , zs2_rhsv ) ! MV DEBUG 599 615 600 616 ! a priori, no lbclnk, because rhsu is only used in the inner domain 601 617 602 ! --- Stress contributions at f-points618 ! --- Stress contributions at F-points 603 619 ! MV NOTE: I applied zfmask on zds, by mimetism on EVP, but without deep understanding of what I was doing 604 620 ! My guess is that this is the way to enforce boundary conditions on strain rate tensor … … 610 626 611 627 ! sig12 contribution to RHS of U equation at F-points 612 zs12_rhsu(ji,jj) = - zef(ji,jj) * ( r1_e2v(ji+1,jj) * zv_c(ji+1,jj) -r1_e2v(ji,jj) * zv_c(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) * zfmask(ji,jj)628 zs12_rhsu(ji,jj) = zef(ji,jj) * ( r1_e2v(ji+1,jj) * zv_c(ji+1,jj) + r1_e2v(ji,jj) * zv_c(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) * zfmask(ji,jj) 613 629 614 630 ! sig12 contribution to RHS of V equation at F-points 615 zs12_rhsv(ji,jj) = zef(ji,jj) * ( r1_e1u(ji,jj+1) * zu_c(ji,jj+1) -r1_e1u(ji,jj) * zu_c(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) * zfmask(ji,jj)631 zs12_rhsv(ji,jj) = zef(ji,jj) * ( r1_e1u(ji,jj+1) * zu_c(ji,jj+1) + r1_e1u(ji,jj) * zu_c(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) * zfmask(ji,jj) 616 632 617 633 END DO 618 634 END DO 635 ! 636 ! TEST 637 ! zs12_rhsv(:,:) = 0._wp 619 638 620 639 CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsu, 'F', 1. ) 621 640 CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsv, 'F', 1. ) 622 641 623 CALL iom_put( 'zs12_rhsu' , zs12_rhsu ) ! MV DEBUG624 CALL iom_put( 'zs12_rhsv' , zs12_rhsv ) ! MV DEBUG642 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs12_rhsu' , zs12_rhsu ) ! MV DEBUG 643 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs12_rhsv' , zs12_rhsv ) ! MV DEBUG 625 644 626 645 ! a priori, no lbclnk, because rhsu are only used in the inner domain … … 631 650 DO jj = 2, jpj - 1 632 651 DO ji = 2, jpi - 1 652 633 653 ! --- U component of internal force contribution to RHS at U points 634 654 zf_rhsu(ji,jj) = 0.5_wp * r1_e1e2u(ji,jj) * & … … 640 660 zf_rhsv(ji,jj) = 0.5_wp * r1_e1e2v(ji,jj) * & 641 661 & ( e1v(ji,jj) * ( zs1_rhsv(ji,jj+1) - zs1_rhsv(ji,jj) ) & 642 & + r1_e1v(ji,jj) * ( e1t(ji,jj+1) * e1t(ji,jj+1) * zs2_rhsv(ji,jj+1) - e1t(ji,jj) * e1t(ji,jj)* zs2_rhsv(ji,jj) ) &662 & - r1_e1v(ji,jj) * ( e1t(ji,jj+1) * e1t(ji,jj+1) * zs2_rhsv(ji,jj+1) - e1t(ji,jj) * e1t(ji,jj) * zs2_rhsv(ji,jj) ) & 643 663 & + 2._wp * r1_e2v(ji,jj) * ( e2f(ji,jj) * e2f(ji,jj) * zs12_rhsv(ji,jj) - e2f(ji-1,jj) * e2f(ji-1,jj) * zs12_rhsv(ji-1,jj) ) ) 644 664 645 665 END DO 646 666 END DO 647 667 648 CALL iom_put( 'zf_rhsu' , zf_rhsu ) ! MV DEBUG649 CALL iom_put( 'zf_rhsv' , zf_rhsv ) ! MV DEBUG668 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zf_rhsu' , zf_rhsu ) ! MV DEBUG 669 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zf_rhsv' , zf_rhsv ) ! MV DEBUG 650 670 651 671 !--------------------------- … … 657 677 DO ji = 2, jpi - 1 658 678 659 ! still miss ice ocean stress and acceleration contribution660 679 zrhsu(ji,jj) = zmU_t(ji,jj) + ztaux_ai(ji,jj) + ztaux_oi_rhsu(ji,jj) + zspgU(ji,jj) + zCorU(ji,jj) + zf_rhsu(ji,jj) 661 zrhsv(ji,jj) = zmV_t(ji,jj) + ztauy_ai(ji,jj) + ztauy_oi_rhsv(ji,jj) + zspgV(ji,jj) + zCorV(ji,jj) + zf_rhs u(ji,jj)680 zrhsv(ji,jj) = zmV_t(ji,jj) + ztauy_ai(ji,jj) + ztauy_oi_rhsv(ji,jj) + zspgV(ji,jj) + zCorV(ji,jj) + zf_rhsv(ji,jj) 662 681 663 682 END DO … … 668 687 CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi_rhsu, 'U', -1., ztauy_oi_rhsv, 'V', -1.) 669 688 670 CALL iom_put( 'zmU_t' , zmU_t ) ! MV DEBUG671 CALL iom_put( 'zmV_t' , zmV_t ) ! MV DEBUG672 CALL iom_put( 'ztaux_oi_rhsu' , ztaux_oi_rhsu ) ! MV DEBUG673 CALL iom_put( 'ztauy_oi_rhsv' , ztauy_oi_rhsv ) ! MV DEBUG674 CALL iom_put( 'zrhsu' , zrhsu ) ! MV DEBUG675 CALL iom_put( 'zrhsv' , zrhsv ) ! MV DEBUG689 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zmU_t' , zmU_t ) ! MV DEBUG 690 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zmV_t' , zmV_t ) ! MV DEBUG 691 IF ( i_out == nn_nout_vp ) CALL iom_put( 'ztaux_oi_rhsu' , ztaux_oi_rhsu ) ! MV DEBUG 692 IF ( i_out == nn_nout_vp ) CALL iom_put( 'ztauy_oi_rhsv' , ztauy_oi_rhsv ) ! MV DEBUG 693 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zrhsu' , zrhsu ) ! MV DEBUG 694 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zrhsv' , zrhsv ) ! MV DEBUG 676 695 677 696 ! inner domain calculations -> no lbclnk … … 712 731 zfac3 = 2._wp * zfac * r1_e1u(ji,jj) 713 732 714 zt12U = - zfac1 * zzt(ji+1,jj)715 733 zt11U = zfac1 * zzt(ji,jj) 716 717 zt22U = - zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj)734 zt12U = zfac1 * zzt(ji+1,jj) 735 718 736 zt21U = zfac2 * zet(ji,jj) * e2t(ji,jj) * e2t(ji,jj) * e2t(ji,jj) * e2t(ji,jj) 719 720 zt122U = - zfac3 * zef(ji,jj) * e1f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) * e1f(ji,jj)737 zt22U = zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) 738 721 739 zt121U = zfac3 * zef(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) 740 zt122U = zfac3 * zef(ji,jj) * e1f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) * e1f(ji,jj) 722 741 723 742 ! 724 743 ! Linear system coefficients 725 744 ! 726 zAU(ji,jj) = - zt11U * e2u(ji-1,jj) - zt21U* r1_e2u(ji-1,jj)727 zBU(ji,jj) = ( zt12U + zt11U ) * e2u(ji,jj) + ( zt22U + zt21U ) * r1_e2u(ji,jj) + ( zt122U + zt121U ) * r1_e1u(ji,jj)728 zCU(ji,jj) = - zt12U * e2u(ji+1,jj) - zt22U* r1_e2u(ji+1,jj)729 730 zDU(ji,jj) = zt121U * r1_e1u(ji,jj-1)731 zEU(ji,jj) = zt122U * r1_e1u(ji,jj+1)745 zAU(ji,jj) = - zt11U * e2u(ji-1,jj) - zt21U * r1_e2u(ji-1,jj) 746 zBU(ji,jj) = ( zt11U + zt12U ) * e2u(ji,jj) + ( zt21U + zt22U ) * r1_e2u(ji,jj) + ( zt121U + zt122U ) * r1_e1u(ji,jj) 747 zCU(ji,jj) = - zt12U * e2u(ji+1,jj) - zt22U * r1_e2u(ji+1,jj) 748 749 zDU(ji,jj) = zt121U * r1_e1u(ji,jj-1) 750 zEU(ji,jj) = zt122U * r1_e1u(ji,jj+1) 732 751 733 752 ! … … 737 756 ! 738 757 zfac = 0.5_wp * r1_e1e2v(ji,jj) 739 zfac1 = zfac * e 2v(ji,jj)758 zfac1 = zfac * e1v(ji,jj) 740 759 zfac2 = zfac * r1_e1v(ji,jj) 741 760 zfac3 = 2._wp * zfac * r1_e2v(ji,jj) 742 761 743 zt12V = - zfac1 * zzt(ji,jj+1)744 762 zt11V = zfac1 * zzt(ji,jj) 745 763 zt12V = zfac1 * zzt(ji,jj+1) 764 765 zt21V = zfac2 * zet(ji,jj) * e1t(ji,jj) * e1t(ji,jj) * e1t(ji,jj) * e1t(ji,jj) 746 766 zt22V = zfac2 * zet(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) 747 zt21V = - zfac2 * zet(ji,jj) * e1t(ji,jj) * e1t(ji,jj) * e1t(ji,jj) * e1t(ji,jj)748 767 768 zt121V = zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) 749 769 zt122V = zfac3 * zef(ji,jj) * e2f(ji,jj) * e2f(ji,jj) * e2f(ji,jj) * e2f(ji,jj) 750 zt121V = - zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) 751 770 752 771 ! 753 772 ! Linear system coefficients 754 773 ! 755 zAV(ji,jj) = - zt11V * e1v(ji,jj-1) + zt21V* r1_e1v(ji,jj-1)756 zBV(ji,jj) = ( zt12V + zt11V ) * e1v(ji,jj) - ( zt22V + zt21V ) * r1_e1v(ji,jj) -( zt122V + zt121V ) * r1_e2v(ji,jj)757 zCV(ji,jj) = - zt12V * e1v(ji,jj+1) + zt22V* r1_e1v(ji,jj+1)758 759 zDV(ji,jj) = - zt121V * r1_e2v(ji-1,jj) ! mistake is in the pdf notes not here760 zEV(ji,jj) = -zt122V * r1_e2v(ji+1,jj)774 zAV(ji,jj) = - zt11V * e1v(ji,jj-1) - zt21V * r1_e1v(ji,jj-1) 775 zBV(ji,jj) = ( zt11V + zt12V ) * e1v(ji,jj) + ( zt21V + zt22V ) * r1_e1v(ji,jj) + ( zt122V + zt121V ) * r1_e2v(ji,jj) 776 zCV(ji,jj) = - zt12V * e1v(ji,jj+1) - zt22V * r1_e1v(ji,jj+1) 777 778 zDV(ji,jj) = zt121V * r1_e2v(ji-1,jj) 779 zEV(ji,jj) = zt122V * r1_e2v(ji+1,jj) 761 780 762 781 !----------------------------------------------------- … … 764 783 !----------------------------------------------------- 765 784 zBU(ji,jj) = zBU(ji,jj) + zCwU(ji,jj) + zmassU_t(ji,jj) 766 zBV(ji,jj) = ZBV(ji,jj) + zCwV(ji,jj) + zmassV_t(ji,jj)785 zBV(ji,jj) = zBV(ji,jj) + zCwV(ji,jj) + zmassV_t(ji,jj) 767 786 768 787 END DO 769 788 END DO 770 789 771 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zAU , 'U', 1., zAV , 'V', 1. ) 790 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zAU , 'U', 1., zAV , 'V', 1. ) ! --> here normal that we use '1' and not '-1' ? 772 791 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zBU , 'U', 1., zBV , 'V', 1. ) 773 792 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCU , 'U', 1., zCV , 'V', 1. ) … … 775 794 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zEU , 'U', 1., zEV , 'V', 1. ) 776 795 777 CALL iom_put( 'zAU' , zAU ) ! MV DEBUG 778 CALL iom_put( 'zBU' , zBU ) ! MV DEBUG 779 CALL iom_put( 'zCU' , zCU ) ! MV DEBUG 780 CALL iom_put( 'zDU' , zDU ) ! MV DEBUG 781 CALL iom_put( 'zEU' , zEU ) ! MV DEBUG 782 CALL iom_put( 'zAV' , zAV ) ! MV DEBUG 783 CALL iom_put( 'zBV' , zBV ) ! MV DEBUG 784 CALL iom_put( 'zCV' , zCV ) ! MV DEBUG 785 CALL iom_put( 'zDV' , zDV ) ! MV DEBUG 786 CALL iom_put( 'zEV' , zEV ) ! MV DEBUG 796 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zAU' , zAU ) ! MV DEBUG 797 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zBU' , zBU ) ! MV DEBUG 798 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zCU' , zCU ) ! MV DEBUG 799 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zDU' , zDU ) ! MV DEBUG 800 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zEU' , zEU ) ! MV DEBUG 801 802 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zAV' , zAV ) ! MV DEBUG 803 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zBV' , zBV ) ! MV DEBUG 804 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zCV' , zCV ) ! MV DEBUG 805 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zDV' , zDV ) ! MV DEBUG 806 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zEV' , zEV ) ! MV DEBUG 787 807 788 808 !------------------------------------------------------------------------------! … … 808 828 zv_b(:,:) = v_ice(:,:) 809 829 810 !zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp811 !zAV(:,:) = 0._wp; zBV(:,:) = 0._wp; zCV(:,:) = 0._wp; zDV(:,:) = 0._wp; zEV(:,:) = 0._wp830 !zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp 831 !zAV(:,:) = 0._wp; zBV(:,:) = 0._wp; zCV(:,:) = 0._wp; zDV(:,:) = 0._wp; zEV(:,:) = 0._wp 812 832 813 833 IF ( ll_u_iterate .OR. ll_v_iterate ) THEN … … 864 884 DO jj = jj_min, jpj - 1, nn_zebra_vp 865 885 886 ! MV BUG OCT 6. Both zBU_prime and zFU_prime were undefined for ji = 2 887 zBU_prime(2,jj) = zBU(2,jj) 888 zFU_prime(2,jj) = zFU(2,jj) 889 866 890 DO ji = 3, jpi - 1 867 891 … … 962 986 DO ji = ji_min, jpi - 1, nn_zebra_vp 963 987 964 DO jj = 3, jpj - 1 988 ! MV BUG OCT 6.. Both zBV_prime and zFV_prime were undefined for jj = 2 989 zBV_prime(ji,2) = zBV(ji,2) 990 zFV_prime(ji,2) = zFV(ji,2) 991 992 DO jj = 3, jpj - 1 965 993 966 994 zfac = SIGN( 1._wp , zBV(ji,jj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj-1) ) - epsi20 ) ) … … 1137 1165 1138 1166 ! MV DEBUG test - replace ice velocity by ocean current to give the model the means to go ahead 1139 DO jj = 2, jpj - 1 1140 DO ji = 2, jpi - 1 1141 1142 u_ice(ji,jj) = zmsk00x(ji,jj) & 1143 & * ( zmsk01x(ji,jj) * u_oce(ji,jj) * 0.01_wp & 1144 + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp ) 1145 1146 v_ice(ji,jj) = zmsk00y(ji,jj) & 1147 & * ( zmsk01y(ji,jj) * v_oce(ji,jj) * 0.01_wp & 1148 + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp ) 1149 1150 END DO 1151 END DO 1152 1153 CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 1154 1155 IF ( lwp ) WRITE(numout,*) ' Velocity replaced ' 1167 ! info: probably wrong because when doing du/dt = div ( sigma ) we have spurious waves instead of u = v = 0 1168 1169 ! DO jj = 2, jpj - 1 1170 ! DO ji = 2, jpi - 1 1171 1172 ! u_ice(ji,jj) = zmsk00x(ji,jj) & 1173 ! & * ( zmsk01x(ji,jj) * u_oce(ji,jj) * 0.01_wp & 1174 ! + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp ) 1175 1176 ! v_ice(ji,jj) = zmsk00y(ji,jj) & 1177 ! & * ( zmsk01y(ji,jj) * v_oce(ji,jj) * 0.01_wp & 1178 ! + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp ) 1179 1180 ! END DO 1181 ! END DO 1182 1183 ! CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 1184 1185 ! IF ( lwp ) WRITE(numout,*) ' Velocity replaced ' 1156 1186 1157 1187 ! END DEBUG
Note: See TracChangeset
for help on using the changeset viewer.
