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

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

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


Revision 1.1 - (hide annotations) (download)
Fri Oct 30 20:23:10 2009 UTC (15 years, 9 months ago) by sannino
Branch: MAIN
External routine driving the two-way nesting

1 sannino 1.1 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
2     C
3     PROGRAM NEST_DRIVER
4     C !DESCRIPTION:
5     C Routine that manages the MPI communication between the CHILD
6     C and PARENT models. It performs also the necessary
7     C interpolations from PARENT2CHILD and CHILD2PARENT.
8     C
9     C ver 1.0 by G. Sannino, V. Ruggiero, A. Carillo, P. Heimbach
10     C
11     C First application described in:
12     C Sannino G.,Herrmann, Carillo, Rupolo, Ruggiero, Artale, Heimbach, 2009:
13     C An eddy-permitting model of the Mediterranean Sea with a two-way grid
14     C refinement at Gibraltar.
15     C Ocean Modelling, 30(1), 56-72, doi: 10.1016/j.ocemod.2009.06.002
16    
17     C
18    
19     C !LOCAL INPUT VARIABLES:
20     c ---------------------------------------------------------------------------------
21     c NST_LEV_TOT :: Total nesting levels (1 for only one nesting)
22     c NST_LEV :: Number of the actual nesting
23     c NCPUs_CHLD :: Number of CPUs used for the CHILD model
24     c at NST_LEV nesting level
25     c NCPUs_PRNT :: Number of CPUs used for the PARENT model
26     c at NST_LEV nesting level
27     c nSxC,nSyF :: Domain decomposition used for CHILD
28     c nSxP,nSyP :: Domain decomposition used for PARENT
29     c OLX,OLY :: Domain dec. overlapping (same for both models)
30     c NrC,NyC,NxC :: Dimension PARENT model
31     c NrF,NyF,NxF :: Dimension CHILD model
32     c WesterB :: Western (i) PARENT index where start the refinement
33     c EasterB :: Eastern (i) PARENT index where finish the refinement
34     c dirNEST :: Directory where are stored the geometry data of both models
35     c ---------------------------------------------------------------------------------
36     CEOP
37     implicit none
38     c--------------------------------------------------------
39     c INPUT VARIABLE DEFINITION
40     c--------------------------------------------------------
41     INTEGER :: NST_LEV_TOT, NST_LEV, NCPUs_PRNT
42     INTEGER :: Count_Lev
43     PARAMETER (NST_LEV_TOT = 1) !Number of Total Nesting Levels
44     PARAMETER (NST_LEV = 1) !Which level am I?
45    
46     INTEGER :: NCPUs_CHLD(NST_LEV_TOT)
47     INTEGER :: MSTR_DRV(NST_LEV_TOT)
48     INTEGER :: MSTR_PRNT(NST_LEV_TOT)
49     INTEGER :: MSTR_CHLD(NST_LEV_TOT)
50    
51     PARAMETER (NCPUs_PRNT = 40)
52    
53     DATA NCPUs_CHLD / 20 /
54     c--------------------------------------------------------
55     INTEGER :: nSxC,nSyC !Domain decomposition CHILD
56     INTEGER :: nSxP,nSyP !Domain decomposition PARENT
57     PARAMETER (nSxC = 4 , nSyC = 5)
58     PARAMETER (nSxP = 8 , nSyP = 5)
59     c--------------------------------------------------------
60     INTEGER :: OLY,OLX !Domain decomposition overlapping
61     c !(the same for both models)
62     PARAMETER (OLX = 3, OLY = 3)
63     c--------------------------------------------------------
64     INTEGER :: NrP,NxP,NyP
65     INTEGER :: NrC,NxC,NyC
66     INTEGER :: IM_C,JM_C
67     INTEGER :: IM_P,JM_P
68     INTEGER :: IndI,IndJ
69     INTEGER :: IndI_P(nSxP*nSyP),IndJ_P(nSxP*nSyP)
70     INTEGER :: IndI_C(nSxC*nSyC),IndJ_C(nSxC*nSyC)
71    
72     INTEGER :: WesternB,EasternB
73     c--------------------------------------------------------
74     PARAMETER (NrP=42, NyP=120,NxP = 400) !PARENT model
75     PARAMETER (NrC=42, NyC=105,NxC = 140) !CHILD model
76     c--------------------------------------------------------
77     PARAMETER (WesternB = 43,EasternB=90)
78     c--------------------------------------------------------
79     CHARACTER :: dirNEST*80
80     c--------------------------------------------------------
81     PARAMETER (dirNEST ="/home/sannino/NESTING/")
82     c--------------------------------------------------------
83     INCLUDE 'mpif.h'
84     INTEGER :: ierr,rank,size,npd
85     INTEGER :: irank,isize
86     INTEGER :: color
87     INTEGER :: istatus,NEST_comm
88     INTEGER :: from,whm,status(MPI_STATUS_SIZE),st_count
89     INTEGER :: I,J,K,II,JJ,Irec,III,JJJ,KK,ICONT
90     INTEGER :: iVar,Indx,Jndx
91     INTEGER :: J1,J2,JJ1,JJ2
92     INTEGER :: I_START,I_END,I_STEP
93     REAL*4 :: XF,YF,XP1,YP1,XP2,YP2,YP3
94     REAL*8 :: TRANSPORT_WEST,TRANSPORT_EAST
95     CHARACTER*10 :: c2i(30)
96     c----------------------------------------------------
97     c Define PARENT Model Geometry
98     c----------------------------------------------------
99     c REAL*4 Xu_P(NxP,NyP)
100     REAL*4 :: Yu_P(NxP,NyP)
101     REAL*4 :: Xv_P(NxP,NyP)
102     REAL*4 :: Yv_P(NxP,NyP)
103     REAL*4 :: Xo_P(NxP,NyP)
104     REAL*4 :: Yo_P(NxP,NyP)
105     REAL*4 :: Xg_P(NxP,NyP)
106     REAL*4 :: Yg_P(NxP,NyP)
107     REAL*4 :: hFacW_P(NxP,NyP,NrP)
108     REAL*4 :: hFacS_P(NxP,NyP,NrP)
109     REAL*4 :: RAC_P(NxP,NyP)
110     REAL*4 :: RAW_P(NxP,NyP)
111     REAL*4 :: RAS_P(NxP,NyP)
112     REAL*4 :: hFacC_P(NxP,NyP,NrP)
113     REAL*4 :: DEEP_P(NxP,NyP,NrP)
114     REAL*4 :: INV_VOL_C_P(NxP,NyP,NrP)
115     REAL*4 :: INV_VOL_S_P(NxP,NyP,NrP)
116     REAL*4 :: INV_VOL_W_P(NxP,NyP,NrP)
117     c----------------------------------------------------
118     c Define CHILD Model Geometry
119     c----------------------------------------------------
120     REAL*4 :: Xu_C(NxC,NyC)
121     REAL*4 :: Yu_C(NxC,NyC)
122     REAL*4 :: Xv_C(NxC,NyC)
123     REAL*4 :: Yv_C(NxC,NyC)
124     REAL*4 :: Xo_C(NxC,NyC)
125     REAL*4 :: Yo_C(NxC,NyC)
126     REAL*4 :: Xg_C(NxC,NyC)
127     REAL*4 :: Yg_C(NxC,NyC)
128     REAL*4 :: hFacW_C(NxC,NyC,NrC)
129     REAL*4 :: hFacS_C(NxC,NyC,NrC)
130     REAL*4 :: RAC_C(NxC,NyC)
131     REAL*4 :: RAW_C(NxC,NyC)
132     REAL*4 :: RAS_C(NxC,NyC)
133     REAL*4 :: hFacC_C(NxC,NyC,NrC)
134     REAL*4 :: DEEP_C(NxC,NyC,NrC)
135     c----------------------------------------------------
136     c Define relative (PARENT-->CHILD) indicies
137     c----------------------------------------------------
138     INTEGER :: P2C_U(NyC)
139     c
140     INTEGER :: P2C_linU(NyC) !Linear interp.
141     INTEGER :: WO3_linU(NyC) !Linear interp. !Which Of 3
142     c
143     INTEGER :: P2C_linV(NyC) !Linear interp.
144     INTEGER :: WO3_linV(NyC) !Linear interp. !Which Of 3
145     c
146     INTEGER :: P2C_V(NyC) !Linear interp.
147     INTEGER :: P2C_o(NyC) !Linear interp.
148     INTEGER :: P2C1_V(NyC) !BiLinear interp.
149     INTEGER :: P2C2_V(NyC) !BiLinear interp.
150     INTEGER :: P2C1_o(NyC) !BiLinear interp.
151     INTEGER :: P2C2_o(NyC) !BiLinear interp.
152     c----------------------------------------------------
153     c Define relative (CHILD-->PARENT) indicies
154     c----------------------------------------------------
155     INTEGER I_C2P(9,NxP,NyP)
156     INTEGER J_C2P(9,NxP,NyP)
157     c----------------------------------------------------
158     c Define CHILD model variable
159     c----------------------------------------------------
160     c _____________ (1) WesternB (2) EasternB
161     c |
162     REAL*8 :: U_C1(NyC,NrC,2)
163     REAL*8 :: V_C1(NyC,NrC,2)
164     REAL*8 :: T_C1(NyC,NrC,2)
165     REAL*8 :: S_C1(NyC,NrC,2)
166     REAL*8 :: ETA_C1(NyC,NrC,2)
167     INTEGER :: MSIZE
168     c
169     REAL*8 :: U_C2(NyC,NrC,2)
170     REAL*8 :: V_C2(NyC,NrC,2)
171     REAL*8 :: T_C2(NyC,NrC,2)
172     REAL*8 :: S_C2(NyC,NrC,2)
173     REAL*8 :: ETA_C2(NyC,NrC,2)
174     c
175     REAL*8,allocatable :: VAR_C1(:,:,:,:)
176     c
177     REAL*8 :: DIFF_U(NyC,NrC,2)
178     REAL*8 :: DIFF_V(NyC,NrC,2)
179     REAL*8 :: DIFF_T(NyC,NrC,2)
180     REAL*8 :: DIFF_S(NyC,NrC,2)
181     REAL*8 :: DIFF_ETA(NyC,NrC,2)
182     c----------------------------------------------------
183     c Define PARENT model variable
184     c----------------------------------------------------
185     REAL*8 :: VAR3D_P(NxP,NyP,NrP,15)
186     REAL*8 :: VAR2D_P(NxP,NyP,4)
187    
188     INTEGER :: ONOFF
189     INTEGER :: index_var3D,index_var2D
190     c---------------------------------------------------------------|
191     c (1) U || (2) V || (3) T || (4) S |
192     c---------------------------------------------------------------|
193     c (5) gU || (6) gV || (7) gT || (8) gS |
194     c---------------------------------------------------------------|
195     c (9) gUNm1 || (10) gVNm1 || (11) gTNm1 || (12) gSNm1 |
196     c---------------------------------------------------------------|
197     c (13) totPhiHyd || (14) IVDConvCount || |
198     c---------------------------------------------------------------|
199     c (15) etaN || (16) etaH || (17) phiHydLow |
200     c---------------------------------------------------------------|
201     c (18) etaNm1 || (19) etaHm1|| |
202     c---------------------------------------------------------------|
203     c---------------------------------------------------------------|
204     c Define Global Variables to Exchange |
205     c---------------------------------------------------------------|
206     REAL*8,allocatable :: globalPA (:,:,:,:) !(6,NyP,NrP,5)
207     REAL*8 :: globalP1(6,NyP,NrP)
208     REAL*8 :: globalP2(6,NyP,NrP)
209     REAL*8 :: globalP3(6,NyP,NrP)
210     REAL*8 :: globalP4(6,NyP,NrP)
211     REAL*8 :: globalP5(6,NyP,NrP)
212     REAL*8 :: globalP6(6,NyP,NrP)
213     REAL*8 :: globalP7(6,NyP,NrP)
214     REAL*8 :: globalP8(6,NyP,NrP)
215     REAL*8 :: globalP9(6,NyP,NrP)
216     REAL*8 :: globalP10(6,NyP,NrP)
217     REAL*8 :: globalP11(6,NyP,NrP)
218     REAL*8 :: globalP12(6,NyP,NrP)
219     REAL*8 :: globalP13(6,NyP,NrP)
220     REAL*8 :: globalP14(6,NyP,NrP)
221     c
222     INTEGER :: index
223     c----------------------------------------------------
224     c Define Global Variables to Exchange
225     c----------------------------------------------------
226     REAL*8 :: globalC3D(NxC,NyC,NrC,15)
227     c |___________ 15 fields
228     c
229    
230     REAL*8 globalC2D(NxC,NyC,4)
231     c |___________ 4 fields
232     c
233     c
234     REAL*8,allocatable :: globalC3D_a(:,:,:,:),globalC2D_a(:,:,:)
235     c
236     INTEGER :: indexF,index1F
237     REAL*4 :: AREA_VOL
238     INTEGER :: vstart,vstop,VCONT,VCONTP(0:3)
239     c----------------------------------------------------
240     c MPI starts here
241     c----------------------------------------------------
242     call MPI_Init(ierr)
243     call MPI_Comm_size(MPI_COMM_WORLD,size,ierr)
244     call MPI_Comm_rank(MPI_COMM_WORLD,rank,ierr)
245     npd=size-(NCPUs_PRNT+NCPUs_CHLD(1))
246     if(rank.lt.npd) color=0
247     call MPI_COMM_SPLIT (MPI_COMM_WORLD, color,0,
248     & NEST_COMM,ierr)
249     c
250     call MPI_Comm_size(NEST_COMM,isize,ierr)
251     call MPI_Comm_rank(NEST_COMM,irank,ierr)
252     c--------------------------------------------------------
253     c COMPUTE MASTER VALUES
254     c--------------------------------------------------------
255     MSTR_DRV(1) = 0
256     c
257     MSTR_PRNT(1) = npd
258     MSTR_CHLD(1) = NCPUs_PRNT + npd
259     c
260     DO Count_Lev = 2, NST_LEV_TOT
261     MSTR_DRV(Count_Lev) = MSTR_CHLD(Count_Lev-1) +
262     & NCPUs_CHLD(Count_Lev - 1)
263     c
264     MSTR_CHLD(Count_Lev) = MSTR_DRV(Count_Lev) + 1
265     MSTR_PRNT(Count_Lev) = MSTR_CHLD(Count_Lev-1)
266     ENDDO
267     c
268     vstart = 1+rank*(nSxP/MSTR_PRNT(1))
269     vstop = (1+rank)*(nSxP/MSTR_PRNT(1))
270     VCONT = (nSxP/npd)*(nSyP)*rank-1
271     VCONTP(rank) = VCONT
272     c--------------------------------------------------------
273     c COMPUTE DOMAIN DECOMPOSITION
274     c--------------------------------------------------------
275     if(rank.eq.0) then
276     c
277     IM_C = int(NxC/nSxC)
278     JM_C = int(NyC/nSyC)
279     c
280     IM_P = int(NxP/nSxP)
281     JM_P = int(NyP/nSyP)
282     c
283     indexF = IM_C*JM_C*NrC*15
284     index1F = IM_C*JM_C*4
285     c
286     indexF = (JM_C+OLY+OLY)*NrC*2*5
287     c
288     index_var3D = IM_P*JM_P*NrP
289     index_var2D = IM_P*JM_P
290     c
291     ICONT = 0
292     DO I = 1,nSxP
293     DO J = 1,nSyP
294     ICONT = ICONT + 1
295     IndI_P(ICONT) = IM_P*(I-1)
296     IndJ_P(ICONT) = JM_P*(J-1)
297     END DO
298     END DO
299     c
300     ICONT = 0
301     DO I = 1,nSxC
302     DO J = 1,nSyC
303     ICONT = ICONT + 1
304     IndI_C(ICONT) = IM_C*(I-1)
305     IndJ_C(ICONT) = JM_C*(J-1)
306     END DO
307     END DO
308     c
309     ALLOCATE (globalPA (6,JM_P,NrP,5))
310     index=6*JM_P*NrP*5
311     c
312     ALLOCATE(globalC3D_a(IM_C,JM_C,NrC,15))
313     ALLOCATE(globalC2D_a(IM_C,JM_C,4))
314     c
315     ALLOCATE (VAR_C1((JM_C+OLY+OLY),NrC,2,5))
316     c--------------------------------------------------------
317     c WARNING
318     c--------------------------------------------------------
319     print*,'*************************************'
320     print*,' have you update geometric files?'
321     print*,' in ./CHILD e ./PARENT'
322     print*,'*************************************'
323     c--------------------------------------------------------
324     c PARENT MODEL
325     c--------------------------------------------------------
326     print*,' [1] Read PARENT model geometry'
327     c----------------------------------------------------
328     c XC & YC
329     c----------------------------------------------------
330     MSIZE = NxP*NyP
331     c
332     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
333     & convert='big_endian',
334     & file=trim(dirNEST)//'/PARENT/XC.data',
335     & form='unformatted')
336     c
337     read (1,REC=1) Xo_P(:,:)
338     close(1)
339     c
340     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
341     & convert='big_endian',
342     & file=trim(dirNEST)//'/PARENT/YC.data',
343     & form='unformatted')
344     c
345     read (1,REC=1) Yo_P(:,:)
346     close(1)
347     c----------------------------------------------------
348     c XG & YG
349     c----------------------------------------------------
350     MSIZE = NxP*NyP
351     c
352     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
353     & convert='big_endian',
354     & file=trim(dirNEST)//'/PARENT/XG.data',
355     & form='unformatted')
356    
357     read (1,REC=1) Xg_P(:,:)
358     close(1)
359     c
360     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
361     & convert='big_endian',
362     & file=trim(dirNEST)//'/PARENT/YG.data',
363     & form='unformatted')
364     c
365     read (1,REC=1) Yg_P(:,:)
366     close(1)
367     c----------------------------------------------------
368     c Yu
369     c----------------------------------------------------
370     DO J = 1,NyP
371     DO I = 1,NxP
372     c Xu_P(I,J) = Xg_P(I,J)
373     Yu_P(I,J) = Yo_P(I,J)
374     ENDDO
375     ENDDO
376     c----------------------------------------------------
377     c Xv & Yv
378     c----------------------------------------------------
379     DO J = 1,NyP
380     DO I = 1,NxP
381     Xv_P(I,J) = Xo_P(I,J)
382     Yv_P(I,J) = Yg_P(I,J)
383     ENDDO
384     ENDDO
385     c----------------------------------------------------
386     c hFacC
387     c----------------------------------------------------
388     MSIZE = NxP*NyP*NrP
389     c
390     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
391     & convert='big_endian',
392     & file=trim(dirNEST)//'/PARENT/hFacC.data',
393     & form='unformatted')
394    
395     read (1,REC=1) hFacC_P(:,:,:)
396     close(1)
397     c----------------------------------------------------
398     c hFacW
399     c----------------------------------------------------
400     MSIZE = NxP*NyP*NrP
401     c
402     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
403     & convert='big_endian',
404     & file=trim(dirNEST)//'/PARENT/hFacW.data',
405     & form='unformatted')
406    
407     read (1,REC=1) hFacW_P(:,:,:)
408     close(1)
409     c----------------------------------------------------
410     c hFacS
411     c----------------------------------------------------
412     MSIZE = NxP*NyP*NrP
413     c
414     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
415     & convert='big_endian',
416     & file=trim(dirNEST)//'/PARENT/hFacS.data',
417     & form='unformatted')
418    
419     read (1,REC=1) hFacS_P(:,:,:)
420     close(1)
421     c----------------------------------------------------
422     c RAC
423     c----------------------------------------------------
424     MSIZE = NxP*NyP
425     c
426     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
427     & convert='big_endian',
428     & file=trim(dirNEST)//'/PARENT/RAC.data',
429     & form='unformatted')
430    
431     read (1,REC=1) RAC_P(:,:)
432     close(1)
433     c----------------------------------------------------
434     c RAW
435     c----------------------------------------------------
436     MSIZE = NxP*NyP
437     c
438     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
439     & convert='big_endian',
440     & file=trim(dirNEST)//'/PARENT/RAW.data',
441     & form='unformatted')
442    
443     read (1,REC=1) RAW_P(:,:)
444     close(1)
445     c----------------------------------------------------
446     c RAS
447     c----------------------------------------------------
448     MSIZE = NxP*NyP
449     c
450     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
451     & convert='big_endian',
452     & file=trim(dirNEST)//'/PARENT/RAS.data',
453     & form='unformatted')
454    
455     read (1,REC=1) RAS_P(:,:)
456     close(1)
457     c----------------------------------------------------
458     c MASK x PARENT
459     c----------------------------------------------------
460     DO K = 1,NrP
461     DO J = 1,NyP
462     DO I = 1,NxP
463     DEEP_P(i,j,k) = 0.
464     IF (hFacC_P(i,j,k).ne.0) then
465     DEEP_P(I,J,K) = 1.
466     ENDIF
467     ENDDO
468     ENDDO
469     ENDDO
470     c----------------------------------------------------
471     c 1/Volume (C)
472     c----------------------------------------------------
473     DO K = 1,NrP
474     DO J = 1,NyP
475     DO I = 1,NxP
476     INV_VOL_C_P(I,J,K) = 1.
477     IF ((RAC_P(I,J)*hFacC_P(I,J,K)).ne.0.) THEN
478     INV_VOL_C_P(I,J,K) = 1./(RAC_P(I,J)*hFacC_P(I,J,K))
479     ENDIF
480     ENDDO
481     ENDDO
482     ENDDO
483     c----------------------------------------------------
484     c 1/Volume (W)
485     c----------------------------------------------------
486     DO K = 1,NrP
487     DO J = 1,NyP
488     DO I = 1,NxP
489     INV_VOL_W_P(I,J,K) = 1.
490     IF ((RAW_P(I,J)*hFacW_P(I,J,K)).ne.0.) THEN
491     INV_VOL_W_P(I,J,K) = 1./(RAW_P(I,J)*hFacW_P(I,J,K))
492     ENDIF
493     ENDDO
494     ENDDO
495     ENDDO
496     c----------------------------------------------------
497     c 1/Volume (S)
498     c----------------------------------------------------
499     DO K = 1,NrP
500     DO J = 1,NyP
501     DO I = 1,NxP
502     INV_VOL_S_P(I,J,K) = 1.
503     IF ((RAS_P(I,J)*hFacS_P(I,J,K)).ne.0.) THEN
504     INV_VOL_S_P(I,J,K) = 1./(RAS_P(I,J)*hFacS_P(I,J,K))
505     ENDIF
506     ENDDO
507     ENDDO
508     ENDDO
509     c--------------------------------------------------------
510     c CHILD MODEL
511     c--------------------------------------------------------
512     print*,' [2] Read CHILD model geometry'
513     c----------------------------------------------------
514     c XC & YC
515     c----------------------------------------------------
516     MSIZE = NxC*NyC
517     c
518     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
519     & convert='big_endian',
520     & file=trim(dirNEST)//'/CHILD/XC.data',
521     & form='unformatted')
522    
523     read (1,REC=1) Xo_C(:,:)
524     close(1)
525     c
526     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
527     & convert='big_endian',
528     & file=trim(dirNEST)//'/CHILD/YC.data',
529     & form='unformatted')
530     c
531     read (1,REC=1) Yo_C(:,:)
532     close(1)
533     c----------------------------------------------------
534     c XG & YG
535     c----------------------------------------------------
536     MSIZE = NxC*NyC
537     c
538     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
539     & convert='big_endian',
540     & file=trim(dirNEST)//'/CHILD/XG.data',
541     & form='unformatted')
542    
543     read (1,REC=1) Xg_C(:,:)
544     close(1)
545     c
546     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
547     & convert='big_endian',
548     & file=trim(dirNEST)//'/CHILD/YG.data',
549     & form='unformatted')
550    
551     read (1,REC=1) Yg_C(:,:)
552     close(1)
553     c----------------------------------------------------
554     c Xu & Yu
555     c----------------------------------------------------
556     DO J = 1,NyC
557     DO I = 1,NxC
558     Xu_C(I,J) = Xg_C(I,J)
559     Yu_C(I,J) = Yo_C(I,J)
560     ENDDO
561     ENDDO
562     c----------------------------------------------------
563     c Xv & Yv
564     c----------------------------------------------------
565     DO J = 1,NyC
566     DO I = 1,NxC
567     Xv_C(I,J) = Xo_C(I,J)
568     Yv_C(I,J) = Yg_C(I,J)
569     ENDDO
570     ENDDO
571     c----------------------------------------------------
572     c hFacC
573     c----------------------------------------------------
574     MSIZE = NxC*NyC*NrC
575     c
576     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
577     & convert='big_endian',
578     & file=trim(dirNEST)//'/CHILD/hFacC.data',
579     & form='unformatted')
580    
581     read (1,REC=1) hFacC_C(:,:,:)
582     close(1)
583     c----------------------------------------------------
584     c hFacW
585     c----------------------------------------------------
586     MSIZE = NxC*NyC*NrC
587     c
588     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
589     & convert='big_endian',
590     & file=trim(dirNEST)//'/CHILD/hFacW.data',
591     & form='unformatted')
592    
593     read (1,REC=1) hFacW_C(:,:,:)
594     close(1)
595     c----------------------------------------------------
596     c hFacC
597     c----------------------------------------------------
598     MSIZE = NxC*NyC*NrC
599     c
600     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
601     & convert='big_endian',
602     & file=trim(dirNEST)//'/CHILD/hFacS.data',
603     & form='unformatted')
604    
605     read (1,REC=1) hFacS_C(:,:,:)
606     close(1)
607     c----------------------------------------------------
608     c RAC
609     c----------------------------------------------------
610     MSIZE = NxC*NyC
611     c
612     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
613     & convert='big_endian',
614     & file=trim(dirNEST)//'/CHILD/RAC.data',
615     & form='unformatted')
616    
617     read (1,REC=1) RAC_C(:,:)
618     close(1)
619     c----------------------------------------------------
620     c RAW
621     c----------------------------------------------------
622     MSIZE = NxC*NyC
623     c
624     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
625     & convert='big_endian',
626     & file=trim(dirNEST)//'/CHILD/RAW.data',
627     & form='unformatted')
628    
629     read (1,REC=1) RAW_C(:,:)
630     close(1)
631     c----------------------------------------------------
632     c RAS
633     c----------------------------------------------------
634     MSIZE = NxC*NyC
635     c
636     open(unit=1,ACCESS='direct',RECL=MSIZE,STATUS='OLD',
637     & convert='big_endian',
638     & file=trim(dirNEST)//'/CHILD/RAS.data',
639     & form='unformatted')
640    
641     read (1,REC=1) RAS_C(:,:)
642     close(1)
643     c----------------------------------------------------
644     c MASK x CHILD
645     c----------------------------------------------------
646     DO K = 1,NrC
647     DO J = 1,NyC
648     DO I = 1,NxC
649     DEEP_C(i,j,k) = 0.
650     IF (hFacC_C(i,j,k).ne.0) then
651     DEEP_C(I,J,K) = 1.
652     ENDIF
653     ENDDO
654     ENDDO
655     ENDDO
656     c----------------------------------------------------
657     c
658     c __/________ ___________
659     c | | | | ||
660     c > o > | | |
661     c |__/__|_____|
662     c | | |
663     c > o > |
664     c |_____|_____|_____|_____|
665     c
666     c
667     c----------------------------------------------------
668     print*,' [3] Compute J index P-->C'
669     c--------------------------------------------------------
670     c Compute J indicies for mapping P->C (I)
671     c--------------------------------------------------------
672     I = 1
673     II = WesternB
674     c
675     DO J = 1,NyC
676     P2C_U(J) = 0.
677     DO JJ = 1,NyP-1
678     YF = Yg_C(I,J)
679     YP1 = Yg_P(II,JJ)
680     YP2 = Yg_P(II,JJ+1)
681     IF (YF.ge.YP1.and.YF.lt.YP2) THEN
682     P2C_U(J) = JJ
683     ENDIF
684     ENDDO
685     ENDDO
686     c--------------------------------------------------------
687     c Compute J indicies for mapping P->C (II)
688     c--------------------------------------------------------
689     I = 1
690     II = WesternB
691     c
692     DO J = 1,NyC
693     P2C_linU(J) = 0.
694     DO JJ = 1,NyP-1
695     YF = Yu_C(I,J)
696     YP1 = Yu_P(II,JJ)
697     YP2 = Yu_P(II,JJ+1)
698     IF (YF.ge.YP1.and.YF.lt.YP2) THEN
699     P2C_linU(J) = JJ
700     ENDIF
701     ENDDO
702     ENDDO
703     c--------------------------------------------------------
704     c Compute J indicies for mapping P->C (III)
705     c--------------------------------------------------------
706     I = 1
707     II = WesternB
708     c
709     DO J = 1,NyC
710     DO JJ = 1,NyP-1
711     YF = Yu_C(I,J)
712     YP1 = Yu_P(II,JJ)
713     IF (YF.eq.YP1) THEN
714     WO3_linU(J) = 0
715     if (J+1.le.NyC) WO3_linU(J+1) = 1
716     if (J+2.le.NyC) WO3_linU(J+2) = 2
717     ENDIF
718     ENDDO
719     ENDDO
720     c--------------------Lower bound
721     DO J = 1,NyC
722     DO JJ = 1,NyP-1
723     YF = Yu_C(I,J)
724     YP1 = Yu_P(II,JJ)
725     IF (YF.eq.YP1) THEN
726     WO3_linU(J) = 0
727     if (J-1.gt.0) WO3_linU(J-1) = 2
728     if (J-2.gt.0) WO3_linU(J-2) = 1
729     GOTO 2345
730     ENDIF
731     ENDDO
732     ENDDO
733     2345 CONTINUE
734     c--------------------Upper bound
735     DO J = NyC,1,-1
736     DO JJ = 1,NyP-1
737     YF = Yu_C(I,J)
738     YP1 = Yu_P(II,JJ)
739     IF (YF.eq.YP1) THEN
740     WO3_linU(J) = 0
741     if (J+1.le.NyC) WO3_linU(J+1) = 1
742     if (J+2.le.NyC) WO3_linU(J+2) = 2
743     GOTO 2346
744     ENDIF
745     ENDDO
746     ENDDO
747     2346 CONTINUE
748     c--------------------------------------------------------
749     c Compute J indicies for mapping P->C (IV)
750     c--------------------------------------------------------
751     I = 1
752     II = WesternB
753     c
754     DO J = 1,NyC
755     P2C_linV(J) = 0.
756     DO JJ = 1,NyP-1
757     YF = Yv_C(I,J)
758     YP1 = Yv_P(II,JJ)
759     YP2 = Yv_P(II,JJ+1)
760     IF (YF.ge.YP1.and.YF.lt.YP2) THEN
761     P2C_linV(J) = JJ
762     ENDIF
763     ENDDO
764     ENDDO
765     c--------------------------------------------------------
766     c Compute J indicies for mapping P->C (V)
767     c--------------------------------------------------------
768     I = 1
769     II = WesternB
770     c
771     DO J = 1,NyC
772     DO JJ = 1,NyP-1
773     YF = Yv_C(I,J)
774     YP1 = Yv_P(II,JJ)
775     IF (YF.eq.YP1) THEN
776     WO3_linV(J) = 0
777     if (J+1.le.NyC) WO3_linV(J+1) = 1
778     if (J+2.le.NyC) WO3_linV(J+2) = 2
779     ENDIF
780     ENDDO
781     ENDDO
782     c--------------------Lower bound
783     DO J = 1,NyC
784     DO JJ = 1,NyP-1
785     YF = Yv_C(I,J)
786     YP1 = Yv_P(II,JJ)
787     IF (YF.eq.YP1) THEN
788     WO3_linV(J) = 0
789     if (J-1.gt.0) WO3_linV(J-1) = 2
790     if (J-2.gt.0) WO3_linV(J-2) = 1
791     GOTO 23451
792     ENDIF
793     ENDDO
794     ENDDO
795     23451 CONTINUE
796     c--------------------Upper bound
797     DO J = NyC,1,-1
798     DO JJ = 1,NyP-1
799     YF = Yv_C(I,J)
800     YP1 = Yv_P(II,JJ)
801     IF (YF.eq.YP1) THEN
802     WO3_linV(J) = 0
803     if (J+1.le.NyC) WO3_linV(J+1) = 1
804     if (J+2.le.NyC) WO3_linV(J+2) = 2
805     GOTO 23461
806     ENDIF
807     ENDDO
808     ENDDO
809     23461 CONTINUE
810     c--------------------------------------------------------
811     c Compute J indicies for mapping P->C (V)
812     c--------------------------------------------------------
813     print*,' [5] Compute J index P-->C for (o)'
814     I = 1
815     II = WesternB
816     c
817     DO J = 1,NyC
818     P2C_o(J) = 0.
819     DO JJ = 1,NyP-1
820     YF = Yo_C(I,J)
821     YP1 = Yg_P(II,JJ)
822     YP2 = Yg_P(II,JJ+1)
823     IF (YF.gt.YP1.and.YF.lt.YP2) THEN
824     P2C_o(J) = JJ
825     ENDIF
826     ENDDO
827     ENDDO
828     c--------------------------------------------------------
829     c Compute J indicies for mapping P->C (VI)
830     c--------------------------------------------------------
831     print*,' [6] Compute J index P-->C for (v bilinear)'
832     I = 1
833     II = WesternB
834     c
835     DO J = 1,NyC
836     DO JJ = 2,NyP-1
837     YF = Yv_C(I,J)
838     YP1 = Yv_P(II,JJ)
839     YP2 = Yv_P(II,JJ+1)
840     YP3 = Yv_P(II,JJ-1)
841     c
842     IF (YF.ge.YP1.and.YF.lt.YP2) THEN
843     P2C1_V(J) = JJ
844     P2C2_V(J) = JJ+1
845     ENDIF
846     ENDDO
847     ENDDO
848     c--------------------------------------------------------
849     c Look for the 9 CHILD indicies in PARENT grid cell
850     c--------------------------------------------------------
851     print*,' [8] Compute I J index C-->P for (o)'
852     c
853     DO J = 1,NyP
854     DO I = 1,NxP
855     I_C2P(:,I,J) = 0
856     J_C2P(:,I,J) = 0
857     c
858     DO JJ = 1,NyC
859     DO II = 1,NxC
860     IF (Xo_C(II,JJ).eq.Xo_P(I,J).and.
861     & Yo_C(II,JJ).eq.Yo_P(I,J)) then
862    
863     KK = 0
864     DO JJJ = JJ-1,JJ+1
865     DO III = II-1,II+1
866     KK = kk +1
867     if (III.lt.1.or.III.gt.NxC) cycle
868     if (JJJ.lt.1.or.JJJ.gt.NyC) cycle
869     I_C2P(KK,I,J) = III
870     J_C2P(KK,I,J) = JJJ
871     ENDDO
872     ENDDO
873     ENDIF
874    
875     ENDDO
876     ENDDO
877     c
878     ENDDO
879     ENDDO
880     ENDIF
881     c--------------------------------------------------------
882     c Broadcast all the above variables
883     c--------------------------------------------------------
884     call MPI_BCAST(I_C2P,9*NxP*NyP,MPI_INTEGER,
885     & 0,NEST_COMM,ierr)
886     call MPI_BCAST(J_C2P,9*NxP*NyP,MPI_INTEGER,
887     & 0,NEST_COMM,ierr)
888    
889     call MPI_BCAST(RAC_C,NxC*NyC,MPI_REAL,
890     & 0,NEST_COMM,ierr)
891     call MPI_BCAST(hFacC_C,NxC*NyC*NrC,MPI_REAL,
892     & 0,NEST_COMM,ierr)
893     call MPI_BCAST(INV_VOL_C_P,NxP*NyP*NrP,MPI_REAL,
894     & 0,NEST_COMM,ierr)
895    
896     call MPI_BCAST(RAW_C,NxC*NyC,MPI_REAL,
897     & 0,NEST_COMM,ierr)
898     call MPI_BCAST(hFacW_C,NxC*NyC*NrC,MPI_REAL,
899     & 0,NEST_COMM,ierr)
900     call MPI_BCAST(INV_VOL_W_P,NxP*NyP*NrP,MPI_REAL,
901     & 0,NEST_COMM,ierr)
902    
903     call MPI_BCAST(RAS_C,NxC*NyC,MPI_REAL,
904     & 0,NEST_COMM,ierr)
905     call MPI_BCAST(hFacS_C,NxC*NyC*NrC,MPI_REAL,
906     & 0,NEST_COMM,ierr)
907     call MPI_BCAST(INV_VOL_S_P,NxP*NyP*NrP,MPI_REAL,
908     & 0,NEST_COMM,ierr)
909     c
910     call MPI_BCAST(DEEP_C,NxC*NyC*NrC,MPI_REAL,
911     & 0,NEST_COMM,ierr)
912     call MPI_BCAST(RAC_P,NxP*NyP,MPI_REAL,
913     & 0,NEST_COMM,ierr)
914     c
915     call MPI_BCAST(IM_P,1,MPI_INTEGER,
916     & 0,NEST_COMM,ierr)
917     call MPI_BCAST(JM_P,1,MPI_INTEGER,
918     & 0,NEST_COMM,ierr)
919     call MPI_BCAST(index_var3D,1,MPI_INTEGER,
920     & 0,NEST_COMM,ierr)
921     call MPI_BCAST(index_var2D,1,MPI_INTEGER,
922     & 0,NEST_COMM,ierr)
923     c
924     call MPI_BCAST(DEEP_P,NxP*NyP*NrP,MPI_REAL,
925     & 0,NEST_COMM,ierr)
926     call MPI_BCAST(hFacS_P,NxP*NyP*NrP,MPI_REAL,
927     & 0,NEST_COMM,ierr)
928     call MPI_BCAST(hFacC_P,NxP*NyP*NrP,MPI_REAL,
929     & 0,NEST_COMM,ierr)
930     call MPI_BCAST(hFacW_P,NxP*NyP*NrP,MPI_REAL,
931     & 0,NEST_COMM,ierr)
932    
933     c--------------------------------------------------------
934     if(rank.eq.0) then
935     c--------------------------------------------------------
936     DO K = 1,NrP
937     DO J = 1,NyP
938     DO I = WesternB+1,EasternB-1
939     cc WesternB side
940    
941     DO II = 1,9
942     IF (I_C2P(II,I,J).eq.0) cycle
943     IF (J_C2P(II,I,J).eq.0) cycle
944     c
945     Indx = I_C2P(II,I,J)
946     Jndx = J_C2P(II,I,J)
947     ENDDO
948     ENDDO
949     ENDDO
950     ENDDO
951     c---------------------------------------------------------
952     ONOFF=0
953     endif
954     c--------------------------------------------------------
955     c BEGIN MAIN LOOP
956     c--------------------------------------------------------
957     do
958     if(rank.eq.0) then
959     c--------------------------------------------------------
960     c (1) READ FROM PARENT MODEL
961     c--------------------------------------------------------
962     ICONT=1
963     DO WHILE(ICONT.le.nSxP*nSyP)
964     from= MPI_ANY_SOURCE
965    
966     call MPI_RECV (globalPA, index, MPI_REAL8,
967     & FROM, 3000,
968     & MPI_COMM_World, status,ierr)
969     c
970     ICONT=ICONT+1
971     c
972     whm=status(MPI_SOURCE)-MSTR_PRNT(NST_LEV)+1
973     c
974     call MPI_GET_COUNT(status,MPI_REAL8,st_count,ierr)
975     c
976     DO II = 1,6
977     IF (globalPA(II,1,1,1).ne.-999.) THEN
978     globalP1(II,1+IndJ_P(whm):JM_P+IndJ_P(whm),:) =
979     & globalPA(II,1:JM_P,:,1)
980     globalP2(II,1+IndJ_P(whm):JM_P+IndJ_P(whm),:) =
981     & globalPA(II,1:JM_P,:,2)
982     globalP3(II,1+IndJ_P(whm):JM_P+IndJ_P(whm),:) =
983     & globalPA(II,1:JM_P,:,3)
984     globalP4(II,1+IndJ_P(whm):JM_P+IndJ_P(whm),:) =
985     & globalPA(II,1:JM_P,:,4)
986     globalP5(II,1+IndJ_P(whm):JM_P+IndJ_P(whm),:) =
987     & globalPA(II,1:JM_P,:,5)
988     ENDIF
989    
990     ENDDO
991     ENDDO
992     c--------------------------------------------------------
993     c Start interpolation for CHILD
994     c--------------------------------------------------------
995     CALL INTERPOLATION_P2C (
996     & globalP1,globalP2,globalP3,globalP4,globalP5,
997     & NxP,NyP,NrP,
998     & NxC,NyC,NrC,
999     $ WesternB,EasternB,
1000     $ P2C_U,P2C_V,P2C_o,P2C1_V,P2C2_V,P2C1_o,P2C2_o,
1001     $ P2C_linU,WO3_linU,P2C_linV,WO3_linV,
1002     $ Xv_C,Yv_C,Xv_P,Yv_P,
1003     $ T_C1,S_C1,U_C1,V_C1,ETA_C1,
1004     $ DEEP_C,DEEP_P
1005     & )
1006     c==============================================================
1007     c Open Files from PARENT MODEL
1008     c==============================================================
1009     ICONT=1
1010    
1011     do while(ICONT.le.nSxP*nSyP)
1012     from= MPI_ANY_SOURCE
1013    
1014     call MPI_RECV (globalPA, index, MPI_REAL8,
1015     & FROM, 3000,
1016     & MPI_COMM_World, status,ierr)
1017    
1018    
1019    
1020     ICONT=ICONT+1
1021    
1022     whm=status(MPI_SOURCE)-MSTR_PRNT(NST_LEV)+1
1023    
1024     call MPI_GET_COUNT(status,MPI_REAL8,st_count,ierr)
1025    
1026     DO II = 1,6
1027     IF (globalPA(II,1,1,1).ne.-999.) THEN
1028     globalP1(II,1+IndJ_P(whm):JM_P+IndJ_P(whm),:) =
1029     & globalPA(II,1:JM_P,:,1)
1030     globalP2(II,1+IndJ_P(whm):JM_P+IndJ_P(whm),:) =
1031     & globalPA(II,1:JM_P,:,2)
1032     globalP3(II,1+IndJ_P(whm):JM_P+IndJ_P(whm),:) =
1033     & globalPA(II,1:JM_P,:,3)
1034     globalP4(II,1+IndJ_P(whm):JM_P+IndJ_P(whm),:) =
1035     & globalPA(II,1:JM_P,:,4)
1036     globalP5(II,1+IndJ_P(whm):JM_P+IndJ_P(whm),:) =
1037     & globalPA(II,1:JM_P,:,5)
1038     ENDIF
1039     ENDDO
1040    
1041     end do
1042     c--------------------------------------------------------
1043     c Start inteprolation for CHILD
1044     c--------------------------------------------------------
1045     CALL INTERPOLATION_P2C (
1046     & globalP1,globalP2,globalP3,globalP4,globalP5,
1047     & NxP,NyP,NrP,
1048     & NxC,NyC,NrC,
1049     $ WesternB,EasternB,
1050     $ P2C_U,P2C_V,P2C_o,P2C1_V,P2C2_V,P2C1_o,P2C2_o,
1051     $ P2C_linU,WO3_linU,P2C_linV,WO3_linV,
1052     $ Xv_C,Yv_C,Xv_P,Yv_P,
1053     $ T_C2,S_C2,U_C2,V_C2,ETA_C2,
1054     $ DEEP_C,DEEP_P
1055     & )
1056    
1057    
1058    
1059     c==============================================================
1060     c Temporal Interpolation OBCs x CHILD MODEL
1061     c==============================================================
1062     c 0 1200
1063     c ---+--.--.--+---- Parent
1064     c
1065     c |--|--|--
1066     c 0 800
1067     c 400
1068     c------------------------------------------------------------
1069     DO I = 1,2 ! WesternB & EasternB
1070     DIFF_T(:,:,I) = (T_C2(:,:,I) - T_C1(:,:,I))/3.
1071     DIFF_S(:,:,I) = (S_C2(:,:,I) - S_C1(:,:,I))/3.
1072     DIFF_U(:,:,I) = (U_C2(:,:,I) - U_C1(:,:,I))/3.
1073     DIFF_V(:,:,I) = (V_C2(:,:,I) - V_C1(:,:,I))/3.
1074     DIFF_eta(:,:,I) = (eta_C2(:,:,I) - eta_C1(:,:,I))/3.
1075     ENDDO
1076     c-------------------------------------------------------------
1077     c
1078     c Step 0 (Rec = 1 ==> WesternB)
1079     c------- (Rec = 2 ==> EasternB)
1080    
1081     DO I = 1,2 !WesternB & EasternB
1082     T_C1(:,:,I) = T_C2(:,:,I) !+ DIFF_T(:,:,I)
1083     S_C1(:,:,I) = S_C2(:,:,I) !+ DIFF_S(:,:,I)
1084     U_C1(:,:,I) = U_C2(:,:,I) !+ DIFF_U(:,:,I)
1085     V_C1(:,:,I) = V_C2(:,:,I) !+ DIFF_V(:,:,I)
1086     ETA_C1(:,:,I) = ETA_C2(:,:,I) !+ DIFF_ETA(:,:,I)
1087     ENDDO
1088    
1089    
1090     if(ONOFF.eq.0) then
1091     c---------------------------------------------------------------------
1092     ICONT = -1
1093     DO I = 1,nSxC
1094     DO J = 1,nSyC
1095     ICONT = ICONT + 1
1096     IndI = IM_C*(I-1)
1097     IndJ = JM_C*(J-1)
1098    
1099     VAR_C1(:,:,:,:) = 0.
1100     c
1101     J1 = 1+IndJ-OLY
1102     J2 = JM_C+IndJ+OLY
1103     c
1104     JJ1 = 1
1105     JJ2 = JM_C+OLY+OLY
1106     c
1107     IF(1 +IndJ-OLY.LT.0) THEN
1108     J1 = 1
1109     JJ1 = 4
1110     ENDIF
1111     c
1112     IF(JM_C+IndJ+OLY.GT.NyC) THEN
1113     J2 = NyC
1114     JJ2 = JM_C
1115     ENDIF
1116     c
1117     VAR_C1(JJ1:JJ2,:,:,1) = U_C1(J1:J2,:,:)
1118     VAR_C1(JJ1:JJ2,:,:,2) = V_C1(J1:J2,:,:)
1119     VAR_C1(JJ1:JJ2,:,:,3) = T_C1(J1:J2,:,:)
1120     VAR_C1(JJ1:JJ2,:,:,4) = S_C1(J1:J2,:,:)
1121     VAR_C1(JJ1:JJ2,:,:,5) = ETA_C1(J1:J2,:,:)
1122     c
1123     call MPI_SEND (VAR_C1, indexF, MPI_REAL8,
1124     & MSTR_CHLD(NST_LEV)+ICONT, 3000,
1125     & MPI_Comm_World,ierr)
1126     c
1127     ENDDO
1128     ENDDO
1129     c----------------------------------------------------------------------
1130     cc write(*,*) 'VIC: MANDO SEGNALE DI OK AL CHILD PER INIZIALIZZARE'
1131     ONOFF=1
1132     ENDIF
1133     cc write(*,*) 'VIC: MANDO SEGNALE DI OK AL CHILD PER IL PASSO 1'
1134     c-----------------------------------------------------------------------
1135     ICONT = -1
1136     DO I = 1,nSxC
1137     DO J = 1,nSyC
1138     ICONT = ICONT + 1
1139     IndI = IM_C*(I-1)
1140     IndJ = JM_C*(J-1)
1141     c
1142     VAR_C1(:,:,:,:) = 0.
1143     c
1144     J1 = 1+IndJ-OLY
1145     J2 = JM_C+IndJ+OLY
1146     c
1147     JJ1 = 1
1148     JJ2 = JM_C+OLY+OLY
1149     c
1150     IF(1 +IndJ-OLY.LT.0) THEN
1151     J1 = 1
1152     JJ1 = 4
1153     ENDIF
1154     c
1155     IF(JM_C+IndJ+OLY.GT.NyC) THEN
1156     J2 = NyC
1157     JJ2 = JM_C
1158     ENDIF
1159     c
1160     VAR_C1(JJ1:JJ2,:,:,1) = U_C1(J1:J2,:,:)
1161     VAR_C1(JJ1:JJ2,:,:,2) = V_C1(J1:J2,:,:)
1162     VAR_C1(JJ1:JJ2,:,:,3) = T_C1(J1:J2,:,:)
1163     VAR_C1(JJ1:JJ2,:,:,4) = S_C1(J1:J2,:,:)
1164     VAR_C1(JJ1:JJ2,:,:,5) = ETA_C1(J1:J2,:,:)
1165     c
1166     call MPI_SEND (VAR_C1, indexF, MPI_REAL8,
1167     & MSTR_CHLD(NST_LEV)+ICONT, 3000,
1168     & MPI_Comm_World,ierr)
1169    
1170     ENDDO
1171     ENDDO
1172     c---------------------------------------------------------------------
1173     goto 8888
1174    
1175     c
1176     c Step 1 (Rec = 3 ==> WesternB)
1177     c------- (Rec = 4 ==> EasternB)
1178    
1179     DO I = 1,2 !WesternB & EasternB
1180     T_C1(:,:,I) = T_C2(:,:,I) !+ DIFF_T(:,:,I)
1181     S_C1(:,:,I) = S_C2(:,:,I) !+ DIFF_S(:,:,I)
1182     U_C1(:,:,I) = U_C2(:,:,I) !+ DIFF_U(:,:,I)
1183     V_C1(:,:,I) = V_C2(:,:,I) !+ DIFF_V(:,:,I)
1184     ETA_C1(:,:,I) = ETA_C2(:,:,I) !+ DIFF_ETA(:,:,I)
1185     ENDDO
1186     c----------------------------------------------------------
1187     ICONT = -1
1188     DO I = 1,nSxC
1189     DO J = 1,nSyC
1190     ICONT = ICONT + 1
1191     IndI = IM_C*(I-1)
1192     IndJ = JM_C*(J-1)
1193    
1194     VAR_C1(:,:,:,:) = 0.
1195     c
1196     J1 = 1+IndJ-OLY
1197     J2 = JM_C+IndJ+OLY
1198     c
1199     JJ1 = 1
1200     JJ2 = JM_C+OLY+OLY
1201     c
1202     IF(1 +IndJ-OLY.LT.0) THEN
1203     J1 = 1
1204     JJ1 = 4
1205     ENDIF
1206     c
1207     IF(JM_C+IndJ+OLY.GT.NyC) THEN
1208     J2 = NyC
1209     JJ2 = JM_C
1210     ENDIF
1211     c
1212     VAR_C1(JJ1:JJ2,:,:,1) = U_C1(J1:J2,:,:)
1213     VAR_C1(JJ1:JJ2,:,:,2) = V_C1(J1:J2,:,:)
1214     VAR_C1(JJ1:JJ2,:,:,3) = T_C1(J1:J2,:,:)
1215     VAR_C1(JJ1:JJ2,:,:,4) = S_C1(J1:J2,:,:)
1216     VAR_C1(JJ1:JJ2,:,:,5) = ETA_C1(J1:J2,:,:)
1217    
1218     call MPI_SEND (VAR_C1, indexF, MPI_REAL8,
1219     & MSTR_CHLD(NST_LEV)+ICONT, 3000,
1220     & MPI_Comm_World,ierr)
1221    
1222    
1223     ENDDO
1224     ENDDO
1225     c-----------------------------------------------------------
1226     c
1227     c Step 2 (Rec = 5 ==> WesternB)
1228     c------- (Rec = 6 ==> EasternB)
1229    
1230     DO I = 1,2 !WesternB & EasternB
1231     T_C1(:,:,I) = T_C2(:,:,I) !+ DIFF_T(:,:,I)
1232     S_C1(:,:,I) = S_C2(:,:,I) !+ DIFF_S(:,:,I)
1233     U_C1(:,:,I) = U_C2(:,:,I) !+ DIFF_U(:,:,I)
1234     V_C1(:,:,I) = V_C2(:,:,I) !+ DIFF_V(:,:,I)
1235     ETA_C1(:,:,I) = ETA_C2(:,:,I) !+ DIFF_ETA(:,:,I)
1236     ENDDO
1237     c----------------------------------------------------------
1238     ICONT = -1
1239     DO I = 1,nSxC
1240     DO J = 1,nSyC
1241     ICONT = ICONT + 1
1242     IndI = IM_C*(I-1)
1243     IndJ = JM_C*(J-1)
1244    
1245     VAR_C1(:,:,:,:) = 0.
1246     c
1247     J1 = 1+IndJ-OLY
1248     J2 = JM_C+IndJ+OLY
1249     c
1250     JJ1 = 1
1251     JJ2 = JM_C+OLY+OLY
1252     c
1253     IF(1 +IndJ-OLY.LT.0) THEN
1254     J1 = 1
1255     JJ1 = 4
1256     ENDIF
1257     c
1258     IF(JM_C+IndJ+OLY.GT.NyC) THEN
1259     J2 = NyC
1260     JJ2 = JM_C
1261     ENDIF
1262     c
1263     VAR_C1(JJ1:JJ2,:,:,1) = U_C1(J1:J2,:,:)
1264     VAR_C1(JJ1:JJ2,:,:,2) = V_C1(J1:J2,:,:)
1265     VAR_C1(JJ1:JJ2,:,:,3) = T_C1(J1:J2,:,:)
1266     VAR_C1(JJ1:JJ2,:,:,4) = S_C1(J1:J2,:,:)
1267     VAR_C1(JJ1:JJ2,:,:,5) = ETA_C1(J1:J2,:,:)
1268    
1269     call MPI_SEND (VAR_C1, indexF, MPI_REAL8,
1270     & MSTR_CHLD(NST_LEV)+ICONT, 3000,
1271     & MPI_Comm_World,ierr)
1272    
1273    
1274     ENDDO
1275     ENDDO
1276     c---------------------------------------------------------------
1277     8888 CONTINUE
1278     c--------------------------------------------------------
1279     c Close OBCs Files x CHILD MODEL
1280     c--------------------------------------------------------
1281     c----------------------------------------------------
1282     c------------- MANDO SEGNALE DI OK AL CHILD
1283     c----------------------------------------------------
1284     c--------------------------------------------------------
1285     c (1) READ FROM CHILD MODEL
1286     c--------------------------------------------------------
1287     call MPI_RECV (TRANSPORT_WEST, 1, MPI_REAL8,
1288     & MSTR_CHLD(NST_LEV), 3000,
1289     & MPI_COMM_World, status,ierr)
1290    
1291    
1292     call MPI_RECV (TRANSPORT_EAST, 1, MPI_REAL8,
1293     & MSTR_CHLD(NST_LEV), 3000,
1294     & MPI_COMM_World, status,ierr)
1295     c---------------------------------------------------------
1296     c---------------------------------------------------------
1297     ICONT=1
1298    
1299     DO WHILE(ICONT.le.nSxC*nSyC)
1300     from= MPI_ANY_SOURCE
1301     call MPI_RECV (globalC3D_a,indexF, MPI_REAL8,
1302     & from, 3000, MPI_COMM_World, status,ierr)
1303    
1304     ICONT=ICONT+1
1305    
1306     whm=status(MPI_SOURCE)-MSTR_CHLD(NST_LEV)+1
1307    
1308     call MPI_GET_COUNT(status,MPI_REAL8,st_count,ierr)
1309    
1310     globalC3D(1+IndI_C(whm):IM_C+IndI_C(whm),
1311     & 1+IndJ_C(whm):JM_C+IndJ_C(whm),:,:)=
1312     & globalC3D_a(:,:,:,:)
1313     END DO
1314     c-----------------------------
1315     ICONT=1
1316    
1317     DO WHILE(ICONT.le.nSxC*nSyC)
1318     from= MPI_ANY_SOURCE
1319     call MPI_RECV (globalC2D_a,index1F, MPI_REAL8,
1320     & from, 3000, MPI_COMM_World, status,ierr)
1321    
1322     ICONT=ICONT+1
1323    
1324     whm=status(MPI_SOURCE)-MSTR_CHLD(NST_LEV)+1
1325    
1326     call MPI_GET_COUNT(status,MPI_REAL8,st_count,ierr)
1327    
1328    
1329     globalC2D(1+IndI_C(whm):IM_C+IndI_C(whm),
1330     & 1+IndJ_C(whm):JM_C+IndJ_C(whm),:)=
1331     & globalC2D_a(:,:,:)
1332     END DO
1333     ENDIF
1334    
1335    
1336    
1337    
1338     call MPI_BCAST(globalC3D,NxC*NyC*NrC*15,MPI_REAL8,
1339     & 0,NEST_COMM,ierr)
1340     call MPI_BCAST(globalC2D,NxC*NyC*4,MPI_REAL8,
1341     & 0,NEST_COMM,ierr)
1342    
1343     2323 CONTINUE
1344     c=======================================================
1345     c (1) READ FROM CHILD MODEL
1346     c=======================================================
1347    
1348     c=======================================================
1349     c (2) INTERPOLATIONS
1350     c=======================================================
1351    
1352     c 3D VAR
1353     c--------
1354     DO iVar = 1,15 ! tipo di variabile
1355     DO K = 1,NrP
1356     DO J = 1,NyP
1357     DO I = WesternB+1,EasternB-1
1358     VAR3D_P(I,J,K,iVar) = 0. ! inizializzo
1359     c WesternB side
1360    
1361     AREA_VOL = 0. !can be area or volume depend on the variable
1362    
1363     SELECT CASE(iVar)
1364     CASE(1,5,9)
1365     I_START = 1
1366     I_END = 9
1367     I_STEP = 1 !3
1368     CASE(2,6,10)
1369     I_START = 1
1370     I_END = 9 !3
1371     I_STEP = 1
1372     CASE DEFAULT
1373     I_START = 1
1374     I_END = 9
1375     I_STEP = 1
1376     END SELECT
1377    
1378    
1379     DO II = I_START,I_END,I_STEP
1380    
1381    
1382     IF (I_C2P(II,I,J).eq.0) cycle
1383     IF (J_C2P(II,I,J).eq.0) cycle
1384     c
1385     Indx = I_C2P(II,I,J)
1386     Jndx = J_C2P(II,I,J)
1387     c
1388     c
1389     SELECT CASE(iVar)
1390    
1391     CASE (1,5,9)
1392     VAR3D_P(I,J,K,ivar) = VAR3D_P(I,J,K,iVar) +
1393     & globalC3D(Indx,Jndx,K,iVar)*
1394     $ RAW_C(Indx,Jndx)*
1395     & hFacW_C(Indx,Jndx,K)
1396    
1397     CASE (2,6,10)
1398     VAR3D_P(I,J,K,ivar) = VAR3D_P(I,J,K,iVar) +
1399     & globalC3D(Indx,Jndx,K,iVar)*
1400     $ RAS_C(Indx,Jndx)*
1401     & hFacS_C(Indx,Jndx,K)
1402    
1403    
1404     CASE DEFAULT
1405     VAR3D_P(I,J,K,ivar) = VAR3D_P(I,J,K,iVar) +
1406     & globalC3D(Indx,Jndx,K,iVar)*
1407     $ RAC_C(Indx,Jndx)*
1408     & hFacC_C(Indx,Jndx,K)
1409    
1410     AREA_VOL = AREA_VOL +
1411     & RAC_C(Indx,Jndx)* hFacC_C(Indx,Jndx,K)
1412    
1413     END SELECT
1414     ENDDO
1415     c-----------------------------------------------
1416     c Make a volume average
1417     c----------------------------------------------
1418     SELECT CASE(iVar)
1419     CASE (1,5,9)
1420     VAR3D_P(I,J,K,ivar) =
1421     & VAR3D_P(I,J,K,iVar)*
1422     & INV_VOL_W_P(I,J,K)
1423     if (hFacW_P(I,J,K).eq.0.) VAR3D_P(I,J,K,ivar)=0.
1424    
1425     CASE (2,6,10)
1426     VAR3D_P(I,J,K,ivar) =
1427     & VAR3D_P(I,J,K,iVar)*
1428     & INV_VOL_S_P(I,J,K)
1429     if (hFacS_P(I,J,K).eq.0.) VAR3D_P(I,J,K,ivar)=0.
1430     CASE DEFAULT
1431     IF (AREA_VOL.ne.0.) then
1432     VAR3D_P(I,J,K,ivar) =
1433     & VAR3D_P(I,J,K,iVar)/AREA_VOL
1434     ENDIF
1435     if (hFacC_P(I,J,K).eq.0.) VAR3D_P(I,J,K,ivar)=0.
1436     END SELECT
1437     ENDDO
1438     ENDDO
1439     ENDDO
1440     ENDDO
1441    
1442    
1443    
1444     c 2D VAR
1445     c--------
1446     DO iVar = 1,4
1447     DO J = 1,NyP
1448     DO I = WesternB+1,EasternB-1
1449     VAR2D_P(I,J,iVar) = 0.
1450     AREA_VOL = 0.
1451     DO II = 1,9
1452     IF (I_C2P(II,I,J).eq.0) cycle
1453     IF (J_C2P(II,I,J).eq.0) cycle
1454     c
1455     Indx = I_C2P(II,I,J)
1456     Jndx = J_C2P(II,I,J)
1457     c
1458     VAR2D_P(I,J,ivar) = VAR2D_P(I,J,iVar) +
1459     & globalC2D(Indx,Jndx,iVar)*
1460     $ RAC_C(Indx,Jndx)*
1461     & DEEP_C(Indx,Jndx,1)
1462    
1463    
1464     AREA_VOL = AREA_VOL +
1465     & RAC_C(Indx,Jndx)* DEEP_C(Indx,Jndx,1)
1466    
1467     ENDDO
1468     c-----------------------------
1469     IF ((RAC_P(I,J)*DEEP_P(I,J,1)).ne.0.) then
1470     c IF (AREA_VOL.ne.0.) then
1471    
1472     VAR2D_P(I,J,ivar) =
1473     & VAR2D_P(I,J,iVar)/
1474     & RAC_P(I,J)
1475     ENDIF
1476     c----------------------------
1477     ENDDO
1478     ENDDO
1479     ENDDO
1480     if(rank.eq.0) then
1481     c--------------------------------------------------------
1482     c Write Files for PARENT MODEL
1483     c--------------------------------------------------------
1484     c print*,' (*) Open Files for PARENT MODEL'
1485    
1486     7575 CONTINUE
1487     c----------------------------------------------------
1488     c------------- MANDO SEGNALE DI OK AL PARENT
1489     c----------------------------------------------------
1490     call MPI_SEND (TRANSPORT_WEST, 1, MPI_REAL8,
1491     & MSTR_PRNT(NST_LEV), 3000,
1492     & MPI_Comm_World,ierr)
1493    
1494     call MPI_SEND (TRANSPORT_EAST, 1, MPI_REAL8,
1495     & MSTR_PRNT(NST_LEV), 3000,
1496     & MPI_Comm_World,ierr)
1497    
1498     ENDIF
1499     c---------------------------------------------------------
1500     !---------------------------------------------------------
1501     VCONT=VCONTP(rank)
1502    
1503     DO I = vstart,vstop
1504     DO J = 1,nSyP
1505     VCONT = VCONT + 1
1506     IndI = IM_P*(I-1)
1507     IndJ = JM_P*(J-1)
1508     c-----------------------------------------------------------
1509     DO iVar=1,15
1510     CALL MPI_SEND (VAR3D_P(1+IndI:IM_P+IndI
1511     & ,1+IndJ:JM_P+IndJ,:,iVar),
1512     & index_var3D,MPI_REAL8,MSTR_PRNT(NST_LEV)+VCONT,
1513     & 3000,MPI_Comm_World,ierr)
1514    
1515     ENDDO
1516    
1517    
1518    
1519     DO iVar=1,4
1520     call MPI_SEND (VAR2D_P(1+IndI:IM_P+IndI
1521     & ,1+IndJ:JM_P+IndJ,iVar),
1522     & index_var2D,MPI_REAL8,MSTR_PRNT(NST_LEV)+VCONT,
1523     & 3000,MPI_Comm_World,ierr)
1524     ENDDO
1525    
1526     c-----------------------------------------------------------
1527     END DO
1528     END DO
1529    
1530     call MPI_BARRIER(NEST_COMM,ierr)
1531     end do
1532     c---------------------------------------------------------
1533     c=======================================================
1534     c END MAIN LOOP
1535     c=======================================================
1536     call MPI_FINALIZE(ierr)
1537     c---------------------------------------------------------
1538     !---------------------------------------------------------
1539    
1540     101 FORMAT (I1)
1541    
1542     STOP
1543     END
1544    
1545    
1546    
1547    
1548    
1549    
1550    
1551    
1552    
1553    
1554    
1555    
1556    
1557    
1558    
1559     SUBROUTINE INTERPOLATION_P2C (
1560     & globalP1,globalP2,globalP3,globalP4,globalP5,
1561     & NxP,NyP,NrP,
1562     & NxC,NyC,NrC,
1563     $ WesternB,EasternB,
1564     $ P2C_U,P2C_V,P2C_o,P2C1_V,P2C2_V,P2C1_o,P2C2_o,
1565     $ P2C_linU,WO3_linU,P2C_linV,WO3_linV,
1566     $ Xv_F,Yv_F,Xv_P,Yv_P,
1567     $ T_F,S_F,U_F,V_F,ETA_F,
1568     $ DEEP_F,DEEP_P
1569     & )
1570     c----------------------------------------------------
1571     implicit none
1572     c----------------------------------------------------
1573     INTEGER :: I,J,K,II,JJ
1574     INTEGER :: WesternB,EasternB
1575     INTEGER :: NrP,NxP,NyP
1576     INTEGER :: NrC,NxC,NyC
1577     c----------------------------------------------------
1578     REAL*8 :: Fp,Fm,Fo,VEL_MEMO
1579     INTEGER :: INDC
1580     c----------------------------------------------------
1581     c Define Global Variables to Exchange
1582     c----------------------------------------------------
1583     REAL*8 :: globalP1(6,NyP,NrP)
1584     REAL*8 :: globalP2(6,NyP,NrP)
1585     REAL*8 :: globalP3(6,NyP,NrP)
1586     REAL*8 :: globalP4(6,NyP,NrP)
1587     REAL*8 :: globalP5(6,NyP,NrP)
1588     c----------------------------------------------------
1589     c Define CHILD Model Geometry
1590     c----------------------------------------------------
1591     REAL*4 :: Xu_F(NxC,NyC)
1592     REAL*4 :: Yu_F(NxC,NyC)
1593     REAL*4 :: Xv_F(NxC,NyC)
1594     REAL*4 :: Yv_F(NxC,NyC)
1595     REAL*4 :: Xo_F(NxC,NyC)
1596     REAL*4 :: Yo_F(NxC,NyC)
1597     REAL*4 :: Xg_F(NxC,NyC)
1598     REAL*4 :: Yg_F(NxC,NyC)
1599     REAL*4 :: DEEP_F(NxC,NyC,NrC)
1600     c----------------------------------------------------
1601     c Define PARENT Model Geometry
1602     c----------------------------------------------------
1603     c REAL*4 :: Xu_P(NxP,NyP)
1604     REAL*4 :: Yu_P(NxP,NyP)
1605     REAL*4 :: Xv_P(NxP,NyP)
1606     REAL*4 :: Yv_P(NxP,NyP)
1607     REAL*4 :: Xo_P(NxP,NyP)
1608     REAL*4 :: Yo_P(NxP,NyP)
1609     REAL*4 :: Xg_P(NxP,NyP)
1610     REAL*4 :: Yg_P(NxP,NyP)
1611     REAL*4 :: DEEP_P(NxP,NyP,NrP)
1612     REAL*4 :: DEPDEP
1613     c-----------------------------------------------------
1614     REAL*4 :: X1,X2,X3,X4,Y1,Y2,Y3,Y4
1615     REAL*4 :: f1,f2,f3,f4,f,x,y
1616     REAL*4 :: gammaT,gammaS,terzo,dueterzi
1617     REAL*4 :: gammaEta
1618     REAL*4 :: gammaV
1619     c----------------------------------------------------
1620     c Define INDICIES MATRIX
1621     c----------------------------------------------------
1622     INTEGER :: P2C_U(NyC) !x imposizione NETTA
1623     INTEGER :: P2C_linU(NyC) !x Lineare
1624     INTEGER :: WO3_linU(NyC) !x Lineare !Which Of 3
1625    
1626     INTEGER :: P2C_linV(NyC) !x Lineare
1627     INTEGER :: WO3_linV(NyC) !x Lineare !Which Of 3
1628    
1629     INTEGER :: P2C_V(NyC) !x Lineare
1630     INTEGER :: P2C_o(NyC) !x Lineare
1631    
1632     INTEGER :: P2C1_V(NyC) !x BiLineare
1633     INTEGER :: P2C2_V(NyC) !x BiLineare
1634     INTEGER :: P2C1_o(NyC) !x BiLineare
1635     INTEGER :: P2C2_o(NyC) !x BiLineare
1636    
1637     REAL*8 :: diff(NrC)
1638     REAL*8 :: DEPDEP_F_WesternB(NrC)
1639     REAL*8 :: DEPDEP_F_EasternB(NrC)
1640     c----------------------------------------------------
1641     c Define CHILD model variable
1642     c----------------------------------------------------
1643     c _____________ (1) WesternB (2) EasternB
1644     c |
1645     REAL*8 :: U_F(NyC,NrC,2)
1646     REAL*8 :: V_F(NyC,NrC,2)
1647     REAL*8 :: T_F(NyC,NrC,2)
1648     REAL*8 :: S_F(NyC,NrC,2)
1649     REAL*8 :: ETA_F(NyC,NrC,2)
1650     c----------------------------------------------------
1651     PARAMETER ( terzo = 1./3.)
1652     PARAMETER (dueterzi = 2./3.)
1653     c=======================================================
1654     c (2) INTERPOLATIONS
1655     c=======================================================
1656     c (2.1) Linear for normal velocity component (u in this case)
1657     c-------------------------------------------------------
1658     DO K = 1,NrP
1659     DO J = 1,NyP-1
1660     IF (globalP4(2,J,K).eq.0.) THEN !uso la salinità come discriminante
1661     globalP1 (2,J,K) = globalP1 (2,J+1,K)
1662     ENDIF
1663     ENDDO
1664    
1665     DO J = NyP,2,-1
1666     IF (globalP4(2,J,K).eq.0.) THEN !uso la salinità come discriminante
1667     globalP1 (2,J,K) = globalP1 (2,J-1,K)
1668     ENDIF
1669     ENDDO
1670     ENDDO
1671     c=======================================================
1672     c (2.1) NOT Linear but simply imposed
1673     c-------------------------------------------------------
1674     DO J = 1,3
1675     U_F(J,:,1) = globalP1(2,P2C_U(J),:)
1676     U_F(J,:,2) = globalP1(6,P2C_U(J),:)
1677     ENDDO
1678    
1679     DO J = NyC-2,NyC
1680     U_F(J,:,1) = globalP1(2,P2C_U(J),:)
1681     U_F(J,:,2) = globalP1(6,P2C_U(J),:)
1682     ENDDO
1683     c=================================================
1684     DO J = 4,NyC-3,3
1685     INDC = P2C_U(J)
1686     DO K = 1,NrC
1687     c-------- WesternB ----------------------
1688     Fp = globalP1(2,INDC+1,K)
1689     Fo = globalP1(2,INDC,K)
1690     Fm = globalP1(2,INDC-1,K)
1691     c
1692     VEL_MEMO = 0.
1693     DO I = -1,1
1694     U_F(J+1+i,K,1) = ((Fp-2.*Fo+Fm)/24.)*
1695     & ((12.*float(i)**2+1.)/9.)+
1696     & ((float(i)/6.)*(Fp-Fm))+(26.*Fo-Fp-Fm)/24.
1697     VEL_MEMO = VEL_MEMO + U_F(J+1+i,K,1)
1698     ENDDO
1699    
1700     VEL_MEMO = ((3.*Fo) - VEL_MEMO)/3.
1701    
1702     DO I = -1,1
1703     U_F(J+1+i,K,1) = U_F(J+1+i,K,1) + VEL_MEMO
1704     ENDDO
1705    
1706     c-------- EasternB ----------------------
1707     Fp = globalP1(6,INDC+1,K)
1708     Fo = globalP1(6,INDC,K)
1709     Fm = globalP1(6,INDC-1,K)
1710    
1711     VEL_MEMO = 0.
1712     DO I = -1,1
1713     U_F(J+1+i,K,2) = ((Fp-2.*Fo+Fm)/24.)*
1714     & ((12.*float(i)**2+1.)/9.)+
1715     & ((float(i)/6.)*(Fp-Fm))+(26.*Fo-Fp-Fm)/24.
1716     VEL_MEMO = VEL_MEMO + U_F(J+1+i,K,2)
1717     ENDDO
1718    
1719     VEL_MEMO = ((3.*Fo) - VEL_MEMO)/3.
1720    
1721     DO I = -1,1
1722     U_F(J+1+i,K,2) = U_F(J+1+i,K,2) + VEL_MEMO
1723     ENDDO
1724     c------------------------------------------------------
1725     ENDDO
1726     ENDDO
1727     c-------------------------------------------------------
1728     c (2.2) BiLinear for tangent velocity component (v in this case)
1729     c-------------------------------------------------------
1730     I = 1
1731     II = WesternB
1732    
1733     V_F(:,:,:) = 0.
1734    
1735     DO K = 1,NrC
1736     DO J = 1,NyC
1737     x1 = Xv_P(II,P2C1_V(J))
1738     x2 = Xv_P(II,P2C2_V(J))
1739    
1740     x3 = Xv_P(II+1,P2C1_V(J))
1741     x4 = Xv_P(II+1,P2C2_V(J))
1742    
1743     y1 = Yv_P(II,P2C1_V(J))
1744     y2 = Yv_P(II,P2C2_V(J))
1745    
1746     y3 = Yv_P(II+1,P2C1_V(J))
1747     y4 = Yv_P(II+1,P2C2_V(J))
1748    
1749     x = Xv_F(I,J)
1750     y = Yv_F(I,J)
1751    
1752     f1 = globalP2(1,P2C1_V(J),K)
1753     f2 = globalP2(1,P2C2_V(J),K)
1754     f3 = globalP2(2,P2C1_V(J),K)
1755     f4 = globalP2(2,P2C2_V(J),K)
1756    
1757     call blint(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2,f3,f4,x,y,f)
1758     V_F(J,K,1) = f
1759     ENDDO
1760     ENDDO
1761     c..............................................
1762     I = NxC
1763     II = EasternB
1764    
1765     DO K = 1,NrC
1766     DO J = 1,NyC
1767     x1 = Xv_P(II,P2C1_V(J))
1768     x2 = Xv_P(II,P2C2_V(J))
1769    
1770     x3 = Xv_P(II-1,P2C1_V(J))
1771     x4 = Xv_P(II-1,P2C2_V(J))
1772    
1773     y1 = Yv_P(II,P2C1_V(J))
1774     y2 = Yv_P(II,P2C2_V(J))
1775    
1776     y3 = Yv_P(II-1,P2C1_V(J))
1777     y4 = Yv_P(II-1,P2C2_V(J))
1778    
1779     x = Xv_F(I,J)
1780     y = Yv_F(I,J)
1781    
1782     f1 = globalP2(6,P2C1_V(J),K)
1783     f2 = globalP2(6,P2C2_V(J),K)
1784     f3 = globalP2(5,P2C1_V(J),K)
1785     f4 = globalP2(5,P2C2_V(J),K)
1786    
1787     call blint(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2,f3,f4,x,y,f)
1788     V_F(J,K,2) = f
1789     ENDDO
1790     ENDDO
1791     c-------------------------------------------------------
1792     c (2.3.1) Linear
1793     c-------------------------------------------------------
1794     c
1795     c WesternB
1796     c
1797     DO K = 1,NrP
1798     DO J = 1,NyP
1799     c
1800     DEPDEP = DEEP_P(WesternB,J,K) * DEEP_P(WesternB+1,J,K)
1801     c
1802     gammaT =(globalP3(2,J,K)-globalP3(1,J,K))*DEPDEP
1803     gammaS =(globalP4(2,J,K)-globalP4(1,J,K))*DEPDEP
1804     gammaeta =(globalP5(2,J,K)-globalP5(1,J,K))*DEPDEP
1805     cgm-------
1806     c gammaV =(globalP2(2,J,K)-globalP2(1,J,K))*DEPDEP
1807     cgm---------
1808     globalP3(1,J,K) = globalP3(1,J,K) + terzo* gammaT
1809     globalP4(1,J,K) = globalP4(1,J,K) + terzo* gammaS
1810     globalP5(1,J,K) = globalP5(1,J,K) + terzo* gammaeta
1811     cgm-------
1812     c globalP2(1,J,K) = globalP2(1,J,K) + terzo* gammaV
1813     cgm---------
1814     ENDDO
1815     c--------------------------------------------------------------
1816     DO J = 1,NyP-1
1817     IF (globalP4(1,J,K).eq.0.) THEN !uso la salinità come discriminante
1818     globalP3(1,J,K) = globalP3(1,J+1,K)
1819     globalP4(1,J,K) = globalP4(1,J+1,K)
1820     globalP5(1,J,K) = globalP5(1,J+1,K)
1821     ENDIF
1822     ENDDO
1823    
1824     DO J = NyP,2,-1
1825     IF (globalP4(1,J,K).eq.0.) THEN !uso la salinità come discriminante
1826     globalP3(1,J,K) = globalP3(1,J-1,K)
1827     globalP4(1,J,K) = globalP4(1,J-1,K)
1828     globalP5(1,J,K) = globalP5(1,J-1,K)
1829     cgm---------
1830     c globalP2(1,J,K) = globalP2(1,J-1,K)
1831     cgm-------------
1832     ENDIF
1833     ENDDO
1834     c---------------------------------------------------------------
1835     ENDDO
1836    
1837    
1838    
1839     c
1840     c EasternB
1841     c
1842     DO K = 1,NrP
1843     DO J = 1,NyP
1844     c
1845     DEPDEP = DEEP_P(EasternB,J,K) * DEEP_P(EasternB-1,J,K)
1846     c
1847     gammaT =(globalP3(5,J,K)-globalP3(6,J,K))*DEPDEP
1848     gammaS =(globalP4(5,J,K)-globalP4(6,J,K))*DEPDEP
1849     gammaeta =(globalP5(5,J,K)-globalP5(6,J,K))*DEPDEP
1850     cgm-------------
1851     c gammaV =(globalP2(5,J,K)-globalP2(6,J,K))*DEPDEP
1852     cgm----------------
1853     globalP3(6,J,K) = globalP3(6,J,K) + terzo* gammaT
1854     globalP4(6,J,K) = globalP4(6,J,K) + terzo* gammaS
1855     globalP5(6,J,K) = globalP5(6,J,K) + terzo* gammaeta
1856     cgm----------------
1857     c globalP2(6,J,K) = globalP2(6,J,K) + terzo* gammaV
1858     cgm----------------
1859     ENDDO
1860     c--------------------------------------------------------------
1861     DO J = 1,NyP-1
1862     IF (globalP4(6,J,K).eq.0.) THEN !uso la salinità come discriminante
1863     globalP3(6,J,K) = globalP3(6,J+1,K)
1864     globalP4(6,J,K) = globalP4(6,J+1,K)
1865     globalP5(6,J,K) = globalP5(6,J+1,K)
1866     cgm----------------
1867     c globalP2(6,J,K) = globalP2(6,J+1,K)
1868     cgm----------------
1869     ENDIF
1870     ENDDO
1871    
1872     DO J = NyP,2,-1
1873     IF (globalP4(6,J,K).eq.0.) THEN !uso la salinità come discriminante
1874     globalP3(6,J,K) = globalP3(6,J-1,K)
1875     globalP4(6,J,K) = globalP4(6,J-1,K)
1876     globalP5(6,J,K) = globalP5(6,J-1,K)
1877     cgm----------------
1878     c globalP2(6,J,K) = globalP2(6,J-1,K)
1879     cgm----------------
1880     ENDIF
1881     ENDDO
1882     c---------------------------------------------------------------
1883     ENDDO
1884     c-------------------------------------------------------
1885     c (2.3.2) Linear
1886     c-------------------------------------------------------
1887     DO J = 1,NyC
1888     c....MASK DEEP_F
1889     IF (J.EQ.NyC) THEN !EVITO ERRORI DI CHECK BOUND x Vittorio
1890     DEPDEP_F_WesternB (:) = DEEP_F(1 ,J,:)*DEEP_F(1 ,J,:)
1891     DEPDEP_F_EasternB(:) = DEEP_F(NxC,J,:)*DEEP_F(NxC,J,:)
1892     ELSE
1893     DEPDEP_F_WesternB (:) = DEEP_F(1 ,J,:)*DEEP_F(1 ,J+1,:)
1894     DEPDEP_F_EasternB(:) = DEEP_F(NxC,J,:)*DEEP_F(NxC,J+1,:)
1895     ENDIF
1896    
1897     c....WesternB...T.....................
1898     diff(:) = globalP3(1,P2C_linU(J)+1,:)-
1899     & globalP3(1,P2C_linU(J) ,:)
1900    
1901     T_F(J,:,1) = globalP3(1,P2C_linU(J) ,:)+
1902     & (diff(:)/3.)*float((WO3_linU(J)))
1903     & *DEPDEP_F_WesternB(:)
1904     c.....EasternB..T.....................
1905     diff(:) = globalP3(6,P2C_linU(J)+1,:)-
1906     & globalP3(6,P2C_linU(J) ,:)
1907    
1908     T_F(J,:,2) = globalP3(6,P2C_linU(J) ,:)+
1909     & (diff(:)/3.)*float((WO3_linU(J)))
1910     & *DEPDEP_F_EasternB(:)
1911     c.....WesternB...S.....................
1912     diff(:) = globalP4(1,P2C_linU(J)+1,:)-
1913     & globalP4 (1,P2C_linU(J) ,:)
1914    
1915     S_F(J,:,1) = globalP4(1,P2C_linU(J) ,:)+
1916     & (diff(:)/3.)*float((WO3_linU(J)))
1917     & *DEPDEP_F_WesternB(:)
1918     c.... EasternB..S.....................
1919     diff(:) = (globalP4(6,P2C_linU(J)+1,:)-
1920     & globalP4 (6,P2C_linU(J) ,:))
1921    
1922     S_F(J,:,2) = globalP4(6,P2C_linU(J) ,:)+
1923     & (diff(:)/3.)*float((WO3_linU(J)))
1924     & *DEPDEP_F_EasternB(:)
1925     c.... WesternB...Eta.....................
1926     diff(:) = globalP5(1,P2C_linU(J)+1,:)-
1927     & globalP5(1,P2C_linU(J) ,:)
1928    
1929     eta_F(J,:,1) = globalP5(1,P2C_linU(J) ,:)+
1930     & (diff(:)/3.)*float((WO3_linU(J)))
1931     & *DEPDEP_F_WesternB(:)
1932     c.... EasternB..Eta.....................
1933     diff(:) = globalP5(6,P2C_linU(J)+1,:)-
1934     & globalP5(6,P2C_linU(J) ,:)
1935    
1936     eta_F(J,:,2) = globalP5(6,P2C_linU(J) ,:)+
1937     & (diff(:)/3.)*float(WO3_linU(J))
1938     & *DEPDEP_F_EasternB(:)
1939    
1940     cgm--------------------------
1941    
1942     c.... WesternB...V.....................
1943     c diff(:) = globalP2(1,P2C_linU(J)+1,:)-
1944     c & globalP2(1,P2C_linU(J) ,:)
1945     c
1946     c V_F(J,:,1) = globalP2(1,P2C_linU(J) ,:)+
1947     c & (diff(:)/3.)*float((WO3_linU(J)))
1948     c & *DEPDEP_F_WesternB(:)
1949     c.... EasternB..V.....................
1950     c diff(:) = globalP2(6,P2C_linU(J)+1,:)-
1951     c & globalP2(6,P2C_linU(J) ,:)
1952     c
1953     c V_F(J,:,2) = globalP2(6,P2C_linU(J) ,:)+
1954     c & (diff(:)/3.)*float((WO3_linU(J)))
1955     c & *DEPDEP_F_EasternB(:)
1956     cgm----------------------------
1957     ENDDO
1958    
1959    
1960     8765 CONTINUE
1961     c=============== END INTERPOLAZIONI x CHILD ====================
1962    
1963     RETURN
1964     END
1965    
1966    
1967    
1968    
1969    
1970    
1971    
1972    
1973    
1974    
1975     c================================================================
1976     SUBROUTINE BLINT(x1,x2,x3,x4,y1,y2,y3,y4,f1,f2,f3,f4,x,y,f)
1977     c================================================================
1978     C
1979     C Bilinear interpolation subroutine.
1980     C (Xi,Yi,fi) = data grid & values surounding model point (x,y)
1981     C f = interpolated value at the model grid point.
1982     IMPLICIT NONE
1983     real*4 a1,a2,a3,a4
1984     real*4 b1,b2,b3,b4
1985     real*4 x1,x2,x3,x4
1986     real*4 y1,y2,y3,y4
1987     real*4 f1,f2,f3,f4
1988     real*4 x,y,A,B,C,t,f,s
1989     C
1990     a1=x1-x2+x3-x4
1991     a2=-x1+x4
1992     a3=-x1+x2
1993     a4=x1-x
1994     b1=y1-y2+y3-y4
1995     b2=-y1+y4
1996     b3=-y1+y2
1997     b4=y1-y
1998     A=a3*b1-a1*b3
1999     B=b2*a3+b1*a4-a1*b4-a2*b3
2000     C=-a2*b4+a4*b2
2001     if(ABS(A*C).gt.0.002*B**2) then
2002     t=(-B-sqrt(B*B-4.*A*C))/(2.*A)
2003     else
2004     t=C/ABS(B)
2005     endif
2006     10 CONTINUE
2007     A=a2*b1-a1*b2
2008     B=b3*a2+b1*a4-a1*b4-a3*b2
2009     C=-a3*b4+a4*b3
2010     if(ABS(A*C).gt.0.002*B**2) then
2011     s=(-B+sqrt(B*B-4.*A*C))/(2.*A)
2012     else
2013     s=-C/ABS(B)
2014     endif
2015     20 CONTINUE
2016     f=f1*(1.-t)*(1.-s)+f2*t*(1.-s)+f3*s*t+f4*(1.-t)*s
2017     return
2018     end

  ViewVC Help
Powered by ViewVC 1.1.22