Changeset 15472 for NEMO/branches/2020
- Timestamp:
- 2021-11-04T17:35:49+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
r15403 r15472 41 41 PUBLIC ice_dyn_rhg_vp ! called by icedyn_rhg.F90 42 42 43 !! 44 INTEGER :: nn_nvp ! total number of VP iterations (n_out_vp*n_inn_vp) 45 43 46 !! for convergence tests 44 47 INTEGER :: ncvgid ! netcdf file id 45 INTEGER :: nvarid_ures 46 INTEGER :: nvarid_vres 47 INTEGER :: nvarid_velres 48 INTEGER :: nvarid_udif 49 INTEGER :: nvarid_vdif 50 INTEGER :: nvarid_veldif 48 INTEGER :: nvarid_ures, nvarid_vres, nvarid_velres 49 INTEGER :: nvarid_uerr_max, nvarid_verr_max, nvarid_velerr_max 50 INTEGER :: nvarid_umad, nvarid_vmad, nvarid_velmad 51 INTEGER :: nvarid_umad_outer, nvarid_vmad_outer, nvarid_velmad_outer 51 52 INTEGER :: nvarid_mke 52 INTEGER :: nvarid_ures_xy, nvarid_vres_xy 53 54 !INTEGER :: nvarid_ures_xy, nvarid_vres_xy 55 56 !!! -> store velocity at last inner iteration zu_end_outer, zv_end_outer 57 !!! -> send it to rhg_cvg 58 !!! -> for 2nd and larger iteration, calculate mad_outer and store it (rhg_cvg) 59 !!! -> store spatially dependent residual (rhg_cvg) 60 !!! -> store spatially dependent mad_outer (rhg_cvg) 53 61 54 62 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmsk00, zmsk15 … … 132 140 ! 133 141 INTEGER :: ji, ji2, jj, jj2, jn ! dummy loop indices 134 INTEGER :: jter, i_out, i_inn!142 INTEGER :: i_out, i_inn, i_inn_tot ! 135 143 INTEGER :: ji_min, jj_min ! 136 144 INTEGER :: nn_zebra_vp ! number of zebra steps 137 145 138 INTEGER :: nn_nvp ! total number of VP iterations (n_out_vp*n_inn_vp)139 146 ! 140 147 REAL(wp) :: zrhoco ! rau0 * rn_cio … … 167 174 REAL(wp), DIMENSION(jpi,jpj) :: zu_c, zv_c ! "current" ice velocity (m/s), average of previous two OL iterates 168 175 REAL(wp), DIMENSION(jpi,jpj) :: zu_b, zv_b ! velocity at previous sub-iterate 176 REAL(wp), DIMENSION(jpi,jpj) :: zu_b_outer, zv_b_outer ! velocity at previous outer-iterate 169 177 REAL(wp), DIMENSION(jpi,jpj) :: zuerr, zverr ! absolute U/Vvelocity difference between current and previous sub-iterates 170 178 ! … … 197 205 REAL(wp), PARAMETER :: zamin = 0.001_wp ! ice concentration below which ice velocity becomes very small 198 206 !! --- diags 199 REAL(wp) :: zsig1, zsig2, zsig12, zdelta, z1_strength, zfac_x, zfac_y207 REAL(wp) :: zsig1, zsig2, zsig12, zdelta, z1_strength, zfac_x, zfac_y 200 208 REAL(wp), DIMENSION(jpi,jpj) :: zs1, zs2, zs12, zs12f ! stress tensor components 201 209 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zsig_I, zsig_II, zsig1_p, zsig2_p … … 204 212 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xmtrp_snw, zdiag_ymtrp_snw ! X/Y-component of snow mass transport (kg/s, SIMIP) 205 213 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdiag_xatrp, zdiag_yatrp ! X/Y-component of area transport (m2/s, SIMIP) 214 215 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zvel_res ! Residual of the linear system at last iteration 216 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zvel_diff ! Absolute velocity difference @last outer iteration 206 217 207 218 !!---------------------------------------------------------------------------------------------------------------------- … … 261 272 ! Initialise convergence checks 262 273 IF( ln_rhg_chkcvg ) THEN 263 264 ! ice area for global mean kinetic energy 265 zglob_area = glob_sum( 'ice_rhg_vp', at_i(:,:) * e1e2t(:,:) ) ! global ice area (km2) 274 275 ! allocate spatially-resolved convergence diagnostics 276 ALLOCATE( zvel_res(jpi,jpj), zvel_diff(jpi,jpj) ) 277 278 ! ice area for proper weighted mean averages 279 zglob_area = glob_sum( 'ice_rhg_vp', at_i(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! global ice area (km2) 266 280 267 281 ENDIF … … 417 431 zv_c(:,:) = v_ice(:,:) 418 432 419 jter= 0433 i_inn_tot = 0 420 434 421 435 DO i_out = 1, nn_nout_vp 422 436 423 437 IF( lwp ) WRITE(numout,*) ' outer loop i_out : ', i_out 424 438 439 zu_b_outer(:,:) = u_ice(:,:) ! velocity at previous outer-iterate 440 zv_b_outer(:,:) = v_ice(:,:) 441 425 442 ! Velocities used in the non linear terms are the average of the past two iterates 426 443 ! u_it = 0.5 * ( u_{it-1} + u_{it-2}) 427 444 ! Also used in Hibler and Ackley (1983); Zhang and Hibler (1997); Lemieux and Tremblay (2009) 428 zu_c(:,:) = 0.5_wp * ( u_ice(:,:) + zu_c(:,:) ) 429 zv_c(:,:) = 0.5_wp * ( v_ice(:,:) + zv_c(:,:) ) 445 446 zu_c(:,:) = 0.5_wp * ( u_ice(:,:) + zu_c(:,:) ) 447 zv_c(:,:) = 0.5_wp * ( v_ice(:,:) + zv_c(:,:) ) 430 448 431 449 !------------------------------------------------------------------------------! … … 525 543 END DO 526 544 527 ! zef(:,:) = 0. ! test528 529 545 CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. ) 530 546 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zef' , zef ) ! MV DEBUG … … 553 569 554 570 ! --- U-component of Coriolis Force (energy conserving formulation) 555 ! Note Lemieux et al 2008 recommend to do that implicitly, but I don't really see how this could be done556 571 zCorU(ji,jj) = 0.25_wp * r1_e1u(ji,jj) * & 557 572 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * zv_c(ji ,jj) + e1v(ji ,jj-1) * zv_c(ji ,jj-1) ) & … … 606 621 END DO 607 622 608 ! TEST609 ! zs2_rhsv(:,:) = 0._wp610 611 623 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs1_rhsu' , zs1_rhsu ) ! MV DEBUG 612 624 IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs2_rhsu' , zs2_rhsu ) ! MV DEBUG … … 633 645 END DO 634 646 END DO 635 !636 ! TEST637 ! zs12_rhsv(:,:) = 0._wp638 647 639 648 CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsu, 'F', 1. ) … … 822 831 !--- mitgcm computes initial value of residual here... 823 832 824 jter = jter + 1 825 ! l_full_nf_update = jter == nn_nvp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 826 827 zu_b(:,:) = u_ice(:,:) ! velocity at previous sub-iterate 828 zv_b(:,:) = v_ice(:,:) 829 830 !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 833 i_inn_tot = i_inn_tot + 1 834 ! l_full_nf_update = i_inn_tot == nn_nvp ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 835 836 zu_b(:,:) = u_ice(:,:) ! velocity at previous inner-iterate 837 zv_b(:,:) = v_ice(:,:) 832 838 833 839 IF ( ll_u_iterate .OR. ll_v_iterate ) THEN … … 1049 1055 ! MV OPT: if the number of iterations to convergence is really variable, and keep the convergence check 1050 1056 ! then we must optimize the use of the mpp_max, which is prohibitive 1051 zuerr_max = 0._wp1057 zuerr_max = 0._wp 1052 1058 1053 1059 IF ( ll_u_iterate .AND. MOD ( i_inn, nn_cvgchk_vp ) == 0 ) THEN … … 1057 1063 DO jj = 2, jpj - 1 1058 1064 DO ji = 2, jpi - 1 1059 zuerr(ji,jj) = ABS ( ( u_ice(ji,jj) - zu_b(ji,jj) ) ) * umask(ji,jj,1) 1065 zuerr(ji,jj) = ABS ( ( u_ice(ji,jj) - zu_b(ji,jj) ) ) * umask(ji,jj,1) 1060 1066 END DO 1061 1067 END DO 1062 1068 zuerr_max = MAXVAL( zuerr ) 1063 1069 CALL mpp_max( 'icedyn_rhg_evp', zuerr_max ) ! max over the global domain - damned! 1064 1065 ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine)1070 1071 ! - Stop if max :error is too large ("safeguard against bad forcing" of original Zhang routine) 1066 1072 IF ( i_inn > 1 .AND. zuerr_max > rn_uerr_max_vp ) THEN 1067 1073 IF ( lwp ) WRITE(numout,*) " VP rheology error was too large : ", zuerr_max, " in outer U-iteration ", i_out, " after ", i_inn, " iterations, we stopped " … … 1091 1097 END DO 1092 1098 END DO 1093 1094 1099 zverr_max = MAXVAL( zverr ) 1095 1100 CALL mpp_max( 'icedyn_rhg_evp', zverr_max ) ! max over the global domain - damned! 1096 1101 1097 1102 ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine) 1098 1103 IF ( i_inn > 1 .AND. zverr_max > rn_uerr_max_vp ) THEN … … 1116 1121 ! 1117 1122 !--------------------------------------------------------------------------------------- 1118 1119 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, & 1120 & zrhsu, zAU, zBU, zCU, zDU, zEU, zrhsv, zAV, zBV, zCV, zDV, zEV ) 1123 IF( ln_rhg_chkcvg .AND. MOD ( i_inn - 1, nn_cvgchk_vp ) == 0 ) THEN 1124 1125 CALL rhg_cvg_vp( kt, i_out, i_inn, i_inn_tot, nn_nout_vp, nn_ninn_vp, nn_nvp, & 1126 & u_ice, v_ice, zu_b, zv_b, zu_b_outer, zv_b_outer, & 1127 & zmt, za_iU, za_iV, zuerr_max, zverr_max, zglob_area, & 1128 & zrhsu, zAU, zBU, zCU, zDU, zEU, zFU, & 1129 & zrhsv, zAV, zBV, zCV, zDV, zEV, zFV, & 1130 zvel_res, zvel_diff ) 1131 1132 ENDIF 1121 1133 1122 1134 IF ( lwp ) WRITE(numout,*) ' Done convergence tests ' … … 1145 1157 CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 1146 1158 1147 IF ( lwp ) WRITE(numout,*) ' We are about to output uice_dbg '1148 IF( iom_use('uice_dbg' ) ) CALL iom_put( 'uice_dbg' , u_ice ) ! ice velocity u after solver1149 IF( iom_use('vice_dbg' ) ) CALL iom_put( 'vice_dbg' , v_ice ) ! ice velocity v after solver1150 1151 !------------------------------------------------------------------------------!1152 !1153 ! --- Convergence diagnostics1154 !1155 !------------------------------------------------------------------------------!1156 1159 1157 1160 IF( ln_rhg_chkcvg ) THEN 1158 1161 1159 IF( iom_use('uice_cvg') ) THEN 1160 CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_b(:,:) ) * umask(:,:,1) , & ! ice velocity difference at last iteration 1161 & ABS( v_ice(:,:) - zv_b(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 1162 ENDIF 1162 IF( iom_use('velo_res') ) CALL iom_put( 'velo_res', zvel_res ) ! linear system residual @last inner&outer iteration 1163 IF( iom_use('velo_ero') ) CALL iom_put( 'velo_ero', zvel_diff ) ! abs velocity difference @last outer iteration 1164 IF( iom_use('uice_eri') ) CALL iom_put( 'uice_eri', zuerr ) ! abs velocity difference @last inner iteration 1165 IF( iom_use('vice_eri') ) CALL iom_put( 'vice_eri', zverr ) ! abs velocity difference @last inner iteration 1166 1167 DEALLOCATE( zvel_res , zvel_diff ) 1163 1168 1164 1169 ENDIF ! ln_rhg_chkcvg … … 1458 1463 DO jj = 2, jpj - 1 1459 1464 DO ji = 2, jpi - 1 1465 1460 1466 zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 1461 1467 & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) & … … 1464 1470 & ) * 2._wp * r1_e1u(ji,jj) & 1465 1471 & ) * r1_e1e2u(ji,jj) 1472 1466 1473 zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 1467 1474 & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) & … … 1470 1477 & ) * 2._wp * r1_e2v(ji,jj) & 1471 1478 & ) * r1_e1e2v(ji,jj) 1479 1472 1480 END DO 1473 1481 END DO … … 1526 1534 1527 1535 1528 1529 SUBROUTINE rhg_cvg_vp( kt, kiter, kitermax, pu, pv, pmt, puerr_max, pverr_max, pglob_area, & 1530 & prhsu, pAU, pBU, pCU, pDU, pEU, prhsv, pAV, pBV, pCV, pDV, pEV ) 1531 1536 SUBROUTINE rhg_cvg_vp( kt, kitout, kitinn, kitinntot, kitoutmax, kitinnmax, kitinntotmax , & 1537 & pu, pv, pub, pvb, pub_outer, pvb_outer , & 1538 & pmt, pat_iu, pat_iv, puerr_max, pverr_max, pglob_area , & 1539 & prhsu, pAU, pBU, pCU, pDU, pEU, pFU , & 1540 & prhsv, pAV, pBV, pCV, pDV, pEV, pFV , & 1541 & pvel_res, pvel_diff ) 1542 !! 1532 1543 !!---------------------------------------------------------------------- 1533 1544 !! *** ROUTINE rhg_cvg_vp *** … … 1544 1555 !! 1545 1556 !! ** Note : for the first sub-iteration, uice_cvg is set to 0 (too large otherwise) 1557 !! 1546 1558 !!---------------------------------------------------------------------- 1547 INTEGER , INTENT(in) :: kt, kiter, kitermax ! ocean time-step index 1548 REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pmt ! now velocity and mass per unit area 1549 REAL(wp), INTENT(in) :: puerr_max, pverr_max ! absolute mean velocity difference 1550 REAL(wp), INTENT(in) :: pglob_area ! global ice area 1551 REAL(wp), DIMENSION(:,:), INTENT(in) :: prhsu, pAU, pBU, pCU, pDU, pEU ! linear system coefficients 1552 REAL(wp), DIMENSION(:,:), INTENT(in) :: prhsv, pAV, pBV, pCV, pDV, pEV 1553 !! 1554 INTEGER :: it, idtime, istatus, ix_dim, iy_dim 1559 !! 1560 INTEGER , INTENT(in) :: kt, kitout, kitinn, kitinntot ! ocean model iterate, outer, inner and total n-iterations 1561 INTEGER , INTENT(in) :: kitoutmax, kitinnmax ! max number of outer & inner iterations 1562 INTEGER , INTENT(in) :: kitinntotmax ! max number of total sub-iterations 1563 REAL(wp), DIMENSION(:,:), INTENT(in) :: pu, pv, pub, pvb ! now & sub-iter-before velocities 1564 REAL(wp), DIMENSION(:,:), INTENT(in) :: pub_outer, pvb_outer ! velocities @before outer iterations 1565 REAL(wp), DIMENSION(:,:), INTENT(in) :: pmt, pat_iu, pat_iv ! mass at T-point, ice concentration at U&V 1566 REAL(wp), INTENT(in) :: puerr_max, pverr_max ! absolute mean velocity difference 1567 REAL(wp), INTENT(in) :: pglob_area ! global ice area 1568 REAL(wp), DIMENSION(:,:), INTENT(in) :: prhsu, pAU, pBU, pCU, pDU, pEU, pFU ! linear system coefficients 1569 REAL(wp), DIMENSION(:,:), INTENT(in) :: prhsv, pAV, pBV, pCV, pDV, pEV, pFV 1570 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pvel_res ! velocity residual @last inner iteration 1571 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pvel_diff ! velocity difference @last outer iteration 1572 !! 1573 1574 INTEGER :: idtime, istatus, ix_dim, iy_dim 1555 1575 INTEGER :: ji, jj ! dummy loop indices 1556 REAL(wp) :: zveldif, zu_res_mean, zv_res_mean, zvelres, zmke, zu, zv ! local scalars 1557 REAL(wp) :: z1_pglob_area 1576 INTEGER :: it_inn_file, it_out_file 1577 REAL(wp) :: zu_res_mean, zv_res_mean, zvel_res_mean ! mean residuals of the linear system 1578 REAL(wp) :: zu_mad, zv_mad, zvel_mad ! mean absolute deviation, sub-iterates 1579 REAL(wp) :: zu_mad_outer, zv_mad_outer, zvel_mad_outer ! mean absolute deviation, outer-iterates 1580 REAL(wp) :: zvel_err_max, zmke, zu, zv ! local scalars 1581 REAL(wp) :: z1_pglob_area ! inverse global ice area 1582 1558 1583 REAL(wp), DIMENSION(jpi,jpj) :: zu_res, zv_res, zvel2 ! local arrays 1584 REAL(wp), DIMENSION(jpi,jpj) :: zu_diff, zv_diff ! local arrays 1559 1585 1560 1586 CHARACTER(len=20) :: clname 1561 1587 !!---------------------------------------------------------------------- 1562 1588 1589 1563 1590 IF( lwp ) THEN 1591 1564 1592 WRITE(numout,*) 1565 1593 WRITE(numout,*) 'rhg_cvg_vp : ice rheology convergence control' 1566 1594 WRITE(numout,*) '~~~~~~~~~~~' 1567 WRITE(numout,*) ' kiter = : ', kiter 1568 WRITE(numout,*) ' kitermax = : ', kitermax 1595 WRITE(numout,*) ' kt = : ', kt 1596 WRITE(numout,*) ' kitout = : ', kitout 1597 WRITE(numout,*) ' kitinn = : ', kitinn 1598 WRITE(numout,*) ' kitinntot = : ', kitinntot 1599 WRITE(numout,*) ' kitoutmax (nn_nout_vp) = ', kitoutmax 1600 WRITE(numout,*) ' kitinnmax (nn_ninn_vp) = ', kitinnmax 1601 WRITE(numout,*) ' kitinntotmax (nn_nvp) = ', kitinntotmax 1602 WRITE(numout,*) 1603 1569 1604 ENDIF 1570 1605 1606 z1_pglob_area = 1._wp / pglob_area ! inverse global ice area 1607 1571 1608 ! create file 1572 IF( kt == nit000 .AND. kit er== 1 ) THEN1609 IF( kt == nit000 .AND. kitinntot == 1 ) THEN 1573 1610 ! 1574 1611 IF( lwm ) THEN … … 1582 1619 istatus = NF90_DEF_DIM( ncvgid, 'y' , jpj, iy_dim ) 1583 1620 1584 ! i suggest vel_dif instead 1585 istatus = NF90_DEF_VAR( ncvgid, 'u_res' , NF90_DOUBLE , (/ idtime /), nvarid_ures ) 1586 istatus = NF90_DEF_VAR( ncvgid, 'v_res' , NF90_DOUBLE , (/ idtime /), nvarid_vres ) 1587 istatus = NF90_DEF_VAR( ncvgid, 'vel_res', NF90_DOUBLE , (/ idtime /), nvarid_velres ) 1588 istatus = NF90_DEF_VAR( ncvgid, 'u_dif' , NF90_DOUBLE , (/ idtime /), nvarid_udif ) 1589 istatus = NF90_DEF_VAR( ncvgid, 'v_dif' , NF90_DOUBLE , (/ idtime /), nvarid_vdif ) 1590 istatus = NF90_DEF_VAR( ncvgid, 'vel_dif', NF90_DOUBLE , (/ idtime /), nvarid_veldif ) 1621 istatus = NF90_DEF_VAR( ncvgid, 'u_res' , NF90_DOUBLE , (/ idtime /), nvarid_ures ) 1622 istatus = NF90_DEF_VAR( ncvgid, 'v_res' , NF90_DOUBLE , (/ idtime /), nvarid_vres ) 1623 istatus = NF90_DEF_VAR( ncvgid, 'vel_res' , NF90_DOUBLE , (/ idtime /), nvarid_velres ) 1624 1625 istatus = NF90_DEF_VAR( ncvgid, 'uerr_max_sub' , NF90_DOUBLE , (/ idtime /), nvarid_uerr_max ) 1626 istatus = NF90_DEF_VAR( ncvgid, 'verr_max_sub' , NF90_DOUBLE , (/ idtime /), nvarid_verr_max ) 1627 istatus = NF90_DEF_VAR( ncvgid, 'velerr_max_sub', NF90_DOUBLE , (/ idtime /), nvarid_velerr_max ) 1628 1629 istatus = NF90_DEF_VAR( ncvgid, 'umad_sub' , NF90_DOUBLE , (/ idtime /), nvarid_umad ) 1630 istatus = NF90_DEF_VAR( ncvgid, 'vmad_sub' , NF90_DOUBLE , (/ idtime /), nvarid_vmad ) 1631 istatus = NF90_DEF_VAR( ncvgid, 'velmad_sub' , NF90_DOUBLE , (/ idtime /), nvarid_velmad ) 1632 1633 istatus = NF90_DEF_VAR( ncvgid, 'umad_outer' , NF90_DOUBLE , (/ idtime /), nvarid_umad_outer ) 1634 istatus = NF90_DEF_VAR( ncvgid, 'vmad_outer' , NF90_DOUBLE , (/ idtime /), nvarid_vmad_outer ) 1635 istatus = NF90_DEF_VAR( ncvgid, 'velmad_outer' , NF90_DOUBLE , (/ idtime /), nvarid_velmad_outer ) 1636 1591 1637 istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE , (/ idtime /), nvarid_mke ) 1592 1593 istatus = NF90_DEF_VAR( ncvgid, 'u_res_xy', NF90_DOUBLE, (/ ix_dim, iy_dim /), nvarid_ures_xy)1594 istatus = NF90_DEF_VAR( ncvgid, 'v_res_xy', NF90_DOUBLE, (/ ix_dim, iy_dim /), nvarid_vres_xy)1595 1638 1596 1639 istatus = NF90_ENDDEF(ncvgid) … … 1600 1643 ENDIF 1601 1644 1602 IF ( lwp ) WRITE(numout,*) ' File created ' 1603 1604 ! --- Max absolute velocity difference with previous iterate (zveldif) 1605 zveldif = MAX( puerr_max, pverr_max ) ! velocity difference with previous iterate, should nearly be equivalent to evp code 1606 ! if puerrmask and pverrmax are masked at 15% (TEST) 1607 1608 ! --- Mean residual and kinetic energy 1609 IF ( kiter == 1 ) THEN 1610 1611 zu_res_mean = 0._wp 1612 zv_res_mean = 0._wp 1613 zvelres = 0._wp 1614 zmke = 0._wp 1645 !------------------------------------------------------------ 1646 ! 1647 ! Max absolute velocity difference with previous sub-iterate 1648 ! ( zvel_err_max ) 1649 ! 1650 !------------------------------------------------------------ 1651 ! 1652 ! This comes from the criterion used to stop the iterative procedure 1653 zvel_err_max = 0.5_wp * ( puerr_max + pverr_max ) ! average of U- and V- maximum error over the whole domain 1654 1655 !---------------------------------------------- 1656 ! 1657 ! Mean-absolute-deviation (sub-iterates) 1658 ! ( zu_mad, zv_mad, zvel_mad) 1659 ! 1660 !---------------------------------------------- 1661 ! 1662 ! U 1663 zu_diff(:,:) = 0._wp 1664 DO jj = 2, jpj - 1 1665 DO ji = 2, jpi - 1 1666 zu_diff(ji,jj) = ABS ( ( pu(ji,jj) - pub(ji,jj) ) ) * e1e2u(ji,jj) * pat_iu(ji,jj) * umask(ji,jj,1) * z1_pglob_area 1667 END DO 1668 END DO 1669 1670 ! V 1671 zv_diff(:,:) = 0._wp 1672 DO jj = 2, jpj - 1 1673 DO ji = 2, jpi - 1 1674 zv_diff(ji,jj) = ABS ( ( pv(ji,jj) - pvb(ji,jj) ) ) * e1e2v(ji,jj) * pat_iv(ji,jj) * vmask(ji,jj,1) * z1_pglob_area 1675 END DO 1676 END DO 1677 1678 1679 ! global sum & U-V average 1680 CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_diff, 'U', 1., zv_diff , 'V', 1. ) 1681 zu_mad = glob_sum( 'icedyn_rhg_vp : ', zu_diff ) 1682 zv_mad = glob_sum( 'icedyn_rhg_vp : ', zv_diff ) 1683 1684 zvel_mad = 0.5_wp * ( zu_mad + zv_mad ) 1685 1686 !----------------------------------------------- 1687 ! 1688 ! Mean-absolute-deviation (outer-iterates) 1689 ! ( zu_mad_outer, zv_mad_outer, zvel_mad_outer) 1690 ! 1691 !----------------------------------------------- 1692 ! 1693 IF ( kitinn == kitinnmax ) THEN ! only work at the end of outer iterates 1694 1695 ! * U 1696 zu_diff(:,:) = 0._wp 1697 DO jj = 2, jpj - 1 1698 DO ji = 2, jpi - 1 1699 zu_diff(ji,jj) = ABS ( ( pu(ji,jj) - pub_outer(ji,jj) ) ) * e1e2u(ji,jj) * pat_iu(ji,jj) * umask(ji,jj,1) * & 1700 & z1_pglob_area 1701 END DO 1702 END DO 1703 1704 ! * V 1705 zv_diff(:,:) = 0._wp 1706 DO jj = 2, jpj - 1 1707 DO ji = 2, jpi - 1 1708 zv_diff(ji,jj) = ABS ( ( pv(ji,jj) - pvb_outer(ji,jj) ) ) * e1e2v(ji,jj) * pat_iv(ji,jj) * vmask(ji,jj,1) * & 1709 & z1_pglob_area 1710 END DO 1711 END DO 1712 1713 ! Global ice-concentration, grid-cell-area weighted mean 1714 CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_diff, 'U', 1., zv_diff , 'V', 1. ) ! abs behaves as a scalar no ? 1715 1716 zu_mad_outer = glob_sum( 'icedyn_rhg_vp : ', zu_diff ) 1717 zv_mad_outer = glob_sum( 'icedyn_rhg_vp : ', zv_diff ) 1718 1719 ! Average of both U & V 1720 zvel_mad_outer = 0.5_wp * ( zu_mad_outer + zv_mad_outer ) 1721 1722 ENDIF 1723 1724 ! --- Spatially-resolved absolute difference to send back to main routine 1725 ! (last iteration only, T-point) 1726 1727 IF ( kitinntot == kitinntotmax) THEN 1728 1729 zu_diff(:,:) = 0._wp 1730 zv_diff(:,:) = 0._wp 1731 1732 DO jj = 2, jpj - 1 1733 DO ji = 2, jpi - 1 1734 1735 zu_diff(ji,jj) = ( ABS ( ( pu(ji-1,jj) - pub_outer(ji-1,jj) ) ) * umask(ji-1,jj,1) & 1736 & + ABS ( ( pu(ji,jj ) - pub_outer(ji,jj) ) ) * umask(ji,jj,1) ) & 1737 & / ( umask(ji-1,jj,1) + umask(ji,jj,1) ) 1738 1739 zv_diff(ji,jj) = ( ABS ( ( pv(ji,jj-1) - pvb_outer(ji,jj-1) ) ) * vmask(ji,jj-1,1) & 1740 & + ABS ( ( pv(ji,jj ) - pvb_outer(ji,jj) ) ) * vmask(ji,jj,1) & 1741 & / ( vmask(ji,jj-1,1) + vmask(ji,jj,1) ) ) 1742 1743 END DO 1744 END DO 1745 1746 CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_diff, 'T', 1., zv_diff , 'T', 1. ) 1747 pvel_diff(:,:) = 0.5_wp * ( zu_diff(:,:) + zv_diff(:,:) ) 1615 1748 1616 1749 ELSE 1617 1750 1618 ! -- Mean residual (N/m^2), zu_res_mean 1619 ! Here we take the residual of the linear system (N/m^2), 1620 ! We define it as in mitgcm: square-root of area-weighted mean square residual 1621 ! Local residual r = Ax - B expresses to which extent the momentum balance is verified 1622 ! i.e., how close we are to a solution 1623 1624 IF ( lwp ) WRITE(numout,*) ' TEST 1 ' 1625 1626 z1_pglob_area = 1._wp / pglob_area 1627 1628 zu_res(:,:) = 0._wp; zv_res(:,:) = 0._wp 1629 1630 DO jj = 2, jpj - 1 1631 DO ji = 2, jpi - 1 1632 zu_res(ji,jj) = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1) & 1633 & - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) ) 1634 1635 zv_res(ji,jj) = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj) & 1636 & - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) ) 1637 1638 zu_res(ji,jj) = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) * e1e2u(ji,jj) * z1_pglob_area 1639 zv_res(ji,jj) = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) * e1e2v(ji,jj) * z1_pglob_area 1640 1641 END DO 1642 END DO 1643 1644 IF ( lwp ) WRITE(numout,*) ' TEST 2 ' 1645 zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) ) 1646 zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) ) 1647 IF ( lwp ) WRITE(numout,*) ' TEST 3 ' 1648 zvelres = 0.5_wp * ( zu_res_mean + zv_res_mean ) 1649 1650 IF ( lwp ) WRITE(numout,*) ' TEST 4 ' 1651 1652 ! -- Global mean kinetic energy per unit area (J/m2) 1653 zvel2(:,:) = 0._wp 1654 DO jj = 2, jpj - 1 1655 DO ji = 2, jpi - 1 1656 zu = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point 1657 zv = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) ) 1658 zvel2(ji,jj) = zu*zu + zv*zv ! square of ice velocity at T-point 1659 END DO 1660 END DO 1661 1662 IF ( lwp ) WRITE(numout,*) ' TEST 5 ' 1663 1664 zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) / pglob_area 1665 1666 IF ( lwp ) WRITE(numout,*) ' TEST 6 ' 1667 1668 ENDIF ! kiter 1751 pvel_diff(:,:) = 0._wp 1752 1753 ENDIF 1754 1755 !--------------------------------------- 1756 ! 1757 ! --- Mean residual & kinetic energy 1758 ! 1759 !--------------------------------------- 1760 1761 IF ( kitinntot == 1 ) THEN 1762 1763 zu_res_mean = 0._wp 1764 zv_res_mean = 0._wp 1765 zvel_res_mean = 0._wp 1766 zmke = 0._wp 1767 1768 ELSE 1769 1770 ! * Mean residual (N/m2) 1771 ! Here we take the residual of the linear system (N/m2), 1772 ! We define it as in mitgcm: global area-weighted mean of square-root residual 1773 ! Local residual r = Ax - B expresses to which extent the momentum balance is verified 1774 ! i.e., how close we are to a solution 1775 !!!! UNITS OF RESIDUAL ARE NOT m/s BUT N/m2 (CHECK) 1776 zu_res(:,:) = 0._wp; zv_res(:,:) = 0._wp 1777 1778 DO jj = 2, jpj - 1 1779 DO ji = 2, jpi - 1 1780 zu_res(ji,jj) = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1) & 1781 & - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) ) 1782 zv_res(ji,jj) = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj) & 1783 & - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) ) 1784 1785 ! zu_res(ji,jj) = pFU(ji,jj) - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) 1786 ! zv_res(ji,jj) = pFV(ji,jj) - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) 1787 1788 zu_res(ji,jj) = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) * pat_iu(ji,jj) * e1e2u(ji,jj) * z1_pglob_area 1789 zv_res(ji,jj) = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) * pat_iv(ji,jj) * e1e2v(ji,jj) * z1_pglob_area 1790 1791 END DO 1792 END DO 1793 1794 ! Global ice-concentration, grid-cell-area weighted mean 1795 CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_res, 'U', 1., zv_res , 'V', 1. ) 1796 1797 zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) ) 1798 zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) ) 1799 zvel_res_mean = 0.5_wp * ( zu_res_mean + zv_res_mean ) 1800 1801 ! --- Global mean kinetic energy per unit area (J/m2) 1802 zvel2(:,:) = 0._wp 1803 DO jj = 2, jpj - 1 1804 DO ji = 2, jpi - 1 1805 zu = 0.5_wp * ( pu(ji-1,jj) + pu(ji,jj) ) ! u-vel at T-point 1806 zv = 0.5_wp * ( pv(ji,jj-1) + pv(ji,jj) ) 1807 zvel2(ji,jj) = zu*zu + zv*zv ! square of ice velocity at T-point 1808 END DO 1809 END DO 1810 1811 zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) / pglob_area 1812 1813 ENDIF ! kitinntot 1814 1815 !--- Spatially-resolved residual at last iteration to send back to main routine (last iteration only) 1816 !--- Calculation @T-point 1817 1818 IF ( kitinntot == kitinntotmax) THEN 1819 1820 DO jj = 2, jpj - 1 ! @ residuals at u and v points 1821 DO ji = 2, jpi - 1 1822 1823 zu_res(ji,jj) = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1) & 1824 & - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) ) 1825 zv_res(ji,jj) = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj) & 1826 & - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) ) 1827 1828 zu_res(ji,jj) = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) 1829 zv_res(ji,jj) = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) 1830 1831 END DO 1832 END DO 1833 1834 CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', zu_res, 'U', 1., zv_res , 'V', 1. ) 1835 1836 DO jj = 2, jpj - 1 ! average @T-point 1837 DO ji = 2, jpi - 1 1838 pvel_res(ji,jj) = 0.25_wp * ( zu_res(ji-1,jj) + zu_res(ji,jj) + zv_res(ji,jj-1) + zv_res(ji,jj) ) 1839 END DO 1840 END DO 1841 1842 CALL lbc_lnk_multi( 'icedyn_rhg_cvg_vp', pvel_res, 'T', 1. ) 1843 1844 ELSE 1845 1846 pvel_res(:,:) = 0._wp 1847 1848 ENDIF 1669 1849 1670 1671 1672 ! time1673 it = ( kt - 1 ) * kitermax + kiter1674 1675 1850 ! ! ==================== ! 1851 1852 it_inn_file = ( kt - 1 ) * kitinntotmax + kitinntot ! time step in the file 1853 it_out_file = ( kt - 1 ) * kitoutmax + kitout 1854 1855 ! write variables 1676 1856 IF( lwm ) THEN 1677 ! write variables 1678 istatus = NF90_PUT_VAR( ncvgid, nvarid_ures, (/zu_res_mean/), (/it/), (/1/) ) ! U-residual of the linear system 1679 istatus = NF90_PUT_VAR( ncvgid, nvarid_vres, (/zv_res_mean/), (/it/), (/1/) ) ! V-residual of the linear system 1680 istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvelres/), (/it/), (/1/) ) ! average of u- and v- residuals 1681 istatus = NF90_PUT_VAR( ncvgid, nvarid_udif, (/puerr_max/), (/it/), (/1/) ) ! max U velocity difference, inner iterations 1682 istatus = NF90_PUT_VAR( ncvgid, nvarid_vdif, (/pverr_max/), (/it/), (/1/) ) ! max V velocity difference, inner iterations 1683 istatus = NF90_PUT_VAR( ncvgid, nvarid_veldif, (/zveldif/), (/it/), (/1/) ) ! max U or V velocity diff between subiterations 1684 istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/it/), (/1/) ) ! mean kinetic energy 1685 1686 ! 1687 IF ( kiter == kitermax ) THEN 1688 WRITE(numout,*) ' Should plot the spatially dependent residual ' 1689 istatus = NF90_PUT_VAR( ncvgid, nvarid_ures_xy, (/zu_res/) ) ! U-residual, spatially dependent 1690 istatus = NF90_PUT_VAR( ncvgid, nvarid_vres_xy, (/zv_res/) ) ! V-residual, spatially dependent 1857 1858 istatus = NF90_PUT_VAR( ncvgid, nvarid_ures , (/zu_res_mean/), (/it_inn_file/), (/1/) ) ! Residuals of the linear system, area weighted mean 1859 istatus = NF90_PUT_VAR( ncvgid, nvarid_vres , (/zv_res_mean/), (/it_inn_file/), (/1/) ) ! 1860 istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvel_res_mean/), (/it_inn_file/), (/1/) ) ! 1861 1862 istatus = NF90_PUT_VAR( ncvgid, nvarid_uerr_max , (/puerr_max/), (/it_inn_file/), (/1/) ) ! Max velocit_inn_filey error, sub-it_inn_fileerates 1863 istatus = NF90_PUT_VAR( ncvgid, nvarid_verr_max , (/pverr_max/), (/it_inn_file/), (/1/) ) ! 1864 istatus = NF90_PUT_VAR( ncvgid, nvarid_velerr_max, (/zvel_err_max/), (/it_inn_file/), (/1/) ) ! 1865 1866 istatus = NF90_PUT_VAR( ncvgid, nvarid_umad , (/zu_mad/) , (/it_inn_file/), (/1/) ) ! velocit_inn_filey MAD, area/sic-weighted, sub-it_inn_fileerates 1867 istatus = NF90_PUT_VAR( ncvgid, nvarid_vmad , (/zv_mad/) , (/it_inn_file/), (/1/) ) ! 1868 istatus = NF90_PUT_VAR( ncvgid, nvarid_velmad , (/zvel_mad/), (/it_inn_file/), (/1/) ) ! 1869 1870 istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/kitinntot/), (/1/) ) ! mean kinetic energy 1871 1872 IF ( kitinn == kitinnmax ) THEN ! only print outer mad at the end of inner loop 1873 1874 istatus = NF90_PUT_VAR( ncvgid, nvarid_umad_outer , (/zu_mad_outer/) , (/it_out_file/), (/1/) ) ! velocity MAD, area/sic-weighted, outer-iterates 1875 istatus = NF90_PUT_VAR( ncvgid, nvarid_vmad_outer , (/zv_mad_outer/) , (/it_out_file/), (/1/) ) ! 1876 istatus = NF90_PUT_VAR( ncvgid, nvarid_velmad_outer , (/zvel_mad_outer/), (/it_out_file/), (/1/) ) ! 1877 1691 1878 ENDIF 1879 1880 ! spatially-dependent things (to be coded) 1881 ! velocity difference ? 1882 ! residual ? 1883 ! ! 1884 ! IF ( kiter == kitermax ) THEN 1885 ! WRITE(numout,*) ' Should plot the spatially dependent residual ' 1886 ! istatus = NF90_PUT_VAR( ncvgid, nvarid_ures_xy, (/zu_res/) ) ! U-residual, spatially dependent 1887 ! istatus = NF90_PUT_VAR( ncvgid, nvarid_vres_xy, (/zv_res/) ) ! V-residual, spatially dependent 1888 ! ENDIF 1692 1889 1693 1890 ! close file
Note: See TracChangeset
for help on using the changeset viewer.
