Changeset 13681 for NEMO/branches
- Timestamp:
- 2020-10-27T13:54:43+01: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_vp.F90
r13629 r13681 43 43 !! for convergence tests 44 44 INTEGER :: ncvgid ! netcdf file id 45 INTEGER :: nvarid_ucvg ! netcdf variable id46 45 INTEGER :: nvarid_ures 47 46 INTEGER :: nvarid_vres … … 49 48 INTEGER :: nvarid_udif 50 49 INTEGER :: nvarid_vdif 50 INTEGER :: nvarid_veldif 51 51 INTEGER :: nvarid_mke 52 INTEGER :: nvarid_u sub, nvarid_vsub52 INTEGER :: nvarid_ures_xy, nvarid_vres_xy 53 53 54 54 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmsk00, zmsk15 … … 136 136 INTEGER :: nn_zebra_vp ! number of zebra steps 137 137 138 INTEGER :: nn_thomas ! tridiagonal system version (1=mitgcm, 2=martin)139 140 138 INTEGER :: nn_nvp ! total number of VP iterations (n_out_vp*n_inn_vp) 141 139 ! … … 144 142 REAL(wp) :: zglob_area ! global ice area for diagnostics 145 143 REAL(wp) :: zkt ! isotropic tensile strength for landfast ice 146 REAL(wp) :: zm 2, zm3, zmassU, zmassV! ice/snow mass and volume144 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV ! ice/snow mass and volume 147 145 REAL(wp) :: zdeltat, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 148 146 REAL(wp) :: zp_deltastar_f ! … … 294 292 END DO 295 293 END DO 296 CALL lbc_lnk( 'icedyn_rhg_vp', zfmask, 'F', 1._wp )297 294 298 295 ! Lateral boundary conditions on velocity (modify zfmask) … … 335 332 zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 336 333 zsshdyn(:,:) = 0._wp ! DEBUG CAREFUL !!! 334 335 zmt(:,:) = rhos * vt_s(:,:) + rhoi * vt_i(:,:) ! Snow and ice mass at T-point 336 zmf(:,:) = zmt(:,:) * ff_t(:,:) ! Coriolis factor at T points (m*f) 337 337 338 338 DO jj = 2, jpj - 1 … … 344 344 345 345 ! Snow and ice mass at U-V points 346 zm t(ji,jj) = ( rhos * vt_s(ji ,jj ) + rhoi * vt_i(ji ,jj ) )347 zm2 = ( rhos * vt_s(ji+1,jj ) + rhoi * vt_i(ji+1,jj ))348 zm3 = ( rhos * vt_s(ji ,jj+1) + rhoi * vt_i(ji ,jj+1))346 zm1 = zmt(ji,jj) 347 zm2 = zmt(ji+1,jj) 348 zm3 = zmt(ji,jj+1) 349 349 350 zmassU = 0.5_wp * ( zm t(ji,jj)* e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)351 zmassV = 0.5_wp * ( zm t(ji,jj)* e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)350 zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 351 zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 352 352 353 353 ! Mass per unit area divided by time step … … 362 362 v_oceU(ji,jj) = 0.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1) 363 363 u_oceV(ji,jj) = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1) 364 365 ! Coriolis factor at T points (m*f)366 zmf(ji,jj) = zmt(ji,jj) * ff_t(ji,jj)367 364 368 365 ! Wind stress … … 379 376 380 377 ! Mask for lots of ice (1) or little ice (0) 381 IF( zmassU <= zmmin .AND. za_iU(ji,jj) <= zamin ) THEN ; zmsk01x(ji,jj) = 0._wp 382 ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF 383 IF( zmassV <= zmmin .AND. za_iV(ji,jj) <= zamin ) THEN ; zmsk01y(ji,jj) = 0._wp 384 ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF 378 IF ( zmassU <= zmmin .AND. za_iU(ji,jj) <= zamin ) THEN ; zmsk01x(ji,jj) = 0._wp 379 ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF 380 IF ( zmassV <= zmmin .AND. za_iV(ji,jj) <= zamin ) THEN ; zmsk01y(ji,jj) = 0._wp 381 ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF 382 383 ! MV TEST DEBUG 384 IF ( ( zmt(ji,jj) <= zmmin .OR. zmt(ji+1,jj) <= zmmin ) .AND. & 385 & ( at_i(ji,jj) <= zamin .OR. at_i(ji+1,jj) <= zamin ) ) THEN ; zmsk01x(ji,jj) = 0._wp 386 ELSE ; zmsk01x(ji,jj) = 1._wp ; ENDIF 387 388 IF ( ( zmt(ji,jj) <= zmmin .OR. zmt(ji,jj+1) <= zmmin ) .AND. & 389 & ( at_i(ji,jj) <= zamin .OR. at_i(ji,jj+1) <= zamin ) ) THEN ; zmsk01y(ji,jj) = 0._wp 390 ELSE ; zmsk01y(ji,jj) = 1._wp ; ENDIF 391 ! END MV TEST DEBUG 385 392 386 393 END DO … … 391 398 CALL iom_put( 'zmsk01x' , zmsk01x ) ! MV DEBUG 392 399 CALL iom_put( 'zmsk01y' , zmsk01y ) ! MV DEBUG 393 CALL iom_put( 'ztauy_ai' , ztauy_ai ) ! MV DEBUG394 400 CALL iom_put( 'ztaux_ai' , ztaux_ai ) ! MV DEBUG 395 401 CALL iom_put( 'ztauy_ai' , ztauy_ai ) ! MV DEBUG … … 449 455 IF( lwp ) WRITE(numout,*) ' outer loop 1a i_out : ', i_out 450 456 451 DO jj = 2, jpj - 1 ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 452 DO ji = 2, jpi - 1 ! 457 !DO jj = 2, jpj - 1 ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 458 ! DO ji = 2, jpi - 1 ! 459 460 ! MV DEBUG 461 DO jj = 2, jpj ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 462 DO ji = 2, jpi ! 463 ! END MV DEBUG 453 464 454 465 ! shear**2 at T points (doc eq. A16) … … 500 511 501 512 ! Temporary zef factor at F-point 502 zef(ji,jj) = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2 513 zef(ji,jj) = zp_deltastar_f * r1_e1e2f(ji,jj) * z1_ecc2 * zfmask(ji,jj) 503 514 504 515 END DO … … 545 556 546 557 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. ) 558 559 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCwU , 'U', -1., zCwV, 'V', -1. ) 560 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCorU, 'U', -1., zCorV, 'V', -1. ) 549 561 550 562 CALL iom_put( 'zCwU' , zCwU ) ! MV DEBUG 551 563 CALL iom_put( 'zCwV' , zCwV ) ! MV DEBUG 552 IF( lwp ) WRITE(numout,*) ' outer loop 1e i_out : ', i_out553 564 CALL iom_put( 'zCorU' , zCorU ) ! MV DEBUG 554 565 CALL iom_put( 'zCorV' , zCorV ) ! MV DEBUG 566 555 567 IF( lwp ) WRITE(numout,*) ' outer loop 1f i_out : ', i_out 556 568 … … 652 664 END DO 653 665 654 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zrhsu, 'U', 1., zrhsv, 'V', 1. ) 655 666 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zrhsu, 'U', -1., zrhsv, 'V', -1.) 667 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zmU_t, 'U', -1., zmV_t, 'V', -1.) 668 CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi_rhsu, 'U', -1., ztauy_oi_rhsv, 'V', -1.) 669 670 CALL iom_put( 'zmU_t' , zmU_t ) ! MV DEBUG 671 CALL iom_put( 'zmV_t' , zmV_t ) ! MV DEBUG 672 CALL iom_put( 'ztaux_oi_rhsu' , ztaux_oi_rhsu ) ! MV DEBUG 673 CALL iom_put( 'ztauy_oi_rhsv' , ztauy_oi_rhsv ) ! MV DEBUG 656 674 CALL iom_put( 'zrhsu' , zrhsu ) ! MV DEBUG 657 675 CALL iom_put( 'zrhsv' , zrhsv ) ! MV DEBUG … … 778 796 ll_v_iterate = .TRUE. 779 797 780 ! nn_thomas = 2 ! Type of Thomas algorithm for tridiagonal system. 1 = martin, 2 = mitgcm ! ---> not working for now781 nn_thomas = 1 ! Type of Thomas algorithm for tridiagonal system. 1 = martin, 2 = mitgcm ! ---> BUG MPP!!!782 783 798 DO i_inn = 1, nn_ninn_vp ! inner loop iterations 784 799 … … 840 855 END DO 841 856 842 !--------------- 843 ! Forward sweep 844 !--------------- 845 846 IF ( nn_thomas == 1 ) THEN ! MV version seemingly more than mitGCM algorithm 847 848 DO ji = 3, jpi - 1 849 850 zfac = SIGN( 1._wp , zBU(ji-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU(ji-1,jj) ) - epsi20 ) ) 851 zw = zfac * zAU(ji,jj) / MAX ( ABS( zBU(ji-1,jj) ) , epsi20 ) 852 zBU_prime(ji,jj) = zBU(ji,jj) - zw * zCU(ji-1,jj) 853 zFU_prime(ji,jj) = zFU(ji,jj) - zw * zFU(ji-1,jj) 854 855 END DO 856 857 ELSE ! nn_thomas == 2 ! mitGCM algorithm 858 859 ! --- ji = 2 860 zw = zBU(2,jj) 861 zfac = SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) ) 862 zw = zfac / MAX ( ABS( zw ) , epsi20 ) 863 zCU_prime(2,jj) = zw * zCU(2,jj) 864 zFU_prime(2,jj) = zw * zFU(2,jj) 865 866 ! --- ji = 3 --> jpi-1 867 DO ji = 3, jpi - 1 868 869 zw = zBU(ji,jj) - zAU(ji,jj) * zCU_prime(ji-1,jj) 870 zfac = SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) ) 871 zw = zfac / MAX ( ABS( zw ) , epsi20 ) 872 zCU_prime(ji,jj) = zw * zCU(ji,jj) 873 zFU_prime(ji,jj) = zw * ( zFU(ji,jj) - zAU(ji,jj) * zFU_prime(ji-1,jj) ) 874 875 END DO 876 877 ENDIF 857 END DO 858 859 CALL lbc_lnk( 'icedyn_rhg_vp', zFU, 'U', 1. ) 860 861 !--------------- 862 ! Forward sweep 863 !--------------- 864 DO jj = jj_min, jpj - 1, nn_zebra_vp 865 866 DO ji = 3, jpi - 1 867 868 zfac = SIGN( 1._wp , zBU(ji-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU(ji-1,jj) ) - epsi20 ) ) 869 zw = zfac * zAU(ji,jj) / MAX ( ABS( zBU(ji-1,jj) ) , epsi20 ) 870 zBU_prime(ji,jj) = zBU(ji,jj) - zw * zCU(ji-1,jj) 871 zFU_prime(ji,jj) = zFU(ji,jj) - zw * zFU(ji-1,jj) 872 873 END DO 874 875 END DO 876 877 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFU_prime, 'U', 1., zBU_prime, 'U', 1. ) 878 878 879 !----------------------------- 880 ! Backward sweep & relaxation 881 !----------------------------- 879 !----------------------------- 880 ! Backward sweep & relaxation 881 !----------------------------- 882 883 DO jj = jj_min, jpj - 1, nn_zebra_vp 882 884 883 885 ! --- Backward sweep 884 IF ( nn_thomas == 1 ) THEN ! MV version 885 886 ! last row 887 zfac = SIGN( 1._wp , zBU_prime(jpi-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(jpi-1,jj) ) - epsi20 ) ) 888 u_ice(jpi-1,jj) = zfac * zFU_prime(jpi-1,jj) / MAX( ABS ( zBU_prime(jpi-1,jj) ) , epsi20 ) & 889 & * umask(jpi-1,jj,1) 890 891 DO ji2 = 2 , jpi-2 ! all other rows 892 !DO ji = jpi-2 , 2, -1 ! all other rows ! ---> original backward loop 893 ji = jpi - ji2 894 ! ji2 = 2 -> ji = jpi - 2 895 ! ji2 = jpi - 2 -> ji = 2 896 zfac = SIGN( 1._wp , zBU_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(ji,jj) ) - epsi20 ) ) 897 u_ice(ji,jj) = zfac * ( zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji+1,jj) ) * umask(ji,jj,1) & ! this 898 ! line is guilty 899 & / MAX ( ABS ( zBU_prime(ji,jj) ) , epsi20 ) 900 END DO 901 902 ELSE ! nn_thomas == 2 ! mitGCM version 903 904 u_ice(jpi-1,jj) = zFU_prime(jpi-1,jj) * umask(jpi-1,jj,1) 905 906 DO ji2 = 2 , jpi-2 ! all other rows ! MV OPT: could be turned into forward loop (by substituting ji) 907 ji = jpi - ji2 908 ! ji2 = 2, ji = jpi - 2 909 ! ji2 = jpi- 2, ji = 2 910 u_ice(ji,jj) = ( zFU_prime(ji,jj) - zCU_prime(ji,jj) * u_ice(ji+1,jj) ) * umask(ji,jj,1) 911 END DO 912 913 ENDIF 886 ! last row 887 zfac = SIGN( 1._wp , zBU_prime(jpi-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(jpi-1,jj) ) - epsi20 ) ) 888 u_ice(jpi-1,jj) = zfac * zFU_prime(jpi-1,jj) / MAX( ABS ( zBU_prime(jpi-1,jj) ) , epsi20 ) & 889 & * umask(jpi-1,jj,1) 890 DO ji = jpi-2 , 2, -1 ! all other rows ! ---> original backward loop 891 zfac = SIGN( 1._wp , zBU_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(ji,jj) ) - epsi20 ) ) 892 u_ice(ji,jj) = zfac * ( zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji+1,jj) ) * umask(ji,jj,1) & 893 & / MAX ( ABS ( zBU_prime(ji,jj) ) , epsi20 ) 894 END DO 914 895 915 896 !--- Relaxation … … 972 953 END DO 973 954 974 !---------------975 ! Forward sweep 976 !---------------955 END DO 956 957 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFV, 'V', 1.) 977 958 978 IF ( nn_thomas == 1 ) THEN ! MV version seemingly more than mitGCM algorithm 959 !--------------- 960 ! Forward sweep 961 !--------------- 962 DO ji = ji_min, jpi - 1, nn_zebra_vp 979 963 980 DO jj = 3, jpj - 1 981 zfac = SIGN( 1._wp , zBV(ji,jj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj-1) ) - epsi20 ) ) 982 zw = zfac * zAV(ji,jj) / MAX ( ABS( zBV(ji,jj-1) ) , epsi20 ) 983 zBV_prime(ji,jj) = zBV(ji,jj) - zw * zCV(ji,jj-1) 984 zFV_prime(ji,jj) = zFV(ji,jj) - zw * zFV(ji,jj-1) 985 END DO 986 987 ELSE ! nn_thomas == 2 ! mitGCM algorithm 964 DO jj = 3, jpj - 1 965 966 zfac = SIGN( 1._wp , zBV(ji,jj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj-1) ) - epsi20 ) ) 967 zw = zfac * zAV(ji,jj) / MAX ( ABS( zBV(ji,jj-1) ) , epsi20 ) 968 zBV_prime(ji,jj) = zBV(ji,jj) - zw * zCV(ji,jj-1) 969 zFV_prime(ji,jj) = zFV(ji,jj) - zw * zFV(ji,jj-1) 970 971 END DO 972 973 END DO 974 975 CALL lbc_lnk_multi( 'icedyn_rhg_vp', zFV_prime, 'V', 1., zBV_prime, 'V', 1. ) 988 976 989 ! jj = 2 990 zw = zBV(ji,2) 991 zfac = SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) ) 992 zw = zfac / MAX ( ABS( zw ) , epsi20 ) 993 zCV_prime(ji,2) = zw * zCV(ji,2) 994 zFV_prime(ji,2) = zw * zFV(ji,2) 995 996 DO jj = 3, jpj - 1 997 zw = zBV(ji,jj) - zAV(ji,jj) * zCV_prime(ji,jj-1) 998 zfac = SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) ) 999 zw = zfac / MAX ( ABS( zw ) , epsi20 ) 1000 zCV_prime(ji,jj) = zw * zCV(ji,jj) 1001 zFV_prime(ji,jj) = zw * ( zFV(ji,jj) - zAV(ji,jj) * zFV_prime(ji,jj-1) ) 1002 END DO 1003 1004 ENDIF 1005 1006 !----------------------------- 1007 ! Backward sweep & relaxation 1008 !----------------------------- 977 !----------------------------- 978 ! Backward sweep & relaxation 979 !----------------------------- 980 DO ji = ji_min, jpi - 1, nn_zebra_vp 1009 981 1010 982 ! --- Backward sweep 1011 IF ( nn_thomas == 1 ) THEN ! nn_thomas = 1 ! --- MV version seemingly more than mitGCM algorithm 1012 1013 ! last row 1014 zfac = SIGN( 1._wp , zBV_prime(ji,jpj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jpj-1) ) - epsi20 ) ) 1015 v_ice(ji,jpj-1) = zfac * zFV_prime(ji,jpj-1) / MAX ( ABS(zBV_prime(ji,jpj-1) ) , epsi20 ) & 1016 & * vmask(ji,jpj-1,1) ! last row 1017 1018 ! other rows 1019 !DO jj = jpj-2, 2, -1 ! original back loop 1020 DO jj2 = 2 , jpj-2 ! modified forward loop 1021 jj = jpj - jj2 ! index swap 1022 zfac = SIGN( 1._wp , zBV_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jj) ) - epsi20 ) ) 1023 v_ice(ji,jj) = zfac * ( zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) ) * vmask(ji,jj,1) & 1024 & / MAX ( ABS( zBV_prime(ji,jj) ) , epsi20 ) 1025 END DO 1026 1027 ELSE ! nn_thomas == 2 ! --- mitGCM algorithm 1028 1029 ! last row 1030 v_ice(ji,jpj-1) = zFV_prime(ji,jpj-1) * vmask(ji,jpj-1,1) 1031 1032 ! other rows 1033 DO jj2 = 2 , jpj-2 ! all other rows 1034 jj = jpj - jj2 1035 v_ice(ji,jj) = ( zFV_prime(ji,jj) - zCV_prime(ji,jj) * v_ice(ji,jj+1) ) * vmask(ji,jj,1) 1036 END DO 1037 1038 ENDIF ! nn_thomas 1039 983 ! last row 984 zfac = SIGN( 1._wp , zBV_prime(ji,jpj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jpj-1) ) - epsi20 ) ) 985 v_ice(ji,jpj-1) = zfac * zFV_prime(ji,jpj-1) / MAX ( ABS(zBV_prime(ji,jpj-1) ) , epsi20 ) & 986 & * vmask(ji,jpj-1,1) ! last row 987 988 ! other rows 989 DO jj = jpj-2, 2, -1 ! original back loop 990 zfac = SIGN( 1._wp , zBV_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jj) ) - epsi20 ) ) 991 v_ice(ji,jj) = zfac * ( zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) ) * vmask(ji,jj,1) & 992 & / MAX ( ABS( zBV_prime(ji,jj) ) , epsi20 ) 993 END DO 994 1040 995 ! --- Relaxation & masking (should it be now or later) 1041 996 DO jj = 2, jpj - 1 … … 1125 1080 1126 1081 ENDIF ! ll_v_iterate 1127 1128 !--------------------------------------------------------------------------------------- 1129 ! 1130 ! --- Calculate extra convergence diagnostics and save them 1131 ! 1132 !--------------------------------------------------------------------------------------- 1133 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, & 1134 & zrhsu, zAU, zBU, zCU, zDU, zEU, zrhsv, zAV, zBV, zCV, zDV, zEV ) 1135 1136 IF ( lwp ) WRITE(numout,*) ' Done convergence tests ' 1137 1082 1138 1083 ENDIF ! --- end ll_u_iterate or ll_v_iterate 1084 1085 !--------------------------------------------------------------------------------------- 1086 ! 1087 ! --- Calculate extra convergence diagnostics and save them 1088 ! 1089 !--------------------------------------------------------------------------------------- 1090 1091 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, & 1092 & zrhsu, zAU, zBU, zCU, zDU, zEU, zrhsv, zAV, zBV, zCV, zDV, zEV ) 1093 1094 IF ( lwp ) WRITE(numout,*) ' Done convergence tests ' 1139 1095 1140 1096 END DO ! i_inn, end of inner loop … … 1596 1552 istatus = NF90_DEF_DIM( ncvgid, 'y' , jpj, iy_dim ) 1597 1553 1598 istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE , (/ idtime /), nvarid_ucvg ) ! the name uice_cvg sucks, no ?1599 1554 ! i suggest vel_dif instead 1600 istatus = NF90_DEF_VAR( ncvgid, 'u_res', NF90_DOUBLE , (/ idtime /), nvarid_ures ) 1601 istatus = NF90_DEF_VAR( ncvgid, 'v_res', NF90_DOUBLE , (/ idtime /), nvarid_vres ) 1602 istatus = NF90_DEF_VAR( ncvgid, 'vel_res', NF90_DOUBLE , (/ idtime /), nvarid_velres ) 1603 istatus = NF90_DEF_VAR( ncvgid, 'u_dif', NF90_DOUBLE , (/ idtime /), nvarid_udif ) 1604 istatus = NF90_DEF_VAR( ncvgid, 'v_dif', NF90_DOUBLE , (/ idtime /), nvarid_vdif ) 1605 istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE , (/ idtime /), nvarid_mke ) 1606 1607 istatus = NF90_DEF_VAR( ncvgid, 'usub' , NF90_DOUBLE, (/ ix_dim, iy_dim, idtime /), nvarid_usub) ! --> u-velocity in sub-iterations 1608 istatus = NF90_DEF_VAR( ncvgid, 'vsub' , NF90_DOUBLE, (/ ix_dim, iy_dim, idtime /), nvarid_vsub) ! --> v-velocity in sub-iterations 1555 istatus = NF90_DEF_VAR( ncvgid, 'u_res' , NF90_DOUBLE , (/ idtime /), nvarid_ures ) 1556 istatus = NF90_DEF_VAR( ncvgid, 'v_res' , NF90_DOUBLE , (/ idtime /), nvarid_vres ) 1557 istatus = NF90_DEF_VAR( ncvgid, 'vel_res', NF90_DOUBLE , (/ idtime /), nvarid_velres ) 1558 istatus = NF90_DEF_VAR( ncvgid, 'u_dif' , NF90_DOUBLE , (/ idtime /), nvarid_udif ) 1559 istatus = NF90_DEF_VAR( ncvgid, 'v_dif' , NF90_DOUBLE , (/ idtime /), nvarid_vdif ) 1560 istatus = NF90_DEF_VAR( ncvgid, 'vel_dif', NF90_DOUBLE , (/ idtime /), nvarid_veldif ) 1561 istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE , (/ idtime /), nvarid_mke ) 1562 1563 istatus = NF90_DEF_VAR( ncvgid, 'u_res_xy', NF90_DOUBLE, (/ ix_dim, iy_dim /), nvarid_ures_xy) 1564 istatus = NF90_DEF_VAR( ncvgid, 'v_res_xy', NF90_DOUBLE, (/ ix_dim, iy_dim /), nvarid_vres_xy) 1609 1565 1610 1566 istatus = NF90_ENDDEF(ncvgid) … … 1620 1576 ! if puerrmask and pverrmax are masked at 15% (TEST) 1621 1577 1622 IF ( lwp ) WRITE(numout,*) ' zfeldif : ', zveldif1623 1624 1578 ! --- Mean residual and kinetic energy 1625 1579 IF ( kiter == 1 ) THEN 1626 1580 1627 1628 1629 1630 1581 zu_res_mean = 0._wp 1582 zv_res_mean = 0._wp 1583 zvelres = 0._wp 1584 zmke = 0._wp 1631 1585 1632 1586 ELSE … … 1641 1595 1642 1596 z1_pglob_area = 1._wp / pglob_area 1597 1598 zu_res(:,:) = 0._wp; zv_res(:,:) = 0._wp 1643 1599 1644 1600 DO jj = 2, jpj - 1 … … 1653 1609 zv_res(ji,jj) = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) * e1e2v(ji,jj) * z1_pglob_area 1654 1610 1655 1656 1657 1658 1659 1660 1661 1662 1663 1664 1611 END DO 1612 END DO 1613 1614 IF ( lwp ) WRITE(numout,*) ' TEST 2 ' 1615 zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) ) 1616 zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) ) 1617 IF ( lwp ) WRITE(numout,*) ' TEST 3 ' 1618 zvelres = 0.5_wp * ( zu_res_mean + zv_res_mean ) 1619 1620 IF ( lwp ) WRITE(numout,*) ' TEST 4 ' 1665 1621 1666 ! -- Global mean kinetic energy per unit area (J/m2) 1667 DO jj = 2, jpj - 1 1668 DO ji = 2, jpi - 1 1669 zu = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point 1670 zv = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) ) 1671 zvel2(:,:) = zu*zu + zv*zv ! square of ice velocity at T-point 1672 END DO 1673 END DO 1622 ! -- Global mean kinetic energy per unit area (J/m2) 1623 zvel2(:,:) = 0._wp 1624 DO jj = 2, jpj - 1 1625 DO ji = 2, jpi - 1 1626 zu = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point 1627 zv = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) ) 1628 zvel2(ji,jj) = zu*zu + zv*zv ! square of ice velocity at T-point 1629 END DO 1630 END DO 1674 1631 1675 IF ( lwp ) WRITE(numout,*) ' TEST 5 ' 1676 zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) / pglob_area 1677 IF ( lwp ) WRITE(numout,*) ' TEST 6 ' 1678 1679 ENDIF 1632 IF ( lwp ) WRITE(numout,*) ' TEST 5 ' 1633 1634 zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) / pglob_area 1635 1636 IF ( lwp ) WRITE(numout,*) ' TEST 6 ' 1637 1638 ENDIF ! kiter 1680 1639 1681 1640 ! ! ==================== ! … … 1687 1646 IF( lwm ) THEN 1688 1647 ! write variables 1689 istatus = NF90_PUT_VAR( ncvgid, nvarid_ucvg, (/zveldif/), (/it/), (/1/) ) ! max U or V velocity diff between subiterations 1690 istatus = NF90_PUT_VAR( ncvgid, nvarid_ures, (/zu_res_mean/), (/it/), (/1/) ) ! U-residual of the linear system 1691 istatus = NF90_PUT_VAR( ncvgid, nvarid_vres, (/zv_res_mean/), (/it/), (/1/) ) ! V-residual of the linear system 1692 istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvelres/), (/it/), (/1/) ) ! average of u- and v- residuals 1693 istatus = NF90_PUT_VAR( ncvgid, nvarid_udif, (/puerr_max/), (/it/), (/1/) ) ! max U velocity difference, inner iterations 1694 istatus = NF90_PUT_VAR( ncvgid, nvarid_vdif, (/pverr_max/), (/it/), (/1/) ) ! max V velocity difference, inner iterations 1695 istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/it/), (/1/) ) ! mean kinetic energy 1696 1697 istatus = NF90_PUT_VAR( ncvgid, nvarid_usub, (/pu/), (/it/) ) ! u-velocity 1698 istatus = NF90_PUT_VAR( ncvgid, nvarid_vsub, (/pv/), (/it/) ) ! v-velocity 1648 istatus = NF90_PUT_VAR( ncvgid, nvarid_ures, (/zu_res_mean/), (/it/), (/1/) ) ! U-residual of the linear system 1649 istatus = NF90_PUT_VAR( ncvgid, nvarid_vres, (/zv_res_mean/), (/it/), (/1/) ) ! V-residual of the linear system 1650 istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvelres/), (/it/), (/1/) ) ! average of u- and v- residuals 1651 istatus = NF90_PUT_VAR( ncvgid, nvarid_udif, (/puerr_max/), (/it/), (/1/) ) ! max U velocity difference, inner iterations 1652 istatus = NF90_PUT_VAR( ncvgid, nvarid_vdif, (/pverr_max/), (/it/), (/1/) ) ! max V velocity difference, inner iterations 1653 istatus = NF90_PUT_VAR( ncvgid, nvarid_veldif, (/zveldif/), (/it/), (/1/) ) ! max U or V velocity diff between subiterations 1654 istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/it/), (/1/) ) ! mean kinetic energy 1655 1656 ! 1657 IF ( kiter == kitermax ) THEN 1658 WRITE(numout,*) ' Should plot the spatially dependent residual ' 1659 istatus = NF90_PUT_VAR( ncvgid, nvarid_ures_xy, (/zu_res/) ) ! U-residual, spatially dependent 1660 istatus = NF90_PUT_VAR( ncvgid, nvarid_vres_xy, (/zv_res/) ) ! V-residual, spatially dependent 1661 ENDIF 1699 1662 1700 1663 ! close file -
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/iceistate.F90
r12988 r13681 165 165 u_ice (:,:) = 0._wp 166 166 v_ice (:,:) = 0._wp 167 168 CALL lbc_lnk_multi( 'ice_istate' , u_ice, 'U', -1., v_ice, 'V', -1. ) 167 169 ! 168 170 !------------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.