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 2497 for branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90 – NEMO

Ignore:
Timestamp:
2010-12-21T07:44:40+01:00 (13 years ago)
Author:
gm
Message:

v3.3beta: #780 zdfgls error in bottom friction formulation

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90

    r2452 r2497  
    2828   USE in_out_manager ! I/O manager 
    2929   USE iom            ! I/O manager library 
    30    USE zdfbfr, ONLY : rn_hbro, wbotu, wbotv ! bottom roughness and bottom stresses 
    3130 
    3231   IMPLICIT NONE 
     
    6160 
    6261   REAL(wp) ::   hsro          =  0.003_wp    ! Minimum surface roughness 
     62   REAL(wp) ::   hbro          =  0.003_wp    ! Bottom roughness (m) 
    6363   REAL(wp) ::   rcm_sf        =  0.73_wp     ! Shear free turbulence parameters 
    6464   REAL(wp) ::   ra_sf         = -2.0_wp      ! Must be negative -2 < ra_sf < -1  
     
    143143      ustars2 = 0._wp   ;   ustarb2 = 0._wp   ;   psi  = 0._wp   ;   zwall_psi = 0._wp 
    144144 
    145       ! Compute wind stress at T-points 
     145      ! Compute surface and bottom friction at T-points 
    146146!CDIR NOVERRCHK 
    147147      DO jj = 2, jpjm1 
    148148!CDIR NOVERRCHK 
    149         DO ji = fs_2, fs_jpim1   ! vector opt. 
    150           !  
    151           ! wind stress 
    152           ! squared surface velocity scale 
    153           ustars2(ji,jj) = rau0r * taum(ji,jj) * tmask(ji,jj,1) 
    154           ! 
    155           ! bottom friction  
    156           ztx2 = ( wbotu(ji-1,jj)*umask(ji-1,jj,1) + wbotu(ji,jj)*umask(ji,jj,1) )   & 
    157             &  / MAX( 1.,         umask(ji-1,jj,1) +              umask(ji,jj,1) ) 
    158           zty2 = ( wbotv(ji,jj-1)*vmask(ji,jj-1,1) + wbotv(ji,jj)*vmask(ji,jj,1) )   & 
    159             &  / MAX( 1.,         vmask(ji,jj-1,1) +              vmask(ji,jj,1) ) 
    160           ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 
    161         END DO 
     149         DO ji = fs_2, fs_jpim1   ! vector opt. 
     150            !  
     151            ! surface friction  
     152            ustars2(ji,jj) = rau0r * taum(ji,jj) * tmask(ji,jj,1) 
     153            ! 
     154            ! bottom friction (explicit before friction) 
     155            ! Note that we chose here not to bound the friction as in dynbfr) 
     156            ztx2 = (  bfrua(ji,jj)  * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) )   & 
     157               & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1) ) 
     158            zty2 = (  bfrva(ji,jj)  * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) )   & 
     159               & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) ) 
     160            ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 
     161         END DO 
    162162      END DO   
    163163 
     
    617617      ! 
    618618      CASE ( 0 )             ! Dirichlet  
    619          !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_hbro 
     619         !                      ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * hbro 
    620620         !                      ! Balance between the production and the dissipation terms 
    621621!CDIR NOVERRCHK 
     
    625625               ibot   = mbkt(ji,jj) + 1      ! k   bottom level of w-point 
    626626               ibotm1 = mbkt(ji,jj)          ! k-1 bottom level of w-point but >=1 
    627                zdep(ji,jj) = vkarmn * rn_hbro 
     627               zdep(ji,jj) = vkarmn * hbro 
    628628               psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 
    629629               z_elem_a(ji,jj,ibot) = 0._wp 
     
    632632               ! 
    633633               ! Just above last level, Dirichlet condition again (GOTM like) 
    634                zdep(ji,jj) = vkarmn * (rn_hbro + fse3t(ji,jj,ibotm1)) 
     634               zdep(ji,jj) = vkarmn * ( hbro + fse3t(ji,jj,ibotm1) ) 
    635635               psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot  )**rmm * zdep(ji,jj)**rnn 
    636636               z_elem_a(ji,jj,ibotm1) = 0._wp 
     
    650650               ! 
    651651               ! Bottom level Dirichlet condition: 
    652                zdep(ji,jj) = vkarmn * rn_hbro 
     652               zdep(ji,jj) = vkarmn * hbro 
    653653               psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 
    654654               ! 
     
    662662               ! 
    663663               ! Set psi vertical flux at the bottom: 
    664                zdep(ji,jj) = rn_hbro + 0.5_wp*fse3t(ji,jj,ibotm1) 
     664               zdep(ji,jj) = hbro + 0.5_wp*fse3t(ji,jj,ibotm1) 
    665665               zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) )   & 
    666666                  &  * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 
     
    906906         WRITE(numout,*) '      minimum value of en                           rn_emin       = ', rn_emin 
    907907         WRITE(numout,*) '      minimum value of eps                          rn_epsmin     = ', rn_epsmin 
    908          WRITE(numout,*) '      Surface roughness (m)                         hsro          = ', hsro 
    909          WRITE(numout,*) '      Bottom roughness (m)                          rn_hbro       = ', rn_hbro 
    910908         WRITE(numout,*) '      Limit dissipation rate under stable stratif.  ln_length_lim = ', ln_length_lim 
    911909         WRITE(numout,*) '      Galperin limit (Standard: 0.53, Holt: 0.26)   rn_clim_galp  = ', rn_clim_galp 
     
    920918         WRITE(numout,*) '      Stability functions                           nn_stab_func   = ', nn_stab_func 
    921919         WRITE(numout,*) '      Type of closure                               nn_clos        = ', nn_clos 
    922          WRITE(numout,*) 
     920         WRITE(numout,*) '   Hard coded parameters' 
     921         WRITE(numout,*) '      Surface roughness (m)                         hsro          = ', hsro 
     922         WRITE(numout,*) '      Bottom roughness (m)                          hbro          = ', hbro 
    923923      ENDIF 
    924924 
     
    11951195            id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. ) 
    11961196            id6 = iom_varid( numror, 'mxln' , ldstop = .FALSE. ) 
    1197             id7 = iom_varid( numror, 'wbotu', ldstop = .FALSE. ) 
    1198             id8 = iom_varid( numror, 'wbotv', ldstop = .FALSE. ) 
    11991197            ! 
    12001198            IF( MIN( id1, id2, id3, id4, id5, id6, id7, id8 ) > 0 ) THEN        ! all required arrays exist 
     
    12051203               CALL iom_get( numror, jpdom_autoglo, 'avmv'  , avmv   ) 
    12061204               CALL iom_get( numror, jpdom_autoglo, 'mxln'  , mxln   ) 
    1207                CALL iom_get( numror, jpdom_autoglo, 'wbotu' , wbotu  ) 
    1208                CALL iom_get( numror, jpdom_autoglo, 'wbotv' , wbotv  ) 
    12091205            ELSE                         
    12101206               IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop' 
    1211                IF(lwp) WRITE(numout,*) ' ===>>>> : The bottom stresses are estimated'  
    12121207               en  (:,:,:) = rn_emin 
    12131208               mxln(:,:,:) = 0.001         
    1214                ! Initialize bottom stresses 
    1215                DO jj = 2, jpjm1 
    1216                   DO ji = fs_2, fs_jpim1   ! vector opt. 
    1217                      ikbu = mbku(ji,jj)         ! bottom ocean level of u-point 
    1218                      ikbv = mbkv(ji,jj) 
    1219                      cbx  = avmu(ji,jj,ikbu+1) / fse3uw(ji,jj,ikbu+1) 
    1220                      cby  = avmv(ji,jj,ikbv+1) / fse3vw(ji,jj,ikbv+1) 
    1221                      wbotu(ji,jj) = -cbx * un(ji,jj,ikbu) * umask(ji,jj,1) 
    1222                      wbotv(ji,jj) = -cby * vn(ji,jj,ikbv) * vmask(ji,jj,1) 
    1223                   END DO 
    1224                END DO 
    12251209               DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_gls( jit )   ;   END DO 
    12261210            ENDIF 
    12271211         ELSE                                   !* Start from rest 
    12281212            IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values' 
    1229             IF(lwp) WRITE(numout,*) ' ===>>>> : The bottom stresses are estimated' 
    12301213            en  (:,:,:) = rn_emin 
    12311214            mxln(:,:,:) = 0.001        
    1232             ! Initialize bottom stresses 
    1233             DO jj = 2, jpjm1 
    1234                DO ji = fs_2, fs_jpim1   ! vector opt. 
    1235                   ikbu = mbku(ji,jj)            ! bottom ocean level of u-point 
    1236                   ikbv = mbkv(ji,jj) 
    1237                   cbx  = avmu(ji,jj,ikbu+1) / fse3uw(ji,jj,ikbu+1) 
    1238                   cby  = avmv(ji,jj,ikbv+1) / fse3vw(ji,jj,ikbv+1) 
    1239                   wbotu(ji,jj) = -cbx * un(ji,jj,ikbu) * umask(ji,jj,1) 
    1240                   wbotv(ji,jj) = -cby * vn(ji,jj,ikbv) * vmask(ji,jj,1) 
    1241                END DO 
    1242             END DO 
    12431215         ENDIF 
    12441216         ! 
     
    12521224         CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv  ) 
    12531225         CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln  ) 
    1254          CALL iom_rstput( kt, nitrst, numrow, 'wbotu' , wbotu ) 
    1255          CALL iom_rstput( kt, nitrst, numrow, 'wbotv' , wbotv ) 
    12561226         ! 
    12571227      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.