Changeset 8861
- Timestamp:
- 2017-11-30T16:58:38+01:00 (7 years ago)
- Location:
- branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC
- Files:
-
- 2 edited
-
agrif_opa_interp.F90 (modified) (9 diffs)
-
agrif_user.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90
r8822 r8861 576 576 !!---------------------------------------------------------------------- 577 577 ! 578 return579 578 580 579 zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp ) … … 602 601 ! 603 602 INTEGER :: ji, jj, jk, jn, iref, jref ! dummy loop indices 604 INTEGER :: imin, imax, jmin, jmax 603 INTEGER :: imin, imax, jmin, jmax, N_in, N_out 605 604 REAL(wp) :: zrhox , zalpha1, zalpha2, zalpha3 606 605 REAL(wp) :: zalpha4, zalpha5, zalpha6, zalpha7 607 606 LOGICAL :: western_side, eastern_side,northern_side,southern_side 608 ! VERTICAL REFINEMENT BEGIN609 607 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child 610 608 REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 611 609 REAL(wp) :: h_in(k1:k2) 612 610 REAL(wp) :: h_out(1:jpk) 613 INTEGER :: N_in, N_out 614 REAL(wp) :: h_diff 615 REAL(wp) :: zrhoxy 616 ! VERTICAL REFINEMENT END 611 REAL(wp) :: h_diff, zrhoxy 617 612 618 613 zrhoxy = Agrif_rhox()*Agrif_rhoy() 619 620 621 622 614 IF (before) THEN 623 IF(Agrif_UseSpecialValue) THEN 624 ! Agrif_SpecialValue = -999._wp 625 Agrif_SpecialValue = 0._wp 626 ELSE 627 Agrif_SpecialValue = 0._wp 628 ENDIF 629 DO jn = n1,n2-1 630 DO jk=k1,k2 631 DO jj=j1,j2 632 DO ji=i1,i2 633 ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk)! * tmask(ji,jj,jk) - & 634 ! & (tmask(ji,jj,jk)-1) * Agrif_SpecialValue 635 END DO 636 END DO 637 END DO 638 END DO 615 DO jn = n1,n2-1 639 616 DO jk=k1,k2 640 617 DO jj=j1,j2 641 DO ji=i1,i2 642 ptab(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 643 END DO 644 END DO 645 END DO 618 DO ji=i1,i2 619 ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 620 END DO 621 END DO 622 END DO 623 END DO 624 DO jk=k1,k2 625 DO jj=j1,j2 626 DO ji=i1,i2 627 ptab(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 628 END DO 629 END DO 630 END DO 646 631 ELSE 647 632 western_side = (nb == 1).AND.(ndir == 1) … … 649 634 southern_side = (nb == 2).AND.(ndir == 1) 650 635 northern_side = (nb == 2).AND.(ndir == 2) 651 Agrif_SpecialValue = 0._wp !reset now interpolation is done652 ! VERTICAL REFINEMENT BEGIN653 636 #ifdef key_vertical 654 637 ptab_child(:,:,:,:) = 0. 655 do jj=j1,j2 656 do ji=i1,i2 657 iref = ji 658 jref = jj 659 if(western_side) iref=2 660 if(eastern_side) iref=nlci-1 661 if(southern_side) jref=2 662 if(northern_side) jref=nlcj-1 663 N_in = 0 664 DO jk=k1,k2 !k2 = jpk of parent grid 665 IF (ptab(ji,jj,jk,n2) == 0) EXIT 666 N_in = N_in + 1 667 tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)/(ptab(ji,jj,jk,n2)) 668 h_in(N_in) = ptab(ji,jj,jk,n2)/(e1e2t(ji,jj)*zrhoxy) 669 END DO 670 N_out = 0 671 DO jk=1,jpk ! jpk of child grid 672 IF (tmask(iref,jref,jk) == 0) EXIT ! TODO: Will not work with ISF. !This doesn't seem to work at the moment in GYRE but is OK in overflow model 673 N_out = N_out + 1 674 h_out(jk) = e3t_n(iref,jref,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above 675 ENDDO 676 IF (N_in > 0) THEN 677 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 678 ! if (h_diff > 0) then 679 ! h_in(N_in+1) = h_diff 680 ! N_in = N_in + 1 681 ! else 682 ! h_out(N_out+1) = -h_diff 683 ! N_out = N_out + 1 684 ! endif 685 ptab(ji,jj,k2,:) = ptab(ji,jj,k2-1,:) !what is this line for????? 686 do jn=1,jpts 687 call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 688 enddo 689 ! if (abs(h_diff) > 1000.) then 690 ! do jn=1,jpts 691 ! do jk=1,N_out 692 ! print *,'AVANT APRES = ',ji,jj,jk,N_out,ptab(ji,jj,jk,jn),ptab_child(ji,jj,jk,jn) 693 ! enddo 694 ! enddo 695 ! endif 696 ENDIF 697 enddo 698 enddo 638 DO jj=j1,j2 639 DO ji=i1,i2 640 iref = ji 641 jref = jj 642 if(western_side) iref=2 643 if(eastern_side) iref=nlci-1 644 if(southern_side) jref=2 645 if(northern_side) jref=nlcj-1 646 N_in = 0 647 DO jk=k1,k2 !k2 = jpk of parent grid 648 IF (ptab(ji,jj,jk,n2) == 0) EXIT 649 N_in = N_in + 1 650 tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)/(ptab(ji,jj,jk,n2)) 651 h_in(N_in) = ptab(ji,jj,jk,n2)/(e1e2t(ji,jj)*zrhoxy) 652 END DO 653 N_out = 0 654 DO jk=1,jpk ! jpk of child grid 655 IF (tmask(iref,jref,jk) == 0) EXIT 656 N_out = N_out + 1 657 h_out(jk) = e3t_n(iref,jref,jk) 658 ENDDO 659 IF (N_in > 0) THEN 660 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 661 DO jn=1,jpts 662 call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 663 ENDDO 664 ENDIF 665 ENDDO 666 ENDDO 699 667 #else 700 do jk=k1,k2 701 do jj=j1,j2 702 do ji=i1,i2 703 ! This will be slow - Is there a way to speed it up and avoid divide by zero? 704 IF (ptab(ji,jj,jk,n2) .NE. 0) THEN 705 ptab_child(ji,jj,jk,n1:n2-1) = ptab(ji,jj,jk,n1:n2-1)/ptab(ji,jj,jk,n2) 706 ELSE 707 ptab_child(ji,jj,jk,n1:n2-1) = 0._wp 708 ENDIF 709 enddo 710 enddo 711 enddo 668 DO jk=k1,k2 669 DO jj=j1,j2 670 DO ji=i1,i2 671 ptab_child(ji,jj,jk,n1:n2-1) = ptab(ji,jj,jk,n1:n2-1)/(zrhoxy*e1e2t(ji,jj)*e3t_n(ji,jj,jk) 672 ENDDO 673 ENDDO 674 ENDDO 712 675 #endif 713 714 715 ! VERTICAL REFINEMENT END716 717 676 ! 718 677 zrhox = Agrif_Rhox() … … 996 955 INTEGER :: ji,jj,jk 997 956 REAL(wp) :: zrhox 998 ! VERTICAL REFINEMENT BEGIN999 957 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child 1000 958 REAL(wp), DIMENSION(k1:k2) :: tabin … … 1006 964 INTEGER :: jref 1007 965 1008 ! VERTICAL REFINEMENT END1009 966 !!--------------------------------------------- 1010 967 ! … … 1029 986 ptab_child(:,:,:) = 0. 1030 987 #ifdef key_vertical 1031 ! VERTICAL REFINEMENT BEGIN1032 988 southern_side = (nb == 2).AND.(ndir == 1) 1033 989 northern_side = (nb == 2).AND.(ndir == 2) … … 1061 1017 CYCLE 1062 1018 ENDIF 1063 1064 ! IF (N_in * N_out > 0) THEN1065 ! h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))1066 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly1067 ! if (h_diff < 0.) then1068 ! print *,'CHECK YOUR BATHY interpvn...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))1069 ! stop1070 ! endif1071 ! ENDIF1072 1019 call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ptab_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 1073 1074 1020 enddo 1075 1021 enddo 1076 ! in the following1077 ! remove division of va by fs e3v, zrhox and e1v (already done)1078 1022 DO jk=1,jpkm1 1079 1023 DO jj=j1,j2 … … 1081 1025 END DO 1082 1026 END DO 1083 ! VERTICAL REFINEMENT END1084 1027 #else 1085 1028 DO jk=1,jpkm1 … … 1448 1391 # if defined key_zdftke 1449 1392 1450 SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, before )1393 SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before ) 1451 1394 !!---------------------------------------------------------------------- 1452 1395 !! *** ROUTINE interavm *** 1453 1396 !!---------------------------------------------------------------------- 1454 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2 1455 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2 ), INTENT(inout) :: ptab1397 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, m1, m2 1398 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 1456 1399 LOGICAL , INTENT(in ) :: before 1400 REAL(wp), DIMENSION(k1:k2) :: tabin 1401 REAL(wp) :: h_in(k1:k2) 1402 REAL(wp) :: h_out(1:jpk) 1403 REAL(wp) :: zrhoxy 1404 INTEGER :: N_in, N_out, ji, jj, jk 1457 1405 !!---------------------------------------------------------------------- 1458 1406 ! 1459 IF( before ) THEN 1460 ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2) 1461 ELSE 1462 avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 1407 zrhoxy = Agrif_rhox()*Agrif_rhoy() 1408 IF (before) THEN 1409 DO jk=k1,k2 1410 DO jj=j1,j2 1411 DO ji=i1,i2 1412 ptab(ji,jj,jk,1) = avm_k(ji,jj,jk) 1413 END DO 1414 END DO 1415 END DO 1416 #ifdef key_vertical 1417 DO jk=k1,k2 1418 DO jj=j1,j2 1419 DO ji=i1,i2 1420 ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e1e2t(ji,jj) * e3w_n(ji,jj,jk) 1421 END DO 1422 END DO 1423 END DO 1424 #else 1425 ptab(i1:i2,j1:j2,k1:k2,2) = 0._wp 1426 #endif 1427 ELSE 1428 #ifdef key_vertical 1429 avm_k(i1:i2,j1:j2,1:jpk) = 0. 1430 DO jj=j1,j2 1431 DO ji=i1,i2 1432 N_in = 0 1433 DO jk=k1,k2 !k2 = jpk of parent grid 1434 IF (ptab(ji,jj,jk,2) == 0) EXIT 1435 N_in = N_in + 1 1436 tabin(jk) = ptab(ji,jj,jk,1) 1437 h_in(N_in) = ptab(ji,jj,jk,2)/(e1e2t(ji,jj)*zrhoxy) 1438 END DO 1439 N_out = 0 1440 DO jk=1,jpk ! jpk of child grid 1441 IF (wmask(ji,jj,jk) == 0) EXIT 1442 N_out = N_out + 1 1443 h_out(jk) = e3t_n(ji,jj,jk) 1444 ENDDO 1445 IF (N_in > 0) THEN 1446 CALL reconstructandremap(tabin(1:N_in),h_in,avm_k(ji,jj,1:N_out),h_out,N_in,N_out) 1447 ENDIF 1448 ENDDO 1449 ENDDO 1450 #else 1451 avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1) 1452 #endif 1463 1453 ENDIF 1464 1454 ! -
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_user.F90
r8835 r8861 374 374 CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/), en_id) 375 375 CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/),avt_id) 376 CALL agrif_declare_variable((/2,2,0 /),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/jpi,jpj,jpk/),avm_id)376 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/jpi,jpj,jpk,2/),avm_id) 377 377 # endif 378 378 … … 757 757 ENDIF 758 758 759 AGRIF_Parent(ln_vertical)=ln_vertical ! Will this work? 760 759 761 CALL Agrif_Update_trc(0) 760 762 ! … … 832 834 INTEGER :: ios ! Local integer output status for namelist read 833 835 INTEGER :: iminspon 834 NAMELIST/namagrif/ nn_cln_update, rn_sponge_tra, rn_sponge_dyn, ln_spc_dyn, ln_chk_bathy 836 NAMELIST/namagrif/ nn_cln_update, rn_sponge_tra, rn_sponge_dyn, ln_spc_dyn, ln_chk_bathy, ln_vertical 835 837 !!-------------------------------------------------------------------------------------- 836 838 ! … … 854 856 WRITE(numout,*) ' use special values for dynamics ln_spc_dyn = ', ln_spc_dyn 855 857 WRITE(numout,*) ' check bathymetry ln_chk_bathy = ', ln_chk_bathy 858 WRITE(numout,*) ' use vertical interpolation ln_vertical = ', ln_vertical 856 859 WRITE(numout,*) 857 860 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
