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 15403 – NEMO

Changeset 15403


Ignore:
Timestamp:
2021-10-19T15:48:28+02:00 (10 months ago)
Author:
vancop
Message:

VP with much less bugs running at least 8 time steps on SAS

File:
1 edited

Legend:

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

    r13681 r15403  
    7979      !! 
    8080      !!  f1(u) = g1(v) 
    81       !!  f2(v) = g2(v) 
     81      !!  f2(v) = g2(u) 
    8282      !! 
    8383      !!  The right-hand side (RHS) is explicit 
     
    143143      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice 
    144144      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                       ! ice/snow mass and volume 
    145       REAL(wp) ::   zdeltat, zds2, zdt, zdt2, zdiv, zdiv2               ! temporary scalars 
     145      REAL(wp) ::   zds2, zdt, zdt2, zdiv, zdiv2                        ! temporary scalars 
    146146      REAL(wp) ::   zp_deltastar_f                                      !  
    147147      REAL(wp) ::   zu_cV, zv_cU                                        !  
     
    156156      REAL(wp), DIMENSION(jpi,jpj) ::   zmassU_t, zmassV_t              ! Mass per unit area divided by time step 
    157157      ! 
    158       REAL(wp), DIMENSION(jpi,jpj) ::   zdeltastar_t                    ! Delta* at T-points 
     158      REAL(wp), DIMENSION(jpi,jpj) ::   zdeltat, zdeltastar_t           ! Delta & Delta* at T-points 
    159159      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! Tension 
    160160      REAL(wp), DIMENSION(jpi,jpj) ::   zp_deltastar_t                  ! P/delta* at T points 
     
    207207      !!---------------------------------------------------------------------------------------------------------------------- 
    208208      ! DEBUG put all forcing terms to zero 
    209          ! air-ice drag 
    210          utau_ice(:,:) = 0._wp 
    211          vtau_ice(:,:) = 0._wp 
    212          ! coriolis 
    213          ff_t(:,:) = 0._wp 
    214          ! ice-ocean drag 
    215          rn_cio = 0._wp 
    216          ! ssh  
    217          ! done line 330 !!! dont forget to act there 
     209      ! air-ice drag 
     210      !utau_ice(:,:) = 0._wp 
     211      !vtau_ice(:,:) = 0._wp 
     212 
     213      ! coriolis 
     214      !ff_t(:,:) = 0._wp 
     215 
     216      ! ice-ocean drag 
     217      !rn_cio = 0._wp 
     218 
     219      ! sea surface height 
     220      !    embedded sea ice: compute representative ice top surface 
     221      !    non-embedded sea ice: use ocean surface for slope calculation 
     222      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) ! DEBUG 
     223      !zsshdyn(:,:) = 0._wp  
     224 
    218225      ! END DEBUG 
     226      !!---------------------------------------------------------------------------------------------------------------------- 
    219227 
    220228      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_vp: VP sea-ice rheology (LSR solver)' 
     
    265273      ENDIF 
    266274 
    267       zs1_rhsu(:,:) = 0._wp; zs2_rhsu(:,:) = 0._wp; zs1_rhsv(:,:) = 0._wp; zs2_rhsv(:,:) = 0._wp 
     275      zs1_rhsu(:,:) = 0._wp; zs2_rhsu(:,:) = 0._wp 
     276      zs1_rhsv(:,:) = 0._wp; zs2_rhsv(:,:) = 0._wp 
    268277      zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp; 
    269278      zAV(:,:) = 0._wp; zBV(:,:) = 0._wp; zCV(:,:) = 0._wp; zDV(:,:) = 0._wp; zEV(:,:) = 0._wp; 
    270       zrhsu(:,:) = 0._wp; zrhsv(:,:) = 0._wp 
    271       zf_rhsu(:,:) = 0._wp; zf_rhsv(:,:) = 0._wp 
     279      zrhsu  (:,:) = 0._wp ; zrhsv  (:,:) = 0._wp 
     280      zf_rhsu(:,:) = 0._wp ; zf_rhsv(:,:) = 0._wp 
    272281 
    273282      !------------------------------------------------------------------------------! 
     
    279288      CALL ice_strength ! strength at T points 
    280289       
    281       !------------------------------ 
    282       ! -- F-mask       (code from EVP) 
    283       !------------------------------ 
     290      !--------------------------- 
     291      ! -- F-mask (code from EVP) 
     292      !--------------------------- 
    284293      ! MartinV:  
    285294      ! In EVP routine, zfmask is applied on shear at F-points, in order to enforce the lateral boundary condition (no-slip, ..., free-slip) 
     
    289298      DO jj = 1, jpj - 1 
    290299         DO ji = 1, jpi - 1 
    291             zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 
     300            zfmask(ji,jj)    = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 
    292301         END DO 
    293302      END DO 
     
    326335      !---------------------------------------------------------------------------------------------------------- 
    327336      ! Compute all terms & factors independent of velocities, or only depending on velocities at previous time step       
    328        
    329       ! sea surface height 
    330       !    embedded sea ice: compute representative ice top surface 
    331       !    non-embedded sea ice: use ocean surface for slope calculation 
    332       zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b) 
    333       zsshdyn(:,:) = 0._wp ! DEBUG CAREFUL !!! 
    334337 
    335338      zmt(:,:) = rhos * vt_s(:,:) + rhoi * vt_i(:,:)       ! Snow and ice mass at T-point 
     
    347350            zm2             = zmt(ji+1,jj) 
    348351            zm3             = zmt(ji,jj+1) 
    349              
    350352            zmassU          = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 
    351353            zmassV          = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
     
    382384 
    383385! MV TEST DEBUG 
    384             IF ( ( zmt(ji,jj)   <= zmmin .OR. zmt(ji+1,jj)  <= zmmin )     .AND.  & 
    385                & ( at_i(ji,jj)  <= zamin .OR. at_i(ji+1,jj) <= zamin ) )    THEN   ;   zmsk01x(ji,jj) = 0._wp 
    386             ELSE                                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF 
    387  
    388             IF ( ( zmt(ji,jj)   <= zmmin .OR. zmt(ji,jj+1)  <= zmmin )     .AND.  & 
    389                & ( at_i(ji,jj)  <= zamin .OR. at_i(ji,jj+1) <= zamin ) )    THEN   ;   zmsk01y(ji,jj) = 0._wp 
    390             ELSE                                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF               
     386!           IF ( ( zmt(ji,jj)   <= zmmin .OR. zmt(ji+1,jj)  <= zmmin )     .AND.  & 
     387!              & ( at_i(ji,jj)  <= zamin .OR. at_i(ji+1,jj) <= zamin ) )    THEN   ;   zmsk01x(ji,jj) = 0._wp 
     388!           ELSE                                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF 
     389 
     390!           IF ( ( zmt(ji,jj)   <= zmmin .OR. zmt(ji,jj+1)  <= zmmin )     .AND.  & 
     391!              & ( at_i(ji,jj)  <= zamin .OR. at_i(ji,jj+1) <= zamin ) )    THEN   ;   zmsk01y(ji,jj) = 0._wp 
     392!           ELSE                                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF               
    391393! END MV TEST DEBUG 
    392394 
    393395         END DO 
    394396      END DO    
     397 
     398      CALL lbc_lnk_multi( 'icedyn_rhg_vp', zspgU   ,  'U', -1., zspgV   , 'V', -1. ) 
     399      CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_ai,  'U', -1., ztauy_ai, 'V', -1. ) 
    395400 
    396401      CALL iom_put( 'zmsk00x'    , zmsk00x  )   ! MV DEBUG 
     
    431436         ! In the outer loop, one needs to update all RHS terms 
    432437         ! with explicit velocity dependencies (viscosities, coriolis, ocean stress) 
    433          ! as a function of uc 
    434          ! 
     438         ! as a function of "current" velocities (uc, vc) 
    435439       
    436440         !------------------------------------------ 
     
    451455 
    452456         CALL lbc_lnk( 'icedyn_rhg_vp', zds, 'F', 1. ) ! MV TEST could be un-necessary according to Gurvan 
    453          CALL iom_put( 'zds'        , zds      )   ! MV DEBUG 
     457 
     458         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zds'          , zds          )   ! MV DEBUG 
    454459 
    455460         IF( lwp )   WRITE(numout,*) ' outer loop  1a i_out : ', i_out 
     
    464469 
    465470               ! shear**2 at T points (doc eq. A16) 
    466                zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  & 
    467                   &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  & 
    468                   &   ) * 0.25_wp * r1_e1e2t(ji,jj) 
     471               zds2  = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  & 
     472                  &    + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  & 
     473                  &    ) * 0.25_wp * r1_e1e2t(ji,jj) 
    469474               
    470475               ! divergence at T points 
     
    475480                
    476481               ! tension at T points 
    477                zdt  = ( ( zu_c(ji,jj) * r1_e2u(ji,jj) - zu_c(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    478                   &   - ( zv_c(ji,jj) * r1_e1v(ji,jj) - zv_c(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    479                   &   ) * r1_e1e2t(ji,jj) 
    480                zdt2 = zdt * zdt 
     482               zdt   = ( ( zu_c(ji,jj) * r1_e2u(ji,jj) - zu_c(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     483                  &    - ( zv_c(ji,jj) * r1_e1v(ji,jj) - zv_c(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     484                  &    ) * r1_e1e2t(ji,jj) 
     485               zdt2  = zdt * zdt 
    481486                
    482487               ! delta at T points 
    483                zdeltat = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
     488               zdeltat(ji,jj)        = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
    484489                
    485490               ! delta* at T points (following Lemieux and Dupont, GMD 2020) 
    486                zdeltastar_t(ji,jj) = zdeltat + rn_creepl 
     491               zdeltastar_t(ji,jj)   = zdeltat(ji,jj) + rn_creepl 
    487492               
    488493               ! P/delta at T-points 
     
    490495                
    491496               ! Temporary zzt and zet factors at T-points 
    492                zzt(ji,jj)     = zp_deltastar_t(ji,jj) * r1_e1e2t(ji,jj) 
    493                zet(ji,jj)     = zzt(ji,jj)     * z1_ecc2  
     497               zzt(ji,jj)            = zp_deltastar_t(ji,jj) * r1_e1e2t(ji,jj) 
     498               zet(ji,jj)            = zzt(ji,jj)     * z1_ecc2  
    494499                           
    495500            END DO 
    496501         END DO 
    497           
    498          CALL lbc_lnk_multi( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. ) 
    499  
    500          CALL iom_put( 'zzt'        , zzt      )   ! MV DEBUG 
    501          CALL iom_put( 'zet'        , zet      )   ! MV DEBUG 
    502          CALL iom_put( 'zp_deltastar_t', zp_deltastar_t ) ! MV DEBUG 
     502 
     503!        zzt(:,:) = 0. ! test 
     504!        zet(:,:) = 0. ! test 
     505          
     506         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zp_deltastar_t , 'T', 1. , zzt , 'T', 1., zet, 'T', 1. , zdeltat, 'T', 1.) 
     507 
     508         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zzt'        , zzt      )   ! MV DEBUG 
     509         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zet'        , zet      )   ! MV DEBUG 
     510         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zp_deltastar_t', zp_deltastar_t ) ! MV DEBUG 
     511         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zdeltat'       , zdeltat        ) ! MV DEBUG 
    503512 
    504513         IF( lwp )   WRITE(numout,*) ' outer loop  1b i_out : ', i_out 
     
    515524            END DO 
    516525         END DO 
     526 
     527!        zef(:,:) = 0. ! test 
    517528          
    518529         CALL lbc_lnk( 'icedyn_rhg_vp', zef, 'F', 1. ) 
    519          CALL iom_put( 'zef'           , zef            ) ! MV DEBUG 
     530         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zef'           , zef            ) ! MV DEBUG 
    520531         IF( lwp )   WRITE(numout,*) ' outer loop  1c i_out : ', i_out 
    521532 
     
    551562                           &             ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * zu_c(ji,jj  ) + e2u(ji-1,jj  ) * zu_c(ji-1,jj  ) )  & 
    552563                           &             + zmf(ji,jj+1) * ( e2u(ji,jj+1) * zu_c(ji,jj+1) + e2u(ji-1,jj+1) * zu_c(ji-1,jj+1) ) ) 
    553           
     564 
    554565             END DO 
    555566         END DO 
     
    560571         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCorU,  'U', -1., zCorV, 'V', -1. ) 
    561572 
    562          CALL iom_put( 'zCwU'          , zCwU           ) ! MV DEBUG 
    563          CALL iom_put( 'zCwV'          , zCwV           ) ! MV DEBUG 
    564          CALL iom_put( 'zCorU'         , zCorU          ) ! MV DEBUG 
    565          CALL iom_put( 'zCorV'         , zCorV          ) ! MV DEBUG 
     573         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zCwU'          , zCwU           ) ! MV DEBUG 
     574         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zCwV'          , zCwV           ) ! MV DEBUG 
     575         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zCorU'         , zCorU          ) ! MV DEBUG 
     576         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zCorV'         , zCorV          ) ! MV DEBUG 
    566577 
    567578         IF( lwp )   WRITE(numout,*) ' outer loop  1f i_out : ', i_out 
     
    575586 
    576587         ! --- Stress contributions at T-points          
    577          DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12 
     588         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1 & zs2 
    578589            DO ji = 2, jpi !  
    579590             
    580591               ! sig1 contribution to RHS of U-equation at T-points 
    581                zs1_rhsu(ji,jj) =   zzt(ji,jj) * ( e1v(ji,jj)    * zv_c(ji,jj) - e1v(ji,jj-1)    * zv_c(ji,jj-1) - 1.0_wp ) 
     592               zs1_rhsu(ji,jj) =   zzt(ji,jj) * ( e1v(ji,jj)    * zv_c(ji,jj) - e1v(ji,jj-1)    * zv_c(ji,jj-1) )   & 
     593                               &                - zp_deltastar_t(ji,jj) * zdeltat(ji,jj) 
    582594                                             
    583595               ! sig2 contribution to RHS of U-equation at T-points             
     
    585597 
    586598               ! sig1 contribution to RHS of V-equation at T-points 
    587                zs1_rhsv(ji,jj) =   zzt(ji,jj) * ( e2u(ji,jj)    * zu_c(ji,jj) - e2u(ji-1,jj)    * zu_c(ji-1,jj) - 1.0_wp ) 
     599               zs1_rhsv(ji,jj) =   zzt(ji,jj) * ( e2u(ji,jj)    * zu_c(ji,jj) - e2u(ji-1,jj)    * zu_c(ji-1,jj) )   &  
     600                               &                - zp_deltastar_t(ji,jj) * zdeltat(ji,jj) 
    588601 
    589602               ! sig2 contribution to RHS of V-equation  at T-points 
     
    593606         END DO 
    594607 
    595          CALL iom_put( 'zs1_rhsu'      , zs1_rhsu       ) ! MV DEBUG 
    596          CALL iom_put( 'zs2_rhsu'      , zs2_rhsu       ) ! MV DEBUG 
    597          CALL iom_put( 'zs1_rhsv'      , zs1_rhsv       ) ! MV DEBUG 
    598          CALL iom_put( 'zs2_rhsv'      , zs2_rhsv       ) ! MV DEBUG 
     608         ! TEST 
     609         ! zs2_rhsv(:,:) = 0._wp 
     610 
     611         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs1_rhsu'      , zs1_rhsu       ) ! MV DEBUG 
     612         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs2_rhsu'      , zs2_rhsu       ) ! MV DEBUG 
     613         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs1_rhsv'      , zs1_rhsv       ) ! MV DEBUG 
     614         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs2_rhsv'      , zs2_rhsv       ) ! MV DEBUG 
    599615          
    600616         ! a priori, no lbclnk, because rhsu is only used in the inner domain 
    601617          
    602          ! --- Stress contributions at f-points          
     618         ! --- Stress contributions at F-points          
    603619         ! MV NOTE: I applied zfmask on zds, by mimetism on EVP, but without deep understanding of what I was doing 
    604620         ! My guess is that this is the way to enforce boundary conditions on strain rate tensor 
     
    610626                
    611627               ! sig12 contribution to RHS of U equation at F-points  
    612                zs12_rhsu(ji,jj) = - zef(ji,jj)  * ( r1_e2v(ji+1,jj) * zv_c(ji+1,jj) - r1_e2v(ji,jj) * zv_c(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) * zfmask(ji,jj) 
     628               zs12_rhsu(ji,jj) =   zef(ji,jj)  * ( r1_e2v(ji+1,jj) * zv_c(ji+1,jj) + r1_e2v(ji,jj) * zv_c(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) * zfmask(ji,jj) 
    613629                
    614630               ! sig12 contribution to RHS of V equation at F-points 
    615                zs12_rhsv(ji,jj) =   zef(ji,jj)  * ( r1_e1u(ji,jj+1) * zu_c(ji,jj+1) - r1_e1u(ji,jj) * zu_c(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) * zfmask(ji,jj) 
     631               zs12_rhsv(ji,jj) =   zef(ji,jj)  * ( r1_e1u(ji,jj+1) * zu_c(ji,jj+1) + r1_e1u(ji,jj) * zu_c(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj) * zfmask(ji,jj) 
    616632 
    617633            END DO 
    618634         END DO 
     635         ! 
     636         ! TEST 
     637         ! zs12_rhsv(:,:) = 0._wp 
    619638 
    620639         CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsu, 'F', 1. ) 
    621640         CALL lbc_lnk( 'icedyn_rhg_vp', zs12_rhsv, 'F', 1. ) 
    622641 
    623          CALL iom_put( 'zs12_rhsu'     , zs12_rhsu      ) ! MV DEBUG 
    624          CALL iom_put( 'zs12_rhsv'     , zs12_rhsv      ) ! MV DEBUG 
     642         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs12_rhsu'     , zs12_rhsu      ) ! MV DEBUG 
     643         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zs12_rhsv'     , zs12_rhsv      ) ! MV DEBUG 
    625644 
    626645         ! a priori, no lbclnk, because rhsu are only used in the inner domain 
     
    631650         DO jj = 2, jpj - 1 
    632651            DO ji = 2, jpi - 1                
     652 
    633653               ! --- U component of internal force contribution to RHS at U points 
    634654               zf_rhsu(ji,jj) = 0.5_wp * r1_e1e2u(ji,jj) * &  
     
    640660               zf_rhsv(ji,jj) = 0.5_wp * r1_e1e2v(ji,jj) * & 
    641661                  &           (    e1v(ji,jj)    * ( zs1_rhsv(ji,jj+1) - zs1_rhsv(ji,jj) )                                                                 & 
    642                   &      +         r1_e1v(ji,jj) * ( e1t(ji,jj+1) * e1t(ji,jj+1) * zs2_rhsv(ji,jj+1) - e1t(ji,jj) * e1t(ji,jj)     * zs2_rhsv(ji,jj) )     & 
     662                  &      -         r1_e1v(ji,jj) * ( e1t(ji,jj+1) * e1t(ji,jj+1) * zs2_rhsv(ji,jj+1) - e1t(ji,jj)   * e1t(ji,jj)   * zs2_rhsv(ji,jj) )     & 
    643663                  &      + 2._wp * r1_e2v(ji,jj) * ( e2f(ji,jj)   * e2f(ji,jj)   * zs12_rhsv(ji,jj)  - e2f(ji-1,jj) * e2f(ji-1,jj) * zs12_rhsv(ji-1,jj) ) ) 
    644                    
     664 
    645665            END DO 
    646666         END DO 
    647667 
    648          CALL iom_put( 'zf_rhsu'       , zf_rhsu        ) ! MV DEBUG 
    649          CALL iom_put( 'zf_rhsv'       , zf_rhsv        ) ! MV DEBUG 
     668         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zf_rhsu'       , zf_rhsu        ) ! MV DEBUG 
     669         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zf_rhsv'       , zf_rhsv        ) ! MV DEBUG 
    650670          
    651671         !--------------------------- 
     
    657677            DO ji = 2, jpi - 1 
    658678             
    659                ! still miss ice ocean stress and acceleration contribution 
    660679               zrhsu(ji,jj) = zmU_t(ji,jj) + ztaux_ai(ji,jj) + ztaux_oi_rhsu(ji,jj) + zspgU(ji,jj) + zCorU(ji,jj) + zf_rhsu(ji,jj) 
    661                zrhsv(ji,jj) = zmV_t(ji,jj) + ztauy_ai(ji,jj) + ztauy_oi_rhsv(ji,jj) + zspgV(ji,jj) + zCorV(ji,jj) + zf_rhsu(ji,jj) 
     680               zrhsv(ji,jj) = zmV_t(ji,jj) + ztauy_ai(ji,jj) + ztauy_oi_rhsv(ji,jj) + zspgV(ji,jj) + zCorV(ji,jj) + zf_rhsv(ji,jj) 
    662681 
    663682            END DO 
     
    668687         CALL lbc_lnk_multi( 'icedyn_rhg_vp', ztaux_oi_rhsu, 'U', -1., ztauy_oi_rhsv, 'V',  -1.) 
    669688 
    670          CALL iom_put( 'zmU_t'         , zmU_t          ) ! MV DEBUG 
    671          CALL iom_put( 'zmV_t'         , zmV_t          ) ! MV DEBUG 
    672          CALL iom_put( 'ztaux_oi_rhsu' , ztaux_oi_rhsu  ) ! MV DEBUG 
    673          CALL iom_put( 'ztauy_oi_rhsv' , ztauy_oi_rhsv  ) ! MV DEBUG 
    674          CALL iom_put( 'zrhsu'         , zrhsu          ) ! MV DEBUG 
    675          CALL iom_put( 'zrhsv'         , zrhsv          ) ! MV DEBUG 
     689         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zmU_t'         , zmU_t          ) ! MV DEBUG 
     690         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zmV_t'         , zmV_t          ) ! MV DEBUG 
     691         IF ( i_out == nn_nout_vp ) CALL iom_put( 'ztaux_oi_rhsu' , ztaux_oi_rhsu  ) ! MV DEBUG 
     692         IF ( i_out == nn_nout_vp ) CALL iom_put( 'ztauy_oi_rhsv' , ztauy_oi_rhsv  ) ! MV DEBUG 
     693         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zrhsu'         , zrhsu          ) ! MV DEBUG 
     694         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zrhsv'         , zrhsv          ) ! MV DEBUG 
    676695          
    677696         ! inner domain calculations -> no lbclnk 
     
    712731               zfac3      = 2._wp * zfac * r1_e1u(ji,jj) 
    713732 
    714                zt12U      = - zfac1 * zzt(ji+1,jj) 
    715733               zt11U      =   zfac1 * zzt(ji,jj) 
    716           
    717                zt22U      = - zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) 
     734               zt12U      =   zfac1 * zzt(ji+1,jj) 
     735          
    718736               zt21U      =   zfac2 * zet(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj)   * e2t(ji,jj) 
    719           
    720                zt122U     = - zfac3 * zef(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj) 
     737               zt22U      =   zfac2 * zet(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) 
     738          
    721739               zt121U     =   zfac3 * zef(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) 
     740               zt122U     =   zfac3 * zef(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj)   * e1f(ji,jj) 
    722741                
    723742               ! 
    724743               ! Linear system coefficients 
    725744               ! 
    726                zAU(ji,jj) = - zt11U * e2u(ji-1,jj) - zt21U * r1_e2u(ji-1,jj) 
    727                zBU(ji,jj) = ( zt12U + zt11U ) * e2u(ji,jj) + ( zt22U + zt21U ) * r1_e2u(ji,jj) + ( zt122U + zt121U ) * r1_e1u(ji,jj) 
    728                zCU(ji,jj) = - zt12U * e2u(ji+1,jj) - zt22U * r1_e2u(ji+1,jj) 
    729           
    730                zDU(ji,jj) =   zt121U * r1_e1u(ji,jj-1) 
    731                zEU(ji,jj) =   zt122U * r1_e1u(ji,jj+1) 
     745               zAU(ji,jj) = -   zt11U           * e2u(ji-1,jj) -   zt21U          * r1_e2u(ji-1,jj) 
     746               zBU(ji,jj) =   ( zt11U + zt12U ) * e2u(ji,jj)   + ( zt21U + zt22U ) * r1_e2u(ji,jj)   + ( zt121U + zt122U ) * r1_e1u(ji,jj) 
     747               zCU(ji,jj) = -   zt12U           * e2u(ji+1,jj) -   zt22U          * r1_e2u(ji+1,jj) 
     748          
     749               zDU(ji,jj) =     zt121U * r1_e1u(ji,jj-1) 
     750               zEU(ji,jj) =     zt122U * r1_e1u(ji,jj+1) 
    732751               
    733752               ! 
     
    737756               ! 
    738757               zfac       = 0.5_wp * r1_e1e2v(ji,jj) 
    739                zfac1      =         zfac * e2v(ji,jj) 
     758               zfac1      =         zfac * e1v(ji,jj) 
    740759               zfac2      =         zfac * r1_e1v(ji,jj) 
    741760               zfac3      = 2._wp * zfac * r1_e2v(ji,jj) 
    742761          
    743                zt12V      = - zfac1 * zzt(ji,jj+1) 
    744762               zt11V      =   zfac1 * zzt(ji,jj) 
    745           
     763               zt12V      =   zfac1 * zzt(ji,jj+1) 
     764          
     765               zt21V      =   zfac2 * zet(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj) 
    746766               zt22V      =   zfac2 * zet(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) 
    747                zt21V      = - zfac2 * zet(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj)   * e1t(ji,jj) 
    748           
     767          
     768               zt121V     =   zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) 
    749769               zt122V     =   zfac3 * zef(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj)   * e2f(ji,jj) 
    750                zt121V     = - zfac3 * zef(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) 
    751           
     770 
    752771               ! 
    753772               ! Linear system coefficients 
    754773               ! 
    755                zAV(ji,jj) = - zt11V * e1v(ji,jj-1) + zt21V * r1_e1v(ji,jj-1) 
    756                zBV(ji,jj) =  ( zt12V + zt11V ) * e1v(ji,jj) - ( zt22V + zt21V ) * r1_e1v(ji,jj) - ( zt122V + zt121V ) * r1_e2v(ji,jj) 
    757                zCV(ji,jj) = - zt12V * e1v(ji,jj+1) + zt22V * r1_e1v(ji,jj+1) 
    758           
    759                zDV(ji,jj) = - zt121V * r1_e2v(ji-1,jj) ! mistake is in the pdf notes not here 
    760                zEV(ji,jj) = - zt122V * r1_e2v(ji+1,jj) 
     774               zAV(ji,jj) = -   zt11V           * e1v(ji,jj-1) -   zt21V          * r1_e1v(ji,jj-1) 
     775               zBV(ji,jj) =   ( zt11V + zt12V ) * e1v(ji,jj)   + ( zt21V + zt22V ) * r1_e1v(ji,jj)   + ( zt122V + zt121V ) * r1_e2v(ji,jj) 
     776               zCV(ji,jj) = -   zt12V           * e1v(ji,jj+1) -   zt22V          * r1_e1v(ji,jj+1) 
     777          
     778               zDV(ji,jj) =     zt121V * r1_e2v(ji-1,jj) 
     779               zEV(ji,jj) =     zt122V * r1_e2v(ji+1,jj) 
    761780                   
    762781               !----------------------------------------------------- 
     
    764783               !----------------------------------------------------- 
    765784               zBU(ji,jj) = zBU(ji,jj) + zCwU(ji,jj) + zmassU_t(ji,jj) 
    766                zBV(ji,jj) = ZBV(ji,jj) + zCwV(ji,jj) + zmassV_t(ji,jj) 
     785               zBV(ji,jj) = zBV(ji,jj) + zCwV(ji,jj) + zmassV_t(ji,jj) 
    767786          
    768787            END DO 
    769788         END DO 
    770789 
    771          CALL lbc_lnk_multi( 'icedyn_rhg_vp', zAU  , 'U', 1., zAV  , 'V',  1. ) 
     790         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zAU  , 'U', 1., zAV  , 'V',  1. ) ! --> here normal that we use '1' and not '-1' ? 
    772791         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zBU  , 'U', 1., zBV  , 'V',  1. ) 
    773792         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zCU  , 'U', 1., zCV  , 'V',  1. ) 
     
    775794         CALL lbc_lnk_multi( 'icedyn_rhg_vp', zEU  , 'U', 1., zEV  , 'V',  1. ) 
    776795                
    777          CALL iom_put( 'zAU'           , zAU            ) ! MV DEBUG 
    778          CALL iom_put( 'zBU'           , zBU            ) ! MV DEBUG 
    779          CALL iom_put( 'zCU'           , zCU            ) ! MV DEBUG 
    780          CALL iom_put( 'zDU'           , zDU            ) ! MV DEBUG 
    781          CALL iom_put( 'zEU'           , zEU            ) ! MV DEBUG 
    782          CALL iom_put( 'zAV'           , zAV            ) ! MV DEBUG 
    783          CALL iom_put( 'zBV'           , zBV            ) ! MV DEBUG 
    784          CALL iom_put( 'zCV'           , zCV            ) ! MV DEBUG 
    785          CALL iom_put( 'zDV'           , zDV            ) ! MV DEBUG 
    786          CALL iom_put( 'zEV'           , zEV            ) ! MV DEBUG 
     796         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zAU'           , zAU            ) ! MV DEBUG 
     797         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zBU'           , zBU            ) ! MV DEBUG 
     798         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zCU'           , zCU            ) ! MV DEBUG 
     799         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zDU'           , zDU            ) ! MV DEBUG 
     800         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zEU'           , zEU            ) ! MV DEBUG 
     801 
     802         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zAV'           , zAV            ) ! MV DEBUG 
     803         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zBV'           , zBV            ) ! MV DEBUG 
     804         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zCV'           , zCV            ) ! MV DEBUG 
     805         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zDV'           , zDV            ) ! MV DEBUG 
     806         IF ( i_out == nn_nout_vp ) CALL iom_put( 'zEV'           , zEV            ) ! MV DEBUG 
    787807 
    788808      !------------------------------------------------------------------------------! 
     
    808828            zv_b(:,:)      = v_ice(:,:) 
    809829 
    810 !           zAU(:,:) = 0._wp; zBU(:,:) = 0._wp; zCU(:,:) = 0._wp; zDU(:,:) = 0._wp; zEU(:,:) = 0._wp 
    811 !           zAV(:,:) = 0._wp; zBV(:,:) = 0._wp; zCV(:,:) = 0._wp; zDV(:,:) = 0._wp; zEV(:,:) = 0._wp 
     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 
    812832 
    813833            IF ( ll_u_iterate .OR. ll_v_iterate )   THEN 
     
    864884                     DO jj = jj_min, jpj - 1, nn_zebra_vp 
    865885       
     886                        ! MV BUG OCT 6. Both zBU_prime and zFU_prime were undefined for ji = 2 
     887                        zBU_prime(2,jj)     = zBU(2,jj) 
     888                        zFU_prime(2,jj)     = zFU(2,jj) 
     889 
    866890                        DO ji = 3, jpi - 1 
    867891 
     
    962986                     DO ji = ji_min, jpi - 1, nn_zebra_vp  
    963987                      
    964                         DO jj = 3, jpj - 1 
     988                        ! MV BUG OCT 6.. Both zBV_prime and zFV_prime were undefined for jj = 2 
     989                        zBV_prime(ji,2)     = zBV(ji,2) 
     990                        zFV_prime(ji,2)     = zFV(ji,2) 
     991 
     992                        DO jj = 3, jpj - 1  
    965993 
    966994                           zfac             = SIGN( 1._wp , zBV(ji,jj-1) ) * MAX( 0._wp , SIGN( 1._wp , ABS( zBV(ji,jj-1) ) - epsi20 ) ) 
     
    11371165 
    11381166      ! MV DEBUG test - replace ice velocity by ocean current to give the model the means to go ahead 
    1139       DO jj = 2, jpj - 1 
    1140          DO ji = 2, jpi - 1    
    1141  
    1142              u_ice(ji,jj) =   zmsk00x(ji,jj)                               &  
    1143       &         * (           zmsk01x(ji,jj)   * u_oce(ji,jj) * 0.01_wp    & 
    1144                   + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp    ) 
    1145  
    1146              v_ice(ji,jj) =   zmsk00y(ji,jj)                               &  
    1147       &         * (           zmsk01y(ji,jj)   * v_oce(ji,jj) * 0.01_wp    & 
    1148                   + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp    ) 
    1149  
    1150          END DO 
    1151       END DO 
    1152  
    1153       CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 
    1154  
    1155       IF ( lwp ) WRITE(numout,*) ' Velocity replaced ' 
     1167      ! info: probably wrong because when doing du/dt = div ( sigma ) we have spurious waves instead of u = v = 0 
     1168 
     1169!     DO jj = 2, jpj - 1 
     1170!        DO ji = 2, jpi - 1    
     1171 
     1172!            u_ice(ji,jj) =   zmsk00x(ji,jj)                               &  
     1173!     &         * (           zmsk01x(ji,jj)   * u_oce(ji,jj) * 0.01_wp    & 
     1174!                 + ( 1._wp - zmsk01x(ji,jj) ) * u_oce(ji,jj) * 0.01_wp    ) 
     1175 
     1176!            v_ice(ji,jj) =   zmsk00y(ji,jj)                               &  
     1177!     &         * (           zmsk01y(ji,jj)   * v_oce(ji,jj) * 0.01_wp    & 
     1178!                 + ( 1._wp - zmsk01y(ji,jj) ) * v_oce(ji,jj) * 0.01_wp    ) 
     1179 
     1180!        END DO 
     1181!     END DO 
     1182 
     1183!     CALL lbc_lnk_multi( 'icedyn_rhg_vp', u_ice, 'U', -1., v_ice, 'V', -1. ) 
     1184 
     1185!     IF ( lwp ) WRITE(numout,*) ' Velocity replaced ' 
    11561186 
    11571187      ! END DEBUG 
Note: See TracChangeset for help on using the changeset viewer.