Changeset 9452
- Timestamp:
- 2018-03-30T18:12:48+02:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90
r9271 r9452 106 106 !! 107 107 !! ** Method : Steps : 108 !! 1) Thickness categories boundaries, ice / o.w. concentrations 109 !! Ridge preparation 110 !! 2) Dynamical inputs (closing rate, divu_adv, opning) 111 !! 3) Ridging iteration 112 !! 4) Ridging diagnostics 113 !! 5) Heat, salt and freshwater fluxes 114 !! 6) Compute increments of tate variables and come back to old values 108 !! 0) Identify grid cells with ice 109 !! 1) Calculate closing rate, divergence and opening 110 !! 2) Identify grid cells with ridging 111 !! 3) Start ridging iterations 112 !! - prep = ridged and rafted ice + closing_gross 113 !! - shift = move ice from one category to another 114 !! 115 !! ** Details 116 !! step1: The net rate of closing is due to convergence and shear, based on Flato and Hibler (1995). 117 !! The energy dissipation rate is equal to the net closing rate times the ice strength. 118 !! 119 !! step3: The gross closing rate is equal to the first two terms (open 120 !! water closing and thin ice ridging) without the third term 121 !! (thick, newly ridged ice). 115 122 !! 116 123 !! References : Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626. … … 123 130 !! 124 131 !! This routine is based on CICE code and authors William H. Lipscomb, 125 !! and Elizabeth C. Hunke, LANL are gratefully acknowledged132 !! and Elizabeth C. Hunke, LANL are gratefully acknowledged 126 133 !!------------------------------------------------------------------- 127 134 INTEGER, INTENT(in) :: kt ! number of iteration … … 153 160 !-------------------------------- 154 161 npti = 0 ; nptidx(:) = 0 162 ipti = 0 ; iptidx(:) = 0 155 163 DO jj = 1, jpj 156 164 DO ji = 1, jpi 157 IF ( SUM(a_i(ji,jj,:)) > 0._wp ) THEN165 IF ( at_i(ji,jj) > 0._wp ) THEN 158 166 npti = npti + 1 159 167 nptidx( npti ) = (jj - 1) * jpi + ji … … 161 169 END DO 162 170 END DO 163 171 172 !-------------------------------------------------------- 173 ! 1) Dynamical inputs (closing rate, divergence, opening) 174 !-------------------------------------------------------- 164 175 IF( npti > 0 ) THEN 165 176 … … 173 184 174 185 DO ji = 1, npti 175 ! -----------------------------------------------------------------------------!176 ! 2) Dynamical inputs (closing rate, divu_adv, opning)177 !-----------------------------------------------------------------------------!186 ! closing_net = rate at which open water area is removed + ice area removed by ridging 187 ! - ice area added in new ridges 188 closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 178 189 ! 179 ! 2.1 closing_net 180 !----------------- 181 ! Compute the net rate of closing due to convergence 182 ! and shear, based on Flato and Hibler (1995). 183 ! 184 ! The energy dissipation rate is equal to the net closing rate 185 ! times the ice strength. 190 ! divergence given by the advection scheme 191 ! (which may not be equal to divu as computed from the velocity field) 192 zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice 186 193 ! 187 ! NOTE: The NET closing rate is equal to the rate that open water 188 ! area is removed, plus the rate at which ice area is removed by 189 ! ridging, minus the rate at which area is added in new ridges. 190 ! The GROSS closing rate is equal to the first two terms (open 191 ! water closing and thin ice ridging) without the third term 192 ! (thick, newly ridged ice). 193 194 closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 195 196 ! 2.2 divu_adv 197 !-------------- 198 ! Compute divu_adv, the divergence rate given by the transport/ 199 ! advection scheme, which may not be equal to divu as computed 200 ! from the velocity field. 201 ! 202 ! If divu_adv < 0, make sure the closing rate is large enough 203 ! to give asum = 1.0 after ridging. 204 zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice 205 206 IF( zdivu_adv(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) ) 207 208 ! 2.3 opning 209 !------------ 210 ! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging. 194 IF( zdivu_adv(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) ) ! make sure the closing rate is large enough 195 ! ! to give asum = 1.0 after ridging 196 ! Opening rate (non-negative) that will give asum = 1.0 after ridging. 211 197 opning(ji) = closing_net(ji) + zdivu_adv(ji) 212 198 END DO 213 214 199 ! 215 !------------------------------------ ------------216 ! 2 .1) Identify grid cells with nonzeroridging217 !------------------------------------ ------------218 CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d )219 ipti = 0 ; iptidx(:) = 0 200 !------------------------------------ 201 ! 2) Identify grid cells with ridging 202 !------------------------------------ 203 CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net ) 204 220 205 DO ji = 1, npti 221 IF( SUM( ABS(apartf(ji,1:jpl)) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN206 IF( SUM( apartf(ji,1:jpl) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 222 207 ipti = ipti + 1 223 208 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 209 ! adjust to new indices 210 a_i_2d (ipti,:) = a_i_2d (ji,:) 211 v_i_2d (ipti,:) = v_i_2d (ji,:) 212 ato_i_1d (ipti) = ato_i_1d (ji) 213 closing_net(ipti) = closing_net(ji) 214 zdivu_adv (ipti) = zdivu_adv (ji) 215 opning (ipti) = opning (ji) 230 216 ENDIF 231 217 END DO 232 nptidx(:) = iptidx(:)233 npti = ipti234 218 235 219 ENDIF 236 220 237 !-------------- 238 ! Start ridging 239 !-------------- 221 ! grid cells with ridging 222 nptidx(:) = iptidx(:) 223 npti = ipti 224 225 !----------------- 226 ! 3) Start ridging 227 !----------------- 240 228 IF( npti > 0 ) THEN 241 229 242 !----------- 243 ! 2D <==> 1D 244 !----------- 245 ! fields used but not modified 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(:,:) ) 248 ! the following fields are modified in this routine 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(:,:,:) ) 255 CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 256 CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 257 DO jl = 1, jpl 258 DO jk = 1, nlay_s 259 CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 260 END DO 261 DO jk = 1, nlay_i 262 CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 263 END DO 264 END DO 265 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_dyn_1d (1:npti), sfx_dyn (:,:) ) 266 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri (:,:) ) 267 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_dyn_1d (1:npti), wfx_dyn (:,:) ) 268 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d (1:npti), hfx_dyn (:,:) ) 269 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 270 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 271 272 !-----------------------------------------------------------------------------! 273 ! 3) Ridging iteration 274 !-----------------------------------------------------------------------------! 230 CALL ice_dyn_1d2d( 1 ) ! --- Move to 1D arrays --- ! 231 275 232 iter = 1 276 233 iterate_ridging = 1 277 DO WHILE( iterate_ridging > 0 .AND. iter < jp_itermax ) 278 279 CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d ) 280 281 ! 3.2 Redistribute area, volume, and energy. 282 !-----------------------------------------------------------------------------! 234 ! !----------------------! 235 DO WHILE( iterate_ridging > 0 .AND. iter < jp_itermax ) ! ridging iterations ! 236 ! !----------------------! 237 ! Calculate participation function (apartf) 238 ! and transfer function 239 ! and closing_gross (+correction on opening) 240 CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net ) 241 242 ! Redistribute area, volume, and energy between categories 283 243 CALL rdgrft_shift 284 244 285 ! 3.4Do we keep on iterating?286 !------------------------- ----------------------------------------------------!287 ! Check whether a sum = 1. If not (because the closing and opening288 ! rates were reduced above), ridge again with new rates.245 ! Do we keep on iterating? 246 !------------------------- 247 ! Check whether a_i + ato_i = 0 248 ! If not, because the closing and opening rates were reduced above, ridge again with new rates 289 249 iterate_ridging = 0 290 250 DO ji = 1, npti … … 307 267 END DO 308 268 309 !----------- 310 ! 1D <==> 2D 311 !----------- 312 CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) ) 313 CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) ) 314 CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) ) 315 CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) ) 316 CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 317 CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) ) 318 CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 319 CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 320 DO jl = 1, jpl 321 DO jk = 1, nlay_s 322 CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 323 END DO 324 DO jk = 1, nlay_i 325 CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 326 END DO 327 END DO 328 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d (1:npti), sfx_dyn (:,:) ) 329 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri (:,:) ) 330 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d (1:npti), wfx_dyn (:,:) ) 331 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d (1:npti), hfx_dyn (:,:) ) 332 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 333 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 334 335 ENDIF ! npti > 0 269 CALL ice_dyn_1d2d( 2 ) ! --- Move to 2D arrays --- ! 270 271 ENDIF 336 272 337 273 CALL ice_var_agg( 1 ) … … 345 281 346 282 347 SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i )283 SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i, pclosing_net ) 348 284 !!------------------------------------------------------------------- 349 285 !! *** ROUTINE rdgrft_prep *** … … 352 288 !! 353 289 !! ** Method : Compute the thickness distribution of the ice and open water 354 !! participating in ridging and of the resulting ridges.355 !!------------------------------------------------------------------- 356 REAL(wp), DIMENSION(:) , INTENT(in) :: pato_i 290 !! participating in ridging and of the resulting ridges. 291 !!------------------------------------------------------------------- 292 REAL(wp), DIMENSION(:) , INTENT(in) :: pato_i, pclosing_net 357 293 REAL(wp), DIMENSION(:,:), INTENT(in) :: pa_i, pv_i 358 294 !! … … 372 308 END WHERE 373 309 310 ! 1) Participation function (apartf): a(h) = b(h).g(h) 374 311 !----------------------------------------------------------------- 375 ! 1) Participation function: a(h) = b(h).g(h) (apartf) 376 !----------------------------------------------------------------- 377 ! Compute the participation function apartf; this is analogous to 378 ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975). 379 ! area lost from category n due to ridging/closing 380 ! apartf(n) = total area lost due to ridging/closing 381 ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar). 382 ! 383 ! The expressions for apartf are found by integrating b(h)g(h) between 384 ! the category boundaries. 312 ! Compute the participation function = total area lost due to ridging/closing 313 ! This is analogous to 314 ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975). 315 ! assuming b(h) = (2/Gstar) * (1 - G(h)/Gstar). 316 ! 317 ! apartf = integrating b(h)g(h) between the category boundaries 385 318 ! apartf is always >= 0 and SUM(apartf(0:jpl))=1 386 319 !----------------------------------------------------------------- … … 393 326 ELSEWHERE ; z1_asum(1:npti) = 0._wp 394 327 END WHERE 328 ! 395 329 ! Compute cumulative thickness distribution function 396 330 ! Compute the cumulative thickness distribution function zGsum, 397 331 ! where zGsum(n) is the fractional area in categories 0 to n. 398 ! initial value (in h = 0) equalsopen water area332 ! initial value (in h = 0) = open water area 399 333 zGsum(1:npti,-1) = 0._wp 400 334 zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti) … … 465 399 ENDIF 466 400 467 !-----------------------------------------------------------------468 401 ! 2) Transfer function 469 402 !----------------------------------------------------------------- … … 479 412 ! constrained by hrmin <= (zhmean + hi)/2. 480 413 ! 481 ! The maximum ridging thickness, hrmax, is determined by 482 ! zhmean and hrmin. 414 ! The maximum ridging thickness, hrmax, is determined by zhmean and hrmin. 483 415 ! 484 416 ! These modifications have the effect of reducing the ice strength 485 ! (relative to the Hibler formulation) when very thick ice is 486 ! ridging. 417 ! (relative to the Hibler formulation) when very thick ice is ridging. 487 418 ! 488 419 ! zaksum = net area removed/ total area removed … … 490 421 ! net area removed = total area removed - area of new ridges 491 422 !----------------------------------------------------------------- 492 493 423 zfac = 1._wp / hi_hrft 494 424 zaksum(1:npti) = apartf(1:npti,0) 495 ! Transfer function496 DO jl = 1, jpl !all categories have a specific transfer function425 ! 426 DO jl = 1, jpl 497 427 DO ji = 1, npti 498 428 IF ( apartf(ji,jl) > 0._wp ) THEN … … 515 445 END DO 516 446 ! 517 ! 3.1 closing_gross 518 !-----------------------------------------------------------------------------! 519 ! Based on the ITD of ridging and ridged ice, convert the net 520 ! closing rate to a gross closing rate. 447 ! 3) closing_gross 448 !----------------- 449 ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate. 521 450 ! NOTE: 0 < aksum <= 1 522 WHERE( zaksum(1:npti) > 0._wp ) ; closing_gross(1:npti) = closing_net(1:npti) / zaksum(1:npti)451 WHERE( zaksum(1:npti) > 0._wp ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 523 452 ELSEWHERE ; closing_gross(1:npti) = 0._wp 524 453 END WHERE 525 454 455 ! correction to closing rate if excessive ice removal 456 !---------------------------------------------------- 457 ! Reduce the closing rate if more than 100% of any ice category would be removed 458 ! Reduce the opening rate in proportion 459 DO jl = 1, jpl 460 DO ji = 1, npti 461 zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice 462 IF( zfac > pa_i(ji,jl) ) THEN 463 closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice 464 ENDIF 465 END DO 466 END DO 467 468 ! 4) correction to opening if excessive open water removal 469 !--------------------------------------------------------- 470 ! Reduce the closing rate if more than 100% of the open water would be removed 471 ! Reduce the opening rate in proportion 526 472 DO ji = 1, npti 527 ! correction to closing rate and opening if closing rate is excessive 528 !--------------------------------------------------------------------- 529 ! Reduce the closing rate if more than 100% of the open water 530 ! would be removed. Reduce the opening rate proportionately. 531 zfac = ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice 532 IF ( zfac < 0._wp .AND. zfac > - pato_i(ji) ) THEN ! would lead to negative ato_i 473 zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice 474 IF( zfac < 0._wp ) THEN ! would lead to negative ato_i 533 475 opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_rdtice 534 ELSEIF( zfac > 0._wp .AND. zfac > ( zasum(ji) - pato_i(ji) ) ) THEN! would lead to ato_i > asum476 ELSEIF( zfac > zasum(ji) ) THEN ! would lead to ato_i > asum 535 477 opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_rdtice 536 478 ENDIF 537 479 END DO 538 539 ! correction to closing rate / opening if excessive ice removal540 !---------------------------------------------------------------541 ! Reduce the closing rate if more than 100% of any ice category542 ! would be removed. Reduce the opening rate proportionately.543 DO jl = 1, jpl544 DO ji = 1, npti545 zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice546 !! IF( zfac > pa_i(ji,jl) .AND. zfac > 0._wp ) THEN547 IF( zfac > pa_i(ji,jl) ) THEN548 closing_gross(ji) = closing_gross(ji) * pa_i(ji,jl) / zfac549 ENDIF550 END DO551 END DO552 480 ! 553 481 END SUBROUTINE rdgrft_prep … … 561 489 !! 562 490 !! ** Method : Remove area, volume, and energy from each ridging category 563 !! and add to thicker ice categories.491 !! and add to thicker ice categories. 564 492 !!------------------------------------------------------------------- 565 493 ! … … 581 509 REAL(wp), DIMENSION(jpij,nlay_s) :: esrdg ! enth*volume of new ridges 582 510 REAL(wp), DIMENSION(jpij,nlay_i) :: eirdg ! enth*volume of new ridges 583 !!------------------------------------------------------------------- 584 585 !------------------------------------------------------------------------------- 586 ! 1) Compute change in open water area due to closing and opening. 587 !------------------------------------------------------------------------------- 511 ! 512 INTEGER , DIMENSION(jpij) :: itest_rdg, itest_rft ! test for conservation 513 !!------------------------------------------------------------------- 514 515 ! 1) Change in open water area due to closing and opening 516 !-------------------------------------------------------- 588 517 DO ji = 1, npti 589 518 ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice ) 590 519 END DO 591 520 592 !----------------------------------------------------------------- 593 ! 2) compute categories participating to ridging/rafting (jl1) 594 !----------------------------------------------------------------- 521 ! 2) compute categories in which ice is removed (jl1) 522 !---------------------------------------------------- 595 523 DO jl1 = 1, jpl 596 524 … … 599 527 DO ji = 1, npti 600 528 601 !------------------------------------------------ 602 ! 2.1) Identify grid cells with nonzero ridging 603 !------------------------------------------------ 604 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 529 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN ! only if ice is ridging 605 530 606 531 z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 607 532 608 !-------------------------------------------------------------------- 609 ! 2.2) Compute area of ridging ice (airdg1) and of new ridge (airdg2) 610 !-------------------------------------------------------------------- 533 ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2) 611 534 airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice 612 535 airft1 = araft (ji,jl1) * closing_gross(ji) * rdt_ice … … 615 538 airft2(ji) = airft1 * hi_hrft 616 539 617 !--------------------------------------------------------------- 618 ! 2.3) Compute ridging /rafting fractions, make sure afrdg <=1 619 !--------------------------------------------------------------- 540 ! ridging /rafting fractions 620 541 afrdg = airdg1 * z1_ai(ji) 621 542 afrft = airft1 * z1_ai(ji) … … 625 546 ersw(ji) = -rhoic * vsw * rcp * sst_1d(ji) ! clem: if sst>0, then ersw <0 (is that possible?) 626 547 627 !--------------------------------------------------------------- 628 ! 2.4) Compute ridging ice and new ridges (vi, vs, sm, oi, es, ei) 629 !--------------------------------------------------------------- 548 ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei) 630 549 virdg1 = v_i_2d (ji,jl1) * afrdg 631 550 virdg2(ji) = v_i_2d (ji,jl1) * afrdg * ( 1. + rn_porordg ) … … 651 570 ENDIF 652 571 653 !----------------------------------------------------------------- 654 ! 2.5) Ice-ocean exchanges 655 !----------------------------------------------------------------- 572 ! Ice-ocean exchanges associated with ice porosity 656 573 wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoic * r1_rdtice ! increase in ice volume due to seawater frozen in voids 657 574 sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoic * r1_rdtice … … 659 576 660 577 ! Put the snow lost by ridging into the ocean 661 !------------------------------------------------662 578 ! Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow. 663 579 wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhosn * vsrdg(ji) * ( 1._wp - rn_fsnwrdg ) & ! fresh water source for ocean … … 665 581 666 582 ! Put the melt pond water into the ocean 667 !------------------------------------------668 583 ! clem: I think the following lines must be commented since there 669 584 ! is no net mass flux between melt ponds and the ocean (see icethd_pnd.F90 for ex.) … … 674 589 675 590 ! virtual salt flux to keep salinity constant 676 !------------------------------------------677 591 IF( nn_icesal /= 2 ) THEN 678 592 sirdg2(ji) = sirdg2(ji) - vsw * ( sss_1d(ji) - s_i_1d(ji) ) ! ridge salinity = s_i … … 681 595 ENDIF 682 596 683 684 !------------------------------------------ 685 ! 2.6) update jl1 (removing ridged/rafted area) 686 !------------------------------------------ 597 ! Remove area, volume of new ridge to each category jl1 598 !------------------------------------------------------ 687 599 a_i_2d (ji,jl1) = a_i_2d (ji,jl1) - airdg1 - airft1 688 600 v_i_2d (ji,jl1) = v_i_2d (ji,jl1) - virdg1 - virft(ji) … … 705 617 afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji) 706 618 afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji) 707 ! Compute ridging ice and new ridges for es619 ! Compute ridging /rafting ice and new ridges for es 708 620 esrdg(ji,jk) = ze_s_2d (ji,jk,jl1) * afrdg 709 621 esrft(ji,jk) = ze_s_2d (ji,jk,jl1) * afrft … … 711 623 hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji,jk) * ( 1._wp - rn_fsnwrdg ) & ! heat sink for ocean (<0, W.m-2) 712 624 & - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice 713 ! Update jl1 625 ! 626 ! Remove energy of new ridge to each category jl1 627 !------------------------------------------------- 714 628 ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 715 629 ENDIF … … 727 641 eirdg(ji,jk) = ze_i_2d (ji,jk,jl1) * afrdg + ersw(ji) * r1_nlay_i 728 642 eirft(ji,jk) = ze_i_2d (ji,jk,jl1) * afrft 729 ! Update jl1 643 ! 644 ! Remove energy of new ridge to each category jl1 645 !------------------------------------------------- 730 646 ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 731 647 ENDIF … … 733 649 END DO 734 650 735 736 !-------------------------------------------------- -----------------------------737 ! 3) Add area, volume, and energy of new ridge to each category jl2738 !-------------------------------------------------------------------------------651 ! 3) compute categories in which ice is added (jl2) 652 !-------------------------------------------------- 653 itest_rdg(1:npti) = 0 654 itest_rft(1:npti) = 0 739 655 DO jl2 = 1, jpl 740 656 ! … … 743 659 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 744 660 745 ! Compute the fraction of ridged ice area and volume going to thickness category jl2 .661 ! Compute the fraction of ridged ice area and volume going to thickness category jl2 746 662 IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN 747 !-----------------------------------------------------------------748 ! 2.8 Compute quantities used to apportion ice among categories749 ! in the n2 loop below750 !-----------------------------------------------------------------751 663 hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) ) 752 664 hR = MIN( hrmax(ji,jl1), hi_max(jl2) ) 753 665 farea = ( hR - hL ) / ( hrmax(ji,jl1) - hrmin(ji,jl1) ) 754 666 fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) ) 667 ! 668 itest_rdg(ji) = 1 ! test for conservation 755 669 ELSE 756 670 farea = 0._wp … … 759 673 760 674 ! Compute the fraction of rafted ice area and volume going to thickness category jl2 761 !!gm see above IF( hraft(ji) <= hi_max(jl2) .AND. hraft(ji) > hi_max(jl2-1) ) THEN 762 IF( hi_max(jl2-1) < hraft(ji,jl1) .AND. hraft(ji,jl1) <= hi_max(jl2) ) THEN ; zswitch(ji) = 1._wp 763 ELSE ; zswitch(ji) = 0._wp 675 IF( hraft(ji,jl1) <= hi_max(jl2) .AND. hraft(ji,jl1) > hi_max(jl2-1) ) THEN 676 zswitch(ji) = 1._wp 677 ! 678 itest_rft(ji) = 1 ! test for conservation 679 ELSE 680 zswitch(ji) = 0._wp 764 681 ENDIF 765 682 ! 683 ! Patch to ensure perfect conservation if ice thickness goes mad 684 ! Sometimes thickness is larger than hi_max(jpl) because of advection scheme (for very small areas) 685 ! Then ice volume is removed from one category but the ridging/rafting scheme 686 ! does not know where to move it, leading to a conservation issue. 687 IF( itest_rdg(ji) == 0 .AND. jl2 == jpl ) THEN ; farea = 1._wp ; fvol(ji) = 1._wp ; ENDIF 688 IF( itest_rft(ji) == 0 .AND. jl2 == jpl ) zswitch(ji) = 1._wp 689 ! 690 ! Add area, volume of new ridge to category jl2 691 !---------------------------------------------- 766 692 a_i_2d (ji,jl2) = a_i_2d (ji,jl2) + ( airdg2(ji) * farea + airft2(ji) * zswitch(ji) ) 767 693 oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ( oirdg2(ji) * farea + oirft2(ji) * zswitch(ji) ) … … 780 706 781 707 END DO 782 ! for snow enthalpy 708 ! Add snow energy of new ridge to category jl2 709 !--------------------------------------------- 783 710 DO jk = 1, nlay_s 784 711 DO ji = 1, npti … … 788 715 END DO 789 716 END DO 790 ! for ice enthalpy 717 ! Add ice energy of new ridge to category jl2 718 !-------------------------------------------- 791 719 DO jk = 1, nlay_i 792 720 DO ji = 1, npti … … 886 814 END SUBROUTINE ice_strength 887 815 816 817 SUBROUTINE ice_dyn_1d2d( kn ) 818 !!----------------------------------------------------------------------- 819 !! *** ROUTINE ice_dyn_1d2d *** 820 !! 821 !! ** Purpose : move arrays from 1d to 2d and the reverse 822 !!----------------------------------------------------------------------- 823 INTEGER, INTENT(in) :: kn ! 1= from 2D to 1D ; 2= from 1D to 2D 824 ! 825 INTEGER :: jl, jk ! dummy loop indices 826 !!----------------------------------------------------------------------- 827 ! 828 SELECT CASE( kn ) 829 ! !---------------------! 830 CASE( 1 ) !== from 2D to 1D ==! 831 ! !---------------------! 832 ! fields used but not modified 833 CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) ) 834 CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) ) 835 ! the following fields are modified in this routine 836 !!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) ) 837 !!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) ) 838 !!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) ) 839 CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) ) 840 CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 841 CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) ) 842 CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 843 CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 844 DO jl = 1, jpl 845 DO jk = 1, nlay_s 846 CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 847 END DO 848 DO jk = 1, nlay_i 849 CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 850 END DO 851 END DO 852 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_dyn_1d (1:npti), sfx_dyn (:,:) ) 853 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri (:,:) ) 854 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_dyn_1d (1:npti), wfx_dyn (:,:) ) 855 CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d (1:npti), hfx_dyn (:,:) ) 856 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 857 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 858 859 ! !---------------------! 860 CASE( 2 ) !== from 1D to 2D ==! 861 ! !---------------------! 862 CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) ) 863 CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) ) 864 CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) ) 865 CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) ) 866 CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 867 CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) ) 868 CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 869 CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 870 DO jl = 1, jpl 871 DO jk = 1, nlay_s 872 CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 873 END DO 874 DO jk = 1, nlay_i 875 CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 876 END DO 877 END DO 878 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d (1:npti), sfx_dyn (:,:) ) 879 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d (1:npti), sfx_bri (:,:) ) 880 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d (1:npti), wfx_dyn (:,:) ) 881 CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d (1:npti), hfx_dyn (:,:) ) 882 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 883 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 884 ! 885 END SELECT 886 ! 887 END SUBROUTINE ice_dyn_1d2d 888 888 889 889 890 SUBROUTINE ice_dyn_rdgrft_init
Note: See TracChangeset
for help on using the changeset viewer.
