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 3970 for branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

Ignore:
Timestamp:
2013-07-11T15:59:14+02:00 (11 years ago)
Author:
cbricaud
Message:

Time splitting update, see ticket #1079

Location:
branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r3851 r3970  
    2929   USE iom             ! IOM library 
    3030   USE in_out_manager  ! I/O logical units 
     31   ! bg jchanut tschanges 
     32   USE dynspg_oce, ONLY: lk_dynspg_ts ! Split-explicit free surface flag 
     33   ! end jchanut tschanges 
     34 
    3135#if defined key_lim2 
    3236   USE ice_2 
     
    314318                     END DO 
    315319                  ENDIF 
    316                   IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
    317                      CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy),  & 
    318                                         & td=tides(ib_bdy), time_offset=time_offset ) 
    319                   ENDIF 
     320                  ! bg jchanut tschanges 
     321                  !IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
     322                  !   CALL bdytide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy),  & 
     323                  !                      & td=tides(ib_bdy), time_offset=time_offset ) 
     324                  !ENDIF 
     325                  ! end jchanut tschanges 
    320326               ENDIF 
    321327            ENDIF 
     
    323329         END IF ! nn_dta(ib_bdy) = 1 
    324330      END DO  ! ib_bdy 
     331 
     332      ! bg jchanut tschanges 
     333#if defined key_tide 
     334      ! Add tides if not split-explicit free surface else this is done in ts loop 
     335      IF (.NOT.lk_dynspg_ts) CALL bdy_dta_tides( kt=kt, time_offset=time_offset ) 
     336#endif 
     337      ! end jchanut tschanges 
    325338 
    326339      IF ( ln_apr_obc ) THEN 
     
    476489            IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN  
    477490 
    478                IF( nn_tra(ib_bdy) .ne. 4 ) THEN ! runoff condition : no ssh reading 
     491               IF( nn_dyn2d(ib_bdy) .ne. jp_frs .and. nn_tra(ib_bdy) .ne. 4 ) THEN ! runoff condition : no ssh reading 
    479492                  jfld = jfld + 1 
    480493                  blf_i(jfld) = bn_ssh 
     
    572585            ! Recalculate field counts 
    573586            !------------------------- 
    574             nb_bdy_fld_sum = 0 
    575587            IF( ib_bdy .eq. 1 ) THEN  
     588               nb_bdy_fld_sum = 0 
    576589               nb_bdy_fld(ib_bdy) = jfld 
    577590               nb_bdy_fld_sum     = jfld               
     
    616629               ALLOCATE( dta_bdy(ib_bdy)%u2d(ilen0(2)) ) 
    617630               ALLOCATE( dta_bdy(ib_bdy)%v2d(ilen0(3)) ) 
    618                IF (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3) THEN 
     631               IF ( nn_dyn2d(ib_bdy) .ne. jp_frs .and. (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3) ) THEN 
    619632                  jfld = jfld + 1 
    620633                  dta_bdy(ib_bdy)%ssh => bf(jfld)%fnow(:,1,1) 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn.F90

    r3294 r3970  
    3030   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    3131   USE in_out_manager  ! 
     32   USE domvvl          ! variable volume 
    3233 
    3334   IMPLICIT NONE 
     
    8485      pu2d(:,:) = 0.e0 
    8586      pv2d(:,:) = 0.e0 
    86       DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
    87           pu2d(:,:) = pu2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
    88           pv2d(:,:) = pv2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
    89       END DO 
    90       pu2d(:,:) = pu2d(:,:) * phur(:,:) 
    91       pv2d(:,:) = pv2d(:,:) * phvr(:,:) 
     87      ! bg jchanut tschanges (not specifically related to ts; this is a bug) 
     88      IF (lk_vvl) THEN 
     89         DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
     90            pu2d(:,:) = pu2d(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
     91            pv2d(:,:) = pv2d(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
     92         END DO 
     93         pu2d(:,:) = pu2d(:,:) / ( hu_0(:,:) + sshu_a(:,:) + 1._wp - umask(:,:,1) ) 
     94         pv2d(:,:) = pv2d(:,:) / ( hv_0(:,:) + sshv_a(:,:) + 1._wp - vmask(:,:,1) ) 
     95      ! end jchanut tschanges 
     96      ELSE 
     97         DO jk = 1, jpkm1   !! Vertically integrated momentum trends 
     98            pu2d(:,:) = pu2d(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 
     99            pv2d(:,:) = pv2d(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 
     100         END DO 
     101         pu2d(:,:) = pu2d(:,:) * phur(:,:) 
     102         pv2d(:,:) = pv2d(:,:) * phvr(:,:) 
     103      ENDIF 
     104 
    92105      DO jk = 1 , jpkm1 
    93          ua(:,:,jk) = ua(:,:,jk) - pu2d(:,:) 
    94          va(:,:,jk) = va(:,:,jk) - pv2d(:,:) 
     106         ua(:,:,jk) = ua(:,:,jk) - pu2d(:,:) * umask(:,:,jk) 
     107         va(:,:,jk) = va(:,:,jk) - pv2d(:,:) * vmask(:,:,jk) 
    95108      END DO 
    96109 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn2d.F90

    r3680 r3970  
    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  
     
    2728 
    2829   PUBLIC   bdy_dyn2d     ! routine called in dynspg_ts and bdy_dyn 
     30   PUBLIC   bdy_ssh       ! routine called in dynspg_ts or sshwzv 
    2931 
    3032   !!---------------------------------------------------------------------- 
     
    135137      REAL(wp) ::   zcorr                            ! Flather correction 
    136138      REAL(wp) ::   zforc                            ! temporary scalar 
     139      REAL(wp) ::   zflag, z1_2                      !    "        " 
    137140      !!---------------------------------------------------------------------- 
    138141 
    139142      IF( nn_timing == 1 ) CALL timing_start('bdy_dyn2d_fla') 
     143 
     144      z1_2 = 0.5_wp 
    140145 
    141146      ! ---------------------------------! 
     
    164169         ! 
    165170         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)  
     171         ! bg jchanut tschanges: Set zflag to 0 below to revert to Flather scheme 
     172!!         zforc = dta%u2d(jb) 
     173         zflag = ABS(idx%flagu(jb)) 
     174         iim1 = ii + idx%flagu(jb) 
     175         zforc = dta%u2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pu2d(iim1,ij) 
     176         pu2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * umask(ii,ij,1)  
     177         ! end jchanut tschanges 
    168178      END DO 
    169179      ! 
     
    177187         ! 
    178188         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) 
     189         ! bg jchanut tschanges: Set zflag to 0 below to revert to std Flather scheme 
     190!!         zforc = dta%v2d(jb) 
     191         zflag = ABS(idx%flagv(jb)) 
     192         ijm1 = ij + idx%flagv(jb) 
     193         zforc = dta%v2d(jb) * (1._wp - z1_2*zflag) + z1_2 * zflag * pv2d(ii,ijm1) 
     194         pv2d(ii,ij) = zforc + (1._wp - z1_2*zflag) * zcorr * vmask(ii,ij,1) 
     195         ! end jchanut tschanges 
    181196      END DO 
    182197      CALL lbc_bdy_lnk( pu2d, 'U', -1., ib_bdy )   ! Boundary points should be updated 
     
    186201      ! 
    187202   END SUBROUTINE bdy_dyn2d_fla 
     203 
     204   SUBROUTINE bdy_ssh( zssh ) 
     205      !!---------------------------------------------------------------------- 
     206      !!                  ***  SUBROUTINE bdy_ssh  *** 
     207      !! 
     208      !! ** Purpose : Duplicate sea level across open boundaries 
     209      !! 
     210      !!---------------------------------------------------------------------- 
     211      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   zssh ! Sea level 
     212      !! 
     213      INTEGER  ::   ib_bdy, ib, igrd                        ! local integers 
     214      INTEGER  ::   ii, ij, zcoef, zcoef1, zcoef2, ip, jp   !   "       " 
     215 
     216      igrd = 1                       ! Everything is at T-points here 
     217 
     218      DO ib_bdy = 1, nb_bdy 
     219         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
     220            ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     221            ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 
     222            ! Set gradient direction: 
     223            zcoef1 = bdytmask(ii-1,ij  ) +  bdytmask(ii+1,ij  ) 
     224            zcoef2 = bdytmask(ii  ,ij-1) +  bdytmask(ii  ,ij+1) 
     225            IF ( zcoef1+zcoef2 == 0 ) THEN 
     226               ! corner 
     227!               zcoef = tmask(ii-1,ij,1) + tmask(ii+1,ij,1) +  tmask(ii,ij-1,1) +  tmask(ii,ij+1,1) 
     228!               zssh(ii,ij) = zssh(ii-1,ij  ) * tmask(ii-1,ij  ,1) + & 
     229!                 &           zssh(ii+1,ij  ) * tmask(ii+1,ij  ,1) + & 
     230!                 &           zssh(ii  ,ij-1) * tmask(ii  ,ij-1,1) + & 
     231!                 &           zssh(ii  ,ij+1) * tmask(ii  ,ij+1,1) 
     232               zcoef = bdytmask(ii-1,ij) + bdytmask(ii+1,ij) +  bdytmask(ii,ij-1) +  bdytmask(ii,ij+1) 
     233               zssh(ii,ij) = zssh(ii-1,ij  ) * bdytmask(ii-1,ij  ) + & 
     234                 &           zssh(ii+1,ij  ) * bdytmask(ii+1,ij  ) + & 
     235                 &           zssh(ii  ,ij-1) * bdytmask(ii  ,ij-1) + & 
     236                 &           zssh(ii  ,ij+1) * bdytmask(ii  ,ij+1) 
     237               zssh(ii,ij) = ( zssh(ii,ij) / MAX( 1, zcoef) ) * tmask(ii,ij,1) 
     238            ELSE 
     239               ip = bdytmask(ii+1,ij  ) - bdytmask(ii-1,ij  ) 
     240               jp = bdytmask(ii  ,ij+1) - bdytmask(ii  ,ij-1) 
     241               zssh(ii,ij) = zssh(ii+ip,ij+jp) * tmask(ii+ip,ij+jp,1) 
     242            ENDIF 
     243         END DO 
     244 
     245         ! Boundary points should be updated 
     246         CALL lbc_bdy_lnk( zssh(:,:), 'T', 1., ib_bdy ) 
     247      END DO 
     248 
     249   END SUBROUTINE bdy_ssh 
    188250#else 
    189251   !!---------------------------------------------------------------------- 
     
    192254CONTAINS 
    193255   SUBROUTINE bdy_dyn2d( kt )      ! Empty routine 
    194       WRITE(*,*) 'bdy_dyn_frs: You should not have seen this print! error?', kt 
     256      INTEGER, intent(in) :: kt 
     257      WRITE(*,*) 'bdy_dyn2: You should not have seen this print! error?', kt 
    195258   END SUBROUTINE bdy_dyn2d 
     259 
    196260#endif 
    197261 
    198262   !!====================================================================== 
    199263END MODULE bdydyn2d 
     264 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90

    r3703 r3970  
    161161 
    162162      DO ib_bdy = 1,nb_bdy 
     163        icount = 0 ! flag to set max rimwidth to 1 if no relaxation 
    163164        IF(lwp) WRITE(numout,*) ' '  
    164165        IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------'  
     
    175176          CASE(jp_none)         ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    176177          CASE(jp_frs)          ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     178          icount = icount + 1 
    177179          CASE(jp_flather)      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition' 
    178180          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' ) 
     
    196198          CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    197199          CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     200          icount = icount + 1 
    198201          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value' 
    199202          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)' 
     
    215218              CALL ctl_stop( 'Use FRS OR relaxation' ) 
    216219           ELSE 
     220              icount = icount + 1 
    217221              IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone' 
    218222              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
     
    228232          CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    229233          CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     234          icount = icount + 1 
    230235          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value' 
    231236          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Neumann conditions' 
     
    248253              CALL ctl_stop( 'Use FRS OR relaxation' ) 
    249254           ELSE 
     255              icount = icount + 1 
    250256              IF(lwp) WRITE(numout,*) '      + T/S relaxation zone' 
    251257              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days' 
     
    262268          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'         
    263269          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme' 
     270          icount = icount + 1 
    264271          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' ) 
    265272        END SELECT 
     
    273280        IF(lwp) WRITE(numout,*) 
    274281#endif 
    275  
    276         IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy) 
    277         IF(lwp) WRITE(numout,*) 
     282        IF ( icount>0 ) THEN 
     283           IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy) 
     284           IF(lwp) WRITE(numout,*) 
     285        ELSE 
     286           nn_rimwidth(ib_bdy) = 1 ! no relaxation 
     287        ENDIF 
    278288 
    279289      ENDDO 
     
    391401            ENDDO 
    392402            CALL iom_close( inum ) 
    393  
    394403         ENDIF  
    395404 
     
    398407      IF (nb_bdy>0) THEN 
    399408         jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy)) 
    400  
    401409         ! Allocate arrays 
    402410         !--------------- 
     
    446454         ENDIF  
    447455 
    448       ENDDO       
     456      ENDDO      
    449457     
    450458      ! 2. Now fill indices corresponding to straight open boundary arrays: 
     
    752760               ! check if point is in local domain 
    753761               IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   & 
    754                   & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in       ) THEN 
     762                  & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND.   & 
     763                  & nbrdta(ib,igrd,ib_bdy) <= nn_rimwidth(ib_bdy)     ) THEN       
    755764                  ! 
    756765                  icount = icount  + 1 
     
    765774         ! Allocate index arrays for this boundary set 
    766775         !-------------------------------------------- 
    767          ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:)) 
     776 
     777         ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(1:jpbgrd)) 
     778         ilen1 = MAX(1,ilen1) 
    768779         ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) ) 
    769780         ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) ) 
     
    773784         ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) ) 
    774785         ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1) ) 
    775          ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) ) 
     786         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) )       
    776787 
    777788         ! Dispatch mapping indices and discrete distances on each processor 
     
    10081019      ! bdytmask = 1  on the computational domain AND on open boundaries 
    10091020      !          = 0  elsewhere    
    1010   
     1021      bdytmask(:,:) = 1.e0 
     1022      bdyumask(:,:) = 1.e0 
     1023      bdyvmask(:,:) = 1.e0 
     1024 
    10111025      IF( ln_mask_file ) THEN 
    10121026         CALL iom_open( cn_mask_file, inum ) 
  • branches/2013/dev_r3867_MERCATOR1_DYN/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90

    r3651 r3970  
    99   !!            3.3  !  2010-09  (D.Storkey and E.O'Dea)  bug fixes 
    1010   !!            3.4  !  2012-09  (G. Reffray and J. Chanut) New inputs + mods 
     11   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes 
    1112   !!---------------------------------------------------------------------- 
    1213#if defined key_bdy 
     
    3233!   USE tide_mod       ! Useless ?? 
    3334   USE fldread, ONLY: fld_map 
     35   USE dynspg_oce, ONLY: lk_dynspg_ts 
    3436 
    3537   IMPLICIT NONE 
     
    3840   PUBLIC   bdytide_init     ! routine called in bdy_init 
    3941   PUBLIC   bdytide_update   ! routine called in bdy_dta 
     42   PUBLIC   bdy_dta_tides    ! routine called in dyn_spg_ts 
    4043 
    4144   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data 
     
    4952 
    5053   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data 
     54   TYPE(OBC_DATA)  , PRIVATE, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component) 
    5155 
    5256   !!---------------------------------------------------------------------- 
     
    252256            ENDIF 
    253257            ! 
     258            IF ( lk_dynspg_ts ) THEN ! Allocate arrays to save slowly varying boundary data during 
     259                                     ! time splitting integration 
     260               ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) ) 
     261               ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) ) 
     262               ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) ) 
     263               dta_bdy_s(ib_bdy)%ssh(:) = 0.e0 
     264               dta_bdy_s(ib_bdy)%u2d(:) = 0.e0 
     265               dta_bdy_s(ib_bdy)%v2d(:) = 0.e0 
     266            ENDIF 
     267            ! 
    254268         ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2 
    255269         ! 
     
    318332          
    319333      IF( PRESENT(jit) ) THEN   
    320          z_arg = ( ((kt-kt_tide)-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) ) 
     334         z_arg = ((kt-kt_tide) * rdt + (jit+0.5*(time_add-1)) * rdt / REAL(nn_baro,wp) ) 
    321335      ELSE                               
    322336         z_arg = ((kt-kt_tide)+time_add) * rdt 
     
    324338 
    325339      ! Linear ramp on tidal component at open boundaries  
    326       zramp = 1. 
    327       IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0.),1.) 
     340      zramp = 1._wp 
     341      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0._wp),1._wp) 
    328342 
    329343      DO itide = 1, nb_harmo 
     
    351365      ! 
    352366   END SUBROUTINE bdytide_update 
     367 
     368   SUBROUTINE bdy_dta_tides( kt, kit, time_offset ) 
     369      !!---------------------------------------------------------------------- 
     370      !!                 ***  SUBROUTINE bdy_dta_tides  *** 
     371      !!                 
     372      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.  
     373      !!                 
     374      !!---------------------------------------------------------------------- 
     375      INTEGER, INTENT( in )            ::   kt          ! Main timestep counter 
     376      INTEGER, INTENT( in ),OPTIONAL   ::   kit         ! Barotropic timestep counter (for timesplitting option) 
     377      INTEGER, INTENT( in ),OPTIONAL   ::   time_offset ! time offset in units of timesteps. NB. if kit 
     378                                                        ! is present then units = subcycle timesteps. 
     379                                                        ! time_offset = 0  => get data at "now"    time level 
     380                                                        ! time_offset = -1 => get data at "before" time level 
     381                                                        ! time_offset = +1 => get data at "after"  time level 
     382                                                        ! etc. 
     383      !! 
     384      LOGICAL  :: lk_first_btstp  ! =.TRUE. if time splitting and first barotropic step 
     385      INTEGER,          DIMENSION(jpbgrd) :: ilen0  
     386      INTEGER, DIMENSION(1:jpbgrd) :: nblen, nblenrim  ! short cuts 
     387      INTEGER  :: itide, ib_bdy, ib, igrd                     ! loop indices 
     388      INTEGER  :: time_add                                    ! time offset in units of timesteps 
     389      REAL(wp) :: z_arg, z_sarg, zramp, zoff, z_cost, z_sist       
     390      !!---------------------------------------------------------------------- 
     391 
     392      IF( nn_timing == 1 ) CALL timing_start('bdy_dta_tides') 
     393 
     394      lk_first_btstp=.TRUE. 
     395      IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF 
     396 
     397      time_add = 0 
     398      IF( PRESENT(time_offset) ) THEN 
     399         time_add = time_offset 
     400      ENDIF 
     401       
     402      ! Absolute time from model initialization:    
     403      IF( PRESENT(kit) ) THEN   
     404         z_arg = ( kt + (kit+0.5*(time_add-1)) / REAL(nn_baro,wp) ) * rdt 
     405      ELSE                               
     406         z_arg = ( kt + time_add ) * rdt 
     407      ENDIF 
     408 
     409      ! Linear ramp on tidal component at open boundaries  
     410      zramp = 1. 
     411      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rdttideramp*rday),0.),1.) 
     412 
     413      DO ib_bdy = 1,nb_bdy 
     414 
     415         ! line below should be simplified (runoff case) 
     416         IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(nn_tra(ib_bdy).NE.4)) THEN 
     417 
     418            nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd) 
     419            nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd) 
     420 
     421            IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 
     422               ilen0(:)=nblen(:) 
     423            ELSE 
     424               ilen0(:)=nblenrim(:) 
     425            ENDIF      
     426 
     427            ! We refresh nodal factors every day below 
     428            ! This should be done somewhere else 
     429            IF ( nsec_day == NINT(0.5 * rdttra(1)) .AND. lk_first_btstp ) THEN 
     430               ! 
     431               kt_tide = kt                
     432               ! 
     433               IF(lwp) THEN 
     434               WRITE(numout,*) 
     435               WRITE(numout,*) 'bdy_tide_dta : Refresh nodal factors for tidal open bdy data at kt=',kt 
     436               WRITE(numout,*) '~~~~~~~~~~~~~~ ' 
     437               ENDIF 
     438               ! 
     439               CALL tide_init_elevation ( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) ) 
     440               CALL tide_init_velocities( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) ) 
     441               ! 
     442            ENDIF 
     443            zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time 
     444            ! 
     445            ! If time splitting, save data at first barotropic iteration 
     446            IF ( PRESENT(kit) ) THEN 
     447               IF ( lk_first_btstp ) THEN ! Save slow varying open boundary data: 
     448                  dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy(ib_bdy)%ssh(1:ilen0(1)) 
     449                  dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy(ib_bdy)%u2d(1:ilen0(2)) 
     450                  dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy(ib_bdy)%v2d(1:ilen0(3)) 
     451 
     452               ELSE ! Initialize arrays from slow varying open boundary data:             
     453                  dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) 
     454                  dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) 
     455                  dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) 
     456               ENDIF 
     457            ENDIF 
     458            ! 
     459            ! Update open boundary data arrays: 
     460            DO itide = 1, nb_harmo 
     461               ! 
     462               z_sarg = (z_arg + zoff) * omega_tide(itide) 
     463               z_cost = zramp * COS( z_sarg ) 
     464               z_sist = zramp * SIN( z_sarg ) 
     465               ! 
     466               igrd=1                              ! SSH on tracer grid 
     467               DO ib = 1, ilen0(igrd) 
     468                  dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + & 
     469                     &                      ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + & 
     470                     &                        tides(ib_bdy)%ssh(ib,itide,2)*z_sist ) 
     471               END DO 
     472               ! 
     473               igrd=2                              ! U grid 
     474               DO ib = 1, ilen0(igrd) 
     475                  dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + & 
     476                     &                      ( tides(ib_bdy)%u(ib,itide,1)*z_cost + & 
     477                     &                        tides(ib_bdy)%u(ib,itide,2)*z_sist ) 
     478               END DO 
     479               ! 
     480               igrd=3                              ! V grid 
     481               DO ib = 1, ilen0(igrd)  
     482                  dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + & 
     483                     &                      ( tides(ib_bdy)%v(ib,itide,1)*z_cost + & 
     484                     &                        tides(ib_bdy)%v(ib,itide,2)*z_sist ) 
     485               END DO 
     486            END DO 
     487         END IF 
     488      END DO 
     489      ! 
     490      IF( nn_timing == 1 ) CALL timing_stop('bdy_dta_tides') 
     491      ! 
     492   END SUBROUTINE bdy_dta_tides 
    353493 
    354494   SUBROUTINE tide_init_elevation( idx, td ) 
     
    457597      WRITE(*,*) 'bdytide_update: You should not have seen this print! error?', kt, jit 
    458598   END SUBROUTINE bdytide_update 
     599   SUBROUTINE bdy_dta_tides( kt, jit )   ! Empty routine 
     600      WRITE(*,*) 'bdy_dta_tides: You should not have seen this print! error?', kt, jit 
     601   END SUBROUTINE bdy_dta_tides 
    459602#endif 
    460603 
    461604   !!====================================================================== 
    462605END MODULE bdytides 
     606 
Note: See TracChangeset for help on using the changeset viewer.