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_bilapg.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_bilapg.F90

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

The Dr Hook changes from my perl code.
File size: 24.9 KB
Line 
1MODULE dynldf_bilapg
2   !!======================================================================
3   !!                       ***  MODULE  dynldf_bilapg  ***
4   !! Ocean dynamics:  lateral viscosity trend
5   !!======================================================================
6   !! History :  OPA  !  1997-07  (G. Madec)  Original code
7   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module
8   !!            2.0  !  2004-08  (C. Talandier) New trends organization
9   !!----------------------------------------------------------------------
10#if defined key_ldfslp   ||   defined key_esopa
11   !!----------------------------------------------------------------------
12   !!   'key_ldfslp'                              Rotation of mixing tensor
13   !!----------------------------------------------------------------------
14   !!   dyn_ldf_bilapg : update the momentum trend with the horizontal part
15   !!                    of the horizontal s-coord. bilaplacian diffusion
16   !!   ldfguv         :
17   !!----------------------------------------------------------------------
18   USE oce             ! ocean dynamics and tracers
19   USE dom_oce         ! ocean space and time domain
20   USE ldfdyn_oce      ! ocean dynamics lateral physics
21   USE zdf_oce         ! ocean vertical physics
22   USE ldfslp          ! iso-neutral slopes available
23   USE ldftra_oce, ONLY: ln_traldf_iso
24   !
25   USE in_out_manager  ! I/O manager
26   USE lib_mpp         ! MPP library
27   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
28   USE prtctl          ! Print control
29   USE wrk_nemo        ! Memory Allocation
30   USE timing          ! Timing
31
32   USE yomhook, ONLY: lhook, dr_hook
33   USE parkind1, ONLY: jprb, jpim
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   dyn_ldf_bilapg       ! called by step.F90
39
40   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zfuw, zfvw , zdiu, zdiv   ! 2D workspace (ldfguv)
41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::  zdju, zdj1u, zdjv, zdj1v  ! 2D workspace (ldfguv)
42
43   !! * Substitutions
44#  include "domzgr_substitute.h90"
45#  include "ldfdyn_substitute.h90"
46   !!----------------------------------------------------------------------
47   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
48   !! $Id$
49   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51CONTAINS
52
53   INTEGER FUNCTION dyn_ldf_bilapg_alloc()
54   INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
55   INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
56   REAL(KIND=jprb)               :: zhook_handle
57
58   CHARACTER(LEN=*), PARAMETER :: RoutineName='DYN_LDF_BILAPG_ALLOC'
59
60   IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
61
62      !!----------------------------------------------------------------------
63      !!               ***  ROUTINE dyn_ldf_bilapg_alloc  ***
64      !!----------------------------------------------------------------------
65      ALLOCATE( zfuw(jpi,jpk) , zfvw (jpi,jpk) , zdiu(jpi,jpk) , zdiv (jpi,jpk) ,     &
66         &      zdju(jpi,jpk) , zdj1u(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_bilapg_alloc )
67         !
68      IF( dyn_ldf_bilapg_alloc /= 0 )   CALL ctl_warn('dyn_ldf_bilapg_alloc: failed to allocate arrays')
69   IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
70   END FUNCTION dyn_ldf_bilapg_alloc
71
72
73   SUBROUTINE dyn_ldf_bilapg( kt )
74      !!----------------------------------------------------------------------
75      !!                   ***  ROUTINE dyn_ldf_bilapg  ***
76      !!                     
77      !! ** Purpose :   Compute the before trend of the horizontal momentum
78      !!      diffusion and add it to the general trend of momentum equation.
79      !!
80      !! ** Method  :   The lateral momentum diffusive trends is provided by a
81      !!      a 4th order operator rotated along geopotential surfaces. It is
82      !!      computed using before fields (forward in time) and geopotential
83      !!      slopes computed in routine inildf.
84      !!         -1- compute the geopotential harmonic operator applied to
85      !!      (ub,vb) and multiply it by the eddy diffusivity coefficient
86      !!      (done by a call to ldfgpu and ldfgpv routines) The result is in
87      !!      (zwk1,zwk2) arrays. Applied the domain lateral boundary conditions
88      !!      by call to lbc_lnk.
89      !!         -2- applied to (zwk1,zwk2) the geopotential harmonic operator
90      !!      by a second call to ldfgpu and ldfgpv routines respectively. The
91      !!      result is in (zwk3,zwk4) arrays.
92      !!         -3- Add this trend to the general trend (ta,sa):
93      !!            (ua,va) = (ua,va) + (zwk3,zwk4)
94      !!
95      !! ** Action  : - Update (ua,va) arrays with the before geopotential
96      !!                biharmonic mixing trend.
97      !!----------------------------------------------------------------------
98      INTEGER, INTENT( in ) ::   kt           ! ocean time-step index
99      !
100      INTEGER ::   ji, jj, jk                 ! dummy loop indices
101      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zwk1, zwk2, zwk3, zwk4
102      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
103      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
104      REAL(KIND=jprb)               :: zhook_handle
105
106      CHARACTER(LEN=*), PARAMETER :: RoutineName='DYN_LDF_BILAPG'
107
108      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
109
110      !!----------------------------------------------------------------------
111      !
112      IF( nn_timing == 1 )  CALL timing_start('dyn_ldf_bilapg')
113      !
114      CALL wrk_alloc( jpi, jpj, jpk, zwk1, zwk2, zwk3, zwk4 ) 
115      !
116      IF( kt == nit000 ) THEN
117         IF(lwp) WRITE(numout,*)
118         IF(lwp) WRITE(numout,*) 'dyn_ldf_bilapg : horizontal biharmonic operator in s-coordinate'
119         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
120         !                                      ! allocate dyn_ldf_bilapg arrays
121         IF( dyn_ldf_bilapg_alloc() /= 0 )   CALL ctl_stop('STOP', 'dyn_ldf_bilapg: failed to allocate arrays')
122      ENDIF
123
124      ! s-coordinate: Iso-level diffusion on tracer, but geopotential level diffusion on momentum
125      IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
126         !
127         DO jk = 1, jpk         ! set the slopes of iso-level
128            DO jj = 2, jpjm1
129               DO ji = 2, jpim1
130                  uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk)
131                  vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk)
132                  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
133                  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
134               END DO
135            END DO
136         END DO
137         ! Lateral boundary conditions on the slopes
138         CALL lbc_lnk( uslp , 'U', -1. )      ;      CALL lbc_lnk( vslp , 'V', -1. )
139         CALL lbc_lnk( wslpi, 'W', -1. )      ;      CALL lbc_lnk( wslpj, 'W', -1. )
140 
141!!bug
142         IF( kt == nit000 ) then
143            IF(lwp) WRITE(numout,*) ' max slop: u', SQRT( MAXVAL(uslp*uslp)), ' v ', SQRT(MAXVAL(vslp)),  &
144               &                             ' wi', sqrt(MAXVAL(wslpi))     , ' wj', sqrt(MAXVAL(wslpj))
145         endif
146!!end
147      ENDIF
148
149      zwk1(:,:,:) = 0.e0   ;   zwk3(:,:,:) = 0.e0
150      zwk2(:,:,:) = 0.e0   ;   zwk4(:,:,:) = 0.e0
151
152      ! Laplacian of (ub,vb) multiplied by ahm
153      ! -------------------------------------- 
154      CALL ldfguv( ub, vb, zwk1, zwk2, 1 )      ! rotated harmonic operator applied to (ub,vb)
155      !                                         ! and multiply by ahmu, ahmv (output in (zwk1,zwk2) )
156      CALL lbc_lnk( zwk1, 'U', -1. )   ;   CALL lbc_lnk( zwk2, 'V', -1. )     ! Lateral boundary conditions
157
158      ! Bilaplacian of (ub,vb)
159      ! ----------------------
160      CALL ldfguv( zwk1, zwk2, zwk3, zwk4, 2 )  ! rotated harmonic operator applied to (zwk1,zwk2)
161      !                                         ! (output in (zwk3,zwk4) )
162
163      ! Update the momentum trends
164      ! --------------------------
165      DO jj = 2, jpjm1               ! add the diffusive trend to the general momentum trends
166         DO jk = 1, jpkm1
167            DO ji = 2, jpim1
168               ua(ji,jj,jk) = ua(ji,jj,jk) + zwk3(ji,jj,jk)
169               va(ji,jj,jk) = va(ji,jj,jk) + zwk4(ji,jj,jk)
170            END DO
171         END DO
172      END DO
173      !
174      CALL wrk_dealloc( jpi, jpj, jpk, zwk1, zwk2, zwk3, zwk4 ) 
175      !
176      IF( nn_timing == 1 )  CALL timing_stop('dyn_ldf_bilapg')
177      !
178      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
179   END SUBROUTINE dyn_ldf_bilapg
180
181
182   SUBROUTINE ldfguv( pu, pv, plu, plv, kahm )
183      !!----------------------------------------------------------------------
184      !!                  ***  ROUTINE ldfguv  ***
185      !!     
186      !! ** Purpose :   Apply a geopotential harmonic operator to (pu,pv)
187      !!      (defined at u- and v-points) and multiply it by the eddy
188      !!      viscosity coefficient (if kahm=1).
189      !!
190      !! ** Method  :   The harmonic operator rotated along geopotential
191      !!      surfaces is applied to (pu,pv) using the slopes of geopotential
192      !!      surfaces computed in inildf routine. The result is provided in
193      !!      (plu,plv) arrays. It is computed in 2 stepv:
194      !!
195      !!      First step: horizontal part of the operator. It is computed on
196      !!      ==========  pu as follows (idem on pv)
197      !!      horizontal fluxes :
198      !!         zftu = e2u*e3u/e1u di[ pu ] - e2u*uslp dk[ mi(mk(pu)) ]
199      !!         zftv = e1v*e3v/e2v dj[ pu ] - e1v*vslp dk[ mj(mk(pu)) ]
200      !!      take the horizontal divergence of the fluxes (no divided by
201      !!      the volume element :
202      !!         plu  = di-1[ zftu ] +  dj-1[ zftv ]
203      !!
204      !!      Second step: vertical part of the operator. It is computed on
205      !!      ===========  pu as follows (idem on pv)
206      !!      vertical fluxes :
207      !!         zftw = e1t*e2t/e3w * (wslpi^2+wslpj^2)  dk-1[ pu ]
208      !!              -     e2t     *       wslpi        di[ mi(mk(pu)) ]
209      !!              -     e1t     *       wslpj        dj[ mj(mk(pu)) ]
210      !!      take the vertical divergence of the fluxes add it to the hori-
211      !!      zontal component, divide the result by the volume element and
212      !!      if kahm=1, multiply by the eddy diffusivity coefficient:
213      !!         plu  = aht / (e1t*e2t*e3t) { plu + dk[ zftw ] }
214      !!      else:
215      !!         plu  =  1  / (e1t*e2t*e3t) { plu + dk[ zftw ] }
216      !!
217      !! ** Action :
218      !!        plu, plv        : partial harmonic operator applied to
219      !!                          pu and pv (all the components except
220      !!                          second order vertical derivative term)
221      !!----------------------------------------------------------------------
222      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::   pu , pv    ! 1st call: before horizontal velocity
223      !                                                               ! 2nd call: ahm x these fields
224      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(  out) ::   plu, plv   ! partial harmonic operator applied to
225      !                                                               ! pu and pv (all the components except
226      !                                                               ! second order vertical derivative term)
227      INTEGER                         , INTENT(in   ) ::   kahm       ! =1 1st call ; =2 2nd call
228      !
229      INTEGER  ::   ji, jj, jk   ! dummy loop indices
230      REAL(wp) ::   zabe1 , zabe2 , zcof1 , zcof2        ! local scalar
231      REAL(wp) ::   zcoef0, zcoef3, zcoef4               !   -      -
232      REAL(wp) ::   zbur, zbvr, zmkt, zmkf, zuav, zvav   !   -      -
233      REAL(wp) ::   zuwslpi, zuwslpj, zvwslpi, zvwslpj   !   -      -
234      !
235      REAL(wp), POINTER, DIMENSION(:,:) :: ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v
236      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
237      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
238      REAL(KIND=jprb)               :: zhook_handle
239
240      CHARACTER(LEN=*), PARAMETER :: RoutineName='LDFGUV'
241
242      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
243
244      !!----------------------------------------------------------------------
245      !
246      IF( nn_timing == 1 )  CALL timing_start('ldfguv')
247      !
248      CALL wrk_alloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v ) 
249      !
250      !                               ! ********** !   ! ===============
251      DO jk = 1, jpkm1                ! First step !   ! Horizontal slab
252         !                            ! ********** !   ! ===============
253
254         ! I.1 Vertical gradient of pu and pv at level jk and jk+1
255         ! -------------------------------------------------------
256         ! surface boundary condition: zdku(jk=1)=zdku(jk=2)
257         !                             zdkv(jk=1)=zdkv(jk=2)
258
259         zdk1u(:,:) = ( pu(:,:,jk) - pu(:,:,jk+1) ) * umask(:,:,jk+1)
260         zdk1v(:,:) = ( pv(:,:,jk) - pv(:,:,jk+1) ) * vmask(:,:,jk+1)
261
262         IF( jk == 1 ) THEN
263            zdku(:,:) = zdk1u(:,:)
264            zdkv(:,:) = zdk1v(:,:)
265         ELSE
266            zdku(:,:) = ( pu(:,:,jk-1) - pu(:,:,jk) ) * umask(:,:,jk)
267            zdkv(:,:) = ( pv(:,:,jk-1) - pv(:,:,jk) ) * vmask(:,:,jk)
268         ENDIF
269
270         !                                -----f-----
271         ! I.2 Horizontal fluxes on U          |
272         ! ------------------------===     t   u   t
273         !                                     |
274         ! i-flux at t-point              -----f-----
275         DO jj = 1, jpjm1
276            DO ji = 2, jpi
277               zabe1 = e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
278
279               zmkt  = 1./MAX(  umask(ji-1,jj,jk  )+umask(ji,jj,jk+1)   &
280                              + umask(ji-1,jj,jk+1)+umask(ji,jj,jk  ), 1. )
281
282               zcof1 = -e2t(ji,jj) * zmkt   &
283                     * 0.5  * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
284
285               ziut(ji,jj) = tmask(ji,jj,jk) *   &
286                           (  zabe1 * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) )   &
287                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj)     &
288                                       +zdk1u(ji,jj) + zdku (ji-1,jj) )  )
289            END DO
290         END DO
291
292         ! j-flux at f-point
293         DO jj = 1, jpjm1
294            DO ji = 1, jpim1
295               zabe2 = e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
296
297               zmkf  = 1./MAX(  umask(ji,jj+1,jk  )+umask(ji,jj,jk+1)   &
298                              + umask(ji,jj+1,jk+1)+umask(ji,jj,jk  ), 1. )
299
300               zcof2 = -e1f(ji,jj) * zmkf   &
301                     * 0.5  * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
302
303               zjuf(ji,jj) = fmask(ji,jj,jk) *   &
304                           (  zabe2 * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) )   &
305                            + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj)     &
306                                       +zdk1u(ji,jj+1) + zdku (ji,jj) )  )
307            END DO
308         END DO
309
310         !                                 |   t   |
311         ! I.3 Horizontal fluxes on V      |       |
312         ! ------------------------===     f---v---f
313         !                                 |       |
314         ! i-flux at f-point               |   t   |
315         DO jj = 1, jpjm1
316            DO ji = 1, jpim1
317               zabe1 = e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
318
319               zmkf  = 1./MAX(  vmask(ji+1,jj,jk  )+vmask(ji,jj,jk+1)   &
320                              + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk  ), 1. )
321
322               zcof1 = -e2f(ji,jj) * zmkf   &
323                     * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
324
325               zivf(ji,jj) = fmask(ji,jj,jk) *   &
326                           (  zabe1 * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) )   &
327                            + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj)     &
328                                       +zdk1u(ji,jj) + zdku (ji+1,jj) )  )
329            END DO
330         END DO
331
332         ! j-flux at t-point
333         DO jj = 2, jpj
334            DO ji = 1, jpim1
335               zabe2 = e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
336
337               zmkt  = 1./MAX(  vmask(ji,jj-1,jk  )+vmask(ji,jj,jk+1)   &
338                              + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk  ), 1. )
339
340               zcof2 = -e1t(ji,jj) * zmkt   &
341                     * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
342
343               zjvt(ji,jj) = tmask(ji,jj,jk) *   &
344                           (  zabe2 * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) )   &
345                            + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj)     &
346                                       +zdk1u(ji,jj-1) + zdku (ji,jj) )  )
347            END DO
348         END DO
349
350
351         ! I.4 Second derivative (divergence) (not divided by the volume)
352         ! ---------------------
353
354         DO jj = 2, jpjm1
355            DO ji = 2, jpim1
356               plu(ji,jj,jk) = ziut (ji+1,jj) - ziut (ji,jj  )   &
357                             + zjuf (ji  ,jj) - zjuf (ji,jj-1)
358               plv(ji,jj,jk) = zivf (ji,jj  ) - zivf (ji-1,jj)   &
359                             + zjvt (ji,jj+1) - zjvt (ji,jj  ) 
360            END DO
361         END DO
362
363         !                                             ! ===============
364      END DO                                           !   End of slab
365      !                                                ! ===============
366
367      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
368
369      !                             ! ************ !   ! ===============
370      DO jj = 2, jpjm1              !  Second step !   ! Horizontal slab
371         !                          ! ************ !   ! ===============
372
373         ! II.1 horizontal (pu,pv) gradients
374         ! ---------------------------------
375
376         DO jk = 1, jpk
377            DO ji = 2, jpi
378               ! i-gradient of u at jj
379               zdiu (ji,jk) = tmask(ji,jj  ,jk) * ( pu(ji,jj  ,jk) - pu(ji-1,jj  ,jk) )
380               ! j-gradient of u and v at jj
381               zdju (ji,jk) = fmask(ji,jj  ,jk) * ( pu(ji,jj+1,jk) - pu(ji  ,jj  ,jk) )
382               zdjv (ji,jk) = tmask(ji,jj  ,jk) * ( pv(ji,jj  ,jk) - pv(ji  ,jj-1,jk) )
383               ! j-gradient of u and v at jj+1
384               zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( pu(ji,jj  ,jk) - pu(ji  ,jj-1,jk) )
385               zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pv(ji,jj+1,jk) - pv(ji  ,jj  ,jk) )
386            END DO
387         END DO
388         DO jk = 1, jpk
389            DO ji = 1, jpim1
390               ! i-gradient of v at jj
391               zdiv (ji,jk) = fmask(ji,jj  ,jk) * ( pv(ji+1,jj,jk) - pv(ji  ,jj  ,jk) )
392            END DO
393         END DO
394
395
396         ! II.2 Vertical fluxes
397         ! --------------------
398
399         ! Surface and bottom vertical fluxes set to zero
400
401         zfuw(:, 1 ) = 0.e0
402         zfvw(:, 1 ) = 0.e0
403         zfuw(:,jpk) = 0.e0
404         zfvw(:,jpk) = 0.e0
405
406         ! interior (2=<jk=<jpk-1) on pu field
407
408         DO jk = 2, jpkm1
409            DO ji = 2, jpim1
410               ! i- and j-slopes at uw-point
411               zuwslpi = 0.5 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
412               zuwslpj = 0.5 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
413               ! coef. for the vertical dirative
414               zcoef0 = e1u(ji,jj) * e2u(ji,jj) / fse3u(ji,jj,jk)   &
415                      * ( zuwslpi * zuwslpi + zuwslpj * zuwslpj )
416               ! weights for the i-k, j-k averaging at t- and f-points, resp.
417               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1)   &
418                             + tmask(ji,jj,jk  )+tmask(ji+1,jj,jk  ), 1. )
419               zmkf = 1./MAX(  fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1)   &
420                             + fmask(ji,jj-1,jk  )+fmask(ji,jj,jk  ), 1. )
421               ! coef. for the horitontal derivative
422               zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
423               zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
424               ! vertical flux on u field
425               zfuw(ji,jk) = umask(ji,jj,jk) *   &
426                           (  zcoef0 * ( pu  (ji,jj,jk-1) - pu  (ji,jj,jk) )   &
427                            + zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1)     &
428                                        +zdiu (ji,jk  ) + zdiu (ji+1,jk  ) )   &
429                            + zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji  ,jk-1)     &
430                                        +zdj1u(ji,jk  ) + zdju (ji  ,jk  ) ) )
431            END DO
432         END DO
433
434         ! interior (2=<jk=<jpk-1) on pv field
435
436         DO jk = 2, jpkm1
437            DO ji = 2, jpim1
438               ! i- and j-slopes at vw-point
439               zvwslpi = 0.5 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
440               zvwslpj = 0.5 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
441               ! coef. for the vertical derivative
442               zcoef0 = e1v(ji,jj) * e2v(ji,jj) / fse3v(ji,jj,jk)   &
443                      * ( zvwslpi * zvwslpi + zvwslpj * zvwslpj )
444               ! weights for the i-k, j-k averaging at f- and t-points, resp.
445               zmkf = 1./MAX(  fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1)   &
446                             + fmask(ji-1,jj,jk  )+fmask(ji,jj,jk  ), 1. )
447               zmkt = 1./MAX(  tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1)   &
448                             + tmask(ji,jj,jk  )+tmask(ji,jj+1,jk  ), 1. )
449               ! coef. for the horizontal derivatives
450               zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
451               zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
452               ! vertical flux on pv field
453               zfvw(ji,jk) = vmask(ji,jj,jk) *   &
454                           (  zcoef0 * ( pv  (ji,jj,jk-1) - pv  (ji,jj,jk) )   &
455                            + zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1)     &
456                                        +zdiv (ji,jk  ) + zdiv (ji-1,jk  ) )   &
457                            + zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji  ,jk-1)     &
458                                        +zdjv (ji,jk  ) + zdj1v(ji  ,jk  ) )  )
459            END DO
460         END DO
461
462
463         ! II.3 Divergence of vertical fluxes added to the horizontal divergence
464         ! ---------------------------------------------------------------------
465         IF( (kahm -nkahm_smag) ==1 ) THEN
466            ! multiply the laplacian by the eddy viscosity coefficient
467            DO jk = 1, jpkm1
468               DO ji = 2, jpim1
469                  ! eddy coef. divided by the volume element
470                  zbur = fsahmu(ji,jj,jk) / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) )
471                  zbvr = fsahmv(ji,jj,jk) / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) )
472                  ! vertical divergence
473                  zuav = zfuw(ji,jk) - zfuw(ji,jk+1)
474                  zvav = zfvw(ji,jk) - zfvw(ji,jk+1)
475                  ! harmonic operator applied to (pu,pv) and multiply by ahm
476                  plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur
477                  plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr
478               END DO
479            END DO
480         ELSEIF( (kahm +nkahm_smag ) == 2 ) THEN
481            ! second call, no multiplication
482            DO jk = 1, jpkm1
483               DO ji = 2, jpim1
484                  ! inverse of the volume element
485                  zbur = 1. / ( e1u(ji,jj)*e2u(ji,jj)*fse3u(ji,jj,jk) )
486                  zbvr = 1. / ( e1v(ji,jj)*e2v(ji,jj)*fse3v(ji,jj,jk) )
487                  ! vertical divergence
488                  zuav = zfuw(ji,jk) - zfuw(ji,jk+1)
489                  zvav = zfvw(ji,jk) - zfvw(ji,jk+1)
490                  ! harmonic operator applied to (pu,pv)
491                  plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur
492                  plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr
493               END DO
494            END DO
495         ELSE
496           
497            WRITE(numout,*) ' ldfguv: kahm= 1 or 2, here =', kahm
498            WRITE(numout,*) '         We stop'
499            CALL ctl_stop('STOP', 'ldfguv: Unexpected kahm value')
500
501         ENDIF
502         !                                             ! ===============
503      END DO                                           !   End of slab
504      !                                                ! ===============
505
506      CALL wrk_dealloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v ) 
507      !
508      IF( nn_timing == 1 )  CALL timing_stop('ldfguv')
509      !
510      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
511   END SUBROUTINE ldfguv
512
513#else
514   !!----------------------------------------------------------------------
515   !!   Dummy module :                         NO rotation of mixing tensor
516   !!----------------------------------------------------------------------
517CONTAINS
518   SUBROUTINE dyn_ldf_bilapg( kt )               ! Dummy routine
519      INTEGER, INTENT(in) :: kt
520      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
521      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
522      REAL(KIND=jprb)               :: zhook_handle
523
524      CHARACTER(LEN=*), PARAMETER :: RoutineName='DYN_LDF_BILAPG'
525
526      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
527
528      WRITE(*,*) 'dyn_ldf_bilapg: You should not have seen this print! error?', kt
529      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
530   END SUBROUTINE dyn_ldf_bilapg
531#endif
532
533   !!======================================================================
534END MODULE dynldf_bilapg
Note: See TracBrowser for help on using the repository browser.