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 4333 for branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90 – NEMO

Ignore:
Timestamp:
2013-12-11T18:34:00+01:00 (10 years ago)
Author:
clem
Message:

remove remaining bugs in LIM3, so that it can run in both regional and global config

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r4269 r4333  
    5050   PUBLIC   lim_rhg        ! routine called by lim_dyn (or lim_dyn_2) 
    5151 
     52   REAL(wp) ::   epsi10 = 1.e-10_wp   ! 
    5253   REAL(wp) ::   rzero   = 0._wp   ! constant values 
    5354   REAL(wp) ::   rone    = 1._wp   ! constant values 
     
    426427               !-Calculate stress tensor components zs1 and zs2  
    427428               !-at centre of grid cells (see section 3.5 of CICE user's guide). 
    428                !zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) +   & 
    429                !   &          ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) )  &        
    430                !   &          / ( 1._wp + alphaevp * dtotel ) 
    431  
    432                !zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) -   & 
    433                !              zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )   & 
    434                !   &          / ( 1._wp + alphaevp * ecc2 * dtotel ) 
     429               zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) +   & 
     430                  &          ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) )  &        
     431                  &          / ( 1._wp + alphaevp * dtotel ) 
     432 
     433               zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) -   & 
     434                             zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )   & 
     435                  &          / ( 1._wp + alphaevp * ecc2 * dtotel ) 
    435436 
    436437               ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp) 
    437                zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) )  & 
    438                   &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel ) 
    439                zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )  & 
    440                   &         / ( 1._wp + dtotel ) 
     438               !zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) )  & 
     439               !   &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel ) 
     440               !zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )  & 
     441               !   &         / ( 1._wp + dtotel ) 
    441442 
    442443            END DO 
     
    472473 
    473474               !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide). 
    474                !zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) /  & 
    475                !   &          ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  & 
    476                !   &          / ( 1._wp + alphaevp * ecc2 * dtotel )  
     475               zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) /  & 
     476                  &          ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  & 
     477                  &          / ( 1._wp + alphaevp * ecc2 * dtotel )  
    477478 
    478479               ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp) 
    479                zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  & 
    480                   &          ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  & 
    481                   &          / ( 1.0 + dtotel )  
     480               !zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  & 
     481               !   &          ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  & 
     482               !   &          / ( 1.0 + dtotel )  
    482483 
    483484            END DO ! ji 
     
    485486 
    486487         CALL lbc_lnk( zs12(:,:), 'F', 1. ) 
    487  
    488 !#if defined key_bdy 
    489 !         ! clem: change zs1, zs2, zs12 at the boundary for each iteration 
    490 !         CALL bdy_ice_lim_dyn( 2, zs1, zs2, zs12 ) 
    491 !         CALL lbc_lnk( zs1 (:,:), 'T', 1. ) 
    492 !         CALL lbc_lnk( zs2 (:,:), 'T', 1. ) 
    493 !         CALL lbc_lnk( zs12(:,:), 'F', 1. ) 
    494 !#endif          
    495488 
    496489         ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 
     
    544537            CALL agrif_rhg_lim2( jter, nevp, 'U' ) 
    545538#endif 
     539#if defined key_bdy 
     540         ! clem: change u_ice and v_ice at the boundary for each iteration 
     541         CALL bdy_ice_lim_dyn( 'U' ) 
     542#endif          
    546543 
    547544!CDIR NOVERRCHK 
     
    572569            CALL agrif_rhg_lim2( jter, nevp, 'V' ) 
    573570#endif 
     571#if defined key_bdy 
     572         ! clem: change u_ice and v_ice at the boundary for each iteration 
     573         CALL bdy_ice_lim_dyn( 'V' ) 
     574#endif          
    574575 
    575576         ELSE  
     
    599600            CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
    600601#if defined key_agrif && defined key_lim2 
    601             CALL agrif_rhg_lim2( jter, nevp , 'V' ) 
    602 #endif 
     602            CALL agrif_rhg_lim2( jter, nevp, 'V' ) 
     603#endif 
     604#if defined key_bdy 
     605         ! clem: change u_ice and v_ice at the boundary for each iteration 
     606         CALL bdy_ice_lim_dyn( 'V' ) 
     607#endif          
    603608 
    604609!CDIR NOVERRCHK 
     
    632637            CALL agrif_rhg_lim2( jter, nevp, 'U' ) 
    633638#endif 
     639#if defined key_bdy 
     640         ! clem: change u_ice and v_ice at the boundary for each iteration 
     641         CALL bdy_ice_lim_dyn( 'U' ) 
     642#endif          
    634643 
    635644         ENDIF 
    636645          
    637 !#if defined key_bdy 
    638 !         ! clem: change u_ice and v_ice at the boundary for each iteration 
    639 !         CALL bdy_ice_lim_dyn( 1 ) 
    640 !#endif          
    641  
    642646         IF(ln_ctl) THEN 
    643647            !---  Convergence test. 
     
    666670!CDIR NOVERRCHK 
    667671         DO ji = fs_2, fs_jpim1 
    668             zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) )  
    669             !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 
     672            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) )  
     673            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 ) 
    670674            zdummy = vt_i(ji,jj) 
    671675            IF ( zdummy .LE. hminrhg ) THEN 
     
    684688#if defined key_bdy 
    685689      ! clem: change u_ice and v_ice at the boundary 
    686       CALL bdy_ice_lim_dyn( 1 ) 
     690      CALL bdy_ice_lim_dyn( 'U' ) 
     691      CALL bdy_ice_lim_dyn( 'V' ) 
    687692#endif          
    688693 
    689694      DO jj = k_j1+1, k_jpj-1  
    690695         DO ji = fs_2, fs_jpim1 
    691             zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) )  
    692             !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 
     696            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) )  
     697            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 ) 
    693698            zdummy = vt_i(ji,jj) 
    694699            IF ( zdummy .LE. hminrhg ) THEN 
     
    714719            !- zdd(:,:), zdt(:,:): divergence and tension at centre  
    715720            !- zds(:,:): shear on northeast corner of grid cells 
    716             zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) )  
    717             !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 
     721            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - epsi10 ) )  
     722            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , epsi10 ) 
    718723            zdummy = vt_i(ji,jj) 
    719724            IF ( zdummy .LE. hminrhg ) THEN 
Note: See TracChangeset for help on using the changeset viewer.