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/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/LDF – NEMO

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90 @ 2633

Last change on this file since 2633 was 2633, checked in by trackstand2, 13 years ago

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