/[MITgcm]/MITgcm_contrib/nesting_sannino/nest_driver/interpolation_p2c.F
ViewVC logotype

Annotation of /MITgcm_contrib/nesting_sannino/nest_driver/interpolation_p2c.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (hide annotations) (download)
Sun Nov 28 03:26:49 2010 UTC (14 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: HEAD
plit driver_nesting.F into "main.F" and "interpolation_p2c.F".

1 jmc 1.1 C $Header: $
2     C $Name: $
3    
4     C-- File interpolation_p2c.F: Routines for interpolation
5     C-- Contents
6     C-- o INTERPOLATION_P2C :: Interpolate from Parent to Child
7     C-- o BLINT :: Bilinear interpolation subroutine.
8    
9     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
10    
11     SUBROUTINE INTERPOLATION_P2C (
12     & globalP1,globalP2,globalP3,globalP4,globalP5,
13     & NxP,NyP,NrP,
14     & NxC,NyC,NrC,
15     $ WesternB,EasternB,
16     $ P2C_U,P2C_V,P2C_o,P2C1_V,P2C2_V,P2C1_o,P2C2_o,
17     $ P2C_linU,WO3_linU,P2C_linV,WO3_linV,
18     $ Xv_F,Yv_F,Xv_P,Yv_P,
19     $ T_F,S_F,U_F,V_F,ETA_F,
20     $ DEEP_F,DEEP_P
21     & )
22     C----------------------------------------------------
23     IMPLICIT NONE
24     C----------------------------------------------------
25     INTEGER :: I,J,K,II,JJ
26     INTEGER :: WesternB,EasternB
27     INTEGER :: NrP,NxP,NyP
28     INTEGER :: NrC,NxC,NyC
29     C----------------------------------------------------
30     REAL*8 :: Fp,Fm,Fo,VEL_MEMO
31     INTEGER :: INDC
32     C----------------------------------------------------
33     C Define Global Variables to Exchange
34     C----------------------------------------------------
35     REAL*8 :: globalP1(6,NyP,NrP)
36     REAL*8 :: globalP2(6,NyP,NrP)
37     REAL*8 :: globalP3(6,NyP,NrP)
38     REAL*8 :: globalP4(6,NyP,NrP)
39     REAL*8 :: globalP5(6,NyP,NrP)
40     C----------------------------------------------------
41     C Define CHILD Model Geometry
42     C----------------------------------------------------
43     REAL*4 :: Xu_F(NxC,NyC)
44     REAL*4 :: Yu_F(NxC,NyC)
45     REAL*4 :: Xv_F(NxC,NyC)
46     REAL*4 :: Yv_F(NxC,NyC)
47     REAL*4 :: Xo_F(NxC,NyC)
48     REAL*4 :: Yo_F(NxC,NyC)
49     REAL*4 :: Xg_F(NxC,NyC)
50     REAL*4 :: Yg_F(NxC,NyC)
51     REAL*4 :: DEEP_F(NxC,NyC,NrC)
52     C----------------------------------------------------
53     C Define PARENT Model Geometry
54     C----------------------------------------------------
55     c REAL*4 :: Xu_P(NxP,NyP)
56     REAL*4 :: Yu_P(NxP,NyP)
57     REAL*4 :: Xv_P(NxP,NyP)
58     REAL*4 :: Yv_P(NxP,NyP)
59     REAL*4 :: Xo_P(NxP,NyP)
60     REAL*4 :: Yo_P(NxP,NyP)
61     REAL*4 :: Xg_P(NxP,NyP)
62     REAL*4 :: Yg_P(NxP,NyP)
63     REAL*4 :: DEEP_P(NxP,NyP,NrP)
64     REAL*4 :: DEPDEP
65     C-----------------------------------------------------
66     REAL*4 :: X1,X2,X3,X4,Y1,Y2,Y3,Y4
67     REAL*4 :: f1,f2,f3,f4,f,x,y
68     REAL*4 :: gammaT,gammaS,terzo,dueterzi
69     REAL*4 :: gammaEta
70     REAL*4 :: gammaV
71     C----------------------------------------------------
72     C Define INDICIES MATRIX
73     C----------------------------------------------------
74     INTEGER :: P2C_U(NyC) !x imposizione NETTA
75     INTEGER :: P2C_linU(NyC) !x Lineare
76     INTEGER :: WO3_linU(NyC) !x Lineare !Which Of 3
77    
78     INTEGER :: P2C_linV(NyC) !x Lineare
79     INTEGER :: WO3_linV(NyC) !x Lineare !Which Of 3
80    
81     INTEGER :: P2C_V(NyC) !x Lineare
82     INTEGER :: P2C_o(NyC) !x Lineare
83    
84     INTEGER :: P2C1_V(NyC) !x BiLineare
85     INTEGER :: P2C2_V(NyC) !x BiLineare
86     INTEGER :: P2C1_o(NyC) !x BiLineare
87     INTEGER :: P2C2_o(NyC) !x BiLineare
88    
89     REAL*8 :: diff(NrC)
90     REAL*8 :: DEPDEP_F_WesternB(NrC)
91     REAL*8 :: DEPDEP_F_EasternB(NrC)
92     C----------------------------------------------------
93     C Define CHILD model variable
94     C----------------------------------------------------
95     C _____________ (1) WesternB (2) EasternB
96     C |
97     REAL*8 :: U_F(NyC,NrC,2)
98     REAL*8 :: V_F(NyC,NrC,2)
99     REAL*8 :: T_F(NyC,NrC,2)
100     REAL*8 :: S_F(NyC,NrC,2)
101     REAL*8 :: ETA_F(NyC,NrC,2)
102     C----------------------------------------------------
103     PARAMETER ( terzo = 1./3.)
104     PARAMETER (dueterzi = 2./3.)
105     C=======================================================
106     C (2) INTERPOLATIONS
107     C=======================================================
108     C (2.1) Linear for normal velocity component (u in this case)
109     C-------------------------------------------------------
110     DO K = 1,NrP
111     DO J = 1,NyP-1
112     IF (globalP4(2,J,K).eq.0.) THEN !uso la salinite come discriminante
113     globalP1 (2,J,K) = globalP1 (2,J+1,K)
114     ENDIF
115     ENDDO
116    
117     DO J = NyP,2,-1
118     IF (globalP4(2,J,K).eq.0.) THEN !uso la salinite come discriminante
119     globalP1 (2,J,K) = globalP1 (2,J-1,K)
120     ENDIF
121     ENDDO
122     ENDDO
123     C=======================================================
124     C (2.1) NOT Linear but simply imposed
125     C-------------------------------------------------------
126     DO J = 1,3
127     U_F(J,:,1) = globalP1(2,P2C_U(J),:)
128     U_F(J,:,2) = globalP1(6,P2C_U(J),:)
129     ENDDO
130    
131     DO J = NyC-2,NyC
132     U_F(J,:,1) = globalP1(2,P2C_U(J),:)
133     U_F(J,:,2) = globalP1(6,P2C_U(J),:)
134     ENDDO
135     C=================================================
136     DO J = 4,NyC-3,3
137     INDC = P2C_U(J)
138     DO K = 1,NrC
139     C-------- WesternB ----------------------
140     Fp = globalP1(2,INDC+1,K)
141     Fo = globalP1(2,INDC,K)
142     Fm = globalP1(2,INDC-1,K)
143    
144     VEL_MEMO = 0.
145     DO I = -1,1
146     U_F(J+1+i,K,1) = ((Fp-2.*Fo+Fm)/24.)*
147     & ((12.*float(i)**2+1.)/9.)+
148     & ((float(i)/6.)*(Fp-Fm))+(26.*Fo-Fp-Fm)/24.
149     VEL_MEMO = VEL_MEMO + U_F(J+1+i,K,1)
150     ENDDO
151    
152     VEL_MEMO = ((3.*Fo) - VEL_MEMO)/3.
153    
154     DO I = -1,1
155     U_F(J+1+i,K,1) = U_F(J+1+i,K,1) + VEL_MEMO
156     ENDDO
157    
158     C-------- EasternB ----------------------
159     Fp = globalP1(6,INDC+1,K)
160     Fo = globalP1(6,INDC,K)
161     Fm = globalP1(6,INDC-1,K)
162    
163     VEL_MEMO = 0.
164     DO I = -1,1
165     U_F(J+1+i,K,2) = ((Fp-2.*Fo+Fm)/24.)*
166     & ((12.*float(i)**2+1.)/9.)+
167     & ((float(i)/6.)*(Fp-Fm))+(26.*Fo-Fp-Fm)/24.
168     VEL_MEMO = VEL_MEMO + U_F(J+1+i,K,2)
169     ENDDO
170    
171     VEL_MEMO = ((3.*Fo) - VEL_MEMO)/3.
172    
173     DO I = -1,1
174     U_F(J+1+i,K,2) = U_F(J+1+i,K,2) + VEL_MEMO
175     ENDDO
176     C------------------------------------------------------
177     ENDDO
178     ENDDO
179     C-------------------------------------------------------
180     C (2.2) BiLinear for tangent velocity component (v in this case)
181     C-------------------------------------------------------
182     I = 1
183     II = WesternB
184    
185     V_F(:,:,:) = 0.
186    
187     DO K = 1,NrC
188     DO J = 1,NyC
189     x1 = Xv_P(II,P2C1_V(J))
190     x2 = Xv_P(II,P2C2_V(J))
191    
192     x3 = Xv_P(II+1,P2C1_V(J))
193     x4 = Xv_P(II+1,P2C2_V(J))
194    
195     y1 = Yv_P(II,P2C1_V(J))
196     y2 = Yv_P(II,P2C2_V(J))
197    
198     y3 = Yv_P(II+1,P2C1_V(J))
199     y4 = Yv_P(II+1,P2C2_V(J))
200    
201     x = Xv_F(I,J)
202     y = Yv_F(I,J)
203    
204     f1 = globalP2(1,P2C1_V(J),K)
205     f2 = globalP2(1,P2C2_V(J),K)
206     f3 = globalP2(2,P2C1_V(J),K)
207     f4 = globalP2(2,P2C2_V(J),K)
208    
209     call blint(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2,f3,f4,x,y,f)
210     V_F(J,K,1) = f
211     ENDDO
212     ENDDO
213     C..............................................
214     I = NxC
215     II = EasternB
216    
217     DO K = 1,NrC
218     DO J = 1,NyC
219     x1 = Xv_P(II,P2C1_V(J))
220     x2 = Xv_P(II,P2C2_V(J))
221    
222     x3 = Xv_P(II-1,P2C1_V(J))
223     x4 = Xv_P(II-1,P2C2_V(J))
224    
225     y1 = Yv_P(II,P2C1_V(J))
226     y2 = Yv_P(II,P2C2_V(J))
227    
228     y3 = Yv_P(II-1,P2C1_V(J))
229     y4 = Yv_P(II-1,P2C2_V(J))
230    
231     x = Xv_F(I,J)
232     y = Yv_F(I,J)
233    
234     f1 = globalP2(6,P2C1_V(J),K)
235     f2 = globalP2(6,P2C2_V(J),K)
236     f3 = globalP2(5,P2C1_V(J),K)
237     f4 = globalP2(5,P2C2_V(J),K)
238    
239     call blint(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2,f3,f4,x,y,f)
240     V_F(J,K,2) = f
241     ENDDO
242     ENDDO
243     C-------------------------------------------------------
244     C (2.3.1) Linear
245     C-------------------------------------------------------
246    
247     C-- WesternB
248    
249     DO K = 1,NrP
250     DO J = 1,NyP
251    
252     DEPDEP = DEEP_P(WesternB,J,K) * DEEP_P(WesternB+1,J,K)
253    
254     gammaT =(globalP3(2,J,K)-globalP3(1,J,K))*DEPDEP
255     gammaS =(globalP4(2,J,K)-globalP4(1,J,K))*DEPDEP
256     gammaeta =(globalP5(2,J,K)-globalP5(1,J,K))*DEPDEP
257     Cgm-------
258     c gammaV =(globalP2(2,J,K)-globalP2(1,J,K))*DEPDEP
259     Cgm---------
260     globalP3(1,J,K) = globalP3(1,J,K) + terzo* gammaT
261     globalP4(1,J,K) = globalP4(1,J,K) + terzo* gammaS
262     globalP5(1,J,K) = globalP5(1,J,K) + terzo* gammaeta
263     Cgm-------
264     c globalP2(1,J,K) = globalP2(1,J,K) + terzo* gammaV
265     Cgm---------
266     ENDDO
267     C--------------------------------------------------------------
268     DO J = 1,NyP-1
269     IF (globalP4(1,J,K).eq.0.) THEN !uso la salinità come discriminante
270     globalP3(1,J,K) = globalP3(1,J+1,K)
271     globalP4(1,J,K) = globalP4(1,J+1,K)
272     globalP5(1,J,K) = globalP5(1,J+1,K)
273     ENDIF
274     ENDDO
275    
276     DO J = NyP,2,-1
277     IF (globalP4(1,J,K).eq.0.) THEN !uso la salinità come discriminante
278     globalP3(1,J,K) = globalP3(1,J-1,K)
279     globalP4(1,J,K) = globalP4(1,J-1,K)
280     globalP5(1,J,K) = globalP5(1,J-1,K)
281     Cgm---------
282     c globalP2(1,J,K) = globalP2(1,J-1,K)
283     Cgm-------------
284     ENDIF
285     ENDDO
286     C---------------------------------------------------------------
287     ENDDO
288    
289     C-- EasternB
290    
291     DO K = 1,NrP
292     DO J = 1,NyP
293    
294     DEPDEP = DEEP_P(EasternB,J,K) * DEEP_P(EasternB-1,J,K)
295    
296     gammaT =(globalP3(5,J,K)-globalP3(6,J,K))*DEPDEP
297     gammaS =(globalP4(5,J,K)-globalP4(6,J,K))*DEPDEP
298     gammaeta =(globalP5(5,J,K)-globalP5(6,J,K))*DEPDEP
299     Cgm-------------
300     c gammaV =(globalP2(5,J,K)-globalP2(6,J,K))*DEPDEP
301     Cgm----------------
302     globalP3(6,J,K) = globalP3(6,J,K) + terzo* gammaT
303     globalP4(6,J,K) = globalP4(6,J,K) + terzo* gammaS
304     globalP5(6,J,K) = globalP5(6,J,K) + terzo* gammaeta
305     Cgm----------------
306     c globalP2(6,J,K) = globalP2(6,J,K) + terzo* gammaV
307     Cgm----------------
308     ENDDO
309     C--------------------------------------------------------------
310     DO J = 1,NyP-1
311     IF (globalP4(6,J,K).eq.0.) THEN !uso la salinità come discriminante
312     globalP3(6,J,K) = globalP3(6,J+1,K)
313     globalP4(6,J,K) = globalP4(6,J+1,K)
314     globalP5(6,J,K) = globalP5(6,J+1,K)
315     Cgm----------------
316     c globalP2(6,J,K) = globalP2(6,J+1,K)
317     Cgm----------------
318     ENDIF
319     ENDDO
320    
321     DO J = NyP,2,-1
322     IF (globalP4(6,J,K).eq.0.) THEN !uso la salinità come discriminante
323     globalP3(6,J,K) = globalP3(6,J-1,K)
324     globalP4(6,J,K) = globalP4(6,J-1,K)
325     globalP5(6,J,K) = globalP5(6,J-1,K)
326     Cgm----------------
327     c globalP2(6,J,K) = globalP2(6,J-1,K)
328     Cgm----------------
329     ENDIF
330     ENDDO
331     C---------------------------------------------------------------
332     ENDDO
333     C-------------------------------------------------------
334     C (2.3.2) Linear
335     C-------------------------------------------------------
336     DO J = 1,NyC
337     C....MASK DEEP_F
338     IF (J.EQ.NyC) THEN !EVITO ERRORI DI CHECK BOUND x Vittorio
339     DEPDEP_F_WesternB (:) = DEEP_F(1 ,J,:)*DEEP_F(1 ,J,:)
340     DEPDEP_F_EasternB(:) = DEEP_F(NxC,J,:)*DEEP_F(NxC,J,:)
341     ELSE
342     DEPDEP_F_WesternB (:) = DEEP_F(1 ,J,:)*DEEP_F(1 ,J+1,:)
343     DEPDEP_F_EasternB(:) = DEEP_F(NxC,J,:)*DEEP_F(NxC,J+1,:)
344     ENDIF
345    
346     C....WesternB...T.....................
347     diff(:) = globalP3(1,P2C_linU(J)+1,:)-
348     & globalP3(1,P2C_linU(J) ,:)
349    
350     T_F(J,:,1) = globalP3(1,P2C_linU(J) ,:)+
351     & (diff(:)/3.)*float((WO3_linU(J)))
352     & *DEPDEP_F_WesternB(:)
353     C.....EasternB..T.....................
354     diff(:) = globalP3(6,P2C_linU(J)+1,:)-
355     & globalP3(6,P2C_linU(J) ,:)
356    
357     T_F(J,:,2) = globalP3(6,P2C_linU(J) ,:)+
358     & (diff(:)/3.)*float((WO3_linU(J)))
359     & *DEPDEP_F_EasternB(:)
360     C.....WesternB...S.....................
361     diff(:) = globalP4(1,P2C_linU(J)+1,:)-
362     & globalP4 (1,P2C_linU(J) ,:)
363    
364     S_F(J,:,1) = globalP4(1,P2C_linU(J) ,:)+
365     & (diff(:)/3.)*float((WO3_linU(J)))
366     & *DEPDEP_F_WesternB(:)
367     C.... EasternB..S.....................
368     diff(:) = (globalP4(6,P2C_linU(J)+1,:)-
369     & globalP4 (6,P2C_linU(J) ,:))
370    
371     S_F(J,:,2) = globalP4(6,P2C_linU(J) ,:)+
372     & (diff(:)/3.)*float((WO3_linU(J)))
373     & *DEPDEP_F_EasternB(:)
374     C.... WesternB...Eta.....................
375     diff(:) = globalP5(1,P2C_linU(J)+1,:)-
376     & globalP5(1,P2C_linU(J) ,:)
377    
378     eta_F(J,:,1) = globalP5(1,P2C_linU(J) ,:)+
379     & (diff(:)/3.)*float((WO3_linU(J)))
380     & *DEPDEP_F_WesternB(:)
381     C.... EasternB..Eta.....................
382     diff(:) = globalP5(6,P2C_linU(J)+1,:)-
383     & globalP5(6,P2C_linU(J) ,:)
384    
385     eta_F(J,:,2) = globalP5(6,P2C_linU(J) ,:)+
386     & (diff(:)/3.)*float(WO3_linU(J))
387     & *DEPDEP_F_EasternB(:)
388    
389     Cgm--------------------------
390    
391     C.... WesternB...V.....................
392     c diff(:) = globalP2(1,P2C_linU(J)+1,:)-
393     c & globalP2(1,P2C_linU(J) ,:)
394     C V_F(J,:,1) = globalP2(1,P2C_linU(J) ,:)+
395     C & (diff(:)/3.)*float((WO3_linU(J)))
396     C & *DEPDEP_F_WesternB(:)
397     C.... EasternB..V.....................
398     C diff(:) = globalP2(6,P2C_linU(J)+1,:)-
399     C & globalP2(6,P2C_linU(J) ,:)
400     C V_F(J,:,2) = globalP2(6,P2C_linU(J) ,:)+
401     C & (diff(:)/3.)*float((WO3_linU(J)))
402     C & *DEPDEP_F_EasternB(:)
403     Cgm----------------------------
404     ENDDO
405    
406     8765 CONTINUE
407     C=============== END INTERPOLAZIONI x CHILD ====================
408    
409     RETURN
410     END
411    
412     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
413    
414     C================================================================
415     SUBROUTINE BLINT(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2,f3,f4,x,y,f)
416     C================================================================
417     C
418     C Bilinear interpolation subroutine.
419     C (Xi,Yi,fi) = data grid & values surounding model point (x,y)
420     C f = interpolated value at the model grid point.
421     IMPLICIT NONE
422     real*4 a1,a2,a3,a4
423     real*4 b1,b2,b3,b4
424     real*4 x1,x2,x3,x4
425     real*4 y1,y2,y3,y4
426     real*4 f1,f2,f3,f4
427     real*4 x,y,A,B,C,t,f,s
428     C
429     a1=x1-x2+x3-x4
430     a2=-x1+x4
431     a3=-x1+x2
432     a4=x1-x
433     b1=y1-y2+y3-y4
434     b2=-y1+y4
435     b3=-y1+y2
436     b4=y1-y
437     A=a3*b1-a1*b3
438     B=b2*a3+b1*a4-a1*b4-a2*b3
439     C=-a2*b4+a4*b2
440     if(ABS(A*C).gt.0.002*B**2) then
441     t=(-B-sqrt(B*B-4.*A*C))/(2.*A)
442     else
443     t=C/ABS(B)
444     endif
445     10 CONTINUE
446     A=a2*b1-a1*b2
447     B=b3*a2+b1*a4-a1*b4-a3*b2
448     C=-a3*b4+a4*b3
449     if(ABS(A*C).gt.0.002*B**2) then
450     s=(-B+sqrt(B*B-4.*A*C))/(2.*A)
451     else
452     s=-C/ABS(B)
453     endif
454     20 CONTINUE
455     f=f1*(1.-t)*(1.-s)+f2*t*(1.-s)+f3*s*t+f4*(1.-t)*s
456     return
457     end

  ViewVC Help
Powered by ViewVC 1.1.22