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 trunk/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 7646

Last change on this file since 7646 was 7646, checked in by timgraham, 7 years ago

Merge of dev_merge_2016 into trunk. UPDATE TO ARCHFILES NEEDED for XIOS2.
LIM_SRC_s/limrhg.F90 to follow in next commit due to change of kind (I'm unable to do it in this commit).
Merged using the following steps:

1) svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk .
2) Resolve minor conflicts in sette.sh and namelist_cfg for ORCA2LIM3 (due to a change in trunk after branch was created)
3) svn commit
4) svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
5) svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2016/dev_merge_2016 .
6) At this stage I checked out a clean copy of the branch to compare against what is about to be committed to the trunk.
6) svn commit #Commit code to the trunk

In this commit I have also reverted a change to Fcheck_archfile.sh which was causing problems on the Paris machine.

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