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 branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_lap_blp.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago

Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.

  • 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.