Changeset 15527
- Timestamp:
- 2021-11-21T11:10:07+01:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90
r15523 r15527 152 152 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV ! ice/snow mass and volume 153 153 REAL(wp) :: zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 154 REAL(wp) :: zp_del tastar_f !154 REAL(wp) :: zp_delstar_f ! 155 155 REAL(wp) :: zu_cV, zv_cU ! 156 156 REAL(wp) :: zfac, zfac1, zfac2, zfac3 … … 164 164 REAL(wp), DIMENSION(jpi,jpj) :: zmassU_t, zmassV_t ! Mass per unit area divided by time step 165 165 ! 166 REAL(wp), DIMENSION(jpi,jpj) :: zdeltat, zdel tastar_t ! Delta & Delta* at T-points166 REAL(wp), DIMENSION(jpi,jpj) :: zdeltat, zdelstar_t ! Delta & Delta* at T-points 167 167 REAL(wp), DIMENSION(jpi,jpj) :: ztens, zshear ! Tension, shear 168 REAL(wp), DIMENSION(jpi,jpj) :: zp_del tastar_t ! P/delta* at T points168 REAL(wp), DIMENSION(jpi,jpj) :: zp_delstar_t ! P/delta* at T points 169 169 REAL(wp), DIMENSION(jpi,jpj) :: zzt, zet ! Viscosity pre-factors at T points 170 170 REAL(wp), DIMENSION(jpi,jpj) :: zef ! Viscosity pre-factor at F point … … 194 194 REAL(wp), DIMENSION(jpi,jpj) :: zFU, zFU_prime, zBU_prime ! Rearranged linear system coefficients, U equation 195 195 REAL(wp), DIMENSION(jpi,jpj) :: zFV, zFV_prime, zBV_prime ! Rearranged linear system coefficients, V equation 196 REAL(wp), DIMENSION(jpi,jpj) :: zCU_prime, zCV_prime ! Rearranged linear system coefficients, V equation197 196 !!! REAL(wp), DIMENSION(jpi,jpj) :: ztaux_bi, ztauy_bi ! ice-OceanBottom stress at U-V points (landfast) 198 197 !!! REAL(wp), DIMENSION(jpi,jpj) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast) … … 302 301 END DO 303 302 END DO 304 305 !CALL lbc_lnk( 'icedyn_rhg_vp', zfmask, 'F', 1._wp ) ! OK to remove306 303 307 304 ! Lateral boundary conditions on velocity (modify zfmask) … … 434 431 END DO 435 432 436 CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! necessary 433 CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! necessary, zds2 uses jpi/jpj values for zds 437 434 438 435 DO jj = 2, jpj ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 … … 460 457 461 458 ! delta* at T points (following Lemieux and Dupont, GMD 2020) 462 zdel tastar_t(ji,jj) = zdeltat(ji,jj) + rn_creepl459 zdelstar_t(ji,jj) = zdeltat(ji,jj) + rn_creepl ! OPT zdelstar_t can be totally removed and put into next line directly. Could change results 463 460 464 ! P/delta at T-points465 zp_del tastar_t(ji,jj) = strength(ji,jj) / zdeltastar_t(ji,jj)461 ! P/delta* at T-points 462 zp_delstar_t(ji,jj) = strength(ji,jj) / zdelstar_t(ji,jj) 466 463 467 464 ! Temporary zzt and zet factors at T-points 468 zzt(ji,jj) = zp_del tastar_t(ji,jj) * r1_e1e2t(ji,jj)465 zzt(ji,jj) = zp_delstar_t(ji,jj) * r1_e1e2t(ji,jj) 469 466 zet(ji,jj) = zzt(ji,jj) * z1_ecc2 470 467 … … 472 469 END DO 473 470 474 ! CALL lbc_lnk_multi( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. , zdeltat, 'T', 1.) ! OK to 475 ! remove zet and zzt 476 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1., zdeltat, 'T', 1.) 471 CALL lbc_lnk( 'icedyn_rhg_vp', zp_delstar_t , 'T', 1. ) ! necessary, used for ji = 1 and jj = 1 477 472 478 473 DO jj = 1, jpj - 1 … … 480 475 481 476 ! P/delta* at F points 482 zp_del tastar_f = 0.25_wp * ( zp_deltastar_t(ji,jj) + zp_deltastar_t(ji+1,jj) + zp_deltastar_t(ji,jj+1) + zp_deltastar_t(ji+1,jj+1) )477 zp_delstar_f = 0.25_wp * ( zp_delstar_t(ji,jj) + zp_delstar_t(ji+1,jj) + zp_delstar_t(ji,jj+1) + zp_delstar_t(ji+1,jj+1) ) 483 478 484 479 ! Temporary zef factor at F-point 485 zef(ji,jj) = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2 * zfmask(ji,jj) * 0.5_wp ! BUG missing un facteur 2 Nov 12 486 487 END DO 488 END DO 489 490 ! CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. ) ! Ok to remove 480 zef(ji,jj) = zp_delstar_f * r1_e1e2f(ji,jj) * z1_ecc2 * zfmask(ji,jj) * 0.5_wp 481 482 END DO 483 END DO 491 484 492 485 !--------------------------------------------------- … … 524 517 END DO 525 518 526 ! CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCwU , 'U', -1., zCwV, 'V', -1. )! OK to remove527 ! CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCorU, 'U', -1., zCorV, 'V', -1. )! OK to remove528 529 ! a priori, Coriolis and drag terms only affect diagonal or independent term of the linear system,530 ! so there is no need for lbclnk on drag and coriolis531 532 519 !------------------------------------- 533 520 ! -- Internal stress RHS contribution … … 540 527 ! sig1 contribution to RHS of U-equation at T-points 541 528 zs1_rhsu(ji,jj) = zzt(ji,jj) * ( e1v(ji,jj) * zv_c(ji,jj) - e1v(ji,jj-1) * zv_c(ji,jj-1) ) & 542 & - zp_del tastar_t(ji,jj) * zdeltat(ji,jj)529 & - zp_delstar_t(ji,jj) * zdeltat(ji,jj) 543 530 544 531 ! sig2 contribution to RHS of U-equation at T-points … … 547 534 ! sig1 contribution to RHS of V-equation at T-points 548 535 zs1_rhsv(ji,jj) = zzt(ji,jj) * ( e2u(ji,jj) * zu_c(ji,jj) - e2u(ji-1,jj) * zu_c(ji-1,jj) ) & 549 & - zp_del tastar_t(ji,jj) * zdeltat(ji,jj)536 & - zp_delstar_t(ji,jj) * zdeltat(ji,jj) 550 537 551 538 ! sig2 contribution to RHS of V-equation at T-points … … 554 541 END DO 555 542 END DO 556 557 ! a priori, no lbclnk, because rhsu is only used in the inner domain558 543 559 544 ! --- Stress contributions at F-points … … 572 557 END DO 573 558 END DO 574 575 ! CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsu, 'F', 1. ) ! Ok to remove576 ! CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsv, 'F', 1. ) ! Ok to remove577 578 ! a priori, no lbclnk, because rhsu are only used in the inner domain579 559 580 560 ! --- Internal force contributions to RHS, taken as divergence of stresses (Appendix C of Hunke and Dukowicz, 2002) … … 613 593 END DO 614 594 615 ! CALL lbc_lnk_multi( 'icedyn_rhg_vp', zrhsu, 'U', -1., zrhsv, 'V', -1.)! OK to remove616 ! CALL lbc_lnk_multi( 'icedyn_rhg_vp', zmU_t, 'U', -1., zmV_t, 'V', -1.)! OK to remove617 ! CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi_rhsu, 'U', -1., ztauy_oi_rhsv, 'V', -1.)! OK to remove618 619 595 !---------------------------------------------------------------------------------------! 620 596 ! … … 707 683 END DO 708 684 709 ! OK to remove710 !CALL lbc_lnk_multi( 'icedyn_rhg_vp', zAU , 'U', 1., zAV , 'V', 1. ) ! --> here normal that we use '1' and not '-1' ?711 !CALL lbc_lnk_multi( 'icedyn_rhg_vp', zBU , 'U', 1., zBV , 'V', 1. )712 !CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCU , 'U', 1., zCV , 'V', 1. )713 !CALL lbc_lnk_multi( 'icedyn_rhg_vp', zDU , 'U', 1., zDV , 'V', 1. )714 !CALL lbc_lnk_multi( 'icedyn_rhg_vp', zEU , 'U', 1., zEV , 'V', 1. )715 716 685 !------------------------------------------------------------------------------! 717 686 ! … … 745 714 ! A*u(i-1,j)+B*u(i,j)+C*u(i+1,j) = F 746 715 747 zFU(:,:) = 0._wp ; zFU_prime(:,:) = 0._wp ; zBU_prime(:,:) = 0._wp; zCU_prime(:,:) = 0._wp716 zFU(:,:) = 0._wp ; zFU_prime(:,:) = 0._wp ; zBU_prime(:,:) = 0._wp; 748 717 749 718 DO jn = 1, nn_zebra_vp ! "zebra" loop (! red-black-sor!!! ) … … 761 730 !------------------------ 762 731 DO ji = 2, jpi - 1 763 764 ! boundary condition substitution 732 ! note: these are key lines linking information between processors 733 ! u_ice/v_ice need to be lbc_linked 734 735 ! sub-domain boundary condition substitution 765 736 ! see Zhang and Hibler, 1997, Appendix B 766 737 zAA3 = 0._wp … … 778 749 END DO 779 750 780 CALL lbc_lnk( 'icedyn_rhg_vp', zFU, 'U', -1. ) ! -1, zFU is a vector781 782 751 !--------------- 783 752 ! Forward sweep … … 785 754 DO jj = jj_min, jpj - 1, nn_zebra_vp 786 755 787 ! MV BUG OCT 6. Both zBU_prime and zFU_prime were undefined for ji = 2788 756 zBU_prime(2,jj) = zBU(2,jj) 789 757 zFU_prime(2,jj) = zFU(2,jj) … … 799 767 800 768 END DO 801 ! CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFU_prime, 'U', -1., zBU_prime, 'U', 1. ) --> can be removed802 769 803 770 !----------------------------- … … 820 787 END DO 821 788 822 !--- Relaxation 823 ! and velocity masking for little-ice and no-ice cases 789 !--- Relaxation and masking (for low-ice/no-ice cases) 824 790 DO ji = 2, jpi - 1 825 791 … … 848 814 !!! ZH97 explain it is critical for convergence speed 849 815 850 zFV(:,:) = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp; zCV_prime(:,:) = 0._wp816 zFV(:,:) = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp; 851 817 852 818 DO jn = 1, nn_zebra_vp ! "zebra" loop … … 863 829 DO jj = 2, jpj - 1 864 830 865 ! boundary condition substitution (check it is correctly applied !!!)831 ! subdomain boundary condition substitution (check it is correctly applied !!!) 866 832 ! see Zhang and Hibler, 1997, Appendix B 867 833 zAA3 = 0._wp … … 879 845 END DO 880 846 881 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFV, 'V', -1.) ! MV Nov 5 bug test to change boundary condition to -1882 ! since zFV is a vector!!!883 884 847 !--------------- 885 848 ! Forward sweep … … 887 850 DO ji = ji_min, jpi - 1, nn_zebra_vp 888 851 889 ! MV BUG OCT 6.. Both zBV_prime and zFV_prime were undefined for jj = 2890 852 zBV_prime(ji,2) = zBV(ji,2) 891 853 zFV_prime(ji,2) = zFV(ji,2) … … 902 864 END DO 903 865 904 ! new test905 ! CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFV_prime, 'V', -1., zBV_prime, 'V', 1. ) ! MV Nov 5 bug test906 907 866 !----------------------------- 908 867 ! Backward sweep & relaxation … … 923 882 END DO 924 883 925 ! --- Relaxation & masking (should it be now or later)884 ! --- Relaxation & masking 926 885 DO jj = 2, jpj - 1 927 886 … … 940 899 ENDIF ! ll_v_iterate 941 900 901 ! I suspect the communication should go into the zebra loop if we want reproducibility 942 902 CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 943 903 … … 1031 991 1032 992 END DO ! End of outer loop (i_out) ============================================================================================= 1033 1034 CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )1035 993 1036 994 IF( ln_rhg_chkcvg ) THEN … … 1123 1081 DO jj = 2, jpj - 1 1124 1082 DO ji = 2, jpi - 1 1125 zp_del tastar_t(ji,jj)= strength(ji,jj) / pdelta_i(ji,jj)1126 zfac = zp_del tastar_t(ji,jj)1083 zp_delstar_t(ji,jj) = strength(ji,jj) / pdelta_i(ji,jj) 1084 zfac = zp_delstar_t(ji,jj) 1127 1085 zs1(ji,jj) = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 1128 1086 zs2(ji,jj) = zfac * z1_ecc2 * ztens(ji,jj) … … 1142 1100 1143 1101 ! P/delta* at F points 1144 zp_del tastar_f = 0.25_wp * ( zp_deltastar_t(ji,jj) + zp_deltastar_t(ji+1,jj) + zp_deltastar_t(ji,jj+1) + zp_deltastar_t(ji+1,jj+1) )1102 zp_delstar_f = 0.25_wp * ( zp_delstar_t(ji,jj) + zp_delstar_t(ji+1,jj) + zp_delstar_t(ji,jj+1) + zp_delstar_t(ji+1,jj+1) ) 1145 1103 1146 1104 ! s12 at F-points 1147 zs12f(ji,jj) = zp_del tastar_f * z1_ecc2 * zds(ji,jj)1105 zs12f(ji,jj) = zp_delstar_f * z1_ecc2 * zds(ji,jj) 1148 1106 1149 1107 END DO … … 1251 1209 ! and **deformations** at current iterates 1252 1210 ! following Lemieux & Dupont (2020) 1253 zfac = zp_del tastar_t(ji,jj)1211 zfac = zp_delstar_t(ji,jj) 1254 1212 zsig1 = zfac * ( pdivu_i(ji,jj) - zdeltat(ji,jj) ) 1255 1213 zsig2 = zfac * z1_ecc2 * ztens(ji,jj) … … 1618 1576 ! Local residual r = Ax - B expresses to which extent the momentum balance is verified 1619 1577 ! i.e., how close we are to a solution 1620 !!!! UNITS OF RESIDUAL ARE NOT m/s BUT N/m2 (CHECK)1621 1578 zu_res(:,:) = 0._wp; zv_res(:,:) = 0._wp 1622 1579
Note: See TracChangeset
for help on using the changeset viewer.
