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 NEMO/branches/UKMO/r8395_coupling_sequence/NEMOGCM/NEMO/NST_SRC – NEMO

source: NEMO/branches/UKMO/r8395_coupling_sequence/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 10762

Last change on this file since 10762 was 10762, checked in by jcastill, 5 years ago

Revert previous changes as the removal of keywords was not uncoupled of the actual changes
File size: 16.4 KB
Line 
1#define SPONGE && define SPONGE_TOP
2
3MODULE agrif_opa_sponge
4   !!======================================================================
5   !!                ***  MODULE agrif_opa_update  ***
6   !! AGRIF :   
7   !!======================================================================
8   !! History : 
9   !!----------------------------------------------------------------------
10#if defined key_agrif
11   USE par_oce
12   USE oce
13   USE dom_oce
14   USE in_out_manager
15   USE agrif_oce
16   USE wrk_nemo 
17   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
18
19   IMPLICIT NONE
20   PRIVATE
21
22   PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn
23   PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge
24
25   !!----------------------------------------------------------------------
26   !! NEMO/NST 3.7 , NEMO Consortium (2015)
27   !! $Id: agrif_opa_sponge.F90 7646 2017-02-06 09:25:03Z timgraham $
28   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
29   !!----------------------------------------------------------------------
30CONTAINS
31
32   SUBROUTINE Agrif_Sponge_Tra
33      !!---------------------------------------------
34      !!   *** ROUTINE Agrif_Sponge_Tra ***
35      !!---------------------------------------------
36      REAL(wp) :: timecoeff
37      !!---------------------------------------------
38      !
39#if defined SPONGE
40      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
41
42      CALL Agrif_Sponge
43      Agrif_SpecialValue=0.
44      Agrif_UseSpecialValue = .TRUE.
45      tabspongedone_tsn = .FALSE.
46
47      CALL Agrif_Bc_Variable(tsn_sponge_id,calledweight=timecoeff,procname=interptsn_sponge)
48
49      Agrif_UseSpecialValue = .FALSE.
50#endif
51      !
52   END SUBROUTINE Agrif_Sponge_Tra
53
54
55   SUBROUTINE Agrif_Sponge_dyn
56      !!---------------------------------------------
57      !!   *** ROUTINE Agrif_Sponge_dyn ***
58      !!---------------------------------------------
59      REAL(wp) :: timecoeff
60      !!---------------------------------------------
61
62#if defined SPONGE
63      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
64
65      Agrif_SpecialValue=0.
66      Agrif_UseSpecialValue = ln_spc_dyn
67
68      tabspongedone_u = .FALSE.
69      tabspongedone_v = .FALSE.         
70      CALL Agrif_Bc_Variable(un_sponge_id,calledweight=timecoeff,procname=interpun_sponge)
71
72      tabspongedone_u = .FALSE.
73      tabspongedone_v = .FALSE.
74      CALL Agrif_Bc_Variable(vn_sponge_id,calledweight=timecoeff,procname=interpvn_sponge)
75
76      Agrif_UseSpecialValue = .FALSE.
77#endif
78      !
79   END SUBROUTINE Agrif_Sponge_dyn
80
81
82   SUBROUTINE Agrif_Sponge
83      !!---------------------------------------------
84      !!   *** ROUTINE  Agrif_Sponge ***
85      !!---------------------------------------------
86      INTEGER  :: ji,jj,jk
87      INTEGER  :: ispongearea, ilci, ilcj
88      LOGICAL  :: ll_spdone
89      REAL(wp) :: z1spongearea, zramp
90      REAL(wp), POINTER, DIMENSION(:,:) :: ztabramp
91
92#if defined SPONGE || defined SPONGE_TOP
93      ll_spdone=.TRUE.
94      IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN
95         ! Define ramp from boundaries towards domain interior
96         ! at T-points
97         ! Store it in ztabramp
98         ll_spdone=.FALSE.
99
100         CALL wrk_alloc( jpi, jpj, ztabramp )
101
102         ispongearea  = 2 + nn_sponge_len * Agrif_irhox()
103         ilci = nlci - ispongearea
104         ilcj = nlcj - ispongearea 
105         z1spongearea = 1._wp / REAL( ispongearea - 2 )
106
107         ztabramp(:,:) = 0._wp
108
109         IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN
110            DO jj = 1, jpj
111               IF ( umask(2,jj,1) == 1._wp ) THEN
112                 DO ji = 2, ispongearea                 
113                    ztabramp(ji,jj) = ( ispongearea-ji ) * z1spongearea
114                 END DO
115               ENDIF
116            ENDDO
117         ENDIF
118
119         IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN
120            DO jj = 1, jpj
121               IF ( umask(nlci-2,jj,1) == 1._wp ) THEN
122                  DO ji = ilci+1,nlci-1
123                     zramp = (ji - (ilci+1) ) * z1spongearea
124                     ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
125                  ENDDO
126               ENDIF
127            ENDDO
128         ENDIF
129
130         IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN
131            DO ji = 1, jpi
132               IF ( vmask(ji,2,1) == 1._wp ) THEN
133                  DO jj = 2, ispongearea
134                     zramp = ( ispongearea-jj ) * z1spongearea
135                     ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
136                  END DO
137               ENDIF
138            ENDDO
139         ENDIF
140
141         IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN
142            DO ji = 1, jpi
143               IF ( vmask(ji,nlcj-2,1) == 1._wp ) THEN
144                  DO jj = ilcj+1,nlcj-1
145                     zramp = (jj - (ilcj+1) ) * z1spongearea
146                     ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
147                  END DO
148               ENDIF
149            ENDDO
150         ENDIF
151
152      ENDIF
153
154      ! Tracers
155      IF( .NOT. spongedoneT ) THEN
156         fsaht_spu(:,:) = 0._wp
157         fsaht_spv(:,:) = 0._wp
158         DO jj = 2, jpjm1
159            DO ji = 2, jpim1   ! vector opt.
160               fsaht_spu(ji,jj) = 0.5_wp * visc_tra * (ztabramp(ji,jj) + ztabramp(ji+1,jj  ))
161               fsaht_spv(ji,jj) = 0.5_wp * visc_tra * (ztabramp(ji,jj) + ztabramp(ji  ,jj+1))
162            END DO
163         END DO
164
165         CALL lbc_lnk( fsaht_spu, 'U', 1. )   ! Lateral boundary conditions
166         CALL lbc_lnk( fsaht_spv, 'V', 1. )
167         spongedoneT = .TRUE.
168      ENDIF
169
170      ! Dynamics
171      IF( .NOT. spongedoneU ) THEN
172         fsahm_spt(:,:) = 0._wp
173         fsahm_spf(:,:) = 0._wp
174         DO jj = 2, jpjm1
175            DO ji = 2, jpim1   ! vector opt.
176               fsahm_spt(ji,jj) = visc_dyn * ztabramp(ji,jj)
177               fsahm_spf(ji,jj) = 0.25_wp * visc_dyn * ( ztabramp(ji,jj) + ztabramp(ji  ,jj+1) &
178                                                     &  +ztabramp(ji,jj) + ztabramp(ji+1,jj  ) )
179            END DO
180         END DO
181
182         CALL lbc_lnk( fsahm_spt, 'T', 1. )   ! Lateral boundary conditions
183         CALL lbc_lnk( fsahm_spf, 'F', 1. )
184         spongedoneU = .TRUE.
185      ENDIF
186      !
187      IF (.NOT.ll_spdone) CALL wrk_dealloc( jpi, jpj, ztabramp )
188      !
189#endif
190      !
191   END SUBROUTINE Agrif_Sponge
192
193
194   SUBROUTINE interptsn_sponge(tabres,i1,i2,j1,j2,k1,k2,n1,n2,before)
195      !!---------------------------------------------
196      !!   *** ROUTINE interptsn_sponge ***
197      !!---------------------------------------------
198      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
199      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
200      LOGICAL, INTENT(in) :: before
201      !
202      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
203      INTEGER  ::   iku, ikv
204      REAL(wp) :: ztsa, zabe1, zabe2, zbtr
205      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ztu, ztv
206      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2) ::tsbdiff
207      !
208      IF( before ) THEN
209         tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
210      ELSE   
211         !
212         tsbdiff(:,:,:,:) = tsb(i1:i2,j1:j2,:,:) - tabres(:,:,:,:)   
213         DO jn = 1, jpts           
214            DO jk = 1, jpkm1
215               DO jj = j1,j2-1
216                  DO ji = i1,i2-1
217                     zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)
218                     zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)
219                     ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj  ,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 
220                     ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji  ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
221                  END DO
222               END DO
223               !
224               IF( ln_zps ) THEN      ! set gradient at partial step level
225                  DO jj = j1,j2-1
226                     DO ji = i1,i2-1
227                        ! last level
228                        iku = mbku(ji,jj)
229                        ikv = mbkv(ji,jj)
230                        IF( iku == jk )   ztu(ji,jj,jk) = 0._wp
231                        IF( ikv == jk )   ztv(ji,jj,jk) = 0._wp
232                     END DO
233                  END DO
234               ENDIF
235            END DO
236            !
237            DO jk = 1, jpkm1
238               DO jj = j1+1,j2-1
239                  DO ji = i1+1,i2-1
240                     IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
241                        zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
242                        ! horizontal diffusive trends
243                        ztsa = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji  ,jj-1,jk)  )
244                        ! add it to the general tracer trends
245                        tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa
246                     ENDIF
247                  END DO
248               END DO
249            END DO
250            !
251         END DO
252         !
253         tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
254         !
255      ENDIF
256      !
257   END SUBROUTINE interptsn_sponge
258
259
260   SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2, before)
261      !!---------------------------------------------
262      !!   *** ROUTINE interpun_sponge ***
263      !!---------------------------------------------   
264      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
265      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
266      LOGICAL, INTENT(in) :: before
267
268      INTEGER :: ji,jj,jk
269
270      ! sponge parameters
271      REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
272      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ubdiff
273      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff
274      INTEGER :: jmax
275      !!---------------------------------------------   
276      !
277      IF( before ) THEN
278         tabres = un(i1:i2,j1:j2,:)
279      ELSE
280         ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(:,:,:))*umask(i1:i2,j1:j2,:)
281         !
282         DO jk = 1, jpkm1                                 ! Horizontal slab
283            !                                             ! ===============
284
285            !                                             ! --------
286            ! Horizontal divergence                       !   div
287            !                                             ! --------
288            DO jj = j1,j2
289               DO ji = i1+1,i2   ! vector opt.
290                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
291                  hdivdiff(ji,jj,jk) = (  e2u(ji  ,jj)*e3u_n(ji  ,jj,jk) * ubdiff(ji  ,jj,jk) &
292                                     &   -e2u(ji-1,jj)*e3u_n(ji-1,jj,jk) * ubdiff(ji-1,jj,jk) ) * zbtr
293               END DO
294            END DO
295
296            DO jj = j1,j2-1
297               DO ji = i1,i2   ! vector opt.
298                  zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
299                  rotdiff(ji,jj,jk) = (-e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) &
300                                       +e1u(ji,jj  ) * ubdiff(ji,jj  ,jk) & 
301                                    & ) * fmask(ji,jj,jk) * zbtr 
302               END DO
303            END DO
304         END DO
305         !
306         DO jj = j1+1, j2-1
307            DO ji = i1+1, i2-1   ! vector opt.
308
309               IF (.NOT. tabspongedone_u(ji,jj)) THEN
310                  DO jk = 1, jpkm1                                 ! Horizontal slab
311                     ze2u = rotdiff (ji,jj,jk)
312                     ze1v = hdivdiff(ji,jj,jk)
313                     ! horizontal diffusive trends
314                     zua = - ( ze2u - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) )   &
315                           + ( hdivdiff(ji+1,jj,jk) - ze1v  ) / e1u(ji,jj)
316
317                     ! add it to the general momentum trends
318                     ua(ji,jj,jk) = ua(ji,jj,jk) + zua
319
320                  END DO
321               ENDIF
322
323            END DO
324         END DO
325
326         tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
327
328         jmax = j2-1
329         IF ((nbondj == 1).OR.(nbondj == 2)) jmax = MIN(jmax,nlcj-3)
330
331         DO jj = j1+1, jmax
332            DO ji = i1+1, i2   ! vector opt.
333
334               IF (.NOT. tabspongedone_v(ji,jj)) THEN
335                  DO jk = 1, jpkm1                                 ! Horizontal slab
336                     ze2u = rotdiff (ji,jj,jk)
337                     ze1v = hdivdiff(ji,jj,jk)
338
339                     ! horizontal diffusive trends
340                     zva = + ( ze2u - rotdiff (ji-1,jj,jk)) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) )   &
341                           + ( hdivdiff(ji,jj+1,jk) - ze1v  ) / e2v(ji,jj)
342
343                     ! add it to the general momentum trends
344                     va(ji,jj,jk) = va(ji,jj,jk) + zva
345                  END DO
346               ENDIF
347               !
348            END DO
349         END DO
350         !
351         tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
352         !
353      ENDIF
354      !
355   END SUBROUTINE interpun_sponge
356
357
358   SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2, before,nb,ndir)
359      !!---------------------------------------------
360      !!   *** ROUTINE interpvn_sponge ***
361      !!---------------------------------------------
362      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
363      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
364      LOGICAL, INTENT(in) :: before
365      INTEGER, INTENT(in) :: nb , ndir
366      !
367      INTEGER  ::   ji, jj, jk
368      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr
369      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: vbdiff
370      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff
371      INTEGER :: imax
372      !!---------------------------------------------
373
374      IF( before ) THEN
375         tabres = vn(i1:i2,j1:j2,:)
376      ELSE
377         !
378         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(:,:,:))*vmask(i1:i2,j1:j2,:)
379         !
380         DO jk = 1, jpkm1                                 ! Horizontal slab
381            !                                             ! ===============
382
383            !                                             ! --------
384            ! Horizontal divergence                       !   div
385            !                                             ! --------
386            DO jj = j1+1,j2
387               DO ji = i1,i2   ! vector opt.
388                  zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
389                  hdivdiff(ji,jj,jk) = ( e1v(ji,jj  ) * e3v_n(ji,jj  ,jk) * vbdiff(ji,jj  ,jk)  &
390                                     &  -e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vbdiff(ji,jj-1,jk)  ) * zbtr
391               END DO
392            END DO
393            DO jj = j1,j2
394               DO ji = i1,i2-1   ! vector opt.
395                  zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
396                  rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 
397                                    &  -e2v(ji  ,jj) * vbdiff(ji  ,jj,jk)  ) * fmask(ji,jj,jk) * zbtr
398               END DO
399            END DO
400         END DO
401
402         !                                                ! ===============
403         !                                               
404
405         imax = i2-1
406         IF ((nbondi == 1).OR.(nbondi == 2))   imax = MIN(imax,nlci-3)
407
408         DO jj = j1+1, j2
409            DO ji = i1+1, imax   ! vector opt.
410               IF( .NOT. tabspongedone_u(ji,jj) ) THEN
411                  DO jk = 1, jpkm1
412                     ua(ji,jj,jk) = ua(ji,jj,jk)                                                               &
413                        & - ( rotdiff (ji  ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) )  &
414                        & + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj  ,jk)) * r1_e1u(ji,jj)
415                  END DO
416               ENDIF
417            END DO
418         END DO
419         !
420         tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
421         !
422         DO jj = j1+1, j2-1
423            DO ji = i1+1, i2-1   ! vector opt.
424               IF( .NOT. tabspongedone_v(ji,jj) ) THEN
425                  DO jk = 1, jpkm1
426                     va(ji,jj,jk) = va(ji,jj,jk)                                                                  &
427                        &  + ( rotdiff (ji,jj  ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) )   &
428                        &  + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji  ,jj,jk) ) * r1_e2v(ji,jj)
429                  END DO
430               ENDIF
431            END DO
432         END DO
433         tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
434      ENDIF
435      !
436   END SUBROUTINE interpvn_sponge
437
438#else
439CONTAINS
440   SUBROUTINE agrif_opa_sponge_empty
441      !!---------------------------------------------
442      !!   *** ROUTINE agrif_OPA_sponge_empty ***
443      !!---------------------------------------------
444      WRITE(*,*)  'agrif_opa_sponge : You should not have seen this print! error?'
445   END SUBROUTINE agrif_opa_sponge_empty
446#endif
447
448   !!======================================================================
449END MODULE agrif_opa_sponge
Note: See TracBrowser for help on using the repository browser.