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/trunk/src/OCE/DYN – NEMO

source: NEMO/trunk/src/OCE/DYN/dynldf_lap_blp.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

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