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/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilap.F90 @ 2633

Last change on this file since 2633 was 2633, checked in by trackstand2, 13 years ago

Renamed wrk_use => wrk_in_use and wrk_release => wrk_not_released
  • Property svn:keywords set to Id
File size: 9.5 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   !! * Routine accessibility
30   PUBLIC dyn_ldf_bilap  ! called by step.F90
31
32   !! * Substitutions
33#  include "domzgr_substitute.h90"
34#  include "ldfdyn_substitute.h90"
35#  include "vectopt_loop_substitute.h90"
36   !!----------------------------------------------------------------------
37   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
38   !! $Id$
39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41CONTAINS
42
43   SUBROUTINE dyn_ldf_bilap( kt )
44      !!----------------------------------------------------------------------
45      !!                     ***  ROUTINE dyn_ldf_bilap  ***
46      !!
47      !! ** Purpose :   Compute the before trend of the lateral momentum
48      !!      diffusion and add it to the general trend of momentum equation.
49      !!
50      !! ** Method  :   The before horizontal momentum diffusion trend is a
51      !!      bi-harmonic operator (bilaplacian type) which separates the
52      !!      divergent and rotational parts of the flow.
53      !!      Its horizontal components are computed as follow:
54      !!      laplacian:
55      !!          zlu = 1/e1u di[ hdivb ] - 1/(e2u*e3u) dj-1[ e3f rotb ]
56      !!          zlv = 1/e2v dj[ hdivb ] + 1/(e1v*e3v) di-1[ e3f rotb ]
57      !!      third derivative:
58      !!       * multiply by the eddy viscosity coef. at u-, v-point, resp.
59      !!          zlu = ahmu * zlu
60      !!          zlv = ahmv * zlv
61      !!       * curl and divergence of the laplacian
62      !!          zuf = 1/(e1f*e2f) ( di[e2v zlv] - dj[e1u zlu] )
63      !!          zut = 1/(e1t*e2t*e3t) ( di[e2u*e3u zlu] + dj[e1v*e3v zlv] )
64      !!      bilaplacian:
65      !!              diffu = 1/e1u di[ zut ] - 1/(e2u*e3u) dj-1[ e3f zuf ]
66      !!              diffv = 1/e2v dj[ zut ] + 1/(e1v*e3v) di-1[ e3f zuf ]
67      !!      If ln_sco=F and ln_zps=F, the vertical scale factors in the
68      !!      rotational part of the diffusion are simplified
69      !!      Add this before trend to the general trend (ua,va):
70      !!            (ua,va) = (ua,va) + (diffu,diffv)
71      !!      'key_trddyn' defined: the two components of the horizontal
72      !!                               diffusion trend are saved.
73      !!
74      !! ** Action : - Update (ua,va) with the before iso-level biharmonic
75      !!               mixing trend.
76      !!----------------------------------------------------------------------
77      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
78      USE wrk_nemo, ONLY: zcu => wrk_2d_1, zcv => wrk_2d_2   ! 3D workspace
79      USE wrk_nemo, ONLY: zuf => wrk_3d_1, zut => wrk_3d_2   ! 3D workspace
80      USE wrk_nemo, ONLY: zlu => wrk_3d_3, zlv => wrk_3d_4
81      !
82      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
83      !
84      INTEGER  ::   ji, jj, jk                ! dummy loop indices
85      REAL(wp) ::   zua, zva, zbt, ze2u, ze2v ! temporary scalar
86      !!----------------------------------------------------------------------
87
88      IF( wrk_in_use(2, 1,2) .OR. wrk_in_use(3, 1,2,3,4) ) THEN
89         CALL ctl_stop('dyn_ldf_bilap : requested workspace arrays unavailable')   ;   RETURN
90      ENDIF
91
92      IF( kt == nit000 .AND. lwp ) THEN
93         WRITE(numout,*)
94         WRITE(numout,*) 'dyn_ldf_bilap : iso-level bilaplacian operator'
95         WRITE(numout,*) '~~~~~~~~~~~~~'
96      ENDIF
97
98!!bug gm this should be enough
99!!$      zuf(:,:,jpk) = 0.e0
100!!$      zut(:,:,jpk) = 0.e0
101!!$      zlu(:,:,jpk) = 0.e0
102!!$      zlv(:,:,jpk) = 0.e0
103      zuf(:,:,:) = 0._wp
104      zut(:,:,:) = 0._wp
105      zlu(:,:,:) = 0._wp
106      zlv(:,:,:) = 0._wp
107
108      !                                                ! ===============
109      DO jk = 1, jpkm1                                 ! Horizontal slab
110         !                                             ! ===============
111         ! Laplacian
112         ! ---------
113
114         IF( ln_sco .OR. ln_zps ) THEN   ! s-coordinate or z-coordinate with partial steps
115            zuf(:,:,jk) = rotb(:,:,jk) * fse3f(:,:,jk)
116            DO jj = 2, jpjm1
117               DO ji = fs_2, fs_jpim1   ! vector opt.
118                  zlu(ji,jj,jk) = - ( zuf(ji,jj,jk) - zuf(ji,jj-1,jk) ) / ( e2u(ji,jj) * fse3u(ji,jj,jk) )   &
119                     &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj,jk) ) / e1u(ji,jj)
120   
121                  zlv(ji,jj,jk) = + ( zuf(ji,jj,jk) - zuf(ji-1,jj,jk) ) / ( e1v(ji,jj) * fse3v(ji,jj,jk) )   &
122                     &         + ( hdivb(ji,jj+1,jk) - hdivb(ji,jj,jk) ) / e2v(ji,jj)
123               END DO
124            END DO
125         ELSE                            ! z-coordinate - full step
126            DO jj = 2, jpjm1
127               DO ji = fs_2, fs_jpim1   ! vector opt.
128                  zlu(ji,jj,jk) = - ( rotb (ji  ,jj,jk) - rotb (ji,jj-1,jk) ) / e2u(ji,jj)   &
129                     &         + ( hdivb(ji+1,jj,jk) - hdivb(ji,jj  ,jk) ) / e1u(ji,jj)
130   
131                  zlv(ji,jj,jk) = + ( rotb (ji,jj  ,jk) - rotb (ji-1,jj,jk) ) / e1v(ji,jj)   &
132                     &         + ( hdivb(ji,jj+1,jk) - hdivb(ji  ,jj,jk) ) / e2v(ji,jj)
133               END DO 
134            END DO 
135         ENDIF
136      END DO
137      CALL lbc_lnk( zlu, 'U', -1. )   ;   CALL lbc_lnk( zlv, 'V', -1. )   ! Boundary conditions
138
139         
140      DO jk = 1, jpkm1
141   
142         ! Third derivative
143         ! ----------------
144         
145         ! Multiply by the eddy viscosity coef. (at u- and v-points)
146         zlu(:,:,jk) = zlu(:,:,jk) * fsahmu(:,:,jk)
147         zlv(:,:,jk) = zlv(:,:,jk) * fsahmv(:,:,jk)
148         
149         ! Contravariant "laplacian"
150         zcu(:,:) = e1u(:,:) * zlu(:,:,jk)
151         zcv(:,:) = e2v(:,:) * zlv(:,:,jk)
152         
153         ! Laplacian curl ( * e3f if s-coordinates or z-coordinate with partial steps)
154         DO jj = 1, jpjm1
155            DO ji = 1, fs_jpim1   ! vector opt.
156               zuf(ji,jj,jk) = fmask(ji,jj,jk) * (  zcv(ji+1,jj  ) - zcv(ji,jj)      &
157                  &                            - zcu(ji  ,jj+1) + zcu(ji,jj)  )   &
158                  &       * fse3f(ji,jj,jk) / ( e1f(ji,jj)*e2f(ji,jj) )
159            END DO 
160         END DO 
161
162         ! Laplacian Horizontal fluxes
163         DO jj = 1, jpjm1
164            DO ji = 1, fs_jpim1   ! vector opt.
165               zlu(ji,jj,jk) = e2u(ji,jj) * fse3u(ji,jj,jk) * zlu(ji,jj,jk)
166               zlv(ji,jj,jk) = e1v(ji,jj) * fse3v(ji,jj,jk) * zlv(ji,jj,jk)
167            END DO
168         END DO
169
170         ! Laplacian divergence
171         DO jj = 2, jpj
172            DO ji = fs_2, jpi   ! vector opt.
173               zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
174               zut(ji,jj,jk) = (  zlu(ji,jj,jk) - zlu(ji-1,jj  ,jk)   &
175                  &             + zlv(ji,jj,jk) - zlv(ji  ,jj-1,jk) ) / zbt
176            END DO
177         END DO
178      END DO
179
180
181      ! boundary conditions on the laplacian curl and div (zuf,zut)
182!!bug gm no need to do this 2 following lbc...
183      CALL lbc_lnk( zuf, 'F', 1. )
184      CALL lbc_lnk( zut, 'T', 1. )
185
186      DO jk = 1, jpkm1     
187   
188         ! Bilaplacian
189         ! -----------
190
191         DO jj = 2, jpjm1
192            DO ji = fs_2, fs_jpim1   ! vector opt.
193               ze2u = e2u(ji,jj) * fse3u(ji,jj,jk)
194               ze2v = e1v(ji,jj) * fse3v(ji,jj,jk)
195               ! horizontal biharmonic diffusive trends
196               zua = - ( zuf(ji  ,jj,jk) - zuf(ji,jj-1,jk) ) / ze2u   &
197                  &  + ( zut(ji+1,jj,jk) - zut(ji,jj  ,jk) ) / e1u(ji,jj)
198
199               zva = + ( zuf(ji,jj  ,jk) - zuf(ji-1,jj,jk) ) / ze2v   &
200                  &  + ( zut(ji,jj+1,jk) - zut(ji  ,jj,jk) ) / e2v(ji,jj)
201               ! add it to the general momentum trends
202               ua(ji,jj,jk) = ua(ji,jj,jk) + zua
203               va(ji,jj,jk) = va(ji,jj,jk) + zva
204            END DO
205         END DO
206
207         !                                             ! ===============
208      END DO                                           !   End of slab
209      !                                                ! ===============
210      IF( wrk_not_released(2, 1,2)      .OR.   &
211          wrk_not_released(3, 1,2,3,4) )   CALL ctl_stop('dyn_ldf_bilap : failed to release workspace arrays')
212      !
213   END SUBROUTINE dyn_ldf_bilap
214
215   !!======================================================================
216END MODULE dynldf_bilap
Note: See TracBrowser for help on using the repository browser.