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 4007 for branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90 – NEMO

Ignore:
Timestamp:
2013-08-28T10:10:35+02:00 (11 years ago)
Author:
davestorkey
Message:
  1. Bug fixes for flagu/flagv calculation in bdyini.F90.
  2. Introduce masking of derivatives in radiation velocity calculation in bdylib.F90
  3. Change relaxation term in Orlanski implementation to explicit timestepping in bdylib.F90.
  4. Remove bdyfmask (redundant).
File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2013/dev_r3987_METO1_MERCATOR6_OBC_BDY_merge/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90

    r3994 r4007  
    3434CONTAINS 
    3535 
    36    SUBROUTINE bdy_orlanski_2d( idx, igrd, phib, phia, phi_ext, mask, ll_npo ) 
     36   SUBROUTINE bdy_orlanski_2d( idx, igrd, phib, phia, phi_ext, ll_npo ) 
    3737      !!---------------------------------------------------------------------- 
    3838      !!                 ***  SUBROUTINE bdy_orlanski_2d  *** 
    … …  
    5050      REAL(wp), DIMENSION(:,:),   INTENT(inout)  ::   phia     ! model after 2D field (to be updated) 
    5151      REAL(wp), DIMENSION(:),     INTENT(in)     ::   phi_ext  ! external forcing data 
    52       REAL(wp), DIMENSION(:,:),   INTENT(in)     ::   mask     ! land/sea mask 
    5352      LOGICAL,                    INTENT(in)     ::   ll_npo   ! switch for NPO version 
    5453 
    … …  
    5756      INTEGER  ::   iijm1, iijp1, ijjm1, ijjp1             ! 2D addresses 
    5857      INTEGER  ::   iibm1jp1, iibm1jm1, ijbm1jp1, ijbm1jm1 ! 2D addresses 
     58      INTEGER  ::   ii_offset, ij_offset                   ! offsets for mask indices 
    5959      INTEGER  ::   flagu, flagv                           ! short cuts 
    6060      REAL(wp) ::   zdt, zdx, zdy, znor2, zcx, zcy         ! intermediate calculations 
    61       REAL(wp) ::   zout, zwgt, zdy_centred, zsign_ups 
     61      REAL(wp) ::   zout, zwgt, zdy_centred 
     62      REAL(wp) ::   zdy_left, zdy_right, zsign_ups 
     63      REAL(wp), POINTER, DIMENSION(:,:)          :: pmask      ! land/sea mask for field 
     64      REAL(wp), POINTER, DIMENSION(:,:)          :: pmask_xdiv ! land/sea mask for x-derivatives 
     65      REAL(wp), POINTER, DIMENSION(:,:)          :: pmask_ydiv ! land/sea mask for y-derivatives 
    6266      !!---------------------------------------------------------------------- 
    6367 
    … …  
    6872      ! ----------------------------------!  
    6973      
     74      SELECT CASE(igrd) 
     75         CASE(1) 
     76            pmask => tmask(:,:,1) 
     77            pmask_xdiv => umask(:,:,1) 
     78            ii_offset = 0 
     79            pmask_ydiv => vmask(:,:,1) 
     80            ij_offset = 0 
     81         CASE(2) 
     82            pmask => umask(:,:,1) 
     83            pmask_xdiv => tmask(:,:,1) 
     84            ii_offset = 1 
     85            pmask_ydiv => fmask(:,:,1) 
     86            ij_offset = 0 
     87         CASE(3) 
     88            pmask => vmask(:,:,1) 
     89            pmask_xdiv => fmask(:,:,1) 
     90            ii_offset = 0 
     91            pmask_ydiv => tmask(:,:,1) 
     92            ij_offset = 1 
     93         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_orlanksi_2d' ) 
     94      END SELECT 
    7095      ! 
    7196      DO jb = 1, idx%nblenrim(igrd) 
    … …  
    86111         ijbm1jm1 = ij + flagv - abs(flagu) ; ijbm1jp1 = ij + flagv + abs(flagu)  
    87112         ! 
    88          ! calculate normal (zcx) and tangential (zcy) components of radiation velocities: 
     113         ! Calculate normal (zcx) and tangential (zcy) components of radiation velocities. 
     114         ! Mask derivatives to ensure correct land boundary conditions for each variable. 
     115         ! Centred derivative is calculated as average of "left" and "right" derivatives for  
     116         ! this reason.  
    89117         zdt = phia(iibm1,ijbm1) - phib(iibm1,ijbm1) 
    90          zdx = phia(iibm1,ijbm1) - phia(iibm2,ijbm2) 
    91          zdy_centred = phib(iibm1jp1,ijbm1jp1) - phib(iibm1jm1,ijbm1jm1) 
     118         zdx = ( phia(iibm1,ijbm1) - phia(iibm2,ijbm2) )                          & 
     119        &    * ( abs(iibm1-iibm2) * pmask_xdiv(iibm2+ii_offset,ijbm2          )   & 
     120        &      + abs(ijbm1-ijbm2) * pmask_ydiv(iibm2          ,ijbm2+ij_offset) )  
     121         zdy_left  = phib(iibm1   ,ijbm1   ) - phib(iibm1jm1,ijbm1jm1)                  & 
     122        &    * ( (iibm1-iibm1jm1) * pmask_xdiv(iibm1jm1+ii_offset,ijbm1jm1          )   &  
     123        &      + (ijbm1-ijbm1jm1) * pmask_ydiv(iibm1jm1          ,ijbm1jm1+ij_offset) )  
     124         zdy_right = phib(iibm1jp1,ijbm1jp1) - phib(iibm1   ,ijbm1)                     & 
     125        &    * ( (iibm1jp1-iibm1) * pmask_xdiv(iibm1+ii_offset,ijbm1)                   & 
     126        &      + (ijbm1jp1-ijbm1) * pmask_ydiv(iibm1          ,ijbm1+ij_offset) )  
     127         zdy_centred = 0.5 * ( zdy_left + zdy_right ) 
     128!!$         zdy_centred = phib(iibm1jp1,ijbm1jp1) - phib(iibm1jm1,ijbm1jm1) 
    92129         ! upstream differencing for tangential derivatives 
    93130         zsign_ups = sign( 1., zdt * zdy_centred ) 
    94131         zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) ) 
    95          zdy = zsign_ups        * ( phib(iibm1   ,ijbm1   ) - phib(iibm1jm1,ijbm1jm1) ) & 
    96        &        + (1. - zsign_ups) * ( phib(iibm1jp1,ijbm1jp1) - phib(iibm1   ,ijbm1   ) ) 
     132         zdy = zsign_ups * zdy_left + (1. - zsign_ups) * zdy_right 
    97133         znor2 = zdx * zdx + zdy * zdy 
    98134         znor2 = max(znor2,rsmall) 
    … …  
    106142         ! only apply radiation on outflow points  
    107143         if( ll_npo ) then     !! NPO version !! 
    108             phia(ii,ij) = (1.-zout) * phia(ii,ij)                                                   & 
     144            phia(ii,ij) = (1.-zout) * phib(ii,ij)                                                   & 
    109145           &            + zout      * ( phib(ii,ij) + zcx*phia(iibm1,ijbm1) ) / ( 1. + zcx )  
    110146         else                  !! full oblique radiation !! 
    111147            zsign_ups = sign( 1., zcy ) 
    112148            zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) ) 
    113             phia(ii,ij) = (1.-zout) * phia(ii,ij)                                                   & 
     149            phia(ii,ij) = (1.-zout) * phib(ii,ij)                                                   & 
    114150           &            + zout      * ( phib(ii,ij) + zcx*phia(iibm1,ijbm1)                         & 
    115151           &                    - zsign_ups      * zcy * ( phib(ii   ,ij   ) - phib(iijm1,ijjm1 ) ) & 
    116152           &                    - (1.-zsign_ups) * zcy * ( phib(iijp1,ijjp1) - phib(ii   ,ij    ) ) ) / ( 1. + zcx )  
    117153         end if 
    118          phia(ii,ij) = phia(ii,ij) + zwgt * ( phi_ext(jb) - phia(ii,ij) )  
    119          phia(ii,ij) = phia(ii,ij) * mask(ii,ij) 
     154!!$         phia(ii,ij) = phia(ii,ij) + zwgt * ( phi_ext(jb) - phib(ii,ij) )  
     155         phia(ii,ij) = phia(ii,ij) * pmask(ii,ij) 
    120156      END DO 
    121157      ! 
    … …  
    125161 
    126162 
    127    SUBROUTINE bdy_orlanski_3d( idx, igrd, phib, phia, phi_ext, mask, ll_npo ) 
     163   SUBROUTINE bdy_orlanski_3d( idx, igrd, phib, phia, phi_ext, ll_npo ) 
    128164      !!---------------------------------------------------------------------- 
    129165      !!                 ***  SUBROUTINE bdy_orlanski_3d  *** 
    … …  
    141177      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   phia     ! model after 3D field (to be updated) 
    142178      REAL(wp), DIMENSION(:,:),   INTENT(in)     ::   phi_ext  ! external forcing data 
    143       REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   mask     ! land/sea mask 
    144179      LOGICAL,                    INTENT(in)     ::   ll_npo   ! switch for NPO version 
    145180 
    … …  
    148183      INTEGER  ::   iijm1, iijp1, ijjm1, ijjp1             ! 2D addresses 
    149184      INTEGER  ::   iibm1jp1, iibm1jm1, ijbm1jp1, ijbm1jm1 ! 2D addresses 
     185      INTEGER  ::   ii_offset, ij_offset                   ! offsets for mask indices 
    150186      INTEGER  ::   flagu, flagv                           ! short cuts 
    151187      REAL(wp) ::   zdt, zdx, zdy, znor2, zcx, zcy         ! intermediate calculations 
    152       REAL(wp) ::   zout, zwgt, zdy_centred, zsign_ups 
     188      REAL(wp) ::   zout, zwgt, zdy_centred 
     189      REAL(wp) ::   zdy_left, zdy_right,  zsign_ups 
     190      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask      ! land/sea mask for field 
     191      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask_xdiv ! land/sea mask for x-derivatives 
     192      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask_ydiv ! land/sea mask for y-derivatives 
    153193      !!---------------------------------------------------------------------- 
    154194 
    … …  
    159199      ! ----------------------------------!  
    160200      
     201      SELECT CASE(igrd) 
     202         CASE(1) 
     203            pmask => tmask(:,:,:) 
     204            pmask_xdiv => umask(:,:,:) 
     205            ii_offset = 0 
     206            pmask_ydiv => vmask(:,:,:) 
     207            ij_offset = 0 
     208         CASE(2) 
     209            pmask => umask(:,:,:) 
     210            pmask_xdiv => tmask(:,:,:) 
     211            ii_offset = 1 
     212            pmask_ydiv => fmask(:,:,:) 
     213            ij_offset = 0 
     214         CASE(3) 
     215            pmask => vmask(:,:,:) 
     216            pmask_xdiv => fmask(:,:,:) 
     217            ii_offset = 0 
     218            pmask_ydiv => tmask(:,:,:) 
     219            ij_offset = 1 
     220         CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_orlanksi_2d' ) 
     221      END SELECT 
     222 
    161223      DO jk = 1, jpk 
    162224         !             
    … …  
    178240            ijbm1jm1 = ij + flagv - abs(flagu) ; ijbm1jp1 = ij + flagv + abs(flagu)  
    179241            ! 
    180             ! calculate normal (zcx) and tangential (zcy) components of radiation velocities: 
     242            ! Calculate normal (zcx) and tangential (zcy) components of radiation velocities. 
     243            ! Mask derivatives to ensure correct land boundary conditions for each variable. 
     244            ! Centred derivative is calculated as average of "left" and "right" derivatives for  
     245            ! this reason.  
    181246            zdt = phia(iibm1,ijbm1,jk) - phib(iibm1,ijbm1,jk) 
    182             zdx = phia(iibm1,ijbm1,jk) - phia(iibm2,ijbm2,jk) 
    183             zdy_centred = phib(iibm1jp1,ijbm1jp1,jk) - phib(iibm1jm1,ijbm1jm1,jk) 
     247            zdx = phia(iibm1,ijbm1,jk) - phia(iibm2,ijbm2,jk)                     & 
     248        &    * ( (iibm1-iibm2) * pmask_xdiv(iibm2+ii_offset,ijbm2          ,jk)   & 
     249        &      + (ijbm1-ijbm2) * pmask_ydiv(iibm2          ,ijbm2+ij_offset,jk) )  
     250            zdy_left  = phib(iibm1   ,ijbm1   ,jk) - phib(iibm1jm1,ijbm1jm1,jk)            & 
     251        &    * ( (iibm1-iibm1jm1) * pmask_xdiv(iibm1jm1+ii_offset,ijbm1jm1          ,jk)   &  
     252        &      + (ijbm1-ijbm1jm1) * pmask_ydiv(iibm1jm1          ,ijbm1jm1+ij_offset,jk) )  
     253            zdy_right = phib(iibm1jp1,ijbm1jp1,jk) - phib(iibm1   ,ijbm1   ,jk)      & 
     254        &    * ( (iibm1jp1-iibm1) * pmask_xdiv(iibm1+ii_offset,ijbm1          ,jk)   & 
     255        &      + (ijbm1jp1-ijbm1) * pmask_ydiv(iibm1          ,ijbm1+ij_offset,jk) )  
     256            zdy_centred = 0.5 * ( zdy_left + zdy_right ) 
     257!!$            zdy_centred = phib(iibm1jp1,ijbm1jp1,jk) - phib(iibm1jm1,ijbm1jm1,jk) 
    184258            ! upstream differencing for tangential derivatives 
    185259            zsign_ups = sign( 1., zdt * zdy_centred ) 
    186260            zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) ) 
    187             zdy = zsign_ups        * ( phib(iibm1   ,ijbm1   ,jk) - phib(iibm1jm1,ijbm1jm1,jk) ) & 
    188            &    + (1. - zsign_ups) * ( phib(iibm1jp1,ijbm1jp1,jk) - phib(iibm1   ,ijbm1   ,jk) ) 
     261            zdy = zsign_ups * zdy_left + (1. - zsign_ups) * zdy_right 
    189262            znor2 = zdx * zdx + zdy * zdy 
    190263            znor2 = max(znor2,rsmall) 
    … …  
    198271            ! only apply radiation on outflow points  
    199272            if( ll_npo ) then     !! NPO version !! 
    200                phia(ii,ij,jk) = (1.-zout) * phia(ii,ij,jk)                                                   & 
     273               phia(ii,ij,jk) = (1.-zout) * phib(ii,ij,jk)                                                   & 
    201274              &               + zout      * ( phib(ii,ij,jk) + zcx*phia(iibm1,ijbm1,jk) ) / ( 1. + zcx )  
    202275            else                  !! full oblique radiation !! 
    203276               zsign_ups = sign( 1., zcy ) 
    204277               zsign_ups = 0.5*( zsign_ups + abs(zsign_ups) ) 
    205                phia(ii,ij,jk) = (1.-zout) * phia(ii,ij,jk)                                                   & 
     278               phia(ii,ij,jk) = (1.-zout) * phib(ii,ij,jk)                                                   & 
    206279              &               + zout      * ( phib(ii,ij,jk) + zcx*phia(iibm1,ijbm1,jk)                      & 
    207280              &                    - zsign_ups      * zcy * ( phib(ii   ,ij   ,jk) - phib(iijm1,ijjm1,jk ) ) & 
    208281              &                    - (1.-zsign_ups) * zcy * ( phib(iijp1,ijjp1,jk) - phib(ii   ,ij   ,jk ) ) ) / ( 1. + zcx )  
    209282            end if 
    210             phia(ii,ij,jk) = phia(ii,ij,jk) + zwgt * ( phi_ext(jb,jk) - phia(ii,ij,jk) )  
    211             phia(ii,ij,jk) = phia(ii,ij,jk) * mask(ii,ij,jk) 
     283!!$            phia(ii,ij,jk) = phia(ii,ij,jk) + zwgt * ( phi_ext(jb,jk) - phib(ii,ij,jk) )  
     284            phia(ii,ij,jk) = phia(ii,ij,jk) * pmask(ii,ij,jk) 
    212285         END DO 
    213286         ! 
Note: See TracChangeset for help on using the changeset viewer.