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