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/src/SWE – NEMO

source: NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynldf_lap_blp.F90 @ 12822

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

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