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
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
[10989]75                       ptab(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kmm)
[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
[10989]85                 ptab(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm) 
[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
[10989]115                  h_out(jk) = e3t(iref,jref,jk,Kmm)
[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
[10989]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) 
[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
[10989]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)
[9788]154                  DO jk = 1, jpkm1
155                     DO jj = jmin,jmax
156                        IF( umask(ibdy-1,jj,jk) == 0._wp ) THEN
[10989]157                           tr(ibdy,jj,jk,jn,Krhs) = tr(ibdy+1,jj,jk,jn,Krhs) * tmask(ibdy,jj,jk)
[9788]158                        ELSE
[10989]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)
[9788]163                           ENDIF
[5656]164                        ENDIF
[9019]165                     END DO
[5656]166                  END DO
[9788]167                  ! Restore ghost points:
[10989]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)
[9031]169               END DO
[9788]170            ENDIF
171            !
172            IF( northern_side ) THEN
[9806]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               !
[9788]180               jbdy = nlcj-nbghostcells         
181               DO jn = 1, jptra
[10989]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)
[9788]183                  DO jk = 1, jpkm1
184                     DO ji = imin,imax
185                        IF( vmask(ji,jbdy-1,jk) == 0._wp ) THEN
[10989]186                           tr(ji,jbdy,jk,jn,Krhs) = tr(ji,jbdy+1,jk,jn,Krhs) * tmask(ji,jbdy,jk)
[9788]187                        ELSE
[10989]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)
[9788]192                           ENDIF
[5656]193                        ENDIF
[9019]194                     END DO
[5656]195                  END DO
[9788]196                  ! Restore ghost points:
[10989]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)
[636]198               END DO
[9788]199            ENDIF
200            !
[9806]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               !   
[9788]209               ibdy = 1+nbghostcells       
210               DO jn = 1, jptra
[10989]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)
[9788]212                  DO jk = 1, jpkm1
213                     DO jj = jmin,jmax
214                        IF( umask(ibdy,jj,jk) == 0._wp ) THEN
[10989]215                           tr(ibdy,jj,jk,jn,Krhs) = tr(ibdy-1,jj,jk,jn,Krhs) * tmask(ibdy,jj,jk)
[9788]216                        ELSE
[10989]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)
[9788]221                           ENDIF
[5656]222                        ENDIF
[9019]223                     END DO
[5656]224                  END DO
[9788]225                  ! Restore ghost points:
[10989]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)
[636]227               END DO
[9788]228            ENDIF
229            !
[9806]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               ! 
[9788]238               jbdy=1+nbghostcells       
239               DO jn = 1, jptra
[10989]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)
[9806]241                  DO jk = 1, jpkm1     
242                     DO ji = imin,imax
[9788]243                        IF( vmask(ji,jbdy,jk) == 0._wp ) THEN
[10989]244                           tr(ji,jbdy,jk,jn,Krhs)=tr(ji,jbdy-1,jk,jn,Krhs) * tmask(ji,jbdy,jk)
[9788]245                        ELSE
[10989]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)
[9788]250                           ENDIF
[5656]251                        ENDIF
[9019]252                     END DO
[5656]253                  END DO
[9788]254                  ! Restore ghost points:
[10989]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)
[636]256               END DO
[9019]257            ENDIF
258            !
[5656]259         ENDIF
[9806]260
[628]261      ENDIF
[3680]262      !
[5656]263   END SUBROUTINE interptrn
[2715]264
[628]265#else
[9019]266   !!----------------------------------------------------------------------
267   !!   Empty module                                           no TOP AGRIF
268   !!----------------------------------------------------------------------
[636]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
[628]276#endif
[9019]277
278   !!======================================================================
[636]279END MODULE agrif_top_interp
Note: See TracBrowser for help on using the repository browser.