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

source: trunk/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 390

Last change on this file since 390 was 390, checked in by opalod, 18 years ago

RB:nemo_v1_update_038: first integration of Agrif :

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