- Timestamp:
- 2012-11-26T11:46:39+01:00 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2012/dev_NOC_MERCATOR_2012/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90
r3294 r3651 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 ! 201 1 (D. Storkey) rewrite in preparation for OBC-BDY merge10 !! 3.4 ! 2012-09 (G. Reffray and J. Chanut) New inputs + mods 11 11 !!---------------------------------------------------------------------- 12 12 #if defined key_bdy … … 15 15 !!---------------------------------------------------------------------- 16 16 !! PUBLIC 17 !! tide_init : read of namelist and initialisation of tidal harmonics data17 !! bdytide_init : read of namelist and initialisation of tidal harmonics data 18 18 !! 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 19 !!---------------------------------------------------------------------- 26 20 USE timing ! Timing … … 34 28 USE bdy_oce ! ocean open boundary conditions 35 29 USE daymod ! calendar 30 USE wrk_nemo ! Memory allocation 31 USE tideini 32 ! USE tide_mod ! Useless ?? 36 33 USE fldread, ONLY: fld_map 37 34 … … 39 36 PRIVATE 40 37 41 PUBLIC tide_init ! routine called in nemo_init42 PUBLIC tide_update ! routine called in bdydyn38 PUBLIC bdytide_init ! routine called in bdy_init 39 PUBLIC bdytide_update ! routine called in bdy_dta 43 40 44 41 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 42 REAL(wp), POINTER, DIMENSION(:,:,:) :: ssh0 !: Tidal constituents : SSH0 (read in file) 43 REAL(wp), POINTER, DIMENSION(:,:,:) :: u0 !: Tidal constituents : U0 (read in file) 44 REAL(wp), POINTER, DIMENSION(:,:,:) :: v0 !: Tidal constituents : V0 (read in file) 45 REAL(wp), POINTER, DIMENSION(:,:,:) :: ssh !: Tidal constituents : SSH (after nodal cor.) 46 REAL(wp), POINTER, DIMENSION(:,:,:) :: u !: Tidal constituents : U (after nodal cor.) 47 REAL(wp), POINTER, DIMENSION(:,:,:) :: v !: Tidal constituents : V (after nodal cor.) 50 48 END TYPE TIDES_DATA 51 49 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 55 50 TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides !: External tidal harmonics data 51 56 52 !!---------------------------------------------------------------------- 57 53 !! NEMO/OPA 3.3 , NEMO Consortium (2010) … … 61 57 CONTAINS 62 58 63 SUBROUTINE tide_init64 !!---------------------------------------------------------------------- 65 !! *** SUBROUTINE tide_init ***59 SUBROUTINE bdytide_init 60 !!---------------------------------------------------------------------- 61 !! *** SUBROUTINE bdytide_init *** 66 62 !! 67 63 !! ** Purpose : - Read in namelist for tides and initialise external … … 71 67 !! namelist variables 72 68 !!------------------- 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) 69 CHARACTER(len=80) :: filtide !: Filename root for tidal input files 70 LOGICAL :: ln_bdytide_2ddta !: If true, read 2d harmonic data 71 LOGICAL :: ln_bdytide_conj !: If true, assume complex conjugate tidal data 77 72 !! 78 INTEGER , DIMENSION(jptides_max) :: nindx !: index of pre-set tidal components79 INTEGER :: i b_bdy, itide, ib!: dummy loop indices73 INTEGER :: ib_bdy, itide, ib !: dummy loop indices 74 INTEGER :: ii, ij !: dummy loop indices 80 75 INTEGER :: inum, igrd 81 INTEGER, DIMENSION(3) :: ilen0 82 CHARACTER(len=80) :: clfile !: full file name for tidal input file83 REAL(wp) :: z_arg, z_atde, z_btde, z1t, z2t84 REAL(wp), DIMENSION(jptides_max) :: z_vplu, z_ftc85 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: dta_read !: work space to read in tidal harmonics data76 INTEGER, DIMENSION(3) :: ilen0 !: length of boundary data (from OBC arrays) 77 INTEGER, POINTER, DIMENSION(:) :: nblen, nblenrim ! short cuts 78 CHARACTER(len=80) :: clfile !: full file name for tidal input file 79 REAL(wp),ALLOCATABLE, DIMENSION(:,:,:) :: dta_read !: work space to read in tidal harmonics data 80 REAL(wp), POINTER, DIMENSION(:,:) :: ztr, zti !: " " " " " " " " 86 81 !! 87 TYPE(TIDES_DATA), POINTER :: td !: local short cut82 TYPE(TIDES_DATA), POINTER :: td !: local short cut 88 83 !! 89 NAMELIST/nambdy_tide/ln_tide_date, filtide, tide_cpt, tide_speed 90 !!---------------------------------------------------------------------- 91 92 IF( nn_timing == 1 ) CALL timing_start('tide_init') 93 94 IF(lwp) WRITE(numout,*) 95 IF(lwp) WRITE(numout,*) 'tide_init : initialization of tidal harmonic forcing at open boundaries' 96 IF(lwp) WRITE(numout,*) '~~~~~~~~~' 84 NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj 85 !!---------------------------------------------------------------------- 86 87 IF( nn_timing == 1 ) CALL timing_start('bdytide_init') 88 89 IF (nb_bdy>0) THEN 90 IF(lwp) WRITE(numout,*) 91 IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries' 92 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 93 ENDIF 94 95 ln_bdytide_2ddta = .FALSE. 96 ln_bdytide_conj = .FALSE. 97 97 98 98 REWIND(numnam) … … 101 101 102 102 td => tides(ib_bdy) 103 nblen => idx_bdy(ib_bdy)%nblen 104 nblenrim => idx_bdy(ib_bdy)%nblenrim 103 105 104 106 ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries 105 ln_tide_date = .false.106 107 filtide(:) = '' 107 tide_speed(:) = 0.0108 tide_cpt(:) = ''109 108 110 109 ! Don't REWIND here - may need to read more than one of these namelists. 111 110 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 111 ! ! Parameter control and print 148 IF( td%ncpt < 1 ) THEN 149 CALL ctl_stop( ' Did not find any tidal components in namelist nambdy_tide' ) 112 IF(lwp) WRITE(numout,*) ' ' 113 IF(lwp) WRITE(numout,*) ' Namelist nambdy_tide : tidal harmonic forcing at open boundaries' 114 IF(lwp) WRITE(numout,*) ' read tidal data in 2d files: ', ln_bdytide_2ddta 115 IF(lwp) WRITE(numout,*) ' assume complex conjugate : ', ln_bdytide_conj 116 IF(lwp) WRITE(numout,*) ' Number of tidal components to read: ', nb_harmo 117 IF(lwp) THEN 118 WRITE(numout,*) ' Tidal cpt name - Phase speed (deg/hr)' 119 DO itide = 1, nb_harmo 120 WRITE(numout,*) ' ', Wave(ntide(itide))%cname_tide, omega_tide(itide) 121 END DO 122 ENDIF 123 IF(lwp) WRITE(numout,*) ' ' 124 125 ! Allocate space for tidal harmonics data - get size from OBC data arrays 126 ! ----------------------------------------------------------------------- 127 128 ! JC: If FRS scheme is used, we assume that tidal is needed over the whole 129 ! relaxation area 130 IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 131 ilen0(:)=nblen(:) 150 132 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) 133 ilen0(:)=nblenrim(:) 156 134 ENDIF 157 135 158 159 ! Allocate space for tidal harmonics data - 160 ! get size from OBC data arrays 161 ! --------------------------------------- 162 163 ilen0(1) = SIZE( dta_bdy(ib_bdy)%ssh ) 164 ALLOCATE( td%ssh( ilen0(1), td%ncpt, 2 ) ) 165 166 ilen0(2) = SIZE( dta_bdy(ib_bdy)%u2d ) 167 ALLOCATE( td%u( ilen0(2), td%ncpt, 2 ) ) 168 169 ilen0(3) = SIZE( dta_bdy(ib_bdy)%v2d ) 170 ALLOCATE( td%v( ilen0(3), td%ncpt, 2 ) ) 171 172 ALLOCATE( dta_read( MAXVAL(ilen0), 1, 1 ) ) 173 174 175 ! Open files and read in tidal forcing data 176 ! ----------------------------------------- 177 178 DO itide = 1, td%ncpt 179 ! ! SSH fields 180 clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc' 181 IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 182 CALL iom_open( clfile, inum ) 183 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) 185 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) 136 ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) ) 137 ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) ) 138 139 ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) ) 140 ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) ) 141 142 ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) ) 143 ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) ) 144 145 td%ssh0(:,:,:) = 0.e0 146 td%ssh(:,:,:) = 0.e0 147 td%u0(:,:,:) = 0.e0 148 td%u(:,:,:) = 0.e0 149 td%v0(:,:,:) = 0.e0 150 td%v(:,:,:) = 0.e0 151 152 IF (ln_bdytide_2ddta) THEN 153 ! It is assumed that each data file contains all complex harmonic amplitudes 154 ! given on the data domain (ie global, jpidta x jpjdta) 155 ! 156 CALL wrk_alloc( jpi, jpj, zti, ztr ) 157 ! 158 ! SSH fields 159 clfile = TRIM(filtide)//'_grid_T.nc' 160 CALL iom_open (clfile , inum ) 161 igrd = 1 ! Everything is at T-points here 162 DO itide = 1, nb_harmo 163 CALL iom_get ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) ) 164 CALL iom_get ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) ) 165 DO ib = 1, ilen0(igrd) 166 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 167 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 168 td%ssh0(ib,itide,1) = ztr(ii,ij) 169 td%ssh0(ib,itide,2) = zti(ii,ij) 170 END DO 171 END DO 187 172 CALL iom_close( inum ) 188 ! ! U fields 189 clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc' 190 IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 191 CALL iom_open( clfile, inum ) 192 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) 194 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) 173 ! 174 ! U fields 175 clfile = TRIM(filtide)//'_grid_U.nc' 176 CALL iom_open (clfile , inum ) 177 igrd = 2 ! Everything is at U-points here 178 DO itide = 1, nb_harmo 179 CALL iom_get ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) ) 180 CALL iom_get ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) ) 181 DO ib = 1, ilen0(igrd) 182 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 183 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 184 td%u0(ib,itide,1) = ztr(ii,ij) 185 td%u0(ib,itide,2) = zti(ii,ij) 186 END DO 187 END DO 196 188 CALL iom_close( inum ) 197 ! ! V fields 198 clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc' 199 IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 200 CALL iom_open( clfile, inum ) 201 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) 203 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) 189 ! 190 ! V fields 191 clfile = TRIM(filtide)//'_grid_V.nc' 192 CALL iom_open (clfile , inum ) 193 igrd = 3 ! Everything is at V-points here 194 DO itide = 1, nb_harmo 195 CALL iom_get ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) ) 196 CALL iom_get ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) ) 197 DO ib = 1, ilen0(igrd) 198 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 199 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 200 td%v0(ib,itide,1) = ztr(ii,ij) 201 td%v0(ib,itide,2) = zti(ii,ij) 202 END DO 203 END DO 205 204 CALL iom_close( inum ) 206 205 ! 207 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 206 CALL wrk_dealloc( jpi, jpj, ztr, zti ) 207 ! 208 ELSE 209 ! 210 ! Read tidal data only on bdy segments 211 ! 212 ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) ) 213 214 ! Open files and read in tidal forcing data 215 ! ----------------------------------------- 216 217 DO itide = 1, nb_harmo 218 ! ! SSH fields 219 clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc' 220 CALL iom_open( clfile, inum ) 221 CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 222 td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1) 223 CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) ) 224 td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1) 225 CALL iom_close( inum ) 226 ! ! U fields 227 clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc' 228 CALL iom_open( clfile, inum ) 229 CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 230 td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1) 231 CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) ) 232 td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1) 233 CALL iom_close( inum ) 234 ! ! V fields 235 clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc' 236 CALL iom_open( clfile, inum ) 237 CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 238 td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1) 239 CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) ) 240 td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1) 241 CALL iom_close( inum ) 225 242 ! 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 243 END DO ! end loop on tidal components 244 ! 245 DEALLOCATE( dta_read ) 246 ENDIF ! ln_bdytide_2ddta=.true. 247 ! 248 IF ( ln_bdytide_conj ) THEN ! assume complex conjugate in data files 249 td%ssh0(:,:,2) = - td%ssh0(:,:,2) 250 td%u0 (:,:,2) = - td%u0 (:,:,2) 251 td%v0 (:,:,2) = - td%v0 (:,:,2) 252 ENDIF 265 253 ! 266 254 ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2 … … 268 256 END DO ! loop on ib_bdy 269 257 270 IF( nn_timing == 1 ) CALL timing_stop('tide_init') 271 272 END SUBROUTINE tide_init 273 274 275 SUBROUTINE tide_update ( kt, idx, dta, td, jit, time_offset ) 276 !!---------------------------------------------------------------------- 277 !! *** SUBROUTINE tide_update *** 258 IF( nn_timing == 1 ) CALL timing_stop('bdytide_init') 259 260 END SUBROUTINE bdytide_init 261 262 SUBROUTINE bdytide_update ( kt, idx, dta, td, jit, time_offset ) 263 !!---------------------------------------------------------------------- 264 !! *** SUBROUTINE bdytide_update *** 278 265 !! 279 266 !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays. 280 267 !! 281 268 !!---------------------------------------------------------------------- 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. 269 INTEGER, INTENT( in ) :: kt ! Main timestep counter 270 TYPE(OBC_INDEX), INTENT( in ) :: idx ! OBC indices 271 TYPE(OBC_DATA), INTENT(inout) :: dta ! OBC external data 272 TYPE(TIDES_DATA),INTENT( inout ) :: td ! tidal harmonics data 273 INTEGER,INTENT(in),OPTIONAL :: jit ! Barotropic timestep counter (for timesplitting option) 274 INTEGER,INTENT( in ), OPTIONAL :: time_offset ! time offset in units of timesteps. NB. if jit 275 ! is present then units = subcycle timesteps. 276 ! time_offset = 0 => get data at "now" time level 277 ! time_offset = -1 => get data at "before" time level 278 ! time_offset = +1 => get data at "after" time level 279 ! etc. 294 280 !! 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 299 !!---------------------------------------------------------------------- 300 301 IF( nn_timing == 1 ) CALL timing_start('tide_update') 281 INTEGER, DIMENSION(3) :: ilen0 !: length of boundary data (from OBC arrays) 282 INTEGER :: itide, igrd, ib ! dummy loop indices 283 INTEGER :: time_add ! time offset in units of timesteps 284 REAL(wp) :: z_arg, z_sarg, zflag, zramp 285 REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost 286 !!---------------------------------------------------------------------- 287 288 IF( nn_timing == 1 ) CALL timing_start('bdytide_update') 289 290 ilen0(1) = SIZE(td%ssh(:,1,1)) 291 ilen0(2) = SIZE(td%u(:,1,1)) 292 ilen0(3) = SIZE(td%v(:,1,1)) 293 294 zflag=1 295 IF ( PRESENT(jit) ) THEN 296 IF ( jit /= 1 ) zflag=0 297 ENDIF 298 299 IF ( nsec_day == NINT(0.5 * rdttra(1)) .AND. zflag==1 ) THEN 300 ! 301 kt_tide = kt 302 ! 303 IF(lwp) THEN 304 WRITE(numout,*) 305 WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt 306 WRITE(numout,*) '~~~~~~~~~~~~~~ ' 307 ENDIF 308 ! 309 CALL tide_init_elevation ( idx, td ) 310 CALL tide_init_velocities( idx, td ) 311 ! 312 ENDIF 302 313 303 314 time_add = 0 … … 306 317 ENDIF 307 318 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 319 IF( PRESENT(jit) ) THEN 311 z_arg = ( ( kt-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) ) * rad / 3600.0320 z_arg = ( ((kt-kt_tide)-1) * rdt + (jit+time_add) * rdt / REAL(nn_baro,wp) ) 312 321 ELSE 313 z_arg = ( kt+time_add) * rdt * rad / 3600.0322 z_arg = ((kt-kt_tide)+time_add) * rdt 314 323 ENDIF 315 324 316 DO itide = 1, td%ncpt 317 z_sarg = z_arg * td%speed(itide) 325 ! Linear ramp on tidal component at open boundaries 326 zramp = 1. 327 IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0.),1.) 328 329 DO itide = 1, nb_harmo 330 z_sarg = z_arg * omega_tide(itide) 318 331 z_cost(itide) = COS( z_sarg ) 319 332 z_sist(itide) = SIN( z_sarg ) 320 333 END DO 321 334 322 DO itide = 1, td%ncpt 323 igrd=1 ! SSH on tracer grid. 324 DO ib = 1, idx%nblenrim(igrd) 325 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) 335 DO itide = 1, nb_harmo 336 igrd=1 ! SSH on tracer grid 337 DO ib = 1, ilen0(igrd) 338 dta%ssh(ib) = dta%ssh(ib) + zramp*(td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide)) 327 339 END DO 328 340 igrd=2 ! U grid 329 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) 341 DO ib = 1, ilen0(igrd) 342 dta%u2d(ib) = dta%u2d(ib) + zramp*(td%u (ib,itide,1)*z_cost(itide) + td%u (ib,itide,2)*z_sist(itide)) 332 343 END DO 333 344 igrd=3 ! V grid 334 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) 345 DO ib = 1, ilen0(igrd) 346 dta%v2d(ib) = dta%v2d(ib) + zramp*(td%v (ib,itide,1)*z_cost(itide) + td%v (ib,itide,2)*z_sist(itide)) 337 347 END DO 338 348 END DO 339 349 ! 340 IF( nn_timing == 1 ) CALL timing_stop(' tide_update')350 IF( nn_timing == 1 ) CALL timing_stop('bdytide_update') 341 351 ! 342 END SUBROUTINE tide_update 343 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) 352 END SUBROUTINE bdytide_update 353 354 SUBROUTINE tide_init_elevation( idx, td ) 355 !!---------------------------------------------------------------------- 356 !! *** ROUTINE tide_init_elevation *** 357 !!---------------------------------------------------------------------- 358 TYPE(OBC_INDEX), INTENT( in ) :: idx ! OBC indices 359 TYPE(TIDES_DATA),INTENT( inout ) :: td ! tidal harmonics data 360 !! * Local declarations 361 INTEGER, DIMENSION(1) :: ilen0 !: length of boundary data (from OBC arrays) 362 REAL(wp),ALLOCATABLE, DIMENSION(:) :: mod_tide, phi_tide 363 INTEGER :: itide, igrd, ib ! dummy loop indices 364 365 igrd=1 366 ! SSH on tracer grid. 367 368 ilen0(1) = SIZE(td%ssh0(:,1,1)) 369 370 ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd))) 371 372 DO itide = 1, nb_harmo 373 DO ib = 1, ilen0(igrd) 374 mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.) 375 phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1)) 376 END DO 377 DO ib = 1 , ilen0(igrd) 378 mod_tide(ib)=mod_tide(ib)*ftide(itide) 379 phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 380 ENDDO 381 DO ib = 1 , ilen0(igrd) 382 td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 383 td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 384 ENDDO 396 385 END DO 397 386 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 387 DEALLOCATE(mod_tide,phi_tide) 388 389 END SUBROUTINE tide_init_elevation 390 391 SUBROUTINE tide_init_velocities( idx, td ) 392 !!---------------------------------------------------------------------- 393 !! *** ROUTINE tide_init_elevation *** 394 !!---------------------------------------------------------------------- 395 TYPE(OBC_INDEX), INTENT( in ) :: idx ! OBC indices 396 TYPE(TIDES_DATA),INTENT( inout ) :: td ! tidal harmonics data 397 !! * Local declarations 398 INTEGER, DIMENSION(3) :: ilen0 !: length of boundary data (from OBC arrays) 399 REAL(wp),ALLOCATABLE, DIMENSION(:) :: mod_tide, phi_tide 400 INTEGER :: itide, igrd, ib ! dummy loop indices 401 402 ilen0(2) = SIZE(td%u0(:,1,1)) 403 ilen0(3) = SIZE(td%v0(:,1,1)) 404 405 igrd=2 ! U grid. 406 407 ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd))) 408 409 DO itide = 1, nb_harmo 410 DO ib = 1, ilen0(igrd) 411 mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.) 412 phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1)) 413 END DO 414 DO ib = 1, ilen0(igrd) 415 mod_tide(ib)=mod_tide(ib)*ftide(itide) 416 phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 417 ENDDO 418 DO ib = 1, ilen0(igrd) 419 td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 420 td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 421 ENDDO 422 END DO 423 424 DEALLOCATE(mod_tide,phi_tide) 425 426 igrd=3 ! V grid. 427 428 ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd))) 429 430 DO itide = 1, nb_harmo 431 DO ib = 1, ilen0(igrd) 432 mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.) 433 phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1)) 434 END DO 435 DO ib = 1, ilen0(igrd) 436 mod_tide(ib)=mod_tide(ib)*ftide(itide) 437 phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide) 438 ENDDO 439 DO ib = 1, ilen0(igrd) 440 td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib)) 441 td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib)) 442 ENDDO 443 END DO 444 445 DEALLOCATE(mod_tide,phi_tide) 446 447 END SUBROUTINE tide_init_velocities 717 448 #else 718 449 !!---------------------------------------------------------------------- 719 450 !! Dummy module NO Unstruct Open Boundary Conditions for tides 720 451 !!---------------------------------------------------------------------- 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 452 CONTAINS 726 SUBROUTINE tide_init ! Empty routine 727 END SUBROUTINE tide_init 728 SUBROUTINE tide_data ! Empty routine 729 END SUBROUTINE tide_data 730 SUBROUTINE tide_update( kt, kit ) ! Empty routine 731 WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, kit 732 END SUBROUTINE tide_update 453 SUBROUTINE bdytide_init ! Empty routine 454 WRITE(*,*) 'bdytide_init: You should not have seen this print! error?' 455 END SUBROUTINE bdytide_init 456 SUBROUTINE bdytide_update( kt, jit ) ! Empty routine 457 WRITE(*,*) 'bdytide_update: You should not have seen this print! error?', kt, jit 458 END SUBROUTINE bdytide_update 733 459 #endif 734 460
Note: See TracChangeset
for help on using the changeset viewer.