source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap_blp.F90 @ 7753

Last change on this file since 7753 was 7753, checked in by mocavero, 4 years ago

Reverting trunk to remove OpenMP

  • Property svn:keywords set to Id
File size: 8.4 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   USE wrk_nemo       ! Memory Allocation
22   USE timing         ! Timing
23
24   IMPLICIT NONE
25   PRIVATE
26
27   PUBLIC dyn_ldf_lap  ! called by dynldf.F90
28   PUBLIC dyn_ldf_blp  ! called by dynldf.F90
29
30   !! * Substitutions
31#  include "vectopt_loop_substitute.h90"
32   !!----------------------------------------------------------------------
33   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
34   !! $Id$
35   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
36   !!----------------------------------------------------------------------
37CONTAINS
38
39   SUBROUTINE dyn_ldf_lap( kt, pub, pvb, pua, pva, kpass )
40      !!----------------------------------------------------------------------
41      !!                     ***  ROUTINE dyn_ldf_lap  ***
42      !!                       
43      !! ** Purpose :   Compute the before horizontal momentum diffusive
44      !!      trend and add it to the general trend of momentum equation.
45      !!
46      !! ** Method  :   The Laplacian operator apply on horizontal velocity is
47      !!      writen as :   grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) )
48      !!
49      !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb.
50      !!----------------------------------------------------------------------
51      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
52      INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
53      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity  [m/s]
54      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! 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), POINTER, DIMENSION(:,:) ::  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( nn_timing == 1 )   CALL timing_start('dyn_ldf_lap')
69      !
70      CALL wrk_alloc( jpi, jpj, zcur, zdiv ) 
71      !
72      IF( kpass == 1 ) THEN   ;   zsign =  1._wp      ! bilaplacian operator require a minus sign
73      ELSE                    ;   zsign = -1._wp      !  (eddy viscosity coef. >0)
74      ENDIF
75      !
76      !                                                ! ===============
77      DO jk = 1, jpkm1                                 ! Horizontal slab
78         !                                             ! ===============
79         DO jj = 2, jpj
80            DO ji = fs_2, jpi   ! vector opt.
81               !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1)
82!!gm open question here : e3f  at before or now ?    probably now...
83!!gm note that ahmf has already been multiplied by fmask
84               zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f_n(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &
85                  &     * (  e2v(ji  ,jj-1) * pvb(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk)  &
86                  &        - e1u(ji-1,jj  ) * pub(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk)  )
87               !                                      ! ahm * div        (computed from 2 to jpi/jpj)
88!!gm note that ahmt has already been multiplied by tmask
89               zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk)                                         &
90                  &     * (  e2u(ji,jj)*e3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk)  &
91                  &        + e1v(ji,jj)*e3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*e3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk)  )
92            END DO 
93         END DO 
94         !
95         DO jj = 2, jpjm1                             ! - curl( curl) + grad( div )
96            DO ji = fs_2, fs_jpim1   ! vector opt.
97               pua(ji,jj,jk) = pua(ji,jj,jk) + zsign * (                                                 &
98                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk)   &
99                  &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     )
100                  !
101               pva(ji,jj,jk) = pva(ji,jj,jk) + zsign * (                                                 &
102                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk)   &
103                  &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     )
104            END DO
105         END DO
106         !                                             ! ===============
107      END DO                                           !   End of slab
108      !                                                ! ===============
109      CALL wrk_dealloc( jpi, jpj, zcur, zdiv ) 
110      !
111      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_lap')
112      !
113   END SUBROUTINE dyn_ldf_lap
114
115
116   SUBROUTINE dyn_ldf_blp( kt, pub, pvb, pua, pva )
117      !!----------------------------------------------------------------------
118      !!                 ***  ROUTINE dyn_ldf_blp  ***
119      !!                   
120      !! ** Purpose :   Compute the before lateral momentum viscous trend
121      !!              and add it to the general trend of momentum equation.
122      !!
123      !! ** Method  :   The lateral viscous trends is provided by a bilaplacian
124      !!      operator applied to before field (forward in time).
125      !!      It is computed by two successive calls to dyn_ldf_lap routine
126      !!
127      !! ** Action :   pta   updated with the before rotated bilaplacian diffusion
128      !!----------------------------------------------------------------------
129      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
130      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity fields
131      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! momentum trend
132      !
133      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zulap, zvlap   ! laplacian at u- and v-point
134      !!----------------------------------------------------------------------
135      !
136      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_blp')
137      !
138      CALL wrk_alloc( jpi, jpj, jpk, zulap, zvlap ) 
139      !
140      IF( kt == nit000 )  THEN
141         IF(lwp) WRITE(numout,*)
142         IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum '
143         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
144      ENDIF
145      !
146      zulap(:,:,:) = 0._wp
147      zvlap(:,:,:) = 0._wp
148      !
149      CALL dyn_ldf_lap( kt, pub, pvb, zulap, zvlap, 1 )   ! rotated laplacian applied to ptb (output in zlap)
150      !
151      CALL lbc_lnk( zulap(:,:,:) , 'U', -1. )             ! Lateral boundary conditions
152      CALL lbc_lnk( zvlap(:,:,:) , 'V', -1. )
153      !
154      CALL dyn_ldf_lap( kt, zulap, zvlap, pua, pva, 2 )   ! rotated laplacian applied to zlap (output in pta)
155      !
156      CALL wrk_dealloc( jpi, jpj, jpk, zulap, zvlap ) 
157      !
158      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_blp')
159      !
160   END SUBROUTINE dyn_ldf_blp
161
162   !!======================================================================
163END MODULE dynldf_lap_blp
Note: See TracBrowser for help on using the repository browser.