Changeset 14789 for NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/ICE/icedyn_rdgrft.F90
- Timestamp:
- 2021-05-05T13:18:04+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icedyn_rdgrft.F90 (modified) (48 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev _r12970_AGRIF_CMEMSext/AGRIF5 ^/vendors/AGRIF/dev@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 ^/vendors/PPR@HEAD ext/PPR 8 9 9 10 # SETTE 10 ^/utils/CI/sette@1 3559sette11 ^/utils/CI/sette@14244 sette
-
- Property svn:externals
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/ICE/icedyn_rdgrft.F90
r13618 r14789 2 2 !!====================================================================== 3 3 !! *** MODULE icedyn_rdgrft *** 4 !! sea-ice : Mechanical impact on ice thickness distribution 4 !! sea-ice : Mechanical impact on ice thickness distribution 5 5 !!====================================================================== 6 !! History : ! 2006-02 (M. Vancoppenolle) Original code 6 !! History : ! 2006-02 (M. Vancoppenolle) Original code 7 7 !! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube] 8 8 !!---------------------------------------------------------------------- … … 16 16 !!---------------------------------------------------------------------- 17 17 USE dom_oce ! ocean domain 18 USE phycst ! physical constants (ocean directory) 18 USE phycst ! physical constants (ocean directory) 19 19 USE sbc_oce , ONLY : sss_m, sst_m ! surface boundary condition: ocean fields 20 20 USE ice1D ! sea-ice: thermodynamics … … 59 59 LOGICAL :: ln_str_H79 ! ice strength parameterization (Hibler79) 60 60 REAL(wp) :: rn_pstar ! determines ice strength, Hibler JPO79 61 REAL(wp) :: rn_csrdg ! fraction of shearing energy contributing to ridging 61 REAL(wp) :: rn_csrdg ! fraction of shearing energy contributing to ridging 62 62 LOGICAL :: ln_partf_lin ! participation function linear (Thorndike et al. (1975)) 63 63 REAL(wp) :: rn_gstar ! fractional area of young ice contributing to ridging 64 64 LOGICAL :: ln_partf_exp ! participation function exponential (Lipscomb et al. (2007)) 65 65 REAL(wp) :: rn_astar ! equivalent of G* for an exponential participation function 66 LOGICAL :: ln_ridging ! ridging of ice or not 66 LOGICAL :: ln_ridging ! ridging of ice or not 67 67 REAL(wp) :: rn_hstar ! thickness that determines the maximal thickness of ridged ice 68 68 REAL(wp) :: rn_porordg ! initial porosity of ridges (0.3 regular value) 69 69 REAL(wp) :: rn_fsnwrdg ! fractional snow loss to the ocean during ridging 70 70 REAL(wp) :: rn_fpndrdg ! fractional pond loss to the ocean during ridging 71 LOGICAL :: ln_rafting ! rafting of ice or not 72 REAL(wp) :: rn_hraft ! threshold thickness (m) for rafting / ridging 71 LOGICAL :: ln_rafting ! rafting of ice or not 72 REAL(wp) :: rn_hraft ! threshold thickness (m) for rafting / ridging 73 73 REAL(wp) :: rn_craft ! coefficient for smoothness of the hyperbolic tangent in rafting 74 74 REAL(wp) :: rn_fsnwrft ! fractional snow loss to the ocean during rafting … … 124 124 !! Hibler, W. D. III, 1980, MWR, 108, 1943-1973, 1980. 125 125 !! Rothrock, D. A., 1975: JGR, 80, 4514-4519. 126 !! Thorndike et al., 1975, JGR, 80, 4501-4513. 126 !! Thorndike et al., 1975, JGR, 80, 4501-4513. 127 127 !! Bitz et al., JGR, 2001 128 128 !! Amundrud and Melling, JGR 2005 129 !! Babko et al., JGR 2002 129 !! Babko et al., JGR 2002 130 130 !! 131 131 !! This routine is based on CICE code and authors William H. Lipscomb, … … 135 135 !! 136 136 INTEGER :: ji, jj, jk, jl ! dummy loop index 137 INTEGER :: iter, iterate_ridging ! local integer 137 INTEGER :: iter, iterate_ridging ! local integer 138 138 INTEGER :: ipti ! local integer 139 139 REAL(wp) :: zfac ! local scalar 140 140 INTEGER , DIMENSION(jpij) :: iptidx ! compute ridge/raft or not 141 141 REAL(wp), DIMENSION(jpij) :: zdivu, zdelt ! 1D divu_i & delta_i 142 ! 143 INTEGER, PARAMETER :: jp_itermax = 20 142 REAL(wp), DIMENSION(jpij) :: zconv ! 1D rdg_conv (if EAP rheology) 143 ! 144 INTEGER, PARAMETER :: jp_itermax = 20 144 145 !!------------------------------------------------------------------- 145 146 ! controls … … 152 153 IF(lwp) WRITE(numout,*)'ice_dyn_rdgrft: ice ridging and rafting' 153 154 IF(lwp) WRITE(numout,*)'~~~~~~~~~~~~~~' 154 ENDIF 155 ENDIF 155 156 156 157 !-------------------------------- … … 167 168 ENDIF 168 169 END_2D 169 170 170 171 !-------------------------------------------------------- 171 172 ! 1) Dynamical inputs (closing rate, divergence, opening) 172 173 !-------------------------------------------------------- 173 174 IF( npti > 0 ) THEN 174 175 175 176 ! just needed here 176 177 CALL tab_2d_1d( npti, nptidx(1:npti), zdelt (1:npti) , delta_i ) 178 CALL tab_2d_1d( npti, nptidx(1:npti), zconv (1:npti) , rdg_conv ) 177 179 ! needed here and in the iteration loop 178 180 CALL tab_2d_1d( npti, nptidx(1:npti), zdivu (1:npti) , divu_i) ! zdivu is used as a work array here (no change in divu_i) … … 182 184 183 185 DO ji = 1, npti 184 ! closing_net = rate at which open water area is removed + ice area removed by ridging 186 ! closing_net = rate at which open water area is removed + ice area removed by ridging 185 187 ! - ice area added in new ridges 186 closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 188 IF( ln_rhg_EVP .OR. ln_rhg_VP ) & 189 & closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 190 IF( ln_rhg_EAP ) closing_net(ji) = zconv(ji) 187 191 ! 188 192 IF( zdivu(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) ) ! make sure the closing rate is large enough … … 221 225 !----------------- 222 226 IF( npti > 0 ) THEN 223 227 224 228 CALL ice_dyn_1d2d( 1 ) ! --- Move to 1D arrays --- ! 225 229 226 230 iter = 1 227 iterate_ridging = 1 231 iterate_ridging = 1 228 232 ! !----------------------! 229 233 DO WHILE( iterate_ridging > 0 .AND. iter < jp_itermax ) ! ridging iterations ! … … 264 268 265 269 ENDIF 266 267 CALL ice_var_agg( 1 ) 270 271 CALL ice_var_agg( 1 ) 268 272 269 273 ! controls … … 283 287 !! ** Purpose : preparation for ridging calculations 284 288 !! 285 !! ** Method : Compute the thickness distribution of the ice and open water 289 !! ** Method : Compute the thickness distribution of the ice and open water 286 290 !! participating in ridging and of the resulting ridges. 287 291 !!------------------------------------------------------------------- 288 REAL(wp), DIMENSION(:) , INTENT(in) :: pato_i, pclosing_net 289 REAL(wp), DIMENSION(:,:), INTENT(in) :: pa_i, pv_i 292 REAL(wp), DIMENSION(:) , INTENT(in) :: pato_i, pclosing_net 293 REAL(wp), DIMENSION(:,:), INTENT(in) :: pa_i, pv_i 290 294 !! 291 295 INTEGER :: ji, jl ! dummy loop indices 292 296 REAL(wp) :: z1_gstar, z1_astar, zhmean, zfac ! local scalar 293 REAL(wp), DIMENSION(jpij) :: zasum, z1_asum, zaksum ! sum of a_i+ato_i and reverse 297 REAL(wp), DIMENSION(jpij) :: zasum, z1_asum, zaksum ! sum of a_i+ato_i and reverse 294 298 REAL(wp), DIMENSION(jpij,jpl) :: zhi ! ice thickness 295 299 REAL(wp), DIMENSION(jpij,-1:jpl) :: zGsum ! zGsum(n) = sum of areas in categories 0 to n … … 317 321 ! This is analogous to 318 322 ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975). 319 ! assuming b(h) = (2/Gstar) * (1 - G(h)/Gstar). 323 ! assuming b(h) = (2/Gstar) * (1 - G(h)/Gstar). 320 324 ! 321 325 ! apartf = integrating b(h)g(h) between the category boundaries … … 342 346 ! 343 347 IF( ln_partf_lin ) THEN !--- Linear formulation (Thorndike et al., 1975) 344 DO jl = 0, jpl 348 DO jl = 0, jpl 345 349 DO ji = 1, npti 346 350 IF ( zGsum(ji,jl) < rn_gstar ) THEN … … 357 361 ! 358 362 ELSEIF( ln_partf_exp ) THEN !--- Exponential, more stable formulation (Lipscomb et al, 2007) 359 ! 363 ! 360 364 zfac = 1._wp / ( 1._wp - EXP(-z1_astar) ) 361 365 DO jl = -1, jpl … … 387 391 END DO 388 392 END DO 389 ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN !- rafting alone 393 ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN !- rafting alone 390 394 DO jl = 1, jpl 391 395 DO ji = 1, npti … … 398 402 DO ji = 1, npti 399 403 aridge(ji,jl) = 0._wp 400 araft (ji,jl) = 0._wp 404 araft (ji,jl) = 0._wp 401 405 END DO 402 406 END DO … … 407 411 ! Compute max and min ridged ice thickness for each ridging category. 408 412 ! Assume ridged ice is uniformly distributed between hrmin and hrmax. 409 ! 413 ! 410 414 ! This parameterization is a modified version of Hibler (1980). 411 415 ! The mean ridging thickness, zhmean, is proportional to hi^(0.5) 412 416 ! and for very thick ridging ice must be >= hrdg_hi_min*hi 413 417 ! 414 ! The minimum ridging thickness, hrmin, is equal to 2*hi 418 ! The minimum ridging thickness, hrmin, is equal to 2*hi 415 419 ! (i.e., rafting) and for very thick ridging ice is 416 420 ! constrained by hrmin <= (zhmean + hi)/2. 417 ! 421 ! 418 422 ! The maximum ridging thickness, hrmax, is determined by zhmean and hrmin. 419 423 ! … … 441 445 & + araft (ji,jl) * ( 1._wp - hi_hrft ) 442 446 ELSE 443 hrmin (ji,jl) = 0._wp 444 hrmax (ji,jl) = 0._wp 445 hraft (ji,jl) = 0._wp 447 hrmin (ji,jl) = 0._wp 448 hrmax (ji,jl) = 0._wp 449 hraft (ji,jl) = 0._wp 446 450 hi_hrdg(ji,jl) = 1._wp 447 451 ENDIF … … 451 455 ! 3) closing_gross 452 456 !----------------- 453 ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate. 457 ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate. 454 458 ! NOTE: 0 < aksum <= 1 455 459 WHERE( zaksum(1:npti) > epsi10 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 456 460 ELSEWHERE ; closing_gross(1:npti) = 0._wp 457 461 END WHERE 458 462 459 463 ! correction to closing rate if excessive ice removal 460 464 !---------------------------------------------------- … … 468 472 ENDIF 469 473 END DO 470 END DO 474 END DO 471 475 472 476 ! 4) correction to opening if excessive open water removal … … 474 478 ! Reduce the closing rate if more than 100% of the open water would be removed 475 479 ! Reduce the opening rate in proportion 476 DO ji = 1, npti 480 DO ji = 1, npti 477 481 zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rDt_ice 478 482 IF( zfac < 0._wp ) THEN ! would lead to negative ato_i 479 opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_Dt_ice 483 opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_Dt_ice 480 484 ELSEIF( zfac > zasum(ji) ) THEN ! would lead to ato_i > asum 481 opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_Dt_ice 485 opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_Dt_ice 482 486 ENDIF 483 487 END DO … … 499 503 REAL(wp) :: hL, hR, farea ! left and right limits of integration and new area going to jl2 500 504 REAL(wp) :: vsw ! vol of water trapped into ridges 501 REAL(wp) :: afrdg, afrft ! fraction of category area ridged/rafted 505 REAL(wp) :: afrdg, afrft ! fraction of category area ridged/rafted 502 506 REAL(wp) :: airdg1, oirdg1, aprdg1, virdg1, sirdg1 503 507 REAL(wp) :: airft1, oirft1, aprft1 … … 512 516 REAL(wp), DIMENSION(jpij,nlay_s) :: esrft ! snow energy of rafting ice 513 517 REAL(wp), DIMENSION(jpij,nlay_i) :: eirft ! ice energy of rafting ice 514 REAL(wp), DIMENSION(jpij,nlay_s) :: esrdg ! enth*volume of new ridges 518 REAL(wp), DIMENSION(jpij,nlay_s) :: esrdg ! enth*volume of new ridges 515 519 REAL(wp), DIMENSION(jpij,nlay_i) :: eirdg ! enth*volume of new ridges 516 520 ! … … 525 529 ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rDt_ice ) 526 530 END DO 527 528 ! 2) compute categories in which ice is removed (jl1) 531 532 ! 2) compute categories in which ice is removed (jl1) 529 533 !---------------------------------------------------- 530 534 DO jl1 = 1, jpl 531 535 532 IF( nn_icesal /= 2 ) THEN 536 IF( nn_icesal /= 2 ) THEN 533 537 CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,jl1) ) 534 538 ENDIF … … 541 545 ELSE ; z1_ai(ji) = 0._wp 542 546 ENDIF 543 547 544 548 ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2) 545 549 airdg1 = aridge(ji,jl1) * closing_gross(ji) * rDt_ice … … 567 571 sirdg2(ji) = sv_i_2d(ji,jl1) * afrdg + vsw * sss_1d(ji) 568 572 oirdg1 = oa_i_2d(ji,jl1) * afrdg 569 oirdg2(ji) = oa_i_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1) 573 oirdg2(ji) = oa_i_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1) 570 574 571 575 virft(ji) = v_i_2d (ji,jl1) * afrft 572 576 vsrft(ji) = v_s_2d (ji,jl1) * afrft 573 sirft(ji) = sv_i_2d(ji,jl1) * afrft 574 oirft1 = oa_i_2d(ji,jl1) * afrft 575 oirft2(ji) = oa_i_2d(ji,jl1) * afrft * hi_hrft 576 577 IF ( ln_pnd_LEV ) THEN577 sirft(ji) = sv_i_2d(ji,jl1) * afrft 578 oirft1 = oa_i_2d(ji,jl1) * afrft 579 oirft2(ji) = oa_i_2d(ji,jl1) * afrft * hi_hrft 580 581 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 578 582 aprdg1 = a_ip_2d(ji,jl1) * afrdg 579 583 aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1) … … 591 595 wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoi * r1_Dt_ice ! increase in ice volume due to seawater frozen in voids 592 596 sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoi * r1_Dt_ice 593 hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_Dt_ice ! > 0 [W.m-2] 597 hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_Dt_ice ! > 0 [W.m-2] 594 598 595 599 ! Put the snow lost by ridging into the ocean … … 602 606 sirdg2(ji) = sirdg2(ji) - vsw * ( sss_1d(ji) - s_i_1d(ji) ) ! ridge salinity = s_i 603 607 sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_Dt_ice & ! put back sss_m into the ocean 604 & - s_i_1d(ji) * vsw * rhoi * r1_Dt_ice ! and get s_i from the ocean 608 & - s_i_1d(ji) * vsw * rhoi * r1_Dt_ice ! and get s_i from the ocean 605 609 ENDIF 606 610 … … 612 616 sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - sirdg1 - sirft(ji) 613 617 oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1 - oirft1 614 IF ( ln_pnd_LEV ) THEN618 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 615 619 a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1 - aprft1 616 620 v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji) … … 639 643 ! Remove energy of new ridge to each category jl1 640 644 !------------------------------------------------- 641 ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 645 ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 642 646 ENDIF 643 647 END DO 644 648 END DO 645 649 646 650 ! special loop for e_i because of layers jk 647 651 DO jk = 1, nlay_i … … 657 661 ! Remove energy of new ridge to each category jl1 658 662 !------------------------------------------------- 659 ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 663 ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 660 664 ENDIF 661 665 END DO 662 666 END DO 663 664 ! 3) compute categories in which ice is added (jl2) 667 668 ! 3) compute categories in which ice is added (jl2) 665 669 !-------------------------------------------------- 666 670 itest_rdg(1:npti) = 0 667 671 itest_rft(1:npti) = 0 668 DO jl2 = 1, jpl 672 DO jl2 = 1, jpl 669 673 ! 670 674 DO ji = 1, npti … … 681 685 itest_rdg(ji) = 1 ! test for conservation 682 686 ELSE 683 farea = 0._wp 684 fvol(ji) = 0._wp 687 farea = 0._wp 688 fvol(ji) = 0._wp 685 689 ENDIF 686 690 … … 697 701 ! Sometimes thickness is larger than hi_max(jpl) because of advection scheme (for very small areas) 698 702 ! Then ice volume is removed from one category but the ridging/rafting scheme 699 ! does not know where to move it, leading to a conservation issue. 703 ! does not know where to move it, leading to a conservation issue. 700 704 IF( itest_rdg(ji) == 0 .AND. jl2 == jpl ) THEN ; farea = 1._wp ; fvol(ji) = 1._wp ; ENDIF 701 705 IF( itest_rft(ji) == 0 .AND. jl2 == jpl ) zswitch(ji) = 1._wp … … 709 713 v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji) + & 710 714 & vsrft (ji) * rn_fsnwrft * zswitch(ji) ) 711 IF ( ln_pnd_LEV ) THEN715 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 712 716 v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + ( vprdg (ji) * rn_fpndrdg * fvol (ji) & 713 717 & + vprft (ji) * rn_fpndrft * zswitch(ji) ) 714 a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + ( aprdg2(ji) * rn_fpndrdg * farea & 718 a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + ( aprdg2(ji) * rn_fpndrdg * farea & 715 719 & + aprft2(ji) * rn_fpndrft * zswitch(ji) ) 716 720 IF ( ln_pnd_lids ) THEN … … 719 723 ENDIF 720 724 ENDIF 721 725 722 726 ENDIF 723 727 … … 737 741 DO ji = 1, npti 738 742 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) & 739 & ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji) 743 & ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji) 740 744 END DO 741 745 END DO … … 759 763 !! ** Purpose : computes ice strength used in dynamics routines of ice thickness 760 764 !! 761 !! ** Method : Compute the strength of the ice pack, defined as the energy (J m-2) 765 !! ** Method : Compute the strength of the ice pack, defined as the energy (J m-2) 762 766 !! dissipated per unit area removed from the ice pack under compression, 763 767 !! and assumed proportional to the change in potential energy caused … … 776 780 ! !--------------------------------------------------! 777 781 strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 778 ismooth = 1 782 ismooth = 1 ! original code 783 ! ismooth = 0 ! try for EAP stability 779 784 ! !--------------------------------------------------! 780 785 ELSE ! Zero strength ! … … 788 793 CASE( 1 ) !--- Spatial smoothing 789 794 DO_2D( 0, 0, 0, 0 ) 790 IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN 795 IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN 791 796 zworka(ji,jj) = ( 4.0 * strength(ji,jj) & 792 & + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 797 & + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 793 798 & + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) & 794 799 & ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) ) … … 797 802 ENDIF 798 803 END_2D 799 804 800 805 DO_2D( 0, 0, 0, 0 ) 801 806 strength(ji,jj) = zworka(ji,jj) … … 810 815 ! 811 816 DO_2D( 0, 0, 0, 0 ) 812 IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN 817 IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN 813 818 itframe = 1 ! number of time steps for the running mean 814 819 IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1 … … 826 831 END SUBROUTINE ice_strength 827 832 828 833 829 834 SUBROUTINE ice_dyn_1d2d( kn ) 830 835 !!----------------------------------------------------------------------- 831 !! *** ROUTINE ice_dyn_1d2d *** 832 !! 836 !! *** ROUTINE ice_dyn_1d2d *** 837 !! 833 838 !! ** Purpose : move arrays from 1d to 2d and the reverse 834 839 !!----------------------------------------------------------------------- … … 900 905 ! 901 906 END SUBROUTINE ice_dyn_1d2d 902 907 903 908 904 909 SUBROUTINE ice_dyn_rdgrft_init … … 906 911 !! *** ROUTINE ice_dyn_rdgrft_init *** 907 912 !! 908 !! ** Purpose : Physical constants and parameters linked 913 !! ** Purpose : Physical constants and parameters linked 909 914 !! to the mechanical ice redistribution 910 915 !! 911 !! ** Method : Read the namdyn_rdgrft namelist 912 !! and check the parameters values 916 !! ** Method : Read the namdyn_rdgrft namelist 917 !! and check the parameters values 913 918 !! called at the first timestep (nit000) 914 919 !! … … 920 925 & rn_csrdg , & 921 926 & ln_partf_lin, rn_gstar, & 922 & ln_partf_exp, rn_astar, & 923 & ln_ridging, rn_hstar, rn_porordg, rn_fsnwrdg, rn_fpndrdg, & 927 & ln_partf_exp, rn_astar, & 928 & ln_ridging, rn_hstar, rn_porordg, rn_fsnwrdg, rn_fpndrdg, & 924 929 & ln_rafting, rn_hraft, rn_craft , rn_fsnwrft, rn_fpndrft 925 930 !!------------------------------------------------------------------- … … 936 941 WRITE(numout,*) '~~~~~~~~~~~~~~~~~~' 937 942 WRITE(numout,*) ' Namelist namdyn_rdgrft:' 938 WRITE(numout,*) ' ice strength parameterization Hibler (1979) ln_str_H79 = ', ln_str_H79 943 WRITE(numout,*) ' ice strength parameterization Hibler (1979) ln_str_H79 = ', ln_str_H79 939 944 WRITE(numout,*) ' 1st bulk-rheology parameter rn_pstar = ', rn_pstar 940 945 WRITE(numout,*) ' 2nd bulk-rhelogy parameter rn_crhg = ', rn_crhg 941 WRITE(numout,*) ' Fraction of shear energy contributing to ridging rn_csrdg = ', rn_csrdg 946 WRITE(numout,*) ' Fraction of shear energy contributing to ridging rn_csrdg = ', rn_csrdg 942 947 WRITE(numout,*) ' linear ridging participation function ln_partf_lin = ', ln_partf_lin 943 948 WRITE(numout,*) ' Fraction of ice coverage contributing to ridging rn_gstar = ', rn_gstar … … 947 952 WRITE(numout,*) ' max ridged ice thickness rn_hstar = ', rn_hstar 948 953 WRITE(numout,*) ' Initial porosity of ridges rn_porordg = ', rn_porordg 949 WRITE(numout,*) ' Fraction of snow volume conserved during ridging rn_fsnwrdg = ', rn_fsnwrdg 950 WRITE(numout,*) ' Fraction of pond volume conserved during ridging rn_fpndrdg = ', rn_fpndrdg 954 WRITE(numout,*) ' Fraction of snow volume conserved during ridging rn_fsnwrdg = ', rn_fsnwrdg 955 WRITE(numout,*) ' Fraction of pond volume conserved during ridging rn_fpndrdg = ', rn_fpndrdg 951 956 WRITE(numout,*) ' Rafting of ice sheets or not ln_rafting = ', ln_rafting 952 957 WRITE(numout,*) ' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft 953 WRITE(numout,*) ' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft 954 WRITE(numout,*) ' Fraction of snow volume conserved during rafting rn_fsnwrft = ', rn_fsnwrft 955 WRITE(numout,*) ' Fraction of pond volume conserved during rafting rn_fpndrft = ', rn_fpndrft 958 WRITE(numout,*) ' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft 959 WRITE(numout,*) ' Fraction of snow volume conserved during rafting rn_fsnwrft = ', rn_fsnwrft 960 WRITE(numout,*) ' Fraction of pond volume conserved during rafting rn_fpndrft = ', rn_fpndrft 956 961 ENDIF 957 962 ! … … 967 972 WRITE(numout,*) ' ==> only ice dynamics is activated, thus some parameters must be changed' 968 973 WRITE(numout,*) ' rn_porordg = ', rn_porordg 969 WRITE(numout,*) ' rn_fsnwrdg = ', rn_fsnwrdg 970 WRITE(numout,*) ' rn_fpndrdg = ', rn_fpndrdg 971 WRITE(numout,*) ' rn_fsnwrft = ', rn_fsnwrft 972 WRITE(numout,*) ' rn_fpndrft = ', rn_fpndrft 974 WRITE(numout,*) ' rn_fsnwrdg = ', rn_fsnwrdg 975 WRITE(numout,*) ' rn_fpndrdg = ', rn_fpndrdg 976 WRITE(numout,*) ' rn_fsnwrft = ', rn_fsnwrft 977 WRITE(numout,*) ' rn_fpndrft = ', rn_fpndrft 973 978 ENDIF 974 979 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
