source: NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/DYN/dynldf_lap_blp.F90 @ 13257

Last change on this file since 13257 was 13257, checked in by orioltp, 3 months ago

Updated with trunk at r13245 and small change allocating variables in icb_oce.F90.

  • Property svn:keywords set to Id
File size: 7.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   !!----------------------------------------------------------------------
8
9   !!----------------------------------------------------------------------
10   !!   dyn_ldf_lap   : update the momentum trend with the lateral viscosity using an iso-level   laplacian operator
11   !!   dyn_ldf_blp   : update the momentum trend with the lateral viscosity using an iso-level bilaplacian operator
12   !!----------------------------------------------------------------------
13   USE oce            ! ocean dynamics and tracers
14   USE dom_oce        ! ocean space and time domain
15   USE ldfdyn         ! lateral diffusion: eddy viscosity coef.
16   USE ldfslp         ! iso-neutral slopes
17   USE zdf_oce        ! ocean vertical physics
18   !
19   USE in_out_manager ! I/O manager
20   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
21
22   IMPLICIT NONE
23   PRIVATE
24
25   PUBLIC dyn_ldf_lap  ! called by dynldf.F90
26   PUBLIC dyn_ldf_blp  ! called by dynldf.F90
27
28   !! * Substitutions
29#  include "do_loop_substitute.h90"
30#  include "domzgr_substitute.h90"
31   !!----------------------------------------------------------------------
32   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
33   !! $Id$
34   !! Software governed by the CeCILL license (see ./LICENSE)
35   !!----------------------------------------------------------------------
36CONTAINS
37
38   SUBROUTINE dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs, kpass )
39      !!----------------------------------------------------------------------
40      !!                     ***  ROUTINE dyn_ldf_lap  ***
41      !!                       
42      !! ** Purpose :   Compute the before horizontal momentum diffusive
43      !!      trend and add it to the general trend of momentum equation.
44      !!
45      !! ** Method  :   The Laplacian operator apply on horizontal velocity is
46      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )
47      !!
48      !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv.
49      !!----------------------------------------------------------------------
50      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
51      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm   ! ocean time level indices
52      INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
53      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv     ! before velocity  [m/s]
54      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! velocity trend   [m/s2]
55      !
56      INTEGER  ::   ji, jj, jk   ! dummy loop indices
57      REAL(wp) ::   zsign        ! local scalars
58      REAL(wp) ::   zua, zva     ! local scalars
59      REAL(wp), DIMENSION(jpi,jpj) ::   zcur, zdiv
60      !!----------------------------------------------------------------------
61      !
62      IF( kt == nit000 .AND. lwp ) THEN
63         WRITE(numout,*)
64         WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass
65         WRITE(numout,*) '~~~~~~~ '
66      ENDIF
67      !
68      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign
69      ELSE                    ;   zsign = -1._wp      !  (eddy viscosity coef. >0)
70      ENDIF
71      !
72      !                                                ! ===============
73      DO jk = 1, jpkm1                                 ! Horizontal slab
74         !                                             ! ===============
75         DO_2D_01_01
76            !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1)
77            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
78               &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  &
79               &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  )
80            !                                      ! ahm * div        (computed from 2 to jpi/jpj)
81            zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb)               &   ! ahmt already * by tmask
82               &     * (  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)  &
83               &        + 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)  )
84         END_2D
85         !
86         DO_2D_00_00
87            pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * umask(ji,jj,jk) * (    &    ! * by umask is mandatory for dyn_ldf_blp use
88               &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm)   &
89               &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                      )
90               !
91            pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * vmask(ji,jj,jk) * (    &    ! * by vmask is mandatory for dyn_ldf_blp use
92               &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm)   &
93               &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                      )
94         END_2D
95         !                                             ! ===============
96      END DO                                           !   End of slab
97      !                                                ! ===============
98      !
99   END SUBROUTINE dyn_ldf_lap
100
101
102   SUBROUTINE dyn_ldf_blp( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs )
103      !!----------------------------------------------------------------------
104      !!                 ***  ROUTINE dyn_ldf_blp  ***
105      !!                   
106      !! ** Purpose :   Compute the before lateral momentum viscous trend
107      !!              and add it to the general trend of momentum equation.
108      !!
109      !! ** Method  :   The lateral viscous trends is provided by a bilaplacian
110      !!      operator applied to before field (forward in time).
111      !!      It is computed by two successive calls to dyn_ldf_lap routine
112      !!
113      !! ** Action :   pt(:,:,:,:,Krhs)   updated with the before rotated bilaplacian diffusion
114      !!----------------------------------------------------------------------
115      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
116      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm   ! ocean time level indices
117      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv     ! before velocity fields
118      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pv_rhs   ! momentum trend
119      !
120      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point
121      !!----------------------------------------------------------------------
122      !
123      IF( kt == nit000 )  THEN
124         IF(lwp) WRITE(numout,*)
125         IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum '
126         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
127      ENDIF
128      !
129      zulap(:,:,:) = 0._wp
130      zvlap(:,:,:) = 0._wp
131      !
132      CALL dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, zulap, zvlap, 1 )   ! rotated laplacian applied to pt (output in zlap,Kbb)
133      !
134      CALL lbc_lnk_multi( 'dynldf_lap_blp', zulap, 'U', -1.0_wp, zvlap, 'V', -1.0_wp )             ! Lateral boundary conditions
135      !
136      CALL dyn_ldf_lap( kt, Kbb, Kmm, zulap, zvlap, pu_rhs, pv_rhs, 2 )   ! rotated laplacian applied to zlap (output in pt(:,:,:,:,Krhs))
137      !
138   END SUBROUTINE dyn_ldf_blp
139
140   !!======================================================================
141END MODULE dynldf_lap_blp
Note: See TracBrowser for help on using the repository browser.