C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/nesting_sannino/nest_driver/interpolation_p2c.F,v 1.1 2010/11/28 03:26:49 jmc Exp $ C $Name: $ C-- File interpolation_p2c.F: Routines for interpolation C-- Contents C-- o INTERPOLATION_P2C :: Interpolate from Parent to Child C-- o BLINT :: Bilinear interpolation subroutine. C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| SUBROUTINE INTERPOLATION_P2C ( & globalP1,globalP2,globalP3,globalP4,globalP5, & NxP,NyP,NrP, & NxC,NyC,NrC, $ WesternB,EasternB, $ P2C_U,P2C_V,P2C_o,P2C1_V,P2C2_V,P2C1_o,P2C2_o, $ P2C_linU,WO3_linU,P2C_linV,WO3_linV, $ Xv_F,Yv_F,Xv_P,Yv_P, $ T_F,S_F,U_F,V_F,ETA_F, $ DEEP_F,DEEP_P & ) C---------------------------------------------------- IMPLICIT NONE C---------------------------------------------------- INTEGER :: I,J,K,II,JJ INTEGER :: WesternB,EasternB INTEGER :: NrP,NxP,NyP INTEGER :: NrC,NxC,NyC C---------------------------------------------------- REAL*8 :: Fp,Fm,Fo,VEL_MEMO INTEGER :: INDC C---------------------------------------------------- C Define Global Variables to Exchange C---------------------------------------------------- REAL*8 :: globalP1(6,NyP,NrP) REAL*8 :: globalP2(6,NyP,NrP) REAL*8 :: globalP3(6,NyP,NrP) REAL*8 :: globalP4(6,NyP,NrP) REAL*8 :: globalP5(6,NyP,NrP) C---------------------------------------------------- C Define CHILD Model Geometry C---------------------------------------------------- REAL*4 :: Xu_F(NxC,NyC) REAL*4 :: Yu_F(NxC,NyC) REAL*4 :: Xv_F(NxC,NyC) REAL*4 :: Yv_F(NxC,NyC) REAL*4 :: Xo_F(NxC,NyC) REAL*4 :: Yo_F(NxC,NyC) REAL*4 :: Xg_F(NxC,NyC) REAL*4 :: Yg_F(NxC,NyC) REAL*4 :: DEEP_F(NxC,NyC,NrC) C---------------------------------------------------- C Define PARENT Model Geometry C---------------------------------------------------- c REAL*4 :: Xu_P(NxP,NyP) REAL*4 :: Yu_P(NxP,NyP) REAL*4 :: Xv_P(NxP,NyP) REAL*4 :: Yv_P(NxP,NyP) REAL*4 :: Xo_P(NxP,NyP) REAL*4 :: Yo_P(NxP,NyP) REAL*4 :: Xg_P(NxP,NyP) REAL*4 :: Yg_P(NxP,NyP) REAL*4 :: DEEP_P(NxP,NyP,NrP) REAL*4 :: DEPDEP C----------------------------------------------------- REAL*4 :: X1,X2,X3,X4,Y1,Y2,Y3,Y4 REAL*4 :: f1,f2,f3,f4,f,x,y REAL*4 :: gammaT,gammaS,terzo,dueterzi REAL*4 :: gammaEta REAL*4 :: gammaV C---------------------------------------------------- C Define INDICIES MATRIX C---------------------------------------------------- INTEGER :: P2C_U(NyC) !x imposizione NETTA INTEGER :: P2C_linU(NyC) !x Lineare INTEGER :: WO3_linU(NyC) !x Lineare !Which Of 3 INTEGER :: P2C_linV(NyC) !x Lineare INTEGER :: WO3_linV(NyC) !x Lineare !Which Of 3 INTEGER :: P2C_V(NyC) !x Lineare INTEGER :: P2C_o(NyC) !x Lineare INTEGER :: P2C1_V(NyC) !x BiLineare INTEGER :: P2C2_V(NyC) !x BiLineare INTEGER :: P2C1_o(NyC) !x BiLineare INTEGER :: P2C2_o(NyC) !x BiLineare REAL*8 :: diff(NrC) REAL*8 :: DEPDEP_F_WesternB(NrC) REAL*8 :: DEPDEP_F_EasternB(NrC) C---------------------------------------------------- C Define CHILD model variable C---------------------------------------------------- C _____________ (1) WesternB (2) EasternB C | REAL*8 :: U_F(NyC,NrC,2) REAL*8 :: V_F(NyC,NrC,2) REAL*8 :: T_F(NyC,NrC,2) REAL*8 :: S_F(NyC,NrC,2) REAL*8 :: ETA_F(NyC,NrC,2) C---------------------------------------------------- PARAMETER ( terzo = 1./3.) PARAMETER (dueterzi = 2./3.) C======================================================= C (2) INTERPOLATIONS C======================================================= C (2.1) Linear for normal velocity component (u in this case) C------------------------------------------------------- DO K = 1,NrP DO J = 1,NyP-1 IF (globalP4(2,J,K).eq.0.) THEN !uso la salinite come discriminante globalP1 (2,J,K) = globalP1 (2,J+1,K) ENDIF ENDDO DO J = NyP,2,-1 IF (globalP4(2,J,K).eq.0.) THEN !uso la salinite come discriminante globalP1 (2,J,K) = globalP1 (2,J-1,K) ENDIF ENDDO ENDDO C======================================================= C (2.1) NOT Linear but simply imposed C------------------------------------------------------- DO J = 1,3 U_F(J,:,1) = globalP1(2,P2C_U(J),:) U_F(J,:,2) = globalP1(6,P2C_U(J),:) ENDDO DO J = NyC-2,NyC U_F(J,:,1) = globalP1(2,P2C_U(J),:) U_F(J,:,2) = globalP1(6,P2C_U(J),:) ENDDO C================================================= DO J = 4,NyC-3,3 INDC = P2C_U(J) DO K = 1,NrC C-------- WesternB ---------------------- Fp = globalP1(2,INDC+1,K) Fo = globalP1(2,INDC,K) Fm = globalP1(2,INDC-1,K) VEL_MEMO = 0. DO I = -1,1 U_F(J+1+i,K,1) = ((Fp-2.*Fo+Fm)/24.)* & ((12.*float(i)**2+1.)/9.)+ & ((float(i)/6.)*(Fp-Fm))+(26.*Fo-Fp-Fm)/24. VEL_MEMO = VEL_MEMO + U_F(J+1+i,K,1) ENDDO VEL_MEMO = ((3.*Fo) - VEL_MEMO)/3. DO I = -1,1 U_F(J+1+i,K,1) = U_F(J+1+i,K,1) + VEL_MEMO ENDDO C-------- EasternB ---------------------- Fp = globalP1(6,INDC+1,K) Fo = globalP1(6,INDC,K) Fm = globalP1(6,INDC-1,K) VEL_MEMO = 0. DO I = -1,1 U_F(J+1+i,K,2) = ((Fp-2.*Fo+Fm)/24.)* & ((12.*float(i)**2+1.)/9.)+ & ((float(i)/6.)*(Fp-Fm))+(26.*Fo-Fp-Fm)/24. VEL_MEMO = VEL_MEMO + U_F(J+1+i,K,2) ENDDO VEL_MEMO = ((3.*Fo) - VEL_MEMO)/3. DO I = -1,1 U_F(J+1+i,K,2) = U_F(J+1+i,K,2) + VEL_MEMO ENDDO C------------------------------------------------------ ENDDO ENDDO C------------------------------------------------------- C (2.2) BiLinear for tangent velocity component (v in this case) C------------------------------------------------------- I = 1 II = WesternB V_F(:,:,:) = 0. DO K = 1,NrC DO J = 1,NyC x1 = Xv_P(II,P2C1_V(J)) x2 = Xv_P(II,P2C2_V(J)) x3 = Xv_P(II+1,P2C1_V(J)) x4 = Xv_P(II+1,P2C2_V(J)) y1 = Yv_P(II,P2C1_V(J)) y2 = Yv_P(II,P2C2_V(J)) y3 = Yv_P(II+1,P2C1_V(J)) y4 = Yv_P(II+1,P2C2_V(J)) x = Xv_F(I,J) y = Yv_F(I,J) f1 = globalP2(1,P2C1_V(J),K) f2 = globalP2(1,P2C2_V(J),K) f3 = globalP2(2,P2C1_V(J),K) f4 = globalP2(2,P2C2_V(J),K) call blint(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2,f3,f4,x,y,f) V_F(J,K,1) = f ENDDO ENDDO C.............................................. I = NxC II = EasternB DO K = 1,NrC DO J = 1,NyC x1 = Xv_P(II,P2C1_V(J)) x2 = Xv_P(II,P2C2_V(J)) x3 = Xv_P(II-1,P2C1_V(J)) x4 = Xv_P(II-1,P2C2_V(J)) y1 = Yv_P(II,P2C1_V(J)) y2 = Yv_P(II,P2C2_V(J)) y3 = Yv_P(II-1,P2C1_V(J)) y4 = Yv_P(II-1,P2C2_V(J)) x = Xv_F(I,J) y = Yv_F(I,J) f1 = globalP2(6,P2C1_V(J),K) f2 = globalP2(6,P2C2_V(J),K) f3 = globalP2(5,P2C1_V(J),K) f4 = globalP2(5,P2C2_V(J),K) call blint(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2,f3,f4,x,y,f) V_F(J,K,2) = f ENDDO ENDDO C------------------------------------------------------- C (2.3.1) Linear C------------------------------------------------------- C-- WesternB DO K = 1,NrP DO J = 1,NyP DEPDEP = DEEP_P(WesternB,J,K) * DEEP_P(WesternB+1,J,K) gammaT =(globalP3(2,J,K)-globalP3(1,J,K))*DEPDEP gammaS =(globalP4(2,J,K)-globalP4(1,J,K))*DEPDEP gammaeta =(globalP5(2,J,K)-globalP5(1,J,K))*DEPDEP Cgm------- c gammaV =(globalP2(2,J,K)-globalP2(1,J,K))*DEPDEP Cgm--------- globalP3(1,J,K) = globalP3(1,J,K) + terzo* gammaT globalP4(1,J,K) = globalP4(1,J,K) + terzo* gammaS globalP5(1,J,K) = globalP5(1,J,K) + terzo* gammaeta Cgm------- c globalP2(1,J,K) = globalP2(1,J,K) + terzo* gammaV Cgm--------- ENDDO C-------------------------------------------------------------- DO J = 1,NyP-1 IF (globalP4(1,J,K).eq.0.) THEN !uso la salinità come discriminante globalP3(1,J,K) = globalP3(1,J+1,K) globalP4(1,J,K) = globalP4(1,J+1,K) globalP5(1,J,K) = globalP5(1,J+1,K) ENDIF ENDDO DO J = NyP,2,-1 IF (globalP4(1,J,K).eq.0.) THEN !uso la salinità come discriminante globalP3(1,J,K) = globalP3(1,J-1,K) globalP4(1,J,K) = globalP4(1,J-1,K) globalP5(1,J,K) = globalP5(1,J-1,K) Cgm--------- c globalP2(1,J,K) = globalP2(1,J-1,K) Cgm------------- ENDIF ENDDO C--------------------------------------------------------------- ENDDO C-- EasternB DO K = 1,NrP DO J = 1,NyP DEPDEP = DEEP_P(EasternB,J,K) * DEEP_P(EasternB-1,J,K) gammaT =(globalP3(5,J,K)-globalP3(6,J,K))*DEPDEP gammaS =(globalP4(5,J,K)-globalP4(6,J,K))*DEPDEP gammaeta =(globalP5(5,J,K)-globalP5(6,J,K))*DEPDEP Cgm------------- c gammaV =(globalP2(5,J,K)-globalP2(6,J,K))*DEPDEP Cgm---------------- globalP3(6,J,K) = globalP3(6,J,K) + terzo* gammaT globalP4(6,J,K) = globalP4(6,J,K) + terzo* gammaS globalP5(6,J,K) = globalP5(6,J,K) + terzo* gammaeta Cgm---------------- c globalP2(6,J,K) = globalP2(6,J,K) + terzo* gammaV Cgm---------------- ENDDO C-------------------------------------------------------------- DO J = 1,NyP-1 IF (globalP4(6,J,K).eq.0.) THEN !uso la salinità come discriminante globalP3(6,J,K) = globalP3(6,J+1,K) globalP4(6,J,K) = globalP4(6,J+1,K) globalP5(6,J,K) = globalP5(6,J+1,K) Cgm---------------- c globalP2(6,J,K) = globalP2(6,J+1,K) Cgm---------------- ENDIF ENDDO DO J = NyP,2,-1 IF (globalP4(6,J,K).eq.0.) THEN !uso la salinità come discriminante globalP3(6,J,K) = globalP3(6,J-1,K) globalP4(6,J,K) = globalP4(6,J-1,K) globalP5(6,J,K) = globalP5(6,J-1,K) Cgm---------------- c globalP2(6,J,K) = globalP2(6,J-1,K) Cgm---------------- ENDIF ENDDO C--------------------------------------------------------------- ENDDO C------------------------------------------------------- C (2.3.2) Linear C------------------------------------------------------- DO J = 1,NyC C....MASK DEEP_F IF (J.EQ.NyC) THEN !EVITO ERRORI DI CHECK BOUND x Vittorio DEPDEP_F_WesternB (:) = DEEP_F(1 ,J,:)*DEEP_F(1 ,J,:) DEPDEP_F_EasternB(:) = DEEP_F(NxC,J,:)*DEEP_F(NxC,J,:) ELSE DEPDEP_F_WesternB (:) = DEEP_F(1 ,J,:)*DEEP_F(1 ,J+1,:) DEPDEP_F_EasternB(:) = DEEP_F(NxC,J,:)*DEEP_F(NxC,J+1,:) ENDIF C....WesternB...T..................... diff(:) = globalP3(1,P2C_linU(J)+1,:)- & globalP3(1,P2C_linU(J) ,:) T_F(J,:,1) = globalP3(1,P2C_linU(J) ,:)+ & (diff(:)/3.)*float((WO3_linU(J))) & *DEPDEP_F_WesternB(:) C.....EasternB..T..................... diff(:) = globalP3(6,P2C_linU(J)+1,:)- & globalP3(6,P2C_linU(J) ,:) T_F(J,:,2) = globalP3(6,P2C_linU(J) ,:)+ & (diff(:)/3.)*float((WO3_linU(J))) & *DEPDEP_F_EasternB(:) C.....WesternB...S..................... diff(:) = globalP4(1,P2C_linU(J)+1,:)- & globalP4 (1,P2C_linU(J) ,:) S_F(J,:,1) = globalP4(1,P2C_linU(J) ,:)+ & (diff(:)/3.)*float((WO3_linU(J))) & *DEPDEP_F_WesternB(:) C.... EasternB..S..................... diff(:) = (globalP4(6,P2C_linU(J)+1,:)- & globalP4 (6,P2C_linU(J) ,:)) S_F(J,:,2) = globalP4(6,P2C_linU(J) ,:)+ & (diff(:)/3.)*float((WO3_linU(J))) & *DEPDEP_F_EasternB(:) C.... WesternB...Eta..................... diff(:) = globalP5(1,P2C_linU(J)+1,:)- & globalP5(1,P2C_linU(J) ,:) eta_F(J,:,1) = globalP5(1,P2C_linU(J) ,:)+ & (diff(:)/3.)*float((WO3_linU(J))) & *DEPDEP_F_WesternB(:) C.... EasternB..Eta..................... diff(:) = globalP5(6,P2C_linU(J)+1,:)- & globalP5(6,P2C_linU(J) ,:) eta_F(J,:,2) = globalP5(6,P2C_linU(J) ,:)+ & (diff(:)/3.)*float(WO3_linU(J)) & *DEPDEP_F_EasternB(:) Cgm-------------------------- C.... WesternB...V..................... c diff(:) = globalP2(1,P2C_linU(J)+1,:)- c & globalP2(1,P2C_linU(J) ,:) C V_F(J,:,1) = globalP2(1,P2C_linU(J) ,:)+ C & (diff(:)/3.)*float((WO3_linU(J))) C & *DEPDEP_F_WesternB(:) C.... EasternB..V..................... C diff(:) = globalP2(6,P2C_linU(J)+1,:)- C & globalP2(6,P2C_linU(J) ,:) C V_F(J,:,2) = globalP2(6,P2C_linU(J) ,:)+ C & (diff(:)/3.)*float((WO3_linU(J))) C & *DEPDEP_F_EasternB(:) Cgm---------------------------- ENDDO 8765 CONTINUE C=============== END INTERPOLAZIONI x CHILD ==================== RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| C================================================================ SUBROUTINE BLINT(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2,f3,f4,x,y,f) C================================================================ C C Bilinear interpolation subroutine. C (Xi,Yi,fi) = data grid & values surounding model point (x,y) C f = interpolated value at the model grid point. IMPLICIT NONE real*4 a1,a2,a3,a4 real*4 b1,b2,b3,b4 real*4 x1,x2,x3,x4 real*4 y1,y2,y3,y4 real*4 f1,f2,f3,f4 real*4 x,y,A,B,C,t,f,s C a1=x1-x2+x3-x4 a2=-x1+x4 a3=-x1+x2 a4=x1-x b1=y1-y2+y3-y4 b2=-y1+y4 b3=-y1+y2 b4=y1-y A=a3*b1-a1*b3 B=b2*a3+b1*a4-a1*b4-a2*b3 C=-a2*b4+a4*b2 if(ABS(A*C).gt.0.002*B**2) then t=(-B-sqrt(B*B-4.*A*C))/(2.*A) else t=C/ABS(B) endif 10 CONTINUE A=a2*b1-a1*b2 B=b3*a2+b1*a4-a1*b4-a3*b2 C=-a3*b4+a4*b3 if(ABS(A*C).gt.0.002*B**2) then s=(-B+sqrt(B*B-4.*A*C))/(2.*A) else s=-C/ABS(B) endif 20 CONTINUE f=f1*(1.-t)*(1.-s)+f2*t*(1.-s)+f3*s*t+f4*(1.-t)*s return end