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/2020/dev_r12527_Gurvan_ShallowWater/cfgs/penzAM98/MY_SRC – NEMO

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/cfgs/penzAM98/MY_SRC/dynldf_lap_blp.F90 @ 13562

Last change on this file since 13562 was 13562, checked in by gm, 4 years ago

zgr_pse created only for NS coast
File size: 15.3 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-04  (A. Nasser, G. Madec)  Add symmetric mixing tensor
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   USE lib_mpp        ! MPP library
23   !
24   USE usrdef_nam , ONLY : nn_dynldf_lap_typ    ! use laplacian parameter
25   !
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC dyn_ldf_lap  ! called by dynldf.F90
30   PUBLIC dyn_ldf_blp  ! called by dynldf.F90
31!!anSYM
32   INTEGER, PUBLIC, PARAMETER ::   np_dynldf_lap_rot  = 1         ! div-rot   laplacian
33   INTEGER, PUBLIC, PARAMETER ::   np_dynldf_lap_sym  = 2         ! symmetric laplacian (Griffies&Hallberg 2000)
34   INTEGER, PUBLIC, PARAMETER ::   np_dynldf_lap_symN = 3         ! symmetric laplacian (cartesian)
35
36   !INTEGER, PUBLIC, PARAMETER ::   nn_dynldf_lap_typ = 1         ! choose type of laplacian (ideally from namelist)
37!!anSYM
38   !! * Substitutions
39#  include "do_loop_substitute.h90"
40!!st21
41#  include "domzgr_substitute.h90"
42   !!----------------------------------------------------------------------
43   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
44   !! $Id: dynldf_lap_blp.F90 13416 2020-08-20 10:10:55Z gm $
45   !! Software governed by the CeCILL license (see ./LICENSE)
46   !!----------------------------------------------------------------------
47CONTAINS
48
49   SUBROUTINE dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs, kpass )
50      !!----------------------------------------------------------------------
51      !!                     ***  ROUTINE dyn_ldf_lap  ***
52      !!
53      !! ** Purpose :   Compute the before horizontal momentum diffusive
54      !!      trend and add it to the general trend of momentum equation.
55      !!
56      !! ** Method  :   The Laplacian operator apply on horizontal velocity is
57      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )
58      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )
59      !!
60      !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv.
61      !!
62      !! Reference : S.Griffies, R.Hallberg 2000 Mon.Wea.Rev., DOI:/
63      !!----------------------------------------------------------------------
64      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
65      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm   ! ocean time level indices
66      INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
67      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv     ! before velocity  [m/s]
68      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! velocity trend   [m/s2]
69      !
70      INTEGER  ::   ji, jj, jk   ! dummy loop indices
71      REAL(wp) ::   zsign        ! local scalars
72      REAL(wp) ::   zua, zva     ! local scalars
73      REAL(wp), DIMENSION(jpi,jpj) ::   zcur, zdiv
74      REAL(wp), DIMENSION(jpi,jpj) ::   zten, zshe   ! tension (diagonal) and shearing (anti-diagonal) terms
75      !
76      REAL(wp), DIMENSION(jpi,jpj) ::   ze1e2f, zr1_e1e2f     ! not penalised scale factors
77      !!----------------------------------------------------------------------
78      !
79!!anSYM TO BE ADDED : reading of laplacian operator (ln_dynldf_lap_typ -> to be written nn_) shall be added in dyn_ldf_init
80!!                 as the writing
81!!                 and an integer as np_dynldf_lap for instance taken as argument by dyn_ldf_lap call in dyn_ldf
82      IF( kt == nit000 .AND. lwp ) THEN
83         WRITE(numout,*)
84         WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass
85         WRITE(numout,*) '~~~~~~~ '
86         WRITE(numout,*) '                  nn_dynldf_lap_typ = ', nn_dynldf_lap_typ
87         SELECT CASE( nn_dynldf_lap_typ )             ! print the choice of operator
88         CASE( np_dynldf_lap_rot )   ;   WRITE(numout,*) '   ==>>>   div-rot   laplacian'
89         CASE( np_dynldf_lap_sym )   ;   WRITE(numout,*) '   ==>>>   symmetric laplacian (covariant form)'
90         CASE( np_dynldf_lap_symN)   ;   WRITE(numout,*) '   ==>>>   symmetric laplacian (cartesian form)'
91         END SELECT
92      ENDIF
93      !
94      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign
95      ELSE                    ;   zsign = -1._wp      !  (eddy viscosity coef. >0)
96      ENDIF
97      !
98      SELECT CASE( nn_dynldf_lap_typ )
99         !
100         CASE ( np_dynldf_lap_rot )       !==  Vorticity-Divergence form  ==!
101           !
102           zr1_e1e2f(:,:) = r1_e1e2f(:,:)
103# if defined key_bvp
104           zr1_e1e2f(:,:) = r1_rpof(:,:,1) * zr1_e1e2f(:,:)
105# endif
106
107            DO jk = 1, jpkm1                                 ! Horizontal slab
108               !
109               DO_2D_01_01
110               !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1)
111!!gm note that ahmf has already been multiplied by fmask
112            zcur(ji-1,jj-1) =  &
113               &      ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * zr1_e1e2f(ji-1,jj-1)      &
114               &  * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  &
115               &     - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  )
116            !                                      ! ahm * div        (computed from 2 to jpi/jpj)
117!!gm note that ahmt has already been multiplied by tmask
118                  zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)                                         &
119                     &     * (  e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk)  &
120                     &        + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk)  )
121               END_2D
122               !
123               DO_2D_00_00
124                  pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * (                                             &
125                     &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   &
126                     &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                       )
127                     !
128                  pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * (                                             &
129                     &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   &
130                     &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                       )
131               END_2D
132               !
133            END DO                                           !   End of slab
134            !
135         CASE ( np_dynldf_lap_sym )       !==  Symmetric form  ==!   (Griffies&Hallberg 2000)
136            !
137            DO jk = 1, jpkm1                                 ! Horizontal slab
138               !
139               DO_2D_01_01
140                  !                                      ! shearing stress component (F-point)   NB : ahmf has already been multiplied by fmask
141                  zshe(ji-1,jj-1) = ahmf(ji-1,jj-1,jk)                                                              &
142                     &     * (    e1f(ji-1,jj-1)    * r1_e2f(ji-1,jj-1)                                             &
143                     &         * ( pu(ji-1,jj  ,jk) * r1_e1u(ji-1,jj  )  - pu(ji-1,jj-1,jk) * r1_e1u(ji-1,jj-1) )   &
144                     &         +  e2f(ji-1,jj-1)    * r1_e1f(ji-1,jj-1)                                             &
145                     &         * ( pv(ji  ,jj-1,jk) * r1_e2v(ji  ,jj-1)  - pv(ji-1,jj-1,jk) * r1_e2v(ji-1,jj-1) )   )
146                  !                                      ! tension stress component (T-point)   NB : ahmt has already been multiplied by tmask
147                  zten(ji,jj)    = ahmt(ji,jj,jk)                                                       &
148                     &     * (    e2t(ji,jj)    * r1_e1t(ji,jj)                                         &
149                     &         * ( pu(ji,jj,jk) * r1_e2u(ji,jj)  - pu(ji-1,jj,jk) * r1_e2u(ji-1,jj) )   &
150                     &         -  e1t(ji,jj)    * r1_e2t(ji,jj)                                         &
151                     &         * ( pv(ji,jj,jk) * r1_e1v(ji,jj)  - pv(ji,jj-1,jk) * r1_e1v(ji,jj-1) )   )
152               END_2D
153               !
154               DO_2D_00_00
155                  pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                               &
156                     &    * (   (   zten(ji+1,jj  ) * e2t(ji+1,jj  )*e2t(ji+1,jj  ) * e3t(ji+1,jj  ,jk,Kmm)                       &
157                     &            - zten(ji  ,jj  ) * e2t(ji  ,jj  )*e2t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm) ) * r1_e2u(ji,jj)     &
158                     &        + (   zshe(ji  ,jj  ) * e1f(ji  ,jj  )*e1f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                           &
159                     &            - zshe(ji  ,jj-1) * e1f(ji  ,jj-1)*e1f(ji  ,jj-1) * e3f(ji  ,jj-1,jk)     ) * r1_e1u(ji,jj) )
160                  !
161                  pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                               &
162                     &    * (   (   zshe(ji  ,jj  ) * e2f(ji  ,jj  )*e2f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                           &
163                     &            - zshe(ji-1,jj  ) * e2f(ji-1,jj  )*e2f(ji-1,jj  ) * e3f(ji-1,jj  ,jk)     ) * r1_e2v(ji,jj)     &
164                     &        - (   zten(ji  ,jj+1) * e1t(ji  ,jj+1)*e1t(ji  ,jj+1) * e3t(ji  ,jj+1,jk,Kmm)                       &
165                     &            - zten(ji  ,jj  ) * e1t(ji  ,jj  )*e1t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm) ) * r1_e1v(ji,jj) )
166                   !
167               END_2D
168               !
169            END DO                                           !   End of slab
170            !
171         CASE ( np_dynldf_lap_symN )       !==  Symmetric form  ==!   (naive way)
172            !
173            DO jk = 1, jpkm1                                 ! Horizontal slab
174               !
175               DO_2D_01_01
176                  !                                      ! shearing stress component (F-point)   NB : ahmf has already been multiplied by fmask
177                  zshe(ji-1,jj-1) = ahmf(ji-1,jj-1,jk)                                           &
178                     &     * (   r1_e2f(ji-1,jj-1) * ( pu(ji-1,jj  ,jk) - pu(ji-1,jj-1,jk)  )   &
179                     &         + r1_e1f(ji-1,jj-1) * ( pv(ji  ,jj-1,jk) - pv(ji-1,jj-1,jk)  )   )
180                  !                                      ! tension stress component (T-point)   NB : ahmt has already been multiplied by tmask
181                  zten(ji,jj)    = ahmt(ji,jj,jk)                                       &
182                     &     * (   r1_e1t(ji,jj) * ( pu(ji,jj,jk) - pu(ji-1,jj  ,jk)  )   &
183                     &         - r1_e2t(ji,jj) * ( pv(ji,jj,jk) - pv(ji  ,jj-1,jk)  )   )
184               END_2D
185               !
186               DO_2D_00_00
187                  pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   &
188                     &    * (   zten(ji+1,jj  ) * e2t(ji+1,jj  ) * e3t(ji+1,jj  ,jk,Kmm)              &
189                     &        - zten(ji  ,jj  ) * e2t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm)              &
190                     &        + zshe(ji  ,jj  ) * e1f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                  &
191                     &        - zshe(ji  ,jj-1) * e1f(ji  ,jj-1) * e3f(ji  ,jj-1,jk)                  )
192                  !
193                  pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)   &
194                     &    * (   zshe(ji  ,jj  ) * e2f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                  &
195                     &        - zshe(ji-1,jj  ) * e2f(ji-1,jj  ) * e3f(ji-1,jj  ,jk)                  &
196                     &        - zten(ji  ,jj+1) * e1t(ji  ,jj+1) * e3t(ji  ,jj+1,jk,Kmm)              &
197                     &        + zten(ji  ,jj  ) * e1t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm)              )
198                   !
199               END_2D
200               !
201            END DO                                           !   End of slab
202            !
203         CASE DEFAULT                                     ! error
204            CALL ctl_stop('STOP','dyn_ldf_lap: wrong value for nn_dynldf_lap_typ'  )
205         END SELECT
206         !
207      !
208   END SUBROUTINE dyn_ldf_lap
209
210
211   SUBROUTINE dyn_ldf_blp( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs )
212      !!----------------------------------------------------------------------
213      !!                 ***  ROUTINE dyn_ldf_blp  ***
214      !!
215      !! ** Purpose :   Compute the before lateral momentum viscous trend
216      !!              and add it to the general trend of momentum equation.
217      !!
218      !! ** Method  :   The lateral viscous trends is provided by a bilaplacian
219      !!      operator applied to before field (forward in time).
220      !!      It is computed by two successive calls to dyn_ldf_lap routine
221      !!
222      !! ** Action :   pt(:,:,:,:,Krhs)   updated with the before rotated bilaplacian diffusion
223      !!----------------------------------------------------------------------
224      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
225      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm   ! ocean time level indices
226      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv     ! before velocity fields
227      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! momentum trend
228      !
229      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point
230      !!----------------------------------------------------------------------
231      !
232      IF( kt == nit000 )  THEN
233         IF(lwp) WRITE(numout,*)
234         IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum '
235         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
236      ENDIF
237      !
238      zulap(:,:,:) = 0._wp
239      zvlap(:,:,:) = 0._wp
240      !
241      CALL dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, zulap, zvlap, 1 )   ! rotated laplacian applied to pt (output in zlap,Kbb)
242      !
243      CALL lbc_lnk_multi( 'dynldf_lap_blp', zulap, 'U', -1., zvlap, 'V', -1. )             ! Lateral boundary conditions
244      !
245      CALL dyn_ldf_lap( kt, Kbb, Kmm, zulap, zvlap, pu_rhs, pv_rhs, 2 )   ! rotated laplacian applied to zlap (output in pt(:,:,:,:,Krhs))
246      !
247   END SUBROUTINE dyn_ldf_blp
248
249   !!======================================================================
250END MODULE dynldf_lap_blp
Note: See TracBrowser for help on using the repository browser.