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

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

UKMO/NEMO_4.0.1_biharmonic_GM : New calculation for zmu with new rn_bgm_msc namelist parameter.
File size: 13.1 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
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, pahmt, pahmf, 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   ) ::   pahmt, pahmf ! viscosity coefficients
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 pahmf has already been multiplied by fmask
82               zcur(ji-1,jj-1) = pahmf(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 pahmt has already been multiplied by tmask
87               zdiv(ji,jj)     = pahmt(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   SUBROUTINE dyn_ldf_bgm( kt, pub, pvb, pua, pva )
111      !!----------------------------------------------------------------------
112      !!                 ***  ROUTINE dyn_bgm  ***
113      !!                   
114      !! ** Purpose :   Compute the before lateral momentum trend due to the bi-Laplacian GM parameterisation
115      !!              and add it to the general trend of momentum equation.
116      !!
117      !! ** Method  :   The bi-Laplacian implementation of GM is via a -d/dz(diffusivity x d/dz(Laplacian of velocity))
118      !!      operator applied at the 'now' time level.  The existing code for the Laplacian contains the 'before' time also in zdiv.
119      !!      It is computed by a call to dyn_ldf_lap routine and vertical differentiation applied twice.
120      !!
121      !! ** Action :   pua, pva increased with the before bi-Laplacian GM momentum trend calculated from pub, pvb.
122      !!----------------------------------------------------------------------
123      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
124      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity fields
125      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! momentum trend
126      !   
127      INTEGER  ::   iku, ikv                         ! local integers
128      !CW
129      INTEGER  ::   ji, jj, jk   ! dummy loop indices
130      !
131      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap     ! laplacian at u- and v-point
132      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulapdz, zvlapdz ! -1*bhm * d/dz(del^2 u) at u- and v-point
133      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmu              ! = bhm / avm
134      !!----------------------------------------------------------------------
135      !
136      IF( kt == nit000 )  THEN
137         IF(lwp) WRITE(numout,*)
138         IF(lwp) WRITE(numout,*) 'dyn_bgm : bi-Laplacian GM operator via momentum '
139         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
140      ENDIF
141      !
142      ! Calculate (del2 u) by calling dyn_ldf_lap with "viscosity" set to 1.0.
143      ! Normally we pass ahmt and ahmf to dyn_ldf_lap, which have been multiplied by tmask and fmask respectively.
144      ! So here we just pass tmask and fmask.
145      zulap(:,:,:) = 0.0_wp ; zvlap(:,:,:) = 0.0_wp
146      CALL dyn_ldf_lap( kt, tmask, fmask, pub, pvb, zulap, zvlap, 1 )   
147      !
148      ! Calculate zmu for stabilising correction terms
149      ! and add to avm to be included in the vertical diffusion calculation later.
150      zmu(:,:,:) = 0.0_wp
151      DO jk = 2, jpkm1 
152         DO jj = 2, jpjm1
153            DO ji = fs_2, jpim1   ! vector opt.
154               zmu(ji,jj,jk) = ( rn_bgm_msc * bhm(ji,jj,jk) / e1e2t(ji,jj) ) * wmask(ji,jj,jk)
155               avm(ji,jj,jk) = avm(ji,jj,jk) + zmu(ji,jj,jk)
156            ENDDO
157         ENDDO
158      ENDDO
159      CALL lbc_lnk_multi( 'dyn_ldf_bgm', zmu, 'W', 1. , avm, 'W', 1. ) 
160     
161      ! ===============
162      !CW: calculate -bhm * d/dz(del^2 u) 
163      ! Use of wumask and wvmask to ensure diffusive fluxes at topography are zero.
164      DO jk = 2, jpkm1 
165         DO jj = 2, jpjm1
166            DO ji = fs_2, jpim1   ! vector opt.
167                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)) &
168     &                                -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) 
169                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)) &
170     &                                -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) 
171            ENDDO
172         ENDDO
173      ENDDO
174
175!CW: set boundary conditions: d/dz(del^2 u) = 0 at top and bottom, so that eddy-induced velocity, w*=0
176!     Surface
177      zulapdz(:,:,1)   = 0._wp   ;   zvlapdz(:,:,1)   = 0._wp 
178!     Flat bottom case
179      zulapdz(:,:,jpk) = 0._wp   ;   zvlapdz(:,:,jpk) = 0._wp
180!     Variable bathymetry case, including z-partial steps, as in dynhpg.F90, subroutine hpg_zps
181        DO jj = 2, jpjm1
182            DO ji = 2, jpim1
183                iku = mbku(ji,jj)+1
184                ikv = mbkv(ji,jj)+1
185                zulapdz(:,:,iku) = 0._wp
186                zvlapdz(:,:,ikv) = 0._wp
187            ENDDO
188        ENDDO
189
190!!    calculate d/dz(-bhm * d/dz(del^2 u))
191                                                       ! ===============
192      DO jk = 1, jpkm1                                 ! Horizontal slab
193         !                                             ! ===============
194         DO jj = 2, jpjm1
195            DO ji = fs_2, jpim1   ! vector opt.
196            pua(ji,jj,jk) = pua(ji,jj,jk) + (zulapdz(ji,jj,jk) - zulapdz(ji,jj,jk+1)) / e3u_n(ji,jj,jk)
197            pva(ji,jj,jk) = pva(ji,jj,jk) + (zvlapdz(ji,jj,jk) - zvlapdz(ji,jj,jk+1)) / e3v_n(ji,jj,jk)
198            ENDDO
199         ENDDO
200      ENDDO
201
202      ! -----
203
204   
205   END SUBROUTINE dyn_ldf_bgm
206   
207
208   SUBROUTINE dyn_ldf_blp( kt, pub, pvb, pua, pva )
209      !!----------------------------------------------------------------------
210      !!                 ***  ROUTINE dyn_ldf_blp  ***
211      !!                   
212      !! ** Purpose :   Compute the before lateral momentum viscous trend
213      !!              and add it to the general trend of momentum equation.
214      !!
215      !! ** Method  :   The lateral viscous trends is provided by a bilaplacian
216      !!      operator applied to before field (forward in time).
217      !!      It is computed by two successive calls to dyn_ldf_lap routine
218      !!
219      !! ** Action :   pta   updated with the before rotated bilaplacian diffusion
220      !!----------------------------------------------------------------------
221      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
222      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity fields
223      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! momentum trend
224      !
225      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point
226      !!----------------------------------------------------------------------
227      !
228      IF( kt == nit000 )  THEN
229         IF(lwp) WRITE(numout,*)
230         IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum '
231         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
232      ENDIF
233      !
234      zulap(:,:,:) = 0._wp
235      zvlap(:,:,:) = 0._wp
236      !
237      CALL dyn_ldf_lap( kt, ahmt, ahmf, pub, pvb, zulap, zvlap, 1 )   ! rotated laplacian applied to ptb (output in zlap)
238      !
239      CALL lbc_lnk_multi( 'dynldf_lap_blp', zulap, 'U', -1., zvlap, 'V', -1. )             ! Lateral boundary conditions
240      !
241      CALL dyn_ldf_lap( kt, ahmt, ahmf, zulap, zvlap, pua, pva, 2 )   ! rotated laplacian applied to zlap (output in pta)
242      !
243   END SUBROUTINE dyn_ldf_blp
244
245   !!======================================================================
246END MODULE dynldf_lap_blp
Note: See TracBrowser for help on using the repository browser.