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

source: trunk/NEMO/OPA_SRC/LDF/ldfeiv.F90 @ 719

Last change on this file since 719 was 719, checked in by ctlod, 17 years ago

get back to the nemo_v2_3 version for trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.1 KB
Line 
1MODULE ldfeiv
2   !!======================================================================
3   !!                     ***  MODULE  ldfeiv  ***
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   !!                  Same results but not same routine if 'key_mpp_omp'
13   !!                  is defined or not
14   !!----------------------------------------------------------------------
15   !! * Modules used
16   USE oce             ! ocean dynamics and tracers
17   USE dom_oce         ! ocean space and time domain
18   USE ldftra_oce      ! ocean tracer   lateral physics
19   USE phycst          ! physical constants
20   USE ldfslp          ! iso-neutral slopes
21   USE flxrnf          !
22   USE in_out_manager  ! I/O manager
23   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
24   USE prtctl          ! Print control
25
26   IMPLICIT NONE
27   PRIVATE
28   
29   !! * Routine accessibility
30   PUBLIC ldf_eiv               ! routine called by step.F90
31   !!----------------------------------------------------------------------
32   !!  OPA 9.0 , LOCEAN-IPSL (2005)
33   !! $Header$
34   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
35   !!----------------------------------------------------------------------
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
40
41CONTAINS
42
43# if defined key_mpp_omp
44   !!----------------------------------------------------------------------
45   !!   'key_mpp_omp' :                  OpenMP /  NEC autotasking (j-slab)
46   !!----------------------------------------------------------------------
47
48   SUBROUTINE ldf_eiv( kt )
49      !!----------------------------------------------------------------------
50      !!                  ***  ROUTINE ldf_eiv  ***
51      !!
52      !! ** Purpose :   Compute the eddy induced velocity coefficient from the
53      !!      growth rate of baroclinic instability.
54      !!
55      !! ** Method :
56      !!
57      !! ** Action :   uslp(),   : i- and j-slopes of neutral surfaces
58      !!               vslp()      at u- and v-points, resp.
59      !!               wslpi(),  : i- and j-slopes of neutral surfaces
60      !!               wslpj()     at w-points.
61      !!
62      !! History :
63      !!   8.1  !  99-03  (G. Madec, A. Jouzeau)  Original code
64      !!   8.5  !  02-06  (G. Madec)  Free form, F90
65      !!----------------------------------------------------------------------
66      !! * Arguments
67      INTEGER, INTENT( in ) ::   kt     ! ocean time-step inedx
68     
69      !! * Local declarations
70      INTEGER ::   ji, jj, jk           ! dummy loop indices
71      REAL(wp) ::   &
72         zfw, ze3w, zn2, zf20,       &  ! temporary scalars
73         zaht, zaht_min
74      REAL(wp), DIMENSION(jpi,jpj) ::   &
75         zn, zah, zhw, zross            ! workspace
76      !!----------------------------------------------------------------------
77
78      IF( kt == nit000 ) THEN
79         IF(lwp) WRITE(numout,*)
80         IF(lwp) WRITE(numout,*) 'ldf_eiv : eddy induced velocity coefficients'
81         IF(lwp) WRITE(numout,*) '~~~~~~~   NEC autotasking / OpenMP : j-slab'
82      ENDIF
83     
84      !                                                ! ===============
85      DO jj = 2, jpjm1                                 !  Vertical slab
86         !                                             ! ===============
87         
88         ! 0. Local initialization
89         ! -----------------------
90         zn   (:,jj) = 0.e0
91         zhw  (:,jj) = 5.e0
92         zah  (:,jj) = 0.e0
93         zross(:,jj) = 0.e0
94         
95         ! 1. Compute lateral diffusive coefficient
96         ! ----------------------------------------
97
98!CDIR NOVERRCHK
99         DO jk = 1, jpk
100!CDIR NOVERRCHK
101            DO ji = 2, jpim1
102               ! Take the max of N^2 and zero then take the vertical sum
103               ! of the square root of the resulting N^2 ( required to compute
104               ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
105               zn2 = MAX( rn2(ji,jj,jk), 0.e0 )
106               ze3w = fse3w(ji,jj,jk) * tmask(ji,jj,jk)
107               zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk)
108               ! Compute elements required for the inverse time scale of baroclinic
109               ! eddies using the isopycnal slopes calculated in ldfslp.F :
110               ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
111               zah(ji,jj) = zah(ji,jj) + zn2   &
112                              * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)    &
113                                + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )   &
114                              * ze3w
115               zhw(ji,jj) = zhw(ji,jj) + ze3w
116            END DO
117         END DO 
118 
119!CDIR NOVERRCHK
120         DO ji = 2, jpim1
121            zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 )
122            ! Rossby radius at w-point taken < 40km and  > 2km
123            zross(ji,jj) = MAX( MIN( .4 * zn(ji,jj) / zfw, 40.e3 ), 2.e3 )
124            ! Compute aeiw by multiplying Ro^2 and T^-1
125            aeiw(ji,jj) = zross(ji,jj) * zross(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1)
126            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA R02
127               ! Take the minimum between aeiw and aeiv0 for depth levels
128               ! lower than 20 (21 in w- point)
129               IF( mbathy(ji,jj) <= 21. ) aeiw(ji,jj) = MIN( aeiw(ji,jj), 1000. )
130            ENDIF
131         END DO
132
133         ! Decrease the coefficient in the tropics (20N-20S)
134         zf20 = 2. * omega * sin( rad * 20. )
135         DO ji = 2, jpim1
136            aeiw(ji,jj) = MIN( 1., ABS( ff(ji,jj) / zf20 ) ) * aeiw(ji,jj)
137         END DO
138 
139         ! ORCA R05: Take the minimum between aeiw and aeiv0
140         IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05
141            DO ji = 2, jpim1
142               aeiw(ji,jj) = MIN( aeiw(ji,jj), aeiv0 )
143            END DO
144         ENDIF
145         !                                             ! ===============
146      END DO                                           !   End of slab
147      !                                                ! ===============
148
149      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
150
151      ! lateral boundary condition on aeiw
152      CALL lbc_lnk( aeiw, 'W', 1. )
153
154      ! Average the diffusive coefficient at u- v- points
155      DO jj = 2, jpjm1
156         DO ji = fs_2, fs_jpim1   ! vector opt.
157            aeiu(ji,jj) = .5 * (aeiw(ji,jj) + aeiw(ji+1,jj  ))
158            aeiv(ji,jj) = .5 * (aeiw(ji,jj) + aeiw(ji  ,jj+1))
159         END DO
160      END DO 
161      !,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
162
163      ! lateral boundary condition on aeiu, aeiv
164      CALL lbc_lnk( aeiu, 'U', 1. )
165      CALL lbc_lnk( aeiv, 'V', 1. )
166
167      IF(ln_ctl)   THEN
168         CALL prt_ctl(tab2d_1=aeiu, clinfo1=' eiv  - u: ', ovlap=1)
169         CALL prt_ctl(tab2d_1=aeiv, clinfo1=' eiv  - v: ', ovlap=1)
170      ENDIF
171     
172      ! ORCA R05: add a space variation on aht (=aeiv except at the equator and river mouth)
173      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN
174         zf20     = 2. * omega * SIN( rad * 20. )
175         zaht_min = 100.                              ! minimum value for aht
176         DO jj = 1, jpj
177            DO ji = 1, jpi
178               zaht      = ( 1. -  MIN( 1., ABS( ff(ji,jj) / zf20 ) ) ) * ( aht0 - zaht_min )  &
179                  &      + aht0 * upsrnfh(ji,jj)                          ! enhanced near river mouths
180               ahtu(ji,jj) = MAX( MAX( zaht_min, aeiu(ji,jj) ) + zaht, aht0 )
181               ahtv(ji,jj) = MAX( MAX( zaht_min, aeiv(ji,jj) ) + zaht, aht0 )
182               ahtw(ji,jj) = MAX( MAX( zaht_min, aeiw(ji,jj) ) + zaht, aht0 )
183            END DO
184         END DO
185         IF(ln_ctl) THEN
186            CALL prt_ctl(tab2d_1=ahtu, clinfo1=' aht  - u: ', ovlap=1)
187            CALL prt_ctl(tab2d_1=ahtv, clinfo1=' aht  - v: ', ovlap=1)
188            CALL prt_ctl(tab2d_1=ahtw, clinfo1=' aht  - w: ', ovlap=1)
189         ENDIF
190      ENDIF
191
192      IF( aeiv0 == 0.e0 ) THEN
193         aeiu(:,:) = 0.e0
194         aeiv(:,:) = 0.e0
195         aeiw(:,:) = 0.e0
196      ENDIF
197
198   END SUBROUTINE ldf_eiv
199
200# else
201   !!----------------------------------------------------------------------
202   !!   Default key                                             k-j-i loops
203   !!----------------------------------------------------------------------
204
205   SUBROUTINE ldf_eiv( kt )
206      !!----------------------------------------------------------------------
207      !!                  ***  ROUTINE ldf_eiv  ***
208      !!
209      !! ** Purpose :   Compute the eddy induced velocity coefficient from the
210      !!      growth rate of baroclinic instability.
211      !!
212      !! ** Method :
213      !!
214      !! ** Action : - uslp(),  : i- and j-slopes of neutral surfaces
215      !!             - vslp()      at u- and v-points, resp.
216      !!             - wslpi(),  : i- and j-slopes of neutral surfaces
217      !!             - wslpj()     at w-points.
218      !!
219      !! History :
220      !!   8.1  !  99-03  (G. Madec, A. Jouzeau)  Original code
221      !!   8.5  !  02-06  (G. Madec)  Free form, F90
222      !!----------------------------------------------------------------------
223      !! * Arguments
224      INTEGER, INTENT( in ) ::   kt     ! ocean time-step inedx
225     
226      !! * Local declarations
227      INTEGER ::   ji, jj, jk           ! dummy loop indices
228      REAL(wp) ::   &
229         zfw, ze3w, zn2, zf20,       &  ! temporary scalars
230         zaht, zaht_min
231      REAL(wp), DIMENSION(jpi,jpj) ::   &
232         zn, zah, zhw, zross            ! workspace
233      !!----------------------------------------------------------------------
234     
235      IF( kt == nit000 ) THEN
236         IF(lwp) WRITE(numout,*)
237         IF(lwp) WRITE(numout,*) 'ldf_eiv : eddy induced velocity coefficients'
238         IF(lwp) WRITE(numout,*) '~~~~~~~'
239      ENDIF
240     
241      ! 0. Local initialization
242      ! -----------------------
243      zn   (:,:) = 0.e0
244      zhw  (:,:) = 5.e0
245      zah  (:,:) = 0.e0
246      zross(:,:) = 0.e0
247
248
249      ! 1. Compute lateral diffusive coefficient
250      ! ----------------------------------------
251
252      DO jk = 1, jpk
253#  if defined key_vectopt_loop  &&  ! defined key_mpp_omp
254!CDIR NOVERRCHK
255         DO ji = 1, jpij   ! vector opt.
256            ! Take the max of N^2 and zero then take the vertical sum
257            ! of the square root of the resulting N^2 ( required to compute
258            ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
259            zn2 = MAX( rn2(ji,1,jk), 0.e0 )
260            zn(ji,1) = zn(ji,1) + SQRT( zn2 ) * fse3w(ji,1,jk)
261            ! Compute elements required for the inverse time scale of baroclinic
262            ! eddies using the isopycnal slopes calculated in ldfslp.F :
263            ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
264            ze3w = fse3w(ji,1,jk) * tmask(ji,1,jk)
265               zah(ji,1) = zah(ji,1) + zn2   &
266                              * ( wslpi(ji,1,jk) * wslpi(ji,1,jk)    &
267                                + wslpj(ji,1,jk) * wslpj(ji,1,jk) )   &
268                              * ze3w
269            zhw(ji,1) = zhw(ji,1) + ze3w
270         END DO
271#  else
272         DO jj = 2, jpjm1
273!CDIR NOVERRCHK
274            DO ji = 2, jpim1
275               ! Take the max of N^2 and zero then take the vertical sum
276               ! of the square root of the resulting N^2 ( required to compute
277               ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
278               zn2 = MAX( rn2(ji,jj,jk), 0.e0 )
279               zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk)
280               ! Compute elements required for the inverse time scale of baroclinic
281               ! eddies using the isopycnal slopes calculated in ldfslp.F :
282               ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
283               ze3w = fse3w(ji,jj,jk) * tmask(ji,jj,jk)
284               zah(ji,jj) = zah(ji,jj) + zn2   &
285                              * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)    &
286                                + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )  &
287                              * ze3w
288               zhw(ji,jj) = zhw(ji,jj) + ze3w
289            END DO
290         END DO 
291#  endif
292      END DO
293
294      DO jj = 2, jpjm1
295!CDIR NOVERRCHK
296         DO ji = fs_2, fs_jpim1   ! vector opt.
297            zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 )
298            ! Rossby radius at w-point taken < 40km and  > 2km
299            zross(ji,jj) = MAX( MIN( .4 * zn(ji,jj) / zfw, 40.e3 ), 2.e3 )
300            ! Compute aeiw by multiplying Ro^2 and T^-1
301            aeiw(ji,jj) = zross(ji,jj) * zross(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1)
302            IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA R02
303               ! Take the minimum between aeiw and aeiv0 for depth levels
304               ! lower than 20 (21 in w- point)
305               IF( mbathy(ji,jj) <= 21. ) aeiw(ji,jj) = MIN( aeiw(ji,jj), 1000. )
306            ENDIF
307         END DO
308      END DO
309
310      ! Decrease the coefficient in the tropics (20N-20S)
311      zf20 = 2. * omega * sin( rad * 20. )
312      DO jj = 2, jpjm1
313         DO ji = fs_2, fs_jpim1   ! vector opt.
314            aeiw(ji,jj) = MIN( 1., ABS( ff(ji,jj) / zf20 ) ) * aeiw(ji,jj)
315         END DO
316      END DO
317
318      ! ORCA R05: Take the minimum between aeiw  and aeiv0
319      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN
320         DO jj = 2, jpjm1
321            DO ji = fs_2, fs_jpim1   ! vector opt.
322               aeiw(ji,jj) = MIN( aeiw(ji,jj), aeiv0 )
323            END DO
324         END DO
325      ENDIF
326
327      ! lateral boundary condition on aeiw
328      CALL lbc_lnk( aeiw, 'W', 1. )
329
330      ! Average the diffusive coefficient at u- v- points
331      DO jj = 2, jpjm1
332         DO ji = fs_2, fs_jpim1   ! vector opt.
333            aeiu(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji+1,jj  ) )
334            aeiv(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji  ,jj+1) )
335         END DO
336      END DO 
337
338      ! lateral boundary condition on aeiu, aeiv
339      CALL lbc_lnk( aeiu, 'U', 1. )
340      CALL lbc_lnk( aeiv, 'V', 1. )
341
342      IF(ln_ctl)   THEN
343         CALL prt_ctl(tab2d_1=aeiu, clinfo1=' eiv  - u: ', ovlap=1)
344         CALL prt_ctl(tab2d_1=aeiv, clinfo1=' eiv  - v: ', ovlap=1)
345      ENDIF
346
347      ! ORCA R05: add a space variation on aht (=aeiv except at the equator and river mouth)
348      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN
349         zf20     = 2. * omega * SIN( rad * 20. )
350         zaht_min = 100.                              ! minimum value for aht
351         DO jj = 1, jpj
352            DO ji = 1, jpi
353               zaht      = ( 1. -  MIN( 1., ABS( ff(ji,jj) / zf20 ) ) ) * ( aht0 - zaht_min )  &
354                  &      + aht0 * upsrnfh(ji,jj)                          ! enhanced near river mouths
355               ahtu(ji,jj) = MAX( MAX( zaht_min, aeiu(ji,jj) ) + zaht, aht0 )
356               ahtv(ji,jj) = MAX( MAX( zaht_min, aeiv(ji,jj) ) + zaht, aht0 )
357               ahtw(ji,jj) = MAX( MAX( zaht_min, aeiw(ji,jj) ) + zaht, aht0 )
358            END DO
359         END DO
360         IF(ln_ctl) THEN
361            CALL prt_ctl(tab2d_1=ahtu, clinfo1=' aht  - u: ', ovlap=1)
362            CALL prt_ctl(tab2d_1=ahtv, clinfo1=' aht  - v: ', ovlap=1)
363            CALL prt_ctl(tab2d_1=ahtw, clinfo1=' aht  - w: ', ovlap=1)
364         ENDIF
365      ENDIF
366
367      IF( aeiv0 == 0.e0 ) THEN
368         aeiu(:,:) = 0.e0
369         aeiv(:,:) = 0.e0
370         aeiw(:,:) = 0.e0
371      ENDIF
372
373   END SUBROUTINE ldf_eiv
374
375# endif
376
377#else
378   !!----------------------------------------------------------------------
379   !!   Default option                                         Dummy module
380   !!----------------------------------------------------------------------
381CONTAINS
382   SUBROUTINE ldf_eiv( kt )       ! Empty routine
383      WRITE(*,*) 'ldf_eiv: You should not have seen this print! error?', kt
384   END SUBROUTINE ldf_eiv
385#endif
386
387   !!======================================================================
388END MODULE ldfeiv
Note: See TracBrowser for help on using the repository browser.