Ignore:
Timestamp:
2011-11-08T13:39:35+01:00 (9 years ago)
Author:
hliu
Message:

update semi-implicit bottom friction branch, Document has been added in NEMO_book Chapter 10.4.4

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r2872 r3057  
    120120      INTEGER  ::   icycle           ! local scalar 
    121121      INTEGER  ::   ikbu, ikbv       ! local scalar 
    122       REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b     ! local scalars 
    123       REAL(wp) ::   z1_8, zx1, zy1                   !   -      - 
    124       REAL(wp) ::   z1_4, zx2, zy2                   !   -      - 
    125       REAL(wp) ::   zu_spg, zu_cor, zu_sld, zu_asp   !   -      - 
    126       REAL(wp) ::   zv_spg, zv_cor, zv_sld, zv_asp   !   -      - 
    127       REAL(wp) ::   ua_btm, va_btm                   !   -      - 
     122      REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b, z2dt_bf     ! local scalars 
     123      REAL(wp) ::   z1_8, zx1, zy1                            !   -      - 
     124      REAL(wp) ::   z1_4, zx2, zy2                            !   -      - 
     125      REAL(wp) ::   zu_spg, zu_cor, zu_sld, zu_asp            !   -      - 
     126      REAL(wp) ::   zv_spg, zv_cor, zv_sld, zv_asp            !   -      - 
     127      REAL(wp) ::   ua_btm, va_btm                            !   -      - 
    128128      !!---------------------------------------------------------------------- 
    129129 
     
    171171      zu_sld = 0._wp   ;   zu_asp = 0._wp                   ! tides trends (lk_tide=F) 
    172172      zv_sld = 0._wp   ;   zv_asp = 0._wp 
     173 
     174      IF( kt == nit000 .AND. neuler == 0) THEN              ! for implicit bottom friction 
     175        z2dt_bf = rdt 
     176      ELSE 
     177        z2dt_bf = 2.0_wp * rdt 
     178      ENDIF 
    173179 
    174180      ! ----------------------------------------------------------------------------- 
     
    207213      END DO 
    208214 
     215      !                                   !* baroclinic momentum trend (remove the vertical mean trend) 
     216      DO jk = 1, jpkm1                    ! -------------------------- 
     217         DO jj = 2, jpjm1 
     218            DO ji = fs_2, fs_jpim1   ! vector opt. 
     219               ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj) 
     220               va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj) 
     221            END DO 
     222         END DO 
     223      END DO 
    209224 
    210225      !                                   !* barotropic Coriolis trends * H (vorticity scheme dependent) 
     
    253268         DO jj = 2, jpjm1  
    254269            DO ji = fs_2, fs_jpim1   ! vector opt. 
    255       !     remove the rhd term according to J. Chanut's suggestion             
    256270               zcu(ji,jj) = zcu(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * hu(ji,jj) / e1u(ji,jj) 
    257271               zcv(ji,jj) = zcv(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * hv(ji,jj) / e2v(ji,jj) 
     
    260274      ENDIF 
    261275 
    262      IF(ln_bfrimp.AND.lk_vvl) THEN 
    263 !    Remove the bottom stress trend from 3-D sea surface level gradient and Coriolis force            
    264 !    H. Liu, Sept. 2011 
    265       DO jj = 2, jpjm1          
    266          DO ji = fs_2, fs_jpim1 
    267             ikbu = mbku(ji,jj) 
    268             ikbv = mbkv(ji,jj) 
    269             ua_btm = zcu(ji,jj) * z2dt_b * hur(ji,jj) * umask (ji,jj,ikbu) 
    270             va_btm = zcv(ji,jj) * z2dt_b * hvr(ji,jj) * vmask (ji,jj,ikbv) 
    271  
    272             ua(ji,jj,ikbu) = ua(ji,jj,ikbu) - bfrua(ji,jj) * ua_btm / fse3u(ji,jj,ikbu) 
    273             va(ji,jj,ikbv) = va(ji,jj,ikbv) - bfrva(ji,jj) * va_btm / fse3v(ji,jj,ikbv) 
    274              
    275             zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm 
    276             zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm 
    277          END DO 
    278       END DO 
    279      END IF 
    280       !                                   !* baroclinic momentum trend (remove the vertical mean trend) 
    281       DO jk = 1, jpkm1                    ! -------------------------- 
    282          DO jj = 2, jpjm1 
    283             DO ji = fs_2, fs_jpim1   ! vector opt. 
    284                ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj) 
    285                va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj) 
    286             END DO 
    287          END DO 
    288       END DO 
    289  
    290276      DO jj = 2, jpjm1                             ! Remove coriolis term (and possibly spg) from barotropic trend 
    291277         DO ji = fs_2, fs_jpim1 
    292             zua(ji,jj) = zua(ji,jj) - zcu(ji,jj) 
    293             zva(ji,jj) = zva(ji,jj) - zcv(ji,jj) 
    294          END DO 
     278             zua(ji,jj) = zua(ji,jj) - zcu(ji,jj) 
     279             zva(ji,jj) = zva(ji,jj) - zcv(ji,jj) 
     280          END DO 
    295281      END DO 
    296282 
     
    299285      !                                             ! from the barotropic transport trend 
    300286      zcoef = -1. / z2dt_b 
    301       IF(.NOT.ln_bfrimp) THEN 
     287 
     288      IF(ln_bfrimp) THEN 
     289      !                                   ! Remove the bottom stress trend from 3-D sea surface level gradient 
     290      !                                   ! and Coriolis forcing in case of 3D semi-implicit bottom friction  
     291        DO jj = 2, jpjm1          
     292           DO ji = fs_2, fs_jpim1 
     293              ikbu = mbku(ji,jj) 
     294              ikbv = mbkv(ji,jj) 
     295              ua_btm = zcu(ji,jj) * z2dt_bf * hur(ji,jj) * umask (ji,jj,ikbu) 
     296              va_btm = zcv(ji,jj) * z2dt_bf * hvr(ji,jj) * vmask (ji,jj,ikbv) 
     297 
     298              zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm 
     299              zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm 
     300           END DO 
     301        END DO 
     302 
     303      ELSE 
     304 
    302305# if defined key_vectopt_loop 
    303       DO jj = 1, 1 
    304          DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
     306        DO jj = 1, 1 
     307           DO ji = 1, jpij-jpi   ! vector opt. (forced unrolling) 
    305308# else 
    306       DO jj = 2, jpjm1 
    307          DO ji = 2, jpim1 
     309        DO jj = 2, jpjm1 
     310           DO ji = 2, jpim1 
    308311# endif 
    309312            ! Apply stability criteria for bottom friction 
    310313            !RBbug for vvl and external mode we may need to use varying fse3 
    311314            !!gm  Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx. 
    312             zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  ) 
    313             zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  ) 
    314          END DO 
    315       END DO 
    316  
    317       IF( lk_vvl ) THEN 
    318          DO jj = 2, jpjm1 
    319             DO ji = fs_2, fs_jpim1   ! vector opt. 
    320                zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   & 
    321                   &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) ) 
    322                zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj)   & 
    323                   &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) ) 
    324             END DO 
    325          END DO 
    326       ELSE 
    327          DO jj = 2, jpjm1 
    328             DO ji = fs_2, fs_jpim1   ! vector opt. 
    329                zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj) 
    330                zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj) 
    331             END DO 
    332          END DO 
    333       ENDIF 
    334       END IF 
    335  
     315              zbfru(ji,jj) = MAX(  bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef  ) 
     316              zbfrv(ji,jj) = MAX(  bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef  ) 
     317           END DO 
     318        END DO 
     319 
     320        IF( lk_vvl ) THEN 
     321           DO jj = 2, jpjm1 
     322              DO ji = fs_2, fs_jpim1   ! vector opt. 
     323                 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj)   & 
     324                    &       / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) ) 
     325                 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj)   & 
     326                    &       / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) ) 
     327              END DO 
     328           END DO 
     329        ELSE 
     330           DO jj = 2, jpjm1 
     331              DO ji = fs_2, fs_jpim1   ! vector opt. 
     332                 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj) 
     333                 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj) 
     334              END DO 
     335           END DO 
     336        ENDIF 
     337      END IF    ! end (ln_bfrimp) 
     338 
     339                     
    336340      !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity) 
    337341      !                                   ! --------------------------  
     
    443447                  ! after velocities with implicit bottom friction 
    444448 
    445                   IF(ln_bfrimp) THEN 
    446                   !   A new method to implement the implicit bottom friction.  
    447                   !   H. Liu 
    448                   !   May 2010 
    449                   ua_e(ji,jj) =  zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )  * umask(ji,jj,1)   & 
     449                  IF(ln_bfrimp) THEN      !implicit bottom friction 
     450                  ua_e(ji,jj) = ((1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj)) * zub_e(ji,jj) +   & 
     451                     &          z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )) * umask(ji,jj,1)   & 
    450452                     &         / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    451453                  ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)    
    452                   va_e(ji,jj) =  zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )  * vmask(ji,jj,1)   & 
     454                  va_e(ji,jj) = ((1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj)) * zvb_e(ji,jj) +   & 
     455                     &          z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )) * vmask(ji,jj,1)   & 
    453456                     &         / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
    454457                  va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
     
    478481                  !   H. Liu 
    479482                  !   May 2010 
    480                   ua_e(ji,jj) =  zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )  * umask(ji,jj,1)   & 
     483                  ua_e(ji,jj) = ((1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj)) * zub_e(ji,jj) +   &  
     484                     &          z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )) * umask(ji,jj,1)   & 
    481485                     &         / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    482486                  ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)    
    483                   va_e(ji,jj) =  zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )  * vmask(ji,jj,1)   & 
     487                  va_e(ji,jj) = ((1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj)) * zvb_e(ji,jj) +   & 
     488                     &          z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )) * vmask(ji,jj,1)   & 
    484489                     &         / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
    485490                  va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
     
    509514                  !   H. Liu 
    510515                  !   May 2010 
    511                   ua_e(ji,jj) =  zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )  * umask(ji,jj,1)   & 
     516                  ua_e(ji,jj) = ((1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj)) * zub_e(ji,jj) +   & 
     517                     &          z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp )) * umask(ji,jj,1)   & 
    512518                     &         / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    513519                  ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e *   zua(ji,jj)  ) * umask(ji,jj,1)    
    514                   va_e(ji,jj) =  zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )  * vmask(ji,jj,1)   & 
     520                  va_e(ji,jj) = ((1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj)) * zvb_e(ji,jj) +   & 
     521                     &          z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp )) * vmask(ji,jj,1)   & 
    515522                     &         / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
    516523                  va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
Note: See TracChangeset for help on using the changeset viewer.