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_bilap.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90 @ 4460

Last change on this file since 4460 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 9.8 KB
RevLine 
[3]1MODULE dynldf_bilap
2   !!======================================================================
3   !!                     ***  MODULE  dynldf_bilap  ***
4   !! Ocean dynamics:  lateral viscosity trend
5   !!======================================================================
[2715]6   !! History :  OPA  ! 1990-09  (G. Madec)  Original code
7   !!            4.0  ! 1993-03  (M. Guyon)  symetrical conditions (M. Guyon)
8   !!            6.0  ! 1996-01  (G. Madec)  statement function for e3
9   !!            8.0  ! 1997-07  (G. Madec)  lbc calls
10   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90: Free form and module
11   !!            2.0  ! 2004-08  (C. Talandier) New trends organization
12   !!----------------------------------------------------------------------
[3]13
14   !!----------------------------------------------------------------------
15   !!   dyn_ldf_bilap : update the momentum trend with the lateral diffusion
16   !!                   using an iso-level bilaplacian operator
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and tracers
19   USE dom_oce         ! ocean space and time domain
20   USE ldfdyn_oce      ! ocean dynamics: lateral physics
21   USE in_out_manager  ! I/O manager
[216]22   USE trdmod          ! ocean dynamics trends
23   USE trdmod_oce      ! ocean variables trends
[3]24   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
25
26   IMPLICIT NONE
27   PRIVATE
28
[2715]29   PUBLIC   dyn_ldf_bilap   ! called by step.F90
[3]30
[3211]31   !! * Control permutation of array indices
32#  include "oce_ftrans.h90"
33#  include "dom_oce_ftrans.h90"
34#  include "ldfdyn_oce_ftrans.h90"
35
[3]36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "ldfdyn_substitute.h90"
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
[2528]41   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1152]42   !! $Id$
[2715]43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
[3]44   !!----------------------------------------------------------------------
45CONTAINS
46
47   SUBROUTINE dyn_ldf_bilap( kt )
48      !!----------------------------------------------------------------------
49      !!                     ***  ROUTINE dyn_ldf_bilap  ***
50      !!
51      !! ** Purpose :   Compute the before trend of the lateral momentum
52      !!      diffusion and add it to the general trend of momentum equation.
53      !!
54      !! ** Method  :   The before horizontal momentum diffusion trend is a
55      !!      bi-harmonic operator (bilaplacian type) which separates the
56      !!      divergent and rotational parts of the flow.
57      !!      Its horizontal components are computed as follow:
58      !!      laplacian:
59      !!          zlu = 1/e1u di[ hdivb ] - 1/(e2u*e3u) dj-1[ e3f rotb ]
60      !!          zlv = 1/e2v dj[ hdivb ] + 1/(e1v*e3v) di-1[ e3f rotb ]
61      !!      third derivative:
62      !!       * multiply by the eddy viscosity coef. at u-, v-point, resp.
63      !!          zlu = ahmu * zlu
64      !!          zlv = ahmv * zlv
65      !!       * curl and divergence of the laplacian
66      !!          zuf = 1/(e1f*e2f) ( di[e2v zlv] - dj[e1u zlu] )
67      !!          zut = 1/(e1t*e2t*e3t) ( di[e2u*e3u zlu] + dj[e1v*e3v zlv] )
68      !!      bilaplacian:
69      !!              diffu = 1/e1u di[ zut ] - 1/(e2u*e3u) dj-1[ e3f zuf ]
70      !!              diffv = 1/e2v dj[ zut ] + 1/(e1v*e3v) di-1[ e3f zuf ]
[455]71      !!      If ln_sco=F and ln_zps=F, the vertical scale factors in the
[3]72      !!      rotational part of the diffusion are simplified
73      !!      Add this before trend to the general trend (ua,va):
74      !!            (ua,va) = (ua,va) + (diffu,diffv)
75      !!      'key_trddyn' defined: the two components of the horizontal
76      !!                               diffusion trend are saved.
77      !!
78      !! ** Action : - Update (ua,va) with the before iso-level biharmonic
79      !!               mixing trend.
80      !!----------------------------------------------------------------------
[2715]81      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
[3211]82      USE wrk_nemo, ONLY:   zcu => wrk_2d_1 , zcv => wrk_2d_2   ! 2D workspace
[2715]83      USE wrk_nemo, ONLY:   zuf => wrk_3d_3 , zut => wrk_3d_4   ! 3D workspace
84      USE wrk_nemo, ONLY:   zlu => wrk_3d_5 , zlv => wrk_3d_6
[3211]85      !! DCSE_NEMO: need additional directives for renamed module variables
86!FTRANS zuf :I :I :z
87!FTRANS zut :I :I :z
88!FTRANS zlu :I :I :z
89!FTRANS zlv :I :I :z
[2715]90      !
91      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
92      !
93      INTEGER  ::   ji, jj, jk                  ! dummy loop indices
94      REAL(wp) ::   zua, zva, zbt, ze2u, ze2v   ! temporary scalar
[3]95      !!----------------------------------------------------------------------
96
[2715]97      IF( wrk_in_use(2, 1,2) .OR. wrk_in_use(3, 3,4,5,6) ) THEN
98         CALL ctl_stop('dyn_ldf_bilap : requested workspace arrays unavailable')   ;   RETURN
[3]99      ENDIF
[216]100
[2715]101      IF( kt == nit000 .AND. lwp ) THEN
102         WRITE(numout,*)
103         WRITE(numout,*) 'dyn_ldf_bilap : iso-level bilaplacian operator'
104         WRITE(numout,*) '~~~~~~~~~~~~~'
105      ENDIF
106
[474]107!!bug gm this should be enough
108!!$      zuf(:,:,jpk) = 0.e0
109!!$      zut(:,:,jpk) = 0.e0
110!!$      zlu(:,:,jpk) = 0.e0
111!!$      zlv(:,:,jpk) = 0.e0
[2715]112      zuf(:,:,:) = 0._wp
113      zut(:,:,:) = 0._wp
114      zlu(:,:,:) = 0._wp
115      zlv(:,:,:) = 0._wp
[474]116
[3]117      !                                                ! ===============
118      DO jk = 1, jpkm1                                 ! Horizontal slab
119         !                                             ! ===============
120         ! Laplacian
121         ! ---------
122
[455]123         IF( ln_sco .OR. ln_zps ) THEN   ! s-coordinate or z-coordinate with partial steps
[474]124            zuf(:,:,jk) = rotb(:,:,jk) * fse3f(:,:,jk)
[3]125            DO jj = 2, jpjm1
126               DO ji = fs_2, fs_jpim1   ! vector opt.
[474]127                  zlu(ji,jj,jk) = - ( zuf(ji,jj,jk) - zuf(ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   &
[3]128                     &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj,jk) ) / e1u(ji,jj)
129   
[474]130                  zlv(ji,jj,jk) = + ( zuf(ji,jj,jk) - zuf(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   &
[3]131                     &         + ( hdivb(ji,jj+1,jk) - hdivb(ji,jj,jk) ) / e2v(ji,jj)
132               END DO
133            END DO
[455]134         ELSE                            ! z-coordinate - full step
[3]135            DO jj = 2, jpjm1
136               DO ji = fs_2, fs_jpim1   ! vector opt.
[474]137                  zlu(ji,jj,jk) = - ( rotb (ji  ,jj,jk) - rotb (ji,jj-1,jk) ) / e2u(ji,jj)   &
[3]138                     &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj  ,jk) ) / e1u(ji,jj)
139   
[474]140                  zlv(ji,jj,jk) = + ( rotb (ji,jj  ,jk) - rotb (ji-1,jj,jk) ) / e1v(ji,jj)   &
[3]141                     &         + ( hdivb(ji,jj+1,jk) - hdivb(ji  ,jj,jk) ) / e2v(ji,jj)
142               END DO 
143            END DO 
144         ENDIF
[2715]145      END DO
146      CALL lbc_lnk( zlu, 'U', -1. )   ;   CALL lbc_lnk( zlv, 'V', -1. )   ! Boundary conditions
[3]147
148         
[474]149      DO jk = 1, jpkm1
150   
[3]151         ! Third derivative
152         ! ----------------
153         
154         ! Multiply by the eddy viscosity coef. (at u- and v-points)
[474]155         zlu(:,:,jk) = zlu(:,:,jk) * fsahmu(:,:,jk)
156         zlv(:,:,jk) = zlv(:,:,jk) * fsahmv(:,:,jk)
[3]157         
158         ! Contravariant "laplacian"
[474]159         zcu(:,:) = e1u(:,:) * zlu(:,:,jk)
160         zcv(:,:) = e2v(:,:) * zlv(:,:,jk)
[3]161         
162         ! Laplacian curl ( * e3f if s-coordinates or z-coordinate with partial steps)
163         DO jj = 1, jpjm1
164            DO ji = 1, fs_jpim1   ! vector opt.
[474]165               zuf(ji,jj,jk) = fmask(ji,jj,jk) * (  zcv(ji+1,jj  ) - zcv(ji,jj)      &
[3]166                  &                            - zcu(ji  ,jj+1) + zcu(ji,jj)  )   &
167                  &       * fse3f(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) )
168            END DO 
169         END DO 
170
171         ! Laplacian Horizontal fluxes
172         DO jj = 1, jpjm1
173            DO ji = 1, fs_jpim1   ! vector opt.
[474]174               zlu(ji,jj,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * zlu(ji,jj,jk)
175               zlv(ji,jj,jk) = e1v(ji,jj) * fse3v(ji,jj,jk) * zlv(ji,jj,jk)
[3]176            END DO
177         END DO
178
179         ! Laplacian divergence
180         DO jj = 2, jpj
181            DO ji = fs_2, jpi   ! vector opt.
182               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
[474]183               zut(ji,jj,jk) = (  zlu(ji,jj,jk) - zlu(ji-1,jj  ,jk)   &
184                  &             + zlv(ji,jj,jk) - zlv(ji  ,jj-1,jk) ) / zbt
[3]185            END DO
186         END DO
[474]187      END DO
[3]188
189
[235]190      ! boundary conditions on the laplacian curl and div (zuf,zut)
[474]191!!bug gm no need to do this 2 following lbc...
[235]192      CALL lbc_lnk( zuf, 'F', 1. )
193      CALL lbc_lnk( zut, 'T', 1. )
[3]194
[474]195      DO jk = 1, jpkm1     
196   
[3]197         ! Bilaplacian
198         ! -----------
199
200         DO jj = 2, jpjm1
201            DO ji = fs_2, fs_jpim1   ! vector opt.
202               ze2u = e2u(ji,jj) * fse3u(ji,jj,jk)
203               ze2v = e1v(ji,jj) * fse3v(ji,jj,jk)
204               ! horizontal biharmonic diffusive trends
[474]205               zua = - ( zuf(ji  ,jj,jk) - zuf(ji,jj-1,jk) ) / ze2u   &
206                  &  + ( zut(ji+1,jj,jk) - zut(ji,jj  ,jk) ) / e1u(ji,jj)
[3]207
[474]208               zva = + ( zuf(ji,jj  ,jk) - zuf(ji-1,jj,jk) ) / ze2v   &
209                  &  + ( zut(ji,jj+1,jk) - zut(ji  ,jj,jk) ) / e2v(ji,jj)
[3]210               ! add it to the general momentum trends
211               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
212               va(ji,jj,jk) = va(ji,jj,jk) + zva
213            END DO
214         END DO
215
216         !                                             ! ===============
217      END DO                                           !   End of slab
218      !                                                ! ===============
[2715]219      IF( wrk_not_released(2, 1,2)     .OR.   &
220          wrk_not_released(3, 3,4,5,6) )   CALL ctl_stop('dyn_ldf_bilap: failed to release workspace arrays')
221      !
[3]222   END SUBROUTINE dyn_ldf_bilap
223
224   !!======================================================================
225END MODULE dynldf_bilap
Note: See TracBrowser for help on using the repository browser.