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 @ 717

Last change on this file since 717 was 717, checked in by smasson, 17 years ago

finalize the first set of modifications related to ticket:3

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