Changeset 3367 for branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO
- Timestamp:
- 2012-04-25T12:21:21+02:00 (12 years ago)
- Location:
- branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 1 added
- 12 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdy_oce.F90
r3294 r3367 73 73 INTEGER, DIMENSION(jp_bdy) :: nn_tra_dta !: = 0 use the initial state as bdy dta ; 74 74 !: = 1 read it in a NetCDF file 75 LOGICAL, DIMENSION(jp_bdy) :: ln_tra_dmp 76 LOGICAL, DIMENSION(jp_bdy) :: ln_dyn3d_dmp 77 REAL, DIMENSION(jp_bdy) :: rn_time_dmp 78 LOGICAL :: ln_tra_dmp_tot 79 LOGICAL :: ln_dyn3d_dmp_tot 80 75 81 #if defined key_lim2 76 82 INTEGER, DIMENSION(jp_bdy) :: nn_ice_lim2 ! Choice of boundary condition for sea ice variables -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
r3294 r3367 32 32 USE ice_2 33 33 #endif 34 USE sbcapr 34 35 35 36 IMPLICIT NONE … … 108 109 109 110 IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 0 ) THEN 110 IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 111 ilen1(:) = nblen(:) 112 ELSE 113 ilen1(:) = nblenrim(:) 114 ENDIF 111 ilen1(:) = nblen(:) 115 112 igrd = 1 116 113 DO ib = 1, ilen1(igrd) … … 134 131 135 132 IF( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN 136 IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 137 ilen1(:) = nblen(:) 138 ELSE 139 ilen1(:) = nblenrim(:) 140 ENDIF 133 ilen1(:) = nblen(:) 141 134 igrd = 2 142 135 DO ib = 1, ilen1(igrd) … … 158 151 159 152 IF( nn_tra(ib_bdy) .gt. 0 .and. nn_tra_dta(ib_bdy) .eq. 0 ) THEN 160 IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 161 ilen1(:) = nblen(:) 162 ELSE 163 ilen1(:) = nblenrim(:) 164 ENDIF 153 ilen1(:) = nblen(:) 165 154 igrd = 1 ! Everything is at T-points here 166 155 DO ib = 1, ilen1(igrd) … … 176 165 #if defined key_lim2 177 166 IF( nn_ice_lim2(ib_bdy) .gt. 0 .and. nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 178 IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 179 ilen1(:) = nblen(:) 180 ELSE 181 ilen1(:) = nblenrim(:) 182 ENDIF 167 ilen1(:) = nblen(:) 183 168 igrd = 1 ! Everything is at T-points here 184 169 DO ib = 1, ilen1(igrd) … … 214 199 dta_bdy(ib_bdy)%v2d(:) = 0.0 215 200 ENDIF 216 IF( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) THEN ! update external data 217 jend = jstart + 2 218 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), & 219 & jit=jit, time_offset=time_offset ) 220 ENDIF 221 IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 222 CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), & 223 & jit=jit, time_offset=time_offset ) 201 IF (nn_tra(ib_bdy).ne.4) THEN 202 IF (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3.or.(ln_full_vel_array(ib_bdy).and.nn_dyn3d_dta(ib_bdy).eq.1)) THEN 203 ! For the runoff case, no need to update the forcing (already done in the baroclinic part) 204 jend = nb_bdy_fld(ib_bdy) 205 IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend - 2 206 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), jit=jit, time_offset=time_offset ) 207 IF ( nn_tra(ib_bdy) .GT. 0 .AND. nn_tra_dta(ib_bdy) .GE. 1 ) jend = jend + 2 208 ! If full velocities in boundary data then split into barotropic and baroclinic data 209 IF( ln_full_vel_array(ib_bdy) .and. & 210 & ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN 211 igrd = 2 ! zonal velocity 212 dta_bdy(ib_bdy)%u2d(:) = 0.0 213 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 214 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 215 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 216 DO ik = 1, jpkm1 217 dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) & 218 & + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik) 219 END DO 220 dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij) 221 DO ik = 1, jpkm1 222 dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib) 223 END DO 224 END DO 225 igrd = 3 ! meridional velocity 226 dta_bdy(ib_bdy)%v2d(:) = 0.0 227 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 228 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 229 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 230 DO ik = 1, jpkm1 231 dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) & 232 & + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik) 233 END DO 234 dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij) 235 DO ik = 1, jpkm1 236 dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib) 237 END DO 238 END DO 239 ENDIF 240 ENDIF 241 IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 242 CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), & 243 & jit=jit, time_offset=time_offset ) 244 ENDIF 224 245 ENDIF 225 246 ENDIF 226 247 ELSE 227 IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 228 dta_bdy(ib_bdy)%ssh(:) = 0.0 229 dta_bdy(ib_bdy)%u2d(:) = 0.0 230 dta_bdy(ib_bdy)%v2d(:) = 0.0 248 IF (nn_tra(ib_bdy).eq.4) then ! runoff condition 249 jend = nb_bdy_fld(ib_bdy) 250 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), time_offset=time_offset ) 251 ! 252 igrd = 2 ! zonal velocity 253 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 254 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 255 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 256 dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) / ( e2u(ii,ij) * hu_0(ii,ij) ) 257 END DO 258 ! 259 igrd = 3 ! meridional velocity 260 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 261 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 262 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 263 dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) / ( e1v(ii,ij) * hv_0(ii,ij) ) 264 END DO 265 ELSE 266 IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .eq. 2 ) THEN ! tidal harmonic forcing ONLY: initialise arrays 267 dta_bdy(ib_bdy)%ssh(:) = 0.0 268 dta_bdy(ib_bdy)%u2d(:) = 0.0 269 dta_bdy(ib_bdy)%v2d(:) = 0.0 270 ENDIF 271 IF( nb_bdy_fld(ib_bdy) .gt. 0 ) THEN ! update external data 272 jend = nb_bdy_fld(ib_bdy) 273 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), time_offset=time_offset ) 274 ENDIF 275 ! If full velocities in boundary data then split into barotropic and baroclinic data 276 IF( ln_full_vel_array(ib_bdy) .and. & 277 & ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN 278 igrd = 2 ! zonal velocity 279 dta_bdy(ib_bdy)%u2d(:) = 0.0 280 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 281 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 282 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 283 DO ik = 1, jpkm1 284 dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) & 285 & + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik) 286 END DO 287 dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij) 288 DO ik = 1, jpkm1 289 dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib) 290 END DO 291 END DO 292 igrd = 3 ! meridional velocity 293 dta_bdy(ib_bdy)%v2d(:) = 0.0 294 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 295 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 296 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 297 DO ik = 1, jpkm1 298 dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) & 299 & + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik) 300 END DO 301 dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij) 302 DO ik = 1, jpkm1 303 dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib) 304 END DO 305 END DO 306 ENDIF 307 IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing 308 CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), time_offset=time_offset ) 309 ENDIF 231 310 ENDIF 232 IF( nb_bdy_fld(ib_bdy) .gt. 0 ) THEN ! update external data233 jend = jstart + nb_bdy_fld(ib_bdy) - 1234 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), map=nbmap_ptr(jstart:jend), time_offset=time_offset )235 ENDIF236 IF( nn_dyn2d(ib_bdy) .gt. 0 .and. nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN ! update tidal harmonic forcing237 CALL tide_update( kt=kt, idx=idx_bdy(ib_bdy), dta=dta_bdy(ib_bdy), td=tides(ib_bdy), time_offset=time_offset )238 ENDIF239 311 ENDIF 240 312 jstart = jend+1 241 242 ! If full velocities in boundary data then split into barotropic and baroclinic data 243 ! (Note that we have already made sure that you can't use ln_full_vel = .true. at the same 244 ! time as the dynspg_ts option). 245 246 IF( ln_full_vel_array(ib_bdy) .and. & 247 & ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 .or. nn_dyn3d_dta(ib_bdy) .eq. 1 ) ) THEN 248 249 igrd = 2 ! zonal velocity 250 dta_bdy(ib_bdy)%u2d(:) = 0.0 251 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 313 END IF ! nn_dta(ib_bdy) = 1 314 END DO ! ib_bdy 315 316 IF ( ln_apr_obc ) THEN 317 DO ib_bdy = 1, nb_bdy 318 IF (nn_tra(ib_bdy).NE.4)THEN 319 igrd = 1 ! meridional velocity 320 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 252 321 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 253 322 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 254 DO ik = 1, jpkm1 255 dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) & 256 & + fse3u(ii,ij,ik) * umask(ii,ij,ik) * dta_bdy(ib_bdy)%u3d(ib,ik) 257 END DO 258 dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) * hur(ii,ij) 259 DO ik = 1, jpkm1 260 dta_bdy(ib_bdy)%u3d(ib,ik) = dta_bdy(ib_bdy)%u3d(ib,ik) - dta_bdy(ib_bdy)%u2d(ib) 261 END DO 262 END DO 263 264 igrd = 3 ! meridional velocity 265 dta_bdy(ib_bdy)%v2d(:) = 0.0 266 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 267 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 268 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 269 DO ik = 1, jpkm1 270 dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) & 271 & + fse3v(ii,ij,ik) * vmask(ii,ij,ik) * dta_bdy(ib_bdy)%v3d(ib,ik) 272 END DO 273 dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) * hvr(ii,ij) 274 DO ik = 1, jpkm1 275 dta_bdy(ib_bdy)%v3d(ib,ik) = dta_bdy(ib_bdy)%v3d(ib,ik) - dta_bdy(ib_bdy)%v2d(ib) 276 END DO 277 END DO 278 279 ENDIF 280 281 END IF ! nn_dta(ib_bdy) = 1 282 END DO ! ib_bdy 323 dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + ssh_ib(ii,ij) 324 ENDDO 325 ENDIF 326 ENDDO 327 ENDIF 283 328 284 329 IF( nn_timing == 1 ) CALL timing_stop('bdy_dta') … … 326 371 IF( nn_timing == 1 ) CALL timing_start('bdy_dta_init') 327 372 373 IF(lwp) WRITE(numout,*) 374 IF(lwp) WRITE(numout,*) 'bdy_dta_ini : initialization of data at the open boundaries' 375 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 376 IF(lwp) WRITE(numout,*) '' 377 328 378 ! Set nn_dta 329 379 DO ib_bdy = 1, nb_bdy 330 nn_dta(ib_bdy) = MAX( nn_dyn2d_dta(ib_bdy) & 331 ,nn_dyn3d_dta(ib_bdy) & 332 ,nn_tra_dta(ib_bdy) & 333 #if defined key_lim2 334 ,nn_ice_lim2_dta(ib_bdy) & 335 #endif 336 ) 337 IF(nn_dta(ib_bdy) .gt. 1) nn_dta(ib_bdy) = 1 380 IF (nn_dyn2d_dta(ib_bdy).eq.1.OR.nn_dyn2d_dta(ib_bdy).eq.3) nn_dta(ib_bdy) = 1 381 IF (nn_dyn3d_dta(ib_bdy).ge.1) nn_dta(ib_bdy) = 1 382 IF (nn_tra_dta (ib_bdy).ge.1) nn_dta(ib_bdy) = 1 383 #if defined key_lim2 384 IF (nn_ice_lim2_dta(ib_bdy).eq.1) nn_dta(ib_bdy) = 1 385 #endif 386 IF (lwp) THEN 387 WRITE(numout,*) 'Bdy package number ',ib_bdy 388 IF (nn_dta(ib_bdy).eq.1) THEN 389 WRITE(numout,*) 'We use external data' 390 ELSE 391 WRITE(numout,*) 'We do not use external data' 392 ENDIF 393 ENDIF 338 394 END DO 339 395 … … 357 413 ENDIF 358 414 #endif 415 IF(lwp) WRITE(numout,*) 'Maximum number of files to open =',nb_bdy_fld(ib_bdy) 359 416 ENDDO 360 417 … … 408 465 ln_full_vel_array(ib_bdy) = ln_full_vel 409 466 410 IF( ln_full_vel_array(ib_bdy) .and. lk_dynspg_ts ) THEN411 CALL ctl_stop( 'bdy_dta_init: ERROR, cannot specify full velocities in boundary data',&412 & 'with dynspg_ts option' ) ; RETURN413 ENDIF414 415 467 nblen => idx_bdy(ib_bdy)%nblen 416 468 nblenrim => idx_bdy(ib_bdy)%nblenrim … … 420 472 IF( nn_dyn2d(ib_bdy) .gt. 0 .and. ( nn_dyn2d_dta(ib_bdy) .eq. 1 .or. nn_dyn2d_dta(ib_bdy) .eq. 3 ) ) THEN 421 473 422 IF( nn_ dyn2d(ib_bdy) .ne. jp_frs ) THEN474 IF( nn_tra(ib_bdy) .ne. 4 ) THEN ! runoff condition : no ssh reading 423 475 jfld = jfld + 1 424 476 blf_i(jfld) = bn_ssh 425 477 ibdy(jfld) = ib_bdy 426 478 igrid(jfld) = 1 427 ilen1(jfld) = nblen rim(igrid(jfld))479 ilen1(jfld) = nblen(igrid(jfld)) 428 480 ilen3(jfld) = 1 429 481 ENDIF 430 482 431 483 IF( .not. ln_full_vel_array(ib_bdy) ) THEN 432 433 484 jfld = jfld + 1 434 485 blf_i(jfld) = bn_u2d 435 486 ibdy(jfld) = ib_bdy 436 487 igrid(jfld) = 2 437 IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 438 ilen1(jfld) = nblen(igrid(jfld)) 439 ELSE 440 ilen1(jfld) = nblenrim(igrid(jfld)) 441 ENDIF 488 ilen1(jfld) = nblen(igrid(jfld)) 442 489 ilen3(jfld) = 1 443 490 … … 446 493 ibdy(jfld) = ib_bdy 447 494 igrid(jfld) = 3 448 IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 449 ilen1(jfld) = nblen(igrid(jfld)) 450 ELSE 451 ilen1(jfld) = nblenrim(igrid(jfld)) 452 ENDIF 495 ilen1(jfld) = nblen(igrid(jfld)) 453 496 ilen3(jfld) = 1 454 455 497 ENDIF 456 498 … … 466 508 ibdy(jfld) = ib_bdy 467 509 igrid(jfld) = 2 468 IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 469 ilen1(jfld) = nblen(igrid(jfld)) 470 ELSE 471 ilen1(jfld) = nblenrim(igrid(jfld)) 472 ENDIF 510 ilen1(jfld) = nblen(igrid(jfld)) 473 511 ilen3(jfld) = jpk 474 512 … … 477 515 ibdy(jfld) = ib_bdy 478 516 igrid(jfld) = 3 479 IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 480 ilen1(jfld) = nblen(igrid(jfld)) 481 ELSE 482 ilen1(jfld) = nblenrim(igrid(jfld)) 483 ENDIF 517 ilen1(jfld) = nblen(igrid(jfld)) 484 518 ilen3(jfld) = jpk 485 519 … … 493 527 ibdy(jfld) = ib_bdy 494 528 igrid(jfld) = 1 495 IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 496 ilen1(jfld) = nblen(igrid(jfld)) 497 ELSE 498 ilen1(jfld) = nblenrim(igrid(jfld)) 499 ENDIF 529 ilen1(jfld) = nblen(igrid(jfld)) 500 530 ilen3(jfld) = jpk 501 531 … … 504 534 ibdy(jfld) = ib_bdy 505 535 igrid(jfld) = 1 506 IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 507 ilen1(jfld) = nblen(igrid(jfld)) 508 ELSE 509 ilen1(jfld) = nblenrim(igrid(jfld)) 510 ENDIF 536 ilen1(jfld) = nblen(igrid(jfld)) 511 537 ilen3(jfld) = jpk 512 538 … … 521 547 ibdy(jfld) = ib_bdy 522 548 igrid(jfld) = 1 523 IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 524 ilen1(jfld) = nblen(igrid(jfld)) 525 ELSE 526 ilen1(jfld) = nblenrim(igrid(jfld)) 527 ENDIF 549 ilen1(jfld) = nblen(igrid(jfld)) 528 550 ilen3(jfld) = 1 529 551 … … 532 554 ibdy(jfld) = ib_bdy 533 555 igrid(jfld) = 1 534 IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 535 ilen1(jfld) = nblen(igrid(jfld)) 536 ELSE 537 ilen1(jfld) = nblenrim(igrid(jfld)) 538 ENDIF 556 ilen1(jfld) = nblen(igrid(jfld)) 539 557 ilen3(jfld) = 1 540 558 … … 543 561 ibdy(jfld) = ib_bdy 544 562 igrid(jfld) = 1 545 IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 546 ilen1(jfld) = nblen(igrid(jfld)) 547 ELSE 548 ilen1(jfld) = nblenrim(igrid(jfld)) 549 ENDIF 563 ilen1(jfld) = nblen(igrid(jfld)) 550 564 ilen3(jfld) = 1 551 565 … … 566 580 ENDDO ! ib_bdy 567 581 568 569 582 DO jfld = 1, nb_bdy_fld_sum 570 583 ALLOCATE( bf(jfld)%fnow(ilen1(jfld),1,ilen3(jfld)) ) … … 577 590 jstart = 1 578 591 DO ib_bdy = 1, nb_bdy 579 jend = jstart + nb_bdy_fld(ib_bdy) - 1592 jend = nb_bdy_fld(ib_bdy) 580 593 CALL fld_fill( bf(jstart:jend), blf_i(jstart:jend), cn_dir_array(ib_bdy), 'bdy_dta', & 581 594 & 'open boundary conditions', 'nambdy_dta' ) … … 596 609 IF (nn_dyn2d(ib_bdy) .gt. 0) THEN 597 610 IF( nn_dyn2d_dta(ib_bdy) .eq. 0 .or. nn_dyn2d_dta(ib_bdy) .eq. 2 .or. ln_full_vel_array(ib_bdy) ) THEN 598 IF( nn_dyn2d(ib_bdy) .eq. jp_frs ) THEN 599 ilen0(1:3) = nblen(1:3) 600 ELSE 601 ilen0(1:3) = nblenrim(1:3) 602 ENDIF 603 ALLOCATE( dta_bdy(ib_bdy)%ssh(ilen0(1)) ) 611 ilen0(1:3) = nblen(1:3) 604 612 ALLOCATE( dta_bdy(ib_bdy)%u2d(ilen0(2)) ) 605 613 ALLOCATE( dta_bdy(ib_bdy)%v2d(ilen0(3)) ) 614 IF (nn_dyn2d_dta(ib_bdy).eq.1.or.nn_dyn2d_dta(ib_bdy).eq.3) THEN 615 jfld = jfld + 1 616 dta_bdy(ib_bdy)%ssh => bf(jfld)%fnow(:,1,1) 617 ELSE 618 ALLOCATE( dta_bdy(ib_bdy)%ssh(nblen(1)) ) 619 ENDIF 606 620 ELSE 607 621 IF( nn_dyn2d(ib_bdy) .ne. jp_frs ) THEN … … 617 631 618 632 IF ( nn_dyn3d(ib_bdy) .gt. 0 .and. nn_dyn3d_dta(ib_bdy) .eq. 0 ) THEN 619 IF( nn_dyn3d(ib_bdy) .eq. jp_frs ) THEN 620 ilen0(1:3) = nblen(1:3) 621 ELSE 622 ilen0(1:3) = nblenrim(1:3) 623 ENDIF 633 ilen0(1:3) = nblen(1:3) 624 634 ALLOCATE( dta_bdy(ib_bdy)%u3d(ilen0(2),jpk) ) 625 635 ALLOCATE( dta_bdy(ib_bdy)%v3d(ilen0(3),jpk) ) … … 636 646 IF (nn_tra(ib_bdy) .gt. 0) THEN 637 647 IF( nn_tra_dta(ib_bdy) .eq. 0 ) THEN 638 IF( nn_tra(ib_bdy) .eq. jp_frs ) THEN 639 ilen0(1:3) = nblen(1:3) 640 ELSE 641 ilen0(1:3) = nblenrim(1:3) 642 ENDIF 648 ilen0(1:3) = nblen(1:3) 643 649 ALLOCATE( dta_bdy(ib_bdy)%tem(ilen0(1),jpk) ) 644 650 ALLOCATE( dta_bdy(ib_bdy)%sal(ilen0(1),jpk) ) … … 654 660 IF (nn_ice_lim2(ib_bdy) .gt. 0) THEN 655 661 IF( nn_ice_lim2_dta(ib_bdy) .eq. 0 ) THEN 656 IF( nn_ice_lim2(ib_bdy) .eq. jp_frs ) THEN 657 ilen0(1:3) = nblen(1:3) 658 ELSE 659 ilen0(1:3) = nblenrim(1:3) 660 ENDIF 662 ilen0(1:3) = nblen(1:3) 661 663 ALLOCATE( dta_bdy(ib_bdy)%frld(ilen0(1)) ) 662 664 ALLOCATE( dta_bdy(ib_bdy)%hicif(ilen0(1)) ) -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdydyn3d.F90
r3294 r3367 19 19 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 20 20 USE in_out_manager ! 21 Use phycst 21 22 22 23 IMPLICIT NONE … … 24 25 25 26 PUBLIC bdy_dyn3d ! routine called by bdy_dyn 27 PUBLIC bdy_dyn3d_dmp ! routine called by step 26 28 27 29 !!---------------------------------------------------------------------- … … 55 57 CASE(jp_frs) 56 58 CALL bdy_dyn3d_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 59 CASE(2) 60 CALL bdy_dyn3d_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 61 CASE(3) 62 CALL bdy_dyn3d_zro( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 57 63 CASE DEFAULT 58 64 CALL ctl_stop( 'bdy_dyn3d : unrecognised option for open boundaries for baroclinic velocities' ) … … 62 68 END SUBROUTINE bdy_dyn3d 63 69 70 SUBROUTINE bdy_dyn3d_spe( idx, dta, kt ) 71 !!---------------------------------------------------------------------- 72 !! *** SUBROUTINE bdy_dyn3d_spe *** 73 !! 74 !! ** Purpose : - Apply a specified value for baroclinic velocities 75 !! at open boundaries. 76 !! 77 !!---------------------------------------------------------------------- 78 INTEGER :: kt 79 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 80 TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data 81 !! 82 INTEGER :: jb, jk ! dummy loop indices 83 INTEGER :: ii, ij, igrd ! local integers 84 REAL(wp) :: zwgt ! boundary weight 85 !!---------------------------------------------------------------------- 86 ! 87 IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_spe') 88 ! 89 igrd = 2 ! Relaxation of zonal velocity 90 DO jb = 1, idx%nblenrim(igrd) 91 DO jk = 1, jpkm1 92 ii = idx%nbi(jb,igrd) 93 ij = idx%nbj(jb,igrd) 94 ua(ii,ij,jk) = dta%u3d(jb,jk) * umask(ii,ij,jk) 95 END DO 96 END DO 97 ! 98 igrd = 3 ! Relaxation of meridional velocity 99 DO jb = 1, idx%nblenrim(igrd) 100 DO jk = 1, jpkm1 101 ii = idx%nbi(jb,igrd) 102 ij = idx%nbj(jb,igrd) 103 va(ii,ij,jk) = dta%v3d(jb,jk) * vmask(ii,ij,jk) 104 END DO 105 END DO 106 CALL lbc_lnk( ua, 'U', -1. ) ; CALL lbc_lnk( va, 'V', -1. ) ! Boundary points should be updated 107 ! 108 IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 109 110 IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_spe') 111 112 END SUBROUTINE bdy_dyn3d_spe 113 114 SUBROUTINE bdy_dyn3d_zro( idx, dta, kt ) 115 !!---------------------------------------------------------------------- 116 !! *** SUBROUTINE bdy_dyn3d_zro *** 117 !! 118 !! ** Purpose : - baroclinic velocities = 0. at open boundaries. 119 !! 120 !!---------------------------------------------------------------------- 121 INTEGER :: kt 122 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 123 TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data 124 !! 125 INTEGER :: ib, ik ! dummy loop indices 126 INTEGER :: ii, ij, igrd, zcoef ! local integers 127 REAL(wp) :: zwgt ! boundary weight 128 !!---------------------------------------------------------------------- 129 ! 130 IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_zro') 131 ! 132 igrd = 2 ! Everything is at T-points here 133 DO ib = 1, idx%nblenrim(igrd) 134 ii = idx%nbi(ib,igrd) 135 ij = idx%nbj(ib,igrd) 136 DO ik = 1, jpkm1 137 ua(ii,ij,ik) = 0._wp 138 END DO 139 END DO 140 141 igrd = 3 ! Everything is at T-points here 142 DO ib = 1, idx%nblenrim(igrd) 143 ii = idx%nbi(ib,igrd) 144 ij = idx%nbj(ib,igrd) 145 DO ik = 1, jpkm1 146 va(ii,ij,ik) = 0._wp 147 END DO 148 END DO 149 ! 150 CALL lbc_lnk( ua, 'U', -1. ) ; CALL lbc_lnk( va, 'V', -1. ) ! Boundary points should be updated 151 ! 152 IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 153 154 IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_zro') 155 156 END SUBROUTINE bdy_dyn3d_zro 157 64 158 SUBROUTINE bdy_dyn3d_frs( idx, dta, kt ) 65 159 !!---------------------------------------------------------------------- … … 111 205 END SUBROUTINE bdy_dyn3d_frs 112 206 207 SUBROUTINE bdy_dyn3d_dmp( kt ) 208 !!---------------------------------------------------------------------- 209 !! *** SUBROUTINE bdy_dyn3d_dmp *** 210 !! 211 !! ** Purpose : Apply damping for baroclinic velocities at open boundaries. 212 !! 213 !!---------------------------------------------------------------------- 214 INTEGER :: kt 215 !! 216 INTEGER :: jb, jk ! dummy loop indices 217 INTEGER :: ii, ij, igrd ! local integers 218 REAL(wp) :: zwgt ! boundary weight 219 INTEGER :: ib_bdy ! loop index 220 !!---------------------------------------------------------------------- 221 ! 222 IF( nn_timing == 1 ) CALL timing_start('bdy_dyn3d_dmp') 223 ! 224 DO ib_bdy=1, nb_bdy 225 IF ( ln_dyn3d_dmp(ib_bdy).and.nn_dyn3d(ib_bdy).gt.0 ) THEN 226 igrd = 2 ! Relaxation of zonal velocity 227 DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd) 228 ii = idx_bdy(ib_bdy)%nbi(jb,igrd) 229 ij = idx_bdy(ib_bdy)%nbj(jb,igrd) 230 zwgt = idx_bdy(ib_bdy)%nbw(jb,igrd) / ( rn_time_dmp(ib_bdy) * rday ) 231 DO jk = 1, jpkm1 232 ua(ii,ij,jk) = ( ua(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%u3d(jb,jk) - ub(ii,ij,jk) ) ) * umask(ii,ij,jk) 233 END DO 234 END DO 235 ! 236 igrd = 3 ! Relaxation of meridional velocity 237 DO jb = 1, idx_bdy(ib_bdy)%nblen(igrd) 238 ii = idx_bdy(ib_bdy)%nbi(jb,igrd) 239 ij = idx_bdy(ib_bdy)%nbj(jb,igrd) 240 zwgt = idx_bdy(ib_bdy)%nbw(jb,igrd) / ( rn_time_dmp(ib_bdy) * rday ) 241 DO jk = 1, jpkm1 242 va(ii,ij,jk) = ( va(ii,ij,jk) + zwgt * ( dta_bdy(ib_bdy)%v3d(jb,jk) - vb(ii,ij,jk) ) ) * vmask(ii,ij,jk) 243 END DO 244 END DO 245 ENDIF 246 ENDDO 247 ! 248 CALL lbc_lnk( ua, 'U', -1. ) ; CALL lbc_lnk( va, 'V', -1. ) ! Boundary points should be updated 249 ! 250 IF( nn_timing == 1 ) CALL timing_stop('bdy_dyn3d_dmp') 251 252 END SUBROUTINE bdy_dyn3d_dmp 113 253 114 254 #else -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90
r3298 r3367 80 80 & ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn2d_dta, & 81 81 & nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta, & 82 & ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, & 82 83 #if defined key_lim2 83 84 & nn_ice_lim2, nn_ice_lim2_dta, & … … 121 122 nn_tra(:) = 0 122 123 nn_tra_dta(:) = -1 ! uninitialised flag 124 ln_tra_dmp(:) = .false. 125 ln_dyn3d_dmp(:) = .false. 126 rn_time_dmp(:) = 1. 123 127 #if defined key_lim2 124 128 nn_ice_lim2(:) = 0 … … 135 139 ! Check and write out namelist parameters 136 140 ! ----------------------------------------- 137 141 142 ln_tra_dmp_tot=.false. 143 ln_dyn3d_dmp_tot=.false. 138 144 ! ! control prints 139 145 IF(lwp) WRITE(numout,*) ' nambdy' … … 178 184 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 179 185 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 186 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Specified value' 187 CASE( 3 ) ; IF(lwp) WRITE(numout,*) ' Zero baroclinic velocities (runoff case)' 180 188 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_dyn3d' ) 181 189 END SELECT … … 189 197 IF(lwp) WRITE(numout,*) 190 198 199 IF ( ln_dyn3d_dmp(ib_bdy) ) THEN 200 IF ( nn_dyn3d(ib_bdy).EQ.0 ) THEN 201 IF(lwp) WRITE(numout,*) 'No open boundary condition for the baroclinic velocitues : we force ln_dyn3d_dmp=.false.' 202 ln_dyn3d_dmp(ib_bdy)=.false. 203 ELSE 204 IF(lwp) WRITE(numout,*) ' Damping of the baroclinic velocities along the boundaries' 205 IF(lwp) WRITE(numout,*) ' Time damping :',rn_time_dmp(ib_bdy),' days' 206 ln_dyn3d_dmp_tot=.true. 207 ENDIF 208 ENDIF 209 IF(lwp) WRITE(numout,*) 210 191 211 IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity: ' 192 212 SELECT CASE( nn_tra(ib_bdy) ) 193 213 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' no open boundary condition' 194 214 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Flow Relaxation Scheme' 215 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Specified value' 216 CASE( 3 ) ; IF(lwp) WRITE(numout,*) ' Neumann conditions' 217 CASE( 4 ) ; IF(lwp) WRITE(numout,*) ' Runoff conditions : Neumann for T and specified to 0.1 for salinity' 195 218 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_tra' ) 196 219 END SELECT … … 201 224 CASE DEFAULT ; CALL ctl_stop( 'nn_tra_dta must be 0 or 1' ) 202 225 END SELECT 226 ENDIF 227 IF(lwp) WRITE(numout,*) 228 229 IF ( ln_tra_dmp(ib_bdy) ) THEN 230 IF ( nn_tra(ib_bdy).EQ.0 ) THEN 231 IF(lwp) WRITE(numout,*) 'No open boundary condition for the tracer : we force ln_tra_dmp=.false.' 232 ln_tra_dmp(ib_bdy)=.false. 233 ELSE 234 IF(lwp) WRITE(numout,*) ' Damping of the scalars along the boundaries' 235 IF(lwp) WRITE(numout,*) ' Time dumping :',rn_time_dmp(ib_bdy),' days' 236 ln_tra_dmp_tot=.true. 237 ENDIF 203 238 ENDIF 204 239 IF(lwp) WRITE(numout,*) … … 247 282 ! --------------------------------------------- 248 283 REWIND( numnam ) 284 jpbdta = 1 285 nblendta(:,:) = 0 249 286 DO ib_bdy = 1, nb_bdy 250 287 251 jpbdta = 1252 288 IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters 253 289 … … 282 318 ENDIF 283 319 284 nblendta(:,ib_bdy) = 0285 320 DO iseg = 1, nbdysege 286 321 igrd = 1 … … 317 352 318 353 nblendta(:,ib_bdy) = nblendta(:,ib_bdy) * nn_rimwidth(ib_bdy) 319 jpbdta = MAXVAL(nblendta(:,ib_bdy))320 321 354 322 355 ELSE ! Read size of arrays in boundary coordinates file. … … 324 357 325 358 CALL iom_open( cn_coords_file(ib_bdy), inum ) 326 jpbdta = 1327 359 DO igrd = 1, jpbgrd 328 360 id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz ) 329 361 nblendta(igrd,ib_bdy) = kdimsz(1) 330 jpbdta = MAX(jpbdta, kdimsz(1))331 362 ENDDO 332 363 … … 334 365 335 366 ENDDO ! ib_bdy 367 368 jpbdta = MAXVAL(nblendta(:,:)) 336 369 337 370 ! Allocate arrays … … 610 643 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 611 644 nbr => idx_bdy(ib_bdy)%nbr(ib,igrd) 612 idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 ) ! tanh formulation613 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2! quadratic614 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth) ! linear645 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 ) ! tanh formulation 646 idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2. ! quadratic 647 ! idx_bdy(ib_bdy)%nbw(ib,igrd) = FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth) ! linear 615 648 END DO 616 649 END DO … … 753 786 END IF 754 787 END DO 755 788 756 789 IF( icount /= 0 ) THEN 757 790 IF(lwp) WRITE(numout,*) -
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 -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/BDY/bdytra.F90
r3294 r3367 22 22 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 23 23 USE in_out_manager ! I/O manager 24 USE phycst 24 25 25 26 IMPLICIT NONE … … 27 28 28 29 PUBLIC bdy_tra ! routine called in tranxt.F90 30 PUBLIC bdy_tra_dmp ! routine called in tranxt.F90 29 31 30 32 !!---------------------------------------------------------------------- … … 53 55 CASE(jp_frs) 54 56 CALL bdy_tra_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 57 CASE(2) 58 CALL bdy_tra_spe( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 59 CASE(3) 60 CALL bdy_tra_nmn( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 61 CASE(4) 62 CALL bdy_tra_rnf( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt ) 55 63 CASE DEFAULT 56 64 CALL ctl_stop( 'bdy_tra : unrecognised option for open boundaries for T and S' ) … … 97 105 ! 98 106 END SUBROUTINE bdy_tra_frs 99 107 108 SUBROUTINE bdy_tra_spe( idx, dta, kt ) 109 !!---------------------------------------------------------------------- 110 !! *** SUBROUTINE bdy_tra_frs *** 111 !! 112 !! ** Purpose : Apply a specified value for tracers at open boundaries. 113 !! 114 !!---------------------------------------------------------------------- 115 INTEGER, INTENT(in) :: kt 116 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 117 TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data 118 !! 119 REAL(wp) :: zwgt ! boundary weight 120 INTEGER :: ib, ik, igrd ! dummy loop indices 121 INTEGER :: ii, ij ! 2D addresses 122 !!---------------------------------------------------------------------- 123 ! 124 IF( nn_timing == 1 ) CALL timing_start('bdy_tra_spe') 125 ! 126 igrd = 1 ! Everything is at T-points here 127 DO ib = 1, idx%nblenrim(igrd) 128 ii = idx%nbi(ib,igrd) 129 ij = idx%nbj(ib,igrd) 130 DO ik = 1, jpkm1 131 tsa(ii,ij,ik,jp_tem) = dta%tem(ib,ik) * tmask(ii,ij,ik) 132 tsa(ii,ij,ik,jp_sal) = dta%sal(ib,ik) * tmask(ii,ij,ik) 133 END DO 134 END DO 135 ! 136 CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) ! Boundary points should be updated 137 ! 138 IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 139 ! 140 IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_spe') 141 ! 142 END SUBROUTINE bdy_tra_spe 143 144 SUBROUTINE bdy_tra_nmn( idx, dta, kt ) 145 !!---------------------------------------------------------------------- 146 !! *** SUBROUTINE bdy_tra_nmn *** 147 !! 148 !! ** Purpose : Duplicate the value for tracers at open boundaries. 149 !! 150 !!---------------------------------------------------------------------- 151 INTEGER, INTENT(in) :: kt 152 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 153 TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data 154 !! 155 REAL(wp) :: zwgt ! boundary weight 156 INTEGER :: ib, ik, igrd ! dummy loop indices 157 INTEGER :: ii, ij,zcoef, zcoef1,zcoef2, ip, jp ! 2D addresses 158 !!---------------------------------------------------------------------- 159 ! 160 IF( nn_timing == 1 ) CALL timing_start('bdy_tra_nmn') 161 ! 162 igrd = 1 ! Everything is at T-points here 163 DO ib = 1, idx%nblenrim(igrd) 164 ii = idx%nbi(ib,igrd) 165 ij = idx%nbj(ib,igrd) 166 DO ik = 1, jpkm1 167 ! search the sense of the gradient 168 zcoef1 = bdytmask(ii-1,ij ) + bdytmask(ii+1,ij ) 169 zcoef2 = bdytmask(ii ,ij-1) + bdytmask(ii ,ij+1) 170 IF ( zcoef1+zcoef2 == 0) THEN 171 ! corner 172 zcoef = tmask(ii-1,ij,ik) + tmask(ii+1,ij,ik) + tmask(ii,ij-1,ik) + tmask(ii,ij+1,ik) 173 tsa(ii,ij,ik,jp_tem) = tsa(ii-1,ij ,ik,jp_tem) * tmask(ii-1,ij ,ik) + & 174 & tsa(ii+1,ij ,ik,jp_tem) * tmask(ii+1,ij ,ik) + & 175 & tsa(ii ,ij-1,ik,jp_tem) * tmask(ii ,ij-1,ik) + & 176 & tsa(ii ,ij+1,ik,jp_tem) * tmask(ii ,ij+1,ik) 177 tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) / MAX( 1, zcoef) ) * tmask(ii,ij,ik) 178 tsa(ii,ij,ik,jp_sal) = tsa(ii-1,ij ,ik,jp_sal) * tmask(ii-1,ij ,ik) + & 179 & tsa(ii+1,ij ,ik,jp_sal) * tmask(ii+1,ij ,ik) + & 180 & tsa(ii ,ij-1,ik,jp_sal) * tmask(ii ,ij-1,ik) + & 181 & tsa(ii ,ij+1,ik,jp_sal) * tmask(ii ,ij+1,ik) 182 tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) / MAX( 1, zcoef) ) * tmask(ii,ij,ik) 183 ELSE 184 ip = bdytmask(ii+1,ij ) - bdytmask(ii-1,ij ) 185 jp = bdytmask(ii ,ij+1) - bdytmask(ii ,ij-1) 186 tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii+ip,ij+jp,ik) 187 tsa(ii,ij,ik,jp_sal) = tsa(ii+ip,ij+jp,ik,jp_sal) * tmask(ii+ip,ij+jp,ik) 188 ENDIF 189 END DO 190 END DO 191 ! 192 CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) ! Boundary points should be updated 193 ! 194 IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 195 ! 196 IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_nmn') 197 ! 198 END SUBROUTINE bdy_tra_nmn 199 200 SUBROUTINE bdy_tra_rnf( idx, dta, kt ) 201 !!---------------------------------------------------------------------- 202 !! *** SUBROUTINE bdy_tra_rnf *** 203 !! 204 !! ** Purpose : Apply the runoff values for tracers at open boundaries: 205 !! - specified to 0.1 PSU for the salinity 206 !! - duplicate the value for the temperature 207 !! 208 !!---------------------------------------------------------------------- 209 INTEGER, INTENT(in) :: kt 210 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 211 TYPE(OBC_DATA), INTENT(in) :: dta ! OBC external data 212 !! 213 REAL(wp) :: zwgt ! boundary weight 214 INTEGER :: ib, ik, igrd ! dummy loop indices 215 INTEGER :: ii, ij, ip, jp ! 2D addresses 216 !!---------------------------------------------------------------------- 217 ! 218 IF( nn_timing == 1 ) CALL timing_start('bdy_tra_rnf') 219 ! 220 igrd = 1 ! Everything is at T-points here 221 DO ib = 1, idx%nblenrim(igrd) 222 ii = idx%nbi(ib,igrd) 223 ij = idx%nbj(ib,igrd) 224 DO ik = 1, jpkm1 225 ip = bdytmask(ii+1,ij ) - bdytmask(ii-1,ij ) 226 jp = bdytmask(ii ,ij+1) - bdytmask(ii ,ij-1) 227 tsa(ii,ij,ik,jp_tem) = tsa(ii+ip,ij+jp,ik,jp_tem) * tmask(ii,ij,ik) 228 tsa(ii,ij,ik,jp_sal) = 0.1 * tmask(ii,ij,ik) 229 END DO 230 END DO 231 ! 232 CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) ! Boundary points should be updated 233 ! 234 IF( kt .eq. nit000 ) CLOSE( unit = 102 ) 235 ! 236 IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_rnf') 237 ! 238 END SUBROUTINE bdy_tra_rnf 239 240 SUBROUTINE bdy_tra_dmp( kt ) 241 !!---------------------------------------------------------------------- 242 !! *** SUBROUTINE bdy_tra_dmp *** 243 !! 244 !! ** Purpose : Apply damping for tracers at open boundaries. 245 !! 246 !!---------------------------------------------------------------------- 247 INTEGER, INTENT(in) :: kt 248 !! 249 REAL(wp) :: zwgt ! boundary weight 250 REAL(wp) :: zta, zsa, ztime 251 INTEGER :: ib, ik, igrd ! dummy loop indices 252 INTEGER :: ii, ij ! 2D addresses 253 INTEGER :: ib_bdy ! Loop index 254 !!---------------------------------------------------------------------- 255 ! 256 IF( nn_timing == 1 ) CALL timing_start('bdy_tra_dmp') 257 ! 258 DO ib_bdy=1, nb_bdy 259 IF ( ln_tra_dmp(ib_bdy) ) THEN 260 igrd = 1 ! Everything is at T-points here 261 DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd) 262 ii = idx_bdy(ib_bdy)%nbi(ib,igrd) 263 ij = idx_bdy(ib_bdy)%nbj(ib,igrd) 264 zwgt = idx_bdy(ib_bdy)%nbw(ib,igrd) / ( rn_time_dmp(ib_bdy) * rday ) 265 DO ik = 1, jpkm1 266 tsa(ii,ij,ik,jp_tem) = ( tsa(ii,ij,ik,jp_tem) + zwgt * ( dta_bdy(ib_bdy)%tem(ib,ik) - tsb(ii,ij,ik,jp_tem) ) ) * tmask(ii,ij,ik) 267 tsa(ii,ij,ik,jp_sal) = ( tsa(ii,ij,ik,jp_sal) + zwgt * ( dta_bdy(ib_bdy)%sal(ib,ik) - tsb(ii,ij,ik,jp_sal) ) ) * tmask(ii,ij,ik) 268 END DO 269 END DO 270 ENDIF 271 ENDDO 272 ! 273 CALL lbc_lnk( tsa(:,:,:,jp_tem), 'T', 1. ) ; CALL lbc_lnk( tsa(:,:,:,jp_sal), 'T', 1. ) ! Boundary points should be updated 274 ! 275 IF( nn_timing == 1 ) CALL timing_stop('bdy_tra_dmp') 276 ! 277 END SUBROUTINE bdy_tra_dmp 278 100 279 #else 101 280 !!---------------------------------------------------------------------- -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r3294 r3367 403 403 IF( lk_obc ) CALL obc_dta_bt ( kt, jn ) 404 404 IF( lk_bdy ) CALL bdy_dta ( kt, jit=jn, time_offset=+1 ) 405 IF ( ln_tide_pot ) CALL upd_tide( kt, jn )405 IF ( ln_tide_pot .AND. lk_tide) CALL upd_tide( kt, jn ) 406 406 407 407 ! !* after ssh_e … … 453 453 ENDIF 454 454 ! add tidal astronomical forcing 455 IF ( ln_tide_pot ) THEN455 IF ( ln_tide_pot .AND. lk_tide ) THEN 456 456 zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 457 457 zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) … … 503 503 ENDIF 504 504 ! add tidal astronomical forcing 505 IF ( ln_tide_pot ) THEN505 IF ( ln_tide_pot .AND. lk_tide ) THEN 506 506 zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 507 507 zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) … … 550 550 ENDIF 551 551 ! add tidal astronomical forcing 552 IF ( ln_tide_pot ) THEN552 IF ( ln_tide_pot .AND. lk_tide ) THEN 553 553 zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 554 554 zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90
r3294 r3367 66 66 !! 67 67 INTEGER :: ierror ! local integer 68 REAL(wp) :: zpref ! local scalar69 68 !! 70 69 CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files 71 70 TYPE(FLD_N) :: sn_apr ! informations about the fields to be read 72 71 !! 73 NAMELIST/namsbc_apr/ cn_dir, sn_apr, ln_ref_apr 72 NAMELIST/namsbc_apr/ cn_dir, sn_apr, ln_ref_apr, rpref, ln_apr_obc 74 73 !!---------------------------------------------------------------------- 75 74 ! … … 113 112 ! 114 113 ! !* control check 115 IF( ln_apr_obc ) & 116 CALL ctl_stop( 'sbc_apr: inverse barometer added to OBC ssh data not yet implemented ' ) 117 IF( ln_apr_obc .AND. .NOT. lk_obc ) & 118 CALL ctl_stop( 'sbc_apr: add inverse barometer to OBC requires to use key_obc' ) 114 IF ( ln_apr_obc ) THEN 115 IF(lwp) WRITE(numout,*) ' Inverse barometer added to OBC ssh data' 116 ENDIF 119 117 IF( ( ln_apr_obc ) .AND. .NOT. lk_dynspg_ts ) & 120 118 CALL ctl_stop( 'sbc_apr: use inverse barometer ssh at open boundary ONLY possible with time-splitting' ) -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/SBC/sbctide.F90
r3294 r3367 14 14 USE daymod 15 15 USE dynspg_oce 16 USE tide _mod16 USE tideini 17 17 USE iom 18 18 … … 21 21 22 22 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: pot_astro 23 LOGICAL, PUBLIC :: ln_tide_pot = .false. 23 24 24 #if defined key_tide 25 25 26 26 LOGICAL, PUBLIC, PARAMETER :: lk_tide = .TRUE. 27 28 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: omega_tide29 30 REAL(wp), ALLOCATABLE, DIMENSION(:) :: &31 v0tide, &32 utide, &33 ftide34 35 27 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: amp_pot,phi_pot 36 37 INTEGER, PUBLIC :: nb_harmo38 INTEGER, PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntide39 INTEGER, PUBLIC :: nn_tide, kt_tide40 41 28 !!--------------------------------------------------------------------------------- 42 29 !! OPA 9.0 , LODYC-IPSL (2003) … … 51 38 !! * Arguments 52 39 INTEGER, INTENT( in ) :: kt ! ocean time-step 53 !! * Local declarations54 INTEGER :: jk,ji55 CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname56 40 !!---------------------------------------------------------------------- 57 41 58 NAMELIST/nam_tide/ln_tide_pot,nb_harmo,clname,nn_tide42 IF ( kt == nit000 .AND. .NOT. lk_dynspg_ts ) CALL ctl_stop( 'STOP', 'sbc_tide : tidal potential use only with time splitting' ) 59 43 60 IF ( kt == nit000 ) THEN 44 IF ( nsec_day == NINT(0.5 * rdttra(1)) ) THEN 45 ! 46 kt_tide = kt 61 47 62 IF( .NOT. lk_dynspg_ts ) CALL ctl_stop( 'STOP', 'sbc_tide : tidal potential use only with time splitting' ) 48 IF(lwp) THEN 49 WRITE(numout,*) 50 WRITE(numout,*) 'sbc_tide : (re)Initialization of the tidal potential at kt=',kt 51 WRITE(numout,*) '~~~~~~~ ' 52 ENDIF 63 53 64 ! Read Namelist nam_tide 54 IF(lwp) THEN 55 IF ( kt == nit000 ) WRITE(numout,*) 'Apply astronomical potential : ln_tide_pot =', ln_tide_pot 56 CALL flush(numout) 57 ENDIF 65 58 66 nn_tide=INT(rday/rdt) 59 IF ( kt == nit000 ) ALLOCATE(amp_pot(jpi,jpj,nb_harmo)) 60 IF ( kt == nit000 ) ALLOCATE(phi_pot(jpi,jpj,nb_harmo)) 61 IF ( kt == nit000 ) ALLOCATE(pot_astro(jpi,jpj)) 67 62 68 CALL tide_init_Wave 63 amp_pot(:,:,:) = 0.e0 64 phi_pot(:,:,:) = 0.e0 65 pot_astro(:,:) = 0.e0 69 66 70 REWIND ( numnam ) 71 READ ( numnam, nam_tide ) 72 73 IF(lwp) THEN 74 WRITE(numout,*) 75 WRITE(numout,*) 'sbc_tide : Initialization of the tidal components' 76 WRITE(numout,*) '~~~~~~~ ' 67 IF ( ln_tide_pot ) CALL tide_init_potential 68 ! 77 69 ENDIF 78 79 IF(lwp) THEN80 WRITE(numout,*) ' Namelist nam_tide'81 WRITE(numout,*) ' Apply astronomical potential : ln_tide_pot =', ln_tide_pot82 WRITE(numout,*) ' nb_harmo = ', nb_harmo83 CALL flush(numout)84 ENDIF85 86 ALLOCATE(ntide (nb_harmo))87 DO jk=1,nb_harmo88 DO ji=1,jpmax_harmo89 IF (TRIM(clname(jk)) .eq. Wave(ji)%cname_tide) THEN90 ntide(jk) = ji91 EXIT92 END IF93 END DO94 END DO95 ALLOCATE(omega_tide(nb_harmo))96 ALLOCATE(v0tide (nb_harmo))97 ALLOCATE(utide (nb_harmo))98 ALLOCATE(ftide (nb_harmo))99 ALLOCATE(amp_pot(jpi,jpj,nb_harmo))100 ALLOCATE(phi_pot(jpi,jpj,nb_harmo))101 ALLOCATE(pot_astro(jpi,jpj))102 ENDIF103 104 IF ( MOD( kt - 1, nn_tide ) == 0 ) THEN105 kt_tide = kt106 CALL tide_harmo(omega_tide, v0tide, utide, ftide, ntide, nb_harmo)107 ENDIF108 109 amp_pot(:,:,:) = 0.e0110 phi_pot(:,:,:) = 0.e0111 pot_astro(:,:) = 0.e0112 113 IF (ln_tide_pot ) CALL tide_init_potential114 70 115 71 END SUBROUTINE sbc_tide -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
r3294 r3367 46 46 USE mppini ! shared/distributed memory setting (mpp_init routine) 47 47 USE domain ! domain initialization (dom_init routine) 48 USE tideini ! tidal components initialization (tide_ini routine) 48 49 USE obcini ! open boundary cond. initialization (obc_ini routine) 49 50 USE bdyini ! open boundary cond. initialization (bdy_init routine) 50 51 USE bdydta ! open boundary cond. initialization (bdy_dta_init routine) 51 USE bdytides ! open boundary cond. initialization ( tide_init routine)52 USE bdytides ! open boundary cond. initialization (bdytide_init routine) 52 53 USE istate ! initial state setting (istate_init routine) 53 54 USE ldfdyn ! lateral viscosity setting (ldfdyn_init routine) … … 314 315 315 316 IF( lk_obc ) CALL obc_init ! Open boundaries 317 318 CALL istate_init ! ocean initial state (Dynamics and tracers) 319 320 CALL tide_init( nit000 ) ! Initialisation of the tidal harmonics 321 316 322 IF( lk_bdy ) CALL bdy_init ! Open boundaries initialisation 317 323 IF( lk_bdy ) CALL bdy_dta_init ! Open boundaries initialisation of external data arrays 318 IF( lk_bdy ) CALL tide_init! Open boundaries initialisation of tidal harmonic forcing324 IF( lk_bdy ) CALL bdytide_init ! Open boundaries initialisation of tidal harmonic forcing 319 325 320 326 CALL flush(numout) 321 327 CALL dyn_nept_init ! simplified form of Neptune effect 322 328 CALL flush(numout) 323 324 CALL istate_init ! ocean initial state (Dynamics and tracers)325 329 326 330 ! ! Ocean physics -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/step.F90
r3294 r3367 95 95 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 96 96 CALL sbc ( kstp ) ! Sea Boundary Condition (including sea-ice) 97 IF( kstp /= nit000 ) CALL tide_init ( kstp ) 97 98 IF( lk_tide ) CALL sbc_tide( kstp ) 98 99 IF( lk_obc ) CALL obc_dta( kstp ) ! update dynamic and tracer data at open boundaries … … 187 188 IF( lk_trabbl ) CALL tra_bbl ( kstp ) ! advective (and/or diffusive) bottom boundary layer scheme 188 189 IF( ln_tradmp ) CALL tra_dmp ( kstp ) ! internal damping trends 190 IF( lk_bdy .AND. & 191 & ln_tra_dmp_tot) CALL bdy_tra_dmp( kstp ) ! bdy damping trends 189 192 CALL tra_adv ( kstp ) ! horizontal & vertical advection 190 193 IF( lk_zdfkpp ) CALL tra_kpp ( kstp ) ! KPP non-local tracer fluxes … … 219 222 & ln_dyninc ) CALL dyn_asm_inc( kstp ) ! apply dynamics assimilation increment 220 223 IF( ln_neptsimp ) CALL dyn_nept_cor( kstp ) ! subtract Neptune velocities (simplified) 224 IF( lk_bdy .AND. & 225 & ln_dyn3d_dmp_tot) CALL bdy_dyn3d_dmp(kstp ) ! bdy damping trends 221 226 CALL dyn_adv( kstp ) ! advection (vector or flux form) 222 227 CALL dyn_vor( kstp ) ! vorticity term including Coriolis -
branches/2012/dev_r3327_MERCATOR1_BDY/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
r3294 r3367 55 55 56 56 USE bdy_par ! for lk_bdy 57 USE bdy_oce ! for dmp logical 57 58 USE bdydta ! open boundary condition data (bdy_dta routine) 59 USE bdytra ! bdy cond. for tracers (bdy_tra routine) 60 USE bdydyn3d ! bdy cond. for baroclinic vel. (bdy_dyn3d routine) 58 61 59 62 USE sshwzv ! vertical velocity and ssh (ssh_wzv routine)
Note: See TracChangeset
for help on using the changeset viewer.