Changeset 13591


Ignore:
Timestamp:
2020-10-14T16:38:22+02:00 (3 months ago)
Author:
vancop
Message:

some debug

Location:
NEMO/branches/2020/SI3-03_VP_rheology/src/ICE
Files:
2 edited

Legend:

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

    r13544 r13591  
    138138      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points 
    139139      REAL(wp), DIMENSION(jpi,jpj) ::   zdeltastar_t                    ! Delta* at T-points 
     140      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! Tension 
    140141      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017 
    141142      ! 
     
    170171      !! --- diags 
    171172      REAL(wp) ::   zsig1, zsig2, zsig12, zfac, z1_strength 
    172       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p, zten_i          
     173      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p 
    173174      !! --- SIMIP diags 
    174175      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 
     
    809810      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence 
    810811      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear 
     812      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta 
    811813      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength 
    812814 
     
    830832               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 
    831833               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                ! 1st stress invariant, aka average normal stress, aka negative pressure 
    832                zsig_II(ji,jj)   =   - SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 )   ! 2nd  ''       '', aka maximum shear stress 
     834               zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 )   ! 2nd  ''       '', aka maximum shear stress 
    833835                
    834836            END DO 
     
    852854      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 
    853855         ! 
    854          ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) )          
     856         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 
    855857         !          
    856858         DO jj = 2, jpj - 1 
     
    867869               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 
    868870               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                ! 1st stress invariant, aka average normal stress, aka negative pressure 
    869                zsig_II(ji,jj)   =   - SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 )   ! 2nd  ''       '', aka maximum shear stress 
     871               zsig_II(ji,jj)   =   SQRT ( zsig2 * zsig2 * 0.25_wp + zsig12 )     ! 2nd  ''       '', aka maximum shear stress 
    870872 
    871873               ! Normalized  principal stresses (used to display the ellipse) 
    872                z1_strength      =   1._wp / strength(ji,jj) 
    873                zsig1_p(ji,jj)   =   ( zsig_I(ji,jj)  - zsig_II(ji,jj) ) * z1_strength 
    874                zsig2_p(ji,jj)   =   ( zsig_II(ji,jj) + zsig_II(ji,jj) ) * z1_strength 
     874               z1_strength      =   1._wp / MAX ( 1.0, strength(ji,jj) ) 
     875               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 
     876               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 
    875877            END DO 
    876878         END DO 
     
    882884         CALL iom_put( 'sig2_pnorm' , zsig2_p )  
    883885 
    884          DEALLOCATE( zsig1_p , zsig2_p ) 
     886         DEALLOCATE( zsig1_p , zsig2_p , zsig_I , zsig_II ) 
    885887          
    886888      ENDIF 
  • NEMO/branches/2020/SI3-03_VP_rheology/src/ICE/icedyn_rhg_vp.F90

    r13544 r13591  
    5050   INTEGER ::   nvarid_vdif 
    5151   INTEGER ::   nvarid_mke 
     52   INTEGER ::   nvarid_usub, nvarid_vsub 
     53 
    5254   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15 
    5355 
     
    129131      LOGICAL ::   ll_u_iterate, ll_v_iterate   ! continue iteration or not 
    130132      ! 
    131       INTEGER ::   ji, jj, jn          ! dummy loop indices 
     133      INTEGER ::   ji, ji2, jj, jn          ! dummy loop indices 
    132134      INTEGER ::   jter, i_out, i_inn  !  
    133135      INTEGER ::   ji_min, jj_min      ! 
    134136      INTEGER ::   nn_zebra_vp         ! number of zebra steps 
     137 
     138      INTEGER ::   nn_thomas           ! tridiagonal system version (1=mitgcm, 2=martin) 
     139 
    135140      INTEGER ::   nn_nvp              ! total number of VP iterations (n_out_vp*n_inn_vp)       
    136141      ! 
     
    146151      REAL(wp) ::   zt12U, zt11U, zt22U, zt21U, zt122U, zt121U 
    147152      REAL(wp) ::   zt12V, zt11V, zt22V, zt21V, zt122V, zt121V 
    148       REAL(wp) ::   zAA3, zw, zuerr_max, zverr_max 
     153      REAL(wp) ::   zAA3, zw, ztau, zuerr_max, zverr_max 
    149154      ! 
    150155      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                          ! mask at F points for the ice 
     
    154159      ! 
    155160      REAL(wp), DIMENSION(jpi,jpj) ::   zdeltastar_t                    ! Delta* at T-points 
     161      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! Tension 
    156162      REAL(wp), DIMENSION(jpi,jpj) ::   zp_deltastar_t                  ! P/delta* at T points 
    157163      REAL(wp), DIMENSION(jpi,jpj) ::   zzt, zet                        ! Viscosity pre-factors at T points 
     
    182188      REAL(wp), DIMENSION(jpi,jpj) ::   zFU, zFU_prime, zBU_prime       ! Rearranged linear system coefficients, U equation 
    183189      REAL(wp), DIMENSION(jpi,jpj) ::   zFV, zFV_prime, zBV_prime       ! Rearranged linear system coefficients, V equation 
     190      REAL(wp), DIMENSION(jpi,jpj) ::   zCU_prime, zCV_prime            ! Rearranged linear system coefficients, V equation 
    184191!!!      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_bi, ztauy_bi              ! ice-OceanBottom stress at U-V points (landfast) 
    185192!!!      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast) 
     
    193200      !! --- diags 
    194201      REAL(wp) ::   zsig1, zsig2, zsig12, zdelta, z1_strength, zfac_x, zfac_y 
    195       REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12, zs12f, zten_i   ! stress tensor components 
     202      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12, zs12f           ! stress tensor components 
    196203      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p 
    197204      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   ztaux_oi, ztauy_oi 
     
    201208                         
    202209      !!---------------------------------------------------------------------------------------------------------------------- 
     210      ! DEBUG put all forcing terms to zero 
     211         ! air-ice drag 
     212         utau_ice(:,:) = 0._wp 
     213         vtau_ice(:,:) = 0._wp 
     214         ! coriolis 
     215         ff_t(:,:) = 0._wp 
     216         ! ice-ocean drag 
     217         rn_cio = 0._wp 
     218         ! ssh  
     219         ! done line 330 !!! dont forget to act there 
     220      ! END DEBUG 
    203221 
    204222      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_vp: VP sea-ice rheology (LSR solver)' 
     223      IF( lwp )   WRITE(numout,*) '-- ice_dyn_rhg_vp: VP sea-ice rheology (LSR solver)' 
    205224             
    206225      !------------------------------------------------------------------------------! 
     
    219238      END DO 
    220239       
    221       IF ( ln_zebra_vp ) THEN; nn_zebra_vp = 2; ELSE; nn_zebra_vp = 1; ENDIF  
    222        
     240      IF ( ln_zebra_vp ) THEN; nn_zebra_vp = 2 
     241                         ELSE; nn_zebra_vp = 1; ENDIF  
     242 
    223243      nn_nvp = nn_nout_vp * nn_ninn_vp ! maximum number of iterations 
     244 
     245      IF( lwp )   WRITE(numout,*) ' ln_zebra_vp : ', ln_zebra_vp 
     246      IF( lwp )   WRITE(numout,*) ' nn_zebra_vp : ', nn_zebra_vp 
     247      IF( lwp )   WRITE(numout,*) ' nn_nvp      : ', nn_nvp 
    224248       
    225249      zrhoco = rau0 * rn_cio  
     
    242266      ELSE                         ;   zkt = 0._wp 
    243267      ENDIF 
     268 
     269      zs1_rhsu(:,:) = 0._wp; zs2_rhsu(:,:) = 0._wp; zs1_rhsv(:,:) = 0._wp; zs2_rhsv(:,:) = 0._wp 
     270      zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp; 
     271      zAV(:,:) = 0._wp; zBV(:,:) = 0._wp; zCV(:,:) = 0._wp; zDV(:,:) = 0._wp; zEV(:,:) = 0._wp; 
     272      zrhsu(:,:) = 0._wp; zrhsv(:,:) = 0._wp 
     273      zf_rhsu(:,:) = 0._wp; zf_rhsv(:,:) = 0._wp 
    244274 
    245275      !------------------------------------------------------------------------------! 
     
    304334      !    non-embedded sea ice: use ocean surface for slope calculation 
    305335      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 
     336      zsshdyn(:,:) = 0._wp ! DEBUG CAREFUL !!! 
    306337       
    307338      DO jj = 2, jpj - 1 
     
    349380            ! Mask for lots of ice (1) or little ice (0) 
    350381            IF( zmassU <= zmmin .AND. za_iU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp 
    351             ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF 
     382            ELSE                                                     ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF 
    352383            IF( zmassV <= zmmin .AND. za_iV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp 
    353             ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF               
     384            ELSE                                                     ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF               
    354385 
    355386         END DO 
    356387      END DO    
     388 
     389      CALL iom_put( 'zmsk00x'    , zmsk00x  )   ! MV DEBUG 
     390      CALL iom_put( 'zmsk00y'    , zmsk00y  )   ! MV DEBUG 
     391      CALL iom_put( 'zmsk01x'    , zmsk01x  )   ! MV DEBUG 
     392      CALL iom_put( 'zmsk01y'    , zmsk01y  )   ! MV DEBUG 
     393      CALL iom_put( 'ztauy_ai'   , ztauy_ai )   ! MV DEBUG 
     394      CALL iom_put( 'ztaux_ai'   , ztaux_ai )   ! MV DEBUG 
     395      CALL iom_put( 'ztauy_ai'   , ztauy_ai )   ! MV DEBUG 
     396      CALL iom_put( 'zspgU'      , zspgU    )   ! MV DEBUG 
     397      CALL iom_put( 'zspgV'      , zspgV    )   ! MV DEBUG 
    357398             
    358399      !------------------------------------------------------------------------------! 
     
    368409 
    369410      DO i_out = 1, nn_nout_vp 
     411 
     412         IF( lwp )   WRITE(numout,*) ' outer loop  i_out : ', i_out 
    370413                
    371414         ! Velocities used in the non linear terms are the average of the past two iterates 
     
    400443            END DO 
    401444         END DO 
    402           
     445 
    403446         CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! MV TEST could be un-necessary according to Gurvan 
     447         CALL iom_put( 'zds'        , zds      )   ! MV DEBUG 
     448 
     449         IF( lwp )   WRITE(numout,*) ' outer loop  1a i_out : ', i_out 
    404450 
    405451         DO jj = 2, jpj - 1    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 
     
    441487         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. ) 
    442488 
     489         CALL iom_put( 'zzt'        , zzt      )   ! MV DEBUG 
     490         CALL iom_put( 'zet'        , zet      )   ! MV DEBUG 
     491         CALL iom_put( 'zp_deltastar_t', zp_deltastar_t ) ! MV DEBUG 
     492 
     493         IF( lwp )   WRITE(numout,*) ' outer loop  1b i_out : ', i_out 
     494 
    443495         DO jj = 1, jpj - 1 
    444496            DO ji = 1, jpi - 1 
     
    454506          
    455507         CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. ) 
     508         CALL iom_put( 'zef'           , zef            ) ! MV DEBUG 
     509         IF( lwp )   WRITE(numout,*) ' outer loop  1c i_out : ', i_out 
    456510 
    457511         !--------------------------------------------------- 
    458512         ! -- Ocean-ice drag and Coriolis RHS contributions 
    459513         !--------------------------------------------------- 
    460           
     514 
    461515         DO jj = 2, jpj - 1 
    462516             DO ji = 2, jpi - 1 
     
    489543             END DO 
    490544         END DO 
     545 
     546         IF( lwp )   WRITE(numout,*) ' outer loop  1d i_out : ', i_out 
     547         CALL lbc_lnk( 'icedyn_rhg_vp', zCwU, 'U', 1. ) 
     548         CALL lbc_lnk( 'icedyn_rhg_vp', zCwV, 'V', 1. ) 
     549 
     550         CALL iom_put( 'zCwU'          , zCwU           ) ! MV DEBUG 
     551         CALL iom_put( 'zCwV'          , zCwV           ) ! MV DEBUG 
     552         IF( lwp )   WRITE(numout,*) ' outer loop  1e i_out : ', i_out 
     553         CALL iom_put( 'zCorU'         , zCorU          ) ! MV DEBUG 
     554         CALL iom_put( 'zCorV'         , zCorV          ) ! MV DEBUG 
     555         IF( lwp )   WRITE(numout,*) ' outer loop  1f i_out : ', i_out 
    491556          
    492557         ! a priori, Coriolis and drag terms only affect diagonal or independent term of the linear system,  
     
    515580            END DO 
    516581         END DO 
     582 
     583         CALL iom_put( 'zs1_rhsu'      , zs1_rhsu       ) ! MV DEBUG 
     584         CALL iom_put( 'zs2_rhsu'      , zs2_rhsu       ) ! MV DEBUG 
     585         CALL iom_put( 'zs1_rhsv'      , zs1_rhsv       ) ! MV DEBUG 
     586         CALL iom_put( 'zs2_rhsv'      , zs2_rhsv       ) ! MV DEBUG 
    517587          
    518588         ! a priori, no lbclnk, because rhsu is only used in the inner domain 
     
    522592         ! My guess is that this is the way to enforce boundary conditions on strain rate tensor 
    523593 
     594         IF( lwp )   WRITE(numout,*) ' outer loop 2 i_out : ', i_out 
     595 
    524596         DO jj = 1, jpj - 1 
    525597            DO ji = 1, jpi - 1 
     
    533605            END DO 
    534606         END DO 
    535           
     607 
     608         CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsu, 'F', 1. ) 
     609         CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsv, 'F', 1. ) 
     610 
     611         CALL iom_put( 'zs12_rhsu'     , zs12_rhsu      ) ! MV DEBUG 
     612         CALL iom_put( 'zs12_rhsv'     , zs12_rhsv      ) ! MV DEBUG 
     613 
    536614         ! a priori, no lbclnk, because rhsu are only used in the inner domain 
    537615 
     
    555633            END DO 
    556634         END DO 
     635 
     636         CALL iom_put( 'zf_rhsu'       , zf_rhsu        ) ! MV DEBUG 
     637         CALL iom_put( 'zf_rhsv'       , zf_rhsv        ) ! MV DEBUG 
    557638          
    558639         !--------------------------- 
     
    571652         END DO 
    572653          
     654         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zrhsu, 'U', 1., zrhsv, 'V',  1. ) 
     655 
     656         CALL iom_put( 'zrhsu'         , zrhsu          ) ! MV DEBUG 
     657         CALL iom_put( 'zrhsv'         , zrhsv          ) ! MV DEBUG 
     658          
    573659         ! inner domain calculations -> no lbclnk 
     660 
     661         IF( lwp )   WRITE(numout,*) ' outer loop 4 i_out : ', i_out 
    574662      
    575663         !---------------------------------------------------------------------------------------! 
     
    662750            END DO 
    663751         END DO 
    664                 
     752 
     753         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zAU  , 'U', 1., zAV  , 'V',  1. ) 
     754         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zBU  , 'U', 1., zBV  , 'V',  1. ) 
     755         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCU  , 'U', 1., zCV  , 'V',  1. ) 
     756         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zDU  , 'U', 1., zDV  , 'V',  1. ) 
     757         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zEU  , 'U', 1., zEV  , 'V',  1. ) 
     758                
     759         CALL iom_put( 'zAU'           , zAU            ) ! MV DEBUG 
     760         CALL iom_put( 'zBU'           , zBU            ) ! MV DEBUG 
     761         CALL iom_put( 'zCU'           , zCU            ) ! MV DEBUG 
     762         CALL iom_put( 'zDU'           , zDU            ) ! MV DEBUG 
     763         CALL iom_put( 'zEU'           , zEU            ) ! MV DEBUG 
     764         CALL iom_put( 'zAV'           , zAV            ) ! MV DEBUG 
     765         CALL iom_put( 'zBV'           , zBV            ) ! MV DEBUG 
     766         CALL iom_put( 'zCV'           , zCV            ) ! MV DEBUG 
     767         CALL iom_put( 'zDV'           , zDV            ) ! MV DEBUG 
     768         CALL iom_put( 'zEV'           , zEV            ) ! MV DEBUG 
     769 
    665770      !------------------------------------------------------------------------------! 
    666771      ! 
     
    673778         ll_v_iterate = .TRUE. 
    674779 
     780         ! zBU 
     781         ! it is probably the way zBU is used that is problematic... 
     782         ! because zBU itself does not have the zipper!!! 
     783         ! BUG only occurs when zBU is different from zero 
     784         ! H1 - the jj loop does not work well with jpj = 149 ? 
     785         ! zBU(:,:) = 0._wp; zBV(:,:) = 0._wp 
     786 
    675787         DO i_inn = 1, nn_ninn_vp ! inner loop iterations 
    676           
    677             ! mitgcm computes initial value of residual 
    678             jter = jter + 1 
    679             l_full_nf_update = jter == nn_nvp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 
    680  
    681             IF ( ll_u_iterate .OR. ll_v_iterate ) THEN 
    682              
    683  
    684                zu_b(:,:) = u_ice(:,:) ! velocity at previous sub-iterate 
    685                zv_b(:,:) = v_ice(:,:) 
     788 
     789            IF( lwp )   WRITE(numout,*) ' inner loop 1 i_inn : ', i_inn 
     790          
     791            !--- mitgcm computes initial value of residual here... 
     792 
     793            jter             = jter + 1 
     794            ! l_full_nf_update = jter == nn_nvp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 
     795 
     796            zu_b(:,:) = u_ice(:,:) ! velocity at previous sub-iterate 
     797            zv_b(:,:) = v_ice(:,:) 
     798 
     799            zFU(:,:)       = 0._wp ; zFU_prime(:,:) = 0._wp ; zBU_prime(:,:) = 0._wp 
     800            zFV(:,:)       = 0._wp ; zFV_prime(:,:) = 0._wp ; zBV_prime(:,:) = 0._wp 
     801 
     802            IF ( ll_u_iterate .OR. ll_v_iterate )   THEN 
    686803 
    687804                                           ! ---------------------------- ! 
     
    689806                                           ! ---------------------------- ! 
    690807       
     808                  ! Thomas Algorithm  for tridiagonal solver 
    691809                  ! what follows could be subroutinized... 
    692810                   
     
    698816                     ELSE                  ;   jj_min = 3 
    699817                     ENDIF 
    700           
    701                      zFU(:,:)       = 0._wp 
    702                      zFU_prime(:,:) = 0._wp 
    703                      zBU_prime(:,:) = 0._wp 
    704           
     818 
     819                     IF ( lwp ) WRITE(numout,*) ' Into the U-zebra loop at step jn = ', jn, ', with jj_min = ', jj_min 
     820 
     821                     !---------------------------------------------------------- 
     822                     ! -- Boundary condition (link with neighbouring processor) 
     823                     !---------------------------------------------------------- 
     824 
    705825                     DO jj = jj_min, jpj - 1, nn_zebra_vp 
    706826 
    707                         !------------------------------------------- 
    708                         ! -- Tridiagonal system solver for each row 
    709                         !------------------------------------------- 
    710827                        ! 
    711828                        ! MV - I am in doubts whether the way I coded it is reproducible - ask Gurvan 
     
    722839                           IF ( ji == 2 )          zAA3 = zAA3 - zAU(ji,jj) * u_ice(ji-1,jj) 
    723840                           IF ( ji == jpi - 1 )    zAA3 = zAA3 - zCU(ji,jj) * u_ice(ji+1,jj) 
    724       
     841 
    725842                           ! right hand side 
    726                            zFU(ji,jj) = ( zrhsu(ji,jj) &                                   ! right-hand side terms 
    727                                &      +   zAA3                                           & ! boundary condition translation 
    728                                &      +   zDU(ji,jj) * u_ice(ji,jj-1)                    & ! internal force, j-1 
    729                                &      +   zEU(ji,jj) * u_ice(ji,jj+1) ) * umask(ji,jj,1)   ! internal force, j+1 
    730                      
     843                           zFU(ji,jj) = ( zrhsu(ji,jj) &                                      ! right-hand side terms 
     844                               &      +   zAA3         &                                      ! boundary condition translation 
     845                               &      +   zDU(ji,jj) * u_ice(ji,jj-1)   &                    ! internal force, j-1 
     846                               &      +   zEU(ji,jj) * u_ice(ji,jj+1) ) * umask(ji,jj,1)      ! internal force, j+1 
     847 
    731848                        END DO 
     849 
     850                     END DO 
     851 
     852                     ! CALL lbc_lnk( 'icedyn_rhg_vp', zFU, 'U', 1. ) 
    732853                         
    733                         ! - Thomas Algorithm  
    734                         ! (MV: I chose a seemingly more efficient form of the algorithm than in mitGCM - not sure) 
    735                         ! Forward sweep 
    736                         DO ji = 3, jpi - 1 
    737                            zw         = zAU(ji,jj) / zBU(ji-1,jj) 
    738                            zBU_prime(ji,jj) = zBU(ji,jj) - zw * zCU(ji-1,jj) 
    739                            zFU_prime(ji,jj) = zFU(ji,jj) - zw * zFU(ji-1,jj) 
     854                     !--------------- 
     855                     ! Forward sweep 
     856                     !--------------- 
     857                     nn_thomas = 2 ! 1 = martin, 2 = mitgcm 
     858 
     859                     IF ( nn_thomas == 1 ) THEN ! MV version seemingly more than mitGCM algorithm 
     860 
     861                        DO jj = jj_min, jpj - 1, nn_zebra_vp 
     862 
     863                           DO ji = 3, jpi - 1 
     864 
     865                              zfac             = SIGN( 1._wp , zBU(ji-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU(ji-1,jj) ) - epsi20 ) ) 
     866                              zw               = zfac * zAU(ji,jj) / MAX ( ABS( zBU(ji-1,jj) ) , epsi20 )  
     867                              zBU_prime(ji,jj) = zBU(ji,jj) - zw * zCU(ji-1,jj) 
     868                              zFU_prime(ji,jj) = zFU(ji,jj) - zw * zFU(ji-1,jj) 
     869 
     870                           END DO 
     871                         
    740872                        END DO 
    741              
    742                         ! Backward sweep 
    743                         ! MV I don't see how this will be reproducible 
    744                         u_ice(jpi-1,jj)     = zFU_prime(jpi-1,jj) / zBU_prime(jpi-1,jj) * umask(jpi-1,jj,1)                     ! do last row first 
    745                         DO ji = jpi-2 , 2, -1 ! all other rows    !   MV OPT: could be turned into forward loop (by substituting ji) 
    746                            u_ice(ji,jj)    = zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji,jj+1) * umask(ji,jj,1) / zBU_prime(ji,jj) ! 
    747                         END DO             
    748              
    749                         !---------------        
    750                         ! -- Relaxation 
    751                         !--------------- 
     873 
     874                     ELSE ! nn_thomas == 2      ! mitGCM algorithm 
     875 
     876                        DO jj = jj_min, jpj - 1, nn_zebra_vp 
     877 
     878                           ! ji = 2 
     879                           zw                   =   zBU(2,jj) 
     880                           zfac                 =   SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) ) 
     881                           zw                   =   zfac * 1.0_wp /  MAX ( ABS( zw ) , epsi20 ) 
     882                           zCU_prime(2,jj)      =   zw * zCU(2,jj) 
     883                           zFU_prime(2,jj)      =   zw * zFU(2,jj) 
     884 
     885                           DO ji = 3, jpi - 1 
     886 
     887                              zw                =   zBU(ji,jj) - zAU(ji,jj) * zCU_prime(ji-1,jj) 
     888                              zfac              =   SIGN( 1._wp , zw ) * MAX( 0._wp , SIGN( 1._wp , ABS(zw) - epsi20 ) ) 
     889                              zw                =   zfac * 1.0_wp /  MAX ( ABS( zw ) , epsi20 ) 
     890                              zCU_prime(ji,jj)  =   zw * zCU(ji,jj)  
     891                              zFU_prime(ji,jj)  =   zw * ( zFU(ji,jj) - zAU(ji,jj) * zFU_prime(ji-1,jj) ) 
     892 
     893                           END DO 
     894 
     895                        END DO 
     896 
     897                     ENDIF 
     898 
     899                     ! CALL lbc_lnk_multi( 'icedyn_rhg_vp', zBU_prime , 'U', 1., zFU_prime, 'U', 1. ) 
     900 
     901                     !---------------- 
     902                     ! Backward sweep 
     903                     !---------------- 
     904                     IF ( nn_thomas == 1 ) THEN ! MV version 
     905 
     906                        DO jj = jj_min, jpj - 1, nn_zebra_vp 
     907 
     908                           ! last row  
     909                           zfac = SIGN( 1._wp , zBU_prime(jpi-1,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(jpi-1,jj) ) - epsi20 ) ) 
     910                           u_ice(jpi-1,jj)     = zfac * zFU_prime(jpi-1,jj) / MAX( ABS ( zBU_prime(jpi-1,jj) ) , epsi20 ) &  
     911                                  &            * umask(jpi-1,jj,1) 
     912 
     913                           DO ji2 = 2 , jpi-2 ! all other rows        
     914                          !DO ji = jpi-2 , 2, -1 ! all other rows    !  ---> original backward loop 
     915                              ji =  jpi - ji2 
     916                              ! ji2 = 2       -> ji = jpi - 2 
     917                              ! ji2 = jpi - 2 -> ji = 2 
     918                              zfac = SIGN( 1._wp , zBU_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBU_prime(ji,jj) ) - epsi20 ) ) 
     919                              u_ice(ji,jj)    = zfac * ( zFU_prime(ji,jj) - zCU(ji,jj) * u_ice(ji+1,jj) * umask(ji,jj,1) ) & ! this 
     920                                   ! line is guilty 
     921                                              &                  / MAX ( ABS ( zBU_prime(ji,jj) ) , epsi20 )  
     922                           END DO             
     923 
     924                        END DO 
     925 
     926                     ELSE ! nn_thomas == 2 ! mitGCM version 
     927 
     928                        DO jj = jj_min, jpj - 1, nn_zebra_vp 
     929 
     930                           u_ice(jpi-1,jpj)   =   zFU_prime(jpi-1,jpj) 
     931 
     932                           DO ji2 = 2 , jpi-2 ! all other rows       !   MV OPT: could be turned into forward loop (by substituting ji) 
     933                              ji =  jpi - ji2 
     934                              u_ice(ji,jj)    =   zFU_prime(ji,jj) - zCU_prime(ji,jj) * u_ice(ji+1,jj) 
     935                           END DO 
     936                      
     937                        END DO 
     938 
     939                     ENDIF 
     940 
     941                     !---------------        
     942                     ! -- Relaxation 
     943                     !--------------- 
     944 
     945                     DO jj = jj_min, jpj - 1, nn_zebra_vp 
     946             
    752947                        DO ji = 2, jpi - 1     
    753948                         
     
    755950                            
    756951                           ! velocity masking for little-ice and no-ice cases 
    757                            ! put 0.01 of ocean current, appropriate choice to avoid too much spreading of the ice edge 
    758952                           u_ice(ji,jj) =   zmsk00x(ji,jj)                                        &  
    759953                              &         * (           zmsk01x(ji,jj)   * u_ice(ji,jj)             & 
    760                                           + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp    ) 
     954                              &           + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp    ) * umask(ji,jj,1) 
    761955                            
    762956                        END DO 
    763957 
    764958                     END DO ! jj 
    765  
     959                      
    766960                  END DO ! zebra loop 
    767                    
     961 
    768962               ENDIF !   ll_u_iterate 
    769                 
     963 
    770964               !                           ! ---------------------------- ! 
    771965               IF ( ll_v_iterate ) THEN    ! --- Solve for V-velocity --- ! 
    772             !                              ! ---------------------------- ! 
     966               !                           ! ---------------------------- ! 
    773967                                           
    774968                  ! MV OPT: what follows could be subroutinized... 
     
    779973                     ELSE                  ;   ji_min = 3 
    780974                     ENDIF 
    781           
    782                      zFV(:,:)       = 0._wp 
    783                      zFV_prime(:,:) = 0._wp 
    784                      zBV_prime(:,:) = 0._wp 
    785  
     975 
     976                     IF ( lwp ) WRITE(numout,*) ' Into the V-zebra loop at step jn = ', jn, ', with ji_min = ', ji_min 
     977          
    786978                     DO ji = ji_min, jpi - 1, nn_zebra_vp  
    787979                         
     
    800992                           ! see Zhang and Hibler, 1997, Appendix B 
    801993                           zAA3 = 0._wp 
    802                            IF ( jj .EQ. 2 )       zAA3 = zAA3 - zAV(ji,jj) * v_ice(ji,jj-1) 
    803                            IF ( jj .EQ. jpj - 1 ) zAA3 = zAA3 - zCV(ji,jj) * v_ice(ji,jj+1) 
     994                           IF ( jj == 2 )       zAA3 = zAA3 - zAV(ji,jj) * v_ice(ji,jj-1) 
     995                           IF ( jj == jpj - 1 ) zAA3 = zAA3 - zCV(ji,jj) * v_ice(ji,jj+1) 
    804996      
    805997                           ! right hand side 
     
    8081000                               &        + zDV(ji,jj) * v_ice(ji-1,jj)                    & ! internal force, j-1 
    8091001                               &        + zEV(ji,jj) * v_ice(ji+1,jj) ) * vmask(ji,jj,1)   ! internal force, j+1 
    810                      
     1002 
    8111003                        END DO 
    812         
     1004 
    8131005                        ! --- Thomas Algorithm  
    8141006                        ! (MV: I chose a seemingly more efficient form of the algorithm than in mitGCM - not sure) 
    8151007                        ! Forward sweep 
    8161008                        DO jj = 3, jpj - 1 
    817                            zw         = zAV(ji,jj) / zBV(ji,jj-1) 
     1009                           zfac             = SIGN( 1._wp , zBV(ji,jj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj-1) ) - epsi20 ) ) 
     1010                           zw               = zfac * zAV(ji,jj) / MAX ( ABS( zBV(ji,jj-1) ) , epsi20 ) 
    8181011                           zBV_prime(ji,jj) = zBV(ji,jj) - zw * zCV(ji,jj-1) 
    8191012                           zFV_prime(ji,jj) = zFV(ji,jj) - zw * zFV(ji,jj-1)  
     
    8211014             
    8221015                        ! Backward sweep 
    823                         v_ice(ji,jpj-1)  = zFV_prime(ji,jpj-1) / zBV_prime(ji,jpj-1) * vmask(ji,jj,jpj-1)  ! last row 
     1016                        zfac = SIGN( 1._wp , zBV_prime(ji,jpj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jpj-1) ) - epsi20 ) ) 
     1017                        v_ice(ji,jpj-1)  = zfac * zFV_prime(ji,jpj-1) / MAX ( ABS(zBV_prime(ji,jpj-1) ) , epsi20 ) &  
     1018                               &         * vmask(ji,jpj-1,1)  ! last row 
     1019 
    8241020                        DO jj = jpj-2, 2, -1 ! can be turned into forward row by substituting jj if needed 
    825                            v_ice(ji,jj)   = zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) * vmask(ji,jj,1) / zBV_prime(ji,jj) 
     1021                           zfac = SIGN( 1._wp , zBV_prime(ji,jj) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV_prime(ji,jj) ) - epsi20 ) ) 
     1022                           v_ice(ji,jj)   = zfac * ( zFV_prime(ji,jj) - zCV(ji,jj) * v_ice(ji,jj+1) * vmask(ji,jj,1) ) & 
     1023                               &                             / MAX ( ABS( zBV_prime(ji,jj) ) , epsi20 ) 
    8261024                        END DO             
    827              
     1025 
    8281026                        !---------------        
    8291027                        ! -- Relaxation 
     
    8351033                           v_ice(ji,jj) =   zmsk00y(ji,jj)                                        &  
    8361034                              &         * (           zmsk01y(ji,jj)   * v_ice(ji,jj)             & 
    837                                           + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp    ) 
     1035                              &           + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp    ) * vmask(ji,jj,1) 
    8381036 
    8391037                        END DO 
     
    8441042                                     
    8451043               ENDIF !   ll_v_iterate 
    846                 
    847                IF ( ( ll_u_iterate .OR. ll_v_iterate ) .OR. jter == nn_nvp ) CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 
     1044 
     1045               CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 
    8481046                               
    8491047               !-------------------------------------------------------------------------------------- 
     
    8561054               ! MV OPT: if the number of iterations to convergence is really variable, and keep the convergence check 
    8571055               ! then we must optimize the use of the mpp_max, which is prohibitive                             
     1056               zuerr_max = 0._wp 
    8581057                                
    8591058               IF ( ll_u_iterate .AND. MOD ( i_inn, nn_cvgchk_vp ) == 0 ) THEN 
     
    8861085               ! on V 
    8871086               !------ 
     1087               zverr_max = 0._wp 
    8881088                
    8891089               IF ( ll_v_iterate .AND. MOD ( i_inn, nn_cvgchk_vp ) == 0 ) THEN 
     
    9191119               ! 
    9201120               !--------------------------------------------------------------------------------------- 
    921                IF( ln_rhg_chkcvg ) CALL rhg_cvg_vp( kt, jter, nn_nvp, u_ice, v_ice, zmt, zuerr_max, zverr_max, zglob_area, & 
     1121               IF( ln_rhg_chkcvg .AND. MOD ( i_inn - 1, nn_cvgchk_vp ) == 0 ) CALL rhg_cvg_vp( kt, jter, nn_nvp, u_ice, v_ice, zmt, zuerr_max, zverr_max, zglob_area, & 
    9221122                          &                         zrhsu, zAU, zBU, zCU, zDU, zEU, zrhsv, zAV, zBV, zCV, zDV, zEV ) 
    923                 
    924              
     1123 
     1124               IF ( lwp ) WRITE(numout,*) ' Done convergence tests ' 
     1125                
    9251126            ENDIF ! ---    end ll_u_iterate or ll_v_iterate 
    926              
     1127 
    9271128         END DO ! i_inn, end of inner loop 
    9281129 
     1130      END DO ! End of outer loop (i_out) ============================================================================================= 
     1131 
     1132      IF ( lwp ) WRITE(numout,*) ' We are out of outer loop ' 
     1133 
     1134      CALL iom_put( 'zFU'           , zFU            ) ! MV DEBUG 
     1135      CALL iom_put( 'zBU_prime'     , zBU_prime      ) ! MV DEBUG 
     1136      CALL iom_put( 'zFU_prime'     , zFU_prime      ) ! MV DEBUG 
     1137 
     1138      CALL iom_put( 'zFV'           , zFV            ) ! MV DEBUG 
     1139      CALL iom_put( 'zBV_prime'     , zBV_prime      ) ! MV DEBUG 
     1140      CALL iom_put( 'zFV_prime'     , zFV_prime      ) ! MV DEBUG 
     1141 
     1142      CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 
     1143 
     1144      IF ( lwp ) WRITE(numout,*) ' We are about to output uice_dbg ' 
     1145      IF( iom_use('uice_dbg' ) )   CALL iom_put( 'uice_dbg'   , u_ice    )                              ! ice velocity u after solver 
     1146      IF( iom_use('vice_dbg' ) )   CALL iom_put( 'vice_dbg'   , v_ice    )                              ! ice velocity v after solver 
     1147             
    9291148      !------------------------------------------------------------------------------! 
    9301149      ! 
    931       ! 6) Mask final velocities 
     1150      ! --- Convergence diagnostics  
    9321151      ! 
    9331152      !------------------------------------------------------------------------------! 
    9341153 
    935       END DO ! End of outer loop (i_out) ============================================================================================= 
     1154      IF( ln_rhg_chkcvg ) THEN 
     1155           
     1156         IF( iom_use('uice_cvg')  ) THEN 
     1157            CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_b(:,:) ) * umask(:,:,1) , &                ! ice velocity difference at last iteration 
     1158                  &                        ABS( v_ice(:,:) - zv_b(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 
     1159         ENDIF    
     1160         
     1161      ENDIF ! ln_rhg_chkcvg 
     1162 
     1163      ! MV DEBUG test - replace ice velocity by ocean current to give the model the means to go ahead 
     1164      DO jj = 2, jpj - 1 
     1165         DO ji = 2, jpi - 1    
     1166 
     1167             u_ice(ji,jj) =   zmsk00x(ji,jj)                               &  
     1168      &         * (           zmsk01x(ji,jj)   * u_oce(ji,jj) * 0.01_wp    & 
     1169                  + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp    ) 
     1170 
     1171             v_ice(ji,jj) =   zmsk00y(ji,jj)                               &  
     1172      &         * (           zmsk01y(ji,jj)   * v_oce(ji,jj) * 0.01_wp    & 
     1173                  + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp    ) 
     1174 
     1175         END DO 
     1176      END DO 
     1177 
     1178      CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 
     1179 
     1180      ! END DEBUG 
    9361181 
    9371182      !------------------------------------------------------------------------------! 
     
    9421187      ! 
    9431188      ! MV OPT: subroutinize ? 
    944              
     1189 
    9451190      DO jj = 1, jpj - 1 
    9461191         DO ji = 1, jpi - 1 
     
    10871332      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence 
    10881333      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear 
     1334      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta 
    10891335      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength 
    10901336 
     
    11231369      ! Recommendation 2 : Need to use deformations at PREVIOUS iterate for viscosities (see p. 1765) 
    11241370      ! R2 means we need to recompute stresses 
     1371 
    11251372      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN 
    11261373         ! 
    1127          ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) )          
     1374         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) ) 
    11281375         !          
    11291376         DO jj = 2, jpj - 1 
     
    11421389 
    11431390               ! Normalized  principal stresses (used to display the ellipse) 
    1144                z1_strength      =   1._wp / strength(ji,jj) 
    1145                zsig1_p(ji,jj)   =   ( zsig_I(ji,jj)  + zsig_II(ji,jj) ) * z1_strength 
    1146                zsig2_p(ji,jj)   =   ( zsig_II(ji,jj) - zsig_II(ji,jj) ) * z1_strength 
     1391               z1_strength      =   1._wp / MAX ( 1._wp , strength(ji,jj) ) 
     1392               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength 
     1393               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength 
    11471394            END DO 
    11481395         END DO 
     
    11541401         CALL iom_put( 'sig2_pnorm' , zsig2_p )  
    11551402 
    1156          DEALLOCATE( zsig1_p , zsig2_p ) 
     1403         DEALLOCATE( zsig1_p , zsig2_p , zsig_I , zsig_II ) 
    11571404          
    11581405      ENDIF 
     
    12531500      ENDIF 
    12541501 
    1255       ! --- Convergence diagnostics --- ! 
    1256       IF( ln_rhg_chkcvg ) THEN 
    1257            
    1258          IF( iom_use('uice_cvg') ) THEN 
    1259             CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_b(:,:) ) * umask(:,:,1) , & 
    1260                   &                        ABS( v_ice(:,:) - zv_b(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 
    1261          ENDIF    
    1262          
    1263       ENDIF ! ln_rhg_chkcvg 
    1264  
    1265       ! 
    12661502      DEALLOCATE( zmsk00, zmsk15 ) 
    12671503 
     
    12851521      !!                   - residuals in U, V and UV-mean taken as square-root of area-weighted mean square residual (u_res, v_res, vel_res, N/m2) 
    12861522      !!                   - mean kinetic energy (mke_ice, J/m2) 
    1287       !!                    
    12881523      !! 
    12891524      !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)    
     
    12961531      REAL(wp), DIMENSION(:,:), INTENT(in) ::   prhsv, pAV, pBV, pCV, pDV, pEV 
    12971532      !! 
    1298       INTEGER           ::   it, idtime, istatus 
     1533      INTEGER           ::   it, idtime, istatus, ix_dim, iy_dim 
    12991534      INTEGER           ::   ji, jj          ! dummy loop indices 
    13001535      REAL(wp)          ::   zveldif, zu_res_mean, zv_res_mean, zvelres, zmke, zu, zv ! local scalars 
     1536      REAL(wp)          ::   z1_pglob_area 
    13011537      REAL(wp), DIMENSION(jpi,jpj) ::   zu_res, zv_res, zvel2                         ! local arrays 
    13021538                                                                              
     
    13041540      !!---------------------------------------------------------------------- 
    13051541 
     1542      IF( lwp ) THEN 
     1543         WRITE(numout,*) 
     1544         WRITE(numout,*) 'rhg_cvg_vp : ice rheology convergence control' 
     1545         WRITE(numout,*) '~~~~~~~~~~~' 
     1546         WRITE(numout,*) ' kiter    =  : ', kiter 
     1547         WRITE(numout,*) ' kitermax =  : ', kitermax 
     1548      ENDIF 
     1549 
    13061550      ! create file 
    13071551      IF( kt == nit000 .AND. kiter == 1 ) THEN 
    13081552         ! 
    1309          IF( lwp ) THEN 
    1310             WRITE(numout,*) 
    1311             WRITE(numout,*) 'rhg_cvg_vp : ice rheology convergence control' 
    1312             WRITE(numout,*) '~~~~~~~~~~~' 
    1313          ENDIF 
    1314          ! 
    13151553         IF( lwm ) THEN 
     1554 
    13161555            clname = 'ice_cvg.nc' 
    13171556            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 
    13181557            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) 
     1558 
    13191559            istatus = NF90_DEF_DIM( ncvgid, 'time'  , NF90_UNLIMITED, idtime ) 
    1320             istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE   , (/ idtime /), nvarid_ucvg ) 
     1560            istatus = NF90_DEF_DIM( ncvgid, 'x'     , jpi, ix_dim ) 
     1561            istatus = NF90_DEF_DIM( ncvgid, 'y'     , jpj, iy_dim ) 
     1562 
     1563            istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE   , (/ idtime /), nvarid_ucvg ) ! the name uice_cvg sucks, no ? 
     1564            ! i suggest vel_dif instead 
    13211565            istatus = NF90_DEF_VAR( ncvgid, 'u_res', NF90_DOUBLE   , (/ idtime /), nvarid_ures ) 
    13221566            istatus = NF90_DEF_VAR( ncvgid, 'v_res', NF90_DOUBLE   , (/ idtime /), nvarid_vres ) 
     
    13251569            istatus = NF90_DEF_VAR( ncvgid, 'v_dif', NF90_DOUBLE   , (/ idtime /), nvarid_vdif ) 
    13261570            istatus = NF90_DEF_VAR( ncvgid, 'mke_ice', NF90_DOUBLE   , (/ idtime /), nvarid_mke ) 
     1571 
     1572            istatus = NF90_DEF_VAR( ncvgid, 'usub'    , NF90_DOUBLE, (/ ix_dim, iy_dim, idtime /), nvarid_usub) ! --> u-velocity in sub-iterations 
     1573            istatus = NF90_DEF_VAR( ncvgid, 'vsub'    , NF90_DOUBLE, (/ ix_dim, iy_dim, idtime /), nvarid_vsub) ! --> v-velocity in sub-iterations 
     1574 
    13271575            istatus = NF90_ENDDEF(ncvgid) 
     1576 
    13281577         ENDIF 
    13291578         ! 
    13301579      ENDIF 
    13311580 
    1332       ! time 
    1333       it = ( kt - 1 ) * kitermax + kiter 
     1581      IF ( lwp ) WRITE(numout,*) ' File created ' 
    13341582 
    13351583      ! --- Max absolute velocity difference with previous iterate (zveldif) 
    1336       ! EVP code       
    1337 !      IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) 
    1338 !         zveldif = 0._wp 
    1339 !      ELSE 
    1340 !         DO jj = 1, jpj 
    1341 !            DO ji = 1, jpi 
    1342 !               zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 
    1343 !                  &               ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 
    1344 !            END DO 
    1345 !         END DO 
    1346 !         zveldif = MAXVAL( zres ) 
    1347 !         CALL mpp_max( 'icedyn_rhg_vp', zveldif )   ! max over the global domain 
    1348 !      ENDIF 
    1349       ! VP code 
    13501584      zveldif = MAX( puerr_max, pverr_max ) ! velocity difference with previous iterate, should nearly be equivalent to evp code  
    13511585                                            ! if puerrmask and pverrmax are masked at 15% (TEST) 
    1352        
     1586 
     1587      IF ( lwp ) WRITE(numout,*) ' zfeldif : ', zveldif 
     1588       
     1589      ! ---  Mean residual and kinetic energy 
     1590      IF ( kiter == 1 ) THEN 
     1591 
     1592              zu_res_mean = 0._wp 
     1593              zv_res_mean = 0._wp 
     1594              zvelres     = 0._wp 
     1595              zmke        = 0._wp 
     1596 
     1597      ELSE 
     1598 
    13531599      ! -- Mean residual (N/m^2), zu_res_mean 
    13541600      ! Here we take the residual of the linear system (N/m^2),  
     
    13561602      ! Local residual r = Ax - B expresses to which extent the momentum balance is verified  
    13571603      ! i.e., how close we are to a solution 
     1604 
     1605      IF ( lwp ) WRITE(numout,*) ' TEST 1 '  
     1606 
     1607      z1_pglob_area = 1._wp / pglob_area 
     1608 
    13581609      DO jj = 2, jpj - 1 
    13591610         DO ji = 2, jpi - 1                                       
    1360             zu_res(ji,jj)  = prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1) & 
    1361                &           - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) 
     1611            zu_res(ji,jj)  = ( prhsu(ji,jj) + pDU(ji,jj) * pu(ji,jj-1) + pEU(ji,jj) * pu(ji,jj+1)              & 
     1612               &             - pAU(ji,jj) * pu(ji-1,jj) - pBU(ji,jj) * pu(ji,jj) - pCU(ji,jj) * pu(ji+1,jj) ) 
    13621613                            
    1363             zv_res(ji,jj)  = prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj) & 
    1364                &           - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1)                            
     1614            zv_res(ji,jj)  = ( prhsv(ji,jj) + pDV(ji,jj) * pv(ji-1,jj) + pEV(ji,jj) * pv(ji+1,jj)               & 
     1615               &             - pAV(ji,jj) * pv(ji,jj-1) - pBV(ji,jj) * pv(ji,jj) - pCV(ji,jj) * pv(ji,jj+1) ) 
     1616 
     1617            zu_res(ji,jj)  = SQRT( zu_res(ji,jj) * zu_res(ji,jj) ) * umask(ji,jj,1) * e1e2u(ji,jj) * z1_pglob_area 
     1618            zv_res(ji,jj)  = SQRT( zv_res(ji,jj) * zv_res(ji,jj) ) * vmask(ji,jj,1) * e1e2v(ji,jj) * z1_pglob_area 
     1619 
    13651620          END DO 
    13661621       END DO                   
    1367        zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) * zu_res(:,:) * e1e2u(:,:) * umask(:,:,1) ) 
    1368        zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) * zv_res(:,:) * e1e2v(:,:) * vmask(:,:,1) ) 
    1369        zu_res_mean = SQRT( zu_res_mean / pglob_area ) 
    1370        zv_res_mean = SQRT( zv_res_mean / pglob_area ) 
     1622 
     1623       IF ( lwp ) WRITE(numout,*) ' TEST 2 '  
     1624       zu_res_mean = glob_sum( 'ice_rhg_vp', zu_res(:,:) ) 
     1625       zv_res_mean = glob_sum( 'ice_rhg_vp', zv_res(:,:) ) 
     1626       IF ( lwp ) WRITE(numout,*) ' TEST 3 '  
    13711627       zvelres     = 0.5_wp * ( zu_res_mean + zv_res_mean ) 
     1628 
     1629       IF ( lwp ) WRITE(numout,*) ' TEST 4 '  
    13721630                                           
    13731631       ! -- Global mean kinetic energy per unit area (J/m2)   
     
    13801638       END DO 
    13811639        
    1382        zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) 
     1640         IF ( lwp ) WRITE(numout,*) ' TEST 5 '  
     1641       zmke = 0.5_wp * glob_sum( 'ice_rhg_vp', pmt(:,:) * e1e2t(:,:) * zvel2(:,:) ) / pglob_area 
     1642         IF ( lwp ) WRITE(numout,*) ' TEST 6 '  
     1643 
     1644      ENDIF 
    13831645                   
    13841646         !                                                ! ==================== ! 
    13851647 
     1648      ! time 
     1649      it = ( kt - 1 ) * kitermax + kiter 
     1650 
     1651 
    13861652      IF( lwm ) THEN 
    13871653         ! write variables 
    1388             istatus = NF90_PUT_VAR( ncvgid, nvarid_ucvg, (/zveldif/), (/it/), (/1/) ) 
    1389             istatus = NF90_PUT_VAR( ncvgid, nvarid_ures, (/zu_res_mean/), (/it/), (/1/) ) 
    1390             istatus = NF90_PUT_VAR( ncvgid, nvarid_vres, (/zv_res_mean/), (/it/), (/1/) ) 
    1391             istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvelres/), (/it/), (/1/) ) 
    1392             istatus = NF90_PUT_VAR( ncvgid, nvarid_udif, (/puerr_max/), (/it/), (/1/) ) 
    1393             istatus = NF90_PUT_VAR( ncvgid, nvarid_vdif, (/pverr_max/), (/it/), (/1/) ) 
    1394             istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/it/), (/1/) ) 
     1654            istatus = NF90_PUT_VAR( ncvgid, nvarid_ucvg, (/zveldif/), (/it/), (/1/) )     ! max U or V velocity diff between subiterations 
     1655            istatus = NF90_PUT_VAR( ncvgid, nvarid_ures, (/zu_res_mean/), (/it/), (/1/) ) ! U-residual of the linear system 
     1656            istatus = NF90_PUT_VAR( ncvgid, nvarid_vres, (/zv_res_mean/), (/it/), (/1/) ) ! V-residual of the linear system 
     1657            istatus = NF90_PUT_VAR( ncvgid, nvarid_velres, (/zvelres/), (/it/), (/1/) )   ! average of u- and v- residuals 
     1658            istatus = NF90_PUT_VAR( ncvgid, nvarid_udif, (/puerr_max/), (/it/), (/1/) )   ! max U velocity difference, inner iterations 
     1659            istatus = NF90_PUT_VAR( ncvgid, nvarid_vdif, (/pverr_max/), (/it/), (/1/) )   ! max V velocity difference, inner iterations 
     1660            istatus = NF90_PUT_VAR( ncvgid, nvarid_mke, (/zmke/), (/it/), (/1/) )         ! mean kinetic energy 
     1661 
     1662            istatus = NF90_PUT_VAR( ncvgid, nvarid_usub, (/pu/), (/it/) )          ! u-velocity 
     1663            istatus = NF90_PUT_VAR( ncvgid, nvarid_vsub, (/pv/), (/it/) )          ! v-velocity 
     1664 
    13951665         ! close file 
    13961666         IF( kt == nitend )   istatus = NF90_CLOSE( ncvgid ) 
Note: See TracChangeset for help on using the changeset viewer.