Changeset 7087


Ignore:
Timestamp:
2016-10-25T15:46:01+02:00 (4 years ago)
Author:
jgraham
Message:

Bring in options for 3D U/V boundaries, following changes made for CO6.
Options added are 'zerograd' or 'neumann' boundary conditions.

Location:
branches/UKMO/AMM15_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/AMM15_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90

    r6928 r7087  
    6161         CASE('zero') 
    6262            CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
     63         CASE('zerograd')  
     64            CALL bdy_dyn3d_zgrad( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )        
     65         CASE('neumann')  
     66            CALL bdy_dyn3d_nmn( idx_bdy(ib_bdy), ib_bdy ) 
    6367         CASE('orlanski') 
    6468            CALL bdy_dyn3d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. ) 
     
    117121 
    118122   END SUBROUTINE bdy_dyn3d_spe 
     123 
     124   SUBROUTINE bdy_dyn3d_zgrad( idx, dta, kt , ib_bdy )  
     125      !!----------------------------------------------------------------------  
     126      !!                  ***  SUBROUTINE bdy_dyn3d_zgrad  ***  
     127      !!  
     128      !! ** Purpose : - Enforce a zero gradient of normal velocity  
     129      !!  
     130      !!----------------------------------------------------------------------  
     131      INTEGER                     ::   kt  
     132      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices  
     133      TYPE(OBC_DATA),  INTENT(in) ::   dta  ! OBC external data  
     134      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index  
     135      !!  
     136      INTEGER  ::   jb, jk         ! dummy loop indices  
     137      INTEGER  ::   ii, ij, igrd   ! local integers  
     138      REAL(wp) ::   zwgt           ! boundary weight  
     139      INTEGER  ::   fu, fv  
     140      !!----------------------------------------------------------------------  
     141      !  
     142      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_zgrad')  
     143      !  
     144      igrd = 2                      ! Copying tangential velocity into bdy points  
     145      DO jb = 1, idx%nblenrim(igrd)  
     146         DO jk = 1, jpkm1  
     147            ii   = idx%nbi(jb,igrd)  
     148            ij   = idx%nbj(jb,igrd)  
     149            fu   = ABS( ABS (NINT( idx%flagu(jb,igrd) ) ) - 1 )  
     150            ua(ii,ij,jk) = ua(ii,ij,jk) * REAL( 1 - fu) + ( ua(ii,ij+fu,jk) * umask(ii,ij+fu,jk) &  
     151                        &+ ua(ii,ij-fu,jk) * umask(ii,ij-fu,jk) ) * umask(ii,ij,jk) * REAL( fu )  
     152         END DO  
     153      END DO  
     154      !  
     155      igrd = 3                      ! Copying tangential velocity into bdy points  
     156      DO jb = 1, idx%nblenrim(igrd)  
     157         DO jk = 1, jpkm1  
     158            ii   = idx%nbi(jb,igrd)  
     159            ij   = idx%nbj(jb,igrd)  
     160            fv   = ABS( ABS (NINT( idx%flagv(jb,igrd) ) ) - 1 )  
     161            va(ii,ij,jk) = va(ii,ij,jk) * REAL( 1 - fv ) + ( va(ii+fv,ij,jk) * vmask(ii+fv,ij,jk) &  
     162                        &+ va(ii-fv,ij,jk) * vmask(ii-fv,ij,jk) ) * vmask(ii,ij,jk) * REAL( fv )  
     163         END DO  
     164      END DO  
     165      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )   ! Boundary points should be updated    
     166      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )     
     167      !  
     168      IF( kt .eq. nit000 ) CLOSE( unit = 102 )  
     169  
     170      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zgrad')  
     171  
     172   END SUBROUTINE bdy_dyn3d_zgrad  
    119173 
    120174   SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy ) 
     
    303357   END SUBROUTINE bdy_dyn3d_dmp 
    304358 
     359   SUBROUTINE bdy_dyn3d_nmn( idx, ib_bdy )  
     360      !!----------------------------------------------------------------------  
     361      !!                 ***  SUBROUTINE bdy_dyn3d_nmn  ***  
     362      !!               
     363      !!              - Apply Neumann condition to baroclinic velocities.   
     364      !!              - Wrapper routine for bdy_nmn  
     365      !!   
     366      !!  
     367      !!----------------------------------------------------------------------  
     368      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices  
     369      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index  
     370  
     371      INTEGER  ::   jb, igrd                               ! dummy loop indices  
     372      !!----------------------------------------------------------------------  
     373  
     374      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_nmn')  
     375      !  
     376      !! Note that at this stage the ub and ua arrays contain the baroclinic velocities.   
     377      !  
     378      igrd = 2      ! Neumann bc on u-velocity;   
     379      !              
     380      CALL bdy_nmn( idx, igrd, ua )  
     381  
     382      igrd = 3      ! Neumann bc on v-velocity  
     383      !    
     384      CALL bdy_nmn( idx, igrd, va )  
     385      !  
     386      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )    ! Boundary points should be updated  
     387      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy )  
     388      !  
     389      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_nmn')  
     390      !  
     391   END SUBROUTINE bdy_dyn3d_nmn  
     392  
     393 
    305394#else 
    306395   !!---------------------------------------------------------------------- 
  • branches/UKMO/AMM15_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r7062 r7087  
    213213             dta_bdy(ib_bdy)%ll_u3d = .true. 
    214214             dta_bdy(ib_bdy)%ll_v3d = .true. 
     215          CASE('neumann')  
     216             IF(lwp) WRITE(numout,*) '      Neumann conditions'  
     217             dta_bdy(ib_bdy)%ll_u3d = .false.  
     218             dta_bdy(ib_bdy)%ll_v3d = .false.  
     219          CASE('zerograd')  
     220             IF(lwp) WRITE(numout,*) '      Zero gradient for baroclinic velocities'  
     221             dta_bdy(ib_bdy)%ll_u3d = .false.  
     222             dta_bdy(ib_bdy)%ll_v3d = .false. 
    215223          CASE('zero') 
    216224             IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)' 
  • branches/UKMO/AMM15_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90

    r6928 r7087  
    2626   PUBLIC   bdy_orlanski_2d     ! routine called where? 
    2727   PUBLIC   bdy_orlanski_3d     ! routine called where? 
     28   PUBLIC   bdy_nmn     ! routine called where?  
    2829 
    2930   !!---------------------------------------------------------------------- 
     
    354355   END SUBROUTINE bdy_orlanski_3d 
    355356 
     357   SUBROUTINE bdy_nmn( idx, igrd, phia )  
     358      !!----------------------------------------------------------------------  
     359      !!                 ***  SUBROUTINE bdy_nmn  ***  
     360      !!                      
     361      !! ** Purpose : Duplicate the value at open boundaries, zero gradient.  
     362      !!   
     363      !!----------------------------------------------------------------------  
     364      INTEGER,                    INTENT(in)     ::   igrd     ! grid index  
     365      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   phia     ! model after 3D field (to be updated)  
     366      TYPE(OBC_INDEX), INTENT(in) ::   idx  ! OBC indices  
     367      !!   
     368      REAL(wp) ::   zcoef, zcoef1, zcoef2  
     369      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask      ! land/sea mask for field  
     370      REAL(wp), POINTER, DIMENSION(:,:)        :: bdypmask      ! land/sea mask for field  
     371      INTEGER  ::   ib, ik   ! dummy loop indices  
     372      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses  
     373      !!----------------------------------------------------------------------  
     374      !  
     375      IF( nn_timing == 1 ) CALL timing_start('bdy_nmn')  
     376      !  
     377      SELECT CASE(igrd)  
     378         CASE(1)  
     379            pmask => tmask(:,:,:)  
     380            bdypmask => bdytmask(:,:)  
     381         CASE(2)  
     382            pmask => umask(:,:,:)  
     383            bdypmask => bdyumask(:,:)  
     384         CASE(3)  
     385            pmask => vmask(:,:,:)  
     386            bdypmask => bdyvmask(:,:)  
     387         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_nmn' )  
     388      END SELECT  
     389      DO ib = 1, idx%nblenrim(igrd)  
     390         ii = idx%nbi(ib,igrd)  
     391         ij = idx%nbj(ib,igrd)  
     392         DO ik = 1, jpkm1  
     393            ! search the sense of the gradient  
     394            zcoef1 = bdypmask(ii-1,ij  )*pmask(ii-1,ij,ik) +  bdypmask(ii+1,ij  )*pmask(ii+1,ij,ik)  
     395            zcoef2 = bdypmask(ii  ,ij-1)*pmask(ii,ij-1,ik) +  bdypmask(ii  ,ij+1)*pmask(ii,ij+1,ik)  
     396            IF ( nint(zcoef1+zcoef2) == 0) THEN  
     397               ! corner **** we probably only want to set the tangentail component for the dynamics here  
     398               zcoef = pmask(ii-1,ij,ik) + pmask(ii+1,ij,ik) +  pmask(ii,ij-1,ik) +  pmask(ii,ij+1,ik)  
     399               IF (zcoef > .5_wp) THEN ! Only set none isolated points.  
     400                 phia(ii,ij,ik) = phia(ii-1,ij  ,ik) * pmask(ii-1,ij  ,ik) + &  
     401                   &              phia(ii+1,ij  ,ik) * pmask(ii+1,ij  ,ik) + &  
     402                   &              phia(ii  ,ij-1,ik) * pmask(ii  ,ij-1,ik) + &  
     403                   &              phia(ii  ,ij+1,ik) * pmask(ii  ,ij+1,ik)  
     404                 phia(ii,ij,ik) = ( phia(ii,ij,ik) / zcoef ) * pmask(ii,ij,ik)  
     405               ELSE  
     406                 phia(ii,ij,ik) = phia(ii,ij  ,ik) * pmask(ii,ij  ,ik)  
     407               ENDIF  
     408            ELSEIF ( nint(zcoef1+zcoef2) == 2) THEN  
     409               ! oblique corner **** we probably only want to set the normal component for the dynamics here  
     410               zcoef = pmask(ii-1,ij,ik)*bdypmask(ii-1,ij  ) + pmask(ii+1,ij,ik)*bdypmask(ii+1,ij  ) + &  
     411                   &   pmask(ii,ij-1,ik)*bdypmask(ii,ij -1 ) +  pmask(ii,ij+1,ik)*bdypmask(ii,ij+1  )  
     412               phia(ii,ij,ik) = phia(ii-1,ij  ,ik) * pmask(ii-1,ij  ,ik)*bdypmask(ii-1,ij  ) + &  
     413                   &            phia(ii+1,ij  ,ik) * pmask(ii+1,ij  ,ik)*bdypmask(ii+1,ij  )  + &  
     414                   &            phia(ii  ,ij-1,ik) * pmask(ii  ,ij-1,ik)*bdypmask(ii,ij -1 ) + &  
     415                   &            phia(ii  ,ij+1,ik) * pmask(ii  ,ij+1,ik)*bdypmask(ii,ij+1  )  
     416    
     417               phia(ii,ij,ik) = ( phia(ii,ij,ik) / MAX(1._wp, zcoef) ) * pmask(ii,ij,ik)  
     418            ELSE  
     419               ip = nint(bdypmask(ii+1,ij  )*pmask(ii+1,ij,ik) - bdypmask(ii-1,ij  )*pmask(ii-1,ij,ik))  
     420               jp = nint(bdypmask(ii  ,ij+1)*pmask(ii,ij+1,ik) - bdypmask(ii  ,ij-1)*pmask(ii,ij-1,ik))  
     421               phia(ii,ij,ik) = phia(ii+ip,ij+jp,ik) * pmask(ii+ip,ij+jp,ik)  
     422            ENDIF  
     423         END DO  
     424      END DO  
     425      !  
     426      IF( nn_timing == 1 ) CALL timing_stop('bdy_nmn')  
     427      !  
     428   END SUBROUTINE bdy_nmn  
     429 
    356430 
    357431#else 
     
    366440      WRITE(*,*) 'bdy_orlanski_3d: You should not have seen this print! error?', kt 
    367441   END SUBROUTINE bdy_orlanski_3d 
     442   SUBROUTINE bdy_nmn( idx, igrd, phia )      ! Empty routine  
     443      WRITE(*,*) 'bdy_nmn: You should not have seen this print! error?', kt  
     444   END SUBROUTINE bdy_nmn  
    368445#endif 
    369446 
  • branches/UKMO/AMM15_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90

    r6928 r7087  
    9191      !!  
    9292      REAL(wp) ::   zwgt           ! boundary weight 
    93       INTEGER  ::   ib, ik, igrd   ! dummy loop indices 
    94       INTEGER  ::   ii, ij         ! 2D addresses 
     93      REAL(wp) ::   zcoef, zcoef1,zcoef2  
     94      INTEGER  ::   ib, ik, igrd   ! dummy loop indices  
     95      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses  
    9596      !!---------------------------------------------------------------------- 
    9697      ! 
     
    160161      !!  
    161162      REAL(wp) ::   zwgt           ! boundary weight 
    162       INTEGER  ::   ib, ik, igrd   ! dummy loop indices 
    163       INTEGER  ::   ii, ij,zcoef, zcoef1,zcoef2, ip, jp   ! 2D addresses 
     163      REAL(wp) ::   zcoef, zcoef1,zcoef2  
     164      INTEGER  ::   ib, ik, igrd   ! dummy loop indices  
     165      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses  
    164166      !!---------------------------------------------------------------------- 
    165167      ! 
     
    174176            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  ) 
    175177            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1) 
    176             IF ( zcoef1+zcoef2 == 0) THEN 
     178            IF ( NINT(zcoef1+zcoef2) == 0) THEN  
    177179               ! corner 
    178180               zcoef = tmask(ii-1,ij,ik) + tmask(ii+1,ij,ik) +  tmask(ii,ij-1,ik) +  tmask(ii,ij+1,ik) 
     
    181183                 &                    tsa(ii  ,ij-1,ik,jp_tem) * tmask(ii  ,ij-1,ik) + & 
    182184                 &                    tsa(ii  ,ij+1,ik,jp_tem) * tmask(ii  ,ij+1,ik) 
    183                tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1, zcoef) ) * tmask(ii,ij,ik) 
     185               tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1._wp, zcoef) ) * tmask(ii,ij,ik)  
    184186               tsa(ii,ij,ik,jp_sal) = tsa(ii-1,ij  ,ik,jp_sal) * tmask(ii-1,ij  ,ik) + & 
    185187                 &                    tsa(ii+1,ij  ,ik,jp_sal) * tmask(ii+1,ij  ,ik) + & 
    186188                 &                    tsa(ii  ,ij-1,ik,jp_sal) * tmask(ii  ,ij-1,ik) + & 
    187189                 &                    tsa(ii  ,ij+1,ik,jp_sal) * tmask(ii  ,ij+1,ik) 
    188                tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1, zcoef) ) * tmask(ii,ij,ik) 
     190               tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1._wp, zcoef) ) * tmask(ii,ij,ik) 
    189191            ELSE 
    190                ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  ) 
    191                jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1) 
     192               ip = NINT(bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  ))  
     193               jp = NINT(bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1))  
    192194               tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii+ip,ij+jp,ik) 
    193195               tsa(ii,ij,ik,jp_sal) = tsa(ii+ip,ij+jp,ik,jp_sal) * tmask(ii+ip,ij+jp,ik) 
  • branches/UKMO/AMM15_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90

    r6928 r7087  
    2424   USE wrk_nemo        ! Memory Allocation 
    2525   USE timing          ! Timing 
     26#if defined key_bdy  
     27   USE bdy_oce         ! ocean open boundary conditions  
     28#endif  
    2629 
    2730   IMPLICIT NONE 
     
    7881      REAL(wp), POINTER, DIMENSION(:,:,:) :: zhke 
    7982      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv  
     83#if defined key_bdy  
     84      INTEGER  ::   jb                 ! dummy loop indices  
     85      INTEGER  ::   ii, ij, igrd, ib_bdy   ! local integers  
     86      INTEGER  ::   fu, fv  
     87#endif  
    8088      !!---------------------------------------------------------------------- 
    8189      ! 
     
    97105       
    98106      zhke(:,:,jpk) = 0._wp 
    99        
     107 
     108#if defined key_bdy  
     109      ! Maria Luneva & Fred Wobus: July-2016  
     110      ! compensate for lack of turbulent kinetic energy on liquid bdy points  
     111      DO ib_bdy = 1, nb_bdy  
     112         IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN  
     113            igrd = 2           ! Copying normal velocity into points outside bdy  
     114            DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)  
     115               DO jk = 1, jpkm1  
     116                  ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)  
     117                  ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)  
     118                  fu   = NINT( idx_bdy(ib_bdy)%flagu(jb,igrd) )  
     119                  un(ii-fu,ij,jk) = un(ii,ij,jk) * umask(ii,ij,jk)  
     120               END DO  
     121            END DO  
     122            !  
     123            igrd = 3           ! Copying normal velocity into points outside bdy  
     124            DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)  
     125               DO jk = 1, jpkm1  
     126                  ii   = idx_bdy(ib_bdy)%nbi(jb,igrd)  
     127                  ij   = idx_bdy(ib_bdy)%nbj(jb,igrd)  
     128                  fv   = NINT( idx_bdy(ib_bdy)%flagv(jb,igrd) )  
     129                  vn(ii,ij-fv,jk) = vn(ii,ij,jk) * vmask(ii,ij,jk)  
     130               END DO  
     131            END DO  
     132         ENDIF  
     133      ENDDO   
     134#endif        
     135 
    100136      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==! 
    101137      ! 
     
    134170      END SELECT 
    135171      ! 
     172 
     173#if defined key_bdy  
     174      ! restore velocity masks at points outside boundary  
     175      un(:,:,:) = un(:,:,:) * umask(:,:,:)  
     176      vn(:,:,:) = vn(:,:,:) * vmask(:,:,:)  
     177#endif   
     178 
    136179      DO jk = 1, jpkm1                    !==  grad( KE ) added to the general momentum trends  ==! 
    137180         DO jj = 2, jpjm1 
Note: See TracChangeset for help on using the changeset viewer.