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.
dynldf_lap_blp.F90 in NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.1_biharmonic_GM/src/OCE/DYN/dynldf_lap_blp.F90 @ 12775

Last change on this file since 12775 was 12775, checked in by davestorkey, 4 years ago

UKMO/NEMO_4.0.1_biharmonic_GM : bug fixes.
File size: 17.7 KB
Line 
1MODULE dynldf_lap_blp
2   !!======================================================================
3   !!                   ***  MODULE  dynldf_lap_blp  ***
4   !! Ocean dynamics:  lateral viscosity trend (laplacian and bilaplacian)
5   !!======================================================================
6   !! History : 3.7  ! 2014-01  (G. Madec, S. Masson)  Original code, re-entrant laplacian
7   !!           4.0  ! 2020-02  (C. Wilson, ...) add bhm coefficient for bi-Laplacian GM implementation via momentum     
8   !!----------------------------------------------------------------------
9
10   !!----------------------------------------------------------------------
11   !!   dyn_ldf_lap   : update the momentum trend with the lateral viscosity using an iso-level   laplacian operator
12   !!   dyn_ldf_blp   : update the momentum trend with the lateral viscosity using an iso-level bilaplacian operator
13   !!----------------------------------------------------------------------
14   USE oce            ! ocean dynamics and tracers
15   USE dom_oce        ! ocean space and time domain
16   USE ldfdyn         ! lateral diffusion: eddy viscosity coef.
17   USE ldfslp         ! iso-neutral slopes
18   USE zdf_oce        ! ocean vertical physics
19   USE phycst
20   !
21   USE in_out_manager ! I/O manager
22   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC dyn_ldf_lap  ! called by dynldf.F90
28   PUBLIC dyn_ldf_blp  ! called by dynldf.F90
29   PUBLIC dyn_ldf_bgm  ! called by dynldf.F90
30   PRIVATE dyn_ldf_lap_no_ahm !needed by dyn_ldf_bgm
31
32   !! * Substitutions
33#  include "vectopt_loop_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
36   !! $Id: dynldf_lap_blp.F90 10425 2018-12-19 21:54:16Z smasson $
37   !! Software governed by the CeCILL license (see ./LICENSE)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE dyn_ldf_lap( kt, pub, pvb, pua, pva, kpass )
42      !!----------------------------------------------------------------------
43      !!                     ***  ROUTINE dyn_ldf_lap  ***
44      !!                       
45      !! ** Purpose :   Compute the before horizontal momentum diffusive
46      !!      trend and add it to the general trend of momentum equation.
47      !!
48      !! ** Method  :   The Laplacian operator apply on horizontal velocity is
49      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )
50      !!
51      !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb.
52      !!----------------------------------------------------------------------
53      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
54      INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
55      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity  [m/s]
56      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! velocity trend   [m/s2]
57      !
58      INTEGER  ::   ji, jj, jk   ! dummy loop indices
59      REAL(wp) ::   zsign        ! local scalars
60      REAL(wp) ::   zua, zva     ! local scalars
61      REAL(wp), DIMENSION(jpi,jpj) ::   zcur, zdiv
62      !!----------------------------------------------------------------------
63      !
64      IF( kt == nit000 .AND. lwp ) THEN
65         WRITE(numout,*)
66         WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass
67         WRITE(numout,*) '~~~~~~~ '
68      ENDIF
69      !
70      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign
71      ELSE                    ;   zsign = -1._wp      !  (eddy viscosity coef. >0)
72      ENDIF
73      !
74      !                                                ! ===============
75      DO jk = 1, jpkm1                                 ! Horizontal slab
76         !                                             ! ===============
77         DO jj = 2, jpj
78            DO ji = fs_2, jpi   ! vector opt.
79               !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1)
80!!gm open question here : e3f  at before or now ?    probably now...
81!!gm note that ahmf has already been multiplied by fmask
82               zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &
83                  &     * (  e2v(ji  ,jj-1) * pvb(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk)  &
84                  &        - e1u(ji-1,jj  ) * pub(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk)  )
85               !                                      ! ahm * div        (computed from 2 to jpi/jpj)
86!!gm note that ahmt has already been multiplied by tmask
87               zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk)                                         &
88                  &     * (  e2u(ji,jj)*e3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk)  &
89                  &        + e1v(ji,jj)*e3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk)  )
90            END DO 
91         END DO 
92         !
93         DO jj = 2, jpjm1                             ! - curl( curl) + grad( div )
94            DO ji = fs_2, fs_jpim1   ! vector opt.
95               pua(ji,jj,jk) = pua(ji,jj,jk) + zsign * (                                                 &
96                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk)   &
97                  &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     )
98                  !
99               pva(ji,jj,jk) = pva(ji,jj,jk) + zsign * (                                                 &
100                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk)   &
101                  &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     )
102            END DO
103         END DO
104         !                                             ! ===============
105      END DO                                           !   End of slab
106      !                                                ! ===============
107      !
108   END SUBROUTINE dyn_ldf_lap
109
110!CW: new subroutine for the Laplacian of velocity, copied from dyn_ldf_lap above, but without eddy viscosity ahmf, ahmt
111
112   SUBROUTINE dyn_ldf_lap_no_ahm( kt, pub, pvb, pulap, pvlap, kpass )
113      !!----------------------------------------------------------------------
114      !!                     ***  ROUTINE dyn_ldf_lap_no_ahm  ***
115      !!                       
116      !! ** Purpose :   Companion subroutine to dyn_ldf_bgm to calculate
117      !!      the before horizontal momentum Laplacian
118      !!      trend and add it to the general trend of momentum equation. 
119      !!      Note: mimics dyn_ldf_lap but without eddy viscosity ahmf, ahmt,
120      !!      because biharmonic GM eddy diffusivity is applied in dyn_bgm.
121      !!
122      !! ** Method  :   The Laplacian operator apply on horizontal velocity is
123      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )
124      !!
125      !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb.
126      !!----------------------------------------------------------------------
127      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
128      INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
129      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity  [m/s]
130      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(out) ::   pulap, pvlap ! velocity trend   [m/s2]
131      !
132      INTEGER  ::   ji, jj, jk   ! dummy loop indices
133      REAL(wp) ::   zsign        ! local scalars
134      REAL(wp) ::   zua, zva     ! local scalars
135      REAL(wp), DIMENSION(jpi,jpj) ::   zcur, zdiv
136      !!----------------------------------------------------------------------
137      !
138      IF( kt == nit000 .AND. lwp ) THEN
139         WRITE(numout,*)
140         WRITE(numout,*) 'dyn_ldf_lap_no_ahm : companion iso-level harmonic (laplacian) operator to dyn_bgm, pass=', kpass
141         WRITE(numout,*) '~~~~~~~ '
142      ENDIF
143      !
144      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign
145      ELSE                    ;   zsign = -1._wp      !  (eddy viscosity coef. >0)
146      ENDIF
147      !
148      !                                                ! ===============
149      DO jk = 1, jpkm1                                 ! Horizontal slab
150         !                                             ! ===============
151         DO jj = 2, jpj
152            DO ji = fs_2, jpi   ! vector opt.
153               !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1)
154!!gm open question here : e3f  at before or now ?    probably now...
155
156!!CW: note that the removed ahmf had already been multiplied by fmask, so we multiply by fmask explicitly here,
157!!    with the i,j,k of fmask aligning with those of ahmf, as defined in ldfdyn.F90.
158               zcur(ji-1,jj-1) = fmask(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)      &
159                  &     * (  e2v(ji  ,jj-1) * pvb(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk)  &
160                  &        - e1u(ji-1,jj  ) * pub(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk)  )
161               !                                      ! ahm * div        (computed from 2 to jpi/jpj)
162!!CW: note that the removed ahmt had already been multiplied by tmask, so we multiply by tmask explicitly here,
163!!    with the i,j,k of tmask aligning with those of ahmt, as defined in ldfdyn.F90
164               zdiv(ji,jj)     = tmask(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk)                                         &
165                  &     * (  e2u(ji,jj)*e3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk)  &
166                  &        + e1v(ji,jj)*e3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk)  )
167            END DO 
168         END DO 
169         !
170!CW: multiply pulap, pvlap by umask, vmask
171         DO jj = 2, jpjm1                             ! - curl( curl) + grad( div )
172            DO ji = fs_2, fs_jpim1   ! vector opt.
173               pulap(ji,jj,jk) = umask(ji,jj,jk) * zsign * (                                            &
174                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk)   &
175                  &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     )
176                  !
177               pvlap(ji,jj,jk) = vmask(ji,jj,jk) * zsign * (                                             &
178                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk)   &
179                  &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     )
180            END DO
181         END DO
182         !                                             ! ===============
183      END DO                                           !   End of slab
184      !                                                ! ===============
185      !
186   END SUBROUTINE dyn_ldf_lap_no_ahm
187
188
189
190
191!CW: new subroutine for biharmonic GM, copied from SUBROUTINE dyn_ldf_blp below
192
193   SUBROUTINE dyn_ldf_bgm( kt, pub, pvb, pua, pva )
194      !!----------------------------------------------------------------------
195      !!                 ***  ROUTINE dyn_bgm  ***
196      !!                   
197      !! ** Purpose :   Compute the before lateral momentum trend due to the bi-Laplacian GM parameterisation
198      !!              and add it to the general trend of momentum equation.
199      !!
200      !! ** Method  :   The bi-Laplacian implementation of GM is via a -d/dz(diffusivity x d/dz(Laplacian of velocity))
201      !!      operator applied at the 'now' time level.  The existing code for the Laplacian contains the 'before' time also in zdiv.
202      !!      It is computed by a call to dyn_ldf_lap routine and vertical differentiation applied twice.
203      !!
204      !! ** Action :   pua, pva increased with the before bi-Laplacian GM momentum trend calculated from pub, pvb.
205      !!----------------------------------------------------------------------
206      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
207      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity fields
208      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! momentum trend
209      !   
210      INTEGER  ::   iku, ikv                         ! local integers
211      !CW
212      INTEGER  ::   ji, jj, jk   ! dummy loop indices
213      !
214      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap     ! laplacian at u- and v-point
215      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulapdz, zvlapdz ! -1*bhm * d/dz(del^2 u) at u- and v-point
216      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmu              ! = bhm / avm
217      !!----------------------------------------------------------------------
218      !
219      IF( kt == nit000 )  THEN
220         IF(lwp) WRITE(numout,*)
221         IF(lwp) WRITE(numout,*) 'dyn_bgm : bi-Laplacian GM operator via momentum '
222         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
223      ENDIF
224      !
225 
226      ! Calculate ratio of bhm and avm for stabilising correction terms
227      ! and add to avm to be included in the vertical diffusion calculation later.
228      DO jk = 2, jpkm1 
229         DO jj = 2, jpjm1
230            DO ji = fs_2, jpim1   ! vector opt.
231               zmu(ji,jj,jk) = ( bhm(ji,jj,jk) / (avm(ji,jj,jk) + rsmall) ) * wmask(ji,jj,jk)
232               avm(ji,jj,jk) = avm(ji,jj,jk) + zmu(ji,jj,jk)
233            ENDDO
234         ENDDO
235      ENDDO
236      CALL lbc_lnk_multi( 'dyn_ldf_bgm', zmu, 'W', 1. , avm, 'W', 1. ) 
237     
238      ! Calculate (del2 u)
239      CALL dyn_ldf_lap_no_ahm( kt, pub, pvb, zulap, zvlap, 1 )   
240
241      ! ===============
242      !CW: calculate -bhm * d/dz(del^2 u) 
243      ! Use of wumask and wvmask to ensure diffusive fluxes at topography are zero.
244      DO jk = 2, jpkm1 
245         DO jj = 2, jpjm1
246            DO ji = fs_2, jpim1   ! vector opt.
247                zulapdz(ji,jj,jk) = ( -0.5_wp*(bhm(ji+1,jj  ,jk)+bhm(ji,jj,jk))*(zulap(ji,jj,jk-1) - zulap(ji,jj,jk)) &
248     &                                -0.5_wp*(zmu(ji+1,jj  ,jk)+zmu(ji,jj,jk))*(pub  (ji,jj,jk-1) - pub  (ji,jj,jk)) ) * wumask(ji,jj,jk) / e3uw_n(ji,jj,jk) 
249                zvlapdz(ji,jj,jk) = ( -0.5_wp*(bhm(ji  ,jj+1,jk)+bhm(ji,jj,jk))*(zvlap(ji,jj,jk-1) - zvlap(ji,jj,jk)) &
250     &                                -0.5_wp*(zmu(ji  ,jj+1,jk)+zmu(ji,jj,jk))*(pvb  (ji,jj,jk-1) - pvb  (ji,jj,jk)) ) * wvmask(ji,jj,jk) / e3vw_n(ji,jj,jk) 
251            ENDDO
252         ENDDO
253      ENDDO
254
255!CW: set boundary conditions: d/dz(del^2 u) = 0 at top and bottom, so that eddy-induced velocity, w*=0
256!     Surface
257      zulapdz(:,:,1)   = 0._wp   ;   zvlapdz(:,:,1)   = 0._wp 
258!     Flat bottom case
259      zulapdz(:,:,jpk) = 0._wp   ;   zvlapdz(:,:,jpk) = 0._wp
260!     Variable bathymetry case, including z-partial steps, as in dynhpg.F90, subroutine hpg_zps
261        DO jj = 2, jpjm1
262            DO ji = 2, jpim1
263                iku = mbku(ji,jj)+1
264                ikv = mbkv(ji,jj)+1
265                zulapdz(:,:,iku) = 0._wp
266                zvlapdz(:,:,ikv) = 0._wp
267            ENDDO
268        ENDDO
269
270!!    calculate d/dz(-bhm * d/dz(del^2 u))
271                                                       ! ===============
272      DO jk = 1, jpkm1                                 ! Horizontal slab
273         !                                             ! ===============
274         DO jj = 2, jpjm1
275            DO ji = fs_2, jpim1   ! vector opt.
276            pua(ji,jj,jk) = pua(ji,jj,jk) + (zulapdz(ji,jj,jk) - zulapdz(ji,jj,jk+1)) / e3u_n(ji,jj,jk)
277            pva(ji,jj,jk) = pva(ji,jj,jk) + (zvlapdz(ji,jj,jk) - zvlapdz(ji,jj,jk+1)) / e3v_n(ji,jj,jk)
278            ENDDO
279         ENDDO
280      ENDDO
281
282      ! -----
283
284   
285   END SUBROUTINE dyn_ldf_bgm
286   
287
288   SUBROUTINE dyn_ldf_blp( kt, pub, pvb, pua, pva )
289      !!----------------------------------------------------------------------
290      !!                 ***  ROUTINE dyn_ldf_blp  ***
291      !!                   
292      !! ** Purpose :   Compute the before lateral momentum viscous trend
293      !!              and add it to the general trend of momentum equation.
294      !!
295      !! ** Method  :   The lateral viscous trends is provided by a bilaplacian
296      !!      operator applied to before field (forward in time).
297      !!      It is computed by two successive calls to dyn_ldf_lap routine
298      !!
299      !! ** Action :   pta   updated with the before rotated bilaplacian diffusion
300      !!----------------------------------------------------------------------
301      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
302      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity fields
303      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! momentum trend
304      !
305      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point
306      !!----------------------------------------------------------------------
307      !
308      IF( kt == nit000 )  THEN
309         IF(lwp) WRITE(numout,*)
310         IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum '
311         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
312      ENDIF
313      !
314      zulap(:,:,:) = 0._wp
315      zvlap(:,:,:) = 0._wp
316      !
317      CALL dyn_ldf_lap( kt, pub, pvb, zulap, zvlap, 1 )   ! rotated laplacian applied to ptb (output in zlap)
318      !
319      CALL lbc_lnk_multi( 'dynldf_lap_blp', zulap, 'U', -1., zvlap, 'V', -1. )             ! Lateral boundary conditions
320      !
321      CALL dyn_ldf_lap( kt, zulap, zvlap, pua, pva, 2 )   ! rotated laplacian applied to zlap (output in pta)
322      !
323   END SUBROUTINE dyn_ldf_blp
324
325   !!======================================================================
326END MODULE dynldf_lap_blp
Note: See TracBrowser for help on using the repository browser.