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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90 @ 2371

Last change on this file since 2371 was 2371, checked in by acc, 13 years ago

nemo_v3_3_beta. Improvement of the Griffies operator code (#680). Some aspects are still undergoing scientific assessment but the code structure is ready for release. Dynamical allocation of the large triad arrays has been introduced to reduce memory when the new operator is not in use.
  • Property svn:keywords set to Id
File size: 11.3 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   !!----------------------------------------------------------------------
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   USE iom
25
26   IMPLICIT NONE
27   PRIVATE
28   
29   !! * Routine accessibility
30   PUBLIC ldf_eiv               ! routine called by step.F90
31   !!----------------------------------------------------------------------
32   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
33   !! $Id$
34   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
35   !!----------------------------------------------------------------------
36   !! * Substitutions
37#  include "domzgr_substitute.h90"
38#  include "vectopt_loop_substitute.h90"
39   !!----------------------------------------------------------------------
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(),  : i- and j-slopes of neutral surfaces
53      !!             - vslp()      at u- and v-points, resp.
54      !!             - wslpi(),  : i- and j-slopes of neutral surfaces
55      !!             - wslpj()     at w-points.
56      !!
57      !! History :
58      !!   8.1  !  99-03  (G. Madec, A. Jouzeau)  Original code
59      !!   8.5  !  02-06  (G. Madec)  Free form, F90
60      !!----------------------------------------------------------------------
61      !! * Arguments
62      INTEGER, INTENT( in ) ::   kt     ! ocean time-step inedx
63     
64      !! * Local declarations
65      INTEGER ::   ji, jj, jk           ! dummy loop indices
66      REAL(wp) ::   &
67         zfw, ze3w, zn2, zf20,       &  ! temporary scalars
68         zaht, zaht_min
69      REAL(wp), DIMENSION(jpi,jpj) ::   &
70         zn, zah, zhw, zross            ! workspace
71      !!----------------------------------------------------------------------
72     
73      IF( kt == nit000 ) THEN
74         IF(lwp) WRITE(numout,*)
75         IF(lwp) WRITE(numout,*) 'ldf_eiv : eddy induced velocity coefficients'
76         IF(lwp) WRITE(numout,*) '~~~~~~~'
77      ENDIF
78     
79      ! 0. Local initialization
80      ! -----------------------
81      zn   (:,:) = 0.e0
82      zhw  (:,:) = 5.e0
83      zah  (:,:) = 0.e0
84      zross(:,:) = 0.e0
85
86
87      ! 1. Compute lateral diffusive coefficient
88      ! ----------------------------------------
89      IF (ln_traldf_grif) THEN
90         DO jk = 1, jpk
91#  if defined key_vectopt_loop 
92            !CDIR NOVERRCHK
93            DO ji = 1, jpij   ! vector opt.
94               ! Take the max of N^2 and zero then take the vertical sum
95               ! of the square root of the resulting N^2 ( required to compute
96               ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
97               zn2 = MAX( rn2b(ji,1,jk), 0.e0 )
98               zn(ji,1) = zn(ji,1) + SQRT( zn2 ) * fse3w(ji,1,jk)
99               ! Compute elements required for the inverse time scale of baroclinic
100               ! eddies using the isopycnal slopes calculated in ldfslp.F :
101               ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
102               ze3w = fse3w(ji,1,jk) * tmask(ji,1,jk)
103               zah(ji,1) = zah(ji,1) + zn2 * wslp2(ji,1,jk) * ze3w
104               zhw(ji,1) = zhw(ji,1) + ze3w
105            END DO
106#  else
107            DO jj = 2, jpjm1
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( rn2b(ji,jj,jk), 0.e0 )
114                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk)
115                  ! Compute elements required for the inverse time scale of baroclinic
116                  ! eddies using the isopycnal slopes calculated in ldfslp.F :
117                  ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
118                  ze3w = fse3w(ji,jj,jk) * tmask(ji,jj,jk)
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#  endif
124         END DO
125      ELSE
126         DO jk = 1, jpk
127#  if defined key_vectopt_loop 
128            !CDIR NOVERRCHK
129            DO ji = 1, jpij   ! vector opt.
130               ! Take the max of N^2 and zero then take the vertical sum
131               ! of the square root of the resulting N^2 ( required to compute
132               ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
133               zn2 = MAX( rn2b(ji,1,jk), 0.e0 )
134               zn(ji,1) = zn(ji,1) + SQRT( zn2 ) * fse3w(ji,1,jk)
135               ! Compute elements required for the inverse time scale of baroclinic
136               ! eddies using the isopycnal slopes calculated in ldfslp.F :
137               ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
138               ze3w = fse3w(ji,1,jk) * tmask(ji,1,jk)
139               zah(ji,1) = zah(ji,1) + zn2   &
140                    * ( wslpi(ji,1,jk) * wslpi(ji,1,jk)    &
141                    + wslpj(ji,1,jk) * wslpj(ji,1,jk) )   &
142                    * ze3w
143               zhw(ji,1) = zhw(ji,1) + ze3w
144            END DO
145#  else
146            DO jj = 2, jpjm1
147               !CDIR NOVERRCHK
148               DO ji = 2, jpim1
149                  ! Take the max of N^2 and zero then take the vertical sum
150                  ! of the square root of the resulting N^2 ( required to compute
151                  ! internal Rossby radius Ro = .5 * sum_jpk(N) / f
152                  zn2 = MAX( rn2b(ji,jj,jk), 0.e0 )
153                  zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * fse3w(ji,jj,jk)
154                  ! Compute elements required for the inverse time scale of baroclinic
155                  ! eddies using the isopycnal slopes calculated in ldfslp.F :
156                  ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w))
157                  ze3w = fse3w(ji,jj,jk) * tmask(ji,jj,jk)
158                  zah(ji,jj) = zah(ji,jj) + zn2   &
159                       * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)    &
160                       + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )  &
161                       * ze3w
162                  zhw(ji,jj) = zhw(ji,jj) + ze3w
163               END DO
164            END DO
165#  endif
166         END DO
167      END IF
168
169      DO jj = 2, jpjm1
170!CDIR NOVERRCHK
171         DO ji = fs_2, fs_jpim1   ! vector opt.
172            zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 )
173            ! Rossby radius at w-point taken < 40km and  > 2km
174            zross(ji,jj) = MAX( MIN( .4 * zn(ji,jj) / zfw, 40.e3 ), 2.e3 )
175            ! Compute aeiw by multiplying Ro^2 and T^-1
176            aeiw(ji,jj) = zross(ji,jj) * zross(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1)
177         END DO
178      END DO
179
180      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN   ! ORCA R02
181         DO jj = 2, jpjm1
182!CDIR NOVERRCHK
183            DO ji = fs_2, fs_jpim1   ! vector opt.
184               ! Take the minimum between aeiw and aeiv0 for depth levels
185               ! lower than 20 (21 in w- point)
186               IF( mbathy(ji,jj) <= 21. ) aeiw(ji,jj) = MIN( aeiw(ji,jj), 1000. )
187            END DO
188         END DO
189      ENDIF
190
191      ! Decrease the coefficient in the tropics (20N-20S)
192      zf20 = 2. * omega * sin( rad * 20. )
193      DO jj = 2, jpjm1
194         DO ji = fs_2, fs_jpim1   ! vector opt.
195            aeiw(ji,jj) = MIN( 1., ABS( ff(ji,jj) / zf20 ) ) * aeiw(ji,jj)
196         END DO
197      END DO
198
199      ! ORCA R05: Take the minimum between aeiw  and aeiv0
200      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN
201         DO jj = 2, jpjm1
202            DO ji = fs_2, fs_jpim1   ! vector opt.
203               aeiw(ji,jj) = MIN( aeiw(ji,jj), aeiv0 )
204            END DO
205         END DO
206      ENDIF
207
208      ! lateral boundary condition on aeiw
209      CALL lbc_lnk( aeiw, 'W', 1. )
210
211      ! Average the diffusive coefficient at u- v- points
212      DO jj = 2, jpjm1
213         DO ji = fs_2, fs_jpim1   ! vector opt.
214            aeiu(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji+1,jj  ) )
215            aeiv(ji,jj) = .5 * ( aeiw(ji,jj) + aeiw(ji  ,jj+1) )
216         END DO
217      END DO 
218
219      ! lateral boundary condition on aeiu, aeiv
220      CALL lbc_lnk( aeiu, 'U', 1. )
221      CALL lbc_lnk( aeiv, 'V', 1. )
222
223      IF(ln_ctl)   THEN
224         CALL prt_ctl(tab2d_1=aeiu, clinfo1=' eiv  - u: ', ovlap=1)
225         CALL prt_ctl(tab2d_1=aeiv, clinfo1=' eiv  - v: ', ovlap=1)
226      ENDIF
227
228      ! ORCA R05: add a space variation on aht (=aeiv except at the equator and river mouth)
229      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN
230         zf20     = 2. * omega * SIN( rad * 20. )
231         zaht_min = 100.                              ! minimum value for aht
232         DO jj = 1, jpj
233            DO ji = 1, jpi
234               zaht      = ( 1. -  MIN( 1., ABS( ff(ji,jj) / zf20 ) ) ) * ( aht0 - zaht_min )  &
235                  &      + aht0 * rnfmsk(ji,jj)                          ! enhanced near river mouths
236               ahtu(ji,jj) = MAX( MAX( zaht_min, aeiu(ji,jj) ) + zaht, aht0 )
237               ahtv(ji,jj) = MAX( MAX( zaht_min, aeiv(ji,jj) ) + zaht, aht0 )
238               ahtw(ji,jj) = MAX( MAX( zaht_min, aeiw(ji,jj) ) + zaht, aht0 )
239            END DO
240         END DO
241         IF(ln_ctl) THEN
242            CALL prt_ctl(tab2d_1=ahtu, clinfo1=' aht  - u: ', ovlap=1)
243            CALL prt_ctl(tab2d_1=ahtv, clinfo1=' aht  - v: ', ovlap=1)
244            CALL prt_ctl(tab2d_1=ahtw, clinfo1=' aht  - w: ', ovlap=1)
245         ENDIF
246      ENDIF
247
248      IF( aeiv0 == 0.e0 ) THEN
249         aeiu(:,:) = 0.e0
250         aeiv(:,:) = 0.e0
251         aeiw(:,:) = 0.e0
252      ENDIF
253
254      CALL iom_put( "aht2d"    , ahtw )   ! lateral eddy diffusivity
255      CALL iom_put( "aht2d_eiv", aeiw )   ! EIV lateral eddy diffusivity
256
257   END SUBROUTINE ldf_eiv
258
259#else
260   !!----------------------------------------------------------------------
261   !!   Default option                                         Dummy module
262   !!----------------------------------------------------------------------
263CONTAINS
264   SUBROUTINE ldf_eiv( kt )       ! Empty routine
265      WRITE(*,*) 'ldf_eiv: You should not have seen this print! error?', kt
266   END SUBROUTINE ldf_eiv
267#endif
268
269   !!======================================================================
270END MODULE ldfeiv
Note: See TracBrowser for help on using the repository browser.