Changeset 11074


Ignore:
Timestamp:
2019-06-04T17:46:46+02:00 (17 months ago)
Author:
girrmann
Message:

dev_r10984_HPC-13 : bug fix and cleaning of flather, see #2287 and #2285

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydyn2d.F90

    r11071 r11074  
    1414   !!   bdy_ssh            : Duplicate sea level across open boundaries 
    1515   !!---------------------------------------------------------------------- 
    16    USE oce             ! ocean dynamics and tracers  
     16   USE oce, only : sshdta => spgu             ! ocean dynamics and tracers  
    1717   USE dom_oce         ! ocean space and time domain 
    1818   USE bdy_oce         ! ocean open boundary conditions 
     
    169169 
    170170      INTEGER  ::   jb, igrd                         ! dummy loop indices 
    171       INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses 
     171      INTEGER  ::   ii, ij                           ! 2D addresses 
     172      INTEGER  ::   iiTrim, ijTrim                   ! T pts i/j-indice on the rim 
     173      INTEGER  ::   iiToce, ijToce, iiUoce, ijVoce   ! T, U and V pts i/j-indice of the ocean next to the rim 
    172174      REAL(wp), POINTER :: flagu, flagv              ! short cuts 
    173       REAL(wp) ::   zcorr                            ! Flather correction 
    174       REAL(wp) ::   zforc                            ! temporary scalar 
    175       REAL(wp) ::   zflag, z1_2                      !    "        " 
     175      REAL(wp) ::   zfla                             ! Flather correction 
     176      REAL(wp) ::   z1_2                             !  
    176177      !!---------------------------------------------------------------------- 
    177178 
     
    180181      ! ---------------------------------! 
    181182      ! Flather boundary conditions     :! 
    182       ! ---------------------------------!  
    183       
    184 !!! REPLACE spgu with nemo_wrk work space 
    185  
    186       ! Fill temporary array with ssh data (here spgu): 
     183      ! ---------------------------------! 
     184 
     185      ! Fill temporary array with ssh data (here we use spgu with the alias sshdta): 
    187186      igrd = 1 
    188       spgu(:,:) = 0.0 
     187      sshdta(:,:) = 0.0 
    189188      DO jb = 1, idx%nblenrim(igrd) 
    190189         ii = idx%nbi(jb,igrd) 
    191190         ij = idx%nbj(jb,igrd) 
    192191         IF( ll_wd ) THEN 
    193             spgu(ii, ij) = dta%ssh(jb)  - ssh_ref  
     192            sshdta(ii, ij) = dta%ssh(jb)  - ssh_ref  
    194193         ELSE 
    195             spgu(ii, ij) = dta%ssh(jb) 
     194            sshdta(ii, ij) = dta%ssh(jb) 
    196195         ENDIF 
    197196      END DO 
     
    201200      !             ! I think we should rather use after ssh ? 
    202201      DO jb = 1, idx%nblenrim(igrd) 
    203          ii  = idx%nbi(jb,igrd) 
    204          ij  = idx%nbj(jb,igrd) 
     202         ii = idx%nbi(jb,igrd) 
     203         ij = idx%nbj(jb,igrd) 
    205204         flagu => idx%flagu(jb,igrd) 
    206          iim1 = ii + MAX( 0, INT( flagu ) )   ! T pts i-indice inside the boundary 
    207          iip1 = ii - MIN( 0, INT( flagu ) )   ! T pts i-indice outside the boundary  
    208          IF( iim1 > jpi .OR. iip1 > jpi )   CYCLE 
    209          ! 
    210          zcorr = - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) ) 
    211  
    212          ! jchanut tschanges: Set zflag to 0 below to revert to Flather scheme 
    213          ! Use characteristics method instead 
    214          zflag = ABS(flagu) 
    215          zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pua2d(iim1,ij) 
    216          pua2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1)  
     205         IF( flagu == 0. ) THEN 
     206            pua2d(ii,ij) = dta%u2d(jb) 
     207         ELSE      ! T pts j-indice       on the rim          on the ocean next to the rim on T and U points 
     208            IF( flagu ==  1. ) THEN   ;   iiTrim = ii     ;   iiToce = ii+1   ;   iiUoce = ii+1   ;   ENDIF 
     209            IF( flagu == -1. ) THEN   ;   iiTrim = ii+1   ;   iiToce = ii     ;   iiUoce = ii-1   ;   ENDIF 
     210            ! 
     211            ! Rare case : rim is parallel to the mpi subdomain border and on the halo : point will be received 
     212            IF( iiTrim > jpi .OR. iiToce > jpi .OR. iiUoce > jpi .OR. iiUoce < 1 )   CYCLE    
     213            ! 
     214            zfla = pua2d(iiUoce,ij) - flagu * SQRT( grav * phur(ii, ij) ) * ( pssh(iiToce,ij) - sshdta(iiTrim,ij) ) 
     215            ! 
     216            ! jchanut tschanges, use characteristics method (Blayo et Debreu, 2005) : 
     217            ! mix Flather scheme with velocity of the ocean next to the rim 
     218            pua2d(ii,ij) =  z1_2 * ( dta%u2d(jb) + zfla ) 
     219         END IF 
    217220      END DO 
    218221      ! 
     
    220223      !             ! remember that flagv=-1 if normal velocity direction is outward 
    221224      DO jb = 1, idx%nblenrim(igrd) 
    222          ii  = idx%nbi(jb,igrd) 
    223          ij  = idx%nbj(jb,igrd) 
     225         ii = idx%nbi(jb,igrd) 
     226         ij = idx%nbj(jb,igrd) 
    224227         flagv => idx%flagv(jb,igrd) 
    225          ijm1 = ij + MAX( 0, INT( flagv ) )   ! T pts j-indice inside the boundary 
    226          ijp1 = ij - MIN( 0, INT( flagv ) )   ! T pts j-indice outside the boundary  
    227          IF( ijm1 > jpj .OR. ijp1 > jpj )   CYCLE 
    228          ! 
    229          zcorr = - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) ) 
    230  
    231          ! jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme 
    232          ! Use characteristics method instead 
    233          zflag = ABS(flagv) 
    234          zforc = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pva2d(ii,ijm1) 
    235          pva2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1) 
     228         IF( flagv == 0. ) THEN 
     229            pva2d(ii,ij) = dta%v2d(jb) 
     230         ELSE      ! T pts j-indice       on the rim          on the ocean next to the rim on T and V points 
     231            IF( flagv ==  1. ) THEN   ;   ijTrim = ij     ;   ijToce = ij+1   ;   ijVoce = ij+1   ;   ENDIF 
     232            IF( flagv == -1. ) THEN   ;   ijTrim = ij+1   ;   ijToce = ij     ;   ijVoce = ij-1   ;   ENDIF 
     233            ! 
     234            ! Rare case : rim is parallel to the mpi subdomain border and on the halo : point will be received 
     235            IF( ijTrim > jpj .OR. ijToce > jpj .OR. ijVoce > jpj .OR. ijVoce < 1 )   CYCLE 
     236            ! 
     237            zfla = pva2d(ii,ijVoce) - flagv * SQRT( grav * phvr(ii, ij) ) * ( pssh(ii,ijToce) - sshdta(ii,ijTrim) ) 
     238            ! 
     239            ! jchanut tschanges, use characteristics method (Blayo et Debreu, 2005) : 
     240            ! mix Flather scheme with velocity of the ocean next to the rim 
     241            pva2d(ii,ij) =  z1_2 * ( dta%v2d(jb) + zfla ) 
     242         END IF 
    236243      END DO 
    237244      ! 
Note: See TracChangeset for help on using the changeset viewer.