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

Last change on this file since 12532 was 12532, checked in by davestorkey, 5 years ago

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