Ignore:
Timestamp:
2013-11-20T17:28:04+01:00 (7 years ago)
Author:
cetlod
Message:

dev_MERGE_2013 : 1st step of the merge, see ticket #1185

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90

    r3680 r4292  
    66   !! History :  3.4  !  2011     (D. Storkey) new module as part of BDY rewrite 
    77   !!            3.5  !  2012     (S. Mocavero, I. Epicoco) Optimization of BDY communications 
     8   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_bdy  
     
    1112   !!   'key_bdy' :                    Unstructured Open Boundary Condition 
    1213   !!---------------------------------------------------------------------- 
    13    !!   bdy_dyn2d      : Apply open boundary conditions to barotropic variables. 
    14    !!   bdy_dyn2d_fla    : Apply Flather condition 
     14   !!   bdy_dyn2d          : Apply open boundary conditions to barotropic variables. 
     15   !!   bdy_dyn2d_frs      : Apply Flow Relaxation Scheme  
     16   !!   bdy_dyn2d_fla      : Apply Flather condition 
     17   !!   bdy_dyn2d_orlanski : Orlanski Radiation 
     18   !!   bdy_ssh            : Duplicate sea level across open boundaries 
    1519   !!---------------------------------------------------------------------- 
    1620   USE timing          ! Timing 
     
    1822   USE dom_oce         ! ocean space and time domain 
    1923   USE bdy_oce         ! ocean open boundary conditions 
     24   USE bdylib          ! BDY library routines 
    2025   USE dynspg_oce      ! for barotropic variables 
    2126   USE phycst          ! physical constants 
     
    2631   PRIVATE 
    2732 
    28    PUBLIC   bdy_dyn2d     ! routine called in dynspg_ts and bdy_dyn 
     33   PUBLIC   bdy_dyn2d   ! routine called in dynspg_ts and bdy_dyn 
     34   PUBLIC   bdy_ssh       ! routine called in dynspg_ts or sshwzv 
    2935 
    3036   !!---------------------------------------------------------------------- 
     
    4854      DO ib_bdy=1, nb_bdy 
    4955 
    50          SELECT CASE( nn_dyn2d(ib_bdy) ) 
    51          CASE(jp_none) 
     56         SELECT CASE( cn_dyn2d(ib_bdy) ) 
     57         CASE('none') 
    5258            CYCLE 
    53          CASE(jp_frs) 
     59         CASE('frs') 
    5460            CALL bdy_dyn2d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy ) 
    55          CASE(jp_flather) 
     61         CASE('flather') 
    5662            CALL bdy_dyn2d_fla( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy ) 
     63         CASE('orlanski') 
     64            CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.false. ) 
     65         CASE('orlanski_npo') 
     66            CALL bdy_dyn2d_orlanski( idx_bdy(ib_bdy), dta_bdy(ib_bdy), ib_bdy, ll_npo=.true. ) 
    5767         CASE DEFAULT 
    5868            CALL ctl_stop( 'bdy_dyn2d : unrecognised option for open boundaries for barotropic variables' ) 
     
    8999         ij   = idx%nbj(jb,igrd) 
    90100         zwgt = idx%nbw(jb,igrd) 
    91          pu2d(ii,ij) = ( pu2d(ii,ij) + zwgt * ( dta%u2d(jb) - pu2d(ii,ij) ) ) * umask(ii,ij,1) 
     101         pua2d(ii,ij) = ( pua2d(ii,ij) + zwgt * ( dta%u2d(jb) - pua2d(ii,ij) ) ) * umask(ii,ij,1) 
    92102      END DO 
    93103      ! 
     
    97107         ij   = idx%nbj(jb,igrd) 
    98108         zwgt = idx%nbw(jb,igrd) 
    99          pv2d(ii,ij) = ( pv2d(ii,ij) + zwgt * ( dta%v2d(jb) - pv2d(ii,ij) ) ) * vmask(ii,ij,1) 
     109         pva2d(ii,ij) = ( pva2d(ii,ij) + zwgt * ( dta%v2d(jb) - pva2d(ii,ij) ) ) * vmask(ii,ij,1) 
    100110      END DO  
    101       CALL lbc_bdy_lnk( pu2d, 'U', -1., ib_bdy )  
    102       CALL lbc_bdy_lnk( pv2d, 'V', -1., ib_bdy)   ! Boundary points should be updated 
     111      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy )  
     112      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy)   ! Boundary points should be updated 
    103113      ! 
    104114      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_frs') 
     
    133143      INTEGER  ::   jb, igrd                         ! dummy loop indices 
    134144      INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses 
     145      REAL(wp), POINTER :: flagu, flagv              ! short cuts 
    135146      REAL(wp) ::   zcorr                            ! Flather correction 
    136147      REAL(wp) ::   zforc                            ! temporary scalar 
     148      REAL(wp) ::   zflag, z1_2                      !    "        " 
    137149      !!---------------------------------------------------------------------- 
    138150 
    139151      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_fla') 
     152 
     153      z1_2 = 0.5_wp 
    140154 
    141155      ! ---------------------------------! 
     
    160174         ii  = idx%nbi(jb,igrd) 
    161175         ij  = idx%nbj(jb,igrd)  
    162          iim1 = ii + MAX( 0, INT( idx%flagu(jb) ) )   ! T pts i-indice inside the boundary 
    163          iip1 = ii - MIN( 0, INT( idx%flagu(jb) ) )   ! T pts i-indice outside the boundary  
     176         flagu => idx%flagu(jb,igrd) 
     177         iim1 = ii + MAX( 0, INT( flagu ) )   ! T pts i-indice inside the boundary 
     178         iip1 = ii - MIN( 0, INT( flagu ) )   ! T pts i-indice outside the boundary  
    164179         ! 
    165          zcorr = - idx%flagu(jb) * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) ) 
    166          zforc = dta%u2d(jb) 
    167          pu2d(ii,ij) = zforc + zcorr * umask(ii,ij,1)  
     180         zcorr = - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) ) 
     181 
     182         ! jchanut tschanges: Set zflag to 0 below to revert to Flather scheme 
     183         ! Use characteristics method instead 
     184         zflag = ABS(flagu) 
     185         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(iim1,ij) 
     186         pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1)  
    168187      END DO 
    169188      ! 
     
    173192         ii  = idx%nbi(jb,igrd) 
    174193         ij  = idx%nbj(jb,igrd)  
    175          ijm1 = ij + MAX( 0, INT( idx%flagv(jb) ) )   ! T pts j-indice inside the boundary 
    176          ijp1 = ij - MIN( 0, INT( idx%flagv(jb) ) )   ! T pts j-indice outside the boundary  
     194         flagv => idx%flagv(jb,igrd) 
     195         ijm1 = ij + MAX( 0, INT( flagv ) )   ! T pts j-indice inside the boundary 
     196         ijp1 = ij - MIN( 0, INT( flagv ) )   ! T pts j-indice outside the boundary  
    177197         ! 
    178          zcorr = - idx%flagv(jb) * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) ) 
    179          zforc = dta%v2d(jb) 
    180          pv2d(ii,ij) = zforc + zcorr * vmask(ii,ij,1) 
    181       END DO 
    182       CALL lbc_bdy_lnk( pu2d, 'U', -1., ib_bdy )   ! Boundary points should be updated 
    183       CALL lbc_bdy_lnk( pv2d, 'V', -1., ib_bdy )   ! 
     198         zcorr = - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) ) 
     199 
     200         ! jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme 
     201         ! Use characteristics method instead 
     202         zflag = ABS(flagv) 
     203         zforc  = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ijm1) 
     204         pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1) 
     205      END DO 
     206      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated 
     207      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy )   ! 
    184208      ! 
    185209      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_fla') 
    186210      ! 
    187211   END SUBROUTINE bdy_dyn2d_fla 
     212 
     213 
     214   SUBROUTINE bdy_dyn2d_orlanski( idx, dta, ib_bdy, ll_npo ) 
     215      !!---------------------------------------------------------------------- 
     216      !!                 ***  SUBROUTINE bdy_dyn2d_orlanski  *** 
     217      !!              
     218      !!              - Apply Orlanski radiation condition adaptively: 
     219      !!                  - radiation plus weak nudging at outflow points 
     220      !!                  - no radiation and strong nudging at inflow points 
     221      !!  
     222      !! 
     223      !! References:  Marchesiello, McWilliams and Shchepetkin, Ocean Modelling vol. 3 (2001)     
     224      !!---------------------------------------------------------------------- 
     225      TYPE(OBC_INDEX),              INTENT(in) ::   idx  ! OBC indices 
     226      TYPE(OBC_DATA),               INTENT(in) ::   dta  ! OBC external data 
     227      INTEGER,                      INTENT(in) ::   ib_bdy  ! number of current open boundary set 
     228      LOGICAL,                      INTENT(in) ::   ll_npo  ! flag for NPO version 
     229 
     230      INTEGER  ::   ib, igrd                               ! dummy loop indices 
     231      INTEGER  ::   ii, ij, iibm1, ijbm1                   ! indices 
     232      !!---------------------------------------------------------------------- 
     233 
     234      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_orlanski') 
     235      ! 
     236      igrd = 2      ! Orlanski bc on u-velocity;  
     237      !             
     238      CALL bdy_orlanski_2d( idx, igrd, pub2d, pua2d, dta%u2d, ll_npo ) 
     239 
     240      igrd = 3      ! Orlanski bc on v-velocity 
     241      !   
     242      CALL bdy_orlanski_2d( idx, igrd, pvb2d, pva2d, dta%v2d, ll_npo ) 
     243      ! 
     244      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_orlanski') 
     245      ! 
     246      CALL lbc_bdy_lnk( pua2d, 'U', -1., ib_bdy )   ! Boundary points should be updated 
     247      CALL lbc_bdy_lnk( pva2d, 'V', -1., ib_bdy )   ! 
     248      ! 
     249      IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn2d_orlanski') 
     250      ! 
     251   END SUBROUTINE bdy_dyn2d_orlanski 
     252 
     253   SUBROUTINE bdy_ssh( zssh ) 
     254      !!---------------------------------------------------------------------- 
     255      !!                  ***  SUBROUTINE bdy_ssh  *** 
     256      !! 
     257      !! ** Purpose : Duplicate sea level across open boundaries 
     258      !! 
     259      !!---------------------------------------------------------------------- 
     260      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zssh ! Sea level 
     261      !! 
     262      INTEGER  ::   ib_bdy, ib, igrd                        ! local integers 
     263      INTEGER  ::   ii, ij, zcoef, zcoef1, zcoef2, ip, jp   !   "       " 
     264 
     265      igrd = 1                       ! Everything is at T-points here 
     266 
     267      DO ib_bdy = 1, nb_bdy 
     268         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     269            ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     270            ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     271            ! Set gradient direction: 
     272            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  ) 
     273            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1) 
     274            IF ( zcoef1+zcoef2 == 0 ) THEN 
     275               ! corner 
     276!               zcoef = tmask(ii-1,ij,1) + tmask(ii+1,ij,1) +  tmask(ii,ij-1,1) +  tmask(ii,ij+1,1) 
     277!               zssh(ii,ij) = zssh(ii-1,ij  ) * tmask(ii-1,ij  ,1) + & 
     278!                 &           zssh(ii+1,ij  ) * tmask(ii+1,ij  ,1) + & 
     279!                 &           zssh(ii  ,ij-1) * tmask(ii  ,ij-1,1) + & 
     280!                 &           zssh(ii  ,ij+1) * tmask(ii  ,ij+1,1) 
     281               zcoef = bdytmask(ii-1,ij) + bdytmask(ii+1,ij) +  bdytmask(ii,ij-1) +  bdytmask(ii,ij+1) 
     282               zssh(ii,ij) = zssh(ii-1,ij  ) * bdytmask(ii-1,ij  ) + & 
     283                 &           zssh(ii+1,ij  ) * bdytmask(ii+1,ij  ) + & 
     284                 &           zssh(ii  ,ij-1) * bdytmask(ii  ,ij-1) + & 
     285                 &           zssh(ii  ,ij+1) * bdytmask(ii  ,ij+1) 
     286               zssh(ii,ij) = ( zssh(ii,ij) / MAX( 1, zcoef) ) * tmask(ii,ij,1) 
     287            ELSE 
     288               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  ) 
     289               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1) 
     290               zssh(ii,ij) = zssh(ii+ip,ij+jp) * tmask(ii+ip,ij+jp,1) 
     291            ENDIF 
     292         END DO 
     293 
     294         ! Boundary points should be updated 
     295         CALL lbc_bdy_lnk( zssh(:,:), 'T', 1., ib_bdy ) 
     296      END DO 
     297 
     298   END SUBROUTINE bdy_ssh 
     299 
    188300#else 
    189301   !!---------------------------------------------------------------------- 
     
    192304CONTAINS 
    193305   SUBROUTINE bdy_dyn2d( kt )      ! Empty routine 
    194       WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt 
     306      INTEGER, intent(in) :: kt 
     307      WRITE(*,*) 'bdy_dyn2d: You should not have seen this print! error?', kt 
    195308   END SUBROUTINE bdy_dyn2d 
     309 
    196310#endif 
    197311 
    198312   !!====================================================================== 
    199313END MODULE bdydyn2d 
     314 
Note: See TracChangeset for help on using the changeset viewer.