Changeset 1125 for trunk/NEMO/OPA_SRC/BDY/bdytides.F90
- Timestamp:
- 2008-06-23T11:05:02+02:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/BDY/bdytides.F90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/BDY/bdytides.F90
r911 r1125 1 1 MODULE bdytides 2 !!====================================================================== ===========2 !!====================================================================== 3 3 !! *** MODULE bdytides *** 4 4 !! Ocean dynamics: Tidal forcing at open boundaries 5 !!================================================================================= 6 #if defined key_bdy_tides 7 !!--------------------------------------------------------------------------------- 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 !!---------------------------------------------------------------------- 10 #if defined key_bdy 11 !!---------------------------------------------------------------------- 12 !! 'key_bdy' Unstructured Open Boundary Condition 13 !!---------------------------------------------------------------------- 8 14 !! PUBLIC 9 !! tide_init: read of namelist10 !! tide_data: read in and initialisation of tidal constituents at boundary11 !! tide_update: calculation of tidal forcing at each timestep15 !! tide_init : read of namelist 16 !! tide_data : read in and initialisation of tidal constituents at boundary 17 !! tide_update : calculation of tidal forcing at each timestep 12 18 !! PRIVATE 13 !! uvset :\ 14 !! vday : | Routines to correct tidal harmonics forcing for 15 !! shpen : | start time of integration 16 !! ufset : | 17 !! vset :/ 18 !!--------------------------------------------------------------------------------- 19 20 !!--------------------------------------------------------------------------------- 21 !! * Modules used 19 !! uvset :\ 20 !! vday : | Routines to correct tidal harmonics forcing for 21 !! shpen : | start time of integration 22 !! ufset : | 23 !! vset :/ 24 !!---------------------------------------------------------------------- 22 25 USE oce ! ocean dynamics and tracers 23 26 USE dom_oce ! ocean space and time domain … … 29 32 USE bdy_oce ! ocean open boundary conditions 30 33 USE daymod ! calendar 34 31 35 IMPLICIT NONE 32 36 PRIVATE 33 37 34 !! * Accessibility 35 PUBLIC tide_init ! routine called in bdyini 36 PUBLIC tide_data ! routine called in bdyini 37 PUBLIC tide_update ! routine called in bdydyn 38 39 !! * Module variables 40 LOGICAL, PUBLIC, PARAMETER :: lk_bdy_tides = .TRUE. !: tidal forcing at boundaries. 41 42 CHARACTER(len=80), PUBLIC :: & 43 filtide !: Filename root for tidal input files 44 45 INTEGER, PUBLIC, PARAMETER :: ntide_max = 15 ! Max number of tidal contituents 46 INTEGER, PUBLIC :: ntide ! Actual number of tidal constituents 47 48 CHARACTER(len=4), PUBLIC, DIMENSION(ntide_max) :: & 49 tide_cpt !: Names of tidal components used. 50 51 logical, PUBLIC :: ln_tide_date ! if true correct tide phases and amplitude for model start date 52 53 REAL(wp), PUBLIC, DIMENSION(ntide_max) :: tide_speed ! Phase speed of tidal constituent (deg/hr) 54 INTEGER, DIMENSION(ntide_max) :: indx 55 REAL(wp), DIMENSION(jpbdim,ntide_max) :: & 56 ssh1, ssh2, & !: Tidal constituents : SSH 57 u1 , u2 , & !: Tidal constituents : U 58 v1 , v2 !: Tidal constituents : V 38 PUBLIC tide_init ! routine called in bdyini 39 PUBLIC tide_data ! routine called in bdyini 40 PUBLIC tide_update ! routine called in bdydyn 41 42 LOGICAL, PUBLIC :: ln_tide_date !: =T correct tide phases and amplitude for model start date 43 44 INTEGER, PARAMETER :: jptides_max = 15 !: Max number of tidal contituents 45 INTEGER :: ntide !: Actual number of tidal constituents 46 47 CHARACTER(len=80), PUBLIC :: filtide !: Filename root for tidal input files 48 CHARACTER(len= 4), PUBLIC, DIMENSION(jptides_max) :: tide_cpt !: Names of tidal components used. 49 50 INTEGER , DIMENSION(jptides_max) :: nindx !: ??? 51 REAL(wp), DIMENSION(jptides_max) :: tide_speed !: Phase speed of tidal constituent (deg/hr) 59 52 53 REAL(wp), DIMENSION(jpbdim,jptides_max) :: ssh1, ssh2 !: Tidal constituents : SSH 54 REAL(wp), DIMENSION(jpbdim,jptides_max) :: u1 , u2 !: Tidal constituents : U 55 REAL(wp), DIMENSION(jpbdim,jptides_max) :: v1 , v2 !: Tidal constituents : V 60 56 61 !!--------------------------------------------------------------------------------- 57 !!---------------------------------------------------------------------- 58 !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 59 !! $Id: $ 60 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 61 !!---------------------------------------------------------------------- 62 62 63 63 CONTAINS 64 64 65 65 SUBROUTINE tide_init 66 !!---------------------------------------------------------------------- --------67 !! SUBROUTINE tide_init68 !! ***********************66 !!---------------------------------------------------------------------- 67 !! *** SUBROUTINE tide_init *** 68 !! 69 69 !! ** Purpose : - Read in namelist for tides 70 70 !! 71 !! History : 72 !! NEMO v2.0 ! 07-01 (D.Storkey) Original 73 !!------------------------------------------------------------------------------ 74 !! * Local declarations 75 INTEGER :: jtide ! dummy loop index 76 ! different from nblendta!) 77 78 !!------------------------------------------------------------------------------ 79 !! OPA 9.0, LODYC-IPSL (2007) 80 !!------------------------------------------------------------------------------ 81 82 NAMELIST/namtide/filtide, tide_cpt, tide_speed, ln_tide_date 83 71 !!---------------------------------------------------------------------- 72 INTEGER :: itide ! dummy loop index 73 !! 74 NAMELIST/namtide/ln_tide_date, filtide, tide_cpt, tide_speed 84 75 !!---------------------------------------------------------------------- 85 76 … … 88 79 IF(lwp) WRITE(numout,*) '~~~~~~~~~' 89 80 90 ! 0. Read namelist parameters 91 ! --------------------------- 92 93 do jtide = 1, ntide_max 94 tide_cpt(jtide) = '' 95 enddo 96 97 REWIND( numnam ) 81 ! Namelist namtide : tidal harmonic forcing at open boundaries 82 ln_tide_date = .false. 83 filtide(:) = '' 84 tide_speed(:) = 0.0 85 tide_cpt(:) = '' 86 REWIND( numnam ) ! Read namelist parameters 98 87 READ ( numnam, namtide ) 99 100 ! control prints 101 IF(lwp) WRITE(numout,*) ' namtide' 102 103 ! Count number of components specified: 104 ntide=ntide_max 105 do jtide = 1, ntide_max 106 if ( tide_cpt(jtide) == '' ) then 107 ntide = jtide-1 108 exit 109 endif 110 111 enddo 112 ! find constients in standard list 113 do jtide = 1, ntide 114 indx(jtide) = 0 115 if ( TRIM(tide_cpt(jtide)) == 'Q1' ) indx(jtide) = 1 116 if ( TRIM(tide_cpt(jtide)) == 'O1' ) indx(jtide) = 2 117 if ( TRIM(tide_cpt(jtide)) == 'P1' ) indx(jtide) = 3 118 if ( TRIM(tide_cpt(jtide)) == 'S1' ) indx(jtide) = 4 119 if ( TRIM(tide_cpt(jtide)) == 'K1' ) indx(jtide) = 5 120 if ( TRIM(tide_cpt(jtide)) == '2N2' ) indx(jtide) = 6 121 if ( TRIM(tide_cpt(jtide)) == 'MU2' ) indx(jtide) = 7 122 if ( TRIM(tide_cpt(jtide)) == 'N2' ) indx(jtide) = 8 123 if ( TRIM(tide_cpt(jtide)) == 'NU2' ) indx(jtide) = 9 124 if ( TRIM(tide_cpt(jtide)) == 'M2' ) indx(jtide) = 10 125 if ( TRIM(tide_cpt(jtide)) == 'L2' ) indx(jtide) = 11 126 if ( TRIM(tide_cpt(jtide)) == 'T2' ) indx(jtide) = 12 127 if ( TRIM(tide_cpt(jtide)) == 'S2' ) indx(jtide) = 13 128 if ( TRIM(tide_cpt(jtide)) == 'K2' ) indx(jtide) = 14 129 if ( TRIM(tide_cpt(jtide)) == 'M4' ) indx(jtide) = 15 130 if (indx(jtide) == 0 ) then 131 if(lwp) write(numout,*) 'WARNING: constitunent', jtide,':', tide_cpt(jtide) & 132 , 'not in standard list' 133 endif 134 enddo 135 ! 136 if ( ntide < 1 ) then 137 if(lwp) write(numout,*) 138 if(lwp) write(numout,*) ' E R R O R : Did not find any tidal components in namelist.' 139 if(lwp) write(numout,*) ' ========== ' 140 if(lwp) write(numout,*) 141 nstop = nstop + 1 142 else 143 if(lwp) write(numout,*) 144 if(lwp) write(numout,*) ntide,' tidal components specified : ' 145 if(lwp) write(numout,*) tide_cpt(1:ntide) 146 if(lwp) write(numout,*) ntide,' phase speeds (deg/hr) : ' 147 if(lwp) write(numout,*) tide_speed(1:ntide) 148 endif 88 ! ! Count number of components specified 89 ntide = jptides_max 90 itide = 1 91 DO WHILE( tide_cpt(itide) /= '' ) 92 ntide = itide 93 itide = itide + 1 94 END DO 95 ! ! find constituents in standard list 96 DO itide = 1, ntide 97 nindx(itide) = 0 98 IF( TRIM( tide_cpt(itide) ) == 'Q1' ) nindx(itide) = 1 99 IF( TRIM( tide_cpt(itide) ) == 'O1' ) nindx(itide) = 2 100 IF( TRIM( tide_cpt(itide) ) == 'P1' ) nindx(itide) = 3 101 IF( TRIM( tide_cpt(itide) ) == 'S1' ) nindx(itide) = 4 102 IF( TRIM( tide_cpt(itide) ) == 'K1' ) nindx(itide) = 5 103 IF( TRIM( tide_cpt(itide) ) == '2N2' ) nindx(itide) = 6 104 IF( TRIM( tide_cpt(itide) ) == 'MU2' ) nindx(itide) = 7 105 IF( TRIM( tide_cpt(itide) ) == 'N2' ) nindx(itide) = 8 106 IF( TRIM( tide_cpt(itide) ) == 'NU2' ) nindx(itide) = 9 107 IF( TRIM( tide_cpt(itide) ) == 'M2' ) nindx(itide) = 10 108 IF( TRIM( tide_cpt(itide) ) == 'L2' ) nindx(itide) = 11 109 IF( TRIM( tide_cpt(itide) ) == 'T2' ) nindx(itide) = 12 110 IF( TRIM( tide_cpt(itide) ) == 'S2' ) nindx(itide) = 13 111 IF( TRIM( tide_cpt(itide) ) == 'K2' ) nindx(itide) = 14 112 IF( TRIM( tide_cpt(itide) ) == 'M4' ) nindx(itide) = 15 113 IF( nindx(itide) == 0 .AND. lwp ) THEN 114 WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list' 115 CALL ctl_warn( ctmp1 ) 116 ENDIF 117 END DO 118 ! ! Parameter control and print 119 IF( ntide < 1 ) THEN 120 CALL ctl_stop( ' Did not find any tidal components in namelist namtide' ) 121 ELSE 122 IF(lwp) WRITE(numout,*) ' Namelist namtide : tidal harmonic forcing at open boundaries' 123 IF(lwp) WRITE(numout,*) ' tidal components specified ', ntide 124 IF(lwp) WRITE(numout,*) ' ', tide_cpt(1:ntide) 125 IF(lwp) WRITE(numout,*) ' associated phase speeds (deg/hr) : ' 126 IF(lwp) WRITE(numout,*) ' ', tide_speed(1:ntide) 127 ENDIF 149 128 150 129 ! Initialisation of tidal harmonics arrays 151 152 130 sshtide(:) = 0.e0 153 utide (:) = 0.e0154 vtide (:) = 0.e0155 131 utide (:) = 0.e0 132 vtide (:) = 0.e0 133 ! 156 134 END SUBROUTINE tide_init 157 135 136 158 137 SUBROUTINE tide_data 138 !!---------------------------------------------------------------------- 139 !! *** SUBROUTINE tide_data *** 140 !! 141 !! ** Purpose : - Read in tidal harmonics data and adjust for the start 142 !! time of the model run. 143 !! 144 !!---------------------------------------------------------------------- 145 INTEGER :: itide, igrd, ib ! dummy loop indices 146 CHARACTER(len=80) :: clfile ! full file name for tidal input file 147 INTEGER :: ipi, ipj, inum, idvar ! temporary integers (netcdf read) 148 INTEGER, DIMENSION(3) :: lendta ! length of data in the file (note may be different from nblendta!) 149 REAL(wp) :: z_arg, z_atde, z_btde, z1t, z2t 150 REAL(wp), DIMENSION(jpbdta,1) :: zdta ! temporary array for data fields 151 REAL(wp), DIMENSION(jptides_max) :: z_vplu, z_ftc 159 152 !!------------------------------------------------------------------------------ 160 !! SUBROUTINE tide_data 161 !! *********************** 162 !! ** Purpose : - Read in tidal harmonics data and adjust for the start time of 163 !! the model run. 164 !! 165 !! History : 166 !! NEMO v2.0 ! 07-01 (D.Storkey) Original 167 !! ! 08-01 (J.Holt) Add date correction. 168 !!------------------------------------------------------------------------------ 169 !! * Local declarations 170 CHARACTER(len=80) :: tidefile ! full file name for tidal input file 171 INTEGER :: jtide, jgrd, jb ! dummy loop indices 172 INTEGER :: ipi, ipj, inum, idvar ! temporary integers (netcdf read) 173 INTEGER, DIMENSION(3) :: lendta ! length of data in the file (note may be 174 ! different from nblendta!) 175 INTEGER :: iday,imonth,iyear 176 REAL(wp) :: arg, atde, btde,z1t,z2t 177 REAL(wp), DIMENSION(jpbdta,1) :: & 178 pdta ! temporary array for data fields 179 REAL(WP), DIMENSION(ntide_max) :: vplu, ftc 180 181 !!------------------------------------------------------------------------------ 182 !! OPA 9.0, LODYC-IPSL (2007) 183 !!------------------------------------------------------------------------------ 184 185 ! 1. Open files and read in tidal forcing data 186 ! -------------------------------------------- 153 154 ! Open files and read in tidal forcing data 155 ! ----------------------------------------- 187 156 188 157 ipj = 1 189 158 190 do jtide = 1 ,ntide 191 192 ! SSH fields 193 ! ---------- 194 195 tidefile = TRIM(filtide)//TRIM(tide_cpt(jtide))//'_grid_T.nc' 196 if(lwp) write(numout,*) 'Reading data from file ',tidefile 197 CALL iom_open( tidefile, inum ) 198 199 jgrd = 1 200 IF ( nblendta(jgrd) .le. 0 ) THEN 201 idvar = iom_varid( inum,'z1' ) 202 if(lwp) write(numout,*) 'iom_file(1)%ndims(idvar) : ',iom_file%ndims(idvar) 203 nblendta(jgrd) = iom_file(inum)%dimsz(1,idvar) 204 WRITE(numout,*) 'Dim size for z1 is ',nblendta(jgrd) 159 DO itide = 1, ntide 160 ! ! SSH fields 161 clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc' 162 IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 163 CALL iom_open( clfile, inum ) 164 igrd = 1 165 IF( nblendta(igrd) <= 0 ) THEN 166 idvar = iom_varid( inum,'z1' ) 167 IF(lwp) WRITE(numout,*) 'iom_file(1)%ndims(idvar) : ',iom_file%ndims(idvar) 168 nblendta(igrd) = iom_file(inum)%dimsz(1,idvar) 169 WRITE(numout,*) 'Dim size for z1 is ', nblendta(igrd) 205 170 ENDIF 206 ipi=nblendta(jgrd) 207 208 CALL iom_get ( inum, jpdom_unknown, 'z1', pdta(1:ipi,1:ipj) ) 209 DO jb=1, nblenrim(jgrd) 210 ssh1(jb,jtide) = pdta(nbmap(jb,jgrd),1) 211 END DO 212 213 CALL iom_get ( inum, jpdom_unknown, 'z2', pdta(1:ipi,1:ipj) ) 214 DO jb=1, nblenrim(jgrd) 215 ssh2(jb,jtide) = pdta(nbmap(jb,jgrd),1) 216 END DO 217 171 ipi = nblendta(igrd) 172 CALL iom_get( inum, jpdom_unknown, 'z1', zdta(1:ipi,1:ipj) ) 173 DO ib = 1, nblenrim(igrd) 174 ssh1(ib,itide) = zdta(nbmap(ib,igrd),1) 175 END DO 176 CALL iom_get( inum, jpdom_unknown, 'z2', zdta(1:ipi,1:ipj) ) 177 DO ib = 1, nblenrim(igrd) 178 ssh2(ib,itide) = zdta(nbmap(ib,igrd),1) 179 END DO 218 180 CALL iom_close( inum ) 219 220 ! U fields 221 ! -------- 222 223 tidefile = TRIM(filtide)//TRIM(tide_cpt(jtide))//'_grid_U.nc' 224 if(lwp) write(numout,*) 'Reading data from file ',tidefile 225 CALL iom_open( tidefile, inum ) 226 227 jgrd = 2 228 IF ( lendta(jgrd) .le. 0 ) THEN 229 idvar = iom_varid( inum,'u1' ) 230 lendta(jgrd) = iom_file(inum)%dimsz(1,idvar) 231 WRITE(numout,*) 'Dim size for u1 is ',lendta(jgrd) 181 ! 182 ! ! U fields 183 clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc' 184 IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile 185 CALL iom_open( clfile, inum ) 186 igrd = 2 187 IF( lendta(igrd) <= 0 ) THEN 188 idvar = iom_varid( inum,'u1' ) 189 lendta(igrd) = iom_file(inum)%dimsz(1,idvar) 190 WRITE(numout,*) 'Dim size for u1 is ',lendta(igrd) 232 191 ENDIF 233 ipi=lendta(jgrd) 234 235 CALL iom_get ( inum, jpdom_unknown, 'u1', pdta(1:ipi,1:ipj) ) 236 DO jb=1, nblenrim(jgrd) 237 u1(jb,jtide) = pdta(nbmap(jb,jgrd),1) 238 END DO 239 240 CALL iom_get ( inum, jpdom_unknown, 'u2', pdta(1:ipi,1:ipj) ) 241 DO jb=1, nblenrim(jgrd) 242 u2(jb,jtide) = pdta(nbmap(jb,jgrd),1) 243 END DO 244 192 ipi = lendta(igrd) 193 CALL iom_get( inum, jpdom_unknown, 'u1', zdta(1:ipi,1:ipj) ) 194 DO ib = 1, nblenrim(igrd) 195 u1(ib,itide) = zdta(nbmap(ib,igrd),1) 196 END DO 197 CALL iom_get( inum, jpdom_unknown, 'u2', zdta(1:ipi,1:ipj) ) 198 DO ib = 1, nblenrim(igrd) 199 u2(ib,itide) = zdta(nbmap(ib,igrd),1) 200 END DO 245 201 CALL iom_close( inum ) 246 247 ! V fields 248 ! -------- 249 250 tidefile = TRIM(filtide)//TRIM(tide_cpt(jtide))//'_grid_V.nc' 251 if(lwp) write(numout,*) 'Reading data from file ',tidefile 252 CALL iom_open( tidefile, inum ) 253 254 jgrd = 3 255 IF ( lendta(jgrd) .le. 0 ) THEN 256 idvar = iom_varid( inum,'v1' ) 257 lendta(jgrd) = iom_file(inum)%dimsz(1,idvar) 258 WRITE(numout,*) 'Dim size for v1 is ',lendta(jgrd) 202 ! 203 ! ! V fields 204 clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc' 205 if(lwp) write(numout,*) 'Reading data from file ', clfile 206 CALL iom_open( clfile, inum ) 207 igrd = 3 208 IF( lendta(igrd) <= 0 ) THEN 209 idvar = iom_varid( inum,'v1' ) 210 lendta(igrd) = iom_file(inum)%dimsz(1,idvar) 211 WRITE(numout,*) 'Dim size for v1 is ', lendta(igrd) 259 212 ENDIF 260 ipi=lendta(jgrd) 261 262 CALL iom_get ( inum, jpdom_unknown, 'v1', pdta(1:ipi,1:ipj) ) 263 DO jb=1, nblenrim(jgrd) 264 v1(jb,jtide) = pdta(nbmap(jb,jgrd),1) 265 END DO 266 267 CALL iom_get ( inum, jpdom_unknown, 'v2', pdta(1:ipi,1:ipj) ) 268 DO jb=1, nblenrim(jgrd) 269 v2(jb,jtide) = pdta(nbmap(jb,jgrd),1) 270 END DO 271 213 ipi = lendta(igrd) 214 CALL iom_get( inum, jpdom_unknown, 'v1', zdta(1:ipi,1:ipj) ) 215 DO ib = 1, nblenrim(igrd) 216 v1(ib,itide) = zdta(nbmap(ib,igrd),1) 217 END DO 218 CALL iom_get( inum, jpdom_unknown, 'v2', zdta(1:ipi,1:ipj) ) 219 DO ib=1, nblenrim(igrd) 220 v2(ib,itide) = zdta(nbmap(ib,igrd),1) 221 END DO 272 222 CALL iom_close( inum ) 273 enddo ! end loop on tidal components 274 ! 275 ! correct for date factors 276 ! 277 if( ln_tide_date ) then 278 279 280 ! Calculate date corrects for 15 standard consituents 281 iyear = int(ndate0 / 10000 ) ! initial year 282 imonth = int((ndate0 - iyear * 10000 ) / 100 ) ! initial month 283 iday = int(ndate0 - iyear * 10000 - imonth * 100) ! initial day betwen 1 and 30 284 285 call uvset(0,iday,imonth,iyear,ftc,vplu) 286 ! 287 if(lwp) write(numout,*) 'Correcting tide for date:',iday,imonth,iyear 288 289 do jtide = 1, ntide 290 ! 291 if(indx(jtide) .ne. 0) then 292 arg=3.14159265d0*vplu(indx(jtide))/180.0d0 293 atde=ftc(indx(jtide))*cos(arg) 294 btde=ftc(indx(jtide))*sin(arg) 295 if(lwp) then 296 write(numout,'(2i5,8f10.6)') jtide,indx(jtide),tide_speed(jtide), & 297 ftc(indx(jtide)),vplu(indx(jtide)) 298 endif 299 else 300 atde = 1.0_wp 301 btde = 0.0_wp 302 endif 303 ! elevation 304 jgrd = 1 305 do jb=1, nblenrim(jgrd) 306 z1t=atde*ssh1(jb,jtide)+btde*ssh2(jb,jtide) 307 z2t=atde*ssh2(jb,jtide)-btde*ssh1(jb,jtide) 308 ssh1(jb,jtide) = z1t 309 ssh2(jb,jtide) = z2t 310 end do 311 ! u 312 jgrd = 2 313 do jb=1, nblenrim(jgrd) 314 z1t=atde*u1(jb,jtide)+btde*u2(jb,jtide) 315 z2t=atde*u2(jb,jtide)-btde*u1(jb,jtide) 316 u1(jb,jtide) = z1t 317 u2(jb,jtide) = z2t 318 end do 319 ! v 320 jgrd = 3 321 do jb=1, nblenrim(jgrd) 322 z1t=atde*v1(jb,jtide)+btde*v2(jb,jtide) 323 z2t=atde*v2(jb,jtide)-btde*v1(jb,jtide) 324 v1(jb,jtide) = z1t 325 v2(jb,jtide) = z2t 326 end do 327 328 enddo ! end loop on tidal components 329 endif ! date correction 330 223 ! 224 END DO ! end loop on tidal components 225 226 IF( ln_tide_date ) THEN ! correct for date factors 227 228 !! used nmonth, nyear and nday from daymod.... 229 ! Calculate date corrects for 15 standard consituents 230 ! This is the initialisation step, so nday, nmonth, nyear are the 231 ! initial date/time of the integration. 232 233 CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu ) 234 235 IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear 236 237 DO itide = 1, ntide ! loop on tidal components 238 ! 239 IF( nindx(itide) /= 0 ) THEN 240 !!gm use rpi and rad global variable 241 z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0 242 z_atde=z_ftc(nindx(itide))*cos(z_arg) 243 z_btde=z_ftc(nindx(itide))*sin(z_arg) 244 IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), tide_speed(itide), & 245 & z_ftc(nindx(itide)), z_vplu(nindx(itide)) 246 ELSE 247 z_atde = 1.0_wp 248 z_btde = 0.0_wp 249 ENDIF 250 ! ! elevation 251 igrd = 1 252 DO ib = 1, nblenrim(igrd) 253 z1t = z_atde * ssh1(ib,itide) + z_btde * ssh2(ib,itide) 254 z2t = z_atde * ssh2(ib,itide) - z_btde * ssh1(ib,itide) 255 ssh1(ib,itide) = z1t 256 ssh2(ib,itide) = z2t 257 END DO 258 ! ! u 259 igrd = 2 260 DO ib = 1, nblenrim(igrd) 261 z1t = z_atde * u1(ib,itide) + z_btde * u2(ib,itide) 262 z2t = z_atde * u2(ib,itide) - z_btde * u1(ib,itide) 263 u1(ib,itide) = z1t 264 u2(ib,itide) = z2t 265 END DO 266 ! ! v 267 igrd = 3 268 DO ib = 1, nblenrim(igrd) 269 z1t = z_atde * v1(ib,itide) + z_btde * v2(ib,itide) 270 z2t = z_atde * v2(ib,itide) - z_btde * v1(ib,itide) 271 v1(ib,itide) = z1t 272 v2(ib,itide) = z2t 273 END DO 274 ! 275 END DO ! end loop on tidal components 276 ! 277 ENDIF ! date correction 278 ! 331 279 END SUBROUTINE tide_data 332 280 281 333 282 SUBROUTINE tide_update ( kt, jit ) 334 !!---------------------------------------------------------------------- --------335 !! SUBROUTINE tide_update336 !! ************************283 !!---------------------------------------------------------------------- 284 !! *** SUBROUTINE tide_update *** 285 !! 337 286 !! ** Purpose : - Add tidal forcing to sshbdy, ubtbdy and vbtbdy arrays. 338 287 !! 339 !! 340 !! History : 341 !! NEMO v2.0 ! 06-12 (D.Storkey) Original 342 !!------------------------------------------------------------------------------ 343 !! * Arguments 288 !!---------------------------------------------------------------------- 344 289 INTEGER, INTENT( in ) :: kt ! Main timestep counter 290 !!gm doctor jit ==> kit 345 291 INTEGER, INTENT( in ) :: jit ! Barotropic timestep counter (for timesplitting option) 346 347 !! * Local declarations 348 INTEGER :: jtide, jgrd, jb ! dummy loop indices 349 REAL(wp) :: arg, sarg 350 REAL(wp), DIMENSION(ntide_max) :: sist, cost 351 352 !!------------------------------------------------------------------------------ 353 !! OPA 9.0, LODYC-IPSL (2003) 354 !!------------------------------------------------------------------------------ 292 !! 293 INTEGER :: itide, igrd, ib ! dummy loop indices 294 REAL(wp) :: z_arg, z_sarg ! 295 REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost 296 !!---------------------------------------------------------------------- 355 297 356 298 ! Note tide phase speeds are in deg/hour, so we need to convert the 357 299 ! elapsed time in seconds to hours by dividing by 3600.0 358 if ( jit .eq. 0 ) then 359 arg = kt*rdt*rad/3600.0 360 else ! we are in a barotropic subcycle (for timesplitting option) 361 arg = ( (kt-1)*rdt + jit*rdtbt ) * rad/3600.0 362 endif 363 364 do jtide = 1,ntide 365 sarg = arg*tide_speed(jtide) 366 cost(jtide) = cos(sarg) 367 sist(jtide) = sin(sarg) 368 enddo 369 370 !! summing of tidal constituents into BDY arrays 371 300 IF( jit == 0 ) THEN 301 z_arg = kt * rdt * rad / 3600.0 302 ELSE ! we are in a barotropic subcycle (for timesplitting option) 303 z_arg = ( (kt-1) * rdt + jit * rdtbt ) * rad / 3600.0 304 ENDIF 305 306 DO itide = 1, ntide 307 z_sarg = z_arg * tide_speed(itide) 308 z_cost(itide) = COS( z_sarg ) 309 z_sist(itide) = SIN( z_sarg ) 310 END DO 311 312 ! summing of tidal constituents into BDY arrays 372 313 sshtide(:) = 0.0 373 utide(:) = 0.0 374 vtide(:) = 0.0 375 376 do jtide = 1 ,ntide 377 jgrd=1 !: SSH on tracer grid. 378 do jb = 1, nblenrim(jgrd) 379 sshtide(jb) =sshtide(jb)+ ssh1(jb,jtide)*cost(jtide) + ssh2(jb,jtide)*sist(jtide) 380 ! if(lwp) write(numout,*) 'z',jb,jtide,sshtide(jb), ssh1(jb,jtide),ssh2(jb,jtide) 381 enddo 382 jgrd=2 !: U grid 383 do jb=1, nblenrim(jgrd) 384 utide(jb) = utide(jb)+ u1(jb,jtide)*cost(jtide) + u2(jb,jtide)*sist(jtide) 385 ! if(lwp) write(numout,*) 'u',jb,jtide,utide(jb), u1(jb,jtide),u2(jb,jtide) 386 387 enddo 388 jgrd=3 !: V grid 389 do jb=1, nblenrim(jgrd) 390 vtide(jb) = vtide(jb)+ v1(jb,jtide)*cost(jtide) + v2(jb,jtide)*sist(jtide) 391 ! if(lwp) write(numout,*) 'v',jb,jtide,vtide(jb), v1(jb,jtide),v2(jb,jtide) 392 393 enddo 394 enddo 395 314 utide (:) = 0.0 315 vtide (:) = 0.0 316 ! 317 DO itide = 1, ntide 318 igrd=1 ! SSH on tracer grid. 319 DO ib = 1, nblenrim(igrd) 320 sshtide(ib) =sshtide(ib)+ ssh1(ib,itide)*z_cost(itide) + ssh2(ib,itide)*z_sist(itide) 321 ! if(lwp) write(numout,*) 'z',ib,itide,sshtide(ib), ssh1(ib,itide),ssh2(ib,itide) 322 END DO 323 igrd=2 ! U grid 324 DO ib=1, nblenrim(igrd) 325 utide(ib) = utide(ib)+ u1(ib,itide)*z_cost(itide) + u2(ib,itide)*z_sist(itide) 326 ! if(lwp) write(numout,*) 'u',ib,itide,utide(ib), u1(ib,itide),u2(ib,itide) 327 END DO 328 igrd=3 ! V grid 329 DO ib=1, nblenrim(igrd) 330 vtide(ib) = vtide(ib)+ v1(ib,itide)*z_cost(itide) + v2(ib,itide)*z_sist(itide) 331 ! if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), v1(ib,itide),v2(ib,itide) 332 END DO 333 END DO 334 ! 396 335 END SUBROUTINE tide_update 397 336 398 ! 399 ! 400 SUBROUTINE uvset (ihs,iday,imnth,iyr,f,vplu)401 !!---------------------------------------------------------------------- --------402 !! SUBROUTINE uvset403 !! ************************337 !!gm doctor naming of dummy argument variables!!! and all local variables 338 339 SUBROUTINE uvset( ihs, iday, imnth, iyr, f, z_vplu ) 340 !!---------------------------------------------------------------------- 341 !! *** SUBROUTINE uvset *** 342 !! 404 343 !! ** Purpose : - adjust tidal forcing for date factors 405 344 !! 406 !! 407 !! History : 408 !! 409 !! Origins POLCOMS v6.3 2007 410 !! NEMO v2.3 ! Jason Holt 411 !!------------------------------------------------------------------------------ 345 !!---------------------------------------------------------------------- 412 346 implicit none 413 !! * Arguments414 347 INTEGER, INTENT( in ) :: ihs ! Start time hours 415 348 INTEGER, INTENT( in ) :: iday ! start time days 416 INTEGER, INTENT( in ) :: imnth ! start time month 417 INTEGER, INTENT( in ) :: iyr ! start time year 418 INTEGER, PARAMETER :: nc =15 ! maximum number of constituents 419 420 REAL(WP) :: f(nc) ! nodal correction 421 REAL(WP) :: vplu(nc) ! phase correction 422 ! 423 INTEGER :: year,vd,ivdy,ndc,i,k 424 REAL(WP) :: u(nc),v(nc),zig(nc),rtd 425 REAL(WP) :: ss,h,p,en,p1 349 INTEGER, INTENT( in ) :: imnth ! start time month 350 INTEGER, INTENT( in ) :: iyr ! start time year 351 !! 352 !!gm nc is jptides_max ???? 353 INTEGER , PARAMETER :: nc =15 ! maximum number of constituents 426 354 CHARACTER(len=8), DIMENSION(nc) :: cname 427 ! 428 !!------------------------------------------------------------------------------ 429 !! OPA 9.0, LODYC-IPSL (2007) 430 !!------------------------------------------------------------------------------ 431 432 data cname/ 'q1', 'o1', 'p1', 's1', 'k1', & 433 '2n2','mu2', 'n2','nu2', 'm2', 'l2', 't2', 's2', 'k2', & 434 'm4'/ 435 data zig/.2338507481,.2433518789,.2610826055,.2617993878 & 436 , .2625161701 & 437 , .4868657873,.4881373225,.4963669182,.4976384533 & 438 , .5058680490,.5153691799,.5228820265,.5235987756 & 439 , .5250323419 & 440 , 1.011736098/ 355 INTEGER :: year, vd, ivdy, ndc, i, k 356 REAL(wp) :: ss, h, p, en, p1, rtd 357 REAL(wp), DIMENSION(nc) :: f ! nodal correction 358 REAL(wp), DIMENSION(nc) :: z_vplu ! phase correction 359 REAL(wp), DIMENSION(nc) :: u, v, zig 360 !! 361 DATA cname/ 'q1' , 'o1' , 'p1' , 's1' , 'k1' , & 362 & '2n2' , 'mu2' , 'n2' , 'nu2' , 'm2' , & 363 & 'l2' , 't2' , 's2' , 'k2' , 'm4' / 364 DATA zig/ .2338507481, .2433518789, .2610826055, .2617993878, .2625161701, & 365 & .4868657873, .4881373225, .4963669182, .4976384533, .5058680490, & 366 & .5153691799, .5228820265, .5235987756, .5250323419, 1.011736098 / 367 !!---------------------------------------------------------------------- 441 368 ! 442 369 ! ihs - start time gmt on ... 443 370 ! iday/imnth/iyr - date e.g. 12/10/87 444 371 ! 445 callvday(iday,imnth,iyr,ivdy)446 vd =ivdy372 CALL vday(iday,imnth,iyr,ivdy) 373 vd = ivdy 447 374 ! 448 375 !rp note change of year number for d. blackman shpen … … 453 380 !.....year = year of required data 454 381 !.....vd = day of required data..set up for 0000gmt day year 455 !456 382 ndc = nc 457 383 !.....ndc = number of constituents allowed 458 ! 459 rtd = 360.0/6.2831852 460 do i = 1,ndc 461 zig(i) = zig(i)*rtd 462 ! sigo(i)= zig(i) 463 enddo 464 ! 465 if(year == 0) then 466 return 467 endif 384 !!gm use rpi ? 385 rtd = 360.0 / 6.2831852 386 DO i = 1, ndc 387 zig(i) = zig(i)*rtd 388 ! sigo(i)= zig(i) 389 END DO 390 391 !!gm try to avoid RETURN in F90 392 IF( year == 0 ) RETURN 468 393 469 call shpen (year,vd,ss,h,p,en,p1) 470 call ufset (p,en,u,f) 471 call vset (ss,h,p,en,p1,v) 472 ! 473 do k=1,nc 474 475 vplu(k)=v(k)+u(k) 476 vplu(k)=vplu(k)+dble(ihs)*zig(k) 477 ! 478 do while ( vplu(k) < 0 ) 479 vplu(k) = vplu(k) + 360.0 480 enddo 481 ! 482 do while (vplu(k) > 360. ) 483 vplu(k) = vplu(k) - 360.0 484 enddo 485 ! 486 enddo 487 ! 394 CALL shpen( year, vd, ss, h , p , en, p1 ) 395 CALL ufset( p , en, u , f ) 396 CALL vset ( ss , h , p , en, p1, v ) 397 ! 398 DO k = 1, nc 399 z_vplu(k) = v(k) + u(k) 400 z_vplu(k) = z_vplu(k) + dble(ihs) * zig(k) 401 DO WHILE( z_vplu(k) < 0 ) 402 z_vplu(k) = z_vplu(k) + 360.0 403 END DO 404 DO WHILE( z_vplu(k) > 360. ) 405 z_vplu(k) = z_vplu(k) - 360.0 406 END DO 407 END DO 408 ! 488 409 END SUBROUTINE uvset 489 410 490 411 491 SUBROUTINE vday( iday,imnth,iy,ivdy)492 !!---------------------------------------------------------------------- --------493 !! SUBROUTINE vday494 !! ****************412 SUBROUTINE vday( iday, imnth, iy, ivdy ) 413 !!---------------------------------------------------------------------- 414 !! *** SUBROUTINE vday *** 415 !! 495 416 !! ** Purpose : - adjust tidal forcing for date factors 496 417 !! 497 !! 498 !! History :499 !!500 !! Origins POLCOMS v6.3 2007501 !! NEMO v2.3 ! Jason Holt418 !!---------------------------------------------------------------------- 419 INTEGER, INTENT(in ) :: iday, imnth, iy ! ???? 420 INTEGER, INTENT( out) :: ivdy ! ??? 421 !! 422 INTEGER :: iyr 502 423 !!------------------------------------------------------------------------------ 503 implicit none 504 ! 505 ! subroutine arguments 506 integer :: iday,imnth,iy,ivdy 507 ! 508 ! local variables 509 integer iyr 510 !!------------------------------------------------------------------------------ 511 !! NEMO 2.3, LODYC-IPSL (2008) 512 !!------------------------------------------------------------------------------ 513 514 ! ================================================================= 515 ! 516 ! calculate day number in year from day/month/year 517 ! 518 ! ================================================================= 424 425 !!gm nday_year in day mode is the variable compiuted here, no? 426 !!gm nday_year , & !: curent day counted from jan 1st of the current year 427 428 !calculate day number in year from day/month/year 519 429 if(imnth.eq.1) ivdy=iday 520 430 if(imnth.eq.2) ivdy=iday+31 … … 529 439 if(imnth.eq.11) ivdy=iday+304 530 440 if(imnth.eq.12) ivdy=iday+334 531 iyr=iy441 iyr=iy 532 442 if(mod(iyr,4).eq.0.and.imnth.gt.2) ivdy=ivdy+1 533 443 if(mod(iyr,100).eq.0.and.imnth.gt.2) ivdy=ivdy-1 534 444 if(mod(iyr,400).eq.0.and.imnth.gt.2) ivdy=ivdy+1 535 536 RETURN 537 END SUBROUTINE vday 538 539 SUBROUTINE shpen (year,vd,s,h,p,en,p1) 540 !!------------------------------------------------------------------------------ 541 !! SUBROUTINE shpen 542 !! ***************** 445 ! 446 END SUBROUTINE vday 447 448 !!doctor norme for dummy arguments 449 450 SUBROUTINE shpen( year, vd, s, h, p, en, p1 ) 451 !!---------------------------------------------------------------------- 452 !! *** SUBROUTINE shpen *** 453 !! 543 454 !! ** Purpose : - calculate astronomical arguments for tides 544 455 !! this version from d. blackman 30 nove 1990 545 !! 546 !! History : 547 !! 548 !! Origins POLCOMS v6.3 2007 549 !! NEMO v2.3 ! Jason Holt 550 !!------------------------------------------------------------------------------ 551 implicit none 552 553 ! subroutine arguments 554 integer ::year,vd 555 real(wp) :: s,h,p,en,p1 556 557 ! local variables 558 559 integer yr,ilc,icent,it,iday,ild,ipos,nn,iyd 456 !! 457 !!---------------------------------------------------------------------- 458 !!gm add INTENT in, out or inout.... DOCTOR name.... 459 !!gm please do not use variable name with a single letter: impossible to search in a code 460 INTEGER :: year,vd 461 REAL(wp) :: s,h,p,en,p1 462 !! 463 INTEGER :: yr,ilc,icent,it,iday,ild,ipos,nn,iyd 560 464 REAL(wp) :: cycle,t,td,delt(84),delta,deltat 561 562 !!------------------------------------------------------------------------------ 563 !! NEMO 2.3, LODYC-IPSL (2008) 564 !!------------------------------------------------------------------------------ 565 566 data delt /-5.04,-3.90,-2.87,-0.58,0.71,1.80, & 567 3.08, 4.63, 5.86, 7.21, 8.58,10.50,12.10, & 568 12.49,14.41,15.59,15.81,17.52,19.01,18.39, & 569 19.55,20.36,21.01,21.81,21.76,22.35,22.68, & 570 22.94,22.93,22.69,22.94,23.20,23.31,23.63, & 571 23.47,23.68,23.62,23.53,23.59,23.99,23.80, & 572 24.20,24.99,24.97,25.72,26.21,26.37,26.89, & 573 27.68,28.13,28.94,29.42,29.66,30.29,30.96, & 574 31.09,31.59,31.52,31.92,32.45,32.91,33.39, & 575 33.80,34.23,34.73,35.40,36.14,36.99,37.87, & 576 38.75,39.70,40.70,41.68,42.82,43.96,45.00, & 577 45.98,47.00,48.03,49.10,50.10,50.97,51.81, & 578 52.57/ 579 ! 465 !! 466 DATA delt /-5.04, -3.90, -2.87, -0.58, 0.71, 1.80, & 467 & 3.08, 4.63, 5.86, 7.21, 8.58, 10.50, 12.10, & 468 & 12.49, 14.41, 15.59, 15.81, 17.52, 19.01, 18.39, & 469 & 19.55, 20.36, 21.01, 21.81, 21.76, 22.35, 22.68, & 470 & 22.94, 22.93, 22.69, 22.94, 23.20, 23.31, 23.63, & 471 & 23.47, 23.68, 23.62, 23.53, 23.59, 23.99, 23.80, & 472 & 24.20, 24.99, 24.97, 25.72, 26.21, 26.37, 26.89, & 473 & 27.68, 28.13, 28.94, 29.42, 29.66, 30.29, 30.96, & 474 & 31.09, 31.59, 31.52, 31.92, 32.45, 32.91, 33.39, & 475 & 33.80, 34.23, 34.73, 35.40, 36.14, 36.99, 37.87, & 476 & 38.75, 39.70, 40.70, 41.68, 42.82, 43.96, 45.00, & 477 & 45.98, 47.00, 48.03, 49.10, 50.10, 50.97, 51.81, & 478 & 52.57 / 479 !!---------------------------------------------------------------------- 480 580 481 cycle = 360.0 581 482 ilc = 0 582 icent = year /100583 yr = year - icent *100483 icent = year / 100 484 yr = year - icent * 100 584 485 t = icent - 20 585 486 ! … … 588 489 ! see notes by cartwright 589 490 ! 491 !!gm old coding style, use CASE instead and avoid GOTO (obsolescence in fortran 90) 492 !!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 590 493 it = icent - 20 591 494 if (it) 1,2,2 … … 602 505 td = 0.0 603 506 ! 507 !!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 604 508 if (yr) 4,5,4 605 509 ! … … 623 527 7 deltat = delta * 1.0e-6 624 528 ! 529 !!gm precision of the computation : example for s it should be replace by: 530 !!gm s = 218.3165 + (481267.8813 - 0.0016*t)*t + 152.0*deltat ==> more precise modify the last digits results 625 531 s = 218.3165 + 481267.8813*t - 0.0016*t*t + 152.0*deltat 626 h = 280.4661 + 36000.7698 *t + 0.0003*t*t +11.0*deltat627 p = 83.3535 + 4069.0139 *t - 0.0103*t*t +deltat628 en = 234.9555 + 1934.1363 *t - 0.0021*t*t +deltat629 p1 = 282.9384 + 1.7195 *t + 0.0005*t*t630 !631 nn = s /cycle632 s = s - nn *cycle633 if ( s .lt. 0.0) s = s+cycle634 !635 nn = h /cycle636 h = h -cycle*nn637 if (h .lt. 0.0) h = h+cycle638 !639 nn = p /cycle640 p = p - cycle*nn641 if (p .lt. 0.0) p = p+cycle642 !643 nn = en /cycle644 en = en -cycle*nn645 if(en .lt. 0.0)en = en + cycle532 h = 280.4661 + 36000.7698 *t + 0.0003*t*t + 11.0*deltat 533 p = 83.3535 + 4069.0139 *t - 0.0103*t*t + deltat 534 en = 234.9555 + 1934.1363 *t - 0.0021*t*t + deltat 535 p1 = 282.9384 + 1.7195 *t + 0.0005*t*t 536 ! 537 nn = s / cycle 538 s = s - nn * cycle 539 IF( s < 0.e0 ) s = s + cycle 540 ! 541 nn = h / cycle 542 h = h - cycle * nn 543 IF( h < 0.e0 ) h = h + cycle 544 ! 545 nn = p / cycle 546 p = p - cycle * nn 547 IF( p < 0.e0) p = p + cycle 548 ! 549 nn = en / cycle 550 en = en - cycle * nn 551 IF( en < 0.e0 ) en = en + cycle 646 552 en = cycle - en 647 ! 648 nn = p1/cycle 649 p1 = p1 - nn*cycle 650 ! 651 RETURN 652 ! 653 END SUBROUTINE shpen 654 655 656 SUBROUTINE ufset (p,cn,b,a) 657 !!------------------------------------------------------------------------------ 658 !! SUBROUTINE ufset 659 !! ***************** 553 ! 554 nn = p1 / cycle 555 p1 = p1 - nn * cycle 556 ! 557 END SUBROUTINE shpen 558 559 560 SUBROUTINE ufset( p, cn, b, a ) 561 !!---------------------------------------------------------------------- 562 !! *** SUBROUTINE ufset *** 563 !! 660 564 !! ** Purpose : - calculate nodal parameters for the tides 661 565 !! 662 !! 663 !! History : 664 !! 665 !! Origins POLCOMS v6.3 2007 666 !! NEMO v2.3 ! Jason Holt 667 !!------------------------------------------------------------------------------ 668 669 implicit none 566 !!---------------------------------------------------------------------- 567 !!gm doctor naming of dummy argument variables!!! and all local variables 568 !!gm nc is jptides_max ???? 670 569 integer nc 671 570 parameter (nc=15) 672 ! subroutine arguments 673 real(wp) p,cn 674 ! 675 ! 676 ! local variables 677 real(wp) :: w1,w2,w3,w4,w5,w6,w7,w8,nw,pw,rad 678 real(wp) :: a(nc),b(nc) 679 integer ndc,k 680 ! 681 !!------------------------------------------------------------------------------ 682 !! NEMO 2.3, LODYC-IPSL (2008) 683 !!------------------------------------------------------------------------------ 684 685 ndc=nc 686 ! 571 REAL(wp) p,cn 572 !! 573 574 !!gm rad is already a public variable defined in phycst.F90 .... ==> doctor norme local real start with "z" 575 REAL(wp) :: w1, w2, w3, w4, w5, w6, w7, w8, nw, pw, rad 576 REAL(wp) :: a(nc), b(nc) 577 INTEGER :: ndc, k 578 !!---------------------------------------------------------------------- 579 580 ndc = nc 581 687 582 ! a=f , b =u 688 583 ! t is zero as compared to tifa. 584 !! use rad defined in phycst (i.e. add a USE phycst at the begining of the module 689 585 rad = 6.2831852d0/360.0 690 pw = p *rad691 nw = cn *rad692 w1 = cos( nw)693 w2 = cos( 2*nw)694 w3 = cos( 3*nw)695 w4 = sin( nw)696 w5 = sin( 2*nw)697 w6 = sin( 3*nw)698 w7 = 1 -0.2505*cos(2*pw) -0.1102*cos(2*pw-nw)&699 -0.156*cos(2*pw-2*nw) -0.037*cos(nw)700 w8 = -0.2505*sin(2*pw) -0.1102*sin(2*pw-nw)&701 -0.0156*sin(2*pw-2*nw) -0.037*sin(nw)702 !703 a(1) = 1.0089 +0.1871*w1-0.0147*w2+0.0014*w3704 b(1) = 0.1885*w4 - 0.0234*w5+.0033*w6705 ! q1586 pw = p * rad 587 nw = cn * rad 588 w1 = cos( nw ) 589 w2 = cos( 2*nw ) 590 w3 = cos( 3*nw ) 591 w4 = sin( nw ) 592 w5 = sin( 2*nw ) 593 w6 = sin( 3*nw ) 594 w7 = 1. - 0.2505 * COS( 2*pw ) - 0.1102 * COS(2*pw-nw ) & 595 & - 0.156 * COS( 2*pw-2*nw ) - 0.037 * COS( nw ) 596 w8 = - 0.2505 * SIN( 2*pw ) - 0.1102 * SIN(2*pw-nw ) & 597 & - 0.0156 * SIN( 2*pw-2*nw ) - 0.037 * SIN( nw ) 598 ! 599 a(1) = 1.0089 + 0.1871 * w1 - 0.0147 * w2 + 0.0014 * w3 600 b(1) = 0.1885 * w4 - 0.0234 * w5 + 0.0033 * w6 601 ! q1 706 602 a(2) = a(1) 707 603 b(2) = b(1) 708 ! o1604 ! o1 709 605 a(3) = 1.0 710 606 b(3) = 0.0 711 ! p1607 ! p1 712 608 a(4) = 1.0 713 609 b(4) = 0.0 714 ! s1610 ! s1 715 611 a(5) = 1.0060+0.1150*w1- 0.0088*w2 +0.0006*w3 716 612 b(5) = -0.1546*w4 + 0.0119*w5 -0.0012*w6 717 ! k1613 ! k1 718 614 a(6) =1.0004 -0.0373*w1+ 0.0002*w2 719 615 b(6) = -0.0374*w4 720 ! 2n2616 ! 2n2 721 617 a(7) = a(6) 722 618 b(7) = b(6) 723 ! mu2619 ! mu2 724 620 a(8) = a(6) 725 621 b(8) = b(6) 726 ! n2622 ! n2 727 623 a(9) = a(6) 728 624 b(9) = b(6) 729 ! nu2625 ! nu2 730 626 a(10) = a(6) 731 627 b(10) = b(6) 732 ! m2 733 a(11) = sqrt(w7*w7+w8*w8) 734 b(11) = atan(w8/w7) 735 if(w7.lt.0) b(11) = b(11) + 3.141992 736 ! l2 628 ! m2 629 a(11) = SQRT( w7 * w7 + w8 * w8 ) 630 b(11) = ATAN( w8 / w7 ) 631 !!gmuse rpi instead of 3.141992 ??? true pi is rpi=3.141592653589793_wp ..... ???? 632 IF( w7 < 0.e0 ) b(11) = b(11) + 3.141992 633 ! l2 737 634 a(12) = 1.0 738 635 b(12) = 0.0 739 ! t2636 ! t2 740 637 a(13)= a(12) 741 638 b(13)= b(12) 742 ! s2639 ! s2 743 640 a(14) = 1.0241+0.2863*w1+0.0083*w2 -0.0015*w3 744 641 b(14) = -0.3096*w4 + 0.0119*w5 - 0.0007*w6 745 ! k2642 ! k2 746 643 a(15) = a(6)*a(6) 747 644 b(15) = 2*b(6) 748 ! m4 749 ! 750 do 40 k = 1,ndc 751 b(k) = b(k)/rad 752 32 if (b(k)) 34,35,35 753 34 b(k) = b(k) + 360.0 754 go to 32 755 35 if (b(k)-360.0) 40,37,37 756 37 b(k) = b(k)-360.0 757 go to 35 645 ! m4 646 !!gm old coding, remove GOTO and label of lines 647 !!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 648 DO 40 k = 1,ndc 649 b(k) = b(k)/rad 650 32 if (b(k)) 34,35,35 651 34 b(k) = b(k) + 360.0 652 go to 32 653 35 if (b(k)-360.0) 40,37,37 654 37 b(k) = b(k)-360.0 655 go to 35 758 656 40 continue 759 RETURN 760 END SUBROUTINE ufset 761 762 SUBROUTINE vset( s,h,p,en,p1,v) 763 !!------------------------------------------------------------------------------ 764 !! SUBROUTINE vset 765 !! **************** 657 ! 658 END SUBROUTINE ufset 659 660 661 SUBROUTINE vset( s,h,p,en,p1,v) 662 !!---------------------------------------------------------------------- 663 !! *** SUBROUTINE vset *** 664 !! 766 665 !! ** Purpose : - calculate tidal phases for 0000gmt on start day of run 767 666 !! 768 !! 769 !! History : 770 !! 771 !! Origins POLCOMS v6.3 2007 772 !! NEMO v2.3 ! Jason Holt 773 !!------------------------------------------------------------------------------ 774 implicit none 667 !!---------------------------------------------------------------------- 668 !!gm doctor naming of dummy argument variables!!! and all local variables 669 !!gm nc is jptides_max ???? 670 !!gm en argument is not used: suppress it ? 775 671 integer nc 776 672 parameter (nc=15) 777 ! subroutine arguments778 673 real(wp) s,h,p,en,p1 779 674 real(wp) v(nc) 780 ! 781 integer ndc,k 782 !!------------------------------------------------------------------------------ 783 !! NEMO 2.3, LODYC-IPSL (2008) 784 !!------------------------------------------------------------------------------ 675 !! 676 integer ndc, k 677 !!---------------------------------------------------------------------- 785 678 786 679 ndc = nc 787 ! v s are computed here. 788 v(1) =-3*s +h +p +270 789 ! q1 790 v(2) =-2*s +h +270.0 791 ! o1 792 v(3) =-h +270 793 ! p1 794 v(4) =180 795 ! s1 796 v(5) =h +90.0 797 ! k1 798 v(6) =-4*s +2*h +2*p 799 ! 2n2 800 v(7) =-4*(s-h) 801 ! mu2 802 v(8) =-3*s +2*h +p 803 ! n2 804 v(9) =-3*s +4*h -p 805 ! mu2 806 v(10) =-2*s +2*h 807 ! m2 808 v(11) =-s +2*h -p +180 809 ! l2 810 v(12) =-h +p1 811 ! t2 812 v(13) =0 813 ! s2 814 v(14) =h+h 815 ! k2 816 v(15) =2*v(10) 817 ! m4 818 ! 680 ! v s are computed here. 681 v(1) =-3*s +h +p +270 ! Q1 682 v(2) =-2*s +h +270.0 ! O1 683 v(3) =-h +270 ! P1 684 v(4) =180 ! S1 685 v(5) =h +90.0 ! K1 686 v(6) =-4*s +2*h +2*p ! 2N2 687 v(7) =-4*(s-h) ! MU2 688 v(8) =-3*s +2*h +p ! N2 689 v(9) =-3*s +4*h -p ! MU2 690 v(10) =-2*s +2*h ! M2 691 v(11) =-s +2*h -p +180 ! L2 692 v(12) =-h +p1 ! T2 693 v(13) =0 ! S2 694 v(14) =h+h ! K2 695 v(15) =2*v(10) ! M4 696 ! 697 !!gm old coding, remove GOTO and label of lines 698 !!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90 819 699 do 72 k = 1, ndc 820 69 if (v(k) )70,71,71821 70 v(k)=v(k)+360.0822 go to 69823 71 if ( v(k) - 360.0)72,73,73824 73 v(k)=v(k)-360.0825 go to 71700 69 if( v(k) ) 70,71,71 701 70 v(k) = v(k)+360.0 702 go to 69 703 71 if( v(k) - 360.0 ) 72,73,73 704 73 v(k) = v(k)-360.0 705 go to 71 826 706 72 continue 827 RETURN828 707 ! 708 END SUBROUTINE vset 829 709 830 710 #else 831 !!================================================================================= 832 !! *** MODULE bdytides *** 833 !!================================================================================= 834 835 LOGICAL, PUBLIC, PARAMETER :: lk_bdy_tides = .FALSE. !: tidal forcing at boundaries. 836 837 CHARACTER(len=80), PUBLIC :: & 838 filtide !: Filename root for tidal input files 839 840 CHARACTER(len=4), PUBLIC, DIMENSION(1) :: & 841 tide_cpt !: Names of tidal components used. 711 !!---------------------------------------------------------------------- 712 !! Dummy module NO Unstruct Open Boundary Conditions for tides 713 !!---------------------------------------------------------------------- 714 !!gm are you sure we need to define filtide and tide_cpt ? 715 CHARACTER(len=80), PUBLIC :: filtide !: Filename root for tidal input files 716 CHARACTER(len=4 ), PUBLIC, DIMENSION(1) :: tide_cpt !: Names of tidal components used. 842 717 843 718 CONTAINS 844 845 SUBROUTINE tide_init 846 ! No tidal forcing at boundaries ==> empty routine 719 SUBROUTINE tide_init ! Empty routine 847 720 END SUBROUTINE tide_init 848 849 SUBROUTINE tide_data 850 ! No tidal forcing at boundaries ==> empty routine 721 SUBROUTINE tide_data ! Empty routine 851 722 END SUBROUTINE tide_data 852 853 SUBROUTINE tide_update( kt, jit ) 854 INTEGER :: kt, jit 855 WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, jit 856 857 ! No tidal forcing at boundaries ==> empty routine 723 SUBROUTINE tide_update( kt, kit ) ! Empty routine 724 WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, kit 858 725 END SUBROUTINE tide_update 859 860 726 #endif 861 727 728 !!====================================================================== 862 729 END MODULE bdytides
Note: See TracChangeset
for help on using the changeset viewer.
