Changeset 6400 for branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO
- Timestamp:
- 2016-03-22T18:27:27+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r5202 r6400 45 45 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: asum ! sum of total ice and open water area 46 46 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: aksum ! ratio of area removed to area ridged 47 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: athorn ! participation function; fraction of ridging/ 48 ! ! closing associated w/ category n 47 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: athorn ! participation function; fraction of ridging/closing associated w/ category n 49 48 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmin ! minimum ridge thickness 50 49 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hrmax ! maximum ridge thickness 51 50 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: hraft ! thickness of rafted ice 52 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: krdg ! mean ridge thickness/thickness of ridging ice51 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: krdg ! thickness of ridging ice / mean ridge thickness 53 52 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: aridge ! participating ice ridging 54 53 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: araft ! participating ice rafting 55 54 56 55 REAL(wp), PARAMETER :: krdgmin = 1.1_wp ! min ridge thickness multiplier 57 REAL(wp), PARAMETER :: kraft = 2.0_wp ! rafting multipliyer 58 REAL(wp), PARAMETER :: kamax = 1.0_wp ! max of ice area authorized (clem: scheme is not stable if kamax <= 0.99) 56 REAL(wp), PARAMETER :: kraft = 0.5_wp ! rafting multipliyer 59 57 60 58 REAL(wp) :: Cp ! 61 59 ! 62 !-----------------------------------------------------------------------63 ! Ridging diagnostic arrays for history files64 !-----------------------------------------------------------------------65 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dardg1dt ! rate of fractional area loss by ridging ice (1/s)66 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dardg2dt ! rate of fractional area gain by new ridges (1/s)67 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: dvirdgdt ! rate of ice volume ridged (m/s)68 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: opening ! rate of opening due to divergence/shear (1/s)69 60 ! 70 61 !!---------------------------------------------------------------------- … … 83 74 & asum (jpi,jpj) , athorn(jpi,jpj,0:jpl) , & 84 75 & aksum(jpi,jpj) , & 85 !86 76 & hrmin(jpi,jpj,jpl) , hraft(jpi,jpj,jpl) , aridge(jpi,jpj,jpl) , & 87 & hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) , & 88 ! 89 !* Ridging diagnostic arrays for history files 90 & dardg1dt(jpi,jpj) , dardg2dt(jpi,jpj) , & 91 & dvirdgdt(jpi,jpj) , opening(jpi,jpj) , STAT=lim_itd_me_alloc ) 77 & hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) , STAT=lim_itd_me_alloc ) 92 78 ! 93 79 IF( lim_itd_me_alloc /= 0 ) CALL ctl_warn( 'lim_itd_me_alloc: failed to allocate arrays' ) … … 132 118 REAL(wp), POINTER, DIMENSION(:,:) :: opning ! rate of opening due to divergence/shear 133 119 REAL(wp), POINTER, DIMENSION(:,:) :: closing_gross ! rate at which area removed, not counting area of new ridges 134 REAL(wp), POINTER, DIMENSION(:,:) :: msnow_mlt ! mass of snow added to ocean (kg m-2)135 REAL(wp), POINTER, DIMENSION(:,:) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2)136 REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories137 120 ! 138 121 INTEGER, PARAMETER :: nitermax = 20 … … 142 125 IF( nn_timing == 1 ) CALL timing_start('limitd_me') 143 126 144 CALL wrk_alloc( jpi,jpj, closing_net, divu_adv, opning, closing_gross , msnow_mlt, esnow_mlt, vt_i_init, vt_i_final)127 CALL wrk_alloc( jpi,jpj, closing_net, divu_adv, opning, closing_gross ) 145 128 146 129 IF(ln_ctl) THEN … … 154 137 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 155 138 156 CALL lim_var_zapsmall157 CALL lim_var_glo2eqv ! equivalent variables, requested for rafting158 159 139 !-----------------------------------------------------------------------------! 160 140 ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons … … 164 144 CALL lim_itd_me_ridgeprep ! prepare ridging 165 145 ! 166 IF( con_i) CALL lim_column_sum( jpl, v_i, vt_i_init ) ! conservation check167 146 168 147 DO jj = 1, jpj ! Initialize arrays. 169 148 DO ji = 1, jpi 170 msnow_mlt(ji,jj) = 0._wp171 esnow_mlt(ji,jj) = 0._wp172 dardg1dt (ji,jj) = 0._wp173 dardg2dt (ji,jj) = 0._wp174 dvirdgdt (ji,jj) = 0._wp175 opening (ji,jj) = 0._wp176 149 177 150 !-----------------------------------------------------------------------------! … … 204 177 ! If divu_adv < 0, make sure the closing rate is large enough 205 178 ! to give asum = 1.0 after ridging. 206 207 divu_adv(ji,jj) = ( kamax- asum(ji,jj) ) * r1_rdtice ! asum found in ridgeprep179 180 divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice ! asum found in ridgeprep 208 181 209 182 IF( divu_adv(ji,jj) < 0._wp ) closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) ) … … 224 197 DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax ) 225 198 199 ! 3.2 closing_gross 200 !-----------------------------------------------------------------------------! 201 ! Based on the ITD of ridging and ridged ice, convert the net 202 ! closing rate to a gross closing rate. 203 ! NOTE: 0 < aksum <= 1 204 closing_gross(:,:) = closing_net(:,:) / aksum(:,:) 205 206 ! correction to closing rate and opening if closing rate is excessive 207 !--------------------------------------------------------------------- 208 ! Reduce the closing rate if more than 100% of the open water 209 ! would be removed. Reduce the opening rate proportionately. 226 210 DO jj = 1, jpj 227 211 DO ji = 1, jpi 228 229 ! 3.2 closing_gross 230 !-----------------------------------------------------------------------------! 231 ! Based on the ITD of ridging and ridged ice, convert the net 232 ! closing rate to a gross closing rate. 233 ! NOTE: 0 < aksum <= 1 234 closing_gross(ji,jj) = closing_net(ji,jj) / aksum(ji,jj) 235 236 ! correction to closing rate and opening if closing rate is excessive 237 !--------------------------------------------------------------------- 238 ! Reduce the closing rate if more than 100% of the open water 239 ! would be removed. Reduce the opening rate proportionately. 240 za = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 241 IF( za > epsi20 ) THEN 242 zfac = MIN( 1._wp, ato_i(ji,jj) / za ) 243 closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 244 opning (ji,jj) = opning (ji,jj) * zfac 212 za = ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice 213 IF( za < 0._wp .AND. za > - ato_i(ji,jj) ) THEN ! would lead to negative ato_i 214 zfac = - ato_i(ji,jj) / za 215 opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice 216 ELSEIF( za > 0._wp .AND. za > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN ! would lead to ato_i > asum 217 zfac = ( asum(ji,jj) - ato_i(ji,jj) ) / za 218 opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice 245 219 ENDIF 246 247 220 END DO 248 221 END DO … … 256 229 DO ji = 1, jpi 257 230 za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 258 IF( za > epsi20) THEN259 zfac = MIN( 1._wp, a_i(ji,jj,jl) / za )231 IF( za > a_i(ji,jj,jl) ) THEN 232 zfac = a_i(ji,jj,jl) / za 260 233 closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 261 opning (ji,jj) = opning (ji,jj) * zfac262 234 ENDIF 263 235 END DO … … 268 240 !-----------------------------------------------------------------------------! 269 241 270 CALL lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt ) 271 242 CALL lim_itd_me_ridgeshift( opning, closing_gross ) 243 244 272 245 ! 3.4 Compute total area of ice plus open water after ridging. 273 246 !-----------------------------------------------------------------------------! 274 247 ! This is in general not equal to one because of divergence during transport 275 asum(:,:) = ato_i(:,:) 276 DO jl = 1, jpl 277 asum(:,:) = asum(:,:) + a_i(:,:,jl) 278 END DO 248 asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 ) 279 249 280 250 ! 3.5 Do we keep on iterating ??? … … 284 254 285 255 iterate_ridging = 0 286 287 256 DO jj = 1, jpj 288 257 DO ji = 1, jpi 289 IF (ABS(asum(ji,jj) - kamax ) < epsi10) THEN258 IF( ABS( asum(ji,jj) - 1._wp ) < epsi10 ) THEN 290 259 closing_net(ji,jj) = 0._wp 291 260 opning (ji,jj) = 0._wp 292 261 ELSE 293 262 iterate_ridging = 1 294 divu_adv (ji,jj) = ( kamax- asum(ji,jj) ) * r1_rdtice263 divu_adv (ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice 295 264 closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) ) 296 265 opning (ji,jj) = MAX( 0._wp, divu_adv(ji,jj) ) … … 309 278 310 279 IF( iterate_ridging == 1 ) THEN 280 CALL lim_itd_me_ridgeprep 311 281 IF( niter > nitermax ) THEN 312 282 WRITE(numout,*) ' ALERTE : non-converging ridging scheme ' 313 283 WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging 314 284 ENDIF 315 CALL lim_itd_me_ridgeprep316 285 ENDIF 317 286 318 287 END DO !! on the do while over iter 319 320 !-----------------------------------------------------------------------------!321 ! 4) Ridging diagnostics322 !-----------------------------------------------------------------------------!323 ! Convert ridging rate diagnostics to correct units.324 ! Update fresh water and heat fluxes due to snow melt.325 DO jj = 1, jpj326 DO ji = 1, jpi327 328 dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice329 dardg2dt(ji,jj) = dardg2dt(ji,jj) * r1_rdtice330 dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * r1_rdtice331 opening (ji,jj) = opening (ji,jj) * r1_rdtice332 333 !-----------------------------------------------------------------------------!334 ! 5) Heat, salt and freshwater fluxes335 !-----------------------------------------------------------------------------!336 wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice ! fresh water source for ocean337 hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice ! heat sink for ocean (<0, W.m-2)338 339 END DO340 END DO341 342 ! Check if there is a ridging error343 IF( lwp ) THEN344 DO jj = 1, jpj345 DO ji = 1, jpi346 IF( ABS( asum(ji,jj) - kamax) > epsi10 ) THEN ! there is a bug347 WRITE(numout,*) ' '348 WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj)349 WRITE(numout,*) ' limitd_me '350 WRITE(numout,*) ' POINT : ', ji, jj351 WRITE(numout,*) ' jpl, a_i, athorn '352 WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0)353 DO jl = 1, jpl354 WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl)355 END DO356 ENDIF357 END DO358 END DO359 END IF360 361 ! Conservation check362 IF ( con_i ) THEN363 CALL lim_column_sum (jpl, v_i, vt_i_final)364 fieldid = ' v_i : limitd_me '365 CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)366 ENDIF367 288 368 289 CALL lim_var_agg( 1 ) … … 410 331 ENDIF ! ln_limdyn=.true. 411 332 ! 412 CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross , msnow_mlt, esnow_mlt, vt_i_init, vt_i_final)333 CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross ) 413 334 ! 414 335 IF( nn_timing == 1 ) CALL timing_stop('limitd_me') 415 336 END SUBROUTINE lim_itd_me 416 337 338 SUBROUTINE lim_itd_me_ridgeprep 339 !!---------------------------------------------------------------------! 340 !! *** ROUTINE lim_itd_me_ridgeprep *** 341 !! 342 !! ** Purpose : preparation for ridging and strength calculations 343 !! 344 !! ** Method : Compute the thickness distribution of the ice and open water 345 !! participating in ridging and of the resulting ridges. 346 !!---------------------------------------------------------------------! 347 INTEGER :: ji,jj, jl ! dummy loop indices 348 REAL(wp) :: Gstari, astari, hrmean, zdummy ! local scalar 349 REAL(wp), POINTER, DIMENSION(:,:,:) :: Gsum ! Gsum(n) = sum of areas in categories 0 to n 350 !------------------------------------------------------------------------------! 351 352 CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 ) 353 354 Gstari = 1.0/rn_gstar 355 astari = 1.0/rn_astar 356 aksum(:,:) = 0.0 357 athorn(:,:,:) = 0.0 358 aridge(:,:,:) = 0.0 359 araft (:,:,:) = 0.0 360 361 ! Zero out categories with very small areas 362 CALL lim_var_zapsmall 363 364 ! Ice thickness needed for rafting 365 DO jl = 1, jpl 366 DO jj = 1, jpj 367 DO ji = 1, jpi 368 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 369 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 370 END DO 371 END DO 372 END DO 373 374 !------------------------------------------------------------------------------! 375 ! 1) Participation function 376 !------------------------------------------------------------------------------! 377 378 ! Compute total area of ice plus open water. 379 ! This is in general not equal to one because of divergence during transport 380 asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 ) 381 382 ! Compute cumulative thickness distribution function 383 ! Compute the cumulative thickness distribution function Gsum, 384 ! where Gsum(n) is the fractional area in categories 0 to n. 385 ! initial value (in h = 0) equals open water area 386 Gsum(:,:,-1) = 0._wp 387 Gsum(:,:,0 ) = ato_i(:,:) 388 ! for each value of h, you have to add ice concentration then 389 DO jl = 1, jpl 390 Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl) 391 END DO 392 393 ! Normalize the cumulative distribution to 1 394 DO jl = 0, jpl 395 Gsum(:,:,jl) = Gsum(:,:,jl) / asum(:,:) 396 END DO 397 398 ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn) 399 !-------------------------------------------------------------------------------------------------- 400 ! Compute the participation function athorn; this is analogous to 401 ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975). 402 ! area lost from category n due to ridging/closing 403 ! athorn(n) = total area lost due to ridging/closing 404 ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar). 405 ! 406 ! The expressions for athorn are found by integrating b(h)g(h) between 407 ! the category boundaries. 408 ! athorn is always >= 0 and SUM(athorn(0:jpl))=1 409 !----------------------------------------------------------------- 410 411 IF( nn_partfun == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975) 412 DO jl = 0, jpl 413 DO jj = 1, jpj 414 DO ji = 1, jpi 415 IF ( Gsum(ji,jj,jl) < rn_gstar ) THEN 416 athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * & 417 & ( 2._wp - ( Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari ) 418 ELSEIF( Gsum(ji,jj,jl-1) < rn_gstar ) THEN 419 athorn(ji,jj,jl) = Gstari * ( rn_gstar - Gsum(ji,jj,jl-1) ) * & 420 & ( 2._wp - ( Gsum(ji,jj,jl-1) + rn_gstar ) * Gstari ) 421 ELSE 422 athorn(ji,jj,jl) = 0._wp 423 ENDIF 424 END DO 425 END DO 426 END DO 427 428 ELSE !--- Exponential, more stable formulation (Lipscomb et al, 2007) 429 ! 430 zdummy = 1._wp / ( 1._wp - EXP(-astari) ) ! precompute exponential terms using Gsum as a work array 431 DO jl = -1, jpl 432 Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy 433 END DO 434 DO jl = 0, jpl 435 athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl) 436 END DO 437 ! 438 ENDIF 439 440 IF( ln_rafting ) THEN ! Ridging and rafting ice participation functions 441 ! 442 DO jl = 1, jpl 443 DO jj = 1, jpj 444 DO ji = 1, jpi 445 zdummy = TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) 446 aridge(ji,jj,jl) = ( 1._wp + zdummy ) * 0.5_wp * athorn(ji,jj,jl) 447 araft (ji,jj,jl) = ( 1._wp - zdummy ) * 0.5_wp * athorn(ji,jj,jl) 448 END DO 449 END DO 450 END DO 451 452 ELSE 453 ! 454 DO jl = 1, jpl 455 aridge(:,:,jl) = athorn(:,:,jl) 456 END DO 457 ! 458 ENDIF 459 460 !----------------------------------------------------------------- 461 ! 2) Transfer function 462 !----------------------------------------------------------------- 463 ! Compute max and min ridged ice thickness for each ridging category. 464 ! Assume ridged ice is uniformly distributed between hrmin and hrmax. 465 ! 466 ! This parameterization is a modified version of Hibler (1980). 467 ! The mean ridging thickness, hrmean, is proportional to hi^(0.5) 468 ! and for very thick ridging ice must be >= krdgmin*hi 469 ! 470 ! The minimum ridging thickness, hrmin, is equal to 2*hi 471 ! (i.e., rafting) and for very thick ridging ice is 472 ! constrained by hrmin <= (hrmean + hi)/2. 473 ! 474 ! The maximum ridging thickness, hrmax, is determined by 475 ! hrmean and hrmin. 476 ! 477 ! These modifications have the effect of reducing the ice strength 478 ! (relative to the Hibler formulation) when very thick ice is 479 ! ridging. 480 ! 481 ! aksum = net area removed/ total area removed 482 ! where total area removed = area of ice that ridges 483 ! net area removed = total area removed - area of new ridges 484 !----------------------------------------------------------------- 485 486 aksum(:,:) = athorn(:,:,0) 487 ! Transfer function 488 DO jl = 1, jpl !all categories have a specific transfer function 489 DO jj = 1, jpj 490 DO ji = 1, jpi 491 492 IF( athorn(ji,jj,jl) > 0._wp ) THEN 493 hrmean = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin ) 494 hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( hrmean + ht_i(ji,jj,jl) ) ) 495 hrmax(ji,jj,jl) = 2._wp * hrmean - hrmin(ji,jj,jl) 496 hraft(ji,jj,jl) = ht_i(ji,jj,jl) / kraft 497 krdg(ji,jj,jl) = ht_i(ji,jj,jl) / MAX( hrmean, epsi20 ) 498 499 ! Normalization factor : aksum, ensures mass conservation 500 aksum(ji,jj) = aksum(ji,jj) + aridge(ji,jj,jl) * ( 1._wp - krdg(ji,jj,jl) ) & 501 & + araft (ji,jj,jl) * ( 1._wp - kraft ) 502 503 ELSE 504 hrmin(ji,jj,jl) = 0._wp 505 hrmax(ji,jj,jl) = 0._wp 506 hraft(ji,jj,jl) = 0._wp 507 krdg (ji,jj,jl) = 1._wp 508 ENDIF 509 510 END DO 511 END DO 512 END DO 513 ! 514 CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 ) 515 ! 516 END SUBROUTINE lim_itd_me_ridgeprep 517 518 519 SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross ) 520 !!---------------------------------------------------------------------- 521 !! *** ROUTINE lim_itd_me_icestrength *** 522 !! 523 !! ** Purpose : shift ridging ice among thickness categories of ice thickness 524 !! 525 !! ** Method : Remove area, volume, and energy from each ridging category 526 !! and add to thicker ice categories. 527 !!---------------------------------------------------------------------- 528 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: opning ! rate of opening due to divergence/shear 529 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: closing_gross ! rate at which area removed, excluding area of new ridges 530 ! 531 CHARACTER (len=80) :: fieldid ! field identifier 532 ! 533 INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices 534 INTEGER :: ij ! horizontal index, combines i and j loops 535 INTEGER :: icells ! number of cells with a_i > puny 536 REAL(wp) :: hL, hR, farea ! left and right limits of integration 537 538 INTEGER , POINTER, DIMENSION(:) :: indxi, indxj ! compressed indices 539 REAL(wp), POINTER, DIMENSION(:) :: zswitch, fvol ! new ridge volume going to n2 540 541 REAL(wp), POINTER, DIMENSION(:) :: afrac ! fraction of category area ridged 542 REAL(wp), POINTER, DIMENSION(:) :: ardg1 , ardg2 ! area of ice ridged & new ridges 543 REAL(wp), POINTER, DIMENSION(:) :: vsrdg , esrdg ! snow volume & energy of ridging ice 544 REAL(wp), POINTER, DIMENSION(:) :: dhr , dhr2 ! hrmax - hrmin & hrmax^2 - hrmin^2 545 546 REAL(wp), POINTER, DIMENSION(:) :: vrdg1 ! volume of ice ridged 547 REAL(wp), POINTER, DIMENSION(:) :: vrdg2 ! volume of new ridges 548 REAL(wp), POINTER, DIMENSION(:) :: vsw ! volume of seawater trapped into ridges 549 REAL(wp), POINTER, DIMENSION(:) :: srdg1 ! sal*volume of ice ridged 550 REAL(wp), POINTER, DIMENSION(:) :: srdg2 ! sal*volume of new ridges 551 REAL(wp), POINTER, DIMENSION(:) :: smsw ! sal*volume of water trapped into ridges 552 REAL(wp), POINTER, DIMENSION(:) :: oirdg1, oirdg2 ! ice age of ice ridged 553 554 REAL(wp), POINTER, DIMENSION(:) :: afrft ! fraction of category area rafted 555 REAL(wp), POINTER, DIMENSION(:) :: arft1 , arft2 ! area of ice rafted and new rafted zone 556 REAL(wp), POINTER, DIMENSION(:) :: virft , vsrft ! ice & snow volume of rafting ice 557 REAL(wp), POINTER, DIMENSION(:) :: esrft , smrft ! snow energy & salinity of rafting ice 558 REAL(wp), POINTER, DIMENSION(:) :: oirft1, oirft2 ! ice age of ice rafted 559 560 REAL(wp), POINTER, DIMENSION(:,:) :: eirft ! ice energy of rafting ice 561 REAL(wp), POINTER, DIMENSION(:,:) :: erdg1 ! enth*volume of ice ridged 562 REAL(wp), POINTER, DIMENSION(:,:) :: erdg2 ! enth*volume of new ridges 563 REAL(wp), POINTER, DIMENSION(:,:) :: ersw ! enth of water trapped into ridges 564 !!---------------------------------------------------------------------- 565 566 CALL wrk_alloc( jpij, indxi, indxj ) 567 CALL wrk_alloc( jpij, zswitch, fvol ) 568 CALL wrk_alloc( jpij, afrac, ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 ) 569 CALL wrk_alloc( jpij, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw, oirdg1, oirdg2 ) 570 CALL wrk_alloc( jpij, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 571 CALL wrk_alloc( jpij,nlay_i, eirft, erdg1, erdg2, ersw ) 572 573 !------------------------------------------------------------------------------- 574 ! 1) Compute change in open water area due to closing and opening. 575 !------------------------------------------------------------------------------- 576 DO jj = 1, jpj 577 DO ji = 1, jpi 578 ato_i(ji,jj) = MAX( 0._wp, ato_i(ji,jj) + & 579 & ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice ) 580 END DO 581 END DO 582 583 !----------------------------------------------------------------- 584 ! 3) Pump everything from ice which is being ridged / rafted 585 !----------------------------------------------------------------- 586 ! Compute the area, volume, and energy of ice ridging in each 587 ! category, along with the area of the resulting ridge. 588 589 DO jl1 = 1, jpl !jl1 describes the ridging category 590 591 !------------------------------------------------ 592 ! 3.1) Identify grid cells with nonzero ridging 593 !------------------------------------------------ 594 icells = 0 595 DO jj = 1, jpj 596 DO ji = 1, jpi 597 IF( athorn(ji,jj,jl1) > 0._wp .AND. closing_gross(ji,jj) > 0._wp ) THEN 598 icells = icells + 1 599 indxi(icells) = ji 600 indxj(icells) = jj 601 ENDIF 602 END DO 603 END DO 604 605 DO ij = 1, icells 606 ji = indxi(ij) ; jj = indxj(ij) 607 608 !-------------------------------------------------------------------- 609 ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2) 610 !-------------------------------------------------------------------- 611 ardg1(ij) = aridge(ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice 612 arft1(ij) = araft (ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice 613 614 !--------------------------------------------------------------- 615 ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1 616 !--------------------------------------------------------------- 617 afrac(ij) = ardg1(ij) / a_i(ji,jj,jl1) !ridging 618 afrft(ij) = arft1(ij) / a_i(ji,jj,jl1) !rafting 619 ardg2(ij) = ardg1(ij) * krdg(ji,jj,jl1) 620 arft2(ij) = arft1(ij) * kraft 621 622 !-------------------------------------------------------------------------- 623 ! 3.4) Subtract area, volume, and energy from ridging 624 ! / rafting category n1. 625 !-------------------------------------------------------------------------- 626 vrdg1(ij) = v_i(ji,jj,jl1) * afrac(ij) 627 vrdg2(ij) = vrdg1(ij) * ( 1. + rn_por_rdg ) 628 vsw (ij) = vrdg1(ij) * rn_por_rdg 629 630 vsrdg (ij) = v_s (ji,jj, jl1) * afrac(ij) 631 esrdg (ij) = e_s (ji,jj,1,jl1) * afrac(ij) 632 srdg1 (ij) = smv_i(ji,jj, jl1) * afrac(ij) 633 oirdg1(ij) = oa_i (ji,jj, jl1) * afrac(ij) 634 oirdg2(ij) = oa_i (ji,jj, jl1) * afrac(ij) * krdg(ji,jj,jl1) 635 636 ! rafting volumes, heat contents ... 637 virft (ij) = v_i (ji,jj, jl1) * afrft(ij) 638 vsrft (ij) = v_s (ji,jj, jl1) * afrft(ij) 639 esrft (ij) = e_s (ji,jj,1,jl1) * afrft(ij) 640 smrft (ij) = smv_i(ji,jj, jl1) * afrft(ij) 641 oirft1(ij) = oa_i (ji,jj, jl1) * afrft(ij) 642 oirft2(ij) = oa_i (ji,jj, jl1) * afrft(ij) * kraft 643 644 !----------------------------------------------------------------- 645 ! 3.5) Compute properties of new ridges 646 !----------------------------------------------------------------- 647 smsw(ij) = vsw(ij) * sss_m(ji,jj) ! salt content of seawater frozen in voids 648 srdg2(ij) = srdg1(ij) + smsw(ij) ! salt content of new ridge 649 650 sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ij) * rhoic * r1_rdtice 651 wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ij) * rhoic * r1_rdtice ! increase in ice volume due to seawater frozen in voids 652 653 !------------------------------------------ 654 ! 3.7 Put the snow somewhere in the ocean 655 !------------------------------------------ 656 ! Place part of the snow lost by ridging into the ocean. 657 ! Note that esrdg > 0; the ocean must cool to melt snow. 658 ! If the ocean temp = Tf already, new ice must grow. 659 ! During the next time step, thermo_rates will determine whether 660 ! the ocean cools or new ice grows. 661 wfx_snw(ji,jj) = wfx_snw(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg ) & 662 & + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice ! fresh water source for ocean 663 664 hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg ) & 665 & - esrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice ! heat sink for ocean (<0, W.m-2) 666 667 !----------------------------------------------------------------- 668 ! 3.8 Compute quantities used to apportion ice among categories 669 ! in the n2 loop below 670 !----------------------------------------------------------------- 671 dhr (ij) = 1._wp / ( hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) ) 672 dhr2(ij) = 1._wp / ( hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) ) 673 674 675 ! update jl1 (removing ridged/rafted area) 676 a_i (ji,jj, jl1) = a_i (ji,jj, jl1) - ardg1 (ij) - arft1 (ij) 677 v_i (ji,jj, jl1) = v_i (ji,jj, jl1) - vrdg1 (ij) - virft (ij) 678 v_s (ji,jj, jl1) = v_s (ji,jj, jl1) - vsrdg (ij) - vsrft (ij) 679 e_s (ji,jj,1,jl1) = e_s (ji,jj,1,jl1) - esrdg (ij) - esrft (ij) 680 smv_i(ji,jj, jl1) = smv_i(ji,jj, jl1) - srdg1 (ij) - smrft (ij) 681 oa_i (ji,jj, jl1) = oa_i (ji,jj, jl1) - oirdg1(ij) - oirft1(ij) 682 683 END DO 684 685 !-------------------------------------------------------------------- 686 ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and 687 ! compute ridged ice enthalpy 688 !-------------------------------------------------------------------- 689 DO jk = 1, nlay_i 690 DO ij = 1, icells 691 ji = indxi(ij) ; jj = indxj(ij) 692 ! heat content of ridged ice 693 erdg1(ij,jk) = e_i(ji,jj,jk,jl1) * afrac(ij) 694 eirft(ij,jk) = e_i(ji,jj,jk,jl1) * afrft(ij) 695 696 ! enthalpy of the trapped seawater (J/m2, >0) 697 ! clem: if sst>0, then ersw <0 (is that possible?) 698 ersw(ij,jk) = - rhoic * vsw(ij) * rcp * sst_m(ji,jj) * r1_nlay_i 699 700 ! heat flux to the ocean 701 hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ij,jk) * r1_rdtice ! > 0 [W.m-2] ocean->ice flux 702 703 ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean 704 erdg2(ij,jk) = erdg1(ij,jk) + ersw(ij,jk) 705 706 ! update jl1 707 e_i (ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk) 708 709 END DO 710 END DO 711 712 !------------------------------------------------------------------------------- 713 ! 4) Add area, volume, and energy of new ridge to each category jl2 714 !------------------------------------------------------------------------------- 715 DO jl2 = 1, jpl 716 ! over categories to which ridged/rafted ice is transferred 717 DO ij = 1, icells 718 ji = indxi(ij) ; jj = indxj(ij) 719 720 ! Compute the fraction of ridged ice area and volume going to thickness category jl2. 721 IF( hrmin(ji,jj,jl1) <= hi_max(jl2) .AND. hrmax(ji,jj,jl1) > hi_max(jl2-1) ) THEN 722 hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) ) 723 hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2) ) 724 farea = ( hR - hL ) * dhr(ij) 725 fvol(ij) = ( hR * hR - hL * hL ) * dhr2(ij) 726 ELSE 727 farea = 0._wp 728 fvol(ij) = 0._wp 729 ENDIF 730 731 ! Compute the fraction of rafted ice area and volume going to thickness category jl2 732 IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN 733 zswitch(ij) = 1._wp 734 ELSE 735 zswitch(ij) = 0._wp 736 ENDIF 737 738 a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + ( ardg2 (ij) * farea + arft2 (ij) * zswitch(ij) ) 739 oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + ( oirdg2(ij) * farea + oirft2(ij) * zswitch(ij) ) 740 v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + ( vrdg2 (ij) * fvol(ij) + virft (ij) * zswitch(ij) ) 741 smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + ( srdg2 (ij) * fvol(ij) + smrft (ij) * zswitch(ij) ) 742 v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + ( vsrdg (ij) * rn_fsnowrdg * fvol(ij) + & 743 & vsrft (ij) * rn_fsnowrft * zswitch(ij) ) 744 e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnowrdg * fvol(ij) + & 745 & esrft (ij) * rn_fsnowrft * zswitch(ij) ) 746 747 END DO 748 749 ! Transfer ice energy to category jl2 by ridging 750 DO jk = 1, nlay_i 751 DO ij = 1, icells 752 ji = indxi(ij) ; jj = indxj(ij) 753 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + erdg2(ij,jk) * fvol(ij) + eirft(ij,jk) * zswitch(ij) 754 END DO 755 END DO 756 ! 757 END DO ! jl2 758 759 END DO ! jl1 (deforming categories) 760 761 ! 762 CALL wrk_dealloc( jpij, indxi, indxj ) 763 CALL wrk_dealloc( jpij, zswitch, fvol ) 764 CALL wrk_dealloc( jpij, afrac, ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 ) 765 CALL wrk_dealloc( jpij, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw, oirdg1, oirdg2 ) 766 CALL wrk_dealloc( jpij, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 767 CALL wrk_dealloc( jpij,nlay_i, eirft, erdg1, erdg2, ersw ) 768 ! 769 END SUBROUTINE lim_itd_me_ridgeshift 417 770 418 771 SUBROUTINE lim_itd_me_icestrength( kstrngth ) … … 434 787 INTEGER :: ksmooth ! smoothing the resistance to deformation 435 788 INTEGER :: numts_rm ! number of time steps for the P smoothing 436 REAL(wp) :: z hi, zp, z1_3! local scalars789 REAL(wp) :: zp, z1_3 ! local scalars 437 790 REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here 438 791 !!---------------------------------------------------------------------- … … 459 812 DO ji = 1, jpi 460 813 ! 461 IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN 462 zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 814 IF( athorn(ji,jj,jl) > 0._wp ) THEN 463 815 !---------------------------- 464 816 ! PE loss from deforming ice 465 817 !---------------------------- 466 strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * zhi * zhi818 strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl) 467 819 468 820 !-------------------------- 469 821 ! PE gain from rafting ice 470 822 !-------------------------- 471 strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * zhi * zhi823 strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl) 472 824 473 825 !---------------------------- 474 826 ! PE gain from ridging ice 475 827 !---------------------------- 476 strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) / krdg(ji,jj,jl) & 477 * z1_3 * ( hrmax(ji,jj,jl)**2 + hrmin(ji,jj,jl)**2 + hrmax(ji,jj,jl) * hrmin(ji,jj,jl) ) 828 strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) * krdg(ji,jj,jl) * z1_3 * & 829 & ( hrmax(ji,jj,jl) * hrmax(ji,jj,jl) + & 830 & hrmin(ji,jj,jl) * hrmin(ji,jj,jl) + & 831 & hrmax(ji,jj,jl) * hrmin(ji,jj,jl) ) 478 832 !!(a**3-b**3)/(a-b) = a*a+ab+b*b 479 833 ENDIF … … 497 851 ! 498 852 ENDIF ! kstrngth 499 500 853 ! 501 854 !------------------------------------------------------------------------------! … … 503 856 !------------------------------------------------------------------------------! 504 857 ! CAN BE REMOVED 505 !506 858 IF( ln_icestr_bvf ) THEN 507 508 859 DO jj = 1, jpj 509 860 DO ji = 1, jpi … … 511 862 END DO 512 863 END DO 513 514 864 ENDIF 515 516 865 ! 517 866 !------------------------------------------------------------------------------! … … 558 907 IF ( ksmooth == 2 ) THEN 559 908 560 561 909 CALL lbc_lnk( strength, 'T', 1. ) 562 910 … … 565 913 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN 566 914 numts_rm = 1 ! number of time steps for the running mean 567 IF ( strp1(ji,jj) > 0. 0) numts_rm = numts_rm + 1568 IF ( strp2(ji,jj) > 0. 0) numts_rm = numts_rm + 1915 IF ( strp1(ji,jj) > 0._wp ) numts_rm = numts_rm + 1 916 IF ( strp2(ji,jj) > 0._wp ) numts_rm = numts_rm + 1 569 917 zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm 570 918 strp2(ji,jj) = strp1(ji,jj) … … 583 931 ! 584 932 END SUBROUTINE lim_itd_me_icestrength 585 586 587 SUBROUTINE lim_itd_me_ridgeprep588 !!---------------------------------------------------------------------!589 !! *** ROUTINE lim_itd_me_ridgeprep ***590 !!591 !! ** Purpose : preparation for ridging and strength calculations592 !!593 !! ** Method : Compute the thickness distribution of the ice and open water594 !! participating in ridging and of the resulting ridges.595 !!---------------------------------------------------------------------!596 INTEGER :: ji,jj, jl ! dummy loop indices597 REAL(wp) :: Gstari, astari, zhi, hrmean, zdummy ! local scalar598 REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here599 REAL(wp), POINTER, DIMENSION(:,:,:) :: Gsum ! Gsum(n) = sum of areas in categories 0 to n600 !------------------------------------------------------------------------------!601 602 CALL wrk_alloc( jpi,jpj, zworka )603 CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )604 605 Gstari = 1.0/rn_gstar606 astari = 1.0/rn_astar607 aksum(:,:) = 0.0608 athorn(:,:,:) = 0.0609 aridge(:,:,:) = 0.0610 araft (:,:,:) = 0.0611 hrmin(:,:,:) = 0.0612 hrmax(:,:,:) = 0.0613 hraft(:,:,:) = 0.0614 krdg (:,:,:) = 1.0615 616 ! ! Zero out categories with very small areas617 CALL lim_var_zapsmall618 619 !------------------------------------------------------------------------------!620 ! 1) Participation function621 !------------------------------------------------------------------------------!622 623 ! Compute total area of ice plus open water.624 ! This is in general not equal to one because of divergence during transport625 asum(:,:) = ato_i(:,:)626 DO jl = 1, jpl627 asum(:,:) = asum(:,:) + a_i(:,:,jl)628 END DO629 630 ! Compute cumulative thickness distribution function631 ! Compute the cumulative thickness distribution function Gsum,632 ! where Gsum(n) is the fractional area in categories 0 to n.633 ! initial value (in h = 0) equals open water area634 635 Gsum(:,:,-1) = 0._wp636 Gsum(:,:,0 ) = ato_i(:,:)637 638 ! for each value of h, you have to add ice concentration then639 DO jl = 1, jpl640 Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl)641 END DO642 643 ! Normalize the cumulative distribution to 1644 zworka(:,:) = 1._wp / Gsum(:,:,jpl)645 DO jl = 0, jpl646 Gsum(:,:,jl) = Gsum(:,:,jl) * zworka(:,:)647 END DO648 649 ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)650 !--------------------------------------------------------------------------------------------------651 ! Compute the participation function athorn; this is analogous to652 ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).653 ! area lost from category n due to ridging/closing654 ! athorn(n) = total area lost due to ridging/closing655 ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).656 !657 ! The expressions for athorn are found by integrating b(h)g(h) between658 ! the category boundaries.659 !-----------------------------------------------------------------660 661 IF( nn_partfun == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975)662 DO jl = 0, jpl663 DO jj = 1, jpj664 DO ji = 1, jpi665 IF( Gsum(ji,jj,jl) < rn_gstar) THEN666 athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * &667 & ( 2.0 - (Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari )668 ELSEIF (Gsum(ji,jj,jl-1) < rn_gstar) THEN669 athorn(ji,jj,jl) = Gstari * ( rn_gstar - Gsum(ji,jj,jl-1) ) * &670 & ( 2.0 - ( Gsum(ji,jj,jl-1) + rn_gstar ) * Gstari )671 ELSE672 athorn(ji,jj,jl) = 0.0673 ENDIF674 END DO675 END DO676 END DO677 678 ELSE !--- Exponential, more stable formulation (Lipscomb et al, 2007)679 !680 zdummy = 1._wp / ( 1._wp - EXP(-astari) ) ! precompute exponential terms using Gsum as a work array681 DO jl = -1, jpl682 Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy683 END DO684 DO jl = 0, jpl685 athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)686 END DO687 !688 ENDIF689 690 IF( ln_rafting ) THEN ! Ridging and rafting ice participation functions691 !692 DO jl = 1, jpl693 DO jj = 1, jpj694 DO ji = 1, jpi695 IF ( athorn(ji,jj,jl) > 0._wp ) THEN696 !!gm TANH( -X ) = - TANH( X ) so can be computed only 1 time....697 aridge(ji,jj,jl) = ( TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)698 araft (ji,jj,jl) = ( TANH ( -rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)699 IF ( araft(ji,jj,jl) < epsi06 ) araft(ji,jj,jl) = 0._wp700 aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 )701 ENDIF702 END DO703 END DO704 END DO705 706 ELSE707 !708 DO jl = 1, jpl709 aridge(:,:,jl) = athorn(:,:,jl)710 END DO711 !712 ENDIF713 714 IF( ln_rafting ) THEN715 716 IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 .AND. lwp ) THEN717 DO jl = 1, jpl718 DO jj = 1, jpj719 DO ji = 1, jpi720 IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) > epsi10 ) THEN721 WRITE(numout,*) ' ALERTE 96 : wrong participation function ... '722 WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl723 WRITE(numout,*) ' lat, lon : ', gphit(ji,jj), glamt(ji,jj)724 WRITE(numout,*) ' aridge : ', aridge(ji,jj,1:jpl)725 WRITE(numout,*) ' araft : ', araft(ji,jj,1:jpl)726 WRITE(numout,*) ' athorn : ', athorn(ji,jj,1:jpl)727 ENDIF728 END DO729 END DO730 END DO731 ENDIF732 733 ENDIF734 735 !-----------------------------------------------------------------736 ! 2) Transfer function737 !-----------------------------------------------------------------738 ! Compute max and min ridged ice thickness for each ridging category.739 ! Assume ridged ice is uniformly distributed between hrmin and hrmax.740 !741 ! This parameterization is a modified version of Hibler (1980).742 ! The mean ridging thickness, hrmean, is proportional to hi^(0.5)743 ! and for very thick ridging ice must be >= krdgmin*hi744 !745 ! The minimum ridging thickness, hrmin, is equal to 2*hi746 ! (i.e., rafting) and for very thick ridging ice is747 ! constrained by hrmin <= (hrmean + hi)/2.748 !749 ! The maximum ridging thickness, hrmax, is determined by750 ! hrmean and hrmin.751 !752 ! These modifications have the effect of reducing the ice strength753 ! (relative to the Hibler formulation) when very thick ice is754 ! ridging.755 !756 ! aksum = net area removed/ total area removed757 ! where total area removed = area of ice that ridges758 ! net area removed = total area removed - area of new ridges759 !-----------------------------------------------------------------760 761 ! Transfer function762 DO jl = 1, jpl !all categories have a specific transfer function763 DO jj = 1, jpj764 DO ji = 1, jpi765 766 IF (a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0.0 ) THEN767 zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)768 hrmean = MAX(SQRT(rn_hstar*zhi), zhi*krdgmin)769 hrmin(ji,jj,jl) = MIN(2.0*zhi, 0.5*(hrmean + zhi))770 hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl)771 hraft(ji,jj,jl) = kraft*zhi772 krdg(ji,jj,jl) = hrmean / zhi773 ELSE774 hraft(ji,jj,jl) = 0.0775 hrmin(ji,jj,jl) = 0.0776 hrmax(ji,jj,jl) = 0.0777 krdg (ji,jj,jl) = 1.0778 ENDIF779 780 END DO781 END DO782 END DO783 784 ! Normalization factor : aksum, ensures mass conservation785 aksum(:,:) = athorn(:,:,0)786 DO jl = 1, jpl787 aksum(:,:) = aksum(:,:) + aridge(:,:,jl) * ( 1._wp - 1._wp / krdg(:,:,jl) ) &788 & + araft (:,:,jl) * ( 1._wp - 1._wp / kraft )789 END DO790 !791 CALL wrk_dealloc( jpi,jpj, zworka )792 CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )793 !794 END SUBROUTINE lim_itd_me_ridgeprep795 796 797 SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )798 !!----------------------------------------------------------------------799 !! *** ROUTINE lim_itd_me_icestrength ***800 !!801 !! ** Purpose : shift ridging ice among thickness categories of ice thickness802 !!803 !! ** Method : Remove area, volume, and energy from each ridging category804 !! and add to thicker ice categories.805 !!----------------------------------------------------------------------806 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: opning ! rate of opening due to divergence/shear807 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: closing_gross ! rate at which area removed, excluding area of new ridges808 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: msnow_mlt ! mass of snow added to ocean (kg m-2)809 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2)810 !811 CHARACTER (len=80) :: fieldid ! field identifier812 LOGICAL, PARAMETER :: l_conservation_check = .true. ! if true, check conservation (useful for debugging)813 !814 INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices815 INTEGER :: ij ! horizontal index, combines i and j loops816 INTEGER :: icells ! number of cells with aicen > puny817 REAL(wp) :: hL, hR, farea, ztmelts ! left and right limits of integration818 819 INTEGER , POINTER, DIMENSION(:) :: indxi, indxj ! compressed indices820 821 REAL(wp), POINTER, DIMENSION(:,:) :: vice_init, vice_final ! ice volume summed over categories822 REAL(wp), POINTER, DIMENSION(:,:) :: eice_init, eice_final ! ice energy summed over layers823 824 REAL(wp), POINTER, DIMENSION(:,:,:) :: aicen_init, vicen_init ! ice area & volume before ridging825 REAL(wp), POINTER, DIMENSION(:,:,:) :: vsnwn_init, esnwn_init ! snow volume & energy before ridging826 REAL(wp), POINTER, DIMENSION(:,:,:) :: smv_i_init, oa_i_init ! ice salinity & age before ridging827 828 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: eicen_init ! ice energy before ridging829 830 REAL(wp), POINTER, DIMENSION(:,:) :: afrac , fvol ! fraction of category area ridged & new ridge volume going to n2831 REAL(wp), POINTER, DIMENSION(:,:) :: ardg1 , ardg2 ! area of ice ridged & new ridges832 REAL(wp), POINTER, DIMENSION(:,:) :: vsrdg , esrdg ! snow volume & energy of ridging ice833 REAL(wp), POINTER, DIMENSION(:,:) :: dhr , dhr2 ! hrmax - hrmin & hrmax^2 - hrmin^2834 835 REAL(wp), POINTER, DIMENSION(:,:) :: vrdg1 ! volume of ice ridged836 REAL(wp), POINTER, DIMENSION(:,:) :: vrdg2 ! volume of new ridges837 REAL(wp), POINTER, DIMENSION(:,:) :: vsw ! volume of seawater trapped into ridges838 REAL(wp), POINTER, DIMENSION(:,:) :: srdg1 ! sal*volume of ice ridged839 REAL(wp), POINTER, DIMENSION(:,:) :: srdg2 ! sal*volume of new ridges840 REAL(wp), POINTER, DIMENSION(:,:) :: smsw ! sal*volume of water trapped into ridges841 REAL(wp), POINTER, DIMENSION(:,:) :: oirdg1, oirdg2 ! ice age of ice ridged842 843 REAL(wp), POINTER, DIMENSION(:,:) :: afrft ! fraction of category area rafted844 REAL(wp), POINTER, DIMENSION(:,:) :: arft1 , arft2 ! area of ice rafted and new rafted zone845 REAL(wp), POINTER, DIMENSION(:,:) :: virft , vsrft ! ice & snow volume of rafting ice846 REAL(wp), POINTER, DIMENSION(:,:) :: esrft , smrft ! snow energy & salinity of rafting ice847 REAL(wp), POINTER, DIMENSION(:,:) :: oirft1, oirft2 ! ice age of ice rafted848 849 REAL(wp), POINTER, DIMENSION(:,:,:) :: eirft ! ice energy of rafting ice850 REAL(wp), POINTER, DIMENSION(:,:,:) :: erdg1 ! enth*volume of ice ridged851 REAL(wp), POINTER, DIMENSION(:,:,:) :: erdg2 ! enth*volume of new ridges852 REAL(wp), POINTER, DIMENSION(:,:,:) :: ersw ! enth of water trapped into ridges853 !!----------------------------------------------------------------------854 855 CALL wrk_alloc( (jpi+1)*(jpj+1), indxi, indxj )856 CALL wrk_alloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final )857 CALL wrk_alloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 )858 CALL wrk_alloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw, oirdg1, oirdg2 )859 CALL wrk_alloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )860 CALL wrk_alloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )861 CALL wrk_alloc( jpi, jpj, nlay_i, eirft, erdg1, erdg2, ersw )862 CALL wrk_alloc( jpi, jpj, nlay_i, jpl, eicen_init )863 864 ! Conservation check865 eice_init(:,:) = 0._wp866 867 IF( con_i ) THEN868 CALL lim_column_sum (jpl, v_i, vice_init )869 CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_init )870 DO ji = mi0(iiceprt), mi1(iiceprt)871 DO jj = mj0(jiceprt), mj1(jiceprt)872 WRITE(numout,*) ' vice_init : ', vice_init(ji,jj)873 WRITE(numout,*) ' eice_init : ', eice_init(ji,jj)874 END DO875 END DO876 ENDIF877 878 !-------------------------------------------------------------------------------879 ! 1) Compute change in open water area due to closing and opening.880 !-------------------------------------------------------------------------------881 DO jj = 1, jpj882 DO ji = 1, jpi883 ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice &884 & + opning(ji,jj) * rdt_ice885 IF ( ato_i(ji,jj) < -epsi10 ) THEN ! there is a bug886 IF(lwp) WRITE(numout,*) 'Ridging error: ato_i < 0 -- ato_i : ',ato_i(ji,jj)887 ELSEIF( ato_i(ji,jj) < 0._wp ) THEN ! roundoff error888 ato_i(ji,jj) = 0._wp889 ENDIF890 END DO891 END DO892 893 !-----------------------------------------------------------------894 ! 2) Save initial state variables895 !-----------------------------------------------------------------896 aicen_init(:,:,:) = a_i (:,:,:)897 vicen_init(:,:,:) = v_i (:,:,:)898 vsnwn_init(:,:,:) = v_s (:,:,:)899 smv_i_init(:,:,:) = smv_i(:,:,:)900 esnwn_init(:,:,:) = e_s (:,:,1,:)901 eicen_init(:,:,:,:) = e_i (:,:,:,:)902 oa_i_init (:,:,:) = oa_i (:,:,:)903 904 !905 !-----------------------------------------------------------------906 ! 3) Pump everything from ice which is being ridged / rafted907 !-----------------------------------------------------------------908 ! Compute the area, volume, and energy of ice ridging in each909 ! category, along with the area of the resulting ridge.910 911 DO jl1 = 1, jpl !jl1 describes the ridging category912 913 !------------------------------------------------914 ! 3.1) Identify grid cells with nonzero ridging915 !------------------------------------------------916 917 icells = 0918 DO jj = 1, jpj919 DO ji = 1, jpi920 IF( aicen_init(ji,jj,jl1) > epsi10 .AND. athorn(ji,jj,jl1) > 0._wp &921 & .AND. closing_gross(ji,jj) > 0._wp ) THEN922 icells = icells + 1923 indxi(icells) = ji924 indxj(icells) = jj925 ENDIF926 END DO927 END DO928 929 DO ij = 1, icells930 ji = indxi(ij)931 jj = indxj(ij)932 933 !--------------------------------------------------------------------934 ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)935 !--------------------------------------------------------------------936 937 ardg1(ji,jj) = aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice938 arft1(ji,jj) = araft (ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice939 ardg2(ji,jj) = ardg1(ji,jj) / krdg(ji,jj,jl1)940 arft2(ji,jj) = arft1(ji,jj) / kraft941 942 !---------------------------------------------------------------943 ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1944 !---------------------------------------------------------------945 946 afrac(ji,jj) = ardg1(ji,jj) / aicen_init(ji,jj,jl1) !ridging947 afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting948 949 IF( afrac(ji,jj) > kamax + epsi10 ) THEN ! there is a bug950 IF(lwp) WRITE(numout,*) ' ardg > a_i -- ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1)951 ELSEIF( afrac(ji,jj) > kamax ) THEN ! roundoff error952 afrac(ji,jj) = kamax953 ENDIF954 955 IF( afrft(ji,jj) > kamax + epsi10 ) THEN ! there is a bug956 IF(lwp) WRITE(numout,*) ' arft > a_i -- arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1)957 ELSEIF( afrft(ji,jj) > kamax) THEN ! roundoff error958 afrft(ji,jj) = kamax959 ENDIF960 961 !--------------------------------------------------------------------------962 ! 3.4) Subtract area, volume, and energy from ridging963 ! / rafting category n1.964 !--------------------------------------------------------------------------965 vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj)966 vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + rn_por_rdg )967 vsw (ji,jj) = vrdg1(ji,jj) * rn_por_rdg968 969 vsrdg (ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj)970 esrdg (ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj)971 srdg1 (ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj)972 oirdg1(ji,jj) = oa_i_init (ji,jj,jl1) * afrac(ji,jj)973 oirdg2(ji,jj) = oa_i_init (ji,jj,jl1) * afrac(ji,jj) / krdg(ji,jj,jl1)974 975 ! rafting volumes, heat contents ...976 virft (ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj)977 vsrft (ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj)978 esrft (ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj)979 smrft (ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj)980 oirft1(ji,jj) = oa_i_init (ji,jj,jl1) * afrft(ji,jj)981 oirft2(ji,jj) = oa_i_init (ji,jj,jl1) * afrft(ji,jj) / kraft982 983 ! substract everything984 a_i(ji,jj,jl1) = a_i(ji,jj,jl1) - ardg1 (ji,jj) - arft1 (ji,jj)985 v_i(ji,jj,jl1) = v_i(ji,jj,jl1) - vrdg1 (ji,jj) - virft (ji,jj)986 v_s(ji,jj,jl1) = v_s(ji,jj,jl1) - vsrdg (ji,jj) - vsrft (ji,jj)987 e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg (ji,jj) - esrft (ji,jj)988 smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1 (ji,jj) - smrft (ji,jj)989 oa_i(ji,jj,jl1) = oa_i(ji,jj,jl1) - oirdg1(ji,jj) - oirft1(ji,jj)990 991 !-----------------------------------------------------------------992 ! 3.5) Compute properties of new ridges993 !-----------------------------------------------------------------994 !---------995 ! Salinity996 !---------997 smsw(ji,jj) = vsw(ji,jj) * sss_m(ji,jj) ! salt content of seawater frozen in voids !! MV HC2014998 srdg2(ji,jj) = srdg1(ji,jj) + smsw(ji,jj) ! salt content of new ridge999 1000 !srdg2(ji,jj) = MIN( rn_simax * vrdg2(ji,jj) , zsrdg2 ) ! impose a maximum salinity1001 1002 sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice1003 wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice ! increase in ice volume du to seawater frozen in voids1004 1005 !------------------------------------1006 ! 3.6 Increment ridging diagnostics1007 !------------------------------------1008 1009 ! jl1 looping 1-jpl1010 ! ij looping 1-icells1011 1012 dardg1dt (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj)1013 dardg2dt (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj)1014 opening (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice1015 1016 IF( con_i ) vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj)1017 1018 !------------------------------------------1019 ! 3.7 Put the snow somewhere in the ocean1020 !------------------------------------------1021 ! Place part of the snow lost by ridging into the ocean.1022 ! Note that esnow_mlt < 0; the ocean must cool to melt snow.1023 ! If the ocean temp = Tf already, new ice must grow.1024 ! During the next time step, thermo_rates will determine whether1025 ! the ocean cools or new ice grows.1026 ! jl1 looping 1-jpl1027 ! ij looping 1-icells1028 1029 msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosn*vsrdg(ji,jj)*(1.0-rn_fsnowrdg) & ! rafting included1030 & + rhosn*vsrft(ji,jj)*(1.0-rn_fsnowrft)1031 1032 ! in J/m2 (same as e_s)1033 esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-rn_fsnowrdg) & !rafting included1034 & - esrft(ji,jj)*(1.0-rn_fsnowrft)1035 1036 !-----------------------------------------------------------------1037 ! 3.8 Compute quantities used to apportion ice among categories1038 ! in the n2 loop below1039 !-----------------------------------------------------------------1040 1041 ! jl1 looping 1-jpl1042 ! ij looping 1-icells1043 1044 dhr (ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1)1045 dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1)1046 1047 END DO1048 1049 !--------------------------------------------------------------------1050 ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and1051 ! compute ridged ice enthalpy1052 !--------------------------------------------------------------------1053 DO jk = 1, nlay_i1054 DO ij = 1, icells1055 ji = indxi(ij)1056 jj = indxj(ij)1057 ! heat content of ridged ice1058 erdg1(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj)1059 eirft(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj)1060 e_i (ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ji,jj,jk) - eirft(ji,jj,jk)1061 1062 1063 ! enthalpy of the trapped seawater (J/m2, >0)1064 ! clem: if sst>0, then ersw <0 (is that possible?)1065 ersw(ji,jj,jk) = - rhoic * vsw(ji,jj) * rcp * sst_m(ji,jj) * r1_nlay_i1066 1067 ! heat flux to the ocean1068 hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice ! > 0 [W.m-2] ocean->ice flux1069 1070 ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean1071 erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk)1072 1073 END DO1074 END DO1075 1076 1077 IF( con_i ) THEN1078 DO jk = 1, nlay_i1079 DO ij = 1, icells1080 ji = indxi(ij)1081 jj = indxj(ij)1082 eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - erdg1(ji,jj,jk)1083 END DO1084 END DO1085 ENDIF1086 1087 !-------------------------------------------------------------------------------1088 ! 4) Add area, volume, and energy of new ridge to each category jl21089 !-------------------------------------------------------------------------------1090 ! jl1 looping 1-jpl1091 DO jl2 = 1, jpl1092 ! over categories to which ridged ice is transferred1093 DO ij = 1, icells1094 ji = indxi(ij)1095 jj = indxj(ij)1096 1097 ! Compute the fraction of ridged ice area and volume going to1098 ! thickness category jl2.1099 ! Transfer area, volume, and energy accordingly.1100 1101 IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR. hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN1102 hL = 0._wp1103 hR = 0._wp1104 ELSE1105 hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )1106 hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2) )1107 ENDIF1108 1109 ! fraction of ridged ice area and volume going to n21110 farea = ( hR - hL ) / dhr(ji,jj)1111 fvol(ji,jj) = ( hR*hR - hL*hL ) / dhr2(ji,jj)1112 1113 a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + ardg2 (ji,jj) * farea1114 v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + vrdg2 (ji,jj) * fvol(ji,jj)1115 v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg1116 e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg1117 smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + srdg2 (ji,jj) * fvol(ji,jj)1118 oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirdg2(ji,jj) * farea1119 1120 END DO1121 1122 ! Transfer ice energy to category jl2 by ridging1123 DO jk = 1, nlay_i1124 DO ij = 1, icells1125 ji = indxi(ij)1126 jj = indxj(ij)1127 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj) * erdg2(ji,jj,jk)1128 END DO1129 END DO1130 !1131 END DO ! jl2 (new ridges)1132 1133 DO jl2 = 1, jpl1134 1135 DO ij = 1, icells1136 ji = indxi(ij)1137 jj = indxj(ij)1138 ! Compute the fraction of rafted ice area and volume going to1139 ! thickness category jl2, transfer area, volume, and energy accordingly.1140 !1141 IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN1142 a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + arft2 (ji,jj)1143 v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + virft (ji,jj)1144 v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrft (ji,jj) * rn_fsnowrft1145 e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrft (ji,jj) * rn_fsnowrft1146 smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + smrft (ji,jj)1147 oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirft2(ji,jj)1148 ENDIF1149 !1150 END DO1151 1152 ! Transfer rafted ice energy to category jl21153 DO jk = 1, nlay_i1154 DO ij = 1, icells1155 ji = indxi(ij)1156 jj = indxj(ij)1157 IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN1158 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk)1159 ENDIF1160 END DO1161 END DO1162 1163 END DO1164 1165 END DO ! jl1 (deforming categories)1166 1167 ! Conservation check1168 IF ( con_i ) THEN1169 CALL lim_column_sum (jpl, v_i, vice_final)1170 fieldid = ' v_i : limitd_me '1171 CALL lim_cons_check (vice_init, vice_final, 1.0e-6, fieldid)1172 1173 CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_final )1174 fieldid = ' e_i : limitd_me '1175 CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid)1176 1177 DO ji = mi0(iiceprt), mi1(iiceprt)1178 DO jj = mj0(jiceprt), mj1(jiceprt)1179 WRITE(numout,*) ' vice_init : ', vice_init (ji,jj)1180 WRITE(numout,*) ' vice_final : ', vice_final(ji,jj)1181 WRITE(numout,*) ' eice_init : ', eice_init (ji,jj)1182 WRITE(numout,*) ' eice_final : ', eice_final(ji,jj)1183 END DO1184 END DO1185 ENDIF1186 !1187 CALL wrk_dealloc( (jpi+1)*(jpj+1), indxi, indxj )1188 CALL wrk_dealloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final )1189 CALL wrk_dealloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 )1190 CALL wrk_dealloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw, oirdg1, oirdg2 )1191 CALL wrk_dealloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )1192 CALL wrk_dealloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )1193 CALL wrk_dealloc( jpi, jpj, nlay_i, eirft, erdg1, erdg2, ersw )1194 CALL wrk_dealloc( jpi, jpj, nlay_i, jpl, eicen_init )1195 !1196 END SUBROUTINE lim_itd_me_ridgeshift1197 933 1198 934 SUBROUTINE lim_itd_me_init
Note: See TracChangeset
for help on using the changeset viewer.