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

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

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