source: branches/UKMO/r6232_tracer_advection/NEMOGCM/TOOLS/NESTING/src/agrif_modutil.f90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

File size: 11.4 KB
Line 
1!
2MODULE agrif_modutil
3  !
4CONTAINS
5  !
6  !************************************************************************
7  !                           *
8  ! MODULE  AGRIF_MODUTIL                    *
9  !                           *
10  ! module containing subroutine used for :           *
11  !   - unrolling 2D arrays to 1D arrays (required for SCRIP package use)  *
12  !   - convert 1D arrays to 2D arrays (required for SCRIP package use) *
13  !   - remapping process (use SCRIP remapping matrix)         *
14  !                           *
15  !************************************************************************
16  !
17  !***********************************************************   
18  SUBROUTINE ssort (x, nb)
19    !***********************************************************
20    !
21    IMPLICIT NONE
22
23    INTEGER :: nb
24    REAL*8, DIMENSION(:) :: x
25    REAL*8 :: temp
26    INTEGER ji,jj,jmax,itemp
27    !     
28    jmax=nb-1
29    !
30    !
31    DO ji=1,nb-1
32       temp=HUGE(1)
33
34       DO jj=1,jmax
35
36          IF(X(jj).LE.X(jj+1)) THEN
37             temp=X(jj)
38             X(jj)=X(jj+1)
39             X(jj+1)=temp
40          ENDIF
41
42       ENDDO
43
44       IF(temp.EQ.HUGE(1)) RETURN
45       jmax=jmax-1
46    ENDDO
47
48    RETURN
49  END SUBROUTINE ssort
50  !
51!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
52  !************************************************************************
53  !   SUBROUTINE 1Dto2D
54  !************************************************************************
55!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
56
57  SUBROUTINE tab1Dto2D(tab1D,tab2D,nx,ny)
58    !
59    IMPLICIT NONE   
60    !
61    REAL*8,DIMENSION(:) :: tab1D
62    REAL*8,DIMENSION(:,:) :: tab2D
63    !         
64    INTEGER :: xpos,ypos
65    INTEGER :: nx,ny   
66    INTEGER :: i     
67    !     
68    xpos=0
69    ypos=1
70    !     
71    DO i=1,nx*ny
72       xpos=xpos+1
73       IF(xpos.GT.nx) THEN
74          xpos=1
75          ypos=ypos+1
76       ENDIF
77       tab2D(ypos,xpos)=tab1D(i)
78    END DO
79    !
80  END SUBROUTINE tab1Dto2D
81
82!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
83  !************************************************************************
84  !   SUBROUTINE tab2Dto1D
85  !************************************************************************
86!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
87
88  SUBROUTINE tab2Dto1D(tab2D,tab1D)
89    !
90    IMPLICIT NONE 
91    !       
92    REAL*8,DIMENSION(:,:) :: tab2D
93    REAL*8,DIMENSION(:) :: tab1D
94    !         
95    INTEGER :: xpos,ypos
96    INTEGER :: nx,ny
97    INTEGER :: i
98    !
99    nx = SIZE(tab2D,2)
100    ny = SIZE(tab2D,1)
101    !     
102    xpos = 0
103    ypos = 1
104    DO i = 1,nx*ny
105       xpos = xpos + 1
106       IF(xpos.GT.nx) THEN
107          xpos = 1
108          ypos = ypos + 1
109       END IF
110       tab1D(i) = tab2D(ypos,xpos)
111    END DO
112    !
113  END SUBROUTINE tab2Dto1D
114
115
116!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
117  !************************************************************************
118  !   SUBROUTINE tab2Dto1D logical
119  !************************************************************************
120!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
121
122  SUBROUTINE logtab2Dto1D(tab2D,tab1D)
123    !
124    IMPLICIT NONE 
125    !       
126    LOGICAL,DIMENSION(:,:) :: tab2D
127    LOGICAL,DIMENSION(:) :: tab1D
128    !         
129    INTEGER :: xpos,ypos
130    INTEGER :: nx,ny
131    INTEGER :: i
132    !
133    nx = SIZE(tab2D,2)
134    ny = SIZE(tab2D,1)
135    !     
136    xpos = 0
137    ypos = 1
138    DO i = 1,nx*ny
139       xpos = xpos + 1
140       IF(xpos.GT.nx) THEN
141          xpos = 1
142          ypos = ypos + 1
143       END IF
144       tab1D(i) = tab2D(ypos,xpos)
145    END DO
146    !
147  END SUBROUTINE logtab2Dto1D
148  !
149!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
150  !************************************************************************
151  !   SUBROUTINE 1Dto2D
152  !************************************************************************
153!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
154
155  SUBROUTINE logtab1Dto2D(tab1D,tab2D,nx,ny)
156    !
157    IMPLICIT NONE   
158    !
159    LOGICAL,DIMENSION(:) :: tab1D
160    LOGICAL,DIMENSION(:,:) :: tab2D
161    !         
162    INTEGER :: xpos,ypos
163    INTEGER :: nx,ny   
164    INTEGER :: i     
165    !     
166    xpos=0
167    ypos=1
168    !     
169    DO i=1,nx*ny
170       xpos=xpos+1
171       IF(xpos.GT.nx) THEN
172          xpos=1
173          ypos=ypos+1
174       ENDIF
175       tab2D(ypos,xpos)=tab1D(i)
176    END DO
177    !
178  END SUBROUTINE logtab1Dto2D
179
180  !**************************************************************
181  !   subroutine make_remap
182  !**************************************************************           
183  !
184  SUBROUTINE make_remap(tabin,tabout,nxfin,nyfin,matrix,src_add,dst_add) 
185    !     
186    IMPLICIT NONE
187    !     
188    REAL*8, DIMENSION(:,:)  :: tabin 
189    REAL*8, DIMENSION(:,:) :: tabout 
190    REAL*8, POINTER, DIMENSION(:,:) :: tabtemp   
191    INTEGER,DIMENSION(:) :: src_add,dst_add
192    INTEGER :: nxfin,nyfin     
193    REAL*8, POINTER, DIMENSION(:) :: var1D,var_interp1D
194    REAL*8,DIMENSION(:,:) :: matrix 
195    INTEGER :: num_links,i
196    !
197    ALLOCATE(var1D(SIZE(tabin,1)*SIZE(tabin,2)))     
198    CALL tab2Dto1D(tabin,var1D)
199    !     
200    ALLOCATE(var_interp1D(nxfin*nyfin))
201    var_interp1D = 0.0
202    num_links = SIZE(dst_add)
203    !
204    DO i = 1,num_links
205       var_interp1D(dst_add(i)) = var_interp1D(dst_add(i)) &
206            + matrix(1,i)*var1D(src_add(i))
207    END DO
208    !
209    ALLOCATE(tabtemp(SIZE(tabout,1),SIZE(tabout,2)))
210    !
211    CALL tab1Dto2D(var_interp1D,tabtemp,nyfin,nxfin)
212    !
213    tabout = tabtemp
214    !
215    DEALLOCATE(var_interp1D,var1D,tabtemp) 
216    !
217  END SUBROUTINE make_remap
218  !         
219  !
220  !**************************************************************
221  !   end subroutine make_remap
222  !**************************************************************         
223  !
224  !
225  !**************************************************************
226  !   subroutine make_bicubic_remap
227  !**************************************************************           
228  !
229  SUBROUTINE make_bicubic_remap(tabin,masksrc,tabout,nxfin,nyfin,matrix,src_add,dst_add) 
230    !     
231    IMPLICIT NONE
232    !     
233    REAL*8, DIMENSION(:,:)  :: tabin
234    LOGICAL, DIMENSION(:,:)  :: masksrc
235    LOGICAL, POINTER, DIMENSION(:)  :: grid1_mask
236    REAL*8, DIMENSION(:,:) :: tabout     
237    INTEGER,DIMENSION(:) :: src_add,dst_add
238    INTEGER :: nxfin,nyfin     
239    REAL*8, POINTER, DIMENSION(:) :: var1D,var_interp1D,gradi,gradj,gradij,deriv1,deriv2
240    REAL*8,DIMENSION(:,:) :: matrix 
241    INTEGER :: num_links,i,j,nx,ny,n,ip1,im1,jp1,jm1
242    INTEGER :: in,is,ie,iw,ine,inw,ise,isw
243    REAL*8 :: delew,delns
244    !
245    nx = SIZE(tabin,1)
246    ny = SIZE(tabin,2)
247    ALLOCATE(gradi(nx*ny),gradj(nx*ny),gradij(nx*ny),deriv1(nx*ny),deriv2(nx*ny))
248    ALLOCATE(var1D(nx*ny),grid1_mask(nx*ny))
249    !           
250    CALL tab2Dto1D(tabin,var1D)
251    CALL logtab2Dto1D(masksrc,grid1_mask)     
252    !
253    gradi  = 0.0
254    gradj  = 0.0
255    gradij = 0.0
256    !
257    DO n = 1,nx*ny
258
259       IF( grid1_mask(n) ) THEN                       
260          !                     
261          delew = 0.5     
262          delns = 0.5
263
264          j = (n-1)/ny + 1
265          i = n - (j-1)*ny
266          !
267          ip1 = i+1
268          im1 = i-1
269          jp1 = j+1
270          jm1 = j-1     
271          !
272          IF (ip1 > ny) ip1 = ip1 - ny
273
274          IF (im1 < 1 ) im1 = ny
275
276          IF (jp1 > nx) THEN
277             jp1 = j
278             delns = 1.
279          ENDIF
280
281          IF (jm1 < 1 ) THEN
282             jm1 = j
283             delns = 1.
284          ENDIF
285          !
286          in  = (jp1-1)*ny + i
287          is  = (jm1-1)*ny + i
288          ie  = (j  -1)*ny + ip1
289          iw  = (j  -1)*ny + im1
290          !
291          ine = (jp1-1)*ny + ip1
292          inw = (jp1-1)*ny + im1
293          ise = (jm1-1)*ny + ip1
294          isw = (jm1-1)*ny + im1
295          !
296          !*** compute i-gradient
297
298          IF (.NOT. grid1_mask(ie)) THEN
299             ie = n
300             delew = 1.
301          ENDIF
302          !           
303          IF (.NOT. grid1_mask(iw)) THEN
304             iw = n
305             delew = 1.
306          ENDIF
307          !
308          gradi(n) = delew*(var1D(ie) - var1D(iw))
309          !
310          !*** compute j-gradient
311
312          IF (.NOT. grid1_mask(in)) THEN
313             in = n
314             delns = 1.
315          ENDIF
316          !                     
317          IF (.NOT. grid1_mask(is)) THEN
318             is = n
319             delns = 1.
320          ENDIF
321          !
322          gradj(n) = delns*(var1D(in) - var1D(is))                   
323          !
324          !*** compute ij-gradient
325
326          delew = 0.5
327
328          IF (jp1 == j .OR. jm1 == j) THEN
329             delns = 1.
330          ELSE
331             delns = 0.5
332          ENDIF
333          !
334          IF (.NOT. grid1_mask(ine)) THEN
335             IF (in /= n) THEN
336                ine = in
337                delew = 1.
338             ELSE IF (ie /= n) THEN
339                ine = ie
340                inw = iw
341                IF (inw == n) delew = 1.
342                delns = 1.
343             ELSE
344                ine = n
345                inw = iw
346                delew = 1
347                delns = 1
348             ENDIF
349          ENDIF
350          !
351          IF (.NOT. grid1_mask(inw)) THEN
352             IF (in /= n) THEN
353                inw = in
354                delew = 1.
355             ELSE IF (iw /= n) THEN
356                inw = iw
357                ine = ie
358                IF (ie == n) delew = 1.
359                delns = 1.
360             ELSE
361                inw = n
362                ine = ie
363                delew = 1.
364                delns = 1.
365             ENDIF
366          ENDIF
367          !
368          deriv1(n) = delew*(var1D(ine)-var1D(inw))                   
369          !
370          IF (.NOT. grid1_mask(ise)) THEN
371             IF (is /= n) THEN
372                ise = is
373                delew = 1.
374             ELSE IF (ie /= n) THEN
375                ise = ie
376                isw = iw
377                IF (isw == n) delew = 1.
378                delns = 1.
379             ELSE
380                ise = n
381                isw = iw
382                delew = 1.
383                delns = 1.
384             ENDIF
385          ENDIF
386          !
387          IF (.NOT. grid1_mask(isw)) THEN
388             IF (is /= n) THEN
389                isw = is
390                delew = 1.
391             ELSE IF (iw /= n) THEN
392                isw = iw
393                ise = ie
394                IF (ie == n) delew = 1.
395                delns = 1.
396             ELSE
397                isw = n
398                ise = ie
399                delew = 1.
400                delns = 1.
401             ENDIF
402          ENDIF
403
404          deriv2(n) = delew*(var1D(ise) - var1D(isw))
405          gradij(n) = delns*(deriv1(n) - deriv2(n))
406       ENDIF
407    END DO
408    !
409    DEALLOCATE(deriv1,deriv2,grid1_mask)
410
411    !       
412    ALLOCATE(var_interp1D(nxfin*nyfin))
413    var_interp1D = 0.0
414    num_links = SIZE(dst_add)
415    !     
416    DO i = 1,num_links
417       !       
418       var_interp1D(dst_add(i)) = var_interp1D(dst_add(i))        +      &
419            matrix(1,i)*var1D(src_add(i))   +      &
420            matrix(2,i)*gradi(src_add(i))   +      &
421            matrix(3,i)*gradj(src_add(i))   +      & 
422            matrix(4,i)*gradij(src_add(i))
423    END DO
424    !
425    DEALLOCATE(gradi,gradj,gradij,var1D)
426    !
427    CALL tab1Dto2D(var_interp1D,tabout,nyfin,nxfin)
428    !     
429    DEALLOCATE(var_interp1D)       
430    !
431  END SUBROUTINE make_bicubic_remap
432  !         
433  !
434  !**************************************************************
435  !   end subroutine make_bicubic_remap
436  !**************************************************************         
437  !
438  !
439END MODULE agrif_modutil
Note: See TracBrowser for help on using the repository browser.