Changeset 1322


Ignore:
Timestamp:
2009-02-20T09:50:55+01:00 (12 years ago)
Author:
ctlod
Message:

correction of the wind module computation + few bugs, see ticket: #346

File:
1 edited

Legend:

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

    r1146 r1322  
    4141   REAL(wp)  ::   rone   = 1.e0 
    4242 
    43    REAL(wp), DIMENSION(jpi,jpj) ::   utau_oce   !: atm.-ocean surface i-stress (ocean referential)     [N/m2] 
    44    REAL(wp), DIMENSION(jpi,jpj) ::   vtau_oce   !: atm.-ocean surface j-stress (ocean referential)     [N/m2] 
     43   REAL(wp), DIMENSION(jpi,jpj) ::   utau_oce, vtau_oce   !: air-ocean surface i- & j-stress            [N/m2] 
     44   REAL(wp), DIMENSION(jpi,jpj) ::   tmod_io              !: modulus of the ice-ocean relative velocity   [m/s] 
     45   REAL(wp), DIMENSION(jpi,jpj) ::   ssu_mb  , ssv_mb     !: before mean ocean surface currents          [m/s] 
    4546   !! * Substitutions 
    4647#  include "vectopt_loop_substitute.h90" 
     
    7980      INTEGER  ::   ji, jj           ! dummy loop indices 
    8081      REAL(wp) ::   zfrldu, zfrldv   ! lead fraction at U- & V-points 
    81       REAL(wp) ::   ztglx , ztgly 
    82       REAL(wp) ::   zat_u, zu_ico, zutaui, zu_u, zv_u, zmodu, zmod 
    83       REAL(wp) ::   zat_v, zv_ico, zvtaui, zu_v, zv_v, zmodv, zsang 
     82      REAL(wp) ::   zat_u, zu_ico, zutaui, zu_u 
     83      REAL(wp) ::   zat_v, zv_ico, zvtaui, zv_v, zsang 
     84      REAL(wp) ::   zu_ij, zu_im1j, zv_ij, zv_ijm1 
    8485 
    8586#if defined key_coupled     
     
    9798 
    9899      SELECT CASE( kcpl ) 
    99          !                                           !--------------------------------! 
     100         !                                        !--------------------------------! 
    100101      CASE( 0 )                                   !  LIM 3 old stress computation  !  (at ice timestep only) 
    101102         !                                        !--------------------------------!  
     103         DO jj = 2, jpjm1                      ! ... modulus of the ice-ocean velocity 
     104            DO ji = fs_2, fs_jpim1 
     105               zu_ij   = u_ice(ji  ,jj) - ssu_m(ji  ,jj)               ! (i  ,j) 
     106               zu_im1j = u_ice(ji-1,jj) - ssu_m(ji-1,jj)               ! (i-1,j) 
     107               zv_ij   = v_ice(ji,jj  ) - ssv_m(ji,jj  )               ! (i,j  ) 
     108               zv_ijm1 = v_ice(ji,jj-1) - ssv_m(ji,jj-1)               ! (i,j-1) 
     109               tmod_io(ji,jj) = SQRT(  0.5 * ( zu_ij * zu_ij + zu_im1j * zu_im1j   & 
     110                  &                          + zv_ij * zv_ij + zv_ijm1 * zv_ijm1  ) ) 
     111            END DO 
     112         END DO 
     113         CALL lbc_lnk( tmod_io, 'T', 1. ) 
    102114         ! ... ice stress over ocean with a ice-ocean rotation angle 
    103115         DO jj = 1, jpj 
     
    106118            DO ji = 1, jpi 
    107119               ! ... ice velocity relative to the ocean 
    108                zu_ico = u_ice(ji,jj) - u_oce(ji,jj) 
    109                zv_ico = v_ice(ji,jj) - v_oce(ji,jj) 
    110                zmod  = SQRT( zu_ico * zu_ico + zv_ico * zv_ico )  
     120               zu_u = u_ice(ji,jj) - u_oce(ji,jj) 
     121               zv_v = v_ice(ji,jj) - v_oce(ji,jj) 
    111122               ! quadratic drag formulation 
    112                ztglx = rhoco * zmod * ( cangvg * zu_ico - zsang * zv_ico ) 
    113                ztgly = rhoco * zmod * ( cangvg * zv_ico + zsang * zu_ico ) 
     123!!gm still an error in the rotation, but usually the angle is zero (zsang=0, cangvg=1) 
     124               zutaui   = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * rhoco * ( cangvg * zu_u - zsang * zv_v )  
     125               zvtaui   = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * rhoco * ( cangvg * zv_v + zsang * zu_u )  
    114126               ! IMPORTANT 
    115127               ! these lines are bound to prevent numerical oscillations 
     
    117129               ! They are physically ill-based. There is a cleaner solution 
    118130               ! to try (remember discussion in Paris Gurvan) 
    119                ztio_u(ji,jj) = ztglx * exp( - zmod / 0.5 )  
    120                ztio_v(ji,jj) = ztgly * exp( - zmod / 0.5 )  
     131               ztio_u(ji,jj) = zutaui * exp( - ( tmod_io(ji,jj) + tmod_io(ji+1,jj) )  ) 
     132               ztio_v(ji,jj) = zvtaui * exp( - ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) )  
    121133               ! 
    122134            END DO 
     
    141153      CASE( 1 )                                   !  Use the now velocity          !  (at each ocean timestep) 
    142154         !                                        !--------------------------------!  
    143          ! ... save the air-ocean stresses at ice time-step 
    144155         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
    145             utau_oce(:,:) = utau(:,:) 
     156            utau_oce(:,:) = utau(:,:)             ! ... save the air-ocean stresses at ice time-step 
    146157            vtau_oce(:,:) = vtau(:,:) 
     158            DO jj = 2, jpjm1                      ! ... modulus of the ice-ocean velocity 
     159               DO ji = fs_2, fs_jpim1 
     160                  zu_ij   = u_ice(ji  ,jj) - ssu_m(ji  ,jj)               ! (i  ,j) 
     161                  zu_im1j = u_ice(ji-1,jj) - ssu_m(ji-1,jj)               ! (i-1,j) 
     162                  zv_ij   = v_ice(ji,jj  ) - ssv_m(ji,jj  )               ! (i,j  ) 
     163                  zv_ijm1 = v_ice(ji,jj-1) - ssv_m(ji,jj-1)               ! (i,j-1) 
     164                  tmod_io(ji,jj) = SQRT( 0.5 * (  zu_ij * zu_ij + zu_im1j * zu_im1j   & 
     165                     &                          + zv_ij * zv_ij + zv_ijm1 * zv_ijm1 ) ) 
     166               END DO 
     167            END DO 
     168         CALL lbc_lnk( tmod_io, 'T', 1. ) 
    147169         ENDIF 
    148170         ! ... ice stress over ocean with a ice-ocean rotation angle 
     
    154176               zat_v  = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5 
    155177 
    156                zu_u   = u_ice(ji,jj) - un(ji,jj,1)                  ! u ice-ocean velocity at U-point 
    157                zv_v   = v_ice(ji,jj) - vn(ji,jj,1)                  ! v ice-ocean velocity at V-point 
     178               zu_u   = u_ice(ji,jj) - ub(ji,jj,1)                  ! u ice-ocean velocity at U-point 
     179               zv_v   = v_ice(ji,jj) - vb(ji,jj,1)                  ! v ice-ocean velocity at V-point 
     180               ! 
     181!!gm still an error in the rotation, but usually the angle is zero (zsang=0, cangvg=1) 
     182               zutaui   = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * rhoco * ( cangvg * zu_u - zsang * zv_v )  
     183               zvtaui   = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * rhoco * ( cangvg * zv_v + zsang * zu_u )  
     184 
     185               utau(ji,jj) = ( 1.- zat_u ) * utau_oce(ji,jj) + zat_u * zutaui    ! stress at the ocean surface 
     186               vtau(ji,jj) = ( 1.- zat_v ) * vtau_oce(ji,jj) + zat_v * zvtaui 
    158187               !                                                    ! u ice-ocean velocity at V-point 
    159                zu_v   = 0.25 * (  u_ice(ji,jj  ) - un(ji,jj  ,1)  +  u_ice(ji+1,jj  ) - un(ji+1,jj  ,1)    & 
    160                   &             + u_ice(ji,jj-1) - un(ji,jj-1,1)  +  u_ice(ji+1,jj-1) - un(ji+1,jj-1,1)  )  
    161                !                                                    ! v ice-ocean velocity at U-point 
    162                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)    & 
    163                   &             + v_ice(ji-1,jj  ) - vn(ji-1,jj  ,1)  +  v_ice(ji,jj  ) - vn(ji,jj  ,1)  ) 
    164                zmodu    = SQRT( zu_u * zu_u + zv_u * zv_u )         ! module of the ice-ocean velocity at U-point 
    165                zmodv    = SQRT( zu_v * zu_v + zv_v * zv_v )         ! module of the ice-ocean velocity at V-point 
    166                zutaui   = rhoco * zmodu * ( cangvg * zu_u - zsang * zv_u )  
    167                zvtaui   = rhoco * zmodv * ( cangvg * zv_v + zsang * zu_v )  
    168  
    169                utau(ji,jj) = ( 1.-zat_u ) * utau_oce(ji,jj) + zat_u * zutaui    ! stress at the ocean surface 
    170                vtau(ji,jj) = ( 1.-zat_v ) * vtau_oce(ji,jj) + zat_v * zvtaui 
    171188            END DO 
    172189         END DO 
     
    177194      CASE( 2 )                                   !  mixed 0 and 2 cases           !  (at each ocean timestep) 
    178195         !                                        !--------------------------------!  
    179          ! ... save the air-ocean stresses at ice time-step 
    180196         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN 
    181             utau_oce(:,:) = utau(:,:) 
    182             vtau_oce(:,:) = vtau(:,:) 
     197            utau_oce(:,:) = utau (:,:)            ! ... save the air-ocean stresses at ice time-step 
     198            vtau_oce(:,:) = vtau (:,:) 
     199            ssu_mb  (:,:) = ssu_m(:,:)            ! ... save the ice-ocean velocity at ice time-step 
     200            ssv_mb  (:,:) = ssv_m(:,:) 
     201            DO jj = 2, jpjm1                      ! ... modulus of the ice-ocean velocity 
     202               DO ji = fs_2, fs_jpim1 
     203                  zu_ij   = u_ice(ji  ,jj) - ssu_m(ji  ,jj)               ! (i  ,j) 
     204                  zu_im1j = u_ice(ji-1,jj) - ssu_m(ji-1,jj)               ! (i-1,j) 
     205                  zv_ij   = v_ice(ji,jj  ) - ssv_m(ji,jj  )               ! (i,j  ) 
     206                  zv_ijm1 = v_ice(ji,jj-1) - ssv_m(ji,jj-1)               ! (i,j-1) 
     207                  tmod_io(ji,jj) = SQRT( 0.5 * ( zu_ij * zu_ij + zu_im1j * zu_im1j   & 
     208                     &                         + zv_ij * zv_ij + zv_ijm1 * zv_ijm1 ) ) 
     209               END DO 
     210            END DO 
     211         CALL lbc_lnk( tmod_io, 'T', 1. ) 
    183212         ENDIF 
    184213         ! ... ice stress over ocean with a ice-ocean rotation angle 
     
    187216            DO ji = fs_2, fs_jpim1 
    188217               ! computation of wind stress over ocean in X and Y direction 
    189                zat_u = at_i(ji,jj) + at_i(ji+1,jj) * 0.5     ! ice area at u and V-points 
    190                zat_v = at_i(ji,jj) + at_i(ji,jj+1) * 0.5  
     218               zat_u = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5     ! ice area at u and V-points 
     219               zat_v = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5  
    191220 
    192221               !!gm bug mixing U and V points value below     ====>>> to be corrected 
    193                zu_ico   = u_ice(ji,jj) - 0.5 * ( un(ji,jj,1) - ssu_m(ji,jj) )   ! ice-oce velocity using un and ssu_m 
    194                zv_ico   = v_ice(ji,jj) - 0.5 * ( vn(ji,jj,1) - ssu_m(ji,jj) ) 
    195                !                                        ! module of the ice-ocean velocity 
    196                zmod     = SQRT( zu_ico * zu_ico + zv_ico * zv_ico ) 
     222               zu_ico = u_ice(ji,jj) - 0.5 * ( ub(ji,jj,1) - ssu_mb(ji,jj) )   ! ice-oce velocity using un and ssu_mb 
     223               zv_ico = v_ice(ji,jj) - 0.5 * ( vb(ji,jj,1) - ssv_mb(ji,jj) ) 
    197224               !                                        ! quadratic drag formulation with rotation 
    198                zutaui   = rhoco * zmod * ( cangvg * zu_ico - zsang * zv_ico ) 
    199                zvtaui   = rhoco * zmod * ( cangvg * zv_ico + zsang * zu_ico ) 
     225!!gm still an error in the rotation, but usually the angle is zero (zsang=0, cangvg=1) 
     226               zutaui = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * rhoco * ( cangvg * zu_ico - zsang * zv_ico ) 
     227               zvtaui = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * rhoco * ( cangvg * zv_ico + zsang * zu_ico ) 
    200228               ! 
    201229               utau(ji,jj) = ( 1.-zat_u ) * utau_oce(ji,jj) + zat_u * zutaui    ! stress at the ocean surface 
Note: See TracChangeset for help on using the changeset viewer.