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 1662 – NEMO

Changeset 1662


Ignore:
Timestamp:
2009-10-14T18:51:06+02:00 (15 years ago)
Author:
rblod
Message:

Correction of major bug on bottom friction, see ticket #233

Location:
trunk/NEMO/OPA_SRC
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r1502 r1662  
    2222   USE phycst          ! physical constants 
    2323   USE domvvl          ! variable volume 
     24   USE zdfbfr          ! bottom friction 
    2425   USE obcdta          ! open boundary condition data      
    2526   USE obcfla          ! Flather open boundary condition   
    … …  
    9192      INTEGER, INTENT(in)  ::   kt   ! ocean time-step index 
    9293      !! 
    93       INTEGER  ::   ji, jj, jk, jn             ! dummy loop indices 
    94       INTEGER  ::   icycle                      ! temporary scalar 
    95       REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b,   &  ! temporary scalars 
    96          z1_8, zspgu, zcubt, zx1, zy1,    &  !     "        " 
    97          z1_4, zspgv, zcvbt, zx2, zy2        !     "        " 
     94      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
     95      INTEGER  ::   icycle           ! temporary scalar 
     96 
     97      REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b     ! temporary scalars 
     98      REAL(wp) ::   z1_8, zx1, zy1                   !    -         - 
     99      REAL(wp) ::   z1_4, zx2, zy2                   !     -         - 
     100      REAL(wp) ::   zu_spg, zu_cor, zu_sld, zu_asp   !     -         - 
     101      REAL(wp) ::   zv_spg, zv_cor, zv_sld, zv_asp   !     -         - 
     102      !! 
    98103      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv, zsshb_e 
    99       !!  
     104      !! 
    100105      REAL(wp), DIMENSION(jpi,jpj) ::   zsshun_e, zsshvn_e   ! 2D workspace 
    101106      !! 
    … …  
    144149      ! 
    145150      zhdiv(:,:) = 0.e0                                     ! barotropic divergence 
     151      zu_sld = 0.e0   ;   zu_asp = 0.e0                     ! tides trends (lk_tide=F) 
     152      zv_sld = 0.e0   ;   zv_asp = 0.e0 
    146153 
    147154      ! ----------------------------------------------------------------------------- 
    … …  
    245252      END DO 
    246253 
     254      !                                            ! Remove barotropic contribution of bottom friction  
     255                                                   ! from the barotropic transport trend 
     256      IF( lk_vvl ) THEN 
     257       DO jj = 2, jpjm1 
     258          DO ji = fs_2, fs_jpim1   ! vector opt. 
     259             zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * zub(ji,jj) / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) ) 
     260             zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * zvb(ji,jj) / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) ) 
     261          END DO 
     262       END DO 
     263      ELSE 
     264       DO jj = 2, jpjm1 
     265          DO ji = fs_2, fs_jpim1   ! vector opt. 
     266             zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * zub(ji,jj) * hur(ji,jj) 
     267             zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * zvb(ji,jj) * hvr(ji,jj) 
     268          END DO 
     269       END DO 
     270      ENDIF 
     271 
    247272      !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity) 
    248273      !                                   ! --------------------------  
    … …  
    272297      hur_e  (:,:) = hur  (:,:)            ! ocean depth inverted at u- and v-points 
    273298      hvr_e  (:,:) = hvr  (:,:) 
    274       zsshb_e(:,:) = sshn (:,:)            ! sea surface height (before and now) 
     299!RBbug     zsshb_e(:,:) = sshn (:,:)   
     300      zsshb_e(:,:) = sshn_b(:,:)           ! sea surface height (before and now) 
    275301      sshn_e (:,:) = sshn (:,:) 
    276302       
    … …  
    343369                  ! surface pressure gradient 
    344370                  IF( lk_vvl) THEN 
    345                      zspgu = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
    346                         &             - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
    347                      zspgv = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
    348                         &             - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
     371                     zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
     372                        &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
     373                     zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
     374                        &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
    349375                  ELSE 
    350                      zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
    351                      zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
     376                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
     377                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    352378                  ENDIF 
    353379                  ! energy conserving formulation for planetary vorticity term 
    … …  
    356382                  zx1 = ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) ) / e2v(ji,jj) 
    357383                  zx2 = ( zwx(ji  ,jj  ) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    358                   zcubt = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj) 
    359                   zcvbt =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) 
    360                   ! after barotropic velocity 
    361                   ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 
    362                   va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 
     384                  zu_cor = z1_4 * ( ff(ji  ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj) 
     385                  zv_cor =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) 
     386                  ! after velocities with implicit bottom friction 
     387                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1)   & 
     388                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
     389                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1)   & 
     390                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
    363391               END DO 
    364392            END DO 
    … …  
    369397                   ! surface pressure gradient 
    370398                  IF( lk_vvl) THEN 
    371                      zspgu = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
    372                         &             - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
    373                      zspgv = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
    374                         &             - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
     399                     zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
     400                        &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
     401                     zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
     402                        &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
    375403                  ELSE 
    376                      zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
    377                      zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
     404                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
     405                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    378406                  ENDIF 
    379407                  ! enstrophy conserving formulation for planetary vorticity term 
    380408                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
    381409                  zx1 = - z1_8 * ( zwx(ji-1,jj  ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji  ,jj+1) ) / e2v(ji,jj) 
    382                   zcubt  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj) 
    383                   zcvbt  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj) 
    384                   ! after barotropic velocity 
    385                   ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 
    386                   va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 
     410                  zu_cor  = zy1 * ( ff(ji  ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj) 
     411                  zv_cor  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj) 
     412                  ! after velocities with implicit bottom friction 
     413                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1)   & 
     414                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
     415                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1)   & 
     416                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
    387417               END DO 
    388418            END DO 
    … …  
    393423                  ! surface pressure gradient 
    394424                  IF( lk_vvl) THEN 
    395                      zspgu = -grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
    396                         &             - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
    397                      zspgv = -grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
    398                         &             - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
     425                     zu_spg = -grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
     426                        &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
     427                     zv_spg = -grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
     428                        &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
    399429                  ELSE 
    400                      zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
    401                      zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
     430                     zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
     431                     zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    402432                  ENDIF 
    403433                  ! energy/enstrophy conserving formulation for planetary vorticity term 
    404                   zcubt = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   & 
     434                  zu_cor = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   & 
    405435                     &                           + ftse(ji,jj  ) * zwy(ji  ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj) 
    406                   zcvbt = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   & 
     436                  zv_cor = - z1_4 / e2v(ji,jj) * (  ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji  ,jj+1)   & 
    407437                     &                           + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) * hvr_e(ji,jj) 
    408                   ! after barotropic velocity 
    409                   ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 
    410                   va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 
     438                  ! after velocities with implicit bottom friction 
     439                  ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1)   & 
     440                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
     441                  va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1)   & 
     442                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
    411443               END DO 
    412444            END DO 
    … …  
    499531      IF( lrst_oce )   CALL ts_rst( kt, 'WRITE' ) 
    500532      ! 
     533      ! 
    501534   END SUBROUTINE dyn_spg_ts 
    502535 

  • trunk/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r1146 r1662  
    7373      INTEGER ::   ji, jj, jk                          ! dummy loop indices 
    7474      REAL(wp) ::   zrau0r, z2dtf, zcoef, zzws, zrhs   ! temporary scalars 
     75      REAL(wp) ::   zzwi                               ! temporary scalars 
    7576      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwi        ! temporary workspace arrays 
    7677      !!---------------------------------------------------------------------- 
    … …  
    8990      ! --------------------------- 
    9091      ! Matrix and second member construction 
    91       ! bottom boundary condition: only zws must be masked as avmu can take 
     92      ! bottom boundary condition: both zwi and zws must be masked as avmu can take 
    9293      ! non zero value at the ocean bottom depending on the bottom friction 
    93       ! used (see zdfmix.F) 
     94      ! used but the bottom velocities have already been updated with the bottom 
     95      ! friction velocity in dyn_bfr using values from the previous timestep. There 
     96      ! is no need to include these in the implicit calculation. 
    9497      DO jk = 1, jpkm1 
    9598         DO jj = 2, jpjm1  
    9699            DO ji = fs_2, fs_jpim1   ! vector opt. 
    97100               zcoef = - p2dt / fse3u(ji,jj,jk) 
    98                zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
     101               zzwi          = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
     102               zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk) 
    99103               zzws          = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
    100104               zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1) 
    … …  
    183187      ! --------------------------- 
    184188      ! Matrix and second member construction 
    185       ! bottom boundary condition: only zws must be masked as avmv can take 
     189      ! bottom boundary condition: both zwi and zws must be masked as avmv can take 
    186190      ! non zero value at the ocean bottom depending on the bottom friction 
    187       ! used (see zdfmix.F) 
     191      ! used but the bottom velocities have already been updated with the bottom 
     192      ! friction velocity in dyn_bfr using values from the previous timestep. There 
     193      ! is no need to include these in the implicit calculation. 
    188194      DO jk = 1, jpkm1 
    189195         DO jj = 2, jpjm1    
    190196            DO ji = fs_2, fs_jpim1   ! vector opt. 
    191197               zcoef = -p2dt / fse3v(ji,jj,jk) 
    192                zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  ) 
     198               zzwi          = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  ) 
     199               zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk) 
    193200               zzws       = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 
    194201               zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1) 

  • trunk/NEMO/OPA_SRC/TRD/trdmod.F90

    r1601 r1662  
    127127                     ztswu(ji,jj) = utau(ji,jj) / ( fse3u(ji,jj,1)*rau0 ) 
    128128                     ztswv(ji,jj) = vtau(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) 
    129                      ! save bottom friction momentum fluxes 
    130                      ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) 
    131                      ikbv   = MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) ) 
    132                      ikbum1 = MAX( ikbu-1, 1 ) 
    133                      ikbvm1 = MAX( ikbv-1, 1 ) 
    134                      zua = ua(ji,jj,ikbum1) * r2dt + ub(ji,jj,ikbum1) 
    135                      zva = va(ji,jj,ikbvm1) * r2dt + vb(ji,jj,ikbvm1) 
    136                      ztbfu(ji,jj) = - avmu(ji,jj,ikbu) * zua / ( fse3u(ji,jj,ikbum1)*fse3uw(ji,jj,ikbu) ) 
    137                      ztbfv(ji,jj) = - avmv(ji,jj,ikbv) * zva / ( fse3v(ji,jj,ikbvm1)*fse3vw(ji,jj,ikbv) ) 
     129                     ! bottom friction contribution now handled explicitly 
    138130                     ! 
    139131                     ptrdx(ji,jj,1     ) = ptrdx(ji,jj,1     ) - ztswu(ji,jj) 
    … …  
    146138               CALL trd_icp( ptrdx, ptrdy, jpicpd_zdf, ctype )    
    147139               CALL trd_icp( ztswu, ztswv, jpicpd_swf, ctype )                               ! wind stress forcing term 
    148                CALL trd_icp( ztbfu, ztbfv, jpicpd_bfr, ctype )                               ! bottom friction term 
     140               ! bottom friction contribution now handled explicitly 
     141            CASE ( jpdyn_trd_bfr )   ;   CALL trd_icp( ptrdx, ptrdy, jpicpd_bfr, ctype )     ! bottom friction term 
    149142            END SELECT 
    150143            ! 

  • trunk/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r1601 r1662  
    44   !! Ocean physics: Bottom friction 
    55   !!====================================================================== 
    6    !! History :  8.0  ! 1997-06  (G. Madec, A.-M. Treguier)  Original code 
     6   !! History :  OPA  ! 1997-06  (G. Madec, A.-M. Treguier)  Original code 
    77   !!   NEMO     1.0  ! 2002-06  (G. Madec)  F90: Free form and module 
     8   !!            3.2  ! 2001-09  (A.C.Coward)  Correction to include barotropic contribution 
    89   !!---------------------------------------------------------------------- 
    910 
    … …  
    1112   !!   zdf_bfr      : update momentum Kz at the ocean bottom due to the type of bottom friction chosen 
    1213   !!   zdf_bfr_init : read in namelist and control the bottom friction parameters. 
     14   !!   zdf_bfr_2d   : read in namelist and control the bottom friction 
     15   !!                  parameters. 
    1316   !!---------------------------------------------------------------------- 
    1417   USE oce             ! ocean dynamics and tracers variables 
    … …  
    2932   REAL(wp) ::   rn_bfri2 = 1.0e-3_wp   ! bottom drag coefficient (non linear case) 
    3033   REAL(wp) ::   rn_bfeb2 = 2.5e-3_wp   ! background bottom turbulent kinetic energy  [m2/s2] 
     34   REAL(wp) ::   rn_bfrien = 30_wp      ! local factor to enhance coefficient bfri 
     35   REAL(wp), DIMENSION(jpi,jpj) :: & 
     36                 bfrcoef2d = 1.e-3_wp   ! 2D bottom drag coefficient 
     37   LOGICAL  ::   ln_bfr2d = .false.     ! logical switch for 2D enhancement 
     38   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: &  
     39      &          bfrua , bfrva          ! Bottom friction coefficients set in zdfbfr 
    3140 
    3241   !! * Substitutions 
    … …  
    4756      !!              Kz at the ocean bottom. 
    4857      !! 
    49       !! ** Method  :   Update the value of avmu and avmv at the ocean bottom  
    50       !!       level following the chosen friction type (no-slip, free-slip,  
    51       !!       linear, or quadratic) 
     58      !! ** Method  :   Calculate and store part of the momentum trend due     
     59      !!              to bottom friction following the chosen friction type  
     60      !!              (free-slip, linear, or quadratic). The component 
     61      !!              calculated here is multiplied by the bottom velocity in 
     62      !!              dyn_bfr to provide the trend term. 
     63      !! 
     64      !! ** Action  :   bfrua , bfrva   bottom friction coefficients 
    5265      !!---------------------------------------------------------------------- 
    5366      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    5467      !! 
    55       INTEGER  ::   ji, jj   ! dummy loop indexes 
    56       INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1   ! temporary integers 
    57       REAL(wp) ::   zvu, zuv, zecu, zecv         ! temporary scalars 
    58       !!---------------------------------------------------------------------- 
    59  
    60       IF( kt == nit000 )   CALL zdf_bfr_init 
    61  
    62       !                               ! -------------------------------------- 
    63       SELECT CASE (nn_bfr)            ! Compute avmu, avmv at the ocean bottom 
    64       !                               ! -------------------------------------- 
    65       ! 
    66       CASE( 0 )                            !==  no-slip boundary condition  ==! 
     68      INTEGER  ::   ji, jj         ! dummy loop indexes 
     69      INTEGER  ::   ikbu, ikbum1   ! temporary integers 
     70      INTEGER  ::   ikbv, ikbvm1   ! temporary integers 
     71      REAL(wp) ::   zvu, zuv, zecu, zecv   ! temporary scalars 
     72      !!---------------------------------------------------------------------- 
     73 
     74 
     75      IF( kt == nit000 )   CALL zdf_bfr_init   ! initialisation 
     76 
     77      IF( nn_bfr == 2 ) THEN                 ! quadratic botton friction 
     78         ! Calculate and store the quadratic bottom friction terms bfrua and bfrva 
     79         ! where bfrUa = C_d*SQRT(u_bot^2 + v_bot^2 + e_b) {U=[u,v]} 
     80         ! from these the trend due to bottom friction:  -F_h/e3U  can be calculated 
     81         ! where -F_h/e3U_bot = bfrUa*Ub/e3U_bot {U=[u,v]} 
     82         ! 
    6783# if defined key_vectopt_loop 
    6884         DO jj = 1, 1 
     85!CDIR NOVERRCHK 
    6986            DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    7087# else 
     88!CDIR NOVERRCHK 
    7189         DO jj = 2, jpjm1 
     90!CDIR NOVERRCHK 
    7291            DO ji = 2, jpim1 
    7392# endif 
    … …  
    7695               ikbum1 = MAX( ikbu-1, 1 ) 
    7796               ikbvm1 = MAX( ikbv-1, 1 ) 
    78                avmu(ji,jj,ikbu) = 2. * avmu(ji,jj,ikbum1) 
    79                avmv(ji,jj,ikbv) = 2. * avmv(ji,jj,ikbvm1) 
     97               ! 
     98               zvu  = 0.25 * (  vn(ji,jj  ,ikbum1) + vn(ji+1,jj  ,ikbum1)     & 
     99                  &           + vn(ji,jj-1,ikbum1) + vn(ji+1,jj-1,ikbum1)  ) 
     100               zuv  = 0.25 * (  un(ji,jj  ,ikbvm1) + un(ji-1,jj  ,ikbvm1)     & 
     101                  &           + un(ji,jj+1,ikbvm1) + un(ji-1,jj+1,ikbvm1)  ) 
     102               ! 
     103               zecu = SQRT(  un(ji,jj,ikbum1) * un(ji,jj,ikbum1) + zvu*zvu + rn_bfeb2  ) 
     104               zecv = SQRT(  vn(ji,jj,ikbvm1) * vn(ji,jj,ikbvm1) + zuv*zuv + rn_bfeb2  ) 
     105               ! 
     106               bfrua(ji,jj) = - bfrcoef2d(ji,jj) * zecu  
     107               bfrva(ji,jj) = - bfrcoef2d(ji,jj) * zecv 
    80108            END DO 
    81109         END DO 
    82  
    83       CASE( 1 )                            !==  linear botton friction  ==! 
    84 # if defined key_vectopt_loop 
    85          DO jj = 1, 1 
    86             DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    87 # else 
    88          DO jj = 2, jpjm1 
    89             DO ji = 2, jpim1 
    90 # endif 
    91                ikbu = MIN( mbathy(ji+1,jj), mbathy(ji,jj) ) 
    92                ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) 
    93                avmu(ji,jj,ikbu) = rn_bfri1 * fse3uw(ji,jj,ikbu) 
    94                avmv(ji,jj,ikbv) = rn_bfri1 * fse3vw(ji,jj,ikbv) 
    95             END DO 
    96          END DO 
    97  
    98       CASE( 2 )                            !==  quadratic botton friction  ==! 
    99 # if defined key_vectopt_loop 
    100          DO jj = 1, 1 
    101 !CDIR NOVERRCHK 
    102             DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    103 # else 
    104 !CDIR NOVERRCHK 
    105          DO jj = 2, jpjm1 
    106 !CDIR NOVERRCHK 
    107             DO ji = 2, jpim1 
    108 # endif 
    109                ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) 
    110                ikbv   = MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) ) 
    111                ikbum1 = MAX( ikbu-1, 1 ) 
    112                ikbvm1 = MAX( ikbv-1, 1 ) 
    113                 
    114                zvu  = 0.25 * (  vn(ji,jj  ,ikbum1) + vn(ji+1,jj  ,ikbum1)     & 
    115                               + vn(ji,jj-1,ikbum1) + vn(ji+1,jj-1,ikbum1)  ) 
    116                 
    117                zuv  = 0.25 * (  un(ji,jj  ,ikbvm1) + un(ji-1,jj  ,ikbvm1)     & 
    118                               + un(ji,jj+1,ikbvm1) + un(ji-1,jj+1,ikbvm1)  ) 
    119                 
    120                zecu = SQRT(  un(ji,jj,ikbum1) * un(ji,jj,ikbum1) + zvu*zvu + rn_bfeb2  ) 
    121                zecv = SQRT(  vn(ji,jj,ikbvm1) * vn(ji,jj,ikbvm1) + zuv*zuv + rn_bfeb2  ) 
    122                 
    123                avmu(ji,jj,ikbu) = rn_bfri2 * zecu * fse3uw(ji,jj,ikbu) 
    124                avmv(ji,jj,ikbv) = rn_bfri2 * zecv * fse3vw(ji,jj,ikbv) 
    125             END DO 
    126          END DO 
    127  
    128       CASE( 3 )                            !==  free-slip boundary condition  ==! 
    129 # if defined key_vectopt_loop 
    130          DO jj = 1, 1 
    131             DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    132 # else 
    133          DO jj = 2, jpjm1 
    134             DO ji = 2, jpim1 
    135 # endif 
    136                ikbu = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) 
    137                ikbv = MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) ) 
    138                avmu(ji,jj,ikbu) = 0.e0 
    139                avmv(ji,jj,ikbv) = 0.e0 
    140             END DO 
    141          END DO 
    142          ! 
    143       END SELECT 
    144       CALL lbc_lnk( avmu, 'U', 1. )   ;    CALL lbc_lnk( avmv, 'V', 1. )   ! Lateral boundary condition   (unchanged sign) 
    145  
    146       IF(ln_ctl)   CALL prt_ctl( tab3d_1=avmu, clinfo1=' bfr  - u: ', mask1=umask,             & 
    147          &                       tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask,ovlap=1, kdim=jpk ) 
     110         ! 
     111      ENDIF 
     112      ! 
     113      CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition 
     114 
     115      IF(ln_ctl)   CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr  - u: ', mask1=umask,        & 
     116         &                       tab2d_2=bfrva, clinfo2=       ' v: ', mask2=vmask,ovlap=1 ) 
    148117      ! 
    149118   END SUBROUTINE zdf_bfr 
    … …  
    157126      !! 
    158127      !! ** Method  :   Read the nammbf namelist and check their consistency 
    159       !!---------------------------------------------------------------------- 
    160       NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfeb2 
    161       !!---------------------------------------------------------------------- 
    162  
    163       REWIND( numnam )              ! Read Namelist nambfr : bottom momentum boundary condition 
    164       READ  ( numnam, nambfr ) 
    165  
    166       IF(lwp) WRITE(numout,*)       ! Parameter print 
     128      !!              called at the first timestep (nit000) 
     129      !!---------------------------------------------------------------------- 
     130      USE iom   ! I/O module for ehanced bottom friction file 
     131      !! 
     132      NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfeb2, ln_bfr2d, rn_bfrien 
     133      INTEGER ::   inum       ! logical unit for enhanced bottom friction file 
     134      INTEGER ::   ji, jj     ! dummy loop indexes 
     135      INTEGER ::   ikbu, ikbv ! temporary integers 
     136      INTEGER ::   ictu, ictv ! temporary integers 
     137      REAL(wp) ::  rminbfr, rmaxbfr 
     138                              ! temporary reals 
     139      !!---------------------------------------------------------------------- 
     140 
     141      REWIND ( numnam )               !* Read Namelist nam_bfr : bottom momentum boundary condition 
     142      READ   ( numnam, nambfr ) 
     143 
     144 
     145      !                               !* Parameter control and print 
     146      IF(lwp) WRITE(numout,*) 
    167147      IF(lwp) WRITE(numout,*) 'zdf_bfr : momentum bottom friction' 
    168148      IF(lwp) WRITE(numout,*) '~~~~~~~' 
    169       IF(lwp) WRITE(numout,*) '   Namelist nambfr : set bottom friction parameters' 
    170  
    171       SELECT CASE (nn_bfr)          ! Parameter control 
    172       ! 
     149      IF(lwp) WRITE(numout,*) '          Namelist nam_bfr : set bottom friction parameters' 
     150 
     151      SELECT CASE (nn_bfr) 
     152 
    173153      CASE( 0 ) 
    174          IF(lwp) WRITE(numout,*) '      no-slip ' 
     154         IF(lwp) WRITE(numout,*) '            free-slip ' 
     155         ! 
     156         bfrua(:,:) = 0.e0 
     157         bfrva(:,:) = 0.e0 
    175158         ! 
    176159      CASE( 1 ) 
    177          IF(lwp) WRITE(numout,*) '      linear botton friction      nn_bfr    = ', nn_bfr 
    178          IF(lwp) WRITE(numout,*) '      friction coef.              rn_bfri1  = ', rn_bfri1 
     160         IF(lwp) WRITE(numout,*) '            linear botton friction' 
     161         IF(lwp) WRITE(numout,*) '            friction coef.   rn_bfri1  = ', rn_bfri1 
     162         IF(lwp)  THEN 
     163            IF( ln_bfr2d ) THEN 
     164                 WRITE(numout,*) '            read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d 
     165                 WRITE(numout,*) '            coef rn_bfri2 enhancement factor   rn_bfrien  = ',rn_bfrien 
     166            ENDIF 
     167         ENDIF 
     168         ! 
     169         bfrcoef2d(:,:) = rn_bfri1  ! initialize bfrcoef2d to the namelist variable 
     170 
     171         IF (ln_bfr2d) THEN  
     172            ! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 
     173            CALL iom_open('bfr_coef.nc',inum) 
     174            CALL iom_get (inum, jpdom_data, 'bfr_coef',bfrcoef2d,1) ! bfrcoef2d is used as tmp array 
     175            CALL iom_close(inum) 
     176            bfrcoef2d(:,:)= rn_bfri1*(1+ rn_bfrien*bfrcoef2d(:,:)) 
     177         ENDIF 
     178         bfrua(:,:) = - bfrcoef2d(:,:) 
     179         bfrva(:,:) = - bfrcoef2d(:,:) 
    179180         ! 
    180181      CASE( 2 ) 
    181          IF(lwp) WRITE(numout,*) '      quadratic botton friction   nn_bfr    = ', nn_bfr 
    182          IF(lwp) WRITE(numout,*) '      friction coef.              rn_bfri2  = ', rn_bfri2 
    183          IF(lwp) WRITE(numout,*) '      background KE               rn_bfeb2  = ', rn_bfeb2 
    184          ! 
    185       CASE( 3 ) 
    186          IF(lwp) WRITE(numout,*) '      free-slip ' 
     182         IF(lwp) WRITE(numout,*) '            quadratic botton friction' 
     183         IF(lwp) WRITE(numout,*) '            friction coef.   rn_bfri2  = ', rn_bfri2 
     184         IF(lwp) WRITE(numout,*) '            background tke   rn_bfeb2  = ', rn_bfeb2 
     185         IF(lwp)  THEN 
     186            IF( ln_bfr2d ) THEN 
     187                 WRITE(numout,*) '            read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d 
     188                 WRITE(numout,*) '            coef rn_bfri2 enhancement factor   rn_bfrien  = ',rn_bfrien 
     189            ENDIF 
     190         ENDIF 
     191         bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable 
     192 
     193         IF (ln_bfr2d) THEN  
     194            ! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 
     195            CALL iom_open('bfr_coef.nc',inum) 
     196            CALL iom_get (inum, jpdom_data, 'bfr_coef',bfrcoef2d,1) ! bfrcoef2d is used as tmp array 
     197            CALL iom_close(inum) 
     198            bfrcoef2d(:,:)= rn_bfri2*(1+ rn_bfrien*bfrcoef2d(:,:)) 
     199         ENDIF 
     200 
    187201         ! 
    188202      CASE DEFAULT 
    189          WRITE(ctmp1,*) 'bad flag value for nn_bfr = ', nn_bfr 
     203         WRITE(ctmp1,*) '         bad flag value for nn_bfr = ', nn_bfr 
    190204         CALL ctl_stop( ctmp1 ) 
    191205         ! 
    192206      END SELECT 
     207      ! 
     208      ! Basic stability check on bottom friction coefficient 
     209      ! 
     210      ictu = 0               ! counter for stability criterion breaches at U-pts 
     211      ictv = 0               ! counter for stability criterion breaches at V-pts 
     212      rminbfr =  1.e10_wp    ! initialise tracker for minimum of bottom friction coefficient 
     213      rmaxbfr = -1.e10_wp    ! initialise tracker for maximum of bottom friction coefficient 
     214      ! 
     215#  if defined key_vectopt_loop 
     216      DO jj = 1, 1 
     217!CDIR NOVERRCHK 
     218         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
     219#  else 
     220!CDIR NOVERRCHK 
     221      DO jj = 2, jpjm1 
     222!CDIR NOVERRCHK 
     223         DO ji = 2, jpim1 
     224#  endif 
     225             ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) 
     226             ikbv   = MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) ) 
     227             IF ( bfrcoef2d(ji,jj).gt.2.0*fse3u(ji,jj,ikbu)/rdt ) then 
     228                IF ( ln_ctl ) THEN 
     229                   write(numout,*) 'BFR ',narea,nimpp+ji,njmpp+jj,ikbu 
     230                   write(numout,*) 'BFR ',bfrcoef2d(ji,jj),2.0*fse3u(ji,jj,ikbu)/rdt 
     231                ENDIF 
     232                ictu = ictu + 1 
     233             ENDIF 
     234             IF ( bfrcoef2d(ji,jj).gt.2.0*fse3v(ji,jj,ikbv)/rdt ) then 
     235                 IF ( ln_ctl ) THEN 
     236                     write(numout,*) 'BFR ',narea,nimpp+ji,njmpp+jj,ikbv 
     237                     write(numout,*) 'BFR ',bfrcoef2d(ji,jj),2.0*fse3v(ji,jj,ikbv)/rdt 
     238                 ENDIF 
     239                 ictv = ictv + 1 
     240             ENDIF 
     241             bfrcoef2d(ji,jj) = MIN(2.0*fse3u(ji,jj,ikbu)/rdt, bfrcoef2d(ji,jj)) 
     242             bfrcoef2d(ji,jj) = MIN(2.0*fse3v(ji,jj,ikbv)/rdt, bfrcoef2d(ji,jj)) 
     243             rminbfr = MIN(rminbfr, bfrcoef2d(ji,jj)) 
     244             rmaxbfr = MAX(rmaxbfr, bfrcoef2d(ji,jj)) 
     245         END DO 
     246      END DO 
     247      IF ( lk_mpp ) THEN 
     248         CALL mpp_sum( ictu ) 
     249         CALL mpp_sum( ictv ) 
     250         CALL mpp_min( rminbfr ) 
     251         CALL mpp_max( rmaxbfr ) 
     252      ENDIF 
     253      IF ( lwp .AND. ictu + ictv .GT. 0 ) THEN 
     254         WRITE(numout,*) ' Bottom friction stability check failed at ', ictu, ' U-points ' 
     255         WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points ' 
     256         WRITE(numout,*) ' Bottom friction coefficient reset where necessary' 
     257         WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', rminbfr, ' to ', rmaxbfr 
     258      ENDIF 
    193259      ! 
    194260   END SUBROUTINE zdf_bfr_init 

  • trunk/NEMO/OPA_SRC/ZDF/zdftke.F90

    r1617 r1662  
    4848   USE in_out_manager ! I/O manager 
    4949   USE iom            ! I/O manager library 
     50   USE zdfbfr          ! bottom friction 
    5051 
    5152   IMPLICIT NONE 
    … …  
    182183      REAL(wp) ::   zus   , zwlc  , zind      !    -         - 
    183184      REAL(wp) ::   zzd_up, zzd_lw            !    -         - 
     185      INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1      ! temporary scalar 
     186      INTEGER  ::   ikbt, ikbumm1, ikbvmm1          ! temporary scalar 
     187      REAL(wp) ::   zebot                           ! temporary scalars 
    184188      INTEGER , DIMENSION(jpi,jpj)     ::   imlc    ! 2D workspace 
    185189      REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc    !  -      - 
    … …  
    209213      !                     !  Bottom boundary condition on tke 
    210214      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    211 !!gm to be added soon 
     215      !                     en(bot)   = 0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin) 
     216!CDIR NOVERRCHK 
     217      DO jj = 2, jpjm1 
     218!CDIR NOVERRCHK 
     219         DO ji = fs_2, fs_jpim1   ! vector opt. 
     220            ikbu = MIN( mbathy(ji+1,jj), mbathy(ji,jj) ) 
     221            ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) 
     222            ikbum1 = MAX( ikbu-1, 1 ) 
     223            ikbvm1 = MAX( ikbv-1, 1 ) 
     224            ikbu = MIN( mbathy(ji,jj), mbathy(ji-1,jj) ) 
     225            ikbv = MIN( mbathy(ji,jj), mbathy(ji,jj-1) ) 
     226            ikbumm1 = MAX( ikbu-1, 1 ) 
     227            ikbvmm1 = MAX( ikbv-1, 1 ) 
     228            ikbt = MAX( mbathy(ji,jj), 1 ) 
     229            ztx2 = bfrua(ji-1,jj) * fse3u(ji-1,jj  ,ikbumm1) + & 
     230                   bfrua(ji  ,jj) * fse3u(ji  ,jj  ,ikbum1 ) 
     231            zty2 = bfrva(ji,jj  ) * fse3v(ji  ,jj  ,ikbvm1) + & 
     232                   bfrva(ji,jj-1) * fse3v(ji  ,jj-1,ikbvmm1 ) 
     233            zebot = 0.25_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) 
     234            en (ji,jj,ikbt) = MAX( zebot, rn_emin ) * tmask(ji,jj,1) 
     235         END DO 
     236      END DO 
    212237      ! 
    213238      ! 

  • trunk/NEMO/OPA_SRC/step.F90

    r1635 r1662  
    6464 
    6565   USE dynadv          ! advection                        (dyn_adv routine) 
     66   USE dynbfr          ! Bottom friction terms            (dyn_bfr routine) 
    6667   USE dynvor          ! vorticity term                   (dyn_vor routine) 
    6768   USE dynhpg          ! hydrostatic pressure grad.       (dyn_hpg routine) 
    … …  
    197198      ! 
    198199      !  VERTICAL PHYSICS    
     200                         CALL zdf_bfr( kstp )         ! bottom friction 
    199201      !                                               ! Vertical eddy viscosity and diffusivity coefficients 
    200202      IF( lk_zdfric     )   CALL zdf_ric    ( kstp )  ! Richardson number dependent Kz 
    … …  
    216218      IF( lk_zdfddm .AND. .NOT. lk_zdfkpp )   & 
    217219         &               CALL zdf_ddm( kstp )         ! double diffusive mixing 
    218  
    219                          CALL zdf_bfr( kstp )         ! bottom friction 
    220  
    221220                         CALL zdf_mxl( kstp )         ! mixed layer depth 
    222221 
    … …  
    306305#endif 
    307306                               CALL dyn_hpg( kstp )         ! horizontal gradient of Hydrostatic pressure 
     307                               CALL dyn_bfr( kstp )           ! bottom friction    
    308308                               CALL dyn_zdf( kstp )         ! vertical diffusion 
    309309                               indic=0 
Note: See TracChangeset for help on using the changeset viewer.