1 | MODULE ldfeke |
---|
2 | !!====================================================================== |
---|
3 | !! *** MODULE ldfeke *** |
---|
4 | !! Ocean physics: Eddy induced velocity coefficient according to the |
---|
5 | !! GEOMETRIC parameterization scheme |
---|
6 | !!===================================================================== |
---|
7 | !! History : 4.0 ! 2017-11 (J. Mak, G. Madec) original code |
---|
8 | !!---------------------------------------------------------------------- |
---|
9 | |
---|
10 | !!---------------------------------------------------------------------- |
---|
11 | !! ldf_eke : time step depth-integrated EKE and update aeiw (and |
---|
12 | !! from that aeiu and aeiv) according to the GEOMETRIC |
---|
13 | !! parameterization scheme |
---|
14 | !! ldf_eke_init : initialization, namelist read, and parameters control |
---|
15 | !! eke_rst : read/write eke restart in ocean restart file |
---|
16 | !!---------------------------------------------------------------------- |
---|
17 | USE oce ! ocean: dynamics and active tracers variables |
---|
18 | USE phycst ! physical constants |
---|
19 | USE dom_oce ! domain: ocean |
---|
20 | USE ldfslp ! lateral physics: slope of iso-neutral surfaces |
---|
21 | USE ldftra ! lateral physics: eddy coefficients |
---|
22 | USE dynldf_lap_blp, ONLY : eke_keS ! source term of eke due to KE dissipation |
---|
23 | USE zdfmxl ! vertical physics: mixed layer |
---|
24 | |
---|
25 | USE dynspg_ts , ONLY : un_adv, vn_adv |
---|
26 | |
---|
27 | ! |
---|
28 | USE in_out_manager ! I/O manager |
---|
29 | USE iom ! I/O manager library |
---|
30 | USE lib_mpp ! MPP library |
---|
31 | USE lbclnk ! ocean lateral boundary conditions (or mpp link) |
---|
32 | USE prtctl ! Print control |
---|
33 | USE timing ! Timing |
---|
34 | USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) |
---|
35 | |
---|
36 | IMPLICIT NONE |
---|
37 | PRIVATE |
---|
38 | |
---|
39 | PUBLIC ldf_eke ! routine called in step module |
---|
40 | PUBLIC ldf_eke_init ! routine called in opa module |
---|
41 | |
---|
42 | |
---|
43 | ! !!** Namelist namldf_eke ** |
---|
44 | REAL(wp) :: rn_ekedis ! dissipation time scale of EKE [days] |
---|
45 | REAL(wp) :: rn_geom ! geometric parameterization master coefficient [-] |
---|
46 | REAL(wp) :: rn_eke_lap ! diffusion of EKE [m2/s] |
---|
47 | REAL(wp) :: rn_eke_init ! initial value of total EKE [m2/s2] |
---|
48 | REAL(wp) :: rn_eke_min ! background value of total EKE [m2/s2] |
---|
49 | REAL(wp) :: rn_ross_min ! tapering based of minimum Rossby radius [m] |
---|
50 | REAL(wp) :: rn_aeiv_min, rn_aeiv_max ! min and max bounds of aeiv coefficient [m2/s] |
---|
51 | REAL(wp) :: rn_SFmin, rn_SFmax ! min and max bounds of Structure Function [-] |
---|
52 | REAL(wp) :: zf0, zbeta ! f0 and beta for computing Rossby speed |
---|
53 | ! |
---|
54 | INTEGER :: nn_eke_opt ! option for what term to include in eddy energy budget |
---|
55 | INTEGER :: nn_eke_dis ! option for taking constant or spatially varying linear dissipation |
---|
56 | ! |
---|
57 | LOGICAL :: ln_adv_wav ! option for having advection by Rossby wave or not |
---|
58 | LOGICAL :: ln_beta_plane ! option for computing long Rossby phase speed |
---|
59 | ! |
---|
60 | REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: eke_b, eke_n, eke_a ! vertical sum of total Eddy Kinetic Energy [m3/s2] |
---|
61 | REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_ekedis ! linear dissipation rate (= 1/rn_ekedis) [ /s ] |
---|
62 | REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zc1, zc_ros, zadv_wav ! 1st baroclinic mode and long Rossby speed [m /s ] |
---|
63 | |
---|
64 | !! * Substitutions |
---|
65 | # include "vectopt_loop_substitute.h90" |
---|
66 | !!---------------------------------------------------------------------- |
---|
67 | !! NEMO/OPA 4.0 , NEMO Consortium (2017) |
---|
68 | !! $Id: ldfeke.F90 7813 2017-03-20 16:17:45Z clem $ |
---|
69 | !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) |
---|
70 | !!---------------------------------------------------------------------- |
---|
71 | CONTAINS |
---|
72 | |
---|
73 | SUBROUTINE ldf_eke( kt ) |
---|
74 | !!---------------------------------------------------------------------- |
---|
75 | !! *** ROUTINE ldf_eke *** |
---|
76 | !! |
---|
77 | !! ** Purpose : implements the GEOMETRIC parameterization |
---|
78 | !! |
---|
79 | !! ** Notes : If nn_aei_ijk_t = 32 then eke and aeiw are BOTH updated |
---|
80 | !! If ln_eke_equ = .true. in namtra_ldfeiv but nn_aei_ijk_t |
---|
81 | !! is something else, then ONLY eddy equation is updated |
---|
82 | !! (but the eddy energy is passive and doesn't do anything) |
---|
83 | !! |
---|
84 | !! ** Method : GEOMETRIC calculates the Gent-McWilliams / eddy induced |
---|
85 | !! velocity coefficient according to |
---|
86 | !! |
---|
87 | !! aeiw = alpha * (\int E dz) / (\int S M^2 / N dz), |
---|
88 | !! |
---|
89 | !! where (\int E dz) is the depth-integrated eddy energy |
---|
90 | !! (at the previous time level), informed by a parameterized |
---|
91 | !! depth-integrated eddy energy, where |
---|
92 | !! |
---|
93 | !! nn_eke_opt = 0 => default: just PE->EKE growth and linear dissipation |
---|
94 | !! ! |
---|
95 | !! = 1 => default + advection by depth-averaged flow |
---|
96 | !! ! |
---|
97 | !! = 2 => default + advection + contribution to EKE from |
---|
98 | !! momentum dissipation |
---|
99 | !! ! |
---|
100 | !! = 88 => ONLY advection |
---|
101 | !! ! |
---|
102 | !! = 99 => ONLY Laplacian diffusion |
---|
103 | !! |
---|
104 | !! S is a structure function, and M and N are horizontal and |
---|
105 | !! vertical buoyancy frequencies |
---|
106 | !! |
---|
107 | !! linear dissipation may be specified by |
---|
108 | !! |
---|
109 | !! nn_eke_dus = 0 => constant |
---|
110 | !! ! |
---|
111 | !! =-20 => read in a geom_diss_2D.nc field (give it in days) |
---|
112 | !! |
---|
113 | !! ** Action : * time step depth-integrated eddy energy budget |
---|
114 | !! * calculate aeiw |
---|
115 | !! |
---|
116 | !! References : Marshall, Maddison, Berloff JPO 2012 ; Mak et al JPO 2018 |
---|
117 | !!---------------------------------------------------------------------- |
---|
118 | INTEGER, INTENT(in) :: kt ! ocean time step |
---|
119 | ! |
---|
120 | INTEGER :: ji, jj, jk ! dummy loop arguments |
---|
121 | INTEGER :: ik ! local integer |
---|
122 | REAL(wp) :: zn2_min = 1.e-8 ! minimum value of N2 used to compute structure function |
---|
123 | REAL(wp) :: z1_f20, zfw, ze3w ! local scalar |
---|
124 | REAL(wp) :: zfp_ui, zfp_vj ! - - |
---|
125 | REAL(wp) :: zfm_ui, zfm_vj ! - - |
---|
126 | REAL(wp) :: zn_slp2, zn2 ! - - |
---|
127 | REAL(wp) :: zmsku, zaeiu_w ! - - |
---|
128 | REAL(wp) :: zmskv, zaeiv_w ! - - |
---|
129 | REAL(wp) :: zen, zed ! - - |
---|
130 | REAL(wp) :: zeke_rhs ! - - |
---|
131 | REAL(wp) :: zck, zwslpi, zwslpj ! - - tapering near coasts |
---|
132 | REAL(wp) :: zc_rosu ! - - |
---|
133 | REAL(wp), DIMENSION(jpi,jpj) :: zeke_peS, zn_slp ! 2D workspace, PE-KE conversion |
---|
134 | REAL(wp), DIMENSION(jpi,jpj) :: zadv_ubt, zwx, zwy ! - - barotropic advection |
---|
135 | REAL(wp), DIMENSION(jpi,jpj) :: zlap, zaheeu, zaheev, ztu, ztv ! - - diffusion |
---|
136 | REAL(wp), DIMENSION(jpi,jpj) :: zdis ! - - linear dissipation |
---|
137 | REAL(wp), DIMENSION(jpi,jpj) :: zn, zross ! - - tapering near coasts |
---|
138 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zaeiw, zzSF ! 3D workspace |
---|
139 | !!-------------------------------------------------------------------- |
---|
140 | ! |
---|
141 | IF( ln_timing ) CALL timing_start('ldf_eke') |
---|
142 | ! |
---|
143 | IF( kt == nit000 .AND. lwp) THEN !* Control print |
---|
144 | WRITE(numout,*) |
---|
145 | WRITE(numout,*) 'ldf_eke : GEOMETRIC parameterization (total EKE time evolution equation)' |
---|
146 | WRITE(numout,*) '~~~~~~~' |
---|
147 | ENDIF |
---|
148 | |
---|
149 | ! !== EIV mean conversion to EKE (ah N^2 slp2) & N slp ==! |
---|
150 | ! |
---|
151 | ! work out the 3D structure function here |
---|
152 | ! |
---|
153 | ! current: Ferreria et al, zzSF = N^2 / N^2_ref, W-points |
---|
154 | ! capped between rn_SFmax and rn_SFmin |
---|
155 | zzSF(:,:,:) = 0._wp |
---|
156 | DO jk = 2, jpkm1 |
---|
157 | DO jj = 1, jpj |
---|
158 | DO ji = 1, jpi |
---|
159 | IF( jk <= nmln(ji,jj) ) THEN ! above and at mixed layer |
---|
160 | zzSF(ji,jj,jk) = rn_SFmax |
---|
161 | ELSE |
---|
162 | ik = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! one level below the mixed layer (MIN in case ML depth is the ocean depth) |
---|
163 | zzSF(ji,jj,jk) = MAX( 0._wp , rn2b(ji,jj,jk) ) / MAX( zn2_min , rn2b(ji,jj,ik) ) ! Structure Function : N^2 / N^2_ref |
---|
164 | zzSF(ji,jj,jk) = MAX( rn_SFmin , MIN( zzSF(ji,jj,jk) , rn_SFmax ) ) |
---|
165 | ENDIF |
---|
166 | END DO |
---|
167 | END DO |
---|
168 | END DO |
---|
169 | ! |
---|
170 | ! !* parametrized PE_EKE conversion due to eddy induced velocity |
---|
171 | zeke_peS(:,:) = 0._wp |
---|
172 | zn_slp (:,:) = 0._wp |
---|
173 | zn (:,:) = 0._wp |
---|
174 | DO jk = 2, jpkm1 ! query: index? |
---|
175 | DO jj = 2, jpjm1 |
---|
176 | DO ji = fs_2, fs_jpim1 ! vector opt. ! NB: ah_slp2 is w-masked |
---|
177 | zmsku = wmask(ji,jj,jk) / MAX( umask(ji ,jj,jk-1) + umask(ji-1,jj,jk) & |
---|
178 | & + umask(ji-1,jj,jk-1) + umask(ji ,jj,jk) , 1._wp ) |
---|
179 | zmskv = wmask(ji,jj,jk) / MAX( vmask(ji,jj ,jk-1) + vmask(ji,jj-1,jk) & |
---|
180 | & + vmask(ji,jj-1,jk-1) + vmask(ji,jj ,jk) , 1._wp ) |
---|
181 | ! |
---|
182 | zaeiu_w = ( aeiu(ji ,jj,jk-1) + aeiu(ji-1,jj,jk) & |
---|
183 | & + aeiu(ji-1,jj,jk-1) + aeiu(ji ,jj,jk) ) * zmsku |
---|
184 | zaeiv_w = ( aeiv(ji,jj ,jk-1) + aeiv(ji,jj-1,jk) & |
---|
185 | & + aeiv(ji,jj-1,jk-1) + aeiv(ji,jj ,jk) ) * zmskv |
---|
186 | ! |
---|
187 | zn_slp2 = ( zaeiu_w * wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & ! (slope ** 2) * aeiv |
---|
188 | & + zaeiv_w * wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) ! JM 28 Jun: undo slope reduction here too? |
---|
189 | zn2 = MAX( 0._wp , rn2b(ji,jj,jk) ) |
---|
190 | ! |
---|
191 | ze3w = e3w_b(ji,jj,jk) * tmask(ji,jj,jk) |
---|
192 | zn(ji,jj) = zn(ji,jj) + SQRT( zn2 ) * ze3w ! for working out taper at small rossby radius regions later |
---|
193 | ! |
---|
194 | zeke_peS(ji,jj) = zeke_peS(ji,jj) + ze3w * zn2 * zn_slp2 ! note this is >=0 |
---|
195 | ! |
---|
196 | zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & ! taken from ldfslp, undo the slope reduction |
---|
197 | & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 ! near topographic features |
---|
198 | zwslpi = wslpi(ji,jj,jk) / MAX( zck, 0.1_wp) ! (just to avoid dividing by zeros) |
---|
199 | zwslpj = wslpj(ji,jj,jk) / MAX( zck, 0.1_wp) |
---|
200 | ! |
---|
201 | zn_slp(ji,jj) = zn_slp(ji,jj) + ze3w * zzSF(ji,jj,jk) & ! note this >=0 and structure function weighted |
---|
202 | & * SQRT( zn2 * ( zwslpi * zwslpi + zwslpj * zwslpj ) ) |
---|
203 | END DO |
---|
204 | END DO |
---|
205 | END DO |
---|
206 | |
---|
207 | ! !* upstream advection with initial mass fluxes & intermediate update |
---|
208 | ! !* upstream tracer flux in the i and j direction |
---|
209 | DO jj = 1, jpjm1 |
---|
210 | DO ji = 1, fs_jpim1 ! vector opt. |
---|
211 | ! upstream scheme |
---|
212 | zfp_ui = un_adv(ji,jj) + ABS( un_adv(ji,jj) ) |
---|
213 | zfm_ui = un_adv(ji,jj) - ABS( un_adv(ji,jj) ) |
---|
214 | zfp_vj = vn_adv(ji,jj) + ABS( vn_adv(ji,jj) ) |
---|
215 | zfm_vj = vn_adv(ji,jj) - ABS( vn_adv(ji,jj) ) |
---|
216 | zwx(ji,jj) = 0.5 * ( zfp_ui * eke_b(ji,jj) + zfm_ui * eke_b(ji+1,jj ) ) |
---|
217 | zwy(ji,jj) = 0.5 * ( zfp_vj * eke_b(ji,jj) + zfm_vj * eke_b(ji ,jj+1) ) |
---|
218 | END DO |
---|
219 | END DO |
---|
220 | ! !* divergence of ubt advective fluxes |
---|
221 | DO jj = 2, jpjm1 |
---|
222 | DO ji = fs_2, fs_jpim1 ! vector opt. |
---|
223 | zadv_ubt(ji,jj) = - ( zwx(ji,jj) - zwx(ji-1,jj ) & |
---|
224 | & + zwy(ji,jj) - zwy(ji ,jj-1) ) * r1_e1e2t(ji,jj) |
---|
225 | END DO |
---|
226 | END DO |
---|
227 | IF( ln_linssh ) THEN !* top value (linear free surf. only as zwz is multiplied by wmask) |
---|
228 | DO jj = 2, jpjm1 |
---|
229 | DO ji = fs_2, fs_jpim1 ! vector opt. |
---|
230 | zadv_ubt(ji,jj) = - wn(ji,jj,1) * eke_b(ji,jj) / ( ht_0(ji,jj) + 1._wp-tmask(ji,jj,1) ) ! jm (03 Mar 18): was eke_n |
---|
231 | END DO |
---|
232 | END DO |
---|
233 | ENDIF |
---|
234 | ! |
---|
235 | ! !* same as above but for advection by Rossby waves |
---|
236 | zadv_wav(:,:) = 0._wp |
---|
237 | IF ( ln_adv_wav ) THEN |
---|
238 | ! |
---|
239 | DO jj = 2, jpjm1 ! first work out the Rossby phase speeds |
---|
240 | DO ji = 2, fs_jpim1 |
---|
241 | ! compute only for deep enough places |
---|
242 | IF ( ht_0(ji,jj) > 300._wp ) THEN ! jm: should use something like ht_b really... |
---|
243 | ! compute vertical mode phase speed on T point |
---|
244 | zc1(ji,jj) = MIN( 10._wp, zn(ji,jj) / rpi ) |
---|
245 | ! compute long Rossby phase speed on T point (minus sign later) |
---|
246 | IF ( ln_beta_plane ) THEN |
---|
247 | zc_ros(ji,jj) = zc1(ji,jj) * zc1(ji,jj) * zbeta / (zf0 * zf0) |
---|
248 | ELSE |
---|
249 | zc_ros(ji,jj) = zc1(ji,jj) * zc1(ji,jj) * COS( rad * gphit(ji,jj) ) & |
---|
250 | / ( ra * ff_t(ji,jj) * SIN( rad * gphit(ji,jj) ) & |
---|
251 | + rsmall ) |
---|
252 | END IF |
---|
253 | ! cap the Rossby phase speeds by the fastest equatorial Rossby wave speed |
---|
254 | ! the minus sign for westward propagation goes here |
---|
255 | zc_ros(ji,jj) = -MIN( zc1(ji,jj) / 3._wp, zc_ros(ji,jj) ) |
---|
256 | END IF |
---|
257 | END DO |
---|
258 | END DO |
---|
259 | ! |
---|
260 | zwx(:,:) = 0._wp ! wipe the advective contribution from above |
---|
261 | ! |
---|
262 | DO jj = 2, jpjm1 ! average Rossby phase speeds onto U points |
---|
263 | DO ji = 1, fs_jpim1 |
---|
264 | zc_rosu = 0.5 * ( zc_ros(ji,jj) + zc_ros(ji+1,jj) ) * umask(ji,jj,1) |
---|
265 | zfp_ui = zc_rosu + ABS( zc_rosu ) |
---|
266 | zfm_ui = zc_rosu - ABS( zc_rosu ) |
---|
267 | zwx(ji,jj) = 0.5 * ( zfp_ui * eke_b(ji,jj) + zfm_ui * eke_b(ji+1,jj) ) |
---|
268 | END DO |
---|
269 | END DO |
---|
270 | ! !* divergence of wav advective fluxes (only e1 here) |
---|
271 | z1_f20 = 1._wp / ( 2._wp * omega * sin( rad * 20._wp ) ) |
---|
272 | DO jj = 2, jpjm1 |
---|
273 | DO ji = fs_2, fs_jpim1 ! vector opt. |
---|
274 | zadv_wav(ji,jj) = - ( zwx(ji,jj) - zwx(ji-1,jj) ) * r1_e1t(ji,jj) |
---|
275 | zadv_wav(ji,jj) = zadv_wav(ji,jj) * MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) ! tropical decrease |
---|
276 | END DO |
---|
277 | END DO |
---|
278 | END IF |
---|
279 | ! |
---|
280 | !* divergence of diffusive fluxes |
---|
281 | ! code adapted from traldf_lap_blp.F90 |
---|
282 | IF ( rn_eke_lap >= 0._wp ) THEN |
---|
283 | DO jj = 1, jpjm1 |
---|
284 | DO ji = 1, fs_jpim1 ! vector opt. |
---|
285 | zaheeu(ji,jj) = rn_eke_lap * e2_e1u(ji,jj) * umask(ji,jj,1) ! rn_eke_lap is constant (for now) and NOT masked |
---|
286 | zaheev(ji,jj) = rn_eke_lap * e1_e2v(ji,jj) * vmask(ji,jj,1) ! before it is pahu and pahv which IS maked |
---|
287 | END DO |
---|
288 | END DO |
---|
289 | DO jj = 1, jpjm1 !== First derivative (gradient) ==! |
---|
290 | DO ji = 1, fs_jpim1 |
---|
291 | ztu(ji,jj) = zaheeu(ji,jj) * ( eke_b(ji+1,jj ) - eke_b(ji,jj) ) |
---|
292 | ztv(ji,jj) = zaheev(ji,jj) * ( eke_b(ji ,jj+1) - eke_b(ji,jj) ) |
---|
293 | END DO |
---|
294 | END DO |
---|
295 | DO jj = 2, jpjm1 !== Second derivative (divergence), form trend ==! |
---|
296 | DO ji = fs_2, fs_jpim1 |
---|
297 | zlap(ji,jj) = ( ztu(ji,jj) - ztu(ji-1,jj ) & |
---|
298 | & + ztv(ji,jj) - ztv(ji ,jj-1) ) & |
---|
299 | & / e1e2t(ji,jj) |
---|
300 | END DO |
---|
301 | END DO |
---|
302 | ELSE |
---|
303 | zlap(:,:) = 0._wp |
---|
304 | END IF |
---|
305 | !* form the trend for linear dissipation |
---|
306 | DO jj = 2, jpjm1 |
---|
307 | DO ji = fs_2, fs_jpim1 |
---|
308 | zdis(ji,jj) = - r1_ekedis(ji,jj) * (eke_b(ji,jj) - rn_eke_min) * tmask(ji,jj,1) |
---|
309 | END DO |
---|
310 | END DO |
---|
311 | ! !== time stepping of EKE Eq. ==! |
---|
312 | ! |
---|
313 | ! note: the rn_eke_min term is a forcing in eddy equation, thus damping in mean equation |
---|
314 | ! added to prevent overshoots and oscillations of energy from exponential growth/decay |
---|
315 | ! maintains a background (depth-integrated) eddy energy level |
---|
316 | ! |
---|
317 | zeke_rhs = 0._wp |
---|
318 | DO jj = 2, jpjm1 |
---|
319 | DO ji = fs_2, fs_jpim1 ! vector opt. |
---|
320 | ! |
---|
321 | SELECT CASE( nn_eke_opt ) ! Specification of what to include in EKE budget |
---|
322 | ! |
---|
323 | CASE( 0 ) ! default: just PE->EKE growth and linear dissipation |
---|
324 | zeke_rhs = zeke_peS(ji,jj) + zdis(ji,jj) + zlap(ji,jj) |
---|
325 | CASE( 1 ) ! as default but with full advection |
---|
326 | zeke_rhs = - zadv_ubt(ji,jj) + zadv_wav(ji,jj) + zeke_peS(ji,jj) + zdis(ji,jj) + zlap(ji,jj) |
---|
327 | CASE( 2 ) ! full thing with additional KE->EKE growth |
---|
328 | zeke_rhs = - zadv_ubt(ji,jj) + zadv_wav(ji,jj) + zeke_peS(ji,jj) + zdis(ji,jj) + zlap(ji,jj) + eke_keS(ji,jj) |
---|
329 | CASE( 88 ) ! ONLY advection by mean flow |
---|
330 | zeke_rhs = - zadv_ubt(ji,jj) |
---|
331 | CASE( 99 ) ! ONLY diffusion |
---|
332 | zeke_rhs = zlap(ji,jj) |
---|
333 | CASE DEFAULT |
---|
334 | CALL ctl_stop('ldf_eke: wrong choice nn_eke_opt, set at least to 0 (default)') |
---|
335 | END SELECT |
---|
336 | ! |
---|
337 | ! leap-frog |
---|
338 | eke_a(ji,jj) = eke_b(ji,jj) + r2dt * zeke_rhs * ssmask(ji,jj) |
---|
339 | ! |
---|
340 | END DO |
---|
341 | END DO |
---|
342 | CALL lbc_lnk( 'ldfeke', eke_a, 'T', 1. ) ! Lateral boundary conditions on zwi (unchanged sign) |
---|
343 | |
---|
344 | IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step (only swap) |
---|
345 | eke_n(:,:) = eke_a(:,:) |
---|
346 | ! |
---|
347 | ELSE ! Leap-Frog + Asselin filter time stepping |
---|
348 | DO jj = 1, jpj |
---|
349 | DO ji = 1, jpi |
---|
350 | zen = eke_n(ji,jj) |
---|
351 | zed = eke_a(ji,jj) - 2._wp * zen + eke_b(ji,jj) ! time laplacian on tracers |
---|
352 | ! |
---|
353 | eke_b(ji,jj) = zen + atfp * zed ! eke_b <-- filtered eke_n |
---|
354 | eke_n(ji,jj) = eke_a(ji,jj) ! ele_n <-- eke_a |
---|
355 | END DO |
---|
356 | END DO |
---|
357 | ENDIF |
---|
358 | |
---|
359 | ! initialise it here so XIOS stops complaining... |
---|
360 | zross(:,:) = 0._wp |
---|
361 | zaeiw(:,:,:) = 0._wp |
---|
362 | IF( l_eke_eiv ) THEN |
---|
363 | ! !== resulting EIV coefficient ==! |
---|
364 | ! |
---|
365 | ! ! surface value |
---|
366 | z1_f20 = 1._wp / ( 2._wp * omega * sin( rad * 20._wp ) ) |
---|
367 | DO jj = 2, jpjm1 |
---|
368 | DO ji = fs_2, fs_jpim1 ! vector opt. |
---|
369 | ! Rossby radius at w-point, Ro = .5 * sum_jpk(N) / f |
---|
370 | zfw = MAX( ABS( ff_t(ji,jj) ) , 1.e-10 ) |
---|
371 | zross(ji,jj) = 0.5 * zn(ji,jj) / zfw |
---|
372 | ! |
---|
373 | zaeiw(ji,jj,1) = rn_geom * eke_b(ji,jj) / MAX( 1.e-10 , zn_slp(ji,jj) ) * tmask(ji,jj,1) ! zn_slp has zzSF multiplied to it |
---|
374 | zaeiw(ji,jj,1) = MIN( rn_aeiv_max, zaeiw(ji,jj,1) ) ! bound aeiv from above |
---|
375 | zaeiw(ji,jj,1) = zaeiw(ji,jj,1) & |
---|
376 | ! tanh taper to deal with some some large values near coast |
---|
377 | & * 0.5 * ( 1._wp - TANH( ( -ht_0(ji,jj) + 800._wp ) / 300._wp ) ) & |
---|
378 | ! tanh taper of aeiv on internal Rossby radius |
---|
379 | & * 0.5 * ( 1._wp + TANH( ( zross(ji,jj) - rn_ross_min ) / 2._wp ) ) |
---|
380 | zaeiw(ji,jj,1) = zaeiw(ji,jj,1) * MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) ! tropical decrease |
---|
381 | zaeiw(ji,jj,1) = MAX( rn_aeiv_min, zaeiw(ji,jj,1) ) ! bound aeiv from below |
---|
382 | END DO |
---|
383 | END DO |
---|
384 | CALL lbc_lnk( 'ldfeke', zaeiw(:,:,1), 'W', 1. ) |
---|
385 | ! |
---|
386 | ! ! inner value |
---|
387 | ! set bottom to zero and use the un-masked zaeiw(ji,jj,1) first... |
---|
388 | zaeiw(:,:,jpk) = 0._wp |
---|
389 | DO jk = 2, jpkm1 |
---|
390 | DO jj = 1, jpj |
---|
391 | DO ji = 1, jpi |
---|
392 | zaeiw(ji,jj,jk) = zaeiw(ji,jj,1) * zzSF(ji,jj,jk) * wmask(ji,jj,jk) ! zzSF has already been capped |
---|
393 | END DO |
---|
394 | END DO |
---|
395 | END DO |
---|
396 | ! ! aei at u- and v-points |
---|
397 | aeiu(:,:,jpk) = 0._wp |
---|
398 | aeiv(:,:,jpk) = 0._wp |
---|
399 | DO jk = 1, jpkm1 |
---|
400 | DO jj = 2, jpjm1 |
---|
401 | DO ji = fs_2, fs_jpim1 ! vector opt. |
---|
402 | aeiu(ji,jj,jk) = ( zaeiw(ji,jj,jk ) + zaeiw(ji+1,jj,jk ) & |
---|
403 | & + zaeiw(ji,jj,jk+1) + zaeiw(ji+1,jj,jk+1) ) & |
---|
404 | & / MAX( wmask(ji,jj,jk ) + wmask(ji+1,jj,jk ) & |
---|
405 | & + wmask(ji,jj,jk+1) + wmask(ji+1,jj,jk+1) , 1._wp ) * umask(ji,jj,jk) |
---|
406 | aeiv(ji,jj,jk) = ( zaeiw(ji,jj,jk ) + zaeiw(ji,jj+1,jk ) & |
---|
407 | & + zaeiw(ji,jj,jk+1) + zaeiw(ji,jj+1,jk+1) ) & |
---|
408 | & / MAX( wmask(ji,jj,jk ) + wmask(ji,jj+1,jk ) & |
---|
409 | & + wmask(ji,jj,jk+1) + wmask(ji,jj+1,jk+1) , 1._wp ) * vmask(ji,jj,jk) |
---|
410 | END DO |
---|
411 | END DO |
---|
412 | END DO |
---|
413 | ! !-- diagnostics --! |
---|
414 | CALL lbc_lnk_multi( 'ldfeke', aeiu , 'U', 1. , aeiv , 'V', 1. ) ! lateral boundary condition |
---|
415 | CALL lbc_lnk( 'ldfeke', zross, 'W', 1. ) |
---|
416 | END IF |
---|
417 | ! |
---|
418 | IF( lrst_oce ) CALL eke_rst( kt, 'WRITE' ) |
---|
419 | |
---|
420 | ! !== output EKE related variables ==! |
---|
421 | CALL lbc_lnk_multi( 'ldfeke', & ! lateral boundary condition |
---|
422 | & zadv_ubt, 'T', 1. , zadv_wav, 'T', 1. , zdis, 'T', 1. , & |
---|
423 | & zlap , 'T', 1. , zeke_peS, 'T', 1. , & |
---|
424 | & zc1 , 'W', 1. , zc_ros , 'W', 1. ) |
---|
425 | ! |
---|
426 | CALL iom_put( "eke" , eke_n ) ! parameterized total EKE (EPE+ EKE) |
---|
427 | CALL iom_put( "trd_eke_adv_ubt", zadv_ubt ) ! ubt advective trend of EKE(LHS) |
---|
428 | CALL iom_put( "trd_eke_adv_wav", zadv_wav ) ! rossby advective trend of EKE(LHS) |
---|
429 | CALL iom_put( "trd_eke_dis" , zdis ) ! dissipative trend of EKE (RHS) |
---|
430 | CALL iom_put( "trd_eke_lap" , zlap ) ! diffusive trend of EKE (RHS) |
---|
431 | CALL iom_put( "trd_eke_peS" , zeke_peS ) ! PE to EKE source trend (RHS) |
---|
432 | CALL iom_put( "trd_eke_keS" , eke_keS ) ! KE to EKE source trend (RHS) |
---|
433 | CALL iom_put( "aeiv_geom" , zaeiw ) ! eddy induced coefficient from GEOMETRIC param |
---|
434 | CALL iom_put( "rossby_rad" , zross ) ! internal Rossby deformation radius |
---|
435 | CALL iom_put( "c1_vert" , zc1 ) ! 1st baroclinic mode phase speed |
---|
436 | CALL iom_put( "c_ros" , zc_ros ) ! long Rossby phase speed |
---|
437 | !! |
---|
438 | !! jm: modified ldf_tra in ldftra.F90 to include the ln_eke_equ flag |
---|
439 | !! into consideration, otherwise iom_put is called twice aeiu and aeiv |
---|
440 | !! |
---|
441 | CALL iom_put( "aeiu_2d", aeiu(:,:,1) ) ! surface u-EIV coeff. |
---|
442 | CALL iom_put( "aeiv_2d", aeiv(:,:,1) ) ! surface v-EIV coeff. |
---|
443 | CALL iom_put( "aeiu_3d", aeiu(:,:,:) ) ! 3D u-EIV coeff. |
---|
444 | CALL iom_put( "aeiv_3d", aeiv(:,:,:) ) ! 3D v-EIV coeff. |
---|
445 | ! |
---|
446 | IF( ln_timing ) CALL timing_stop('ldf_eke') |
---|
447 | ! |
---|
448 | END SUBROUTINE ldf_eke |
---|
449 | |
---|
450 | SUBROUTINE ldf_eke_init |
---|
451 | !!---------------------------------------------------------------------- |
---|
452 | !! *** ROUTINE ldf_eke_init *** |
---|
453 | !! |
---|
454 | !! ** Purpose : Initialization of the depth-integrated eke, vertical |
---|
455 | !! structure function and max/min caps of aeiw when using the |
---|
456 | !! GEOMETRIC parameterization |
---|
457 | !! |
---|
458 | !! ** Method : Read the namldf_eke namelist and check the parameters |
---|
459 | !! called at the first timestep (nit000) |
---|
460 | !! |
---|
461 | !! ** input : Namlist namldf_eke |
---|
462 | !!---------------------------------------------------------------------- |
---|
463 | INTEGER :: ios, inum |
---|
464 | INTEGER :: ji, jj |
---|
465 | !! |
---|
466 | NAMELIST/namldf_eke/ rn_ekedis , nn_eke_dis , & ! GEOM master params (lambda and option, alpha, e0) |
---|
467 | & rn_geom , rn_eke_init, & |
---|
468 | & rn_eke_lap , & ! coeff of Laplacian diffusion of EKE |
---|
469 | & rn_eke_min , & ! background value of (depth-integrated) EKE |
---|
470 | & rn_ross_min, & ! taper aeiv based on Rossby internal radius |
---|
471 | & rn_aeiv_max, rn_aeiv_min, & ! caps and options on result aeiv |
---|
472 | & rn_SFmin, rn_SFmax, & ! caps on vertical structure function |
---|
473 | & nn_eke_opt , & ! option for EKE budget terms |
---|
474 | & ln_adv_wav , & ! option for advection by Rossby wave |
---|
475 | & ln_beta_plane ! beta plane option for computing Rossby wave speed |
---|
476 | !!---------------------------------------------------------------------- |
---|
477 | ! |
---|
478 | ! |
---|
479 | ! note: GEOMETRIC is *** off *** by default (so no ref list) |
---|
480 | ! REWIND( numnam_ref ) ! Namelist namldf_eke in reference namelist : Turbulent Kinetic Energy |
---|
481 | ! READ ( numnam_ref, namldf_eke, IOSTAT = ios, ERR = 901) |
---|
482 | !901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namldf_eke in reference namelist' ) |
---|
483 | |
---|
484 | REWIND( numnam_cfg ) ! Namelist namldf_eke in configuration namelist : GEOMETRIC parameterization |
---|
485 | READ ( numnam_cfg, namldf_eke, IOSTAT = ios, ERR = 902 ) |
---|
486 | 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namldf_eke in configuration namelist' ) |
---|
487 | IF(lwm) WRITE ( numond, namldf_eke ) |
---|
488 | ! |
---|
489 | IF(lwp) THEN !* Control print |
---|
490 | WRITE(numout,*) |
---|
491 | WRITE(numout,*) 'ldf_eke_init : GEOMETRIC parameterization (total EKE time evolution equation)' |
---|
492 | WRITE(numout,*) '~~~~~~~~~~~~' |
---|
493 | WRITE(numout,*) ' Namelist namldf_eke : set the GEOMETRIC parameters' |
---|
494 | WRITE(numout,*) ' aeiw updated according to GEOMETRIC ln_eke_eiv = ', l_eke_eiv |
---|
495 | WRITE(numout,*) ' dissipation time scale of EKE in days rn_ekedis = ', rn_ekedis, ' days' |
---|
496 | WRITE(numout,*) ' option for dissipation field nn_eke_dis = ', nn_eke_dis |
---|
497 | SELECT CASE( nn_eke_dis ) |
---|
498 | ! |
---|
499 | CASE( 0 ) |
---|
500 | WRITE(numout,*) ' spatially constant ' |
---|
501 | CASE( 2 ) |
---|
502 | WRITE(numout,*) ' read in a 2d varying file geom_diss_2D.nc ' |
---|
503 | END SELECT |
---|
504 | ! |
---|
505 | WRITE(numout,*) ' geometric parameterization master coefficient rn_geom = ', rn_geom |
---|
506 | WRITE(numout,*) ' coeff of Lapican diffusion of EKE rn_eke_lap = ', rn_eke_lap |
---|
507 | WRITE(numout,*) ' initial total EKE value rn_eke_init = ', rn_eke_init |
---|
508 | WRITE(numout,*) ' background values of total EKE rn_eke_min = ', rn_eke_min |
---|
509 | WRITE(numout,*) ' taper aeiv subject to min Rossby radius rn_ross_min = ', rn_ross_min |
---|
510 | WRITE(numout,*) ' maximum bound of aeiv coeff rn_aeiv_max = ', rn_aeiv_max |
---|
511 | WRITE(numout,*) ' minimum bound of aeiv coeff rn_aeiv_min = ', rn_aeiv_min |
---|
512 | WRITE(numout,*) ' minimum bound of Structure Function rn_SFmin = ', rn_SFmin |
---|
513 | WRITE(numout,*) ' maximum bound of Structure Function rn_SFmax = ', rn_SFmax |
---|
514 | WRITE(numout,*) ' [set rnSFmin = rnSFmax = 1 to use aeiv = aeiv(x,y,t)] ' |
---|
515 | WRITE(numout,*) ' option for terms in eddy energy budget nn_eke_opt = ', nn_eke_opt |
---|
516 | WRITE(numout,*) ' default: just PE->EKE growth and linear dissipation ' |
---|
517 | SELECT CASE( nn_eke_opt ) |
---|
518 | ! |
---|
519 | CASE( 0 ) |
---|
520 | WRITE(numout,*) ' default ' |
---|
521 | CASE( 1 ) |
---|
522 | WRITE(numout,*) ' default + advection ' |
---|
523 | CASE( 2 ) |
---|
524 | WRITE(numout,*) ' default + advection + KE->EKE ' |
---|
525 | CASE( 88 ) |
---|
526 | WRITE(numout,*) ' ONLY advection by z-avg mean flow (no growth or dissipation)' |
---|
527 | CASE( 99 ) |
---|
528 | WRITE(numout,*) ' ONLY Laplacian diffusion (no growth or dissipation)' |
---|
529 | END SELECT |
---|
530 | WRITE(numout,*) ' advection of energy by Rossby waves ln_adv_wav = ', ln_adv_wav |
---|
531 | ENDIF |
---|
532 | ! ! allocate eke arrays |
---|
533 | ALLOCATE( eke_b(jpi,jpj) , eke_n(jpi,jpj) , eke_a(jpi,jpj) , eke_keS(jpi,jpj) , & |
---|
534 | & r1_ekedis(jpi,jpj) , zc1(jpi,jpj) , zc_ros(jpi,jpj) , zadv_wav(jpi,jpj) ) |
---|
535 | ! |
---|
536 | SELECT CASE( nn_eke_dis ) ! Specification of linear dissipation |
---|
537 | ! |
---|
538 | CASE( 0 ) !== constant ==! |
---|
539 | IF(lwp) WRITE(numout,*) ' linear EKE dissipation coef. = constant = ', rn_ekedis, ' days' |
---|
540 | r1_ekedis(:,:) = 1._wp / (rn_ekedis * rday) |
---|
541 | ! |
---|
542 | CASE( -20 ) !== fixed horizontal shape read in file ==! |
---|
543 | IF(lwp) WRITE(numout,*) ' linear EKE dissipation coef. = F(i,j) read in geom_diss_2D.nc file' |
---|
544 | CALL iom_open ( 'geom_diss_2D.nc', inum ) |
---|
545 | CALL iom_get ( inum, jpdom_data, 'r1_ekedis', r1_ekedis ) ! read in as time-scale in days... |
---|
546 | CALL iom_close( inum ) |
---|
547 | DO jj = 1, jpj |
---|
548 | DO ji = 1, jpi |
---|
549 | r1_ekedis(ji,jj) = 1._wp / (r1_ekedis(ji,jj) * rday) ! ...convert rate in per second |
---|
550 | END DO |
---|
551 | END DO |
---|
552 | ! |
---|
553 | CASE DEFAULT |
---|
554 | CALL ctl_stop('ldf_eke_init: wrong choice for nn_eke_dis, the option for linear dissipation in GEOMETRIC') |
---|
555 | END SELECT |
---|
556 | ! |
---|
557 | CALL eke_rst( nit000, 'READ' ) !* read or initialize all required files |
---|
558 | ! |
---|
559 | ! tidy up the initilisation of zc1 |
---|
560 | zc_ros(:,:) = 0._wp |
---|
561 | zadv_wav(:,:) = 0._wp |
---|
562 | IF ( ln_adv_wav ) THEN |
---|
563 | DO jj = 2, jpjm1 |
---|
564 | DO ji = 1, fs_jpim1 |
---|
565 | ! set the ones that are not deep enough to be zero |
---|
566 | IF ( ht_0(ji,jj) < 300._wp ) THEN ! jm: should use something like ht_b really... |
---|
567 | zc1(ji,jj) = 0._wp |
---|
568 | END IF |
---|
569 | END DO |
---|
570 | END DO |
---|
571 | END IF |
---|
572 | ! if using beta_plane, compute beta and f0 locally |
---|
573 | IF ( ln_beta_plane ) THEN |
---|
574 | zf0 = ff_f(1,1) |
---|
575 | zbeta = (ff_f(1,2) - zf0) / e2t(1,1) |
---|
576 | IF(lwp) WRITE(numout,*) ' beta plane option is ln_beta_plane =', ln_beta_plane |
---|
577 | IF(lwp) WRITE(numout,*) ' f0 = ', zf0, ' s-1' |
---|
578 | IF(lwp) WRITE(numout,*) ' beta = ', zbeta, 'm-1 s-1' |
---|
579 | END IF |
---|
580 | ! |
---|
581 | ! |
---|
582 | END SUBROUTINE ldf_eke_init |
---|
583 | |
---|
584 | |
---|
585 | SUBROUTINE eke_rst( kt, cdrw ) |
---|
586 | !!--------------------------------------------------------------------- |
---|
587 | !! *** ROUTINE eke_rst *** |
---|
588 | !! |
---|
589 | !! ** Purpose : Read or write EKE file (eke) in restart file |
---|
590 | !! |
---|
591 | !! ** Method : use of IOM library |
---|
592 | !! if the restart does not contain EKE, eke is set |
---|
593 | !! according to rn_eke_init, and aeiu = aeiv = 10 m s^-2 |
---|
594 | !!---------------------------------------------------------------------- |
---|
595 | INTEGER , INTENT(in) :: kt ! ocean time-step |
---|
596 | CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag |
---|
597 | ! |
---|
598 | INTEGER :: jit, jk ! dummy loop indices |
---|
599 | INTEGER :: id1, id2, id3, id4, id5 ! local integer |
---|
600 | !!---------------------------------------------------------------------- |
---|
601 | ! |
---|
602 | IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise |
---|
603 | ! ! --------------- |
---|
604 | IF( ln_rstart ) THEN !* Read the restart file |
---|
605 | id1 = iom_varid( numror, 'eke_b' , ldstop = .FALSE. ) |
---|
606 | id2 = iom_varid( numror, 'eke_n' , ldstop = .FALSE. ) |
---|
607 | id3 = iom_varid( numror, 'aeiu' , ldstop = .FALSE. ) |
---|
608 | id4 = iom_varid( numror, 'aeiv' , ldstop = .FALSE. ) |
---|
609 | IF(lwp) write(numout,*) ' ldfeke init restart' |
---|
610 | ! |
---|
611 | IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN ! all required arrays exist |
---|
612 | IF(lwp) write(numout,*) ' all files for ldfeke exist, loading ...' |
---|
613 | CALL iom_get( numror, jpdom_autoglo, 'eke_b', eke_b ) |
---|
614 | CALL iom_get( numror, jpdom_autoglo, 'eke_n', eke_n ) |
---|
615 | CALL iom_get( numror, jpdom_autoglo, 'aeiu' , aeiu ) |
---|
616 | CALL iom_get( numror, jpdom_autoglo, 'aeiv' , aeiv ) |
---|
617 | ELSE ! one at least array is missing |
---|
618 | IF(lwp) write(numout,*) ' not all files for ldfeke exist ' |
---|
619 | IF(lwp) write(numout,*) ' --- initialize from namelist' |
---|
620 | eke_b(:,:) = rn_eke_init * ht_n(:,:) * ssmask(:,:) |
---|
621 | eke_n(:,:) = eke_b(:,:) |
---|
622 | aeiu(:,:,:) = 10._wp * umask(:,:,:) ! bottom eddy coeff set to zero at last level |
---|
623 | aeiv(:,:,:) = 10._wp * vmask(:,:,:) |
---|
624 | ENDIF |
---|
625 | ! ln_wav_adv restarts + initialisations |
---|
626 | IF( ln_adv_wav ) THEN |
---|
627 | ! id5 = iom_varid( numror, 'zc1' , ldstop = .FALSE. ) |
---|
628 | ! IF( id5 > 0 ) THEN |
---|
629 | ! IF(lwp) write(numout,*) ' ln_adv_wav on and found zc1 variable, loading ...' |
---|
630 | ! CALL iom_get( numror, jpdom_autoglo, 'zc1' , zc1 ) |
---|
631 | ! ELSE |
---|
632 | IF(lwp) write(numout,*) ' ln_adv_wav on but not found zc1 variable, initialising ...' |
---|
633 | zc1(:,:) = 4._wp ! set it to something everywhere but some are set to zero in ldf_eke_init |
---|
634 | ! END IF |
---|
635 | END IF |
---|
636 | ELSE !* Start from rest with a non zero value (required) |
---|
637 | IF(lwp) write(numout,*) ' ldfeke init restart from namelist' |
---|
638 | eke_b(:,:) = rn_eke_init * ht_n(:,:) * ssmask(:,:) |
---|
639 | eke_n(:,:) = eke_b(:,:) |
---|
640 | aeiu(:,:,:) = 10._wp * umask(:,:,:) ! bottom eddy coeff set to zero at last level |
---|
641 | aeiv(:,:,:) = 10._wp * vmask(:,:,:) |
---|
642 | zc1(:,:) = 4._wp |
---|
643 | ! |
---|
644 | ENDIF |
---|
645 | ! |
---|
646 | ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file |
---|
647 | ! ! ------------------- |
---|
648 | IF(lwp) WRITE(numout,*) '---- eke-rst ----' |
---|
649 | CALL iom_rstput( kt, nitrst, numrow, 'eke_b', eke_b ) |
---|
650 | CALL iom_rstput( kt, nitrst, numrow, 'eke_n', eke_n ) |
---|
651 | CALL iom_rstput( kt, nitrst, numrow, 'aeiu' , aeiu ) |
---|
652 | CALL iom_rstput( kt, nitrst, numrow, 'aeiv' , aeiv ) |
---|
653 | CALL iom_rstput( kt, nitrst, numrow, 'zc1' , zc1 ) |
---|
654 | ! |
---|
655 | ENDIF |
---|
656 | ! |
---|
657 | END SUBROUTINE eke_rst |
---|
658 | |
---|
659 | !!====================================================================== |
---|
660 | END MODULE ldfeke |
---|