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

source: trunk/NEMO/OPA_SRC/DYN/dynzdf_iso.F90 @ 3

Last change on this file since 3 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: 16.7 KB
Line 
1MODULE dynzdf_iso
2   !!==============================================================================
3   !!                       ***  MODULE  dynzdf_iso  ***
4   !! Ocean dynamics:  vertical component(s) of the momentum mixing trend
5   !!==============================================================================
6#if   defined key_ldfslp   ||   defined key_esopa
7   !!----------------------------------------------------------------------
8   !!   'key_ldfslp'                          rotation of the mixing tensor
9   !!----------------------------------------------------------------------
10   !!   dyn_zdf_iso  : update the momentum trend with the vertical diffusion
11   !!                  (vertical mixing + vertical component of lateral
12   !!                  mixing) (rotated lateral operator case)
13   !!----------------------------------------------------------------------
14   !! * Modules used
15   USE oce             ! ocean dynamics and tracers
16   USE dom_oce         ! ocean space and time domain
17   USE phycst          ! physical constants
18   USE zdf_oce         ! ocean vertical physics
19   USE in_out_manager  ! I/O manager
20   USE taumod          ! surface ocean stress
21   USE trddyn_oce     ! dynamics trends diagnostics variables
22
23   IMPLICIT NONE
24   PRIVATE
25
26   !! * Routine accessibility
27   PUBLIC dyn_zdf_iso    ! called by step.F90
28
29   !! * Substitutions
30#  include "domzgr_substitute.h90"
31#  include "vectopt_loop_substitute.h90"
32   !!----------------------------------------------------------------------
33   !!   OPA 9.0 , LODYC-IPSL  (2003)
34   !!----------------------------------------------------------------------
35
36CONTAINS
37
38   SUBROUTINE dyn_zdf_iso( kt )
39      !!----------------------------------------------------------------------
40      !!                  ***  ROUTINE dyn_zdf_iso  ***
41      !!                   
42      !! ** Purpose :
43      !!         Compute the vertical momentum trend due to both vertical and
44      !!      lateral mixing (only for second order lateral operator, for
45      !!      fourth order it is already computed and add to the general trend
46      !!      in dynldf.F) and the surface forcing, and add it to the general
47      !!      trend of the momentum equations.
48      !!
49      !! ** Method :
50      !!         The vertical component of the lateral diffusive trends is
51      !!      provided by a 2nd order operator rotated along neural or geopo-
52      !!      tential surfaces to which an eddy induced advection can be added
53      !!      It is computed using before fields (forward in time) and isopyc-
54      !!      nal or geopotential slopes computed in routine ldfslp.
55      !!
56      !!      First part: vertical trends associated with the lateral mixing
57      !!      ==========  (excluding the vertical flux proportional to dk[U] )
58      !!      vertical fluxes associated with the rotated lateral mixing:
59      !!         zfuw =-ahm {  e2t*mi(wslpi) di[ mi(mk(ub)) ]
60      !!                     + e1t*mj(wslpj) dj[ mj(mk(ub)) ]  }
61      !!      update and save in zavt the vertical eddy viscosity coefficient:
62      !!         avmu = avmu + mi(wslpi)^2 + mj(wslj)^2
63      !!      take the horizontal divergence of the fluxes:
64      !!         diffu = 1/(e1u*e2u*e3u) dk[ zfuw ]
65      !!      Add this trend to the general trend (ta,sa):
66      !!         ua = ua + difft
67      !!
68      !!      Second part: vertical trend associated with the vertical physics
69      !!      ===========  (including the vertical flux proportional to dk[U]
70      !!                    associated with the lateral mixing, through the
71      !!                    update of avmu)
72      !!      The vertical diffusion of momentum is given by:
73      !!             diffu = dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) )
74      !!      using a backward (implicit) time stepping.
75      !!      Bottom boundary conditions : bottom stress (cf zdfbfr.F)
76      !!      Add this trend to the general trend ua :
77      !!         ua = ua + dz( avmu dz(u) )
78      !!
79      !!      'key_trddyn' defined: trend saved for futher diagnostics.
80      !!
81      !!      macro-tasked on vertical slab (jj-loop)
82      !!
83      !! ** Action :
84      !!            /comaft/ ua, va   : general momentum trend increased
85      !!                                by the after vertical diffusive trend
86      !!            /comtra/ utrd,vtrd: after vertical momentum diffusive
87      !!                                trend ('key_trddyn' defined)
88      !!
89      !! History :
90      !!      original : 90-10 (B. Blanke)
91      !!      addition : 97-05 (G. Madec) vertical component of isopycnal
92      !!---------------------------------------------------------------------
93      !! * Modules used
94      USE ldfslp    , ONLY : wslpi, wslpj
95      USE ldftra_oce, ONLY : aht0
96
97      !! * Arguments
98      INTEGER, INTENT( in ) ::   kt       ! ocean time-step index
99     
100      !! * Local declarations
101      INTEGER ::   ji, jj, jk             ! dummy loop indices
102      INTEGER ::   &
103         ikst, ikenm2, ikstp1             ! temporary integers
104      REAL(wp) ::   &
105         zrau0r, z2dt,                 &  ! temporary scalars
106         z2dtf, zua, zva, zcoef, zzws
107      REAL(wp) ::   &
108         zcoef0, zcoef3, zcoef4, zbu, zbv, zmkt, zmkf,   &
109         zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj
110      REAL(wp), DIMENSION(jpi,jpk) ::   &
111         zwx, zwy, zwz                 &  ! workspace
112       , zwd, zws, zwi, zwt
113      REAL(wp), DIMENSION(jpi,jpk) ::   &
114         zfuw, zdiu, zdju, zdj1u,      &  ! workspace
115         zfvw, zdiv, zdjv, zdj1v
116#if defined key_trddyn
117      INTEGER ::   &
118         ikbu, ikbum1 , ikbv, ikbvm1      ! temporary integers
119#endif
120      !!----------------------------------------------------------------------
121
122      IF( kt == nit000 ) THEN
123         IF(lwp) WRITE(numout,*)
124         IF(lwp) WRITE(numout,*) 'dyn_zdf_iso : vertical momentum diffusion isopycnal operator'
125         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
126      ENDIF
127
128      ! 0. Local constant initialization
129      ! --------------------------------
130     
131      ! inverse of the reference density
132      zrau0r = 1. / rau0
133      ! Leap-frog environnement
134      z2dt = 2. * rdt
135      ! Euler time stepping when starting from rest
136      IF ( neuler == 0 .AND. kt == nit000 ) z2dt = rdt
137
138      !                                                ! ===============
139      DO jj = 2, jpjm1                                 !  Vertical slab
140         !                                             ! ===============
141
142 
143         ! I. vertical trends associated with the lateral mixing
144         ! =====================================================
145         !  (excluding the vertical flux proportional to dk[t]
146
147
148         ! I.1 horizontal momentum gradient
149         ! --------------------------------
150
151         DO jk = 1, jpk
152            DO ji = 2, jpi
153               ! i-gradient of u at jj
154               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( ub(ji,jj  ,jk) - ub(ji-1,jj  ,jk) )
155               ! j-gradient of u and v at jj
156               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( ub(ji,jj+1,jk) - ub(ji  ,jj  ,jk) )
157               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( vb(ji,jj  ,jk) - vb(ji  ,jj-1,jk) )
158               ! j-gradient of u and v at jj+1
159               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj  ,jk) - ub(ji  ,jj-1,jk) )
160               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji  ,jj  ,jk) )
161            END DO
162         END DO
163         DO jk = 1, jpk
164            DO ji = 1, jpim1
165               ! i-gradient of v at jj
166               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( vb(ji+1,jj,jk) - vb(ji  ,jj  ,jk) )
167            END DO
168         END DO
169
170
171         ! I.2 Vertical fluxes
172         ! -------------------
173
174         ! Surface and bottom vertical fluxes set to zero
175         DO ji = 1, jpi
176            zfuw(ji, 1 ) = 0.e0
177            zfvw(ji, 1 ) = 0.e0
178            zfuw(ji,jpk) = 0.e0
179            zfvw(ji,jpk) = 0.e0
180         END DO
181
182         ! interior (2=<jk=<jpk-1) on U field
183         DO jk = 2, jpkm1
184            DO ji = 2, jpim1
185               zcoef0= 0.5 * aht0 * umask(ji,jj,jk)
186
187               zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
188               zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
189
190               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   &
191                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. )
192               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   &
193                             + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. )
194
195               zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
196               zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
197               ! vertical flux on u field
198               zfuw(ji,jk) = zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)     &
199                                       +zdiu (ji,jk  ) + zdiu (ji+1,jk  ) )   &
200                           + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     &
201                                       +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) )
202               ! update avmu (add isopycnal vertical coefficient to avmu)
203               avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi          &
204                                                 + zuwslpj * zuwslpj ) / aht0
205            END DO
206         END DO
207
208         ! interior (2=<jk=<jpk-1) on V field
209         DO jk = 2, jpkm1
210            DO ji = 2, jpim1
211               zcoef0= 0.5 * aht0 * vmask(ji,jj,jk)
212
213               zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
214               zvwslpj = zcoef0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
215
216               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)   &
217                             + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ), 1. )
218               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)   &
219                             + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ), 1. )
220
221               zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
222               zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
223               ! vertical flux on v field
224               zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     &
225                                       +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
226                           + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
227                                       +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )
228               ! update avmv (add isopycnal vertical coefficient to avmv)
229               avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi          &
230                                                 + zvwslpj * zvwslpj ) / aht0
231            END DO
232         END DO
233
234
235         ! I.3 Divergence of vertical fluxes added to the general tracer trend
236         ! -------------------------------------------------------------------
237
238         DO jk = 1, jpkm1
239            DO ji = 2, jpim1
240               ! volume elements
241               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
242               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
243               ! part of the k-component of isopycnal momentum diffusive trends
244               zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu
245               zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv
246               ! add the trends to the general trends
247               ua(ji,jj,jk) = ua(ji,jj,jk) + zuav
248               va(ji,jj,jk) = va(ji,jj,jk) + zvav
249#if defined key_trddyn
250               ! save the trends for diagnostics
251               utrd(ji,jj,jk,5) = utrd(ji,jj,jk,5) + zuav
252               vtrd(ji,jj,jk,5) = vtrd(ji,jj,jk,5) + zvav
253#endif
254            END DO
255         END DO
256
257
258         ! 1. Vertical diffusion on u
259         ! ---------------------------
260
261         ! 1.0 Matrix and second member construction
262         ! bottom boundary condition: only zws must be masked as avmu can take
263         ! non zero value at the ocean bottom depending on the bottom friction
264         ! used (see zdfmix.F)
265         DO jk = 1, jpkm1
266            DO ji = 2, jpim1
267               zcoef = - z2dt / fse3u(ji,jj,jk)
268               zwi(ji,jk) = zcoef * avmu(ji,jj,jk  ) / fse3uw(ji,jj,jk  )
269               zzws       = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1)
270               zws(ji,jk) = zzws * umask(ji,jj,jk+1)
271               zwd(ji,jk) = 1. - zwi(ji,jk) - zzws
272               zwy(ji,jk) = ub(ji,jj,jk) + z2dt * ua(ji,jj,jk)
273            END DO 
274         END DO 
275
276         ! 1.1 Surface boudary conditions
277         DO ji = 2, jpim1
278            z2dtf = z2dt / ( fse3u(ji,jj,1)*rau0 )
279            zwi(ji,1) = 0.
280            zwd(ji,1) = 1. - zws(ji,1)
281            zwy(ji,1) = zwy(ji,1) + z2dtf * taux(ji,jj)
282         END DO 
283
284         ! 1.2 Matrix inversion starting from the first level
285         ikst = 1
286#include "zdf.matrixsolver.h90"
287
288         ! 1.3 Normalization to obtain the general momentum trend ua
289         DO jk = 1, jpkm1
290            DO ji = 2, jpim1
291               zua = ( zwx(ji,jk) - ub(ji,jj,jk) ) / z2dt
292#if defined key_trddyn
293               ! save the vertical diffusive momentum trend
294               utrd(ji,jj,jk,7) = zua - ua(ji,jj,jk)
295#endif
296               ua(ji,jj,jk) = zua
297            END DO 
298         END DO 
299
300#if defined key_trddyn
301
302         ! 1.4 diagnose surface and bottom momentum fluxes
303         DO ji = 2, jpim1
304            ! save the surface forcing momentum fluxes
305            tautrd(ji,jj,1) = taux(ji,jj) / ( fse3u(ji,jj,1)*rau0 )
306            ! save bottom friction momentum fluxes
307            ikbu   = min( mbathy(ji+1,jj), mbathy(ji,jj) )
308            ikbum1 = max( ikbu-1, 1 )
309            tautrd(ji,jj,3) = - avmu(ji,jj,ikbu) * zwx(ji,ikbum1)   &
310                            / ( fse3u(ji,jj,ikbum1)*fse3uw(ji,jj,ikbu) )
311            ! subtract surface forcing and bottom friction trend from vertical
312            ! diffusive momentum trend
313            utrd(ji,jj,1     ,7) = utrd(ji,jj,1     ,7) - tautrd(ji,jj,1)
314            utrd(ji,jj,ikbum1,7) = utrd(ji,jj,ikbum1,7) - tautrd(ji,jj,3)
315         END DO
316
317#endif
318
319         ! 2. Vertical diffusion on v
320         ! ---------------------------
321
322         ! 2.0 Matrix and second member construction
323         ! bottom boundary condition: only zws must be masked as avmv can take
324         ! non zero value at the ocean bottom depending on the bottom friction
325         ! used (see zdfmix.F)
326         DO jk = 1, jpkm1
327            DO ji = 2, jpim1
328               zcoef = -z2dt/fse3v(ji,jj,jk)
329               zwi(ji,jk) = zcoef * avmv(ji,jj,jk  ) / fse3vw(ji,jj,jk  )
330               zzws       = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)
331               zws(ji,jk) =  zzws * vmask(ji,jj,jk+1)
332               zwd(ji,jk) = 1. - zwi(ji,jk) - zzws
333               zwy(ji,jk) = vb(ji,jj,jk) + z2dt * va(ji,jj,jk)
334            END DO 
335         END DO 
336
337         ! 2.1 Surface boudary conditions
338         DO ji = 2, jpim1
339            z2dtf = z2dt / ( fse3v(ji,jj,1)*rau0 )
340            zwi(ji,1) = 0.e0
341            zwd(ji,1) = 1. - zws(ji,1)
342            zwy(ji,1) = zwy(ji,1) + z2dtf * tauy(ji,jj)
343         END DO 
344
345         ! 2.2 Matrix inversion starting from the first level
346         ikst = 1
347#include "zdf.matrixsolver.h90"
348
349         ! 2.3 Normalization to obtain the general momentum trend va
350         DO jk = 1, jpkm1
351            DO ji = 2, jpim1
352               zva = ( zwx(ji,jk) - vb(ji,jj,jk) ) / z2dt
353#if defined key_trddyn
354               ! save the vertical diffusive momentum fluxes
355               vtrd(ji,jj,jk,7) = zva - va(ji,jj,jk)
356#endif
357               va(ji,jj,jk) = zva
358            END DO   
359         END DO   
360
361#if defined key_trddyn
362         ! 2.4 diagnose surface and bottom momentum fluxes
363         DO ji = 2, jpim1
364            ! save the surface forcing momentum fluxes
365            tautrd(ji,jj,2) = tauy(ji,jj) / ( fse3v(ji,jj,1)*rau0 )
366            ! save bottom friction momentum fluxes
367            ikbv   = min( mbathy(ji,jj+1), mbathy(ji,jj) )
368            ikbvm1 = max( ikbv-1, 1 )
369            tautrd(ji,jj,4) = - avmv(ji,jj,ikbv) * zwx(ji,ikbvm1)   &
370                            / ( fse3v(ji,jj,ikbvm1)*fse3vw(ji,jj,ikbv) )
371            ! subtract surface forcing and bottom friction trend from vertical
372            ! diffusive momentum trend
373            vtrd(ji,jj,1     ,7) = vtrd(ji,jj,1     ,7) - tautrd(ji,jj,2)
374            vtrd(ji,jj,ikbvm1,7) = vtrd(ji,jj,ikbvm1,7) - tautrd(ji,jj,4)
375         END DO
376
377#endif
378         !                                             ! ===============
379      END DO                                           !   End of slab
380      !                                                ! ===============
381   END SUBROUTINE dyn_zdf_iso
382
383#else
384   !!----------------------------------------------------------------------
385   !!   Dummy module                       NO rotation of the mixing tensor
386   !!----------------------------------------------------------------------
387CONTAINS
388   SUBROUTINE dyn_zdf_iso( kt )                        ! Dummy routine
389      WRITE(*,*) kt
390   END SUBROUTINE dyn_zdf_iso
391#endif
392
393   !!==============================================================================
394END MODULE dynzdf_iso
Note: See TracBrowser for help on using the repository browser.