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