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 12059 for NEMO/branches – NEMO

Changeset 12059 for NEMO/branches


Ignore:
Timestamp:
2019-12-05T09:34:39+01:00 (4 years ago)
Author:
smueller
Message:

Reversal of changeset [10865] to avoid a conflict in the upcoming sync merge with the trunk; a new version of the modification implemented by changeset [10865] will be applied following the upcoming sync merge with the trunk (ticket #2194)

Location:
NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/BDY
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/BDY/bdydta.F90

    r10865 r12059  
    6565CONTAINS 
    6666 
    67       SUBROUTINE bdy_dta( kt, time_offset ) 
     67      SUBROUTINE bdy_dta( kt, jit, time_offset ) 
    6868      !!---------------------------------------------------------------------- 
    6969      !!                   ***  SUBROUTINE bdy_dta  *** 
    … …  
    7575      !!---------------------------------------------------------------------- 
    7676      INTEGER, INTENT(in)           ::   kt           ! ocean time-step index  
    77       INTEGER, INTENT(in), OPTIONAL ::   time_offset  ! time offset in units of timesteps. 
     77      INTEGER, INTENT(in), OPTIONAL ::   jit          ! subcycle time-step index (for timesplitting option) 
     78      INTEGER, INTENT(in), OPTIONAL ::   time_offset  ! time offset in units of timesteps. NB. if jit 
     79      !                                               ! is present then units = subcycle timesteps. 
    7880      !                                               ! time_offset = 0 => get data at "now" time level 
    7981      !                                               ! time_offset = -1 => get data at "before" time level 
    … …  
    9294      ! Initialise data arrays once for all from initial conditions where required 
    9395      !--------------------------------------------------------------------------- 
    94       IF( kt == nit000 ) THEN 
     96      IF( kt == nit000 .AND. .NOT.PRESENT(jit) ) THEN 
    9597 
    9698         ! Calculate depth-mean currents 
    … …  
    226228         IF( nn_dta(jbdy) == 1 ) THEN ! skip this bit if no external data required 
    227229       
    228             IF (cn_tra(jbdy) == 'runoff') then      ! runoff condition 
    229                jend = nb_bdy_fld(jbdy) 
    230                CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend),  & 
    231                   & map=nbmap_ptr(jstart:jend), kt_offset=time_offset ) 
    232                ! 
    233                igrd = 2                      ! zonal velocity 
    234                DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    235                   ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
    236                   ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    237                   dta%u2d(ib) = dta%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) ) 
    238                END DO 
    239                ! 
    240                igrd = 3                      ! meridional velocity 
    241                DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    242                   ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
    243                   ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    244                   dta%v2d(ib) = dta%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) ) 
    245                END DO 
     230            IF( PRESENT(jit) ) THEN 
     231               ! Update barotropic boundary conditions only 
     232               ! jit is optional argument for fld_read and bdytide_update 
     233               IF( cn_dyn2d(jbdy) /= 'none' ) THEN 
     234                  IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
     235                     IF( dta%ll_ssh ) dta%ssh(:) = 0._wp 
     236                     IF( dta%ll_u2d ) dta%u2d(:) = 0._wp 
     237                     IF( dta%ll_u3d ) dta%v2d(:) = 0._wp 
     238                  ENDIF 
     239                  IF (cn_tra(jbdy) /= 'runoff') THEN 
     240                     IF( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 ) THEN 
     241 
     242                        jend = jstart + dta%nread(2) - 1 
     243                        IF( ln_full_vel_array(jbdy) ) THEN 
     244                           CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  & 
     245                                     & kit=jit, kt_offset=time_offset , jpk_bdy=nb_jpk_bdy,   & 
     246                                     & fvl=ln_full_vel_array(jbdy)  ) 
     247                        ELSE 
     248                           CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend),  & 
     249                                     & kit=jit, kt_offset=time_offset  ) 
     250                        ENDIF 
     251 
     252                        ! If full velocities in boundary data then extract barotropic velocities from 3D fields 
     253                        IF( ln_full_vel_array(jbdy) .AND.                                             & 
     254                          &    ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR.  & 
     255                          &      nn_dyn3d_dta(jbdy) == 1 ) )THEN 
     256 
     257                           igrd = 2                      ! zonal velocity 
     258                           dta%u2d(:) = 0._wp 
     259                           DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
     260                              ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
     261                              ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
     262                              DO ik = 1, jpkm1 
     263                                 dta%u2d(ib) = dta%u2d(ib) & 
     264                       &                          + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik) 
     265                              END DO 
     266                              dta%u2d(ib) =  dta%u2d(ib) * r1_hu_n(ii,ij) 
     267                           END DO 
     268                           igrd = 3                      ! meridional velocity 
     269                           dta%v2d(:) = 0._wp 
     270                           DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
     271                              ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
     272                              ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
     273                              DO ik = 1, jpkm1 
     274                                 dta%v2d(ib) = dta%v2d(ib) & 
     275                       &                       + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik) 
     276                              END DO 
     277                              dta%v2d(ib) =  dta%v2d(ib) * r1_hv_n(ii,ij) 
     278                           END DO 
     279                        ENDIF                     
     280                     ENDIF 
     281                     IF( nn_dyn2d_dta(jbdy) .ge. 2 ) THEN ! update tidal harmonic forcing 
     282                        CALL bdytide_update( kt=kt, idx=idx_bdy(jbdy), dta=dta, td=tides(jbdy),   &  
     283                          &                 jit=jit, time_offset=time_offset ) 
     284                     ENDIF 
     285                  ENDIF 
     286               ENDIF 
    246287            ELSE 
    247                IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
    248                   IF( dta%ll_ssh ) dta%ssh(:) = 0._wp 
    249                   IF( dta%ll_u2d ) dta%u2d(:) = 0._wp 
    250                   IF( dta%ll_v2d ) dta%v2d(:) = 0._wp 
    251                ENDIF 
    252                IF( dta%nread(1) .gt. 0 ) THEN ! update external data 
    253                   jend = jstart + dta%nread(1) - 1 
    254                   CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), & 
    255                      & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy,   & 
    256                      & fvl=ln_full_vel_array(jbdy) ) 
    257                ENDIF 
    258                ! If full velocities in boundary data then split into barotropic and baroclinic data 
    259                IF( ln_full_vel_array(jbdy) .and.                                             & 
    260                   & ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR. & 
    261                   &   nn_dyn3d_dta(jbdy) == 1 ) ) THEN 
     288               IF (cn_tra(jbdy) == 'runoff') then      ! runoff condition 
     289                  jend = nb_bdy_fld(jbdy) 
     290                  CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend),  & 
     291                               & map=nbmap_ptr(jstart:jend), kt_offset=time_offset ) 
     292                  ! 
    262293                  igrd = 2                      ! zonal velocity 
    263                   dta%u2d(:) = 0._wp 
    264294                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    265295                     ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
    266296                     ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    267                      DO ik = 1, jpkm1 
    268                         dta%u2d(ib) = dta%u2d(ib) & 
    269               &                       + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik) 
    270                      END DO 
    271                      dta%u2d(ib) =  dta%u2d(ib) * r1_hu_n(ii,ij) 
    272                      DO ik = 1, jpkm1 
    273                         dta%u3d(ib,ik) = dta%u3d(ib,ik) - dta%u2d(ib) 
    274                      END DO 
     297                     dta%u2d(ib) = dta%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) ) 
    275298                  END DO 
     299                  ! 
    276300                  igrd = 3                      ! meridional velocity 
    277                   dta%v2d(:) = 0._wp 
    278301                  DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
    279302                     ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
    280303                     ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
    281                      DO ik = 1, jpkm1 
    282                         dta%v2d(ib) = dta%v2d(ib) & 
    283               &                       + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik) 
     304                     dta%v2d(ib) = dta%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) ) 
     305                  END DO 
     306               ELSE 
     307                  IF( nn_dyn2d_dta(jbdy) == 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 
     308                     IF( dta%ll_ssh ) dta%ssh(:) = 0._wp 
     309                     IF( dta%ll_u2d ) dta%u2d(:) = 0._wp 
     310                     IF( dta%ll_v2d ) dta%v2d(:) = 0._wp 
     311                  ENDIF 
     312                  IF( dta%nread(1) .gt. 0 ) THEN ! update external data 
     313                     jend = jstart + dta%nread(1) - 1 
     314                     CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), & 
     315                                  & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy,   & 
     316                                  & fvl=ln_full_vel_array(jbdy) ) 
     317                  ENDIF 
     318                  ! If full velocities in boundary data then split into barotropic and baroclinic data 
     319                  IF( ln_full_vel_array(jbdy) .and.                                             & 
     320                    & ( nn_dyn2d_dta(jbdy) == 1 .OR. nn_dyn2d_dta(jbdy) == 3 .OR. & 
     321                    &   nn_dyn3d_dta(jbdy) == 1 ) ) THEN 
     322                     igrd = 2                      ! zonal velocity 
     323                     dta%u2d(:) = 0._wp 
     324                     DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
     325                        ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
     326                        ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
     327                        DO ik = 1, jpkm1 
     328                           dta%u2d(ib) = dta%u2d(ib) & 
     329                 &                       + e3u_n(ii,ij,ik) * umask(ii,ij,ik) * dta%u3d(ib,ik) 
     330                        END DO 
     331                        dta%u2d(ib) =  dta%u2d(ib) * r1_hu_n(ii,ij) 
     332                        DO ik = 1, jpkm1 
     333                           dta%u3d(ib,ik) = dta%u3d(ib,ik) - dta%u2d(ib) 
     334                        END DO 
    284335                     END DO 
    285                      dta%v2d(ib) =  dta%v2d(ib) * r1_hv_n(ii,ij) 
    286                      DO ik = 1, jpkm1 
    287                         dta%v3d(ib,ik) = dta%v3d(ib,ik) - dta%v2d(ib) 
     336                     igrd = 3                      ! meridional velocity 
     337                     dta%v2d(:) = 0._wp 
     338                     DO ib = 1, idx_bdy(jbdy)%nblen(igrd) 
     339                        ii   = idx_bdy(jbdy)%nbi(ib,igrd) 
     340                        ij   = idx_bdy(jbdy)%nbj(ib,igrd) 
     341                        DO ik = 1, jpkm1 
     342                           dta%v2d(ib) = dta%v2d(ib) & 
     343                 &                       + e3v_n(ii,ij,ik) * vmask(ii,ij,ik) * dta%v3d(ib,ik) 
     344                        END DO 
     345                        dta%v2d(ib) =  dta%v2d(ib) * r1_hv_n(ii,ij) 
     346                        DO ik = 1, jpkm1 
     347                           dta%v3d(ib,ik) = dta%v3d(ib,ik) - dta%v2d(ib) 
     348                        END DO 
    288349                     END DO 
    289                   END DO 
    290                ENDIF 
    291  
    292             ENDIF 
    293 #if defined key_si3 
    294             ! convert N-cat fields (input) into jpl-cat (output) 
    295             IF( cn_ice(jbdy) /= 'none' .AND. nn_ice_dta(jbdy) == 1 ) THEN 
    296                jfld_hti = jfld_htit(jbdy) 
    297                jfld_hts = jfld_htst(jbdy) 
    298                jfld_ai  = jfld_ait(jbdy) 
    299                IF    ( jpl /= 1 .AND. nice_cat == 1 ) THEN                       ! case input cat = 1 
    300                   CALL ice_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), & 
    301                      &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    ) 
    302                ELSEIF( jpl /= 1 .AND. nice_cat /= 1 .AND. nice_cat /= jpl ) THEN ! case input cat /=1 and /=jpl 
    303                   CALL ice_var_itd2( bf(jfld_hti)%fnow(:,1,:), bf(jfld_hts)%fnow(:,1,:), bf(jfld_ai)%fnow(:,1,:), & 
    304                      &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    ) 
    305                ENDIF 
    306             ENDIF 
    307 #endif 
     350                  ENDIF 
     351 
     352               ENDIF 
     353#if defined key_si3 
     354               ! convert N-cat fields (input) into jpl-cat (output) 
     355               IF( cn_ice(jbdy) /= 'none' .AND. nn_ice_dta(jbdy) == 1 ) THEN 
     356                  jfld_hti = jfld_htit(jbdy) 
     357                  jfld_hts = jfld_htst(jbdy) 
     358                  jfld_ai  = jfld_ait(jbdy) 
     359                  IF    ( jpl /= 1 .AND. nice_cat == 1 ) THEN                       ! case input cat = 1 
     360                     CALL ice_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), & 
     361                        &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    ) 
     362                  ELSEIF( jpl /= 1 .AND. nice_cat /= 1 .AND. nice_cat /= jpl ) THEN ! case input cat /=1 and /=jpl 
     363                     CALL ice_var_itd2( bf(jfld_hti)%fnow(:,1,:), bf(jfld_hts)%fnow(:,1,:), bf(jfld_ai)%fnow(:,1,:), & 
     364                        &               dta_bdy(jbdy)%h_i     , dta_bdy(jbdy)%h_s     , dta_bdy(jbdy)%a_i    ) 
     365                  ENDIF 
     366               ENDIF 
     367#endif 
     368            ENDIF 
    308369            jstart = jstart + dta%nread(1) 
    309370         ENDIF    ! nn_dta(jbdy) = 1 

  • NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/BDY/bdytides.F90

    r10865 r12059  
    3030 
    3131   PUBLIC   bdytide_init     ! routine called in bdy_init 
     32   PUBLIC   bdytide_update   ! routine called in bdy_dta 
    3233   PUBLIC   bdy_dta_tides    ! routine called in dyn_spg_ts 
    3334 
    … …  
    262263 
    263264 
     265   SUBROUTINE bdytide_update( kt, idx, dta, td, jit, time_offset ) 
     266      !!---------------------------------------------------------------------- 
     267      !!                 ***  SUBROUTINE bdytide_update  *** 
     268      !!                 
     269      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.  
     270      !!                 
     271      !!---------------------------------------------------------------------- 
     272      INTEGER          , INTENT(in   ) ::   kt          ! Main timestep counter 
     273      TYPE(OBC_INDEX)  , INTENT(in   ) ::   idx         ! OBC indices 
     274      TYPE(OBC_DATA)   , INTENT(inout) ::   dta         ! OBC external data 
     275      TYPE(TIDES_DATA) , INTENT(inout) ::   td          ! tidal harmonics data 
     276      INTEGER, OPTIONAL, INTENT(in   ) ::   jit         ! Barotropic timestep counter (for timesplitting option) 
     277      INTEGER, OPTIONAL, INTENT(in   ) ::   time_offset ! time offset in units of timesteps. NB. if jit 
     278      !                                                 ! is present then units = subcycle timesteps. 
     279      !                                                 ! time_offset = 0  => get data at "now"    time level 
     280      !                                                 ! time_offset = -1 => get data at "before" time level 
     281      !                                                 ! time_offset = +1 => get data at "after"  time level 
     282      !                                                 ! etc. 
     283      ! 
     284      INTEGER  ::   itide, igrd, ib       ! dummy loop indices 
     285      INTEGER  ::   time_add              ! time offset in units of timesteps 
     286      INTEGER, DIMENSION(3) ::   ilen0    ! length of boundary data (from OBC arrays) 
     287      REAL(wp) ::   z_arg, z_sarg, zflag, zramp   ! local scalars     
     288      REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost 
     289      !!---------------------------------------------------------------------- 
     290      ! 
     291      ilen0(1) =  SIZE(td%ssh(:,1,1)) 
     292      ilen0(2) =  SIZE(td%u(:,1,1)) 
     293      ilen0(3) =  SIZE(td%v(:,1,1)) 
     294 
     295      zflag=1 
     296      IF ( PRESENT(jit) ) THEN 
     297        IF ( jit /= 1 ) zflag=0 
     298      ENDIF 
     299 
     300      IF ( (nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000) .AND. zflag==1 ) THEN 
     301        ! 
     302        kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt 
     303        ! 
     304        IF(lwp) THEN 
     305           WRITE(numout,*) 
     306           WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt 
     307           WRITE(numout,*) '~~~~~~~~~~~~~~ ' 
     308        ENDIF 
     309        ! 
     310        CALL tide_init_elevation ( idx, td ) 
     311        CALL tide_init_velocities( idx, td ) 
     312        ! 
     313      ENDIF  
     314 
     315      time_add = 0 
     316      IF( PRESENT(time_offset) ) THEN 
     317         time_add = time_offset 
     318      ENDIF 
     319          
     320      IF( PRESENT(jit) ) THEN   
     321         z_arg = ((kt-kt_tide) * rdt + (jit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) ) 
     322      ELSE                               
     323         z_arg = ((kt-kt_tide)+time_add) * rdt 
     324      ENDIF 
     325 
     326      ! Linear ramp on tidal component at open boundaries  
     327      zramp = 1._wp 
     328      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rn_tide_ramp_dt*rday),0._wp),1._wp) 
     329 
     330      DO itide = 1, nb_harmo 
     331         z_sarg = z_arg * tide_harmonics(itide)%omega 
     332         z_cost(itide) = COS( z_sarg ) 
     333         z_sist(itide) = SIN( z_sarg ) 
     334      END DO 
     335 
     336      DO itide = 1, nb_harmo 
     337         igrd=1                              ! SSH on tracer grid 
     338         DO ib = 1, ilen0(igrd) 
     339            dta%ssh(ib) = dta%ssh(ib) + zramp*(td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide)) 
     340         END DO 
     341         igrd=2                              ! U grid 
     342         DO ib = 1, ilen0(igrd) 
     343            dta%u2d(ib) = dta%u2d(ib) + zramp*(td%u  (ib,itide,1)*z_cost(itide) + td%u  (ib,itide,2)*z_sist(itide)) 
     344         END DO 
     345         igrd=3                              ! V grid 
     346         DO ib = 1, ilen0(igrd)  
     347            dta%v2d(ib) = dta%v2d(ib) + zramp*(td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide)) 
     348         END DO 
     349      END DO 
     350      ! 
     351   END SUBROUTINE bdytide_update 
     352 
     353 
    264354   SUBROUTINE bdy_dta_tides( kt, kit, time_offset ) 
    265355      !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.