- Timestamp:
- 2020-06-24T15:31:32+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynkeg.F90
r13151 r13154 29 29 30 30 PUBLIC dyn_keg ! routine called by step module 31 PUBLIC dyn_kegAD ! routine called by step module 31 32 32 33 INTEGER, PARAMETER, PUBLIC :: nkeg_C2 = 0 !: 2nd order centered scheme (standard scheme) … … 155 156 ! 156 157 END SUBROUTINE dyn_keg 157 158 159 160 SUBROUTINE dyn_kegAD( kt, kscheme, puu, pvv, pu_rhs, pv_rhs ) 161 !!---------------------------------------------------------------------- 162 !! *** ROUTINE dyn_kegAD *** 163 !! 164 !! ** Purpose : Compute the now momentum trend due to the horizontal 165 !! gradient of the horizontal kinetic energy and add it to the 166 !! general momentum trend. 167 !! 168 !! ** Method : * kscheme = nkeg_C2 : 2nd order centered scheme that 169 !! conserve kinetic energy. Compute the now horizontal kinetic energy 170 !! zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ] 171 !! * kscheme = nkeg_HW : Hollingsworth correction following 172 !! Arakawa (2001). The now horizontal kinetic energy is given by: 173 !! zhke = 1/6 [ mi-1( 2 * un^2 + ((u(j+1)+u(j-1))/2)^2 ) 174 !! + mj-1( 2 * vn^2 + ((v(i+1)+v(i-1))/2)^2 ) ] 175 !! 176 !! Take its horizontal gradient and add it to the general momentum 177 !! trend. 178 !! u(rhs) = u(rhs) - 1/e1u di[ zhke ] 179 !! v(rhs) = v(rhs) - 1/e2v dj[ zhke ] 180 !! 181 !! ** Action : - Update the (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) with the hor. ke gradient trend 182 !! - send this trends to trd_dyn (l_trddyn=T) for post-processing 183 !! 184 !! ** References : Arakawa, A., International Geophysics 2001. 185 !! Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983. 186 !!---------------------------------------------------------------------- 187 ! 188 INTEGER , INTENT( in ) :: kt ! ocean time-step index 189 INTEGER , INTENT( in ) :: kscheme ! =0/1/2 type of KEG scheme 190 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: puu, pvv ! ocean velocities at Kmm 191 REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) :: pu_rhs, pv_rhs ! RHS 192 ! 193 INTEGER :: ji, jj, jk ! dummy loop indices 194 REAL(wp) :: zu, zv ! local scalars 195 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhke 196 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv 197 !!---------------------------------------------------------------------- 198 ! 199 IF( ln_timing ) CALL timing_start('dyn_keg') 200 ! 201 IF( kt == nit000 ) THEN 202 IF(lwp) WRITE(numout,*) 203 IF(lwp) WRITE(numout,*) 'dyn_kegAD : kinetic energy gradient trend, scheme number=', kscheme 204 IF(lwp) WRITE(numout,*) '~~~~~~~~~' 205 ENDIF 206 207 zhke(:,:,jpk) = 0._wp 208 209 SELECT CASE ( kscheme ) !== Horizontal kinetic energy at T-point ==! 210 !!an45 to be ADDED : que cas C2 - "wet points only" il suffit de x2 le terme quadratic a la coast (nn_dynkeg_adv = 2) 211 CASE ( nkeg_C2_wpo ) !-- Standard scheme --! 212 DO_3D_01_01( 1, jpkm1 ) 213 zu = ( puu(ji-1,jj ,jk) * puu(ji-1,jj ,jk) & 214 & + puu(ji ,jj ,jk) * puu(ji ,jj ,jk) ) * ( 2._wp - umask(ji-1,jj,jk) * umask(ji,jj,jk) ) 215 zv = ( pvv(ji ,jj-1,jk) * pvv(ji ,jj-1,jk) & 216 & + pvv(ji ,jj ,jk) * pvv(ji ,jj ,jk) ) * ( 2._wp - vmask(ji,jj-1,jk) * vmask(ji,jj,jk) ) 217 zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 218 END_3D 219 !!an45 220 ! 221 CASE ( nkeg_C2 ) !-- Standard scheme --! 222 DO_3D_01_01( 1, jpkm1 ) 223 zu = puu(ji-1,jj ,jk) * puu(ji-1,jj ,jk) & 224 & + puu(ji ,jj ,jk) * puu(ji ,jj ,jk) 225 zv = pvv(ji ,jj-1,jk) * pvv(ji ,jj-1,jk) & 226 & + pvv(ji ,jj ,jk) * pvv(ji ,jj ,jk) 227 zhke(ji,jj,jk) = 0.25_wp * ( zv + zu ) 228 END_3D 229 !!an 00_00 ? 230 CASE ( nkeg_HW ) !-- Hollingsworth scheme --! 231 DO_3D_00_00( 1, jpkm1 ) 232 zu = 8._wp * ( puu(ji-1,jj ,jk) * puu(ji-1,jj ,jk) & 233 & + puu(ji ,jj ,jk) * puu(ji ,jj ,jk) ) & 234 & + ( puu(ji-1,jj-1,jk) + puu(ji-1,jj+1,jk) ) * ( puu(ji-1,jj-1,jk) + puu(ji-1,jj+1,jk) ) & 235 & + ( puu(ji ,jj-1,jk) + puu(ji ,jj+1,jk) ) * ( puu(ji ,jj-1,jk) + puu(ji ,jj+1,jk) ) 236 ! 237 zv = 8._wp * ( pvv(ji ,jj-1,jk) * pvv(ji ,jj-1,jk) & 238 & + pvv(ji ,jj ,jk) * pvv(ji ,jj ,jk) ) & 239 & + ( pvv(ji-1,jj-1,jk) + pvv(ji+1,jj-1,jk) ) * ( pvv(ji-1,jj-1,jk) + pvv(ji+1,jj-1,jk) ) & 240 & + ( pvv(ji-1,jj ,jk) + pvv(ji+1,jj ,jk) ) * ( pvv(ji-1,jj ,jk) + pvv(ji+1,jj ,jk) ) 241 zhke(ji,jj,jk) = r1_48 * ( zv + zu ) 242 END_3D 243 CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. ) 244 ! 245 END SELECT 246 ! 247 IF( PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** NO alternating direction ***! 248 ! 249 DO_3D_00_00( 1, jpkm1 ) 250 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 251 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) 252 END_3D 253 ! 254 ELSEIF( PRESENT( pu_rhs ) .AND. .NOT. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : i-component ***! 255 ! 256 DO_3D_00_00( 1, jpkm1 ) 257 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - ( zhke(ji+1,jj ,jk) - zhke(ji,jj,jk) ) / e1u(ji,jj) 258 END_3D 259 ! 260 ELSEIF( .NOT. PRESENT( pu_rhs ) .AND. PRESENT( pv_rhs ) ) THEN !*** Alternating direction : j-component ***! 261 ! 262 DO_3D_00_00( 1, jpkm1 ) 263 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - ( zhke(ji ,jj+1,jk) - zhke(ji,jj,jk) ) / e2v(ji,jj) 264 END_3D 265 ! 266 ENDIF 267 IF( ln_timing ) CALL timing_stop('dyn_kegAD') 268 ! 269 END SUBROUTINE dyn_kegAD 158 270 !!====================================================================== 159 271 END MODULE dynkeg
Note: See TracChangeset
for help on using the changeset viewer.
