- Timestamp:
- 2017-09-27T12:09:10+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90
r8564 r8565 61 61 REAL(wp) :: rn_pstar ! determines ice strength, Hibler JPO79 62 62 REAL(wp) :: rn_crhg ! determines changes in ice strength 63 LOGICAL :: ln_str_R75 ! ice strength parameterization (Rothrock75)64 REAL(wp) :: rn_perdg ! ridging work divided by pot. energy change in ridging65 63 REAL(wp) :: rn_csrdg ! fraction of shearing energy contributing to ridging 66 64 LOGICAL :: ln_partf_lin ! participation function linear (Thorndike et al. (1975)) … … 130 128 !! 131 129 INTEGER :: ji, jj, jk, jl ! dummy loop index 132 INTEGER :: niter, iterate_ridging! local integer133 INTEGER :: nidx2! local integer130 INTEGER :: iter, iterate_ridging ! local integer 131 INTEGER :: ipti ! local integer 134 132 REAL(wp) :: zfac ! local scalar 135 INTEGER , DIMENSION(jpij) :: i dxice2! compute ridge/raft or not133 INTEGER , DIMENSION(jpij) :: iptidx ! compute ridge/raft or not 136 134 REAL(wp), DIMENSION(jpij) :: zdivu_adv ! divu as implied by transport scheme (1/s) 137 135 REAL(wp), DIMENSION(jpij) :: zdivu, zdelt ! 1D divu_i & delta_i 138 136 ! 139 INTEGER, PARAMETER :: nitermax = 20137 INTEGER, PARAMETER :: jp_itermax = 20 140 138 !!------------------------------------------------------------------- 141 139 ! controls … … 154 152 ! 0) Identify grid cells with ice 155 153 !-------------------------------- 156 n idx = 0 ; idxice(:) = 0154 npti = 0 ; nptidx(:) = 0 157 155 DO jj = 1, jpj 158 156 DO ji = 1, jpi 159 157 IF ( SUM(a_i(ji,jj,:)) > 0._wp ) THEN 160 n idx = nidx+ 1161 idxice( nidx) = (jj - 1) * jpi + ji158 npti = npti + 1 159 nptidx( npti ) = (jj - 1) * jpi + ji 162 160 ENDIF 163 161 END DO 164 162 END DO 165 163 166 IF( n idx> 0 ) THEN164 IF( npti > 0 ) THEN 167 165 168 166 ! just needed here 169 CALL tab_2d_1d( n idx, idxice(1:nidx), zdivu(1:nidx), divu_i(:,:) )170 CALL tab_2d_1d( n idx, idxice(1:nidx), zdelt(1:nidx), delta_i(:,:) )167 CALL tab_2d_1d( npti, nptidx(1:npti), zdivu(1:npti), divu_i(:,:) ) 168 CALL tab_2d_1d( npti, nptidx(1:npti), zdelt(1:npti), delta_i(:,:) ) 171 169 ! needed here and in the iteration loop 172 CALL tab_3d_2d( n idx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i(:,:,:) )173 CALL tab_3d_2d( n idx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i(:,:,:) )174 CALL tab_2d_1d( n idx, idxice(1:nidx), ato_i_1d(1:nidx) , ato_i(:,:) )175 176 DO ji = 1, n idx170 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i(:,:,:) ) 171 CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i(:,:,:) ) 172 CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i(:,:) ) 173 174 DO ji = 1, npti 177 175 !-----------------------------------------------------------------------------! 178 176 ! 2) Dynamical inputs (closing rate, divu_adv, opning) … … 219 217 !------------------------------------------------ 220 218 CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d ) 221 nidx2 = 0 ; idxice2(:) = 0222 DO ji = 1, n idx219 ipti = 0 ; iptidx(:) = 0 220 DO ji = 1, npti 223 221 IF( SUM( ABS(apartf(ji,1:jpl)) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 224 nidx2 = nidx2+ 1225 i dxice2 (nidx2) = idxice(ji)226 a_i_2d ( nidx2,:) = a_i_2d (ji,:) ! adjust to new indices227 v_i_2d ( nidx2,:) = v_i_2d (ji,:) ! adjust to new indices228 ato_i_1d ( nidx2) = ato_i_1d (ji) ! adjust to new indices229 closing_net( nidx2) = closing_net(ji) ! adjust to new indices230 zdivu_adv ( nidx2) = zdivu_adv (ji) ! adjust to new indices231 opning ( nidx2) = opning (ji) ! adjust to new indices222 ipti = ipti + 1 223 iptidx (ipti) = nptidx (ji) 224 a_i_2d (ipti,:) = a_i_2d (ji,:) ! adjust to new indices 225 v_i_2d (ipti,:) = v_i_2d (ji,:) ! adjust to new indices 226 ato_i_1d (ipti) = ato_i_1d (ji) ! adjust to new indices 227 closing_net(ipti) = closing_net(ji) ! adjust to new indices 228 zdivu_adv (ipti) = zdivu_adv (ji) ! adjust to new indices 229 opning (ipti) = opning (ji) ! adjust to new indices 232 230 ENDIF 233 231 END DO 234 idxice(:) = idxice2(:)235 n idx = nidx2232 nptidx(:) = iptidx(:) 233 npti = ipti 236 234 237 235 ENDIF … … 240 238 ! Start ridging 241 239 !-------------- 242 IF( n idx> 0 ) THEN240 IF( npti > 0 ) THEN 243 241 244 242 !----------- … … 246 244 !----------- 247 245 ! fields used but not modified 248 CALL tab_2d_1d( n idx, idxice(1:nidx), sss_1d(1:nidx), sss_m(:,:) )249 CALL tab_2d_1d( n idx, idxice(1:nidx), sst_1d(1:nidx), sst_m(:,:) )246 CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) ) 247 CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) ) 250 248 ! the following fields are modified in this routine 251 !!CALL tab_2d_1d( n idx, idxice(1:nidx), ato_i_1d(1:nidx), ato_i(:,:) )252 !!CALL tab_3d_2d( n idx, idxice(1:nidx), a_i_2d(1:nidx,1:jpl), a_i(:,:,:) )253 !!CALL tab_3d_2d( n idx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i (:,:,:) )254 CALL tab_3d_2d( n idx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s (:,:,:) )255 CALL tab_3d_2d( n idx, idxice(1:nidx), sv_i_2d(1:nidx,1:jpl), sv_i(:,:,:) )256 CALL tab_3d_2d( n idx, idxice(1:nidx), oa_i_2d(1:nidx,1:jpl), oa_i(:,:,:) )249 !!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) ) 250 !!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) ) 251 !!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) ) 252 CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) ) 253 CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 254 CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) ) 257 255 IF ( nn_pnd_scheme > 0 ) THEN 258 CALL tab_3d_2d( n idx, idxice(1:nidx), a_ip_2d(1:nidx,1:jpl), a_ip(:,:,:) )259 CALL tab_3d_2d( n idx, idxice(1:nidx), v_ip_2d(1:nidx,1:jpl), v_ip(:,:,:) )256 CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 257 CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 260 258 ENDIF 261 259 DO jl = 1, jpl 262 260 DO jk = 1, nlay_s 263 CALL tab_2d_1d( n idx, idxice(1:nidx), ze_s_2d(1:nidx,jk,jl), e_s(:,:,jk,jl) )261 CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 264 262 END DO 265 263 DO jk = 1, nlay_i 266 CALL tab_2d_1d( n idx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) )267 END DO 268 END DO 269 CALL tab_2d_1d( n idx, idxice(1:nidx), sfx_dyn_1d (1:nidx), sfx_dyn (:,:) )270 CALL tab_2d_1d( n idx, idxice(1:nidx), sfx_bri_1d (1:nidx), sfx_bri (:,:) )271 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_dyn_1d (1:nidx), wfx_dyn (:,:) )272 CALL tab_2d_1d( n idx, idxice(1:nidx), hfx_dyn_1d (1:nidx), hfx_dyn (:,:) )273 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_snw_dyn_1d(1:nidx), wfx_snw_dyn(:,:) )264 CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 265 END DO 266 END DO 267 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_dyn_1d (1:npti), sfx_dyn (:,:) ) 268 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri (:,:) ) 269 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_dyn_1d (1:npti), wfx_dyn (:,:) ) 270 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d (1:npti), hfx_dyn (:,:) ) 271 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 274 272 IF ( nn_pnd_scheme > 0 ) THEN 275 CALL tab_2d_1d( n idx, idxice(1:nidx), wfx_pnd_1d(1:nidx), wfx_pnd(:,:) )273 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 276 274 ENDIF 277 275 … … 279 277 ! 3) Ridging iteration 280 278 !-----------------------------------------------------------------------------! 281 niter= 1279 iter = 1 282 280 iterate_ridging = 1 283 DO WHILE( iterate_ridging > 0 .AND. niter < nitermax )281 DO WHILE( iterate_ridging > 0 .AND. iter < jp_itermax ) 284 282 285 283 CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d ) … … 294 292 ! rates were reduced above), ridge again with new rates. 295 293 iterate_ridging = 0 296 DO ji = 1, n idx294 DO ji = 1, npti 297 295 zfac = 1._wp - ( ato_i_1d(ji) + SUM( a_i_2d(ji,:) ) ) 298 296 IF( ABS( zfac ) < epsi10 ) THEN … … 308 306 END DO 309 307 ! 310 niter = niter + 1311 IF( niter > nitermax ) CALL ctl_warn( 'icedyn_rdgrft: non-converging ridging scheme' )308 iter = iter + 1 309 IF( iter > jp_itermax ) CALL ctl_warn( 'icedyn_rdgrft: non-converging ridging scheme' ) 312 310 ! 313 311 END DO … … 316 314 ! 1D <==> 2D 317 315 !----------- 318 CALL tab_1d_2d( n idx, idxice(1:nidx), ato_i_1d(1:nidx), ato_i(:,:) )319 CALL tab_2d_3d( n idx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i (:,:,:) )320 CALL tab_2d_3d( n idx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i (:,:,:) )321 CALL tab_2d_3d( n idx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s (:,:,:) )322 CALL tab_2d_3d( n idx, idxice(1:nidx), sv_i_2d(1:nidx,1:jpl), sv_i(:,:,:) )323 CALL tab_2d_3d( n idx, idxice(1:nidx), oa_i_2d(1:nidx,1:jpl), oa_i(:,:,:) )316 CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) ) 317 CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) ) 318 CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) ) 319 CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) ) 320 CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 321 CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) ) 324 322 IF ( nn_pnd_scheme > 0 ) THEN 325 CALL tab_2d_3d( n idx, idxice(1:nidx), a_ip_2d(1:nidx,1:jpl), a_ip(:,:,:) )326 CALL tab_2d_3d( n idx, idxice(1:nidx), v_ip_2d(1:nidx,1:jpl), v_ip(:,:,:) )323 CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 324 CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 327 325 ENDIF 328 326 DO jl = 1, jpl 329 327 DO jk = 1, nlay_s 330 CALL tab_1d_2d( n idx, idxice(1:nidx), ze_s_2d(1:nidx,jk,jl), e_s(:,:,jk,jl) )328 CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 331 329 END DO 332 330 DO jk = 1, nlay_i 333 CALL tab_1d_2d( n idx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) )334 END DO 335 END DO 336 CALL tab_1d_2d( n idx, idxice(1:nidx), sfx_dyn_1d (1:nidx), sfx_dyn (:,:) )337 CALL tab_1d_2d( n idx, idxice(1:nidx), sfx_bri_1d (1:nidx), sfx_bri (:,:) )338 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_dyn_1d (1:nidx), wfx_dyn (:,:) )339 CALL tab_1d_2d( n idx, idxice(1:nidx), hfx_dyn_1d (1:nidx), hfx_dyn (:,:) )340 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_snw_dyn_1d(1:nidx), wfx_snw_dyn(:,:) )331 CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 332 END DO 333 END DO 334 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d (1:npti), sfx_dyn (:,:) ) 335 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri (:,:) ) 336 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d (1:npti), wfx_dyn (:,:) ) 337 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d (1:npti), hfx_dyn (:,:) ) 338 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 341 339 IF ( nn_pnd_scheme > 0 ) THEN 342 CALL tab_1d_2d( n idx, idxice(1:nidx), wfx_pnd_1d(1:nidx), wfx_pnd(:,:) )340 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 343 341 ENDIF 344 342 345 ENDIF ! n idx> 0343 ENDIF ! npti > 0 346 344 347 345 CALL ice_var_agg( 1 ) … … 378 376 379 377 ! ! Ice thickness needed for rafting 380 WHERE( pa_i(1:n idx,:) > 0._wp ) ; zhi(1:nidx,:) = pv_i(1:nidx,:) / pa_i(1:nidx,:)381 ELSEWHERE ; zhi(1:n idx,:) = 0._wp378 WHERE( pa_i(1:npti,:) > 0._wp ) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 379 ELSEWHERE ; zhi(1:npti,:) = 0._wp 382 380 END WHERE 383 381 … … 398 396 ! Compute total area of ice plus open water. 399 397 ! This is in general not equal to one because of divergence during transport 400 zasum(1:n idx) = pato_i(1:nidx) + SUM( pa_i(1:nidx,:), dim=2 )401 ! 402 WHERE( zasum(1:n idx) > 0._wp ) ; z1_asum(1:nidx) = 1._wp / zasum(1:nidx)403 ELSEWHERE ; z1_asum(1:n idx) = 0._wp398 zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 399 ! 400 WHERE( zasum(1:npti) > 0._wp ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti) 401 ELSEWHERE ; z1_asum(1:npti) = 0._wp 404 402 END WHERE 405 403 ! Compute cumulative thickness distribution function … … 407 405 ! where zGsum(n) is the fractional area in categories 0 to n. 408 406 ! initial value (in h = 0) equals open water area 409 zGsum(1:n idx,-1) = 0._wp410 zGsum(1:n idx,0 ) = pato_i(1:nidx) * z1_asum(1:nidx)407 zGsum(1:npti,-1) = 0._wp 408 zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti) 411 409 DO jl = 1, jpl 412 zGsum(1:n idx,jl) = ( pato_i(1:nidx) + SUM( pa_i(1:nidx,1:jl), dim=2 ) ) * z1_asum(1:nidx) ! sum(1:jl) is ok (and not jpl)410 zGsum(1:npti,jl) = ( pato_i(1:npti) + SUM( pa_i(1:npti,1:jl), dim=2 ) ) * z1_asum(1:npti) ! sum(1:jl) is ok (and not jpl) 413 411 END DO 414 412 ! 415 413 IF( ln_partf_lin ) THEN !--- Linear formulation (Thorndike et al., 1975) 416 414 DO jl = 0, jpl 417 DO ji = 1, n idx415 DO ji = 1, npti 418 416 IF ( zGsum(ji,jl) < rn_gstar ) THEN 419 417 apartf(ji,jl) = z1_gstar * ( zGsum(ji,jl) - zGsum(ji,jl-1) ) * & … … 432 430 zfac = 1._wp / ( 1._wp - EXP(-z1_astar) ) 433 431 DO jl = -1, jpl 434 DO ji = 1, n idx432 DO ji = 1, npti 435 433 zGsum(ji,jl) = EXP( -zGsum(ji,jl) * z1_astar ) * zfac 436 434 END DO 437 435 END DO 438 436 DO jl = 0, jpl 439 DO ji = 1, n idx437 DO ji = 1, npti 440 438 apartf(ji,jl) = zGsum(ji,jl-1) - zGsum(ji,jl) 441 439 END DO … … 447 445 IF( ln_rafting .AND. ln_ridging ) THEN !- ridging & rafting 448 446 DO jl = 1, jpl 449 DO ji = 1, n idx447 DO ji = 1, npti 450 448 aridge(ji,jl) = ( 1._wp + TANH ( rn_craft * ( zhi(ji,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jl) 451 449 araft (ji,jl) = apartf(ji,jl) - aridge(ji,jl) … … 454 452 ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN !- ridging alone 455 453 DO jl = 1, jpl 456 DO ji = 1, n idx454 DO ji = 1, npti 457 455 aridge(ji,jl) = apartf(ji,jl) 458 456 araft (ji,jl) = 0._wp … … 461 459 ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN !- rafting alone 462 460 DO jl = 1, jpl 463 DO ji = 1, n idx461 DO ji = 1, npti 464 462 aridge(ji,jl) = 0._wp 465 463 araft (ji,jl) = apartf(ji,jl) … … 468 466 ELSE !- no ridging & no rafting 469 467 DO jl = 1, jpl 470 DO ji = 1, n idx468 DO ji = 1, npti 471 469 aridge(ji,jl) = 0._wp 472 470 araft (ji,jl) = 0._wp … … 502 500 503 501 zfac = 1._wp / hi_hrft 504 zaksum(1:n idx) = apartf(1:nidx,0)502 zaksum(1:npti) = apartf(1:npti,0) 505 503 ! Transfer function 506 504 DO jl = 1, jpl !all categories have a specific transfer function 507 DO ji = 1, n idx505 DO ji = 1, npti 508 506 IF ( apartf(ji,jl) > 0._wp ) THEN 509 507 zhmean = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min ) … … 530 528 ! closing rate to a gross closing rate. 531 529 ! NOTE: 0 < aksum <= 1 532 WHERE( zaksum(1:n idx) > 0._wp ) ; closing_gross(1:nidx) = closing_net(1:nidx) / zaksum(1:nidx)533 ELSEWHERE ; closing_gross(1:n idx) = 0._wp530 WHERE( zaksum(1:npti) > 0._wp ) ; closing_gross(1:npti) = closing_net(1:npti) / zaksum(1:npti) 531 ELSEWHERE ; closing_gross(1:npti) = 0._wp 534 532 END WHERE 535 533 536 DO ji = 1, n idx534 DO ji = 1, npti 537 535 ! correction to closing rate and opening if closing rate is excessive 538 536 !--------------------------------------------------------------------- … … 552 550 ! would be removed. Reduce the opening rate proportionately. 553 551 DO jl = 1, jpl 554 DO ji = 1, n idx552 DO ji = 1, npti 555 553 zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice 556 554 !! IF( zfac > pa_i(ji,jl) .AND. zfac > 0._wp ) THEN … … 594 592 ! 1) Compute change in open water area due to closing and opening. 595 593 !------------------------------------------------------------------------------- 596 DO ji = 1, n idx594 DO ji = 1, npti 597 595 ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice ) 598 596 END DO … … 603 601 DO jl1 = 1, jpl 604 602 605 CALL tab_2d_1d( n idx, idxice(1:nidx), s_i_1d(1:nidx), s_i(:,:,jl1) )606 607 DO ji = 1, n idx603 CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,jl1) ) 604 605 DO ji = 1, npti 608 606 609 607 !------------------------------------------------ … … 716 714 ! special loop for e_i because of layers jk 717 715 DO jk = 1, nlay_i 718 DO ji = 1, n idx716 DO ji = 1, npti 719 717 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 720 718 ! Compute ridging /rafting fractions … … 736 734 DO jl2 = 1, jpl 737 735 ! 738 DO ji = 1, n idx736 DO ji = 1, npti 739 737 740 738 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN … … 783 781 784 782 DO jk = 1, nlay_i 785 DO ji = 1, n idx783 DO ji = 1, npti 786 784 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) & 787 785 & ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji) … … 794 792 ! 795 793 ! In case ridging/rafting lead to very small negative values (sometimes it happens) 796 WHERE( a_i_2d(1:n idx,:) < 0._wp ) a_i_2d(1:nidx,:) = 0._wp797 WHERE( v_i_2d(1:n idx,:) < 0._wp ) v_i_2d(1:nidx,:) = 0._wp794 WHERE( a_i_2d(1:npti,:) < 0._wp ) a_i_2d(1:npti,:) = 0._wp 795 WHERE( v_i_2d(1:npti,:) < 0._wp ) v_i_2d(1:npti,:) = 0._wp 798 796 ! 799 797 END SUBROUTINE rdgrft_shift … … 820 818 !!---------------------------------------------------------------------- 821 819 ! !--------------------------------------------------! 822 IF( ln_str_R75 ) THEN ! Ice strength => Rothrock (1975) method ! 823 ! !--------------------------------------------------! 824 ! !--------------------------------------------------! 825 ELSEIF( ln_str_H79 ) THEN ! Ice strength => Hibler (1979) method ! 820 IF( ln_str_H79 ) THEN ! Ice strength => Hibler (1979) method ! 826 821 ! !--------------------------------------------------! 827 822 strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 828 !829 823 ismooth = 1 830 ! 824 ! !--------------------------------------------------! 825 ELSE ! Zero strength ! 826 ! !--------------------------------------------------! 827 strength(:,:) = 0._wp 828 ismooth = 0 831 829 ENDIF 832 830 ! !--------------------------------------------------! … … 896 894 !! 897 895 NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, & 898 & ln_str_R75, rn_perdg, &899 896 & rn_csrdg , & 900 897 & ln_partf_lin, rn_gstar, & … … 921 918 WRITE(numout,*) ' 1st bulk-rheology parameter rn_pstar = ', rn_pstar 922 919 WRITE(numout,*) ' 2nd bulk-rhelogy parameter rn_crhg = ', rn_crhg 923 WRITE(numout,*) ' ice strength parameterization Rothrock (1975) ln_str_R75 = ', ln_str_R75924 WRITE(numout,*) ' Ratio of ridging work to PotEner change in ridging rn_perdg = ', rn_perdg925 920 WRITE(numout,*) ' Fraction of shear energy contributing to ridging rn_csrdg = ', rn_csrdg 926 921 WRITE(numout,*) ' linear ridging participation function ln_partf_lin = ', ln_partf_lin … … 940 935 ENDIF 941 936 ! 942 IF ( ( ln_str_H79 .AND. ln_str_R75 ) .OR. ( .NOT.ln_str_H79 .AND. .NOT.ln_str_R75 ) ) THEN943 CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one formulation for ice strength (ln_str_H79 or ln_str_R75)' )944 ENDIF945 !946 937 IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN 947 938 CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' )
Note: See TracChangeset
for help on using the changeset viewer.
