Changeset 14789 for NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/ICB/icbdyn.F90
- Timestamp:
- 2021-05-05T13:18:04+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev _r12970_AGRIF_CMEMSext/AGRIF5 ^/vendors/AGRIF/dev@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 ^/vendors/PPR@HEAD ext/PPR 8 9 9 10 # SETTE 10 ^/utils/CI/sette@1 3559sette11 ^/utils/CI/sette@14244 sette
-
- Property svn:externals
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/ICB/icbdyn.F90
r13281 r14789 14 14 USE dom_oce ! NEMO ocean domain 15 15 USE phycst ! NEMO physical constants 16 USE in_out_manager ! IO parameters 16 17 ! 17 18 USE icb_oce ! define iceberg arrays … … 84 85 85 86 ! !** A1 = A(X1,V1) 86 CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1, &87 & zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 )87 CALL icb_accel( kt, berg , zxi1, ze1, zuvel1, zuvel1, zax1, & 88 & zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2, 0.5_wp ) 88 89 ! 89 90 zu1 = zuvel1 / ze1 !** V1 in d(i,j)/dt … … 97 98 zyj2 = zyj1 + zdt_2 * zv1 ; zvvel2 = zvvel1 + zdt_2 * zay1 98 99 ! 99 CALL icb_ground( zxi2, zxi1, zu1, &100 & zyj2, zyj1, zv1, ll_bounced )100 CALL icb_ground( berg, zxi2, zxi1, zu1, & 101 & zyj2, zyj1, zv1, ll_bounced ) 101 102 102 103 ! !** A2 = A(X2,V2) 103 CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2, &104 & zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 )104 CALL icb_accel( kt, berg , zxi2, ze1, zuvel2, zuvel1, zax2, & 105 & zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2, 0.5_wp ) 105 106 ! 106 107 zu2 = zuvel2 / ze1 !** V2 in d(i,j)/dt … … 113 114 zyj3 = zyj1 + zdt_2 * zv2 ; zvvel3 = zvvel1 + zdt_2 * zay2 114 115 ! 115 CALL icb_ground( zxi3, zxi1, zu3, &116 & zyj3, zyj1, zv3, ll_bounced )116 CALL icb_ground( berg, zxi3, zxi1, zu2, & 117 & zyj3, zyj1, zv2, ll_bounced ) 117 118 118 119 ! !** A3 = A(X3,V3) 119 CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3, &120 & zyj3, ze2, zvvel3, zvvel1, zay3, zdt )120 CALL icb_accel( kt, berg , zxi3, ze1, zuvel3, zuvel1, zax3, & 121 & zyj3, ze2, zvvel3, zvvel1, zay3, zdt, 1._wp ) 121 122 ! 122 123 zu3 = zuvel3 / ze1 !** V3 in d(i,j)/dt … … 129 130 zyj4 = zyj1 + zdt * zv3 ; zvvel4 = zvvel1 + zdt * zay3 130 131 131 CALL icb_ground( zxi4, zxi1, zu4, &132 & zyj4, zyj1, zv4, ll_bounced )132 CALL icb_ground( berg, zxi4, zxi1, zu3, & 133 & zyj4, zyj1, zv3, ll_bounced ) 133 134 134 135 ! !** A4 = A(X4,V4) 135 CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4, &136 & zyj4, ze2, zvvel4, zvvel1, zay4, zdt )136 CALL icb_accel( kt, berg , zxi4, ze1, zuvel4, zuvel1, zax4, & 137 & zyj4, ze2, zvvel4, zvvel1, zay4, zdt, 1._wp ) 137 138 138 139 zu4 = zuvel4 / ze1 !** V4 in d(i,j)/dt … … 148 149 zvvel_n = pt%vvel + zdt_6 * ( zay1 + 2.*(zay2 + zay3) + zay4 ) 149 150 150 CALL icb_ground( zxi_n, zxi1, zuvel_n, &151 & zyj_n, zyj1, zvvel_n, ll_bounced )151 CALL icb_ground( berg, zxi_n, zxi1, zuvel_n, & 152 & zyj_n, zyj1, zvvel_n, ll_bounced ) 152 153 153 154 pt%uvel = zuvel_n !** save in berg structure … … 156 157 pt%yj = zyj_n 157 158 158 ! update actual position159 pt%lon = icb_utl_bilin_x(glamt, pt%xi, pt%yj )160 pt%lat = icb_utl_bilin(gphit, pt%xi, pt%yj, 'T' )161 162 159 berg => berg%next ! switch to the next berg 163 160 ! … … 167 164 168 165 169 SUBROUTINE icb_ground( pi, pi0, pu, &170 & pj, pj0, pv, ld_bounced )166 SUBROUTINE icb_ground( berg, pi, pi0, pu, & 167 & pj, pj0, pv, ld_bounced ) 171 168 !!---------------------------------------------------------------------- 172 169 !! *** ROUTINE icb_ground *** … … 177 174 !! NB two possibilities available one of which is hard-coded here 178 175 !!---------------------------------------------------------------------- 176 TYPE(iceberg ), POINTER, INTENT(in ) :: berg ! berg 177 ! 179 178 REAL(wp), INTENT(inout) :: pi , pj ! current iceberg position 180 179 REAL(wp), INTENT(in ) :: pi0, pj0 ! previous iceberg position … … 184 183 INTEGER :: ii, ii0 185 184 INTEGER :: ij, ij0 185 INTEGER :: ikb 186 186 INTEGER :: ibounce_method 187 ! 188 REAL(wp) :: zD 189 REAL(wp), DIMENSION(jpk) :: ze3t 187 190 !!---------------------------------------------------------------------- 188 191 ! … … 200 203 ij = mj1( ij ) 201 204 ! 202 IF( tmask(ii,ij,1) /= 0._wp ) RETURN ! berg reach a new t-cell, but an ocean one 205 ! assume icb is grounded if tmask(ii,ij,1) or tmask(ii,ij,ikb), depending of the option is not 0 206 IF ( ln_M2016 .AND. ln_icb_grd ) THEN 207 ! 208 ! draught (keel depth) 209 zD = rho_berg_1_oce * berg%current_point%thickness 210 ! 211 ! interpol needed data 212 CALL icb_utl_interp( pi, pj, pe3t=ze3t ) 213 ! 214 !compute bottom level 215 CALL icb_utl_getkb( ikb, ze3t, zD ) 216 ! 217 ! berg reach a new t-cell, but an ocean one 218 ! .AND. needed in case berg hit an isf (tmask(ii,ij,1) == 0 and tmask(ii,ij,ikb) /= 0) 219 IF( tmask(ii,ij,ikb) /= 0._wp .AND. tmask(ii,ij,1) /= 0._wp ) RETURN 220 ! 221 ELSE 222 IF( tmask(ii,ij,1) /= 0._wp ) RETURN ! berg reach a new t-cell, but an ocean one 223 END IF 203 224 ! 204 225 ! From here, berg have reach land: treat grounding/bouncing … … 234 255 235 256 236 SUBROUTINE icb_accel( berg , pxi, pe1, puvel, puvel0, pax,&237 & pyj, pe2, pvvel, pvvel0, pay, pdt)257 SUBROUTINE icb_accel( kt, berg , pxi, pe1, puvel, puvel0, pax, & 258 & pyj, pe2, pvvel, pvvel0, pay, pdt, pcfl_scale ) 238 259 !!---------------------------------------------------------------------- 239 260 !! *** ROUTINE icb_accel *** … … 244 265 !!---------------------------------------------------------------------- 245 266 TYPE(iceberg ), POINTER, INTENT(in ) :: berg ! berg 267 INTEGER , INTENT(in ) :: kt ! time step 268 REAL(wp) , INTENT(in ) :: pcfl_scale 246 269 REAL(wp) , INTENT(in ) :: pxi , pyj ! berg position in (i,j) referential 247 270 REAL(wp) , INTENT(in ) :: puvel , pvvel ! berg velocity [m/s] … … 257 280 REAL(wp), PARAMETER :: pp_Cr0 = 0.06_wp ! 258 281 ! 259 INTEGER :: itloop 260 REAL(wp) :: zuo, z ui, zua, zuwave, zssh_x, zsst, zcn, zhi, zsss261 REAL(wp) :: zvo, z vi, zva, zvwave, zssh_y282 INTEGER :: itloop, ikb, jk 283 REAL(wp) :: zuo, zssu, zui, zua, zuwave, zssh_x, zcn, zhi 284 REAL(wp) :: zvo, zssv, zvi, zva, zvwave, zssh_y 262 285 REAL(wp) :: zff, zT, zD, zW, zL, zM, zF 263 286 REAL(wp) :: zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad 264 REAL(wp) :: z_ocn, z_atm, z_ice 287 REAL(wp) :: z_ocn, z_atm, z_ice, zdep 265 288 REAL(wp) :: zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop 266 289 REAL(wp) :: zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi 267 290 REAL(wp) :: zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new 291 REAL(wp), DIMENSION(jpk) :: zuoce, zvoce, ze3t, zdepw 268 292 !!---------------------------------------------------------------------- 269 293 270 294 ! Interpolate gridded fields to berg 271 295 nknberg = berg%number(1) 272 CALL icb_utl_interp( pxi, pe1, zuo, zui, zua, zssh_x, & 273 & pyj, pe2, zvo, zvi, zva, zssh_y, zsst, zcn, zhi, zff, zsss ) 296 CALL icb_utl_interp( pxi, pyj, pe1=pe1, pe2=pe2, & ! scale factor 297 & pssu=zssu, pui=zui, pua=zua, & ! oce/ice/atm velocities 298 & pssv=zssv, pvi=zvi, pva=zva, & ! oce/ice/atm velocities 299 & pssh_i=zssh_x, pssh_j=zssh_y, & ! ssh gradient 300 & phi=zhi, pff=zff) ! ice thickness and coriolis 274 301 275 302 zM = berg%current_point%mass 276 303 zT = berg%current_point%thickness ! total thickness 277 zD = ( rn_rho_bergs / pp_rho_seawater ) * zT! draught (keel depth)304 zD = rho_berg_1_oce * zT ! draught (keel depth) 278 305 zF = zT - zD ! freeboard 279 306 zW = berg%current_point%width … … 282 309 zhi = MIN( zhi , zD ) 283 310 zD_hi = MAX( 0._wp, zD-zhi ) 284 285 286 zuwave = zua - z uo ; zvwave = zva - zvo! Use wind speed rel. to ocean for wave model311 312 ! Wave radiation 313 zuwave = zua - zssu ; zvwave = zva - zssv ! Use wind speed rel. to ocean for wave model 287 314 zwmod = zuwave*zuwave + zvwave*zvwave ! The wave amplitude and length depend on the current; 288 315 ! ! wind speed relative to the ocean. Actually wmod is wmod**2 here. … … 309 336 IF( abs(zui) + abs(zvi) == 0._wp ) z_ice = 0._wp 310 337 338 ! lateral velocities 339 ! default ssu and ssv 340 ! ln_M2016: mean velocity along the profile 341 IF ( ln_M2016 ) THEN 342 ! interpol needed data 343 CALL icb_utl_interp( pxi, pyj, puoce=zuoce, pvoce=zvoce, pe3t=ze3t ) ! 3d velocities 344 345 !compute bottom level 346 CALL icb_utl_getkb( ikb, ze3t, zD ) 347 348 ! compute mean velocity 349 CALL icb_utl_zavg(zuo, zuoce, ze3t, zD, ikb) 350 CALL icb_utl_zavg(zvo, zvoce, ze3t, zD, ikb) 351 ELSE 352 zuo = zssu 353 zvo = zssv 354 END IF 355 311 356 zuveln = puvel ; zvveln = pvvel ! Copy starting uvel, vvel 312 357 ! … … 321 366 ! Explicit accelerations 322 367 !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave & 323 ! -zdrag_ocn*(puvel-z uo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui)368 ! -zdrag_ocn*(puvel-zssu) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui) 324 369 !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave & 325 ! -zdrag_ocn*(pvvel-z vo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi)370 ! -zdrag_ocn*(pvvel-zssv) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi) 326 371 zaxe = -grav * zssh_x + zwave_rad * zuwave 327 372 zaye = -grav * zssh_y + zwave_rad * zvwave … … 361 406 zspeed = SQRT( zuveln*zuveln + zvveln*zvveln ) ! Speed of berg 362 407 IF( zspeed > 0._wp ) THEN 363 zloc_dx = MIN( pe1, pe2 ) ! minimum grid spacing 364 zspeed_new = zloc_dx / pdt * rn_speed_limit ! Speed limit as a factor of dx / dt 408 zloc_dx = MIN( pe1, pe2 ) ! minimum grid spacing 409 ! cfl scale is function of the RK4 step 410 zspeed_new = zloc_dx / pdt * rn_speed_limit * pcfl_scale ! Speed limit as a factor of dx / dt 365 411 IF( zspeed_new < zspeed ) THEN 366 zuveln = zuveln * ( zspeed_new / zspeed ) ! Scale velocity to reduce speed 367 zvveln = zvveln * ( zspeed_new / zspeed ) ! without changing the direction 412 zuveln = zuveln * ( zspeed_new / zspeed ) ! Scale velocity to reduce speed 413 zvveln = zvveln * ( zspeed_new / zspeed ) ! without changing the direction 414 pax = (zuveln - puvel0)/pdt 415 pay = (zvveln - pvvel0)/pdt 416 ! 417 ! print speeding ticket 418 IF (nn_verbose_level > 0) THEN 419 WRITE(numicb, 9200) 'icb speeding : ',kt, nknberg, zspeed, & 420 & pxi, pyj, zuo, zvo, zua, zva, zui, zvi 421 9200 FORMAT(a,i9,i6,f9.2,1x,4(1x,2f9.2)) 422 END IF 423 ! 368 424 CALL icb_dia_speed() 369 425 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.