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.
ldfeiv_vis.F90 in branches/DEV_r1924_nocs_latphys/NEMO/OPA_SRC/LDF – NEMO

source: branches/DEV_r1924_nocs_latphys/NEMO/OPA_SRC/LDF/ldfeiv_vis.F90 @ 2060

Last change on this file since 2060 was 2060, checked in by sga, 14 years ago

Commits to NEMO branch DEV_r1924_nocs_latphys.

Added Griffies skew flux as formulated by Chris Harris
and modified by George Nurser at NOC. Not tested yet.

Unfinished, expecting to put in new metric terms.

File size: 21.5 KB
Line 
1MODULE ldfeiv_vis
2   !!======================================================================
3   !!                     ***  MODULE  ldfeiv_vis  ***
4   !! Ocean physics:  variable eddy induced velocity coefficients
5   !!======================================================================
6#if   defined key_traldf_eiv   &&   defined key_traldf_c2d
7   !!----------------------------------------------------------------------
8   !!   'key_traldf_eiv'      and                     eddy induced velocity
9   !!   'key_traldf_c2d'                    2D tracer lateral  mixing coef.
10   !!----------------------------------------------------------------------
11   !!   ldf_eiv      : compute the eddy induced velocity coefficients
12   !!----------------------------------------------------------------------
13   !! * Modules used
14   USE oce             ! ocean dynamics and tracers
15   USE dom_oce         ! ocean space and time domain
16   USE sbc_oce         ! surface boundary condition: ocean
17   USE sbcrnf          ! river runoffs
18   USE ldftra_oce      ! ocean tracer   lateral physics
19   USE phycst          ! physical constants
20   USE ldfslp          ! iso-neutral slopes
21   USE in_out_manager  ! I/O manager
22   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
23   USE prtctl          ! Print control
24
25   IMPLICIT NONE
26   PRIVATE
27
28   !! * Routine accessibility
29   PUBLIC ldf_eiv_vis          ! routine called by step.F90
30   !!----------------------------------------------------------------------
31   !!  OPA 9.0 , LOCEAN-IPSL (2005)
32   !! $Id: ldfeiv.F90 1146 2008-06-25 11:42:56Z rblod $
33   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
34   !!----------------------------------------------------------------------
35   !! * Substitutions
36#  include "domzgr_substitute.h90"
37#  include "vectopt_loop_substitute.h90"
38   !!----------------------------------------------------------------------
39
40CONTAINS
41
42# if defined key_mpp_omp
43   !!----------------------------------------------------------------------
44   !!   'key_mpp_omp' :                  OpenMP /  NEC autotasking (j-slab)
45   !!----------------------------------------------------------------------
46
47   SUBROUTINE ldf_eiv_vis( kt )
48      !!----------------------------------------------------------------------
49      !!                  ***  ROUTINE ldf_eiv  ***
50      !!
51      !! ** Purpose :   Compute the eddy induced velocity coefficient from the
52      !!      growth rate of baroclinic instability.
53      !!
54      !! ** Method :
55      !!
56      !! ** Action :   uslp(),   : i- and j-slopes of neutral surfaces
57      !!               vslp()      at u- and v-points, resp.
58      !!               wslpi(),  : i- and j-slopes of neutral surfaces
59      !!               wslpj()     at w-points.
60      !!
61      !! History :
62      !!   8.1  !  99-03  (G. Madec, A. Jouzeau)  Original code
63      !!   8.5  !  02-06  (G. Madec)  Free form, F90
64      !!----------------------------------------------------------------------
65      !! * Arguments
66      INTEGER, INTENT( in ) ::   kt     ! ocean time-step inedx
67
68      !! * Local declarations
69      INTEGER ::   ji, jj, jk, zgap, znsend     ! dummy loop indices
70      REAL(wp) ::   &
71         zfw, ze3w, zn2, zf20,       &  ! temporary scalars
72         zaht, zaht_min, &
73         ztmin1max, zan,zas,zaw,zae,  &
74         zlscale
75      REAL(wp), DIMENSION(jpi,jpj) ::   &
76         zn, zah, zhw, zross, zana, zasa, zawa, zaea, ztmin1, &
77         zindn, zinds, zinde, zindw, znstot, zewtot   ! workspace
78      !!----------------------------------------------------------------------
79
80      IF( kt == nit000 ) THEN
81         IF(lwp) WRITE(numout,*)
82         IF(lwp) WRITE(numout,*) 'ldf_eiv_vis : eddy induced velocity coefficients'
83         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   NEC autotasking / OpenMP : j-slab'
84      ENDIF
85
86      ztmin1max=3.5e-06
87
88      IF ( .NOT. ln_traldf_grif) THEN
89         wslp2(:,:,:)=wslpi(:,:,:) * wslpi(:,:,:) + wslpj(:,:,:) * wslpj(:,:,:)
90      END IF
91
92      !                                                ! ===============
93      DO jj = 2, jpjm1                                 !  Vertical slab
94         !                                             ! ===============
95
96         ! 0. Local initialization
97         ! -----------------------
98         zn   (:,jj) = 0.e0
99         zhw  (:,jj) = 5.e0
100         zah  (:,jj) = 0.e0
101         zross(:,jj) = 0.e0
102
103         ! 1. Compute lateral diffusive coefficient
104         ! ----------------------------------------
105
106!CDIR NOVERRCHK
107         DO jk = 11,27
108!CDIR NOVERRCHK
109            DO ji = 2, jpim1
110               ! Take the max of N^2 and zero then take the vertical sum
111               ! of the square root of the resulting N^2 ( required to compute
112               ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
113               zn2 = MAX( rn2(ji,jj,jk), 0.e0 )
114               ze3w = fse3w(ji,jj,jk) * tmask(ji,jj,jk)
115               zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk)
116               ! Compute elements required for the inverse time scale of baroclinic
117               ! eddies using the isopycnal slopes calculated in ldfslp.F :
118               ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
119               zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w
120               zhw(ji,jj) = zhw(ji,jj) + ze3w
121            END DO
122         END DO
123
124!CDIR NOVERRCHK
125         DO ji = 2, jpim1
126              IF (zhw(ji,jj) > 0) THEN
127               ztmin1(ji,jj)=SQRT( zah(ji,jj) / zhw(ji,jj) )
128             ELSE
129               ztmin1(ji,jj)=0.0
130             ENDIF
131         ENDDO
132      ENDDO
133
134      zindn(:,:)=0.0
135      zinds(:,:)=0.0
136      zindw(:,:)=0.0
137      zinde(:,:)=0.0
138      zana(:,:)=0.0
139      zasa(:,:)=0.0
140      zawa(:,:)=0.0
141      zaea(:,:)=0.0
142
143      DO jj = 2, jpjm1
144         DO ji = 2, jpim1
145             IF ( ztmin1(ji,jj) > ztmin1max ) THEN
146               zan=e2t(ji,jj)/2
147               DO zgap=1,MIN(15,nlcj-1-jj)
148                 IF (ztmin1(ji,jj+zgap) < ztmin1max ) GOTO 100
149                 zan=zan+e2t(ji,jj+zgap)
150               ENDDO
151               IF (zgap == nlcj-jj) zindn(ji,jj)=1
152       100     zas=e2t(ji,jj)/2
153               DO zgap=1,MIN(15,jj-2)
154                 IF (ztmin1(ji,jj-zgap) < ztmin1max ) GOTO 200
155                 zas=zas+e2t(ji,jj-zgap)
156               ENDDO
157               IF (zgap == jj-1) zinds(ji,jj)=1
158       200     zaw=e1t(ji,jj)/2
159               DO zgap=1,MIN(15,ji-2)
160                 IF (ztmin1(ji-zgap,jj) < ztmin1max ) GOTO 300
161                 zaw=zaw+e1t(ji-zgap,jj)
162               ENDDO
163               IF (zgap == ji-1) zindw(ji,jj)=1
164       300     zae=e1t(ji,jj)/2
165               DO zgap=1,MIN(15,nlci-1-ji)
166                 IF (ztmin1(ji+zgap,jj) < ztmin1max ) GOTO 400
167                 zae=zae+e1t(ji+zgap,jj)
168               ENDDO
169               IF (zgap == nlci-ji) zinde(ji,jj)=1
170       400     zana(ji,jj)=zan
171               zasa(ji,jj)=zas
172               zawa(ji,jj)=zaw
173               zaea(ji,jj)=zae
174             ENDIF
175         ENDDO
176       ENDDO
177
178       znsend=nlcj-1
179       ! This code deals with the case of T-pivot at the North fold to ensure the
180       ! correct elements of znstot are set before calling lbc_lnk
181       IF (jperio == 3 .OR. jperio == 4) THEN
182          IF (nproc >= jpnij-jpni) THEN
183             znsend=nlcj-2
184          ENDIF
185       ENDIF
186
187       znstot(:,znsend)=zana(:,znsend)+zasa(:,znsend)
188       znstot(:,2)=zana(:,2)+zasa(:,2)
189       zewtot(nlci-1,:)=zawa(nlci-1,:)+zaea(nlci-1,:)
190       zewtot(2,:)=zawa(2,:)+zaea(2,:)
191
192       CALL lbc_lnk( znstot, 'T', 1. )
193       CALL lbc_lnk( zewtot, 'T', 1. )
194
195       DO jj = 2, jpjm1
196!CDIR NOVERRCHK
197          DO ji = 2, jpim1
198               zana(ji,jj)=zana(ji,jj)+zindn(ji,jj)*znstot(ji,nlcj)
199               zasa(ji,jj)=zasa(ji,jj)+zinds(ji,jj)*znstot(ji,1)
200               zawa(ji,jj)=zawa(ji,jj)+zindw(ji,jj)*zewtot(1,jj)
201               zaea(ji,jj)=zaea(ji,jj)+zinde(ji,jj)*zewtot(nlci,jj)
202
203               zana(ji,jj)=MIN(zana(ji,jj),10*e2t(ji,jj))
204               zasa(ji,jj)=MIN(zasa(ji,jj),10*e2t(ji,jj))
205               zawa(ji,jj)=MIN(zawa(ji,jj),10*e1t(ji,jj))
206               zaea(ji,jj)=MIN(zaea(ji,jj),10*e1t(ji,jj))
207
208               znstot(ji,jj)=zana(ji,jj)+zasa(ji,jj)
209               zewtot(ji,jj)=zawa(ji,jj)+zaea(ji,jj)
210
211               IF ( znstot(ji,jj) == 0) THEN
212                 zlscale=0
213               ELSE
214                 IF ( znstot(ji,jj) < zewtot(ji,jj) ) THEN
215                   zlscale=(MIN(zana(ji,jj),zasa(ji,jj))/MAX(zana(ji,jj),zasa(ji,jj)))*znstot(ji,jj)
216                 ELSE
217                   zlscale=(MIN(zaea(ji,jj),zawa(ji,jj))/MAX(zaea(ji,jj),zawa(ji,jj)))*zewtot(ji,jj)
218                 ENDIF
219               ENDIF
220
221               zlscale=MAX(e2t(ji,jj),zlscale)
222               zlscale=MAX(e1t(ji,jj),zlscale)
223
224               aeiw(ji,jj) = 0.015 * (zlscale**2) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1)
225               aeiw(ji,jj) = MIN (aeiw(ji,jj),2000.0)
226               aeiw(ji,jj) = MAX (aeiw(ji,jj),150.0)*tmask(ji,jj,1)
227         END DO
228
229         ! ORCA R05: Take the minimum between aeiw  and 1000m2/s
230         IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05
231            DO ji = 2, jpim1
232               aeiw(ji,jj) = MIN( aeiw(ji,jj), 1000. )
233            END DO
234         ENDIF
235         !                                             ! ===============
236      END DO                                           !   End of slab
237      !                                                ! ===============
238
239      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
240
241      ! lateral boundary condition on aeiw
242      CALL lbc_lnk( aeiw, 'W', 1. )
243
244      ! Average the diffusive coefficient at u- v- points
245      DO jj = 2, jpjm1
246         DO ji = fs_2, fs_jpim1   ! vector opt.
247            aeiu(ji,jj) = .5 * (aeiw(ji,jj) + aeiw(ji+1,jj  ))
248            aeiv(ji,jj) = .5 * (aeiw(ji,jj) + aeiw(ji  ,jj+1))
249         END DO
250      END DO
251      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
252
253      ! lateral boundary condition on aeiu, aeiv
254      CALL lbc_lnk( aeiu, 'U', 1. )
255      CALL lbc_lnk( aeiv, 'V', 1. )
256
257      IF(ln_ctl)   THEN
258         CALL prt_ctl(tab2d_1=aeiu, clinfo1=' eiv  - u: ', ovlap=1)
259         CALL prt_ctl(tab2d_1=aeiv, clinfo1=' eiv  - v: ', ovlap=1)
260      ENDIF
261
262      ! ORCA R05: add a space variation on aht (=aeiv except at the equator and river mouth)
263      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN
264         zf20     = 2. * omega * SIN( rad * 20. )
265         zaht_min = 100.                              ! minimum value for aht
266         DO jj = 1, jpj
267            DO ji = 1, jpi
268               zaht      = ( 1. -  MIN( 1., ABS( ff(ji,jj) / zf20 ) ) ) * ( aht0 - zaht_min )  &
269                  &      + aht0 * upsrnfh(ji,jj)                          ! enhanced near river mouths
270               ahtu(ji,jj) = MAX( zaht_min, aeiu(ji,jj) ) + zaht
271               ahtv(ji,jj) = MAX( zaht_min, aeiv(ji,jj) ) + zaht
272               ahtw(ji,jj) = MAX( zaht_min, aeiw(ji,jj) ) + zaht
273            END DO
274         END DO
275         IF(ln_ctl) THEN
276            CALL prt_ctl(tab2d_1=ahtu, clinfo1=' aht  - u: ', ovlap=1)
277            CALL prt_ctl(tab2d_1=ahtv, clinfo1=' aht  - v: ', ovlap=1)
278            CALL prt_ctl(tab2d_1=ahtw, clinfo1=' aht  - w: ', ovlap=1)
279         ENDIF
280      ENDIF
281
282      IF( aeiv0 == 0.e0 ) THEN
283         aeiu(:,:) = 0.e0
284         aeiv(:,:) = 0.e0
285         aeiw(:,:) = 0.e0
286      ENDIF
287
288   END SUBROUTINE ldf_eiv_vis
289
290# else
291   !!----------------------------------------------------------------------
292   !!   Default key                                             k-j-i loops
293   !!----------------------------------------------------------------------
294
295   SUBROUTINE ldf_eiv_vis( kt )
296      !!----------------------------------------------------------------------
297      !!                  ***  ROUTINE ldf_eiv_vis  ***
298      !!
299      !! ** Purpose :   Compute the eddy induced velocity coefficient from the
300      !!      growth rate of baroclinic instability.
301      !!
302      !! ** Method :
303      !!
304      !! ** Action : - uslp(),  : i- and j-slopes of neutral surfaces
305      !!             - vslp()      at u- and v-points, resp.
306      !!             - wslpi(),  : i- and j-slopes of neutral surfaces
307      !!             - wslpj()     at w-points.
308      !!
309      !! History :
310      !!   8.1  !  99-03  (G. Madec, A. Jouzeau)  Original code
311      !!   8.5  !  02-06  (G. Madec)  Free form, F90
312      !!----------------------------------------------------------------------
313      !! * Arguments
314      INTEGER, INTENT( in ) ::   kt     ! ocean time-step inedx
315
316      !! * Local declarations
317      INTEGER ::   ji, jj, jk, zgap, znsend         ! dummy loop indices
318      REAL(wp) ::   &
319         zfw, ze3w, zn2, zf20,       &  ! temporary scalars
320         zaht, zaht_min, &
321         ztmin1max, zan,zas,zaw,zae,  &
322         zlscale
323      REAL(wp), DIMENSION(jpi,jpj) ::   &
324         zn, zah, zhw, zross, zana, zasa, zawa, zaea, ztmin1, &
325         zindn, zinds, zinde, zindw, znstot, zewtot   ! workspace
326      !!----------------------------------------------------------------------
327
328      IF( kt == nit000 ) THEN
329         IF(lwp) WRITE(numout,*)
330         IF(lwp) WRITE(numout,*) 'ldf_eiv_vis : eddy induced velocity coefficients'
331         IF(lwp) WRITE(numout,*) '~~~~~~~'
332      ENDIF
333
334      ! 0. Local initialization
335      ! -----------------------
336      zn   (:,:) = 0.e0
337      zhw  (:,:) = 5.e0
338      zah  (:,:) = 0.e0
339      zross(:,:) = 0.e0
340
341      ztmin1max=3.5e-06
342
343      IF ( .NOT. ln_traldf_grif) THEN
344         wslp2(:,:,:)=wslpi(:,:,:) * wslpi(:,:,:) + wslpj(:,:,:) * wslpj(:,:,:)
345      END IF
346
347      ! 1. Compute lateral diffusive coefficient
348      ! ----------------------------------------
349
350      DO jk = 11, 27
351#  if defined key_vectopt_loop
352!CDIR NOVERRCHK
353         DO ji = 1, jpij   ! vector opt.
354            ! Take the max of N^2 and zero then take the vertical sum
355            ! of the square root of the resulting N^2 ( required to compute
356            ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
357            zn2 = MAX( rn2(ji,1,jk), 0.e0 )
358            zn(ji,1) = zn(ji,1) + SQRT( zn2 ) * fse3w(ji,1,jk)
359            ! Compute elements required for the inverse time scale of baroclinic
360            ! eddies using the isopycnal slopes calculated in ldfslp.F :
361            ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
362            ze3w = fse3w(ji,1,jk) * tmask(ji,1,jk)
363            zah(ji,1) = zah(ji,1) + zn2 * wslp2(ji,1,jk) * ze3w
364            zhw(ji,1) = zhw(ji,1) + ze3w
365         END DO
366#  else
367         DO jj = 2, jpjm1
368!CDIR NOVERRCHK
369            DO ji = 2, jpim1
370               ! Take the max of N^2 and zero then take the vertical sum
371               ! of the square root of the resulting N^2 ( required to compute
372               ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
373               zn2 = MAX( rn2(ji,jj,jk), 0.e0 )
374               zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk)
375               ! Compute elements required for the inverse time scale of baroclinic
376               ! eddies using the isopycnal slopes calculated in ldfslp.F :
377               ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
378               ze3w = fse3w(ji,jj,jk) * tmask(ji,jj,jk)
379               zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w
380               zhw(ji,jj) = zhw(ji,jj) + ze3w
381            END DO
382         END DO
383#  endif
384      END DO
385
386      zindn(:,:)=0.0
387      zinds(:,:)=0.0
388      zindw(:,:)=0.0
389      zinde(:,:)=0.0
390      zana(:,:)=0.0
391      zasa(:,:)=0.0
392      zawa(:,:)=0.0
393      zaea(:,:)=0.0
394
395      DO jj = 2, jpjm1
396!CDIR NOVERRCHK
397         DO ji = fs_2, fs_jpim1   ! vector opt.
398             ! Need to calculate  zlscale2 by searching E-W and N-S
399             IF (zhw(ji,jj) > 0) THEN
400               ztmin1(ji,jj)=SQRT( zah(ji,jj) / zhw(ji,jj) )
401             ELSE
402               ztmin1(ji,jj)=0.0
403             ENDIF
404         ENDDO
405      ENDDO
406
407      DO jj = 2, jpjm1
408         DO ji = 2, jpim1
409             IF ( ztmin1(ji,jj) > ztmin1max ) THEN
410               zan=e2t(ji,jj)/2
411               DO zgap=1,MIN(15,nlcj-1-jj)
412                 IF (ztmin1(ji,jj+zgap) < ztmin1max ) GOTO 100
413                 zan=zan+e2t(ji,jj+zgap)
414               ENDDO
415               IF (zgap == nlcj-jj) zindn(ji,jj)=1
416       100     zas=e2t(ji,jj)/2
417               DO zgap=1,MIN(15,jj-2)
418                 IF (ztmin1(ji,jj-zgap) < ztmin1max ) GOTO 200
419                 zas=zas+e2t(ji,jj-zgap)
420               ENDDO
421               IF (zgap == jj-1) zinds(ji,jj)=1
422       200     zaw=e1t(ji,jj)/2
423               DO zgap=1,MIN(15,ji-2)
424                 IF (ztmin1(ji-zgap,jj) < ztmin1max ) GOTO 300
425                 zaw=zaw+e1t(ji-zgap,jj)
426               ENDDO
427               IF (zgap == ji-1) zindw(ji,jj)=1
428       300     zae=e1t(ji,jj)/2
429               DO zgap=1,MIN(15,nlci-1-ji)
430                 IF (ztmin1(ji+zgap,jj) < ztmin1max ) GOTO 400
431                 zae=zae+e1t(ji+zgap,jj)
432               ENDDO
433               IF (zgap == nlci-ji) zinde(ji,jj)=1
434       400     zana(ji,jj)=zan
435               zasa(ji,jj)=zas
436               zawa(ji,jj)=zaw
437               zaea(ji,jj)=zae
438             ENDIF
439         ENDDO
440       ENDDO
441
442       znsend=nlcj-1
443       ! This code deals with the case of T-pivot at the North fold to ensure the
444       ! correct elements of znstot are set before calling lbc_lnk
445       IF (jperio == 3 .OR. jperio == 4) THEN
446          IF (nproc >= jpnij-jpni) THEN
447             znsend=nlcj-2
448          ENDIF
449       ENDIF
450
451       znstot(:,znsend)=zana(:,znsend)+zasa(:,znsend)
452       znstot(:,2)=zana(:,2)+zasa(:,2)
453       zewtot(nlci-1,:)=zawa(nlci-1,:)+zaea(nlci-1,:)
454       zewtot(2,:)=zawa(2,:)+zaea(2,:)
455
456       CALL lbc_lnk( znstot, 'T', 1. )
457       CALL lbc_lnk( zewtot, 'T', 1. )
458
459       DO jj = 2, jpjm1
460!CDIR NOVERRCHK
461!          DO ji = fs_2, fs_jpim1   ! vector opt.
462           DO ji = fs_2, jpim1
463               zana(ji,jj)=zana(ji,jj)+zindn(ji,jj)*znstot(ji,nlcj)
464               zasa(ji,jj)=zasa(ji,jj)+zinds(ji,jj)*znstot(ji,1)
465               zawa(ji,jj)=zawa(ji,jj)+zindw(ji,jj)*zewtot(1,jj)
466               zaea(ji,jj)=zaea(ji,jj)+zinde(ji,jj)*zewtot(nlci,jj)
467
468               zana(ji,jj)=MIN(zana(ji,jj),10*e2t(ji,jj))
469               zasa(ji,jj)=MIN(zasa(ji,jj),10*e2t(ji,jj))
470               zawa(ji,jj)=MIN(zawa(ji,jj),10*e1t(ji,jj))
471               zaea(ji,jj)=MIN(zaea(ji,jj),10*e1t(ji,jj))
472
473               znstot(ji,jj)=zana(ji,jj)+zasa(ji,jj)
474               zewtot(ji,jj)=zawa(ji,jj)+zaea(ji,jj)
475
476               IF ( znstot(ji,jj) == 0) THEN
477                 zlscale=0
478               ELSE
479                 IF ( znstot(ji,jj) < zewtot(ji,jj) ) THEN
480                   zlscale=(MIN(zana(ji,jj),zasa(ji,jj))/MAX(zana(ji,jj),zasa(ji,jj)))*znstot(ji,jj)
481                 ELSE
482                   zlscale=(MIN(zaea(ji,jj),zawa(ji,jj))/MAX(zaea(ji,jj),zawa(ji,jj)))*zewtot(ji,jj)
483                 ENDIF
484               ENDIF
485
486               zlscale=MAX(e2t(ji,jj),zlscale)
487               zlscale=MAX(e1t(ji,jj),zlscale)
488
489               aeiw(ji,jj) = 0.015 * (zlscale**2) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1)
490               aeiw(ji,jj) = MIN (aeiw(ji,jj),2000.0)
491               aeiw(ji,jj) = MAX (aeiw(ji,jj),150.0)*tmask(ji,jj,1)
492         END DO
493      END DO
494
495      ! ORCA R05: Take the minimum between aeiw  and aeiv0
496      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN
497         DO jj = 2, jpjm1
498            DO ji = fs_2, fs_jpim1   ! vector opt.
499               aeiw(ji,jj) = MIN( aeiw(ji,jj), aeiv0 )
500            END DO
501         END DO
502      ENDIF
503
504      ! lateral boundary condition on aeiw
505      CALL lbc_lnk( aeiw, 'W', 1. )
506
507      ! Average the diffusive coefficient at u- v- points
508      DO jj = 2, jpjm1
509         DO ji = fs_2, fs_jpim1   ! vector opt.
510            aeiu(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji+1,jj  ) )
511            aeiv(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji  ,jj+1) )
512         END DO
513      END DO
514
515      ! lateral boundary condition on aeiu, aeiv
516      CALL lbc_lnk( aeiu, 'U', 1. )
517      CALL lbc_lnk( aeiv, 'V', 1. )
518
519      IF(ln_ctl)   THEN
520         CALL prt_ctl(tab2d_1=aeiu, clinfo1=' eiv  - u: ', ovlap=1)
521         CALL prt_ctl(tab2d_1=aeiv, clinfo1=' eiv  - v: ', ovlap=1)
522      ENDIF
523
524      ! ORCA R05: add a space variation on aht (=aeiv except at the equator and river mouth)
525      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN
526         zf20     = 2. * omega * SIN( rad * 20. )
527         zaht_min = 100.                              ! minimum value for aht
528         DO jj = 1, jpj
529            DO ji = 1, jpi
530               zaht      = ( 1. -  MIN( 1., ABS( ff(ji,jj) / zf20 ) ) ) * ( aht0 - zaht_min )  &
531                  &      + aht0 * rnfmsk(ji,jj)                          ! enhanced near river mouths
532               ahtu(ji,jj) = MAX( MAX( zaht_min, aeiu(ji,jj) ) + zaht, aht0 )
533               ahtv(ji,jj) = MAX( MAX( zaht_min, aeiv(ji,jj) ) + zaht, aht0 )
534               ahtw(ji,jj) = MAX( MAX( zaht_min, aeiw(ji,jj) ) + zaht, aht0 )
535            END DO
536         END DO
537         IF(ln_ctl) THEN
538            CALL prt_ctl(tab2d_1=ahtu, clinfo1=' aht  - u: ', ovlap=1)
539            CALL prt_ctl(tab2d_1=ahtv, clinfo1=' aht  - v: ', ovlap=1)
540            CALL prt_ctl(tab2d_1=ahtw, clinfo1=' aht  - w: ', ovlap=1)
541         ENDIF
542      ENDIF
543
544      IF( aeiv0 == 0.e0 ) THEN
545         aeiu(:,:) = 0.e0
546         aeiv(:,:) = 0.e0
547         aeiw(:,:) = 0.e0
548      ENDIF
549
550   END SUBROUTINE ldf_eiv_vis
551
552# endif
553
554#else
555   !!----------------------------------------------------------------------
556   !!   Default option                                         Dummy module
557   !!----------------------------------------------------------------------
558CONTAINS
559   SUBROUTINE ldf_eiv_vis( kt )       ! Empty routine
560      WRITE(*,*) 'ldf_eiv_vis: You should not have seen this print! error?', kt
561   END SUBROUTINE ldf_eiv_vis
562#endif
563
564   !!======================================================================
565END MODULE ldfeiv_vis
Note: See TracBrowser for help on using the repository browser.