- Timestamp:
- 2012-04-25T12:21:21+02:00 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90
r3294 r3367 8 8 !! 3.0 ! 2008-04 (NEMO team) add in the reference version 9 9 !! 3.3 ! 2010-09 (D.Storkey and E.O'Dea) bug fixes 10 !! 3.4 ! 2011 (D. Storkey) rewrite in preparation for OBC-BDY merge11 10 !!---------------------------------------------------------------------- 12 11 #if defined key_bdy … … 15 14 !!---------------------------------------------------------------------- 16 15 !! PUBLIC 17 !! tide_init : read of namelist and initialisation of tidal harmonics data16 !! bdytide_init : read of namelist and initialisation of tidal harmonics data 18 17 !! tide_update : calculation of tidal forcing at each timestep 19 !! PRIVATE20 !! uvset :\21 !! vday : | Routines to correct tidal harmonics forcing for22 !! shpen : | start time of integration23 !! ufset : |24 !! vset :/25 18 !!---------------------------------------------------------------------- 26 19 USE timing ! Timing … … 34 27 USE bdy_oce ! ocean open boundary conditions 35 28 USE daymod ! calendar 29 USE tideini 30 USE tide_mod 36 31 USE fldread, ONLY: fld_map 37 32 … … 39 34 PRIVATE 40 35 41 PUBLIC tide_init! routine called in nemo_init36 PUBLIC bdytide_init ! routine called in nemo_init 42 37 PUBLIC tide_update ! routine called in bdydyn 43 38 44 39 TYPE, PUBLIC :: TIDES_DATA !: Storage for external tidal harmonics data 45 INTEGER :: ncpt !: Actual number of tidal components 46 REAL(wp), POINTER, DIMENSION(:) :: speed !: Phase speed of tidal constituent (deg/hr) 47 REAL(wp), POINTER, DIMENSION(:,:,:) :: ssh !: Tidal constituents : SSH 48 REAL(wp), POINTER, DIMENSION(:,:,:) :: u !: Tidal constituents : U 49 REAL(wp), POINTER, DIMENSION(:,:,:) :: v !: Tidal constituents : V 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.) 50 46 END TYPE TIDES_DATA 51 47 52 INTEGER, PUBLIC, PARAMETER :: jptides_max = 15 !: Max number of tidal contituents 53 54 TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides !: External tidal harmonics data 48 TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides !: External tidal harmonics data 55 49 56 50 !!---------------------------------------------------------------------- … … 61 55 CONTAINS 62 56 63 SUBROUTINE tide_init64 !!---------------------------------------------------------------------- 65 !! *** SUBROUTINE tide_init ***57 SUBROUTINE bdytide_init 58 !!---------------------------------------------------------------------- 59 !! *** SUBROUTINE bdytide_init *** 66 60 !! 67 61 !! ** Purpose : - Read in namelist for tides and initialise external … … 71 65 !! namelist variables 72 66 !!------------------- 73 LOGICAL :: ln_tide_date !: =T correct tide phases and amplitude for model start date 74 CHARACTER(len=80) :: filtide !: Filename root for tidal input files 75 CHARACTER(len= 4), DIMENSION(jptides_max) :: tide_cpt !: Names of tidal components used. 76 REAL(wp), DIMENSION(jptides_max) :: tide_speed !: Phase speed of tidal constituent (deg/hr) 77 !! 78 INTEGER, DIMENSION(jptides_max) :: nindx !: index of pre-set tidal components 79 INTEGER :: ib_bdy, itide, ib !: dummy loop indices 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 80 71 INTEGER :: inum, igrd 81 INTEGER, DIMENSION(3) :: ilen0 !: length of boundary data (from OBC arrays) 82 CHARACTER(len=80) :: clfile !: full file name for tidal input file 83 REAL(wp) :: z_arg, z_atde, z_btde, z1t, z2t 84 REAL(wp),DIMENSION(jptides_max) :: z_vplu, z_ftc 85 REAL(wp),ALLOCATABLE, DIMENSION(:,:,:) :: dta_read !: work space to read in tidal harmonics data 86 !! 87 TYPE(TIDES_DATA), POINTER :: td !: local short cut 88 !! 89 NAMELIST/nambdy_tide/ln_tide_date, filtide, tide_cpt, tide_speed 90 !!---------------------------------------------------------------------- 91 92 IF( nn_timing == 1 ) CALL timing_start('tide_init') 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') 93 82 94 83 IF(lwp) WRITE(numout,*) 95 IF(lwp) WRITE(numout,*) ' tide_init : initialization of tidal harmonic forcing at open boundaries'84 IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries' 96 85 IF(lwp) WRITE(numout,*) '~~~~~~~~~' 97 86 … … 103 92 104 93 ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries 105 ln_tide_date = .false.106 94 filtide(:) = '' 107 tide_speed(:) = 0.0108 tide_cpt(:) = ''109 95 110 96 ! Don't REWIND here - may need to read more than one of these namelists. 111 97 READ ( numnam, nambdy_tide ) 112 ! ! Count number of components specified113 td%ncpt = 0114 DO itide = 1, jptides_max115 IF( tide_cpt(itide) /= '' ) THEN116 td%ncpt = td%ncpt + 1117 ENDIF118 END DO119 120 ! Fill in phase speeds from namelist121 ALLOCATE( td%speed(td%ncpt) )122 td%speed = tide_speed(1:td%ncpt)123 124 ! Find constituents in standard list125 DO itide = 1, td%ncpt126 nindx(itide) = 0127 IF( TRIM( tide_cpt(itide) ) == 'Q1' ) nindx(itide) = 1128 IF( TRIM( tide_cpt(itide) ) == 'O1' ) nindx(itide) = 2129 IF( TRIM( tide_cpt(itide) ) == 'P1' ) nindx(itide) = 3130 IF( TRIM( tide_cpt(itide) ) == 'S1' ) nindx(itide) = 4131 IF( TRIM( tide_cpt(itide) ) == 'K1' ) nindx(itide) = 5132 IF( TRIM( tide_cpt(itide) ) == '2N2' ) nindx(itide) = 6133 IF( TRIM( tide_cpt(itide) ) == 'MU2' ) nindx(itide) = 7134 IF( TRIM( tide_cpt(itide) ) == 'N2' ) nindx(itide) = 8135 IF( TRIM( tide_cpt(itide) ) == 'NU2' ) nindx(itide) = 9136 IF( TRIM( tide_cpt(itide) ) == 'M2' ) nindx(itide) = 10137 IF( TRIM( tide_cpt(itide) ) == 'L2' ) nindx(itide) = 11138 IF( TRIM( tide_cpt(itide) ) == 'T2' ) nindx(itide) = 12139 IF( TRIM( tide_cpt(itide) ) == 'S2' ) nindx(itide) = 13140 IF( TRIM( tide_cpt(itide) ) == 'K2' ) nindx(itide) = 14141 IF( TRIM( tide_cpt(itide) ) == 'M4' ) nindx(itide) = 15142 IF( nindx(itide) == 0 .AND. lwp ) THEN143 WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list'144 CALL ctl_warn( ctmp1 )145 ENDIF146 END DO147 98 ! ! Parameter control and print 148 IF( td%ncpt < 1 ) THEN 149 CALL ctl_stop( ' Did not find any tidal components in namelist nambdy_tide' ) 150 ELSE 151 IF(lwp) WRITE(numout,*) ' Namelist nambdy_tide : tidal harmonic forcing at open boundaries' 152 IF(lwp) WRITE(numout,*) ' tidal components specified ', td%ncpt 153 IF(lwp) WRITE(numout,*) ' ', tide_cpt(1:td%ncpt) 154 IF(lwp) WRITE(numout,*) ' associated phase speeds (deg/hr) : ' 155 IF(lwp) WRITE(numout,*) ' ', tide_speed(1:td%ncpt) 156 ENDIF 157 158 159 ! Allocate space for tidal harmonics data - 160 ! get size from OBC data arrays 161 ! --------------------------------------- 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 ! ----------------------------------------------------------------------- 162 107 163 108 ilen0(1) = SIZE( dta_bdy(ib_bdy)%ssh ) 164 ALLOCATE( td%ssh( ilen0(1), td%ncpt, 2 ) ) 109 ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) ) 110 ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) ) 165 111 166 112 ilen0(2) = SIZE( dta_bdy(ib_bdy)%u2d ) 167 ALLOCATE( td%u( ilen0(2), td%ncpt, 2 ) ) 113 ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) ) 114 ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) ) 168 115 169 116 ilen0(3) = SIZE( dta_bdy(ib_bdy)%v2d ) 170 ALLOCATE( td%v( ilen0(3), td%ncpt, 2 ) ) 117 ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) ) 118 ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) ) 171 119 172 120 ALLOCATE( dta_read( MAXVAL(ilen0), 1, 1 ) ) 173 174 121 175 122 ! Open files and read in tidal forcing data 176 123 ! ----------------------------------------- 177 124 178 DO itide = 1, td%ncpt125 DO itide = 1, nb_harmo 179 126 ! ! SSH fields 180 clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc' 181 IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 127 clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc' 182 128 CALL iom_open( clfile, inum ) 183 129 CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 184 td%ssh (:,itide,1) = dta_read(1:ilen0(1),1,1)130 td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1) 185 131 CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 186 td%ssh (:,itide,2) = dta_read(1:ilen0(1),1,1)132 td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1) 187 133 CALL iom_close( inum ) 188 134 ! ! U fields 189 clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc' 190 IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 135 clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc' 191 136 CALL iom_open( clfile, inum ) 192 137 CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 193 td%u (:,itide,1) = dta_read(1:ilen0(2),1,1)138 td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1) 194 139 CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 195 td%u (:,itide,2) = dta_read(1:ilen0(2),1,1)140 td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1) 196 141 CALL iom_close( inum ) 197 142 ! ! V fields 198 clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc' 199 IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 143 clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc' 200 144 CALL iom_open( clfile, inum ) 201 145 CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 202 td%v (:,itide,1) = dta_read(1:ilen0(3),1,1)146 td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1) 203 147 CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 204 td%v (:,itide,2) = dta_read(1:ilen0(3),1,1)148 td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1) 205 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 206 156 ! 207 157 END DO ! end loop on tidal components 208 209 IF( ln_tide_date ) THEN ! correct for date factors 210 211 !! used nmonth, nyear and nday from daymod.... 212 ! Calculate date corrects for 15 standard consituents 213 ! This is the initialisation step, so nday, nmonth, nyear are the 214 ! initial date/time of the integration. 215 print *, nday,nmonth,nyear 216 nyear = int(ndate0 / 10000 ) ! initial year 217 nmonth = int((ndate0 - nyear * 10000 ) / 100 ) ! initial month 218 nday = int(ndate0 - nyear * 10000 - nmonth * 100) 219 220 CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu ) 221 222 IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear 223 224 DO itide = 1, td%ncpt ! loop on tidal components 225 ! 226 IF( nindx(itide) /= 0 ) THEN 227 !!gm use rpi and rad global variable 228 z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0 229 z_atde=z_ftc(nindx(itide))*cos(z_arg) 230 z_btde=z_ftc(nindx(itide))*sin(z_arg) 231 IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), td%speed(itide), & 232 & z_ftc(nindx(itide)), z_vplu(nindx(itide)) 233 ELSE 234 z_atde = 1.0_wp 235 z_btde = 0.0_wp 236 ENDIF 237 ! ! elevation 238 igrd = 1 239 DO ib = 1, ilen0(igrd) 240 z1t = z_atde * td%ssh(ib,itide,1) + z_btde * td%ssh(ib,itide,2) 241 z2t = z_atde * td%ssh(ib,itide,2) - z_btde * td%ssh(ib,itide,1) 242 td%ssh(ib,itide,1) = z1t 243 td%ssh(ib,itide,2) = z2t 244 END DO 245 ! ! u 246 igrd = 2 247 DO ib = 1, ilen0(igrd) 248 z1t = z_atde * td%u(ib,itide,1) + z_btde * td%u(ib,itide,2) 249 z2t = z_atde * td%u(ib,itide,2) - z_btde * td%u(ib,itide,1) 250 td%u(ib,itide,1) = z1t 251 td%u(ib,itide,2) = z2t 252 END DO 253 ! ! v 254 igrd = 3 255 DO ib = 1, ilen0(igrd) 256 z1t = z_atde * td%v(ib,itide,1) + z_btde * td%v(ib,itide,2) 257 z2t = z_atde * td%v(ib,itide,2) - z_btde * td%v(ib,itide,1) 258 td%v(ib,itide,1) = z1t 259 td%v(ib,itide,2) = z2t 260 END DO 261 ! 262 END DO ! end loop on tidal components 263 ! 264 ENDIF ! date correction 158 ! 159 DEALLOCATE( dta_read ) 265 160 ! 266 161 ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2 … … 268 163 END DO ! loop on ib_bdy 269 164 270 IF( nn_timing == 1 ) CALL timing_stop('tide_init') 271 272 END SUBROUTINE tide_init 273 165 IF( nn_timing == 1 ) CALL timing_stop('bdytide_init') 166 167 END SUBROUTINE bdytide_init 274 168 275 169 SUBROUTINE tide_update ( kt, idx, dta, td, jit, time_offset ) … … 280 174 !! 281 175 !!---------------------------------------------------------------------- 282 INTEGER, INTENT( in ) :: kt ! Main timestep counter 283 !!gm doctor jit ==> kit 284 TYPE(OBC_INDEX), INTENT( in ) :: idx ! OBC indices 285 TYPE(OBC_DATA), INTENT(inout) :: dta ! OBC external data 286 TYPE(TIDES_DATA),INTENT( in ) :: td ! tidal harmonics data 287 INTEGER,INTENT(in),OPTIONAL :: jit ! Barotropic timestep counter (for timesplitting option) 288 INTEGER,INTENT( in ), OPTIONAL :: time_offset ! time offset in units of timesteps. NB. if jit 289 ! is present then units = subcycle timesteps. 290 ! time_offset = 0 => get data at "now" time level 291 ! time_offset = -1 => get data at "before" time level 292 ! time_offset = +1 => get data at "after" time level 293 ! etc. 294 !! 295 INTEGER :: itide, igrd, ib ! dummy loop indices 296 INTEGER :: time_add ! time offset in units of timesteps 297 REAL(wp) :: z_arg, z_sarg 298 REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost 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 299 192 !!---------------------------------------------------------------------- 300 193 301 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 302 215 303 216 time_add = 0 … … 306 219 ENDIF 307 220 308 ! Note tide phase speeds are in deg/hour, so we need to convert the309 ! elapsed time in seconds to hours by dividing by 3600.0310 221 IF( PRESENT(jit) ) THEN 311 z_arg = ( ( kt-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) ) * rad / 3600.0222 z_arg = ( ((kt-kt_tide)-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) ) 312 223 ELSE 313 z_arg = ( kt+time_add) * rdt * rad / 3600.0224 z_arg = ((kt-kt_tide)+time_add) * rdt 314 225 ENDIF 315 226 316 DO itide = 1, td%ncpt317 z_sarg = z_arg * td%speed(itide)227 DO itide = 1, nb_harmo 228 z_sarg = z_arg * omega_tide(itide) 318 229 z_cost(itide) = COS( z_sarg ) 319 230 z_sist(itide) = SIN( z_sarg ) 320 231 END DO 321 232 322 DO itide = 1, td%ncpt323 igrd=1 ! SSH on tracer grid .233 DO itide = 1, nb_harmo 234 igrd=1 ! SSH on tracer grid 324 235 DO ib = 1, idx%nblenrim(igrd) 325 236 dta%ssh(ib) = dta%ssh(ib) + td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide) 326 ! if(lwp) write(numout,*) 'z', ib, itide, dta%ssh(ib), td%ssh(ib,itide,1),td%ssh(ib,itide,2)327 237 END DO 328 238 igrd=2 ! U grid 329 239 DO ib=1, idx%nblenrim(igrd) 330 dta%u2d(ib) = dta%u2d(ib) + td%u(ib,itide,1)*z_cost(itide) + td%u(ib,itide,2)*z_sist(itide) 331 ! if(lwp) write(numout,*) 'u',ib,itide,utide(ib), td%u(ib,itide,1),td%u(ib,itide,2) 240 dta%u2d(ib) = dta%u2d(ib) + td%u (ib,itide,1)*z_cost(itide) + td%u (ib,itide,2)*z_sist(itide) 332 241 END DO 333 242 igrd=3 ! V grid 334 243 DO ib=1, idx%nblenrim(igrd) 335 dta%v2d(ib) = dta%v2d(ib) + td%v(ib,itide,1)*z_cost(itide) + td%v(ib,itide,2)*z_sist(itide) 336 ! if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), td%v(ib,itide,1),td%v(ib,itide,2) 244 dta%v2d(ib) = dta%v2d(ib) + td%v (ib,itide,1)*z_cost(itide) + td%v (ib,itide,2)*z_sist(itide) 337 245 END DO 338 246 END DO … … 342 250 END SUBROUTINE tide_update 343 251 344 !!gm doctor naming of dummy argument variables!!! and all local variables 345 346 SUBROUTINE uvset( ihs, iday, imnth, iyr, f, z_vplu ) 347 !!---------------------------------------------------------------------- 348 !! *** SUBROUTINE uvset *** 349 !! 350 !! ** Purpose : - adjust tidal forcing for date factors 351 !! 352 !!---------------------------------------------------------------------- 353 implicit none 354 INTEGER, INTENT( in ) :: ihs ! Start time hours 355 INTEGER, INTENT( in ) :: iday ! start time days 356 INTEGER, INTENT( in ) :: imnth ! start time month 357 INTEGER, INTENT( in ) :: iyr ! start time year 358 !! 359 !!gm nc is jptides_max ???? 360 INTEGER , PARAMETER :: nc =15 ! maximum number of constituents 361 CHARACTER(len=8), DIMENSION(nc) :: cname 362 INTEGER :: year, vd, ivdy, ndc, i, k 363 REAL(wp) :: ss, h, p, en, p1, rtd 364 REAL(wp), DIMENSION(nc) :: f ! nodal correction 365 REAL(wp), DIMENSION(nc) :: z_vplu ! phase correction 366 REAL(wp), DIMENSION(nc) :: u, v, zig 367 !! 368 DATA cname/ 'q1' , 'o1' , 'p1' , 's1' , 'k1' , & 369 & '2n2' , 'mu2' , 'n2' , 'nu2' , 'm2' , & 370 & 'l2' , 't2' , 's2' , 'k2' , 'm4' / 371 DATA zig/ .2338507481, .2433518789, .2610826055, .2617993878, .2625161701, & 372 & .4868657873, .4881373225, .4963669182, .4976384533, .5058680490, & 373 & .5153691799, .5228820265, .5235987756, .5250323419, 1.011736098 / 374 !!---------------------------------------------------------------------- 375 ! 376 ! ihs - start time gmt on ... 377 ! iday/imnth/iyr - date e.g. 12/10/87 378 ! 379 CALL vday(iday,imnth,iyr,ivdy) 380 vd = ivdy 381 ! 382 !rp note change of year number for d. blackman shpen 383 !rp if(iyr.ge.1000) year=iyr-1900 384 !rp if(iyr.lt.1000) year=iyr 385 year = iyr 386 ! 387 !.....year = year of required data 388 !.....vd = day of required data..set up for 0000gmt day year 389 ndc = nc 390 !.....ndc = number of constituents allowed 391 !!gm use rpi ? 392 rtd = 360.0 / 6.2831852 393 DO i = 1, ndc 394 zig(i) = zig(i)*rtd 395 ! sigo(i)= zig(i) 396 END DO 397 398 !!gm try to avoid RETURN in F90 399 IF( year == 0 ) RETURN 400 401 CALL shpen( year, vd, ss, h , p , en, p1 ) 402 CALL ufset( p , en, u , f ) 403 CALL vset ( ss , h , p , en, p1, v ) 404 ! 405 DO k = 1, nc 406 z_vplu(k) = v(k) + u(k) 407 z_vplu(k) = z_vplu(k) + dble(ihs) * zig(k) 408 DO WHILE( z_vplu(k) < 0 ) 409 z_vplu(k) = z_vplu(k) + 360.0 410 END DO 411 DO WHILE( z_vplu(k) > 360. ) 412 z_vplu(k) = z_vplu(k) - 360.0 413 END DO 414 END DO 415 ! 416 END SUBROUTINE uvset 417 418 419 SUBROUTINE vday( iday, imnth, iy, ivdy ) 420 !!---------------------------------------------------------------------- 421 !! *** SUBROUTINE vday *** 422 !! 423 !! ** Purpose : - adjust tidal forcing for date factors 424 !! 425 !!---------------------------------------------------------------------- 426 INTEGER, INTENT(in ) :: iday, imnth, iy ! ???? 427 INTEGER, INTENT( out) :: ivdy ! ??? 428 !! 429 INTEGER :: iyr 430 !!------------------------------------------------------------------------------ 431 432 !!gm nday_year in day mode is the variable compiuted here, no? 433 !!gm nday_year , & !: curent day counted from jan 1st of the current year 434 435 !calculate day number in year from day/month/year 436 if(imnth.eq.1) ivdy=iday 437 if(imnth.eq.2) ivdy=iday+31 438 if(imnth.eq.3) ivdy=iday+59 439 if(imnth.eq.4) ivdy=iday+90 440 if(imnth.eq.5) ivdy=iday+120 441 if(imnth.eq.6) ivdy=iday+151 442 if(imnth.eq.7) ivdy=iday+181 443 if(imnth.eq.8) ivdy=iday+212 444 if(imnth.eq.9) ivdy=iday+243 445 if(imnth.eq.10) ivdy=iday+273 446 if(imnth.eq.11) ivdy=iday+304 447 if(imnth.eq.12) ivdy=iday+334 448 iyr=iy 449 if(mod(iyr,4).eq.0.and.imnth.gt.2) ivdy=ivdy+1 450 if(mod(iyr,100).eq.0.and.imnth.gt.2) ivdy=ivdy-1 451 if(mod(iyr,400).eq.0.and.imnth.gt.2) ivdy=ivdy+1 452 ! 453 END SUBROUTINE vday 454 455 !!doctor norme for dummy arguments 456 457 SUBROUTINE shpen( year, vd, s, h, p, en, p1 ) 458 !!---------------------------------------------------------------------- 459 !! *** SUBROUTINE shpen *** 460 !! 461 !! ** Purpose : - calculate astronomical arguments for tides 462 !! this version from d. blackman 30 nove 1990 463 !! 464 !!---------------------------------------------------------------------- 465 !!gm add INTENT in, out or inout.... DOCTOR name.... 466 !!gm please do not use variable name with a single letter: impossible to search in a code 467 INTEGER :: year,vd 468 REAL(wp) :: s,h,p,en,p1 469 !! 470 INTEGER :: yr,ilc,icent,it,iday,ild,ipos,nn,iyd 471 REAL(wp) :: cycle,t,td,delt(84),delta,deltat 472 !! 473 DATA delt /-5.04, -3.90, -2.87, -0.58, 0.71, 1.80, & 474 & 3.08, 4.63, 5.86, 7.21, 8.58, 10.50, 12.10, & 475 & 12.49, 14.41, 15.59, 15.81, 17.52, 19.01, 18.39, & 476 & 19.55, 20.36, 21.01, 21.81, 21.76, 22.35, 22.68, & 477 & 22.94, 22.93, 22.69, 22.94, 23.20, 23.31, 23.63, & 478 & 23.47, 23.68, 23.62, 23.53, 23.59, 23.99, 23.80, & 479 & 24.20, 24.99, 24.97, 25.72, 26.21, 26.37, 26.89, & 480 & 27.68, 28.13, 28.94, 29.42, 29.66, 30.29, 30.96, & 481 & 31.09, 31.59, 31.52, 31.92, 32.45, 32.91, 33.39, & 482 & 33.80, 34.23, 34.73, 35.40, 36.14, 36.99, 37.87, & 483 & 38.75, 39.70, 40.70, 41.68, 42.82, 43.96, 45.00, & 484 & 45.98, 47.00, 48.03, 49.10, 50.10, 50.97, 51.81, & 485 & 52.57 / 486 !!---------------------------------------------------------------------- 487 488 cycle = 360.0 489 ilc = 0 490 icent = year / 100 491 yr = year - icent * 100 492 t = icent - 20 493 ! 494 ! for the following equations 495 ! time origin is fixed at 00 hr of jan 1st,2000. 496 ! see notes by cartwright 497 ! 498 !!gm old coding style, use CASE instead and avoid GOTO (obsolescence in fortran 90) 499 !!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 500 it = icent - 20 501 if (it) 1,2,2 502 1 iday = it/4 -it 503 go to 3 504 2 iday = (it+3)/4 - it 505 ! 506 ! t is in julian century 507 ! correction in gegorian calander where only century year divisible 508 ! by 4 is leap year. 509 ! 510 3 continue 511 ! 512 td = 0.0 513 ! 514 !!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 515 if (yr) 4,5,4 516 ! 517 4 iyd = 365*yr 518 ild = (yr-1)/4 519 if((icent - (icent/4)*4) .eq. 0) ilc = 1 520 td = iyd + ild + ilc 521 ! 522 5 td = td + iday + vd -1.0 - 0.5 523 t = t + (td/36525.0) 524 ! 525 ipos=year-1899 526 if (ipos .lt. 0) go to 7 527 if (ipos .gt. 83) go to 6 528 ! 529 delta = (delt(ipos+1)+delt(ipos))/2.0 530 go to 7 531 ! 532 6 delta= (65.0-50.5)/20.0*(year-1980)+50.5 533 ! 534 7 deltat = delta * 1.0e-6 535 ! 536 !!gm precision of the computation : example for s it should be replace by: 537 !!gm s = 218.3165 + (481267.8813 - 0.0016*t)*t + 152.0*deltat ==> more precise modify the last digits results 538 s = 218.3165 + 481267.8813*t - 0.0016*t*t + 152.0*deltat 539 h = 280.4661 + 36000.7698 *t + 0.0003*t*t + 11.0*deltat 540 p = 83.3535 + 4069.0139 *t - 0.0103*t*t + deltat 541 en = 234.9555 + 1934.1363 *t - 0.0021*t*t + deltat 542 p1 = 282.9384 + 1.7195 *t + 0.0005*t*t 543 ! 544 nn = s / cycle 545 s = s - nn * cycle 546 IF( s < 0.e0 ) s = s + cycle 547 ! 548 nn = h / cycle 549 h = h - cycle * nn 550 IF( h < 0.e0 ) h = h + cycle 551 ! 552 nn = p / cycle 553 p = p - cycle * nn 554 IF( p < 0.e0) p = p + cycle 555 ! 556 nn = en / cycle 557 en = en - cycle * nn 558 IF( en < 0.e0 ) en = en + cycle 559 en = cycle - en 560 ! 561 nn = p1 / cycle 562 p1 = p1 - nn * cycle 563 ! 564 END SUBROUTINE shpen 565 566 567 SUBROUTINE ufset( p, cn, b, a ) 568 !!---------------------------------------------------------------------- 569 !! *** SUBROUTINE ufset *** 570 !! 571 !! ** Purpose : - calculate nodal parameters for the tides 572 !! 573 !!---------------------------------------------------------------------- 574 !!gm doctor naming of dummy argument variables!!! and all local variables 575 !!gm nc is jptides_max ???? 576 integer nc 577 parameter (nc=15) 578 REAL(wp) p,cn 579 !! 580 581 !!gm rad is already a public variable defined in phycst.F90 .... ==> doctor norme local real start with "z" 582 REAL(wp) :: w1, w2, w3, w4, w5, w6, w7, w8, nw, pw, rad 583 REAL(wp) :: a(nc), b(nc) 584 INTEGER :: ndc, k 585 !!---------------------------------------------------------------------- 586 587 ndc = nc 588 589 ! a=f , b =u 590 ! t is zero as compared to tifa. 591 !! use rad defined in phycst (i.e. add a USE phycst at the begining of the module 592 rad = 6.2831852d0/360.0 593 pw = p * rad 594 nw = cn * rad 595 w1 = cos( nw ) 596 w2 = cos( 2*nw ) 597 w3 = cos( 3*nw ) 598 w4 = sin( nw ) 599 w5 = sin( 2*nw ) 600 w6 = sin( 3*nw ) 601 w7 = 1. - 0.2505 * COS( 2*pw ) - 0.1102 * COS(2*pw-nw ) & 602 & - 0.156 * COS( 2*pw-2*nw ) - 0.037 * COS( nw ) 603 w8 = - 0.2505 * SIN( 2*pw ) - 0.1102 * SIN(2*pw-nw ) & 604 & - 0.0156 * SIN( 2*pw-2*nw ) - 0.037 * SIN( nw ) 605 ! 606 a(1) = 1.0089 + 0.1871 * w1 - 0.0147 * w2 + 0.0014 * w3 607 b(1) = 0.1885 * w4 - 0.0234 * w5 + 0.0033 * w6 608 ! q1 609 a(2) = a(1) 610 b(2) = b(1) 611 ! o1 612 a(3) = 1.0 613 b(3) = 0.0 614 ! p1 615 a(4) = 1.0 616 b(4) = 0.0 617 ! s1 618 a(5) = 1.0060+0.1150*w1- 0.0088*w2 +0.0006*w3 619 b(5) = -0.1546*w4 + 0.0119*w5 -0.0012*w6 620 ! k1 621 a(6) =1.0004 -0.0373*w1+ 0.0002*w2 622 b(6) = -0.0374*w4 623 ! 2n2 624 a(7) = a(6) 625 b(7) = b(6) 626 ! mu2 627 a(8) = a(6) 628 b(8) = b(6) 629 ! n2 630 a(9) = a(6) 631 b(9) = b(6) 632 ! nu2 633 a(10) = a(6) 634 b(10) = b(6) 635 ! m2 636 a(11) = SQRT( w7 * w7 + w8 * w8 ) 637 b(11) = ATAN( w8 / w7 ) 638 !!gmuse rpi instead of 3.141992 ??? true pi is rpi=3.141592653589793_wp ..... ???? 639 IF( w7 < 0.e0 ) b(11) = b(11) + 3.141992 640 ! l2 641 a(12) = 1.0 642 b(12) = 0.0 643 ! t2 644 a(13)= a(12) 645 b(13)= b(12) 646 ! s2 647 a(14) = 1.0241+0.2863*w1+0.0083*w2 -0.0015*w3 648 b(14) = -0.3096*w4 + 0.0119*w5 - 0.0007*w6 649 ! k2 650 a(15) = a(6)*a(6) 651 b(15) = 2*b(6) 652 ! m4 653 !!gm old coding, remove GOTO and label of lines 654 !!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 655 DO 40 k = 1,ndc 656 b(k) = b(k)/rad 657 32 if (b(k)) 34,35,35 658 34 b(k) = b(k) + 360.0 659 go to 32 660 35 if (b(k)-360.0) 40,37,37 661 37 b(k) = b(k)-360.0 662 go to 35 663 40 continue 664 ! 665 END SUBROUTINE ufset 666 667 668 SUBROUTINE vset( s,h,p,en,p1,v) 669 !!---------------------------------------------------------------------- 670 !! *** SUBROUTINE vset *** 671 !! 672 !! ** Purpose : - calculate tidal phases for 0000gmt on start day of run 673 !! 674 !!---------------------------------------------------------------------- 675 !!gm doctor naming of dummy argument variables!!! and all local variables 676 !!gm nc is jptides_max ???? 677 !!gm en argument is not used: suppress it ? 678 integer nc 679 parameter (nc=15) 680 real(wp) s,h,p,en,p1 681 real(wp) v(nc) 682 !! 683 integer ndc, k 684 !!---------------------------------------------------------------------- 685 686 ndc = nc 687 ! v s are computed here. 688 v(1) =-3*s +h +p +270 ! Q1 689 v(2) =-2*s +h +270.0 ! O1 690 v(3) =-h +270 ! P1 691 v(4) =180 ! S1 692 v(5) =h +90.0 ! K1 693 v(6) =-4*s +2*h +2*p ! 2N2 694 v(7) =-4*(s-h) ! MU2 695 v(8) =-3*s +2*h +p ! N2 696 v(9) =-3*s +4*h -p ! MU2 697 v(10) =-2*s +2*h ! M2 698 v(11) =-s +2*h -p +180 ! L2 699 v(12) =-h +p1 ! T2 700 v(13) =0 ! S2 701 v(14) =h+h ! K2 702 v(15) =2*v(10) ! M4 703 ! 704 !!gm old coding, remove GOTO and label of lines 705 !!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 706 do 72 k = 1, ndc 707 69 if( v(k) ) 70,71,71 708 70 v(k) = v(k)+360.0 709 go to 69 710 71 if( v(k) - 360.0 ) 72,73,73 711 73 v(k) = v(k)-360.0 712 go to 71 713 72 continue 714 ! 715 END SUBROUTINE vset 716 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 717 338 #else 718 339 !!---------------------------------------------------------------------- 719 340 !! Dummy module NO Unstruct Open Boundary Conditions for tides 720 341 !!---------------------------------------------------------------------- 721 !!gm are you sure we need to define filtide and tide_cpt ?722 CHARACTER(len=80), PUBLIC :: filtide !: Filename root for tidal input files723 CHARACTER(len=4 ), PUBLIC, DIMENSION(1) :: tide_cpt !: Names of tidal components used.724 725 342 CONTAINS 726 SUBROUTINE tide_init ! Empty routine727 END SUBROUTINE tide_init343 SUBROUTINE bdytide_init ! Empty routine 344 END SUBROUTINE bdytide_init 728 345 SUBROUTINE tide_data ! Empty routine 729 346 END SUBROUTINE tide_data 730 SUBROUTINE tide_update( kt, kit ) ! Empty routine731 WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, kit347 SUBROUTINE tide_update( kt, jit ) ! Empty routine 348 WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, jit 732 349 END SUBROUTINE tide_update 733 350 #endif
Note: See TracChangeset
for help on using the changeset viewer.
