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.
sbctide.F90 in NEMO/branches/NERC/NEMO_4.0.4_CO9_package_tides/src/OCE/SBC – NEMO

source: NEMO/branches/NERC/NEMO_4.0.4_CO9_package_tides/src/OCE/SBC/sbctide.F90 @ 14266

Last change on this file since 14266 was 14266, checked in by dbyrne, 3 years ago

Long period tide forcing, variable love number, new nodal equation 20 added
File size: 7.4 KB
Line 
1MODULE sbctide
2   !!======================================================================
3   !!                       ***  MODULE  sbctide  ***
4   !! Initialization of tidal forcing
5   !!======================================================================
6   !! History :  9.0  !  2007  (O. Le Galloudec)  Original code
7   !!----------------------------------------------------------------------
8   USE oce            ! ocean dynamics and tracers variables
9   USE dom_oce        ! ocean space and time domain
10   USE phycst         ! physical constant
11   USE daymod         ! calandar
12   USE tideini        !
13   !
14   USE in_out_manager ! I/O units
15   USE iom            ! xIOs server
16   USE ioipsl         ! NetCDF IPSL library
17   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
18
19   USE bdytides ! davbyr - Access to love number
20
21   IMPLICIT NONE
22   PUBLIC
23
24   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   pot_astro   !
25
26   !!----------------------------------------------------------------------
27   !!   tidal potential
28   !!----------------------------------------------------------------------
29   !!   sbc_tide            :
30   !!   tide_init_potential :
31   !!----------------------------------------------------------------------
32
33   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   amp_pot, phi_pot
34   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   amp_load, phi_load
35 
36   !!----------------------------------------------------------------------
37   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
38   !! $Id$
39   !! Software governed by the CeCILL license (see ./LICENSE)
40   !!----------------------------------------------------------------------
41CONTAINS
42
43   SUBROUTINE sbc_tide( kt )
44      !!----------------------------------------------------------------------
45      !!                 ***  ROUTINE sbc_tide  ***
46      !!----------------------------------------------------------------------     
47      INTEGER, INTENT( in ) ::   kt     ! ocean time-step
48      INTEGER               ::   jk     ! dummy loop index
49      INTEGER               ::   nsec_day_orig     ! Temporary variable
50      !!----------------------------------------------------------------------
51     
52      IF( nsec_day == NINT(0.5_wp * rdt) .OR. kt == nit000 ) THEN      ! start a new day
53         !
54         IF( kt == nit000 )THEN
55            ALLOCATE( amp_pot(jpi,jpj,nb_harmo),                      &
56               &      phi_pot(jpi,jpj,nb_harmo), pot_astro(jpi,jpj)   )
57            IF( ln_read_load )THEN
58               ALLOCATE( amp_load(jpi,jpj,nb_harmo), phi_load(jpi,jpj,nb_harmo) )
59               CALL tide_init_load
60            ENDIF
61         ENDIF
62         !
63         IF( ln_read_load )THEN
64            amp_pot(:,:,:) = amp_load(:,:,:)
65            phi_pot(:,:,:) = phi_load(:,:,:)
66         ELSE
67            amp_pot(:,:,:) = 0._wp
68            phi_pot(:,:,:) = 0._wp
69         ENDIF
70         pot_astro(:,:) = 0._wp
71         !
72         ! If the run does not start from midnight then need to initialise tides
73         ! at the start of the current day (only occurs when kt==nit000)
74         ! Temporarily set nsec_day to beginning of day.
75         nsec_day_orig = nsec_day
76         IF ( nsec_day /= NINT(0.5_wp * rdt) ) THEN
77            kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
78            nsec_day = NINT(0.5_wp * rdt)
79         ELSE
80            kt_tide = kt 
81         ENDIF
82         CALL tide_harmo( omega_tide, v0tide, utide, ftide, ntide, nb_harmo )
83         !
84         !
85         IF(lwp) THEN
86            WRITE(numout,*)
87            WRITE(numout,*) 'sbc_tide : Update of the components and (re)Init. the potential at kt=', kt
88            WRITE(numout,*) '~~~~~~~~ '
89            DO jk = 1, nb_harmo
90               WRITE(numout,*) Wave(ntide(jk))%cname_tide, utide(jk), ftide(jk), v0tide(jk), omega_tide(jk)
91            END DO
92         ENDIF
93         !
94         IF( ln_tide_pot )   CALL tide_init_potential
95         !
96         ! Reset nsec_day
97         nsec_day = nsec_day_orig 
98      ENDIF
99      !
100   END SUBROUTINE sbc_tide
101
102
103   SUBROUTINE tide_init_potential
104      !!----------------------------------------------------------------------
105      !!                 ***  ROUTINE tide_init_potential  ***
106      !!----------------------------------------------------------------------
107      INTEGER  ::   ji, jj, jk   ! dummy loop indices
108      REAL(wp) ::   zcons, ztmp1, ztmp2, zlat, zlon, ztmp, zamp, zcs   ! local scalar
109      !!----------------------------------------------------------------------
110
111      DO jk = 1, nb_harmo
112         ! davbyr - Insert variable Love number where once was 0.7
113         zcons = dn_love_number * Wave(ntide(jk))%equitide * ftide(jk)
114         ! END davbyr
115         DO ji = 1, jpi
116            DO jj = 1, jpj
117               ztmp1 =  ftide(jk) * amp_pot(ji,jj,jk) * COS( phi_pot(ji,jj,jk) + v0tide(jk) + utide(jk) )
118               ztmp2 = -ftide(jk) * amp_pot(ji,jj,jk) * SIN( phi_pot(ji,jj,jk) + v0tide(jk) + utide(jk) )
119               zlat = gphit(ji,jj)*rad !! latitude en radian
120               zlon = glamt(ji,jj)*rad !! longitude en radian
121               ztmp = v0tide(jk) + utide(jk) + Wave(ntide(jk))%nutide * zlon
122               ! le potentiel est composé des effets des astres:
123               IF    ( Wave(ntide(jk))%nutide == 1 )  THEN  ;  zcs = zcons * SIN( 2._wp*zlat )
124               ELSEIF( Wave(ntide(jk))%nutide == 2 )  THEN  ;  zcs = zcons * COS( zlat )**2
125               ! davbyr - Include long period tidal forcing
126               ELSEIF( Wave(ntide(jk))%nutide == 0 )  THEN  ;  zcs = zcons * (0.5_wp-1.5_wp*SIN(zlat)**2._wp)
127               ! END - davbyr
128               ELSE                                         ;  zcs = 0._wp
129               ENDIF
130               ztmp1 = ztmp1 + zcs * COS( ztmp )
131               ztmp2 = ztmp2 - zcs * SIN( ztmp )
132               zamp = SQRT( ztmp1*ztmp1 + ztmp2*ztmp2 )
133               amp_pot(ji,jj,jk) = zamp
134               phi_pot(ji,jj,jk) = ATAN2( -ztmp2 / MAX( 1.e-10_wp , zamp ) ,   &
135                  &                        ztmp1 / MAX( 1.e-10_wp,  zamp )   )
136            END DO
137         END DO
138      END DO
139      !
140   END SUBROUTINE tide_init_potential
141
142   SUBROUTINE tide_init_load
143      !!----------------------------------------------------------------------
144      !!                 ***  ROUTINE tide_init_load  ***
145      !!----------------------------------------------------------------------
146      INTEGER :: inum                 ! Logical unit of input file
147      INTEGER :: ji, jj, itide        ! dummy loop indices
148      REAL(wp), DIMENSION(jpi,jpj) ::   ztr, zti   !: workspace to read in tidal harmonics data
149      !!----------------------------------------------------------------------
150      IF(lwp) THEN
151         WRITE(numout,*)
152         WRITE(numout,*) 'tide_init_load : Initialization of load potential from file'
153         WRITE(numout,*) '~~~~~~~~~~~~~~ '
154      ENDIF
155      !
156      CALL iom_open ( cn_tide_load , inum )
157      !
158      DO itide = 1, nb_harmo
159         CALL iom_get  ( inum, jpdom_data,TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) )
160         CALL iom_get  ( inum, jpdom_data,TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) )
161         !
162         DO ji=1,jpi
163            DO jj=1,jpj
164               amp_load(ji,jj,itide) =  SQRT( ztr(ji,jj)**2. + zti(ji,jj)**2. )
165               phi_load(ji,jj,itide) = ATAN2(-zti(ji,jj), ztr(ji,jj) )
166            END DO
167         END DO
168         !
169      END DO
170      CALL iom_close( inum )
171      !
172   END SUBROUTINE tide_init_load
173
174  !!======================================================================
175END MODULE sbctide
Note: See TracBrowser for help on using the repository browser.