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.
modmask.F in trunk/AGRIF/AGRIF_FILES – NEMO

source: trunk/AGRIF/AGRIF_FILES/modmask.F @ 779

Last change on this file since 779 was 779, checked in by rblod, 16 years ago

Agrif improvment for vectorization, see ticket #41

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 21.8 KB
RevLine 
[396]1!
2! $Id$
3!
4C     AGRIF (Adaptive Grid Refinement In Fortran)
5C
6C     Copyright (C) 2003 Laurent Debreu (Laurent.Debreu@imag.fr)
7C                        Christophe Vouland (Christophe.Vouland@imag.fr)   
8C
9C     This program is free software; you can redistribute it and/or modify
10C     it under the terms of the GNU General Public License as published by
11C     the Free Software Foundation; either version 2 of the License, or
12C     (at your option) any later version.
13C
14C     This program is distributed in the hope that it will be useful,
15C     but WITHOUT ANY WARRANTY; without even the implied warranty of
16C     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17C     GNU General Public License for more details.
18C
19C     You should have received a copy of the GNU General Public License
20C     along with this program; if not, write to the Free Software
21C     Foundation, Inc., 59 Temple Place -  Suite 330, Boston, MA 02111-1307, USA.
22C
23C
24C
25CCC   Module Agrif_Mask
26C
27      Module Agrif_Mask
28C
29CCC   Description:
30CCC   Module for masks
31C
32C     Modules used: 
33C
34      Use Agrif_Types       
35C
36      IMPLICIT NONE
[779]37      Integer, Parameter :: MaxSearch = 3
[396]38C
39      CONTAINS
40C     Define procedures contained in this module
41C
42C     **************************************************************************
43C     Subroutine Agrif_CheckMasknD
44C     **************************************************************************
45C       
46      Subroutine Agrif_CheckMasknD(tempP,parent,pbtab,petab,ppbtab,
47     &               ppetab,noraftab,nbdim)
48C
49CCC   Description:
50CCC   Subroutine called in the procedure Agrif_InterpnD to recalculate the value
51CCC   of the parent grid variable when this one is equal to Agrif_SpecialValue. 
52C
53C     Declarations:
54C
55       
56C
57C     Arrays arguments     
58      INTEGER :: nbdim
59      INTEGER,DIMENSION(nbdim) :: pbtab  ! Limits of the parent grid used 
60      INTEGER,DIMENSION(nbdim) :: petab  ! interpolation of the child grid
61      LOGICAL,DIMENSION(nbdim) :: noraftab
62      INTEGER,DIMENSION(nbdim) :: ppbtab,ppetab
63C
64C     Pointer argument
65      TYPE(AGRIF_PVARIABLE) :: tempP  ! Part of the parent grid used for
66                                      ! the interpolation of the child grid                     
67C     Data TYPE argument                                   
68      TYPE(AGRIF_PVARIABLE) :: parent      ! The parent grid
69C
70C     Local scalar
71      INTEGER                   :: i0,j0,k0,l0,m0,n0
72C     
73C     Local arrays
74C
75C     
76      SELECT CASE (nbdim)
77      CASE (1)
78         do i0 = pbtab(1),petab(1)
79         IF (tempP%var%array1(i0)
80     &                        == Agrif_SpecialValue) Then
81            Call CalculNewValTempP((/i0/),
82     &                             tempP,parent,
83     &                             ppbtab,ppetab,
84     &                             noraftab,nbdim)
85         ENDIF
86         enddo
87      CASE (2)
88         do j0 = pbtab(2),petab(2)
89         do i0 = pbtab(1),petab(1)
90         IF (tempP%var%array2(i0,j0)
91     &                        == Agrif_SpecialValue) Then
92            Call CalculNewValTempP((/i0,j0/),
93     &                             tempP,parent,
94     &                             ppbtab,ppetab,
95     &                             noraftab,nbdim)
96         ENDIF
97         enddo 
98         enddo
99      CASE (3)
100         do k0 = pbtab(3),petab(3)
101         do j0 = pbtab(2),petab(2)
102         do i0 = pbtab(1),petab(1)
103         IF (tempP%var%array3(i0,j0,k0)
104     &                        == Agrif_SpecialValue) Then
[779]105!------CDIR NEXPAND
106            Call CalculNewValTempP3D((/i0,j0,k0/),
107     &      tempP%var%array3(ppbtab(1),ppbtab(2),ppbtab(3)),
108     &      parent%var%array3(ppbtab(1),ppbtab(2),ppbtab(3)),
109     &                           ppbtab,ppetab,
110     &                           noraftab,MaxSearch,Agrif_SpecialValue)
111     
112c            Call CalculNewValTempP((/i0,j0,k0/),
113c     &                             tempP,parent,
114c     &                             ppbtab,ppetab,
115c     &                             noraftab,nbdim)
116         
[396]117         ENDIF
118         enddo
119         enddo 
120         enddo
121      CASE (4)
122         do l0 = pbtab(4),petab(4)
123         do k0 = pbtab(3),petab(3)
124         do j0 = pbtab(2),petab(2)
125         do i0 = pbtab(1),petab(1)
126         IF (tempP%var%array4(i0,j0,k0,l0) 
127     &                        == Agrif_SpecialValue) Then
[779]128            Call CalculNewValTempP4D((/i0,j0,k0,l0/),
129     &      tempP%var%array4(ppbtab(1),ppbtab(2),ppbtab(3),ppbtab(4)),
130     &      parent%var%array4(ppbtab(1),ppbtab(2),ppbtab(3),ppbtab(4)),
131     &                           ppbtab,ppetab,
132     &                           noraftab,MaxSearch,Agrif_SpecialValue)
[396]133         ENDIF
134         enddo
135         enddo
136         enddo 
137         enddo
138      CASE (5)
139         do m0 = pbtab(5),petab(5)
140         do l0 = pbtab(4),petab(4)
141         do k0 = pbtab(3),petab(3)
142         do j0 = pbtab(2),petab(2)
143         do i0 = pbtab(1),petab(1)
144         IF (tempP%var%array5(i0,j0,k0,l0,m0) 
145     &                       == Agrif_SpecialValue) Then
146            Call CalculNewValTempP((/i0,j0,k0,l0,m0/),
147     &                             tempP,parent,
148     &                             ppbtab,ppetab,
149     &                             noraftab,nbdim)
150         ENDIF
151         enddo
152         enddo
153         enddo 
154         enddo
155         enddo
156      CASE (6)
157         do n0 = pbtab(6),petab(6)
158         do m0 = pbtab(5),petab(5)
159         do l0 = pbtab(4),petab(4)
160         do k0 = pbtab(3),petab(3)
161         do j0 = pbtab(2),petab(2)
162         do i0 = pbtab(1),petab(1)
163         IF (tempP%var%array6(i0,j0,k0,l0,m0,n0) 
164     &                       == Agrif_SpecialValue) Then
165            Call CalculNewValTempP((/i0,j0,k0,l0,m0,n0/),
166     &                             tempP,parent,
167     &                             ppbtab,ppetab,
168     &                             noraftab,nbdim)
169         ENDIF
170         enddo
171         enddo
172         enddo
173         enddo 
174         enddo
175         enddo
176      END SELECT
177C
178C     
179C     
180      End Subroutine Agrif_CheckMasknD
181C
182C
183C     **************************************************************************
184C     Subroutine CalculNewValTempP
185C     **************************************************************************
186C       
187      Subroutine CalculNewValTempP(indic,
188     &               tempP,parent,ppbtab,
189     &               ppetab,noraftab,nbdim)
190C
191CCC   Description:
192CCC   Subroutine called in the procedure Agrif_InterpnD to recalculate the value
193CCC   of the parent grid variable when this one is equal to Agrif_SpecialValue. 
194C
195C     Declarations:
196C
197       
198C
199C     Arrays arguments     
200      INTEGER :: nbdim
201      INTEGER,DIMENSION(nbdim) :: indic
202      LOGICAL,DIMENSION(nbdim) :: noraftab
203      INTEGER,DIMENSION(nbdim) :: ppbtab,ppetab
204C
205C     Pointer argument
206      TYPE(AGRIF_PVARIABLE) :: tempP  ! Part of the parent grid used for
207                                      ! the interpolation of the child grid                     
208C     Data TYPE argument                                   
209      TYPE(AGRIF_PVARIABLE) :: parent      ! The parent grid
210C
211C     Local scalar
212      INTEGER                  :: i,ii,iii,jj,kk,ll,mm,nn 
[662]213      INTEGER,DIMENSION(nbdim) :: imin,imax,idecal
[396]214      INTEGER                  :: Nbvals
215      REAL                     :: Res
216      REAL                     :: ValParent
217      INTEGER                  :: ValMax
218      LOGICAL                  :: firsttest
219C     
220C     Local arrays     
221C
222      ValMax = 1
223      do iii = 1 , nbdim
224         IF (.NOT.noraftab(iii)) THEN
225         ValMax = max(ValMax,ppetab(iii)-indic(iii))
226         ValMax = max(ValMax,indic(iii)-ppbtab(iii))
227         ENDIF
228      enddo
229C
230      Valmax = min(Valmax,MaxSearch)
231C
[779]232!CDIR NOVECTOR
[396]233      imin = indic
[779]234!CDIR NOVECTOR
[396]235      imax = indic
236C
237         i = 1
238         firsttest = .TRUE.
239C
240         do While(i <= ValMax)
241C
242         IF ((i == 1).AND.(firsttest)) i = Valmax
243
244            do iii = 1 , nbdim
245               if (.NOT.noraftab(iii)) then
246                  imin(iii) = max(indic(iii) - i,ppbtab(iii))
247                  imax(iii) = min(indic(iii) + i,ppetab(iii))
[662]248                  if (firsttest) then               
249                  if (indic(iii).GT.ppbtab(iii)) then
250
[779]251!CDIR NOVECTOR
[662]252                     idecal = indic
253                     idecal(iii) = idecal(iii)-1
254                     SELECT CASE(nbdim)
255                     CASE (1)
256                        if (tempP%var%array1(idecal(1)
257     &                            ) == Agrif_SpecialValue) then
258                           imin(iii) = imax(iii)
259                        endif                     
260                     CASE (2)
261                        if (tempP%var%array2(idecal(1),
262     &            idecal(2)) == Agrif_SpecialValue) then
263                           imin(iii) = imax(iii)
264                        endif
265                     CASE (3)
266                        if (tempP%var%array3(idecal(1),
267     &            idecal(2),idecal(3)) 
268     &                               == Agrif_SpecialValue) then
269                           imin(iii) = imax(iii)
270                        endif 
271                     CASE (4)
272                        if (tempP%var%array4(idecal(1),
273     &            idecal(2),idecal(3),idecal(4)) 
274     &                               == Agrif_SpecialValue) then
275                           imin(iii) = imax(iii)
276                        endif     
277                     CASE (5)
278                        if (tempP%var%array5(idecal(1),
279     &            idecal(2),idecal(3),idecal(4),idecal(5)) 
280     &                               == Agrif_SpecialValue) then
281                           imin(iii) = imax(iii)
282                        endif 
283                     CASE (6)
284                        if (tempP%var%array6(idecal(1),
285     &            idecal(2),idecal(3),idecal(4),idecal(5),idecal(6)) 
286     &                               == Agrif_SpecialValue) then
287                           imin(iii) = imax(iii)
288                        endif                                                                                           
289                     END SELECT
290                  endif
291                  endif
[396]292               endif           
293            enddo
294C
295            Res = 0.
296            Nbvals = 0
297C
[662]298            SELECT CASE(nbdim)
299            CASE (1)
[779]300!CDIR ALTCODE
301!CDIR SHORTLOOP
[396]302               do ii = imin(1),imax(1)
303                    ValParent = parent%var%array1(ii)
304                    if ( ValParent .NE. Agrif_SpecialValue) then
305                        Res = Res + ValParent
306                        Nbvals = Nbvals + 1
307                    endif
308               enddo
309C
[662]310            CASE (2)
[396]311               do jj = imin(2),imax(2)
[779]312!CDIR ALTCODE
313!CDIR SHORTLOOP
[396]314               do ii = imin(1),imax(1)
315                    ValParent = parent%var%array2(ii,jj)
316                    if ( ValParent .NE. Agrif_SpecialValue) then
317                        Res = Res + ValParent
318                        Nbvals = Nbvals + 1
319                    endif
320               enddo 
321               enddo
[662]322               
323            CASE (3)
[396]324               do kk = imin(3),imax(3)
325               do jj = imin(2),imax(2)
[779]326!CDIR ALTCODE
327!CDIR SHORTLOOP
[396]328               do ii = imin(1),imax(1)
329                    ValParent = parent%var%array3(ii,jj,kk)
330                    if ( ValParent .NE. Agrif_SpecialValue) then
331                        Res = Res + ValParent
332                        Nbvals = Nbvals + 1
333                    endif
334                        enddo
335                  enddo 
336               enddo
[662]337
338            CASE (4)
[396]339               do ll = imin(4),imax(4)
340               do kk = imin(3),imax(3)
341               do jj = imin(2),imax(2)
[779]342!CDIR ALTCODE
343!CDIR SHORTLOOP
[396]344               do ii = imin(1),imax(1)
345                    ValParent = parent%var%array4(ii,jj,kk,ll)
346                    if ( ValParent .NE. Agrif_SpecialValue) then
347                        Res = Res + ValParent
348                        Nbvals = Nbvals + 1
349                    endif
350                              enddo
351                        enddo
352                  enddo 
353               enddo
[662]354
355            CASE (5)
[396]356               do mm = imin(5),imax(5)
357               do ll = imin(4),imax(4)
358               do kk = imin(3),imax(3)
359               do jj = imin(2),imax(2)
[779]360!CDIR ALTCODE
361!CDIR SHORTLOOP
[396]362               do ii = imin(1),imax(1)
363                    ValParent = parent%var%array5(ii,jj,kk,ll,mm)
364                    if ( ValParent .NE. Agrif_SpecialValue) then
365                        Res = Res + ValParent
366                        Nbvals = Nbvals + 1
367                    endif
368                                    enddo
369                              enddo
370                        enddo
371                  enddo 
372               enddo
[662]373
374            CASE (6)
[396]375               do nn = imin(6),imax(6)
376               do mm = imin(5),imax(5)
377               do ll = imin(4),imax(4)
378               do kk = imin(3),imax(3)
379               do jj = imin(2),imax(2)
[779]380!CDIR ALTCODE
381!CDIR SHORTLOOP
[396]382               do ii = imin(1),imax(1)
383                    ValParent = parent%var%array6(ii,jj,kk,ll,mm,nn)
384                    if ( ValParent .NE. Agrif_SpecialValue) then
385                        Res = Res + ValParent
386                        Nbvals = Nbvals + 1
387                    endif
388                                          enddo
389                                    enddo
390                              enddo
391                        enddo
392                  enddo 
393               enddo
[662]394
395            END SELECT
[396]396C
397C
398           
399            if (Nbvals.GT.0) then
400              if (firsttest) then
401                   firsttest = .FALSE.
[662]402                   i=1
[396]403                   cycle
404              endif
[662]405            SELECT CASE(nbdim)
406            CASE (1)             
407              tempP%var%array1(indic(1)) 
[396]408     &           = Res/Nbvals
[662]409            CASE (2)
410              tempP%var%array2(indic(1),
[396]411     &                            indic(2)) = Res/Nbvals
[662]412            CASE (3)
413              tempP%var%array3(indic(1),
[396]414     &                            indic(2),indic(3)) = Res/Nbvals
[662]415            CASE (4)
416              tempP%var%array4(indic(1),
[396]417     &                            indic(2),indic(3),indic(4))
418     &                = Res/Nbvals
[662]419            CASE (5)
420              tempP%var%array5(indic(1),
[396]421     &                            indic(2),indic(3),indic(4),
422     &                   indic(5)) = Res/Nbvals
[662]423            CASE (6)
424              tempP%var%array6(indic(1),
[396]425     &                            indic(2),indic(3),indic(4),
426     &                           indic(5),indic(6)) = Res/Nbvals
[662]427            END SELECT
[396]428              exit
429            else
430               if (firsttest) exit
431               i = i + 1                     
432            endif
433          enddo           
434C     
435      End Subroutine CalculNewValTempP
436C
437C
[779]438      End Module Agrif_Mask     
439     
440      Subroutine CalculNewValTempP3D(indic,
441     &               tempP,parent,ppbtab,
442     &               ppetab,noraftab,MaxSearch,Agrif_SpecialValue)
443C
444CCC   Description:
445CCC   Subroutine called in the procedure Agrif_InterpnD to recalculate the value
446CCC   of the parent grid variable when this one is equal to Agrif_SpecialValue. 
447C
448C     Declarations:
449C
450       
451C
452C     Arrays arguments     
453      INTEGER, PARAMETER :: nbdim = 3
454      INTEGER,DIMENSION(nbdim) :: indic
455      LOGICAL,DIMENSION(nbdim) :: noraftab
456      INTEGER,DIMENSION(nbdim) :: ppbtab,ppetab
457      REAL :: Agrif_SpecialValue
458C
459C     Pointer argument
460      REAL,DIMENSION(ppbtab(1):ppetab(1),ppbtab(2):ppetab(2),
461     &    ppbtab(3):ppetab(3)) :: tempP, parent  ! Part of the parent grid used for
462                                      ! the interpolation of the child grid                     
463C
464C     Local scalar
465      INTEGER                  :: i,ii,iii,jj,kk,ll,mm,nn 
466      INTEGER,DIMENSION(nbdim) :: imin,imax,idecal
467      INTEGER                  :: Nbvals
468      REAL                     :: Res
469      REAL                     :: ValParent
470      INTEGER                  :: ValMax
471      INTEGER :: MaxSearch
472      LOGICAL :: Existunmasked
473C     
474C     Local arrays     
475C
476      ValMax = 1
477!CDIR NOVECTOR
478      do iii = 1 , nbdim
479         IF (.NOT.noraftab(iii)) THEN
480         ValMax = max(ValMax,ppetab(iii)-indic(iii))
481         ValMax = max(ValMax,indic(iii)-ppbtab(iii))
482         ENDIF
483      enddo
484C
485      Valmax = min(Valmax,MaxSearch)
486C
487!CDIR NOVECTOR
488      imin = indic
489!CDIR NOVECTOR
490      imax = indic
491     
492!CDIR NOVECTOR
493         idecal = indic 
494C
495         i = Valmax
496
497            do iii = 1 , nbdim
498               if (.NOT.noraftab(iii)) then
499                  imin(iii) = max(indic(iii) - i,ppbtab(iii))
500                  imax(iii) = min(indic(iii) + i,ppetab(iii))
501         
502                  if (indic(iii).GT.ppbtab(iii)) then
503
504                     idecal(iii) = idecal(iii)-1
505
506                        if (tempP(idecal(1),idecal(2),idecal(3)) 
507     &                               == Agrif_SpecialValue) then
508                           imin(iii) = imax(iii)
509                        endif
510
511                     idecal(iii) = idecal(iii)+1
512                  endif
513
514               endif           
515            enddo
516C
517            Existunmasked = .FALSE.
518C
519               do kk = imin(3),imax(3)
520               do jj = imin(2),imax(2)
521!CDIR NOVECTOR
522               do ii = imin(1),imax(1)
523                    if ( parent(ii,jj,kk) .NE. Agrif_SpecialValue) then
524                        Existunmasked = .TRUE.
525                       exit
526                    endif
527                        enddo
528                  enddo 
529               enddo
530C
531C
532          if (.Not.Existunmasked) return
533C
534         i = 1
535C
536         do While(i <= ValMax)
537C
538
539            do iii = 1 , nbdim
540               if (.NOT.noraftab(iii)) then
541                  imin(iii) = max(indic(iii) - i,ppbtab(iii))
542                  imax(iii) = min(indic(iii) + i,ppetab(iii))
543               endif           
544            enddo
545C
546            Res = 0.
547            Nbvals = 0
548C
549               do kk = imin(3),imax(3)
550               do jj = imin(2),imax(2)
551!CDIR NOVECTOR
552               do ii = imin(1),imax(1)
553                    ValParent = parent(ii,jj,kk)
554                    if ( ValParent .NE. Agrif_SpecialValue) then
555                        Res = Res + ValParent
556                        Nbvals = Nbvals + 1
557                    endif
558                        enddo
559                  enddo 
560               enddo
561C
562C
563           
564            if (Nbvals.GT.0) then
565              tempP(indic(1),indic(2),indic(3)) = Res/Nbvals
566              exit
567            else
568               i = i + 1                     
569            endif
570          enddo           
571C     
572      End Subroutine CalculNewValTempP3D       
573
574      Subroutine CalculNewValTempP4D(indic,
575     &               tempP,parent,ppbtab,
576     &               ppetab,noraftab,MaxSearch,Agrif_SpecialValue)
577C
578CCC   Description:
579CCC   Subroutine called in the procedure Agrif_InterpnD to recalculate the value
580CCC   of the parent grid variable when this one is equal to Agrif_SpecialValue. 
581C
582C     Declarations:
583C
584       
585C
586C     Arrays arguments     
587      INTEGER, PARAMETER :: nbdim = 4
588      INTEGER,DIMENSION(nbdim) :: indic
589      LOGICAL,DIMENSION(nbdim) :: noraftab
590      INTEGER,DIMENSION(nbdim) :: ppbtab,ppetab
591      INTEGER :: MaxSearch     
592      REAL :: Agrif_SpecialValue     
593C
594C     Pointer argument
595      REAL,DIMENSION(ppbtab(1):ppetab(1),ppbtab(2):ppetab(2),
596     &    ppbtab(3):ppetab(3),
597     &    ppbtab(4):ppetab(4)) :: tempP, parent  ! Part of the parent grid used for
598                                      ! the interpolation of the child grid                     
599C
600C     Local scalar
601      INTEGER                  :: i,ii,iii,jj,kk,ll,mm,nn 
602      INTEGER,DIMENSION(nbdim) :: imin,imax,idecal
603      INTEGER                  :: Nbvals
604      REAL                     :: Res
605      REAL                     :: ValParent
606      INTEGER                  :: ValMax
607      LOGICAL                  :: firsttest
608C     
609C     Local arrays     
610C
611      ValMax = 1
612      do iii = 1 , nbdim
613         IF (.NOT.noraftab(iii)) THEN
614         ValMax = max(ValMax,ppetab(iii)-indic(iii))
615         ValMax = max(ValMax,indic(iii)-ppbtab(iii))
616         ENDIF
617      enddo
618C
619      Valmax = min(Valmax,MaxSearch)
620C
621      imin = indic
622      imax = indic
623C
624         i = 1
625         firsttest = .TRUE.
626         idecal = indic 
627   
628C
629         do While(i <= ValMax)
630C
631         IF ((i == 1).AND.(firsttest)) i = Valmax
632
633            do iii = 1 , nbdim
634               if (.NOT.noraftab(iii)) then
635                  imin(iii) = max(indic(iii) - i,ppbtab(iii))
636                  imax(iii) = min(indic(iii) + i,ppetab(iii))
637                  if (firsttest) then               
638                  if (indic(iii).GT.ppbtab(iii)) then
639
640
641                     idecal(iii) = idecal(iii)-1
642           
643                     if (tempP(idecal(1),idecal(2),idecal(3),idecal(4))
644     &                               == Agrif_SpecialValue) then
645                           imin(iii) = imax(iii)
646                     endif
647           
648           idecal(iii) = idecal(iii)+1
649                  endif
650                  endif
651               endif           
652            enddo
653C
654            Res = 0.
655            Nbvals = 0
656C
657               do ll = imin(4),imax(4)
658               do kk = imin(3),imax(3)
659               do jj = imin(2),imax(2)
660               do ii = imin(1),imax(1)
661                    ValParent = parent(ii,jj,kk,ll)
662                    if ( ValParent .NE. Agrif_SpecialValue) then
663                        Res = Res + ValParent
664                        Nbvals = Nbvals + 1
665                    endif
666                              enddo
667                        enddo
668                  enddo 
669               enddo
670C
671C
672           
673            if (Nbvals.GT.0) then
674              if (firsttest) then
675                   firsttest = .FALSE.
676                   i=1
677                   cycle
678              endif
679         
680              tempP(indic(1),indic(2),indic(3),indic(4)) = Res/Nbvals
681       
682              exit
683            else
684               if (firsttest) exit
685               i = i + 1                     
686            endif
687          enddo           
688C     
689      End Subroutine CalculNewValTempP4D       
Note: See TracBrowser for help on using the repository browser.