Changeset 2872


Ignore:
Timestamp:
2011-09-28T10:18:52+02:00 (9 years ago)
Author:
hliu
Message:

1) added the semi-implicit bottom friction code (3 in DYN, 1 in ZDF), 2) added par_AMM7/12.h90, 3) added two arch files for NOCL mobius and ubuntu linux

Location:
branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM
Files:
4 added
5 edited

Legend:

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

    r2715 r2872  
    1313   USE dom_oce         ! ocean space and time domain variables  
    1414   USE zdf_oce         ! ocean vertical physics variables 
     15   USE zdfbfr          ! ocean bottom friction variables 
    1516   USE trdmod          ! ocean active dynamics and tracers trends  
    1617   USE trdmod_oce      ! ocean variables trends 
     
    5152      !!--------------------------------------------------------------------- 
    5253      ! 
    53       zm1_2dt = - 1._wp / ( 2._wp * rdt ) 
     54      IF( .not. ln_bfrimp) THEN     ! only for explicit bottom friction form 
     55                                    ! implicit bfr is implemented in dynzdf_imp 
     56                                    ! H. Liu,  Sept. 2011 
    5457 
    55       IF( l_trddyn )   THEN                      ! temporary save of ua and va trends 
    56          ztrduv(:,:,:,1) = ua(:,:,:) 
    57          ztrduv(:,:,:,2) = va(:,:,:) 
    58       ENDIF 
     58        zm1_2dt = - 1._wp / ( 2._wp * rdt ) 
     59 
     60        IF( l_trddyn )   THEN                      ! temporary save of ua and va trends 
     61           ztrduv(:,:,:,1) = ua(:,:,:) 
     62           ztrduv(:,:,:,2) = va(:,:,:) 
     63        ENDIF 
     64 
    5965 
    6066# if defined key_vectopt_loop 
    61       DO jj = 1, 1 
    62          DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
     67        DO jj = 1, 1 
     68           DO ji = jpi+2, jpij-jpi-1   ! vector opt. (forced unrolling) 
    6369# else 
    64       DO jj = 2, jpjm1 
    65          DO ji = 2, jpim1 
     70        DO jj = 2, jpjm1 
     71           DO ji = 2, jpim1 
    6672# endif 
    67             ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels 
    68             ikbv = mbkv(ji,jj) 
    69             ! 
    70             ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 
    71             ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  bfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu) 
    72             va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  bfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv) 
    73          END DO 
    74       END DO 
     73              ikbu = mbku(ji,jj)          ! deepest ocean u- & v-levels 
     74              ikbv = mbkv(ji,jj) 
     75              ! 
     76              ! Apply stability criteria on absolute value  : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 
     77              ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX(  bfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt  ) * ub(ji,jj,ikbu) 
     78              va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX(  bfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt  ) * vb(ji,jj,ikbv) 
     79           END DO 
     80        END DO 
    7581 
    76       ! 
    77       IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics 
    78          ztrduv(:,:,:,1) = ua(:,:,:) - ztrduv(:,:,:,1) 
    79          ztrduv(:,:,:,2) = va(:,:,:) - ztrduv(:,:,:,2) 
    80          CALL trd_mod( ztrduv(:,:,:,1), ztrduv(:,:,:,2), jpdyn_trd_bfr, 'DYN', kt ) 
    81       ENDIF 
    82       !                                          ! print mean trends (used for debugging) 
    83       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' bfr  - Ua: ', mask1=umask,               & 
    84          &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
    85       ! 
     82        ! 
     83        IF( l_trddyn )   THEN                      ! save the vertical diffusive trends for further diagnostics 
     84           ztrduv(:,:,:,1) = ua(:,:,:) - ztrduv(:,:,:,1) 
     85           ztrduv(:,:,:,2) = va(:,:,:) - ztrduv(:,:,:,2) 
     86           CALL trd_mod( ztrduv(:,:,:,1), ztrduv(:,:,:,2), jpdyn_trd_bfr, 'DYN', kt ) 
     87        ENDIF 
     88        !                                          ! print mean trends (used for debugging) 
     89        IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' bfr  - Ua: ', mask1=umask,               & 
     90           &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' ) 
     91        ! 
     92      ENDIF     ! end explicit bottom friction 
    8693   END SUBROUTINE dyn_bfr 
    8794 
  • branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90

    r2724 r2872  
    119119      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices 
    120120      INTEGER  ::   icycle           ! local scalar 
     121      INTEGER  ::   ikbu, ikbv       ! local scalar 
    121122      REAL(wp) ::   zraur, zcoef, z2dt_e, z2dt_b     ! local scalars 
    122123      REAL(wp) ::   z1_8, zx1, zy1                   !   -      - 
     
    124125      REAL(wp) ::   zu_spg, zu_cor, zu_sld, zu_asp   !   -      - 
    125126      REAL(wp) ::   zv_spg, zv_cor, zv_sld, zv_asp   !   -      - 
     127      REAL(wp) ::   ua_btm, va_btm                   !   -      - 
    126128      !!---------------------------------------------------------------------- 
    127129 
     
    147149         hvr_e (:,:) = hvr  (:,:) 
    148150         IF( ln_dynvor_een ) THEN 
    149             ftne(1,:) = 0.e0   ;   ftnw(1,:) = 0.e0   ;   ftse(1,:) = 0.e0   ;   ftsw(1,:) = 0.e0 
     151            ftne(1,:) = 0._wp   ;   ftnw(1,:) = 0._wp   ;   ftse(1,:) = 0._wp   ;   ftsw(1,:) = 0._wp 
    150152            DO jj = 2, jpj 
    151153               DO ji = fs_2, jpi   ! vector opt. 
    152                   ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3. 
    153                   ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3. 
    154                   ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3. 
    155                   ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3. 
     154                  ftne(ji,jj) = ( ff(ji-1,jj  ) + ff(ji  ,jj  ) + ff(ji  ,jj-1) ) / 3._wp 
     155                  ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj  ) + ff(ji  ,jj  ) ) / 3._wp 
     156                  ftse(ji,jj) = ( ff(ji  ,jj  ) + ff(ji  ,jj-1) + ff(ji-1,jj-1) ) / 3._wp 
     157                  ftsw(ji,jj) = ( ff(ji  ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj  ) ) / 3._wp 
    156158               END DO 
    157159            END DO 
     
    161163 
    162164      !                                   !* Local constant initialization 
    163       z2dt_b = 2.0 * rdt                                    ! baroclinic time step 
    164       z1_8 = 0.5 * 0.25                                     ! coefficient for vorticity estimates 
    165       z1_4 = 0.5 * 0.5 
    166       zraur  = 1. / rau0                                    ! 1 / volumic mass 
    167       ! 
    168       zhdiv(:,:) = 0.e0                                     ! barotropic divergence 
    169       zu_sld = 0.e0   ;   zu_asp = 0.e0                     ! tides trends (lk_tide=F) 
    170       zv_sld = 0.e0   ;   zv_asp = 0.e0 
     165      z2dt_b = 2.0_wp * rdt                                 ! baroclinic time step 
     166      z1_8   = 0.5_wp * 0.25_wp                             ! coefficient for vorticity estimates 
     167      z1_4   = 0.5_wp * 0.5_wp 
     168      zraur  = 1._wp / rau0                                 ! 1 / volumic mass 
     169      ! 
     170      zhdiv(:,:) = 0._wp                                    ! barotropic divergence 
     171      zu_sld = 0._wp   ;   zu_asp = 0._wp                   ! tides trends (lk_tide=F) 
     172      zv_sld = 0._wp   ;   zv_asp = 0._wp 
    171173 
    172174      ! ----------------------------------------------------------------------------- 
     
    176178      !                                   !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated) 
    177179      !                                   ! -------------------------- 
    178       zua(:,:) = 0.e0   ;   zun(:,:) = 0.e0   ;   ub_b(:,:) = 0.e0 
    179       zva(:,:) = 0.e0   ;   zvn(:,:) = 0.e0   ;   vb_b(:,:) = 0.e0 
     180      zua(:,:) = 0._wp   ;   zun(:,:) = 0._wp   ;   ub_b(:,:) = 0._wp 
     181      zva(:,:) = 0._wp   ;   zvn(:,:) = 0._wp   ;   vb_b(:,:) = 0._wp 
    180182      ! 
    181183      DO jk = 1, jpkm1 
     
    205207      END DO 
    206208 
    207       !                                   !* baroclinic momentum trend (remove the vertical mean trend) 
    208       DO jk = 1, jpkm1                    ! -------------------------- 
    209          DO jj = 2, jpjm1 
    210             DO ji = fs_2, fs_jpim1   ! vector opt. 
    211                ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj) 
    212                va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj) 
    213             END DO 
    214          END DO 
    215       END DO 
    216209 
    217210      !                                   !* barotropic Coriolis trends * H (vorticity scheme dependent) 
     
    260253         DO jj = 2, jpjm1  
    261254            DO ji = fs_2, fs_jpim1   ! vector opt. 
    262                zcu(ji,jj) = zcu(ji,jj) - grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn(ji+1,jj  )    & 
    263                   &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hu(ji,jj) / e1u(ji,jj) 
    264                zcv(ji,jj) = zcv(ji,jj) - grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn(ji  ,jj+1)    & 
    265                   &                              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn(ji  ,jj  )  ) * hv(ji,jj) / e2v(ji,jj) 
     255      !     remove the rhd term according to J. Chanut's suggestion             
     256               zcu(ji,jj) = zcu(ji,jj) - grav * (  sshn(ji+1,jj  ) - sshn(ji  ,jj  )  ) * hu(ji,jj) / e1u(ji,jj) 
     257               zcv(ji,jj) = zcv(ji,jj) - grav * (  sshn(ji  ,jj+1) - sshn(ji  ,jj  )  ) * hv(ji,jj) / e2v(ji,jj) 
    266258            END DO 
    267259         END DO 
    268260      ENDIF 
     261 
     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 
    269289 
    270290      DO jj = 2, jpjm1                             ! Remove coriolis term (and possibly spg) from barotropic trend 
     
    279299      !                                             ! from the barotropic transport trend 
    280300      zcoef = -1. / z2dt_b 
     301      IF(.NOT.ln_bfrimp) THEN 
    281302# if defined key_vectopt_loop 
    282303      DO jj = 1, 1 
     
    311332         END DO 
    312333      ENDIF 
     334      END IF 
    313335 
    314336      !                                   !* d/dt(Ua), Ub, Vn (Vertical mean velocity) 
     
    410432               DO ji = fs_2, fs_jpim1   ! vector opt. 
    411433                  ! surface pressure gradient 
    412                   IF( lk_vvl) THEN 
    413                      zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
    414                         &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
    415                      zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
    416                         &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
    417                   ELSE 
    418                      zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
    419                      zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    420                   ENDIF 
     434                  zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
     435                  zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    421436                  ! energy conserving formulation for planetary vorticity term 
    422437                  zy1 = ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) 
     
    427442                  zv_cor =-z1_4 * ( ff(ji-1,jj  ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) 
    428443                  ! after velocities with implicit bottom friction 
     444 
     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)   & 
     450                     &         / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
     451                  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)   & 
     453                     &         / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
     454                  va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
     455                  ELSE 
    429456                  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)   & 
    430457                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    431458                  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)   & 
    432459                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
     460                  ENDIF 
    433461               END DO 
    434462            END DO 
     
    438466               DO ji = fs_2, fs_jpim1   ! vector opt. 
    439467                   ! surface pressure gradient 
    440                   IF( lk_vvl) THEN 
    441                      zu_spg = -grav * (  ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
    442                         &              - ( rhd(ji  ,jj ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
    443                      zv_spg = -grav * (  ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
    444                         &              - ( rhd(ji ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
    445                   ELSE 
    446                      zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
    447                      zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    448                   ENDIF 
     468                  zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
     469                  zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    449470                  ! enstrophy conserving formulation for planetary vorticity term 
    450471                  zy1 =   z1_8 * ( zwy(ji  ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj  ) ) / e1u(ji,jj) 
     
    453474                  zv_cor  = zx1 * ( ff(ji-1,jj  ) + ff(ji,jj) ) * hvr_e(ji,jj) 
    454475                  ! after velocities with implicit bottom friction 
     476                  IF(ln_bfrimp) THEN 
     477                  !   A new method to implement the implicit bottom friction.  
     478                  !   H. Liu 
     479                  !   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)   & 
     481                     &         / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
     482                  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)   & 
     484                     &         / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
     485                  va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
     486                  ELSE 
    455487                  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)   & 
    456488                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    457489                  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)   & 
    458490                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
     491                  ENDIF 
    459492               END DO 
    460493            END DO 
     
    464497               DO ji = fs_2, fs_jpim1   ! vector opt. 
    465498                  ! surface pressure gradient 
    466                   IF( lk_vvl) THEN 
    467                      zu_spg = -grav * (  ( rhd(ji+1,jj  ,1) + 1 ) * sshn_e(ji+1,jj  )    & 
    468                         &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e1u(ji,jj) 
    469                      zv_spg = -grav * (  ( rhd(ji  ,jj+1,1) + 1 ) * sshn_e(ji  ,jj+1)    & 
    470                         &              - ( rhd(ji  ,jj  ,1) + 1 ) * sshn_e(ji  ,jj  )  ) / e2v(ji,jj) 
    471                   ELSE 
    472                      zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
    473                      zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    474                   ENDIF 
     499                  zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 
     500                  zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 
    475501                  ! energy/enstrophy conserving formulation for planetary vorticity term 
    476502                  zu_cor = + z1_4 / e1u(ji,jj) * (  ftne(ji,jj  ) * zwy(ji  ,jj  ) + ftnw(ji+1,jj) * zwy(ji+1,jj  )   & 
     
    479505                     &                           + ftnw(ji,jj  ) * zwx(ji-1,jj  ) + ftne(ji,jj  ) * zwx(ji  ,jj  ) ) * hvr_e(ji,jj) 
    480506                  ! after velocities with implicit bottom friction 
     507                  IF(ln_bfrimp) THEN 
     508                  !   A new method to implement the implicit bottom friction.  
     509                  !   H. Liu 
     510                  !   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)   & 
     512                     &         / ( 1._wp      - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
     513                  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)   & 
     515                     &         / ( 1._wp      - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
     516                  va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e *   zva(ji,jj)  ) * vmask(ji,jj,1)    
     517                  ELSE 
    481518                  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)   & 
    482519                     &         / ( 1.e0         - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 
    483520                  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)   & 
    484521                     &         / ( 1.e0         - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 
     522                  ENDIF 
    485523               END DO 
    486524            END DO 
  • branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90

    r2715 r2872  
    2020   USE in_out_manager  ! I/O manager 
    2121   USE lib_mpp         ! MPP library 
     22   USE zdfbfr          ! bottom friction setup 
    2223 
    2324   IMPLICIT NONE 
     
    6162      REAL(wp), INTENT(in) ::  p2dt   ! vertical profile of tracer time-step 
    6263      !! 
    63       INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    64       REAL(wp) ::   z1_p2dt, zcoef, zzwi, zzws, zrhs   ! local scalars 
     64      INTEGER  ::   ji, jj, jk     ! dummy loop indices 
     65      INTEGER  ::   ikbum1, ikbvm1 ! local variable 
     66      REAL(wp) ::   z1_p2dt, z2dtf, zcoef, zzwi, zzws, zrhs   ! local scalars 
     67 
     68      !! * Local variables for implicit bottom friction.    H. Liu 
     69      REAL(wp) ::   zbfru, zbfrv  
     70      REAL(wp) ::   bfr_imp = 1._wp 
     71      !!---------------------------------------------------------------------- 
    6572      !!---------------------------------------------------------------------- 
    6673 
     
    7380         IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' 
    7481         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 
     82         IF(.NOT.ln_bfrimp) bfr_imp = 0._wp 
    7583      ENDIF 
    7684 
     
    8088 
    8189      ! 1. Vertical diffusion on u 
     90 
     91      ! Vertical diffusion on u&v 
    8292      ! --------------------------- 
    8393      ! Matrix and second member construction 
    84       ! bottom boundary condition: both zwi and zws must be masked as avmu can take 
    85       ! non zero value at the ocean bottom depending on the bottom friction 
    86       ! used but the bottom velocities have already been updated with the bottom 
    87       ! friction velocity in dyn_bfr using values from the previous timestep. There 
    88       ! is no need to include these in the implicit calculation. 
    89       ! 
    90       DO jk = 1, jpkm1        ! Matrix 
    91          DO jj = 2, jpjm1  
    92             DO ji = fs_2, fs_jpim1   ! vector opt. 
     94      !! bottom boundary condition: both zwi and zws must be masked as avmu can take 
     95      !! non zero value at the ocean bottom depending on the bottom friction 
     96      !! used but the bottom velocities have already been updated with the bottom 
     97      !! friction velocity in dyn_bfr using values from the previous timestep. There 
     98      !! is no need to include these in the implicit calculation. 
     99 
     100      ! The code has been modified here to implicitly implement bottom 
     101      ! friction: u(v)mask is not necessary here anymore.  
     102      ! H. Liu, April 2010. 
     103 
     104      ! 1. Vertical diffusion on u 
     105      DO jj = 2, jpjm1  
     106         DO ji = fs_2, fs_jpim1   ! vector opt. 
     107            ikbum1 = mbku(ji,jj) 
     108            ! Apply stability criteria on absolute value  : Min abs(bfr) => Max (bfr) 
     109!               zbfru = MAX( bfrua(ji,jj), fse3u(ji,jj,ikbum1)*zinv ) 
     110               zbfru = bfrua(ji,jj) 
     111 
     112            DO jk = 1, ikbum1 
    93113               zcoef = - p2dt / fse3u(ji,jj,jk) 
    94                zzwi          = zcoef * avmu (ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
    95                zwi(ji,jj,jk) = zzwi  * umask(ji,jj,jk) 
    96                zzws          = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
    97                zws(ji,jj,jk) = zzws  * umask(ji,jj,jk+1) 
    98                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 
    99             END DO 
    100          END DO 
    101       END DO 
    102       DO jj = 2, jpjm1        ! Surface boudary conditions 
    103          DO ji = fs_2, fs_jpim1   ! vector opt. 
     114               zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  ) 
     115               zws(ji,jj,jk) = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 
     116               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 
     117            END DO 
     118 
     119      ! Surface boudary conditions 
    104120            zwi(ji,jj,1) = 0._wp 
    105121            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
    106          END DO 
    107       END DO 
     122 
     123      ! Bottom boudary conditions  ! H. Liu, May, 2010 
     124!            !commented out to be consisent with v3.2, h.liu 
     125!            z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 2.0 * bfr_imp 
     126            z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 1.0_wp * bfr_imp 
     127            zws(ji,jj,ikbum1) = 0._wp 
     128            zwd(ji,jj,ikbum1) = 1._wp - zwi(ji,jj,ikbum1) - z2dtf  
    108129 
    109130      ! Matrix inversion starting from the first level 
     
    121142      !   The solution (the after velocity) is in ua 
    122143      !----------------------------------------------------------------------- 
    123       ! 
    124       DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    125          DO jj = 2, jpjm1    
    126             DO ji = fs_2, fs_jpim1   ! vector opt. 
     144 
     145      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k) 
     146            DO jk = 2, ikbum1 
    127147               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
    128148            END DO 
    129          END DO 
    130       END DO 
    131       ! 
    132       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    133          DO ji = fs_2, fs_jpim1   ! vector opt. 
    134             ua(ji,jj,1) = ub(ji,jj,1) + p2dt * (  ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) )   & 
    135                &                                                       / ( fse3u(ji,jj,1) * rau0       )  ) 
    136          END DO 
    137       END DO 
    138       DO jk = 2, jpkm1 
    139          DO jj = 2, jpjm1    
    140             DO ji = fs_2, fs_jpim1   ! vector opt. 
     149 
     150      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1 
     151            z2dtf = 0.5_wp * p2dt / ( fse3u(ji,jj,1) * rau0 ) 
     152            ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ua(ji,jj,1) + z2dtf * (utau_b(ji,jj) + utau(ji,jj)) 
     153           DO jk = 2, ikbum1    
    141154               zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk)   ! zrhs=right hand side 
    142155               ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 
    143156            END DO 
    144          END DO 
    145       END DO 
    146       ! 
    147       DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk  == 
    148          DO ji = fs_2, fs_jpim1   ! vector opt. 
    149             ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    150          END DO 
    151       END DO 
    152       DO jk = jpk-2, 1, -1 
    153          DO jj = 2, jpjm1    
    154             DO ji = fs_2, fs_jpim1   ! vector opt. 
    155                ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
     157 
     158 
     159      ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk 
     160            ua(ji,jj,ikbum1) = ua(ji,jj,ikbum1) / zwd(ji,jj,ikbum1) 
     161            DO jk = ikbum1-1, 1, -1 
     162               ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    156163            END DO 
    157164         END DO 
     
    170177      ! 2. Vertical diffusion on v 
    171178      ! --------------------------- 
    172       ! Matrix and second member construction 
    173       ! bottom boundary condition: both zwi and zws must be masked as avmv can take 
    174       ! non zero value at the ocean bottom depending on the bottom friction 
    175       ! used but the bottom velocities have already been updated with the bottom 
    176       ! friction velocity in dyn_bfr using values from the previous timestep. There 
    177       ! is no need to include these in the implicit calculation. 
    178       ! 
    179       DO jk = 1, jpkm1        ! Matrix 
     179 
     180      DO ji = fs_2, fs_jpim1   ! vector opt. 
    180181         DO jj = 2, jpjm1    
    181             DO ji = fs_2, fs_jpim1   ! vector opt. 
     182            ikbvm1 = mbkv(ji,jj) 
     183            ! Apply stability criteria on absolute value  : Min abs(bfr) => Max (bfr) 
     184!               zbfrv = MAX( bfrva(ji,jj), fse3v(ji,jj,ikbvm1)*zinv ) 
     185               zbfrv = bfrva(ji,jj) 
     186 
     187            DO jk = 1, ikbvm1 
    182188               zcoef = -p2dt / fse3v(ji,jj,jk) 
    183                zzwi          = zcoef * avmv (ji,jj,jk  ) / fse3vw(ji,jj,jk  ) 
    184                zwi(ji,jj,jk) =  zzwi * vmask(ji,jj,jk) 
    185                zzws          = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 
    186                zws(ji,jj,jk) =  zzws * vmask(ji,jj,jk+1) 
    187                zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 
    188             END DO 
    189          END DO 
    190       END DO 
    191       DO jj = 2, jpjm1        ! Surface boudary conditions 
    192          DO ji = fs_2, fs_jpim1   ! vector opt. 
     189               zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  ) 
     190               zws(ji,jj,jk) = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 
     191               zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 
     192            END DO 
     193 
     194      ! Surface boudary conditions 
    193195            zwi(ji,jj,1) = 0._wp 
    194196            zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 
    195          END DO 
    196       END DO 
     197 
     198      ! Bottom boudary conditions   
     199      ! Bottom boudary conditions  ! H. Liu, May, 2010 
     200!            !commented out to be consisent with v3.2, h.liu 
     201!            z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 2.0 * bfr_imp 
     202            z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 1.0_wp * bfr_imp 
     203            zws(ji,jj,ikbvm1) = 0._wp 
     204            zwd(ji,jj,ikbvm1) = 1._wp - zwi(ji,jj,ikbvm1) - z2dtf  
    197205 
    198206      ! Matrix inversion 
     
    206214      !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk ) 
    207215      ! 
    208       !   m is decomposed in the product of an upper and lower triangular matrix 
     216      !   m is decomposed in the product of an upper and lower triangular 
     217      !   matrix 
    209218      !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 
    210219      !   The solution (after velocity) is in 2d array va 
    211220      !----------------------------------------------------------------------- 
    212       ! 
    213       DO jk = 2, jpkm1        !==  First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k)  == 
    214          DO jj = 2, jpjm1    
    215             DO ji = fs_2, fs_jpim1   ! vector opt. 
     221 
     222      ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1   (increasing k) 
     223            DO jk = 2, ikbvm1 
    216224               zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 
    217225            END DO 
    218          END DO 
    219       END DO 
    220       ! 
    221       DO jj = 2, jpjm1        !==  second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1  == 
    222          DO ji = fs_2, fs_jpim1   ! vector opt. 
    223             va(ji,jj,1) = vb(ji,jj,1) + p2dt * (  va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) )   & 
    224                &                                                       / ( fse3v(ji,jj,1) * rau0       )  ) 
    225          END DO 
    226       END DO 
    227       DO jk = 2, jpkm1 
    228          DO jj = 2, jpjm1 
    229             DO ji = fs_2, fs_jpim1   ! vector opt. 
     226 
     227      ! second recurrence:    SOLk = RHSk - Lk / Dk-1  Lk-1 
     228            z2dtf = 0.5_wp * p2dt / ( fse3v(ji,jj,1)*rau0 ) 
     229            va(ji,jj,1) = vb(ji,jj,1) + p2dt * va(ji,jj,1) + z2dtf * (vtau_b(ji,jj) + vtau(ji,jj)) 
     230            DO jk = 2, ikbvm1 
    230231               zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk)   ! zrhs=right hand side 
    231232               va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 
    232233            END DO 
    233          END DO 
    234       END DO 
    235       ! 
    236       DO jj = 2, jpjm1        !==  thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk  == 
    237          DO ji = fs_2, fs_jpim1   ! vector opt. 
    238             va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
    239          END DO 
    240       END DO 
    241       DO jk = jpk-2, 1, -1 
    242          DO jj = 2, jpjm1    
    243             DO ji = fs_2, fs_jpim1   ! vector opt. 
    244                va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
    245             END DO 
     234 
     235      ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk 
     236            va(ji,jj,ikbvm1) = va(ji,jj,ikbvm1) / zwd(ji,jj,ikbvm1) 
     237 
     238            DO jk = ikbvm1-1, 1, -1 
     239               va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 
     240            END DO 
     241 
    246242         END DO 
    247243      END DO 
     
    258254      IF( wrk_not_released(3, 3) )   CALL ctl_stop('dyn_zdf_imp: failed to release workspace array') 
    259255      ! 
     256 
    260257   END SUBROUTINE dyn_zdf_imp 
    261258 
  • branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90

    r2715 r2872  
    3636   REAL(wp) ::   rn_bfrien = 30._wp      ! local factor to enhance coefficient bfri 
    3737   LOGICAL  ::   ln_bfr2d  = .false.     ! logical switch for 2D enhancement 
     38   LOGICAL, PUBLIC  ::   ln_bfrimp = .false.     ! logical switch for implicit bottom friction 
    3839    
    3940   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  bfrcoef2d   ! 2D bottom drag coefficient 
     
    142143      REAL(wp) ::  zfru, zfrv         !    -         - 
    143144      !! 
    144       NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfeb2, ln_bfr2d, rn_bfrien 
     145      NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfeb2, ln_bfr2d, rn_bfrien, ln_bfrimp 
    145146      !!---------------------------------------------------------------------- 
    146147 
  • branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/par_oce.F90

    r2715 r2872  
    8181   !!--------------------------------------------------------------------- 
    8282#             include "par_POMME_R025.h90" 
     83#elif defined key_amm12  
     84   !!---------------------------------------------------------------------  
     85   !!   'key_amm12'         :   Atlantic Margin Model (~12km)     :      AMM  
     86   !!---------------------------------------------------------------------  
     87#             include "par_AMM12.h90"  
     88#elif defined key_amm7  
     89   !!---------------------------------------------------------------------  
     90   !!   'key_amm7'         :   Atlantic Margin Model (~7km)     :      AMM  
     91   !!---------------------------------------------------------------------  
     92#             include "par_AMM7.h90"  
    8393#else 
    8494   !!--------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.