New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
agrif_modutil.f90 in utils/tools/DOMAINcfg/src – NEMO

source: utils/tools/DOMAINcfg/src/agrif_modutil.f90 @ 13204

Last change on this file since 13204 was 13204, checked in by smasson, 4 years ago

tools: update with tools_dev_r12970_AGRIF_CMEMS

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