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
Line 
1MODULE dynldf_bilap
2   !!======================================================================
3   !!                     ***  MODULE  dynldf_bilap  ***
4   !! Ocean dynamics:  lateral viscosity trend
5   !!======================================================================
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   !!----------------------------------------------------------------------
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
22   USE trdmod          ! ocean dynamics trends
23   USE trdmod_oce      ! ocean variables trends
24   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC   dyn_ldf_bilap   ! called by step.F90
30
31   !! * Control permutation of array indices
32#  include "oce_ftrans.h90"
33#  include "dom_oce_ftrans.h90"
34#  include "ldfdyn_oce_ftrans.h90"
35
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "ldfdyn_substitute.h90"
39#  include "vectopt_loop_substitute.h90"
40   !!----------------------------------------------------------------------
41   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
42   !! $Id$
43   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
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 ]
71      !!      If ln_sco=F and ln_zps=F, the vertical scale factors in the
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      !!----------------------------------------------------------------------
81      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
82      USE wrk_nemo, ONLY:   zcu => wrk_2d_1 , zcv => wrk_2d_2   ! 2D workspace
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
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
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
95      !!----------------------------------------------------------------------
96
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
99      ENDIF
100
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
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
112      zuf(:,:,:) = 0._wp
113      zut(:,:,:) = 0._wp
114      zlu(:,:,:) = 0._wp
115      zlv(:,:,:) = 0._wp
116
117      !                                                ! ===============
118      DO jk = 1, jpkm1                                 ! Horizontal slab
119         !                                             ! ===============
120         ! Laplacian
121         ! ---------
122
123         IF( ln_sco .OR. ln_zps ) THEN   ! s-coordinate or z-coordinate with partial steps
124            zuf(:,:,jk) = rotb(:,:,jk) * fse3f(:,:,jk)
125            DO jj = 2, jpjm1
126               DO ji = fs_2, fs_jpim1   ! vector opt.
127                  zlu(ji,jj,jk) = - ( zuf(ji,jj,jk) - zuf(ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   &
128                     &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj,jk) ) / e1u(ji,jj)
129   
130                  zlv(ji,jj,jk) = + ( zuf(ji,jj,jk) - zuf(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   &
131                     &         + ( hdivb(ji,jj+1,jk) - hdivb(ji,jj,jk) ) / e2v(ji,jj)
132               END DO
133            END DO
134         ELSE                            ! z-coordinate - full step
135            DO jj = 2, jpjm1
136               DO ji = fs_2, fs_jpim1   ! vector opt.
137                  zlu(ji,jj,jk) = - ( rotb (ji  ,jj,jk) - rotb (ji,jj-1,jk) ) / e2u(ji,jj)   &
138                     &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj  ,jk) ) / e1u(ji,jj)
139   
140                  zlv(ji,jj,jk) = + ( rotb (ji,jj  ,jk) - rotb (ji-1,jj,jk) ) / e1v(ji,jj)   &
141                     &         + ( hdivb(ji,jj+1,jk) - hdivb(ji  ,jj,jk) ) / e2v(ji,jj)
142               END DO 
143            END DO 
144         ENDIF
145      END DO
146      CALL lbc_lnk( zlu, 'U', -1. )   ;   CALL lbc_lnk( zlv, 'V', -1. )   ! Boundary conditions
147
148         
149      DO jk = 1, jpkm1
150   
151         ! Third derivative
152         ! ----------------
153         
154         ! Multiply by the eddy viscosity coef. (at u- and v-points)
155         zlu(:,:,jk) = zlu(:,:,jk) * fsahmu(:,:,jk)
156         zlv(:,:,jk) = zlv(:,:,jk) * fsahmv(:,:,jk)
157         
158         ! Contravariant "laplacian"
159         zcu(:,:) = e1u(:,:) * zlu(:,:,jk)
160         zcv(:,:) = e2v(:,:) * zlv(:,:,jk)
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.
165               zuf(ji,jj,jk) = fmask(ji,jj,jk) * (  zcv(ji+1,jj  ) - zcv(ji,jj)      &
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.
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)
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)
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
185            END DO
186         END DO
187      END DO
188
189
190      ! boundary conditions on the laplacian curl and div (zuf,zut)
191!!bug gm no need to do this 2 following lbc...
192      CALL lbc_lnk( zuf, 'F', 1. )
193      CALL lbc_lnk( zut, 'T', 1. )
194
195      DO jk = 1, jpkm1     
196   
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
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)
207
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)
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      !                                                ! ===============
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      !
222   END SUBROUTINE dyn_ldf_bilap
223
224   !!======================================================================
225END MODULE dynldf_bilap
Note: See TracBrowser for help on using the repository browser.