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_top_sponge.F90 in NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST – NEMO

source: NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_top_sponge.F90 @ 14013

Last change on this file since 14013 was 14013, checked in by jchanut, 3 years ago

Manually merge with Christian's branch agrif top sponge/interp #2129

  • Property svn:keywords set to Id
File size: 11.3 KB
Line 
1#define SPONGE_TOP
2
3MODULE agrif_top_sponge
4   !!======================================================================
5   !!                ***  MODULE agrif_top_sponge  ***
6   !! AGRIF :   sponge layer pakage for passive tracers (TOP)
7   !!======================================================================
8   !! History :  2.0  ! 2006-08  (R. Benshila, L. Debreu)  Original code
9   !!----------------------------------------------------------------------
10#if defined key_agrif && defined key_top
11   !!----------------------------------------------------------------------
12   !!   Agrif_Sponge_trc :
13   !!   interptrn_sponge : 
14   !!----------------------------------------------------------------------
15   USE par_oce
16   USE par_trc
17   USE oce
18   USE trc
19   USE dom_oce
20   USE agrif_oce
21   USE agrif_oce_sponge
22   USE vremap
23   !
24   USE in_out_manager
25   USE lib_mpp
26
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC Agrif_Sponge_trc, interptrn_sponge
31
32   !!----------------------------------------------------------------------
33   !! NEMO/NST 4.0 , NEMO Consortium (2018)
34   !! $Id$
35   !! Software governed by the CeCILL license (see ./LICENSE)
36   !!----------------------------------------------------------------------
37CONTAINS
38
39   SUBROUTINE Agrif_Sponge_trc
40      !!----------------------------------------------------------------------
41      !!                   *** ROUTINE Agrif_Sponge_Trc ***
42      !!----------------------------------------------------------------------
43      REAL(wp) ::   zcoef   ! local scalar
44      !!----------------------------------------------------------------------
45      !
46#if defined SPONGE_TOP
47      !! Assume persistence:
48      zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
49
50      Agrif_SpecialValue    = 0._wp
51      Agrif_UseSpecialValue = .TRUE.
52      l_vremap              = ln_vert_remap
53      tabspongedone_trn     = .FALSE.
54      !
55      CALL Agrif_Bc_Variable( trn_sponge_id, calledweight=zcoef, procname=interptrn_sponge )
56      !
57      Agrif_UseSpecialValue = .FALSE.
58      l_vremap              = .FALSE.
59#endif
60      !
61   END SUBROUTINE Agrif_Sponge_Trc
62
63
64   SUBROUTINE interptrn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before) 
65      !!----------------------------------------------------------------------
66      !!                 *** ROUTINE interptrn_sponge ***
67      !!----------------------------------------------------------------------
68      INTEGER                                     , INTENT(in   ) ::   i1, i2, j1, j2, k1, k2, n1, n2
69      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) ::   tabres
70      LOGICAL                                     , INTENT(in   ) ::   before
71      !
72      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
73      INTEGER  ::   iku, ikv
74      REAL(wp) :: ztra, zabe1, zabe2, zbtr, zhtot
75      REAL(wp), DIMENSION(i1-1:i2,j1-1:j2,jpk) :: ztu, ztv
76      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::trbdiff
77      ! vertical interpolation:
78      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
79      REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin, tabin_i
80      REAL(wp), DIMENSION(k1:k2) :: z_in, z_in_i, h_in_i
81      REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
82      INTEGER :: N_in, N_out
83      !!----------------------------------------------------------------------
84      !
85      IF( before ) THEN
86         DO jn = 1, jptra
87            DO jk=k1,k2
88               DO jj=j1,j2
89                  DO ji=i1,i2
90                     tabres(ji,jj,jk,jn) = tr(ji,jj,jk,jn,Kbb_a)
91                  END DO
92               END DO
93            END DO
94         END DO
95
96         IF ( l_vremap.OR.ln_zps ) THEN
97
98            ! Fill cell depths (i.e. gdept) to be interpolated
99            ! Warning: these are masked, hence extrapolated prior interpolation.
100            DO jj=j1,j2
101               DO ji=i1,i2
102                  tabres(ji,jj,k1,jptra+1) = 0.5_wp * tmask(ji,jj,k1) * e3t(ji,jj,k1,Kbb_a)
103                  DO jk=k1+1,k2
104                     tabres(ji,jj,jk,jptra+1) = tmask(ji,jj,jk) * &
105                        & ( tabres(ji,jj,jk-1,jptra+1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a)+e3t(ji,jj,jk,Kbb_a)) )
106                  END DO
107               END DO
108            END DO
109
110            ! Save ssh at last level:
111            IF ( .NOT.ln_linssh ) THEN
112               tabres(i1:i2,j1:j2,k2,jptra+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1) 
113            END IF 
114   
115         END IF
116
117      ELSE 
118         !
119         IF ( l_vremap ) THEN
120
121            IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp
122
123            DO jj=j1,j2
124               DO ji=i1,i2
125
126                  tabres_child(ji,jj,:,:) = 0._wp 
127                  ! Build vertical grids:
128                  N_in = mbkt_parent(ji,jj)
129                  ! Input grid (account for partial cells if any):
130                  DO jk=1,N_in
131                     z_in(jk) = tabres(ji,jj,jk,n2) - tabres(ji,jj,k2,n2)
132                     tabin(jk,1:jptra) = tabres(ji,jj,jk,1:jptra)
133                  END DO
134                 
135                  ! Intermediate grid:
136                  DO jk = 1, N_in
137                     h_in_i(jk) = e3t0_parent(ji,jj,jk) * & 
138                       &       (1._wp + tabres(ji,jj,k2,n2)/(ht0_parent(ji,jj)*ssmask(ji,jj) + 1._wp - ssmask(ji,jj)))
139                  END DO
140                  z_in_i(1) = 0.5_wp * h_in_i(1)
141                  DO jk=2,N_in
142                     z_in_i(jk) = z_in_i(jk-1) + 0.5_wp * ( h_in_i(jk) + h_in_i(jk-1) )
143                  END DO
144                  z_in_i(1:N_in) = z_in_i(1:N_in)  - tabres(ji,jj,k2,n2)
145
146                  ! Output (Child) grid:
147                  N_out = mbkt(ji,jj)
148                  DO jk=1,N_out
149                     h_out(jk) = e3t(ji,jj,jk,Kbb_a)
150                  END DO
151                  z_out(1) = 0.5_wp * h_out(1)
152                  DO jk=2,N_out
153                     z_out(jk) = z_out(jk-1) + 0.5_wp * ( h_out(jk)+h_out(jk-1) )
154                  END DO
155                  IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out)  - ssh(ji,jj,Kbb_a)
156
157                  ! Account for small differences in the free-surface
158                  IF ( sum(h_out(1:N_out)) > sum(h_in_i(1:N_in) )) THEN
159                     h_out(1) = h_out(1)  - ( sum(h_out(1:N_out))-sum(h_in_i(1:N_in)) )
160                  ELSE
161                     h_in_i(1)= h_in_i(1) - ( sum(h_in_i(1:N_in))-sum(h_out(1:N_out)) )
162                  END IF
163                  IF (N_in*N_out > 0) THEN
164                     CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),tabin_i(1:N_in,1:jptra),z_in_i(1:N_in),N_in,N_in,jptra)
165                     CALL reconstructandremap(tabin_i(1:N_in,1:jptra),h_in_i(1:N_in),tabres_child(ji,jj,1:N_out,1:jptra),h_out(1:N_out),N_in,N_out,jptra)
166!                     CALL remap_linear(tabin(1:N_in,1:jptra),z_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jptra),z_out(1:N_in),N_in,N_out,jptra) 
167                  ENDIF
168               END DO
169            END DO
170
171            DO jj=j1,j2
172               DO ji=i1,i2
173                  DO jk=1,jpkm1
174                     trbdiff(ji,jj,jk,1:jptra) = (tr(ji,jj,jk,1:jptra,Kbb_a) - tabres_child(ji,jj,jk,1:jptra)) * tmask(ji,jj,jk)
175                  END DO
176               END DO
177            END DO
178
179         ELSE
180
181            IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells
182
183               DO jj=j1,j2
184                  DO ji=i1,i2
185                     !
186                     N_in  = mbkt(ji,jj) 
187                     N_out = mbkt(ji,jj) 
188                     z_in(1) = tabres(ji,jj,1,n2)
189                     tabin(1,1:jptra) = tabres(ji,jj,1,1:jptra)
190                     DO jk=2, N_in
191                        z_in(jk) = tabres(ji,jj,jk,n2)
192                        tabin(jk,1:jptra) = tabres(ji,jj,jk,1:jptra)
193                     END DO
194                     IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - tabres(ji,jj,k2,n2)
195
196                     z_out(1) = 0.5_wp * e3t(ji,jj,1,Kbb_a)
197                     DO jk=2, N_out
198                        z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk,Kbb_a)) 
199                     END DO
200                     IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Kbb_a)
201
202                     CALL remap_linear(tabin(1:N_in,1:jptra), z_in(1:N_in), tabres(ji,jj,1:N_out,1:jptra), &
203                                         &   z_out(1:N_out), N_in, N_out, jptra)
204                  END DO
205               END DO
206            ENDIF
207
208            DO jj=j1,j2
209               DO ji=i1,i2
210                  DO jk=1,jpkm1
211                     trbdiff(ji,jj,jk,1:jptra) = (tr(ji,jj,jk,1:jptra,Kbb_a) - tabres(ji,jj,jk,1:jptra))*tmask(ji,jj,jk)
212                  END DO
213               END DO
214            END DO
215
216         END IF
217
218         DO jn = 1, jptra           
219            DO jk = 1, jpkm1
220               ztu(i1-1:i2,j1-1:j2,jk) = 0._wp
221               DO jj = j1,j2
222                  DO ji = i1,i2-1
223                     zabe1 = rn_sponge_tra * r1_Dt * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
224                     ztu(ji,jj,jk) = zabe1 * fspu(ji,jj) * ( trbdiff(ji+1,jj  ,jk,jn) - trbdiff(ji,jj,jk,jn) ) 
225                  END DO
226               END DO
227               ztv(i1-1:i2,j1-1:j2,jk) = 0._wp
228               DO ji = i1,i2
229                  DO jj = j1,j2-1
230                     zabe2 = rn_sponge_tra * r1_Dt * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
231                     ztv(ji,jj,jk) = zabe2 * fspv(ji,jj) * ( trbdiff(ji  ,jj+1,jk,jn) - trbdiff(ji,jj,jk,jn) )
232                  END DO
233               END DO
234               !
235               IF( ln_zps ) THEN      ! set gradient at partial step level
236                  DO jj = j1,j2
237                     DO ji = i1,i2
238                        ! last level
239                        iku = mbku(ji,jj)
240                        ikv = mbkv(ji,jj)
241                        IF( iku == jk )   ztu(ji,jj,jk) = 0._wp
242                        IF( ikv == jk )   ztv(ji,jj,jk) = 0._wp
243                     END DO
244                  END DO
245               ENDIF
246            END DO
247            !
248! JC: there is something wrong with the Laplacian in corners
249            DO jk = 1, jpkm1
250               DO jj = j1,j2
251                  DO ji = i1,i2
252                     IF (.NOT. tabspongedone_trn(ji,jj)) THEN
253                        zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)
254                        ! horizontal diffusive trends
255                        ztra = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk)   & 
256                          &           + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) &
257                          &   - rn_trelax_tra * r1_Dt * fspt(ji,jj) * trbdiff(ji,jj,jk,jn)
258                        ! add it to the general tracer trends
259                        tr(ji,jj,jk,jn,Krhs_a) = tr(ji,jj,jk,jn,Krhs_a) + ztra
260                     ENDIF
261                  END DO
262               END DO
263
264            END DO
265            !
266         END DO
267         !
268         tabspongedone_trn(i1:i2,j1:j2) = .TRUE.
269         !
270      ENDIF
271      !
272   END SUBROUTINE interptrn_sponge
273
274#else
275   !!----------------------------------------------------------------------
276   !!   Empty module                                           no TOP AGRIF
277   !!----------------------------------------------------------------------
278CONTAINS
279   SUBROUTINE agrif_top_sponge_empty
280      WRITE(*,*)  'agrif_top_sponge : You should not have seen this print! error?'
281   END SUBROUTINE agrif_top_sponge_empty
282#endif
283
284   !!======================================================================
285END MODULE agrif_top_sponge
Note: See TracBrowser for help on using the repository browser.