New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 15472 – NEMO

Changeset 15472


Ignore:
Timestamp:
2021-11-04T17:35:49+01:00 (7 months ago)
Author:
vancop
Message:

Further fixes to VP rheology

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90

    r15403 r15472  
    4141   PUBLIC   ice_dyn_rhg_vp   ! called by icedyn_rhg.F90 
    4242 
     43   !! 
     44   INTEGER ::   nn_nvp              ! total number of VP iterations (n_out_vp*n_inn_vp)       
     45 
    4346   !! for convergence tests 
    4447   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 
    5152   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) 
    5361 
    5462   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15 
     
    132140      ! 
    133141      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  !  
    135143      INTEGER ::   ji_min, jj_min      ! 
    136144      INTEGER ::   nn_zebra_vp         ! number of zebra steps 
    137145 
    138       INTEGER ::   nn_nvp              ! total number of VP iterations (n_out_vp*n_inn_vp)       
    139146      ! 
    140147      REAL(wp) ::   zrhoco                                              ! rau0 * rn_cio 
     
    167174      REAL(wp), DIMENSION(jpi,jpj) ::   zu_c, zv_c                      ! "current" ice velocity (m/s), average of previous two OL iterates 
    168175      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 
    169177      REAL(wp), DIMENSION(jpi,jpj) ::   zuerr, zverr                    ! absolute U/Vvelocity difference between current and previous sub-iterates 
    170178      ! 
     
    197205      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small 
    198206      !! --- diags 
    199       REAL(wp) ::   zsig1, zsig2, zsig12, zdelta, z1_strength, zfac_x, zfac_y 
     207      REAL(wp)                     ::   zsig1, zsig2, zsig12, zdelta, z1_strength, zfac_x, zfac_y 
    200208      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12, zs12f           ! stress tensor components 
    201209      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p 
     
    204212      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw, zdiag_ymtrp_snw ! X/Y-component of snow mass transport (kg/s, SIMIP) 
    205213      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 
    206217                         
    207218      !!---------------------------------------------------------------------------------------------------------------------- 
     
    261272      ! Initialise convergence checks 
    262273      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) 
    266280          
    267281      ENDIF 
     
    417431      zv_c(:,:) = v_ice(:,:) 
    418432       
    419       jter = 0 
     433      i_inn_tot = 0 
    420434 
    421435      DO i_out = 1, nn_nout_vp 
    422436 
    423437         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 
    425442         ! Velocities used in the non linear terms are the average of the past two iterates 
    426443         ! u_it = 0.5 * ( u_{it-1} + u_{it-2}) 
    427444         ! 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(:,:) ) 
    430448          
    431449         !------------------------------------------------------------------------------! 
     
    525543         END DO 
    526544 
    527 !        zef(:,:) = 0. ! test 
    528           
    529545         CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. ) 
    530546         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zef'           , zef            ) ! MV DEBUG 
     
    553569                 
    554570                ! --- 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 done 
    556571                zCorU(ji,jj)         =   0.25_wp * r1_e1u(ji,jj) *  & 
    557572                           &             ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * zv_c(ji  ,jj) + e1v(ji  ,jj-1) * zv_c(ji  ,jj-1) )  & 
     
    606621         END DO 
    607622 
    608          ! TEST 
    609          ! zs2_rhsv(:,:) = 0._wp 
    610  
    611623         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs1_rhsu'      , zs1_rhsu       ) ! MV DEBUG 
    612624         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs2_rhsu'      , zs2_rhsu       ) ! MV DEBUG 
     
    633645            END DO 
    634646         END DO 
    635          ! 
    636          ! TEST 
    637          ! zs12_rhsv(:,:) = 0._wp 
    638647 
    639648         CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsu, 'F', 1. ) 
     
    822831            !--- mitgcm computes initial value of residual here... 
    823832 
    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(:,:) 
    832838 
    833839            IF ( ll_u_iterate .OR. ll_v_iterate )   THEN 
     
    10491055               ! MV OPT: if the number of iterations to convergence is really variable, and keep the convergence check 
    10501056               ! then we must optimize the use of the mpp_max, which is prohibitive                             
    1051                zuerr_max = 0._wp 
     1057               zuerr_max  = 0._wp 
    10521058                                
    10531059               IF ( ll_u_iterate .AND. MOD ( i_inn, nn_cvgchk_vp ) == 0 ) THEN 
     
    10571063                  DO jj = 2, jpj - 1 
    10581064                     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)  
    10601066                     END DO 
    10611067                  END DO 
    10621068                  zuerr_max = MAXVAL( zuerr ) 
    10631069                  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) 
    10661072                  IF ( i_inn > 1 .AND. zuerr_max > rn_uerr_max_vp ) THEN 
    10671073                      IF ( lwp ) WRITE(numout,*) " VP rheology error was too large : ", zuerr_max, " in outer U-iteration ", i_out, " after ", i_inn, " iterations, we stopped " 
     
    10911097                     END DO 
    10921098                  END DO 
    1093                    
    10941099                  zverr_max = MAXVAL( zverr ) 
    10951100                  CALL mpp_max( 'icedyn_rhg_evp', zverr_max )   ! max over the global domain - damned! 
    1096                    
     1101 
    10971102                  ! - Stop if error is too large ("safeguard against bad forcing" of original Zhang routine) 
    10981103                  IF ( i_inn > 1 .AND. zverr_max > rn_uerr_max_vp ) THEN 
     
    11161121            ! 
    11171122            !--------------------------------------------------------------------------------------- 
    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 
    11211133 
    11221134            IF ( lwp ) WRITE(numout,*) ' Done convergence tests ' 
     
    11451157      CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 
    11461158 
    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 solver 
    1149       IF( iom_use('vice_dbg' ) )   CALL iom_put( 'vice_dbg'   , v_ice    )                              ! ice velocity v after solver 
    1150              
    1151       !------------------------------------------------------------------------------! 
    1152       ! 
    1153       ! --- Convergence diagnostics  
    1154       ! 
    1155       !------------------------------------------------------------------------------! 
    11561159 
    11571160      IF( ln_rhg_chkcvg ) THEN 
    11581161           
    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 ) 
    11631168         
    11641169      ENDIF ! ln_rhg_chkcvg 
     
    14581463         DO jj = 2, jpj - 1 
    14591464            DO ji = 2, jpi - 1 
     1465 
    14601466               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             & 
    14611467                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    & 
     
    14641470                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              & 
    14651471                  &                  ) * r1_e1e2u(ji,jj) 
     1472 
    14661473               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             & 
    14671474                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    & 
     
    14701477                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              & 
    14711478                  &                  ) * r1_e1e2v(ji,jj) 
     1479 
    14721480            END DO 
    14731481         END DO 
     
    15261534    
    15271535    
    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      !! 
    15321543      !!---------------------------------------------------------------------- 
    15331544      !!                    ***  ROUTINE rhg_cvg_vp  *** 
     
    15441555      !! 
    15451556      !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)    
     1557      !! 
    15461558      !!---------------------------------------------------------------------- 
    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 
    15551575      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 
    15581583      REAL(wp), DIMENSION(jpi,jpj) ::   zu_res, zv_res, zvel2                         ! local arrays 
     1584      REAL(wp), DIMENSION(jpi,jpj) ::   zu_diff, zv_diff                              ! local arrays 
    15591585                                                                              
    15601586      CHARACTER(len=20) ::   clname 
    15611587      !!---------------------------------------------------------------------- 
    15621588 
     1589 
    15631590      IF( lwp ) THEN 
     1591 
    15641592         WRITE(numout,*) 
    15651593         WRITE(numout,*) 'rhg_cvg_vp : ice rheology convergence control' 
    15661594         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 
    15691604      ENDIF 
    15701605 
     1606      z1_pglob_area = 1._wp / pglob_area      ! inverse global ice area 
     1607 
    15711608      ! create file 
    1572       IF( kt == nit000 .AND. kiter == 1 ) THEN 
     1609      IF( kt == nit000 .AND. kitinntot == 1 ) THEN 
    15731610         ! 
    15741611         IF( lwm ) THEN 
     
    15821619            istatus = NF90_DEF_DIM( ncvgid, 'y'     , jpj, iy_dim ) 
    15831620 
    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 
    15911637            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) 
    15951638 
    15961639            istatus = NF90_ENDDEF(ncvgid) 
     
    16001643      ENDIF 
    16011644 
    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(:,:) ) 
    16151748 
    16161749      ELSE 
    16171750 
    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 
    16691849                   
    1670          !                                                ! ==================== ! 
    1671  
    1672       ! time 
    1673       it = ( kt - 1 ) * kitermax + kiter 
    1674  
    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 
    16761856      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 
    16911878         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 
    16921889 
    16931890         ! close file 
Note: See TracChangeset for help on using the changeset viewer.