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

Changeset 1708


Ignore:
Timestamp:
2009-11-04T14:24:34+01:00 (14 years ago)
Author:
rblod
Message:

Stability for explicit bottom friction, see ticket #588

Location:
trunk/NEMO/OPA_SRC
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r1671 r1708  
    44   !! Ocean dynamics :  bottom friction component of the momentum mixing trend 
    55   !!============================================================================== 
    6    !! History :  9.0  !  08-11  (A. C. Coward)  Original code 
     6   !! History :  9.0  !  2008-11  (A. C. Coward)  Original code 
    77   !!---------------------------------------------------------------------- 
    88 
    … …  
    4545      !! ** Action  :   (ua,va)   momentum trend increased by bottom friction trend 
    4646      !!--------------------------------------------------------------------- 
     47      USE oce, ONLY :   ztrdu => ta   ! use ta as 3D workspace 
     48      USE oce, ONLY :   ztrdv => sa   ! use sa as 3D workspace 
     49      !! 
    4750      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
    4851      !!  
    49       INTEGER ::   ji, jj       ! dummy loop indexes 
    50       INTEGER ::   ikbu, ikbv   ! temporary integers 
    51       INTEGER ::   ikbum1, ikbvm1 
    52 !! + Enforce stability 
    53 !     REAL(wp):: frdt           ! temporary real 
    54 !! - Enforce stability 
    55       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztrdu, ztrdv   ! 3D workspace 
     52      INTEGER  ::   ji, jj          ! dummy loop indexes 
     53      INTEGER  ::   ikbu  , ikbv    ! temporary integers 
     54      INTEGER  ::   ikbum1, ikbvm1  !    -          - 
     55      REAL(wp) ::   zinv, zbfru, zbfrv   ! temporary scalar 
    5656      !!--------------------------------------------------------------------- 
    5757      ! 
    58 !! + Enforce stability 
    59 !     frdt = 0.125/rdt 
    60 !! - Enforce stability 
    61       IF( l_trddyn )   THEN                      ! temporary save of ta and sa trends 
     58      zinv = -1. / ( 2.*rdt ) 
     59 
     60      IF( l_trddyn )   THEN                      ! temporary save of ua and va trends 
    6261         ztrdu(:,:,:) = ua(:,:,:) 
    6362         ztrdv(:,:,:) = va(:,:,:) 
    6463      ENDIF 
     64 
    6565# if defined key_vectopt_loop 
    6666      DO jj = 1, 1 
    … …  
    7575            ikbvm1 = MAX( ikbv-1, 1 ) 
    7676            ! 
    77 !! + Enforce stability 
    7877            ! Apply stability criteria 
    79 !            bfrua(ji,jj) = MIN(bfrua(ji,jj), fse3u(ji,jj,ikbu)*frdt) 
    80 !            bfrva(ji,jj) = MIN(bfrva(ji,jj), fse3v(ji,jj,ikbv)*frdt) 
     78            zbfru = MAX( bfrua(ji,jj), fse3u(ji,jj,ikbu)*zinv ) 
     79            zbfrv = MAX( bfrva(ji,jj), fse3v(ji,jj,ikbv)*zinv ) 
    8180            ! 
    82 !! - Enforce stability 
    83             ua(ji,jj,ikbum1) = ua(ji,jj,ikbum1) + bfrua(ji,jj) * ub(ji,jj,ikbum1) / fse3u(ji,jj,ikbu) 
    84             va(ji,jj,ikbvm1) = va(ji,jj,ikbvm1) + bfrva(ji,jj) * vb(ji,jj,ikbvm1) / fse3v(ji,jj,ikbu) 
     81            ua(ji,jj,ikbum1) = ua(ji,jj,ikbum1) + zbfru * ub(ji,jj,ikbum1) / fse3u(ji,jj,ikbu) 
     82            va(ji,jj,ikbvm1) = va(ji,jj,ikbvm1) + zbfrv * vb(ji,jj,ikbvm1) / fse3v(ji,jj,ikbv) 
     83            ! 
    8584         END DO 
    8685      END DO 
     86 
    8787      ! 
    8888      IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics 

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

    r1694 r1708  
    9494      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    9595      INTEGER  ::   icycle           ! temporary scalar 
     96      INTEGER  ::   ikbu, ikbv        ! temporary scalar 
    9697 
    9798      REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b     ! temporary scalars 
    … …  
    102103      !! 
    103104      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv, zsshb_e 
     105      !! 
     106      REAL(wp), DIMENSION(jpi,jpj) ::   zbfru  , zbfrv     ! 2D workspace 
    104107      !! 
    105108      REAL(wp), DIMENSION(jpi,jpj) ::   zsshun_e, zsshvn_e   ! 2D workspace 
    … …  
    119122         IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ',  2*nn_baro 
    120123         ! 
    121          CALL ts_rst( nit000, 'READ' )   ! read or initialize the following fields: 
    122          !                               ! un_b, vn_b   
     124         CALL ts_rst( nit000, 'READ' )   ! read or initialize the following fields: un_b, vn_b   
    123125         ! 
    124126         ua_e  (:,:) = un_b (:,:) 
    … …  
    131133            ftne(1,:) = 0.e0   ;   ftnw(1,:) = 0.e0   ;   ftse(1,:) = 0.e0   ;   ftsw(1,:) = 0.e0 
    132134            DO jj = 2, jpj 
    133                DO ji = 2, jpi   ! NO vector opt. 
     135               DO ji = fs_2, jpi   ! vector opt. 
    134136                  ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3. 
    135137                  ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3. 
    … …  
    252254      END DO 
    253255 
    254       !                                            ! Remove barotropic contribution of bottom friction  
    255                                                    ! from the barotropic transport trend 
     256                     
     257      !                                             ! Remove barotropic contribution of bottom friction  
     258      !                                             ! from the barotropic transport trend 
     259      zcoef = -1. / z2dt_b 
     260# if defined key_vectopt_loop 
     261      DO jj = 1, 1 
     262         DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
     263# else 
     264      DO jj = 2, jpjm1 
     265         DO ji = 2, jpim1 
     266# endif 
     267            ikbu = MIN( mbathy(ji+1,jj), mbathy(ji,jj) ) 
     268            ikbv = MIN( mbathy(ji,jj+1), mbathy(ji,jj) ) 
     269            ! 
     270            ! Apply stability criteria for bottom friction 
     271            !RBbug for vvl and external mode we may need to 
     272            ! use varying fse3 
     273            zbfru  (ji,jj) = MAX( bfrua(ji,jj), fse3u(ji,jj,ikbu)*zcoef ) 
     274            zbfrv  (ji,jj) = MAX( bfrva(ji,jj), fse3v(ji,jj,ikbv)*zcoef ) 
     275         END DO 
     276      END DO 
     277 
    256278      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 
     279         DO jj = 2, jpjm1 
     280            DO ji = fs_2, fs_jpim1   ! vector opt. 
     281               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * zub(ji,jj)   & 
     282                  &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) ) 
     283               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * zvb(ji,jj)   & 
     284                  &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) ) 
     285            END DO 
     286         END DO 
    263287      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 
     288         DO jj = 2, jpjm1 
     289            DO ji = fs_2, fs_jpim1   ! vector opt. 
     290               zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * zub(ji,jj) * hur(ji,jj) 
     291               zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * zvb(ji,jj) * hvr(ji,jj) 
     292            END DO 
     293         END DO 
    270294      ENDIF 
    271295 

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

    r1601 r1708  
    117117                  umo(ktrd) = umo(ktrd) + ptrd2dx(ji,jj) * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1) 
    118118                  vmo(ktrd) = vmo(ktrd) + ptrd2dy(ji,jj) * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1) 
    119                END DO 
    120             END DO 
    121             ! 
    122          CASE( jpdyn_trd_bfr )         ! bottom friction fluxes 
    123             DO jj = 1, jpj 
    124                DO ji = 1, jpi 
    125                   umo(ktrd) = umo(ktrd) + ptrd2dx(ji,jj) 
    126                   vmo(ktrd) = vmo(ktrd) + ptrd2dy(ji,jj) 
    127119               END DO 
    128120            END DO 
    … …  
    520512            WRITE (numout,9530) hke(jpicpd_swf) / tvolt 
    521513            WRITE (numout,9531) hke(jpicpd_dat) / tvolt 
    522             WRITE (numout,9532) 
    523             WRITE (numout,9533)   & 
     514            WRITE (numout,9532) hke(jpicpd_bfr) / tvolt 
     515            WRITE (numout,9533) 
     516            WRITE (numout,9534)   & 
    524517            &     (  hke(jpicpd_hpg) + hke(jpicpd_keg) + hke(jpicpd_rvo) + hke(jpicpd_pvo) + hke(jpicpd_ldf)   & 
    525518            &      + hke(jpicpd_had) + hke(jpicpd_zad) + hke(jpicpd_zdf) + hke(jpicpd_spg) + hke(jpicpd_dat)   & 
    526             &      + hke(jpicpd_swf) ) / tvolt 
     519            &      + hke(jpicpd_swf) + hke(jpicpd_bfr) ) / tvolt 
    527520         ENDIF 
    528521 
    … …  
    539532 9530    FORMAT(' surface wind forcing      u2= ', e20.13) 
    540533 9531    FORMAT(' dampimg term              u2= ', e20.13) 
    541  9532    FORMAT(' --------------------------------------------------') 
    542  9533    FORMAT(' total trend               u2= ', e20.13) 
     534 9532    FORMAT(' bottom flux               u2= ', e20.13) 
     535 9533    FORMAT(' --------------------------------------------------') 
     536 9534    FORMAT(' total trend               u2= ', e20.13) 
    543537 
    544538         IF(lwp) THEN 

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

    r1662 r1708  
    5757      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::   ptrdy ! Salinity    or V trend 
    5858      !! 
    59       INTEGER ::   ji, ikbu, ikbum1 
    60       INTEGER ::   jj, ikbv, ikbvm1 
     59      INTEGER ::   ji, jj 
    6160      CHARACTER(len=3) ::   ccpas                                  ! number of passage 
    62       REAL(wp) ::   zua, zva                                       ! scalars 
    6361      REAL(wp), DIMENSION(jpi,jpj) ::   ztswu, ztswv               ! 2D workspace 
    6462      REAL(wp), DIMENSION(jpi,jpj) ::   ztbfu, ztbfv               ! 2D workspace 
    … …  
    131129                     ptrdx(ji,jj,1     ) = ptrdx(ji,jj,1     ) - ztswu(ji,jj) 
    132130                     ptrdy(ji,jj,1     ) = ptrdy(ji,jj,1     ) - ztswv(ji,jj) 
    133                      ptrdx(ji,jj,ikbum1) = ptrdx(ji,jj,ikbum1) - ztbfu(ji,jj) 
    134                      ptrdy(ji,jj,ikbvm1) = ptrdy(ji,jj,ikbvm1) - ztbfv(ji,jj) 
    135131                  END DO 
    136132               END DO 
    … …  
    172168                  ztswu(ji,jj) = utau(ji,jj) / ( fse3u(ji,jj,1)*rau0 ) 
    173169                  ztswv(ji,jj) = vtau(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) 
    174                   ! save bottom friction momentum fluxes 
    175                   ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) 
    176                   ikbv   = MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) ) 
    177                   ikbum1 = MAX( ikbu-1, 1 ) 
    178                   ikbvm1 = MAX( ikbv-1, 1 ) 
    179                   zua = ua(ji,jj,ikbum1) * r2dt + ub(ji,jj,ikbum1) 
    180                   zva = va(ji,jj,ikbvm1) * r2dt + vb(ji,jj,ikbvm1) 
    181                   ztbfu(ji,jj) = - avmu(ji,jj,ikbu) * zua / ( fse3u(ji,jj,ikbum1)*fse3uw(ji,jj,ikbu) ) 
    182                   ztbfv(ji,jj) = - avmv(ji,jj,ikbv) * zva / ( fse3v(ji,jj,ikbvm1)*fse3vw(ji,jj,ikbv) ) 
    183170                  ! 
    184171                  ptrdx(ji,jj,1     ) = ptrdx(ji,jj,1     ) - ztswu(ji,jj) 
    185                   ptrdx(ji,jj,ikbum1) = ptrdx(ji,jj,ikbum1) - ztbfu(ji,jj) 
    186172                  ptrdy(ji,jj,1     ) = ptrdy(ji,jj,1     ) - ztswv(ji,jj) 
    187                   ptrdy(ji,jj,ikbvm1) = ptrdy(ji,jj,ikbvm1) - ztbfv(ji,jj) 
    188173               END DO 
    189174            END DO 
    … …  
    191176            CALL trd_vor_zint( ptrdx, ptrdy, jpvor_zdf )    
    192177            CALL trd_vor_zint( ztswu, ztswv, jpvor_swf )                               ! Wind stress forcing term 
    193             CALL trd_vor_zint( ztbfu, ztbfv, jpvor_bfr )                               ! Bottom friction term 
     178         CASE ( jpdyn_trd_bfr ) 
     179            CALL trd_vor_zint( ptrdx, ptrdy, jpvor_bfr )                               ! Bottom friction term 
    194180         END SELECT 
    195181         ! 

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

    r1671 r1708  
    66   !! 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 
     8   !!            3.2  ! 2009-09  (A.C.Coward)  Correction to include barotropic contribution 
    99   !!---------------------------------------------------------------------- 
    1010 
    … …  
    2727 
    2828   PUBLIC   zdf_bfr    ! called by step.F90 
    29  
    30    !                                   !!* Namelist nambfr: bottom friction namelist * 
    31    INTEGER  ::   nn_bfr   = 0           ! = 0/1/2/3 type of bottom friction  
    32    REAL(wp) ::   rn_bfri1 = 4.0e-4_wp   ! bottom drag coefficient (linear case)  
    33    REAL(wp) ::   rn_bfri2 = 1.0e-3_wp   ! bottom drag coefficient (non linear case) 
    34    REAL(wp) ::   rn_bfeb2 = 2.5e-3_wp   ! background bottom turbulent kinetic energy  [m2/s2] 
    35    REAL(wp) ::   rn_bfrien = 30_wp      ! local factor to enhance coefficient bfri 
    36    REAL(wp), DIMENSION(jpi,jpj) :: & 
    37                  bfrcoef2d = 1.e-3_wp   ! 2D bottom drag coefficient 
    38    LOGICAL  ::   ln_bfr2d = .false.     ! logical switch for 2D enhancement 
    39    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: &  
    40       &          bfrua , bfrva          ! Bottom friction coefficients set in zdfbfr 
     29    
     30   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   bfrua , bfrva   !: Bottom friction coefficients set in zdfbfr 
     31 
     32   !                                    !!* Namelist nambfr: bottom friction namelist * 
     33   INTEGER  ::   nn_bfr    = 0           ! = 0/1/2/3 type of bottom friction  
     34   REAL(wp) ::   rn_bfri1  = 4.0e-4_wp   ! bottom drag coefficient (linear case)  
     35   REAL(wp) ::   rn_bfri2  = 1.0e-3_wp   ! bottom drag coefficient (non linear case) 
     36   REAL(wp) ::   rn_bfeb2  = 2.5e-3_wp   ! background bottom turbulent kinetic energy  [m2/s2] 
     37   REAL(wp) ::   rn_bfrien = 30_wp       ! local factor to enhance coefficient bfri 
     38   LOGICAL  ::   ln_bfr2d  = .false.     ! logical switch for 2D enhancement 
     39    
     40   REAL(wp), DIMENSION(jpi,jpj) ::   bfrcoef2d = 1.e-3_wp   ! 2D bottom drag coefficient 
    4141 
    4242   !! * Substitutions 
    … …  
    5454      !!                   ***  ROUTINE zdf_bfr  *** 
    5555      !!                  
    56       !! ** Purpose :   Applied the bottom friction through a specification of  
    57       !!              Kz at the ocean bottom. 
     56      !! ** Purpose :   compute the bottom friction coefficient. 
    5857      !! 
    5958      !! ** Method  :   Calculate and store part of the momentum trend due     
    … …  
    6261      !!              calculated here is multiplied by the bottom velocity in 
    6362      !!              dyn_bfr to provide the trend term. 
     63      !!                The coefficients are updated at each time step only 
     64      !!              in the quadratic case. 
    6465      !! 
    6566      !! ** Action  :   bfrua , bfrva   bottom friction coefficients 
    … …  
    6768      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index 
    6869      !! 
    69       INTEGER  ::   ji, jj         ! dummy loop indexes 
     70      INTEGER  ::   ji, jj         ! dummy loop indices 
    7071      INTEGER  ::   ikbu, ikbum1   ! temporary integers 
    71       INTEGER  ::   ikbv, ikbvm1   ! temporary integers 
     72      INTEGER  ::   ikbv, ikbvm1   !    -          - 
    7273      REAL(wp) ::   zvu, zuv, zecu, zecv   ! temporary scalars 
    7374      !!---------------------------------------------------------------------- 
    … …  
    7778 
    7879      IF( nn_bfr == 2 ) THEN                 ! quadratic botton friction 
    79          ! Calculate and store the quadratic bottom friction terms bfrua and bfrva 
     80         ! Calculate and store the quadratic bottom friction coefficient bfrua and bfrva 
    8081         ! where bfrUa = C_d*SQRT(u_bot^2 + v_bot^2 + e_b) {U=[u,v]} 
    8182         ! from these the trend due to bottom friction:  -F_h/e3U  can be calculated 
    … …  
    105106               zecv = SQRT(  vn(ji,jj,ikbvm1) * vn(ji,jj,ikbvm1) + zuv*zuv + rn_bfeb2  ) 
    106107               ! 
    107                bfrua(ji,jj) = - bfrcoef2d(ji,jj) * zecu  
    108                bfrva(ji,jj) = - bfrcoef2d(ji,jj) * zecv 
     108               bfrua(ji,jj) = - 0.5 * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji+1,jj  ) ) * zecu  
     109               bfrva(ji,jj) = - 0.5 * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji  ,jj+1) ) * zecv 
    109110            END DO 
    110111         END DO 
    111112         ! 
     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 ) 
    112117      ENDIF 
    113       ! 
    114       CALL lbc_lnk( bfrua, 'U', 1. )   ;   CALL lbc_lnk( bfrva, 'V', 1. )      ! Lateral boundary condition 
    115  
    116       IF(ln_ctl)   CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr  - u: ', mask1=umask,        & 
    117          &                       tab2d_2=bfrva, clinfo2=       ' v: ', mask2=vmask,ovlap=1 ) 
    118118      ! 
    119119   END SUBROUTINE zdf_bfr 
    … …  
    131131      USE iom   ! I/O module for ehanced bottom friction file 
    132132      !! 
     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   !    -          - 
     137      REAL(wp) ::  zminbfr, zmaxbfr   ! temporary scalars 
     138      REAL(wp) ::  zfru, zfrv         !    -         - 
     139      !! 
    133140      NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfeb2, ln_bfr2d, rn_bfrien 
    134       INTEGER ::   inum       ! logical unit for enhanced bottom friction file 
    135       INTEGER ::   ji, jj     ! dummy loop indexes 
    136       INTEGER ::   ikbu, ikbv ! temporary integers 
    137       INTEGER ::   ictu, ictv ! temporary integers 
    138       REAL(wp) ::  rminbfr, rmaxbfr 
    139                               ! temporary reals 
    140141      !!---------------------------------------------------------------------- 
    141142 
    142143      REWIND ( numnam )               !* Read Namelist nam_bfr : bottom momentum boundary condition 
    143144      READ   ( numnam, nambfr ) 
    144  
    145145 
    146146      !                               !* Parameter control and print 
    … …  
    148148      IF(lwp) WRITE(numout,*) 'zdf_bfr : momentum bottom friction' 
    149149      IF(lwp) WRITE(numout,*) '~~~~~~~' 
    150       IF(lwp) WRITE(numout,*) '          Namelist nam_bfr : set bottom friction parameters' 
     150      IF(lwp) WRITE(numout,*) '   Namelist nam_bfr : set bottom friction parameters' 
    151151 
    152152      SELECT CASE (nn_bfr) 
    153153 
    154154      CASE( 0 ) 
    155          IF(lwp) WRITE(numout,*) '            free-slip ' 
    156          ! 
     155         IF(lwp) WRITE(numout,*) '      free-slip ' 
    157156         bfrua(:,:) = 0.e0 
    158157         bfrva(:,:) = 0.e0 
    159158         ! 
    160159      CASE( 1 ) 
    161          IF(lwp) WRITE(numout,*) '            linear botton friction' 
    162          IF(lwp) WRITE(numout,*) '            friction coef.   rn_bfri1  = ', rn_bfri1 
    163          IF(lwp)  THEN 
    164             IF( ln_bfr2d ) THEN 
    165                  WRITE(numout,*) '            read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d 
    166                  WRITE(numout,*) '            coef rn_bfri2 enhancement factor   rn_bfrien  = ',rn_bfrien 
    167             ENDIF 
     160         IF(lwp) WRITE(numout,*) '      linear botton friction' 
     161         IF(lwp) WRITE(numout,*) '      friction coef.   rn_bfri1  = ', rn_bfri1 
     162         IF( ln_bfr2d ) THEN 
     163            IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d 
     164            IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien 
    168165         ENDIF 
    169166         ! 
    170167         bfrcoef2d(:,:) = rn_bfri1  ! initialize bfrcoef2d to the namelist variable 
    171  
    172          IF (ln_bfr2d) THEN  
     168         ! 
     169         IF(ln_bfr2d) THEN  
    173170            ! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 
    174171            CALL iom_open('bfr_coef.nc',inum) 
    175172            CALL iom_get (inum, jpdom_data, 'bfr_coef',bfrcoef2d,1) ! bfrcoef2d is used as tmp array 
    176173            CALL iom_close(inum) 
    177             bfrcoef2d(:,:)= rn_bfri1*(1+ rn_bfrien*bfrcoef2d(:,:)) 
     174            bfrcoef2d(:,:)= rn_bfri1 * ( 1 + rn_bfrien * bfrcoef2d(:,:) ) 
    178175         ENDIF 
    179176         bfrua(:,:) = - bfrcoef2d(:,:) 
    … …  
    181178         ! 
    182179      CASE( 2 ) 
    183          IF(lwp) WRITE(numout,*) '            quadratic botton friction' 
    184          IF(lwp) WRITE(numout,*) '            friction coef.   rn_bfri2  = ', rn_bfri2 
    185          IF(lwp) WRITE(numout,*) '            background tke   rn_bfeb2  = ', rn_bfeb2 
    186          IF(lwp)  THEN 
    187             IF( ln_bfr2d ) THEN 
    188                  WRITE(numout,*) '            read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d 
    189                  WRITE(numout,*) '            coef rn_bfri2 enhancement factor   rn_bfrien  = ',rn_bfrien 
    190             ENDIF 
     180         IF(lwp) WRITE(numout,*) '      quadratic botton friction' 
     181         IF(lwp) WRITE(numout,*) '      friction coef.   rn_bfri2  = ', rn_bfri2 
     182         IF(lwp) WRITE(numout,*) '      background tke   rn_bfeb2  = ', rn_bfeb2 
     183         IF( ln_bfr2d ) THEN 
     184            IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d 
     185            IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien 
    191186         ENDIF 
    192187         bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable 
    193  
    194          IF (ln_bfr2d) THEN  
     188         ! 
     189         IF(ln_bfr2d) THEN  
    195190            ! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 
    196191            CALL iom_open('bfr_coef.nc',inum) 
    197192            CALL iom_get (inum, jpdom_data, 'bfr_coef',bfrcoef2d,1) ! bfrcoef2d is used as tmp array 
    198193            CALL iom_close(inum) 
    199             bfrcoef2d(:,:)= rn_bfri2*(1+ rn_bfrien*bfrcoef2d(:,:)) 
    200          ENDIF 
    201  
     194            bfrcoef2d(:,:)= rn_bfri2 * ( 1 + rn_bfrien * bfrcoef2d(:,:) ) 
     195         ENDIF 
    202196         ! 
    203197      CASE DEFAULT 
    204          WRITE(ctmp1,*) '         bad flag value for nn_bfr = ', nn_bfr 
     198         IF(lwp) WRITE(ctmp1,*) '         bad flag value for nn_bfr = ', nn_bfr 
    205199         CALL ctl_stop( ctmp1 ) 
    206200         ! 
    … …  
    211205      ictu = 0               ! counter for stability criterion breaches at U-pts 
    212206      ictv = 0               ! counter for stability criterion breaches at V-pts 
    213       rminbfr =  1.e10_wp    ! initialise tracker for minimum of bottom friction coefficient 
    214       rmaxbfr = -1.e10_wp    ! initialise tracker for maximum of bottom friction coefficient 
     207      zminbfr =  1.e10_wp    ! initialise tracker for minimum of bottom friction coefficient 
     208      zmaxbfr = -1.e10_wp    ! initialise tracker for maximum of bottom friction coefficient 
    215209      ! 
    216210#  if defined key_vectopt_loop 
    … …  
    224218         DO ji = 2, jpim1 
    225219#  endif 
    226              ikbu   = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) 
    227              ikbv   = MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) ) 
    228              IF ( bfrcoef2d(ji,jj).gt.2.0*fse3u(ji,jj,ikbu)/rdt ) then 
    229                 IF ( ln_ctl ) THEN 
    230                    write(numout,*) 'BFR ',narea,nimpp+ji,njmpp+jj,ikbu 
    231                    write(numout,*) 'BFR ',bfrcoef2d(ji,jj),2.0*fse3u(ji,jj,ikbu)/rdt 
     220             ikbu = MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) 
     221             ikbv = MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) ) 
     222             zfru = 0.5 * fse3u(ji,jj,ikbu) / rdt 
     223             zfrv = 0.5 * fse3v(ji,jj,ikbv) / rdt 
     224             IF( ABS( bfrcoef2d(ji,jj) ) > zfru ) THEN 
     225                IF( ln_ctl ) THEN 
     226                   WRITE(numout,*) 'BFR ',narea,nimpp+ji,njmpp+jj,ikbu 
     227                   WRITE(numout,*) 'BFR ',ABS( bfrcoef2d(ji,jj) ), zfru 
    232228                ENDIF 
    233229                ictu = ictu + 1 
    234230             ENDIF 
    235              IF ( bfrcoef2d(ji,jj).gt.2.0*fse3v(ji,jj,ikbv)/rdt ) then 
    236                  IF ( ln_ctl ) THEN 
    237                      write(numout,*) 'BFR ',narea,nimpp+ji,njmpp+jj,ikbv 
    238                      write(numout,*) 'BFR ',bfrcoef2d(ji,jj),2.0*fse3v(ji,jj,ikbv)/rdt 
     231             IF( ABS( bfrcoef2d(ji,jj) ) > zfrv ) THEN 
     232                 IF( ln_ctl ) THEN 
     233                     WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbv 
     234                     WRITE(numout,*) 'BFR ', bfrcoef2d(ji,jj), zfrv 
    239235                 ENDIF 
    240236                 ictv = ictv + 1 
    241237             ENDIF 
    242              bfrcoef2d(ji,jj) = MIN(2.0*fse3u(ji,jj,ikbu)/rdt, bfrcoef2d(ji,jj)) 
    243              bfrcoef2d(ji,jj) = MIN(2.0*fse3v(ji,jj,ikbv)/rdt, bfrcoef2d(ji,jj)) 
    244              rminbfr = MIN(rminbfr, bfrcoef2d(ji,jj)) 
    245              rmaxbfr = MAX(rmaxbfr, bfrcoef2d(ji,jj)) 
     238             zminbfr = MIN(  zminbfr, MIN( zfru, ABS( bfrcoef2d(ji,jj) ) )  ) 
     239             zmaxbfr = MAX(  zmaxbfr, MIN( zfrv, ABS( bfrcoef2d(ji,jj) ) )  ) 
    246240         END DO 
    247241      END DO 
    248       IF ( lk_mpp ) THEN 
     242      IF( lk_mpp ) THEN 
    249243         CALL mpp_sum( ictu ) 
    250244         CALL mpp_sum( ictv ) 
    251          CALL mpp_min( rminbfr ) 
    252          CALL mpp_max( rmaxbfr ) 
     245         CALL mpp_min( zminbfr ) 
     246         CALL mpp_max( zmaxbfr ) 
    253247      ENDIF 
    254       IF ( lwp .AND. ictu + ictv .GT. 0 ) THEN 
     248      IF( lwp .AND. ictu + ictv .GT. 0 ) THEN 
    255249         WRITE(numout,*) ' Bottom friction stability check failed at ', ictu, ' U-points ' 
    256250         WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points ' 
    257          WRITE(numout,*) ' Bottom friction coefficient reset where necessary' 
    258          WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', rminbfr, ' to ', rmaxbfr 
     251         WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zminbfr, ' to ', zmaxbfr 
     252         WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary' 
    259253      ENDIF 
    260254      ! 
Note: See TracChangeset for help on using the changeset viewer.