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

Last change on this file since 106 was 106, checked in by opalod, 20 years ago

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