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

source: NEMO/trunk/src/SWE/dynldf_lap_blp.F90 @ 13295

Last change on this file since 13295 was 13295, checked in by acc, 4 years ago

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