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 branches/2012/dev_r3387_LOCEAN6_AGRIF_LIM/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2012/dev_r3387_LOCEAN6_AGRIF_LIM/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 4736

Last change on this file since 4736 was 3566, checked in by cetlod, 12 years ago

branch dev_r3387_LOCEAN6_AGRIF_LIM: add some corrections to make AGRIF compatible with TOP

  • Property svn:keywords set to Id
File size: 13.9 KB
RevLine 
[3566]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
[3566]15   PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn, interptsn, interpun, interpvn
[390]16
[3566]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
36      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztab
37      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: tsbdiff
[390]38
39#if defined SPONGE
[3566]40      CALL wrk_alloc( jpi, jpj, ztu, ztv )
[3294]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.
[636]47      ztab = 0.e0
[3294]48      CALL Agrif_Bc_Variable(ztab, tsa_id,calledweight=timecoeff,procname=interptsn)
[390]49      Agrif_UseSpecialValue = .FALSE.
50
[3294]51      tsbdiff(:,:,:,:) = tsb(:,:,:,:) - ztab(:,:,:,:)
[390]52
[3566]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
[3566]81      CALL wrk_dealloc( jpi, jpj, ztu, ztv )
[3294]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
97      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztab
[390]98
99#if defined SPONGE
[3294]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
[636]106      ztab = 0.e0
[2715]107      CALL Agrif_Bc_Variable(ztab, ua_id,calledweight=timecoeff,procname=interpun)
[390]108      Agrif_UseSpecialValue = .FALSE.
109
[3566]110      ubdiff(:,:,:) = ( ub(:,:,:) - ztab(:,:,:) ) * umask(:,:,:)
[390]111
[636]112      ztab = 0.e0
[390]113      Agrif_SpecialValue=0.
[782]114      Agrif_UseSpecialValue = ln_spc_dyn
[2715]115      CALL Agrif_Bc_Variable(ztab, va_id,calledweight=timecoeff,procname=interpvn)
[390]116      Agrif_UseSpecialValue = .FALSE.
117
[3566]118      vbdiff(:,:,:) = ( vb(:,:,:) - ztab(:,:,:) ) * vmask(:,:,:)
[390]119
[3566]120      CALL Agrif_Sponge
[390]121
[3566]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)
[3566]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)    &
[3566]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
[3566]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
[3566]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      !                                                ! ===============
[3294]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
[3566]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
210            spe1ur(2:ispongearea-1,:      ) = 0.5 * ( zlocalviscsponge(2:ispongearea-1,:      ) + zlocalviscsponge(3:ispongearea,:    ) ) &
211               &                         * e2u(2:ispongearea-1,:      ) / e1u(2:ispongearea-1,:      )
212            spe2vr(2:ispongearea  ,1:jpjm1) = 0.5 * ( zlocalviscsponge(2:ispongearea  ,1:jpjm1) + zlocalviscsponge(2:ispongearea,2:jpj) ) &
213               &                         * e1v(2:ispongearea  ,1:jpjm1) / e2v(2:ispongearea  ,1:jpjm1)
214         ENDIF
215
216         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN
217            DO ji = ilci+1,nlci-1
218               zlocalviscsponge(ji,:) = visc_tra * (ji - (ilci+1) ) * z1spongearea
219            ENDDO
220 
221            spe1ur(ilci+1:nlci-2,:      ) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-2,:) + zlocalviscsponge(ilci+2:nlci-1,:) )  &
222               &                                   * e2u(ilci+1:nlci-2,:) / e1u(ilci+1:nlci-2,:)
223
224            spe2vr(ilci+1:nlci-1,1:jpjm1) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-1,1:jpjm1) + zlocalviscsponge(ilci+1:nlci-1,2:jpj  )  ) & 
225               &                                   * e1v(ilci+1:nlci-1,1:jpjm1) / e2v(ilci+1:nlci-1,1:jpjm1)
226         ENDIF
227
228         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN
229            DO jj = 2, ispongearea
230               zlocalviscsponge(:,jj) = visc_tra * ( ispongearea-jj ) * z1spongearea
231            ENDDO
232            spe1ur(1:jpim1,2:ispongearea  ) = 0.5 * ( zlocalviscsponge(1:jpim1,2:ispongearea) + zlocalviscsponge(2:jpi,2:ispongearea) ) &
233               &                                * e2u(1:jpim1,2:ispongearea) / e1u(1:jpim1,2:ispongearea)
234   
235            spe2vr(:      ,2:ispongearea-1) = 0.5 * ( zlocalviscsponge(:,2:ispongearea-1)     + zlocalviscsponge(:,3:ispongearea)     ) &
236               &                                  * e1v(:,2:ispongearea-1) / e2v(:,2:ispongearea-1)
237         ENDIF
238
239         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN
240            DO jj = ilcj+1,nlcj-1
241               zlocalviscsponge(:,jj) = visc_tra * (jj - (ilcj+1) ) * z1spongearea
242            ENDDO
243            spe1ur(1:jpim1,ilcj+1:nlcj-1) = 0.5 * ( zlocalviscsponge(1:jpim1,ilcj+1:nlcj-1) + zlocalviscsponge(2:jpi,ilcj+1:nlcj-1) ) &
244               &                                * e2u(1:jpim1,ilcj+1:nlcj-1) / e1u(1:jpim1,ilcj+1:nlcj-1)
245            spe2vr(:      ,ilcj+1:nlcj-2) = 0.5 * ( zlocalviscsponge(:,ilcj+1:nlcj-2      ) + zlocalviscsponge(:,ilcj+2:nlcj-1)     ) &
246               &                                * e1v(:,ilcj+1:nlcj-2) / e2v(:,ilcj+1:nlcj-2)
247         ENDIF
248         spongedoneT = .TRUE.
249      ENDIF
250
251      ! Dynamics
252      IF( .NOT. spongedoneU ) THEN
253         zlocalviscsponge(:,:) = 0.
254         spe1ur2(:,:) = 0.
255         spe2vr2(:,:) = 0.
256
257         IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN
258            DO ji = 2, ispongearea
259               zlocalviscsponge(ji,:) = visc_dyn * ( ispongearea-ji ) * z1spongearea
260            ENDDO
261            spe1ur2(2:ispongearea-1,:      ) = 0.5 * ( zlocalviscsponge(2:ispongearea-1,:      ) + zlocalviscsponge(3:ispongearea,:    ) )
262            spe2vr2(2:ispongearea  ,1:jpjm1) = 0.5 * ( zlocalviscsponge(2:ispongearea  ,1:jpjm1) + zlocalviscsponge(2:ispongearea,2:jpj) ) 
263         ENDIF
264
265         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN
266            DO ji = ilci+1,nlci-1
267               zlocalviscsponge(ji,:) = visc_dyn * (ji - (ilci+1) ) * z1spongearea
268            ENDDO
269            spe1ur2(ilci+1:nlci-2,:      ) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-2,:) + zlocalviscsponge(ilci+2:nlci-1,:) ) 
270            spe2vr2(ilci+1:nlci-1,1:jpjm1) = 0.5 * (  zlocalviscsponge(ilci+1:nlci-1,1:jpjm1) + zlocalviscsponge(ilci+1:nlci-1,2:jpj  )  ) 
271         ENDIF
272
273         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN
274            DO jj = 2, ispongearea
275               zlocalviscsponge(:,jj) = visc_dyn * ( ispongearea-jj ) * z1spongearea
276            ENDDO
277            spe1ur2(1:jpim1,2:ispongearea  ) = 0.5 * ( zlocalviscsponge(1:jpim1,2:ispongearea) + zlocalviscsponge(2:jpi,2:ispongearea) ) 
278            spe2vr2(:      ,2:ispongearea-1) = 0.5 * ( zlocalviscsponge(:,2:ispongearea-1)     + zlocalviscsponge(:,3:ispongearea)     )
279         ENDIF
280
281         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN
282            DO jj = ilcj+1,nlcj-1
283               zlocalviscsponge(:,jj) = visc_dyn * (jj - (ilcj+1) ) * z1spongearea
284            ENDDO
285            spe1ur2(1:jpim1,ilcj+1:nlcj-1) = 0.5 * ( zlocalviscsponge(1:jpim1,ilcj+1:nlcj-1) + zlocalviscsponge(2:jpi,ilcj+1:nlcj-1) ) 
286            spe2vr2(:      ,ilcj+1:nlcj-2) = 0.5 * ( zlocalviscsponge(:,ilcj+1:nlcj-2      ) + zlocalviscsponge(:,ilcj+2:nlcj-1)     )
287         ENDIF
288         spongedoneU = .TRUE.
289         spbtr3(:,:) = 1. / ( e1f(:,:) * e2f(:,:) )
290      ENDIF
291      !
292      CALL wrk_dealloc( jpi, jpj, zlocalviscsponge )
293      !
294#endif
295
296   END SUBROUTINE Agrif_Sponge
297
[3294]298   SUBROUTINE interptsn(tabres,i1,i2,j1,j2,k1,k2,n1,n2)
[636]299      !!---------------------------------------------
[3294]300      !!   *** ROUTINE interptsn ***
[636]301      !!---------------------------------------------
[3294]302      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
303      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
[390]304
[3294]305      tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
[390]306
[3294]307   END SUBROUTINE interptsn
[636]308
309   SUBROUTINE interpun(tabres,i1,i2,j1,j2,k1,k2)
310      !!---------------------------------------------
311      !!   *** ROUTINE interpun ***
312      !!---------------------------------------------
313      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
314      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
[390]315
[636]316      tabres(i1:i2,j1:j2,k1:k2) = un(i1:i2,j1:j2,k1:k2)
[390]317
[636]318   END SUBROUTINE interpun
319
320   SUBROUTINE interpvn(tabres,i1,i2,j1,j2,k1,k2)
321      !!---------------------------------------------
322      !!   *** ROUTINE interpvn ***
323      !!---------------------------------------------
324      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
325      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
[390]326
[636]327      tabres(i1:i2,j1:j2,k1:k2) = vn(i1:i2,j1:j2,k1:k2)
[390]328
[636]329   END SUBROUTINE interpvn
[390]330
331#else
[636]332CONTAINS
333
334   SUBROUTINE agrif_opa_sponge_empty
335      !!---------------------------------------------
336      !!   *** ROUTINE agrif_OPA_sponge_empty ***
337      !!---------------------------------------------
338      WRITE(*,*)  'agrif_opa_sponge : You should not have seen this print! error?'
339   END SUBROUTINE agrif_opa_sponge_empty
[390]340#endif
[636]341
342END MODULE agrif_opa_sponge
Note: See TracBrowser for help on using the repository browser.