Ignore:
Timestamp:
2019-12-09T10:46:17+01:00 (10 months ago)
Author:
smueller
Message:

Merging of the developments in /NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides@12096 with respect to /NEMO/trunk@12072 into /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis (tickets #2175 and #2194)

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/BDY/bdytides.F90

    r11536 r12117  
    1818   USE phycst         ! physical constants 
    1919   USE bdy_oce        ! ocean open boundary conditions 
    20    USE tideini        !  
     20   USE tide_mod       !  
    2121   USE daymod         ! calendar 
    2222   ! 
     
    3030 
    3131   PUBLIC   bdytide_init     ! routine called in bdy_init 
    32    PUBLIC   bdytide_update   ! routine called in bdy_dta 
    3332   PUBLIC   bdy_dta_tides    ! routine called in dyn_spg_ts 
    3433 
     
    4544   TYPE(OBC_DATA)  , PUBLIC, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component) 
    4645 
     46   INTEGER ::   kt_tide 
     47 
    4748   !!---------------------------------------------------------------------- 
    4849   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
     
    6465      CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files 
    6566      LOGICAL                                   ::   ln_bdytide_2ddta    !: If true, read 2d harmonic data 
    66       LOGICAL                                   ::   ln_bdytide_conj     !: If true, assume complex conjugate tidal data 
    6767      !! 
    6868      INTEGER                                   ::   ib_bdy, itide, ib   !: dummy loop indices 
     
    7777      TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut    
    7878      !! 
    79       NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj 
     79      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta 
    8080      !!---------------------------------------------------------------------- 
    8181      ! 
     
    105105            IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries' 
    106106            IF(lwp) WRITE(numout,*) '             read tidal data in 2d files: ', ln_bdytide_2ddta 
    107             IF(lwp) WRITE(numout,*) '             assume complex conjugate   : ', ln_bdytide_conj 
    108107            IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo 
    109108            IF(lwp) THEN  
    110109                    WRITE(numout,*) '             Tidal components: '  
    111110               DO itide = 1, nb_harmo 
    112                   WRITE(numout,*)  '                 ', Wave(ntide(itide))%cname_tide  
     111                  WRITE(numout,*)  '                 ', tide_harmonics(itide)%cname_tide  
    113112               END DO 
    114113            ENDIF  
     
    151150               igrd = 1                       ! Everything is at T-points here 
    152151               DO itide = 1, nb_harmo 
    153                   CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) ) 
    154                   CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) )  
     152                  CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z1', ztr(:,:) ) 
     153                  CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z2', zti(:,:) )  
    155154                  DO ib = 1, ilen0(igrd) 
    156155                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     
    168167               igrd = 2                       ! Everything is at U-points here 
    169168               DO itide = 1, nb_harmo 
    170                   CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) ) 
    171                   CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) ) 
     169                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u1', ztr(:,:) ) 
     170                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u2', zti(:,:) ) 
    172171                  DO ib = 1, ilen0(igrd) 
    173172                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     
    185184               igrd = 3                       ! Everything is at V-points here 
    186185               DO itide = 1, nb_harmo 
    187                   CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) ) 
    188                   CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) ) 
     186                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v1', ztr(:,:) ) 
     187                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v2', zti(:,:) ) 
    189188                  DO ib = 1, ilen0(igrd) 
    190189                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 
     
    210209               DO itide = 1, nb_harmo 
    211210                  !                                                              ! SSH fields 
    212                   clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc' 
     211                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_T.nc' 
    213212                  CALL iom_open( clfile, inum ) 
    214213                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 
     
    218217                  CALL iom_close( inum ) 
    219218                  !                                                              ! U fields 
    220                   clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc' 
     219                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_U.nc' 
    221220                  CALL iom_open( clfile, inum ) 
    222221                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 
     
    226225                  CALL iom_close( inum ) 
    227226                  !                                                              ! V fields 
    228                   clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc' 
     227                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_V.nc' 
    229228                  CALL iom_open( clfile, inum ) 
    230229                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 
     
    240239            ENDIF ! ln_bdytide_2ddta=.true. 
    241240            ! 
    242             IF( ln_bdytide_conj ) THEN    ! assume complex conjugate in data files 
    243                td%ssh0(:,:,2) = - td%ssh0(:,:,2) 
    244                td%u0  (:,:,2) = - td%u0  (:,:,2) 
    245                td%v0  (:,:,2) = - td%v0  (:,:,2) 
    246             ENDIF 
    247             ! 
    248241            ! Allocate slow varying data in the case of time splitting: 
    249242            ! Do it anyway because at this stage knowledge of free surface scheme is unknown 
     
    260253      ! 
    261254   END SUBROUTINE bdytide_init 
    262  
    263  
    264    SUBROUTINE bdytide_update( kt, idx, dta, td, kit, kt_offset ) 
    265       !!---------------------------------------------------------------------- 
    266       !!                 ***  SUBROUTINE bdytide_update  *** 
    267       !!                 
    268       !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.  
    269       !!                 
    270       !!---------------------------------------------------------------------- 
    271       INTEGER          , INTENT(in   ) ::   kt          ! Main timestep counter 
    272       TYPE(OBC_INDEX)  , INTENT(in   ) ::   idx         ! OBC indices 
    273       TYPE(OBC_DATA)   , INTENT(inout) ::   dta         ! OBC external data 
    274       TYPE(TIDES_DATA) , INTENT(inout) ::   td          ! tidal harmonics data 
    275       INTEGER, OPTIONAL, INTENT(in   ) ::   kit         ! Barotropic timestep counter (for timesplitting option) 
    276       INTEGER, OPTIONAL, INTENT(in   ) ::   kt_offset   ! time offset in units of timesteps. NB. if kit 
    277       !                                                 ! is present then units = subcycle timesteps. 
    278       !                                                 ! kt_offset = 0  => get data at "now"    time level 
    279       !                                                 ! kt_offset = -1 => get data at "before" time level 
    280       !                                                 ! kt_offset = +1 => get data at "after"  time level 
    281       !                                                 ! etc. 
    282       ! 
    283       INTEGER  ::   itide, igrd, ib       ! dummy loop indices 
    284       INTEGER  ::   time_add              ! time offset in units of timesteps 
    285       INTEGER, DIMENSION(3) ::   ilen0    ! length of boundary data (from OBC arrays) 
    286       REAL(wp) ::   z_arg, z_sarg, zflag, zramp   ! local scalars     
    287       REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost 
    288       !!---------------------------------------------------------------------- 
    289       ! 
    290       ilen0(1) =  SIZE(td%ssh(:,1,1)) 
    291       ilen0(2) =  SIZE(td%u(:,1,1)) 
    292       ilen0(3) =  SIZE(td%v(:,1,1)) 
    293  
    294       zflag=1 
    295       IF ( PRESENT(kit) ) THEN 
    296         IF ( kit /= 1 ) zflag=0 
    297       ENDIF 
    298  
    299       IF ( (nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000) .AND. zflag==1 ) THEN 
    300         ! 
    301         kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt 
    302         ! 
    303         IF(lwp) THEN 
    304            WRITE(numout,*) 
    305            WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt 
    306            WRITE(numout,*) '~~~~~~~~~~~~~~ ' 
    307         ENDIF 
    308         ! 
    309         CALL tide_init_elevation ( idx, td ) 
    310         CALL tide_init_velocities( idx, td ) 
    311         ! 
    312       ENDIF  
    313  
    314       time_add = 0 
    315       IF( PRESENT(kt_offset) ) THEN 
    316          time_add = kt_offset 
    317       ENDIF 
    318           
    319       IF( PRESENT(kit) ) THEN   
    320          z_arg = ((kt-kt_tide) * rdt + (kit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) ) 
    321       ELSE                               
    322          z_arg = ((kt-kt_tide)+time_add) * rdt 
    323       ENDIF 
    324  
    325       ! Linear ramp on tidal component at open boundaries  
    326       zramp = 1._wp 
    327       IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0._wp),1._wp) 
    328  
    329       DO itide = 1, nb_harmo 
    330          z_sarg = z_arg * omega_tide(itide) 
    331          z_cost(itide) = COS( z_sarg ) 
    332          z_sist(itide) = SIN( z_sarg ) 
    333       END DO 
    334  
    335       DO itide = 1, nb_harmo 
    336          igrd=1                              ! SSH on tracer grid 
    337          DO ib = 1, ilen0(igrd) 
    338             dta%ssh(ib) = dta%ssh(ib) + zramp*(td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide)) 
    339          END DO 
    340          igrd=2                              ! U grid 
    341          DO ib = 1, ilen0(igrd) 
    342             dta%u2d(ib) = dta%u2d(ib) + zramp*(td%u  (ib,itide,1)*z_cost(itide) + td%u  (ib,itide,2)*z_sist(itide)) 
    343          END DO 
    344          igrd=3                              ! V grid 
    345          DO ib = 1, ilen0(igrd)  
    346             dta%v2d(ib) = dta%v2d(ib) + zramp*(td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide)) 
    347          END DO 
    348       END DO 
    349       ! 
    350    END SUBROUTINE bdytide_update 
    351255 
    352256 
     
    392296      ! Linear ramp on tidal component at open boundaries  
    393297      zramp = 1. 
    394       IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rdttideramp*rday),0.),1.) 
     298      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rn_tide_ramp_dt*rday),0.),1.) 
    395299 
    396300      DO ib_bdy = 1,nb_bdy 
     
    433337            DO itide = 1, nb_harmo 
    434338               ! 
    435                z_sarg = (z_arg + zoff) * omega_tide(itide) 
     339               z_sarg = (z_arg + zoff) * tide_harmonics(itide)%omega 
    436340               z_cost = zramp * COS( z_sarg ) 
    437341               z_sist = zramp * SIN( z_sarg ) 
     
    491395         END DO 
    492396         DO ib = 1 , ilen0(igrd) 
    493             mod_tide(ib)=mod_tide(ib)*ftide(itide) 
    494             phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 
     397            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f 
     398            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0+tide_harmonics(itide)%u 
    495399         ENDDO 
    496400         DO ib = 1 , ilen0(igrd) 
     
    530434         END DO 
    531435         DO ib = 1, ilen0(igrd) 
    532             mod_tide(ib)=mod_tide(ib)*ftide(itide) 
    533             phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 
     436            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f 
     437            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u 
    534438         ENDDO 
    535439         DO ib = 1, ilen0(igrd) 
     
    551455         END DO 
    552456         DO ib = 1, ilen0(igrd) 
    553             mod_tide(ib)=mod_tide(ib)*ftide(itide) 
    554             phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 
     457            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f 
     458            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u 
    555459         ENDDO 
    556460         DO ib = 1, ilen0(igrd) 
Note: See TracChangeset for help on using the changeset viewer.