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

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/LDF/ldfeiv.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 5 years ago

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