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_opa_sponge.F90 in trunk/NEMOGCM/NEMO/NST_SRC – NEMO

source: trunk/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 3941

Last change on this file since 3941 was 3918, checked in by smasson, 11 years ago

trunk: better fortran and error in the conv of agrif, see ticket #1111 and #1112
  • Property svn:keywords set to Id
File size: 14.7 KB
RevLine 
[3680]1#define SPONGE && define SPONGE_TOP
[390]2
[636]3Module agrif_opa_sponge
[2528]4#if defined key_agrif  && ! defined key_offline
[636]5   USE par_oce
6   USE oce
7   USE dom_oce
8   USE in_out_manager
[782]9   USE agrif_oce
[3294]10   USE wrk_nemo 
[390]11
[636]12   IMPLICIT NONE
13   PRIVATE
[390]14
[3680]15   PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn, interptsn, interpun, interpvn
[390]16
[3680]17  !! * Substitutions
18#  include "domzgr_substitute.h90"
[1156]19   !!----------------------------------------------------------------------
[2528]20   !! NEMO/NST 3.3 , NEMO Consortium (2010)
[1156]21   !! $Id$
[2528]22   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[1156]23   !!----------------------------------------------------------------------
[390]24
[636]25   CONTAINS
26
[782]27   SUBROUTINE Agrif_Sponge_Tra
[636]28      !!---------------------------------------------
29      !!   *** ROUTINE Agrif_Sponge_Tra ***
30      !!---------------------------------------------
[2715]31      !!
[3294]32      INTEGER :: ji,jj,jk,jn
[636]33      REAL(wp) :: timecoeff
[3294]34      REAL(wp) :: ztsa, zabe1, zabe2, zbtr
35      REAL(wp), POINTER, DIMENSION(:,:    ) :: ztu, ztv
[3918]36      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztab
[3294]37      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: tsbdiff
[390]38
39#if defined SPONGE
[3680]40      CALL wrk_alloc( jpi, jpj, ztu, ztv )
[3918]41      CALL wrk_alloc( jpi, jpj, jpk, jpts, ztab, tsbdiff  )
[390]42
[636]43      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
44
[390]45      Agrif_SpecialValue=0.
46      Agrif_UseSpecialValue = .TRUE.
[3918]47      ztab = 0.e0
48      CALL Agrif_Bc_Variable(ztab, tsa_id,calledweight=timecoeff,procname=interptsn)
[390]49      Agrif_UseSpecialValue = .FALSE.
50
[3918]51      tsbdiff(:,:,:,:) = tsb(:,:,:,:) - ztab(:,:,:,:)
[390]52
[3680]53      CALL Agrif_Sponge
[390]54
[3294]55      DO jn = 1, jpts
56         DO jk = 1, jpkm1
57            !
58            DO jj = 1, jpjm1
59               DO ji = 1, jpim1
60                  zabe1 = umask(ji,jj,jk) * spe1ur(ji,jj) * fse3u(ji,jj,jk)
61                  zabe2 = vmask(ji,jj,jk) * spe2vr(ji,jj) * fse3v(ji,jj,jk)
62                  ztu(ji,jj) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) )
63                  ztv(ji,jj) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
64               ENDDO
[390]65            ENDDO
66
[3294]67            DO jj = 2, jpjm1
68               DO ji = 2, jpim1
69                  zbtr = spbtr2(ji,jj) / fse3t(ji,jj,jk)
70                  ! horizontal diffusive trends
71                  ztsa = zbtr * (  ztu(ji,jj) - ztu(ji-1,jj  )   &
72                  &              + ztv(ji,jj) - ztv(ji  ,jj-1)  )
73                  ! add it to the general tracer trends
74                  tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa
75               END DO
[390]76            END DO
[3294]77            !
78         ENDDO
[636]79      ENDDO
[390]80
[3680]81      CALL wrk_dealloc( jpi, jpj, ztu, ztv )
[3918]82      CALL wrk_dealloc( jpi, jpj, jpk, jpts, ztab, tsbdiff  )
[390]83#endif
84
[636]85   END SUBROUTINE Agrif_Sponge_Tra
[390]86
[782]87   SUBROUTINE Agrif_Sponge_dyn
[636]88      !!---------------------------------------------
89      !!   *** ROUTINE Agrif_Sponge_dyn ***
90      !!---------------------------------------------
[2715]91      !!
[390]92      INTEGER :: ji,jj,jk
[636]93      REAL(wp) :: timecoeff
[782]94      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
[2715]95      REAL(wp), POINTER, DIMENSION(:,:,:) :: ubdiff, vbdiff
96      REAL(wp), POINTER, DIMENSION(:,:,:) :: rotdiff, hdivdiff
[3918]97      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztab
[390]98
99#if defined SPONGE
[3918]100      CALL wrk_alloc( jpi, jpj, jpk, ztab, ubdiff, vbdiff, rotdiff, hdivdiff )
[390]101
[636]102      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
[390]103
104      Agrif_SpecialValue=0.
[782]105      Agrif_UseSpecialValue = ln_spc_dyn
[3918]106      ztab = 0.e0
107      CALL Agrif_Bc_Variable(ztab, ua_id,calledweight=timecoeff,procname=interpun)
[390]108      Agrif_UseSpecialValue = .FALSE.
109
[3918]110      ubdiff(:,:,:) = ( ub(:,:,:) - ztab(:,:,:) ) * umask(:,:,:)
[390]111
[3918]112      ztab = 0.e0
[390]113      Agrif_SpecialValue=0.
[782]114      Agrif_UseSpecialValue = ln_spc_dyn
[3918]115      CALL Agrif_Bc_Variable(ztab, va_id,calledweight=timecoeff,procname=interpvn)
[390]116      Agrif_UseSpecialValue = .FALSE.
117
[3918]118      vbdiff(:,:,:) = ( vb(:,:,:) - ztab(:,:,:) ) * vmask(:,:,:)
[390]119
[3680]120      CALL Agrif_Sponge
[390]121
[3680]122      DO jk = 1,jpkm1
123         ubdiff(:,:,jk) = ubdiff(:,:,jk) * spe1ur2(:,:)
124         vbdiff(:,:,jk) = vbdiff(:,:,jk) * spe2vr2(:,:)
[782]125      ENDDO
126     
[390]127      hdivdiff = 0.
128      rotdiff = 0.
129
130      DO jk = 1, jpkm1                                 ! Horizontal slab
131         !                                             ! ===============
132
133         !                                             ! --------
134         ! Horizontal divergence                       !   div
135         !                                             ! --------
136         DO jj = 2, jpjm1
137            DO ji = 2, jpim1   ! vector opt.
[782]138               zbtr = spbtr2(ji,jj) / fse3t(ji,jj,jk)
[3680]139               hdivdiff(ji,jj,jk) =  (  e2u(ji  ,jj  ) * fse3u(ji  ,jj  ,jk) * ubdiff(ji  ,jj  ,jk)     &
140                  &                   - e2u(ji-1,jj  ) * fse3u(ji-1,jj  ,jk) * ubdiff(ji-1,jj  ,jk)     &
141                  &                   + e1v(ji  ,jj  ) * fse3v(ji  ,jj  ,jk) * vbdiff(ji  ,jj  ,jk)     &
142                  &                   - e1v(ji  ,jj-1) * fse3v(ji  ,jj-1,jk) * vbdiff(ji  ,jj-1,jk)  ) * zbtr
[390]143            END DO
144         END DO
[636]145
[390]146         DO jj = 1, jpjm1
147            DO ji = 1, jpim1   ! vector opt.
[782]148               zbtr = spbtr3(ji,jj) * fse3f(ji,jj,jk)
149               rotdiff(ji,jj,jk) = (  e2v(ji+1,jj  ) * vbdiff(ji+1,jj  ,jk) - e2v(ji,jj) * vbdiff(ji,jj,jk)    &
[3680]150                  &                 - e1u(ji  ,jj+1) * ubdiff(ji  ,jj+1,jk) + e1u(ji,jj) * ubdiff(ji,jj,jk)  ) &
151                  &               * fmask(ji,jj,jk) * zbtr
[390]152            END DO
153         END DO
[636]154
155      ENDDO
156
[390]157      !                                                ! ===============
158      DO jk = 1, jpkm1                                 ! Horizontal slab
159         !                                             ! ===============
160         DO jj = 2, jpjm1
161            DO ji = 2, jpim1   ! vector opt.
162               ! horizontal diffusive trends
[3680]163               zua = - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   &
164                     + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj  ,jk) ) / e1u(ji,jj)
[390]165
[3680]166               zva = + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   &
167                     + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) / e2v(ji,jj)
[390]168               ! add it to the general momentum trends
169               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
170               va(ji,jj,jk) = va(ji,jj,jk) + zva
171            END DO
172         END DO
173         !                                             ! ===============
174      END DO                                           !   End of slab
175      !                                                ! ===============
[3918]176      CALL wrk_dealloc( jpi, jpj, jpk, ztab, ubdiff, vbdiff, rotdiff, hdivdiff )
[390]177#endif
178
[636]179   END SUBROUTINE Agrif_Sponge_dyn
[390]180
[3680]181   SUBROUTINE Agrif_Sponge
182      !!---------------------------------------------
183      !!   *** ROUTINE  Agrif_Sponge ***
184      !!---------------------------------------------
185      INTEGER  :: ji,jj,jk
186      INTEGER  :: ispongearea, ilci, ilcj
187      REAL(wp) :: z1spongearea
188      REAL(wp), POINTER, DIMENSION(:,:) :: zlocalviscsponge
189
190#if defined SPONGE || defined SPONGE_TOP
191
192      CALL wrk_alloc( jpi, jpj, zlocalviscsponge )
193
194      ispongearea  = 2 + 2 * Agrif_irhox()
195      ilci = nlci - ispongearea
196      ilcj = nlcj - ispongearea 
197      z1spongearea = 1._wp / REAL( ispongearea - 2 )
198      spbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
199
200      ! Tracers
201      IF( .NOT. spongedoneT ) THEN
202         zlocalviscsponge(:,:) = 0.
203         spe1ur(:,:) = 0.
204         spe2vr(:,:) = 0.
205
206         IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN
207            DO ji = 2, ispongearea
208               zlocalviscsponge(ji,:) = visc_tra * ( ispongearea-ji ) * z1spongearea
209            ENDDO
[3698]210            spe1ur(2:ispongearea-1,:      ) = 0.5 * ( zlocalviscsponge(2:ispongearea-1,:      )   &
211               &                         +            zlocalviscsponge(3:ispongearea  ,:      ) ) &
[3680]212               &                         * e2u(2:ispongearea-1,:      ) / e1u(2:ispongearea-1,:      )
[3698]213            spe2vr(2:ispongearea  ,1:jpjm1) = 0.5 * ( zlocalviscsponge(2:ispongearea  ,1:jpjm1)   &
214               &                         +            zlocalviscsponge(2:ispongearea,2  :jpj  ) ) &
[3680]215               &                         * e1v(2:ispongearea  ,1:jpjm1) / e2v(2:ispongearea  ,1:jpjm1)
216         ENDIF
217
218         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN
219            DO ji = ilci+1,nlci-1
220               zlocalviscsponge(ji,:) = visc_tra * (ji - (ilci+1) ) * z1spongearea
221            ENDDO
222 
[3698]223            spe1ur(ilci+1:nlci-2,:      ) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-2,:)    & 
224               &                          +          zlocalviscsponge(ilci+2:nlci-1,:) )  &
225               &                          * e2u(ilci+1:nlci-2,:) / e1u(ilci+1:nlci-2,:)
[3680]226
[3698]227            spe2vr(ilci+1:nlci-1,1:jpjm1) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-1,1:jpjm1)    & 
228               &                            +        zlocalviscsponge(ilci+1:nlci-1,2:jpj  )  ) & 
[3680]229               &                                   * e1v(ilci+1:nlci-1,1:jpjm1) / e2v(ilci+1:nlci-1,1:jpjm1)
230         ENDIF
231
232         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN
233            DO jj = 2, ispongearea
234               zlocalviscsponge(:,jj) = visc_tra * ( ispongearea-jj ) * z1spongearea
235            ENDDO
[3698]236            spe1ur(1:jpim1,2:ispongearea  ) = 0.5 * ( zlocalviscsponge(1:jpim1,2:ispongearea  ) & 
237               &                            +         zlocalviscsponge(2:jpi  ,2:ispongearea) ) &
238               &                            * e2u(1:jpim1,2:ispongearea) / e1u(1:jpim1,2:ispongearea)
[3680]239   
[3698]240            spe2vr(:      ,2:ispongearea-1) = 0.5 * ( zlocalviscsponge(:,2:ispongearea-1)       &
241               &                            +         zlocalviscsponge(:,3:ispongearea  )     ) &
242               &                            * e1v(:,2:ispongearea-1) / e2v(:,2:ispongearea-1)
[3680]243         ENDIF
244
245         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN
246            DO jj = ilcj+1,nlcj-1
247               zlocalviscsponge(:,jj) = visc_tra * (jj - (ilcj+1) ) * z1spongearea
248            ENDDO
[3698]249            spe1ur(1:jpim1,ilcj+1:nlcj-1) = 0.5 * ( zlocalviscsponge(1:jpim1,ilcj+1:nlcj-1)   &
250               &                          +         zlocalviscsponge(2:jpi  ,ilcj+1:nlcj-1) ) &
[3680]251               &                                * e2u(1:jpim1,ilcj+1:nlcj-1) / e1u(1:jpim1,ilcj+1:nlcj-1)
[3698]252            spe2vr(:      ,ilcj+1:nlcj-2) = 0.5 * ( zlocalviscsponge(:,ilcj+1:nlcj-2      )   &
253               &                          +         zlocalviscsponge(:,ilcj+2:nlcj-1)     )   &
[3680]254               &                                * e1v(:,ilcj+1:nlcj-2) / e2v(:,ilcj+1:nlcj-2)
255         ENDIF
256         spongedoneT = .TRUE.
257      ENDIF
258
259      ! Dynamics
260      IF( .NOT. spongedoneU ) THEN
261         zlocalviscsponge(:,:) = 0.
262         spe1ur2(:,:) = 0.
263         spe2vr2(:,:) = 0.
264
265         IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN
266            DO ji = 2, ispongearea
267               zlocalviscsponge(ji,:) = visc_dyn * ( ispongearea-ji ) * z1spongearea
268            ENDDO
[3698]269            spe1ur2(2:ispongearea-1,:      ) = 0.5 * ( zlocalviscsponge(2:ispongearea-1,:      ) &
270                                             &     +   zlocalviscsponge(3:ispongearea,:    ) )
271            spe2vr2(2:ispongearea  ,1:jpjm1) = 0.5 * ( zlocalviscsponge(2:ispongearea  ,1:jpjm1) &
272                                             &     +   zlocalviscsponge(2:ispongearea,2:jpj) ) 
[3680]273         ENDIF
274
275         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN
276            DO ji = ilci+1,nlci-1
277               zlocalviscsponge(ji,:) = visc_dyn * (ji - (ilci+1) ) * z1spongearea
278            ENDDO
[3698]279            spe1ur2(ilci+1:nlci-2,:      ) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-2,:) &
280                                           &        + zlocalviscsponge(ilci+2:nlci-1,:) ) 
281            spe2vr2(ilci+1:nlci-1,1:jpjm1) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-1,1:jpjm1) &
282                                           &        + zlocalviscsponge(ilci+1:nlci-1,2:jpj  )  ) 
[3680]283         ENDIF
284
285         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN
286            DO jj = 2, ispongearea
287               zlocalviscsponge(:,jj) = visc_dyn * ( ispongearea-jj ) * z1spongearea
288            ENDDO
[3698]289            spe1ur2(1:jpim1,2:ispongearea  ) = 0.5 * ( zlocalviscsponge(1:jpim1,2:ispongearea) &
290                                             &      + zlocalviscsponge(2:jpi,2:ispongearea) ) 
291            spe2vr2(:      ,2:ispongearea-1) = 0.5 * ( zlocalviscsponge(:,2:ispongearea-1)     &
292                                             &      + zlocalviscsponge(:,3:ispongearea)     )
[3680]293         ENDIF
294
295         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN
296            DO jj = ilcj+1,nlcj-1
297               zlocalviscsponge(:,jj) = visc_dyn * (jj - (ilcj+1) ) * z1spongearea
298            ENDDO
[3698]299            spe1ur2(1:jpim1,ilcj+1:nlcj-1) = 0.5 * ( zlocalviscsponge(1:jpim1,ilcj+1:nlcj-1) &
300                                           &         + zlocalviscsponge(2:jpi,ilcj+1:nlcj-1) ) 
301            spe2vr2(:      ,ilcj+1:nlcj-2) = 0.5 * ( zlocalviscsponge(:,ilcj+1:nlcj-2      ) &
302                                           &         + zlocalviscsponge(:,ilcj+2:nlcj-1)     )
[3680]303         ENDIF
304         spongedoneU = .TRUE.
305         spbtr3(:,:) = 1. / ( e1f(:,:) * e2f(:,:) )
306      ENDIF
307      !
308      CALL wrk_dealloc( jpi, jpj, zlocalviscsponge )
309      !
310#endif
311
312   END SUBROUTINE Agrif_Sponge
313
[3294]314   SUBROUTINE interptsn(tabres,i1,i2,j1,j2,k1,k2,n1,n2)
[636]315      !!---------------------------------------------
[3294]316      !!   *** ROUTINE interptsn ***
[636]317      !!---------------------------------------------
[3294]318      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
319      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
[390]320
[3294]321      tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
[390]322
[3294]323   END SUBROUTINE interptsn
[636]324
325   SUBROUTINE interpun(tabres,i1,i2,j1,j2,k1,k2)
326      !!---------------------------------------------
327      !!   *** ROUTINE interpun ***
328      !!---------------------------------------------
329      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
330      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
[390]331
[636]332      tabres(i1:i2,j1:j2,k1:k2) = un(i1:i2,j1:j2,k1:k2)
[390]333
[636]334   END SUBROUTINE interpun
335
336   SUBROUTINE interpvn(tabres,i1,i2,j1,j2,k1,k2)
337      !!---------------------------------------------
338      !!   *** ROUTINE interpvn ***
339      !!---------------------------------------------
340      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
341      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
[390]342
[636]343      tabres(i1:i2,j1:j2,k1:k2) = vn(i1:i2,j1:j2,k1:k2)
[390]344
[636]345   END SUBROUTINE interpvn
[390]346
347#else
[636]348CONTAINS
349
350   SUBROUTINE agrif_opa_sponge_empty
351      !!---------------------------------------------
352      !!   *** ROUTINE agrif_OPA_sponge_empty ***
353      !!---------------------------------------------
354      WRITE(*,*)  'agrif_opa_sponge : You should not have seen this print! error?'
355   END SUBROUTINE agrif_opa_sponge_empty
[390]356#endif
[636]357
358END MODULE agrif_opa_sponge
Note: See TracBrowser for help on using the repository browser.