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 tags/nemo_v3_2/nemo_v3_2/NEMO/NST_SRC – NEMO

source: tags/nemo_v3_2/nemo_v3_2/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 1878

Last change on this file since 1878 was 1878, checked in by flavoni, 14 years ago

initial test for nemogcm
File size: 15.4 KB
Line 
1#define SPONGE
2
3Module agrif_opa_sponge
4#if defined key_agrif  && ! defined key_off_tra
5   USE par_oce
6   USE oce
7   USE dom_oce
8   USE in_out_manager
9   USE agrif_oce
10
11   IMPLICIT NONE
12   PRIVATE
13
14   PUBLIC Agrif_Sponge_Tra, Agrif_Sponge_Dyn, interptn, interpsn, interpun, interpvn
15
16   !!----------------------------------------------------------------------
17   !!   OPA 9.0 , LOCEAN-IPSL (2006)
18   !! $Id: agrif_opa_sponge.F90 1300 2009-02-09 15:36:04Z rblod $
19   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
20   !!----------------------------------------------------------------------
21
22   CONTAINS
23
24   SUBROUTINE Agrif_Sponge_Tra
25      !!---------------------------------------------
26      !!   *** ROUTINE Agrif_Sponge_Tra ***
27      !!---------------------------------------------
28#include "domzgr_substitute.h90"
29
30      INTEGER :: ji,jj,jk
31      INTEGER :: spongearea
32      REAL(wp) :: timecoeff
33      REAL(wp) :: zta, zsa, zabe1, zabe2, zbtr
34      REAL(wp), DIMENSION(jpi,jpj) :: localviscsponge
35      REAL(wp), DIMENSION(jpi,jpj,jpk) :: tbdiff, sbdiff
36      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztu ,ztv, zsu ,zsv
37      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  ztab
38
39#if defined SPONGE
40
41      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
42
43      Agrif_SpecialValue=0.
44      Agrif_UseSpecialValue = .TRUE.
45      ztab = 0.e0
46      CALL Agrif_Bc_Variable(ztab, ta,calledweight=timecoeff,procname=interptn)
47      Agrif_UseSpecialValue = .FALSE.
48
49      tbdiff(:,:,:) = tb(:,:,:) - ztab(:,:,:)
50
51      ztab = 0.e0
52      Agrif_SpecialValue=0.
53      Agrif_UseSpecialValue = .TRUE.
54      CALL Agrif_Bc_Variable(ztab, sa,calledweight=timecoeff,procname=interpsn)
55      Agrif_UseSpecialValue = .FALSE.
56
57      sbdiff(:,:,:) = sb(:,:,:) - ztab(:,:,:)
58
59
60      spongearea = 2 + 2 * Agrif_irhox()
61
62      localviscsponge = 0.
63     
64      IF (.NOT. spongedoneT) THEN
65         spe1ur(:,:) = 0.
66         spe2vr(:,:) = 0.
67
68      IF ((nbondi == -1).OR.(nbondi == 2)) THEN
69         DO ji = 2, spongearea
70            localviscsponge(ji,:) = visc_tra * (spongearea-ji)/real(spongearea-2)
71         ENDDO
72   
73    spe1ur(2:spongearea-1,:)=0.5 * (localviscsponge(2:spongearea-1,:) + localviscsponge(3:spongearea,:)) &
74          * e2u(2:spongearea-1,:) / e1u(2:spongearea-1,:)
75
76         spe2vr(2:spongearea,1:jpjm1) = 0.5 * (localviscsponge(2:spongearea,1:jpjm1) + &
77             localviscsponge(2:spongearea,2:jpj)) &
78           * e1v(2:spongearea,1:jpjm1) / e2v(2:spongearea,1:jpjm1)
79      ENDIF
80
81      IF ((nbondi == 1).OR.(nbondi == 2)) THEN
82         DO ji = nlci-spongearea + 1,nlci-1
83            localviscsponge(ji,:) = visc_tra * (ji - (nlci-spongearea+1))/real(spongearea-2)
84         ENDDO
85   
86    spe1ur(nlci-spongearea + 1:nlci-2,:)=0.5 * (localviscsponge(nlci-spongearea + 1:nlci-2,:) + &
87           localviscsponge(nlci-spongearea + 2:nlci-1,:)) &
88          * e2u(nlci-spongearea + 1:nlci-2,:) / e1u(nlci-spongearea + 1:nlci-2,:)
89
90         spe2vr(nlci-spongearea + 1:nlci-1,1:jpjm1) = 0.5 * (localviscsponge(nlci-spongearea + 1:nlci-1,1:jpjm1) &
91              + localviscsponge(nlci-spongearea + 1:nlci-1,2:jpj)) &
92           * e1v(nlci-spongearea + 1:nlci-1,1:jpjm1) / e2v(nlci-spongearea + 1:nlci-1,1:jpjm1)
93      ENDIF
94
95
96      IF ((nbondj == -1).OR.(nbondj == 2)) THEN
97         DO jj = 2, spongearea
98            localviscsponge(:,jj) = visc_tra * (spongearea-jj)/real(spongearea-2)
99         ENDDO
100   
101    spe1ur(1:jpim1,2:spongearea)=0.5 * (localviscsponge(1:jpim1,2:spongearea) + &
102           localviscsponge(2:jpi,2:spongearea)) &
103          * e2u(1:jpim1,2:spongearea) / e1u(1:jpim1,2:spongearea)
104
105         spe2vr(:,2:spongearea-1) = 0.5 * (localviscsponge(:,2:spongearea-1) + &
106             localviscsponge(:,3:spongearea)) &
107           * e1v(:,2:spongearea-1) / e2v(:,2:spongearea-1)
108      ENDIF
109
110      IF ((nbondj == 1).OR.(nbondj == 2)) THEN
111         DO jj = nlcj-spongearea + 1,nlcj-1
112            localviscsponge(:,jj) = visc_tra * (jj - (nlcj-spongearea+1))/real(spongearea-2)
113         ENDDO
114   
115    spe1ur(1:jpim1,nlcj-spongearea + 1:nlcj-1)=0.5 * (localviscsponge(1:jpim1,nlcj-spongearea + 1:nlcj-1) + &
116            localviscsponge(2:jpi,nlcj-spongearea + 1:nlcj-1)) &
117          * e2u(1:jpim1,nlcj-spongearea + 1:nlcj-1) / e1u(1:jpim1,nlcj-spongearea + 1:nlcj-1)
118
119         spe2vr(:,nlcj-spongearea + 1:nlcj-2) = 0.5 * (localviscsponge(:,nlcj-spongearea + 1:nlcj-2) + &
120            localviscsponge(:,nlcj-spongearea + 2:nlcj-1)) &
121           * e1v(:,nlcj-spongearea + 1:nlcj-2) / e2v(:,nlcj-spongearea + 1:nlcj-2)
122      ENDIF
123     
124         spbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:))
125
126         spongedoneT = .TRUE.
127      ENDIF
128
129      DO jk = 1, jpkm1
130         DO jj = 1, jpjm1
131            DO ji = 1, jpim1
132#if defined key_zco
133               zabe1 = umask(ji,jj,jk) * spe1ur(ji,jj)
134               zabe2 = vmask(ji,jj,jk) * spe2vr(ji,jj)
135#else
136               zabe1 = umask(ji,jj,jk) * spe1ur(ji,jj) * fse3u(ji,jj,jk)
137               zabe2 = vmask(ji,jj,jk) * spe2vr(ji,jj) * fse3v(ji,jj,jk)
138#endif
139               ztu(ji,jj,jk) = zabe1 * ( tbdiff(ji+1,jj  ,jk) - tbdiff(ji,jj,jk) )
140               zsu(ji,jj,jk) = zabe1 * ( sbdiff(ji+1,jj  ,jk) - sbdiff(ji,jj,jk) )
141               ztv(ji,jj,jk) = zabe2 * ( tbdiff(ji  ,jj+1,jk) - tbdiff(ji,jj,jk) )
142               zsv(ji,jj,jk) = zabe2 * ( sbdiff(ji  ,jj+1,jk) - sbdiff(ji,jj,jk) )
143            ENDDO
144         ENDDO
145
146         DO jj = 2,jpjm1
147            DO ji = 2,jpim1
148#if defined key_zco
149               zbtr = spbtr2(ji,jj)
150#else
151               zbtr = spbtr2(ji,jj) / fse3t(ji,jj,jk)
152#endif
153               ! horizontal diffusive trends
154               zta = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   &
155                  &          + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  )
156               zsa = zbtr * (  zsu(ji,jj,jk) - zsu(ji-1,jj,jk)   &
157                  &          + zsv(ji,jj,jk) - zsv(ji,jj-1,jk)  )
158               ! add it to the general tracer trends
159               ta(ji,jj,jk) = (ta(ji,jj,jk) + zta)
160               sa(ji,jj,jk) = (sa(ji,jj,jk) + zsa)
161            END DO
162         END DO
163
164      ENDDO
165
166#endif
167
168   END SUBROUTINE Agrif_Sponge_Tra
169
170   SUBROUTINE Agrif_Sponge_dyn
171      !!---------------------------------------------
172      !!   *** ROUTINE Agrif_Sponge_dyn ***
173      !!---------------------------------------------
174#include "domzgr_substitute.h90"
175
176      INTEGER :: ji,jj,jk
177      INTEGER :: spongearea
178      REAL(wp) :: timecoeff
179      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
180      REAL(wp), DIMENSION(jpi,jpj) :: localviscsponge
181      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztab, ubdiff, vbdiff,rotdiff,hdivdiff
182
183#if defined SPONGE
184
185      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
186
187      Agrif_SpecialValue=0.
188      Agrif_UseSpecialValue = ln_spc_dyn
189      ztab = 0.e0
190      CALL Agrif_Bc_Variable(ztab, ua,calledweight=timecoeff,procname=interpun)
191      Agrif_UseSpecialValue = .FALSE.
192
193      ubdiff(:,:,:) = (ub(:,:,:) - ztab(:,:,:))*umask(:,:,:)
194
195      ztab = 0.e0
196      Agrif_SpecialValue=0.
197      Agrif_UseSpecialValue = ln_spc_dyn
198      CALL Agrif_Bc_Variable(ztab, va,calledweight=timecoeff,procname=interpvn)
199      Agrif_UseSpecialValue = .FALSE.
200
201      vbdiff(:,:,:) = (vb(:,:,:) - ztab(:,:,:))*vmask(:,:,:)
202
203      spongearea = 2 + 2 * Agrif_irhox()
204
205      localviscsponge = 0.
206
207      IF (.NOT. spongedoneU) THEN
208         spe1ur2(:,:) = 0.
209         spe2vr2(:,:) = 0.
210
211      IF ((nbondi == -1).OR.(nbondi == 2)) THEN
212         DO ji = 2, spongearea
213            localviscsponge(ji,:) = visc_dyn * (spongearea-ji)/real(spongearea-2)
214         ENDDO
215   
216    spe1ur2(2:spongearea-1,:)=0.5 * (localviscsponge(2:spongearea-1,:) + localviscsponge(3:spongearea,:))
217
218         spe2vr2(2:spongearea,1:jpjm1) = 0.5 * (localviscsponge(2:spongearea,1:jpjm1) + &
219             localviscsponge(2:spongearea,2:jpj))
220      ENDIF
221
222      IF ((nbondi == 1).OR.(nbondi == 2)) THEN
223         DO ji = nlci-spongearea + 1,nlci-1
224            localviscsponge(ji,:) = visc_dyn * (ji - (nlci-spongearea+1))/real(spongearea-2)
225         ENDDO
226   
227    spe1ur2(nlci-spongearea + 1:nlci-2,:)=0.5 * (localviscsponge(nlci-spongearea + 1:nlci-2,:) + &
228           localviscsponge(nlci-spongearea + 2:nlci-1,:))
229
230         spe2vr2(nlci-spongearea + 1:nlci-1,1:jpjm1) = 0.5 * (localviscsponge(nlci-spongearea + 1:nlci-1,1:jpjm1) &
231              + localviscsponge(nlci-spongearea + 1:nlci-1,2:jpj))
232      ENDIF
233
234
235      IF ((nbondj == -1).OR.(nbondj == 2)) THEN
236         DO jj = 2, spongearea
237            localviscsponge(:,jj) = visc_dyn * (spongearea-jj)/real(spongearea-2)
238         ENDDO
239   
240    spe1ur2(1:jpim1,2:spongearea)=0.5 * (localviscsponge(1:jpim1,2:spongearea) + &
241           localviscsponge(2:jpi,2:spongearea))
242
243         spe2vr2(:,2:spongearea-1) = 0.5 * (localviscsponge(:,2:spongearea-1) + &
244             localviscsponge(:,3:spongearea))
245      ENDIF
246
247      IF ((nbondj == 1).OR.(nbondj == 2)) THEN
248         DO jj = nlcj-spongearea + 1,nlcj-1
249            localviscsponge(:,jj) = visc_dyn * (jj - (nlcj-spongearea+1))/real(spongearea-2)
250         ENDDO
251   
252    spe1ur2(1:jpim1,nlcj-spongearea + 1:nlcj-1)=0.5 * (localviscsponge(1:jpim1,nlcj-spongearea + 1:nlcj-1) + &
253            localviscsponge(2:jpi,nlcj-spongearea + 1:nlcj-1))
254
255         spe2vr2(:,nlcj-spongearea + 1:nlcj-2) = 0.5 * (localviscsponge(:,nlcj-spongearea + 1:nlcj-2) + &
256            localviscsponge(:,nlcj-spongearea + 2:nlcj-1))
257      ENDIF
258
259         spongedoneU = .TRUE.
260   
261    spbtr3(:,:) = 1./( e1f(:,:) * e2f(:,:))
262      ENDIF
263     
264      IF (.NOT. spongedoneT) THEN
265        spbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:))     
266      ENDIF
267     
268      DO jk=1,jpkm1
269      ubdiff(:,:,jk) = ubdiff(:,:,jk) * spe1ur2(:,:)
270      vbdiff(:,:,jk) = vbdiff(:,:,jk) * spe2vr2(:,:)
271      ENDDO
272     
273      hdivdiff = 0.
274      rotdiff = 0.
275
276      DO jk = 1, jpkm1                                 ! Horizontal slab
277         !                                             ! ===============
278
279         !                                             ! --------
280         ! Horizontal divergence                       !   div
281         !                                             ! --------
282         DO jj = 2, jpjm1
283            DO ji = 2, jpim1   ! vector opt.
284#if defined key_zco
285               zbtr = spbtr2(ji,jj)
286               hdivdiff(ji,jj,jk) = (  e2u(ji,jj) * ubdiff(ji,jj,jk) &
287                  - e2u(ji-1,jj  ) * ubdiff(ji-1,jj  ,jk)      &
288                  &               + e1v(ji,jj) * vbdiff(ji,jj,jk) - &
289                  &              e1v(ji  ,jj-1) * vbdiff(ji  ,jj-1,jk)  ) * zbtr
290#else
291               zbtr = spbtr2(ji,jj) / fse3t(ji,jj,jk)
292               hdivdiff(ji,jj,jk) =   &
293                  (  e2u(ji,jj)*fse3u(ji,jj,jk) * & 
294                  ubdiff(ji,jj,jk) - e2u(ji-1,jj  )* &
295                  fse3u(ji-1,jj  ,jk)  * ubdiff(ji-1,jj  ,jk)       &
296                  + e1v(ji,jj)*fse3v(ji,jj,jk) * &
297                  vbdiff(ji,jj,jk) - e1v(ji  ,jj-1)* &
298                  fse3v(ji  ,jj-1,jk)  * vbdiff(ji  ,jj-1,jk)  ) * zbtr
299#endif
300            END DO
301         END DO
302
303         DO jj = 1, jpjm1
304            DO ji = 1, jpim1   ! vector opt.
305#if defined key_zco     
306               zbtr = spbtr3(ji,jj)
307               rotdiff(ji,jj,jk) = (  e2v(ji+1,jj  ) * vbdiff(ji+1,jj  ,jk) - e2v(ji,jj) * vbdiff(ji,jj,jk)    &
308                  &              - e1u(ji  ,jj+1) * ubdiff(ji  ,jj+1,jk) + e1u(ji,jj) * ubdiff(ji,jj,jk)  ) &
309                  &           * fmask(ji,jj,jk) * zbtr
310#else
311               zbtr = spbtr3(ji,jj) * fse3f(ji,jj,jk)
312               rotdiff(ji,jj,jk) = (  e2v(ji+1,jj  ) * vbdiff(ji+1,jj  ,jk) - e2v(ji,jj) * vbdiff(ji,jj,jk)    &
313                  &              - e1u(ji  ,jj+1) * ubdiff(ji  ,jj+1,jk) + e1u(ji,jj) * ubdiff(ji,jj,jk)  ) &
314                  &           * fmask(ji,jj,jk) * zbtr
315#endif       
316            END DO
317         END DO
318
319      ENDDO
320
321      !                                                ! ===============
322      DO jk = 1, jpkm1                                 ! Horizontal slab
323         !                                             ! ===============
324         DO jj = 2, jpjm1
325            DO ji = 2, jpim1   ! vector opt.
326#if defined key_zco
327               ! horizontal diffusive trends
328               ze2u = rotdiff (ji,jj,jk)
329               ze1v = hdivdiff(ji,jj,jk)
330               zua = - (                ze2u                  - &
331                  rotdiff (ji,jj-1,jk) ) / e2u(ji,jj)   &
332                  + ( hdivdiff(ji+1,jj,jk) -     &
333                  ze1v                  ) / e1u(ji,jj)
334
335               zva = + (                ze2u                  - &
336                  rotdiff (ji-1,jj,jk) ) / e1v(ji,jj)   &
337                  + ( hdivdiff(ji,jj+1,jk) -       &
338                  ze1v                  ) / e2v(ji,jj)
339#else
340               ze2u = rotdiff (ji,jj,jk)
341               ze1v = hdivdiff(ji,jj,jk)
342               ! horizontal diffusive trends
343               zua = - ( ze2u - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   &
344                  + ( hdivdiff(ji+1,jj,jk) - ze1v      &
345                  ) / e1u(ji,jj)
346
347               zva = + ( ze2u - rotdiff (ji-1,jj,jk)) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   &
348                  + ( hdivdiff(ji,jj+1,jk) - ze1v    &
349                  ) / e2v(ji,jj)
350#endif
351
352               ! add it to the general momentum trends
353               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
354               va(ji,jj,jk) = va(ji,jj,jk) + zva
355            END DO
356         END DO
357         !                                             ! ===============
358      END DO                                           !   End of slab
359      !                                                ! ===============
360
361#endif
362
363   END SUBROUTINE Agrif_Sponge_dyn
364
365   SUBROUTINE interptn(tabres,i1,i2,j1,j2,k1,k2)
366      !!---------------------------------------------
367      !!   *** ROUTINE interptn ***
368      !!---------------------------------------------
369#  include "domzgr_substitute.h90"       
370     
371      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
372      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
373
374      tabres(i1:i2,j1:j2,k1:k2) = tn(i1:i2,j1:j2,k1:k2)
375
376   END SUBROUTINE interptn
377
378   SUBROUTINE interpsn(tabres,i1,i2,j1,j2,k1,k2)
379      !!---------------------------------------------
380      !!   *** ROUTINE interpsn ***
381      !!---------------------------------------------
382#  include "domzgr_substitute.h90"       
383     
384      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
385      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
386
387      tabres(i1:i2,j1:j2,k1:k2) = sn(i1:i2,j1:j2,k1:k2)
388
389   END SUBROUTINE interpsn
390
391   SUBROUTINE interpun(tabres,i1,i2,j1,j2,k1,k2)
392      !!---------------------------------------------
393      !!   *** ROUTINE interpun ***
394      !!---------------------------------------------
395#  include "domzgr_substitute.h90"       
396     
397      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
398      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
399
400      tabres(i1:i2,j1:j2,k1:k2) = un(i1:i2,j1:j2,k1:k2)
401
402   END SUBROUTINE interpun
403
404   SUBROUTINE interpvn(tabres,i1,i2,j1,j2,k1,k2)
405      !!---------------------------------------------
406      !!   *** ROUTINE interpvn ***
407      !!---------------------------------------------
408#  include "domzgr_substitute.h90"       
409     
410      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
411      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
412
413      tabres(i1:i2,j1:j2,k1:k2) = vn(i1:i2,j1:j2,k1:k2)
414
415   END SUBROUTINE interpvn
416
417#else
418CONTAINS
419
420   SUBROUTINE agrif_opa_sponge_empty
421      !!---------------------------------------------
422      !!   *** ROUTINE agrif_OPA_sponge_empty ***
423      !!---------------------------------------------
424      WRITE(*,*)  'agrif_opa_sponge : You should not have seen this print! error?'
425   END SUBROUTINE agrif_opa_sponge_empty
426#endif
427
428END MODULE agrif_opa_sponge
Note: See TracBrowser for help on using the repository browser.