Changeset 15548 for NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/DYN/dynldf_lap_blp.F90
- Timestamp:
- 2021-11-28T18:59:49+01:00 (3 years ago)
- Location:
- NEMO/branches/2021/ticket2632_r14588_theta_sbcblk
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/ticket2632_r14588_theta_sbcblk
- Property svn:externals
-
old new 9 9 10 10 # SETTE 11 ^/utils/CI/sette@14244 sette 11 ^/utils/CI/sette@HEAD sette 12
-
- Property svn:externals
-
NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/DYN/dynldf_lap_blp.F90
r14433 r15548 14 14 USE oce ! ocean dynamics and tracers 15 15 USE dom_oce ! ocean space and time domain 16 USE domutl, ONLY : is_tile 16 17 USE ldfdyn ! lateral diffusion: eddy viscosity coef. 17 18 USE ldfslp ! iso-neutral slopes … … 21 22 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 22 23 USE lib_mpp 23 24 #if defined key_loop_fusion 25 USE dynldf_lap_blp_lf 26 #endif 27 24 28 IMPLICIT NONE 25 29 PRIVATE … … 39 43 40 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 ) 41 62 !!---------------------------------------------------------------------- 42 63 !! *** ROUTINE dyn_ldf_lap *** … … 52 73 !! Reference : S.Griffies, R.Hallberg 2000 Mon.Wea.Rev., DOI:/ 53 74 !!---------------------------------------------------------------------- 54 INTEGER , INTENT(in ) :: kt ! ocean time-step index 55 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 56 INTEGER , INTENT(in ) :: kpass ! =1/2 first or second passage 57 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu, pv ! before velocity [m/s] 58 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! velocity trend [m/s2] 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] 59 81 ! 60 82 INTEGER :: ji, jj, jk ! dummy loop indices 83 INTEGER :: iij 61 84 REAL(wp) :: zsign ! local scalars 62 85 REAL(wp) :: zua, zva ! local scalars … … 65 88 !!---------------------------------------------------------------------- 66 89 ! 67 IF( kt == nit000 .AND. lwp ) THEN 68 WRITE(numout,*) 69 WRITE(numout,*) 'dyn_ldf : iso-level harmonic (laplacian) operator, pass=', kpass 70 WRITE(numout,*) '~~~~~~~ ' 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 71 101 ENDIF 72 102 ! … … 79 109 CASE ( np_typ_rot ) !== Vorticity-Divergence operator ==! 80 110 ! 81 ALLOCATE( zcur( jpi,jpj) , zdiv(jpi,jpj) )111 ALLOCATE( zcur(A2D(nn_hls)) , zdiv(A2D(nn_hls)) ) 82 112 ! 83 113 DO jk = 1, jpkm1 ! Horizontal slab 84 114 ! 85 DO_2D( 0, 1, 0, 1)86 ! ! ahm * e3 * curl ( computed from 1 to jpim1/jpjm1)115 DO_2D( iij-1, iij, iij-1, iij ) 116 ! ! ahm * e3 * curl (warning: computed for ji-1,jj-1) 87 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 88 118 & * ( e2v(ji ,jj-1) * pv(ji ,jj-1,jk) - e2v(ji-1,jj-1) * pv(ji-1,jj-1,jk) & 89 119 & - e1u(ji-1,jj ) * pu(ji-1,jj ,jk) + e1u(ji-1,jj-1) * pu(ji-1,jj-1,jk) ) 90 ! ! ahm * div ( computed from 2 to jpi/jpj)120 ! ! ahm * div (warning: computed for ji,jj) 91 121 zdiv(ji,jj) = ahmt(ji,jj,jk) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb) & ! ahmt already * by tmask 92 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) & … … 94 124 END_2D 95 125 ! 96 DO_2D( 0, 0, 0, 0 )! - curl( curl) + grad( div )126 DO_2D( iij-1, iij-1, iij-1, iij-1 ) ! - curl( curl) + grad( div ) 97 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 98 128 & - ( zcur(ji ,jj) - zcur(ji,jj-1) ) * r1_e2u(ji,jj) / e3u(ji,jj,jk,Kmm) & … … 110 140 CASE ( np_typ_sym ) !== Symmetric operator ==! 111 141 ! 112 ALLOCATE( zten( jpi,jpj) , zshe(jpi,jpj) )142 ALLOCATE( zten(A2D(nn_hls)) , zshe(A2D(nn_hls)) ) 113 143 ! 114 144 DO jk = 1, jpkm1 ! Horizontal slab 115 145 ! 116 DO_2D( 0, 1, 0, 1)146 DO_2D( iij-1, iij, iij-1, iij ) 117 147 ! ! shearing stress component (F-point) NB : ahmf has already been multiplied by fmask 118 148 zshe(ji-1,jj-1) = ahmf(ji-1,jj-1,jk) & … … 129 159 END_2D 130 160 ! 131 DO_2D( 0, 0, 0, 0)161 DO_2D( iij-1, iij-1, iij-1, iij-1 ) 132 162 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zsign * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) & 133 163 & * ( ( zten(ji+1,jj ) * e2t(ji+1,jj )*e2t(ji+1,jj ) * e3t(ji+1,jj ,jk,Kmm) & … … 150 180 END SELECT 151 181 ! 152 END SUBROUTINE dyn_ldf_lap 182 END SUBROUTINE dyn_ldf_lap_t 153 183 154 184 … … 171 201 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! momentum trend 172 202 ! 173 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zulap, zvlap ! laplacian at u- and v-point 174 !!---------------------------------------------------------------------- 175 ! 176 IF( kt == nit000 ) THEN 177 IF(lwp) WRITE(numout,*) 178 IF(lwp) WRITE(numout,*) 'dyn_ldf_blp : bilaplacian operator momentum ' 179 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 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 180 215 ENDIF 181 216 ! … … 185 220 CALL dyn_ldf_lap( kt, Kbb, Kmm, pu, pv, zulap, zvlap, 1 ) ! rotated laplacian applied to pt (output in zlap,Kbb) 186 221 ! 187 CALL lbc_lnk( 'dynldf_lap_blp', zulap, 'U', -1.0_wp, zvlap, 'V', -1.0_wp ) ! Lateral boundary conditions222 IF (nn_hls==1) CALL lbc_lnk( 'dynldf_lap_blp', zulap, 'U', -1.0_wp, zvlap, 'V', -1.0_wp ) ! Lateral boundary conditions 188 223 ! 189 224 CALL dyn_ldf_lap( kt, Kbb, Kmm, zulap, zvlap, pu_rhs, pv_rhs, 2 ) ! rotated laplacian applied to zlap (output in pt(:,:,:,:,Krhs)) 190 225 ! 226 #endif 191 227 END SUBROUTINE dyn_ldf_blp 192 228
Note: See TracChangeset
for help on using the changeset viewer.