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

Last change on this file since 5861 was 5861, checked in by gm, 5 years ago

#1612 : correct a bug in dynldf_lap_blp in vvl case: forgot en e3t_b

  • Property svn:keywords set to Id
File size: 8.9 KB
Line 
1MODULE dynldf_lap_blp
2   !!======================================================================
3   !!                   ***  MODULE  dynldf_lap_blp  ***
4   !! Ocean dynamics:  lateral viscosity trend (laplacian and bilaplacian)
5   !!======================================================================
6   !! History :  OPA  ! 1990-09 (G. Madec) Original code
7   !!            4.0  ! 1991-11 (G. Madec)
8   !!            6.0  ! 1996-01 (G. Madec) statement function for e3 and ahm
9   !!   NEMO     1.0  ! 2002-06 (G. Madec)  F90: Free form and module
10   !!             -   ! 2004-08 (C. Talandier) New trends organization
11   !!            3.7  ! 2014-01  (F. Lemarie, G. Madec)  restructuration/simplification of ahm specification,
12   !!                 !                                  add velocity dependent coefficient and optional read in file
13   !!----------------------------------------------------------------------
14
15   !!----------------------------------------------------------------------
16   !!   dyn_ldf_lap   : update the momentum trend with the lateral viscosity using an iso-level   laplacian operator
17   !!   dyn_ldf_blp   : update the momentum trend with the lateral viscosity using an iso-level bilaplacian operator
18   !!----------------------------------------------------------------------
19   USE oce            ! ocean dynamics and tracers
20   USE dom_oce        ! ocean space and time domain
21   USE ldfdyn         ! lateral diffusion: eddy viscosity coef.
22   USE ldfslp         ! iso-neutral slopes
23   USE zdf_oce        ! ocean vertical physics
24   !
25   USE in_out_manager ! I/O manager
26   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
27   USE wrk_nemo       ! Memory Allocation
28   USE timing         ! Timing
29
30   IMPLICIT NONE
31   PRIVATE
32
33   PUBLIC dyn_ldf_lap  ! called by dynldf.F90
34   PUBLIC dyn_ldf_blp  ! called by dynldf.F90
35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40   !! NEMO/OPA 3.7 , NEMO Consortium (2014)
41   !! $Id$
42   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
43   !!----------------------------------------------------------------------
44CONTAINS
45
46   SUBROUTINE dyn_ldf_lap( kt, pub, pvb, pua, pva, 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      !!
56      !! ** Action : - pua, pva increased by the harmonic operator applied on pub, pvb.
57      !!----------------------------------------------------------------------
58      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
59      INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
60      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity  [m/s]
61      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! 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), POINTER, DIMENSION(:,:) ::  zcur, zdiv
67      !!----------------------------------------------------------------------
68      !
69      IF( kt == nit000 .AND. lwp ) THEN
70         WRITE(numout,*)
71         WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass
72         WRITE(numout,*) '~~~~~~~ '
73      ENDIF
74      !
75      IF( nn_timing == 1 )   CALL timing_start('dyn_ldf_lap')
76      !
77      CALL wrk_alloc( jpi, jpj, zcur, zdiv ) 
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      !                                                ! ===============
84      DO jk = 1, jpkm1                                 ! Horizontal slab
85         !                                             ! ===============
86         DO jj = 2, jpj
87            DO ji = fs_2, jpi   ! vector opt.
88               !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1)
89!!gm open question here : fse3f  at before or now ?    probably now...
90!!gm note that ahmf has already been multiplied by fmask
91               zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * fse3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &
92                  &     * (  e2v(ji  ,jj-1) * pvb(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pvb(ji-1,jj-1,jk)  &
93                  &        - e1u(ji-1,jj  ) * pub(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pub(ji-1,jj-1,jk)  ) 
94               !                                      ! ahm * div        (computed from 2 to jpi/jpj)
95               zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / fse3t_b(ji,jj,jk) * tmask(ji,jj,jk)                           &
96                  &     * (  e2u(ji,jj)*fse3u_b(ji,jj,jk) * pub(ji,jj,jk) - e2u(ji-1,jj)*fse3u_b(ji-1,jj,jk) * pub(ji-1,jj,jk)  &
97                  &        + e1v(ji,jj)*fse3v_b(ji,jj,jk) * pvb(ji,jj,jk) - e1v(ji,jj-1)*fse3v_b(ji,jj-1,jk) * pvb(ji,jj-1,jk)  )
98            END DO 
99         END DO 
100         !
101         DO jj = 2, jpjm1                             ! - curl( curl) + grad( div )
102            DO ji = fs_2, fs_jpim1   ! vector opt.
103               pua(ji,jj,jk) = pua(ji,jj,jk) + zsign * (                                                   &
104                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) /  ( e2u(ji,jj) * fse3u(ji,jj,jk) )   &
105                  &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     )
106                  !
107               pva(ji,jj,jk) = pva(ji,jj,jk) + zsign * (                                                   &
108                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) /  ( e1v(ji,jj) * fse3v(ji,jj,jk) )   &
109                  &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     )
110            END DO
111         END DO
112         !                                             ! ===============
113      END DO                                           !   End of slab
114      !                                                ! ===============
115      CALL wrk_dealloc( jpi, jpj, zcur, zdiv ) 
116      !
117      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_lap')
118      !
119   END SUBROUTINE dyn_ldf_lap
120
121
122   SUBROUTINE dyn_ldf_blp( kt, pub, pvb, pua, pva )
123      !!----------------------------------------------------------------------
124      !!                 ***  ROUTINE dyn_ldf_blp  ***
125      !!                   
126      !! ** Purpose :   Compute the before lateral momentum viscous trend
127      !!              and add it to the general trend of momentum equation.
128      !!
129      !! ** Method  :   The lateral viscous trends is provided by a bilaplacian
130      !!      operator applied to before field (forward in time).
131      !!      It is computed by two successive calls to dyn_ldf_lap routine
132      !!
133      !! ** Action :   pta   updated with the before rotated bilaplacian diffusion
134      !!----------------------------------------------------------------------
135      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
136      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pub, pvb   ! before velocity fields
137      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pua, pva   ! momentum trend
138      !
139      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zulap, zvlap   ! laplacian at u- and v-point
140      !!----------------------------------------------------------------------
141      !
142      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_blp')
143      !
144      CALL wrk_alloc( jpi, jpj, jpk, zulap, zvlap ) 
145      !
146      IF( kt == nit000 )  THEN
147         IF(lwp) WRITE(numout,*)
148         IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum '
149         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
150      ENDIF
151      !
152      zulap(:,:,:) = 0._wp
153      zvlap(:,:,:) = 0._wp
154      !
155      CALL dyn_ldf_lap( kt, pub, pvb, zulap, zvlap, 1 )   ! rotated laplacian applied to ptb (output in zlap)
156      !
157      CALL lbc_lnk( zulap(:,:,:) , 'U', -1. )             ! Lateral boundary conditions
158      CALL lbc_lnk( zvlap(:,:,:) , 'V', -1. )
159      !
160      CALL dyn_ldf_lap( kt, zulap, zvlap, pua, pva, 2 )   ! rotated laplacian applied to zlap (output in pta)
161      !
162      CALL wrk_dealloc( jpi, jpj, jpk, zulap, zvlap ) 
163      !
164      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_blp')
165      !
166   END SUBROUTINE dyn_ldf_blp
167
168   !!======================================================================
169END MODULE dynldf_lap_blp
Note: See TracBrowser for help on using the repository browser.