Changeset 6991


Ignore:
Timestamp:
2016-10-05T12:08:37+02:00 (4 years ago)
Author:
kingr
Message:

Adding Momme's bdy changes for zerograd and neumann conditions.

Location:
branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC
Files:
1 deleted
6 edited

Legend:

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

    r6331 r6991  
    5959         CASE('specified') 
    6060            CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
     61         CASE('zerograd') 
     62            CALL bdy_dyn3d_zgrad( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy ) 
    6163         CASE('zero') 
    6264            CALL bdy_dyn3d_zro( 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      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zgrad') 
     170   END SUBROUTINE bdy_dyn3d_zgrad 
    119171 
    120172   SUBROUTINE bdy_dyn3d_zro( idx, dta, kt, ib_bdy ) 
     
    303355   END SUBROUTINE bdy_dyn3d_dmp 
    304356 
     357   SUBROUTINE bdy_dyn3d_nmn( idx, ib_bdy ) 
     358      !!---------------------------------------------------------------------- 
     359      !!                 ***  SUBROUTINE bdy_dyn3d_nmn  *** 
     360      !!              
     361      !!              - Apply Neumann condition to baroclinic velocities. 
     362      !!              - Wrapper routine for bdy_nmn 
     363      !! 
     364      !! 
     365      !!---------------------------------------------------------------------- 
     366      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices 
     367      INTEGER,                      INTENT(in) ::   ib_bdy  ! BDY set index 
     368      INTEGER  ::   jb, igrd                               ! dummy loop indices 
     369      !!---------------------------------------------------------------------- 
     370      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_nmn') 
     371      ! 
     372      !! Note that at this stage the ub and ua arrays contain the baroclinic velocities. 
     373      ! 
     374      igrd = 2      ! Neumann bc on u-velocity; 
     375      !            
     376      CALL bdy_nmn( idx, igrd, ua ) 
     377      igrd = 3      ! Neumann bc on v-velocity 
     378      !  
     379      CALL bdy_nmn( idx, igrd, va ) 
     380      ! 
     381      CALL lbc_bdy_lnk( ua, 'U', -1., ib_bdy )    ! Boundary points should be updated 
     382      CALL lbc_bdy_lnk( va, 'V', -1., ib_bdy ) 
     383      ! 
     384      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_nmn') 
     385      ! 
     386   END SUBROUTINE bdy_dyn3d_nmn 
     387    
    305388#else 
    306389   !!---------------------------------------------------------------------- 
  • branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r6331 r6991  
    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/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90

    r6331 r6991  
    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) * pmask(ii,ij,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 
    356429 
    357430#else 
     
    366439      WRITE(*,*) 'bdy_orlanski_3d: You should not have seen this print! error?', kt 
    367440   END SUBROUTINE bdy_orlanski_3d 
     441   SUBROUTINE bdy_nmn( idx, igrd, phia )      ! Empty routine 
     442      WRITE(*,*) 'bdy_nmn: You should not have seen this print! error?', kt 
     443   END SUBROUTINE bdy_nmn 
    368444#endif 
    369445 
  • branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90

    r6331 r6991  
    160160      !!  
    161161      REAL(wp) ::   zwgt           ! boundary weight 
    162       INTEGER  ::   ib, ik, igrd   ! dummy loop indices 
    163       INTEGER  ::   ii, ij,zcoef, zcoef1,zcoef2, ip, jp   ! 2D addresses 
     162      REAL(wp) ::   zcoef, zcoef1,zcoef2 
     163      INTEGER  ::   ib, ik, igrd   ! dummy loop indices 
     164      INTEGER  ::   ii, ij, ip, jp   ! 2D addresses 
    164165      !!---------------------------------------------------------------------- 
    165166      ! 
     
    174175            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  ) 
    175176            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1) 
    176             IF ( zcoef1+zcoef2 == 0) THEN 
     177            IF ( NINT(zcoef1+zcoef2) == 0) THEN 
    177178               ! corner 
    178179               zcoef = tmask(ii-1,ij,ik) + tmask(ii+1,ij,ik) +  tmask(ii,ij-1,ik) +  tmask(ii,ij+1,ik) 
     
    181182                 &                    tsa(ii  ,ij-1,ik,jp_tem) * tmask(ii  ,ij-1,ik) + & 
    182183                 &                    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) 
     184               tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1._wp, zcoef) ) * tmask(ii,ij,ik) 
    184185               tsa(ii,ij,ik,jp_sal) = tsa(ii-1,ij  ,ik,jp_sal) * tmask(ii-1,ij  ,ik) + & 
    185186                 &                    tsa(ii+1,ij  ,ik,jp_sal) * tmask(ii+1,ij  ,ik) + & 
    186187                 &                    tsa(ii  ,ij-1,ik,jp_sal) * tmask(ii  ,ij-1,ik) + & 
    187188                 &                    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) 
     189               tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1._wp, zcoef) ) * tmask(ii,ij,ik) 
    189190            ELSE 
    190                ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  ) 
    191                jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1) 
    192                tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii+ip,ij+jp,ik) 
    193                tsa(ii,ij,ik,jp_sal) = tsa(ii+ip,ij+jp,ik,jp_sal) * tmask(ii+ip,ij+jp,ik) 
     191               ip = NINT(bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  )) 
     192               jp = NINT(bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1)) 
     193               tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii+ip,ij+jp,ik) * tmask(ii,ij,ik) 
     194               tsa(ii,ij,ik,jp_sal) = tsa(ii+ip,ij+jp,ik,jp_sal) * tmask(ii+ip,ij+jp,ik) * tmask(ii,ij,ik) 
    194195            ENDIF 
    195196         END DO 
  • branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

    r6332 r6991  
    189189               zztmpy = (  bfrva(ji,  jj) * vn(ji,jj  ,mbkv(ji,jj  ))  & 
    190190                      &  + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1))  )  
    191                z2d(ji,jj) = rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1)  
     191               z2d(ji,jj) = 0.5_wp * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1)  
    192192               ! 
    193193            ENDDO 
  • branches/UKMO/CO6_KD490_amm7oper/NEMOGCM/NEMO/OPA_SRC/DYN/dynkeg.F90

    r6331 r6991  
    2424   USE wrk_nemo        ! Memory Allocation 
    2525   USE timing          ! Timing 
    26  
     26#if defined key_bdy 
     27   USE bdy_oce         ! ocean open boundary conditions 
     28#endif 
     29    
    2730   IMPLICIT NONE 
    2831   PRIVATE 
     
    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 
     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  
    99135       
    100136      SELECT CASE ( kscheme )             !== Horizontal kinetic energy at T-point  ==! 
     
    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.