Changeset 14030 for NEMO/trunk/src/OCE/ICB/icbdyn.F90
- Timestamp:
- 2020-12-03T10:26:33+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/OCE/ICB/icbdyn.F90
r13281 r14030 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 … … 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) … … 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, zu3, & 117 & zyj3, zyj1, zv3, ll_bounced ) 117 118 118 119 ! !** A3 = A(X3,V3) … … 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, zu4, & 133 & zyj4, zyj1, zv4, ll_bounced ) 133 134 134 135 ! !** A4 = A(X4,V4) … … 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 … … 257 278 REAL(wp), PARAMETER :: pp_Cr0 = 0.06_wp ! 258 279 ! 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_y280 INTEGER :: itloop, ikb, jk 281 REAL(wp) :: zuo, zssu, zui, zua, zuwave, zssh_x, zcn, zhi 282 REAL(wp) :: zvo, zssv, zvi, zva, zvwave, zssh_y 262 283 REAL(wp) :: zff, zT, zD, zW, zL, zM, zF 263 284 REAL(wp) :: zdrag_ocn, zdrag_atm, zdrag_ice, zwave_rad 264 REAL(wp) :: z_ocn, z_atm, z_ice 285 REAL(wp) :: z_ocn, z_atm, z_ice, zdep 265 286 REAL(wp) :: zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop 266 287 REAL(wp) :: zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi 267 288 REAL(wp) :: zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new 289 REAL(wp), DIMENSION(jpk) :: zuoce, zvoce, ze3t, zdepw 268 290 !!---------------------------------------------------------------------- 269 291 270 292 ! Interpolate gridded fields to berg 271 293 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 ) 294 CALL icb_utl_interp( pxi, pyj, pe1=pe1, pe2=pe2, & ! scale factor 295 & pssu=zssu, pui=zui, pua=zua, & ! oce/ice/atm velocities 296 & pssv=zssv, pvi=zvi, pva=zva, & ! oce/ice/atm velocities 297 & pssh_i=zssh_x, pssh_j=zssh_y, & ! ssh gradient 298 & phi=zhi, pff=zff) ! ice thickness and coriolis 274 299 275 300 zM = berg%current_point%mass 276 301 zT = berg%current_point%thickness ! total thickness 277 zD = ( rn_rho_bergs / pp_rho_seawater ) * zT! draught (keel depth)302 zD = rho_berg_1_oce * zT ! draught (keel depth) 278 303 zF = zT - zD ! freeboard 279 304 zW = berg%current_point%width … … 282 307 zhi = MIN( zhi , zD ) 283 308 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 model309 310 ! Wave radiation 311 zuwave = zua - zssu ; zvwave = zva - zssv ! Use wind speed rel. to ocean for wave model 287 312 zwmod = zuwave*zuwave + zvwave*zvwave ! The wave amplitude and length depend on the current; 288 313 ! ! wind speed relative to the ocean. Actually wmod is wmod**2 here. … … 309 334 IF( abs(zui) + abs(zvi) == 0._wp ) z_ice = 0._wp 310 335 336 ! lateral velocities 337 ! default ssu and ssv 338 ! ln_M2016: mean velocity along the profile 339 IF ( ln_M2016 ) THEN 340 ! interpol needed data 341 CALL icb_utl_interp( pxi, pyj, puoce=zuoce, pvoce=zvoce, pe3t=ze3t ) ! 3d velocities 342 343 !compute bottom level 344 CALL icb_utl_getkb( ikb, ze3t, zD ) 345 346 ! compute mean velocity 347 CALL icb_utl_zavg(zuo, zuoce, ze3t, zD, ikb) 348 CALL icb_utl_zavg(zvo, zvoce, ze3t, zD, ikb) 349 ELSE 350 zuo = zssu 351 zvo = zssv 352 END IF 353 311 354 zuveln = puvel ; zvveln = pvvel ! Copy starting uvel, vvel 312 355 ! … … 321 364 ! Explicit accelerations 322 365 !zaxe= zff*pvvel -grav*zssh_x +zwave_rad*zuwave & 323 ! -zdrag_ocn*(puvel-z uo) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui)366 ! -zdrag_ocn*(puvel-zssu) -zdrag_atm*(puvel-zua) -zdrag_ice*(puvel-zui) 324 367 !zaye=-zff*puvel -grav*zssh_y +zwave_rad*zvwave & 325 ! -zdrag_ocn*(pvvel-z vo) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi)368 ! -zdrag_ocn*(pvvel-zssv) -zdrag_atm*(pvvel-zva) -zdrag_ice*(pvvel-zvi) 326 369 zaxe = -grav * zssh_x + zwave_rad * zuwave 327 370 zaye = -grav * zssh_y + zwave_rad * zvwave
Note: See TracChangeset
for help on using the changeset viewer.