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 branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 4 years ago

The Dr Hook changes from my perl code.
File size: 22.6 KB
Line 
1MODULE dynldf_iso
2   !!======================================================================
3   !!                     ***  MODULE  dynldf_iso  ***
4   !! Ocean dynamics:  lateral viscosity trend
5   !!======================================================================
6   !! History :  OPA  !  97-07  (G. Madec)  Original code
7   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module
8   !!             -   !  2004-08  (C. Talandier) New trends organization
9   !!            2.0  !  2005-11  (G. Madec)  s-coordinate: horizontal diffusion
10   !!----------------------------------------------------------------------
11#if defined key_ldfslp   ||   defined key_esopa
12   !!----------------------------------------------------------------------
13   !!   'key_ldfslp'                      slopes of the direction of mixing
14   !!----------------------------------------------------------------------
15   !!   dyn_ldf_iso  : update the momentum trend with the horizontal part
16   !!                  of the lateral diffusion using isopycnal or horizon-
17   !!                  tal s-coordinate laplacian operator.
18   !!----------------------------------------------------------------------
19   USE oce             ! ocean dynamics and tracers
20   USE dom_oce         ! ocean space and time domain
21   USE ldfdyn_oce      ! ocean dynamics lateral physics
22   USE ldftra_oce      ! ocean tracer   lateral physics
23   USE zdf_oce         ! ocean vertical physics
24   USE ldfslp          ! iso-neutral slopes
25   !
26   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
27   USE in_out_manager  ! I/O manager
28   USE lib_mpp         ! MPP library
29   USE prtctl          ! Print control
30   USE wrk_nemo        ! Memory Allocation
31   USE timing          ! Timing
32
33   USE yomhook, ONLY: lhook, dr_hook
34   USE parkind1, ONLY: jprb, jpim
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC   dyn_ldf_iso           ! called by step.F90
40   PUBLIC   dyn_ldf_iso_alloc     ! called by nemogcm.F90
41
42   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u   ! 2D workspace (dyn_ldf_iso)
43   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v   !  -      -
44
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47#  include "ldfdyn_substitute.h90"
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/OPA 3.3 , NEMO Consortium (2011)
51   !! $Id$
52   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   INTEGER FUNCTION dyn_ldf_iso_alloc()
57   INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
58   INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
59   REAL(KIND=jprb)               :: zhook_handle
60
61   CHARACTER(LEN=*), PARAMETER :: RoutineName='DYN_LDF_ISO_ALLOC'
62
63   IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
64
65      !!----------------------------------------------------------------------
66      !!                  ***  ROUTINE dyn_ldf_iso_alloc  ***
67      !!----------------------------------------------------------------------
68      ALLOCATE( zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) ,     & 
69         &      zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc )
70         !
71      IF( dyn_ldf_iso_alloc /= 0 )   CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.')
72   IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
73   END FUNCTION dyn_ldf_iso_alloc
74
75
76   SUBROUTINE dyn_ldf_iso( kt )
77      !!----------------------------------------------------------------------
78      !!                     ***  ROUTINE dyn_ldf_iso  ***
79      !!                       
80      !! ** Purpose :   Compute the before trend of the rotated laplacian
81      !!      operator of lateral momentum diffusion except the diagonal
82      !!      vertical term that will be computed in dynzdf module. Add it
83      !!      to the general trend of momentum equation.
84      !!
85      !! ** Method :
86      !!        The momentum lateral diffusive trend is provided by a 2nd
87      !!      order operator rotated along neutral or geopotential surfaces
88      !!      (in s-coordinates).
89      !!      It is computed using before fields (forward in time) and isopyc-
90      !!      nal or geopotential slopes computed in routine ldfslp.
91      !!      Here, u and v components are considered as 2 independent scalar
92      !!      fields. Therefore, the property of splitting divergent and rota-
93      !!      tional part of the flow of the standard, z-coordinate laplacian
94      !!      momentum diffusion is lost.
95      !!      horizontal fluxes associated with the rotated lateral mixing:
96      !!      u-component:
97      !!         ziut = ( ahmt + ahmb0 ) e2t * e3t / e1t  di[ ub ]
98      !!               -      ahmt       e2t * mi-1(uslp) dk[ mi(mk(ub)) ]
99      !!         zjuf = ( ahmf + ahmb0 ) e1f * e3f / e2f  dj[ ub ]
100      !!               -      ahmf       e1f * mi(vslp)   dk[ mj(mk(ub)) ]
101      !!      v-component:
102      !!         zivf = ( ahmf + ahmb0 ) e2t * e3t / e1t  di[ vb ]
103      !!               -      ahmf       e2t * mj(uslp)   dk[ mi(mk(vb)) ]
104      !!         zjvt = ( ahmt + ahmb0 ) e1f * e3f / e2f  dj[ ub ]
105      !!               -      ahmt       e1f * mj-1(vslp) dk[ mj(mk(vb)) ]
106      !!      take the horizontal divergence of the fluxes:
107      !!         diffu = 1/(e1u*e2u*e3u) {  di  [ ziut ] + dj-1[ zjuf ]  }
108      !!         diffv = 1/(e1v*e2v*e3v) {  di-1[ zivf ] + dj  [ zjvt ]  }
109      !!      Add this trend to the general trend (ua,va):
110      !!         ua = ua + diffu
111      !!      CAUTION: here the isopycnal part is with a coeff. of aht. This
112      !!      should be modified for applications others than orca_r2 (!!bug)
113      !!
114      !! ** Action :
115      !!        Update (ua,va) arrays with the before geopotential biharmonic
116      !!      mixing trend.
117      !!        Update (avmu,avmv) to accompt for the diagonal vertical component
118      !!      of the rotated operator in dynzdf module
119      !!----------------------------------------------------------------------
120      !
121      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
122      !
123      INTEGER  ::   ji, jj, jk   ! dummy loop indices
124      REAL(wp) ::   zabe1, zabe2, zcof1, zcof2                       ! local scalars
125      REAL(wp) ::   zmskt, zmskf, zbu, zbv, zuah, zvah               !   -      -
126      REAL(wp) ::   zcoef0, zcoef3, zcoef4, zmkt, zmkf               !   -      -
127      REAL(wp) ::   zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      -
128      !
129      REAL(wp), POINTER, DIMENSION(:,:) :: ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v
130      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
131      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
132      REAL(KIND=jprb)               :: zhook_handle
133
134      CHARACTER(LEN=*), PARAMETER :: RoutineName='DYN_LDF_ISO'
135
136      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
137
138      !!----------------------------------------------------------------------
139      !
140      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_iso')
141      !
142      CALL wrk_alloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v ) 
143      !
144      IF( kt == nit000 ) THEN
145         IF(lwp) WRITE(numout,*)
146         IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
147         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate horizontal diffusive operator'
148         !                                      ! allocate dyn_ldf_bilap arrays
149         IF( dyn_ldf_iso_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
150      ENDIF
151
152      ! s-coordinate: Iso-level diffusion on tracer, but geopotential level diffusion on momentum
153      IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
154         !
155         DO jk = 1, jpk         ! set the slopes of iso-level
156            DO jj = 2, jpjm1
157               DO ji = 2, jpim1
158                  uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk)
159                  vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk)
160                  wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
161                  wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
162               END DO
163            END DO
164         END DO
165         ! Lateral boundary conditions on the slopes
166         CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
167         CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
168 
169!!bug
170         IF( kt == nit000 ) then
171            IF(lwp) WRITE(numout,*) ' max slop: u', SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)),  &
172               &                             ' wi', sqrt(MAXVAL(wslpi))     , ' wj', sqrt(MAXVAL(wslpj))
173         endif
174!!end
175      ENDIF
176
177      !                                                ! ===============
178      DO jk = 1, jpkm1                                 ! Horizontal slab
179         !                                             ! ===============
180
181         ! Vertical u- and v-shears at level jk and jk+1
182         ! ---------------------------------------------
183         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
184         !                             zdkv(jk=1)=zdkv(jk=2)
185
186         zdk1u(:,:) = ( ub(:,:,jk) -ub(:,:,jk+1) ) * umask(:,:,jk+1)
187         zdk1v(:,:) = ( vb(:,:,jk) -vb(:,:,jk+1) ) * vmask(:,:,jk+1)
188
189         IF( jk == 1 ) THEN
190            zdku(:,:) = zdk1u(:,:)
191            zdkv(:,:) = zdk1v(:,:)
192         ELSE
193            zdku(:,:) = ( ub(:,:,jk-1) - ub(:,:,jk) ) * umask(:,:,jk)
194            zdkv(:,:) = ( vb(:,:,jk-1) - vb(:,:,jk) ) * vmask(:,:,jk)
195         ENDIF
196
197         !                               -----f-----
198         ! Horizontal fluxes on U             | 
199         ! --------------------===        t   u   t
200         !                                    | 
201         ! i-flux at t-point             -----f-----
202
203         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
204            DO jj = 2, jpjm1
205               DO ji = fs_2, jpi   ! vector opt.
206                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj)
207
208                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
209                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
210
211                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
212   
213                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
214                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
215                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
216               END DO
217            END DO
218         ELSE                   ! other coordinate system (zco or sco) : e3t
219            DO jj = 2, jpjm1
220               DO ji = fs_2, jpi   ! vector opt.
221                  zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
222
223                  zmskt = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
224                     &           + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
225
226                  zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
227
228                  ziut(ji,jj) = (  zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) )   &
229                     &           + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
230                     &                      +zdk1u(ji,jj) + zdku (ji-1,jj) )  ) * tmask(ji,jj,jk)
231               END DO
232            END DO
233         ENDIF
234
235         ! j-flux at f-point
236         DO jj = 1, jpjm1
237            DO ji = 1, fs_jpim1   ! vector opt.
238               zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
239
240               zmskf = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
241                  &           + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )
242
243               zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
244
245               zjuf(ji,jj) = (  zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) )   &
246                  &           + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
247                  &                      +zdk1u(ji,jj+1) + zdku (ji,jj) )  ) * fmask(ji,jj,jk)
248            END DO
249         END DO
250
251         !                                |   t   |
252         ! Horizontal fluxes on V         |       |
253         ! --------------------===        f---v---f
254         !                                |       |
255         ! i-flux at f-point              |   t   |
256
257         DO jj = 2, jpjm1
258            DO ji = 1, fs_jpim1   ! vector opt.
259               zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
260
261               zmskf = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
262                  &           + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. )
263
264               zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
265
266               zivf(ji,jj) = (  zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) )   &
267                  &           + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj)     &
268                  &                      +zdk1v(ji,jj) + zdkv (ji+1,jj) )  ) * fmask(ji,jj,jk)
269            END DO
270         END DO
271
272         ! j-flux at t-point
273         IF( ln_zps ) THEN      ! z-coordinate - partial steps : min(e3u)
274            DO jj = 2, jpj
275               DO ji = 1, fs_jpim1   ! vector opt.
276                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj)
277
278                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
279                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
280
281                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
282
283                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
284                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
285                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
286               END DO
287            END DO
288         ELSE                   ! other coordinate system (zco or sco) : e3t
289            DO jj = 2, jpj
290               DO ji = 1, fs_jpim1   ! vector opt.
291                  zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
292
293                  zmskt = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
294                     &           + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
295
296                  zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
297
298                  zjvt(ji,jj) = (  zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) )   &
299                     &           + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj)     &
300                     &                      +zdk1v(ji,jj-1) + zdkv (ji,jj) )  ) * tmask(ji,jj,jk)
301               END DO
302            END DO
303         ENDIF
304
305
306         ! Second derivative (divergence) and add to the general trend
307         ! -----------------------------------------------------------
308
309         DO jj = 2, jpjm1
310            DO ji = 2, jpim1          !! Question vectop possible??? !!bug
311               ! volume elements
312               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
313               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
314               ! horizontal component of isopycnal momentum diffusive trends
315               zuah =( ziut (ji+1,jj) - ziut (ji,jj  ) +   &
316                  &    zjuf (ji  ,jj) - zjuf (ji,jj-1)  ) / zbu
317               zvah =( zivf (ji,jj  ) - zivf (ji-1,jj) +   &
318                  &    zjvt (ji,jj+1) - zjvt (ji,jj  )  ) / zbv
319               ! add the trends to the general trends
320               ua (ji,jj,jk) = ua (ji,jj,jk) + zuah
321               va (ji,jj,jk) = va (ji,jj,jk) + zvah
322            END DO
323         END DO
324         !                                             ! ===============
325      END DO                                           !   End of slab
326      !                                                ! ===============
327
328      ! print sum trends (used for debugging)
329      IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' ldfh - Ua: ', mask1=umask, &
330         &                       tab3d_2=va, clinfo2=       ' Va: ', mask2=vmask, clinfo3='dyn' )
331
332
333      !                                                ! ===============
334      DO jj = 2, jpjm1                                 !  Vertical slab
335         !                                             ! ===============
336
337 
338         ! I. vertical trends associated with the lateral mixing
339         ! =====================================================
340         !  (excluding the vertical flux proportional to dk[t]
341
342
343         ! I.1 horizontal momentum gradient
344         ! --------------------------------
345
346         DO jk = 1, jpk
347            DO ji = 2, jpi
348               ! i-gradient of u at jj
349               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( ub(ji,jj  ,jk) - ub(ji-1,jj  ,jk) )
350               ! j-gradient of u and v at jj
351               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( ub(ji,jj+1,jk) - ub(ji  ,jj  ,jk) )
352               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( vb(ji,jj  ,jk) - vb(ji  ,jj-1,jk) )
353               ! j-gradient of u and v at jj+1
354               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj  ,jk) - ub(ji  ,jj-1,jk) )
355               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji  ,jj  ,jk) )
356            END DO
357         END DO
358         DO jk = 1, jpk
359            DO ji = 1, jpim1
360               ! i-gradient of v at jj
361               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( vb(ji+1,jj,jk) - vb(ji  ,jj  ,jk) )
362            END DO
363         END DO
364
365
366         ! I.2 Vertical fluxes
367         ! -------------------
368
369         ! Surface and bottom vertical fluxes set to zero
370         DO ji = 1, jpi
371            zfuw(ji, 1 ) = 0.e0
372            zfvw(ji, 1 ) = 0.e0
373            zfuw(ji,jpk) = 0.e0
374            zfvw(ji,jpk) = 0.e0
375         END DO
376
377         ! interior (2=<jk=<jpk-1) on U field
378         DO jk = 2, jpkm1
379            DO ji = 2, jpim1
380               zcoef0= 0.5 * aht0 * umask(ji,jj,jk)
381
382               zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
383               zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
384
385               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   &
386                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. )
387               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   &
388                             + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. )
389
390               zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
391               zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
392               ! vertical flux on u field
393               zfuw(ji,jk) = zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)     &
394                                       +zdiu (ji,jk  ) + zdiu (ji+1,jk  ) )   &
395                           + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     &
396                                       +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) )
397               ! update avmu (add isopycnal vertical coefficient to avmu)
398               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
399               avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0
400            END DO
401         END DO
402
403         ! interior (2=<jk=<jpk-1) on V field
404         DO jk = 2, jpkm1
405            DO ji = 2, jpim1
406               zcoef0= 0.5 * aht0 * vmask(ji,jj,jk)
407
408               zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
409               zvwslpj = zcoef0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
410
411               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)   &
412                             + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ), 1. )
413               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)   &
414                             + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ), 1. )
415
416               zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
417               zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
418               ! vertical flux on v field
419               zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     &
420                                       +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
421                           + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
422                                       +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )
423               ! update avmv (add isopycnal vertical coefficient to avmv)
424               ! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
425               avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0
426            END DO
427         END DO
428
429
430         ! I.3 Divergence of vertical fluxes added to the general tracer trend
431         ! -------------------------------------------------------------------
432         DO jk = 1, jpkm1
433            DO ji = 2, jpim1
434               ! volume elements
435               zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
436               zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
437               ! part of the k-component of isopycnal momentum diffusive trends
438               zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu
439               zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv
440               ! add the trends to the general trends
441               ua(ji,jj,jk) = ua(ji,jj,jk) + zuav
442               va(ji,jj,jk) = va(ji,jj,jk) + zvav
443            END DO
444         END DO
445         !                                             ! ===============
446      END DO                                           !   End of slab
447      !                                                ! ===============
448      CALL wrk_dealloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v ) 
449      !
450      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_iso')
451      !
452      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
453   END SUBROUTINE dyn_ldf_iso
454
455# else
456   !!----------------------------------------------------------------------
457   !!   Dummy module                           NO rotation of mixing tensor
458   !!----------------------------------------------------------------------
459CONTAINS
460   SUBROUTINE dyn_ldf_iso( kt )               ! Empty routine
461      INTEGER, INTENT(in) :: kt
462      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
463      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
464      REAL(KIND=jprb)               :: zhook_handle
465
466      CHARACTER(LEN=*), PARAMETER :: RoutineName='DYN_LDF_ISO'
467
468      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
469
470      WRITE(*,*) 'dyn_ldf_iso: You should not have seen this print! error?', kt
471      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
472   END SUBROUTINE dyn_ldf_iso
473#endif
474
475   !!======================================================================
476END MODULE dynldf_iso
Note: See TracBrowser for help on using the repository browser.