Changeset 13591
- Timestamp:
- 2020-10-14T16:38:22+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/SI3-03_VP_rheology/src/ICE
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_evp.F90
r13544 r13591 138 138 REAL(wp), DIMENSION(jpi,jpj) :: zp_delt ! P/delta at T points 139 139 REAL(wp), DIMENSION(jpi,jpj) :: zdeltastar_t ! Delta* at T-points 140 REAL(wp), DIMENSION(jpi,jpj) :: zten_i ! Tension 140 141 REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017 141 142 ! … … 170 171 !! --- diags 171 172 REAL(wp) :: zsig1, zsig2, zsig12, zfac, z1_strength 172 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p , zten_i173 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p 173 174 !! --- SIMIP diags 174 175 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) … … 809 810 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence 810 811 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear 812 IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * zmsk00 ) ! delta 811 813 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 812 814 … … 830 832 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 831 833 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 832 zsig_II(ji,jj) = -SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 ) ! 2nd '' '', aka maximum shear stress834 zsig_II(ji,jj) = SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 ) ! 2nd '' '', aka maximum shear stress 833 835 834 836 END DO … … 852 854 IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 853 855 ! 854 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) )856 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 855 857 ! 856 858 DO jj = 2, jpj - 1 … … 867 869 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 868 870 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 869 zsig_II(ji,jj) = - SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 )! 2nd '' '', aka maximum shear stress871 zsig_II(ji,jj) = SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 ) ! 2nd '' '', aka maximum shear stress 870 872 871 873 ! Normalized principal stresses (used to display the ellipse) 872 z1_strength = 1._wp / strength(ji,jj)873 zsig1_p(ji,jj) = ( zsig_I(ji,jj) -zsig_II(ji,jj) ) * z1_strength874 zsig2_p(ji,jj) = ( zsig_I I(ji,jj) +zsig_II(ji,jj) ) * z1_strength874 z1_strength = 1._wp / MAX ( 1.0, strength(ji,jj) ) 875 zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 876 zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 875 877 END DO 876 878 END DO … … 882 884 CALL iom_put( 'sig2_pnorm' , zsig2_p ) 883 885 884 DEALLOCATE( zsig1_p , zsig2_p )886 DEALLOCATE( zsig1_p , zsig2_p , zsig_I , zsig_II ) 885 887 886 888 ENDIF -
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90
r13544 r13591 50 50 INTEGER :: nvarid_vdif 51 51 INTEGER :: nvarid_mke 52 INTEGER :: nvarid_usub, nvarid_vsub 53 52 54 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmsk00, zmsk15 53 55 … … 129 131 LOGICAL :: ll_u_iterate, ll_v_iterate ! continue iteration or not 130 132 ! 131 INTEGER :: ji, j j, jn ! dummy loop indices133 INTEGER :: ji, ji2, jj, jn ! dummy loop indices 132 134 INTEGER :: jter, i_out, i_inn ! 133 135 INTEGER :: ji_min, jj_min ! 134 136 INTEGER :: nn_zebra_vp ! number of zebra steps 137 138 INTEGER :: nn_thomas ! tridiagonal system version (1=mitgcm, 2=martin) 139 135 140 INTEGER :: nn_nvp ! total number of VP iterations (n_out_vp*n_inn_vp) 136 141 ! … … 146 151 REAL(wp) :: zt12U, zt11U, zt22U, zt21U, zt122U, zt121U 147 152 REAL(wp) :: zt12V, zt11V, zt22V, zt21V, zt122V, zt121V 148 REAL(wp) :: zAA3, zw, z uerr_max, zverr_max153 REAL(wp) :: zAA3, zw, ztau, zuerr_max, zverr_max 149 154 ! 150 155 REAL(wp), DIMENSION(jpi,jpj) :: zfmask ! mask at F points for the ice … … 154 159 ! 155 160 REAL(wp), DIMENSION(jpi,jpj) :: zdeltastar_t ! Delta* at T-points 161 REAL(wp), DIMENSION(jpi,jpj) :: zten_i ! Tension 156 162 REAL(wp), DIMENSION(jpi,jpj) :: zp_deltastar_t ! P/delta* at T points 157 163 REAL(wp), DIMENSION(jpi,jpj) :: zzt, zet ! Viscosity pre-factors at T points … … 182 188 REAL(wp), DIMENSION(jpi,jpj) :: zFU, zFU_prime, zBU_prime ! Rearranged linear system coefficients, U equation 183 189 REAL(wp), DIMENSION(jpi,jpj) :: zFV, zFV_prime, zBV_prime ! Rearranged linear system coefficients, V equation 190 REAL(wp), DIMENSION(jpi,jpj) :: zCU_prime, zCV_prime ! Rearranged linear system coefficients, V equation 184 191 !!! REAL(wp), DIMENSION(jpi,jpj) :: ztaux_bi, ztauy_bi ! ice-OceanBottom stress at U-V points (landfast) 185 192 !!! REAL(wp), DIMENSION(jpi,jpj) :: ztaux_base, ztauy_base ! ice-bottom stress at U-V points (landfast) … … 193 200 !! --- diags 194 201 REAL(wp) :: zsig1, zsig2, zsig12, zdelta, z1_strength, zfac_x, zfac_y 195 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12, zs12f , zten_i! stress tensor components202 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12, zs12f ! stress tensor components 196 203 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p 197 204 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztaux_oi, ztauy_oi … … 201 208 202 209 !!---------------------------------------------------------------------------------------------------------------------- 210 ! DEBUG put all forcing terms to zero 211 ! air-ice drag 212 utau_ice(:,:) = 0._wp 213 vtau_ice(:,:) = 0._wp 214 ! coriolis 215 ff_t(:,:) = 0._wp 216 ! ice-ocean drag 217 rn_cio = 0._wp 218 ! ssh 219 ! done line 330 !!! dont forget to act there 220 ! END DEBUG 203 221 204 222 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_dyn_rhg_vp: VP sea-ice rheology (LSR solver)' 223 IF( lwp ) WRITE(numout,*) '-- ice_dyn_rhg_vp: VP sea-ice rheology (LSR solver)' 205 224 206 225 !------------------------------------------------------------------------------! … … 219 238 END DO 220 239 221 IF ( ln_zebra_vp ) THEN; nn_zebra_vp = 2; ELSE; nn_zebra_vp = 1; ENDIF 222 240 IF ( ln_zebra_vp ) THEN; nn_zebra_vp = 2 241 ELSE; nn_zebra_vp = 1; ENDIF 242 223 243 nn_nvp = nn_nout_vp * nn_ninn_vp ! maximum number of iterations 244 245 IF( lwp ) WRITE(numout,*) ' ln_zebra_vp : ', ln_zebra_vp 246 IF( lwp ) WRITE(numout,*) ' nn_zebra_vp : ', nn_zebra_vp 247 IF( lwp ) WRITE(numout,*) ' nn_nvp : ', nn_nvp 224 248 225 249 zrhoco = rau0 * rn_cio … … 242 266 ELSE ; zkt = 0._wp 243 267 ENDIF 268 269 zs1_rhsu(:,:) = 0._wp; zs2_rhsu(:,:) = 0._wp; zs1_rhsv(:,:) = 0._wp; zs2_rhsv(:,:) = 0._wp 270 zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp; 271 zAV(:,:) = 0._wp; zBV(:,:) = 0._wp; zCV(:,:) = 0._wp; zDV(:,:) = 0._wp; zEV(:,:) = 0._wp; 272 zrhsu(:,:) = 0._wp; zrhsv(:,:) = 0._wp 273 zf_rhsu(:,:) = 0._wp; zf_rhsv(:,:) = 0._wp 244 274 245 275 !------------------------------------------------------------------------------! … … 304 334 ! non-embedded sea ice: use ocean surface for slope calculation 305 335 zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 336 zsshdyn(:,:) = 0._wp ! DEBUG CAREFUL !!! 306 337 307 338 DO jj = 2, jpj - 1 … … 349 380 ! Mask for lots of ice (1) or little ice (0) 350 381 IF( zmassU <= zmmin .AND. za_iU(ji,jj) <= zamin ) THEN ; zmsk01x(ji,jj) = 0._wp 351 ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF382 ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF 352 383 IF( zmassV <= zmmin .AND. za_iV(ji,jj) <= zamin ) THEN ; zmsk01y(ji,jj) = 0._wp 353 ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF384 ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF 354 385 355 386 END DO 356 387 END DO 388 389 CALL iom_put( 'zmsk00x' , zmsk00x ) ! MV DEBUG 390 CALL iom_put( 'zmsk00y' , zmsk00y ) ! MV DEBUG 391 CALL iom_put( 'zmsk01x' , zmsk01x ) ! MV DEBUG 392 CALL iom_put( 'zmsk01y' , zmsk01y ) ! MV DEBUG 393 CALL iom_put( 'ztauy_ai' , ztauy_ai ) ! MV DEBUG 394 CALL iom_put( 'ztaux_ai' , ztaux_ai ) ! MV DEBUG 395 CALL iom_put( 'ztauy_ai' , ztauy_ai ) ! MV DEBUG 396 CALL iom_put( 'zspgU' , zspgU ) ! MV DEBUG 397 CALL iom_put( 'zspgV' , zspgV ) ! MV DEBUG 357 398 358 399 !------------------------------------------------------------------------------! … … 368 409 369 410 DO i_out = 1, nn_nout_vp 411 412 IF( lwp ) WRITE(numout,*) ' outer loop i_out : ', i_out 370 413 371 414 ! Velocities used in the non linear terms are the average of the past two iterates … … 400 443 END DO 401 444 END DO 402 445 403 446 CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! MV TEST could be un-necessary according to Gurvan 447 CALL iom_put( 'zds' , zds ) ! MV DEBUG 448 449 IF( lwp ) WRITE(numout,*) ' outer loop 1a i_out : ', i_out 404 450 405 451 DO jj = 2, jpj - 1 ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 … … 441 487 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. ) 442 488 489 CALL iom_put( 'zzt' , zzt ) ! MV DEBUG 490 CALL iom_put( 'zet' , zet ) ! MV DEBUG 491 CALL iom_put( 'zp_deltastar_t', zp_deltastar_t ) ! MV DEBUG 492 493 IF( lwp ) WRITE(numout,*) ' outer loop 1b i_out : ', i_out 494 443 495 DO jj = 1, jpj - 1 444 496 DO ji = 1, jpi - 1 … … 454 506 455 507 CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. ) 508 CALL iom_put( 'zef' , zef ) ! MV DEBUG 509 IF( lwp ) WRITE(numout,*) ' outer loop 1c i_out : ', i_out 456 510 457 511 !--------------------------------------------------- 458 512 ! -- Ocean-ice drag and Coriolis RHS contributions 459 513 !--------------------------------------------------- 460 514 461 515 DO jj = 2, jpj - 1 462 516 DO ji = 2, jpi - 1 … … 489 543 END DO 490 544 END DO 545 546 IF( lwp ) WRITE(numout,*) ' outer loop 1d i_out : ', i_out 547 CALL lbc_lnk( 'icedyn_rhg_vp', zCwU, 'U', 1. ) 548 CALL lbc_lnk( 'icedyn_rhg_vp', zCwV, 'V', 1. ) 549 550 CALL iom_put( 'zCwU' , zCwU ) ! MV DEBUG 551 CALL iom_put( 'zCwV' , zCwV ) ! MV DEBUG 552 IF( lwp ) WRITE(numout,*) ' outer loop 1e i_out : ', i_out 553 CALL iom_put( 'zCorU' , zCorU ) ! MV DEBUG 554 CALL iom_put( 'zCorV' , zCorV ) ! MV DEBUG 555 IF( lwp ) WRITE(numout,*) ' outer loop 1f i_out : ', i_out 491 556 492 557 ! a priori, Coriolis and drag terms only affect diagonal or independent term of the linear system, … … 515 580 END DO 516 581 END DO 582 583 CALL iom_put( 'zs1_rhsu' , zs1_rhsu ) ! MV DEBUG 584 CALL iom_put( 'zs2_rhsu' , zs2_rhsu ) ! MV DEBUG 585 CALL iom_put( 'zs1_rhsv' , zs1_rhsv ) ! MV DEBUG 586 CALL iom_put( 'zs2_rhsv' , zs2_rhsv ) ! MV DEBUG 517 587 518 588 ! a priori, no lbclnk, because rhsu is only used in the inner domain … … 522 592 ! My guess is that this is the way to enforce boundary conditions on strain rate tensor 523 593 594 IF( lwp ) WRITE(numout,*) ' outer loop 2 i_out : ', i_out 595 524 596 DO jj = 1, jpj - 1 525 597 DO ji = 1, jpi - 1 … … 533 605 END DO 534 606 END DO 535 607 608 CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsu, 'F', 1. ) 609 CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsv, 'F', 1. ) 610 611 CALL iom_put( 'zs12_rhsu' , zs12_rhsu ) ! MV DEBUG 612 CALL iom_put( 'zs12_rhsv' , zs12_rhsv ) ! MV DEBUG 613 536 614 ! a priori, no lbclnk, because rhsu are only used in the inner domain 537 615 … … 555 633 END DO 556 634 END DO 635 636 CALL iom_put( 'zf_rhsu' , zf_rhsu ) ! MV DEBUG 637 CALL iom_put( 'zf_rhsv' , zf_rhsv ) ! MV DEBUG 557 638 558 639 !--------------------------- … … 571 652 END DO 572 653 654 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zrhsu, 'U', 1., zrhsv, 'V', 1. ) 655 656 CALL iom_put( 'zrhsu' , zrhsu ) ! MV DEBUG 657 CALL iom_put( 'zrhsv' , zrhsv ) ! MV DEBUG 658 573 659 ! inner domain calculations -> no lbclnk 660 661 IF( lwp ) WRITE(numout,*) ' outer loop 4 i_out : ', i_out 574 662 575 663 !---------------------------------------------------------------------------------------! … … 662 750 END DO 663 751 END DO 664 752 753 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zAU , 'U', 1., zAV , 'V', 1. ) 754 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zBU , 'U', 1., zBV , 'V', 1. ) 755 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCU , 'U', 1., zCV , 'V', 1. ) 756 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zDU , 'U', 1., zDV , 'V', 1. ) 757 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zEU , 'U', 1., zEV , 'V', 1. ) 758 759 CALL iom_put( 'zAU' , zAU ) ! MV DEBUG 760 CALL iom_put( 'zBU' , zBU ) ! MV DEBUG 761 CALL iom_put( 'zCU' , zCU ) ! MV DEBUG 762 CALL iom_put( 'zDU' , zDU ) ! MV DEBUG 763 CALL iom_put( 'zEU' , zEU ) ! MV DEBUG 764 CALL iom_put( 'zAV' , zAV ) ! MV DEBUG 765 CALL iom_put( 'zBV' , zBV ) ! MV DEBUG 766 CALL iom_put( 'zCV' , zCV ) ! MV DEBUG 767 CALL iom_put( 'zDV' , zDV ) ! MV DEBUG 768 CALL iom_put( 'zEV' , zEV ) ! MV DEBUG 769 665 770 !------------------------------------------------------------------------------! 666 771 ! … … 673 778 ll_v_iterate = .TRUE. 674 779 780 ! zBU 781 ! it is probably the way zBU is used that is problematic... 782 ! because zBU itself does not have the zipper!!! 783 ! BUG only occurs when zBU is different from zero 784 ! H1 - the jj loop does not work well with jpj = 149 ? 785 ! zBU(:,:) = 0._wp; zBV(:,:) = 0._wp 786 675 787 DO i_inn = 1, nn_ninn_vp ! inner loop iterations 676 677 ! mitgcm computes initial value of residual 678 jter = jter + 1 679 l_full_nf_update = jter == nn_nvp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 680 681 IF ( ll_u_iterate .OR. ll_v_iterate ) THEN 682 683 684 zu_b(:,:) = u_ice(:,:) ! velocity at previous sub-iterate 685 zv_b(:,:) = v_ice(:,:) 788 789 IF( lwp ) WRITE(numout,*) ' inner loop 1 i_inn : ', i_inn 790 791 !--- mitgcm computes initial value of residual here... 792 793 jter = jter + 1 794 ! l_full_nf_update = jter == nn_nvp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 795 796 zu_b(:,:) = u_ice(:,:) ! velocity at previous sub-iterate 797 zv_b(:,:) = v_ice(:,:) 798 799 zFU(:,:) = 0._wp ; zFU_prime(:,:) = 0._wp ; zBU_prime(:,:) = 0._wp 800 zFV(:,:) = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp 801 802 IF ( ll_u_iterate .OR. ll_v_iterate ) THEN 686 803 687 804 ! ---------------------------- ! … … 689 806 ! ---------------------------- ! 690 807 808 ! Thomas Algorithm for tridiagonal solver 691 809 ! what follows could be subroutinized... 692 810 … … 698 816 ELSE ; jj_min = 3 699 817 ENDIF 700 701 zFU(:,:) = 0._wp 702 zFU_prime(:,:) = 0._wp 703 zBU_prime(:,:) = 0._wp 704 818 819 IF ( lwp ) WRITE(numout,*) ' Into the U-zebra loop at step jn = ', jn, ', with jj_min = ', jj_min 820 821 !---------------------------------------------------------- 822 ! -- Boundary condition (link with neighbouring processor) 823 !---------------------------------------------------------- 824 705 825 DO jj = jj_min, jpj - 1, nn_zebra_vp 706 826 707 !-------------------------------------------708 ! -- Tridiagonal system solver for each row709 !-------------------------------------------710 827 ! 711 828 ! MV - I am in doubts whether the way I coded it is reproducible - ask Gurvan … … 722 839 IF ( ji == 2 ) zAA3 = zAA3 - zAU(ji,jj) * u_ice(ji-1,jj) 723 840 IF ( ji == jpi - 1 ) zAA3 = zAA3 - zCU(ji,jj) * u_ice(ji+1,jj) 724 841 725 842 ! right hand side 726 zFU(ji,jj) = ( zrhsu(ji,jj) & ! right-hand side terms727 & + zAA3 &! boundary condition translation728 & + zDU(ji,jj) * u_ice(ji,jj-1) &! internal force, j-1729 & + zEU(ji,jj) * u_ice(ji,jj+1) ) * umask(ji,jj,1) ! internal force, j+1730 843 zFU(ji,jj) = ( zrhsu(ji,jj) & ! right-hand side terms 844 & + zAA3 & ! boundary condition translation 845 & + zDU(ji,jj) * u_ice(ji,jj-1) & ! internal force, j-1 846 & + zEU(ji,jj) * u_ice(ji,jj+1) ) * umask(ji,jj,1) ! internal force, j+1 847 731 848 END DO 849 850 END DO 851 852 ! CALL lbc_lnk( 'icedyn_rhg_vp', zFU, 'U', 1. ) 732 853 733 ! - Thomas Algorithm 734 ! (MV: I chose a seemingly more efficient form of the algorithm than in mitGCM - not sure) 735 ! Forward sweep 736 DO ji = 3, jpi - 1 737 zw = zAU(ji,jj) / zBU(ji-1,jj) 738 zBU_prime(ji,jj) = zBU(ji,jj) - zw * zCU(ji-1,jj) 739 zFU_prime(ji,jj) = zFU(ji,jj) - zw * zFU(ji-1,jj) 854 !--------------- 855 ! Forward sweep 856 !--------------- 857 nn_thomas = 2 ! 1 = martin, 2 = mitgcm 858 859 IF ( nn_thomas == 1 ) THEN ! MV version seemingly more than mitGCM algorithm 860 861 DO jj = jj_min, jpj - 1, nn_zebra_vp 862 863 DO ji = 3, jpi - 1 864 865 zfac = SIGN( 1._wp , zBU(ji-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU(ji-1,jj) ) - epsi20 ) ) 866 zw = zfac * zAU(ji,jj) / MAX ( ABS( zBU(ji-1,jj) ) , epsi20 ) 867 zBU_prime(ji,jj) = zBU(ji,jj) - zw * zCU(ji-1,jj) 868 zFU_prime(ji,jj) = zFU(ji,jj) - zw * zFU(ji-1,jj) 869 870 END DO 871 740 872 END DO 741 742 ! Backward sweep 743 ! MV I don't see how this will be reproducible 744 u_ice(jpi-1,jj) = zFU_prime(jpi-1,jj) / zBU_prime(jpi-1,jj) * umask(jpi-1,jj,1) ! do last row first 745 DO ji = jpi-2 , 2, -1 ! all other rows ! MV OPT: could be turned into forward loop (by substituting ji) 746 u_ice(ji,jj) = zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji,jj+1) * umask(ji,jj,1) / zBU_prime(ji,jj) ! 747 END DO 748 749 !--------------- 750 ! -- Relaxation 751 !--------------- 873 874 ELSE ! nn_thomas == 2 ! mitGCM algorithm 875 876 DO jj = jj_min, jpj - 1, nn_zebra_vp 877 878 ! ji = 2 879 zw = zBU(2,jj) 880 zfac = SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) ) 881 zw = zfac * 1.0_wp / MAX ( ABS( zw ) , epsi20 ) 882 zCU_prime(2,jj) = zw * zCU(2,jj) 883 zFU_prime(2,jj) = zw * zFU(2,jj) 884 885 DO ji = 3, jpi - 1 886 887 zw = zBU(ji,jj) - zAU(ji,jj) * zCU_prime(ji-1,jj) 888 zfac = SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) ) 889 zw = zfac * 1.0_wp / MAX ( ABS( zw ) , epsi20 ) 890 zCU_prime(ji,jj) = zw * zCU(ji,jj) 891 zFU_prime(ji,jj) = zw * ( zFU(ji,jj) - zAU(ji,jj) * zFU_prime(ji-1,jj) ) 892 893 END DO 894 895 END DO 896 897 ENDIF 898 899 ! CALL lbc_lnk_multi( 'icedyn_rhg_vp', zBU_prime , 'U', 1., zFU_prime, 'U', 1. ) 900 901 !---------------- 902 ! Backward sweep 903 !---------------- 904 IF ( nn_thomas == 1 ) THEN ! MV version 905 906 DO jj = jj_min, jpj - 1, nn_zebra_vp 907 908 ! last row 909 zfac = SIGN( 1._wp , zBU_prime(jpi-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(jpi-1,jj) ) - epsi20 ) ) 910 u_ice(jpi-1,jj) = zfac * zFU_prime(jpi-1,jj) / MAX( ABS ( zBU_prime(jpi-1,jj) ) , epsi20 ) & 911 & * umask(jpi-1,jj,1) 912 913 DO ji2 = 2 , jpi-2 ! all other rows 914 !DO ji = jpi-2 , 2, -1 ! all other rows ! ---> original backward loop 915 ji = jpi - ji2 916 ! ji2 = 2 -> ji = jpi - 2 917 ! ji2 = jpi - 2 -> ji = 2 918 zfac = SIGN( 1._wp , zBU_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(ji,jj) ) - epsi20 ) ) 919 u_ice(ji,jj) = zfac * ( zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji+1,jj) * umask(ji,jj,1) ) & ! this 920 ! line is guilty 921 & / MAX ( ABS ( zBU_prime(ji,jj) ) , epsi20 ) 922 END DO 923 924 END DO 925 926 ELSE ! nn_thomas == 2 ! mitGCM version 927 928 DO jj = jj_min, jpj - 1, nn_zebra_vp 929 930 u_ice(jpi-1,jpj) = zFU_prime(jpi-1,jpj) 931 932 DO ji2 = 2 , jpi-2 ! all other rows ! MV OPT: could be turned into forward loop (by substituting ji) 933 ji = jpi - ji2 934 u_ice(ji,jj) = zFU_prime(ji,jj) - zCU_prime(ji,jj) * u_ice(ji+1,jj) 935 END DO 936 937 END DO 938 939 ENDIF 940 941 !--------------- 942 ! -- Relaxation 943 !--------------- 944 945 DO jj = jj_min, jpj - 1, nn_zebra_vp 946 752 947 DO ji = 2, jpi - 1 753 948 … … 755 950 756 951 ! velocity masking for little-ice and no-ice cases 757 ! put 0.01 of ocean current, appropriate choice to avoid too much spreading of the ice edge758 952 u_ice(ji,jj) = zmsk00x(ji,jj) & 759 953 & * ( zmsk01x(ji,jj) * u_ice(ji,jj) & 760 + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp)954 & + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp ) * umask(ji,jj,1) 761 955 762 956 END DO 763 957 764 958 END DO ! jj 765 959 766 960 END DO ! zebra loop 767 961 768 962 ENDIF ! ll_u_iterate 769 963 770 964 ! ! ---------------------------- ! 771 965 IF ( ll_v_iterate ) THEN ! --- Solve for V-velocity --- ! 772 !! ---------------------------- !966 ! ! ---------------------------- ! 773 967 774 968 ! MV OPT: what follows could be subroutinized... … … 779 973 ELSE ; ji_min = 3 780 974 ENDIF 781 782 zFV(:,:) = 0._wp 783 zFV_prime(:,:) = 0._wp 784 zBV_prime(:,:) = 0._wp 785 975 976 IF ( lwp ) WRITE(numout,*) ' Into the V-zebra loop at step jn = ', jn, ', with ji_min = ', ji_min 977 786 978 DO ji = ji_min, jpi - 1, nn_zebra_vp 787 979 … … 800 992 ! see Zhang and Hibler, 1997, Appendix B 801 993 zAA3 = 0._wp 802 IF ( jj .EQ.2 ) zAA3 = zAA3 - zAV(ji,jj) * v_ice(ji,jj-1)803 IF ( jj .EQ.jpj - 1 ) zAA3 = zAA3 - zCV(ji,jj) * v_ice(ji,jj+1)994 IF ( jj == 2 ) zAA3 = zAA3 - zAV(ji,jj) * v_ice(ji,jj-1) 995 IF ( jj == jpj - 1 ) zAA3 = zAA3 - zCV(ji,jj) * v_ice(ji,jj+1) 804 996 805 997 ! right hand side … … 808 1000 & + zDV(ji,jj) * v_ice(ji-1,jj) & ! internal force, j-1 809 1001 & + zEV(ji,jj) * v_ice(ji+1,jj) ) * vmask(ji,jj,1) ! internal force, j+1 810 1002 811 1003 END DO 812 1004 813 1005 ! --- Thomas Algorithm 814 1006 ! (MV: I chose a seemingly more efficient form of the algorithm than in mitGCM - not sure) 815 1007 ! Forward sweep 816 1008 DO jj = 3, jpj - 1 817 zw = zAV(ji,jj) / zBV(ji,jj-1) 1009 zfac = SIGN( 1._wp , zBV(ji,jj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj-1) ) - epsi20 ) ) 1010 zw = zfac * zAV(ji,jj) / MAX ( ABS( zBV(ji,jj-1) ) , epsi20 ) 818 1011 zBV_prime(ji,jj) = zBV(ji,jj) - zw * zCV(ji,jj-1) 819 1012 zFV_prime(ji,jj) = zFV(ji,jj) - zw * zFV(ji,jj-1) … … 821 1014 822 1015 ! Backward sweep 823 v_ice(ji,jpj-1) = zFV_prime(ji,jpj-1) / zBV_prime(ji,jpj-1) * vmask(ji,jj,jpj-1) ! last row 1016 zfac = SIGN( 1._wp , zBV_prime(ji,jpj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jpj-1) ) - epsi20 ) ) 1017 v_ice(ji,jpj-1) = zfac * zFV_prime(ji,jpj-1) / MAX ( ABS(zBV_prime(ji,jpj-1) ) , epsi20 ) & 1018 & * vmask(ji,jpj-1,1) ! last row 1019 824 1020 DO jj = jpj-2, 2, -1 ! can be turned into forward row by substituting jj if needed 825 v_ice(ji,jj) = zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) * vmask(ji,jj,1) / zBV_prime(ji,jj) 1021 zfac = SIGN( 1._wp , zBV_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jj) ) - epsi20 ) ) 1022 v_ice(ji,jj) = zfac * ( zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) * vmask(ji,jj,1) ) & 1023 & / MAX ( ABS( zBV_prime(ji,jj) ) , epsi20 ) 826 1024 END DO 827 1025 828 1026 !--------------- 829 1027 ! -- Relaxation … … 835 1033 v_ice(ji,jj) = zmsk00y(ji,jj) & 836 1034 & * ( zmsk01y(ji,jj) * v_ice(ji,jj) & 837 + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp)1035 & + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp ) * vmask(ji,jj,1) 838 1036 839 1037 END DO … … 844 1042 845 1043 ENDIF ! ll_v_iterate 846 847 IF ( ( ll_u_iterate .OR. ll_v_iterate ) .OR. jter == nn_nvp )CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. )1044 1045 CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 848 1046 849 1047 !-------------------------------------------------------------------------------------- … … 856 1054 ! MV OPT: if the number of iterations to convergence is really variable, and keep the convergence check 857 1055 ! then we must optimize the use of the mpp_max, which is prohibitive 1056 zuerr_max = 0._wp 858 1057 859 1058 IF ( ll_u_iterate .AND. MOD ( i_inn, nn_cvgchk_vp ) == 0 ) THEN … … 886 1085 ! on V 887 1086 !------ 1087 zverr_max = 0._wp 888 1088 889 1089 IF ( ll_v_iterate .AND. MOD ( i_inn, nn_cvgchk_vp ) == 0 ) THEN … … 919 1119 ! 920 1120 !--------------------------------------------------------------------------------------- 921 IF( ln_rhg_chkcvg ) CALL rhg_cvg_vp( kt, jter, nn_nvp, u_ice, v_ice, zmt, zuerr_max, zverr_max, zglob_area, &1121 IF( ln_rhg_chkcvg .AND. MOD ( i_inn - 1, nn_cvgchk_vp ) == 0 ) CALL rhg_cvg_vp( kt, jter, nn_nvp, u_ice, v_ice, zmt, zuerr_max, zverr_max, zglob_area, & 922 1122 & zrhsu, zAU, zBU, zCU, zDU, zEU, zrhsv, zAV, zBV, zCV, zDV, zEV ) 923 924 1123 1124 IF ( lwp ) WRITE(numout,*) ' Done convergence tests ' 1125 925 1126 ENDIF ! --- end ll_u_iterate or ll_v_iterate 926 1127 927 1128 END DO ! i_inn, end of inner loop 928 1129 1130 END DO ! End of outer loop (i_out) ============================================================================================= 1131 1132 IF ( lwp ) WRITE(numout,*) ' We are out of outer loop ' 1133 1134 CALL iom_put( 'zFU' , zFU ) ! MV DEBUG 1135 CALL iom_put( 'zBU_prime' , zBU_prime ) ! MV DEBUG 1136 CALL iom_put( 'zFU_prime' , zFU_prime ) ! MV DEBUG 1137 1138 CALL iom_put( 'zFV' , zFV ) ! MV DEBUG 1139 CALL iom_put( 'zBV_prime' , zBV_prime ) ! MV DEBUG 1140 CALL iom_put( 'zFV_prime' , zFV_prime ) ! MV DEBUG 1141 1142 CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 1143 1144 IF ( lwp ) WRITE(numout,*) ' We are about to output uice_dbg ' 1145 IF( iom_use('uice_dbg' ) ) CALL iom_put( 'uice_dbg' , u_ice ) ! ice velocity u after solver 1146 IF( iom_use('vice_dbg' ) ) CALL iom_put( 'vice_dbg' , v_ice ) ! ice velocity v after solver 1147 929 1148 !------------------------------------------------------------------------------! 930 1149 ! 931 ! 6) Mask final velocities1150 ! --- Convergence diagnostics 932 1151 ! 933 1152 !------------------------------------------------------------------------------! 934 1153 935 END DO ! End of outer loop (i_out) ============================================================================================= 1154 IF( ln_rhg_chkcvg ) THEN 1155 1156 IF( iom_use('uice_cvg') ) THEN 1157 CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_b(:,:) ) * umask(:,:,1) , & ! ice velocity difference at last iteration 1158 & ABS( v_ice(:,:) - zv_b(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 1159 ENDIF 1160 1161 ENDIF ! ln_rhg_chkcvg 1162 1163 ! MV DEBUG test - replace ice velocity by ocean current to give the model the means to go ahead 1164 DO jj = 2, jpj - 1 1165 DO ji = 2, jpi - 1 1166 1167 u_ice(ji,jj) = zmsk00x(ji,jj) & 1168 & * ( zmsk01x(ji,jj) * u_oce(ji,jj) * 0.01_wp & 1169 + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp ) 1170 1171 v_ice(ji,jj) = zmsk00y(ji,jj) & 1172 & * ( zmsk01y(ji,jj) * v_oce(ji,jj) * 0.01_wp & 1173 + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp ) 1174 1175 END DO 1176 END DO 1177 1178 CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 1179 1180 ! END DEBUG 936 1181 937 1182 !------------------------------------------------------------------------------! … … 942 1187 ! 943 1188 ! MV OPT: subroutinize ? 944 1189 945 1190 DO jj = 1, jpj - 1 946 1191 DO ji = 1, jpi - 1 … … 1087 1332 IF( iom_use('icediv') ) CALL iom_put( 'icediv' , pdivu_i * zmsk00 ) ! divergence 1088 1333 IF( iom_use('iceshe') ) CALL iom_put( 'iceshe' , pshear_i * zmsk00 ) ! shear 1334 IF( iom_use('icedlt') ) CALL iom_put( 'icedlt' , pdelta_i * zmsk00 ) ! delta 1089 1335 IF( iom_use('icestr') ) CALL iom_put( 'icestr' , strength * zmsk00 ) ! strength 1090 1336 … … 1123 1369 ! Recommendation 2 : Need to use deformations at PREVIOUS iterate for viscosities (see p. 1765) 1124 1370 ! R2 means we need to recompute stresses 1371 1125 1372 IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 1126 1373 ! 1127 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) )1374 ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 1128 1375 ! 1129 1376 DO jj = 2, jpj - 1 … … 1142 1389 1143 1390 ! Normalized principal stresses (used to display the ellipse) 1144 z1_strength = 1._wp / strength(ji,jj)1145 zsig1_p(ji,jj) = ( zsig_I(ji,jj) 1146 zsig2_p(ji,jj) = ( zsig_I I(ji,jj) - zsig_II(ji,jj) ) * z1_strength1391 z1_strength = 1._wp / MAX ( 1._wp , strength(ji,jj) ) 1392 zsig1_p(ji,jj) = ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 1393 zsig2_p(ji,jj) = ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 1147 1394 END DO 1148 1395 END DO … … 1154 1401 CALL iom_put( 'sig2_pnorm' , zsig2_p ) 1155 1402 1156 DEALLOCATE( zsig1_p , zsig2_p )1403 DEALLOCATE( zsig1_p , zsig2_p , zsig_I , zsig_II ) 1157 1404 1158 1405 ENDIF … … 1253 1500 ENDIF 1254 1501 1255 ! --- Convergence diagnostics --- !1256 IF( ln_rhg_chkcvg ) THEN1257 1258 IF( iom_use('uice_cvg') ) THEN1259 CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_b(:,:) ) * umask(:,:,1) , &1260 & ABS( v_ice(:,:) - zv_b(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )1261 ENDIF1262 1263 ENDIF ! ln_rhg_chkcvg1264 1265 !1266 1502 DEALLOCATE( zmsk00, zmsk15 ) 1267 1503 … … 1285 1521 !! - residuals in U, V and UV-mean taken as square-root of area-weighted mean square residual (u_res, v_res, vel_res, N/m2) 1286 1522 !! - mean kinetic energy (mke_ice, J/m2) 1287 !!1288 1523 !! 1289 1524 !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) … … 1296 1531 REAL(wp), DIMENSION(:,:), INTENT(in) :: prhsv, pAV, pBV, pCV, pDV, pEV 1297 1532 !! 1298 INTEGER :: it, idtime, istatus 1533 INTEGER :: it, idtime, istatus, ix_dim, iy_dim 1299 1534 INTEGER :: ji, jj ! dummy loop indices 1300 1535 REAL(wp) :: zveldif, zu_res_mean, zv_res_mean, zvelres, zmke, zu, zv ! local scalars 1536 REAL(wp) :: z1_pglob_area 1301 1537 REAL(wp), DIMENSION(jpi,jpj) :: zu_res, zv_res, zvel2 ! local arrays 1302 1538 … … 1304 1540 !!---------------------------------------------------------------------- 1305 1541 1542 IF( lwp ) THEN 1543 WRITE(numout,*) 1544 WRITE(numout,*) 'rhg_cvg_vp : ice rheology convergence control' 1545 WRITE(numout,*) '~~~~~~~~~~~' 1546 WRITE(numout,*) ' kiter = : ', kiter 1547 WRITE(numout,*) ' kitermax = : ', kitermax 1548 ENDIF 1549 1306 1550 ! create file 1307 1551 IF( kt == nit000 .AND. kiter == 1 ) THEN 1308 1552 ! 1309 IF( lwp ) THEN1310 WRITE(numout,*)1311 WRITE(numout,*) 'rhg_cvg_vp : ice rheology convergence control'1312 WRITE(numout,*) '~~~~~~~~~~~'1313 ENDIF1314 !1315 1553 IF( lwm ) THEN 1554 1316 1555 clname = 'ice_cvg.nc' 1317 1556 IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 1318 1557 istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) 1558 1319 1559 istatus = NF90_DEF_DIM( ncvgid, 'time' , NF90_UNLIMITED, idtime ) 1320 istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid_ucvg ) 1560 istatus = NF90_DEF_DIM( ncvgid, 'x' , jpi, ix_dim ) 1561 istatus = NF90_DEF_DIM( ncvgid, 'y' , jpj, iy_dim ) 1562 1563 istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid_ucvg ) ! the name uice_cvg sucks, no ? 1564 ! i suggest vel_dif instead 1321 1565 istatus = NF90_DEF_VAR( ncvgid, 'u_res', NF90_DOUBLE , (/ idtime /), nvarid_ures ) 1322 1566 istatus = NF90_DEF_VAR( ncvgid, 'v_res', NF90_DOUBLE , (/ idtime /), nvarid_vres ) … … 1325 1569 istatus = NF90_DEF_VAR( ncvgid, 'v_dif', NF90_DOUBLE , (/ idtime /), nvarid_vdif ) 1326 1570 istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE , (/ idtime /), nvarid_mke ) 1571 1572 istatus = NF90_DEF_VAR( ncvgid, 'usub' , NF90_DOUBLE, (/ ix_dim, iy_dim, idtime /), nvarid_usub) ! --> u-velocity in sub-iterations 1573 istatus = NF90_DEF_VAR( ncvgid, 'vsub' , NF90_DOUBLE, (/ ix_dim, iy_dim, idtime /), nvarid_vsub) ! --> v-velocity in sub-iterations 1574 1327 1575 istatus = NF90_ENDDEF(ncvgid) 1576 1328 1577 ENDIF 1329 1578 ! 1330 1579 ENDIF 1331 1580 1332 ! time 1333 it = ( kt - 1 ) * kitermax + kiter 1581 IF ( lwp ) WRITE(numout,*) ' File created ' 1334 1582 1335 1583 ! --- Max absolute velocity difference with previous iterate (zveldif) 1336 ! EVP code1337 ! IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large)1338 ! zveldif = 0._wp1339 ! ELSE1340 ! DO jj = 1, jpj1341 ! DO ji = 1, jpi1342 ! zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &1343 ! & ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj)1344 ! END DO1345 ! END DO1346 ! zveldif = MAXVAL( zres )1347 ! CALL mpp_max( 'icedyn_rhg_vp', zveldif ) ! max over the global domain1348 ! ENDIF1349 ! VP code1350 1584 zveldif = MAX( puerr_max, pverr_max ) ! velocity difference with previous iterate, should nearly be equivalent to evp code 1351 1585 ! if puerrmask and pverrmax are masked at 15% (TEST) 1352 1586 1587 IF ( lwp ) WRITE(numout,*) ' zfeldif : ', zveldif 1588 1589 ! --- Mean residual and kinetic energy 1590 IF ( kiter == 1 ) THEN 1591 1592 zu_res_mean = 0._wp 1593 zv_res_mean = 0._wp 1594 zvelres = 0._wp 1595 zmke = 0._wp 1596 1597 ELSE 1598 1353 1599 ! -- Mean residual (N/m^2), zu_res_mean 1354 1600 ! Here we take the residual of the linear system (N/m^2), … … 1356 1602 ! Local residual r = Ax - B expresses to which extent the momentum balance is verified 1357 1603 ! i.e., how close we are to a solution 1604 1605 IF ( lwp ) WRITE(numout,*) ' TEST 1 ' 1606 1607 z1_pglob_area = 1._wp / pglob_area 1608 1358 1609 DO jj = 2, jpj - 1 1359 1610 DO ji = 2, jpi - 1 1360 zu_res(ji,jj) = prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1)&1361 & - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj)1611 zu_res(ji,jj) = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1) & 1612 & - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) ) 1362 1613 1363 zv_res(ji,jj) = prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj) & 1364 & - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) 1614 zv_res(ji,jj) = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj) & 1615 & - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) ) 1616 1617 zu_res(ji,jj) = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) * e1e2u(ji,jj) * z1_pglob_area 1618 zv_res(ji,jj) = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) * e1e2v(ji,jj) * z1_pglob_area 1619 1365 1620 END DO 1366 1621 END DO 1367 zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) * zu_res(:,:) * e1e2u(:,:) * umask(:,:,1) ) 1368 zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) * zv_res(:,:) * e1e2v(:,:) * vmask(:,:,1) ) 1369 zu_res_mean = SQRT( zu_res_mean / pglob_area ) 1370 zv_res_mean = SQRT( zv_res_mean / pglob_area ) 1622 1623 IF ( lwp ) WRITE(numout,*) ' TEST 2 ' 1624 zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) ) 1625 zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) ) 1626 IF ( lwp ) WRITE(numout,*) ' TEST 3 ' 1371 1627 zvelres = 0.5_wp * ( zu_res_mean + zv_res_mean ) 1628 1629 IF ( lwp ) WRITE(numout,*) ' TEST 4 ' 1372 1630 1373 1631 ! -- Global mean kinetic energy per unit area (J/m2) … … 1380 1638 END DO 1381 1639 1382 zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) 1640 IF ( lwp ) WRITE(numout,*) ' TEST 5 ' 1641 zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) / pglob_area 1642 IF ( lwp ) WRITE(numout,*) ' TEST 6 ' 1643 1644 ENDIF 1383 1645 1384 1646 ! ! ==================== ! 1385 1647 1648 ! time 1649 it = ( kt - 1 ) * kitermax + kiter 1650 1651 1386 1652 IF( lwm ) THEN 1387 1653 ! write variables 1388 istatus = NF90_PUT_VAR( ncvgid, nvarid_ucvg, (/zveldif/), (/it/), (/1/) ) 1389 istatus = NF90_PUT_VAR( ncvgid, nvarid_ures, (/zu_res_mean/), (/it/), (/1/) ) 1390 istatus = NF90_PUT_VAR( ncvgid, nvarid_vres, (/zv_res_mean/), (/it/), (/1/) ) 1391 istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvelres/), (/it/), (/1/) ) 1392 istatus = NF90_PUT_VAR( ncvgid, nvarid_udif, (/puerr_max/), (/it/), (/1/) ) 1393 istatus = NF90_PUT_VAR( ncvgid, nvarid_vdif, (/pverr_max/), (/it/), (/1/) ) 1394 istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/it/), (/1/) ) 1654 istatus = NF90_PUT_VAR( ncvgid, nvarid_ucvg, (/zveldif/), (/it/), (/1/) ) ! max U or V velocity diff between subiterations 1655 istatus = NF90_PUT_VAR( ncvgid, nvarid_ures, (/zu_res_mean/), (/it/), (/1/) ) ! U-residual of the linear system 1656 istatus = NF90_PUT_VAR( ncvgid, nvarid_vres, (/zv_res_mean/), (/it/), (/1/) ) ! V-residual of the linear system 1657 istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvelres/), (/it/), (/1/) ) ! average of u- and v- residuals 1658 istatus = NF90_PUT_VAR( ncvgid, nvarid_udif, (/puerr_max/), (/it/), (/1/) ) ! max U velocity difference, inner iterations 1659 istatus = NF90_PUT_VAR( ncvgid, nvarid_vdif, (/pverr_max/), (/it/), (/1/) ) ! max V velocity difference, inner iterations 1660 istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/it/), (/1/) ) ! mean kinetic energy 1661 1662 istatus = NF90_PUT_VAR( ncvgid, nvarid_usub, (/pu/), (/it/) ) ! u-velocity 1663 istatus = NF90_PUT_VAR( ncvgid, nvarid_vsub, (/pv/), (/it/) ) ! v-velocity 1664 1395 1665 ! close file 1396 1666 IF( kt == nitend ) istatus = NF90_CLOSE( ncvgid )
Note: See TracChangeset
for help on using the changeset viewer.