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 918 for trunk/NEMO/LIM_SRC_3/limsbc.F90 – NEMO

Ignore:
Timestamp:
2008-05-09T12:00:09+02:00 (16 years ago)
Author:
rblod
Message:

Change ice_ocean stress computation, step II, see ticket #128

File:
1 edited

Legend:

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

    r917 r918  
    44   !!           computation of the flux at the sea ice/ocean interface 
    55   !!====================================================================== 
    6    !! History : 00-01 (H. Goosse) Original code 
    7    !!           02-07 (C. Ethe, G. Madec) re-writing F90 
    8    !!           06-07 (G. Madec) surface module 
     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 
    99   !!---------------------------------------------------------------------- 
    1010#if defined key_lim3 
     
    1212   !!   'key_lim3'                                    LIM 3.0 sea-ice model 
    1313   !!---------------------------------------------------------------------- 
    14    !!---------------------------------------------------------------------- 
    1514   !!   lim_sbc  : flux at the ice / ocean interface 
    1615   !!---------------------------------------------------------------------- 
     16   USE oce 
    1717   USE par_oce          ! ocean parameters 
    1818   USE par_ice          ! ice parameters 
     
    3535   PRIVATE 
    3636 
    37    PUBLIC lim_sbc       ! called by sbc_ice_lim 
     37   PUBLIC   lim_sbc_flx   ! called by sbc_ice_lim 
     38   PUBLIC   lim_sbc_tau   ! called by sbc_ice_lim 
    3839 
    3940   REAL(wp)  ::   epsi16 = 1.e-16  ! constant values 
     
    4142   REAL(wp)  ::   rone   = 1.e0 
    4243 
     44   REAL(wp), DIMENSION(jpi,jpj) ::   utau_oce   !: atm.-ocean surface i-stress (ocean referential)     [N/m2] 
     45   REAL(wp), DIMENSION(jpi,jpj) ::   vtau_oce   !: atm.-ocean surface j-stress (ocean referential)     [N/m2] 
    4346   !! * Substitutions 
    4447#  include "vectopt_loop_substitute.h90" 
    4548   !!---------------------------------------------------------------------- 
    46    !!   LIM 3.0,  UCL-LOCEAN-IPSL (2006)  
    47    !! $ Id: $ 
     49   !! NEMO/LIM  3.0 ,  UCL-LOCEAN-IPSL (2008)  
     50   !! $Id: $ 
    4851   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
    4952   !!---------------------------------------------------------------------- 
     
    5154CONTAINS 
    5255 
    53    SUBROUTINE lim_sbc( kt ) 
     56   SUBROUTINE lim_sbc_tau( kt, kcpl ) 
    5457      !!------------------------------------------------------------------- 
    55       !!                ***  ROUTINE lim_sbc *** 
     58      !!                ***  ROUTINE lim_sbc_tau *** 
    5659      !!   
    57       !! ** Purpose : Update surface ocean boundary condition over areas 
    58       !!      that are at least partially covered by sea-ice 
     60      !! ** Purpose : Update the ocean surface stresses due to the ice 
    5961      !!          
    60       !! ** Action  : - comput. of the momentum, heat and freshwater/salt 
    61       !!      fluxes at the ice-ocean interface. 
    62       !!              - Update  
     62      !! ** Action  : - compute the ice-ocean stress depending on kcpl: 
     63      !!              fluxes at the ice-ocean interface. 
     64      !!     Case 0  :  old LIM-3 way, computed at ice time-step only 
     65      !!     Case 1  :  at each ocean time step the stresses are computed 
     66      !!                using the current ocean velocity (now) 
     67      !!     Case 2  :  combination of half case 0 + half case 1 
     68      !!      
     69      !! ** Outputs : - utau   : sea surface i-stress (ocean referential) 
     70      !!              - vtau   : sea surface j-stress (ocean referential) 
     71      !! 
     72      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90. 
     73      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108. 
     74      !!--------------------------------------------------------------------- 
     75      INTEGER ::   kt     ! number of ocean iteration 
     76      INTEGER ::   kcpl   ! ice-ocean coupling indicator: =0  LIM-3 old case 
     77      !                   !                               =1  stresses computed using now ocean velocity 
     78      !                   !                               =2  combination of 0 and 1 cases 
     79      !! 
     80      INTEGER  ::   ji, jj           ! dummy loop indices 
     81      REAL(wp) ::   zfrldu, zfrldv   ! lead fraction at U- & V-points 
     82      REAL(wp) ::   ztglx , ztgly 
     83      REAL(wp) ::   zat_u, zu_ico, zutaui, zu_u, zv_u, zmodu, zmod 
     84      REAL(wp) ::   zat_v, zv_ico, zvtaui, zu_v, zv_v, zmodv, zsang 
     85       
     86#if defined key_coupled     
     87      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb     ! albedo of ice under overcast sky 
     88      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalbp    ! albedo of ice under clear sky 
     89#endif 
     90     REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! ocean stress below sea-ice 
     91      !!--------------------------------------------------------------------- 
     92      
     93      IF( kt == nit000 ) THEN 
     94         IF(lwp) WRITE(numout,*) 
     95         IF(lwp) WRITE(numout,*) 'lim_sbc_tau : LIM 3.0 sea-ice - surface ocean momentum fluxes' 
     96         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
     97      ENDIF 
     98 
     99      SELECT CASE( kcpl ) 
     100      !                                           !--------------------------------! 
     101      CASE( 0 )                                   !  LIM 3 old stress computation  !  (at ice timestep only) 
     102         !                                        !--------------------------------!  
     103         ! ... ice stress over ocean with a ice-ocean rotation angle 
     104         DO jj = 1, jpj 
     105            ! ... change the sinus angle sign in the south hemisphere 
     106            zsang  = SIGN(1.e0, gphif(1,jj) ) * sangvg 
     107            DO ji = 1, jpi 
     108               ! ... ice velocity relative to the ocean 
     109               zu_ico = u_ice(ji,jj) - u_oce(ji,jj) 
     110               zv_ico = v_ice(ji,jj) - v_oce(ji,jj) 
     111               zmod  = SQRT( zu_ico * zu_ico + zv_ico * zv_ico )  
     112               ! quadratic drag formulation 
     113               ztglx = rhoco * zmod * ( cangvg * zu_ico - zsang * zv_ico ) 
     114               ztgly = rhoco * zmod * ( cangvg * zv_ico + zsang * zu_ico ) 
     115               ! IMPORTANT 
     116               ! these lines are bound to prevent numerical oscillations 
     117               ! in the ice-ocean stress 
     118               ! They are physically ill-based. There is a cleaner solution 
     119               ! to try (remember discussion in Paris Gurvan) 
     120               ztio_u(ji,jj) = ztglx * exp( - zmod / 0.5 )  
     121               ztio_v(ji,jj) = ztgly * exp( - zmod / 0.5 )  
     122               ! 
     123            END DO 
     124         END DO 
     125 
     126         DO jj = 2, jpjm1 
     127            DO ji = fs_2, fs_jpim1   ! vertor opt. 
     128               ! ... open-ocean (lead) fraction at U- & V-points (from T-point values) 
     129               zfrldu = 0.5 * ( ato_i(ji,jj) + ato_i(ji+1,jj  ) ) 
     130               zfrldv = 0.5 * ( ato_i(ji,jj) + ato_i(ji  ,jj+1) ) 
     131               ! update surface ocean stress 
     132               utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * ztio_u(ji,jj) 
     133               vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * ztio_v(ji,jj) 
     134               ! 
     135            END DO 
     136         END DO 
     137 
     138         ! boundary condition on the stress (utau,vtau) 
     139         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. ) 
     140         ! 
     141         !                                        !--------------------------------! 
     142      CASE( 1 )                                   !  Use the now velocity          !  (at each ocean timestep) 
     143         !                                        !--------------------------------!  
     144         ! ... save the air-ocean stresses at ice time-step 
     145         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
     146            utau_oce(:,:) = utau(:,:) 
     147            vtau_oce(:,:) = vtau(:,:) 
     148         ENDIF 
     149         ! ... ice stress over ocean with a ice-ocean rotation angle 
     150         DO jj = 2, jpjm1 
     151            zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg 
     152            DO ji = fs_2, fs_jpim1 
     153               ! computation of wind stress over ocean in X and Y direction 
     154               zat_u  = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5      ! ice area at u and V-points 
     155               zat_v  = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5 
     156 
     157               zu_u   = u_ice(ji,jj) - un(ji,jj,1)                  ! u ice-ocean velocity at U-point 
     158               zv_v   = v_ice(ji,jj) - vn(ji,jj,1)                  ! v ice-ocean velocity at V-point 
     159               !                                                    ! u ice-ocean velocity at V-point 
     160               zu_v   = 0.25 * (  u_ice(ji,jj  ) - un(ji,jj  ,1)  +  u_ice(ji+1,jj  ) - un(ji+1,jj  ,1)    & 
     161                  &             + u_ice(ji,jj-1) - un(ji,jj-1,1)  +  u_ice(ji+1,jj-1) - un(ji+1,jj-1,1)  )  
     162               !                                                    ! v ice-ocean velocity at U-point 
     163               zv_u   = 0.25 * (  v_ice(ji-1,jj+1) - vn(ji-1,jj+1,1)  +  v_ice(ji,jj+1) - vn(ji,jj+1,1)    & 
     164                  &             + v_ice(ji-1,jj  ) - vn(ji-1,jj  ,1)  +  v_ice(ji,jj  ) - vn(ji,jj  ,1)  ) 
     165               zmodu    = SQRT( zu_u * zu_u + zv_u * zv_u )         ! module of the ice-ocean velocity at U-point 
     166               zmodv    = SQRT( zu_v * zu_v + zv_v * zv_v )         ! module of the ice-ocean velocity at V-point 
     167               zutaui   = rhoco * zmodu * ( cangvg * zu_u - zsang * zv_u )  
     168               zvtaui   = rhoco * zmodv * ( cangvg * zv_v + zsang * zu_v )  
     169 
     170               utau(ji,jj) = ( 1.-zat_u ) * utau_oce(ji,jj) + zat_u * zutaui    ! stress at the ocean surface 
     171               vtau(ji,jj) = ( 1.-zat_v ) * vtau_oce(ji,jj) + zat_v * zvtaui 
     172            END DO 
     173         END DO 
     174         ! boundary condition on the stress (utau,vtau) 
     175         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. ) 
     176         ! 
     177         !                                        !--------------------------------! 
     178      CASE( 2 )                                   !  mixed 0 and 2 cases           !  (at each ocean timestep) 
     179         !                                        !--------------------------------!  
     180         ! ... save the air-ocean stresses at ice time-step 
     181         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
     182            utau_oce(:,:) = utau(:,:) 
     183            vtau_oce(:,:) = vtau(:,:) 
     184         ENDIF 
     185         ! ... ice stress over ocean with a ice-ocean rotation angle 
     186         DO jj = 2, jpjm1 
     187            zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg 
     188            DO ji = fs_2, fs_jpim1 
     189               ! computation of wind stress over ocean in X and Y direction 
     190               zat_u = at_i(ji,jj) + at_i(ji+1,jj) * 0.5     ! ice area at u and V-points 
     191               zat_v = at_i(ji,jj) + at_i(ji,jj+1) * 0.5  
     192 
     193!!gm bug mixing U and V points value below     ====>>> to be corrected 
     194               zu_ico   = u_ice(ji,jj) - 0.5 * ( un(ji,jj,1) - ssu_m(ji,jj) )   ! ice-oce velocity using un and ssu_m 
     195               zv_ico   = v_ice(ji,jj) - 0.5 * ( vn(ji,jj,1) - ssu_m(ji,jj) ) 
     196               !                                        ! module of the ice-ocean velocity 
     197               zmod     = SQRT( zu_ico * zu_ico + zv_ico * zv_ico ) 
     198               !                                        ! quadratic drag formulation with rotation 
     199               zutaui   = rhoco * zmod * ( cangvg * zu_ico - zsang * zv_ico ) 
     200               zvtaui   = rhoco * zmod * ( cangvg * zv_ico + zsang * zu_ico ) 
     201! 
     202               utau(ji,jj) = ( 1.-zat_u ) * utau_oce(ji,jj) + zat_u * zutaui    ! stress at the ocean surface 
     203               vtau(ji,jj) = ( 1.-zat_v ) * vtau_oce(ji,jj) + zat_v * zvtaui 
     204            END DO 
     205         END DO 
     206         ! boundary condition on the stress (utau,vtau) 
     207         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. ) 
     208         ! 
     209      END SELECT 
     210      IF(ln_ctl)   CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau   : ', mask1=umask,   & 
     211         &                       tab2d_2=vtau, clinfo2=' vtau    : '        , mask2=vmask ) 
     212      !   
     213   END SUBROUTINE lim_sbc_tau 
     214 
     215 
     216   SUBROUTINE lim_sbc_flx( kt ) 
     217      !!------------------------------------------------------------------- 
     218      !!                ***  ROUTINE lim_sbc_flx *** 
     219      !!   
     220      !! ** Purpose :   Update the surface ocean boundary condition for heat  
     221      !!              salt and mass over areas where sea-ice is non-zero 
     222      !!          
     223      !! ** Action  : - computes the heat and freshwater/salt fluxes 
     224      !!              at the ice-ocean interface. 
     225      !!              - Update the ocean sbc 
    63226      !!      
    64227      !! ** Outputs : - qsr    : sea heat flux:     solar  
     
    66229      !!              - emp    : freshwater budget: volume flux  
    67230      !!              - emps   : freshwater budget: concentration/dillution  
    68       !!              - utau   : sea surface i-stress (ocean referential) 
    69       !!              - vtau   : sea surface j-stress (ocean referential) 
    70231      !! 
    71232      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90. 
     
    80241      REAL(wp) ::   zfons            ! salt exchanges at the ice/ocean interface 
    81242      REAL(wp) ::   zpme             ! freshwater exchanges at the ice/ocean interface 
    82       REAL(wp) ::   zfrldu, zfrldv   ! lead fraction at U- & V-points 
    83       REAL(wp) ::   zu_io , zv_io    ! 2 components of the ice-ocean velocity 
    84       REAL(wp) ::   ztglx , ztgly    ! temporary scalar 
    85        
     243      REAL(wp), DIMENSION(jpi,jpj) ::   zfcm1 , zfcm2    ! solar/non solar heat fluxes 
    86244#if defined key_coupled     
    87245      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb     ! albedo of ice under overcast sky 
    88246      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalbp    ! albedo of ice under clear sky 
    89247#endif 
    90       REAL(wp) ::   zsang, zmod 
    91       REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! ocean stress below sea-ice 
    92       REAL(wp), DIMENSION(jpi,jpj) ::   zfcm1 , zfcm2    ! solar/non solar heat fluxes 
    93  
    94248      !!--------------------------------------------------------------------- 
    95249      
    96250      IF( kt == nit000 ) THEN 
    97251         IF(lwp) WRITE(numout,*) 
    98          IF(lwp) WRITE(numout,*) 'lim_sbc : LIM 3.0 sea-ice - surface boundary condition' 
    99          IF(lwp) WRITE(numout,*) '~~~~~~~  ' 
     252         IF(lwp) WRITE(numout,*) 'lim_sbc_flx : LIM 3.0 sea-ice - heat salt and mass ocean surface fluxes' 
     253         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
    100254      ENDIF 
    101255 
     
    169323!                           ! fdtcn : turbulent oceanic heat flux 
    170324 
    171  
     325!!gm   this IF prevents the vertorisation of the whole loop 
    172326            IF ( ( ji .EQ. jiindx ) .AND. ( jj .EQ. jjindx) ) THEN 
    173327               WRITE(numout,*) ' lim_sbc : heat fluxes ' 
     
    198352               WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx) 
    199353            ENDIF 
     354!!gm   end 
    200355         END DO 
    201356      END DO 
     
    205360      !------------------------------------------! 
    206361 
     362!!gm   optimisation: this loop have to be merged with the previous one 
    207363      DO jj = 1, jpj 
    208364         DO ji = 1, jpi 
     
    254410      END DO 
    255411 
    256       emps(:,:) = fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:) 
    257       IF (num_sal.eq.2) THEN 
    258          !In case of variable salinity the salt flux has to be accounted for differently 
    259          ! Brine drainage has to be added 
     412      IF( num_sal == 2 ) THEN      ! variable ice salinity: brine drainage included in the salt flux 
    260413         emps(:,:) = fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:) 
     414      ELSE                         ! constant ice salinity: 
     415         emps(:,:) =              fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:) 
    261416      ENDIF 
    262417       
    263418      IF( lk_dynspg_rl )    emp (:,:) = emps(:,:)      ! rigid-lid formulation : emp = emps 
    264  
    265       !------------------------------------------! 
    266       !    momentum flux at the ocean surface    ! 
    267       !------------------------------------------! 
    268  
    269       IF ( ln_limdyn ) THEN                        ! Update the stress over ice-over area (only in ice-dynamic case) 
    270          !                                         ! otherwise the atmosphere-ocean stress is used everywhere 
    271  
    272          ! ... ice stress over ocean with a ice-ocean rotation angle 
    273          DO jj = 1, jpj 
    274             ! ... change the sinus angle sign in the south hemisphere 
    275             zsang  = SIGN(1.e0, gphif(1,jj) ) * sangvg 
    276             DO ji = 1, jpi 
    277                ! ... ice velocity relative to the ocean 
    278                zu_io = u_ice(ji,jj) - u_oce(ji,jj) 
    279                zv_io = v_ice(ji,jj) - v_oce(ji,jj) 
    280                zmod  = SQRT( zu_io * zu_io + zv_io * zv_io )  
    281                ! quadratic drag formulation 
    282                ztglx = rhoco * zmod * ( cangvg * zu_io - zsang * zv_io ) 
    283                ztgly = rhoco * zmod * ( cangvg * zv_io + zsang * zu_io ) 
    284                ! IMPORTANT 
    285                ! these lines are bound to prevent numerical oscillations 
    286                ! in the ice-ocean stress 
    287                ! They are physically ill-based. There is a cleaner solution 
    288                ! to try (remember discussion in Paris Gurvan) 
    289                ztio_u(ji,jj) = ztglx * exp( - zmod / 0.5 )  
    290                ztio_v(ji,jj) = ztgly * exp( - zmod / 0.5 )  
    291                ! 
    292             END DO 
    293          END DO 
    294  
    295          DO jj = 2, jpjm1 
    296             DO ji = fs_2, fs_jpim1   ! vertor opt. 
    297                ! ... open-ocean (lead) fraction at U- & V-points (from T-point values) 
    298                zfrldu = 0.5 * ( ato_i(ji,jj) + ato_i(ji+1,jj  ) ) 
    299                zfrldv = 0.5 * ( ato_i(ji,jj) + ato_i(ji  ,jj+1) ) 
    300                ! update surface ocean stress 
    301                utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * ztio_u(ji,jj) 
    302                vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * ztio_v(ji,jj) 
    303                ! 
    304             END DO 
    305          END DO 
    306  
    307          ! boundary condition on the stress (utau,vtau) 
    308          CALL lbc_lnk( utau, 'U', -1. ) 
    309          CALL lbc_lnk( vtau, 'V', -1. ) 
    310  
    311       ENDIF 
    312419 
    313420      !-----------------------------------------------! 
     
    331438 
    332439      IF(ln_ctl) THEN 
    333          CALL prt_ctl(tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns   , clinfo2=' qns     : ') 
    334          CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps  , clinfo2=' emps    : ') 
    335          CALL prt_ctl(tab2d_1=utau  , clinfo1=' lim_sbc: utau   : ', mask1=umask,   & 
    336             &         tab2d_2=vtau  , clinfo2=' vtau    : '        , mask2=vmask ) 
    337          CALL prt_ctl(tab2d_1=freeze, clinfo1=' lim_sbc: freeze : ') 
    338          CALL prt_ctl(tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl) 
     440         CALL prt_ctl( tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns , clinfo2=' qns     : ' ) 
     441         CALL prt_ctl( tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps, clinfo2=' emps    : ' ) 
     442         CALL prt_ctl( tab2d_1=freeze, clinfo1=' lim_sbc: freeze : ' ) 
     443         CALL prt_ctl( tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl ) 
    339444      ENDIF  
    340     
    341     END SUBROUTINE lim_sbc 
     445      !  
     446   END SUBROUTINE lim_sbc_flx 
     447 
    342448 
    343449#else 
Note: See TracChangeset for help on using the changeset viewer.