source: NEMO/trunk/src/NST/agrif_top_interp.F90 @ 10068

Last change on this file since 10068 was 10068, checked in by nicolasmartin, 2 years ago

First part of modifications to have a common default header : fix typos and SVN keywords properties

  • Property svn:keywords set to Id
File size: 12.3 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) = trn(ji,jj,jk,jn)
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_n(ji,jj,jk) 
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_n(iref,jref,jk)
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            tra(i1:i2,j1:j2,1:jpk,jn)=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                  tra(ibdy+1,jmin:jmax,1:jpkm1,jn) = 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                           tra(ibdy,jj,jk,jn) = tra(ibdy+1,jj,jk,jn) * tmask(ibdy,jj,jk)
158                        ELSE
159                           tra(ibdy,jj,jk,jn)=(z4*tra(ibdy+1,jj,jk,jn)+z3*tra(ibdy-1,jj,jk,jn))*tmask(ibdy,jj,jk)
160                           IF( un(ibdy-1,jj,jk) > 0._wp ) THEN
161                              tra(ibdy,jj,jk,jn)=( z6*tra(ibdy-1,jj,jk,jn)+z5*tra(ibdy+1,jj,jk,jn) & 
162                                                 + z7*tra(ibdy-2,jj,jk,jn) ) * tmask(ibdy,jj,jk)
163                           ENDIF
164                        ENDIF
165                     END DO
166                  END DO
167                  ! Restore ghost points:
168                  tra(ibdy+1,jmin:jmax,1:jpkm1,jn) = 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                  tra(imin:imax,jbdy+1,1:jpkm1,jn) = 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                           tra(ji,jbdy,jk,jn) = tra(ji,jbdy+1,jk,jn) * tmask(ji,jbdy,jk)
187                        ELSE
188                           tra(ji,jbdy,jk,jn)=(z4*tra(ji,jbdy+1,jk,jn)+z3*tra(ji,jbdy-1,jk,jn))*tmask(ji,jbdy,jk)       
189                           IF (vn(ji,jbdy-1,jk) > 0._wp ) THEN
190                              tra(ji,jbdy,jk,jn)=( z6*tra(ji,jbdy-1,jk,jn)+z5*tra(ji,jbdy+1,jk,jn)  &
191                                                 + z7*tra(ji,jbdy-2,jk,jn) ) * tmask(ji,jbdy,jk)
192                           ENDIF
193                        ENDIF
194                     END DO
195                  END DO
196                  ! Restore ghost points:
197                  tra(imin:imax,jbdy+1,1:jpkm1,jn) = 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                  tra(ibdy-1,jmin:jmax,1:jpkm1,jn) = 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                           tra(ibdy,jj,jk,jn) = tra(ibdy-1,jj,jk,jn) * tmask(ibdy,jj,jk)
216                        ELSE
217                           tra(ibdy,jj,jk,jn)=(z4*tra(ibdy-1,jj,jk,jn)+z3*tra(ibdy+1,jj,jk,jn))*tmask(ibdy,jj,jk)       
218                           IF( un(ibdy,jj,jk) < 0._wp ) THEN
219                              tra(ibdy,jj,jk,jn)=( z6*tra(ibdy+1,jj,jk,jn)+z5*tra(ibdy-1,jj,jk,jn) &
220                                                 + z7*tra(ibdy+2,jj,jk,jn) ) * tmask(ibdy,jj,jk)
221                           ENDIF
222                        ENDIF
223                     END DO
224                  END DO
225                  ! Restore ghost points:
226                  tra(ibdy-1,jmin:jmax,1:jpkm1,jn) = 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                  tra(imin:imax,jbdy-1,1:jpkm1,jn) = 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                           tra(ji,jbdy,jk,jn)=tra(ji,jbdy-1,jk,jn) * tmask(ji,jbdy,jk)
245                        ELSE
246                           tra(ji,jbdy,jk,jn)=(z4*tra(ji,jbdy-1,jk,jn)+z3*tra(ji,jbdy+1,jk,jn))*tmask(ji,jbdy,jk)
247                           IF( vn(ji,jbdy,jk) < 0._wp ) THEN
248                              tra(ji,jbdy,jk,jn)=( z6*tra(ji,jbdy+1,jk,jn)+z5*tra(ji,jbdy-1,jk,jn) & 
249                                                 + z7*tra(ji,jbdy+2,jk,jn) ) * tmask(ji,jbdy,jk)
250                           ENDIF
251                        ENDIF
252                     END DO
253                  END DO
254                  ! Restore ghost points:
255                  tra(imin:imax,jbdy-1,1:jpkm1,jn) = 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.