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 @ 31

Last change on this file since 31 was 3, checked in by opalod, 20 years ago

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