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

source: trunk/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 6140

Last change on this file since 6140 was 6140, checked in by timgraham, 8 years ago

Merge of branches/2015/dev_merge_2015 back into trunk. Merge excludes NEMOGCM/TOOLS/OBSTOOLS/ for now due to issues with the change of file type. Will sort these manually with further commits.

Branch merged as follows:
In the working copy of branch ran:
svn merge svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk@HEAD
Small conflicts due to bug fixes applied to trunk since the dev_merge_2015 was copied. Bug fixes were applied to the branch as well so these were easy to resolve.
Branch committed at this stage

In working copy run:
svn switch svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/trunk
to switch working copy

Run:
svn merge --reintegrate svn+ssh://forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/2015/dev_merge_2015
to merge the branch into the trunk and then commit - no conflicts at this stage.

  • Property svn:keywords set to Id
File size: 16.4 KB
RevLine 
[3680]1#define SPONGE && define SPONGE_TOP
[390]2
[5656]3MODULE agrif_opa_sponge
[6140]4   !!======================================================================
5   !!                ***  MODULE agrif_opa_update  ***
6   !! AGRIF :   
7   !!======================================================================
8   !! History : 
9   !!----------------------------------------------------------------------
[2528]10#if defined key_agrif  && ! defined key_offline
[636]11   USE par_oce
12   USE oce
13   USE dom_oce
14   USE in_out_manager
[782]15   USE agrif_oce
[3294]16   USE wrk_nemo 
[5656]17   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
[390]18
[636]19   IMPLICIT NONE
20   PRIVATE
[390]21
[5656]22   PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn
23   PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge
[390]24
[1156]25   !!----------------------------------------------------------------------
[6140]26   !! NEMO/NST 3.7 , NEMO Consortium (2015)
[1156]27   !! $Id$
[2528]28   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[1156]29   !!----------------------------------------------------------------------
[5656]30CONTAINS
[636]31
[782]32   SUBROUTINE Agrif_Sponge_Tra
[636]33      !!---------------------------------------------
34      !!   *** ROUTINE Agrif_Sponge_Tra ***
35      !!---------------------------------------------
36      REAL(wp) :: timecoeff
[6140]37      !!---------------------------------------------
38      !
[390]39#if defined SPONGE
[636]40      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
41
[5656]42      CALL Agrif_Sponge
[390]43      Agrif_SpecialValue=0.
44      Agrif_UseSpecialValue = .TRUE.
[5656]45      tabspongedone_tsn = .FALSE.
[390]46
[5656]47      CALL Agrif_Bc_Variable(tsn_sponge_id,calledweight=timecoeff,procname=interptsn_sponge)
[390]48
[5656]49      Agrif_UseSpecialValue = .FALSE.
[390]50#endif
[6140]51      !
[636]52   END SUBROUTINE Agrif_Sponge_Tra
[390]53
[6140]54
[782]55   SUBROUTINE Agrif_Sponge_dyn
[636]56      !!---------------------------------------------
57      !!   *** ROUTINE Agrif_Sponge_dyn ***
58      !!---------------------------------------------
59      REAL(wp) :: timecoeff
[6140]60      !!---------------------------------------------
[390]61
62#if defined SPONGE
[636]63      timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
[390]64
65      Agrif_SpecialValue=0.
[782]66      Agrif_UseSpecialValue = ln_spc_dyn
[390]67
[5656]68      tabspongedone_u = .FALSE.
69      tabspongedone_v = .FALSE.         
70      CALL Agrif_Bc_Variable(un_sponge_id,calledweight=timecoeff,procname=interpun_sponge)
[390]71
[5656]72      tabspongedone_u = .FALSE.
73      tabspongedone_v = .FALSE.
74      CALL Agrif_Bc_Variable(vn_sponge_id,calledweight=timecoeff,procname=interpvn_sponge)
75
[390]76      Agrif_UseSpecialValue = .FALSE.
77#endif
[6140]78      !
[636]79   END SUBROUTINE Agrif_Sponge_dyn
[390]80
[6140]81
[3680]82   SUBROUTINE Agrif_Sponge
83      !!---------------------------------------------
84      !!   *** ROUTINE  Agrif_Sponge ***
85      !!---------------------------------------------
86      INTEGER  :: ji,jj,jk
87      INTEGER  :: ispongearea, ilci, ilcj
[4153]88      LOGICAL  :: ll_spdone
89      REAL(wp) :: z1spongearea, zramp
90      REAL(wp), POINTER, DIMENSION(:,:) :: ztabramp
[3680]91
92#if defined SPONGE || defined SPONGE_TOP
[4153]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.
[3680]99
[4153]100         CALL wrk_alloc( jpi, jpj, ztabramp )
[3680]101
[5656]102         ispongearea  = 2 + nn_sponge_len * Agrif_irhox()
[4153]103         ilci = nlci - ispongearea
104         ilcj = nlcj - ispongearea 
105         z1spongearea = 1._wp / REAL( ispongearea - 2 )
[3680]106
[5656]107         ztabramp(:,:) = 0._wp
[4153]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
[3680]154      ! Tracers
155      IF( .NOT. spongedoneT ) THEN
[5656]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
[3680]164
[5656]165         CALL lbc_lnk( fsaht_spu, 'U', 1. )   ! Lateral boundary conditions
166         CALL lbc_lnk( fsaht_spv, 'V', 1. )
[3680]167         spongedoneT = .TRUE.
168      ENDIF
169
170      ! Dynamics
171      IF( .NOT. spongedoneU ) THEN
[5656]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
[3680]181
[5656]182         CALL lbc_lnk( fsahm_spt, 'T', 1. )   ! Lateral boundary conditions
183         CALL lbc_lnk( fsahm_spf, 'F', 1. )
[3680]184         spongedoneU = .TRUE.
185      ENDIF
186      !
[4153]187      IF (.NOT.ll_spdone) CALL wrk_dealloc( jpi, jpj, ztabramp )
[3680]188      !
189#endif
[6140]190      !
[3680]191   END SUBROUTINE Agrif_Sponge
192
[6140]193
[5656]194   SUBROUTINE interptsn_sponge(tabres,i1,i2,j1,j2,k1,k2,n1,n2,before)
[636]195      !!---------------------------------------------
[5656]196      !!   *** ROUTINE interptsn_sponge ***
[636]197      !!---------------------------------------------
[3294]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
[5656]200      LOGICAL, INTENT(in) :: before
[6140]201      !
[5656]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      !
[6140]208      IF( before ) THEN
[5656]209         tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2)
210      ELSE   
[6140]211         !
[5656]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
[6140]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)
[5656]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) )
[6140]221                  END DO
222               END DO
223               !
[5656]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)
[6140]230                        IF( iku == jk )   ztu(ji,jj,jk) = 0._wp
231                        IF( ikv == jk )   ztv(ji,jj,jk) = 0._wp
[5656]232                     END DO
233                  END DO
234               ENDIF
[6140]235            END DO
236            !
[5656]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
[6140]241                        zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
[5656]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
[6140]247                  END DO
248               END DO
249            END DO
250            !
251         END DO
252         !
[5656]253         tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
[6140]254         !
[5656]255      ENDIF
[6140]256      !
[5656]257   END SUBROUTINE interptsn_sponge
258
[6140]259
[5656]260   SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2, before)
[636]261      !!---------------------------------------------
[5656]262      !!   *** ROUTINE interpun_sponge ***
263      !!---------------------------------------------   
[636]264      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
265      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
[5656]266      LOGICAL, INTENT(in) :: before
[390]267
[5656]268      INTEGER :: ji,jj,jk
[390]269
[5656]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
[6140]275      !!---------------------------------------------   
[5656]276      !
[6140]277      IF( before ) THEN
[5656]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,:)
[6140]281         !
[5656]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.
[6140]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
[5656]293               END DO
294            END DO
295
296            DO jj = j1,j2-1
297               DO ji = i1,i2   ! vector opt.
[6140]298                  zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
[5656]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
[6140]304         END DO
[5656]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
[6140]314                     zua = - ( ze2u - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) )   &
[5656]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
[6140]340                     zva = + ( ze2u - rotdiff (ji-1,jj,jk)) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) )   &
[5656]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
[6140]347               !
[5656]348            END DO
349         END DO
[6140]350         !
[5656]351         tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
[6140]352         !
[5656]353      ENDIF
[6140]354      !
[5656]355   END SUBROUTINE interpun_sponge
356
357
358   SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2, before,nb,ndir)
[636]359      !!---------------------------------------------
[5656]360      !!   *** ROUTINE interpvn_sponge ***
361      !!---------------------------------------------
[636]362      INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
363      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
[5656]364      LOGICAL, INTENT(in) :: before
365      INTEGER, INTENT(in) :: nb , ndir
[6140]366      !
367      INTEGER  ::   ji, jj, jk
368      REAL(wp) ::   ze2u, ze1v, zua, zva, zbtr
[5656]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
[6140]372      !!---------------------------------------------
[5656]373
[6140]374      IF( before ) THEN
[5656]375         tabres = vn(i1:i2,j1:j2,:)
376      ELSE
[6140]377         !
[5656]378         vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(:,:,:))*vmask(i1:i2,j1:j2,:)
[6140]379         !
[5656]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.
[6140]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
[5656]391               END DO
392            END DO
393            DO jj = j1,j2
394               DO ji = i1,i2-1   ! vector opt.
[6140]395                  zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
[5656]396                  rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 
[6140]397                                    &  -e2v(ji  ,jj) * vbdiff(ji  ,jj,jk)  ) * fmask(ji,jj,jk) * zbtr
[5656]398               END DO
399            END DO
[6140]400         END DO
[5656]401
402         !                                                ! ===============
403         !                                               
404
405         imax = i2-1
[6140]406         IF ((nbondi == 1).OR.(nbondi == 2))   imax = MIN(imax,nlci-3)
[5656]407
408         DO jj = j1+1, j2
409            DO ji = i1+1, imax   ! vector opt.
[6140]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)
[5656]415                  END DO
416               ENDIF
417            END DO
418         END DO
[6140]419         !
[5656]420         tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
[6140]421         !
[5656]422         DO jj = j1+1, j2-1
423            DO ji = i1+1, i2-1   ! vector opt.
[6140]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)
[5656]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
[6140]435      !
[5656]436   END SUBROUTINE interpvn_sponge
437
[390]438#else
[636]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
[390]446#endif
[636]447
[6140]448   !!======================================================================
[636]449END MODULE agrif_opa_sponge
Note: See TracBrowser for help on using the repository browser.