Changeset 10455


Ignore:
Timestamp:
2019-01-04T19:34:14+01:00 (2 years ago)
Author:
mathiot
Message:

suggestion for ticket #2200 (bdyvol)

Location:
NEMO/branches/UKMO/dev_r10448_bdyvol/src/OCE
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/dev_r10448_bdyvol/src/OCE/BDY/bdyini.F90

    r10425 r10455  
    133133      INTEGER  ::   i_offset, j_offset                     !   -       - 
    134134      INTEGER , POINTER  ::  nbi, nbj, nbr                 ! short cuts 
    135       REAL(wp), POINTER  ::  flagu, flagv                  !    -   - 
    136135      REAL(wp), POINTER, DIMENSION(:,:)       ::   pmask    ! pointer to 2D mask fields 
    137136      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars 
     
    12441243         ! 
    12451244      END DO 
    1246  
    1247       ! Compute total lateral surface for volume correction: 
    1248       ! ---------------------------------------------------- 
    1249       ! JC: this must be done at each time step with non-linear free surface 
    1250       bdysurftot = 0._wp  
    1251       IF( ln_vol ) THEN   
    1252          igrd = 2      ! Lateral surface at U-points 
    1253          DO ib_bdy = 1, nb_bdy 
    1254             DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    1255                nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
    1256                nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
    1257                flagu => idx_bdy(ib_bdy)%flagu(ib,igrd) 
    1258                bdysurftot = bdysurftot + hu_n   (nbi  , nbj)                           & 
    1259                   &                    * e2u    (nbi  , nbj) * ABS( flagu ) & 
    1260                   &                    * tmask_i(nbi  , nbj)                           & 
    1261                   &                    * tmask_i(nbi+1, nbj)                    
    1262             END DO 
    1263          END DO 
    1264  
    1265          igrd=3 ! Add lateral surface at V-points 
    1266          DO ib_bdy = 1, nb_bdy 
    1267             DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
    1268                nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
    1269                nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
    1270                flagv => idx_bdy(ib_bdy)%flagv(ib,igrd) 
    1271                bdysurftot = bdysurftot + hv_n   (nbi, nbj  )                           & 
    1272                   &                    * e1v    (nbi, nbj  ) * ABS( flagv ) & 
    1273                   &                    * tmask_i(nbi, nbj  )                           & 
    1274                   &                    * tmask_i(nbi, nbj+1) 
    1275             END DO 
    1276          END DO 
    1277          ! 
    1278          CALL mpp_sum( 'bdyini', bdysurftot )      ! sum over the global domain 
    1279       END IF    
    12801245      ! 
    12811246      ! Tidy up 
     
    12841249      ! 
    12851250   END SUBROUTINE bdy_segs 
    1286  
    12871251 
    12881252   SUBROUTINE bdy_ctl_seg 
  • NEMO/branches/UKMO/dev_r10448_bdyvol/src/OCE/BDY/bdyvol.F90

    r10425 r10455  
    99   !!            3.0  !  2008-04  (NEMO team)  add in the reference version 
    1010   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge 
     11   !!            4.0  !  2019-01  (P. Mathiot) adapted to time splitting      
    1112   !!---------------------------------------------------------------------- 
    1213   USE oce            ! ocean dynamics and tracers  
     
    2425   PRIVATE 
    2526 
    26    PUBLIC   bdy_vol    ! called by ??? 
     27   PUBLIC   bdy_vol2d    ! called by dynspg_ts 
    2728 
    2829   !!---------------------------------------------------------------------- 
     
    3334CONTAINS 
    3435 
    35    SUBROUTINE bdy_vol( kt ) 
     36   SUBROUTINE bdy_vol2d( kt, kc, pua2d, pva2d, phu, phv ) 
    3637      !!---------------------------------------------------------------------- 
    3738      !!                      ***  ROUTINE bdyvol  *** 
     
    4041      !!      A correction velocity is calculated to correct the total transport  
    4142      !!      through the unstructured OBC.  
    42       !!      The total depth used is constant (H0) to be consistent with the  
    43       !!      linear free surface coded in OPA 8.2    <<<=== !!gm  ???? true ???? 
    4443      !! 
    4544      !! ** Method  :   The correction velocity (zubtpecor here) is defined calculating 
     
    6564      !!            (set nn_volctl to 1 in tne namelist for this option) 
    6665      !!---------------------------------------------------------------------- 
    67       INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     66      INTEGER, INTENT(in) ::   kt, kc   ! ocean time-step index, cycle time step 
    6867      ! 
    6968      INTEGER  ::   ji, jj, jk, jb, jgrd 
    7069      INTEGER  ::   ib_bdy, ii, ij 
    71       REAL(wp) ::   zubtpecor, z_cflxemp, ztranst 
     70      REAL(wp) ::   zubtpecor, ztranst, zeh 
     71      REAL(wp), SAVE :: z_cflxemp    
     72      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pua2d, pva2d 
     73      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phu, phv 
    7274      TYPE(OBC_INDEX), POINTER :: idx 
    7375      !!----------------------------------------------------------------------------- 
    7476      ! 
    75       IF( ln_vol ) THEN 
    76       ! 
    77       IF( kt == nit000 ) THEN  
    78          IF(lwp) WRITE(numout,*) 
    79          IF(lwp) WRITE(numout,*)'bdy_vol : Correction of velocities along unstructured OBC' 
    80          IF(lwp) WRITE(numout,*)'~~~~~~~' 
    81       END IF  
    82  
    8377      ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain 
    8478      ! ----------------------------------------------------------------------- 
    85 !!gm replace these lines : 
    86       z_cflxemp = SUM ( ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:) ) / rau0 
    87       CALL mpp_sum( 'bdyvol', z_cflxemp )     ! sum over the global domain 
    88 !!gm   by : 
    89 !!gm      z_cflxemp = glob_sum(  'bdyvol', ( emp(:,:)-rnf(:,:)+fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rau0 
    90 !!gm 
     79      IF ( kc == 1 ) z_cflxemp = glob_sum( 'bdyvol', ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) * bdytmask(:,:) * e1e2t(:,:)  ) / rau0 
     80 
     81      ! Compute bdy surface each cycle if non linear free surface 
     82      ! --------------------------------------------------------- 
     83      IF ( .NOT. ln_linssh ) THEN 
     84         ! compute area each time step 
     85         bdysurftot = bdy_segs_surf( phu, phv ) 
     86      ELSE 
     87         ! compute area only the first time 
     88         IF ( ( kt == nit000 ) .AND. ( kc == 1 ) ) bdysurftot = bdy_segs_surf( phu, phv ) 
     89      END IF 
    9190 
    9291      ! Transport through the unstructured open boundary 
     
    9897         jgrd = 2                               ! cumulate u component contribution first  
    9998         DO jb = 1, idx%nblenrim(jgrd) 
    100             DO jk = 1, jpkm1 
     99            ii = idx%nbi(jb,jgrd) 
     100            ij = idx%nbj(jb,jgrd) 
     101            zubtpecor = zubtpecor + idx%flagu(jb,jgrd) * pua2d(ii,ij) * e2u(ii,ij) * phu(ii,ij) * tmask_i(ii,ij) * tmask_i(ii+1,ij) 
     102         END DO 
     103         jgrd = 3                               ! then add v component contribution 
     104         DO jb = 1, idx%nblenrim(jgrd) 
     105            ii = idx%nbi(jb,jgrd) 
     106            ij = idx%nbj(jb,jgrd) 
     107            zubtpecor = zubtpecor + idx%flagv(jb,jgrd) * pva2d(ii,ij) * e1v(ii,ij) * phv(ii,ij) * tmask_i(ii,ij) * tmask_i(ii,ij+1) 
     108         END DO 
     109         ! 
     110      END DO 
     111      IF( lk_mpp )   CALL mpp_sum( 'bdyvol', zubtpecor )   ! sum over the global domain 
     112 
     113      ! The normal velocity correction 
     114      ! ------------------------------ 
     115      IF( nn_volctl==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp ) / bdysurftot  ! maybe should be apply only once at the end 
     116      ELSE                      ;   zubtpecor =   zubtpecor               / bdysurftot 
     117      END IF 
     118 
     119      ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation 
     120      ! ------------------------------------------------------------- 
     121      DO ib_bdy = 1, nb_bdy 
     122         idx => idx_bdy(ib_bdy) 
     123         ! 
     124         jgrd = 2                               ! correct u component 
     125         DO jb = 1, idx%nblenrim(jgrd) 
    101126               ii = idx%nbi(jb,jgrd) 
    102127               ij = idx%nbj(jb,jgrd) 
    103                zubtpecor = zubtpecor + idx%flagu(jb,jgrd) * ua(ii,ij, jk) * e2u(ii,ij) * e3u_n(ii,ij,jk) 
    104             END DO 
    105          END DO 
    106          jgrd = 3                               ! then add v component contribution 
    107          DO jb = 1, idx%nblenrim(jgrd) 
    108             DO jk = 1, jpkm1 
     128               pua2d(ii,ij) = pua2d(ii,ij) - idx%flagu(jb,jgrd) * zubtpecor * tmask_i(ii,ij) * tmask_i(ii+1,ij) 
     129         END DO 
     130         jgrd = 3                              ! correct v component 
     131         DO jb = 1, idx%nblenrim(jgrd) 
    109132               ii = idx%nbi(jb,jgrd) 
    110133               ij = idx%nbj(jb,jgrd) 
    111                zubtpecor = zubtpecor + idx%flagv(jb,jgrd) * va(ii,ij, jk) * e1v(ii,ij) * e3v_n(ii,ij,jk)  
    112             END DO 
    113          END DO 
    114          ! 
    115       END DO 
    116       CALL mpp_sum( 'bdyvol', zubtpecor )   ! sum over the global domain 
    117  
    118       ! The normal velocity correction 
    119       ! ------------------------------ 
    120       IF( nn_volctl==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp ) / bdysurftot  
    121       ELSE                      ;   zubtpecor =   zubtpecor               / bdysurftot 
    122       END IF 
    123  
    124       ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation 
    125       ! ------------------------------------------------------------- 
    126       ztranst = 0._wp 
    127       DO ib_bdy = 1, nb_bdy 
    128          idx => idx_bdy(ib_bdy) 
    129          ! 
    130          jgrd = 2                               ! correct u component 
    131          DO jb = 1, idx%nblenrim(jgrd) 
    132             DO jk = 1, jpkm1 
    133                ii = idx%nbi(jb,jgrd) 
    134                ij = idx%nbj(jb,jgrd) 
    135                ua(ii,ij,jk) = ua(ii,ij,jk) - idx%flagu(jb,jgrd) * zubtpecor * umask(ii,ij,jk) 
    136                ztranst = ztranst + idx%flagu(jb,jgrd) * ua(ii,ij,jk) * e2u(ii,ij) * e3u_n(ii,ij,jk) 
    137             END DO 
    138          END DO 
    139          jgrd = 3                              ! correct v component 
    140          DO jb = 1, idx%nblenrim(jgrd) 
    141             DO jk = 1, jpkm1 
    142                ii = idx%nbi(jb,jgrd) 
    143                ij = idx%nbj(jb,jgrd) 
    144                va(ii,ij,jk) = va(ii,ij,jk) -idx%flagv(jb,jgrd) * zubtpecor * vmask(ii,ij,jk) 
    145                ztranst = ztranst + idx%flagv(jb,jgrd) * va(ii,ij,jk) * e1v(ii,ij) * e3v_n(ii,ij,jk) 
    146             END DO 
    147          END DO 
    148          ! 
    149       END DO 
    150       CALL mpp_sum( 'bdyvol', ztranst )   ! sum over the global domain 
    151   
     134               pva2d(ii,ij) = pva2d(ii,ij) - idx%flagv(jb,jgrd) * zubtpecor * tmask_i(ii,ij) * tmask_i(ii,ij+1) 
     135         END DO 
     136         ! 
     137      END DO 
     138      !  
    152139      ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected 
    153140      ! ------------------------------------------------------ 
    154       IF( lwp .AND. MOD( kt, nwrite ) == 0 ) THEN 
     141      IF( MOD( kt, nwrite ) == 0 .AND. ( kc == 1 ) ) THEN 
     142         ! 
     143         ! compute residual transport across boundary 
     144         ztranst = 0._wp 
     145         DO ib_bdy = 1, nb_bdy 
     146            idx => idx_bdy(ib_bdy) 
     147            ! 
     148            jgrd = 2                               ! correct u component 
     149            DO jb = 1, idx%nblenrim(jgrd) 
     150                  ii = idx%nbi(jb,jgrd) 
     151                  ij = idx%nbj(jb,jgrd) 
     152                  ztranst = ztranst + idx%flagu(jb,jgrd) * pua2d(ii,ij) * e2u(ii,ij) * phu(ii,ij) * tmask_i(ii,ij) * tmask_i(ii+1,ij) 
     153            END DO 
     154            jgrd = 3                              ! correct v component 
     155            DO jb = 1, idx%nblenrim(jgrd) 
     156                  ii = idx%nbi(jb,jgrd) 
     157                  ij = idx%nbj(jb,jgrd) 
     158                  ztranst = ztranst + idx%flagv(jb,jgrd) * pva2d(ii,ij) * e1v(ii,ij) * phv(ii,ij) * tmask_i(ii,ij) * tmask_i(ii,ij+1) 
     159            END DO 
     160            ! 
     161         END DO 
     162         IF( lk_mpp )   CALL mpp_sum('bdyvol', ztranst )   ! sum over the global domain 
     163 
     164 
    155165         IF(lwp) WRITE(numout,*) 
    156166         IF(lwp) WRITE(numout,*)'bdy_vol : time step :', kt 
     
    162172      END IF  
    163173      ! 
    164       END IF ! ln_vol 
    165       ! 
    166    END SUBROUTINE bdy_vol 
    167  
     174   END SUBROUTINE bdy_vol2d 
     175   ! 
     176   REAL(wp) FUNCTION bdy_segs_surf(phu, phv) 
     177      !!---------------------------------------------------------------------- 
     178      !!                 ***  ROUTINE bdy_ctl_seg  *** 
     179      !! 
     180      !! ** Purpose :   Compute total lateral surface for volume correction 
     181      !! 
     182      !!---------------------------------------------------------------------- 
     183 
     184      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: phu, phv ! water column thickness at U and V points 
     185      INTEGER            ::  igrd, ib_bdy, ib                ! loop indexes 
     186      INTEGER , POINTER  ::  nbi, nbj                        ! short cuts 
     187      REAL(wp), POINTER  ::  zflagu, zflagv                  !    -   - 
     188 
     189      ! Compute total lateral surface for volume correction: 
     190      ! ---------------------------------------------------- 
     191      bdy_segs_surf = 0._wp 
     192      igrd = 2      ! Lateral surface at U-points 
     193      DO ib_bdy = 1, nb_bdy 
     194         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     195            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     196            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
     197            zflagu => idx_bdy(ib_bdy)%flagu(ib,igrd) 
     198            bdy_segs_surf = bdy_segs_surf + phu(nbi, nbj)                              & 
     199               &                            * e2u(nbi, nbj) * ABS( zflagu )            & 
     200               &                            * tmask_i(nbi, nbj) * tmask_i(nbi+1, nbj) 
     201         END DO 
     202      END DO 
     203 
     204      igrd=3 ! Add lateral surface at V-points 
     205      DO ib_bdy = 1, nb_bdy 
     206         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     207            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd) 
     208            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd) 
     209            zflagv => idx_bdy(ib_bdy)%flagv(ib,igrd) 
     210            bdy_segs_surf = bdy_segs_surf + phv(nbi, nbj)                              & 
     211               &                            * e1v(nbi, nbj) * ABS( zflagv )            & 
     212               &                            * tmask_i(nbi, nbj) * tmask_i(nbi, nbj+1) 
     213         END DO 
     214      END DO 
     215      ! 
     216      ! redirect the time to bdyvol as this variable is only used by bdyvol 
     217      IF( lk_mpp )   CALL mpp_sum( 'bdyvol', bdy_segs_surf ) ! sum over the global domain 
     218      ! 
     219   END FUNCTION bdy_segs_surf 
    168220   !!====================================================================== 
    169221END MODULE bdyvol 
  • NEMO/branches/UKMO/dev_r10448_bdyvol/src/OCE/DYN/dynspg_ts.F90

    r10425 r10455  
    4141   USE wet_dry         ! wetting/drying flux limter 
    4242   USE bdy_oce         ! open boundary 
     43   USE bdyvol          ! open boundary volume conservation 
    4344   USE bdytides        ! open boundary condition data 
    4445   USE bdydyn2d        ! open boundary conditions on barotropic variables 
     
    792793         !                                                !* after ssh 
    793794         !                                                !  ----------- 
    794          ! One should enforce volume conservation at open boundaries here 
    795          ! considering fluxes below: 
     795         ! 
     796         ! Enforce volume conservation at open boundaries: 
     797         IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e ) 
    796798         ! 
    797799         zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:)         ! fluxes at jn+0.5 
Note: See TracChangeset for help on using the changeset viewer.