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.
bdytides.F90 in NEMO/trunk/src/OCE/BDY – NEMO

source: NEMO/trunk/src/OCE/BDY/bdytides.F90 @ 10068

Last change on this file since 10068 was 10068, checked in by nicolasmartin, 6 years ago

First part of modifications to have a common default header : fix typos and SVN keywords properties

  • Property svn:keywords set to Id
File size: 27.4 KB
RevLine 
[911]1MODULE bdytides
[1125]2   !!======================================================================
[911]3   !!                       ***  MODULE  bdytides  ***
4   !! Ocean dynamics:   Tidal forcing at open boundaries
[1125]5   !!======================================================================
6   !! History :  2.0  !  2007-01  (D.Storkey)  Original code
7   !!            2.3  !  2008-01  (J.Holt)  Add date correction. Origins POLCOMS v6.3 2007
8   !!            3.0  !  2008-04  (NEMO team)  add in the reference version
[2528]9   !!            3.3  !  2010-09  (D.Storkey and E.O'Dea)  bug fixes
[3651]10   !!            3.4  !  2012-09  (G. Reffray and J. Chanut) New inputs + mods
[4292]11   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
[1125]12   !!----------------------------------------------------------------------
[6140]13   !!   bdytide_init  : read of namelist and initialisation of tidal harmonics data
14   !!   tide_update   : calculation of tidal forcing at each timestep
[1125]15   !!----------------------------------------------------------------------
[6140]16   USE oce            ! ocean dynamics and tracers
17   USE dom_oce        ! ocean space and time domain
18   USE phycst         ! physical constants
19   USE bdy_oce        ! ocean open boundary conditions
20   USE tideini        !
21   USE daymod         ! calendar
22   !
23   USE in_out_manager ! I/O units
24   USE iom            ! xIO server
25   USE fldread        !
26   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
[1125]27
[911]28   IMPLICIT NONE
29   PRIVATE
30
[3651]31   PUBLIC   bdytide_init     ! routine called in bdy_init
32   PUBLIC   bdytide_update   ! routine called in bdy_dta
[4292]33   PUBLIC   bdy_dta_tides    ! routine called in dyn_spg_ts
[911]34
[3294]35   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data
[6140]36      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh0     !: Tidal constituents : SSH0   (read in file)
37      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u0, v0   !: Tidal constituents : U0, V0 (read in file)
38      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh      !: Tidal constituents : SSH    (after nodal cor.)
39      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u , v    !: Tidal constituents : U , V  (after nodal cor.)
[3294]40   END TYPE TIDES_DATA
[911]41
[4354]42!$AGRIF_DO_NOT_TREAT
[3651]43   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data
[4354]44!$AGRIF_END_DO_NOT_TREAT
[5930]45   TYPE(OBC_DATA)  , PUBLIC, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component)
[911]46
[1125]47   !!----------------------------------------------------------------------
[9598]48   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[6140]49   !! $Id$
[10068]50   !! Software governed by the CeCILL license (see ./LICENSE)
[1125]51   !!----------------------------------------------------------------------
[911]52CONTAINS
53
[3651]54   SUBROUTINE bdytide_init
[1125]55      !!----------------------------------------------------------------------
[3651]56      !!                    ***  SUBROUTINE bdytide_init  ***
[1125]57      !!                     
[3294]58      !! ** Purpose : - Read in namelist for tides and initialise external
59      !!                tidal harmonics data
[911]60      !!
61      !!----------------------------------------------------------------------
[3294]62      !! namelist variables
63      !!-------------------
[3651]64      CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files
65      LOGICAL                                   ::   ln_bdytide_2ddta    !: If true, read 2d harmonic data
66      LOGICAL                                   ::   ln_bdytide_conj     !: If true, assume complex conjugate tidal data
[1125]67      !!
[3651]68      INTEGER                                   ::   ib_bdy, itide, ib   !: dummy loop indices
69      INTEGER                                   ::   ii, ij              !: dummy loop indices
[3294]70      INTEGER                                   ::   inum, igrd
[3651]71      INTEGER, DIMENSION(3)                     ::   ilen0       !: length of boundary data (from OBC arrays)
72      INTEGER, POINTER, DIMENSION(:)            ::   nblen, nblenrim     ! short cuts
[4147]73      INTEGER                                   ::   ios                 ! Local integer output status for namelist read
[3651]74      CHARACTER(len=80)                         ::   clfile              !: full file name for tidal input file
75      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            !: work space to read in tidal harmonics data
[9125]76      REAL(wp),ALLOCATABLE, DIMENSION(:,:)      ::   ztr, zti            !:  "     "    "   "   "   "        "      "
[3294]77      !!
[3651]78      TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut   
[5132]79      TYPE(MAP_POINTER), DIMENSION(jpbgrd)      ::   ibmap_ptr           !: array of pointers to nbmap
[3294]80      !!
[3651]81      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj
[1125]82      !!----------------------------------------------------------------------
[6140]83      !
[3651]84      IF (nb_bdy>0) THEN
85         IF(lwp) WRITE(numout,*)
86         IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries'
87         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
88      ENDIF
[911]89
[4147]90      REWIND(numnam_cfg)
[3651]91
[3294]92      DO ib_bdy = 1, nb_bdy
[7646]93         IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
94            !
[3294]95            td => tides(ib_bdy)
[3651]96            nblen => idx_bdy(ib_bdy)%nblen
97            nblenrim => idx_bdy(ib_bdy)%nblenrim
[911]98
[3294]99            ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
100            filtide(:) = ''
[911]101
[4147]102            ! Don't REWIND here - may need to read more than one of these namelists.
103            READ  ( numnam_ref, nambdy_tide, IOSTAT = ios, ERR = 901)
[9168]104901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_tide in reference namelist', lwp )
[4147]105            READ  ( numnam_cfg, nambdy_tide, IOSTAT = ios, ERR = 902 )
[9168]106902         IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy_tide in configuration namelist', lwp )
[4624]107            IF(lwm) WRITE ( numond, nambdy_tide )
[3651]108            !                                               ! Parameter control and print
109            IF(lwp) WRITE(numout,*) '  '
110            IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
111            IF(lwp) WRITE(numout,*) '             read tidal data in 2d files: ', ln_bdytide_2ddta
112            IF(lwp) WRITE(numout,*) '             assume complex conjugate   : ', ln_bdytide_conj
113            IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo
114            IF(lwp) THEN
[5084]115                    WRITE(numout,*) '             Tidal components: ' 
[3651]116               DO itide = 1, nb_harmo
[5084]117                  WRITE(numout,*)  '                 ', Wave(ntide(itide))%cname_tide 
[3651]118               END DO
119            ENDIF
120            IF(lwp) WRITE(numout,*) ' '
[1125]121
[3651]122            ! Allocate space for tidal harmonics data - get size from OBC data arrays
123            ! -----------------------------------------------------------------------
[911]124
[3651]125            ! JC: If FRS scheme is used, we assume that tidal is needed over the whole
126            ! relaxation area     
[7646]127            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = nblen   (:)
128            ELSE                                   ;   ilen0(:) = nblenrim(:)
[3294]129            ENDIF
[911]130
[3651]131            ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) )
132            ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) )
[911]133
[3651]134            ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) )
135            ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) )
[911]136
[3651]137            ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) )
138            ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) )
[911]139
[4292]140            td%ssh0(:,:,:) = 0._wp
141            td%ssh (:,:,:) = 0._wp
142            td%u0  (:,:,:) = 0._wp
143            td%u   (:,:,:) = 0._wp
144            td%v0  (:,:,:) = 0._wp
145            td%v   (:,:,:) = 0._wp
[3294]146
[7646]147            IF( ln_bdytide_2ddta ) THEN
[3651]148               ! It is assumed that each data file contains all complex harmonic amplitudes
[7646]149               ! given on the global domain (ie global, jpiglo x jpjglo)
[3651]150               !
[9125]151               ALLOCATE( zti(jpi,jpj), ztr(jpi,jpj) )
[3651]152               !
153               ! SSH fields
154               clfile = TRIM(filtide)//'_grid_T.nc'
[7646]155               CALL iom_open( clfile , inum ) 
[3651]156               igrd = 1                       ! Everything is at T-points here
157               DO itide = 1, nb_harmo
[7646]158                  CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) )
159                  CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) ) 
[3651]160                  DO ib = 1, ilen0(igrd)
161                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
162                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
163                     td%ssh0(ib,itide,1) = ztr(ii,ij)
164                     td%ssh0(ib,itide,2) = zti(ii,ij)
165                  END DO
166               END DO
[3294]167               CALL iom_close( inum )
[3651]168               !
169               ! U fields
170               clfile = TRIM(filtide)//'_grid_U.nc'
[7646]171               CALL iom_open( clfile , inum ) 
[3651]172               igrd = 2                       ! Everything is at U-points here
173               DO itide = 1, nb_harmo
[7646]174                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) )
175                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) )
[3651]176                  DO ib = 1, ilen0(igrd)
177                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
178                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
179                     td%u0(ib,itide,1) = ztr(ii,ij)
180                     td%u0(ib,itide,2) = zti(ii,ij)
181                  END DO
182               END DO
[3294]183               CALL iom_close( inum )
[3651]184               !
185               ! V fields
186               clfile = TRIM(filtide)//'_grid_V.nc'
[7646]187               CALL iom_open( clfile , inum ) 
[3651]188               igrd = 3                       ! Everything is at V-points here
189               DO itide = 1, nb_harmo
[7646]190                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) )
191                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) )
[3651]192                  DO ib = 1, ilen0(igrd)
193                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
194                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
195                     td%v0(ib,itide,1) = ztr(ii,ij)
196                     td%v0(ib,itide,2) = zti(ii,ij)
197                  END DO
198               END DO 
[3294]199               CALL iom_close( inum )
200               !
[9125]201               DEALLOCATE( ztr, zti ) 
[3651]202               !
203            ELSE           
204               !
205               ! Read tidal data only on bdy segments
206               !
207               ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) )
[5132]208               !
209               ! Set map structure
[7646]210               ibmap_ptr(1)%ptr => idx_bdy(ib_bdy)%nbmap(:,1)   ;   ibmap_ptr(1)%ll_unstruc = ln_coords_file(ib_bdy)
211               ibmap_ptr(2)%ptr => idx_bdy(ib_bdy)%nbmap(:,2)   ;   ibmap_ptr(2)%ll_unstruc = ln_coords_file(ib_bdy)
212               ibmap_ptr(3)%ptr => idx_bdy(ib_bdy)%nbmap(:,3)   ;   ibmap_ptr(3)%ll_unstruc = ln_coords_file(ib_bdy)
[3294]213
[3651]214               ! Open files and read in tidal forcing data
215               ! -----------------------------------------
[3294]216
[3651]217               DO itide = 1, nb_harmo
218                  !                                                              ! SSH fields
219                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc'
220                  CALL iom_open( clfile, inum )
[5132]221                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1,  ibmap_ptr(1) )
[3651]222                  td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1)
[5132]223                  CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1,  ibmap_ptr(1) )
[3651]224                  td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1)
225                  CALL iom_close( inum )
226                  !                                                              ! U fields
227                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc'
228                  CALL iom_open( clfile, inum )
[5132]229                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, ibmap_ptr(2) )
[3651]230                  td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1)
[5132]231                  CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, ibmap_ptr(2) )
[3651]232                  td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1)
233                  CALL iom_close( inum )
234                  !                                                              ! V fields
235                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc'
236                  CALL iom_open( clfile, inum )
[5132]237                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, ibmap_ptr(3) )
[3651]238                  td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1)
[5132]239                  CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, ibmap_ptr(3) )
[3651]240                  td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1)
241                  CALL iom_close( inum )
[3294]242                  !
[3651]243               END DO ! end loop on tidal components
[3294]244               !
[3651]245               DEALLOCATE( dta_read )
[7646]246               !
[3651]247            ENDIF ! ln_bdytide_2ddta=.true.
[1125]248            !
[6140]249            IF( ln_bdytide_conj ) THEN    ! assume complex conjugate in data files
[3651]250               td%ssh0(:,:,2) = - td%ssh0(:,:,2)
251               td%u0  (:,:,2) = - td%u0  (:,:,2)
252               td%v0  (:,:,2) = - td%v0  (:,:,2)
253            ENDIF
254            !
[5930]255            ! Allocate slow varying data in the case of time splitting:
256            ! Do it anyway because at this stage knowledge of free surface scheme is unknown
257            ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) )
258            ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) )
259            ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) )
[6140]260            dta_bdy_s(ib_bdy)%ssh(:) = 0._wp
261            dta_bdy_s(ib_bdy)%u2d(:) = 0._wp
262            dta_bdy_s(ib_bdy)%v2d(:) = 0._wp
[4292]263            !
[7646]264         ENDIF ! nn_dyn2d_dta(ib_bdy) >= 2
[1125]265         !
[3294]266      END DO ! loop on ib_bdy
[6140]267      !
268   END SUBROUTINE bdytide_init
[911]269
[1125]270
[6140]271   SUBROUTINE bdytide_update( kt, idx, dta, td, jit, time_offset )
[1125]272      !!----------------------------------------------------------------------
[3651]273      !!                 ***  SUBROUTINE bdytide_update  ***
[1125]274      !!               
[3294]275      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
[911]276      !!               
[1125]277      !!----------------------------------------------------------------------
[6140]278      INTEGER          , INTENT(in   ) ::   kt          ! Main timestep counter
279      TYPE(OBC_INDEX)  , INTENT(in   ) ::   idx         ! OBC indices
280      TYPE(OBC_DATA)   , INTENT(inout) ::   dta         ! OBC external data
281      TYPE(TIDES_DATA) , INTENT(inout) ::   td          ! tidal harmonics data
282      INTEGER, OPTIONAL, INTENT(in   ) ::   jit         ! Barotropic timestep counter (for timesplitting option)
283      INTEGER, OPTIONAL, INTENT(in   ) ::   time_offset ! time offset in units of timesteps. NB. if jit
284      !                                                 ! is present then units = subcycle timesteps.
285      !                                                 ! time_offset = 0  => get data at "now"    time level
286      !                                                 ! time_offset = -1 => get data at "before" time level
287      !                                                 ! time_offset = +1 => get data at "after"  time level
288      !                                                 ! etc.
289      !
290      INTEGER  ::   itide, igrd, ib       ! dummy loop indices
291      INTEGER  ::   time_add              ! time offset in units of timesteps
292      INTEGER, DIMENSION(3) ::   ilen0    ! length of boundary data (from OBC arrays)
293      REAL(wp) ::   z_arg, z_sarg, zflag, zramp   ! local scalars   
[3651]294      REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost
[1125]295      !!----------------------------------------------------------------------
[6140]296      !
[3651]297      ilen0(1) =  SIZE(td%ssh(:,1,1))
298      ilen0(2) =  SIZE(td%u(:,1,1))
299      ilen0(3) =  SIZE(td%v(:,1,1))
300
301      zflag=1
302      IF ( PRESENT(jit) ) THEN
303        IF ( jit /= 1 ) zflag=0
304      ENDIF
305
[6140]306      IF ( (nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000) .AND. zflag==1 ) THEN
[3651]307        !
[6140]308        kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
[3651]309        !
310        IF(lwp) THEN
311           WRITE(numout,*)
312           WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt
313           WRITE(numout,*) '~~~~~~~~~~~~~~ '
314        ENDIF
315        !
316        CALL tide_init_elevation ( idx, td )
317        CALL tide_init_velocities( idx, td )
318        !
319      ENDIF
320
[3294]321      time_add = 0
322      IF( PRESENT(time_offset) ) THEN
323         time_add = time_offset
324      ENDIF
325         
326      IF( PRESENT(jit) ) THEN 
[4292]327         z_arg = ((kt-kt_tide) * rdt + (jit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) )
[3294]328      ELSE                             
[3651]329         z_arg = ((kt-kt_tide)+time_add) * rdt
[1125]330      ENDIF
[911]331
[3651]332      ! Linear ramp on tidal component at open boundaries
[4292]333      zramp = 1._wp
334      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0._wp),1._wp)
[3651]335
336      DO itide = 1, nb_harmo
337         z_sarg = z_arg * omega_tide(itide)
[1125]338         z_cost(itide) = COS( z_sarg )
339         z_sist(itide) = SIN( z_sarg )
340      END DO
[911]341
[3651]342      DO itide = 1, nb_harmo
343         igrd=1                              ! SSH on tracer grid
344         DO ib = 1, ilen0(igrd)
345            dta%ssh(ib) = dta%ssh(ib) + zramp*(td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide))
[1125]346         END DO
[3294]347         igrd=2                              ! U grid
[3651]348         DO ib = 1, ilen0(igrd)
349            dta%u2d(ib) = dta%u2d(ib) + zramp*(td%u  (ib,itide,1)*z_cost(itide) + td%u  (ib,itide,2)*z_sist(itide))
[1125]350         END DO
[3294]351         igrd=3                              ! V grid
[3651]352         DO ib = 1, ilen0(igrd) 
353            dta%v2d(ib) = dta%v2d(ib) + zramp*(td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide))
[1125]354         END DO
355      END DO
356      !
[3651]357   END SUBROUTINE bdytide_update
[911]358
[7646]359
[4292]360   SUBROUTINE bdy_dta_tides( kt, kit, time_offset )
361      !!----------------------------------------------------------------------
362      !!                 ***  SUBROUTINE bdy_dta_tides  ***
363      !!               
364      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
365      !!               
366      !!----------------------------------------------------------------------
[6140]367      INTEGER,           INTENT(in) ::   kt          ! Main timestep counter
368      INTEGER, OPTIONAL, INTENT(in) ::   kit         ! Barotropic timestep counter (for timesplitting option)
369      INTEGER, OPTIONAL, INTENT(in) ::   time_offset ! time offset in units of timesteps. NB. if kit
370      !                                              ! is present then units = subcycle timesteps.
371      !                                              ! time_offset = 0  => get data at "now"    time level
372      !                                              ! time_offset = -1 => get data at "before" time level
373      !                                              ! time_offset = +1 => get data at "after"  time level
374      !                                              ! etc.
375      !
376      LOGICAL  ::   lk_first_btstp            ! =.TRUE. if time splitting and first barotropic step
377      INTEGER  ::   itide, ib_bdy, ib, igrd   ! loop indices
378      INTEGER  ::   time_add                  ! time offset in units of timesteps
379      INTEGER, DIMENSION(jpbgrd)   ::   ilen0 
380      INTEGER, DIMENSION(1:jpbgrd) ::   nblen, nblenrim  ! short cuts
381      REAL(wp) ::   z_arg, z_sarg, zramp, zoff, z_cost, z_sist     
[4292]382      !!----------------------------------------------------------------------
[6140]383      !
[4292]384      lk_first_btstp=.TRUE.
385      IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF
386
387      time_add = 0
388      IF( PRESENT(time_offset) ) THEN
389         time_add = time_offset
390      ENDIF
391     
392      ! Absolute time from model initialization:   
393      IF( PRESENT(kit) ) THEN 
[5913]394         z_arg = ( kt + (kit+time_add-1) / REAL(nn_baro,wp) ) * rdt
[4292]395      ELSE                             
396         z_arg = ( kt + time_add ) * rdt
397      ENDIF
398
399      ! Linear ramp on tidal component at open boundaries
400      zramp = 1.
401      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rdttideramp*rday),0.),1.)
402
403      DO ib_bdy = 1,nb_bdy
[7646]404         !
405         IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
406            !
[4292]407            nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd)
408            nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd)
[7646]409            !
410            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = nblen   (:)
411            ELSE                                   ;   ilen0(:) = nblenrim(:)
[4292]412            ENDIF     
[7646]413            !
[4292]414            ! We refresh nodal factors every day below
415            ! This should be done somewhere else
[6140]416            IF ( ( nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000 ) .AND. lk_first_btstp ) THEN
[4292]417               !
[6140]418               kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
[4292]419               !
420               IF(lwp) THEN
421               WRITE(numout,*)
422               WRITE(numout,*) 'bdy_tide_dta : Refresh nodal factors for tidal open bdy data at kt=',kt
423               WRITE(numout,*) '~~~~~~~~~~~~~~ '
424               ENDIF
425               !
426               CALL tide_init_elevation ( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
427               CALL tide_init_velocities( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
428               !
429            ENDIF
430            zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time
431            !
[5930]432            ! If time splitting, initialize arrays from slow varying open boundary data:
433            IF ( PRESENT(kit) ) THEN           
434               IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1))
435               IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2))
436               IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3))
[4292]437            ENDIF
438            !
439            ! Update open boundary data arrays:
440            DO itide = 1, nb_harmo
441               !
442               z_sarg = (z_arg + zoff) * omega_tide(itide)
443               z_cost = zramp * COS( z_sarg )
444               z_sist = zramp * SIN( z_sarg )
445               !
[4758]446               IF ( dta_bdy(ib_bdy)%ll_ssh ) THEN
447                  igrd=1                              ! SSH on tracer grid
448                  DO ib = 1, ilen0(igrd)
449                     dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + &
450                        &                      ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + &
451                        &                        tides(ib_bdy)%ssh(ib,itide,2)*z_sist )
452                  END DO
453               ENDIF
[4292]454               !
[4758]455               IF ( dta_bdy(ib_bdy)%ll_u2d ) THEN
456                  igrd=2                              ! U grid
457                  DO ib = 1, ilen0(igrd)
458                     dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + &
459                        &                      ( tides(ib_bdy)%u(ib,itide,1)*z_cost + &
460                        &                        tides(ib_bdy)%u(ib,itide,2)*z_sist )
461                  END DO
462               ENDIF
[4292]463               !
[4758]464               IF ( dta_bdy(ib_bdy)%ll_v2d ) THEN
465                  igrd=3                              ! V grid
466                  DO ib = 1, ilen0(igrd) 
467                     dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + &
468                        &                      ( tides(ib_bdy)%v(ib,itide,1)*z_cost + &
469                        &                        tides(ib_bdy)%v(ib,itide,2)*z_sist )
470                  END DO
471               ENDIF
472            END DO             
[4292]473         END IF
474      END DO
475      !
476   END SUBROUTINE bdy_dta_tides
477
[6140]478
[3651]479   SUBROUTINE tide_init_elevation( idx, td )
[1125]480      !!----------------------------------------------------------------------
[3651]481      !!                 ***  ROUTINE tide_init_elevation  ***
[1125]482      !!----------------------------------------------------------------------
[6140]483      TYPE(OBC_INDEX) , INTENT(in   ) ::   idx   ! OBC indices
484      TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
485      !
486      INTEGER ::   itide, igrd, ib       ! dummy loop indices
487      INTEGER, DIMENSION(1) ::   ilen0   ! length of boundary data (from OBC arrays)
[3651]488      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
[6140]489      !!----------------------------------------------------------------------
490      !
[3651]491      igrd=1   
492                              ! SSH on tracer grid.
493      ilen0(1) =  SIZE(td%ssh0(:,1,1))
[6140]494      !
495      ALLOCATE( mod_tide(ilen0(igrd)), phi_tide(ilen0(igrd)) )
496      !
[3651]497      DO itide = 1, nb_harmo
498         DO ib = 1, ilen0(igrd)
499            mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.)
500            phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
[1125]501         END DO
[3651]502         DO ib = 1 , ilen0(igrd)
503            mod_tide(ib)=mod_tide(ib)*ftide(itide)
504            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
505         ENDDO
506         DO ib = 1 , ilen0(igrd)
507            td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
508            td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
509         ENDDO
510      END DO
[6140]511      !
512      DEALLOCATE( mod_tide, phi_tide )
513      !
514   END SUBROUTINE tide_init_elevation
[911]515
516
[3651]517   SUBROUTINE tide_init_velocities( idx, td )
[1125]518      !!----------------------------------------------------------------------
[3651]519      !!                 ***  ROUTINE tide_init_elevation  ***
[1125]520      !!----------------------------------------------------------------------
[6140]521      TYPE(OBC_INDEX) , INTENT(in   ) ::   idx   ! OBC indices
522      TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
523      !
524      INTEGER ::   itide, igrd, ib       ! dummy loop indices
525      INTEGER, DIMENSION(3) ::   ilen0   ! length of boundary data (from OBC arrays)
[3651]526      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
[6140]527      !!----------------------------------------------------------------------
528      !
[3651]529      ilen0(2) =  SIZE(td%u0(:,1,1))
530      ilen0(3) =  SIZE(td%v0(:,1,1))
[6140]531      !
[3651]532      igrd=2                                 ! U grid.
[6140]533      !
534      ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
535      !
[3651]536      DO itide = 1, nb_harmo
537         DO ib = 1, ilen0(igrd)
538            mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.)
539            phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
540         END DO
541         DO ib = 1, ilen0(igrd)
542            mod_tide(ib)=mod_tide(ib)*ftide(itide)
543            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
544         ENDDO
545         DO ib = 1, ilen0(igrd)
546            td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
547            td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
548         ENDDO
549      END DO
[6140]550      !
551      DEALLOCATE( mod_tide , phi_tide )
552      !
[3651]553      igrd=3                                 ! V grid.
[6140]554      !
555      ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
[911]556
[3651]557      DO itide = 1, nb_harmo
558         DO ib = 1, ilen0(igrd)
559            mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.)
560            phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
561         END DO
562         DO ib = 1, ilen0(igrd)
563            mod_tide(ib)=mod_tide(ib)*ftide(itide)
564            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
565         ENDDO
566         DO ib = 1, ilen0(igrd)
567            td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
568            td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
569         ENDDO
570      END DO
[6140]571      !
572      DEALLOCATE( mod_tide, phi_tide )
573      !
574  END SUBROUTINE tide_init_velocities
[1125]575
576   !!======================================================================
[911]577END MODULE bdytides
[4292]578
Note: See TracBrowser for help on using the repository browser.