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_lf.F90 in NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DYN – NEMO

source: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DYN/dynldf_lap_blp_lf.F90 @ 15540

Last change on this file since 15540 was 15540, checked in by sparonuz, 3 years ago

Mixed precision version, tested up to 30 years on ORCA2.
File size: 15.1 KB
Line 
1MODULE dynldf_lap_blp_lf
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 domutl, ONLY : is_tile
17   USE ldfdyn         ! lateral diffusion: eddy viscosity coef.
18   USE ldfslp         ! iso-neutral slopes
19   USE zdf_oce        ! ocean vertical physics
20   !
21   USE in_out_manager ! I/O manager
22   USE lib_mpp
23   
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC dyn_ldf_lap_lf  ! called by dynldf.F90
28   PUBLIC dyn_ldf_blp_lf  ! called by dynldf.F90
29
30   !! * Substitutions
31#  include "do_loop_substitute.h90"
32#  include "single_precision_substitute.h90"
33#  include "domzgr_substitute.h90"
34   !!----------------------------------------------------------------------
35   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
36   !! $Id: dynldf_lap_blp.F90 14757 2021-04-27 15:33:44Z francesca $
37   !! Software governed by the CeCILL license (see ./LICENSE)
38   !!----------------------------------------------------------------------
39CONTAINS
40
41   SUBROUTINE dyn_ldf_lap_lf( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs, kpass )
42      !!
43      INTEGER                   , INTENT(in   ) ::   kt               ! ocean time-step index
44      INTEGER                   , INTENT(in   ) ::   Kbb, Kmm         ! ocean time level indices
45      INTEGER                   , INTENT(in   ) ::   kpass            ! =1/2 first or second passage
46      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pu, pv           ! before velocity  [m/s]
47      REAL(dp), DIMENSION(:,:,:), INTENT(inout) ::   pu_rhs, pv_rhs   ! velocity trend   [m/s2]
48      !!
49      CALL dyn_ldf_lap_lf_t( kt, Kbb, Kmm, pu, pv, is_tile(pu), pu_rhs, pv_rhs, is_tile(pu_rhs), kpass )
50
51   END SUBROUTINE dyn_ldf_lap_lf
52 
53   SUBROUTINE dyn_ldf_lap_lf_t( kt, Kbb, Kmm, pu, pv, ktuv, pu_rhs, pv_rhs, ktuv_rhs, kpass )
54      !!----------------------------------------------------------------------
55      !!                     ***  ROUTINE dyn_ldf_lap  ***
56      !!                       
57      !! ** Purpose :   Compute the before horizontal momentum diffusive
58      !!      trend and add it to the general trend of momentum equation.
59      !!
60      !! ** Method  :   The Laplacian operator apply on horizontal velocity is
61      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )
62      !!
63      !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv.
64      !!
65      !! Reference : S.Griffies, R.Hallberg 2000 Mon.Wea.Rev., DOI:/
66      !!----------------------------------------------------------------------
67      INTEGER                                 , INTENT(in   ) ::   kt         ! ocean time-step index
68      INTEGER                                 , INTENT(in   ) ::   Kbb, Kmm   ! ocean time level indices
69      INTEGER                                 , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
70      INTEGER                                 , INTENT(in   ) ::   ktuv, ktuv_rhs
71      REAL(wp), DIMENSION(A2D_T(ktuv)    ,JPK), INTENT(in   ) ::   pu, pv ! before velocity  [m/s]
72      REAL(dp), DIMENSION(A2D_T(ktuv_rhs),JPK), INTENT(inout) ::   pu_rhs, pv_rhs   ! velocity trend   [m/s2]
73      !
74      INTEGER  ::   ji, jj, jk   ! dummy loop indices
75      INTEGER  ::   iij
76      REAL(dp) ::   zsign        ! local scalars
77      REAL(dp) ::   zcur, zcur_im1, zcur_jm1     ! local scalars
78      REAL(dp) ::   zdiv, zdiv_ip1, zdiv_jp1     ! local scalars
79      REAL(dp) ::   zten, zten_ip1, zten_jp1, zshe, zshe_im1, zshe_jm1  ! tension (diagonal) and shearing (anti-diagonal) terms
80      !!----------------------------------------------------------------------
81      !
82      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
83         IF( kt == nit000 .AND. lwp ) THEN
84            WRITE(numout,*)
85            WRITE(numout,*) 'dyn_ldf_lf : iso-level harmonic (laplacian) operator, pass=', kpass
86            WRITE(numout,*) '~~~~~~~ '
87         ENDIF
88      ENDIF
89      !
90      ! Define pu_rhs/pv_rhs halo points for multi-point haloes in bilaplacian case
91      IF( nldf_dyn == np_blp .AND. kpass == 1 ) THEN ; iij = nn_hls
92      ELSE                                           ; iij = 1
93      ENDIF
94      !
95      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign
96      ELSE                    ;   zsign = -1._wp      !  (eddy viscosity coef. >0)
97      ENDIF
98      !
99      SELECT CASE( nn_dynldf_typ ) 
100      !             
101      CASE ( np_typ_rot )       !==  Vorticity-Divergence operator  ==!
102         !
103         DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 )                           ! Horizontal slab
104            !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1)
105            zcur     = ahmf(ji,jj,jk) * e3f(ji,jj,jk) * r1_e1e2f(ji,jj)               &   ! ahmf already * by fmask   
106               &       * ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk)  &
107               &       - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) )
108            zcur_jm1 = ahmf(ji  ,jj-1,jk) * e3f(ji,jj-1,jk) * r1_e1e2f(ji,jj-1)               &   ! ahmf already * by fmask
109               &       * ( e2v(ji+1,jj-1) * pv(ji+1,jj-1,jk) - e2v(ji,jj-1) * pv(ji,jj-1,jk)  &
110               &       - e1u(ji,jj) * pu(ji,jj,jk) + e1u(ji,jj-1) * pu(ji,jj-1,jk) )
111            zcur_im1 = ahmf(ji-1,jj,jk) * e3f(ji-1,jj,jk) * r1_e1e2f(ji-1,jj)         &   ! ahmf already * by fmask
112               &       * ( e2v(ji,jj) * pv(ji,jj,jk) - e2v(ji-1,jj) * pv(ji-1,jj,jk)  &
113               &       - e1u(ji-1,jj+1) * pu(ji-1,jj+1,jk) + e1u(ji-1,jj) * pu(ji-1,jj,jk) )
114            !                                      ! ahm * div        (computed from 2 to jpi/jpj)
115            zdiv     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)               &   ! ahmt already * by tmask
116               &     * ( 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)  &
117               &       + 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)  )
118            zdiv_ip1 = ahmt(ji+1,jj,jk) * r1_e1e2t(ji+1,jj) / e3t(ji+1,jj,jk,Kbb)         &   ! ahmt already * by tmask
119               &     * ( e2u(ji+1,jj)*e3u(ji+1,jj,jk,Kbb) * pu(ji+1,jj,jk) - e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk)  &
120               &       + e1v(ji+1,jj)*e3v(ji+1,jj,jk,Kbb) * pv(ji+1,jj,jk) - e1v(ji+1,jj-1)*e3v(ji+1,jj-1,jk,Kbb) * pv(ji+1,jj-1,jk) )
121            zdiv_jp1 = ahmt(ji,jj+1,jk) * r1_e1e2t(ji,jj+1) / e3t(ji,jj+1,jk,Kbb)         &   ! ahmt already * by tmask
122               &     * ( e2u(ji,jj+1)*e3u(ji,jj+1,jk,Kbb) * pu(ji,jj+1,jk) - e2u(ji-1,jj+1)*e3u(ji-1,jj+1,jk,Kbb) * pu(ji-1,jj+1,jk)  &
123               &       + e1v(ji,jj+1)*e3v(ji,jj+1,jk,Kbb) * pv(ji,jj+1,jk) - e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) )
124            !                                      ! - curl( curl) + grad( div )
125            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * umask(ji,jj,jk) * (    &    ! * by umask is mandatory for dyn_ldf_blp use
126               &              - ( zcur - zcur_jm1 ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   &
127               &              + ( zdiv_ip1 - zdiv ) * r1_e1u(ji,jj)                      )
128            !
129            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * vmask(ji,jj,jk) * (    &    ! * by vmask is mandatory for dyn_ldf_blp use
130               &                ( zcur - zcur_im1 ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   &
131               &              + ( zdiv_jp1 - zdiv ) * r1_e2v(ji,jj)                      )
132         END_3D                                                  !   End of slab
133         !
134      CASE ( np_typ_sym )       !==  Symmetric operator  ==!
135         !
136         DO_3D( iij-1, iij-1, iij-1, iij-1, 1, jpkm1 )                           ! Horizontal slab
137            !                                      ! shearing stress component (F-point)   NB : ahmf has already been multiplied by fmask
138            zshe = ahmf(ji,jj,jk)                                                            &
139               &    * (    e1f(ji,jj)    * r1_e2f(ji,jj)                                     &
140               &    * ( pu(ji,jj+1,jk) * r1_e1u(ji,jj+1)  - pu(ji,jj,jk) * r1_e1u(ji,jj) )   &
141               &    +  e2f(ji,jj)    * r1_e1f(ji,jj)                                         &
142               &    * ( pv(ji+1,jj,jk) * r1_e2v(ji+1,jj)  - pv(ji,jj,jk) * r1_e2v(ji,jj) )   ) 
143            zshe_im1 = ahmf(ji-1,jj,jk)                                                                  &
144               &        * (    e1f(ji-1,jj)    * r1_e2f(ji-1,jj)                                         &
145               &        * ( pu(ji-1,jj+1,jk) * r1_e1u(ji-1,jj+1)  - pu(ji-1,jj,jk) * r1_e1u(ji-1,jj) )   &
146               &        +  e2f(ji-1,jj)    * r1_e1f(ji-1,jj)                                             &
147               &        * ( pv(ji  ,jj,jk) * r1_e2v(ji  ,jj)  - pv(ji-1,jj,jk) * r1_e2v(ji-1,jj) )   ) 
148            zshe_jm1 = ahmf(ji,jj-1,jk)                                                                &
149               &        * (    e1f(ji,jj-1)    * r1_e2f(ji,jj-1)                                       &
150               &        * ( pu(ji,jj,jk) * r1_e1u(ji,jj)  - pu(ji,jj-1,jk) * r1_e1u(ji,jj-1) )         &
151               &        +  e2f(ji,jj-1)    * r1_e1f(ji,jj-1)                                           &
152               &        * ( pv(ji+1,jj-1,jk) * r1_e2v(ji+1,jj-1)  - pv(ji,jj-1,jk) * r1_e2v(ji,jj-1) )   ) 
153            !                                      ! tension stress component (T-point)   NB : ahmt has already been multiplied by tmask
154            zten        = ahmt(ji,jj,jk)                                                          &
155               &     * (    e2t(ji,jj)    * r1_e1t(ji,jj)                                         &
156               &         * ( pu(ji,jj,jk) * r1_e2u(ji,jj)  - pu(ji-1,jj,jk) * r1_e2u(ji-1,jj) )   &
157               &         -  e1t(ji,jj)    * r1_e2t(ji,jj)                                         &
158               &         * ( pv(ji,jj,jk) * r1_e1v(ji,jj)  - pv(ji,jj-1,jk) * r1_e1v(ji,jj-1) )   )   
159            zten_ip1    = ahmt(ji+1,jj,jk)                                                              &
160               &     * (    e2t(ji+1,jj)    * r1_e1t(ji+1,jj)                                           &
161               &         * ( pu(ji+1,jj,jk) * r1_e2u(ji+1,jj)  - pu(ji,jj,jk) * r1_e2u(ji,jj) )         &
162               &         -  e1t(ji+1,jj)    * r1_e2t(ji+1,jj)                                           &
163               &         * ( pv(ji+1,jj,jk) * r1_e1v(ji+1,jj)  - pv(ji+1,jj-1,jk) * r1_e1v(ji+1,jj-1) ) )   
164            zten_jp1    = ahmt(ji,jj+1,jk)                                                              &
165               &     * (    e2t(ji,jj+1)    * r1_e1t(ji,jj+1)                                           &
166               &         * ( pu(ji,jj+1,jk) * r1_e2u(ji,jj+1)  - pu(ji-1,jj+1,jk) * r1_e2u(ji-1,jj+1) ) &
167               &         -  e1t(ji,jj+1)    * r1_e2t(ji,jj+1)                                           &
168               &         * ( pv(ji,jj+1,jk) * r1_e1v(ji,jj+1)  - pv(ji,jj,jk) * r1_e1v(ji,jj) )   )   
169            !
170            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                    &
171               &    * (   (   zten_ip1 * e2t(ji+1,jj  )*e2t(ji+1,jj  ) * e3t(ji+1,jj  ,jk,Kmm)                   &
172               &            - zten * e2t(ji  ,jj  )*e2t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm) ) * r1_e2u(ji,jj)     & 
173               &        + (   zshe * e1f(ji  ,jj  )*e1f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                           &
174               &            - zshe_jm1 * e1f(ji  ,jj-1)*e1f(ji  ,jj-1) * e3f(ji  ,jj-1,jk)     ) * r1_e1u(ji,jj) )   
175            !
176            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                    &
177               &    * (   (   zshe * e2f(ji  ,jj  )*e2f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                           &
178               &            - zshe_im1 * e2f(ji-1,jj  )*e2f(ji-1,jj  ) * e3f(ji-1,jj  ,jk)     ) * r1_e2v(ji,jj) &
179               &        - (   zten_jp1 * e1t(ji  ,jj+1)*e1t(ji  ,jj+1) * e3t(ji  ,jj+1,jk,Kmm)                   &
180               &            - zten * e1t(ji  ,jj  )*e1t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm) ) * r1_e1v(ji,jj) )
181            !
182         END_3D
183         !
184      END SELECT
185      !
186   END SUBROUTINE dyn_ldf_lap_lf_t
187
188
189   SUBROUTINE dyn_ldf_blp_lf( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs )
190      !!----------------------------------------------------------------------
191      !!                 ***  ROUTINE dyn_ldf_blp  ***
192      !!                   
193      !! ** Purpose :   Compute the before lateral momentum viscous trend
194      !!              and add it to the general trend of momentum equation.
195      !!
196      !! ** Method  :   The lateral viscous trends is provided by a bilaplacian
197      !!      operator applied to before field (forward in time).
198      !!      It is computed by two successive calls to dyn_ldf_lap routine
199      !!
200      !! ** Action :   pt(:,:,:,:,Krhs)   updated with the before rotated bilaplacian diffusion
201      !!----------------------------------------------------------------------
202      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
203      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm   ! ocean time level indices
204      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv     ! before velocity fields
205      REAL(dp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! momentum trend
206      !
207      REAL(dp), DIMENSION(A2D(nn_hls),jpk) ::   zulap, zvlap   ! laplacian at u- and v-point
208      !!----------------------------------------------------------------------
209      !
210      IF( kt == nit000 )  THEN
211         IF(lwp) WRITE(numout,*)
212         IF(lwp) WRITE(numout,*) 'dyn_ldf_blp_lf : bilaplacian operator momentum '
213         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
214      ENDIF
215      !
216      zulap(:,:,:) = 0._wp
217      zvlap(:,:,:) = 0._wp
218      !
219      CALL dyn_ldf_lap_lf( kt, Kbb, Kmm, pu, pv, zulap, zvlap, 1 )   ! rotated laplacian applied to pt (output in zlap,Kbb)
220      !
221      CALL dyn_ldf_lap_lf( kt, Kbb, Kmm, CASTSP(zulap), CASTSP(zvlap), pu_rhs, pv_rhs, 2 )   ! rotated laplacian applied to zlap (output in pt(:,:,:,:,Krhs))
222      !
223   END SUBROUTINE dyn_ldf_blp_lf
224
225   !!======================================================================
226END MODULE dynldf_lap_blp_lf
Note: See TracBrowser for help on using the repository browser.