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 branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90 @ 4460

Last change on this file since 4460 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.
  • Property svn:keywords set to Id
File size: 11.9 KB
Line 
1MODULE ldfeiv
2   !!======================================================================
3   !!                     ***  MODULE  ldfeiv  ***
4   !! Ocean physics:  variable eddy induced velocity coefficients
5   !!======================================================================
6   !! History :  OPA  ! 1999-03  (G. Madec, A. Jouzeau)  Original code
7   !!   NEMO     1.0  ! 2002-06  (G. Madec)  Free form, F90
8   !!----------------------------------------------------------------------
9#if   defined key_traldf_eiv   &&   defined key_traldf_c2d
10   !!----------------------------------------------------------------------
11   !!   'key_traldf_eiv'      and                     eddy induced velocity
12   !!   'key_traldf_c2d'                    2D tracer lateral  mixing coef.
13   !!----------------------------------------------------------------------
14   !!   ldf_eiv      : compute the eddy induced velocity coefficients
15   !!----------------------------------------------------------------------
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   USE iom             ! I/O library
27
28   IMPLICIT NONE
29   PRIVATE
30   
31   PUBLIC   ldf_eiv    ! routine called by step.F90
32   
33   !! * Control permutation of array indices
34#  include "oce_ftrans.h90"
35#  include "dom_oce_ftrans.h90"
36#  include "sbc_oce_ftrans.h90"
37#  include "ldftra_oce_ftrans.h90"
38#  include "ldfslp_ftrans.h90"
39
40   !! * Substitutions
41#  include "domzgr_substitute.h90"
42#  include "vectopt_loop_substitute.h90"
43   !!----------------------------------------------------------------------
44   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
45   !! $Id$
46   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
47   !!----------------------------------------------------------------------
48CONTAINS
49
50   SUBROUTINE ldf_eiv( kt )
51      !!----------------------------------------------------------------------
52      !!                  ***  ROUTINE ldf_eiv  ***
53      !!
54      !! ** Purpose :   Compute the eddy induced velocity coefficient from the
55      !!              growth rate of baroclinic instability.
56      !!
57      !! ** Method :
58      !!
59      !! ** Action : - uslp , vslp  : i- and j-slopes of neutral surfaces at u- & v-points
60      !!             - wslpi, wslpj : i- and j-slopes of neutral surfaces at w-points.
61      !!----------------------------------------------------------------------
62      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
63      USE wrk_nemo, ONLY:   zn  => wrk_2d_1 , zah   => wrk_2d_2   ! 2D workspace
64      USE wrk_nemo, ONLY:   zhw => wrk_2d_3 , zross => wrk_2d_4
65      !
66      INTEGER, INTENT(in) ::   kt   ! ocean time-step inedx
67      !
68      INTEGER  ::   ji, jj, jk   ! dummy loop indices
69      REAL(wp) ::   zfw, ze3w, zn2, zf20, zaht, zaht_min      ! temporary scalars
70      !!----------------------------------------------------------------------
71     
72      IF( wrk_in_use(2, 1,2,3,4) ) THEN
73         CALL ctl_stop('ldf_eiv: requested workspace arrays are unavailable.')   ;   RETURN
74      ENDIF
75
76      IF( kt == nit000 ) THEN
77         IF(lwp) WRITE(numout,*)
78         IF(lwp) WRITE(numout,*) 'ldf_eiv : eddy induced velocity coefficients'
79         IF(lwp) WRITE(numout,*) '~~~~~~~'
80      ENDIF
81     
82      ! 0. Local initialization
83      ! -----------------------
84      zn   (:,:) = 0._wp
85      zhw  (:,:) = 5._wp
86      zah  (:,:) = 0._wp
87      zross(:,:) = 0._wp
88
89
90      ! 1. Compute lateral diffusive coefficient
91      ! ----------------------------------------
92      IF( ln_traldf_grif ) THEN
93#  if defined key_vectopt_loop 
94         DO jk = 1, jpk
95!CDIR NOVERRCHK
96            DO ji = 1, jpij   ! vector opt.
97               ! Take the max of N^2 and zero then take the vertical sum
98               ! of the square root of the resulting N^2 ( required to compute
99               ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
100               zn2 = MAX( rn2b(ji,1,jk), 0._wp )
101               zn(ji,1) = zn(ji,1) + SQRT( zn2 ) * fse3w(ji,1,jk)
102               ! Compute elements required for the inverse time scale of baroclinic
103               ! eddies using the isopycnal slopes calculated in ldfslp.F :
104               ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
105               ze3w = fse3w(ji,1,jk) * tmask(ji,1,jk)
106               zah(ji,1) = zah(ji,1) + zn2 * wslp2(ji,1,jk) * ze3w
107               zhw(ji,1) = zhw(ji,1) + ze3w
108            END DO
109         END DO
110#  else
111#     if defined key_z_first
112         DO jj = 2, jpjm1
113            DO ji = 2, jpim1
114               DO jk = 1, jpk
115#     else
116         DO jk = 1, jpk
117            DO jj = 2, jpjm1
118!CDIR NOVERRCHK
119               DO ji = 2, jpim1
120#     endif
121                  ! Take the max of N^2 and zero then take the vertical sum
122                  ! of the square root of the resulting N^2 ( required to compute
123                  ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
124                  zn2 = MAX( rn2b(ji,jj,jk), 0._wp )
125                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk)
126                  ! Compute elements required for the inverse time scale of baroclinic
127                  ! eddies using the isopycnal slopes calculated in ldfslp.F :
128                  ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
129                  ze3w = fse3w(ji,jj,jk) * tmask(ji,jj,jk)
130                  zah(ji,jj) = zah(ji,jj) + zn2 * wslp2(ji,jj,jk) * ze3w
131                  zhw(ji,jj) = zhw(ji,jj) + ze3w
132               END DO
133            END DO
134         END DO
135#  endif
136      ELSE
137#  if defined key_vectopt_loop 
138         DO jk = 1, jpk
139!CDIR NOVERRCHK
140            DO ji = 1, jpij   ! vector opt.
141               ! Take the max of N^2 and zero then take the vertical sum
142               ! of the square root of the resulting N^2 ( required to compute
143               ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
144               zn2 = MAX( rn2b(ji,1,jk), 0._wp )
145               zn(ji,1) = zn(ji,1) + SQRT( zn2 ) * fse3w(ji,1,jk)
146               ! Compute elements required for the inverse time scale of baroclinic
147               ! eddies using the isopycnal slopes calculated in ldfslp.F :
148               ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
149               ze3w = fse3w(ji,1,jk) * tmask(ji,1,jk)
150               zah(ji,1) = zah(ji,1) + zn2 * ( wslpi(ji,1,jk) * wslpi(ji,1,jk)   &
151                  &                          + wslpj(ji,1,jk) * wslpj(ji,1,jk) ) * ze3w
152               zhw(ji,1) = zhw(ji,1) + ze3w
153            END DO
154         END DO
155#  else
156#     if defined key_z_first
157         DO jj = 2, jpjm1
158            DO ji = 2, jpim1
159               DO jk = 1, jpk
160#     else
161         DO jk = 1, jpk
162            DO jj = 2, jpjm1
163!CDIR NOVERRCHK
164               DO ji = 2, jpim1
165#     endif
166                  ! Take the max of N^2 and zero then take the vertical sum
167                  ! of the square root of the resulting N^2 ( required to compute
168                  ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
169                  zn2 = MAX( rn2b(ji,jj,jk), 0._wp )
170                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk)
171                  ! Compute elements required for the inverse time scale of baroclinic
172                  ! eddies using the isopycnal slopes calculated in ldfslp.F :
173                  ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
174                  ze3w = fse3w(ji,jj,jk) * tmask(ji,jj,jk)
175                  zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
176                     &                            + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w
177                  zhw(ji,jj) = zhw(ji,jj) + ze3w
178               END DO
179            END DO
180         END DO
181#  endif
182      END IF
183
184      DO jj = 2, jpjm1
185!CDIR NOVERRCHK
186         DO ji = fs_2, fs_jpim1   ! vector opt.
187            zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 )
188            ! Rossby radius at w-point taken < 40km and  > 2km
189            zross(ji,jj) = MAX( MIN( .4 * zn(ji,jj) / zfw, 40.e3 ), 2.e3 )
190            ! Compute aeiw by multiplying Ro^2 and T^-1
191            aeiw(ji,jj) = zross(ji,jj) * zross(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1)
192         END DO
193      END DO
194
195      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA R2
196         DO jj = 2, jpjm1
197!CDIR NOVERRCHK
198            DO ji = fs_2, fs_jpim1   ! vector opt.
199               ! Take the minimum between aeiw and 1000 m2/s over shelves (depth shallower than 650 m)
200               IF( mbkt(ji,jj) <= 20 )   aeiw(ji,jj) = MIN( aeiw(ji,jj), 1000. )
201            END DO
202         END DO
203      ENDIF
204
205      ! Decrease the coefficient in the tropics (20N-20S)
206      zf20 = 2._wp * omega * sin( rad * 20._wp )
207      DO jj = 2, jpjm1
208         DO ji = fs_2, fs_jpim1   ! vector opt.
209            aeiw(ji,jj) = MIN( 1., ABS( ff(ji,jj) / zf20 ) ) * aeiw(ji,jj)
210         END DO
211      END DO
212
213      ! ORCA R05: Take the minimum between aeiw  and aeiv0
214      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN
215         DO jj = 2, jpjm1
216            DO ji = fs_2, fs_jpim1   ! vector opt.
217               aeiw(ji,jj) = MIN( aeiw(ji,jj), aeiv0 )
218            END DO
219         END DO
220      ENDIF
221      CALL lbc_lnk( aeiw, 'W', 1. )      ! lateral boundary condition on aeiw
222
223
224      ! Average the diffusive coefficient at u- v- points
225      DO jj = 2, jpjm1
226         DO ji = fs_2, fs_jpim1   ! vector opt.
227            aeiu(ji,jj) = 0.5_wp * ( aeiw(ji,jj) + aeiw(ji+1,jj  ) )
228            aeiv(ji,jj) = 0.5_wp * ( aeiw(ji,jj) + aeiw(ji  ,jj+1) )
229         END DO
230      END DO
231      CALL lbc_lnk( aeiu, 'U', 1. )   ;   CALL lbc_lnk( aeiv, 'V', 1. )      ! lateral boundary condition
232
233
234      IF(ln_ctl) THEN
235         CALL prt_ctl(tab2d_1=aeiu, clinfo1=' eiv  - u: ', ovlap=1)
236         CALL prt_ctl(tab2d_1=aeiv, clinfo1=' eiv  - v: ', ovlap=1)
237      ENDIF
238
239      ! ORCA R05: add a space variation on aht (=aeiv except at the equator and river mouth)
240      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN
241         zf20     = 2._wp * omega * SIN( rad * 20._wp )
242         zaht_min = 100._wp                           ! minimum value for aht
243         DO jj = 1, jpj
244            DO ji = 1, jpi
245               zaht      = ( 1._wp -  MIN( 1._wp , ABS( ff(ji,jj) / zf20 ) ) ) * ( aht0 - zaht_min )  &
246                  &      + aht0 * rnfmsk(ji,jj)                          ! enhanced near river mouths
247               ahtu(ji,jj) = MAX( MAX( zaht_min, aeiu(ji,jj) ) + zaht, aht0 )
248               ahtv(ji,jj) = MAX( MAX( zaht_min, aeiv(ji,jj) ) + zaht, aht0 )
249               ahtw(ji,jj) = MAX( MAX( zaht_min, aeiw(ji,jj) ) + zaht, aht0 )
250            END DO
251         END DO
252         IF(ln_ctl) THEN
253            CALL prt_ctl(tab2d_1=ahtu, clinfo1=' aht  - u: ', ovlap=1)
254            CALL prt_ctl(tab2d_1=ahtv, clinfo1=' aht  - v: ', ovlap=1)
255            CALL prt_ctl(tab2d_1=ahtw, clinfo1=' aht  - w: ', ovlap=1)
256         ENDIF
257      ENDIF
258
259      IF( aeiv0 == 0._wp ) THEN
260         aeiu(:,:) = 0._wp
261         aeiv(:,:) = 0._wp
262         aeiw(:,:) = 0._wp
263      ENDIF
264
265      CALL iom_put( "aht2d"    , ahtw )   ! lateral eddy diffusivity
266      CALL iom_put( "aht2d_eiv", aeiw )   ! EIV lateral eddy diffusivity
267      ! 
268      IF( wrk_not_released(2, 1,2,3,4) )   CALL ctl_stop('ldf_eiv: failed to release workspace arrays')
269      !
270   END SUBROUTINE ldf_eiv
271
272#else
273   !!----------------------------------------------------------------------
274   !!   Default option                                         Dummy module
275   !!----------------------------------------------------------------------
276CONTAINS
277   SUBROUTINE ldf_eiv( kt )       ! Empty routine
278      INTEGER :: kt
279      WRITE(*,*) 'ldf_eiv: You should not have seen this print! error?', kt
280   END SUBROUTINE ldf_eiv
281#endif
282
283   !!======================================================================
284END MODULE ldfeiv
Note: See TracBrowser for help on using the repository browser.