| 1 | MODULE dynldf_lap_blp |
|---|
| 2 | !!====================================================================== |
|---|
| 3 | !! *** MODULE dynldf_lap_blp *** |
|---|
| 4 | !! Ocean dynamics: lateral viscosity trend (laplacian and bilaplacian) |
|---|
| 5 | !!====================================================================== |
|---|
| 6 | !! History : 3.7 ! 2014-01 (G. Madec, S. Masson) Original code, re-entrant laplacian |
|---|
| 7 | !! 4.0 ! 2020-04 (A. Nasser, G. Madec) Add symmetric mixing tensor |
|---|
| 8 | !!---------------------------------------------------------------------- |
|---|
| 9 | |
|---|
| 10 | !!---------------------------------------------------------------------- |
|---|
| 11 | !! dyn_ldf_lap : update the momentum trend with the lateral viscosity using an iso-level laplacian operator |
|---|
| 12 | !! dyn_ldf_blp : update the momentum trend with the lateral viscosity using an iso-level bilaplacian operator |
|---|
| 13 | !!---------------------------------------------------------------------- |
|---|
| 14 | USE oce ! ocean dynamics and tracers |
|---|
| 15 | USE dom_oce ! ocean space and time domain |
|---|
| 16 | USE domutl, ONLY : is_tile |
|---|
| 17 | USE ldfdyn ! lateral diffusion: eddy viscosity coef. |
|---|
| 18 | USE ldfslp ! iso-neutral slopes |
|---|
| 19 | USE zdf_oce ! ocean vertical physics |
|---|
| 20 | ! |
|---|
| 21 | USE in_out_manager ! I/O manager |
|---|
| 22 | USE lbclnk ! ocean lateral boundary conditions (or mpp link) |
|---|
| 23 | USE lib_mpp |
|---|
| 24 | #if defined key_loop_fusion |
|---|
| 25 | USE dynldf_lap_blp_lf |
|---|
| 26 | #endif |
|---|
| 27 | |
|---|
| 28 | IMPLICIT NONE |
|---|
| 29 | PRIVATE |
|---|
| 30 | |
|---|
| 31 | PUBLIC dyn_ldf_lap ! called by dynldf.F90 |
|---|
| 32 | PUBLIC dyn_ldf_blp ! called by dynldf.F90 |
|---|
| 33 | |
|---|
| 34 | !! * Substitutions |
|---|
| 35 | # include "do_loop_substitute.h90" |
|---|
| 36 | # include "domzgr_substitute.h90" |
|---|
| 37 | !!---------------------------------------------------------------------- |
|---|
| 38 | !! NEMO/OCE 4.0 , NEMO Consortium (2018) |
|---|
| 39 | !! $Id$ |
|---|
| 40 | !! Software governed by the CeCILL license (see ./LICENSE) |
|---|
| 41 | !!---------------------------------------------------------------------- |
|---|
| 42 | CONTAINS |
|---|
| 43 | |
|---|
| 44 | SUBROUTINE dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs, kpass ) |
|---|
| 45 | !! |
|---|
| 46 | INTEGER , INTENT(in ) :: kt ! ocean time-step index |
|---|
| 47 | INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices |
|---|
| 48 | INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage |
|---|
| 49 | REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pu, pv ! before velocity [m/s] |
|---|
| 50 | REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pu_rhs, pv_rhs ! velocity trend [m/s2] |
|---|
| 51 | !! |
|---|
| 52 | #if defined key_loop_fusion |
|---|
| 53 | CALL dyn_ldf_lap_lf( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs, kpass ) |
|---|
| 54 | #else |
|---|
| 55 | CALL dyn_ldf_lap_t( kt, Kbb, Kmm, pu, pv, is_tile(pu), pu_rhs, pv_rhs, is_tile(pu_rhs), kpass ) |
|---|
| 56 | #endif |
|---|
| 57 | |
|---|
| 58 | END SUBROUTINE dyn_ldf_lap |
|---|
| 59 | |
|---|
| 60 | |
|---|
| 61 | SUBROUTINE dyn_ldf_lap_t( kt, Kbb, Kmm, pu, pv, ktuv, pu_rhs, pv_rhs, ktuv_rhs, kpass ) |
|---|
| 62 | !!---------------------------------------------------------------------- |
|---|
| 63 | !! *** ROUTINE dyn_ldf_lap *** |
|---|
| 64 | !! |
|---|
| 65 | !! ** Purpose : Compute the before horizontal momentum diffusive |
|---|
| 66 | !! trend and add it to the general trend of momentum equation. |
|---|
| 67 | !! |
|---|
| 68 | !! ** Method : The Laplacian operator apply on horizontal velocity is |
|---|
| 69 | !! writen as : grad_h( ahmt div_h(U )) - curl_h( ahmf curl_z(U) ) |
|---|
| 70 | !! |
|---|
| 71 | !! ** Action : - pu_rhs, pv_rhs increased by the harmonic operator applied on pu, pv. |
|---|
| 72 | !! |
|---|
| 73 | !! Reference : S.Griffies, R.Hallberg 2000 Mon.Wea.Rev., DOI:/ |
|---|
| 74 | !!---------------------------------------------------------------------- |
|---|
| 75 | INTEGER , INTENT(in ) :: kt ! ocean time-step index |
|---|
| 76 | INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices |
|---|
| 77 | INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage |
|---|
| 78 | INTEGER , INTENT(in ) :: ktuv, ktuv_rhs |
|---|
| 79 | REAL(wp), DIMENSION(A2D_T(ktuv) ,JPK), INTENT(in ) :: pu, pv ! before velocity [m/s] |
|---|
| 80 | REAL(wp), DIMENSION(A2D_T(ktuv_rhs),JPK), INTENT(inout) :: pu_rhs, pv_rhs ! velocity trend [m/s2] |
|---|
| 81 | ! |
|---|
| 82 | INTEGER :: ji, jj, jk ! dummy loop indices |
|---|
| 83 | INTEGER :: iij |
|---|
| 84 | REAL(wp) :: zsign ! local scalars |
|---|
| 85 | REAL(wp) :: zua, zva ! local scalars |
|---|
| 86 | REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcur, zdiv |
|---|
| 87 | REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zten, zshe ! tension (diagonal) and shearing (anti-diagonal) terms |
|---|
| 88 | !!---------------------------------------------------------------------- |
|---|
| 89 | ! |
|---|
| 90 | IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile |
|---|
| 91 | IF( kt == nit000 .AND. lwp ) THEN |
|---|
| 92 | WRITE(numout,*) |
|---|
| 93 | WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass |
|---|
| 94 | WRITE(numout,*) '~~~~~~~ ' |
|---|
| 95 | ENDIF |
|---|
| 96 | ENDIF |
|---|
| 97 | ! |
|---|
| 98 | ! Define pu_rhs/pv_rhs halo points for multi-point haloes in bilaplacian case |
|---|
| 99 | IF( nldf_dyn == np_blp .AND. kpass == 1 ) THEN ; iij = nn_hls |
|---|
| 100 | ELSE ; iij = 1 |
|---|
| 101 | ENDIF |
|---|
| 102 | ! |
|---|
| 103 | IF( kpass == 1 ) THEN ; zsign = 1._wp ! bilaplacian operator require a minus sign |
|---|
| 104 | ELSE ; zsign = -1._wp ! (eddy viscosity coef. >0) |
|---|
| 105 | ENDIF |
|---|
| 106 | ! |
|---|
| 107 | SELECT CASE( nn_dynldf_typ ) |
|---|
| 108 | ! |
|---|
| 109 | CASE ( np_typ_rot ) !== Vorticity-Divergence operator ==! |
|---|
| 110 | ! |
|---|
| 111 | ALLOCATE( zcur(A2D(nn_hls)) , zdiv(A2D(nn_hls)) ) |
|---|
| 112 | ! |
|---|
| 113 | DO jk = 1, jpkm1 ! Horizontal slab |
|---|
| 114 | ! |
|---|
| 115 | DO_2D( iij-1, iij, iij-1, iij ) |
|---|
| 116 | ! ! ahm * e3 * curl (warning: computed for ji-1,jj-1) |
|---|
| 117 | zcur(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) * e3f(ji-1,jj-1,jk) * r1_e1e2f(ji-1,jj-1) & ! ahmf already * by fmask |
|---|
| 118 | & * ( e2v(ji ,jj-1) * pv(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk) & |
|---|
| 119 | & - e1u(ji-1,jj ) * pu(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk) ) |
|---|
| 120 | ! ! ahm * div (warning: computed for ji,jj) |
|---|
| 121 | zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb) & ! ahmt already * by tmask |
|---|
| 122 | & * ( e2u(ji,jj)*e3u(ji,jj,jk,Kbb) * pu(ji,jj,jk) - e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb) * pu(ji-1,jj,jk) & |
|---|
| 123 | & + e1v(ji,jj)*e3v(ji,jj,jk,Kbb) * pv(ji,jj,jk) - e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kbb) * pv(ji,jj-1,jk) ) |
|---|
| 124 | END_2D |
|---|
| 125 | ! |
|---|
| 126 | DO_2D( iij-1, iij-1, iij-1, iij-1 ) ! - curl( curl) + grad( div ) |
|---|
| 127 | pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * umask(ji,jj,jk) * ( & ! * by umask is mandatory for dyn_ldf_blp use |
|---|
| 128 | & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm) & |
|---|
| 129 | & + ( zdiv(ji+1,jj) - zdiv(ji,jj ) ) * r1_e1u(ji,jj) ) |
|---|
| 130 | ! |
|---|
| 131 | pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * vmask(ji,jj,jk) * ( & ! * by vmask is mandatory for dyn_ldf_blp use |
|---|
| 132 | & ( zcur(ji,jj ) - zcur(ji-1,jj) ) * r1_e1v(ji,jj) / e3v(ji,jj,jk,Kmm) & |
|---|
| 133 | & + ( zdiv(ji,jj+1) - zdiv(ji ,jj) ) * r1_e2v(ji,jj) ) |
|---|
| 134 | END_2D |
|---|
| 135 | ! |
|---|
| 136 | END DO ! End of slab |
|---|
| 137 | ! |
|---|
| 138 | DEALLOCATE( zcur , zdiv ) |
|---|
| 139 | ! |
|---|
| 140 | CASE ( np_typ_sym ) !== Symmetric operator ==! |
|---|
| 141 | ! |
|---|
| 142 | ALLOCATE( zten(A2D(nn_hls)) , zshe(A2D(nn_hls)) ) |
|---|
| 143 | ! |
|---|
| 144 | DO jk = 1, jpkm1 ! Horizontal slab |
|---|
| 145 | ! |
|---|
| 146 | DO_2D( iij-1, iij, iij-1, iij ) |
|---|
| 147 | ! ! shearing stress component (F-point) NB : ahmf has already been multiplied by fmask |
|---|
| 148 | zshe(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) & |
|---|
| 149 | & * ( e1f(ji-1,jj-1) * r1_e2f(ji-1,jj-1) & |
|---|
| 150 | & * ( pu(ji-1,jj ,jk) * r1_e1u(ji-1,jj ) - pu(ji-1,jj-1,jk) * r1_e1u(ji-1,jj-1) ) & |
|---|
| 151 | & + e2f(ji-1,jj-1) * r1_e1f(ji-1,jj-1) & |
|---|
| 152 | & * ( pv(ji ,jj-1,jk) * r1_e2v(ji ,jj-1) - pv(ji-1,jj-1,jk) * r1_e2v(ji-1,jj-1) ) ) |
|---|
| 153 | ! ! tension stress component (T-point) NB : ahmt has already been multiplied by tmask |
|---|
| 154 | zten(ji,jj) = ahmt(ji,jj,jk) & |
|---|
| 155 | & * ( e2t(ji,jj) * r1_e1t(ji,jj) & |
|---|
| 156 | & * ( pu(ji,jj,jk) * r1_e2u(ji,jj) - pu(ji-1,jj,jk) * r1_e2u(ji-1,jj) ) & |
|---|
| 157 | & - e1t(ji,jj) * r1_e2t(ji,jj) & |
|---|
| 158 | & * ( pv(ji,jj,jk) * r1_e1v(ji,jj) - pv(ji,jj-1,jk) * r1_e1v(ji,jj-1) ) ) |
|---|
| 159 | END_2D |
|---|
| 160 | ! |
|---|
| 161 | DO_2D( iij-1, iij-1, iij-1, iij-1 ) |
|---|
| 162 | pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) & |
|---|
| 163 | & * ( ( zten(ji+1,jj ) * e2t(ji+1,jj )*e2t(ji+1,jj ) * e3t(ji+1,jj ,jk,Kmm) & |
|---|
| 164 | & - zten(ji ,jj ) * e2t(ji ,jj )*e2t(ji ,jj ) * e3t(ji ,jj ,jk,Kmm) ) * r1_e2u(ji,jj) & |
|---|
| 165 | & + ( zshe(ji ,jj ) * e1f(ji ,jj )*e1f(ji ,jj ) * e3f(ji ,jj ,jk) & |
|---|
| 166 | & - zshe(ji ,jj-1) * e1f(ji ,jj-1)*e1f(ji ,jj-1) * e3f(ji ,jj-1,jk) ) * r1_e1u(ji,jj) ) |
|---|
| 167 | ! |
|---|
| 168 | pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zsign * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) & |
|---|
| 169 | & * ( ( zshe(ji ,jj ) * e2f(ji ,jj )*e2f(ji ,jj ) * e3f(ji ,jj ,jk) & |
|---|
| 170 | & - zshe(ji-1,jj ) * e2f(ji-1,jj )*e2f(ji-1,jj ) * e3f(ji-1,jj ,jk) ) * r1_e2v(ji,jj) & |
|---|
| 171 | & - ( zten(ji ,jj+1) * e1t(ji ,jj+1)*e1t(ji ,jj+1) * e3t(ji ,jj+1,jk,Kmm) & |
|---|
| 172 | & - zten(ji ,jj ) * e1t(ji ,jj )*e1t(ji ,jj ) * e3t(ji ,jj ,jk,Kmm) ) * r1_e1v(ji,jj) ) |
|---|
| 173 | ! |
|---|
| 174 | END_2D |
|---|
| 175 | ! |
|---|
| 176 | END DO |
|---|
| 177 | ! |
|---|
| 178 | DEALLOCATE( zten , zshe ) |
|---|
| 179 | ! |
|---|
| 180 | END SELECT |
|---|
| 181 | ! |
|---|
| 182 | END SUBROUTINE dyn_ldf_lap_t |
|---|
| 183 | |
|---|
| 184 | |
|---|
| 185 | SUBROUTINE dyn_ldf_blp( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs ) |
|---|
| 186 | !!---------------------------------------------------------------------- |
|---|
| 187 | !! *** ROUTINE dyn_ldf_blp *** |
|---|
| 188 | !! |
|---|
| 189 | !! ** Purpose : Compute the before lateral momentum viscous trend |
|---|
| 190 | !! and add it to the general trend of momentum equation. |
|---|
| 191 | !! |
|---|
| 192 | !! ** Method : The lateral viscous trends is provided by a bilaplacian |
|---|
| 193 | !! operator applied to before field (forward in time). |
|---|
| 194 | !! It is computed by two successive calls to dyn_ldf_lap routine |
|---|
| 195 | !! |
|---|
| 196 | !! ** Action : pt(:,:,:,:,Krhs) updated with the before rotated bilaplacian diffusion |
|---|
| 197 | !!---------------------------------------------------------------------- |
|---|
| 198 | INTEGER , INTENT(in ) :: kt ! ocean time-step index |
|---|
| 199 | INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices |
|---|
| 200 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu, pv ! before velocity fields |
|---|
| 201 | REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! momentum trend |
|---|
| 202 | ! |
|---|
| 203 | REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zulap, zvlap ! laplacian at u- and v-point |
|---|
| 204 | !!---------------------------------------------------------------------- |
|---|
| 205 | ! |
|---|
| 206 | #if defined key_loop_fusion |
|---|
| 207 | CALL dyn_ldf_blp_lf( kt, Kbb, Kmm, pu, pv, pu_rhs, pv_rhs ) |
|---|
| 208 | #else |
|---|
| 209 | IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile |
|---|
| 210 | IF( kt == nit000 ) THEN |
|---|
| 211 | IF(lwp) WRITE(numout,*) |
|---|
| 212 | IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum ' |
|---|
| 213 | IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' |
|---|
| 214 | ENDIF |
|---|
| 215 | ENDIF |
|---|
| 216 | ! |
|---|
| 217 | zulap(:,:,:) = 0._wp |
|---|
| 218 | zvlap(:,:,:) = 0._wp |
|---|
| 219 | ! |
|---|
| 220 | CALL dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, zulap, zvlap, 1 ) ! rotated laplacian applied to pt (output in zlap,Kbb) |
|---|
| 221 | ! |
|---|
| 222 | IF (nn_hls==1) CALL lbc_lnk( 'dynldf_lap_blp', zulap, 'U', -1.0_wp, zvlap, 'V', -1.0_wp ) ! Lateral boundary conditions |
|---|
| 223 | ! |
|---|
| 224 | CALL dyn_ldf_lap( kt, Kbb, Kmm, zulap, zvlap, pu_rhs, pv_rhs, 2 ) ! rotated laplacian applied to zlap (output in pt(:,:,:,:,Krhs)) |
|---|
| 225 | ! |
|---|
| 226 | #endif |
|---|
| 227 | END SUBROUTINE dyn_ldf_blp |
|---|
| 228 | |
|---|
| 229 | !!====================================================================== |
|---|
| 230 | END MODULE dynldf_lap_blp |
|---|
