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/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/DYN – NEMO

source: NEMO/branches/2021/dev_r14393_HPC-03_Mele_Comm_Cleanup/src/OCE/DYN/dynldf_lap_blp.F90 @ 14801

Last change on this file since 14801 was 14801, checked in by francesca, 3 years ago

add loop fusion to DYN and TRA modules - ticket #2607

  • Property svn:keywords set to Id
File size: 11.8 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
23#if defined key_loop_fusion
24   USE dynldf_lap_blp_lf
25#endif
26   
27   IMPLICIT NONE
28   PRIVATE
29
30   PUBLIC dyn_ldf_lap  ! called by dynldf.F90
31   PUBLIC dyn_ldf_blp  ! called by dynldf.F90
32
33   !! * Substitutions
34#  include "do_loop_substitute.h90"
35#  include "domzgr_substitute.h90"
36   !!----------------------------------------------------------------------
37   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
38   !! $Id$
39   !! Software governed by the CeCILL license (see ./LICENSE)
40   !!----------------------------------------------------------------------
41CONTAINS
42
43   SUBROUTINE dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs, kpass )
44      !!----------------------------------------------------------------------
45      !!                     ***  ROUTINE dyn_ldf_lap  ***
46      !!                       
47      !! ** Purpose :   Compute the before horizontal momentum diffusive
48      !!      trend and add it to the general trend of momentum equation.
49      !!
50      !! ** Method  :   The Laplacian operator apply on horizontal velocity is
51      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )
52      !!
53      !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv.
54      !!
55      !! Reference : S.Griffies, R.Hallberg 2000 Mon.Wea.Rev., DOI:/
56      !!----------------------------------------------------------------------
57      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
58      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm   ! ocean time level indices
59      INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
60      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv     ! before velocity  [m/s]
61      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! velocity trend   [m/s2]
62      !
63      INTEGER  ::   ji, jj, jk   ! dummy loop indices
64      REAL(wp) ::   zsign        ! local scalars
65      REAL(wp) ::   zua, zva     ! local scalars
66      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zcur, zdiv
67      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zten, zshe   ! tension (diagonal) and shearing (anti-diagonal) terms
68      !!----------------------------------------------------------------------
69      !
70#if defined key_loop_fusion
71      CALL dyn_ldf_lap_lf( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs, kpass ) 
72#else
73      IF( kt == nit000 .AND. lwp ) THEN
74         WRITE(numout,*)
75         WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass
76         WRITE(numout,*) '~~~~~~~ '
77      ENDIF
78      !
79      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign
80      ELSE                    ;   zsign = -1._wp      !  (eddy viscosity coef. >0)
81      ENDIF
82      !
83      SELECT CASE( nn_dynldf_typ ) 
84      !             
85      CASE ( np_typ_rot )       !==  Vorticity-Divergence operator  ==!
86         !
87         ALLOCATE( zcur(jpi,jpj) , zdiv(jpi,jpj) )
88         !
89         DO jk = 1, jpkm1                                 ! Horizontal slab
90            !
91            DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
92               !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1)
93               zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &   ! ahmf already * by fmask
94                  &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  &
95                  &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  )
96               !                                      ! ahm * div        (computed from 2 to jpi/jpj)
97               zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)               &   ! ahmt already * by tmask
98                  &     * (  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)  &
99                  &        + 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)  )
100            END_2D
101            !
102            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                       ! - curl( curl) + grad( div )
103               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * umask(ji,jj,jk) * (    &    ! * by umask is mandatory for dyn_ldf_blp use
104                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   &
105                  &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                      )
106               !
107               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * vmask(ji,jj,jk) * (    &    ! * by vmask is mandatory for dyn_ldf_blp use
108                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   &
109                  &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                      )
110            END_2D
111            !
112         END DO                                           !   End of slab
113         !
114         DEALLOCATE( zcur , zdiv )
115         !
116      CASE ( np_typ_sym )       !==  Symmetric operator  ==!
117         !
118         ALLOCATE( zten(jpi,jpj) , zshe(jpi,jpj) )
119         !
120         DO jk = 1, jpkm1                                 ! Horizontal slab
121            !
122            DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
123               !                                      ! shearing stress component (F-point)   NB : ahmf has already been multiplied by fmask
124               zshe(ji-1,jj-1) = ahmf(ji-1,jj-1,jk)                                                              &
125                  &     * (    e1f(ji-1,jj-1)    * r1_e2f(ji-1,jj-1)                                             &
126                  &         * ( pu(ji-1,jj  ,jk) * r1_e1u(ji-1,jj  )  - pu(ji-1,jj-1,jk) * r1_e1u(ji-1,jj-1) )   &
127                  &         +  e2f(ji-1,jj-1)    * r1_e1f(ji-1,jj-1)                                             &
128                  &         * ( pv(ji  ,jj-1,jk) * r1_e2v(ji  ,jj-1)  - pv(ji-1,jj-1,jk) * r1_e2v(ji-1,jj-1) )   ) 
129               !                                      ! tension stress component (T-point)   NB : ahmt has already been multiplied by tmask
130               zten(ji,jj)    = ahmt(ji,jj,jk)                                                       &
131                  &     * (    e2t(ji,jj)    * r1_e1t(ji,jj)                                         &
132                  &         * ( pu(ji,jj,jk) * r1_e2u(ji,jj)  - pu(ji-1,jj,jk) * r1_e2u(ji-1,jj) )   &
133                  &         -  e1t(ji,jj)    * r1_e2t(ji,jj)                                         &
134                  &         * ( pv(ji,jj,jk) * r1_e1v(ji,jj)  - pv(ji,jj-1,jk) * r1_e1v(ji,jj-1) )   )   
135            END_2D
136            !
137            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
138               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)                               &
139                  &    * (   (   zten(ji+1,jj  ) * e2t(ji+1,jj  )*e2t(ji+1,jj  ) * e3t(ji+1,jj  ,jk,Kmm)                       &
140                  &            - zten(ji  ,jj  ) * e2t(ji  ,jj  )*e2t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm) ) * r1_e2u(ji,jj)     &                                                   
141                  &        + (   zshe(ji  ,jj  ) * e1f(ji  ,jj  )*e1f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                           &
142                  &            - zshe(ji  ,jj-1) * e1f(ji  ,jj-1)*e1f(ji  ,jj-1) * e3f(ji  ,jj-1,jk)     ) * r1_e1u(ji,jj) )   
143               !
144               pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)                               &
145                  &    * (   (   zshe(ji  ,jj  ) * e2f(ji  ,jj  )*e2f(ji  ,jj  ) * e3f(ji  ,jj  ,jk)                           &
146                  &            - zshe(ji-1,jj  ) * e2f(ji-1,jj  )*e2f(ji-1,jj  ) * e3f(ji-1,jj  ,jk)     ) * r1_e2v(ji,jj)     &
147                  &        - (   zten(ji  ,jj+1) * e1t(ji  ,jj+1)*e1t(ji  ,jj+1) * e3t(ji  ,jj+1,jk,Kmm)                       &
148                  &            - zten(ji  ,jj  ) * e1t(ji  ,jj  )*e1t(ji  ,jj  ) * e3t(ji  ,jj  ,jk,Kmm) ) * r1_e1v(ji,jj) )
149               !
150            END_2D
151            !
152         END DO
153         !
154         DEALLOCATE( zten , zshe )
155         !
156      END SELECT
157      !
158#endif
159   END SUBROUTINE dyn_ldf_lap
160
161
162   SUBROUTINE dyn_ldf_blp( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs )
163      !!----------------------------------------------------------------------
164      !!                 ***  ROUTINE dyn_ldf_blp  ***
165      !!                   
166      !! ** Purpose :   Compute the before lateral momentum viscous trend
167      !!              and add it to the general trend of momentum equation.
168      !!
169      !! ** Method  :   The lateral viscous trends is provided by a bilaplacian
170      !!      operator applied to before field (forward in time).
171      !!      It is computed by two successive calls to dyn_ldf_lap routine
172      !!
173      !! ** Action :   pt(:,:,:,:,Krhs)   updated with the before rotated bilaplacian diffusion
174      !!----------------------------------------------------------------------
175      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
176      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm   ! ocean time level indices
177      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv     ! before velocity fields
178      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! momentum trend
179      !
180      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point
181      !!----------------------------------------------------------------------
182      !
183#if defined key_loop_fusion
184      CALL dyn_ldf_blp_lf( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs ) 
185#else
186      IF( kt == nit000 )  THEN
187         IF(lwp) WRITE(numout,*)
188         IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum '
189         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
190      ENDIF
191      !
192      zulap(:,:,:) = 0._wp
193      zvlap(:,:,:) = 0._wp
194      !
195      CALL dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, zulap, zvlap, 1 )   ! rotated laplacian applied to pt (output in zlap,Kbb)
196      !
197      IF (nn_hls==1) CALL lbc_lnk( 'dynldf_lap_blp', zulap, 'U', -1.0_wp, zvlap, 'V', -1.0_wp )             ! Lateral boundary conditions
198      !
199      CALL dyn_ldf_lap( kt, Kbb, Kmm, zulap, zvlap, pu_rhs, pv_rhs, 2 )   ! rotated laplacian applied to zlap (output in pt(:,:,:,:,Krhs))
200      !
201#endif
202   END SUBROUTINE dyn_ldf_blp
203
204   !!======================================================================
205END MODULE dynldf_lap_blp
Note: See TracBrowser for help on using the repository browser.