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_MERGE_2012/NEMOGCM/NEMO/NST_SRC – NEMO

source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 3680

Last change on this file since 3680 was 3680, checked in by rblod, 11 years ago

First commit of the final branch for 2012 (future nemo_3_5), see ticket #1028
  • Property svn:keywords set to Id
File size: 13.9 KB
Line 
1#define SPONGE && define SPONGE_TOP
2
3Module agrif_opa_sponge
4#if defined key_agrif  && ! defined key_offline
5   USE par_oce
6   USE oce
7   USE dom_oce
8   USE in_out_manager
9   USE agrif_oce
10   USE wrk_nemo 
11
12   IMPLICIT NONE
13   PRIVATE
14
15   PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn, interptsn, interpun, interpvn
16
17  !! * Substitutions
18#  include "domzgr_substitute.h90"
19   !!----------------------------------------------------------------------
20   !! NEMO/NST 3.3 , NEMO Consortium (2010)
21   !! $Id$
22   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
23   !!----------------------------------------------------------------------
24
25   CONTAINS
26
27   SUBROUTINE Agrif_Sponge_Tra
28      !!---------------------------------------------
29      !!   *** ROUTINE Agrif_Sponge_Tra ***
30      !!---------------------------------------------
31      !!
32      INTEGER :: ji,jj,jk,jn
33      REAL(wp) :: timecoeff
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
38
39#if defined SPONGE
40      CALL wrk_alloc( jpi, jpj, ztu, ztv )
41      CALL wrk_alloc( jpi, jpj, jpk, jpts, ztab, tsbdiff  )
42
43      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
44
45      Agrif_SpecialValue=0.
46      Agrif_UseSpecialValue = .TRUE.
47      ztab = 0.e0
48      CALL Agrif_Bc_Variable(ztab, tsa_id,calledweight=timecoeff,procname=interptsn)
49      Agrif_UseSpecialValue = .FALSE.
50
51      tsbdiff(:,:,:,:) = tsb(:,:,:,:) - ztab(:,:,:,:)
52
53      CALL Agrif_Sponge
54
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
65            ENDDO
66
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
76            END DO
77            !
78         ENDDO
79      ENDDO
80
81      CALL wrk_dealloc( jpi, jpj, ztu, ztv )
82      CALL wrk_dealloc( jpi, jpj, jpk, jpts, ztab, tsbdiff  )
83#endif
84
85   END SUBROUTINE Agrif_Sponge_Tra
86
87   SUBROUTINE Agrif_Sponge_dyn
88      !!---------------------------------------------
89      !!   *** ROUTINE Agrif_Sponge_dyn ***
90      !!---------------------------------------------
91      !!
92      INTEGER :: ji,jj,jk
93      REAL(wp) :: timecoeff
94      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
95      REAL(wp), POINTER, DIMENSION(:,:,:) :: ubdiff, vbdiff
96      REAL(wp), POINTER, DIMENSION(:,:,:) :: rotdiff, hdivdiff
97      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztab
98
99#if defined SPONGE
100      CALL wrk_alloc( jpi, jpj, jpk, ztab, ubdiff, vbdiff, rotdiff, hdivdiff )
101
102      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
103
104      Agrif_SpecialValue=0.
105      Agrif_UseSpecialValue = ln_spc_dyn
106      ztab = 0.e0
107      CALL Agrif_Bc_Variable(ztab, ua_id,calledweight=timecoeff,procname=interpun)
108      Agrif_UseSpecialValue = .FALSE.
109
110      ubdiff(:,:,:) = ( ub(:,:,:) - ztab(:,:,:) ) * umask(:,:,:)
111
112      ztab = 0.e0
113      Agrif_SpecialValue=0.
114      Agrif_UseSpecialValue = ln_spc_dyn
115      CALL Agrif_Bc_Variable(ztab, va_id,calledweight=timecoeff,procname=interpvn)
116      Agrif_UseSpecialValue = .FALSE.
117
118      vbdiff(:,:,:) = ( vb(:,:,:) - ztab(:,:,:) ) * vmask(:,:,:)
119
120      CALL Agrif_Sponge
121
122      DO jk = 1,jpkm1
123         ubdiff(:,:,jk) = ubdiff(:,:,jk) * spe1ur2(:,:)
124         vbdiff(:,:,jk) = vbdiff(:,:,jk) * spe2vr2(:,:)
125      ENDDO
126     
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.
138               zbtr = spbtr2(ji,jj) / fse3t(ji,jj,jk)
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
143            END DO
144         END DO
145
146         DO jj = 1, jpjm1
147            DO ji = 1, jpim1   ! vector opt.
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)    &
150                  &                 - e1u(ji  ,jj+1) * ubdiff(ji  ,jj+1,jk) + e1u(ji,jj) * ubdiff(ji,jj,jk)  ) &
151                  &               * fmask(ji,jj,jk) * zbtr
152            END DO
153         END DO
154
155      ENDDO
156
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
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)
165
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)
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      !                                                ! ===============
176      CALL wrk_dealloc( jpi, jpj, jpk, ztab, ubdiff, vbdiff, rotdiff, hdivdiff )
177#endif
178
179   END SUBROUTINE Agrif_Sponge_dyn
180
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
298   SUBROUTINE interptsn(tabres,i1,i2,j1,j2,k1,k2,n1,n2)
299      !!---------------------------------------------
300      !!   *** ROUTINE interptsn ***
301      !!---------------------------------------------
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
304
305      tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
306
307   END SUBROUTINE interptsn
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
315
316      tabres(i1:i2,j1:j2,k1:k2) = un(i1:i2,j1:j2,k1:k2)
317
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
326
327      tabres(i1:i2,j1:j2,k1:k2) = vn(i1:i2,j1:j2,k1:k2)
328
329   END SUBROUTINE interpvn
330
331#else
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
340#endif
341
342END MODULE agrif_opa_sponge
Note: See TracBrowser for help on using the repository browser.