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_rewrite_time_filterswap/src/NST – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/NST/agrif_top_interp.F90 @ 11463

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

Branch: dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap. Minor bugfix in step.F90 to enable AGRIF SETTE tests to run. Also merged prettification changes to NST routines from the dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps branch. Neither AGRIF_DEMO nor VORTEX restart perfectly (drifting after 8 and 121 timesteps, respectively).
  • Property svn:keywords set to Id
File size: 14.1 KB
RevLine 
[636]1MODULE agrif_top_interp
[9019]2   !!======================================================================
3   !!                   ***  MODULE  agrif_top_interp  ***
4   !! AGRIF: interpolation package for TOP
5   !!======================================================================
6   !! History :  2.0  !  ???
7   !!----------------------------------------------------------------------
[1206]8#if defined key_agrif && defined key_top
[9019]9   !!----------------------------------------------------------------------
10   !!   'key_agrif'                                              AGRIF zoom
11   !!   'key_top'                                           on-line tracers
12   !!----------------------------------------------------------------------
[636]13   USE par_oce
14   USE oce
15   USE dom_oce     
[782]16   USE agrif_oce
[2715]17   USE agrif_top_sponge
[5656]18   USE par_trc
[1271]19   USE trc
[9019]20   !
21   USE lib_mpp     ! MPP library
[628]22
[636]23   IMPLICIT NONE
24   PRIVATE
[628]25
[5656]26   PUBLIC Agrif_trc, interptrn
[636]27
[2715]28  !!----------------------------------------------------------------------
[9598]29   !! NEMO/NST 4.0 , NEMO Consortium (2018)
[1156]30   !! $Id$
[10068]31   !! Software governed by the CeCILL license (see ./LICENSE)
[1156]32   !!----------------------------------------------------------------------
[6140]33CONTAINS
[1156]34
[1271]35   SUBROUTINE Agrif_trc
[3680]36      !!----------------------------------------------------------------------
[9019]37      !!                   ***  ROUTINE Agrif_trc  ***
[3680]38      !!----------------------------------------------------------------------
39      !
40      IF( Agrif_Root() )   RETURN
[9019]41      !
42      Agrif_SpecialValue    = 0._wp
[636]43      Agrif_UseSpecialValue = .TRUE.
[9019]44      !
[5656]45      CALL Agrif_Bc_variable( trn_id, procname=interptrn )
[636]46      Agrif_UseSpecialValue = .FALSE.
[5656]47      !
48   END SUBROUTINE Agrif_trc
[636]49
[9019]50   SUBROUTINE interptrn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
51      !!----------------------------------------------------------------------
[9788]52      !!                  *** ROUTINE interptrn ***
[9019]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
[5656]58      !
[9788]59      INTEGER  ::   ji, jj, jk, jn, iref, jref, ibdy, jbdy   ! dummy loop indices
[9031]60      INTEGER  ::   imin, imax, jmin, jmax, N_in, N_out
[9806]61      REAL(wp) ::   zrho, z1, z2, z3, z4, z5, z6, z7
[9031]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
[9788]67      REAL(wp), DIMENSION(1:jpk) :: h_out
68      REAL(wp) :: h_diff
[9031]69
[9788]70      IF( before ) THEN         
71         DO jn = 1,jptra
[9031]72            DO jk=k1,k2
73               DO jj=j1,j2
74                 DO ji=i1,i2
[11053]75                       ptab(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm_a)
[9031]76                 END DO
[9788]77              END DO
78           END DO
79        END DO
80
[9031]81# if defined key_vertical
82        DO jk=k1,k2
83           DO jj=j1,j2
84              DO ji=i1,i2
[11053]85                 ptab(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a) 
[9031]86              END DO
87           END DO
88        END DO
89# endif
[9788]90      ELSE
[9031]91
[9788]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)
[9031]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
[11053]115                  h_out(jk) = e3t(iref,jref,jk,Kmm_a)
[9031]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
[9788]127         !
128         DO jn=1, jptra
[11053]129            tr(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab_child(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk) 
[9788]130         END DO
[9031]131
[9788]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
[9806]140            IF((nbondi == -1).OR.(nbondi == 2)) imin = 2 + nbghostcells
[9788]141            IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci - nbghostcells - 1     
142            !
143            IF( eastern_side ) THEN
[9806]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               !
[9788]151               ibdy = nlci-nbghostcells
152               DO jn = 1, jptra
[11463]153                  tr(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs_a) =  z1 * ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn)  &
154                    &                                     + z2 * ptab_child(ibdy  ,jmin:jmax,1:jpkm1,jn)
[9788]155                  DO jk = 1, jpkm1
156                     DO jj = jmin,jmax
157                        IF( umask(ibdy-1,jj,jk) == 0._wp ) THEN
[11053]158                           tr(ibdy,jj,jk,jn,Krhs_a) = tr(ibdy+1,jj,jk,jn,Krhs_a) * tmask(ibdy,jj,jk)
[9788]159                        ELSE
[11463]160                           tr(ibdy,jj,jk,jn,Krhs_a) = (  z4 * tr(ibdy+1,jj,jk,jn,Krhs_a)     &
161                             &                         + z3 * tr(ibdy-1,jj,jk,jn,Krhs_a)     &
162                             &                        ) *  tmask(ibdy  ,jj,jk)
[11053]163                           IF( uu(ibdy-1,jj,jk,Kmm_a) > 0._wp ) THEN
[11463]164                              tr(ibdy,jj,jk,jn,Krhs_a) = (  z6 * tr(ibdy-1,jj,jk,jn,Krhs_a)  &
165                                &                         + z5 * tr(ibdy+1,jj,jk,jn,Krhs_a)  & 
166                                &                         + z7 * tr(ibdy-2,jj,jk,jn,Krhs_a)  &
167                                &                         ) * tmask(ibdy  ,jj,jk)
[9788]168                           ENDIF
[5656]169                        ENDIF
[9019]170                     END DO
[5656]171                  END DO
[9788]172                  ! Restore ghost points:
[11463]173                  tr(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs_a) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) &
174                    &                                     *     tmask(ibdy+1,jmin:jmax,1:jpkm1)
[9031]175               END DO
[9788]176            ENDIF
177            !
178            IF( northern_side ) THEN
[9806]179               zrho = Agrif_Rhoy()
180               z1 = ( zrho - 1._wp ) * 0.5_wp                   
181               z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )         
182               z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )
183               z7 =       - ( zrho - 1._wp ) / ( zrho + 3._wp )
184               z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7
185               !
[9788]186               jbdy = nlcj-nbghostcells         
187               DO jn = 1, jptra
[11463]188                  tr(imin:imax,jbdy+1,1:jpkm1,jn,Krhs_a) =  z1 * ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) &
189                    &                                     + z2 * ptab_child(imin:imax,jbdy  ,1:jpkm1,jn)
[9788]190                  DO jk = 1, jpkm1
191                     DO ji = imin,imax
192                        IF( vmask(ji,jbdy-1,jk) == 0._wp ) THEN
[11053]193                           tr(ji,jbdy,jk,jn,Krhs_a) = tr(ji,jbdy+1,jk,jn,Krhs_a) * tmask(ji,jbdy,jk)
[9788]194                        ELSE
[11463]195                           tr(ji,jbdy,jk,jn,Krhs_a) = ( z4 * tr(ji,jbdy+1,jk,jn,Krhs_a) 
196                             &                        + z3 * tr(ji,jbdy-1,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)       
[11053]197                           IF (vv(ji,jbdy-1,jk,Kmm_a) > 0._wp ) THEN
[11463]198                              tr(ji,jbdy,jk,jn,Krhs_a) = (  z6 * tr(ji,jbdy-1,jk,jn,Krhs_a)                       &
199                                &                         + z5 * tr(ji,jbdy+1,jk,jn,Krhs_a)                       &
200                                &                         + z7 * tr(ji,jbdy-2,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)
[9788]201                           ENDIF
[5656]202                        ENDIF
[9019]203                     END DO
[5656]204                  END DO
[9788]205                  ! Restore ghost points:
[11463]206                  tr(imin:imax,jbdy+1,1:jpkm1,jn,Krhs_a) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn)  &
207                    &                                     *     tmask(imin:imax,jbdy+1,1:jpkm1)
[636]208               END DO
[9788]209            ENDIF
210            !
[9806]211            IF( western_side ) THEN
212               zrho = Agrif_Rhox()
213               z1 = ( zrho - 1._wp ) * 0.5_wp                   
214               z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )         
215               z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )
216               z7 =       - ( zrho - 1._wp ) / ( zrho + 3._wp )
217               z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7
218               !   
[9788]219               ibdy = 1+nbghostcells       
220               DO jn = 1, jptra
[11463]221                  tr(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs_a) =  z1 * ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn)         &
222                    &                                     + z2 * ptab_child(ibdy  ,jmin:jmax,1:jpkm1,jn)
[9788]223                  DO jk = 1, jpkm1
224                     DO jj = jmin,jmax
225                        IF( umask(ibdy,jj,jk) == 0._wp ) THEN
[11053]226                           tr(ibdy,jj,jk,jn,Krhs_a) = tr(ibdy-1,jj,jk,jn,Krhs_a) * tmask(ibdy,jj,jk)
[9788]227                        ELSE
[11463]228                           tr(ibdy,jj,jk,jn,Krhs_a) = (  z4 * tr(ibdy-1,jj,jk,jn,Krhs_a)                          &
229                             &                         + z3 * tr(ibdy+1,jj,jk,jn,Krhs_a) ) * tmask(ibdy,jj,jk)       
[11053]230                           IF( uu(ibdy,jj,jk,Kmm_a) < 0._wp ) THEN
[11463]231                              tr(ibdy,jj,jk,jn,Krhs_a) = (  z6 * tr(ibdy+1,jj,jk,jn,Krhs_a)                       &
232                                &                         + z5 * tr(ibdy-1,jj,jk,jn,Krhs_a)                       &
233                                &                         + z7 * tr(ibdy+2,jj,jk,jn,Krhs_a) ) * tmask(ibdy,jj,jk)
[9788]234                           ENDIF
[5656]235                        ENDIF
[9019]236                     END DO
[5656]237                  END DO
[9788]238                  ! Restore ghost points:
[11463]239                  tr(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs_a) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn)                &
240                    &                                     *     tmask(ibdy-1,jmin:jmax,1:jpkm1)
[636]241               END DO
[9788]242            ENDIF
243            !
[9806]244            IF( southern_side ) THEN
245               zrho = Agrif_Rhoy()
246               z1 = ( zrho - 1._wp ) * 0.5_wp                   
247               z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )         
248               z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )
249               z7 =       - ( zrho - 1._wp ) / ( zrho + 3._wp )
250               z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7
251               ! 
[9788]252               jbdy=1+nbghostcells       
253               DO jn = 1, jptra
[11463]254                  tr(imin:imax,jbdy-1,1:jpkm1,jn,Krhs_a) =  z1 * ptab_child(imin:imax,jbdy-1,1:jpkm1,jn)          &
255                    &                                     + z2 * ptab_child(imin:imax,jbdy  ,1:jpkm1,jn)
[9806]256                  DO jk = 1, jpkm1     
257                     DO ji = imin,imax
[9788]258                        IF( vmask(ji,jbdy,jk) == 0._wp ) THEN
[11463]259                           tr(ji,jbdy,jk,jn,Krhs_a) = tr(ji,jbdy-1,jk,jn,Krhs_a) * tmask(ji,jbdy,jk)
[9788]260                        ELSE
[11463]261                           tr(ji,jbdy,jk,jn,Krhs_a) = (  z4 * tr(ji,jbdy-1,jk,jn,Krhs_a)                          &
262                             &                         + z3 * tr(ji,jbdy+1,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)
[11053]263                           IF( vv(ji,jbdy,jk,Kmm_a) < 0._wp ) THEN
[11463]264                              tr(ji,jbdy,jk,jn,Krhs_a) = (  z6 * tr(ji,jbdy+1,jk,jn,Krhs_a)                       &
265                                &                         + z5 * tr(ji,jbdy-1,jk,jn,Krhs_a)                       & 
266                                &                         + z7 * tr(ji,jbdy+2,jk,jn,Krhs_a) ) * tmask(ji,jbdy,jk)
[9788]267                           ENDIF
[5656]268                        ENDIF
[9019]269                     END DO
[5656]270                  END DO
[9788]271                  ! Restore ghost points:
[11463]272                  tr(imin:imax,jbdy-1,1:jpkm1,jn,Krhs_a) = ptab_child(imin:imax,jbdy-1,1:jpkm1,jn)                &
273                    &                                     *     tmask(imin:imax,jbdy-1,1:jpkm1)
[636]274               END DO
[9019]275            ENDIF
276            !
[5656]277         ENDIF
[9806]278
[628]279      ENDIF
[3680]280      !
[5656]281   END SUBROUTINE interptrn
[2715]282
[628]283#else
[9019]284   !!----------------------------------------------------------------------
285   !!   Empty module                                           no TOP AGRIF
286   !!----------------------------------------------------------------------
[636]287CONTAINS
288   SUBROUTINE Agrif_TOP_Interp_empty
289      !!---------------------------------------------
290      !!   *** ROUTINE agrif_Top_Interp_empty ***
291      !!---------------------------------------------
292      WRITE(*,*)  'agrif_top_interp : You should not have seen this print! error?'
293   END SUBROUTINE Agrif_TOP_Interp_empty
[628]294#endif
[9019]295
296   !!======================================================================
[636]297END MODULE agrif_top_interp
Note: See TracBrowser for help on using the repository browser.