source: vendors/AGRIF/CMEMS_2020/AGRIF_FILES/modmask.F90 @ 10087

Last change on this file since 10087 was 10087, checked in by rblod, 2 years ago

update AGRIF library

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