Changeset 15742
- Timestamp:
- 2022-03-08T16:20:40+01:00 (2 years ago)
- Location:
- NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix
- Files:
-
- 5 edited
-
cfgs/SHARED/field_def_nemo-ice.xml (modified) (2 diffs)
-
src/ICE/icedyn_rdgrft.F90 (modified) (27 diffs)
-
src/ICE/icedyn_rhg_evp.F90 (modified) (8 diffs)
-
src/ICE/icewri.F90 (modified) (1 diff)
-
src/OCE/timing.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/cfgs/SHARED/field_def_nemo-ice.xml
r15570 r15742 91 91 <field id="yield22" long_name="yield surface tensor component 22" standard_name="yield22" unit="N/m" /> 92 92 <field id="yield12" long_name="yield surface tensor component 12" standard_name="yield12" unit="N/m" /> 93 <field id="beta_evp" long_name="Relaxation parameter of ice rheology (beta)" standard_name="relaxation_parameter_of_ice_rheology" unit="" /> 94 93 <field id="beta_evp" long_name="Relaxation parameter of ice rheology (beta)" standard_name="relaxation_parameter_of_ice_rheology" unit="" /> 94 <field id="opning" long_name="Lead area opening rate" standard_name="opning" unit="1/s" /> 95 <field id="dairdg1dt" long_name="Ridging ice area loss rate" standard_name="dairdg1dt" unit="1/s" /> 96 <field id="dairft1dt" long_name="Rafting ice area loss rate" 97 standard_name="dairft1dt" unit="1/s" /> 98 <field id="dairdg2dt" long_name="New ridged ice area gain rate" standard_name="dairdg2dt" unit="1/s" /> 99 <field id="dairft2dt" long_name="New rafted ice area gain rate" standard_name="dairft2dt" unit="1/s" /> 100 95 101 <!-- surface heat fluxes --> 96 102 <field id="qt_ice" long_name="total heat flux at ice surface" standard_name="surface_downward_heat_flux_in_air" unit="W/m2" /> … … 409 415 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 410 416 <field field_ref="icedlt" name="sidelta" /> 417 <field field_ref="opning" name="opning" /> 418 <field field_ref="dairdg1dt" name="dairdg1dt" /> 419 <field field_ref="dairft1dt" name="dairft1dt" /> 420 <field field_ref="dairdg2dt" name="dairdg2dt" /> 421 <field field_ref="dairft2dt" name="dairft2dt" /> 422 411 423 412 424 <!-- heat fluxes --> -
NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rdgrft.F90
r15570 r15742 43 43 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: opning ! rate of opening due to divergence/shear 44 44 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: closing_gross ! rate at which area removed, not counting area of new ridges 45 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: zaksum ! normalisation factor 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: zasum ! sum of ice area plus open water 45 47 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: apartf ! participation function; fraction of ridging/closing associated w/ category n 48 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zhi ! ice thickness in category n 46 49 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrmin ! minimum ridge thickness 47 50 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrmax ! maximum ridge thickness 51 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrexp ! e-folding ridge thickness 48 52 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hraft ! thickness of rafted ice 49 53 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hi_hrdg ! thickness of ridging ice / mean ridge thickness … … 53 57 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ze_s_2d 54 58 ! 59 ! For ridging diagnostics 60 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: airdg1 ! Ridging ice area loss 61 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: airft1 ! Rafting ice area loss 62 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: airdg2 ! New ridged ice area gain 63 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: airft2 ! New rafted ice area gain 64 ! 55 65 REAL(wp), PARAMETER :: hrdg_hi_min = 1.1_wp ! min ridge thickness multiplier: min(hrdg/hi) 56 REAL(wp), PARAMETER :: hi_hrft = 0.5_wp ! rafting multipli yer: (hi/hraft)66 REAL(wp), PARAMETER :: hi_hrft = 0.5_wp ! rafting multiplier: (hi/hraft) 57 67 ! 58 68 ! ** namelist (namdyn_rdgrft) ** 59 69 LOGICAL :: ln_str_H79 ! ice strength parameterization (Hibler79) 60 REAL(wp) :: rn_pstar ! determines ice strength, Hibler JPO79 70 REAL(wp) :: rn_pstar ! determines ice strength, Hibler JPO79 71 LOGICAL :: ln_str_R75 ! ice strength parameterization (Rothrock75) 72 REAL(wp) :: rn_Cf ! ratio of ridging work to PE change in ridging, R75 73 INTEGER :: ismooth ! smoothing the resistance to deformation (strength) 74 LOGICAL :: ln_redist_ufm ! uniform redistribution of ridged ice (Hibler, 1980) 75 LOGICAL :: ln_redist_exp ! exponential redistribution of ridged ice 76 REAL(wp) :: max_strength ! maximum strength cap 77 REAL(wp) :: rn_murdg ! gives e-folding scale of ridged ice (m^.5) for ln_redist_exp 61 78 REAL(wp) :: rn_csrdg ! fraction of shearing energy contributing to ridging 62 79 LOGICAL :: ln_partf_lin ! participation function linear (Thorndike et al. (1975)) … … 86 103 !! *** ROUTINE ice_dyn_rdgrft_alloc *** 87 104 !!------------------------------------------------------------------- 88 ALLOCATE( closing_net(jpij) , opning(jpij) , closing_gross(jpij) , & 89 & apartf(jpij,0:jpl) , hrmin (jpij,jpl) , hraft(jpij,jpl) , aridge(jpij,jpl), & 90 & hrmax (jpij,jpl) , hi_hrdg(jpij,jpl) , araft(jpij,jpl) , & 91 & ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc ) 105 ALLOCATE( closing_net(jpij) , opning(jpij) , closing_gross(jpij) , & 106 & apartf(jpij,0:jpl) , hrmin(jpij,jpl) , hrmax(jpij,jpl) , zhi(jpij,jpl), & 107 & zaksum(jpij) , hraft(jpij,jpl) , hrexp(jpij,jpl) , & 108 & hi_hrdg(jpij,jpl) , araft(jpij,jpl) , aridge(jpij,jpl) , zasum(jpij), & 109 & ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), & 110 & airdg1(jpij), airft1(jpij), airdg2(jpij), airft2(jpij), STAT=ice_dyn_rdgrft_alloc ) 92 111 93 112 CALL mpp_sum ( 'icedyn_rdgrft', ice_dyn_rdgrft_alloc ) … … 139 158 REAL(wp), DIMENSION(jpij) :: zdivu, zdelt ! 1D divu_i & delta_i 140 159 REAL(wp), DIMENSION(jpij) :: zconv ! 1D rdg_conv (if EAP rheology) 160 161 REAL(wp), DIMENSION(jpi,jpj) :: opning_2d ! Lead opening diagnostic 162 REAL(wp), DIMENSION(jpi,jpj) :: airdg1_2d ! Ridging ice area loss diagnostic 163 REAL(wp), DIMENSION(jpi,jpj) :: airft1_2d ! Rafting ice area loss diagnostic 164 REAL(wp), DIMENSION(jpi,jpj) :: airdg2_2d ! New ridged ice area gain diagnostic 165 REAL(wp), DIMENSION(jpi,jpj) :: airft2_2d ! New rafted ice area gain diagnostic 166 REAL(wp), DIMENSION(jpi,jpj) :: dairdg1dt ! Ridging ice area loss rate diagnostic 167 REAL(wp), DIMENSION(jpi,jpj) :: dairft1dt ! Rafting ice area loss rate diagnostic 168 REAL(wp), DIMENSION(jpi,jpj) :: dairdg2dt ! New ridged ice area gain rate diagnostic 169 REAL(wp), DIMENSION(jpi,jpj) :: dairft2dt ! New rafted ice area gain rate diagnostic 170 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00 ! Mask for ice presence 141 171 ! 142 172 INTEGER, PARAMETER :: jp_itermax = 20 … … 151 181 IF(lwp) WRITE(numout,*)'ice_dyn_rdgrft: ice ridging and rafting' 152 182 IF(lwp) WRITE(numout,*)'~~~~~~~~~~~~~~' 153 ENDIF 183 ENDIF 184 185 ! Initialise ridging diagnostics 186 opning_2d(:,:) = 0.0_wp 187 dairdg1dt(:,:) = 0.0_wp 188 dairft1dt(:,:) = 0.0_wp 189 dairdg2dt(:,:) = 0.0_wp 190 dairft2dt(:,:) = 0.0_wp 191 airdg1_2d(:,:) = 0.0_wp 192 airft1_2d(:,:) = 0.0_wp 193 airdg2_2d(:,:) = 0.0_wp 194 airft2_2d(:,:) = 0.0_wp 154 195 155 196 !-------------------------------- … … 162 203 DO jj = 1, jpj 163 204 DO ji = 1, jpi 205 zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice. epsi06 used here rather than epsi10 as below for consistency with ice masking elsewhere 164 206 IF ( at_i(ji,jj) > epsi10 ) THEN 165 207 npti = npti + 1 … … 188 230 IF( ln_rhg_EVP ) closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 189 231 IF( ln_rhg_EAP ) closing_net(ji) = zconv(ji) 190 ! 232 ! 191 233 IF( zdivu(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) ) ! make sure the closing rate is large enough 192 234 ! ! to give asum = 1.0 after ridging 193 235 ! Opening rate (non-negative) that will give asum = 1.0 after ridging. 194 236 opning(ji) = closing_net(ji) + zdivu(ji) 237 ! 195 238 END DO 196 239 ! … … 266 309 CALL ice_dyn_1d2d( 2 ) ! --- Move to 2D arrays --- ! 267 310 311 ! --- Ridging diagnostics --- ! 312 313 IF( iom_use('opning') ) THEN 314 CALL tab_1d_2d( npti, nptidx(1:npti), opning(1:npti), opning_2d(:,:) ) 315 END IF 316 317 IF( iom_use('dairdg1dt') .OR. iom_use('dairft1dt') .OR. iom_use('dairdg2dt') .OR. iom_use('dairft2dt') ) THEN 318 ! 319 CALL tab_1d_2d( npti, nptidx(1:npti), airdg1(1:npti), airdg1_2d(:,:) ) 320 CALL tab_1d_2d( npti, nptidx(1:npti), airft1(1:npti), airft1_2d(:,:) ) 321 CALL tab_1d_2d( npti, nptidx(1:npti), airdg2(1:npti), airdg2_2d(:,:) ) 322 CALL tab_1d_2d( npti, nptidx(1:npti), airft2(1:npti), airft2_2d(:,:) ) 323 ! 324 DO jj = 1, jpj 325 DO ji = 1, jpi 326 ! 327 dairdg1dt(ji,jj) = airdg1_2d(ji,jj) * r1_rdtice 328 dairft1dt(ji,jj) = airft1_2d(ji,jj) * r1_rdtice 329 dairdg2dt(ji,jj) = airdg2_2d(ji,jj) * r1_rdtice 330 dairft2dt(ji,jj) = airft2_2d(ji,jj) * r1_rdtice 331 ! 332 END DO 333 END DO 334 ! 335 ENDIF 336 ! 337 ENDIF ! npti>0 338 339 340 IF( iom_use('opning') ) THEN 341 CALL lbc_lnk( 'icedyn_rdgrft', opning_2d(:,:), 'T', 1. ) 342 CALL iom_put( 'opning', opning_2d(:,:) * zmsk00(:,:) ) ! Lead area opening rate 268 343 ENDIF 269 344 345 IF( iom_use('dairdg1dt') ) THEN 346 CALL lbc_lnk( 'icedyn_rdgrft', dairdg1dt(:,:), 'T', 1. ) 347 CALL iom_put( 'dairdg1dt', dairdg1dt(:,:) * zmsk00(:,:) ) ! Ridging ice area loss rate 348 ENDIF 349 350 IF( iom_use('dairft1dt') ) THEN 351 CALL lbc_lnk( 'icedyn_rdgrft', dairft1dt(:,:), 'T', 1. ) 352 CALL iom_put( 'dairft1dt', dairft1dt(:,:) * zmsk00(:,:) ) ! Rafting ice area loss rate 353 ENDIF 354 355 IF( iom_use('dairdg2dt') ) THEN 356 CALL lbc_lnk( 'icedyn_rdgrft', dairdg2dt(:,:), 'T', 1. ) 357 CALL iom_put( 'dairdg2dt', dairdg2dt(:,:) * zmsk00(:,:) ) ! New ridged ice area gain rate 358 ENDIF 359 360 IF( iom_use('dairft2dt') ) THEN 361 CALL lbc_lnk( 'icedyn_rdgrft', dairft2dt(:,:), 'T', 1. ) 362 CALL iom_put( 'dairft2dt', dairft2dt(:,:) * zmsk00(:,:) ) ! New rafted ice area gain rate 363 ENDIF 364 365 ! ------- ! 366 270 367 CALL ice_var_agg( 1 ) 271 368 … … 280 377 281 378 282 SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i, pclosing_net ) 379 380 SUBROUTINE rdgrft_prep( pa_i_in, pv_i, pato_i, pclosing_net ) 283 381 !!------------------------------------------------------------------- 284 382 !! *** ROUTINE rdgrft_prep *** … … 290 388 !!------------------------------------------------------------------- 291 389 REAL(wp), DIMENSION(:) , INTENT(in) :: pato_i, pclosing_net 292 REAL(wp), DIMENSION(:,:), INTENT(in) :: pa_i, pv_i 390 REAL(wp), DIMENSION(:,:), INTENT(in) :: pa_i_in, pv_i 391 REAL(wp), DIMENSION(jpij,jpl) :: pa_i ! area adjusted for minimum zhi 293 392 !! 294 393 INTEGER :: ji, jl ! dummy loop indices 295 394 REAL(wp) :: z1_gstar, z1_astar, zhmean, zfac ! local scalar 296 REAL(wp), DIMENSION(jpij) :: zasum, z1_asum, zaksum ! sum of a_i+ato_i and reverse 297 REAL(wp), DIMENSION(jpij,jpl) :: zhi ! ice thickness 395 REAL(wp), DIMENSION(jpij) :: z1_asum ! inverse of sum of a_i+ato_i 298 396 REAL(wp), DIMENSION(jpij,-1:jpl) :: zGsum ! zGsum(n) = sum of areas in categories 0 to n 299 397 !-------------------------------------------------------------------- … … 301 399 z1_gstar = 1._wp / rn_gstar 302 400 z1_astar = 1._wp / rn_astar 401 402 pa_i(:,:) = pa_i_in(:,:) 303 403 304 404 ! ! Ice thickness needed for rafting … … 306 406 ELSEWHERE ; zhi(1:npti,:) = 0._wp 307 407 END WHERE 408 ! Make sure ice thickness is not below the minimum (from iceitd.F90; no minimum was causing issues for R75 409 ! strength) 410 WHERE( pa_i(1:npti,:) > epsi10 .AND. zhi(1:npti,:) < rn_himin ) 411 pa_i(1:npti,:) = pa_i(1:npti,:) * zhi(1:npti,:) / rn_himin 412 zhi(1:npti,:) = rn_himin 413 END WHERE 414 415 ! Set tiny values of pa_i to zero 416 ! If this isn't done, can end up with zhi=0 but apartf>0 417 WHERE( pa_i(1:npti,:) <= epsi10 ) ; pa_i(1:npti,:) = 0._wp 418 END WHERE 419 308 420 309 421 ! 1) Participation function (apartf): a(h) = b(h).g(h) … … 320 432 ! Compute total area of ice plus open water. 321 433 ! This is in general not equal to one because of divergence during transport 322 zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti, :), dim=2 )323 ! 434 zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,1:jpl), dim=2 ) 435 324 436 WHERE( zasum(1:npti) > epsi10 ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti) 325 437 ELSEWHERE ; z1_asum(1:npti) = 0._wp 326 438 END WHERE 327 ! 439 328 440 ! Compute cumulative thickness distribution function 329 441 ! Compute the cumulative thickness distribution function zGsum, … … 333 445 zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti) 334 446 DO jl = 1, jpl 335 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)447 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 correct (and not jpl) 336 448 END DO 337 449 ! … … 398 510 ENDIF 399 511 512 513 400 514 ! 2) Transfer function 401 515 !----------------------------------------------------------------- 516 ! If assuming ridged ice is uniformly distributed between hrmin and 517 ! hrmax (ln_redist_ufm): 518 ! 402 519 ! Compute max and min ridged ice thickness for each ridging category. 403 ! Assume ridged ice is uniformly distributed between hrmin and hrmax.404 520 ! 405 521 ! This parameterization is a modified version of Hibler (1980). … … 416 532 ! (relative to the Hibler formulation) when very thick ice is ridging. 417 533 ! 418 ! zaksum = net area removed/ total area removed 419 ! where total area removed = area of ice that ridges 420 ! net area removed = total area removed - area of new ridges 534 !----------------------------------------------------------------- 535 ! If assuming ridged ice ITD is a negative exponential 536 ! (ln_redist_exp) and following CICE implementation: 537 ! 538 ! g(h) ~ exp[-(h-hrmin)/hrexp], h >= hrmin 539 ! 540 ! where hrmin is the minimum thickness of ridging ice and 541 ! hrexp is the e-folding thickness. 542 ! 543 ! Here, assume as above that hrmin = min(2*hi, hi+maxraft). 544 ! That is, the minimum ridge thickness results from rafting, 545 ! unless the ice is thicker than maxraft. 546 ! 547 ! Also, assume that hrexp = mu_rdg*sqrt(hi). 548 ! The parameter mu_rdg is tuned to give e-folding scales mostly 549 ! in the range 2-4 m as observed by upward-looking sonar. 550 ! 551 ! Values of mu_rdg in the right column give ice strengths 552 ! roughly equal to values of Hstar in the left column 553 ! (within ~10 kN/m for typical ITDs): 554 ! 555 ! Hstar mu_rdg 556 ! 557 ! 25 3.0 558 ! 50 4.0 559 ! 75 5.0 560 ! 100 6.0 561 ! 562 ! zaksum = net area removed/ total area participating 563 ! where total area participating = area of ice that ridges 564 ! net area removed = total area participating - area of new ridges 421 565 !----------------------------------------------------------------- 422 566 zfac = 1._wp / hi_hrft … … 428 572 zhmean = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min ) 429 573 hrmin (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) ) 430 hrmax (ji,jl) = 2._wp * zhmean - hrmin(ji,jl)431 574 hraft (ji,jl) = zhi(ji,jl) * zfac 432 hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 ) 575 ! 576 IF ( ln_redist_ufm ) THEN 577 hrmax (ji,jl) = 2._wp * zhmean - hrmin(ji,jl) 578 hi_hrdg(ji,jl) = zhi(ji,jl) / zhmean 579 ELSE IF ( ln_redist_exp ) THEN 580 hrexp (ji,jl) = rn_murdg * SQRT( zhi(ji,jl) ) 581 hi_hrdg(ji,jl) = zhi(ji,jl) / ( hrmin(ji,jl) + hrexp(ji,jl) ) 582 END IF 433 583 ! 434 584 ! Normalization factor : zaksum, ensures mass conservation … … 439 589 hrmax (ji,jl) = 0._wp 440 590 hraft (ji,jl) = 0._wp 591 hrexp (ji,jl) = 0._wp 441 592 hi_hrdg(ji,jl) = 1._wp 593 zaksum (ji) = 0._wp 442 594 ENDIF 443 595 END DO … … 493 645 INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices 494 646 REAL(wp) :: hL, hR, farea ! left and right limits of integration and new area going to jl2 647 REAL(wp) :: expL, expR ! exponentials involving hL, hR 495 648 REAL(wp) :: vsw ! vol of water trapped into ridges 496 649 REAL(wp) :: afrdg, afrft ! fraction of category area ridged/rafted 497 REAL(wp) :: airdg1,oirdg1, aprdg1, virdg1, sirdg1498 REAL(wp) :: airft1,oirft1, aprft1499 REAL(wp), DIMENSION(jpij) :: airdg2,oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, vlrdg ! area etc of new ridges500 REAL(wp), DIMENSION(jpij) :: airft2,oirft2, aprft2, virft , sirft , vsrft, vprft, vlrft ! area etc of rafted ice650 REAL(wp) :: oirdg1, aprdg1, virdg1, sirdg1 651 REAL(wp) :: oirft1, aprft1 652 REAL(wp), DIMENSION(jpij) :: oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, vlrdg ! area etc of new ridges 653 REAL(wp), DIMENSION(jpij) :: oirft2, aprft2, virft , sirft , vsrft, vprft, vlrft ! area etc of rafted ice 501 654 ! 502 655 REAL(wp), DIMENSION(jpij) :: ersw ! enth of water trapped into ridges … … 538 691 539 692 ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2) 540 airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice541 airft1 = araft (ji,jl1) * closing_gross(ji) * rdt_ice542 543 airdg2(ji) = airdg1 * hi_hrdg(ji,jl1)544 airft2(ji) = airft1 * hi_hrft693 airdg1(ji) = aridge(ji,jl1) * closing_gross(ji) * rdt_ice 694 airft1(ji) = araft (ji,jl1) * closing_gross(ji) * rdt_ice 695 696 airdg2(ji) = airdg1(ji) * hi_hrdg(ji,jl1) 697 airft2(ji) = airft1(ji) * hi_hrft 545 698 546 699 ! ridging /rafting fractions 547 afrdg = airdg1 * z1_ai(ji)548 afrft = airft1 * z1_ai(ji)700 afrdg = airdg1(ji) * z1_ai(ji) 701 afrft = airft1(ji) * z1_ai(ji) 549 702 550 703 ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges … … 602 755 ! Remove area, volume of new ridge to each category jl1 603 756 !------------------------------------------------------ 604 a_i_2d (ji,jl1) = a_i_2d (ji,jl1) - airdg1 - airft1757 a_i_2d (ji,jl1) = a_i_2d (ji,jl1) - airdg1(ji)- airft1(ji) 605 758 v_i_2d (ji,jl1) = v_i_2d (ji,jl1) - virdg1 - virft(ji) 606 759 v_s_2d (ji,jl1) = v_s_2d (ji,jl1) - vsrdg(ji) - vsrft(ji) … … 661 814 itest_rdg(1:npti) = 0 662 815 itest_rft(1:npti) = 0 816 ! 663 817 DO jl2 = 1, jpl 664 818 ! 665 819 DO ji = 1, npti 666 820 667 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 668 669 ! Compute the fraction of ridged ice area and volume going to thickness category jl2 821 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 822 ! 823 ! Compute the fraction of ridged ice area and volume going to thickness category jl2 824 ! 825 IF ( ln_redist_ufm ) THEN ! Hibler (1980) uniform formulation 826 ! 670 827 IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN 671 828 hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) ) … … 679 836 fvol(ji) = 0._wp 680 837 ENDIF 838 ! 839 ELSE IF ( ln_redist_exp ) THEN ! Lipscomb et al. (2007) exponential formulation 840 ! 841 IF (jl2 < jpl) THEN 842 ! 843 IF (hrmin(ji,jl1) <= hi_max(jl2)) THEN 844 hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) ) 845 hR = hi_max(jl2) 846 expL = EXP( -( hL - hrmin(ji,jl1) ) / hrexp(ji,jl1) ) 847 expR = EXP( -( hR - hrmin(ji,jl1) ) / hrexp(ji,jl1) ) 848 farea = expL - expR 849 fvol(ji) = ( ( hL + hrexp(ji,jl1) ) * expL & 850 - ( hR + hrexp(ji,jl1) ) * expR ) / ( hrmin(ji,jl1) + hrexp(ji,jl1) ) 851 ELSE 852 farea = 0._wp 853 fvol(ji) = 0._wp 854 END IF 855 ! 856 ELSE ! jl2 = jpl 857 ! 858 hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) ) 859 expL = EXP(-( hL - hrmin(ji,jl1) ) / hrexp(ji,jl1) ) 860 farea = expL 861 fvol(ji) = ( hL + hrexp(ji,jl1) ) * expL / ( hrmin(ji,jl1) + hrexp(ji,jl1) ) 862 ! 863 END IF ! jl2 < jpl 864 ! 865 END IF ! ridge redistribution 681 866 682 867 ! Compute the fraction of rafted ice area and volume going to thickness category jl2 … … 757 942 !! dissipated per unit area removed from the ice pack under compression, 758 943 !! and assumed proportional to the change in potential energy caused 759 !! by ridging. Note that only Hibler's formulation is stable and that760 !! ice strength has to be smoothed944 !! by ridging. Note that ice strength using Hibler's formulation must be 945 !! smoothed. 761 946 !!---------------------------------------------------------------------- 762 947 INTEGER :: ji, jj, jl ! dummy loop indices 763 INTEGER :: ismooth ! smoothing the resistance to deformation764 948 INTEGER :: itframe ! number of time steps for the P smoothing 765 949 REAL(wp) :: zp, z1_3 ! local scalars 950 REAL(wp) :: Cp ! proportional constant for PE 951 REAL(wp) :: h2rdg ! mean value of h^2 for new ridge 952 REAL(wp) :: dh2rdg ! change in mean value of h^2 per unit area consumed by ridging 766 953 REAL(wp), DIMENSION(jpi,jpj) :: zworka ! temporary array used here 767 954 REAL(wp), DIMENSION(jpi,jpj) :: zstrp1, zstrp2 ! strength at previous time steps 955 REAL(wp), DIMENSION(jpij) :: strength_1d ! 1-dimensional strength 956 REAL(wp), DIMENSION(jpij,jpl) :: strength_pe ! PE component of R75 strength 957 REAL(wp), DIMENSION(jpij) :: strength_pe_sum ! PE component of R75 strength 768 958 !!---------------------------------------------------------------------- 769 959 ! !--------------------------------------------------! 770 960 IF( ln_str_H79 ) THEN ! Ice strength => Hibler (1979) method ! 771 961 ! !--------------------------------------------------! 772 strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 773 ismooth = 1 962 ! Recommend using ismooth = 1 (spatial smoothing) 963 strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 964 ! !--------------------------------------------------! 965 ELSE IF( ln_str_R75 ) THEN ! Ice strength => Rothrock (1975) method ! 966 ! !--------------------------------------------------! 967 ! Recommend using ismooth = 0 (no smoothing needed) 968 ! Initialise the strength field at each timestep 969 strength(:,:) = 0._wp 970 ! 971 !-------------------------------- 972 ! 0) Identify grid cells with ice 973 !-------------------------------- 974 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 975 vt_i(:,:) = SUM( v_i(:,:,:), dim=3 ) 976 ato_i(:,:) = 1 - at_i(:,:) 977 ! 978 npti = 0 ; nptidx(:) = 0 979 DO jj = 1, jpj 980 DO ji = 1, jpi 981 IF ( at_i(ji,jj) > epsi10 ) THEN 982 npti = npti + 1 983 nptidx( npti ) = (jj - 1) * jpi + ji 984 ENDIF 985 END DO 986 END DO 987 988 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i ) 989 CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i ) 990 CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i ) 991 992 !----------------------------------- 993 ! Get thickness distribution of ice 994 !----------------------------------- 995 996 Cp = 0.5_wp * grav * (rhow-rhoi) * rhoi / rhow ! EF: This should go in an initialisation routine but unclear where, won't let me set as a parameter above 997 998 CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net ) 999 ! EF: Note that closing rate is not needed here, could split this function up, and refactor as initial rdgrft_prep is 1000 ! called here, before going into ridging iterations 1001 1002 ! Initialise 1003 strength_pe_sum(:) = 0._wp 1004 1005 IF ( npti > 0 ) THEN 1006 1007 DO jl = 1, jpl ! categories 1008 DO ji = 1, npti ! grid cells 1009 IF ( apartf(ji,jl) > 0._wp ) THEN 1010 ! 1011 IF ( ln_redist_ufm ) THEN ! Uniform redistribution of ridged ice 1012 h2rdg = (1._wp/3._wp) * (hrmax(ji,jl)**3._wp - hrmin(ji,jl)**3._wp) & 1013 / (hrmax(ji,jl) - hrmin(ji,jl)) 1014 ! 1015 ELSE IF ( ln_redist_exp ) THEN ! Exponential redistribution of ridged ice 1016 h2rdg = hrmin(ji,jl) * hrmin(ji,jl) & 1017 + 2._wp * hrmin(ji,jl) * hrexp(ji,jl) & 1018 + 2._wp * hrexp(ji,jl) * hrexp(ji,jl) 1019 END IF 1020 1021 dh2rdg = -zhi(ji,jl)*zhi(ji,jl) + h2rdg * hi_hrdg(ji,jl) 1022 1023 strength_pe_sum(ji) = strength_pe_sum(ji) + apartf(ji,jl) * dh2rdg 1024 1025 ELSE 1026 1027 strength_pe_sum(ji) = strength_pe_sum(ji) + 0._wp 1028 1029 END IF 1030 1031 END DO ! grid cells 1032 END DO ! categories 1033 1034 1035 DO ji = 1, npti 1036 IF ( zaksum(ji) > epsi10 ) THEN 1037 strength_1d(ji) = rn_Cf * Cp * strength_pe_sum(ji) / zaksum(ji) 1038 ELSE 1039 strength_1d(ji) = 0._wp 1040 END IF 1041 END DO 1042 1043 1044 ! Enforce a maximum for strength 1045 ! Richter-Menge and Elder (1998) estimate maximum in Beaufort Sea in wintertime of the order 150 kN/m 1046 WHERE( strength_1d(:) > max_strength ) ; strength_1d(:) = max_strength 1047 END WHERE 1048 1049 1050 ! Convert strength_1d to 2d 1051 CALL tab_1d_2d( npti, nptidx(1:npti), strength_1d(1:npti), strength(:,:) ) 1052 1053 END IF ! number of ice points 1054 1055 774 1056 ! !--------------------------------------------------! 775 ELSE ! Zero strength 1057 ELSE ! Zero strength (this just crashes) ! 776 1058 ! !--------------------------------------------------! 777 1059 strength(:,:) = 0._wp 778 ismooth = 0 779 ENDIF 780 ! !--------------------------------------------------! 1060 ! 1061 ENDIF ! Ice strength formulation 1062 ! 1063 !--------------------------------------------------! 781 1064 SELECT CASE( ismooth ) ! Smoothing ice strength ! 782 1065 ! !--------------------------------------------------! 1066 CASE( 0 ) !--- No smoothing 1067 1068 CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. ) 1069 783 1070 CASE( 1 ) !--- Spatial smoothing 784 1071 DO jj = 2, jpjm1 … … 918 1205 INTEGER :: ios ! Local integer output status for namelist read 919 1206 !! 920 NAMELIST/namdyn_rdgrft/ ln_str_H79, rn_pstar, rn_crhg, & 921 & rn_csrdg , & 1207 NAMELIST/namdyn_rdgrft/ ismooth, ln_str_H79, rn_pstar, & 1208 & ln_str_R75, rn_Cf, & 1209 & ln_redist_ufm, ln_redist_exp, & 1210 & max_strength, rn_murdg, & 1211 & rn_crhg, rn_csrdg, & 922 1212 & ln_partf_lin, rn_gstar, & 923 1213 & ln_partf_exp, rn_astar, & … … 939 1229 WRITE(numout,*) '~~~~~~~~~~~~~~~~~~' 940 1230 WRITE(numout,*) ' Namelist namdyn_rdgrft:' 941 WRITE(numout,*) ' ice strength parameterization Hibler (1979) ln_str_H79 = ', ln_str_H79 1231 WRITE(numout,*) ' ice strength parameterization Hibler (1979) ln_str_H79 = ', ln_str_H79 942 1232 WRITE(numout,*) ' 1st bulk-rheology parameter rn_pstar = ', rn_pstar 943 1233 WRITE(numout,*) ' 2nd bulk-rhelogy parameter rn_crhg = ', rn_crhg 1234 WRITE(numout,*) ' ice strength parameterization Rothrock (1975) ln_str_R75 = ', ln_str_R75 1235 WRITE(numout,*) ' ratio of ridging work to PE change in ridging rn_Cf = ', rn_Cf 1236 WRITE(numout,*) ' smoothing the resistance to deformation (strength) ismooth = ', ismooth 1237 WRITE(numout,*) ' maximum strength cap max_strength = ', max_strength 1238 WRITE(numout,*) ' uniform redistribution of ridged ice (Hibler, 1980) ln_redist_ufm = ', ln_redist_ufm 1239 WRITE(numout,*) ' exponential redistribution of ridged ice ln_redist_exp = ', ln_redist_exp 1240 WRITE(numout,*) ' e-folding scale of ridged ice for ln_redist_exp rn_murdg = ', rn_murdg 944 1241 WRITE(numout,*) ' Fraction of shear energy contributing to ridging rn_csrdg = ', rn_csrdg 945 1242 WRITE(numout,*) ' linear ridging participation function ln_partf_lin = ', ln_partf_lin … … 959 1256 ENDIF 960 1257 ! 1258 IF ( ( ln_str_H79 .AND. ln_str_R75 ) .OR. ( .NOT.ln_str_H79 .AND. .NOT.ln_str_R75 ) ) THEN 1259 CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one ice strength formulation (ln_str_H79 or ln_str_R75)' ) 1260 ENDIF 1261 ! 961 1262 IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN 962 1263 CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' ) 963 1264 ENDIF 1265 ! 1266 IF ( ( ln_redist_ufm .AND. ln_redist_exp ) .OR. ( .NOT.ln_redist_ufm .AND. .NOT.ln_redist_exp ) ) THEN 1267 CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one ice redistribution function (ln_redist_ufm or ln_redist_exp)' ) 1268 ENDIF 1269 ! 1270 IF ( ( ismooth .NE. 0 ) .AND. ( ismooth .NE. 1 ) .AND. (ismooth .NE. 2 ) ) THEN 1271 CALL ctl_stop( 'ice_dyn_rdgrft_init: ismooth must be one of 0 (no smoothing), 1 (spatial), or 2 (temporal)' ) 1272 ENDIF 1273 ! 964 1274 ! 965 1275 IF( .NOT. ln_icethd ) THEN -
NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rhg_evp.F90
r15692 r15742 137 137 ! 138 138 REAL(wp), DIMENSION(jpi,jpj) :: zdelta, zp_delt ! delta and P/delta at T points 139 REAL(wp), DIMENSION(jpi,jpj) :: zten_i ! tension139 REAL(wp), DIMENSION(jpi,jpj) :: zten_i, zshear ! tension, shear 140 140 REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017 141 141 ! … … 773 773 & ) * 0.25_wp * r1_e1e2t(ji,jj) 774 774 775 ! shear at T points775 ! maximum shear rate at T points (includes tension, output only) 776 776 pshear_i(ji,jj) = SQRT( zdt2 + zds2 ) 777 778 ! shear at T-points 779 zshear(ji,jj) = SQRT( zds2 ) 777 780 778 781 ! divergence at T points … … 789 792 END DO 790 793 CALL lbc_lnk_multi( 'icedyn_rhg_evp', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1., zten_i, 'T', 1., & 791 & zs1 , 'T', 1., zs2 , 'T', 1., zs12 , 'F', 1. )794 & zs1 , 'T', 1., zs2 , 'T', 1., zs12 , 'F', 1., zshear, 'T', 1. ) 792 795 793 796 ! --- Store the stress tensor for the next time step --- ! … … 835 838 zsig1 = zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) ) 836 839 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 837 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj)838 840 zsig12 = zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp 841 839 842 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008) 840 843 zsig_I (ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 841 zsig_II(ji,jj) = SQRT( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 ) ! 2nd '' '' , aka maximum shear stress844 zsig_II(ji,jj) = SQRT( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 ) ! 2nd '' '' , aka maximum shear stress 842 845 843 846 END DO … … 866 869 ! and **deformations** at current iterates 867 870 ! following Lemieux & Dupont (2020) 868 zfac = zp_delt(ji,jj)869 zsig1 = zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl) )871 zfac = strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 872 zsig1 = zfac * ( pdivu_i(ji,jj) - zdelta(ji,jj) ) 870 873 zsig2 = zfac * z1_ecc2 * zten_i(ji,jj) 871 zsig12 = zfac * z1_ecc2 * pshear_i(ji,jj)872 874 zsig12 = zfac * z1_ecc2 * zshear(ji,jj) * 0.5_wp 875 873 876 ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point 874 877 zsig_I(ji,jj) = zsig1 * 0.5_wp ! 1st stress invariant, aka average normal stress, aka negative pressure 875 zsig_II(ji,jj) = SQRT( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 ) ! 2nd '' '' , aka maximum shear stress878 zsig_II(ji,jj) = SQRT( zsig2 * zsig2 * 0.25_wp + zsig12 * zsig12 ) ! 2nd '' '' , aka maximum shear stress 876 879 877 880 ! Normalized principal stresses (used to display the ellipse) … … 882 885 END DO 883 886 ! 884 IF( iom_use('sig1_pnorm') ) CALL iom_put( 'sig1_pnorm' , zsig1_p(:,:) * zmsk00(:,:) ) 887 IF( iom_use('sig1_pnorm') ) CALL iom_put( 'sig1_pnorm' , zsig1_p(:,:) * zmsk00(:,:) ) 885 888 IF( iom_use('sig2_pnorm') ) CALL iom_put( 'sig2_pnorm' , zsig2_p(:,:) * zmsk00(:,:) ) 886 889 … … 1006 1009 1007 1010 ! time 1008 it = ( kt - 1) * kitermax + kiter1011 it = ( kt - nit000 ) * kitermax + kiter 1009 1012 1010 1013 ! convergence … … 1026 1029 istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 1027 1030 ! close file 1028 IF( kt == nitend - nn_fsbc + 1 ) istatus = NF90_CLOSE(ncvgid)1031 IF( kt == nitend - nn_fsbc + 1 .AND. kiter == kitermax ) istatus = NF90_CLOSE(ncvgid) 1029 1032 ENDIF 1030 1033 -
NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icewri.F90
r15570 r15742 254 254 CALL iom_rstput( 0, 0, kid, 'sivolu', vt_i ) ! Ice volume 255 255 CALL iom_rstput( 0, 0, kid, 'sidive', divu_i*1.0e8 ) ! Ice divergence 256 CALL iom_rstput( 0, 0, kid, 'sishea', shear_i*1.0e8 ) 256 CALL iom_rstput( 0, 0, kid, 'sishea', shear_i*1.0e8 ) ! Ice shear 257 257 CALL iom_rstput( 0, 0, kid, 'si_amp', at_ip ) ! Melt pond fraction 258 258 CALL iom_rstput( 0, 0, kid, 'si_vmp', vt_ip ) ! Melt pond volume -
NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/OCE/timing.F90
r14075 r15742 295 295 296 296 ll_averep = .TRUE. 297 zperc = 0._wp 297 298 298 299 ! total CPU and elapse
Note: See TracChangeset
for help on using the changeset viewer.
