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