Ignore:
Timestamp:
2008-04-29T14:40:00+02:00 (13 years ago)
Author:
ctlod
Message:

Manage wind stress at the sea-ice (LIM 3.0)/ocean interface in a proper way, see ticket: #128

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/LIM_SRC_3/limsbc.F90

    r888 r913  
    7777      INTEGER  ::   ifvt, i1mfr, idfr               ! some switches 
    7878      INTEGER  ::   iflt, ial, iadv, ifral, ifrdv 
    79       REAL(wp) ::   zqsr  , zqns     ! solar & non solar heat flux 
    8079      REAL(wp) ::   zinda            ! switch for testing the values of ice concentration 
    8180      REAL(wp) ::   zfons            ! salt exchanges at the ice/ocean interface 
    8281      REAL(wp) ::   zpme             ! freshwater exchanges at the ice/ocean interface 
    8382      REAL(wp) ::   zfrldu, zfrldv   ! lead fraction at U- & V-points 
    84       REAL(wp) ::   zutau , zvtau    ! lead fraction at U- & V-points 
    8583      REAL(wp) ::   zu_io , zv_io    ! 2 components of the ice-ocean velocity 
    86       REAL(wp) ::   ztairx, ztairy    
    87       REAL(wp) ::   zsfrldmx2, zsfrldmy2 
    88       REAL(wp) ::   zu_ice, zv_ice 
    89       REAL(wp) ::   ztglx , ztgly 
    90       REAL(wp) ::   ztaux , ztauy 
    9184       
    9285#if defined key_coupled     
     
    9487      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalbp    ! albedo of ice under clear sky 
    9588#endif 
    96       REAL(wp) ::   zsang, zmod, zfm 
     89      REAL(wp) ::   zsang, zmod 
    9790      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! ocean stress below sea-ice 
    9891      REAL(wp), DIMENSION(jpi,jpj) ::   zfcm1 , zfcm2    ! solar/non solar heat fluxes 
     
    277270 
    278271         ! ... ice stress over ocean with a ice-ocean rotation angle 
    279          DO jj = 2, jpjm1 
    280             zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg 
    281  
    282             DO ji = fs_2, fs_jpim1 
    283                ! computation of wind stress over ocean in X and Y direction 
    284                ztairx =  ( 2.0 - at_i(ji,jj) - at_i(ji+1,jj) ) * utau(ji,jj) 
    285                ztairy =  ( 2.0 - at_i(ji,jj) - at_i(ji,jj+1) ) * vtau(ji,jj) 
    286  
    287                zsfrldmx2 = at_i(ji,jj) + at_i(ji+1,jj) 
    288                zsfrldmy2 = at_i(ji,jj) + at_i(ji,jj+1) 
    289  
    290                zu_ice   = u_ice(ji,jj) - u_oce(ji,jj) 
    291                zv_ice   = v_ice(ji,jj) - v_oce(ji,jj) 
    292                zmod     = SQRT( zu_ice * zu_ice + zv_ice * zv_ice )  
    293  
     272         DO jj = 1, jpj 
     273            ! ... change the sinus angle sign in the south hemisphere 
     274            zsang  = SIGN(1.e0, gphif(1,jj) ) * sangvg 
     275            DO ji = 1, jpi 
     276               ! ... ice velocity relative to the ocean 
     277               zu_io = u_ice(ji,jj) - u_oce(ji,jj) 
     278               zv_io = v_ice(ji,jj) - v_oce(ji,jj) 
     279               zmod  = SQRT( zu_io * zu_io + zv_io * zv_io )  
    294280               ! quadratic drag formulation 
    295                ztglx   = zsfrldmx2 * rhoco * zmod * ( cangvg * zu_ice - zsang * zv_ice )  
    296                ztgly   = zsfrldmy2 * rhoco * zmod * ( cangvg * zv_ice + zsang * zu_ice )  
    297 ! 
    298 !              ! IMPORTANT 
    299 !              ! these lines are bound to prevent numerical oscillations 
    300 !              ! in the ice-ocean stress 
    301 !              ! They are physically ill-based. There is a cleaner solution 
    302 !              ! to try (remember discussion in Paris Gurvan) 
    303 ! 
    304                ztglx   = ztglx * exp( - zmod / 0.5 ) 
    305                ztgly   = ztglx * exp( - zmod / 0.5 ) 
    306  
    307                tio_u(ji,jj) = ( ztairx + 1.0 * ztglx ) / 2. 
    308                tio_v(ji,jj) = ( ztairy + 1.0 * ztgly ) / 2. 
     281               ztio_u(ji,jj) = rhoco * zmod * ( cangvg * zu_io - zsang * zv_io )  
     282               ztio_v(ji,jj) = rhoco * zmod * ( cangvg * zv_io + zsang * zu_io )  
     283               ! 
    309284            END DO 
    310285         END DO 
    311286 
    312          DO jj = 1, jpjm1 
    313             DO ji = 1, fs_jpim1   ! vertor opt. 
    314                zfrldu = MAX( freezn(ji,jj), freezn(ji,jj+1) )   ! ice/ocean indicator at U- and V-points 
    315                zfrldv = MAX( freezn(ji,jj), freezn(ji+1,jj) ) 
    316                ztaux = tio_u(ji,jj) 
    317                ztauy = tio_v(ji,jj) 
    318                utau(ji,jj) = (1.-zfrldu) * utau(ji,jj) + zfrldu * ztaux    ! stress at the ocean surface 
    319                vtau(ji,jj) = (1.-zfrldv) * vtau(ji,jj) + zfrldv * ztauy 
     287         DO jj = 2, jpjm1 
     288            DO ji = fs_2, fs_jpim1   ! vertor opt. 
     289               ! ... open-ocean (lead) fraction at U- & V-points (from T-point values) 
     290               zfrldu = 0.5 * ( ato_i(ji,jj) + ato_i(ji+1,jj  ) ) 
     291               zfrldv = 0.5 * ( ato_i(ji,jj) + ato_i(ji  ,jj+1) ) 
     292               ! update surface ocean stress 
     293               utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * ztio_u(ji,jj) 
     294               vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * ztio_v(ji,jj) 
     295               ! 
    320296            END DO 
    321297         END DO 
    322  
    323 !new fashion         DO jj = 2, jpjm1 
    324 !new fashion            DO ji = fs_2, fs_jpim1   ! vertor opt. 
    325 !new fashion               ! ... ice-cover wheighted ice-ocean stress at U and V-points  (from I-point values) 
    326 !new fashion               zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) ) 
    327 !new fashion               zvtau  = 0.5 * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) ) 
    328 !new fashion               ! ... open-ocean (lead) fraction at U- & V-points (from T-point values) 
    329 !new fashion               zfrldu = 0.5 * ( frld (ji,jj) + frld (ji+1,jj  ) ) 
    330 !new fashion               zfrldv = 0.5 * ( frld (ji,jj) + frld (ji  ,jj+1) ) 
    331 !new fashion               ! update surface ocean stress 
    332 !new fashion               utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * zutau 
    333 !new fashion               vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * zvtau 
    334 !new fashion               ! 
    335 !new fashion            END DO 
    336 !new fashion         END DO 
    337298 
    338299         ! boundary condition on the stress (utau,vtau) 
Note: See TracChangeset for help on using the changeset viewer.