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 4381 for branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90 – NEMO

Ignore:
Timestamp:
2014-01-29T15:26:21+01:00 (10 years ago)
Author:
jchanut
Message:

bug with log bottom friction and namelist update

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r4147 r4381  
    1313 
    1414   !!---------------------------------------------------------------------- 
    15    !!   zdf_bfr      : update momentum Kz at the ocean bottom due to the type of bottom friction chosen 
     15   !!   zdf_bfr      : update bottom friction coefficient (non-linear bottom friction only) 
    1616   !!   zdf_bfr_init : read in namelist and control the bottom friction parameters. 
    17    !!   zdf_bfr_2d   : read in namelist and control the bottom friction parameters. 
    1817   !!---------------------------------------------------------------------- 
    1918   USE oce             ! ocean dynamics and tracers variables 
     
    2524   USE prtctl          ! Print control 
    2625   USE timing          ! Timing 
    27  
     26   USE wrk_nemo        ! Memory Allocation 
    2827   USE phycst, ONLY: vkarmn 
    2928 
     
    3231 
    3332   PUBLIC   zdf_bfr         ! called by step.F90 
    34    PUBLIC   zdf_bfr_init    ! called by opa.F90 
     33   PUBLIC   zdf_bfr_init    ! called by nemogcm.F90 
    3534 
    3635   !                                 !!* Namelist nambfr: bottom friction namelist * 
     
    3837   REAL(wp), PUBLIC ::   rn_bfri1     ! bottom drag coefficient (linear case)  (PUBLIC for TAM) 
    3938   REAL(wp), PUBLIC ::   rn_bfri2     ! bottom drag coefficient (non linear case) (PUBLIC for TAM) 
     39   REAL(wp), PUBLIC ::   rn_bfri2_max ! Maximum bottom drag coefficient (non linear case and ln_loglayer=T) (PUBLIC for TAM) 
    4040   REAL(wp), PUBLIC ::   rn_bfeb2     ! background bottom turbulent kinetic energy  [m2/s2] (PUBLIC for TAM) 
    4141   REAL(wp), PUBLIC ::   rn_bfrien    ! local factor to enhance coefficient bfri (PUBLIC for TAM) 
     
    4343   LOGICAL , PUBLIC ::   ln_loglayer  ! switch for log layer bfr coeff. (PUBLIC for TAM) 
    4444   REAL(wp), PUBLIC ::   rn_bfrz0     ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM) 
    45    LOGICAL , PUBLIC ::  ln_bfrimp     ! logical switch for implicit bottom friction 
     45   LOGICAL , PUBLIC ::   ln_bfrimp    ! logical switch for implicit bottom friction 
    4646   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::  bfrcoef2d            ! 2D bottom drag coefficient (PUBLIC for TAM) 
    4747 
     
    8686      !! 
    8787      INTEGER  ::   ji, jj                       ! dummy loop indices 
    88       INTEGER  ::   ikbu, ikbv                   ! local integers 
    89       REAL(wp) ::   zvu, zuv, zecu, zecv         ! temporary scalars 
    90       REAL(wp) ::   ztmp, ztmp1                  ! temporary scalars 
     88      INTEGER  ::   ikbt, ikbu, ikbv             ! local integers 
     89      REAL(wp) ::   zvu, zuv, zecu, zecv, ztmp   ! temporary scalars 
     90      REAL(wp), POINTER, DIMENSION(:,:) ::  zbfrt 
    9191      !!---------------------------------------------------------------------- 
    9292      ! 
    9393      IF( nn_timing == 1 )  CALL timing_start('zdf_bfr') 
    9494      ! 
    95       IF( nn_bfr == 2 ) THEN                 ! quadratic botton friction 
    96          ! Calculate and store the quadratic bottom friction coefficient bfrua and bfrva 
    97          ! where bfrUa = C_d*SQRT(u_bot^2 + v_bot^2 + e_b) {U=[u,v]} 
    98          ! from these the trend due to bottom friction:  -F_h/e3U  can be calculated 
    99          ! where -F_h/e3U_bot = bfrUa*Ub/e3U_bot {U=[u,v]} 
    100          ! 
    101  
    102          IF(ln_loglayer) THEN       ! "log layer" bottom friction coefficient 
    103  
    104            ! add 2D-enhancement bottom friction 
    105            ztmp1 = 1._wp 
    106            IF(ABS(rn_bfri2) >= 1.e-10 ) THEN 
    107              ztmp1 = 1._wp / rn_bfri2 
    108            ELSE 
    109              CALL ctl_stop( 'rn_bfri2 must not be less than 1.e-10') 
    110            END IF 
     95      IF( kt == nit000 .AND. lwp ) THEN 
     96         WRITE(numout,*) 
     97         WRITE(numout,*) 'zdf_bfr : Set bottom friction coefficient (non-linear case)' 
     98         WRITE(numout,*) '~~~~~~~~' 
     99      ENDIF 
     100      ! 
     101      IF( nn_bfr == 2 ) THEN                 ! quadratic bottom friction only 
     102         ! 
     103         CALL wrk_alloc( jpi, jpj, zbfrt ) 
     104 
     105         IF ( ln_loglayer.AND.lk_vvl ) THEN ! "log layer" bottom friction coefficient 
    111106 
    112107#  if defined key_vectopt_loop 
    113            DO jj = 1, 1 
    114              DO ji = 1, jpij   ! vector opt. (forced unrolling) 
     108            DO jj = 1, 1 
     109!CDIR NOVERRCHK 
     110               DO ji = 1, jpij   ! vector opt. (forced unrolling) 
    115111#  else 
    116            DO jj = 1, jpj 
    117              DO ji = 1, jpi 
     112!CDIR NOVERRCHK 
     113            DO jj = 1, jpj 
     114!CDIR NOVERRCHK 
     115               DO ji = 1, jpi 
    118116#  endif 
    119                 ztmp = 0.5_wp * fse3t(ji,jj,mbkt(ji,jj)) 
    120                 ztmp = max(ztmp, rn_bfrz0 + 1.e-10) 
    121                 bfrcoef2d(ji,jj) = bfrcoef2d(ji,jj) * ztmp1 * & 
    122                                  & ( log( ztmp / rn_bfrz0 ) / vkarmn ) ** (-2) 
    123              END DO 
    124            END DO 
     117                  ikbt = mbkt(ji,jj) 
     118! JC: possible WAD implementation should modify line below if layers vanish 
     119                  ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp 
     120                  zbfrt(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp) 
     121                  zbfrt(ji,jj) = MIN(zbfrt(ji,jj), rn_bfri2_max) 
     122               END DO 
     123            END DO 
     124         !    
     125         ELSE 
     126            zbfrt(:,:) = bfrcoef2d(:,:) 
    125127         ENDIF 
    126128 
     
    146148               zecv = SQRT(  vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2  ) 
    147149               ! 
    148                bfrua(ji,jj) = - 0.5_wp * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji+1,jj  ) ) * zecu 
    149                bfrva(ji,jj) = - 0.5_wp * ( bfrcoef2d(ji,jj) + bfrcoef2d(ji  ,jj+1) ) * zecv 
     150               bfrua(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji+1,jj  ) ) * zecu 
     151               bfrva(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji  ,jj+1) ) * zecv 
    150152            END DO 
    151153         END DO 
     
    156158         IF(ln_ctl)   CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr  - u: ', mask1=umask,        & 
    157159            &                       tab2d_2=bfrva, clinfo2=       ' v: ', mask2=vmask,ovlap=1 ) 
    158       ENDIF 
    159  
     160         CALL wrk_dealloc( jpi,jpj,zbfrt ) 
     161      ENDIF 
    160162      ! 
    161163      IF( nn_timing == 1 )  CALL timing_stop('zdf_bfr') 
     
    170172      !! ** Purpose :   Initialization of the bottom friction 
    171173      !! 
    172       !! ** Method  :   Read the nammbf namelist and check their consistency 
    173       !!              called at the first timestep (nit000) 
     174      !! ** Method  :   Read the nambfr namelist and check their consistency 
     175      !!                called at the first timestep (nit000) 
    174176      !!---------------------------------------------------------------------- 
    175177      USE iom   ! I/O module for ehanced bottom friction file 
     
    177179      INTEGER   ::   inum         ! logical unit for enhanced bottom friction file 
    178180      INTEGER   ::   ji, jj       ! dummy loop indexes 
    179       INTEGER   ::   ikbu, ikbv   ! temporary integers 
    180       INTEGER   ::   ictu, ictv   !    -          - 
     181      INTEGER   ::   ikbt, ikbu, ikbv   ! temporary integers 
     182      INTEGER   ::   ictu, ictv         !    -          - 
    181183      INTEGER   ::   ios 
    182184      REAL(wp)  ::   zminbfr, zmaxbfr   ! temporary scalars 
    183       REAL(wp)  ::   zfru, zfrv         !    -         - 
    184       !! 
    185       NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfeb2, rn_bfrz0, ln_bfr2d, & 
     185      REAL(wp)  ::   ztmp, zfru, zfrv   !    -         - 
     186      !! 
     187      NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfri2_max, rn_bfeb2, rn_bfrz0, ln_bfr2d, & 
    186188                    &  rn_bfrien, ln_bfrimp, ln_loglayer 
    187189      !!---------------------------------------------------------------------- 
    188190      ! 
    189191      IF( nn_timing == 1 )  CALL timing_start('zdf_bfr_init') 
     192      ! 
     193      !                              !* Allocate zdfbfr arrays 
     194      IF( zdf_bfr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_bfr_init : unable to allocate arrays' ) 
     195      ! 
     196      !                              !* Parameter control and print 
    190197      ! 
    191198      REWIND( numnam_ref )              ! Namelist nambfr in reference namelist : Bottom momentum boundary condition 
     
    197204902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambfr in configuration namelist', lwp ) 
    198205      WRITE ( numond, nambfr ) 
    199  
    200       !                               !* Parameter control and print 
    201206      IF(lwp) WRITE(numout,*) 
    202       IF(lwp) WRITE(numout,*) 'zdf_bfr : momentum bottom friction' 
    203       IF(lwp) WRITE(numout,*) '~~~~~~~' 
     207      IF(lwp) WRITE(numout,*) 'zdf_bfr_init : momentum bottom friction' 
     208      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~' 
    204209      IF(lwp) WRITE(numout,*) '   Namelist nam_bfr : set bottom friction parameters' 
    205  
    206       !                              ! allocate zdfbfr arrays 
    207       IF( zdf_bfr_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_bfr_init : unable to allocate arrays' ) 
    208  
     210      ! 
     211      SELECT CASE (nn_bfr) 
     212      ! 
     213      CASE( 0 ) 
     214         IF(lwp) WRITE(numout,*) '      free-slip ' 
     215         bfrua(:,:) = 0.e0 
     216         bfrva(:,:) = 0.e0 
     217         ! 
     218      CASE( 1 ) 
     219         IF(lwp) WRITE(numout,*) '      linear botton friction' 
     220         IF(lwp) WRITE(numout,*) '      friction coef.   rn_bfri1  = ', rn_bfri1 
     221         IF( ln_bfr2d ) THEN 
     222            IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d 
     223            IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien 
     224         ENDIF 
     225         ! 
     226         IF(ln_bfr2d) THEN 
     227            ! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 
     228            CALL iom_open('bfr_coef.nc',inum) 
     229            CALL iom_get (inum, jpdom_data, 'bfr_coef',bfrcoef2d,1) ! bfrcoef2d is used as tmp array 
     230            CALL iom_close(inum) 
     231            bfrcoef2d(:,:) = rn_bfri1 * ( 1 + rn_bfrien * bfrcoef2d(:,:) ) 
     232         ELSE 
     233            bfrcoef2d(:,:) = rn_bfri1  ! initialize bfrcoef2d to the namelist variable 
     234         ENDIF 
     235         ! 
     236         bfrua(:,:) = - bfrcoef2d(:,:) 
     237         bfrva(:,:) = - bfrcoef2d(:,:) 
     238         ! 
     239      CASE( 2 ) 
     240         IF(lwp) WRITE(numout,*) '      quadratic bottom friction' 
     241         IF(lwp) WRITE(numout,*) '      friction coef.   rn_bfri2  = ', rn_bfri2 
     242         IF(lwp) WRITE(numout,*) '      Max. coef. (log case)   rn_bfri2_max  = ', rn_bfri2_max 
     243         IF(lwp) WRITE(numout,*) '      background tke   rn_bfeb2  = ', rn_bfeb2 
     244         IF(lwp) WRITE(numout,*) '      log formulation   ln_bfr2d = ', ln_loglayer 
     245         IF(lwp) WRITE(numout,*) '      bottom roughness  rn_bfrz0 [m] = ', rn_bfrz0 
     246         IF( rn_bfrz0<=0.e0 ) THEN 
     247            WRITE(ctmp1,*) '      bottom roughness must be strictly positive' 
     248            CALL ctl_stop( ctmp1 ) 
     249         ENDIF 
     250         IF( ln_bfr2d ) THEN 
     251            IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d 
     252            IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien 
     253         ENDIF 
     254         ! 
     255         IF(ln_bfr2d) THEN 
     256            ! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 
     257            CALL iom_open('bfr_coef.nc',inum) 
     258            CALL iom_get (inum, jpdom_data, 'bfr_coef',bfrcoef2d,1) ! bfrcoef2d is used as tmp array 
     259            CALL iom_close(inum) 
     260            ! 
     261            bfrcoef2d(:,:) = rn_bfri2 * ( 1 + rn_bfrien * bfrcoef2d(:,:) ) 
     262         ELSE 
     263            bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable 
     264         ENDIF 
     265         ! 
     266         IF ( ln_loglayer.AND.(.NOT.lk_vvl) ) THEN ! set "log layer" bottom friction once for all 
     267#  if defined key_vectopt_loop 
     268            DO jj = 1, 1 
     269!CDIR NOVERRCHK 
     270               DO ji = 1, jpij   ! vector opt. (forced unrolling) 
     271#  else 
     272!CDIR NOVERRCHK 
     273            DO jj = 1, jpj 
     274!CDIR NOVERRCHK 
     275               DO ji = 1, jpi 
     276#  endif 
     277                  ikbt = mbkt(ji,jj) 
     278                  ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp 
     279                  bfrcoef2d(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp) 
     280                  bfrcoef2d(ji,jj) = MIN(bfrcoef2d(ji,jj), rn_bfri2_max) 
     281               END DO 
     282            END DO 
     283         ENDIF 
     284         ! 
     285      CASE DEFAULT 
     286         IF(lwp) WRITE(ctmp1,*) '         bad flag value for nn_bfr = ', nn_bfr 
     287         CALL ctl_stop( ctmp1 ) 
     288         ! 
     289      END SELECT 
     290      ! 
     291      IF(lwp) WRITE(numout,*) '      implicit bottom friction switch                ln_bfrimp  = ', ln_bfrimp 
     292      ! 
    209293      !                              ! Make sure ln_zdfexp=.false. when use implicit bfr 
    210294      IF( ln_bfrimp .AND. ln_zdfexp ) THEN 
     
    217301         END IF 
    218302      END IF 
    219  
    220       SELECT CASE (nn_bfr) 
    221       ! 
    222       CASE( 0 ) 
    223          IF(lwp) WRITE(numout,*) '      free-slip ' 
    224          bfrua(:,:) = 0.e0 
    225          bfrva(:,:) = 0.e0 
    226          ! 
    227       CASE( 1 ) 
    228          IF(lwp) WRITE(numout,*) '      linear botton friction' 
    229          IF(lwp) WRITE(numout,*) '      friction coef.   rn_bfri1  = ', rn_bfri1 
    230          IF( ln_bfr2d ) THEN 
    231             IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d 
    232             IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien 
    233          ENDIF 
    234          ! 
    235          bfrcoef2d(:,:) = rn_bfri1  ! initialize bfrcoef2d to the namelist variable 
    236          ! 
    237          IF(ln_bfr2d) THEN 
    238             ! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 
    239             CALL iom_open('bfr_coef.nc',inum) 
    240             CALL iom_get (inum, jpdom_data, 'bfr_coef',bfrcoef2d,1) ! bfrcoef2d is used as tmp array 
    241             CALL iom_close(inum) 
    242             bfrcoef2d(:,:) = rn_bfri1 * ( 1 + rn_bfrien * bfrcoef2d(:,:) ) 
    243          ENDIF 
    244          bfrua(:,:) = - bfrcoef2d(:,:) 
    245          bfrva(:,:) = - bfrcoef2d(:,:) 
    246          ! 
    247       CASE( 2 ) 
    248          IF(lwp) WRITE(numout,*) '      quadratic botton friction' 
    249          IF(lwp) WRITE(numout,*) '      friction coef.   rn_bfri2  = ', rn_bfri2 
    250          IF(lwp) WRITE(numout,*) '      background tke   rn_bfeb2  = ', rn_bfeb2 
    251          IF( ln_bfr2d ) THEN 
    252             IF(lwp) WRITE(numout,*) '      read coef. enhancement distribution from file   ln_bfr2d  = ', ln_bfr2d 
    253             IF(lwp) WRITE(numout,*) '      coef rn_bfri2 enhancement factor                rn_bfrien  = ',rn_bfrien 
    254          ENDIF 
    255          bfrcoef2d(:,:) = rn_bfri2  ! initialize bfrcoef2d to the namelist variable 
    256  
    257          ! 
    258          IF(ln_bfr2d) THEN 
    259             ! bfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 
    260             CALL iom_open('bfr_coef.nc',inum) 
    261             CALL iom_get (inum, jpdom_data, 'bfr_coef',bfrcoef2d,1) ! bfrcoef2d is used as tmp array 
    262             CALL iom_close(inum) 
    263             bfrcoef2d(:,:)= rn_bfri2 * ( 1 + rn_bfrien * bfrcoef2d(:,:) ) 
    264          ENDIF 
    265          ! 
    266       CASE DEFAULT 
    267          IF(lwp) WRITE(ctmp1,*) '         bad flag value for nn_bfr = ', nn_bfr 
    268          CALL ctl_stop( ctmp1 ) 
    269          ! 
    270       END SELECT 
    271  
    272       IF( nn_bfr /= 2 .AND. ln_loglayer ) THEN 
    273          IF(lwp) THEN 
    274             WRITE(numout,*) 
    275             WRITE(numout,*) 'Loglayer can only be by applied for quadratic bottom friction'  
    276             WRITE(numout,*) 'but you have set: nn_bfr /= 2 and ln_loglayer=.true.!!!!' 
    277             WRITE(ctmp1,*)  'check nn_bfr and ln_loglayer (should be 2 and true)' 
    278             CALL ctl_stop( ctmp1 ) 
    279          END IF 
    280       END IF 
    281  
    282  
    283  
    284       IF(lwp) WRITE(numout,*) '      implicit bottom friction switch                ln_bfrimp  = ', ln_bfrimp 
    285303      ! 
    286304      ! Basic stability check on bottom friction coefficient 
Note: See TracChangeset for help on using the changeset viewer.