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_iso.F90 in trunk/NEMO/OPA_SRC/DYN – NEMO

source: trunk/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 247

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

CL : Add CVS Header and CeCILL licence information
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.5 KB
Line 
1MODULE dynldf_iso
2   !!======================================================================
3   !!                     ***  MODULE  dynldf_iso  ***
4   !! Ocean dynamics:  lateral viscosity trend
5   !!======================================================================
6#if defined key_ldfslp   ||   defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_ldfslp'                      slopes of the direction of mixing
9   !!----------------------------------------------------------------------
10   !!   dyn_ldf_iso  : update the momentum trend with the horizontal part
11   !!                  of the lateral diffusion using isopycnal or horizon-
12   !!                  tal s-coordinate laplacian operator.
13   !!----------------------------------------------------------------------
14   !! * Modules used
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE ldfdyn_oce      ! ocean dynamics lateral physics
18   USE ldftra_oce      ! ocean tracer   lateral physics
19   USE zdf_oce         ! ocean vertical physics
20   USE trdmod          ! ocean dynamics trends
21   USE trdmod_oce      ! ocean variables trends
22   USE ldfslp          ! iso-neutral slopes
23   USE in_out_manager  ! I/O manager
24
25   IMPLICIT NONE
26   PRIVATE
27
28   !! * Routine accessibility
29   PUBLIC dyn_ldf_iso           ! called by step.F90
30
31   !! * Substitutions
32#  include "domzgr_substitute.h90"
33#  include "ldfdyn_substitute.h90"
34#  include "vectopt_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36   !!   OPA 9.0 , LOCEAN-IPSL (2005)
37   !! $Header$
38   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
39   !!----------------------------------------------------------------------
40
41CONTAINS
42
43   SUBROUTINE dyn_ldf_iso( kt )
44      !!----------------------------------------------------------------------
45      !!                     ***  ROUTINE dyn_ldf_iso  ***
46      !!                       
47      !! ** Purpose :   Compute the before trend of the horizontal part of the
48      !!      lateral momentum diffusion and add it to the general trend of
49      !!      momentum equation.
50      !!
51      !! ** Method :
52      !!        The horizontal component of the lateral diffusive trends on
53      !!      momentum is provided by a 2nd order operator rotated along neu-
54      !!      tral or geopotential surfaces (in s-coordinates).
55      !!      It is computed using before fields (forward in time) and isopyc-
56      !!      nal or geopotential slopes computed in routine ldfslp or inildf.
57      !!      Here, u and v components are considered as 2 independent scalar
58      !!      fields. Therefore, the property of splitting divergent and rota-
59      !!      tional part of the flow of the standard, z-coordinate laplacian
60      !!      momentum diffusion is lost.
61      !!      horizontal fluxes associated with the rotated lateral mixing:
62      !!      u-component:
63      !!         ziut = ( ahmt + ahmb0 ) e2t * e3t / e1t  di[ ub ]
64      !!               -      ahmt       e2t * mi-1(uslp) dk[ mi(mk(ub)) ]
65      !!         zjuf = ( ahmf + ahmb0 ) e1f * e3f / e2f  dj[ ub ]
66      !!               -      ahmf       e1f * mi(vslp)   dk[ mj(mk(ub)) ]
67      !!      v-component:
68      !!         zivf = ( ahmf + ahmb0 ) e2t * e3t / e1t  di[ vb ]
69      !!               -      ahmf       e2t * mj(uslp)   dk[ mi(mk(vb)) ]
70      !!         zjvt = ( ahmt + ahmb0 ) e1f * e3f / e2f  dj[ ub ]
71      !!               -      ahmt       e1f * mj-1(vslp) dk[ mj(mk(vb)) ]
72      !!      take the horizontal divergence of the fluxes:
73      !!         diffu = 1/(e1u*e2u*e3u) {  di  [ ziut ] + dj-1[ zjuf ]  }
74      !!         diffv = 1/(e1v*e2v*e3v) {  di-1[ zivf ] + dj  [ zjvt ]  }
75      !!      Add this trend to the general trend (ua,va):
76      !!         ua = ua + diffu
77      !!      'key_trddyn' defined: the trends are saved for diagnostics.
78      !!
79      !! ** Action :
80      !!        Update (ua,va) arrays with the before geopotential biharmonic
81      !!      mixing trend.
82      !!        Save in (uldftrd,vldftrd) arrays the trends if 'key_trddyn' defined
83      !!
84      !! History :
85      !!   8.0  !  97-07  (G. Madec)  Original code
86      !!   8.5  !  02-08  (G. Madec)  F90: Free form and module
87      !!   9.0  !  04-08  (C. Talandier) New trends organization
88      !!----------------------------------------------------------------------
89      !! * Modules used     
90      USE oce, ONLY :    ztdua => ta,   & ! use ta as 3D workspace   
91                         ztdva => sa      ! use sa as 3D workspace   
92      !! * Arguments
93      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
94
95      !! * Local declarations
96      INTEGER  ::   ji, jj, jk            ! dummy loop indices
97      REAL(wp) ::   &
98         zabe1, zabe2, zcof1, zcof2,   &  ! temporary scalars
99         zmskt, zmskf, zbu, zbv,       &
100         zuah, zvah, zua, zva
101      REAL(wp), DIMENSION(jpi,jpj) ::   &
102         ziut, zjuf, zjvt, zivf,        & ! temporary workspace
103         zdku, zdk1u, zdkv, zdk1v
104      !!----------------------------------------------------------------------
105
106      IF( kt == nit000 ) THEN
107         IF(lwp) WRITE(numout,*)
108         IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
109         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate horizontal diffusive operator'
110      ENDIF
111
112      ! Save ua and va trends
113      IF( l_trddyn )   THEN
114         ztdua(:,:,:) = ua(:,:,:) 
115         ztdva(:,:,:) = va(:,:,:) 
116      ENDIF
117      !                                                ! ===============
118      DO jk = 1, jpkm1                                 ! Horizontal slab
119         !                                             ! ===============
120
121         ! Vertical u- and v-shears at level jk and jk+1
122         ! ---------------------------------------------
123         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
124         !                             zdkv(jk=1)=zdkv(jk=2)
125
126         zdk1u(:,:) = ( ub(:,:,jk) -ub(:,:,jk+1) ) * umask(:,:,jk+1)
127         zdk1v(:,:) = ( vb(:,:,jk) -vb(:,:,jk+1) ) * vmask(:,:,jk+1)
128
129         IF( jk == 1 ) THEN
130            zdku(:,:) = zdk1u(:,:)
131            zdkv(:,:) = zdk1v(:,:)
132         ELSE
133            zdku(:,:) = ( ub(:,:,jk-1) - ub(:,:,jk) ) * umask(:,:,jk)
134            zdkv(:,:) = ( vb(:,:,jk-1) - vb(:,:,jk) ) * vmask(:,:,jk)
135         ENDIF
136
137         !                               -----f-----
138         ! Horizontal fluxes on U             | 
139         ! --------------------===        t   u   t
140         !                                    | 
141         ! i-flux at t-point             -----f-----
142
143         DO jj = 2, jpjm1
144            DO ji = fs_2, jpi   ! vector opt.
145               zabe1 = ( fsahmt(ji,jj,jk) + ahmb0 )   &
146#if defined key_partial_steps
147                     * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1, jj,jk) ) / e1t(ji,jj)
148#else
149                     * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
150#endif
151
152               zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
153                              + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
154
155               zcof1 = - aht0 * e2t(ji,jj) * zmskt   &
156                     * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
157
158               ziut(ji,jj) = tmask(ji,jj,jk) *   &
159                           (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
160                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
161                                       +zdk1u(ji,jj) + zdku (ji-1,jj) )  )
162            END DO
163         END DO
164
165         ! j-flux at f-point
166         DO jj = 1, jpjm1
167            DO ji = 1, fs_jpim1   ! vector opt.
168               zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 )   &
169                     * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
170
171               zmskf = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
172                              + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )
173
174               zcof2 = - aht0 * e1f(ji,jj) * zmskf   &
175                     * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
176
177               zjuf(ji,jj) = fmask(ji,jj,jk) *   &
178                           (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   &
179                            + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
180                                       +zdk1u(ji,jj+1) + zdku (ji,jj) )  )
181            END DO
182         END DO
183
184         !                                |   t   |
185         ! Horizontal fluxes on V         |       |
186         ! --------------------===        f---v---f
187         !                                |       |
188         ! i-flux at f-point              |   t   |
189
190         DO jj = 2, jpjm1
191            DO ji = 1, fs_jpim1   ! vector opt.
192               zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 )   &
193                     * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
194
195               zmskf = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
196                              + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. )
197
198               zcof1 = - aht0 * e2f(ji,jj) * zmskf   &
199                     * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
200
201               zivf(ji,jj) = fmask(ji,jj,jk) *   &
202                           (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )   &
203                            + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj)     &
204                                       +zdk1v(ji,jj) + zdkv (ji+1,jj) )  )
205            END DO
206         END DO
207
208         ! j-flux at t-point
209         DO jj = 2, jpj
210            DO ji = 1, fs_jpim1   ! vector opt.
211               zabe2 = ( fsahmt(ji,jj,jk) + ahmb0 )   &
212#if defined key_partial_steps
213                     * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji, jj-1, jk) ) / e2t(ji,jj)
214#else
215                     * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
216#endif
217
218               zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
219                              + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
220
221               zcof2 = - aht0 * e1t(ji,jj) * zmskt   &
222                     * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
223
224               zjvt(ji,jj) = tmask(ji,jj,jk) *   &
225                           (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
226                            + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
227                                       +zdk1v(ji,jj-1) + zdkv (ji,jj) )  )
228            END DO
229         END DO
230
231
232         ! Second derivative (divergence) and add to the general trend
233         ! -----------------------------------------------------------
234
235         DO jj = 2, jpjm1
236            DO ji = 2, jpim1          !! Question vectop possible??? !!bug
237               ! volume elements
238               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
239               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
240               ! horizontal component of isopycnal momentum diffusive trends
241               zuah =( ziut (ji+1,jj) - ziut (ji,jj  ) +   &
242                       zjuf (ji  ,jj) - zjuf (ji,jj-1)  ) / zbu
243               zvah =( zivf (ji,jj  ) - zivf (ji-1,jj) +   &
244                       zjvt (ji,jj+1) - zjvt (ji,jj  )  ) / zbv
245               ! add the trends to the general trends
246               ua (ji,jj,jk) = ua (ji,jj,jk) + zuah
247               va (ji,jj,jk) = va (ji,jj,jk) + zvah
248            END DO
249         END DO
250         !                                             ! ===============
251      END DO                                           !   End of slab
252      !                                                ! ===============
253
254      ! save the lateral diffusion trends for diagnostic
255      ! momentum trends will be saved in dynzdf_iso.F90
256      IF( l_trddyn )   THEN
257         uldftrd(:,:,:) = ua(:,:,:) - ztdua(:,:,:)
258         vldftrd(:,:,:) = va(:,:,:) - ztdva(:,:,:)
259      ENDIF
260
261      IF(l_ctl) THEN         ! print sum trends (used for debugging)
262         zua = SUM( ua(2:nictl,2:njctl,1:jpkm1) * umask(2:nictl,2:njctl,1:jpkm1) )
263         zva = SUM( va(2:nictl,2:njctl,1:jpkm1) * vmask(2:nictl,2:njctl,1:jpkm1) )
264         WRITE(numout,*) ' ldf  - Ua: ', zua-u_ctl, ' Va: ', zva-v_ctl
265         u_ctl = zua   ;   v_ctl = zva
266      ENDIF
267
268   END SUBROUTINE dyn_ldf_iso
269
270# else
271   !!----------------------------------------------------------------------
272   !!   Dummy module                           NO rotation of mixing tensor
273   !!----------------------------------------------------------------------
274CONTAINS
275   SUBROUTINE dyn_ldf_iso( kt )               ! Empty routine
276      WRITE(*,*) 'dyn_ldf_iso: You should not have seen this print! error?', kt
277   END SUBROUTINE dyn_ldf_iso
278#endif
279
280   !!======================================================================
281END MODULE dynldf_iso
Note: See TracBrowser for help on using the repository browser.