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 @ 2715

Last change on this file since 2715 was 2715, checked in by rblod, 13 years ago

First attempt to put dynamic allocation on the trunk

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