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

Last change on this file since 247 was 247, checked in by opalod, 19 years ago

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