Changeset 9007
- Timestamp:
- 2017-12-13T12:32:38+01:00 (7 years ago)
- Location:
- branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/NST_SRC
- Files:
-
- 8 edited
-
agrif_oce.F90 (modified) (3 diffs)
-
agrif_opa_interp.F90 (modified) (17 diffs)
-
agrif_opa_sponge.F90 (modified) (5 diffs)
-
agrif_opa_update.F90 (modified) (16 diffs)
-
agrif_top_interp.F90 (modified) (7 diffs)
-
agrif_top_sponge.F90 (modified) (1 diff)
-
agrif_top_update.F90 (modified) (2 diffs)
-
agrif_user.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90
r8993 r9007 17 17 18 18 PUBLIC agrif_oce_alloc ! routine called by nemo_init in nemogcm.F90 19 19 #if defined key_vertical 20 PUBLIC reconstructandremap ! remapping routine 21 #endif 20 22 ! !!* Namelist namagrif: AGRIF parameters 21 23 LOGICAL , PUBLIC :: ln_spc_dyn = .FALSE. !: … … 24 26 REAL(wp), PUBLIC :: rn_sponge_dyn = 2800. !: sponge coeff. for dynamics 25 27 LOGICAL , PUBLIC :: ln_chk_bathy = .FALSE. !: check of parent bathymetry 26 LOGICAL , PUBLIC :: lk_agrif_clp = .FALSE. !: Flag to retrieve clamped open boundaries 27 28 LOGICAL , PUBLIC :: lk_agrif_clp = .TRUE. !: Force clamped bcs 28 29 ! !!! OLD namelist names 29 30 REAL(wp), PUBLIC :: visc_tra !: sponge coeff. for tracers … … 101 102 END FUNCTION agrif_oce_alloc 102 103 104 #if defined key_vertical 105 SUBROUTINE reconstructandremap(tabin,hin,tabout,hout,N,Nout) 106 !!---------------------------------------------------------------------- 107 !! *** FUNCTION reconstructandremap *** 108 !!---------------------------------------------------------------------- 109 IMPLICIT NONE 110 INTEGER N, Nout 111 REAL(wp) tabin(N), tabout(Nout) 112 REAL(wp) hin(N), hout(Nout) 113 REAL(wp) coeffremap(N,3),zwork(N,3) 114 REAL(wp) zwork2(N+1,3) 115 INTEGER jk 116 DOUBLE PRECISION, PARAMETER :: dsmll=1.0d-8 117 REAL(wp) q,q01,q02,q001,q002,q0 118 REAL(wp) z_win(1:N+1), z_wout(1:Nout+1) 119 REAL(wp),PARAMETER :: dpthin = 1.D-3 120 INTEGER :: k1, kbox, ktop, ka, kbot 121 REAL(wp) :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop 122 123 z_win(1)=0.; z_wout(1)= 0. 124 DO jk=1,N 125 z_win(jk+1)=z_win(jk)+hin(jk) 126 ENDDO 127 128 DO jk=1,Nout 129 z_wout(jk+1)=z_wout(jk)+hout(jk) 130 ENDDO 131 132 DO jk=2,N 133 zwork(jk,1)=1./(hin(jk-1)+hin(jk)) 134 ENDDO 135 136 DO jk=2,N-1 137 q0 = 1./(hin(jk-1)+hin(jk)+hin(jk+1)) 138 zwork(jk,2)=hin(jk-1)+2.*hin(jk)+hin(jk+1) 139 zwork(jk,3)=q0 140 ENDDO 141 142 DO jk= 2,N 143 zwork2(jk,1)=zwork(jk,1)*(tabin(jk)-tabin(jk-1)) 144 ENDDO 145 146 coeffremap(:,1) = tabin(:) 147 148 DO jk=2,N-1 149 q001 = hin(jk)*zwork2(jk+1,1) 150 q002 = hin(jk)*zwork2(jk,1) 151 IF (q001*q002 < 0) then 152 q001 = 0. 153 q002 = 0. 154 ENDIF 155 q=zwork(jk,2) 156 q01=q*zwork2(jk+1,1) 157 q02=q*zwork2(jk,1) 158 IF (abs(q001) > abs(q02)) q001 = q02 159 IF (abs(q002) > abs(q01)) q002 = q01 160 161 q=(q001-q002)*zwork(jk,3) 162 q001=q001-q*hin(jk+1) 163 q002=q002+q*hin(jk-1) 164 165 coeffremap(jk,3)=coeffremap(jk,1)+q001 166 coeffremap(jk,2)=coeffremap(jk,1)-q002 167 168 zwork2(jk,1)=(2.*q001-q002)**2 169 zwork2(jk,2)=(2.*q002-q001)**2 170 ENDDO 171 172 DO jk=1,N 173 IF(jk.EQ.1 .OR. jk.EQ.N .OR. hin(jk).LE.dpthin) THEN 174 coeffremap(jk,3) = coeffremap(jk,1) 175 coeffremap(jk,2) = coeffremap(jk,1) 176 zwork2(jk,1) = 0. 177 zwork2(jk,2) = 0. 178 ENDIF 179 ENDDO 180 181 DO jk=2,N 182 q002=max(zwork2(jk-1,2),dsmll) 183 q001=max(zwork2(jk,1),dsmll) 184 zwork2(jk,3)=(q001*coeffremap(jk-1,3)+q002*coeffremap(jk,2))/(q001+q002) 185 ENDDO 186 187 zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3) 188 zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3) 189 190 DO jk=1,N 191 q01=zwork2(jk+1,3)-coeffremap(jk,1) 192 q02=coeffremap(jk,1)-zwork2(jk,3) 193 q001=2.*q01 194 q002=2.*q02 195 IF (q01*q02<0) then 196 q01=0. 197 q02=0. 198 ELSEIF (abs(q01)>abs(q002)) then 199 q01=q002 200 ELSEIF (abs(q02)>abs(q001)) then 201 q02=q001 202 ENDIF 203 coeffremap(jk,2)=coeffremap(jk,1)-q02 204 coeffremap(jk,3)=coeffremap(jk,1)+q01 205 ENDDO 206 207 zbot=0.0 208 kbot=1 209 DO jk=1,Nout 210 ztop=zbot !top is bottom of previous layer 211 ktop=kbot 212 IF (ztop.GE.z_win(ktop+1)) then 213 ktop=ktop+1 214 ENDIF 215 216 zbot=z_wout(jk+1) 217 zthk=zbot-ztop 218 219 IF(zthk.GT.dpthin .AND. ztop.LT.z_wout(Nout+1)) THEN 220 221 kbot=ktop 222 DO while (z_win(kbot+1).lt.zbot.and.kbot.lt.N) 223 kbot=kbot+1 224 ENDDO 225 zbox=zbot 226 DO k1= jk+1,Nout 227 IF (z_wout(k1+1)-z_wout(k1).GT.dpthin) THEN 228 exit !thick layer 229 ELSE 230 zbox=z_wout(k1+1) !include thin adjacent layers 231 IF(zbox.EQ.z_wout(Nout+1)) THEN 232 exit !at bottom 233 ENDIF 234 ENDIF 235 ENDDO 236 zthk=zbox-ztop 237 238 kbox=ktop 239 DO while (z_win(kbox+1).lt.zbox.and.kbox.lt.N) 240 kbox=kbox+1 241 ENDDO 242 243 IF(ktop.EQ.kbox) THEN 244 IF(z_wout(jk).NE.z_win(kbox).OR.z_wout(jk+1).NE.z_win(kbox+1)) THEN 245 IF(hin(kbox).GT.dpthin) THEN 246 q001 = (zbox-z_win(kbox))/hin(kbox) 247 q002 = (ztop-z_win(kbox))/hin(kbox) 248 q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002) 249 q02=q01-1.+(q001+q002) 250 q0=1.-q01-q02 251 ELSE 252 q0 = 1.0 253 q01 = 0. 254 q02 = 0. 255 ENDIF 256 tabout(jk)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3) 257 ELSE 258 tabout(jk) = tabin(kbox) 259 ENDIF 260 ELSE 261 IF(ktop.LE.jk .AND. kbox.GE.jk) THEN 262 ka = jk 263 ELSEIF (kbox-ktop.GE.3) THEN 264 ka = (kbox+ktop)/2 265 ELSEIF (hin(ktop).GE.hin(kbox)) THEN 266 ka = ktop 267 ELSE 268 ka = kbox 269 ENDIF !choose ka 270 271 offset=coeffremap(ka,1) 272 273 qtop = z_win(ktop+1)-ztop !partial layer thickness 274 IF(hin(ktop).GT.dpthin) THEN 275 q=(ztop-z_win(ktop))/hin(ktop) 276 q01=q*(q-1.) 277 q02=q01+q 278 q0=1-q01-q02 279 ELSE 280 q0 = 1. 281 q01 = 0. 282 q02 = 0. 283 ENDIF 284 285 tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop 286 287 DO k1= ktop+1,kbox-1 288 tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1) 289 ENDDO !k1 290 291 qbot = zbox-z_win(kbox) !partial layer thickness 292 IF(hin(kbox).GT.dpthin) THEN 293 q=qbot/hin(kbox) 294 q01=(q-1.)**2 295 q02=q01-1.+q 296 q0=1-q01-q02 297 ELSE 298 q0 = 1.0 299 q01 = 0. 300 q02 = 0. 301 ENDIF 302 303 tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot 304 305 rpsum=1.0d0/zthk 306 tabout(jk)=offset+tsum*rpsum 307 308 ENDIF !single or multiple layers 309 ELSE 310 IF (jk==1) THEN 311 write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(jk+1),hout(1) 312 ENDIF 313 tabout(jk) = tabout(jk-1) 314 315 ENDIF !normal:thin layer 316 ENDDO !jk 317 318 return 319 end subroutine reconstructandremap 320 #endif 321 103 322 #endif 104 323 !!====================================================================== -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90
r8993 r9007 118 118 ! 119 119 IF (.NOT.lk_agrif_clp) THEN 120 DO jk=1,jpkm1 ! Smooth120 DO jk=1,jpkm1 ! Smooth 121 121 DO jj=j1,j2 122 122 ua(2,jj,jk) = 0.25_wp*(ua(1,jj,jk)+2._wp*ua(2,jj,jk)+ua(3,jj,jk)) … … 124 124 END DO 125 125 END DO 126 END 126 ENDIF 127 127 ! 128 128 zub(2,:) = 0._wp ! Correct transport … … 189 189 190 190 IF (.NOT.lk_agrif_clp) THEN 191 DO jk = 1, jpkm1 ! Smooth191 DO jk = 1, jpkm1 ! Smooth 192 192 DO jj = j1, j2 193 193 ua(nlci-2,jj,jk) = 0.25_wp * umask(nlci-2,jj,jk) & … … 196 196 END DO 197 197 ENDIF 198 199 198 zub(nlci-2,:) = 0._wp ! Correct transport 200 199 DO jk = 1, jpkm1 … … 331 330 ! 332 331 IF (.NOT.lk_agrif_clp) THEN 333 DO jk = 1, jpkm1 ! Smooth332 DO jk = 1, jpkm1 ! Smooth 334 333 DO ji = i1, i2 335 334 va(ji,nlcj-2,jk) = 0.25_wp * vmask(ji,nlcj-2,jk) & … … 337 336 END DO 338 337 END DO 339 END IF338 END IF 340 339 ! 341 340 zvb(:,nlcj-2) = 0._wp ! Correct transport … … 614 613 INTEGER , INTENT(in ) :: nb , ndir 615 614 ! 616 INTEGER :: ji, jj, jk, jn ! dummy loop indices617 INTEGER :: imin, imax, jmin, jmax 615 INTEGER :: ji, jj, jk, jn, iref, jref ! dummy loop indices 616 INTEGER :: imin, imax, jmin, jmax, N_in, N_out 618 617 REAL(wp) :: zrhox , zalpha1, zalpha2, zalpha3 619 618 REAL(wp) :: zalpha4, zalpha5, zalpha6, zalpha7 620 LOGICAL :: western_side, eastern_side,northern_side,southern_side 621 !!---------------------------------------------------------------------- 622 ! 619 LOGICAL :: western_side, eastern_side,northern_side,southern_side 620 ! vertical interpolation: 621 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child 622 REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 623 REAL(wp), DIMENSION(k1:k2) :: h_in 624 REAL(wp), DIMENSION(1:jpk) :: h_out(1:jpk) 625 REAL(wp) :: h_diff, zrhoxy 626 627 zrhoxy = Agrif_rhox()*Agrif_rhoy() 623 628 IF (before) THEN 624 ptab(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2) 625 ELSE 629 DO jn = 1,jpts 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) 634 END DO 635 END DO 636 END DO 637 END DO 638 639 # if defined key_vertical 640 DO jk=k1,k2 641 DO jj=j1,j2 642 DO ji=i1,i2 643 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 644 END DO 645 END DO 646 END DO 647 # endif 648 ELSE 649 650 western_side = (nb == 1).AND.(ndir == 1) 651 eastern_side = (nb == 1).AND.(ndir == 2) 652 southern_side = (nb == 2).AND.(ndir == 1) 653 northern_side = (nb == 2).AND.(ndir == 2) 654 655 # if defined key_vertical 656 DO jj=j1,j2 657 DO ji=i1,i2 658 iref = ji 659 jref = jj 660 if(western_side) iref=MAX(2,ji) 661 if(eastern_side) iref=MIN(nlci-1,ji) 662 if(southern_side) jref=MAX(2,jj) 663 if(northern_side) jref=MIN(nlcj-1,jj) 664 N_in = 0 665 DO jk=k1,k2 !k2 = jpk of parent grid 666 IF (ptab(ji,jj,jk,n2) == 0) EXIT 667 N_in = N_in + 1 668 tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1) 669 h_in(N_in) = ptab(ji,jj,jk,n2) 670 END DO 671 N_out = 0 672 DO jk=1,jpk ! jpk of child grid 673 IF (tmask(iref,jref,jk) == 0) EXIT 674 N_out = N_out + 1 675 h_out(jk) = e3t_n(iref,jref,jk) 676 ENDDO 677 IF (N_in > 0) THEN 678 DO jn=1,jpts 679 call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 680 ENDDO 681 ENDIF 682 ENDDO 683 ENDDO 684 # else 685 ptab_child(i1:i2,j1:j2,1:jpk,1:jpts) = ptab(i1:i2,j1:j2,1:jpk,1:jpts) 686 # endif 687 ! 626 688 IF (lk_agrif_clp) THEN 627 689 DO jn = 1, jpts … … 629 691 DO ji = i1,i2 630 692 DO jj = j1,j2 631 tsa(ji,jj,jk,jn) = ptab (ji,jj,jk,jn)693 tsa(ji,jj,jk,jn) = ptab_child(ji,jj,jk,jn) 632 694 END DO 633 695 END DO … … 636 698 return 637 699 ENDIF 638 !639 western_side = (nb == 1).AND.(ndir == 1)640 eastern_side = (nb == 1).AND.(ndir == 2)641 southern_side = (nb == 2).AND.(ndir == 1)642 northern_side = (nb == 2).AND.(ndir == 2)643 700 ! 644 701 zrhox = Agrif_Rhox() … … 667 724 IF( eastern_side ) THEN 668 725 DO jn = 1, jpts 669 tsa(nlci,j1:j2, k1:k2,jn) = zalpha1 * ptab(nlci,j1:j2,k1:k2,jn) + zalpha2 * ptab(nlci-1,j1:j2,k1:k2,jn)726 tsa(nlci,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(nlci,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(nlci-1,j1:j2,1:jpk,jn) 670 727 DO jk = 1, jpkm1 671 728 DO jj = jmin,jmax … … 687 744 IF( northern_side ) THEN 688 745 DO jn = 1, jpts 689 tsa(i1:i2,nlcj, k1:k2,jn) = zalpha1 * ptab(i1:i2,nlcj,k1:k2,jn) + zalpha2 * ptab(i1:i2,nlcj-1,k1:k2,jn)746 tsa(i1:i2,nlcj,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,nlcj,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,nlcj-1,1:jpk,jn) 690 747 DO jk = 1, jpkm1 691 748 DO ji = imin,imax … … 707 764 IF( western_side ) THEN 708 765 DO jn = 1, jpts 709 tsa(1,j1:j2, k1:k2,jn) = zalpha1 * ptab(1,j1:j2,k1:k2,jn) + zalpha2 * ptab(2,j1:j2,k1:k2,jn)766 tsa(1,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(1,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(2,j1:j2,1:jpk,jn) 710 767 DO jk = 1, jpkm1 711 768 DO jj = jmin,jmax … … 726 783 IF( southern_side ) THEN 727 784 DO jn = 1, jpts 728 tsa(i1:i2,1, k1:k2,jn) = zalpha1 * ptab(i1:i2,1,k1:k2,jn) + zalpha2 * ptab(i1:i2,2,k1:k2,jn)785 tsa(i1:i2,1,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,1,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,2,1:jpk,jn) 729 786 DO jk = 1, jpk 730 787 DO ji=imin,imax … … 747 804 ! East south 748 805 IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN 749 tsa(nlci-1,2,:,:) = ptab (nlci-1,2,:,:)806 tsa(nlci-1,2,:,:) = ptab_child(nlci-1,2,:,1:jpts) 750 807 ENDIF 751 808 ! East north 752 809 IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN 753 tsa(nlci-1,nlcj-1,:,:) = ptab (nlci-1,nlcj-1,:,:)810 tsa(nlci-1,nlcj-1,:,:) = ptab_child(nlci-1,nlcj-1,:,1:jpts) 754 811 ENDIF 755 812 ! West south 756 813 IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN 757 tsa(2,2,:,:) = ptab (2,2,:,:)814 tsa(2,2,:,:) = ptab_child(2,2,:,1:jpts) 758 815 ENDIF 759 816 ! West north 760 817 IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN 761 tsa(2,nlcj-1,:,:) = ptab (2,nlcj-1,:,:)818 tsa(2,nlcj-1,:,:) = ptab_child(2,nlcj-1,:,1:jpts) 762 819 ENDIF 763 820 ! … … 765 822 ! 766 823 END SUBROUTINE interptsn 767 768 824 769 825 SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before, nb, ndir ) … … 794 850 END SUBROUTINE interpsshn 795 851 796 797 SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, before ) 852 SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir ) 798 853 !!---------------------------------------------------------------------- 799 854 !! *** ROUTINE interpun *** 800 !!---------------------------------------------------------------------- 801 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2 802 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 803 LOGICAL , INTENT(in ) :: before 804 ! 805 INTEGER :: ji, jj, jk 806 REAL(wp) :: zrhoy 807 !!---------------------------------------------------------------------- 808 ! 809 IF( before ) THEN 810 DO jk = 1, jpkm1 811 ptab(i1:i2,j1:j2,jk) = e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk) 855 !!--------------------------------------------- 856 !! 857 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 858 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 859 LOGICAL, INTENT(in) :: before 860 INTEGER, INTENT(in) :: nb , ndir 861 !! 862 INTEGER :: ji,jj,jk 863 REAL(wp) :: zrhoy 864 ! vertical interpolation: 865 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 866 REAL(wp), DIMENSION(1:jpk) :: h_out 867 INTEGER :: N_in, N_out, iref 868 REAL(wp) :: h_diff 869 LOGICAL :: western_side, eastern_side 870 !!--------------------------------------------- 871 ! 872 zrhoy = Agrif_rhoy() 873 IF (before) THEN 874 DO jk=1,jpk 875 DO jj=j1,j2 876 DO ji=i1,i2 877 ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) 878 # if defined key_vertical 879 ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)) 880 # endif 881 END DO 882 END DO 812 883 END DO 813 884 ELSE 814 zrhoy = Agrif_Rhoy() 885 zrhoy = Agrif_rhoy() 886 # if defined key_vertical 887 ! VERTICAL REFINEMENT BEGIN 888 western_side = (nb == 1).AND.(ndir == 1) 889 eastern_side = (nb == 1).AND.(ndir == 2) 890 891 DO ji=i1,i2 892 iref = ji 893 IF (western_side) iref = MAX(2,ji) 894 IF (eastern_side) iref = MIN(nlci-2,ji) 895 DO jj=j1,j2 896 N_in = 0 897 DO jk=k1,k2 898 IF (ptab(ji,jj,jk,2) == 0) EXIT 899 N_in = N_in + 1 900 tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 901 h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 902 ENDDO 903 904 IF (N_in == 0) THEN 905 ua(ji,jj,:) = 0._wp 906 CYCLE 907 ENDIF 908 909 N_out = 0 910 DO jk=1,jpk 911 if (umask(iref,jj,jk) == 0) EXIT 912 N_out = N_out + 1 913 h_out(N_out) = e3u_a(iref,jj,jk) 914 ENDDO 915 916 IF (N_out == 0) THEN 917 ua(ji,jj,:) = 0._wp 918 CYCLE 919 ENDIF 920 921 IF (N_in * N_out > 0) THEN 922 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 923 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 924 if (h_diff < -1.e4) then 925 print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 926 ! stop 927 endif 928 ENDIF 929 call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ua(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 930 ENDDO 931 ENDDO 932 933 # else 815 934 DO jk = 1, jpkm1 816 935 DO jj=j1,j2 817 ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk) / ( zrhoy * e2u(i1:i2,jj) * e3u_a(i1:i2,jj,jk) ) 818 END DO 819 END DO 936 ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u_a(i1:i2,jj,jk) ) 937 END DO 938 END DO 939 # endif 940 820 941 ENDIF 821 942 ! 822 943 END SUBROUTINE interpun 823 944 824 825 SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, before ) 945 SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir ) 826 946 !!---------------------------------------------------------------------- 827 947 !! *** ROUTINE interpvn *** 828 948 !!---------------------------------------------------------------------- 829 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2 830 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 831 LOGICAL , INTENT(in ) :: before 832 ! 833 INTEGER :: ji, jj, jk 834 REAL(wp) :: zrhox 835 !!---------------------------------------------------------------------- 949 ! 950 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 951 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 952 LOGICAL, INTENT(in) :: before 953 INTEGER, INTENT(in) :: nb , ndir 954 ! 955 INTEGER :: ji,jj,jk 956 REAL(wp) :: zrhox 957 ! vertical interpolation: 958 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 959 REAL(wp), DIMENSION(1:jpk) :: h_out 960 INTEGER :: N_in, N_out, jref 961 REAL(wp) :: h_diff 962 LOGICAL :: northern_side,southern_side 963 !!--------------------------------------------- 836 964 ! 837 IF( before ) THEN 965 IF (before) THEN 966 DO jk=k1,k2 967 DO jj=j1,j2 968 DO ji=i1,i2 969 ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)*vmask(ji,jj,jk)) 970 # if defined key_vertical 971 ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk) 972 # endif 973 END DO 974 END DO 975 END DO 976 ELSE 977 zrhox = Agrif_rhox() 978 # if defined key_vertical 979 980 southern_side = (nb == 2).AND.(ndir == 1) 981 northern_side = (nb == 2).AND.(ndir == 2) 982 983 DO jj=j1,j2 984 jref = jj 985 IF (southern_side) jref = MAX(2,jj) 986 IF (northern_side) jref = MIN(nlcj-2,jj) 987 DO ji=i1,i2 988 N_in = 0 989 DO jk=k1,k2 990 if (ptab(ji,jj,jk,2) == 0) EXIT 991 N_in = N_in + 1 992 tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 993 h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 994 END DO 995 IF (N_in == 0) THEN 996 va(ji,jj,:) = 0._wp 997 CYCLE 998 ENDIF 999 1000 N_out = 0 1001 DO jk=1,jpk 1002 if (vmask(ji,jref,jk) == 0) EXIT 1003 N_out = N_out + 1 1004 h_out(N_out) = e3v_a(ji,jref,jk) 1005 END DO 1006 IF (N_out == 0) THEN 1007 va(ji,jj,:) = 0._wp 1008 CYCLE 1009 ENDIF 1010 call reconstructandremap(tabin(1:N_in),h_in(1:N_in),va(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 1011 END DO 1012 END DO 1013 # else 838 1014 DO jk = 1, jpkm1 839 ptab(i1:i2,j1:j2,jk) = e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) * vn(i1:i2,j1:j2,jk) 840 END DO 841 ELSE 842 zrhox= Agrif_Rhox() 843 DO jk = 1, jpkm1 844 va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_a(i1:i2,j1:j2,jk) ) 845 END DO 1015 va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_a(i1:i2,j1:j2,jk) ) 1016 END DO 1017 # endif 846 1018 ENDIF 847 1019 ! 848 1020 END SUBROUTINE interpvn 849 850 1021 851 1022 SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before, nb, ndir ) … … 1212 1383 # if defined key_zdftke || defined key_zdfgls 1213 1384 1214 SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, before )1385 SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before ) 1215 1386 !!---------------------------------------------------------------------- 1216 1387 !! *** ROUTINE interavm *** 1217 1388 !!---------------------------------------------------------------------- 1218 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2 1219 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2 ), INTENT(inout) :: ptab1389 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, m1, m2 1390 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 1220 1391 LOGICAL , INTENT(in ) :: before 1392 REAL(wp), DIMENSION(k1:k2) :: tabin 1393 REAL(wp) :: h_in(k1:k2) 1394 REAL(wp) :: h_out(1:jpk) 1395 REAL(wp) :: zrhoxy 1396 INTEGER :: N_in, N_out, ji, jj, jk 1221 1397 !!---------------------------------------------------------------------- 1222 1398 ! 1223 IF( before ) THEN 1224 ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2) 1225 ELSE 1226 avm (i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2) 1399 zrhoxy = Agrif_rhox()*Agrif_rhoy() 1400 IF (before) THEN 1401 DO jk=k1,k2 1402 DO jj=j1,j2 1403 DO ji=i1,i2 1404 ptab(ji,jj,jk,1) = avm_k(ji,jj,jk) 1405 END DO 1406 END DO 1407 END DO 1408 #ifdef key_vertical 1409 DO jk=k1,k2 1410 DO jj=j1,j2 1411 DO ji=i1,i2 1412 ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e1e2t(ji,jj) * e3w_n(ji,jj,jk) 1413 END DO 1414 END DO 1415 END DO 1416 #else 1417 ptab(i1:i2,j1:j2,k1:k2,2) = 0._wp 1418 #endif 1419 ELSE 1420 #ifdef key_vertical 1421 avm_k(i1:i2,j1:j2,1:jpk) = 0. 1422 DO jj=j1,j2 1423 DO ji=i1,i2 1424 N_in = 0 1425 DO jk=k1,k2 !k2 = jpk of parent grid 1426 IF (ptab(ji,jj,jk,2) == 0) EXIT 1427 N_in = N_in + 1 1428 tabin(jk) = ptab(ji,jj,jk,1) 1429 h_in(N_in) = ptab(ji,jj,jk,2)/(e1e2t(ji,jj)*zrhoxy) 1430 END DO 1431 N_out = 0 1432 DO jk=1,jpk ! jpk of child grid 1433 IF (wmask(ji,jj,jk) == 0) EXIT 1434 N_out = N_out + 1 1435 h_out(jk) = e3t_n(ji,jj,jk) 1436 ENDDO 1437 IF (N_in > 0) THEN 1438 CALL reconstructandremap(tabin(1:N_in),h_in,avm_k(ji,jj,1:N_out),h_out,N_in,N_out) 1439 ENDIF 1440 ENDDO 1441 ENDDO 1442 #else 1443 avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1) 1444 #endif 1227 1445 ENDIF 1228 1446 ! -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90
r8993 r9007 39 39 #if defined SPONGE 40 40 !! timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot() 41 !! Assume persistence:41 !! Assume persistence: 42 42 timecoeff = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot()) 43 43 … … 195 195 END SUBROUTINE Agrif_Sponge 196 196 197 198 197 SUBROUTINE interptsn_sponge(tabres,i1,i2,j1,j2,k1,k2,n1,n2,before) 199 198 !!--------------------------------------------- … … 207 206 INTEGER :: iku, ikv 208 207 REAL(wp) :: ztsa, zabe1, zabe2, zbtr 209 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ztu, ztv 210 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2) ::tsbdiff 208 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv 209 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff 210 ! vertical interpolation: 211 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child 212 REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 213 REAL(wp), DIMENSION(k1:k2) :: h_in 214 REAL(wp), DIMENSION(1:jpk) :: h_out 215 INTEGER :: N_in, N_out 216 REAL(wp) :: h_diff 211 217 ! 212 218 IF( before ) THEN 213 tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsb(i1:i2,j1:j2,k1:k2,n1:n2) 219 DO jn = 1, jpts 220 DO jk=k1,k2 221 DO jj=j1,j2 222 DO ji=i1,i2 223 tabres(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) 224 END DO 225 END DO 226 END DO 227 END DO 228 229 # if defined key_vertical 230 DO jk=k1,k2 231 DO jj=j1,j2 232 DO ji=i1,i2 233 tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 234 END DO 235 END DO 236 END DO 237 # endif 238 214 239 ELSE 215 240 ! 216 tsbdiff(:,:,:,:) = tsb(i1:i2,j1:j2,:,:) - tabres(:,:,:,:) 241 # if defined key_vertical 242 tabres_child(:,:,:,:) = 0. 243 DO jj=j1,j2 244 DO ji=i1,i2 245 N_in = 0 246 DO jk=k1,k2 !k2 = jpk of parent grid 247 IF (tabres(ji,jj,jk,n2) == 0) EXIT 248 N_in = N_in + 1 249 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1) 250 h_in(N_in) = tabres(ji,jj,jk,n2) 251 END DO 252 N_out = 0 253 DO jk=1,jpk ! jpk of child grid 254 IF (tmask(ji,jj,jk) == 0) EXIT 255 N_out = N_out + 1 256 h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above 257 ENDDO 258 IF (N_in > 0) THEN 259 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 260 tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for????? 261 DO jn=1,jpts 262 call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 263 ENDDO 264 ENDIF 265 ENDDO 266 ENDDO 267 # endif 268 269 DO jj=j1,j2 270 DO ji=i1,i2 271 DO jk=1,jpkm1 272 # if defined key_vertical 273 tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres_child(ji,jj,jk,1:jpts) 274 # else 275 tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres(ji,jj,jk,1:jpts) 276 # endif 277 ENDDO 278 ENDDO 279 ENDDO 280 217 281 DO jn = 1, jpts 218 282 DO jk = 1, jpkm1 … … 261 325 END SUBROUTINE interptsn_sponge 262 326 263 264 SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2, before) 327 SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before) 265 328 !!--------------------------------------------- 266 329 !! *** ROUTINE interpun_sponge *** 267 330 !!--------------------------------------------- 268 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 269 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2 ), INTENT(inout) :: tabres331 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 332 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres 270 333 LOGICAL, INTENT(in) :: before 271 334 272 INTEGER :: ji,jj,jk 335 INTEGER :: ji,jj,jk,jmax 273 336 274 337 ! sponge parameters 275 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr 276 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ubdiff 277 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff 278 INTEGER :: jmax 338 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff 339 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff 340 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff 341 ! vertical interpolation: 342 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 343 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 344 REAL(wp), DIMENSION(1:jpk) :: h_out 345 INTEGER ::N_in,N_out 279 346 !!--------------------------------------------- 280 347 ! 281 348 IF( before ) THEN 282 tabres = ub(i1:i2,j1:j2,:) 349 DO jk=1,jpkm1 350 DO jj=j1,j2 351 DO ji=i1,i2 352 tabres(ji,jj,jk,m1) = ub(ji,jj,jk) 353 # if defined key_vertical 354 tabres(ji,jj,jk,m2) = e3u_n(ji,jj,jk)*umask(ji,jj,jk) 355 # endif 356 END DO 357 END DO 358 END DO 359 283 360 ELSE 284 ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(:,:,:))*umask(i1:i2,j1:j2,:) 361 362 # if defined key_vertical 363 tabres_child(:,:,:) = 0._wp 364 DO jj=j1,j2 365 DO ji=i1,i2 366 N_in = 0 367 DO jk=k1,k2 368 IF (tabres(ji,jj,jk,m2) == 0) EXIT 369 N_in = N_in + 1 370 tabin(jk) = tabres(ji,jj,jk,m1) 371 h_in(N_in) = tabres(ji,jj,jk,m2) 372 ENDDO 373 ! 374 IF (N_in == 0) THEN 375 tabres_child(ji,jj,:) = 0. 376 CYCLE 377 ENDIF 378 379 N_out = 0 380 DO jk=1,jpk 381 if (umask(ji,jj,jk) == 0) EXIT 382 N_out = N_out + 1 383 h_out(N_out) = e3u_n(ji,jj,jk) 384 ENDDO 385 386 IF (N_out == 0) THEN 387 tabres_child(ji,jj,:) = 0. 388 CYCLE 389 ENDIF 390 391 IF (N_in * N_out > 0) THEN 392 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 393 if (h_diff < -1.e4) then 394 print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 395 endif 396 ENDIF 397 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) 398 399 ENDDO 400 ENDDO 401 402 ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:) 403 #else 404 ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:) 405 #endif 285 406 ! 286 407 DO jk = 1, jpkm1 ! Horizontal slab … … 359 480 END SUBROUTINE interpun_sponge 360 481 361 362 SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2, before,nb,ndir) 482 SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before,nb,ndir) 363 483 !!--------------------------------------------- 364 484 !! *** ROUTINE interpvn_sponge *** 365 485 !!--------------------------------------------- 366 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 367 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2 ), INTENT(inout) :: tabres486 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 487 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres 368 488 LOGICAL, INTENT(in) :: before 369 489 INTEGER, INTENT(in) :: nb , ndir 370 490 ! 371 INTEGER :: ji, jj, jk 372 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr 373 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: vbdiff 374 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff 375 INTEGER :: imax 491 INTEGER :: ji, jj, jk, imax 492 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff 493 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff 494 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff 495 ! vertical interpolation: 496 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 497 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 498 REAL(wp), DIMENSION(1:jpk) :: h_out 499 INTEGER :: N_in, N_out 376 500 !!--------------------------------------------- 377 501 378 502 IF( before ) THEN 379 tabres = vb(i1:i2,j1:j2,:) 503 DO jk=1,jpkm1 504 DO jj=j1,j2 505 DO ji=i1,i2 506 tabres(ji,jj,jk,m1) = vb(ji,jj,jk) 507 # if defined key_vertical 508 tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_n(ji,jj,jk) 509 # endif 510 END DO 511 END DO 512 END DO 380 513 ELSE 381 ! 382 vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(:,:,:))*vmask(i1:i2,j1:j2,:) 514 515 # if defined key_vertical 516 tabres_child(:,:,:) = 0._wp 517 DO jj=j1,j2 518 DO ji=i1,i2 519 N_in = 0 520 DO jk=k1,k2 521 IF (tabres(ji,jj,jk,m2) == 0) EXIT 522 N_in = N_in + 1 523 tabin(jk) = tabres(ji,jj,jk,m1) 524 h_in(N_in) = tabres(ji,jj,jk,m2) 525 ENDDO 526 527 IF (N_in == 0) THEN 528 tabres_child(ji,jj,:) = 0. 529 CYCLE 530 ENDIF 531 532 N_out = 0 533 DO jk=1,jpk 534 if (vmask(ji,jj,jk) == 0) EXIT 535 N_out = N_out + 1 536 h_out(N_out) = e3v_n(ji,jj,jk) 537 ENDDO 538 539 IF (N_in * N_out > 0) THEN 540 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 541 if (h_diff < -1.e4) then 542 print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 543 endif 544 ENDIF 545 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) 546 ENDDO 547 ENDDO 548 549 vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:) 550 # else 551 vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:) 552 # endif 383 553 ! 384 554 DO jk = 1, jpkm1 ! Horizontal slab -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90
r8993 r9007 1 1 #define TWO_WAY /* TWO WAY NESTING */ 2 # undef DECAL_FEEDBACK/* SEPARATION of INTERFACES*/2 #define DECAL_FEEDBACK /* SEPARATION of INTERFACES*/ 3 3 #undef VOL_REFLUX /* VOLUME REFLUXING*/ 4 4 … … 97 97 ! Account for updated thicknesses at boundary edges 98 98 IF (.NOT.ln_linssh) THEN 99 ! For the time being calls below do not ensure reproducible results100 99 ! CALL Agrif_Update_Variable(un_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_u_bdy) 101 100 ! CALL Agrif_Update_Variable(vn_update_id,locupdate1=(/0,0/),locupdate2=(/0,0/),procname = correct_v_bdy) … … 283 282 END SUBROUTINE dom_vvl_update_UVF 284 283 284 #if defined key_vertical 285 286 SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 287 !!--------------------------------------------- 288 !! *** ROUTINE updateT *** 289 !!--------------------------------------------- 290 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 291 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 292 LOGICAL, INTENT(in) :: before 293 !! 294 INTEGER :: ji,jj,jk,jn 295 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child 296 REAL(wp) :: h_in(k1:k2) 297 REAL(wp) :: h_out(1:jpk) 298 INTEGER :: N_in, N_out 299 REAL(wp) :: h_diff 300 REAL(wp) :: zrho_xy 301 REAL(wp) :: tabin(k1:k2,n1:n2) 302 !!--------------------------------------------- 303 ! 304 IF (before) THEN 305 AGRIF_SpecialValue = -999._wp 306 zrho_xy = Agrif_rhox() * Agrif_rhoy() 307 DO jn = n1,n2-1 308 DO jk=k1,k2 309 DO jj=j1,j2 310 DO ji=i1,i2 311 tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) & 312 * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp 313 END DO 314 END DO 315 END DO 316 END DO 317 DO jk=k1,k2 318 DO jj=j1,j2 319 DO ji=i1,i2 320 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) & 321 + (tmask(ji,jj,jk)-1)*999._wp 322 END DO 323 END DO 324 END DO 325 ELSE 326 tabres_child(:,:,:,:) = 0. 327 AGRIF_SpecialValue = 0._wp 328 DO jj=j1,j2 329 DO ji=i1,i2 330 N_in = 0 331 DO jk=k1,k2 !k2 = jpk of child grid 332 IF (tabres(ji,jj,jk,n2) == 0 ) EXIT 333 N_in = N_in + 1 334 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 335 h_in(N_in) = tabres(ji,jj,jk,n2) 336 ENDDO 337 N_out = 0 338 DO jk=1,jpk ! jpk of parent grid 339 IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF 340 N_out = N_out + 1 341 h_out(N_out) = e3t_n(ji,jj,jk) 342 ENDDO 343 IF (N_in > 0) THEN !Remove this? 344 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 345 IF (h_diff < -1.e-4) THEN 346 print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out)) 347 print *,h_in(1:N_in) 348 print *,h_out(1:N_out) 349 STOP 350 ENDIF 351 DO jn=n1,n2-1 352 CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 353 ENDDO 354 ENDIF 355 ENDDO 356 ENDDO 357 358 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 359 ! Add asselin part 360 DO jn = n1,n2-1 361 DO jk=1,jpk 362 DO jj=j1,j2 363 DO ji=i1,i2 364 IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 365 tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 366 & + atfp * ( tabres_child(ji,jj,jk,jn) & 367 & - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 368 ENDIF 369 ENDDO 370 ENDDO 371 ENDDO 372 ENDDO 373 ENDIF 374 DO jn = n1,n2-1 375 DO jk=1,jpk 376 DO jj=j1,j2 377 DO ji=i1,i2 378 IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 379 tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk) 380 END IF 381 END DO 382 END DO 383 END DO 384 END DO 385 ENDIF 386 ! 387 END SUBROUTINE updateTS 388 389 # else 390 285 391 SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 286 392 !!--------------------------------------------- … … 296 402 ! 297 403 IF (before) THEN 298 DO jn = n1,n2404 DO jn = 1,jpts 299 405 DO jk=k1,k2 300 406 DO jj=j1,j2 … … 310 416 ELSE 311 417 !> jc tmp 312 DO jn = n1,n2418 DO jn = 1,jpts 313 419 tabres(i1:i2,j1:j2,k1:k2,jn) = tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) & 314 420 & * tmask(i1:i2,j1:j2,k1:k2) … … 317 423 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 318 424 ! Add asselin part 319 DO jn = n1,n2425 DO jn = 1,jpts 320 426 DO jk=k1,k2 321 427 DO jj=j1,j2 … … 333 439 ENDDO 334 440 ENDIF 335 DO jn = n1,n2441 DO jn = 1,jpts 336 442 DO jk=k1,k2 337 443 DO jj=j1,j2 … … 346 452 ! 347 453 IF ((neuler==0).AND.(Agrif_Nb_Step()==0) ) THEN 348 tsb(i1:i2,j1:j2,k1:k2, n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2)454 tsb(i1:i2,j1:j2,k1:k2,1:jpts) = tsn(i1:i2,j1:j2,k1:k2,1:jpts) 349 455 ENDIF 350 456 ! … … 353 459 END SUBROUTINE updateTS 354 460 355 356 SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, before ) 461 # endif 462 463 # if defined key_vertical 464 465 SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 357 466 !!--------------------------------------------- 358 467 !! *** ROUTINE updateu *** 359 468 !!--------------------------------------------- 360 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2 361 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres 362 LOGICAL , INTENT(in ) :: before 469 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 470 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 471 LOGICAL , INTENT(in ) :: before 472 ! 473 INTEGER :: ji, jj, jk 474 REAL(wp) :: zrhoy 475 ! VERTICAL REFINEMENT BEGIN 476 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 477 REAL(wp) :: h_in(k1:k2) 478 REAL(wp) :: h_out(1:jpk) 479 INTEGER :: N_in, N_out 480 REAL(wp) :: h_diff, excess, thick 481 REAL(wp) :: tabin(k1:k2) 482 ! VERTICAL REFINEMENT END 483 !!--------------------------------------------- 484 ! 485 IF( before ) THEN 486 zrhoy = Agrif_Rhoy() 487 AGRIF_SpecialValue = -999._wp 488 DO jk=k1,k2 489 DO jj=j1,j2 490 DO ji=i1,i2 491 tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk) & 492 + (umask(ji,jj,jk)-1)*999._wp 493 tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk) & 494 + (umask(ji,jj,jk)-1)*999._wp 495 END DO 496 END DO 497 END DO 498 ELSE 499 tabres_child(:,:,:) = 0. 500 AGRIF_SpecialValue = 0._wp 501 DO jj=j1,j2 502 DO ji=i1,i2 503 N_in = 0 504 h_in(:) = 0._wp 505 tabin(:) = 0._wp 506 DO jk=k1,k2 !k2=jpk of child grid 507 IF( tabres(ji,jj,jk,2) < -900) EXIT 508 N_in = N_in + 1 509 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 510 h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj) 511 ENDDO 512 N_out = 0 513 DO jk=1,jpk 514 IF (umask(ji,jj,jk) == 0) EXIT 515 N_out = N_out + 1 516 h_out(N_out) = e3u_n(ji,jj,jk) 517 ENDDO 518 IF (N_in * N_out > 0) THEN 519 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 520 IF (h_diff < -1.e-4) THEN 521 !Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 522 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 523 excess = 0._wp 524 DO jk=N_in,1,-1 525 thick = MIN(-1*h_diff, h_in(jk)) 526 excess = excess + tabin(jk)*thick*e2u(ji,jj) 527 tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 528 h_diff = h_diff + thick 529 IF ( h_diff == 0) THEN 530 N_in = jk 531 h_in(jk) = h_in(jk) - thick 532 EXIT 533 ENDIF 534 ENDDO 535 ENDIF 536 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) 537 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) 538 ENDIF 539 ENDDO 540 ENDDO 541 542 DO jk=1,jpk 543 DO jj=j1,j2 544 DO ji=i1,i2 545 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 546 ub(ji,jj,jk) = ub(ji,jj,jk) & 547 & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 548 ENDIF 549 ! 550 un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk) 551 END DO 552 END DO 553 END DO 554 ENDIF 555 ! 556 END SUBROUTINE updateu 557 558 #else 559 560 SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 561 !!--------------------------------------------- 562 !! *** ROUTINE updateu *** 563 !!--------------------------------------------- 564 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 565 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 566 LOGICAL , INTENT(in ) :: before 363 567 ! 364 568 INTEGER :: ji, jj, jk … … 369 573 zrhoy = Agrif_Rhoy() 370 574 DO jk = k1, k2 371 tabres(i1:i2,j1:j2,jk ) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk)575 tabres(i1:i2,j1:j2,jk,1) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk) 372 576 END DO 373 577 ELSE … … 375 579 DO jj=j1,j2 376 580 DO ji=i1,i2 377 tabres(ji,jj,jk ) = tabres(ji,jj,jk) * r1_e2u(ji,jj)581 tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e2u(ji,jj) 378 582 ! 379 583 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 380 584 zub = ub(ji,jj,jk) * e3u_b(ji,jj,jk) ! fse3t_b prior update should be used 381 585 zuno = un(ji,jj,jk) * e3u_a(ji,jj,jk) 382 zunu = tabres(ji,jj,jk )586 zunu = tabres(ji,jj,jk,1) 383 587 ub(ji,jj,jk) = ( zub + atfp * ( zunu - zuno) ) & 384 588 & * umask(ji,jj,jk) / e3u_b(ji,jj,jk) 385 589 ENDIF 386 590 ! 387 un(ji,jj,jk) = tabres(ji,jj,jk ) * umask(ji,jj,jk) / e3u_n(ji,jj,jk)591 un(ji,jj,jk) = tabres(ji,jj,jk,1) * umask(ji,jj,jk) / e3u_n(ji,jj,jk) 388 592 END DO 389 593 END DO … … 398 602 END SUBROUTINE updateu 399 603 400 SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, before, nb, ndir ) 604 # endif 605 606 SUBROUTINE correct_u_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 401 607 !!--------------------------------------------- 402 608 !! *** ROUTINE correct_u_bdy *** 403 609 !!--------------------------------------------- 404 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2405 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2 ), INTENT(inout) :: tabres406 LOGICAL , INTENT(in ) :: before407 INTEGER , INTENT(in) :: nb, ndir610 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 611 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 612 LOGICAL , INTENT(in ) :: before 613 INTEGER , INTENT(in) :: nb, ndir 408 614 !! 409 615 LOGICAL :: western_side, eastern_side … … 442 648 END SUBROUTINE correct_u_bdy 443 649 444 445 SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, before) 650 # if defined key_vertical 651 652 SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 446 653 !!--------------------------------------------- 447 654 !! *** ROUTINE updatev *** 448 655 !!--------------------------------------------- 449 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2 450 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres 451 LOGICAL , INTENT(in ) :: before 656 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 657 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 658 LOGICAL , INTENT(in ) :: before 659 ! 660 INTEGER :: ji, jj, jk 661 REAL(wp) :: zrhox 662 ! VERTICAL REFINEMENT BEGIN 663 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 664 REAL(wp) :: h_in(k1:k2) 665 REAL(wp) :: h_out(1:jpk) 666 INTEGER :: N_in, N_out 667 REAL(wp) :: h_diff, excess, thick 668 REAL(wp) :: tabin(k1:k2) 669 ! VERTICAL REFINEMENT END 670 !!--------------------------------------------- 671 ! 672 IF (before) THEN 673 zrhox = Agrif_Rhox() 674 AGRIF_SpecialValue = -999._wp 675 DO jk=k1,k2 676 DO jj=j1,j2 677 DO ji=i1,i2 678 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) & 679 + (vmask(ji,jj,jk)-1)*999._wp 680 tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) & 681 + (vmask(ji,jj,jk)-1)*999._wp 682 END DO 683 END DO 684 END DO 685 ELSE 686 tabres_child(:,:,:) = 0. 687 AGRIF_SpecialValue = 0._wp 688 DO jj=j1,j2 689 DO ji=i1,i2 690 N_in = 0 691 DO jk=k1,k2 692 IF (tabres(ji,jj,jk,2) < -900) EXIT 693 N_in = N_in + 1 694 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 695 h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj) 696 ENDDO 697 N_out = 0 698 DO jk=1,jpk 699 IF (vmask(ji,jj,jk) == 0) EXIT 700 N_out = N_out + 1 701 h_out(N_out) = e3v_n(ji,jj,jk) 702 ENDDO 703 IF (N_in * N_out > 0) THEN 704 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 705 IF (h_diff < -1.e-4) then 706 !Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 707 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 708 excess = 0._wp 709 DO jk=N_in,1,-1 710 thick = MIN(-1*h_diff, h_in(jk)) 711 excess = excess + tabin(jk)*thick*e2u(ji,jj) 712 tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 713 h_diff = h_diff + thick 714 IF ( h_diff == 0) THEN 715 N_in = jk 716 h_in(jk) = h_in(jk) - thick 717 EXIT 718 ENDIF 719 ENDDO 720 ENDIF 721 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) 722 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out)) 723 ENDIF 724 ENDDO 725 ENDDO 726 727 DO jk=1,jpk 728 DO jj=j1,j2 729 DO ji=i1,i2 730 ! 731 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 732 vb(ji,jj,jk) = vb(ji,jj,jk) & 733 & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 734 ENDIF 735 ! 736 vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk) 737 END DO 738 END DO 739 END DO 740 ENDIF 741 ! 742 END SUBROUTINE updatev 743 744 # else 745 746 SUBROUTINE updatev( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before) 747 !!--------------------------------------------- 748 !! *** ROUTINE updatev *** 749 !!--------------------------------------------- 750 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 751 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 752 LOGICAL , INTENT(in ) :: before 452 753 ! 453 754 INTEGER :: ji, jj, jk … … 460 761 DO jj=j1,j2 461 762 DO ji=i1,i2 462 tabres(ji,jj,jk ) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)763 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) 463 764 END DO 464 765 END DO … … 468 769 DO jj=j1,j2 469 770 DO ji=i1,i2 470 tabres(ji,jj,jk ) = tabres(ji,jj,jk) * r1_e1v(ji,jj)771 tabres(ji,jj,jk,1) = tabres(ji,jj,jk,1) * r1_e1v(ji,jj) 471 772 ! 472 773 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 473 774 zvb = vb(ji,jj,jk) * e3v_b(ji,jj,jk) ! fse3t_b prior update should be used 474 775 zvno = vn(ji,jj,jk) * e3v_a(ji,jj,jk) 475 zvnu = tabres(ji,jj,jk )776 zvnu = tabres(ji,jj,jk,1) 476 777 vb(ji,jj,jk) = ( zvb + atfp * ( zvnu - zvno) ) & 477 778 & * vmask(ji,jj,jk) / e3v_b(ji,jj,jk) 478 779 ENDIF 479 780 ! 480 vn(ji,jj,jk) = tabres(ji,jj,jk ) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk)781 vn(ji,jj,jk) = tabres(ji,jj,jk,1) * vmask(ji,jj,jk) / e3v_n(ji,jj,jk) 481 782 END DO 482 783 END DO … … 491 792 END SUBROUTINE updatev 492 793 493 SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, before, nb, ndir ) 794 # endif 795 796 SUBROUTINE correct_v_bdy( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir ) 494 797 !!--------------------------------------------- 495 798 !! *** ROUTINE correct_u_bdy *** 496 799 !!--------------------------------------------- 497 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2498 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2 ), INTENT(inout) :: tabres499 LOGICAL , INTENT(in ) :: before500 INTEGER , INTENT(in) :: nb, ndir800 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 801 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 802 LOGICAL , INTENT(in ) :: before 803 INTEGER , INTENT(in) :: nb, ndir 501 804 !! 502 805 LOGICAL :: southern_side, northern_side -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/NST_SRC/agrif_top_interp.F90
r6140 r9007 50 50 ! 51 51 INTEGER :: ji, jj, jk, jn ! dummy loop indices 52 INTEGER :: imin, imax, jmin, jmax52 INTEGER :: imin, imax, jmin, jmax, N_in, N_out 53 53 REAL(wp) :: zrhox , zalpha1, zalpha2, zalpha3 54 54 REAL(wp) :: zalpha4, zalpha5, zalpha6, zalpha7 55 55 LOGICAL :: western_side, eastern_side,northern_side,southern_side 56 56 ! vertical interpolation: 57 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child 58 REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 59 REAL(wp), DIMENSION(k1:k2) :: h_in 60 REAL(wp), DIMENSION(1:jpk) :: h_out(1:jpk) 61 REAL(wp) :: h_diff, zrhoxy 62 63 zrhoxy = Agrif_rhox()*Agrif_rhoy() 57 64 IF (before) THEN 58 ptab(i1:i2,j1:j2,k1:k2,n1:n2) = trn(i1:i2,j1:j2,k1:k2,n1:n2) 65 DO jn = 1,jpts 66 DO jk=k1,k2 67 DO jj=j1,j2 68 DO ji=i1,i2 69 ptab(ji,jj,jk,jn) = trn(ji,jj,jk,jn) 70 END DO 71 END DO 72 END DO 73 END DO 74 # if defined key_vertical 75 DO jk=k1,k2 76 DO jj=j1,j2 77 DO ji=i1,i2 78 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 79 END DO 80 END DO 81 END DO 82 # endif 83 59 84 ELSE 60 85 ! … … 63 88 southern_side = (nb == 2).AND.(ndir == 1) 64 89 northern_side = (nb == 2).AND.(ndir == 2) 90 91 # if defined key_vertical 92 DO jj=j1,j2 93 DO ji=i1,i2 94 iref = ji 95 jref = jj 96 if(western_side) iref=MAX(2,ji) 97 if(eastern_side) iref=MIN(nlci-1,ji) 98 if(southern_side) jref=MAX(2,jj) 99 if(northern_side) jref=MIN(nlcj-1,jj) 100 N_in = 0 101 DO jk=k1,k2 !k2 = jpk of parent grid 102 IF (ptab(ji,jj,jk,n2) == 0) EXIT 103 N_in = N_in + 1 104 tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1) 105 h_in(N_in) = ptab(ji,jj,jk,n2) 106 END DO 107 N_out = 0 108 DO jk=1,jpk ! jpk of child grid 109 IF (tmask(iref,jref,jk) == 0) EXIT 110 N_out = N_out + 1 111 h_out(jk) = e3t_n(iref,jref,jk) 112 ENDDO 113 IF (N_in > 0) THEN 114 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 115 DO jn=1,jptra 116 call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 117 ENDDO 118 ENDIF 119 ENDDO 120 ENDDO 121 # else 122 ptab_child(i1:i2,j1:j2,1:jpk,1:jptra) = ptab(i1:i2,j1:j2,1:jpk,1:jptra) 123 # endif 124 65 125 ! 66 126 zrhox = Agrif_Rhox() … … 89 149 IF( eastern_side) THEN 90 150 DO jn = 1, jptra 91 tra(nlci,j1:j2, k1:k2,jn) = zalpha1 * ptab(nlci,j1:j2,k1:k2,jn) + zalpha2 * ptab(nlci-1,j1:j2,k1:k2,jn)151 tra(nlci,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(nlci,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(nlci-1,j1:j2,1:jpk,jn) 92 152 DO jk = 1, jpkm1 93 153 DO jj = jmin,jmax … … 108 168 IF( northern_side ) THEN 109 169 DO jn = 1, jptra 110 tra(i1:i2,nlcj, k1:k2,jn) = zalpha1 * ptab(i1:i2,nlcj,k1:k2,jn) + zalpha2 * ptab(i1:i2,nlcj-1,k1:k2,jn)170 tra(i1:i2,nlcj,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,nlcj,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,nlcj-1,1:jpk,jn) 111 171 DO jk = 1, jpkm1 112 172 DO ji = imin,imax … … 127 187 IF( western_side) THEN 128 188 DO jn = 1, jptra 129 tra(1,j1:j2, k1:k2,jn) = zalpha1 * ptab(1,j1:j2,k1:k2,jn) + zalpha2 * ptab(2,j1:j2,k1:k2,jn)189 tra(1,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(1,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(2,j1:j2,1:jpk,jn) 130 190 DO jk = 1, jpkm1 131 191 DO jj = jmin,jmax … … 145 205 IF( southern_side ) THEN 146 206 DO jn = 1, jptra 147 tra(i1:i2,1, k1:k2,jn) = zalpha1 * ptab(i1:i2,1,k1:k2,jn) + zalpha2 * ptab(i1:i2,2,k1:k2,jn)148 DO jk=1,jpk 207 tra(i1:i2,1,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,1,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,2,1:jpk,jn) 208 DO jk=1,jpkm1 149 209 DO ji=imin,imax 150 210 IF( vmask(ji,2,jk) == 0.e0 ) THEN … … 165 225 ! East south 166 226 IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN 167 tra(nlci-1,2,:,:) = ptab (nlci-1,2,:,:)227 tra(nlci-1,2,:,:) = ptab_child(nlci-1,2,:,:) 168 228 ENDIF 169 229 ! East north 170 230 IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN 171 tra(nlci-1,nlcj-1,:,:) = ptab (nlci-1,nlcj-1,:,:)231 tra(nlci-1,nlcj-1,:,:) = ptab_child(nlci-1,nlcj-1,:,:) 172 232 ENDIF 173 233 ! West south 174 234 IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN 175 tra(2,2,:,:) = ptab (2,2,:,:)235 tra(2,2,:,:) = ptab_child(2,2,:,:) 176 236 ENDIF 177 237 ! West north 178 238 IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN 179 tra(2,nlcj-1,:,:) = ptab (2,nlcj-1,:,:)239 tra(2,nlcj-1,:,:) = ptab_child(2,nlcj-1,:,:) 180 240 ENDIF 181 241 ! -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/NST_SRC/agrif_top_sponge.F90
r8993 r9007 72 72 REAL(wp), DIMENSION(i1:i2,j1:j2) :: ztu, ztv 73 73 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2) :: trbdiff 74 ! vertical interpolation: 75 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child 76 REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 77 REAL(wp), DIMENSION(k1:k2) :: h_in 78 REAL(wp), DIMENSION(1:jpk) :: h_out 79 INTEGER :: N_in, N_out 80 REAL(wp) :: h_diff 74 81 !!---------------------------------------------------------------------- 75 82 ! 76 83 IF( before ) THEN 77 tabres(i1:i2,j1:j2,k1:k2,n1:n2) = trb(i1:i2,j1:j2,k1:k2,n1:n2) 84 DO jn = 1, jptra 85 DO jk=k1,k2 86 DO jj=j1,j2 87 DO ji=i1,i2 88 tabres(ji,jj,jk,jn) = trb(ji,jj,jk,jn) 89 END DO 90 END DO 91 END DO 92 END DO 93 94 # if defined key_vertical 95 DO jk=k1,k2 96 DO jj=j1,j2 97 DO ji=i1,i2 98 tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 99 END DO 100 END DO 101 END DO 102 # endif 78 103 ELSE 79 !!gm line below use of :,: versus i1:i2,j1:j2 .... strange, not wrong. ===>> to be corrected 80 trbdiff(:,:,:,:) = trb(i1:i2,j1:j2,:,:) - tabres(:,:,:,:) 104 # if defined key_vertical 105 tabres_child(:,:,:,:) = 0. 106 DO jj=j1,j2 107 DO ji=i1,i2 108 N_in = 0 109 DO jk=k1,k2 !k2 = jpk of parent grid 110 IF (tabres(ji,jj,jk,n2) == 0) EXIT 111 N_in = N_in + 1 112 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1) 113 h_in(N_in) = tabres(ji,jj,jk,n2) 114 END DO 115 N_out = 0 116 DO jk=1,jpk ! jpk of child grid 117 IF (tmask(ji,jj,jk) == 0) EXIT 118 N_out = N_out + 1 119 h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above 120 ENDDO 121 IF (N_in > 0) THEN 122 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 123 tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for????? 124 DO jn=1,jptra 125 call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 126 ENDDO 127 ENDIF 128 ENDDO 129 ENDDO 130 # endif 131 132 DO jj=j1,j2 133 DO ji=i1,i2 134 DO jk=1,jpkm1 135 # if defined key_vertical 136 trbdiff(ji,jj,jk,1:jptra) = trb(ji,jj,jk,1:jptra) - tabres_child(ji,jj,jk,1:jptra) 137 # else 138 trbdiff(ji,jj,jk,1:jptra) = trb(ji,jj,jk,1:jptra) - tabres(ji,jj,jk,1:jptra) 139 # endif 140 ENDDO 141 ENDDO 142 ENDDO 143 81 144 DO jn = 1, jptra 82 145 DO jk = 1, jpkm1 -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/NST_SRC/agrif_top_update.F90
r8993 r9007 56 56 END SUBROUTINE Agrif_Update_Trc 57 57 58 58 #ifdef key_vertical 59 SUBROUTINE updateTRC( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 60 !!--------------------------------------------- 61 !! *** ROUTINE updateT *** 62 !!--------------------------------------------- 63 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 64 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 65 LOGICAL, INTENT(in) :: before 66 !! 67 INTEGER :: ji,jj,jk,jn 68 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child 69 REAL(wp) :: h_in(k1:k2) 70 REAL(wp) :: h_out(1:jpk) 71 INTEGER :: N_in, N_out 72 REAL(wp) :: h_diff 73 REAL(wp) :: zrho_xy 74 REAL(wp) :: tabin(k1:k2,n1:n2) 75 !!--------------------------------------------- 76 ! 77 IF (before) THEN 78 AGRIF_SpecialValue = -999._wp 79 zrho_xy = Agrif_rhox() * Agrif_rhoy() 80 DO jn = n1,n2-1 81 DO jk=k1,k2 82 DO jj=j1,j2 83 DO ji=i1,i2 84 tabres(ji,jj,jk,jn) = (trn(ji,jj,jk,jn) * e3t_n(ji,jj,jk) ) & 85 * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp 86 END DO 87 END DO 88 END DO 89 END DO 90 DO jk=k1,k2 91 DO jj=j1,j2 92 DO ji=i1,i2 93 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) & 94 + (tmask(ji,jj,jk)-1)*999._wp 95 END DO 96 END DO 97 END DO 98 ELSE 99 tabres_child(:,:,:,:) = 0. 100 AGRIF_SpecialValue = 0._wp 101 DO jj=j1,j2 102 DO ji=i1,i2 103 N_in = 0 104 DO jk=k1,k2 !k2 = jpk of child grid 105 IF (tabres(ji,jj,jk,n2) == 0 ) EXIT 106 N_in = N_in + 1 107 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 108 h_in(N_in) = tabres(ji,jj,jk,n2) 109 ENDDO 110 N_out = 0 111 DO jk=1,jpk ! jpk of parent grid 112 IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF 113 N_out = N_out + 1 114 h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above 115 ENDDO 116 IF (N_in > 0) THEN !Remove this? 117 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 118 IF (h_diff < -1.e-4) THEN 119 print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out)) 120 print *,h_in(1:N_in) 121 print *,h_out(1:N_out) 122 STOP 123 ENDIF 124 DO jn=1,jptra 125 CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 126 ENDDO 127 ENDIF 128 ENDDO 129 ENDDO 130 131 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 132 ! Add asselin part 133 DO jn = 1,jptra 134 DO jk=1,jpk 135 DO jj=j1,j2 136 DO ji=i1,i2 137 IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 138 trb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 139 & + atfp * ( tabres_child(ji,jj,jk,jn) & 140 & - trn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 141 ENDIF 142 ENDDO 143 ENDDO 144 ENDDO 145 ENDDO 146 ENDIF 147 DO jn = 1,jptra 148 DO jk=1,jpk 149 DO jj=j1,j2 150 DO ji=i1,i2 151 IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 152 trn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk) 153 END IF 154 END DO 155 END DO 156 END DO 157 END DO 158 ENDIF 159 ! 160 END SUBROUTINE updateTRC 161 162 163 #else 59 164 SUBROUTINE updateTRC( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 60 165 !!---------------------------------------------------------------------- … … 127 232 ! 128 233 END SUBROUTINE updateTRC 234 #endif 129 235 130 236 #else -
branches/2017/dev_METO_MERCATOR_2017/NEMOGCM/NEMO/NST_SRC/agrif_user.F90
r8993 r9007 1 # undefUPD_HIGH /* MIX HIGH UPDATE */1 #define UPD_HIGH /* MIX HIGH UPDATE */ 2 2 #if defined key_agrif 3 3 !!---------------------------------------------------------------------- … … 330 330 IF (Agrif_Root()) RETURN 331 331 ! 332 CALL Agrif_Update_ssh()333 332 IF (.NOT.ln_linssh) CALL Agrif_Update_vvl() 334 333 CALL Agrif_Update_tra() 335 334 #if defined key_top 336 335 CALL Agrif_Update_Trc() 337 336 #endif 338 337 CALL Agrif_Update_dyn() 339 338 # if defined key_zdftke 340 339 ! JC remove update because this precludes from perfect restartability 341 !! 340 !! CALL Agrif_Update_tke(0) 342 341 # endif 343 342 … … 364 363 ! 1. Declaration of the type of variable which have to be interpolated 365 364 !--------------------------------------------------------------------- 365 # if defined key_vertical 366 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_id) 367 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_sponge_id) 368 369 CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_interp_id) 370 CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_interp_id) 371 CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_update_id) 372 CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_update_id) 373 CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_sponge_id) 374 CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_sponge_id) 375 # else 366 376 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_id) 367 377 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_sponge_id) 368 378 369 CALL agrif_declare_variable((/1,2,0/),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_interp_id) 370 CALL agrif_declare_variable((/2,1,0/),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_interp_id) 371 CALL agrif_declare_variable((/1,2,0/),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_update_id) 372 CALL agrif_declare_variable((/2,1,0/),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_update_id) 373 CALL agrif_declare_variable((/1,2,0/),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_sponge_id) 374 CALL agrif_declare_variable((/2,1,0/),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_sponge_id) 379 CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_interp_id) 380 CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_interp_id) 381 CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_update_id) 382 CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_update_id) 383 CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),un_sponge_id) 384 CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),vn_sponge_id) 385 # endif 375 386 376 387 CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),e3t_id) … … 392 403 CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/), en_id) 393 404 CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),avt_id) 394 CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),avm_id) 405 # if defined key_vertical 406 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) 407 # else 408 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,1/),avm_id) 409 # endif 395 410 # endif 396 411 … … 777 792 ! 1. Declaration of the type of variable which have to be interpolated 778 793 !--------------------------------------------------------------------- 794 # if defined key_vertical 795 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra+1/),trn_id) 796 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra+1/),trn_sponge_id) 797 # else 779 798 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra/),trn_id) 780 799 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jptra/),trn_sponge_id) 800 # endif 781 801 782 802 ! 2. Type of interpolation
Note: See TracChangeset
for help on using the changeset viewer.
