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 3116 for branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90 – NEMO

Ignore:
Timestamp:
2011-11-15T21:55:40+01:00 (13 years ago)
Author:
cetlod
Message:

dev_NEMO_MERGE_2011: add in changes dev_NOC_UKMO_MERGE developments

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90

    r2528 r3116  
    1515   !!   'key_bdy' :                    Unstructured Open Boundary Condition 
    1616   !!---------------------------------------------------------------------- 
    17    !!   bdy_dyn_frs    : relaxation of velocities on unstructured open boundary 
    18    !!   bdy_dyn_fla    : Flather condition for barotropic solution 
     17   !!   bdy_dyn3d        : apply open boundary conditions to baroclinic velocities 
     18   !!   bdy_dyn3d_frs    : apply Flow Relaxation Scheme 
    1919   !!---------------------------------------------------------------------- 
    2020   USE oce             ! ocean dynamics and tracers  
    2121   USE dom_oce         ! ocean space and time domain 
     22   USE dynspg_oce       
    2223   USE bdy_oce         ! ocean open boundary conditions 
    23    USE dynspg_oce      ! for barotropic variables 
    24    USE phycst          ! physical constants 
     24   USE bdydyn2d        ! open boundary conditions for barotropic solution 
     25   USE bdydyn3d        ! open boundary conditions for baroclinic velocities 
    2526   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    26    USE bdytides        ! for tidal harmonic forcing at boundary 
    2727   USE in_out_manager  ! 
    2828 
     
    3030   PRIVATE 
    3131 
    32    PUBLIC   bdy_dyn_frs   ! routine called in dynspg_flt (free surface case ONLY) 
    33 # if defined key_dynspg_exp || defined key_dynspg_ts 
    34    PUBLIC   bdy_dyn_fla   ! routine called in dynspg_flt (free surface case ONLY) 
    35 # endif 
     32   PUBLIC   bdy_dyn     ! routine called in dynspg_flt (if lk_dynspg_flt) or  
     33                        ! dyn_nxt (if lk_dynspg_ts or lk_dynspg_exp) 
    3634 
     35#  include "domzgr_substitute.h90" 
    3736   !!---------------------------------------------------------------------- 
    3837   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
     
    4241CONTAINS 
    4342 
    44    SUBROUTINE bdy_dyn_frs( kt ) 
     43   SUBROUTINE bdy_dyn( kt, dyn3d_only ) 
    4544      !!---------------------------------------------------------------------- 
    46       !!                  ***  SUBROUTINE bdy_dyn_frs  *** 
     45      !!                  ***  SUBROUTINE bdy_dyn  *** 
    4746      !! 
    48       !! ** Purpose : - Apply the Flow Relaxation Scheme for dynamic in the   
    49       !!                case of unstructured open boundaries. 
     47      !! ** Purpose : - Wrapper routine for bdy_dyn2d and bdy_dyn3d. 
    5048      !! 
    51       !! References :- Engedahl H., 1995: Use of the flow relaxation scheme in  
    52       !!               a three-dimensional baroclinic ocean model with realistic 
    53       !!               topography. Tellus, 365-382. 
    5449      !!---------------------------------------------------------------------- 
    55       INTEGER, INTENT( in ) ::   kt   ! Main time step counter 
     50      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released 
     51      USE wrk_nemo, ONLY: wrk_2d_7, wrk_2d_8      ! 2D workspace 
    5652      !! 
    57       INTEGER  ::   jb, jk         ! dummy loop indices 
    58       INTEGER  ::   ii, ij, igrd   ! local integers 
    59       REAL(wp) ::   zwgt           ! boundary weight 
    60       !!---------------------------------------------------------------------- 
    61       ! 
    62       IF(ln_dyn_frs) THEN       ! If this is false, then this routine does nothing.  
    63          ! 
    64          IF( kt == nit000 ) THEN 
    65             IF(lwp) WRITE(numout,*) 
    66             IF(lwp) WRITE(numout,*) 'bdy_dyn_frs : Flow Relaxation Scheme on momentum' 
    67             IF(lwp) WRITE(numout,*) '~~~~~~~' 
    68          ENDIF 
    69          ! 
    70          igrd = 2                      ! Relaxation of zonal velocity 
    71          DO jb = 1, nblen(igrd) 
    72             DO jk = 1, jpkm1 
    73                ii   = nbi(jb,igrd) 
    74                ij   = nbj(jb,igrd) 
    75                zwgt = nbw(jb,igrd) 
    76                ua(ii,ij,jk) = ( ua(ii,ij,jk) * ( 1.- zwgt ) + ubdy(jb,jk) * zwgt ) * umask(ii,ij,jk) 
    77             END DO 
    78          END DO 
    79          ! 
    80          igrd = 3                      ! Relaxation of meridional velocity 
    81          DO jb = 1, nblen(igrd) 
    82             DO jk = 1, jpkm1 
    83                ii   = nbi(jb,igrd) 
    84                ij   = nbj(jb,igrd) 
    85                zwgt = nbw(jb,igrd) 
    86                va(ii,ij,jk) = ( va(ii,ij,jk) * ( 1.- zwgt ) + vbdy(jb,jk) * zwgt ) * vmask(ii,ij,jk) 
    87             END DO 
    88          END DO  
    89          CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated 
    90          ! 
    91       ENDIF ! ln_dyn_frs 
    92       ! 
    93    END SUBROUTINE bdy_dyn_frs 
     53      INTEGER, INTENT( in )           :: kt               ! Main time step counter 
     54      LOGICAL, INTENT( in ), OPTIONAL :: dyn3d_only       ! T => only update baroclinic velocities 
     55      !! 
     56      INTEGER               :: jk,ii,ij,ib,igrd     ! Loop counter 
     57      LOGICAL               :: ll_dyn2d, ll_dyn3d   
     58      !! 
    9459 
     60      IF(wrk_in_use(2, 7,8) ) THEN 
     61         CALL ctl_stop('bdy_dyn: ERROR: requested workspace arrays are unavailable.')   ;   RETURN 
     62      END IF 
    9563 
    96 # if defined   key_dynspg_exp   ||   defined key_dynspg_ts 
    97    !!---------------------------------------------------------------------- 
    98    !!   'key_dynspg_exp'        OR              explicit sea surface height 
    99    !!   'key_dynspg_ts '                  split-explicit sea surface height 
    100    !!---------------------------------------------------------------------- 
    101     
    102 !! Option to use Flather with dynspg_flt not coded yet... 
     64      ll_dyn2d = .true. 
     65      ll_dyn3d = .true. 
    10366 
    104    SUBROUTINE bdy_dyn_fla( pssh ) 
    105       !!---------------------------------------------------------------------- 
    106       !!                 ***  SUBROUTINE bdy_dyn_fla  *** 
    107       !!              
    108       !!              - Apply Flather boundary conditions on normal barotropic velocities  
    109       !!                (ln_dyn_fla=.true. or ln_tides=.true.) 
    110       !! 
    111       !! ** WARNINGS about FLATHER implementation: 
    112       !!1. According to Palma and Matano, 1998 "after ssh" is used.  
    113       !!   In ROMS and POM implementations, it is "now ssh". In the current  
    114       !!   implementation (tested only in the EEL-R5 conf.), both cases were unstable.  
    115       !!   So I use "before ssh" in the following. 
    116       !! 
    117       !!2. We assume that the normal ssh gradient at the bdy is zero. As a matter of  
    118       !!   fact, the model ssh just inside the dynamical boundary is used (the outside   
    119       !!   ssh in the code is not updated). 
    120       !! 
    121       !! References:  Flather, R. A., 1976: A tidal model of the northwest European 
    122       !!              continental shelf. Mem. Soc. R. Sci. Liege, Ser. 6,10, 141-164.      
    123       !!---------------------------------------------------------------------- 
    124       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pssh 
     67      IF( PRESENT(dyn3d_only) ) THEN 
     68         IF( dyn3d_only ) ll_dyn2d = .false. 
     69      ENDIF 
    12570 
    126       INTEGER  ::   jb, igrd                         ! dummy loop indices 
    127       INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses 
    128       REAL(wp) ::   zcorr                            ! Flather correction 
    129       REAL(wp) ::   zforc                            ! temporary scalar 
    130       !!---------------------------------------------------------------------- 
     71      !------------------------------------------------------- 
     72      ! Set pointers 
     73      !------------------------------------------------------- 
    13174 
    132       ! ---------------------------------! 
    133       ! Flather boundary conditions     :! 
    134       ! ---------------------------------!  
    135       
    136       IF(ln_dyn_fla .OR. ln_tides) THEN ! If these are both false, then this routine does nothing.  
     75      pssh => sshn 
     76      phur => hur 
     77      phvr => hvr 
     78      pu2d => wrk_2d_7 
     79      pv2d => wrk_2d_8 
    13780 
    138          ! Fill temporary array with ssh data (here spgu): 
    139          igrd = 4 
    140          spgu(:,:) = 0.0 
    141          DO jb = 1, nblenrim(igrd) 
    142             ii = nbi(jb,igrd) 
    143             ij = nbj(jb,igrd) 
    144             IF( ln_dyn_fla ) spgu(ii, ij) = sshbdy(jb) 
    145             IF( ln_tides )   spgu(ii, ij) = spgu(ii, ij) + sshtide(jb) 
    146          END DO 
    147          ! 
    148          igrd = 5      ! Flather bc on u-velocity;  
    149          !             ! remember that flagu=-1 if normal velocity direction is outward 
    150          !             ! I think we should rather use after ssh ? 
    151          DO jb = 1, nblenrim(igrd) 
    152             ii  = nbi(jb,igrd) 
    153             ij  = nbj(jb,igrd)  
    154             iim1 = ii + MAX( 0, INT( flagu(jb) ) )   ! T pts i-indice inside the boundary 
    155             iip1 = ii - MIN( 0, INT( flagu(jb) ) )   ! T pts i-indice outside the boundary  
    156             ! 
    157             zcorr = - flagu(jb) * SQRT( grav * hur_e(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) ) 
    158             zforc = ubtbdy(jb) + utide(jb) 
    159             ua_e(ii,ij) = zforc + zcorr * umask(ii,ij,1)  
    160          END DO 
    161          ! 
    162          igrd = 6      ! Flather bc on v-velocity 
    163          !             ! remember that flagv=-1 if normal velocity direction is outward 
    164          DO jb = 1, nblenrim(igrd) 
    165             ii  = nbi(jb,igrd) 
    166             ij  = nbj(jb,igrd)  
    167             ijm1 = ij + MAX( 0, INT( flagv(jb) ) )   ! T pts j-indice inside the boundary 
    168             ijp1 = ij - MIN( 0, INT( flagv(jb) ) )   ! T pts j-indice outside the boundary  
    169             ! 
    170             zcorr = - flagv(jb) * SQRT( grav * hvr_e(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) ) 
    171             zforc = vbtbdy(jb) + vtide(jb) 
    172             va_e(ii,ij) = zforc + zcorr * vmask(ii,ij,1) 
    173          END DO 
    174          CALL lbc_lnk( ua_e, 'U', -1. )   ! Boundary points should be updated 
    175          CALL lbc_lnk( va_e, 'V', -1. )   ! 
    176          ! 
    177       ENDIF ! ln_dyn_fla .or. ln_tides 
    178       ! 
    179    END SUBROUTINE bdy_dyn_fla 
    180 #endif 
     81      !------------------------------------------------------- 
     82      ! Split velocities into barotropic and baroclinic parts 
     83      !------------------------------------------------------- 
     84 
     85      pu2d(:,:) = 0.e0 
     86      pv2d(:,:) = 0.e0 
     87      DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
     88          pu2d(:,:) = pu2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
     89          pv2d(:,:) = pv2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
     90      END DO 
     91      pu2d(:,:) = pu2d(:,:) * phur(:,:) 
     92      pv2d(:,:) = pv2d(:,:) * phvr(:,:) 
     93      DO jk = 1 , jpkm1 
     94         ua(:,:,jk) = ua(:,:,jk) - pu2d(:,:) 
     95         va(:,:,jk) = va(:,:,jk) - pv2d(:,:) 
     96      END DO 
     97 
     98      !------------------------------------------------------- 
     99      ! Apply boundary conditions to barotropic and baroclinic 
     100      ! parts separately 
     101      !------------------------------------------------------- 
     102 
     103      IF( ll_dyn2d ) CALL bdy_dyn2d( kt ) 
     104 
     105      IF( ll_dyn3d ) CALL bdy_dyn3d( kt ) 
     106 
     107      !------------------------------------------------------- 
     108      ! Recombine velocities 
     109      !------------------------------------------------------- 
     110 
     111      DO jk = 1 , jpkm1 
     112         ua(:,:,jk) = ( ua(:,:,jk) + pu2d(:,:) ) * umask(:,:,jk) 
     113         va(:,:,jk) = ( va(:,:,jk) + pv2d(:,:) ) * vmask(:,:,jk) 
     114      END DO 
     115 
     116      IF(wrk_not_released(2, 7,8) )    CALL ctl_stop('bdy_dyn: ERROR: failed to release workspace arrays.') 
     117 
     118   END SUBROUTINE bdy_dyn 
    181119 
    182120#else 
     
    185123   !!---------------------------------------------------------------------- 
    186124CONTAINS 
    187    SUBROUTINE bdy_dyn_frs( kt )      ! Empty routine 
    188       WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt 
    189    END SUBROUTINE bdy_dyn_frs 
    190    SUBROUTINE bdy_dyn_fla( pssh )    ! Empty routine 
    191       REAL :: pssh(:,:) 
    192       WRITE(*,*) 'bdy_dyn_fla: You should not have seen this print! error?', pssh(1,1) 
    193    END SUBROUTINE bdy_dyn_fla 
     125   SUBROUTINE bdy_dyn( kt )      ! Empty routine 
     126      WRITE(*,*) 'bdy_dyn: You should not have seen this print! error?', kt 
     127   END SUBROUTINE bdy_dyn 
    194128#endif 
    195129 
Note: See TracChangeset for help on using the changeset viewer.