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_top_interp.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_top_interp.F90 @ 10989

Last change on this file since 10989 was 10989, checked in by acc, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert NST routines in preparation for getting AGRIF back up and running. AGRIF conv stage now works but requires some renaming of recently changes DIU modules (included in this commit). AGRIF compile and link stage not yet working (agrif routines need to be passed the time-level indices) but non-AGRIF SETTE tests are all OK
  • Property svn:keywords set to Id
File size: 12.5 KB
Line 
1MODULE agrif_top_interp
2   !!======================================================================
3   !!                   ***  MODULE  agrif_top_interp  ***
4   !! AGRIF: interpolation package for TOP
5   !!======================================================================
6   !! History :  2.0  !  ???
7   !!----------------------------------------------------------------------
8#if defined key_agrif && defined key_top
9   !!----------------------------------------------------------------------
10   !!   'key_agrif'                                              AGRIF zoom
11   !!   'key_top'                                           on-line tracers
12   !!----------------------------------------------------------------------
13   USE par_oce
14   USE oce
15   USE dom_oce     
16   USE agrif_oce
17   USE agrif_top_sponge
18   USE par_trc
19   USE trc
20   !
21   USE lib_mpp     ! MPP library
22
23   IMPLICIT NONE
24   PRIVATE
25
26   PUBLIC Agrif_trc, interptrn
27
28  !!----------------------------------------------------------------------
29   !! NEMO/NST 4.0 , NEMO Consortium (2018)
30   !! $Id$
31   !! Software governed by the CeCILL license (see ./LICENSE)
32   !!----------------------------------------------------------------------
33CONTAINS
34
35   SUBROUTINE Agrif_trc
36      !!----------------------------------------------------------------------
37      !!                   ***  ROUTINE Agrif_trc  ***
38      !!----------------------------------------------------------------------
39      !
40      IF( Agrif_Root() )   RETURN
41      !
42      Agrif_SpecialValue    = 0._wp
43      Agrif_UseSpecialValue = .TRUE.
44      !
45      CALL Agrif_Bc_variable( trn_id, procname=interptrn )
46      Agrif_UseSpecialValue = .FALSE.
47      !
48   END SUBROUTINE Agrif_trc
49
50   SUBROUTINE interptrn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
51      !!----------------------------------------------------------------------
52      !!                  *** ROUTINE interptrn ***
53      !!----------------------------------------------------------------------
54      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   ptab
55      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
56      LOGICAL                                     , INTENT(in   ) ::   before
57      INTEGER                                     , INTENT(in   ) ::   nb , ndir
58      !
59      INTEGER  ::   ji, jj, jk, jn, iref, jref, ibdy, jbdy   ! dummy loop indices
60      INTEGER  ::   imin, imax, jmin, jmax, N_in, N_out
61      REAL(wp) ::   zrho, z1, z2, z3, z4, z5, z6, z7
62      LOGICAL :: western_side, eastern_side,northern_side,southern_side
63      ! vertical interpolation:
64      REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child
65      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
66      REAL(wp), DIMENSION(k1:k2) :: h_in
67      REAL(wp), DIMENSION(1:jpk) :: h_out
68      REAL(wp) :: h_diff
69
70      IF( before ) THEN         
71         DO jn = 1,jptra
72            DO jk=k1,k2
73               DO jj=j1,j2
74                 DO ji=i1,i2
75                       ptab(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm)
76                 END DO
77              END DO
78           END DO
79        END DO
80
81# if defined key_vertical
82        DO jk=k1,k2
83           DO jj=j1,j2
84              DO ji=i1,i2
85                 ptab(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) 
86              END DO
87           END DO
88        END DO
89# endif
90      ELSE
91
92         western_side  = (nb == 1).AND.(ndir == 1)   ;   eastern_side  = (nb == 1).AND.(ndir == 2)
93         southern_side = (nb == 2).AND.(ndir == 1)   ;   northern_side = (nb == 2).AND.(ndir == 2)
94
95# if defined key_vertical             
96         DO jj=j1,j2
97            DO ji=i1,i2
98               iref = ji
99               jref = jj
100               if(western_side) iref=MAX(2,ji)
101               if(eastern_side) iref=MIN(nlci-1,ji)
102               if(southern_side) jref=MAX(2,jj)
103               if(northern_side) jref=MIN(nlcj-1,jj)
104               N_in = 0
105               DO jk=k1,k2 !k2 = jpk of parent grid
106                  IF (ptab(ji,jj,jk,n2) == 0) EXIT
107                  N_in = N_in + 1
108                  tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)
109                  h_in(N_in) = ptab(ji,jj,jk,n2)
110               END DO
111               N_out = 0
112               DO jk=1,jpk ! jpk of child grid
113                  IF (tmask(iref,jref,jk) == 0) EXIT
114                  N_out = N_out + 1
115                  h_out(jk) = e3t(iref,jref,jk,Kmm)
116               ENDDO
117               IF (N_in > 0) THEN
118                  DO jn=1,jptra
119                     call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
120                  ENDDO
121               ENDIF
122            ENDDO
123         ENDDO
124# else
125         ptab_child(i1:i2,j1:j2,1:jpk,1:jptra) = ptab(i1:i2,j1:j2,1:jpk,1:jptra)
126# endif
127         !
128         DO jn=1, jptra
129            tr(i1:i2,j1:j2,1:jpk,jn,Krhs)=ptab_child(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
130         END DO
131
132         IF ( .NOT.lk_agrif_clp ) THEN 
133            !
134            imin = i1 ; imax = i2
135            jmin = j1 ; jmax = j2
136            !
137            ! Remove CORNERS
138            IF((nbondj == -1).OR.(nbondj == 2)) jmin = 2 + nbghostcells
139            IF((nbondj == +1).OR.(nbondj == 2)) jmax = nlcj - nbghostcells - 1
140            IF((nbondi == -1).OR.(nbondi == 2)) imin = 2 + nbghostcells
141            IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci - nbghostcells - 1     
142            !
143            IF( eastern_side ) THEN
144               zrho = Agrif_Rhox()
145               z1 = ( zrho - 1._wp ) * 0.5_wp                   
146               z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )         
147               z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )
148               z7 =       - ( zrho - 1._wp ) / ( zrho + 3._wp )
149               z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7
150               !
151               ibdy = nlci-nbghostcells
152               DO jn = 1, jptra
153                  tr(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs) = z1 * ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
154                  DO jk = 1, jpkm1
155                     DO jj = jmin,jmax
156                        IF( umask(ibdy-1,jj,jk) == 0._wp ) THEN
157                           tr(ibdy,jj,jk,jn,Krhs) = tr(ibdy+1,jj,jk,jn,Krhs) * tmask(ibdy,jj,jk)
158                        ELSE
159                           tr(ibdy,jj,jk,jn,Krhs)=(z4*tr(ibdy+1,jj,jk,jn,Krhs)+z3*tr(ibdy-1,jj,jk,jn,Krhs))*tmask(ibdy,jj,jk)
160                           IF( uu(ibdy-1,jj,jk,Kmm) > 0._wp ) THEN
161                              tr(ibdy,jj,jk,jn,Krhs)=( z6*tr(ibdy-1,jj,jk,jn,Krhs)+z5*tr(ibdy+1,jj,jk,jn,Krhs) & 
162                                                 + z7*tr(ibdy-2,jj,jk,jn,Krhs) ) * tmask(ibdy,jj,jk)
163                           ENDIF
164                        ENDIF
165                     END DO
166                  END DO
167                  ! Restore ghost points:
168                  tr(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy+1,jmin:jmax,1:jpkm1)
169               END DO
170            ENDIF
171            !
172            IF( northern_side ) THEN
173               zrho = Agrif_Rhoy()
174               z1 = ( zrho - 1._wp ) * 0.5_wp                   
175               z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )         
176               z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )
177               z7 =       - ( zrho - 1._wp ) / ( zrho + 3._wp )
178               z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7
179               !
180               jbdy = nlcj-nbghostcells         
181               DO jn = 1, jptra
182                  tr(imin:imax,jbdy+1,1:jpkm1,jn,Krhs) = z1 * ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
183                  DO jk = 1, jpkm1
184                     DO ji = imin,imax
185                        IF( vmask(ji,jbdy-1,jk) == 0._wp ) THEN
186                           tr(ji,jbdy,jk,jn,Krhs) = tr(ji,jbdy+1,jk,jn,Krhs) * tmask(ji,jbdy,jk)
187                        ELSE
188                           tr(ji,jbdy,jk,jn,Krhs)=(z4*tr(ji,jbdy+1,jk,jn,Krhs)+z3*tr(ji,jbdy-1,jk,jn,Krhs))*tmask(ji,jbdy,jk)       
189                           IF (vv(ji,jbdy-1,jk,Kmm) > 0._wp ) THEN
190                              tr(ji,jbdy,jk,jn,Krhs)=( z6*tr(ji,jbdy-1,jk,jn,Krhs)+z5*tr(ji,jbdy+1,jk,jn,Krhs)  &
191                                                 + z7*tr(ji,jbdy-2,jk,jn,Krhs) ) * tmask(ji,jbdy,jk)
192                           ENDIF
193                        ENDIF
194                     END DO
195                  END DO
196                  ! Restore ghost points:
197                  tr(imin:imax,jbdy+1,1:jpkm1,jn,Krhs) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) * tmask(imin:imax,jbdy+1,1:jpkm1)
198               END DO
199            ENDIF
200            !
201            IF( western_side ) THEN
202               zrho = Agrif_Rhox()
203               z1 = ( zrho - 1._wp ) * 0.5_wp                   
204               z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )         
205               z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )
206               z7 =       - ( zrho - 1._wp ) / ( zrho + 3._wp )
207               z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7
208               !   
209               ibdy = 1+nbghostcells       
210               DO jn = 1, jptra
211                  tr(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs) = z1 * ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) + z2 * ptab_child(ibdy,jmin:jmax,1:jpkm1,jn)
212                  DO jk = 1, jpkm1
213                     DO jj = jmin,jmax
214                        IF( umask(ibdy,jj,jk) == 0._wp ) THEN
215                           tr(ibdy,jj,jk,jn,Krhs) = tr(ibdy-1,jj,jk,jn,Krhs) * tmask(ibdy,jj,jk)
216                        ELSE
217                           tr(ibdy,jj,jk,jn,Krhs)=(z4*tr(ibdy-1,jj,jk,jn,Krhs)+z3*tr(ibdy+1,jj,jk,jn,Krhs))*tmask(ibdy,jj,jk)       
218                           IF( uu(ibdy,jj,jk,Kmm) < 0._wp ) THEN
219                              tr(ibdy,jj,jk,jn,Krhs)=( z6*tr(ibdy+1,jj,jk,jn,Krhs)+z5*tr(ibdy-1,jj,jk,jn,Krhs) &
220                                                 + z7*tr(ibdy+2,jj,jk,jn,Krhs) ) * tmask(ibdy,jj,jk)
221                           ENDIF
222                        ENDIF
223                     END DO
224                  END DO
225                  ! Restore ghost points:
226                  tr(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) * tmask(ibdy-1,jmin:jmax,1:jpkm1)
227               END DO
228            ENDIF
229            !
230            IF( southern_side ) THEN
231               zrho = Agrif_Rhoy()
232               z1 = ( zrho - 1._wp ) * 0.5_wp                   
233               z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )         
234               z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )
235               z7 =       - ( zrho - 1._wp ) / ( zrho + 3._wp )
236               z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7
237               ! 
238               jbdy=1+nbghostcells       
239               DO jn = 1, jptra
240                  tr(imin:imax,jbdy-1,1:jpkm1,jn,Krhs) = z1 * ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) + z2 * ptab_child(imin:imax,jbdy,1:jpkm1,jn)
241                  DO jk = 1, jpkm1     
242                     DO ji = imin,imax
243                        IF( vmask(ji,jbdy,jk) == 0._wp ) THEN
244                           tr(ji,jbdy,jk,jn,Krhs)=tr(ji,jbdy-1,jk,jn,Krhs) * tmask(ji,jbdy,jk)
245                        ELSE
246                           tr(ji,jbdy,jk,jn,Krhs)=(z4*tr(ji,jbdy-1,jk,jn,Krhs)+z3*tr(ji,jbdy+1,jk,jn,Krhs))*tmask(ji,jbdy,jk)
247                           IF( vv(ji,jbdy,jk,Kmm) < 0._wp ) THEN
248                              tr(ji,jbdy,jk,jn,Krhs)=( z6*tr(ji,jbdy+1,jk,jn,Krhs)+z5*tr(ji,jbdy-1,jk,jn,Krhs) & 
249                                                 + z7*tr(ji,jbdy+2,jk,jn,Krhs) ) * tmask(ji,jbdy,jk)
250                           ENDIF
251                        ENDIF
252                     END DO
253                  END DO
254                  ! Restore ghost points:
255                  tr(imin:imax,jbdy-1,1:jpkm1,jn,Krhs) = ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) * tmask(imin:imax,jbdy-1,1:jpkm1)
256               END DO
257            ENDIF
258            !
259         ENDIF
260
261      ENDIF
262      !
263   END SUBROUTINE interptrn
264
265#else
266   !!----------------------------------------------------------------------
267   !!   Empty module                                           no TOP AGRIF
268   !!----------------------------------------------------------------------
269CONTAINS
270   SUBROUTINE Agrif_TOP_Interp_empty
271      !!---------------------------------------------
272      !!   *** ROUTINE agrif_Top_Interp_empty ***
273      !!---------------------------------------------
274      WRITE(*,*)  'agrif_top_interp : You should not have seen this print! error?'
275   END SUBROUTINE Agrif_TOP_Interp_empty
276#endif
277
278   !!======================================================================
279END MODULE agrif_top_interp
Note: See TracBrowser for help on using the repository browser.