- Timestamp:
- 2020-06-24T14:38:26+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12527_Gurvan_ShallowWater/src/SWE/dynkeg.F90
r13005 r13151 29 29 30 30 PUBLIC dyn_keg ! routine called by step module 31 PUBLIC dyn_kegAD ! routine called by step module32 31 33 32 INTEGER, PARAMETER, PUBLIC :: nkeg_C2 = 0 !: 2nd order centered scheme (standard scheme) … … 156 155 ! 157 156 END SUBROUTINE dyn_keg 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 horizontal165 !! gradient of the horizontal kinetic energy and add it to the166 !! general momentum trend.167 !!168 !! ** Method : * kscheme = nkeg_C2 : 2nd order centered scheme that169 !! conserve kinetic energy. Compute the now horizontal kinetic energy170 !! zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]171 !! * kscheme = nkeg_HW : Hollingsworth correction following172 !! 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 momentum177 !! 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 trend182 !! - send this trends to trd_dyn (l_trddyn=T) for post-processing183 !!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 index189 INTEGER , INTENT( in ) :: kscheme ! =0/1/2 type of KEG scheme190 REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: puu, pvv ! ocean velocities at Kmm191 REAL(wp), DIMENSION(jpi,jpj,jpk),OPTIONAL, INTENT(inout) :: pu_rhs, pv_rhs ! RHS192 !193 INTEGER :: ji, jj, jk ! dummy loop indices194 REAL(wp) :: zu, zv ! local scalars195 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhke196 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv197 !!----------------------------------------------------------------------198 !199 IF( ln_timing ) CALL timing_start('dyn_keg')200 !201 IF( kt == nit000 ) THEN202 IF(lwp) WRITE(numout,*)203 IF(lwp) WRITE(numout,*) 'dyn_kegAD : kinetic energy gradient trend, scheme number=', kscheme204 IF(lwp) WRITE(numout,*) '~~~~~~~~~'205 ENDIF206 207 zhke(:,:,jpk) = 0._wp208 157 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_3D219 !!an45220 !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_3D229 !!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_3D243 CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. )244 !245 END SELECT246 !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_3D253 !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_3D259 !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_3D265 !266 ENDIF267 IF( ln_timing ) CALL timing_stop('dyn_kegAD')268 !269 END SUBROUTINE dyn_kegAD270 158 !!====================================================================== 271 159 END MODULE dynkeg
Note: See TracChangeset
for help on using the changeset viewer.
