source: utils/tools/NESTING/src/agrif_modutil.f90 @ 10027

Last change on this file since 10027 was 10027, checked in by clem, 2 years ago

solve issues with median averages (greatly increase speed and sort out land issues)

  • Property svn:keywords set to Id
File size: 12.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  !                    --- quicksort ---
53  ! Author: t-nissie
54  ! License: GPLv3
55  ! Gist: https://gist.github.com/t-nissie/479f0f16966925fa29ea
56  !***********************************************************
57  RECURSIVE SUBROUTINE quicksort(var, first, last)
58     IMPLICIT NONE
59
60     REAL*8, DIMENSION(:), INTENT(inout) :: var
61     INTEGER,              INTENT(in)    :: first, last
62     REAL*8  :: x, t
63     INTEGER :: ji, jj
64     
65     x = var( (first+last) / 2 )
66     ji = first
67     jj = last
68     DO
69        DO WHILE (var(ji) < x)
70           ji=ji+1
71        END DO
72        DO WHILE (x < var(jj))
73           jj=jj-1
74        END DO
75        IF (ji >= jj) EXIT
76        t = var(ji);  var(ji) = var(jj);  var(jj) = t
77        ji=ji+1
78        jj=jj-1
79     END DO
80     IF (first < ji-1) CALL quicksort(var, first, ji-1)
81     IF (jj+1 < last)  CALL quicksort(var, jj+1, last)
82  END SUBROUTINE quicksort
83 
84!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
85  !************************************************************************
86  !   SUBROUTINE 1Dto2D
87  !************************************************************************
88!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
89
90  SUBROUTINE tab1Dto2D(tab1D,tab2D,nx,ny)
91    !
92    IMPLICIT NONE   
93    !
94    REAL*8,DIMENSION(:) :: tab1D
95    REAL*8,DIMENSION(:,:) :: tab2D
96    !         
97    INTEGER :: xpos,ypos
98    INTEGER :: nx,ny   
99    INTEGER :: i     
100    !     
101    xpos=0
102    ypos=1
103    !     
104    DO i=1,nx*ny
105       xpos=xpos+1
106       IF(xpos.GT.nx) THEN
107          xpos=1
108          ypos=ypos+1
109       ENDIF
110       tab2D(ypos,xpos)=tab1D(i)
111    END DO
112    !
113  END SUBROUTINE tab1Dto2D
114
115!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
116  !************************************************************************
117  !   SUBROUTINE tab2Dto1D
118  !************************************************************************
119!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
120
121  SUBROUTINE tab2Dto1D(tab2D,tab1D)
122    !
123    IMPLICIT NONE 
124    !       
125    REAL*8,DIMENSION(:,:) :: tab2D
126    REAL*8,DIMENSION(:) :: tab1D
127    !         
128    INTEGER :: xpos,ypos
129    INTEGER :: nx,ny
130    INTEGER :: i
131    !
132    nx = SIZE(tab2D,2)
133    ny = SIZE(tab2D,1)
134    !     
135    xpos = 0
136    ypos = 1
137    DO i = 1,nx*ny
138       xpos = xpos + 1
139       IF(xpos.GT.nx) THEN
140          xpos = 1
141          ypos = ypos + 1
142       END IF
143       tab1D(i) = tab2D(ypos,xpos)
144    END DO
145    !
146  END SUBROUTINE tab2Dto1D
147
148
149!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
150  !************************************************************************
151  !   SUBROUTINE tab2Dto1D logical
152  !************************************************************************
153!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
154
155  SUBROUTINE logtab2Dto1D(tab2D,tab1D)
156    !
157    IMPLICIT NONE 
158    !       
159    LOGICAL,DIMENSION(:,:) :: tab2D
160    LOGICAL,DIMENSION(:) :: tab1D
161    !         
162    INTEGER :: xpos,ypos
163    INTEGER :: nx,ny
164    INTEGER :: i
165    !
166    nx = SIZE(tab2D,2)
167    ny = SIZE(tab2D,1)
168    !     
169    xpos = 0
170    ypos = 1
171    DO i = 1,nx*ny
172       xpos = xpos + 1
173       IF(xpos.GT.nx) THEN
174          xpos = 1
175          ypos = ypos + 1
176       END IF
177       tab1D(i) = tab2D(ypos,xpos)
178    END DO
179    !
180  END SUBROUTINE logtab2Dto1D
181  !
182!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
183  !************************************************************************
184  !   SUBROUTINE 1Dto2D
185  !************************************************************************
186!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
187
188  SUBROUTINE logtab1Dto2D(tab1D,tab2D,nx,ny)
189    !
190    IMPLICIT NONE   
191    !
192    LOGICAL,DIMENSION(:) :: tab1D
193    LOGICAL,DIMENSION(:,:) :: tab2D
194    !         
195    INTEGER :: xpos,ypos
196    INTEGER :: nx,ny   
197    INTEGER :: i     
198    !     
199    xpos=0
200    ypos=1
201    !     
202    DO i=1,nx*ny
203       xpos=xpos+1
204       IF(xpos.GT.nx) THEN
205          xpos=1
206          ypos=ypos+1
207       ENDIF
208       tab2D(ypos,xpos)=tab1D(i)
209    END DO
210    !
211  END SUBROUTINE logtab1Dto2D
212
213  !**************************************************************
214  !   subroutine make_remap
215  !**************************************************************           
216  !
217  SUBROUTINE make_remap(tabin,tabout,nxfin,nyfin,matrix,src_add,dst_add) 
218    !     
219    IMPLICIT NONE
220    !     
221    REAL*8, DIMENSION(:,:)  :: tabin 
222    REAL*8, DIMENSION(:,:) :: tabout 
223    REAL*8, POINTER, DIMENSION(:,:) :: tabtemp   
224    INTEGER,DIMENSION(:) :: src_add,dst_add
225    INTEGER :: nxfin,nyfin     
226    REAL*8, POINTER, DIMENSION(:) :: var1D,var_interp1D
227    REAL*8,DIMENSION(:,:) :: matrix 
228    INTEGER :: num_links,i
229    !
230    ALLOCATE(var1D(SIZE(tabin,1)*SIZE(tabin,2)))     
231    CALL tab2Dto1D(tabin,var1D)
232    !     
233    ALLOCATE(var_interp1D(nxfin*nyfin))
234    var_interp1D = 0.0
235    num_links = SIZE(dst_add)
236    !
237    DO i = 1,num_links
238       var_interp1D(dst_add(i)) = var_interp1D(dst_add(i)) &
239            + matrix(1,i)*var1D(src_add(i))
240    END DO
241    !
242    ALLOCATE(tabtemp(SIZE(tabout,1),SIZE(tabout,2)))
243    !
244    CALL tab1Dto2D(var_interp1D,tabtemp,nyfin,nxfin)
245    !
246    tabout = tabtemp
247    !
248    DEALLOCATE(var_interp1D,var1D,tabtemp) 
249    !
250  END SUBROUTINE make_remap
251  !         
252  !
253  !**************************************************************
254  !   end subroutine make_remap
255  !**************************************************************         
256  !
257  !
258  !**************************************************************
259  !   subroutine make_bicubic_remap
260  !**************************************************************           
261  !
262  SUBROUTINE make_bicubic_remap(tabin,masksrc,tabout,nxfin,nyfin,matrix,src_add,dst_add) 
263    !     
264    IMPLICIT NONE
265    !     
266    REAL*8, DIMENSION(:,:)  :: tabin
267    LOGICAL, DIMENSION(:,:)  :: masksrc
268    LOGICAL, POINTER, DIMENSION(:)  :: grid1_mask
269    REAL*8, DIMENSION(:,:) :: tabout     
270    INTEGER,DIMENSION(:) :: src_add,dst_add
271    INTEGER :: nxfin,nyfin     
272    REAL*8, POINTER, DIMENSION(:) :: var1D,var_interp1D,gradi,gradj,gradij,deriv1,deriv2
273    REAL*8,DIMENSION(:,:) :: matrix 
274    INTEGER :: num_links,i,j,nx,ny,n,ip1,im1,jp1,jm1
275    INTEGER :: in,is,ie,iw,ine,inw,ise,isw
276    REAL*8 :: delew,delns
277    !
278    nx = SIZE(tabin,1)
279    ny = SIZE(tabin,2)
280    ALLOCATE(gradi(nx*ny),gradj(nx*ny),gradij(nx*ny),deriv1(nx*ny),deriv2(nx*ny))
281    ALLOCATE(var1D(nx*ny),grid1_mask(nx*ny))
282    !           
283    CALL tab2Dto1D(tabin,var1D)
284    CALL logtab2Dto1D(masksrc,grid1_mask)     
285    !
286    gradi  = 0.0
287    gradj  = 0.0
288    gradij = 0.0
289    !
290    DO n = 1,nx*ny
291
292       IF( grid1_mask(n) ) THEN                       
293          !                     
294          delew = 0.5     
295          delns = 0.5
296
297          j = (n-1)/ny + 1
298          i = n - (j-1)*ny
299          !
300          ip1 = i+1
301          im1 = i-1
302          jp1 = j+1
303          jm1 = j-1     
304          !
305          IF (ip1 > ny) ip1 = ip1 - ny
306
307          IF (im1 < 1 ) im1 = ny
308
309          IF (jp1 > nx) THEN
310             jp1 = j
311             delns = 1.
312          ENDIF
313
314          IF (jm1 < 1 ) THEN
315             jm1 = j
316             delns = 1.
317          ENDIF
318          !
319          in  = (jp1-1)*ny + i
320          is  = (jm1-1)*ny + i
321          ie  = (j  -1)*ny + ip1
322          iw  = (j  -1)*ny + im1
323          !
324          ine = (jp1-1)*ny + ip1
325          inw = (jp1-1)*ny + im1
326          ise = (jm1-1)*ny + ip1
327          isw = (jm1-1)*ny + im1
328          !
329          !*** compute i-gradient
330
331          IF (.NOT. grid1_mask(ie)) THEN
332             ie = n
333             delew = 1.
334          ENDIF
335          !           
336          IF (.NOT. grid1_mask(iw)) THEN
337             iw = n
338             delew = 1.
339          ENDIF
340          !
341          gradi(n) = delew*(var1D(ie) - var1D(iw))
342          !
343          !*** compute j-gradient
344
345          IF (.NOT. grid1_mask(in)) THEN
346             in = n
347             delns = 1.
348          ENDIF
349          !                     
350          IF (.NOT. grid1_mask(is)) THEN
351             is = n
352             delns = 1.
353          ENDIF
354          !
355          gradj(n) = delns*(var1D(in) - var1D(is))                   
356          !
357          !*** compute ij-gradient
358
359          delew = 0.5
360
361          IF (jp1 == j .OR. jm1 == j) THEN
362             delns = 1.
363          ELSE
364             delns = 0.5
365          ENDIF
366          !
367          IF (.NOT. grid1_mask(ine)) THEN
368             IF (in /= n) THEN
369                ine = in
370                delew = 1.
371             ELSE IF (ie /= n) THEN
372                ine = ie
373                inw = iw
374                IF (inw == n) delew = 1.
375                delns = 1.
376             ELSE
377                ine = n
378                inw = iw
379                delew = 1
380                delns = 1
381             ENDIF
382          ENDIF
383          !
384          IF (.NOT. grid1_mask(inw)) THEN
385             IF (in /= n) THEN
386                inw = in
387                delew = 1.
388             ELSE IF (iw /= n) THEN
389                inw = iw
390                ine = ie
391                IF (ie == n) delew = 1.
392                delns = 1.
393             ELSE
394                inw = n
395                ine = ie
396                delew = 1.
397                delns = 1.
398             ENDIF
399          ENDIF
400          !
401          deriv1(n) = delew*(var1D(ine)-var1D(inw))                   
402          !
403          IF (.NOT. grid1_mask(ise)) THEN
404             IF (is /= n) THEN
405                ise = is
406                delew = 1.
407             ELSE IF (ie /= n) THEN
408                ise = ie
409                isw = iw
410                IF (isw == n) delew = 1.
411                delns = 1.
412             ELSE
413                ise = n
414                isw = iw
415                delew = 1.
416                delns = 1.
417             ENDIF
418          ENDIF
419          !
420          IF (.NOT. grid1_mask(isw)) THEN
421             IF (is /= n) THEN
422                isw = is
423                delew = 1.
424             ELSE IF (iw /= n) THEN
425                isw = iw
426                ise = ie
427                IF (ie == n) delew = 1.
428                delns = 1.
429             ELSE
430                isw = n
431                ise = ie
432                delew = 1.
433                delns = 1.
434             ENDIF
435          ENDIF
436
437          deriv2(n) = delew*(var1D(ise) - var1D(isw))
438          gradij(n) = delns*(deriv1(n) - deriv2(n))
439       ENDIF
440    END DO
441    !
442    DEALLOCATE(deriv1,deriv2,grid1_mask)
443
444    !       
445    ALLOCATE(var_interp1D(nxfin*nyfin))
446    var_interp1D = 0.0
447    num_links = SIZE(dst_add)
448    !     
449    DO i = 1,num_links
450       !       
451       var_interp1D(dst_add(i)) = var_interp1D(dst_add(i))        +      &
452            matrix(1,i)*var1D(src_add(i))   +      &
453            matrix(2,i)*gradi(src_add(i))   +      &
454            matrix(3,i)*gradj(src_add(i))   +      & 
455            matrix(4,i)*gradij(src_add(i))
456    END DO
457    !
458    DEALLOCATE(gradi,gradj,gradij,var1D)
459    !
460    CALL tab1Dto2D(var_interp1D,tabout,nyfin,nxfin)
461    !     
462    DEALLOCATE(var_interp1D)       
463    !
464  END SUBROUTINE make_bicubic_remap
465  !         
466  !
467  !**************************************************************
468  !   end subroutine make_bicubic_remap
469  !**************************************************************         
470  !
471  !
472END MODULE agrif_modutil
Note: See TracBrowser for help on using the repository browser.