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 branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 11134

Last change on this file since 11134 was 11134, checked in by jcastill, 5 years ago

Full set of changes as in the original branch

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