Changeset 15788
- Timestamp:
- 2022-04-22T12:12:35+02:00 (2 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0.4_CO9_wadcpl/src/OCE
- Files:
-
- 14 edited
- 2 copied
-
BDY/bdyini.F90 (modified) (1 diff)
-
DIA/diaharm.F90 (modified) (14 diffs)
-
DIA/diaharm_fast.F90 (copied) (copied from NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/DIA/diaharm_fast.F90)
-
DOM/dtatsd.F90 (modified) (8 diffs)
-
DYN/dynnxt.F90 (modified) (7 diffs)
-
DYN/dynspg.F90 (modified) (3 diffs)
-
SBC/fldread.F90 (modified) (1 diff)
-
SBC/sbctide.F90 (modified) (3 diffs)
-
SBC/tide_FES14.h90 (copied) (copied from NEMO/branches/UKMO/r14075_India_uncoupled/src/OCE/SBC/tide_FES14.h90)
-
SBC/tide_mod.F90 (modified) (3 diffs)
-
SBC/tideini.F90 (modified) (3 diffs)
-
USR/usrdef_istate.F90 (modified) (1 diff)
-
par_oce.F90 (modified) (1 diff)
-
step.F90 (modified) (1 diff)
-
step_oce.F90 (modified) (1 diff)
-
stpctl.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.4_CO9_wadcpl/src/OCE/BDY/bdyini.F90
r15784 r15788 196 196 ! Check and write out namelist parameters 197 197 ! ----------------------------------------- 198 IF( jperio /= 0 ) CALL ctl_stop( 'bdy_segs: Cyclic or symmetric,', &199 & ' and general open boundary condition are not compatible' )198 ! IF( jperio /= 0 ) CALL ctl_stop( 'bdy_segs: Cyclic or symmetric,', & 199 ! & ' and general open boundary condition are not compatible' ) 200 200 201 201 IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy -
NEMO/branches/UKMO/NEMO_4.0.4_CO9_wadcpl/src/OCE/DIA/diaharm.F90
r14075 r15788 5 5 !!====================================================================== 6 6 !! History : 3.1 ! 2007 (O. Le Galloudec, J. Chanut) Original code 7 !! 8 !! NB: 2017-12 : add 3D harmonic analysis of velocities 9 !! integration of Maria Luneva's development 10 !! 'key_3Ddiaharm 7 11 !!---------------------------------------------------------------------- 8 12 USE oce ! ocean dynamics and tracers variables … … 13 17 USE sbctide ! Tidal forcing or not 14 18 ! 19 # if defined key_3Ddiaharm 20 USE zdf_oce 21 #endif 22 ! 15 23 USE in_out_manager ! I/O units 16 24 USE iom ! I/0 library … … 33 41 INTEGER :: nb_ana ! Number of harmonics to analyse 34 42 35 INTEGER , ALLOCATABLE, DIMENSION(:) :: name 36 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ana_temp 37 REAL(wp), ALLOCATABLE, DIMENSION(:) :: ana_freq, ut, vt, ft 38 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: out_eta, out_u, out_v 43 INTEGER , ALLOCATABLE, DIMENSION(:) :: name 44 REAL(wp), ALLOCATABLE, DIMENSION(:) :: ana_freq, ut, vt, ft 45 # if defined key_3Ddiaharm 46 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:,:) :: ana_temp 47 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: out_eta, out_u, out_v, out_w, out_dzi 48 # else 49 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ana_temp 50 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: out_eta, out_u, out_v 51 # endif 39 52 40 53 INTEGER :: ninco, nsparse … … 76 89 WRITE(numout,*) 77 90 WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization' 91 # if defined key_3Ddiaharm 92 WRITE(numout,*) ' - 3D harmonic analysis of currents activated (key_3Ddiaharm)' 93 #endif 78 94 WRITE(numout,*) '~~~~~~~~~~~~~ ' 79 95 ENDIF … … 155 171 ! Initialize temporary arrays: 156 172 ! ---------------------------- 157 ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) ) 158 ana_temp(:,:,:,:) = 0._wp 173 # if defined key_3Ddiaharm 174 ALLOCATE( ana_temp( jpi, jpj, 2*nb_ana, 5, jpk ) ) 175 ana_temp(:,:,:,:,:) = 0._wp 176 # else 177 ALLOCATE( ana_temp( jpi, jpj, 2*nb_ana, 3 ) ) 178 ana_temp(:,:,:,: ) = 0._wp 179 #endif 159 180 160 181 ENDIF … … 175 196 ! 176 197 INTEGER :: ji, jj, jh, jc, nhc 198 # if defined key_3Ddiaharm 199 INTEGER :: jk 200 # endif 177 201 REAL(wp) :: ztime, ztemp 178 202 !!-------------------------------------------------------------------- … … 190 214 & +(1.-MOD(jc,2))* ft(jh) *SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh))) 191 215 ! 216 ! ssh, ub, vb are stored at the last level of 5d array 192 217 DO jj = 2, jpjm1 193 218 DO ji = 2, jpim1 194 ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp * sshn(ji,jj) * ssmask (ji,jj) ! elevation 195 ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp * un_b(ji,jj) * ssumask(ji,jj) ! u-vel 196 ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp * vn_b(ji,jj) * ssvmask(ji,jj) ! v-vel 219 ! Elevation and currents 220 # if defined key_3Ddiaharm 221 ana_temp(ji,jj,nhc,1,jpk) = ana_temp(ji,jj,nhc,1,jpk) + ztemp*sshn(ji,jj)*ssmask (ji,jj) 222 ana_temp(ji,jj,nhc,2,jpk) = ana_temp(ji,jj,nhc,2,jpk) + ztemp*un_b(ji,jj)*ssumask(ji,jj) 223 ana_temp(ji,jj,nhc,3,jpk) = ana_temp(ji,jj,nhc,3,jpk) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj) 224 225 ana_temp(ji,jj,nhc,5,jpk) = ana_temp(ji,jj,nhc,5,jpk) & 226 & + ztemp*bfrva(ji,jj)*vn(ji,jj,mbkv(ji,jj))*ssvmask(ji,jj) 227 ana_temp(ji,jj,nhc,4,jpk) = ana_temp(ji,jj,nhc,4,jpk) & 228 & + ztemp*bfrua(ji,jj)*un(ji,jj,mbku(ji,jj))*ssumask(ji,jj) 229 # else 230 ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztemp*sshn(ji,jj)*ssmask (ji,jj) 231 ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztemp*un_b(ji,jj)*ssumask(ji,jj) 232 ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztemp*vn_b(ji,jj)*ssvmask(ji,jj) 233 # endif 197 234 END DO 198 235 END DO 236 ! 237 # if defined key_3Ddiaharm 238 ! 3d velocity and density: 239 DO jk=1,jpk-1 240 DO jj = 1,jpj 241 DO ji = 1,jpi 242 ! density and velocity 243 ana_temp(ji,jj,nhc,1,jk) = ana_temp(ji,jj,nhc,1,jk) + ztemp*rhd(ji,jj,jk) 244 ana_temp(ji,jj,nhc,2,jk) = ana_temp(ji,jj,nhc,2,jk) + ztemp*(un(ji,jj,jk)-un_b(ji,jj)) & 245 & *umask(ji,jj,jk) 246 ana_temp(ji,jj,nhc,3,jk) = ana_temp(ji,jj,nhc,3,jk) + ztemp*(vn(ji,jj,jk)-vn_b(ji,jj)) & 247 & *vmask(ji,jj,jk) 248 ana_temp(ji,jj,nhc,4,jk) = ana_temp(ji,jj,nhc,4,jk) + ztemp*wn(ji,jj,jk) 249 250 ana_temp(ji,jj,nhc,5,jk) = ana_temp(ji,jj,nhc,5,jk) - 0.5*grav*ztemp*(rhd(ji,jj,jk)+rhd(ji,jj,jk+1))/max(rn2(ji,jj,jk),1.e-8_wp) 251 END DO 252 END DO 253 ENDDO 254 # endif 199 255 END DO 200 256 END DO … … 218 274 !!-------------------------------------------------------------------- 219 275 INTEGER :: ji, jj, jh, jc, jn, nhan 276 # if defined key_3Ddiaharm 277 INTEGER :: jk 278 # endif 220 279 INTEGER :: ksp, kun, keq 221 280 REAL(wp) :: ztime, ztime_ini, ztime_end, z1_han … … 226 285 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 227 286 228 ALLOCATE( out_eta(jpi,jpj,2*nb_ana), out_u(jpi,jpj,2*nb_ana), out_v(jpi,jpj,2*nb_ana) )229 230 287 ztime_ini = nit000_han*rdt ! Initial time in seconds at the beginning of analysis 231 288 ztime_end = nitend_han*rdt ! Final time in seconds at the end of analysis … … 233 290 z1_han = 1._wp / REAL(nhan-1) 234 291 292 # if defined key_3Ddiaharm 293 ALLOCATE( out_eta(jpi,jpj,jpk,2*nb_ana), & 294 & out_u (jpi,jpj,jpk,2*nb_ana), & 295 & out_v (jpi,jpj,jpk,2*nb_ana), & 296 & out_w (jpi,jpj,jpk,2*nb_ana), & 297 & out_dzi(jpi,jpj,jpk,2*nb_ana) ) 298 # else 299 ALLOCATE( out_eta(jpi,jpj,2*nb_ana), & 300 & out_u (jpi,jpj,2*nb_ana), & 301 & out_v (jpi,jpj,2*nb_ana) ) 302 # endif 303 304 IF(lwp) WRITE(numout,*) 'ANA F OLD', ft 305 IF(lwp) WRITE(numout,*) 'ANA U OLD', ut 306 IF(lwp) WRITE(numout,*) 'ANA V OLD', vt 307 235 308 ninco = 2*nb_ana 236 309 … … 260 333 CALL SUR_DETERMINE_INIT 261 334 262 ! Elevation: 335 ! Density and Elevation: 336 # if defined key_3Ddiaharm 337 DO jk=1,jpk 338 # endif 263 339 DO jj = 2, jpjm1 264 340 DO ji = 2, jpim1 265 341 266 342 ! Fill input array 343 # if defined key_3Ddiaharm 344 ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,1,jk) 345 # else 267 346 ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,1) 347 # endif 268 348 CALL SUR_DETERMINE 269 349 270 350 ! Fill output array 271 351 DO jh = 1, nb_ana 272 out_eta(ji,jj,jh ) = ztmp7((jh-1)*2+1) * ssmask(ji,jj) 273 out_eta(ji,jj,jh+nb_ana) = -ztmp7((jh-1)*2+2) * ssmask(ji,jj) 352 # if defined key_3Ddiaharm 353 out_eta(ji,jj,jk,jh ) = ztmp7((jh-1)*2+1) * ssmask(ji,jj) 354 out_eta(ji,jj,jk,jh+nb_ana) = -ztmp7((jh-1)*2+2) * ssmask(ji,jj) 355 # else 356 out_eta(ji,jj, jh ) = ztmp7((jh-1)*2+1) * ssmask(ji,jj) 357 out_eta(ji,jj, jh+nb_ana) = -ztmp7((jh-1)*2+2) * ssmask(ji,jj) 358 # endif 274 359 END DO 275 360 END DO … … 281 366 282 367 ! Fill input array 368 # if defined key_3Ddiaharm 369 ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,2,jk) 370 # else 283 371 ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,2) 372 # endif 284 373 CALL SUR_DETERMINE 285 374 286 375 ! Fill output array 287 376 DO jh = 1, nb_ana 288 out_u(ji,jj, jh) = ztmp7((jh-1)*2+1) * ssumask(ji,jj) 289 out_u(ji,jj,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssumask(ji,jj) 377 # if defined key_3Ddiaharm 378 out_u(ji,jj,jk, jh) = ztmp7((jh-1)*2+1) * ssumask(ji,jj) 379 out_u(ji,jj,jk,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssumask(ji,jj) 380 # else 381 out_u(ji,jj, jh) = ztmp7((jh-1)*2+1) * ssumask(ji,jj) 382 out_u(ji,jj, nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssumask(ji,jj) 383 # endif 384 290 385 END DO 291 386 … … 298 393 299 394 ! Fill input array 395 # if defined key_3Ddiaharm 396 ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,3,jk) 397 # else 300 398 ztmp4(1:nb_ana*2) = ana_temp(ji,jj,1:nb_ana*2,3) 399 # endif 301 400 CALL SUR_DETERMINE 302 401 303 402 ! Fill output array 304 403 DO jh = 1, nb_ana 305 out_v(ji,jj, jh) = ztmp7((jh-1)*2+1) * ssvmask(ji,jj) 306 out_v(ji,jj,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssvmask(ji,jj) 404 # if defined key_3Ddiaharm 405 out_v(ji,jj,jk, jh) = ztmp7((jh-1)*2+1) * ssvmask(ji,jj) 406 out_v(ji,jj,jk,nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssvmask(ji,jj) 407 # else 408 out_v(ji,jj, jh) = ztmp7((jh-1)*2+1) * ssvmask(ji,jj) 409 out_v(ji,jj, nb_ana+jh) = -ztmp7((jh-1)*2+2) * ssvmask(ji,jj) 410 # endif 307 411 END DO 308 412 309 413 END DO 310 414 END DO 415 416 # if defined key_3Ddiaharm 417 ! w- velocity 418 DO jj = 1, jpj 419 DO ji = 1, jpi 420 ! Fill input array 421 kun=0 422 DO jh = 1,nb_ana 423 DO jc = 1,2 424 kun = kun + 1 425 ztmp4(kun)=ana_temp(ji,jj,kun,4,jk) 426 END DO 427 END DO 428 429 CALL SUR_DETERMINE(jj+1) 430 431 ! Fill output array 432 DO jh = 1, nb_ana 433 ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1) 434 ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2) 435 END DO 436 437 END DO 438 END DO 439 440 DO jj = 1, jpj 441 DO ji = 1, jpi 442 DO jh = 1, nb_ana 443 X1=ana_amp(ji,jj,jh,1) 444 X2=-ana_amp(ji,jj,jh,2) 445 out_w(ji,jj,jk, jh)=X1 * tmask_i(ji,jj) 446 out_w(ji,jj,jk,nb_ana+jh)=X2 * tmask_i(ji,jj) 447 END DO 448 END DO 449 END DO 450 451 ! dzi- isopycnal displacements 452 DO jj = 1, jpj 453 DO ji = 1, jpi 454 ! Fill input array 455 kun=0 456 DO jh = 1,nb_ana 457 DO jc = 1,2 458 kun = kun + 1 459 ztmp4(kun)=ana_temp(ji,jj,kun,5,jk) 460 END DO 461 END DO 462 463 CALL SUR_DETERMINE(jj+1) 464 465 ! Fill output array 466 DO jh = 1, nb_ana 467 ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1) 468 ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2) 469 END DO 470 471 END DO 472 END DO 473 474 DO jj = 1, jpj 475 DO ji = 1, jpi 476 DO jh = 1, nb_ana 477 X1=ana_amp(ji,jj,jh,1) 478 X2=-ana_amp(ji,jj,jh,2) 479 out_dzi(ji,jj,jk, jh)=X1 * tmask_i(ji,jj) 480 out_dzi(ji,jj,jk,nb_ana+jh)=X2 * tmask_i(ji,jj) 481 END DO 482 END DO 483 END DO 484 485 ENDDO ! jk 486 # endif 311 487 ! 312 488 ! clem: we could avoid this call if all the loops were from 1:jpi and 1:jpj … … 328 504 !!-------------------------------------------------------------------- 329 505 INTEGER :: jh 330 !!---------------------------------------------------------------------- 506 507 # if defined key_3Ddiaharm 508 CHARACTER(LEN=lc) :: cdfile_name_W ! name of the file created (W-points) 509 INTEGER :: jk 510 REAL(WP), ALLOCATABLE, DIMENSION (:,:,:) :: z3real, z3im 511 REAL(WP), ALLOCATABLE, DIMENSION (:,:) :: z2real, z2im 512 # endif 513 !!---------------------------------------------------------------------- 514 515 #if defined key_dimgout 516 cdfile_name_T = TRIM(cexper)//'_Tidal_harmonics_gridT.dimgproc' 517 cdfile_name_U = TRIM(cexper)//'_Tidal_harmonics_gridU.dimgproc' 518 cdfile_name_V = TRIM(cexper)//'_Tidal_harmonics_gridV.dimgproc' 519 # if defined key_3Ddiaharm 520 cdfile_name_W = TRIM(cexper)//'_Tidal_harmonics_gridW.dimgproc' 521 # endif 522 #endif 331 523 332 524 IF(lwp) WRITE(numout,*) ' ' 333 525 IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results' 526 #if defined key_dimgout 527 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ Output files: ', TRIM(cdfile_name_T) 528 IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_U) 529 IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_V) 530 # if defined key_3Ddiaharm 531 IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_W) 532 # endif 533 #endif 334 534 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 335 535 336 ! A) Elevation 536 # if defined key_3Ddiaharm 537 ALLOCATE(z3real(jpi,jpj,jpk),z3im(jpi,jpj,jpk),z2real(jpi,jpj),z2im(jpi,jpj)) 538 # endif 539 540 ! A) density and elevation 337 541 !///////////// 542 #if defined key_dimgout 543 cltext='density amplitude and phase; elevation is level=jpk ' 544 CALL dia_wri_dimg(TRIM(cdfile_name_T), TRIM(cltext), out_eta, 2*nb_ana, '2') 545 #else 546 # if defined key_3Ddiaharm 547 z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp 548 # endif 338 549 DO jh = 1, nb_ana 550 # if defined key_3Ddiaharm 551 DO jk=1,jpkm1 552 z3real(:,:,jk)=out_eta(:,:,jk,jh) 553 z3im (:,:,jk)=out_eta(:,:,jk,jh+nb_ana) 554 ENDDO 555 z2real(:,:)=out_eta(:,:,jpk,jh); z2im(:,:)=out_eta(:,:,jpk,jh+nb_ana) 556 CALL iom_put( TRIM(tname(jh))//'x_ro', z3real(:,:,:) ) 557 CALL iom_put( TRIM(tname(jh))//'y_ro', z3im (:,:,:) ) 558 CALL iom_put( TRIM(tname(jh))//'x' , z2real(:,: ) ) 559 CALL iom_put( TRIM(tname(jh))//'y' , z2im (:,: ) ) 560 # else 561 WRITE(numout,*) "OUTPUT ORI: ", TRIM(tname(jh))//'x', ' & ', TRIM(tname(jh))//'y', MAXVAL(out_eta(:,:,jh)) 339 562 CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) ) 340 563 CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,jh+nb_ana) ) 341 END DO 342 343 ! B) ubar 564 # endif 565 END DO 566 #endif 567 568 ! B) u 344 569 !///////// 570 #if defined key_dimgout 571 cltext='3d u amplitude and phase; ubar is the last level' 572 CALL dia_wri_dimg(TRIM(cdfile_name_U), TRIM(cltext), out_u, 2*nb_ana, '2') 573 #else 574 # if defined key_3Ddiaharm 575 z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp 576 # endif 345 577 DO jh = 1, nb_ana 346 CALL iom_put( TRIM(tname(jh))//'x_u', out_u(:,:,jh) ) 347 CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:,jh+nb_ana) ) 348 END DO 349 350 ! C) vbar 578 # if defined key_3Ddiaharm 579 DO jk=1,jpkm1 580 z3real(:,:,jk)=out_u(:,:,jk,jh) 581 z3im (:,:,jk)=out_u(:,:,jk,jh+nb_ana) 582 ENDDO 583 z2real(:,:)=out_u(:,:,jpk,jh); z2im(:,:)=out_u(:,:,jpk,jh+nb_ana) 584 CALL iom_put( TRIM(tname(jh))//'x_u3d', z3real(:,:,:) ) 585 CALL iom_put( TRIM(tname(jh))//'y_u3d', z3im (:,:,:) ) 586 CALL iom_put( TRIM(tname(jh))//'x_u2d', z2real(:,:) ) 587 CALL iom_put( TRIM(tname(jh))//'y_u2d', z2im (:,:) ) 588 z2real(:,:)=out_w(:,:,jpk,jh); z2im(:,:)=out_w(:,:,jpk,jh+nb_ana) 589 CALL iom_put( TRIM(tname(jh))//'x_tabx', z2real(:,:) ) 590 CALL iom_put( TRIM(tname(jh))//'y_tabx', z2im (:,:) ) 591 # else 592 CALL iom_put( TRIM(tname(jh))//'x_u2d', out_u(:,:,jh) ) 593 CALL iom_put( TRIM(tname(jh))//'y_u2d', out_u(:,:,nb_ana+jh) ) 594 # endif 595 END DO 596 #endif 597 598 ! C) v 351 599 !///////// 600 #if defined key_dimgout 601 cltext='3d v amplitude and phase; vbar is the last level' 602 CALL dia_wri_dimg(TRIM(cdfile_name_V), TRIM(cltext), out_v, 2*nb_ana, '2') 603 #else 604 # if defined key_3Ddiaharm 605 z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp 606 # endif 352 607 DO jh = 1, nb_ana 353 CALL iom_put( TRIM(tname(jh))//'x_v', out_v(:,:,jh ) ) 354 CALL iom_put( TRIM(tname(jh))//'y_v', out_v(:,:,jh+nb_ana) ) 355 END DO 356 ! 608 # if defined key_3Ddiaharm 609 DO jk=1,jpkm1 610 z3real(:,:,jk)=out_v(:,:,jk,jh) 611 z3im (:,:,jk)=out_v(:,:,jk,jh+nb_ana) 612 ENDDO 613 z2real(:,:)=out_v(:,:,jpk,jh); z2im(:,:)=out_v(:,:,jpk,jh+nb_ana) 614 CALL iom_put( TRIM(tname(jh))//'x_v3d', z3real(:,:,:) ) 615 CALL iom_put( TRIM(tname(jh))//'y_v3d', z3im (:,:,:) ) 616 CALL iom_put( TRIM(tname(jh))//'x_v2d' , z2real(:,:) ) 617 CALL iom_put( TRIM(tname(jh))//'y_v2d' , z2im (:,:) ) 618 z2real(:,:)=out_dzi(:,:,jpk,jh); z2im(:,:)=out_dzi(:,:,jpk,jh+nb_ana) 619 CALL iom_put( TRIM(tname(jh))//'x_taby', z2real(:,:) ) 620 CALL iom_put( TRIM(tname(jh))//'y_taby', z2im (:,:) ) 621 # else 622 CALL iom_put( TRIM(tname(jh))//'x_v2d', out_v(:,:,jh ) ) 623 CALL iom_put( TRIM(tname(jh))//'y_v2d', out_v(:,:,jh+nb_ana) ) 624 # endif 625 END DO 626 627 #endif 628 ! D) w 629 # if defined key_3Ddiaharm 630 # if defined key_dimgout 631 cltext='3d w amplitude and phase; vort_baro is the last level' 632 CALL dia_wri_dimg(TRIM(cdfile_name_W), TRIM(cltext), out_w, 2*nb_ana, '2') 633 # else 634 DO jh = 1, nb_ana 635 DO jk=1,jpkm1 636 z3real(:,:,jk)=out_w(:,:,jk,jh) 637 z3im(:,:,jk)=out_w(:,:,jk,jh+nb_ana) 638 ENDDO 639 CALL iom_put( TRIM(tname(jh))//'x_w3d', z3real(:,:,:) ) 640 CALL iom_put( TRIM(tname(jh))//'y_w3d', z3im(:,:,:) ) 641 END DO 642 # endif 643 644 ! E) dzi + tau_bot 645 # if defined key_dimgout 646 cltext='dzi=g*ro/N2 amplitude and phase' 647 CALL dia_wri_dimg(TRIM(cdfile_name_W), TRIM(cltext), out_w, 2*nb_ana, '2') 648 # else 649 DO jh = 1, nb_ana 650 DO jk=1,jpkm1 651 z3real(:,:,jk)=out_dzi(:,:,jk,jh) 652 z3im(:,:,jk)=out_dzi(:,:,jk,jh+nb_ana) 653 ENDDO 654 CALL iom_put( TRIM(tname(jh))//'x_dzi', z3real(:,:,:) ) 655 CALL iom_put( TRIM(tname(jh))//'y_dzi', z3im(:,:,:) ) 656 END DO 657 # endif 658 # endif 659 660 ! 661 # if defined key_3Ddiaharm 662 DEALLOCATE(z3real, z3im, z2real,z2im) 663 # endif 664 357 665 END SUBROUTINE dia_wri_harm 358 666 -
NEMO/branches/UKMO/NEMO_4.0.4_CO9_wadcpl/src/OCE/DOM/dtatsd.F90
r15784 r15788 21 21 ! 22 22 USE in_out_manager ! I/O manager 23 USE iom ! IOM library 23 24 USE lib_mpp ! MPP library 24 25 … … 31 32 ! !!* namtsd namelist : Temperature & Salinity Data * 32 33 LOGICAL , PUBLIC :: ln_tsd_init !: T & S data flag 34 LOGICAL , PUBLIC :: ln_tsd_interp !: vertical interpolation flag 33 35 LOGICAL , PUBLIC :: ln_tsd_dmp !: internal damping toward input data flag 34 36 35 37 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_tsd ! structure of input SST (file informations, fields read) 38 INTEGER :: jpk_init, inum_dta 39 INTEGER :: id, linum ! local integers 40 INTEGER :: zdim(4) 36 41 37 42 !!---------------------------------------------------------------------- … … 53 58 LOGICAL, INTENT(in), OPTIONAL :: ld_tradmp ! force the initialization when tradp is used 54 59 ! 55 INTEGER :: ios, ierr0, ierr1, ierr2, ierr3 ! local integers60 INTEGER :: ios, ierr0, ierr1, ierr2, ierr3, ierr4, ierr5 ! local integers 56 61 !! 57 62 CHARACTER(len=100) :: cn_dir ! Root directory for location of ssr files 58 TYPE(FLD_N), DIMENSION( jpts) :: slf_i ! array of namelist informations on the fields to read59 TYPE(FLD_N) :: sn_tem, sn_sal 63 TYPE(FLD_N), DIMENSION(jpts+2):: slf_i ! array of namelist informations on the fields to read 64 TYPE(FLD_N) :: sn_tem, sn_sal, sn_dep, sn_msk 60 65 !! 61 NAMELIST/namtsd/ ln_tsd_init, ln_tsd_ dmp, cn_dir, sn_tem, sn_sal66 NAMELIST/namtsd/ ln_tsd_init, ln_tsd_interp, ln_tsd_dmp, cn_dir, sn_tem, sn_sal, sn_dep, sn_msk 62 67 !!---------------------------------------------------------------------- 63 68 ! 64 69 ! Initialisation 65 ierr0 = 0 ; ierr1 = 0 ; ierr2 = 0 ; ierr3 = 0 70 ierr0 = 0 ; ierr1 = 0 ; ierr2 = 0 ; ierr3 = 0 ; ierr4 = 0 ; ierr5 = 0 66 71 ! 67 72 REWIND( numnam_ref ) ! Namelist namtsd in reference namelist : … … 80 85 WRITE(numout,*) '~~~~~~~~~~~~ ' 81 86 WRITE(numout,*) ' Namelist namtsd' 82 WRITE(numout,*) ' Initialisation of ocean T & S with T &S input data ln_tsd_init = ', ln_tsd_init 83 WRITE(numout,*) ' damping of ocean T & S toward T &S input data ln_tsd_dmp = ', ln_tsd_dmp 87 WRITE(numout,*) ' Initialisation of ocean T & S with T & S input data ln_tsd_init = ', ln_tsd_init 88 WRITE(numout,*) ' Interpolation of initial conditions in the vertical ln_tsd_interp = ', ln_tsd_interp 89 WRITE(numout,*) ' damping of ocean T & S toward T & S input data ln_tsd_dmp = ', ln_tsd_dmp 84 90 WRITE(numout,*) 85 91 IF( .NOT.ln_tsd_init .AND. .NOT.ln_tsd_dmp ) THEN … … 94 100 ln_tsd_init = .FALSE. 95 101 ENDIF 102 IF( ln_tsd_interp .AND. ln_tsd_dmp ) THEN 103 CALL ctl_stop( 'dta_tsd_init: Tracer damping and vertical interpolation not yet configured' ) ; RETURN 104 ENDIF 105 IF( ln_tsd_interp .AND. LEN(TRIM(sn_msk%wname)) > 0 ) THEN 106 CALL ctl_stop( 'dta_tsd_init: Using vertical interpolation and weights files not recommended' ) ; RETURN 107 ENDIF 96 108 ! 97 109 ! ! allocate the arrays (if necessary) 98 110 IF( ln_tsd_init .OR. ln_tsd_dmp ) THEN 99 111 ! 100 ALLOCATE( sf_tsd(jpts), STAT=ierr0 ) 112 IF( ln_tsd_interp ) THEN 113 ALLOCATE( sf_tsd(jpts+2), STAT=ierr0 ) ! to carry the addtional depth information 114 ELSE 115 ALLOCATE( sf_tsd(jpts ), STAT=ierr0 ) 116 ENDIF 101 117 IF( ierr0 > 0 ) THEN 102 118 CALL ctl_stop( 'dta_tsd_init: unable to allocate sf_tsd structure' ) ; RETURN 103 119 ENDIF 104 120 ! 105 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk) , STAT=ierr0 ) 106 IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 ) 107 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk) , STAT=ierr2 ) 108 IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 ) 109 ! 110 IF( ierr0 + ierr1 + ierr2 + ierr3 > 0 ) THEN 121 IF( ln_tsd_interp ) THEN 122 CALL iom_open ( trim(cn_dir) // trim(sn_dep%clname), inum_dta ) 123 id = iom_varid( inum_dta, sn_dep%clvar, zdim ) 124 jpk_init = zdim(3) 125 IF(lwp) WRITE(numout,*) 'Dimension of vertical coordinate in ICs: ', jpk_init 126 CALL iom_close( inum_dta ) ! Close the input file 127 ! 128 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk_init ) , STAT=ierr0 ) 129 IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk_init,2) , STAT=ierr1 ) 130 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk_init ) , STAT=ierr2 ) 131 IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk_init,2) , STAT=ierr3 ) 132 ALLOCATE( sf_tsd(jp_dep)%fnow(jpi,jpj,jpk_init ) , STAT=ierr4 ) 133 ALLOCATE( sf_tsd(jp_msk)%fnow(jpi,jpj,jpk_init ) , STAT=ierr5 ) 134 ELSE 135 ALLOCATE( sf_tsd(jp_tem)%fnow(jpi,jpj,jpk) , STAT=ierr0 ) 136 IF( sn_tem%ln_tint ) ALLOCATE( sf_tsd(jp_tem)%fdta(jpi,jpj,jpk,2) , STAT=ierr1 ) 137 ALLOCATE( sf_tsd(jp_sal)%fnow(jpi,jpj,jpk) , STAT=ierr2 ) 138 IF( sn_sal%ln_tint ) ALLOCATE( sf_tsd(jp_sal)%fdta(jpi,jpj,jpk,2) , STAT=ierr3 ) 139 ENDIF ! ln_tsd_interp 140 ! 141 IF( ierr0 + ierr1 + ierr2 + ierr3 + ierr4 + ierr5 > 0 ) THEN 111 142 CALL ctl_stop( 'dta_tsd : unable to allocate T & S data arrays' ) ; RETURN 112 143 ENDIF 113 144 ! ! fill sf_tsd with sn_tem & sn_sal and control print 114 145 slf_i(jp_tem) = sn_tem ; slf_i(jp_sal) = sn_sal 146 IF( ln_tsd_interp ) slf_i(jp_dep) = sn_dep ; slf_i(jp_msk) = sn_msk 115 147 CALL fld_fill( sf_tsd, slf_i, cn_dir, 'dta_tsd', 'Temperature & Salinity data', 'namtsd', no_print ) 116 148 ! … … 138 170 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT( out) :: ptsd ! T & S data 139 171 ! 140 INTEGER :: ji, jj, jk, jl, jk k ! dummy loop indicies141 INTEGER :: ik, il0, il1, ii0, ii1, ij0, ij1 ! local integers172 INTEGER :: ji, jj, jk, jl, jk_init ! dummy loop indices 173 INTEGER :: ik, il0, il1, ii0, ii1, ij0, ij1 ! local integers 142 174 REAL(wp):: zl, zi ! local scalars 143 REAL(wp), DIMENSION(jpk) :: ztp, zsp ! 1D workspace144 175 !!---------------------------------------------------------------------- 145 176 ! … … 176 207 !!gm end 177 208 ! 178 ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:) ! NO mask 179 ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:) 180 ! 181 !CEOD, We think this is incorrect we dont want to do this. 209 IF( kt == nit000 .AND. lwp )THEN 210 WRITE(numout,*) 211 WRITE(numout,*) 'dta_tsd: interpolates T & S data onto current mesh' 212 ENDIF 213 ! 214 IF( ln_tsd_interp ) THEN ! probably should use pointers in the following to make more readable 215 ! 216 DO jk = 1, jpk ! determines the intepolated T-S profiles at each (i,j) points 217 DO jj= 1, jpj 218 DO ji= 1, jpi 219 zl = gdept_0(ji,jj,jk) 220 IF( zl < sf_tsd(jp_dep)%fnow(ji,jj,1) ) THEN ! above the first level of data 221 ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,1) 222 ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,1) 223 ELSEIF( zl > sf_tsd(jp_dep)%fnow(ji,jj,jpk_init) ) THEN ! below the last level of data 224 ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jpk_init) 225 ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jpk_init) 226 ELSE ! inbetween : vertical interpolation between jk_init & jk_init+1 227 DO jk_init = 1, jpk_init-1 ! when gdept(jk_init) < zl < gdept(jk_init+1) 228 IF( sf_tsd(jp_msk)%fnow(ji,jj,jk_init+1) == 0 ) THEN ! if there is no data fill down 229 sf_tsd(jp_tem)%fnow(ji,jj,jk_init+1) = sf_tsd(jp_tem)%fnow(ji,jj,jk_init) 230 sf_tsd(jp_sal)%fnow(ji,jj,jk_init+1) = sf_tsd(jp_sal)%fnow(ji,jj,jk_init) 231 ENDIF 232 IF((zl-sf_tsd(jp_dep)%fnow(ji,jj,jk_init)) * (zl-sf_tsd(jp_dep)%fnow(ji,jj,jk_init+1)) <= 0._wp ) THEN 233 zi = ( zl - sf_tsd(jp_dep)%fnow(ji,jj,jk_init) ) / & 234 & (sf_tsd(jp_dep)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_dep)%fnow(ji,jj,jk_init)) 235 ptsd(ji,jj,jk,jp_tem) = sf_tsd(jp_tem)%fnow(ji,jj,jk_init) + & 236 & (sf_tsd(jp_tem)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_tem)%fnow(ji,jj,jk_init)) * zi 237 ptsd(ji,jj,jk,jp_sal) = sf_tsd(jp_sal)%fnow(ji,jj,jk_init) + & 238 & (sf_tsd(jp_sal)%fnow(ji,jj,jk_init+1)-sf_tsd(jp_sal)%fnow(ji,jj,jk_init)) * zi 239 ENDIF 240 END DO 241 ENDIF 242 END DO 243 END DO 244 END DO 245 ! 246 ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) *tmask(:,:,:) 247 ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) *tmask(:,:,:) 248 ELSE !== z- or zps- coordinate ==! 182 249 ! 183 ptsd(:,:,:,jp_tem) = ptsd(:,:,:,jp_tem) * tmask(:,:,:) ! Mask184 ptsd(:,:,:,jp_sal) = ptsd(:,:,:,jp_sal) * tmask(:,:,:)250 ptsd(:,:,:,jp_tem) = sf_tsd(jp_tem)%fnow(:,:,:) * tmask(:,:,:) ! Mask 251 ptsd(:,:,:,jp_sal) = sf_tsd(jp_sal)%fnow(:,:,:) * tmask(:,:,:) 185 252 ! 186 253 IF( ln_zps ) THEN ! zps-coordinate (partial steps) interpolation at the last ocean level … … 212 279 DEALLOCATE( sf_tsd(jp_sal)%fnow ) ! S arrays in the structure 213 280 IF( sf_tsd(jp_sal)%ln_tint ) DEALLOCATE( sf_tsd(jp_sal)%fdta ) 281 IF( ln_tsd_interp ) DEALLOCATE( sf_tsd(jp_dep)%fnow ) ! T arrays in the structure 282 IF( ln_tsd_interp ) DEALLOCATE( sf_tsd(jp_msk)%fnow ) ! T arrays in the structure 214 283 DEALLOCATE( sf_tsd ) ! the structure itself 215 284 ENDIF -
NEMO/branches/UKMO/NEMO_4.0.4_CO9_wadcpl/src/OCE/DYN/dynnxt.F90
r14075 r15788 33 33 USE dynadv ! dynamics: vector invariant versus flux form 34 34 USE dynspg_ts ! surface pressure gradient: split-explicit scheme 35 USE dynspg 35 36 USE domvvl ! variable volume 36 37 USE bdy_oce , ONLY : ln_bdy … … 58 59 PUBLIC dyn_nxt ! routine called by step.F90 59 60 61 !! Substitution 62 # include "vectopt_loop_substitute.h90" 60 63 !!---------------------------------------------------------------------- 61 64 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 182 185 vn(:,:,jk) = va(:,:,jk) 183 186 END DO 187 ! limit velocities 188 IF (ln_ulimit) THEN 189 call dyn_limit_velocity (kt) 190 ENDIF 184 191 IF( .NOT.ln_linssh ) THEN ! e3._b <-- e3._n 185 192 !!gm BUG ???? I don't understand why it is not : e3._n <-- e3._a … … 214 221 END DO 215 222 END DO 223 ! limit velocities 224 IF (ln_ulimit) THEN 225 call dyn_limit_velocity (kt) 226 ENDIF 216 227 ! ! ================! 217 228 ELSE ! Variable volume ! … … 263 274 END DO 264 275 END DO 276 ! limit velocities 277 IF (ln_ulimit) THEN 278 call dyn_limit_velocity (kt) 279 ENDIF 265 280 ! 266 281 ELSE ! Asselin filter applied on thickness weighted velocity … … 290 305 END DO 291 306 END DO 307 ! limit velocities 308 IF (ln_ulimit) THEN 309 call dyn_limit_velocity (kt) 310 ENDIF 292 311 e3u_b(:,:,1:jpkm1) = ze3u_f(:,:,1:jpkm1) ! e3u_b <-- filtered scale factor 293 312 e3v_b(:,:,1:jpkm1) = ze3v_f(:,:,1:jpkm1) … … 401 420 END SUBROUTINE dyn_nxt 402 421 422 SUBROUTINE dyn_limit_velocity (kt) 423 ! limits maximum values of un and vn by volume courant number 424 INTEGER, INTENT( in ) :: kt ! ocean time-step index 425 ! 426 INTEGER :: ji, jj, jk ! dummy loop indices 427 REAL(wp) :: zzu,zplim,zmlim,isp,ism,zcn,ze3e1,zzcn,zcnn,idivp,idivm 428 429 ! limit fluxes 430 zcn =cn_ulimit !0.9 ! maximum velocity inverse courant number 431 zcnn = cnn_ulimit !0.54 ! how much to reduce cn by in divergen flow 432 433 DO jk = 1, jpkm1 434 DO jj = 2, jpjm1 435 DO ji = fs_2, fs_jpim1 ! vect. opt. 436 ! U direction 437 zzu = un(ji,jj,jk) 438 ze3e1 = e3u_n(ji ,jj,jk) * e2u(ji ,jj) 439 ! ips is 1 if flow is positive othersize zero 440 isp = 0.5 * (sign(1.0_wp,zzu) + 1.0_wp ) 441 ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp ) 442 ! idev is 1 if divergent flow otherwise zero 443 idivp = -isp * 0.5 * (sign(1.0_wp, un(ji-1,jj,jk)) - 1.0_wp ) 444 idivm = ism * 0.5 * (sign(1.0_wp, un(ji+1,jj,jk)) + 1.0_wp ) 445 zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp 446 zzcn = zzcn * zcn 447 zplim = zzcn * (e3t_n(ji ,jj,jk) * e1t(ji ,jj) * e2t(ji ,jj)) / & 448 (2.0*rdt * ze3e1)*umask(ji,jj,jk) 449 zmlim = -zzcn * (e3t_n(ji+1,jj,jk) * e1t(ji+1,jj) * e2t(ji+1,jj)) / & 450 (2.0*rdt * ze3e1)*umask(ji,jj,jk) 451 ! limit currents 452 un(ji,jj,jk) = min ( zzu,zplim) * isp + max(zzu,zmlim) *ism 453 454 ! V direction 455 zzu = vn(ji,jj,jk) 456 ze3e1 = e3v_n(ji ,jj,jk) * e1v(ji ,jj) 457 isp = 0.5 * (sign(1.0_wp,zzu) + 1.0_wp ) 458 ism = -0.5 * (sign(1.0_wp,zzu) - 1.0_wp ) 459 ! idev is 1 if divergent flow otherwise zero 460 idivp = -isp * 0.5 * (sign(1.0_wp, vn(ji,jj-1,jk)) - 1.0_wp ) 461 idivm = ism * 0.5 * (sign(1.0_wp, vn(ji,jj+1,jk)) + 1.0_wp ) 462 zzcn = (idivp+idivm)*(zcnn-1.0_wp)+1.0_wp 463 zzcn = zzcn * zcn 464 zplim = zzcn * (e3t_n(ji,jj ,jk) * e1t(ji,jj ) * e2t(ji,jj )) / & 465 (2.0*rdt * ze3e1)*vmask(ji,jj,jk) 466 zmlim = -zzcn * (e3t_n(ji,jj+1,jk) * e1t(ji,jj+1) * e2t(ji,jj+1)) / & 467 (2.0*rdt * ze3e1)*vmask(ji,jj,jk) 468 ! limit currents 469 vn(ji,jj,jk) = min ( zzu,zplim) * isp + max(zzu,zmlim) *ism 470 ENDDO 471 ENDDO 472 ENDDO 473 CALL lbc_lnk_multi( 'dynnxt', un(:,:,:), 'U', -1., vn(:,:,:), 'V', -1. ) 474 475 END SUBROUTINE dyn_limit_velocity 476 403 477 !!========================================================================= 404 478 END MODULE dynnxt -
NEMO/branches/UKMO/NEMO_4.0.4_CO9_wadcpl/src/OCE/DYN/dynspg.F90
r14075 r15788 39 39 INTEGER :: nspg = 0 ! type of surface pressure gradient scheme defined from lk_dynspg_... 40 40 41 LOGICAL, PUBLIC :: ln_ulimit 42 REAL(wp), PUBLIC :: cn_ulimit,cnn_ulimit 43 41 44 ! ! Parameter to control the surface pressure gradient scheme 42 45 INTEGER, PARAMETER :: np_TS = 1 ! split-explicit time stepping (Time-Splitting) … … 191 194 NAMELIST/namdyn_spg/ ln_dynspg_exp , ln_dynspg_ts, & 192 195 & ln_bt_fw, ln_bt_av , ln_bt_auto , & 193 & nn_baro , rn_bt_cmax, nn_bt_flt, rn_bt_alpha 196 & nn_baro , rn_bt_cmax, nn_bt_flt, rn_bt_alpha, & 197 & ln_ulimit, cn_ulimit, cnn_ulimit 194 198 !!---------------------------------------------------------------------- 195 199 ! … … 213 217 WRITE(numout,*) ' Explicit free surface ln_dynspg_exp = ', ln_dynspg_exp 214 218 WRITE(numout,*) ' Free surface with time splitting ln_dynspg_ts = ', ln_dynspg_ts 219 WRITE(numout,*) ' Limit velocities ln_ulimit = ', ln_ulimit 220 WRITE(numout,*) ' Limit velocities max inverse Courant number = ', cn_ulimit 221 WRITE(numout,*) ' Limit velocities multiplier for divergant flow = ', cnn_ulimit 215 222 ENDIF 216 223 ! ! Control of surface pressure gradient scheme options -
NEMO/branches/UKMO/NEMO_4.0.4_CO9_wadcpl/src/OCE/SBC/fldread.F90
r15786 r15788 909 909 IF( zdepth(jk) <= pdta_read_z(jb,1,1) ) THEN ! above the first level of external data 910 910 pdta(jb,1,jk) = pdta_read(jb,1,1) 911 ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkmax) ) THEN ! below the last level of external data /= FillValue912 pdta(jb,1,jk) = pdta_read(jb,1,ipkmax)911 ELSEIF( zdepth(jk) > pdta_read_z(jb,1,ipkmax) ) THEN ! below the last level of external data 912 pdta(jb,1,jk) = pdta_read(jb,1,ipkmax) 913 913 ELSE ! inbetween: vertical interpolation between jkb & jkb+1 914 914 DO jkb = 1, ipkmax-1 -
NEMO/branches/UKMO/NEMO_4.0.4_CO9_wadcpl/src/OCE/SBC/sbctide.F90
r15784 r15788 16 16 USE ioipsl ! NetCDF IPSL library 17 17 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 18 USE bdytides 18 19 19 20 USE bdytides ! davbyr - Access to love number … … 110 111 111 112 DO jk = 1, nb_harmo 112 ! davbyr - Insert variable Love number where once was 0.7113 ! love number now provides in tide namelist 113 114 zcons = dn_love_number * Wave(ntide(jk))%equitide * ftide(jk) 114 ! END davbyr115 115 DO ji = 1, jpi 116 116 DO jj = 1, jpj … … 123 123 IF ( Wave(ntide(jk))%nutide == 1 ) THEN ; zcs = zcons * SIN( 2._wp*zlat ) 124 124 ELSEIF( Wave(ntide(jk))%nutide == 2 ) THEN ; zcs = zcons * COS( zlat )**2 125 ! davbyr - Include long period tidal forcing125 ! Add tide potential for long period tides 126 126 ELSEIF( Wave(ntide(jk))%nutide == 0 ) THEN ; zcs = zcons * (0.5_wp-1.5_wp*SIN(zlat)**2._wp) 127 ! END - davbyr128 127 ELSE ; zcs = 0._wp 129 128 ENDIF -
NEMO/branches/UKMO/NEMO_4.0.4_CO9_wadcpl/src/OCE/SBC/tide_mod.F90
r15785 r15788 16 16 PUBLIC tide_init_Wave ! called by tideini and diaharm modules 17 17 18 ! davbyr: increase maximum number of harmonics from 19 to 3419 ! RCS: merge in extra tides.20 18 INTEGER, PUBLIC, PARAMETER :: jpmax_harmo = 43 !: maximum number of harmonic 21 19 … … 43 41 44 42 SUBROUTINE tide_init_Wave 43 # if defined key_FES14_tides 44 # include "tide_FES14.h90" 45 # else 45 46 # include "tide.h90" 47 # endif 46 48 END SUBROUTINE tide_init_Wave 47 49 … … 333 335 zf = zf * zf1 * zf1 334 336 ! 335 336 !--- davbyr 11/2017 337 CASE( 20 ) !== formule 20, compound waves ( 78 x 78 x 78 x 78 ) 338 zf1 = nodal_factort(78) 339 zf = zf1 * zf1 * zf1 * zf1 340 !--- END davbyr 337 CASE( 20 ) !== formule 20, compound waves ( 78 x 78 x 78 x 78 ) 338 zf1 = nodal_factort(78) 339 zf = zf1 * zf1 * zf1 * zf1 340 ! 341 341 CASE( 73 ) !== formule 73 342 342 zs = sin(sh_I) -
NEMO/branches/UKMO/NEMO_4.0.4_CO9_wadcpl/src/OCE/SBC/tideini.F90
r15784 r15788 33 33 INTEGER , PUBLIC :: kt_tide !: 34 34 REAL(wp), PUBLIC :: rdttideramp !: 35 REAL(wp), PUBLIC :: dn_love_number !: 35 36 REAL(wp), PUBLIC :: rn_scal_load !: 36 37 ! davbyr - read love number from namelist … … 81 82 WRITE(numout,*) ' Read load potential from file ln_read_load = ', ln_read_load 82 83 WRITE(numout,*) ' Apply ramp on tides at startup ln_tide_ramp = ', ln_tide_ramp 84 WRITE(numout,*) ' dn_love_number = ', dn_love_number 83 85 WRITE(numout,*) ' Fraction of SSH used in scal. approx. rn_scal_load = ', rn_scal_load 84 86 WRITE(numout,*) ' Duration (days) of ramp rdttideramp = ', rdttideramp … … 104 106 END DO 105 107 ! 108 IF (ln_tide .and.lwp) WRITE(numout,*) 'nb_harmo = ', nb_harmo 109 106 110 ! Ensure that tidal components have been set in namelist_cfg 107 111 IF( nb_harmo == 0 ) CALL ctl_stop( 'tide_init : No tidal components set in nam_tide' ) -
NEMO/branches/UKMO/NEMO_4.0.4_CO9_wadcpl/src/OCE/USR/usrdef_istate.F90
r14075 r15788 62 62 DO jj = 1, jpj 63 63 DO ji = 1, jpi 64 pts(ji,jj,jk,jp_tem) = ( ( 16. - 12. * TANH( (pdept(ji,jj,jk) - 400) / 700 ) ) & 65 & * (-TANH( (500. - pdept(ji,jj,jk)) / 150. ) + 1.) / 2. & 66 & + ( 15. * ( 1. - TANH( (pdept(ji,jj,jk)-50.) / 1500.) ) & 67 & - 1.4 * TANH((pdept(ji,jj,jk)-100.) / 100.) & 68 & + 7. * (1500. - pdept(ji,jj,jk) ) / 1500.) & 69 & * (-TANH( (pdept(ji,jj,jk) - 500.) / 150.) + 1.) / 2. ) * ptmask(ji,jj,jk) 70 71 pts(ji,jj,jk,jp_sal) = ( ( 36.25 - 1.13 * TANH( (pdept(ji,jj,jk) - 305) / 460 ) ) & 72 & * (-TANH((500. - pdept(ji,jj,jk)) / 150.) + 1.) / 2 & 73 & + ( 35.55 + 1.25 * (5000. - pdept(ji,jj,jk)) / 5000. & 74 & - 1.62 * TANH( (pdept(ji,jj,jk) - 60. ) / 650. ) & 75 & + 0.2 * TANH( (pdept(ji,jj,jk) - 35. ) / 100. ) & 76 & + 0.2 * TANH( (pdept(ji,jj,jk) - 1000.) / 5000.) ) & 77 & * (-TANH( (pdept(ji,jj,jk) - 500.) / 150.) + 1.) / 2 ) * ptmask(ji,jj,jk) 64 pts(ji,jj,jk,jp_tem) = 20._wp * ptmask(ji,jj,jk) 65 pts(ji,jj,jk,jp_sal) = 36.25_wp * ptmask(ji,jj,jk) 78 66 END DO 79 67 END DO -
NEMO/branches/UKMO/NEMO_4.0.4_CO9_wadcpl/src/OCE/par_oce.F90
r14075 r15788 64 64 INTEGER, PUBLIC, PARAMETER :: jp_tem = 1 !: indice for temperature 65 65 INTEGER, PUBLIC, PARAMETER :: jp_sal = 2 !: indice for salinity 66 INTEGER, PUBLIC, PARAMETER :: jp_dep = 3 !: indice for depth 67 INTEGER, PUBLIC, PARAMETER :: jp_msk = 4 !: indice for mask 66 68 67 69 !!---------------------------------------------------------------------- -
NEMO/branches/UKMO/NEMO_4.0.4_CO9_wadcpl/src/OCE/step.F90
r15783 r15788 206 206 CALL dia_ar5 ( kstp ) ! ar5 diag 207 207 IF( ln_diaptr ) CALL dia_ptr ! Poleward adv/ldf TRansports diagnostics 208 IF( lk_diaharm_fast ) & 209 & CALL dia_harm_fast( kstp ) ! Tidal harmonic analysis - restart and faster version 208 210 IF( ln_diaharm ) CALL dia_harm( kstp ) ! Tidal harmonic analysis 209 211 CALL dia_prod( kstp ) ! ocean model: product diagnostics -
NEMO/branches/UKMO/NEMO_4.0.4_CO9_wadcpl/src/OCE/step_oce.F90
r15783 r15788 79 79 USE diahth ! thermocline depth (dia_hth routine) 80 80 USE diahsb ! heat, salt and volume budgets (dia_hsb routine) 81 USE diaharm_fast ! harmonic analysis of tides (harm_ana routine) 81 82 USE diaharm 82 83 USE diaprod -
NEMO/branches/UKMO/NEMO_4.0.4_CO9_wadcpl/src/OCE/stpctl.F90
r14075 r15788 126 126 zmax(1) = MAXVAL( ABS( sshn(:,:) ) ) ! ssh max 127 127 ENDIF 128 zmax(2) = MAXVAL( ABS( un(:,:,:) ) ) ! velocity max (zonal only)128 zmax(2) = MAXVAL( un(:,:,:)*un(:,:,:) + vn(:,:,:)*vn(:,:,:) ) ! velocity max 129 129 llmsk(:,:,:) = tmask(:,:,:) == 1._wp 130 130 IF( COUNT( llmsk(:,:,:) ) > 0 ) THEN ! avoid huge values sent back for land processors... … … 167 167 ! !== error handling ==! 168 168 IF( zmax(1) > 20._wp .OR. & ! too large sea surface height ( > 20 m ) 169 & zmax(2) > 10._wp .OR. & ! too large velocity ( > 10 m/s)169 & zmax(2) > 100._wp .OR. & ! too large velocity ( > 10 m/s) 170 170 & zmax(3) >= 0._wp .OR. & ! negative or zero sea surface salinity 171 171 & zmax(4) >= 100._wp .OR. & ! too large sea surface salinity ( > 100 ) … … 177 177 IF( lwm .AND. kt /= nitend ) istatus = NF90_CLOSE(idrun) 178 178 CALL mpp_maxloc( 'stpctl', ABS(sshn) , ssmask(:,:) , zzz, ih(1:2) ) ; ih(3) = 0 179 CALL mpp_maxloc( 'stpctl', ABS(un), umask (:,:,:), zzz, iu )179 CALL mpp_maxloc( 'stpctl', un*un + vn*vn , umask (:,:,:), zzz, iu ) 180 180 CALL mpp_minloc( 'stpctl', tsn(:,:,:,jp_sal), tmask (:,:,:), zzz, is1 ) 181 181 CALL mpp_maxloc( 'stpctl', tsn(:,:,:,jp_sal), tmask (:,:,:), zzz, is2 ) … … 196 196 is2(:) = MAXLOC( tsn(:,:,:,jp_sal), mask = tmask(:,:,:) == 1._wp ) + (/ nimpp - 1, njmpp - 1, 0 /) 197 197 iareamin(:) = narea ; iareamax(:) = narea ; iareasum(:) = 1 ! this is local information 198 ! The following line could overwrite one of the lines above but it is here to allow merging another branch 199 iu(:) = MAXLOC( un(:,:,:)*un(:,:,:) + vn(:,:,:)*vn(:,:,:) ) + (/ nimpp - 1, njmpp - 1, 0 /) 198 200 ENDIF 199 201 ! 200 202 WRITE(ctmp1,*) ' stp_ctl: |ssh| > 20 m or |U| > 10 m/s or S <= 0 or S >= 100 or NaN encounter in the tests' 201 203 CALL wrt_line(ctmp2, kt, ' |ssh| max ', zmax(1), ih , iareasum(1), iareamin(1), iareamax(1) ) 202 CALL wrt_line(ctmp3, kt, ' |U|max ', zmax(2), iu , iareasum(2), iareamin(2), iareamax(2) )204 CALL wrt_line(ctmp3, kt, ' Vel2 max ', zmax(2), iu , iareasum(2), iareamin(2), iareamax(2) ) 203 205 CALL wrt_line(ctmp4, kt, ' Sal min ', - zmax(3), is1, iareasum(3), iareamin(3), iareamax(3) ) 204 206 CALL wrt_line(ctmp5, kt, ' Sal max ', zmax(4), is2, iareasum(4), iareamin(4), iareamax(4) ) … … 226 228 ENDIF 227 229 ! 228 9500 FORMAT(' it :', i8, ' |ssh|_max: ', D23.16, ' |U|_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16)230 9500 FORMAT(' it :', i8, ' |ssh|_max: ', D23.16, ' Vel2_max: ', D23.16,' S_min: ', D23.16,' S_max: ', D23.16) 229 231 ! 230 232 END SUBROUTINE stp_ctl
Note: See TracChangeset
for help on using the changeset viewer.
