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

source: branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 3367

Last change on this file since 3367 was 3367, checked in by greffray, 12 years ago

Merge tidal packages
Merge OBC-BDY packages
Add atmospheric pressure at the open boundaries
Modification of the namelist for AMM12 configuration
  • Property svn:keywords set to Id
File size: 16.4 KB
Line 
1MODULE bdytides
2   !!======================================================================
3   !!                       ***  MODULE  bdytides  ***
4   !! Ocean dynamics:   Tidal forcing at open boundaries
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
9   !!            3.3  !  2010-09  (D.Storkey and E.O'Dea)  bug fixes
10   !!----------------------------------------------------------------------
11#if defined key_bdy
12   !!----------------------------------------------------------------------
13   !!   'key_bdy'     Open Boundary Condition
14   !!----------------------------------------------------------------------
15   !!   PUBLIC
16   !!      bdytide_init     : read of namelist and initialisation of tidal harmonics data
17   !!      tide_update   : calculation of tidal forcing at each timestep
18   !!----------------------------------------------------------------------
19   USE timing          ! Timing
20   USE oce             ! ocean dynamics and tracers
21   USE dom_oce         ! ocean space and time domain
22   USE iom
23   USE in_out_manager  ! I/O units
24   USE phycst          ! physical constants
25   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
26   USE bdy_par         ! Unstructured boundary parameters
27   USE bdy_oce         ! ocean open boundary conditions
28   USE daymod          ! calendar
29   USE tideini
30   USE tide_mod
31   USE fldread, ONLY: fld_map
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   bdytide_init  ! routine called in nemo_init
37   PUBLIC   tide_update   ! routine called in bdydyn
38
39   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data
40      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh0       !: Tidal constituents : SSH0 (read in file)
41      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u0         !: Tidal constituents : U0   (read in file)
42      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v0         !: Tidal constituents : V0   (read in file)
43      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh        !: Tidal constituents : SSH  (after nodal cor.)
44      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u          !: Tidal constituents : U    (after nodal cor.)
45      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v          !: Tidal constituents : V    (after nodal cor.)
46   END TYPE TIDES_DATA
47
48   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data
49   
50   !!----------------------------------------------------------------------
51   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
52   !! $Id$
53   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE bdytide_init
58      !!----------------------------------------------------------------------
59      !!                    ***  SUBROUTINE bdytide_init  ***
60      !!                     
61      !! ** Purpose : - Read in namelist for tides and initialise external
62      !!                tidal harmonics data
63      !!
64      !!----------------------------------------------------------------------
65      !! namelist variables
66      !!-------------------
67      CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files
68      LOGICAL                                   ::   ln_conjug = .FALSE. !: F/T : tidal data in complex/complex conjugate
69      !!
70      INTEGER                                   ::   ib_bdy, itide, ib   !: dummy loop indices
71      INTEGER                                   ::   inum, igrd
72      INTEGER, DIMENSION(3)                     ::   ilen0               !: length of boundary data (from OBC arrays)
73      CHARACTER(len=80)                         ::   clfile              !: full file name for tidal input file
74      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            !: work space to read in tidal harmonics data
75      !!
76      TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut   
77      !!
78      NAMELIST/nambdy_tide/filtide, ln_conjug
79      !!----------------------------------------------------------------------
80
81      IF( nn_timing == 1 ) CALL timing_start('bdytide_init')
82
83      IF(lwp) WRITE(numout,*)
84      IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries'
85      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
86
87      REWIND(numnam)
88      DO ib_bdy = 1, nb_bdy
89         IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN
90
91            td => tides(ib_bdy)
92
93            ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
94            filtide(:) = ''
95
96            ! Don't REWIND here - may need to read more than one of these namelists.
97            READ  ( numnam, nambdy_tide )
98            !                                               ! Parameter control and print
99            IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
100            IF(lwp) WRITE(numout,*) '             tidal components specified ', nb_harmo
101            IF(lwp) WRITE(numout,*) '                ', Wave(ntide(1:nb_harmo))%cname_tide
102            IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : '
103            IF(lwp) WRITE(numout,*) '                ', omega_tide(1:nb_harmo)
104
105            ! Allocate space for tidal harmonics data - get size from OBC data arrays
106            ! -----------------------------------------------------------------------
107
108            ilen0(1) = SIZE( dta_bdy(ib_bdy)%ssh ) 
109            ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) )
110            ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) )
111
112            ilen0(2) = SIZE( dta_bdy(ib_bdy)%u2d ) 
113            ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) )
114            ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) )
115
116            ilen0(3) = SIZE( dta_bdy(ib_bdy)%v2d ) 
117            ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) )
118            ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) )
119
120            ALLOCATE( dta_read( MAXVAL(ilen0), 1, 1 ) )
121
122            ! Open files and read in tidal forcing data
123            ! -----------------------------------------
124
125            DO itide = 1, nb_harmo
126               !                                                              ! SSH fields
127               clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc'
128               CALL iom_open( clfile, inum )
129               CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
130               td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1)
131               CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
132               td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1)
133               CALL iom_close( inum )
134               !                                                              ! U fields
135               clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc'
136               CALL iom_open( clfile, inum )
137               CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
138               td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1)
139               CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
140               td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1)
141               CALL iom_close( inum )
142               !                                                              ! V fields
143               clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc'
144               CALL iom_open( clfile, inum )
145               CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
146               td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1)
147               CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
148               td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1)
149               CALL iom_close( inum )
150               IF ( ln_conjug ) THEN
151                  IF(lwp) WRITE(numout,*) '       The tidal input data are written in complex conjugate'
152                  td%ssh0(:,itide,2) = - td%ssh0(:,itide,2)
153                  td%u0  (:,itide,2) = - td%u0  (:,itide,2)
154                  td%v0  (:,itide,2) = - td%v0  (:,itide,2)
155               ENDIF
156               !
157            END DO ! end loop on tidal components
158            !
159            DEALLOCATE( dta_read )
160            !
161         ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2
162         !
163      END DO ! loop on ib_bdy
164
165      IF( nn_timing == 1 ) CALL timing_stop('bdytide_init')
166
167   END SUBROUTINE bdytide_init
168
169   SUBROUTINE tide_update ( kt, idx, dta, td, jit, time_offset )
170      !!----------------------------------------------------------------------
171      !!                 ***  SUBROUTINE tide_update  ***
172      !!               
173      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
174      !!               
175      !!----------------------------------------------------------------------
176      INTEGER, INTENT( in )            ::   kt          ! Main timestep counter
177      TYPE(OBC_INDEX), INTENT( in )    ::   idx         ! OBC indices
178      TYPE(OBC_DATA),  INTENT(inout)   ::   dta         ! OBC external data
179      TYPE(TIDES_DATA),INTENT( inout ) ::   td          ! tidal harmonics data
180      INTEGER,INTENT(in),OPTIONAL      ::   jit         ! Barotropic timestep counter (for timesplitting option)
181      INTEGER,INTENT( in ), OPTIONAL   ::   time_offset ! time offset in units of timesteps. NB. if jit
182                                                        ! is present then units = subcycle timesteps.
183                                                        ! time_offset = 0  => get data at "now"    time level
184                                                        ! time_offset = -1 => get data at "before" time level
185                                                        ! time_offset = +1 => get data at "after"  time level
186                                                        ! etc.
187      !!
188      INTEGER                          :: itide, igrd, ib   ! dummy loop indices
189      INTEGER                          :: time_add          ! time offset in units of timesteps
190      REAL(wp)                         :: z_arg, z_sarg, zflag     
191      REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost
192      !!----------------------------------------------------------------------
193
194      IF( nn_timing == 1 ) CALL timing_start('tide_update')
195
196      zflag=1
197      IF ( PRESENT(jit) ) THEN
198        IF ( jit /= 1 ) zflag=0
199      ENDIF
200
201      IF ( nsec_day == NINT(0.5 * rdttra(1)) .AND. zflag==1 ) THEN
202        !
203        kt_tide = kt
204        !
205        IF(lwp) THEN
206           WRITE(numout,*)
207           WRITE(numout,*) 'bdy_tide : (re)Initialization of the tidal bdy forcing at kt=',kt
208           WRITE(numout,*) '~~~~~~~ '
209        ENDIF
210        !
211        CALL tide_init_elevation ( idx, td )
212        CALL tide_init_velocities( idx, td )
213        !
214      ENDIF
215
216      time_add = 0
217      IF( PRESENT(time_offset) ) THEN
218         time_add = time_offset
219      ENDIF
220         
221      IF( PRESENT(jit) ) THEN 
222         z_arg = ( ((kt-kt_tide)-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) )
223      ELSE                             
224         z_arg = ((kt-kt_tide)+time_add) * rdt
225      ENDIF
226
227      DO itide = 1, nb_harmo
228         z_sarg = z_arg * omega_tide(itide)
229         z_cost(itide) = COS( z_sarg )
230         z_sist(itide) = SIN( z_sarg )
231      END DO
232
233      DO itide = 1, nb_harmo
234         igrd=1                              ! SSH on tracer grid
235         DO ib = 1, idx%nblenrim(igrd)
236            dta%ssh(ib) = dta%ssh(ib) + td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide)
237         END DO
238         igrd=2                              ! U grid
239         DO ib=1, idx%nblenrim(igrd)
240            dta%u2d(ib) = dta%u2d(ib) + td%u  (ib,itide,1)*z_cost(itide) + td%u  (ib,itide,2)*z_sist(itide)
241         END DO
242         igrd=3                              ! V grid
243         DO ib=1, idx%nblenrim(igrd)
244            dta%v2d(ib) = dta%v2d(ib) + td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide)
245         END DO
246      END DO
247      !
248      IF( nn_timing == 1 ) CALL timing_stop('tide_update')
249      !
250   END SUBROUTINE tide_update
251
252   SUBROUTINE tide_init_elevation( idx, td )
253      !!----------------------------------------------------------------------
254      !!                 ***  ROUTINE tide_init_elevation  ***
255      !!----------------------------------------------------------------------
256      TYPE(OBC_INDEX), INTENT( in )      ::   idx     ! OBC indices
257      TYPE(TIDES_DATA),INTENT( inout )   ::   td      ! tidal harmonics data
258      !! * Local declarations
259      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
260      INTEGER                            ::   itide, igrd, ib      ! dummy loop indices
261
262      igrd=1                                 ! SSH on tracer grid.
263
264      ALLOCATE(mod_tide(idx%nblenrim(igrd)),phi_tide(idx%nblenrim(igrd)))
265
266      DO itide = 1, nb_harmo
267         DO ib = 1, idx%nblenrim(igrd)
268            mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.)
269            phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
270         END DO
271         DO ib = 1, idx%nblenrim(igrd)
272            mod_tide(ib)=mod_tide(ib)*ftide(itide)
273            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
274         ENDDO
275         DO ib = 1, idx%nblenrim(igrd)
276            td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
277            td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
278         ENDDO
279      END DO
280
281      DEALLOCATE(mod_tide,phi_tide)
282
283   END SUBROUTINE tide_init_elevation
284
285   SUBROUTINE tide_init_velocities( idx, td )
286      !!----------------------------------------------------------------------
287      !!                 ***  ROUTINE tide_init_elevation  ***
288      !!----------------------------------------------------------------------
289      TYPE(OBC_INDEX), INTENT( in )      ::   idx     ! OBC indices
290      TYPE(TIDES_DATA),INTENT( inout )      ::   td      ! tidal harmonics data
291      !! * Local declarations
292      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
293      INTEGER                            ::   itide, igrd, ib      ! dummy loop indices
294
295      igrd=2                                 ! U grid.
296
297      ALLOCATE(mod_tide(idx%nblenrim(igrd)),phi_tide(idx%nblenrim(igrd)))
298
299      DO itide = 1, nb_harmo
300         DO ib = 1, idx%nblenrim(igrd)
301            mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.)
302            phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
303         END DO
304         DO ib = 1, idx%nblenrim(igrd)
305            mod_tide(ib)=mod_tide(ib)*ftide(itide)
306            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
307         ENDDO
308         DO ib = 1, idx%nblenrim(igrd)
309            td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
310            td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
311         ENDDO
312      END DO
313
314      DEALLOCATE(mod_tide,phi_tide)
315
316      igrd=3                                 ! U grid.
317
318      ALLOCATE(mod_tide(idx%nblenrim(igrd)),phi_tide(idx%nblenrim(igrd)))
319
320      DO itide = 1, nb_harmo
321         DO ib = 1, idx%nblenrim(igrd)
322            mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.)
323            phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
324         END DO
325         DO ib = 1, idx%nblenrim(igrd)
326            mod_tide(ib)=mod_tide(ib)*ftide(itide)
327            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
328         ENDDO
329         DO ib = 1, idx%nblenrim(igrd)
330            td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
331            td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
332         ENDDO
333      END DO
334
335      DEALLOCATE(mod_tide,phi_tide)
336
337  END SUBROUTINE tide_init_velocities
338#else
339   !!----------------------------------------------------------------------
340   !!   Dummy module         NO Unstruct Open Boundary Conditions for tides
341   !!----------------------------------------------------------------------
342CONTAINS
343   SUBROUTINE bdytide_init                ! Empty routine
344   END SUBROUTINE bdytide_init
345   SUBROUTINE tide_data                ! Empty routine
346   END SUBROUTINE tide_data
347   SUBROUTINE tide_update( kt, jit )   ! Empty routine
348      WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, jit
349   END SUBROUTINE tide_update
350#endif
351
352   !!======================================================================
353END MODULE bdytides
Note: See TracBrowser for help on using the repository browser.