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 2370 for branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90 – NEMO

Ignore:
Timestamp:
2010-11-10T08:48:54+01:00 (13 years ago)
Author:
gm
Message:

v3.3beta: ice-ocean stress at kt with VP & EVP (LIM-2 and -3)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r2287 r2370  
    44   !!           computation of the flux at the sea ice/ocean interface 
    55   !!====================================================================== 
    6    !! History :   -   !  2006-07 (M. Vancoppelle)  LIM3 original code 
    7    !!            3.0  !  2008-03 (C. Tallandier)  surface module 
    8    !!             -   !  2008-04 (C. Tallandier)  split in 2 + new ice-ocean coupling 
     6   !! History :   -   ! 2006-07 (M. Vancoppelle)  LIM3 original code 
     7   !!            3.0  ! 2008-03 (C. Tallandier)  surface module 
     8   !!             -   ! 2008-04 (C. Tallandier)  split in 2 + new ice-ocean coupling 
     9   !!            3.3  ! 2010-05 (G. Madec) decrease ocean & ice reference salinities in the Baltic sea 
     10   !!                 !                  + simplification of the ice-ocean stress calculation 
    911   !!---------------------------------------------------------------------- 
    1012#if defined key_lim3 
     
    1214   !!   'key_lim3'                                    LIM 3.0 sea-ice model 
    1315   !!---------------------------------------------------------------------- 
    14    !!   lim_sbc  : flux at the ice / ocean interface 
    15    !!---------------------------------------------------------------------- 
    16    USE oce 
     16   !!   lim_sbc_flx  : updates mass, heat and salt fluxes at the ocean surface 
     17   !!   lim_sbc_tau  : update i- and j-stresses, and its modulus at the ocean surface 
     18   !!---------------------------------------------------------------------- 
    1719   USE par_oce          ! ocean parameters 
    1820   USE par_ice          ! ice parameters 
     
    3537   PUBLIC   lim_sbc_tau   ! called by sbc_ice_lim 
    3638 
    37    REAL(wp)  ::   epsi16 = 1.e-16  ! constant values 
    38    REAL(wp)  ::   rzero  = 0.e0     
    39    REAL(wp)  ::   rone   = 1.e0 
    40  
    41    REAL(wp), DIMENSION(jpi,jpj) ::   utau_oce, vtau_oce   !: air-ocean surface i- & j-stress              [N/m2] 
    42    REAL(wp), DIMENSION(jpi,jpj) ::   tmod_io              !: modulus of the ice-ocean relative velocity   [m/s] 
    43    REAL(wp), DIMENSION(jpi,jpj) ::   ssu_mb  , ssv_mb     !: before mean ocean surface currents           [m/s] 
     39   REAL(wp)  ::   r1_rdtice            ! = 1. / rdt_ice  
     40   REAL(wp)  ::   epsi16 = 1.e-16_wp   ! constant values 
     41   REAL(wp)  ::   rzero  = 0._wp     
     42   REAL(wp)  ::   rone   = 1._wp 
     43 
     44   REAL(wp), DIMENSION(jpi,jpj) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress              [N/m2] 
     45   REAL(wp), DIMENSION(jpi,jpj) ::   tmod_io              ! modulus of the ice-ocean relative velocity   [m/s] 
     46 
     47   REAL(wp), DIMENSION(jpi,jpj) ::   soce_0, sice_0   ! constant SSS and ice salinity used in levitating sea-ice case 
    4448 
    4549   !! * Substitutions 
     
    4852   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 
    4953   !! $Id$ 
    50    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 
    51    !!---------------------------------------------------------------------- 
    52  
     54   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     55   !!---------------------------------------------------------------------- 
    5356CONTAINS 
    54  
    55    SUBROUTINE lim_sbc_tau( kt, kcpl ) 
    56       !!------------------------------------------------------------------- 
    57       !!                ***  ROUTINE lim_sbc_tau *** 
    58       !!   
    59       !! ** Purpose : Update the ocean surface stresses due to the ice 
    60       !!          
    61       !! ** Action  : - compute the ice-ocean stress depending on kcpl: 
    62       !!              fluxes at the ice-ocean interface. 
    63       !!     Case 0  :  old LIM-3 way, computed at ice time-step only 
    64       !!     Case 1  :  at each ocean time step the stresses are computed 
    65       !!                using the current ocean velocity (now) 
    66       !!     Case 2  :  combination of half case 0 + half case 1 
    67       !!      
    68       !! ** Outputs : - utau   : sea surface i-stress (ocean referential) 
    69       !!              - vtau   : sea surface j-stress (ocean referential) 
    70       !! 
    71       !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90. 
    72       !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108. 
    73       !!--------------------------------------------------------------------- 
    74       INTEGER ::   kt     ! number of ocean iteration 
    75       INTEGER ::   kcpl   ! ice-ocean coupling indicator: =0  LIM-3 old case 
    76       !                   !                               =1  stresses computed using now ocean velocity 
    77       !                   !                               =2  combination of 0 and 1 cases 
    78       !! 
    79       INTEGER  ::   ji, jj   ! dummy loop indices 
    80       REAL(wp) ::   zfrldu, zat_u, zu_ico, zutaui, zu_u, zu_ij, zu_im1j   ! temporary scalar 
    81       REAL(wp) ::   zfrldv, zat_v, zv_ico, zvtaui, zv_v, zv_ij, zv_ijm1   !    -         - 
    82       REAL(wp) ::   zsang, zztmp                                          !    -         - 
    83       REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! ocean stress below sea-ice 
    84 #if defined key_coupled     
    85       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb     ! albedo of ice under overcast sky 
    86       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalbp    ! albedo of ice under clear sky 
    87 #endif 
    88      !!--------------------------------------------------------------------- 
    89  
    90       IF( kt == nit000 ) THEN 
    91          IF(lwp) WRITE(numout,*) 
    92          IF(lwp) WRITE(numout,*) 'lim_sbc_tau : LIM 3.0 sea-ice - surface ocean momentum fluxes' 
    93          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
    94       ENDIF 
    95  
    96       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
    97 !CDIR NOVERRCHK 
    98          DO jj = 2, jpjm1                         !* modulus of the ice-ocean velocity 
    99 !CDIR NOVERRCHK 
    100             DO ji = fs_2, fs_jpim1 
    101                zu_ij   = u_ice(ji  ,jj) - ssu_m(ji  ,jj)               ! (i  ,j) 
    102                zu_im1j = u_ice(ji-1,jj) - ssu_m(ji-1,jj)               ! (i-1,j) 
    103                zv_ij   = v_ice(ji,jj  ) - ssv_m(ji,jj  )               ! (i,j  ) 
    104                zv_ijm1 = v_ice(ji,jj-1) - ssv_m(ji,jj-1)               ! (i,j-1) 
    105                zztmp =  0.5 * ( zu_ij * zu_ij + zu_im1j * zu_im1j   &  ! mean of the square values instead 
    106                   &           + zv_ij * zv_ij + zv_ijm1 * zv_ijm1 )    ! of the square of the mean values... 
    107                ! WARNING, here we make an approximation for case 1 and 2: taum is not computed at each time 
    108                ! step but only every nn_fsbc time step. This avoid mutiple avarage to pass from T -> U,V grids 
    109                ! and next from U,V grids to T grid. Taum is used in tke, which should not be too sensitive to 
    110                ! this approximaton... 
    111                taum(ji,jj) = ( 1. - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * rhoco * zztmp 
    112                tmod_io(ji,jj) = SQRT( zztmp )  
    113             END DO 
    114          END DO 
    115          CALL lbc_lnk( taum, 'T', 1. )   ;   CALL lbc_lnk( tmod_io, 'T', 1. ) 
    116       ENDIF 
    117  
    118       SELECT CASE( kcpl ) 
    119          !                                        !--------------------------------! 
    120       CASE( 0 )                                   !  LIM 3 old stress computation  !  (at ice timestep only) 
    121          !                                        !--------------------------------!  
    122          !                                             !* ice stress over ocean with a ice-ocean rotation angle 
    123          DO jj = 1, jpjm1 
    124             zsang  = SIGN( 1.e0, gphif(1,jj) ) * sangvg         ! change the sinus angle sign in the south hemisphere 
    125             DO ji = 1, fs_jpim1 
    126                zu_u = u_ice(ji,jj) - u_oce(ji,jj)               ! ice velocity relative to the ocean 
    127                zv_v = v_ice(ji,jj) - v_oce(ji,jj) 
    128                !                                                ! quadratic drag formulation with rotation 
    129 !!gm still an error in the rotation, but usually the angle is zero (zsang=0, cangvg=1) 
    130                zutaui   = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * rhoco * ( cangvg * zu_u - zsang * zv_v )  
    131                zvtaui   = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * rhoco * ( cangvg * zv_v + zsang * zu_u )  
    132                !                                                ! bound for too large stress values 
    133                ! IMPORTANT: the exponential below prevents numerical oscillations in the ice-ocean stress 
    134                ! This is not physically based. A cleaner solution is offer in CASE kcpl=2 
    135                ztio_u(ji,jj) = zutaui * EXP( - ( tmod_io(ji,jj) + tmod_io(ji+1,jj) )  ) 
    136                ztio_v(ji,jj) = zvtaui * EXP( - ( tmod_io(ji,jj) + tmod_io(ji,jj+1) )  )  
    137             END DO 
    138          END DO 
    139          DO jj = 2, jpjm1                                       ! stress at the surface of the ocean 
    140             DO ji = fs_2, fs_jpim1   ! vertor opt. 
    141                zfrldu = 0.5 * ( ato_i(ji,jj) + ato_i(ji+1,jj  ) )   ! open-ocean fraction at U- & V-points (from T-point values) 
    142                zfrldv = 0.5 * ( ato_i(ji,jj) + ato_i(ji  ,jj+1) ) 
    143                !                                                    ! update surface ocean stress 
    144                utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * ztio_u(ji,jj) 
    145                vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * ztio_v(ji,jj) 
    146             END DO 
    147          END DO 
    148          CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition 
    149  
    150          ! 
    151          !                                        !--------------------------------! 
    152       CASE( 1 )                                   !  Use the now velocity          !  (at each ocean timestep) 
    153          !                                        !--------------------------------!  
    154          IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
    155             utau_oce(:,:) = utau(:,:)                 !* save the air-ocean stresses at ice time-step 
    156             vtau_oce(:,:) = vtau(:,:) 
    157          ENDIF 
    158          ! 
    159         DO jj = 2, jpjm1                              !* ice stress over ocean with a ice-ocean rotation angle 
    160             zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg        ! change the sinus angle sign in the south hemisphere 
    161             DO ji = fs_2, fs_jpim1 
    162                zat_u  = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5   ! ice area at u and V-points 
    163                zat_v  = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5 
    164                !                                                ! (u,v) ice-ocean velocity at (U,V)-point, resp. 
    165                zu_u   = u_ice(ji,jj) - ub(ji,jj,1) 
    166                zv_v   = v_ice(ji,jj) - vb(ji,jj,1) 
    167                !                                                ! quadratic drag formulation with rotation 
    168 !!gm still an error in the rotation, but usually the angle is zero (zsang=0, cangvg=1) 
    169                zutaui   = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * rhoco * ( cangvg * zu_u - zsang * zv_v )  
    170                zvtaui   = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * rhoco * ( cangvg * zv_v + zsang * zu_u )  
    171                !                                                   ! stress at the ocean surface 
    172                utau(ji,jj) = ( 1.- zat_u ) * utau_oce(ji,jj) + zat_u * zutaui 
    173                vtau(ji,jj) = ( 1.- zat_v ) * vtau_oce(ji,jj) + zat_v * zvtaui 
    174             END DO 
    175          END DO 
    176          CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition 
    177  
    178          ! 
    179          !                                        !--------------------------------! 
    180       CASE( 2 )                                   !  mixed 0 and 2 cases           !  (at each ocean timestep) 
    181          !                                        !--------------------------------!  
    182          IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
    183             utau_oce(:,:) = utau (:,:)               !* save the air-ocean stresses at ice time-step 
    184             vtau_oce(:,:) = vtau (:,:) 
    185             ssu_mb  (:,:) = ssu_m(:,:)               !* save the ice-ocean velocity at ice time-step 
    186             ssv_mb  (:,:) = ssv_m(:,:) 
    187          ENDIF 
    188          DO jj = 2, jpjm1                            !* ice stress over ocean with a ice-ocean rotation angle 
    189             zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg 
    190             DO ji = fs_2, fs_jpim1 
    191                zat_u = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5     ! ice area at u and V-points 
    192                zat_v = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5  
    193                ! 
    194                zu_ico = u_ice(ji,jj) - 0.5 * ( ub(ji,jj,1) + ssu_mb(ji,jj) )   ! ice-oce velocity using ub and ssu_mb 
    195                zv_ico = v_ice(ji,jj) - 0.5 * ( vb(ji,jj,1) + ssv_mb(ji,jj) ) 
    196                !                                        ! quadratic drag formulation with rotation 
    197 !!gm still an error in the rotation, but usually the angle is zero (zsang=0, cangvg=1) 
    198                zutaui = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * rhoco * ( cangvg * zu_ico - zsang * zv_ico ) 
    199                zvtaui = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * rhoco * ( cangvg * zv_ico + zsang * zu_ico ) 
    200                ! 
    201                utau(ji,jj) = ( 1.-zat_u ) * utau_oce(ji,jj) + zat_u * zutaui    ! stress at the ocean surface 
    202                vtau(ji,jj) = ( 1.-zat_v ) * vtau_oce(ji,jj) + zat_v * zvtaui 
    203             END DO 
    204          END DO 
    205          CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition 
    206          ! 
    207       END SELECT 
    208  
    209       IF(ln_ctl)   CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau   : ', mask1=umask,   & 
    210          &                       tab2d_2=vtau, clinfo2=' vtau    : '        , mask2=vmask ) 
    211       !   
    212    END SUBROUTINE lim_sbc_tau 
    213  
    21457 
    21558   SUBROUTINE lim_sbc_flx( kt ) 
     
    23578      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108. 
    23679      !!--------------------------------------------------------------------- 
    237       INTEGER ::   kt    ! number of iteration 
     80      INTEGER, INTENT(in) ::   kt    ! number of iteration 
    23881      !! 
    23982      INTEGER  ::   ji, jj           ! dummy loop indices 
     
    25497         IF(lwp) WRITE(numout,*) 'lim_sbc_flx : LIM 3.0 sea-ice - heat salt and mass ocean surface fluxes' 
    25598         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
     99         ! 
     100         r1_rdtice = 1. / rdt_ice 
     101         ! 
     102         soce_0(:,:) = soce 
     103         sice_0(:,:) = sice 
     104         ! 
     105         IF( cp_cfg == "orca" ) THEN           ! decrease ocean & ice reference salinities in the Baltic sea  
     106            WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.   & 
     107               &   54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp         )  
     108               soce_0(:,:) = 4._wp 
     109               sice_0(:,:) = 2._wp 
     110            END WHERE 
     111         ENDIF 
     112         ! 
    256113      ENDIF 
    257114 
     
    298155               ! fscmbq and ffltbif are obsolete 
    299156               !              &           + iflt * ffltbif(ji,jj) !!! only if one category is used 
    300                &           + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) / rdt_ice   & 
    301                &           + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) / rdt_ice                     & 
     157               &           + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice   & 
     158               &           + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) )                   * r1_rdtice   & 
    302159               &           + fhmec(ji,jj)     & ! new contribution due to snow melt in ridging!! 
    303160               &           + fheat_rpo(ji,jj) & ! contribution from ridge formation 
     
    340197               WRITE(numout,*) ' qcmif     : ', qcmif(jiindx,jjindx) 
    341198               WRITE(numout,*) ' qldif     : ', qldif(jiindx,jjindx) 
    342                WRITE(numout,*) ' qcmif / dt: ', qcmif(jiindx,jjindx) / rdt_ice 
    343                WRITE(numout,*) ' qldif / dt: ', qldif(jiindx,jjindx) / rdt_ice 
     199               WRITE(numout,*) ' qcmif / dt: ', qcmif(jiindx,jjindx) * r1_rdtice 
     200               WRITE(numout,*) ' qldif / dt: ', qldif(jiindx,jjindx) * r1_rdtice 
    344201               WRITE(numout,*) ' ifrdv     : ', ifrdv 
    345202               WRITE(numout,*) ' qfvbq     : ', qfvbq(jiindx,jjindx) 
    346203               WRITE(numout,*) ' qdtcn     : ', qdtcn(jiindx,jjindx) 
    347                WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) / rdt_ice 
    348                WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) / rdt_ice 
     204               WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) * r1_rdtice 
     205               WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) * r1_rdtice 
    349206               WRITE(numout,*) ' ' 
    350207               WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx) 
     
    380237               !              &   - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )  &   !  remov. snow precip over ice 
    381238               &   - sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   !  remov. snow precip over ice 
    382                &   - rdmsnif(ji,jj) / rdt_ice                &   !  freshwaterflux due to snow melting  
     239               &   - rdmsnif(ji,jj) * r1_rdtice                       &   !  freshwaterflux due to snow melting  
    383240               ! new contribution from snow falling when ridging 
    384241               &   + fmmec(ji,jj) 
     
    386243            !  computing salt exchanges at the ice/ocean interface 
    387244            !  sice should be the same as computed with the ice model 
    388             zfons =  ( soce - sice ) * ( rdmicif(ji,jj) / rdt_ice )  
     245            zfons =  ( soce_0(ji,jj) - sice_0(ji,jj) ) * rdmicif(ji,jj) * r1_rdtice  
    389246            ! SOCE 
    390             zfons =  ( sss_m(ji,jj) - sice ) * ( rdmicif(ji,jj) / rdt_ice )  
     247            zfons =  ( sss_m (ji,jj) - sice_0(ji,jj) ) * rdmicif(ji,jj) * r1_rdtice 
    391248 
    392249            !CT useless            !  salt flux for constant salinity 
     
    446303   END SUBROUTINE lim_sbc_flx 
    447304 
     305 
     306   SUBROUTINE lim_sbc_tau( kt , pu_oce, pv_oce ) 
     307      !!------------------------------------------------------------------- 
     308      !!                ***  ROUTINE lim_sbc_tau *** 
     309      !!   
     310      !! ** Purpose : Update the ocean surface stresses due to the ice 
     311      !!          
     312      !! ** Action  : * at each ice time step (every nn_fsbc time step): 
     313      !!                - compute the modulus of ice-ocean relative velocity  
     314      !!                  (*rho*Cd) at T-point (C-grid) or I-point (B-grid) 
     315      !!                      tmod_io = rhoco * | U_ice-U_oce | 
     316      !!                - update the modulus of stress at ocean surface 
     317      !!                      taum = frld * taum + (1-frld) * tmod_io * | U_ice-U_oce | 
     318      !!              * at each ocean time step (every kt):  
     319      !!                  compute linearized ice-ocean stresses as 
     320      !!                      Utau = tmod_io * | U_ice - pU_oce | 
     321      !!                using instantaneous current ocean velocity (usually before) 
     322      !! 
     323      !!    NB: - ice-ocean rotation angle no more allowed 
     324      !!        - here we make an approximation: taum is only computed every ice time step 
     325      !!          This avoids mutiple average to pass from T -> U,V grids and next from U,V grids  
     326      !!          to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton... 
     327      !! 
     328      !! ** Outputs : - utau, vtau   : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes 
     329      !!              - taum         : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes 
     330      !!--------------------------------------------------------------------- 
     331      INTEGER ,                     INTENT(in) ::   kt               ! ocean time-step index 
     332      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pu_oce, pv_oce   ! surface ocean currents 
     333      !! 
     334      INTEGER  ::   ji, jj   ! dummy loop indices 
     335      REAL(wp) ::   zat_u, zutau_ice, zu_t, zmodt   ! local scalar 
     336      REAL(wp) ::   zat_v, zvtau_ice, zv_t          !   -      - 
     337     !!--------------------------------------------------------------------- 
     338 
     339      IF( kt == nit000 ) THEN 
     340         IF(lwp) WRITE(numout,*) 
     341         IF(lwp) WRITE(numout,*) 'lim_sbc_tau : LIM-3 sea-ice - surface ocean momentum fluxes' 
     342         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
     343      ENDIF 
     344 
     345      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step) 
     346!CDIR NOVERRCHK 
     347         DO jj = 2, jpjm1                             !* update the modulus of stress at ocean surface (T-point) 
     348!CDIR NOVERRCHK 
     349            DO ji = fs_2, fs_jpim1 
     350               !                                               ! 2*(U_ice-U_oce) at T-point 
     351               zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)    
     352               zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1)  
     353               !                                              ! |U_ice-U_oce|^2 
     354               zmodt =  0.25_wp * (  zu_t * zu_t + zv_t * zv_t  ) 
     355               !                                               ! update the ocean stress modulus 
     356               taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * rhoco * zmodt 
     357               tmod_io(ji,jj) = rhoco * SQRT( zmodt )          ! rhoco * |U_ice-U_oce| at T-point 
     358            END DO 
     359         END DO 
     360         CALL lbc_lnk( taum, 'T', 1. )   ;   CALL lbc_lnk( tmod_io, 'T', 1. ) 
     361         ! 
     362         utau_oce(:,:) = utau(:,:)                    !* save the air-ocean stresses at ice time-step 
     363         vtau_oce(:,:) = vtau(:,:) 
     364         ! 
     365      ENDIF 
     366         ! 
     367         !                                     !==  every ocean time-step  ==! 
     368         ! 
     369      DO jj = 2, jpjm1                                !* update the stress WITHOUT a ice-ocean rotation angle 
     370         DO ji = fs_2, fs_jpim1   ! Vect. Opt. 
     371            zat_u  = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5_wp   ! ice area at u and V-points 
     372            zat_v  = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5_wp 
     373            !                                                   ! linearized quadratic drag formulation 
     374            zutau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) ) 
     375            zvtau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) ) 
     376            !                                                   ! stresses at the ocean surface 
     377            utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice 
     378            vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice 
     379         END DO 
     380      END DO 
     381      CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition 
     382      ! 
     383      IF(ln_ctl)   CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau   : ', mask1=umask,   & 
     384         &                       tab2d_2=vtau, clinfo2=' vtau    : '        , mask2=vmask ) 
     385      !   
     386   END SUBROUTINE lim_sbc_tau 
     387 
    448388#else 
    449389   !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.