Changeset 11741 for NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_sponge.F90
- Timestamp:
- 2019-10-21T12:26:39+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_sponge.F90
r11625 r11741 258 258 CALL lbc_lnk( 'agrif_Sponge', fsaht_spu, 'U', 1. ) ! Lateral boundary conditions 259 259 CALL lbc_lnk( 'agrif_Sponge', fsaht_spv, 'V', 1. ) 260 260 261 261 spongedoneT = .TRUE. 262 262 ENDIF … … 293 293 INTEGER :: ji, jj, jk, jn ! dummy loop indices 294 294 INTEGER :: iku, ikv 295 REAL(wp) :: ztsa, zabe1, zabe2, zbtr 295 REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot 296 296 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv 297 297 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff … … 302 302 REAL(wp), DIMENSION(1:jpk) :: h_out 303 303 INTEGER :: N_in, N_out 304 REAL(wp) :: h_diff305 304 !!---------------------------------------------------------------------- 306 305 ! … … 317 316 318 317 # if defined key_vertical 319 DO jk=k1,k2 320 DO jj=j1,j2 321 DO ji=i1,i2 322 tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_b(ji,jj,jk) 323 END DO 324 END DO 325 END DO 318 ! Interpolate thicknesses 319 ! Warning: these are masked, hence extrapolated prior interpolation. 320 DO jk=k1,k2-1 321 DO jj=j1,j2 322 DO ji=i1,i2 323 tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_b(ji,jj,jk) 324 END DO 325 END DO 326 END DO 327 328 ! Extrapolate thicknesses in partial bottom cells: 329 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 330 IF (ln_zps) THEN 331 DO jj=j1,j2 332 DO ji=i1,i2 333 jk = mbkt(ji,jj) 334 tabres(ji,jj,jk,jpts+1) = 0._wp 335 END DO 336 END DO 337 END IF 338 339 ! Save ssh at last level: 340 IF (.NOT.ln_linssh) THEN 341 tabres(i1:i2,j1:j2,k2,jpts+1) = sshb(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1) 342 ELSE 343 tabres(i1:i2,j1:j2,k2,jpts+1) = 0._wp 344 END IF 326 345 # endif 327 346 … … 329 348 ! 330 349 # if defined key_vertical 350 351 IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp 352 331 353 DO jj=j1,j2 332 354 DO ji=i1,i2 333 355 tabres_child(ji,jj,:,:) = 0._wp 334 N_in = 0 335 DO jk=k1,k2 !k2 = jpk of parent grid 336 IF (tabres(ji,jj,jk,n2) == 0) EXIT 337 N_in = N_in + 1 356 N_in = mbkt_parent(ji,jj) 357 zhtot = 0._wp 358 DO jk=1,N_in !k2 = jpk of parent grid 359 IF (jk==N_in) THEN 360 h_in(jk) = ht0_parent(ji,jj) + tabres(ji,jj,k2,n2) - zhtot 361 ELSE 362 h_in(jk) = tabres(ji,jj,jk,n2) 363 ENDIF 364 zhtot = zhtot + h_in(jk) 338 365 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1) 339 h_in(N_in) = tabres(ji,jj,jk,n2)340 366 END DO 341 367 N_out = 0 … … 345 371 h_out(jk) = e3t_b(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above 346 372 ENDDO 373 374 ! Account for small differences in free-surface 375 IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 376 h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 377 ELSE 378 h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 379 ENDIF 347 380 IF (N_in*N_out > 0) THEN 348 381 CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) … … 356 389 DO jk=1,jpkm1 357 390 # if defined key_vertical 358 tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres_child(ji,jj,jk,1:jpts)391 tsbdiff(ji,jj,jk,1:jpts) = (tsb(ji,jj,jk,1:jpts) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk) 359 392 # else 360 tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres(ji,jj,jk,1:jpts)393 tsbdiff(ji,jj,jk,1:jpts) = (tsb(ji,jj,jk,1:jpts) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk) 361 394 # endif 362 395 ENDDO … … 427 460 428 461 ! sponge parameters 429 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff462 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot 430 463 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff 431 464 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff … … 434 467 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 435 468 REAL(wp), DIMENSION(1:jpk) :: h_out 436 INTEGER ::N_in, N_out469 INTEGER ::N_in, N_out 437 470 !!--------------------------------------------- 438 471 ! … … 449 482 END DO 450 483 484 # if defined key_vertical 485 ! Extrapolate thicknesses in partial bottom cells: 486 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 487 IF (ln_zps) THEN 488 DO jj=j1,j2 489 DO ji=i1,i2 490 jk = mbku(ji,jj) 491 tabres(ji,jj,jk,m2) = 0._wp 492 END DO 493 END DO 494 END IF 495 ! Save ssh at last level: 496 tabres(i1:i2,j1:j2,k2,m2) = 0._wp 497 IF (.NOT.ln_linssh) THEN 498 ! This vertical sum below should be replaced by the sea-level at U-points (optimization): 499 DO jk=1,jpk 500 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3u_b(i1:i2,j1:j2,jk) * umask(i1:i2,j1:j2,jk) 501 END DO 502 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2) 503 END IF 504 # endif 505 451 506 ELSE 452 507 453 508 # if defined key_vertical 454 tabres_child(:,:,:) = 0._wp 509 IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 510 455 511 DO jj=j1,j2 456 512 DO ji=i1,i2 457 N_in = 0 458 DO jk=k1,k2 459 IF (tabres(ji,jj,jk,m2) == 0) EXIT 460 N_in = N_in + 1 513 tabres_child(ji,jj,:) = 0._wp 514 N_in = mbku_parent(ji,jj) 515 zhtot = 0._wp 516 DO jk=1,N_in 517 IF (jk==N_in) THEN 518 h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 519 ELSE 520 h_in(jk) = tabres(ji,jj,jk,m2) 521 ENDIF 522 zhtot = zhtot + h_in(jk) 461 523 tabin(jk) = tabres(ji,jj,jk,m1) 462 h_in(N_in) = tabres(ji,jj,jk,m2) 463 ENDDO 464 ! 465 IF (N_in == 0) THEN 466 tabres_child(ji,jj,:) = 0. 467 CYCLE 468 ENDIF 469 470 N_out = 0 471 DO jk=1,jpk 472 if (umask(ji,jj,jk) == 0) EXIT 473 N_out = N_out + 1 474 h_out(N_out) = e3u_b(ji,jj,jk) 475 ENDDO 476 477 IF (N_out == 0) THEN 478 tabres_child(ji,jj,:) = 0. 479 CYCLE 480 ENDIF 481 482 IF (N_in * N_out > 0) THEN 483 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 484 if (h_diff < -1.e4) then 485 print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 486 endif 487 ENDIF 488 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 489 524 ENDDO 525 ! 526 N_out = 0 527 DO jk=1,jpk 528 IF (umask(ji,jj,jk) == 0) EXIT 529 N_out = N_out + 1 530 h_out(N_out) = e3u_b(ji,jj,jk) 531 ENDDO 532 533 ! Account for small differences in free-surface 534 IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 535 h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 536 ELSE 537 h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 538 ENDIF 539 540 IF (N_in * N_out > 0) THEN 541 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 542 ENDIF 490 543 ENDDO 491 544 ENDDO … … 580 633 ! 581 634 INTEGER :: ji, jj, jk, imax 582 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff635 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot 583 636 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff 584 637 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff … … 601 654 END DO 602 655 END DO 656 657 # if defined key_vertical 658 ! Extrapolate thicknesses in partial bottom cells: 659 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 660 IF (ln_zps) THEN 661 DO jj=j1,j2 662 DO ji=i1,i2 663 jk = mbkv(ji,jj) 664 tabres(ji,jj,jk,m2) = 0._wp 665 END DO 666 END DO 667 END IF 668 ! Save ssh at last level: 669 tabres(i1:i2,j1:j2,k2,m2) = 0._wp 670 IF (.NOT.ln_linssh) THEN 671 ! This vertical sum below should be replaced by the sea-level at V-points (optimization): 672 DO jk=1,jpk 673 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3v_b(i1:i2,j1:j2,jk) * vmask(i1:i2,j1:j2,jk) 674 END DO 675 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2) 676 END IF 677 # endif 678 603 679 ELSE 604 680 605 681 # if defined key_vertical 606 tabres_child(:,:,:) = 0._wp682 IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 607 683 DO jj=j1,j2 608 684 DO ji=i1,i2 609 N_in = 0 610 DO jk=k1,k2 611 IF (tabres(ji,jj,jk,m2) == 0) EXIT 612 N_in = N_in + 1 685 tabres_child(ji,jj,:) = 0._wp 686 N_in = mbkv_parent(ji,jj) 687 zhtot = 0._wp 688 DO jk=1,N_in 689 IF (jk==N_in) THEN 690 h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 691 ELSE 692 h_in(jk) = tabres(ji,jj,jk,m2) 693 ENDIF 694 zhtot = zhtot + h_in(jk) 613 695 tabin(jk) = tabres(ji,jj,jk,m1) 614 h_in(N_in) = tabres(ji,jj,jk,m2) 615 ENDDO 696 ENDDO 697 ! 698 N_out = 0 699 DO jk=1,jpk 700 IF (vmask(ji,jj,jk) == 0) EXIT 701 N_out = N_out + 1 702 h_out(N_out) = e3v_b(ji,jj,jk) 703 ENDDO 704 705 ! Account for small differences in free-surface 706 IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 707 h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 708 ELSE 709 h_in(1) = h_in(1) - ( sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 710 ENDIF 616 711 617 IF (N_in == 0) THEN 618 tabres_child(ji,jj,:) = 0. 619 CYCLE 620 ENDIF 621 622 N_out = 0 623 DO jk=1,jpk 624 if (vmask(ji,jj,jk) == 0) EXIT 625 N_out = N_out + 1 626 h_out(N_out) = e3v_b(ji,jj,jk) 627 ENDDO 628 629 IF (N_in * N_out > 0) THEN 630 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 631 if (h_diff < -1.e4) then 632 print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 633 endif 634 ENDIF 635 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 712 IF (N_in * N_out > 0) THEN 713 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 714 ENDIF 636 715 ENDDO 637 716 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
