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
Line 
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
37      Integer, Parameter :: MaxSearch = 3
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
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         
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
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)
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 
213      INTEGER,DIMENSION(nbdim) :: imin,imax,idecal
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
232!CDIR NOVECTOR
233      imin = indic
234!CDIR NOVECTOR
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))
248                  if (firsttest) then               
249                  if (indic(iii).GT.ppbtab(iii)) then
250
251!CDIR NOVECTOR
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
292               endif           
293            enddo
294C
295            Res = 0.
296            Nbvals = 0
297C
298            SELECT CASE(nbdim)
299            CASE (1)
300!CDIR ALTCODE
301!CDIR SHORTLOOP
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
310            CASE (2)
311               do jj = imin(2),imax(2)
312!CDIR ALTCODE
313!CDIR SHORTLOOP
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
322               
323            CASE (3)
324               do kk = imin(3),imax(3)
325               do jj = imin(2),imax(2)
326!CDIR ALTCODE
327!CDIR SHORTLOOP
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
337
338            CASE (4)
339               do ll = imin(4),imax(4)
340               do kk = imin(3),imax(3)
341               do jj = imin(2),imax(2)
342!CDIR ALTCODE
343!CDIR SHORTLOOP
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
354
355            CASE (5)
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)
360!CDIR ALTCODE
361!CDIR SHORTLOOP
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
373
374            CASE (6)
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)
380!CDIR ALTCODE
381!CDIR SHORTLOOP
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
394
395            END SELECT
396C
397C
398           
399            if (Nbvals.GT.0) then
400              if (firsttest) then
401                   firsttest = .FALSE.
402                   i=1
403                   cycle
404              endif
405            SELECT CASE(nbdim)
406            CASE (1)             
407              tempP%var%array1(indic(1)) 
408     &           = Res/Nbvals
409            CASE (2)
410              tempP%var%array2(indic(1),
411     &                            indic(2)) = Res/Nbvals
412            CASE (3)
413              tempP%var%array3(indic(1),
414     &                            indic(2),indic(3)) = Res/Nbvals
415            CASE (4)
416              tempP%var%array4(indic(1),
417     &                            indic(2),indic(3),indic(4))
418     &                = Res/Nbvals
419            CASE (5)
420              tempP%var%array5(indic(1),
421     &                            indic(2),indic(3),indic(4),
422     &                   indic(5)) = Res/Nbvals
423            CASE (6)
424              tempP%var%array6(indic(1),
425     &                            indic(2),indic(3),indic(4),
426     &                           indic(5),indic(6)) = Res/Nbvals
427            END SELECT
428              exit
429            else
430               if (firsttest) exit
431               i = i + 1                     
432            endif
433          enddo           
434C     
435      End Subroutine CalculNewValTempP
436C
437C
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.