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

Changeset 10996


Ignore:
Timestamp:
2019-05-17T18:02:04+02:00 (6 years ago)
Author:
girrmann
Message:

trunk : fixed reproductibility issue in dynkeg.F90 with bdy, see #2283

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/OCE/DYN/dynkeg.F90

    r10425 r10996  
    7474      INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme  
    7575      ! 
    76       INTEGER  ::   ji, jj, jk, jb    ! dummy loop indices 
    77       INTEGER  ::   ii, ifu, ib_bdy   ! local integers 
    78       INTEGER  ::   ij, ifv, igrd     !   -       - 
    79       REAL(wp) ::   zu, zv            ! local scalars 
     76      INTEGER  ::   ji, jj, jk, jb           ! dummy loop indices 
     77      INTEGER  ::   ifu, ifv, igrd, ib_bdy   ! local integers 
     78      REAL(wp) ::   zu, zv                   ! local scalars 
    8079      REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zhke 
    8180      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv  
     81      REAL(wp)  :: zweightu, zweightv 
    8282      !!---------------------------------------------------------------------- 
    8383      ! 
     
    9797       
    9898      zhke(:,:,jpk) = 0._wp 
    99        
    100       IF (ln_bdy) THEN 
    101          ! Maria Luneva & Fred Wobus: July-2016 
    102          ! compensate for lack of turbulent kinetic energy on liquid bdy points 
    103          DO ib_bdy = 1, nb_bdy 
    104             IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 
    105                igrd = 2           ! Copying normal velocity into points outside bdy 
    106                DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    107                   DO jk = 1, jpkm1 
    108                      ii   = idx_bdy(ib_bdy)%nbi(jb,igrd) 
    109                      ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
    110                      ifu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) ) 
    111                      un(ii-ifu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk) 
    112                   END DO 
    113                END DO 
    114                ! 
    115                igrd = 3           ! Copying normal velocity into points outside bdy 
    116                DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    117                   DO jk = 1, jpkm1 
    118                      ii   = idx_bdy(ib_bdy)%nbi(jb,igrd) 
    119                      ij   = idx_bdy(ib_bdy)%nbj(jb,igrd) 
    120                      ifv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) ) 
    121                      vn(ii,ij-ifv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk) 
    122                   END DO 
    123                END DO 
    124             ENDIF 
    125          ENDDO   
    126       ENDIF  
    12799 
    128100      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==! 
     
    140112            END DO 
    141113         END DO 
     114         ! 
     115         IF (ln_bdy) THEN 
     116            ! Maria Luneva & Fred Wobus: July-2016 
     117            ! compensate for lack of turbulent kinetic energy on liquid bdy points 
     118            DO ib_bdy = 1, nb_bdy 
     119               IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 
     120                  igrd = 1           ! compensating null velocity on the bdy 
     121                  DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     122                     ji   = idx_bdy(ib_bdy)%nbi(jb,igrd)   ! maximum extent : from 2 to jpi-1 
     123                     jj   = idx_bdy(ib_bdy)%nbj(jb,igrd)   ! maximum extent : from 2 to jpj-1 
     124                     DO jk = 1, jpkm1 
     125                        zhke(ji,jj,jk) = 0._wp 
     126                        zweightu = umask(ji-1,jj  ,jk) + umask(ji,jj,jk) 
     127                        zweightv = vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk) 
     128                        zu = un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)  +  un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) 
     129                        zv = vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  +  vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) 
     130                        IF( zweightu > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) + zu / (2._wp * zweightu)  
     131                        IF( zweightv > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) + zv / (2._wp * zweightv)  
     132                     END DO 
     133                  END DO 
     134               END IF 
     135               CALL lbc_bdy_lnk( 'dynkeg', zhke, 'T', 1., ib_bdy )   ! send 2 and recv jpi, jpj used in the computation of the speed tendencies 
     136            END DO 
     137         END IF 
    142138         ! 
    143139      CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --! 
     
    158154            END DO 
    159155         END DO 
     156         IF (ln_bdy) THEN 
     157            ! Maria Luneva & Fred Wobus: July-2016 
     158            ! compensate for lack of turbulent kinetic energy on liquid bdy points 
     159            DO ib_bdy = 1, nb_bdy 
     160               IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN 
     161                  igrd = 1           ! compensation null velocity on land at the bdy 
     162                  DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     163                     ji   = idx_bdy(ib_bdy)%nbi(jb,igrd)   ! maximum extent : from 2 to jpi-1 
     164                     jj   = idx_bdy(ib_bdy)%nbj(jb,igrd)   ! maximum extent : from 2 to jpj-1 
     165                     DO jk = 1, jpkm1 
     166                        zhke(ji,jj,jk) = 0._wp 
     167                        zweightu = 8._wp * ( umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk) ) & 
     168                             &   + 2._wp * ( umask(ji-1,jj-1,jk) + umask(ji-1,jj+1,jk) + umask(ji  ,jj-1,jk) + umask(ji  ,jj+1,jk) )  
     169                        zweightv = 8._wp * ( vmask(ji  ,jj-1,jk) + vmask(ji  ,jj-1,jk) ) & 
     170                             &   + 2._wp * ( vmask(ji-1,jj-1,jk) + vmask(ji+1,jj-1,jk) + vmask(ji-1,jj  ,jk) + vmask(ji+1,jj  ,jk) ) 
     171                        zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    & 
     172                           &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  & 
     173                           &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   & 
     174                           &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) 
     175                        zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    & 
     176                           &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  & 
     177                           &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   & 
     178                           &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) 
     179                        IF( zweightu > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) +  zu / ( 2._wp * zweightu ) 
     180                        IF( zweightv > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) +  zv / ( 2._wp * zweightv ) 
     181                     END DO 
     182                  END DO 
     183               END IF 
     184            END DO 
     185         END IF 
    160186         CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. ) 
    161187         ! 
    162       END SELECT 
    163  
    164       IF (ln_bdy) THEN 
    165          ! restore velocity masks at points outside boundary 
    166          un(:,:,:) = un(:,:,:) * umask(:,:,:) 
    167          vn(:,:,:) = vn(:,:,:) * vmask(:,:,:) 
    168       ENDIF       
    169  
     188      END SELECT  
    170189      ! 
    171190      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==! 
Note: See TracChangeset for help on using the changeset viewer.