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/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynldf_lap_blp.F90 @ 10806

Last change on this file since 10806 was 10806, checked in by davestorkey, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps branch: Latest updates. Make sure all time-dependent 3D variables are converted in previously modified modules.
  • Property svn:keywords set to Id
File size: 8.2 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
29#  include "vectopt_loop_substitute.h90"
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
[10806]37   SUBROUTINE dyn_ldf_lap( kt, ktlev1, ktlev2, pu, pv, pu_rhs, pva_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      !!
[10806]47      !! ** Action : - pu_rhs, pva_rhs increased by the harmonic operator applied on pu, pv.
[3]48      !!----------------------------------------------------------------------
[10806]49      INTEGER                         , INTENT(in   ) ::   kt              ! ocean time-step index
50      INTEGER                         , INTENT(in   ) ::   ktlev1, ktlev2  ! time level index for scale factors
[5777]51      INTEGER                         , INTENT(in   ) ::   kpass      ! =1/2 first or second passage
[10806]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, pva_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         !                                             ! ===============
[5777]74         DO jj = 2, jpj
75            DO ji = fs_2, jpi   ! vector opt.
76               !                                      ! ahm * e3 * curl  (computed from 1 to jpim1/jpjm1)
[6140]77!!gm open question here : e3f  at before or now ?    probably now...
[5860]78!!gm note that ahmf has already been multiplied by fmask
[10806]79               zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1)       &
80                  &     * (  e2v(ji  ,jj-1) * pv(ji  ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk)  &
81                  &        - e1u(ji-1,jj  ) * pu(ji-1,jj  ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk)  )
[5777]82               !                                      ! ahm * div        (computed from 2 to jpi/jpj)
[6140]83!!gm note that ahmt has already been multiplied by tmask
[10806]84               zdiv(ji,jj)     = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev1)                                         &
85                  &     * (  e2u(ji,jj)*e3u(ji,jj,jk,ktlev1) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,ktlev1) * pu(ji-1,jj,jk)  &
86                  &        + e1v(ji,jj)*e3v(ji,jj,jk,ktlev1) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,ktlev1) * pv(ji,jj-1,jk)  )
[5777]87            END DO 
88         END DO 
89         !
90         DO jj = 2, jpjm1                             ! - curl( curl) + grad( div )
[3]91            DO ji = fs_2, fs_jpim1   ! vector opt.
[10806]92               pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * (                                                 &
93                  &              - ( zcur(ji  ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,ktlev2)   &
[5777]94                  &              + ( zdiv(ji+1,jj) - zdiv(ji,jj  ) ) * r1_e1u(ji,jj)                     )
95                  !
[10806]96               pva_rhs(ji,jj,jk) = pva_rhs(ji,jj,jk) + zsign * (                                                 &
97                  &                ( zcur(ji,jj  ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,ktlev2)   &
[5777]98                  &              + ( zdiv(ji,jj+1) - zdiv(ji  ,jj) ) * r1_e2v(ji,jj)                     )
[3]99            END DO
100         END DO
101         !                                             ! ===============
102      END DO                                           !   End of slab
103      !                                                ! ===============
[5777]104      !
[3]105   END SUBROUTINE dyn_ldf_lap
106
[5777]107
[10806]108   SUBROUTINE dyn_ldf_blp( kt, ktlev1, ktlev2, pu, pv, pu_rhs, pva_rhs )
[5777]109      !!----------------------------------------------------------------------
110      !!                 ***  ROUTINE dyn_ldf_blp  ***
111      !!                   
112      !! ** Purpose :   Compute the before lateral momentum viscous trend
113      !!              and add it to the general trend of momentum equation.
114      !!
115      !! ** Method  :   The lateral viscous trends is provided by a bilaplacian
116      !!      operator applied to before field (forward in time).
117      !!      It is computed by two successive calls to dyn_ldf_lap routine
118      !!
[10806]119      !! ** Action :   pt_rhs   updated with the before rotated bilaplacian diffusion
[5777]120      !!----------------------------------------------------------------------
121      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
[10806]122      INTEGER                         , INTENT(in   ) ::   ktlev1, ktlev2  ! time level index for scale factors
123      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu, pv   ! before velocity fields
124      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu_rhs, pva_rhs   ! momentum trend
[5777]125      !
[9019]126      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zulap, zvlap   ! laplacian at u- and v-point
[5777]127      !!----------------------------------------------------------------------
128      !
129      IF( kt == nit000 )  THEN
130         IF(lwp) WRITE(numout,*)
131         IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum '
132         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
133      ENDIF
134      !
[7753]135      zulap(:,:,:) = 0._wp
136      zvlap(:,:,:) = 0._wp
[5777]137      !
[10806]138      CALL dyn_ldf_lap( kt, ktlev1, ktlev2, pu, pv, zulap, zvlap, 1 )   ! rotated laplacian applied to pt (output in zlap)
[5777]139      !
[10425]140      CALL lbc_lnk_multi( 'dynldf_lap_blp', zulap, 'U', -1., zvlap, 'V', -1. )             ! Lateral boundary conditions
[5777]141      !
[10806]142      CALL dyn_ldf_lap( kt, ktlev1, ktlev2, zulap, zvlap, pu_rhs, pva_rhs, 2 )   ! rotated laplacian applied to zlap (output in pt_rhs)
[5777]143      !
144   END SUBROUTINE dyn_ldf_blp
145
[3]146   !!======================================================================
[5777]147END MODULE dynldf_lap_blp
Note: See TracBrowser for help on using the repository browser.