- Timestamp:
- 2017-10-05T16:44:46+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90
r8564 r8597 7 7 !! 1.0 ! 2012 (O. Lecomte) Adaptation for scientific tests (NEMO3.1) 8 8 !! 2.0 ! 2017 (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO3.6 9 !! ! NB: Only lim_mp_cstt and lim_mp_cesm work in this10 !! version11 9 !!---------------------------------------------------------------------- 12 10 #if defined key_lim3 … … 16 14 !! lim_mp_init : some initialization and namelist read 17 15 !! lim_mp : main calling routine 18 !! lim_mp_topo : main melt pond routine for the "topographic" formulation (FloccoFeltham)19 !! lim_mp_area : ??? compute melt pond fraction per category20 !! lim_mp_perm : computes permeability (should be a FUNCTION!)21 !! calc_hpond : computes melt pond depth22 !! permeability_phy : computes permeability23 24 16 !!---------------------------------------------------------------------- 25 17 USE phycst ! physical constants 26 18 USE dom_oce ! ocean space and time domain 27 ! USE sbc_ice ! Surface boundary condition: ice fields28 19 USE ice ! LIM-3 variables 29 ! USE icecons ! conservation tests30 ! USE limctl ! control prints31 ! USE limvar32 20 ! 33 21 USE lbclnk ! lateral boundary conditions - MPP exchanges … … 37 25 USE timing ! Timing 38 26 39 !OLI_CODE USE ice_oce, ONLY: rdt_ice, tatm_ice40 !OLI_CODE USE phycst41 !OLI_CODE USE dom_ice42 !OLI_CODE USE dom_oce43 !OLI_CODE USE sbc_oce44 !OLI_CODE USE sbc_ice45 !OLI_CODE USE par_ice46 !OLI_CODE USE par_oce47 !OLI_CODE USE ice48 !OLI_CODE USE thd_ice49 !OLI_CODE USE in_out_manager50 !OLI_CODE USE lbclnk51 !OLI_CODE USE lib_mpp52 !OLI_CODE53 !OLI_CODE IMPLICIT NONE54 !OLI_CODE PRIVATE55 !OLI_CODE56 !OLI_CODE PUBLIC lim_mp_init57 !OLI_CODE PUBLIC lim_mp58 59 27 IMPLICIT NONE 60 28 PRIVATE … … 62 30 PUBLIC lim_mp_init ! routine called by icestp.F90 63 31 PUBLIC lim_mp ! routine called by icestp.F90 32 33 INTEGER :: nice_pnd ! choice of the type of pond scheme 34 ! ! associated indices: 35 INTEGER, PARAMETER :: np_pndNO = 0 ! No pond scheme 36 INTEGER, PARAMETER :: np_pndCST = 1 ! Constant pond scheme 37 INTEGER, PARAMETER :: np_pndH12 = 2 ! Evolutive pond scheme (Holland et al. 2012) 64 38 65 39 !! * Substitutions … … 71 45 !!---------------------------------------------------------------------- 72 46 CONTAINS 73 74 SUBROUTINE lim_mp_init75 !!-------------------------------------------------------------------76 !! *** ROUTINE lim_mp_init ***77 !!78 !! ** Purpose : Physical constants and parameters linked to melt ponds79 !! over sea ice80 !!81 !! ** Method : Read the nammp namelist and check the melt pond82 !! parameter values called at the first timestep (nit000)83 !!84 !! ** input : Namelist nammp85 !!-------------------------------------------------------------------86 INTEGER :: ios ! Local integer output status for namelist read87 NAMELIST/nammp/ ln_pnd, ln_pnd_rad, ln_pnd_fw, nn_pnd_scheme, rn_apnd, rn_hpnd88 !!-------------------------------------------------------------------89 90 REWIND( numnam_ice_ref ) ! Namelist nammp in reference namelist : Melt Ponds91 READ ( numnam_ice_ref, nammp, IOSTAT = ios, ERR = 901)92 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammp in reference namelist', lwp )93 94 REWIND( numnam_ice_cfg ) ! Namelist nammp in configuration namelist : Melt Ponds95 READ ( numnam_ice_cfg, nammp, IOSTAT = ios, ERR = 902 )96 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammp in configuration namelist', lwp )97 IF(lwm) WRITE ( numoni, nammp )98 99 IF(lwp) THEN ! control print100 WRITE(numout,*)101 WRITE(numout,*) 'lim_mp_init : ice parameters for melt ponds'102 WRITE(numout,*) '~~~~~~~~~~~~'103 WRITE(numout,*) ' Active melt ponds ln_pnd = ', ln_pnd104 WRITE(numout,*) ' Active melt ponds radiative coupling ln_pnd_rad = ', ln_pnd_rad105 WRITE(numout,*) ' Active melt ponds freshwater coupling ln_pnd_fw = ', ln_pnd_fw106 WRITE(numout,*) ' Type of melt pond scheme =0 presc, =1 empirical = 2 topo nn_pnd_scheme = ', nn_pnd_scheme107 WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd108 WRITE(numout,*) ' Prescribed pond depth rn_hpnd = ', rn_hpnd109 ENDIF110 111 IF ( .NOT. ln_pnd ) THEN112 IF(lwp) WRITE(numout,*)113 IF(lwp) WRITE(numout,*) ' Melt ponds are not activated '114 IF(lwp) WRITE(numout,*) ' ln_pnd_rad and ln_pnd_fw set to .FALSE. '115 IF(lwp) WRITE(numout,*) ' nn_pnd_scheme, rn_apnd, rn_hpnd set to zero '116 ln_pnd_rad = .FALSE.117 ln_pnd_fw = .FALSE.118 nn_pnd_scheme = 0119 rn_apnd = 0._wp120 rn_hpnd = 0._wp121 122 IF(lwp) THEN ! control print123 WRITE(numout,*) ' Active melt ponds radiative coupling ln_pnd_rad = ', ln_pnd_rad124 WRITE(numout,*) ' Active melt ponds freshwater coupling ln_pnd_fw = ', ln_pnd_fw125 WRITE(numout,*) ' Type of melt pond scheme =0 presc, =1 empirical = 2 topo nn_pnd_scheme = ', nn_pnd_scheme126 WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd127 WRITE(numout,*) ' Prescribed pond depth rn_hpnd = ', rn_hpnd128 ENDIF129 ENDIF130 131 IF ( ln_pnd .AND. ( nn_pnd_scheme == 0 ) .AND. ( ln_pnd_fw ) ) THEN132 IF(lwp) WRITE(numout,*) ' Prescribed melt ponds do not conserve fresh water mass, hence ln_pnd_fw must be set to false '133 ln_pnd_fw = .FALSE.134 IF(lwp) THEN ! control print135 WRITE(numout,*) ' Active melt ponds freshwater coupling ln_pnd_fw = ', ln_pnd_fw136 ENDIF137 ENDIF138 139 IF ( ln_pnd .AND. ( nn_pnd_scheme == 2 ) .AND. ( jpl == 1 ) ) THEN140 IF(lwp) WRITE(numout,*) ' Topographic melt ponds are incompatible with jpl = 1 '141 IF(lwp) WRITE(numout,*) ' Run aborted '142 CALL ctl_stop( 'STOP', 'lim_mp_init: uncompatible options, reset namelist_ice_ref ' )143 ENDIF144 145 !146 END SUBROUTINE lim_mp_init147 148 47 149 48 SUBROUTINE lim_mp( kt ) … … 158 57 !! - 159 58 !!------------------------------------------------------------------------------------ 160 161 INTEGER, INTENT(in) :: kt ! number of iteration 162 INTEGER :: ji, jj, jl ! dummy loop indices 163 59 INTEGER, INTENT(in) :: kt ! number of iteration 60 INTEGER :: ji, jj, jl ! dummy loop indices 164 61 !!------------------------------------------------------------------- 165 62 166 63 IF( nn_timing == 1 ) CALL timing_start('lim_mp') 167 64 168 SELECT CASE ( nn_pnd_scheme ) 169 170 CASE (0) 171 172 CALL lim_mp_cstt ! staircase melt ponds 173 174 CASE (1) 175 176 CALL lim_mp_cesm ! empirical melt ponds 177 178 CASE (2) 179 180 CALL lim_mp_topo & ! topographic melt ponds 181 & (at_i, a_i, & 182 & vt_i, v_i, v_s, t_i, sz_i, a_ip_frac, & 183 & h_ip, t_su) 184 65 SELECT CASE ( nice_pnd ) 66 67 CASE (np_pndCST) 68 ! !------------------------------! 69 CALL lim_mp_CST ! Constant melt ponds ! 70 ! !------------------------------! 71 CASE (np_pndH12) 72 ! !------------------------------! 73 CALL lim_mp_H12 ! Holland et al 2012 melt ponds! 74 ! !------------------------------! 185 75 END SELECT 186 76 … … 189 79 END SUBROUTINE lim_mp 190 80 191 SUBROUTINE lim_mp_ cstt192 !!------------------------------------------------------------------- 193 !! *** ROUTINE lim_mp_ cstt***81 SUBROUTINE lim_mp_CST 82 !!------------------------------------------------------------------- 83 !! *** ROUTINE lim_mp_CST *** 194 84 !! 195 85 !! ** Purpose : Compute melt pond evolution … … 206 96 !! 207 97 !!------------------------------------------------------------------- 208 INTEGER :: ji, jj, jl209 REAL(wp) :: z1_jpl ! 1/jpl210 !!-------------------------------------------------------------------211 212 z1_jpl = 1. / FLOAT(jpl)213 98 214 99 WHERE ( ( a_i > epsi10 ) .AND. ( t_su >= rt0-epsi06 ) ) 100 !! WHERE ( ( a_i > 0._wp ) .AND. ( t_su >= rt0 ) ) 215 101 a_ip_frac = rn_apnd 216 102 h_ip = rn_hpnd … … 226 112 wfx_pnd(:,:) = 0._wp 227 113 228 END SUBROUTINE lim_mp_ cstt229 230 SUBROUTINE lim_mp_ cesm231 !!------------------------------------------------------------------- 232 !! *** ROUTINE lim_mp_ cesm***114 END SUBROUTINE lim_mp_CST 115 116 SUBROUTINE lim_mp_H12 117 !!------------------------------------------------------------------- 118 !! *** ROUTINE lim_mp_H12 *** 233 119 !! 234 120 !! ** Purpose : Compute melt pond evolution … … 249 135 !!------------------------------------------------------------------- 250 136 251 INTEGER, DIMENSION(jpij) :: indxi ! compressed indices for cells with ice melting252 INTEGER, DIMENSION(jpij) :: indxj !253 254 REAL(wp), DIMENSION(jpi,jpj) :: zwfx_mlw ! available meltwater for melt ponding255 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zrfrac ! fraction of available meltwater retained for melt ponding256 257 REAL(wp), PARAMETER :: zrmin = 0.15_wp ! minimum fraction of available meltwater retained for melt ponding258 REAL(wp), PARAMETER :: zrmax = 0.70_wp ! maximum '' '' '' '' ''259 REAL(wp), PARAMETER :: zrexp = 0.01_wp ! rate constant to refreeze melt ponds260 REAL(wp), PARAMETER :: zpnd_aspect = 0.8_wp ! pond aspect ratio261 262 REAL(wp) :: zhi ! dummy ice thickness263 REAL(wp) :: zhs ! dummy snow depth264 REAL(wp) :: zTp ! reference temperature265 REAL(wp) :: zdTs ! dummy temperature difference266 REAL(wp) :: z1_rhofw ! inverse freshwater density267 REAL(wp) :: z1_zpnd_aspect ! inverse pond aspect ratio268 REAL(wp) :: zvpold ! dummy pond volume269 270 INTEGER :: ji, jj, jl, ij ! loop indices271 INTEGER :: icells ! size of dummy array137 INTEGER, DIMENSION(jpij) :: indxi ! compressed indices for cells with ice melting 138 INTEGER, DIMENSION(jpij) :: indxj ! 139 140 REAL(wp), DIMENSION(jpi,jpj) :: zwfx_mlw ! available meltwater for melt ponding 141 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zrfrac ! fraction of available meltwater retained for melt ponding 142 143 REAL(wp), PARAMETER :: zrmin = 0.15_wp ! minimum fraction of available meltwater retained for melt ponding 144 REAL(wp), PARAMETER :: zrmax = 0.70_wp ! maximum '' '' '' '' '' 145 REAL(wp), PARAMETER :: zrexp = 0.01_wp ! rate constant to refreeze melt ponds 146 REAL(wp), PARAMETER :: zpnd_aspect = 0.8_wp ! pond aspect ratio 147 148 REAL(wp) :: zhi ! dummy ice thickness 149 REAL(wp) :: zhs ! dummy snow depth 150 REAL(wp) :: zTp ! reference temperature 151 REAL(wp) :: zdTs ! dummy temperature difference 152 REAL(wp) :: z1_rhofw ! inverse freshwater density 153 REAL(wp) :: z1_zpnd_aspect ! inverse pond aspect ratio 154 REAL(wp) :: zvpold ! dummy pond volume 155 156 INTEGER :: ji, jj, jl, ij ! loop indices 157 INTEGER :: icells ! size of dummy array 272 158 !!------------------------------------------------------------------- 273 159 z1_rhofw = 1. / rhofw … … 301 187 ! Identify grid cells where ponds should be updated (can probably be improved) 302 188 !------------------------------------------------------------------------------ 303 304 189 indxi(:) = 0 305 190 indxj(:) = 0 306 191 icells = 0 307 308 192 DO jj = 1, jpj 309 193 DO ji = 1, jpi … … 313 197 indxj(icells) = jj 314 198 ENDIF 315 END DO ! ji316 END DO ! jj199 END DO 200 END DO 317 201 318 202 DO ij = 1, icells … … 331 215 h_ip(ji,jj,jl) = 0._wp 332 216 333 IF ( ln_pnd_fw ) & !--- Give freshwater to the ocean217 IF ( ln_pnd_fwb ) & !--- Give freshwater to the ocean 334 218 wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + v_ip(ji,jj,jl) 335 219 … … 352 236 !--- Dump meltwater due to refreezing ( of course this is wrong 353 237 !--- but this parameterization is too simple ) 354 IF ( ln_pnd_fw ) &238 IF ( ln_pnd_fwb ) & 355 239 wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + rhofw * ( v_ip(ji,jj,jl) - zvpold ) * r1_rdtice 356 240 … … 361 245 362 246 a_ip(ji,jj,jl) = a_ip_frac(ji,jj,jl) * a_i(ji,jj,jl) 363 364 !-----------------------------------------------------------365 ! Limit pond depth366 !-----------------------------------------------------------367 ! The original version has pond depth limitation, which I did not368 ! keep here. Maybe we need it later on369 !370 247 371 248 ENDIF … … 377 254 !--- Remove retained meltwater from surface fluxes 378 255 379 IF ( ln_pnd_fw ) THEN 380 381 wfx_snw_sum(:,:) = wfx_snw_sum(:,:) * ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) ) 382 383 wfx_sum(:,:) = wfx_sum(:,:) * ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) ) 256 IF ( ln_pnd_fwb ) THEN 257 258 wfx_snw_sum(:,:) = wfx_snw_sum(:,:) * ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) ) 259 wfx_sum(:,:) = wfx_sum(:,:) * ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) ) 384 260 385 261 ENDIF 386 262 387 END SUBROUTINE lim_mp_cesm 388 389 SUBROUTINE lim_mp_topo (aice, aicen, & 390 vice, vicen, & 391 vsnon, & 392 ticen, salin, & 393 a_ip_frac, h_ip, & 394 Tsfc ) 395 !!------------------------------------------------------------------- 396 !! *** ROUTINE lim_mp_topo *** 397 !! 398 !! ** Purpose : Compute melt pond evolution based on the ice 399 !! topography as inferred from the ice thickness 400 !! distribution. 401 !! 402 !! ** Method : This code is initially based on Flocco and Feltham 403 !! (2007) and Flocco et al. (2010). More to come... 404 !! 405 !! ** Tunable parameters : 406 !! 407 !! ** Note : 408 !! 409 !! ** References 410 !! Flocco, D. and D. L. Feltham, 2007. A continuum model of melt pond 411 !! evolution on Arctic sea ice. J. Geophys. Res. 112, C08016, doi: 412 !! 10.1029/2006JC003836. 413 !! Flocco, D., D. L. Feltham and A. K. Turner, 2010. Incorporation of 414 !! a physically based melt pond scheme into the sea ice component of a 415 !! climate model. J. Geophys. Res. 115, C08012, 416 !! doi: 10.1029/2009JC005568. 417 !! 418 !!------------------------------------------------------------------- 419 420 REAL (wp), DIMENSION (jpi,jpj), & 421 INTENT(IN) :: & 422 aice, & ! total ice area fraction 423 vice ! total ice volume (m) 424 425 REAL (wp), DIMENSION (jpi,jpj,jpl), & 426 INTENT(IN) :: & 427 aicen, & ! ice area fraction, per category 428 vsnon, & ! snow volume, per category (m) 429 vicen ! ice volume, per category (m) 430 431 REAL (wp), DIMENSION (jpi,jpj,nlay_i,jpl), & 432 INTENT(IN) :: & 433 ticen, & ! ice enthalpy, per category 434 salin 435 436 REAL (wp), DIMENSION (jpi,jpj,jpl), & 437 INTENT(INOUT) :: & 438 a_ip_frac , & ! pond area fraction of ice, per ice category 439 h_ip ! pond depth, per ice category (m) 440 441 REAL (wp), DIMENSION (jpi,jpj,jpl), & 442 INTENT(IN) :: & 443 Tsfc ! snow/sea ice surface temperature 444 445 ! local variables 446 REAL (wp), DIMENSION (jpi,jpj,jpl) :: & 447 zTsfcn, & ! ice/snow surface temperature (C) 448 zvolpn, & ! pond volume per unit area, per category (m) 449 zvuin ! water-equivalent volume of ice lid on melt pond ('upper ice', m) 450 451 REAL (wp), DIMENSION (jpi,jpj,jpl) :: & 452 zapondn,& ! pond area fraction, per category 453 zhpondn ! pond depth, per category (m) 454 455 REAL (wp), DIMENSION (jpi,jpj) :: & 456 zvolp ! total volume of pond, per unit area of pond (m) 457 458 REAL (wp) :: & 459 zhi, & ! ice thickness (m) 460 zdHui, & ! change in thickness of ice lid (m) 461 zomega, & ! conduction 462 zdTice, & ! temperature difference across ice lid (C) 463 zdvice, & ! change in ice volume (m) 464 zTavg, & ! mean surface temperature across categories (C) 465 zTp, & ! pond freezing temperature (C) 466 zdvn ! change in melt pond volume for fresh water budget 467 INTEGER, DIMENSION (jpi*jpj) :: & 468 indxi, indxj ! compressed indices for cells with ice melting 469 470 INTEGER :: n,k,i,j,ij,icells,indxij ! loop indices 471 472 INTEGER, DIMENSION (jpl) :: & 473 kcells ! cells where ice lid combines with vice 474 475 INTEGER, DIMENSION (jpi*jpj,jpl) :: & 476 indxii, indxjj ! i,j indices for kcells loop 477 478 REAL (wp), parameter :: & 479 zhicemin = 0.1_wp , & ! minimum ice thickness with ponds (m) 480 zTd = 0.15_wp, & ! temperature difference for freeze-up (C) 481 zr1_rlfus = 1._wp / 0.334e+6 / 917._wp , & ! (J/m^3) 482 zmin_volp = 1.e-4_wp, & ! minimum pond volume (m) 483 z0 = 0._wp, & 484 zTimelt = 0._wp, & 485 z01 = 0.01_wp, & 486 z25 = 0.25_wp, & 487 z5 = 0.5_wp 488 489 !--------------------------------------------------------------- 490 ! Initialization 491 !--------------------------------------------------------------- 492 493 zhpondn(:,:,:) = 0._wp 494 zapondn(:,:,:) = 0._wp 495 indxii(:,:) = 0 496 indxjj(:,:) = 0 497 kcells(:) = 0 498 499 zvolp(:,:) = wfx_sum(:,:) + wfx_snw_sum(:,:) + vt_ip(:,:) ! Total available melt water, to be distributed as melt ponds 500 zTsfcn(:,:,:) = zTsfcn(:,:,:) - rt0 ! Convert in Celsius 501 502 ! The freezing temperature for meltponds is assumed slightly below 0C, 503 ! as if meltponds had a little salt in them. The salt budget is not 504 ! altered for meltponds, but if it were then an actual pond freezing 505 ! temperature could be computed. 506 507 ! zTp = zTimelt - zTd ---> for lids 508 509 !----------------------------------------------------------------- 510 ! Identify grid cells with ponds 511 !----------------------------------------------------------------- 512 513 icells = 0 514 DO j = 1, jpj 515 DO i = 1, jpi 516 zhi = z0 517 IF (aice(i,j) > epsi10 ) zhi = vice(i,j)/aice(i,j) 518 IF ( aice(i,j) > z01 .and. zhi > zhicemin .and. & 519 zvolp(i,j) > zmin_volp*aice(i,j)) THEN 520 icells = icells + 1 521 indxi(icells) = i 522 indxj(icells) = j 523 ELSE ! remove ponds on thin ice 524 !fpond(i,j) = fpond(i,j) - zvolp(i,j) 525 zvolpn(i,j,:) = z0 526 zvuin (i,j,:) = z0 527 zvolp (i,j) = z0 528 END IF 529 END DO ! i 530 END DO ! j 531 532 DO ij = 1, icells 533 i = indxi(ij) 534 j = indxj(ij) 535 536 !-------------------------------------------------------------- 537 ! calculate pond area and depth 538 !-------------------------------------------------------------- 539 CALL lim_mp_area(aice(i,j),vice(i,j), & 540 aicen(i,j,:), vicen(i,j,:), vsnon(i,j,:), & 541 ticen(i,j,:,:), salin(i,j,:,:), & 542 zvolpn(i,j,:), zvolp(i,j), & 543 zapondn(i,j,:),zhpondn(i,j,:), zdvn) 544 ! outputs are 545 ! - zdvn 546 ! - zvolpn 547 ! - zvolp 548 ! - zapondn 549 ! - zhpondn 550 551 wfx_pnd(i,j) = wfx_pnd(i,j) + zdvn ! update flux from ponds to ocean 552 553 ! mean surface temperature MV - why do we need that ? --> for the lid 554 555 ! zTavg = z0 556 ! DO n = 1, jpl 557 ! zTavg = zTavg + zTsfcn(i,j,n)*aicen(i,j,n) 558 ! END DO 559 ! zTavg = zTavg / aice(i,j) 560 561 END DO ! ij 562 563 !--------------------------------------------------------------- 564 ! Update pond volume and fraction 565 !--------------------------------------------------------------- 566 567 a_ip(:,:,:) = zapondn(:,:,:) 568 v_ip(:,:,:) = zapondn(:,:,:) * zhpondn(:,:,:) 569 a_ip_frac(:,:,:) = 0._wp 570 h_ip (:,:,:) = 0._wp 571 572 END SUBROUTINE lim_mp_topo 573 574 SUBROUTINE lim_mp_area(aice,vice, & 575 aicen, vicen, vsnon, ticen, & 576 salin, zvolpn, zvolp, & 577 zapondn,zhpondn,dvolp) 578 579 !!------------------------------------------------------------------- 580 !! *** ROUTINE lim_mp_area *** 581 !! 582 !! ** Purpose : Given the total volume of meltwater, update 583 !! pond fraction (a_ip) and depth (should be volume) 584 !! 585 !! ** 586 !! 587 !!------------------------------------------------------------------ 588 589 REAL (wp), INTENT(IN) :: & 590 aice,vice 591 592 REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 593 aicen, vicen, vsnon 594 595 REAL (wp), DIMENSION(nlay_i,jpl), INTENT(IN) :: & 596 ticen, salin 597 598 REAL (wp), DIMENSION(jpl), INTENT(INOUT) :: & 599 zvolpn 600 601 REAL (wp), INTENT(INOUT) :: & 602 zvolp, dvolp 603 604 REAL (wp), DIMENSION(jpl), INTENT(OUT) :: & 605 zapondn, zhpondn 606 607 INTEGER :: & 608 n, ns, & 609 m_index, & 610 permflag 611 612 REAL (wp), DIMENSION(jpl) :: & 613 hicen, & 614 hsnon, & 615 asnon, & 616 alfan, & 617 betan, & 618 cum_max_vol, & 619 reduced_aicen 620 621 REAL (wp), DIMENSION(0:jpl) :: & 622 cum_max_vol_tmp 623 624 REAL (wp) :: & 625 hpond, & 626 drain, & 627 floe_weight, & 628 pressure_head, & 629 hsl_rel, & 630 deltah, & 631 perm, & 632 msno 633 634 REAL (wp), parameter :: & 635 viscosity = 1.79e-3_wp, & ! kinematic water viscosity in kg/m/s 636 z0 = 0.0_wp , & 637 c1 = 1.0_wp , & 638 p4 = 0.4_wp , & 639 p6 = 0.6_wp , & 640 epsi10 = 1.0e-11_wp 641 642 !-----------| 643 ! | 644 ! |-----------| 645 !___________|___________|______________________________________sea-level 646 ! | | 647 ! | |---^--------| 648 ! | | | | 649 ! | | | |-----------| |------- 650 ! | | |alfan(n)| | | 651 ! | | | | |--------------| 652 ! | | | | | | 653 !---------------------------v------------------------------------------- 654 ! | | ^ | | | 655 ! | | | | |--------------| 656 ! | | |betan(n)| | | 657 ! | | | |-----------| |------- 658 ! | | | | 659 ! | |---v------- | 660 ! | | 661 ! |-----------| 662 ! | 663 !-----------| 664 665 !------------------------------------------------------------------- 666 ! initialize 667 !------------------------------------------------------------------- 668 669 DO n = 1, jpl 670 671 zapondn(n) = z0 672 zhpondn(n) = z0 673 674 !---------------------------------------- 675 ! X) compute the effective snow fraction 676 !---------------------------------------- 677 IF (aicen(n) < epsi10) THEN 678 hicen(n) = z0 679 hsnon(n) = z0 680 reduced_aicen(n) = z0 681 ELSE 682 hicen(n) = vicen(n) / aicen(n) 683 hsnon(n) = vsnon(n) / aicen(n) 684 reduced_aicen(n) = c1 ! n=jpl 685 IF (n < jpl) reduced_aicen(n) = aicen(n) & 686 * (-0.024_wp*hicen(n) + 0.832_wp) 687 asnon(n) = reduced_aicen(n) ! effective snow fraction (empirical) 688 ! MV should check whether this makes sense to have the same effective snow fraction in here 689 END IF 690 691 ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 692 ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming 693 ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all 694 ! categories. alfa and beta partition the ITD - they are areas not thicknesses! 695 ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. 696 ! Here, alfa = 60% of the ice area (and since hice is constant in a category, 697 ! alfan = 60% of the ice volume) in each category lies above the reference line, 698 ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 699 700 ! MV: 701 ! Note that this choice is not in the original FF07 paper and has been adopted in CICE 702 ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe 703 704 ! Where does that choice come from 705 706 alfan(n) = 0.6 * hicen(n) 707 betan(n) = 0.4 * hicen(n) 708 709 cum_max_vol(n) = z0 710 cum_max_vol_tmp(n) = z0 711 712 END DO ! jpl 713 714 cum_max_vol_tmp(0) = z0 715 drain = z0 716 dvolp = z0 717 718 !---------------------------------------------------------- 719 ! x) Drain overflow water, update pond fraction and volume 720 !---------------------------------------------------------- 721 722 !-------------------------------------------------------------------------- 723 ! the maximum amount of water that can be contained up to each ice category 724 !-------------------------------------------------------------------------- 725 726 ! MV 727 ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) 728 ! Then the excess volume cum_max_vol(jl) drains out of the system 729 ! It should be added to wfx_pnd 730 ! END MV 731 732 DO n = 1, jpl-1 ! last category can not hold any volume 733 734 IF (alfan(n+1) >= alfan(n) .and. alfan(n+1) > z0) THEN 735 736 ! total volume in level including snow 737 cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) + & 738 (alfan(n+1) - alfan(n)) * sum(reduced_aicen(1:n)) 739 740 ! subtract snow solid volumes from lower categories in current level 741 DO ns = 1, n 742 cum_max_vol_tmp(n) = cum_max_vol_tmp(n) & 743 - rhosn/rhofw * & ! free air fraction that can be filled by water 744 asnon(ns) * & ! effective areal fraction of snow in that category 745 max(min(hsnon(ns)+alfan(ns)-alfan(n), alfan(n+1)- & 746 alfan(n)), z0) 747 END DO 748 749 ELSE ! assume higher categories unoccupied 750 cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) 751 END IF 752 !IF (cum_max_vol_tmp(n) < z0) THEN 753 ! call abort_ice('negative melt pond volume') 754 !END IF 755 END DO 756 cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume 757 cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl) 758 759 !---------------------------------------------------------------- 760 ! is there more meltwater than can be held in the floe? 761 !---------------------------------------------------------------- 762 IF (zvolp >= cum_max_vol(jpl)) THEN 763 drain = zvolp - cum_max_vol(jpl) + epsi10 764 zvolp = zvolp - drain ! update meltwater volume available 765 dvolp = drain ! this is the drained water 766 IF (zvolp < epsi10) THEN 767 dvolp = dvolp + zvolp 768 zvolp = z0 769 END IF 770 END IF 771 772 ! height and area corresponding to the remaining volume 773 774 ! call calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, & 775 ! zvolp, cum_max_vol, hpond, m_index) 776 777 DO n=1, m_index 778 zhpondn(n) = hpond - alfan(n) + alfan(1) ! here oui choulde update 779 ! volume instead, no ? 780 zapondn(n) = reduced_aicen(n) 781 ! in practise, pond fraction depends on the empirical snow fraction 782 ! so in turn on ice thickness 783 END DO 784 785 !------------------------------------------------------------------------ 786 ! Drainage through brine network (permeability) 787 !------------------------------------------------------------------------ 788 !!! drainage due to ice permeability - Darcy's law 789 790 ! sea water level 791 msno = z0 792 DO n=1,jpl 793 msno = msno + vsnon(n) * rhosn 794 END DO 795 floe_weight = (msno + rhoic*vice + rau0*zvolp) / aice 796 hsl_rel = floe_weight / rau0 & 797 - ((sum(betan(:)*aicen(:))/aice) + alfan(1)) 798 799 deltah = hpond - hsl_rel 800 pressure_head = grav * rau0 * max(deltah, z0) 801 802 ! drain IF ice is permeable 803 permflag = 0 804 IF (pressure_head > z0) THEN 805 DO n = 1, jpl-1 806 IF (hicen(n) /= z0) THEN 807 perm = 0. ! MV ugly dummy patch 808 CALL lim_mp_perm(ticen(:,n), salin(:,n), vicen(n), perm) 809 IF (perm > z0) permflag = 1 810 811 drain = perm*zapondn(n)*pressure_head*rdt_ice / & 812 (viscosity*hicen(n)) 813 dvolp = dvolp + min(drain, zvolp) 814 zvolp = max(zvolp - drain, z0) 815 IF (zvolp < epsi10) THEN 816 dvolp = dvolp + zvolp 817 zvolp = z0 818 END IF 819 END IF 820 END DO 821 822 ! adjust melt pond DIMENSIONs 823 IF (permflag > 0) THEN 824 ! recompute pond depth 825 ! CALL calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, & 826 ! zvolp, cum_max_vol, hpond, m_index) 827 DO n=1, m_index 828 zhpondn(n) = hpond - alfan(n) + alfan(1) 829 zapondn(n) = reduced_aicen(n) 830 END DO 831 END IF 832 END IF ! pressure_head 833 834 !------------------------------- 835 ! X) remove water from the snow 836 !------------------------------- 837 !------------------------------------------------------------------------ 838 ! total melt pond volume in category DOes not include snow volume 839 ! snow in melt ponds is not melted 840 !------------------------------------------------------------------------ 841 842 ! Calculate pond volume for lower categories 843 DO n=1,m_index-1 844 zvolpn(n) = zapondn(n) * zhpondn(n) & ! what is not in the snow 845 - (rhosn/rhofw) * asnon(n) * min(hsnon(n), zhpondn(n)) 846 END DO 847 848 ! Calculate pond volume for highest category = remaining pond volume 849 850 ! The following is completely unclear to Martin at least 851 ! Could we redefine properly and recode in a more readable way ? 852 853 ! m_index = last category with melt pond 854 855 IF (m_index == 1) zvolpn(m_index) = zvolp ! volume of mw in 1st category is the total volume of melt water 856 857 IF (m_index > 1) THEN 858 IF (zvolp > sum(zvolpn(1:m_index-1))) THEN 859 zvolpn(m_index) = zvolp - sum(zvolpn(1:m_index-1)) ! 860 ELSE 861 zvolpn(m_index) = z0 862 zhpondn(m_index) = z0 863 zapondn(m_index) = z0 864 ! If remaining pond volume is negative reduce pond volume of 865 ! lower category 866 IF (zvolp+epsi10 < sum(zvolpn(1:m_index-1))) & 867 zvolpn(m_index-1) = zvolpn(m_index-1)-sum(zvolpn(1:m_index-1))& 868 + zvolp 869 END IF 870 END IF 871 872 DO n=1,m_index 873 IF (zapondn(n) > epsi10) THEN 874 zhpondn(n) = zvolpn(n) / zapondn(n) 875 ELSE 876 dvolp = dvolp + zvolpn(n) 877 zhpondn(n) = z0 878 zvolpn(n) = z0 879 zapondn(n) = z0 880 end IF 881 END DO 882 DO n = m_index+1, jpl 883 zhpondn(n) = z0 884 zapondn(n) = z0 885 zvolpn (n) = z0 886 END DO 887 888 END SUBROUTINE lim_mp_area 889 890 !OLI_CODE 891 !OLI_CODE 892 !OLI_CODE SUBROUTINE calc_hpond(aicen, asnon, hsnon, rhos, alfan, & 893 !OLI_CODE zvolp, cum_max_vol, & 894 !OLI_CODE hpond, m_index) 895 !OLI_CODE !!------------------------------------------------------------------- 896 !OLI_CODE !! *** ROUTINE calc_hpond *** 897 !OLI_CODE !! 898 !OLI_CODE !! ** Purpose : Compute melt pond depth 899 !OLI_CODE !!------------------------------------------------------------------- 900 !OLI_CODE 901 !OLI_CODE REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 902 !OLI_CODE aicen, & 903 !OLI_CODE asnon, & 904 !OLI_CODE hsnon, & 905 !OLI_CODE rhos, & 906 !OLI_CODE alfan, & 907 !OLI_CODE cum_max_vol 908 !OLI_CODE 909 !OLI_CODE REAL (wp), INTENT(IN) :: & 910 !OLI_CODE zvolp 911 !OLI_CODE 912 !OLI_CODE REAL (wp), INTENT(OUT) :: & 913 !OLI_CODE hpond 914 !OLI_CODE 915 !OLI_CODE INTEGER, INTENT(OUT) :: & 916 !OLI_CODE m_index 917 !OLI_CODE 918 !OLI_CODE INTEGER :: n, ns 919 !OLI_CODE 920 !OLI_CODE REAL (wp), DIMENSION(0:jpl+1) :: & 921 !OLI_CODE hitl, & 922 !OLI_CODE aicetl 923 !OLI_CODE 924 !OLI_CODE REAL (wp) :: & 925 !OLI_CODE rem_vol, & 926 !OLI_CODE area, & 927 !OLI_CODE vol, & 928 !OLI_CODE tmp, & 929 !OLI_CODE z0 = 0.0_wp, & 930 !OLI_CODE epsi10 = 1.0e-11_wp 931 !OLI_CODE 932 !OLI_CODE !---------------------------------------------------------------- 933 !OLI_CODE ! hpond is zero if zvolp is zero - have we fully drained? 934 !OLI_CODE !---------------------------------------------------------------- 935 !OLI_CODE 936 !OLI_CODE IF (zvolp < epsi10) THEN 937 !OLI_CODE hpond = z0 938 !OLI_CODE m_index = 0 939 !OLI_CODE ELSE 940 !OLI_CODE 941 !OLI_CODE !---------------------------------------------------------------- 942 !OLI_CODE ! Calculate the category where water fills up to 943 !OLI_CODE !---------------------------------------------------------------- 944 !OLI_CODE 945 !OLI_CODE !----------| 946 !OLI_CODE ! | 947 !OLI_CODE ! | 948 !OLI_CODE ! |----------| -- -- 949 !OLI_CODE !__________|__________|_________________________________________ ^ 950 !OLI_CODE ! | | rem_vol ^ | Semi-filled 951 !OLI_CODE ! | |----------|-- -- -- - ---|-- ---- -- -- --v layer 952 !OLI_CODE ! | | | | 953 !OLI_CODE ! | | | |hpond 954 !OLI_CODE ! | | |----------| | |------- 955 !OLI_CODE ! | | | | | | 956 !OLI_CODE ! | | | |---v-----| 957 !OLI_CODE ! | | m_index | | | 958 !OLI_CODE !------------------------------------------------------------- 959 !OLI_CODE 960 !OLI_CODE m_index = 0 ! 1:m_index categories have water in them 961 !OLI_CODE DO n = 1, jpl 962 !OLI_CODE IF (zvolp <= cum_max_vol(n)) THEN 963 !OLI_CODE m_index = n 964 !OLI_CODE IF (n == 1) THEN 965 !OLI_CODE rem_vol = zvolp 966 !OLI_CODE ELSE 967 !OLI_CODE rem_vol = zvolp - cum_max_vol(n-1) 968 !OLI_CODE END IF 969 !OLI_CODE exit ! to break out of the loop 970 !OLI_CODE END IF 971 !OLI_CODE END DO 972 !OLI_CODE m_index = min(jpl-1, m_index) 973 !OLI_CODE 974 !OLI_CODE !---------------------------------------------------------------- 975 !OLI_CODE ! semi-filled layer may have m_index different snow in it 976 !OLI_CODE !---------------------------------------------------------------- 977 !OLI_CODE 978 !OLI_CODE !----------------------------------------------------------- ^ 979 !OLI_CODE ! | alfan(m_index+1) 980 !OLI_CODE ! | 981 !OLI_CODE !hitl(3)--> |----------| | 982 !OLI_CODE !hitl(2)--> |------------| * * * * *| | 983 !OLI_CODE !hitl(1)--> |----------|* * * * * * |* * * * * | | 984 !OLI_CODE !hitl(0)-->------------------------------------------------- | ^ 985 !OLI_CODE ! various snow from lower categories | |alfa(m_index) 986 !OLI_CODE 987 !OLI_CODE ! hitl - heights of the snow layers from thinner and current categories 988 !OLI_CODE ! aicetl - area of each snow depth in this layer 989 !OLI_CODE 990 !OLI_CODE hitl(:) = z0 991 !OLI_CODE aicetl(:) = z0 992 !OLI_CODE DO n = 1, m_index 993 !OLI_CODE hitl(n) = max(min(hsnon(n) + alfan(n) - alfan(m_index), & 994 !OLI_CODE alfan(m_index+1) - alfan(m_index)), z0) 995 !OLI_CODE aicetl(n) = asnon(n) 996 !OLI_CODE 997 !OLI_CODE aicetl(0) = aicetl(0) + (aicen(n) - asnon(n)) 998 !OLI_CODE END DO 999 !OLI_CODE hitl(m_index+1) = alfan(m_index+1) - alfan(m_index) 1000 !OLI_CODE aicetl(m_index+1) = z0 1001 !OLI_CODE 1002 !OLI_CODE !---------------------------------------------------------------- 1003 !OLI_CODE ! reorder array according to hitl 1004 !OLI_CODE ! snow heights not necessarily in height order 1005 !OLI_CODE !---------------------------------------------------------------- 1006 !OLI_CODE 1007 !OLI_CODE DO ns = 1, m_index+1 1008 !OLI_CODE DO n = 0, m_index - ns + 1 1009 !OLI_CODE IF (hitl(n) > hitl(n+1)) THEN ! swap order 1010 !OLI_CODE tmp = hitl(n) 1011 !OLI_CODE hitl(n) = hitl(n+1) 1012 !OLI_CODE hitl(n+1) = tmp 1013 !OLI_CODE tmp = aicetl(n) 1014 !OLI_CODE aicetl(n) = aicetl(n+1) 1015 !OLI_CODE aicetl(n+1) = tmp 1016 !OLI_CODE END IF 1017 !OLI_CODE END DO 1018 !OLI_CODE END DO 1019 !OLI_CODE 1020 !OLI_CODE !---------------------------------------------------------------- 1021 !OLI_CODE ! divide semi-filled layer into set of sublayers each vertically homogenous 1022 !OLI_CODE !---------------------------------------------------------------- 1023 !OLI_CODE 1024 !OLI_CODE !hitl(3)---------------------------------------------------------------- 1025 !OLI_CODE ! | * * * * * * * * 1026 !OLI_CODE ! |* * * * * * * * * 1027 !OLI_CODE !hitl(2)---------------------------------------------------------------- 1028 !OLI_CODE ! | * * * * * * * * | * * * * * * * * 1029 !OLI_CODE ! |* * * * * * * * * |* * * * * * * * * 1030 !OLI_CODE !hitl(1)---------------------------------------------------------------- 1031 !OLI_CODE ! | * * * * * * * * | * * * * * * * * | * * * * * * * * 1032 !OLI_CODE ! |* * * * * * * * * |* * * * * * * * * |* * * * * * * * * 1033 !OLI_CODE !hitl(0)---------------------------------------------------------------- 1034 !OLI_CODE ! aicetl(0) aicetl(1) aicetl(2) aicetl(3) 1035 !OLI_CODE 1036 !OLI_CODE ! move up over layers incrementing volume 1037 !OLI_CODE DO n = 1, m_index+1 1038 !OLI_CODE 1039 !OLI_CODE area = sum(aicetl(:)) - & ! total area of sub-layer 1040 !OLI_CODE (rhos(n)/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow 1041 !OLI_CODE 1042 !OLI_CODE vol = (hitl(n) - hitl(n-1)) * area ! thickness of sub-layer times area 1043 !OLI_CODE 1044 !OLI_CODE IF (vol >= rem_vol) THEN ! have reached the sub-layer with the depth within 1045 !OLI_CODE hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - & 1046 !OLI_CODE alfan(1) 1047 !OLI_CODE exit 1048 !OLI_CODE ELSE ! still in sub-layer below the sub-layer with the depth 1049 !OLI_CODE rem_vol = rem_vol - vol 1050 !OLI_CODE END IF 1051 !OLI_CODE 1052 !OLI_CODE END DO 1053 !OLI_CODE 1054 !OLI_CODE END IF 1055 !OLI_CODE 1056 !OLI_CODE END SUBROUTINE calc_hpond 1057 !OLI_CODE 1058 !OLI_CODE 1059 SUBROUTINE lim_mp_perm(ticen, salin, vicen, perm) 1060 !!------------------------------------------------------------------- 1061 !! *** ROUTINE lim_mp_perm *** 1062 !! 1063 !! ** Purpose : Determine the liquid fraction of brine in the ice 1064 !! and its permeability 1065 !!------------------------------------------------------------------- 1066 REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: & 1067 ticen, & ! energy of melting for each ice layer (J/m2) 1068 salin 1069 1070 REAL (wp), INTENT(IN) :: & 1071 vicen ! ice volume 1072 1073 REAL (wp), INTENT(OUT) :: & 1074 perm ! permeability 1075 1076 REAL (wp) :: & 1077 Sbr ! brine salinity 1078 1079 REAL (wp), DIMENSION(nlay_i) :: & 1080 Tin, & ! ice temperature 1081 phi ! liquid fraction 1082 1083 INTEGER :: k 1084 1085 REAL (wp) :: & 1086 c2 = 2.0_wp 1087 1088 !----------------------------------------------------------------- 1089 ! Compute ice temperatures from enthalpies using quadratic formula 1090 !----------------------------------------------------------------- 1091 1092 DO k = 1,nlay_i 1093 Tin(k) = ticen(k) 1094 END DO 1095 1096 !----------------------------------------------------------------- 1097 ! brine salinity and liquid fraction 1098 !----------------------------------------------------------------- 1099 1100 IF (maxval(Tin-rtt) <= -c2) THEN 1101 1102 DO k = 1,nlay_i 1103 Sbr = - 1.2_wp & 1104 -21.8_wp * (Tin(k)-rtt) & 1105 - 0.919_wp * (Tin(k)-rtt)**2 & 1106 - 0.01878_wp * (Tin(k)-rtt)**3 1107 phi(k) = salin(k)/Sbr ! liquid fraction 1108 END DO ! k 1109 1110 ELSE 1111 1112 DO k = 1,nlay_i 1113 Sbr = -17.6_wp * (Tin(k)-rtt) & 1114 - 0.389_wp * (Tin(k)-rtt)**2 & 1115 - 0.00362_wp* (Tin(k)-rtt)**3 1116 phi(k) = salin(k)/Sbr ! liquid fraction 1117 END DO 1118 1119 END IF 1120 1121 !----------------------------------------------------------------- 1122 ! permeability 1123 !----------------------------------------------------------------- 1124 1125 perm = 3.0e-08_wp * (minval(phi))**3 ! REFERENCE PLEASE (this fucking 1126 ! bastard of Golden) 1127 1128 END SUBROUTINE lim_mp_perm 1129 !OLI_CODE 1130 !OLI_CODE #else 1131 !OLI_CODE !!---------------------------------------------------------------------- 1132 !OLI_CODE !! Default option Dummy Module No LIM-3 sea-ice model 1133 !OLI_CODE !!---------------------------------------------------------------------- 1134 !OLI_CODE CONTAINS 1135 !OLI_CODE SUBROUTINE lim_mp_init ! Empty routine 1136 !OLI_CODE END SUBROUTINE lim_mp_init 1137 !OLI_CODE SUBROUTINE lim_mp ! Empty routine 1138 !OLI_CODE END SUBROUTINE lim_mp 1139 !OLI_CODE SUBROUTINE compute_mp_topo ! Empty routine 1140 !OLI_CODE END SUBROUTINE compute_mp_topo 1141 !OLI_CODE SUBROUTINE pond_area ! Empty routine 1142 !OLI_CODE END SUBROUTINE pond_area 1143 !OLI_CODE SUBROUTINE calc_hpond ! Empty routine 1144 !OLI_CODE END SUBROUTINE calc_hpond 1145 !OLI_CODE SUBROUTINE permeability_phy ! Empty routine 1146 !OLI_CODE END SUBROUTINE permeability_phy 1147 !OLI_CODE #endif 1148 !OLI_CODE !!====================================================================== 1149 !OLI_CODE END MODULE limmp_topo 1150 1151 1152 !OLI_CODE MODULE limmp_topo 1153 !OLI_CODE !!====================================================================== 1154 !OLI_CODE !! *** MODULE limmp_topo *** 1155 !OLI_CODE !! LIM-3 sea-ice : computation of melt ponds' properties 1156 !OLI_CODE !!====================================================================== 1157 !OLI_CODE !! History : Original code by Daniela Flocco and Adrian Turner 1158 !OLI_CODE !! ! 2012-09 (O. Lecomte) Adaptation for routine inclusion in 1159 !OLI_CODE !! NEMO-LIM3.1 1160 !OLI_CODE !! ! 2016-11 (O. Lecomte, C. Rousset, M. Vancoppenolle) 1161 !OLI_CODE !! Adaptation for merge with NEMO-LIM3.6 1162 !OLI_CODE !!---------------------------------------------------------------------- 1163 !OLI_CODE #if defined key_lim3 1164 !OLI_CODE !!---------------------------------------------------------------------- 1165 !OLI_CODE !! 'key_lim3' LIM-3 sea-ice model 1166 !OLI_CODE !!---------------------------------------------------------------------- 1167 !OLI_CODE !! lim_mp_init : melt pond properties initialization 1168 !OLI_CODE !! lim_mp : melt pond routine caller 1169 !OLI_CODE !! compute_mp_topo : Actual melt pond routine 1170 !OLI_CODE !!---------------------------------------------------------------------- 1171 !OLI_CODE USE ice_oce, ONLY: rdt_ice, tatm_ice 1172 !OLI_CODE USE phycst 1173 !OLI_CODE USE dom_ice 1174 !OLI_CODE USE dom_oce 1175 !OLI_CODE USE sbc_oce 1176 !OLI_CODE USE sbc_ice 1177 !OLI_CODE USE par_ice 1178 !OLI_CODE USE par_oce 1179 !OLI_CODE USE ice 1180 !OLI_CODE USE thd_ice 1181 !OLI_CODE USE in_out_manager 1182 !OLI_CODE USE lbclnk 1183 !OLI_CODE USE lib_mpp 1184 !OLI_CODE 1185 !OLI_CODE IMPLICIT NONE 1186 !OLI_CODE PRIVATE 1187 !OLI_CODE 1188 !OLI_CODE PUBLIC lim_mp_init 1189 !OLI_CODE PUBLIC lim_mp 1190 !OLI_CODE 1191 !OLI_CODE CONTAINS 1192 !OLI_CODE 1193 !OLI_CODE SUBROUTINE lim_mp_init 1194 !OLI_CODE !!------------------------------------------------------------------- 1195 !OLI_CODE !! *** ROUTINE lim_mp_init *** 1196 !OLI_CODE !! 1197 !OLI_CODE !! ** Purpose : Initialize melt ponds 1198 !OLI_CODE !!------------------------------------------------------------------- 1199 !OLI_CODE a_ip_frac(:,:,:) = 0._wp 1200 !OLI_CODE a_ip(:,:,:) = 0._wp 1201 !OLI_CODE h_ip(:,:,:) = 0._wp 1202 !OLI_CODE v_ip(:,:,:) = 0._wp 1203 !OLI_CODE h_il(:,:,:) = 0._wp 1204 !OLI_CODE v_il(:,:,:) = 0._wp 1205 !OLI_CODE 1206 !OLI_CODE END SUBROUTINE lim_mp_init 1207 !OLI_CODE 1208 !OLI_CODE 1209 !OLI_CODE SUBROUTINE lim_mp 1210 !OLI_CODE !!------------------------------------------------------------------- 1211 !OLI_CODE !! *** ROUTINE lim_mp *** 1212 !OLI_CODE !! 1213 !OLI_CODE !! ** Purpose : Compute surface heat flux and call main melt pond 1214 !OLI_CODE !! routine 1215 !OLI_CODE !!------------------------------------------------------------------- 1216 !OLI_CODE 1217 !OLI_CODE INTEGER :: ji, jj, jl ! dummy loop indices 1218 !OLI_CODE 1219 !OLI_CODE fsurf(:,:) = 0.e0 1220 !OLI_CODE DO jl = 1, jpl 1221 !OLI_CODE DO jj = 1, jpj 1222 !OLI_CODE DO ji = 1, jpi 1223 !OLI_CODE fsurf(ji,jj) = fsurf(ji,jj) + a_i(ji,jj,jl) * & 1224 !OLI_CODE (qns_ice(ji,jj,jl) + (1.0 - izero(ji,jj,jl)) & 1225 !OLI_CODE * qsr_ice(ji,jj,jl)) 1226 !OLI_CODE END DO 1227 !OLI_CODE END DO 1228 !OLI_CODE END DO 1229 !OLI_CODE 1230 !OLI_CODE CALL compute_mp_topo(at_i, a_i, & 1231 !OLI_CODE vt_i, v_i, v_s, rhosn_glo, t_i, sz_i, a_ip_frac, & 1232 !OLI_CODE h_ip, h_il, t_su, tatm_ice, diag_sur_me*rdt_ice, & 1233 !OLI_CODE fsurf, fwoc) 1234 !OLI_CODE 1235 !OLI_CODE at_ip(:,:) = 0.0 1236 !OLI_CODE vt_ip(:,:) = 0.0 1237 !OLI_CODE vt_il(:,:) = 0.0 1238 !OLI_CODE DO jl = 1, jpl 1239 !OLI_CODE DO jj = 1, jpj 1240 !OLI_CODE DO ji = 1, jpi 1241 !OLI_CODE a_ip(ji,jj,jl) = MAX(0.0_wp, a_ip_frac(ji,jj,jl) * & 1242 !OLI_CODE a_i(ji,jj,jl)) 1243 !OLI_CODE v_ip(ji,jj,jl) = MAX(0.0_wp, a_ip_frac(ji,jj,jl) * & 1244 !OLI_CODE a_i(ji,jj,jl) * h_ip(ji,jj,jl)) 1245 !OLI_CODE v_il(ji,jj,jl) = MAX(0.0_wp, a_ip_frac(ji,jj,jl) * & 1246 !OLI_CODE a_i(ji,jj,jl) * h_il(ji,jj,jl)) 1247 !OLI_CODE at_ip(ji,jj) = at_ip(ji,jj) + a_ip(ji,jj,jl) 1248 !OLI_CODE vt_ip(ji,jj) = vt_ip(ji,jj) + v_ip(ji,jj,jl) 1249 !OLI_CODE vt_il(ji,jj) = vt_il(ji,jj) + v_il(ji,jj,jl) 1250 !OLI_CODE END DO 1251 !OLI_CODE END DO 1252 !OLI_CODE END DO 1253 !OLI_CODE 1254 !OLI_CODE END SUBROUTINE lim_mp 1255 !OLI_CODE 1256 !OLI_CODE 1257 !OLI_CODE SUBROUTINE compute_mp_topo(aice, aicen, & 1258 !OLI_CODE vice, vicen, & 1259 !OLI_CODE vsnon, rhos, & 1260 !OLI_CODE ticen, salin, & 1261 !OLI_CODE a_ip_frac, h_ip, & 1262 !OLI_CODE h_il, Tsfc, & 1263 !OLI_CODE potT, meltt, & 1264 !OLI_CODE fsurf, fwoc) 1265 !OLI_CODE !!------------------------------------------------------------------- 1266 !OLI_CODE !! *** ROUTINE compute_mp_topo *** 1267 !OLI_CODE !! 1268 !OLI_CODE !! ** Purpose : Compute melt pond evolution based on the ice 1269 !OLI_CODE !! topography as inferred from the ice thickness 1270 !OLI_CODE !! distribution. 1271 !OLI_CODE !! 1272 !OLI_CODE !! ** Method : This code is initially based on Flocco and Feltham 1273 !OLI_CODE !! (2007) and Flocco et al. (2010). More to come... 1274 !OLI_CODE !! 1275 !OLI_CODE !! ** Tunable parameters : 1276 !OLI_CODE !! 1277 !OLI_CODE !! ** Note : 1278 !OLI_CODE !! 1279 !OLI_CODE !! ** References 1280 !OLI_CODE !! Flocco, D. and D. L. Feltham, 2007. A continuum model of melt pond 1281 !OLI_CODE !! evolution on Arctic sea ice. J. Geophys. Res. 112, C08016, doi: 1282 !OLI_CODE !! 10.1029/2006JC003836. 1283 !OLI_CODE !! Flocco, D., D. L. Feltham and A. K. Turner, 2010. Incorporation of 1284 !OLI_CODE !! a physically based melt pond scheme into the sea ice component of a 1285 !OLI_CODE !! climate model. J. Geophys. Res. 115, C08012, 1286 !OLI_CODE !! doi: 10.1029/2009JC005568. 1287 !OLI_CODE !! 1288 !OLI_CODE !!------------------------------------------------------------------- 1289 !OLI_CODE 1290 !OLI_CODE REAL (wp), DIMENSION (jpi,jpj), & 1291 !OLI_CODE INTENT(IN) :: & 1292 !OLI_CODE aice, & ! total ice area fraction 1293 !OLI_CODE vice ! total ice volume (m) 1294 !OLI_CODE 1295 !OLI_CODE REAL (wp), DIMENSION (jpi,jpj,jpl), & 1296 !OLI_CODE INTENT(IN) :: & 1297 !OLI_CODE aicen, & ! ice area fraction, per category 1298 !OLI_CODE vsnon, & ! snow volume, per category (m) 1299 !OLI_CODE rhos, & ! equivalent snow density, per category (kg/m^3) 1300 !OLI_CODE vicen ! ice volume, per category (m) 1301 !OLI_CODE 1302 !OLI_CODE REAL (wp), DIMENSION (jpi,jpj,nlay_i,jpl), & 1303 !OLI_CODE INTENT(IN) :: & 1304 !OLI_CODE ticen, & ! ice enthalpy, per category 1305 !OLI_CODE salin 1306 !OLI_CODE 1307 !OLI_CODE REAL (wp), DIMENSION (jpi,jpj,jpl), & 1308 !OLI_CODE INTENT(INOUT) :: & 1309 !OLI_CODE a_ip_frac , & ! pond area fraction of ice, per ice category 1310 !OLI_CODE h_ip , & ! pond depth, per ice category (m) 1311 !OLI_CODE h_il ! Refrozen ice lid thickness, per ice category (m) 1312 !OLI_CODE 1313 !OLI_CODE REAL (wp), DIMENSION (jpi,jpj), & 1314 !OLI_CODE INTENT(IN) :: & 1315 !OLI_CODE potT, & ! air potential temperature 1316 !OLI_CODE meltt, & ! total surface meltwater flux 1317 !OLI_CODE fsurf ! thermodynamic heat flux at ice/snow surface (W/m^2) 1318 !OLI_CODE 1319 !OLI_CODE REAL (wp), DIMENSION (jpi,jpj), & 1320 !OLI_CODE INTENT(INOUT) :: & 1321 !OLI_CODE fwoc ! fresh water flux to the ocean (from draining and other pond volume adjustments) 1322 !OLI_CODE ! (m) 1323 !OLI_CODE 1324 !OLI_CODE REAL (wp), DIMENSION (jpi,jpj,jpl), & 1325 !OLI_CODE INTENT(IN) :: & 1326 !OLI_CODE Tsfc ! snow/sea ice surface temperature 1327 !OLI_CODE 1328 !OLI_CODE ! local variables 1329 !OLI_CODE REAL (wp), DIMENSION (jpi,jpj,jpl) :: & 1330 !OLI_CODE zTsfcn, & ! ice/snow surface temperature (C) 1331 !OLI_CODE zvolpn, & ! pond volume per unit area, per category (m) 1332 !OLI_CODE zvuin ! water-equivalent volume of ice lid on melt pond ('upper ice', m) 1333 !OLI_CODE 1334 !OLI_CODE REAL (wp), DIMENSION (jpi,jpj,jpl) :: & 1335 !OLI_CODE zapondn,& ! pond area fraction, per category 1336 !OLI_CODE zhpondn ! pond depth, per category (m) 1337 !OLI_CODE 1338 !OLI_CODE REAL (wp), DIMENSION (jpi,jpj) :: & 1339 !OLI_CODE zvolp ! total volume of pond, per unit area of pond (m) 1340 !OLI_CODE 1341 !OLI_CODE REAL (wp) :: & 1342 !OLI_CODE zhi, & ! ice thickness (m) 1343 !OLI_CODE zdHui, & ! change in thickness of ice lid (m) 1344 !OLI_CODE zomega, & ! conduction 1345 !OLI_CODE zdTice, & ! temperature difference across ice lid (C) 1346 !OLI_CODE zdvice, & ! change in ice volume (m) 1347 !OLI_CODE zTavg, & ! mean surface temperature across categories (C) 1348 !OLI_CODE zTp, & ! pond freezing temperature (C) 1349 !OLI_CODE zdvn ! change in melt pond volume for fresh water budget 1350 !OLI_CODE INTEGER, DIMENSION (jpi*jpj) :: & 1351 !OLI_CODE indxi, indxj ! compressed indices for cells with ice melting 1352 !OLI_CODE 1353 !OLI_CODE INTEGER :: n,k,i,j,ij,icells,indxij ! loop indices 1354 !OLI_CODE 1355 !OLI_CODE INTEGER, DIMENSION (jpl) :: & 1356 !OLI_CODE kcells ! cells where ice lid combines with vice 1357 !OLI_CODE 1358 !OLI_CODE INTEGER, DIMENSION (jpi*jpj,jpl) :: & 1359 !OLI_CODE indxii, indxjj ! i,j indices for kcells loop 1360 !OLI_CODE 1361 !OLI_CODE REAL (wp), parameter :: & 1362 !OLI_CODE zhicemin = 0.1_wp , & ! minimum ice thickness with ponds (m) 1363 !OLI_CODE zTd = 0.15_wp, & ! temperature difference for freeze-up (C) 1364 !OLI_CODE zr1_rlfus = 1._wp / 0.334e+6 / 917._wp , & ! (J/m^3) 1365 !OLI_CODE zmin_volp = 1.e-4_wp, & ! minimum pond volume (m) 1366 !OLI_CODE z0 = 0._wp, & 1367 !OLI_CODE zTimelt = 0._wp, & 1368 !OLI_CODE z01 = 0.01_wp, & 1369 !OLI_CODE z25 = 0.25_wp, & 1370 !OLI_CODE z5 = 0.5_wp, & 1371 !OLI_CODE epsi10 = 1.0e-11_wp 1372 !OLI_CODE !--------------------------------------------------------------- 1373 !OLI_CODE ! initialize 1374 !OLI_CODE !--------------------------------------------------------------- 1375 !OLI_CODE 1376 !OLI_CODE DO j = 1, jpj 1377 !OLI_CODE DO i = 1, jpi 1378 !OLI_CODE zvolp(i,j) = z0 1379 !OLI_CODE END DO 1380 !OLI_CODE END DO 1381 !OLI_CODE DO n = 1, jpl 1382 !OLI_CODE DO j = 1, jpj 1383 !OLI_CODE DO i = 1, jpi 1384 !OLI_CODE ! load tracers 1385 !OLI_CODE zvolp(i,j) = zvolp(i,j) + h_ip(i,j,n) & 1386 !OLI_CODE * a_ip_frac(i,j,n) * aicen(i,j,n) 1387 !OLI_CODE zTsfcn(i,j,n) = Tsfc(i,j,n) - rtt ! convert in Celsius - Oli 1388 !OLI_CODE zvuin (i,j,n) = h_il(i,j,n) & 1389 !OLI_CODE * a_ip_frac(i,j,n) * aicen(i,j,n) 1390 !OLI_CODE 1391 !OLI_CODE zhpondn(i,j,n) = z0 ! pond depth, per category 1392 !OLI_CODE zapondn(i,j,n) = z0 ! pond area, per category 1393 !OLI_CODE END DO 1394 !OLI_CODE END DO 1395 !OLI_CODE indxii(:,n) = 0 1396 !OLI_CODE indxjj(:,n) = 0 1397 !OLI_CODE kcells (n) = 0 1398 !OLI_CODE END DO 1399 !OLI_CODE 1400 !OLI_CODE ! The freezing temperature for meltponds is assumed slightly below 0C, 1401 !OLI_CODE ! as if meltponds had a little salt in them. The salt budget is not 1402 !OLI_CODE ! altered for meltponds, but if it were then an actual pond freezing 1403 !OLI_CODE ! temperature could be computed. 1404 !OLI_CODE 1405 !OLI_CODE zTp = zTimelt - zTd 1406 !OLI_CODE 1407 !OLI_CODE !----------------------------------------------------------------- 1408 !OLI_CODE ! Identify grid cells with ponds 1409 !OLI_CODE !----------------------------------------------------------------- 1410 !OLI_CODE 1411 !OLI_CODE icells = 0 1412 !OLI_CODE DO j = 1, jpj 1413 !OLI_CODE DO i = 1, jpi 1414 !OLI_CODE zhi = z0 1415 !OLI_CODE IF (aice(i,j) > epsi10) zhi = vice(i,j)/aice(i,j) 1416 !OLI_CODE IF ( aice(i,j) > z01 .and. zhi > zhicemin .and. & 1417 !OLI_CODE zvolp(i,j) > zmin_volp*aice(i,j)) THEN 1418 !OLI_CODE icells = icells + 1 1419 !OLI_CODE indxi(icells) = i 1420 !OLI_CODE indxj(icells) = j 1421 !OLI_CODE ELSE ! remove ponds on thin ice 1422 !OLI_CODE !fpond(i,j) = fpond(i,j) - zvolp(i,j) 1423 !OLI_CODE zvolpn(i,j,:) = z0 1424 !OLI_CODE zvuin (i,j,:) = z0 1425 !OLI_CODE zvolp (i,j) = z0 1426 !OLI_CODE END IF 1427 !OLI_CODE END DO ! i 1428 !OLI_CODE END DO ! j 1429 !OLI_CODE 1430 !OLI_CODE DO ij = 1, icells 1431 !OLI_CODE i = indxi(ij) 1432 !OLI_CODE j = indxj(ij) 1433 !OLI_CODE 1434 !OLI_CODE !-------------------------------------------------------------- 1435 !OLI_CODE ! calculate pond area and depth 1436 !OLI_CODE !-------------------------------------------------------------- 1437 !OLI_CODE CALL pond_area(aice(i,j),vice(i,j),rhos(i,j,:), & 1438 !OLI_CODE aicen(i,j,:), vicen(i,j,:), vsnon(i,j,:), & 1439 !OLI_CODE ticen(i,j,:,:), salin(i,j,:,:), & 1440 !OLI_CODE zvolpn(i,j,:), zvolp(i,j), & 1441 !OLI_CODE zapondn(i,j,:),zhpondn(i,j,:), zdvn) 1442 !OLI_CODE 1443 !OLI_CODE fwoc(i,j) = fwoc(i,j) + zdvn ! -> Goes to fresh water budget 1444 !OLI_CODE 1445 !OLI_CODE ! mean surface temperature 1446 !OLI_CODE zTavg = z0 1447 !OLI_CODE DO n = 1, jpl 1448 !OLI_CODE zTavg = zTavg + zTsfcn(i,j,n)*aicen(i,j,n) 1449 !OLI_CODE END DO 1450 !OLI_CODE zTavg = zTavg / aice(i,j) 1451 !OLI_CODE 1452 !OLI_CODE DO n = 1, jpl-1 1453 !OLI_CODE 1454 !OLI_CODE IF (zvuin(i,j,n) > epsi10) THEN 1455 !OLI_CODE 1456 !OLI_CODE !---------------------------------------------------------------- 1457 !OLI_CODE ! melting: floating upper ice layer melts in whole or part 1458 !OLI_CODE !---------------------------------------------------------------- 1459 !OLI_CODE ! IF (zTsfcn(i,j,n) > zTp) THEN 1460 !OLI_CODE IF (zTavg > zTp) THEN 1461 !OLI_CODE 1462 !OLI_CODE zdvice = min(meltt(i,j)*zapondn(i,j,n), zvuin(i,j,n)) 1463 !OLI_CODE IF (zdvice > epsi10) THEN 1464 !OLI_CODE zvuin (i,j,n) = zvuin (i,j,n) - zdvice 1465 !OLI_CODE zvolpn(i,j,n) = zvolpn(i,j,n) + zdvice 1466 !OLI_CODE zvolp (i,j) = zvolp (i,j) + zdvice 1467 !OLI_CODE !fwoc(i,j) = fwoc(i,j) + zdvice 1468 !OLI_CODE 1469 !OLI_CODE IF (zvuin(i,j,n) < epsi10 .and. zvolpn(i,j,n) > puny) THEN 1470 !OLI_CODE ! ice lid melted and category is pond covered 1471 !OLI_CODE zvolpn(i,j,n) = zvolpn(i,j,n) + zvuin(i,j,n) 1472 !OLI_CODE !fwoc(i,j) = fwoc(i,j) + zvuin(i,j,n) 1473 !OLI_CODE zvuin(i,j,n) = z0 1474 !OLI_CODE END IF 1475 !OLI_CODE zhpondn(i,j,n) = zvolpn(i,j,n) / zapondn(i,j,n) 1476 !OLI_CODE END IF 1477 !OLI_CODE 1478 !OLI_CODE !---------------------------------------------------------------- 1479 !OLI_CODE ! freezing: existing upper ice layer grows 1480 !OLI_CODE !---------------------------------------------------------------- 1481 !OLI_CODE ELSE IF (zvolpn(i,j,n) > epsi10) THEN ! zTavg <= zTp 1482 !OLI_CODE 1483 !OLI_CODE ! dIFferential growth of base of surface floating ice layer 1484 !OLI_CODE zdTice = max(-zTavg, z0) ! > 0 1485 !OLI_CODE zomega = rcdic*zdTice * zr1_rlfus 1486 !OLI_CODE zdHui = sqrt(zomega*rdt_ice + z25*(zvuin(i,j,n)/ & 1487 !OLI_CODE aicen(i,j,n))**2)- z5 * zvuin(i,j,n)/aicen(i,j,n) 1488 !OLI_CODE 1489 !OLI_CODE zdvice = min(zdHui*zapondn(i,j,n), zvolpn(i,j,n)) 1490 !OLI_CODE IF (zdvice > epsi10) THEN 1491 !OLI_CODE zvuin (i,j,n) = zvuin (i,j,n) + zdvice 1492 !OLI_CODE zvolpn(i,j,n) = zvolpn(i,j,n) - zdvice 1493 !OLI_CODE zvolp (i,j) = zvolp (i,j) - zdvice 1494 !OLI_CODE !fwoc(i,j) = fwoc(i,j) - zdvice 1495 !OLI_CODE zhpondn(i,j,n) = zvolpn(i,j,n) / zapondn(i,j,n) 1496 !OLI_CODE END IF 1497 !OLI_CODE 1498 !OLI_CODE END IF ! zTavg 1499 !OLI_CODE 1500 !OLI_CODE !---------------------------------------------------------------- 1501 !OLI_CODE ! freezing: upper ice layer begins to form 1502 !OLI_CODE ! note: albedo does not change 1503 !OLI_CODE !---------------------------------------------------------------- 1504 !OLI_CODE ELSE ! zvuin < epsi10 1505 !OLI_CODE 1506 !OLI_CODE ! thickness of newly formed ice 1507 !OLI_CODE ! the surface temperature of a meltpond is the same as that 1508 !OLI_CODE ! of the ice underneath (0C), and the thermodynamic surface 1509 !OLI_CODE ! flux is the same 1510 !OLI_CODE zdHui = max(-fsurf(i,j)*rdt_ice*zr1_rlfus, z0) 1511 !OLI_CODE zdvice = min(zdHui*zapondn(i,j,n), zvolpn(i,j,n)) 1512 !OLI_CODE IF (zdvice > epsi10) THEN 1513 !OLI_CODE zvuin (i,j,n) = zdvice 1514 !OLI_CODE zvolpn(i,j,n) = zvolpn(i,j,n) - zdvice 1515 !OLI_CODE zvolp (i,j) = zvolp (i,j) - zdvice 1516 !OLI_CODE !fwoc(i,j) = fwoc(i,j) - zdvice 1517 !OLI_CODE zhpondn(i,j,n)= zvolpn(i,j,n) / zapondn(i,j,n) 1518 !OLI_CODE END IF 1519 !OLI_CODE 1520 !OLI_CODE END IF ! zvuin 1521 !OLI_CODE 1522 !OLI_CODE END DO ! jpl 1523 !OLI_CODE 1524 !OLI_CODE END DO ! ij 1525 !OLI_CODE 1526 !OLI_CODE !--------------------------------------------------------------- 1527 !OLI_CODE ! remove ice lid if there is no liquid pond 1528 !OLI_CODE ! zvuin may be nonzero on category jpl due to dynamics 1529 !OLI_CODE !--------------------------------------------------------------- 1530 !OLI_CODE 1531 !OLI_CODE DO j = 1, jpj 1532 !OLI_CODE DO i = 1, jpi 1533 !OLI_CODE DO n = 1, jpl 1534 !OLI_CODE IF (aicen(i,j,n) > epsi10 .and. zvolpn(i,j,n) < puny & 1535 !OLI_CODE .and. zvuin (i,j,n) > epsi10) THEN 1536 !OLI_CODE kcells(n) = kcells(n) + 1 1537 !OLI_CODE indxij = kcells(n) 1538 !OLI_CODE indxii(indxij,n) = i 1539 !OLI_CODE indxjj(indxij,n) = j 1540 !OLI_CODE END IF 1541 !OLI_CODE END DO 1542 !OLI_CODE END DO ! i 1543 !OLI_CODE END DO ! j 1544 !OLI_CODE 1545 !OLI_CODE DO n = 1, jpl 1546 !OLI_CODE 1547 !OLI_CODE IF (kcells(n) > 0) THEN 1548 !OLI_CODE DO ij = 1, kcells(n) 1549 !OLI_CODE i = indxii(ij,n) 1550 !OLI_CODE j = indxjj(ij,n) 1551 !OLI_CODE fwoc(i,j) = fwoc(i,j) + rhoic/rauw * zvuin(i,j,n) ! Completely refrozen lid goes into ocean (to be changed) 1552 !OLI_CODE zvuin(i,j,n) = z0 1553 !OLI_CODE END DO ! ij 1554 !OLI_CODE END IF 1555 !OLI_CODE 1556 !OLI_CODE ! reload tracers 1557 !OLI_CODE DO j = 1, jpj 1558 !OLI_CODE DO i = 1, jpi 1559 !OLI_CODE IF (zapondn(i,j,n) > epsi10) THEN 1560 !OLI_CODE h_il(i,j,n) = zvuin(i,j,n) / zapondn(i,j,n) 1561 !OLI_CODE ELSE 1562 !OLI_CODE zvuin(i,j,n) = z0 1563 !OLI_CODE h_il(i,j,n) = z0 1564 !OLI_CODE END IF 1565 !OLI_CODE IF (aicen(i,j,n) > epsi10) THEN 1566 !OLI_CODE a_ip_frac(i,j,n) = zapondn(i,j,n) / aicen(i,j,n) * & 1567 !OLI_CODE (1.0_wp - MAX(z0, SIGN(1.0_wp, -zvolpn(i,j,n)))) 1568 !OLI_CODE h_ip(i,j,n) = zhpondn(i,j,n) 1569 !OLI_CODE ELSE 1570 !OLI_CODE a_ip_frac(i,j,n) = z0 1571 !OLI_CODE h_ip(i,j,n) = z0 1572 !OLI_CODE h_il(i,j,n) = z0 1573 !OLI_CODE END IF 1574 !OLI_CODE END DO ! i 1575 !OLI_CODE END DO ! j 1576 !OLI_CODE 1577 !OLI_CODE END DO ! n 1578 !OLI_CODE 1579 !OLI_CODE END SUBROUTINE compute_mp_topo 1580 !OLI_CODE 1581 !OLI_CODE 1582 !OLI_CODE SUBROUTINE pond_area(aice,vice,rhos, & 1583 !OLI_CODE aicen, vicen, vsnon, ticen, & 1584 !OLI_CODE salin, zvolpn, zvolp, & 1585 !OLI_CODE zapondn,zhpondn,dvolp) 1586 !OLI_CODE !!------------------------------------------------------------------- 1587 !OLI_CODE !! *** ROUTINE pond_area *** 1588 !OLI_CODE !! 1589 !OLI_CODE !! ** Purpose : Compute melt pond area, depth and melting rates 1590 !OLI_CODE !!------------------------------------------------------------------ 1591 !OLI_CODE REAL (wp), INTENT(IN) :: & 1592 !OLI_CODE aice,vice 1593 !OLI_CODE 1594 !OLI_CODE REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 1595 !OLI_CODE aicen, vicen, vsnon, rhos 1596 !OLI_CODE 1597 !OLI_CODE REAL (wp), DIMENSION(nlay_i,jpl), INTENT(IN) :: & 1598 !OLI_CODE ticen, salin 1599 !OLI_CODE 1600 !OLI_CODE REAL (wp), DIMENSION(jpl), INTENT(INOUT) :: & 1601 !OLI_CODE zvolpn 1602 !OLI_CODE 1603 !OLI_CODE REAL (wp), INTENT(INOUT) :: & 1604 !OLI_CODE zvolp, dvolp 1605 !OLI_CODE 1606 !OLI_CODE REAL (wp), DIMENSION(jpl), INTENT(OUT) :: & 1607 !OLI_CODE zapondn, zhpondn 1608 !OLI_CODE 1609 !OLI_CODE INTEGER :: & 1610 !OLI_CODE n, ns, & 1611 !OLI_CODE m_index, & 1612 !OLI_CODE permflag 1613 !OLI_CODE 1614 !OLI_CODE REAL (wp), DIMENSION(jpl) :: & 1615 !OLI_CODE hicen, & 1616 !OLI_CODE hsnon, & 1617 !OLI_CODE asnon, & 1618 !OLI_CODE alfan, & 1619 !OLI_CODE betan, & 1620 !OLI_CODE cum_max_vol, & 1621 !OLI_CODE reduced_aicen 1622 !OLI_CODE 1623 !OLI_CODE REAL (wp), DIMENSION(0:jpl) :: & 1624 !OLI_CODE cum_max_vol_tmp 1625 !OLI_CODE 1626 !OLI_CODE REAL (wp) :: & 1627 !OLI_CODE hpond, & 1628 !OLI_CODE drain, & 1629 !OLI_CODE floe_weight, & 1630 !OLI_CODE pressure_head, & 1631 !OLI_CODE hsl_rel, & 1632 !OLI_CODE deltah, & 1633 !OLI_CODE perm, & 1634 !OLI_CODE apond, & 1635 !OLI_CODE msno 1636 !OLI_CODE 1637 !OLI_CODE REAL (wp), parameter :: & 1638 !OLI_CODE viscosity = 1.79e-3_wp, & ! kinematic water viscosity in kg/m/s 1639 !OLI_CODE z0 = 0.0_wp , & 1640 !OLI_CODE c1 = 1.0_wp , & 1641 !OLI_CODE p4 = 0.4_wp , & 1642 !OLI_CODE p6 = 0.6_wp , & 1643 !OLI_CODE epsi10 = 1.0e-11_wp 1644 !OLI_CODE 1645 !OLI_CODE !-----------| 1646 !OLI_CODE ! | 1647 !OLI_CODE ! |-----------| 1648 !OLI_CODE !___________|___________|______________________________________sea-level 1649 !OLI_CODE ! | | 1650 !OLI_CODE ! | |---^--------| 1651 !OLI_CODE ! | | | | 1652 !OLI_CODE ! | | | |-----------| |------- 1653 !OLI_CODE ! | | |alfan(n)| | | 1654 !OLI_CODE ! | | | | |--------------| 1655 !OLI_CODE ! | | | | | | 1656 !OLI_CODE !---------------------------v------------------------------------------- 1657 !OLI_CODE ! | | ^ | | | 1658 !OLI_CODE ! | | | | |--------------| 1659 !OLI_CODE ! | | |betan(n)| | | 1660 !OLI_CODE ! | | | |-----------| |------- 1661 !OLI_CODE ! | | | | 1662 !OLI_CODE ! | |---v------- | 1663 !OLI_CODE ! | | 1664 !OLI_CODE ! |-----------| 1665 !OLI_CODE ! | 1666 !OLI_CODE !-----------| 1667 !OLI_CODE 1668 !OLI_CODE !------------------------------------------------------------------- 1669 !OLI_CODE ! initialize 1670 !OLI_CODE !------------------------------------------------------------------- 1671 !OLI_CODE 1672 !OLI_CODE DO n = 1, jpl 1673 !OLI_CODE 1674 !OLI_CODE zapondn(n) = z0 1675 !OLI_CODE zhpondn(n) = z0 1676 !OLI_CODE 1677 !OLI_CODE IF (aicen(n) < epsi10) THEN 1678 !OLI_CODE hicen(n) = z0 1679 !OLI_CODE hsnon(n) = z0 1680 !OLI_CODE reduced_aicen(n) = z0 1681 !OLI_CODE ELSE 1682 !OLI_CODE hicen(n) = vicen(n) / aicen(n) 1683 !OLI_CODE hsnon(n) = vsnon(n) / aicen(n) 1684 !OLI_CODE reduced_aicen(n) = c1 ! n=jpl 1685 !OLI_CODE IF (n < jpl) reduced_aicen(n) = aicen(n) & 1686 !OLI_CODE * (-0.024_wp*hicen(n) + 0.832_wp) 1687 !OLI_CODE asnon(n) = reduced_aicen(n) 1688 !OLI_CODE END IF 1689 !OLI_CODE 1690 !OLI_CODE ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 1691 !OLI_CODE ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming 1692 !OLI_CODE ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all 1693 !OLI_CODE ! categories. alfa and beta partition the ITD - they are areas not thicknesses! 1694 !OLI_CODE ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. 1695 !OLI_CODE ! Here, alfa = 60% of the ice area (and since hice is constant in a category, 1696 !OLI_CODE ! alfan = 60% of the ice volume) in each category lies above the reference line, 1697 !OLI_CODE ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 1698 !OLI_CODE 1699 !OLI_CODE alfan(n) = p6 * hicen(n) 1700 !OLI_CODE betan(n) = p4 * hicen(n) 1701 !OLI_CODE 1702 !OLI_CODE cum_max_vol(n) = z0 1703 !OLI_CODE cum_max_vol_tmp(n) = z0 1704 !OLI_CODE 1705 !OLI_CODE END DO ! jpl 1706 !OLI_CODE 1707 !OLI_CODE cum_max_vol_tmp(0) = z0 1708 !OLI_CODE drain = z0 1709 !OLI_CODE dvolp = z0 1710 !OLI_CODE 1711 !OLI_CODE !-------------------------------------------------------------------------- 1712 !OLI_CODE ! the maximum amount of water that can be contained up to each ice category 1713 !OLI_CODE !-------------------------------------------------------------------------- 1714 !OLI_CODE 1715 !OLI_CODE DO n = 1, jpl-1 ! last category can not hold any volume 1716 !OLI_CODE 1717 !OLI_CODE IF (alfan(n+1) >= alfan(n) .and. alfan(n+1) > z0) THEN 1718 !OLI_CODE 1719 !OLI_CODE ! total volume in level including snow 1720 !OLI_CODE cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) + & 1721 !OLI_CODE (alfan(n+1) - alfan(n)) * sum(reduced_aicen(1:n)) 1722 !OLI_CODE 1723 !OLI_CODE 1724 !OLI_CODE ! subtract snow solid volumes from lower categories in current level 1725 !OLI_CODE DO ns = 1, n 1726 !OLI_CODE cum_max_vol_tmp(n) = cum_max_vol_tmp(n) & 1727 !OLI_CODE - rhos(ns)/rauw * & ! fraction of snow that is occupied by solid ??rauw 1728 !OLI_CODE asnon(ns) * & ! area of snow from that category 1729 !OLI_CODE max(min(hsnon(ns)+alfan(ns)-alfan(n), alfan(n+1)- & 1730 !OLI_CODE alfan(n)), z0) 1731 !OLI_CODE ! thickness of snow from ns layer in n layer 1732 !OLI_CODE END DO 1733 !OLI_CODE 1734 !OLI_CODE ELSE ! assume higher categories unoccupied 1735 !OLI_CODE cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) 1736 !OLI_CODE END IF 1737 !OLI_CODE !IF (cum_max_vol_tmp(n) < z0) THEN 1738 !OLI_CODE ! call abort_ice('negative melt pond volume') 1739 !OLI_CODE !END IF 1740 !OLI_CODE END DO 1741 !OLI_CODE cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume 1742 !OLI_CODE cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl) 1743 !OLI_CODE 1744 !OLI_CODE !---------------------------------------------------------------- 1745 !OLI_CODE ! is there more meltwater than can be held in the floe? 1746 !OLI_CODE !---------------------------------------------------------------- 1747 !OLI_CODE IF (zvolp >= cum_max_vol(jpl)) THEN 1748 !OLI_CODE drain = zvolp - cum_max_vol(jpl) + epsi10 1749 !OLI_CODE zvolp = zvolp - drain 1750 !OLI_CODE dvolp = drain 1751 !OLI_CODE IF (zvolp < epsi10) THEN 1752 !OLI_CODE dvolp = dvolp + zvolp 1753 !OLI_CODE zvolp = z0 1754 !OLI_CODE END IF 1755 !OLI_CODE END IF 1756 !OLI_CODE 1757 !OLI_CODE ! height and area corresponding to the remaining volume 1758 !OLI_CODE 1759 !OLI_CODE call calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, & 1760 !OLI_CODE zvolp, cum_max_vol, hpond, m_index) 1761 !OLI_CODE 1762 !OLI_CODE DO n=1, m_index 1763 !OLI_CODE zhpondn(n) = hpond - alfan(n) + alfan(1) 1764 !OLI_CODE zapondn(n) = reduced_aicen(n) 1765 !OLI_CODE END DO 1766 !OLI_CODE apond = sum(zapondn(1:m_index)) 1767 !OLI_CODE 1768 !OLI_CODE !------------------------------------------------------------------------ 1769 !OLI_CODE ! drainage due to ice permeability - Darcy's law 1770 !OLI_CODE !------------------------------------------------------------------------ 1771 !OLI_CODE 1772 !OLI_CODE ! sea water level 1773 !OLI_CODE msno = z0 1774 !OLI_CODE DO n=1,jpl 1775 !OLI_CODE msno = msno + vsnon(n) * rhos(n) 1776 !OLI_CODE END DO 1777 !OLI_CODE floe_weight = (msno + rhoic*vice + rau0*zvolp) / aice 1778 !OLI_CODE hsl_rel = floe_weight / rau0 & 1779 !OLI_CODE - ((sum(betan(:)*aicen(:))/aice) + alfan(1)) 1780 !OLI_CODE 1781 !OLI_CODE deltah = hpond - hsl_rel 1782 !OLI_CODE pressure_head = grav * rau0 * max(deltah, z0) 1783 !OLI_CODE 1784 !OLI_CODE ! drain IF ice is permeable 1785 !OLI_CODE permflag = 0 1786 !OLI_CODE IF (pressure_head > z0) THEN 1787 !OLI_CODE DO n = 1, jpl-1 1788 !OLI_CODE IF (hicen(n) /= z0) THEN 1789 !OLI_CODE CALL permeability_phi(ticen(:,n), salin(:,n), vicen(n), perm) 1790 !OLI_CODE IF (perm > z0) permflag = 1 1791 !OLI_CODE drain = perm*zapondn(n)*pressure_head*rdt_ice / & 1792 !OLI_CODE (viscosity*hicen(n)) 1793 !OLI_CODE dvolp = dvolp + min(drain, zvolp) 1794 !OLI_CODE zvolp = max(zvolp - drain, z0) 1795 !OLI_CODE IF (zvolp < epsi10) THEN 1796 !OLI_CODE dvolp = dvolp + zvolp 1797 !OLI_CODE zvolp = z0 1798 !OLI_CODE END IF 1799 !OLI_CODE END IF 1800 !OLI_CODE END DO 1801 !OLI_CODE 1802 !OLI_CODE ! adjust melt pond DIMENSIONs 1803 !OLI_CODE IF (permflag > 0) THEN 1804 !OLI_CODE ! recompute pond depth 1805 !OLI_CODE CALL calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, & 1806 !OLI_CODE zvolp, cum_max_vol, hpond, m_index) 1807 !OLI_CODE DO n=1, m_index 1808 !OLI_CODE zhpondn(n) = hpond - alfan(n) + alfan(1) 1809 !OLI_CODE zapondn(n) = reduced_aicen(n) 1810 !OLI_CODE END DO 1811 !OLI_CODE apond = sum(zapondn(1:m_index)) 1812 !OLI_CODE END IF 1813 !OLI_CODE END IF ! pressure_head 1814 !OLI_CODE 1815 !OLI_CODE !------------------------------------------------------------------------ 1816 !OLI_CODE ! total melt pond volume in category DOes not include snow volume 1817 !OLI_CODE ! snow in melt ponds is not melted 1818 !OLI_CODE !------------------------------------------------------------------------ 1819 !OLI_CODE 1820 !OLI_CODE ! Calculate pond volume for lower categories 1821 !OLI_CODE DO n=1,m_index-1 1822 !OLI_CODE zvolpn(n) = zapondn(n) * zhpondn(n) & 1823 !OLI_CODE - (rhos(n)/rauw) * asnon(n) * min(hsnon(n), zhpondn(n))! 1824 !OLI_CODE END DO 1825 !OLI_CODE 1826 !OLI_CODE ! Calculate pond volume for highest category = remaining pond volume 1827 !OLI_CODE IF (m_index == 1) zvolpn(m_index) = zvolp 1828 !OLI_CODE IF (m_index > 1) THEN 1829 !OLI_CODE IF (zvolp > sum(zvolpn(1:m_index-1))) THEN 1830 !OLI_CODE zvolpn(m_index) = zvolp - sum(zvolpn(1:m_index-1)) 1831 !OLI_CODE ELSE 1832 !OLI_CODE zvolpn(m_index) = z0 1833 !OLI_CODE zhpondn(m_index) = z0 1834 !OLI_CODE zapondn(m_index) = z0 1835 !OLI_CODE ! If remaining pond volume is negative reduce pond volume of 1836 !OLI_CODE ! lower category 1837 !OLI_CODE IF (zvolp+epsi10 < sum(zvolpn(1:m_index-1))) & 1838 !OLI_CODE zvolpn(m_index-1) = zvolpn(m_index-1)-sum(zvolpn(1:m_index-1))& 1839 !OLI_CODE + zvolp 1840 !OLI_CODE END IF 1841 !OLI_CODE END IF 1842 !OLI_CODE 1843 !OLI_CODE DO n=1,m_index 1844 !OLI_CODE IF (zapondn(n) > epsi10) THEN 1845 !OLI_CODE zhpondn(n) = zvolpn(n) / zapondn(n) 1846 !OLI_CODE ELSE 1847 !OLI_CODE dvolp = dvolp + zvolpn(n) 1848 !OLI_CODE zhpondn(n) = z0 1849 !OLI_CODE zvolpn(n) = z0 1850 !OLI_CODE zapondn(n) = z0 1851 !OLI_CODE end IF 1852 !OLI_CODE END DO 1853 !OLI_CODE DO n = m_index+1, jpl 1854 !OLI_CODE zhpondn(n) = z0 1855 !OLI_CODE zapondn(n) = z0 1856 !OLI_CODE zvolpn (n) = z0 1857 !OLI_CODE END DO 1858 !OLI_CODE 1859 !OLI_CODE END SUBROUTINE pond_area 1860 !OLI_CODE 1861 !OLI_CODE 1862 !OLI_CODE SUBROUTINE calc_hpond(aicen, asnon, hsnon, rhos, alfan, & 1863 !OLI_CODE zvolp, cum_max_vol, & 1864 !OLI_CODE hpond, m_index) 1865 !OLI_CODE !!------------------------------------------------------------------- 1866 !OLI_CODE !! *** ROUTINE calc_hpond *** 1867 !OLI_CODE !! 1868 !OLI_CODE !! ** Purpose : Compute melt pond depth 1869 !OLI_CODE !!------------------------------------------------------------------- 1870 !OLI_CODE 1871 !OLI_CODE REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 1872 !OLI_CODE aicen, & 1873 !OLI_CODE asnon, & 1874 !OLI_CODE hsnon, & 1875 !OLI_CODE rhos, & 1876 !OLI_CODE alfan, & 1877 !OLI_CODE cum_max_vol 1878 !OLI_CODE 1879 !OLI_CODE REAL (wp), INTENT(IN) :: & 1880 !OLI_CODE zvolp 1881 !OLI_CODE 1882 !OLI_CODE REAL (wp), INTENT(OUT) :: & 1883 !OLI_CODE hpond 1884 !OLI_CODE 1885 !OLI_CODE INTEGER, INTENT(OUT) :: & 1886 !OLI_CODE m_index 1887 !OLI_CODE 1888 !OLI_CODE INTEGER :: n, ns 1889 !OLI_CODE 1890 !OLI_CODE REAL (wp), DIMENSION(0:jpl+1) :: & 1891 !OLI_CODE hitl, & 1892 !OLI_CODE aicetl 1893 !OLI_CODE 1894 !OLI_CODE REAL (wp) :: & 1895 !OLI_CODE rem_vol, & 1896 !OLI_CODE area, & 1897 !OLI_CODE vol, & 1898 !OLI_CODE tmp, & 1899 !OLI_CODE z0 = 0.0_wp, & 1900 !OLI_CODE epsi10 = 1.0e-11_wp 1901 !OLI_CODE 1902 !OLI_CODE !---------------------------------------------------------------- 1903 !OLI_CODE ! hpond is zero if zvolp is zero - have we fully drained? 1904 !OLI_CODE !---------------------------------------------------------------- 1905 !OLI_CODE 1906 !OLI_CODE IF (zvolp < epsi10) THEN 1907 !OLI_CODE hpond = z0 1908 !OLI_CODE m_index = 0 1909 !OLI_CODE ELSE 1910 !OLI_CODE 1911 !OLI_CODE !---------------------------------------------------------------- 1912 !OLI_CODE ! Calculate the category where water fills up to 1913 !OLI_CODE !---------------------------------------------------------------- 1914 !OLI_CODE 1915 !OLI_CODE !----------| 1916 !OLI_CODE ! | 1917 !OLI_CODE ! | 1918 !OLI_CODE ! |----------| -- -- 1919 !OLI_CODE !__________|__________|_________________________________________ ^ 1920 !OLI_CODE ! | | rem_vol ^ | Semi-filled 1921 !OLI_CODE ! | |----------|-- -- -- - ---|-- ---- -- -- --v layer 1922 !OLI_CODE ! | | | | 1923 !OLI_CODE ! | | | |hpond 1924 !OLI_CODE ! | | |----------| | |------- 1925 !OLI_CODE ! | | | | | | 1926 !OLI_CODE ! | | | |---v-----| 1927 !OLI_CODE ! | | m_index | | | 1928 !OLI_CODE !------------------------------------------------------------- 1929 !OLI_CODE 1930 !OLI_CODE m_index = 0 ! 1:m_index categories have water in them 1931 !OLI_CODE DO n = 1, jpl 1932 !OLI_CODE IF (zvolp <= cum_max_vol(n)) THEN 1933 !OLI_CODE m_index = n 1934 !OLI_CODE IF (n == 1) THEN 1935 !OLI_CODE rem_vol = zvolp 1936 !OLI_CODE ELSE 1937 !OLI_CODE rem_vol = zvolp - cum_max_vol(n-1) 1938 !OLI_CODE END IF 1939 !OLI_CODE exit ! to break out of the loop 1940 !OLI_CODE END IF 1941 !OLI_CODE END DO 1942 !OLI_CODE m_index = min(jpl-1, m_index) 1943 !OLI_CODE 1944 !OLI_CODE !---------------------------------------------------------------- 1945 !OLI_CODE ! semi-filled layer may have m_index different snow in it 1946 !OLI_CODE !---------------------------------------------------------------- 1947 !OLI_CODE 1948 !OLI_CODE !----------------------------------------------------------- ^ 1949 !OLI_CODE ! | alfan(m_index+1) 1950 !OLI_CODE ! | 1951 !OLI_CODE !hitl(3)--> |----------| | 1952 !OLI_CODE !hitl(2)--> |------------| * * * * *| | 1953 !OLI_CODE !hitl(1)--> |----------|* * * * * * |* * * * * | | 1954 !OLI_CODE !hitl(0)-->------------------------------------------------- | ^ 1955 !OLI_CODE ! various snow from lower categories | |alfa(m_index) 1956 !OLI_CODE 1957 !OLI_CODE ! hitl - heights of the snow layers from thinner and current categories 1958 !OLI_CODE ! aicetl - area of each snow depth in this layer 1959 !OLI_CODE 1960 !OLI_CODE hitl(:) = z0 1961 !OLI_CODE aicetl(:) = z0 1962 !OLI_CODE DO n = 1, m_index 1963 !OLI_CODE hitl(n) = max(min(hsnon(n) + alfan(n) - alfan(m_index), & 1964 !OLI_CODE alfan(m_index+1) - alfan(m_index)), z0) 1965 !OLI_CODE aicetl(n) = asnon(n) 1966 !OLI_CODE 1967 !OLI_CODE aicetl(0) = aicetl(0) + (aicen(n) - asnon(n)) 1968 !OLI_CODE END DO 1969 !OLI_CODE hitl(m_index+1) = alfan(m_index+1) - alfan(m_index) 1970 !OLI_CODE aicetl(m_index+1) = z0 1971 !OLI_CODE 1972 !OLI_CODE !---------------------------------------------------------------- 1973 !OLI_CODE ! reorder array according to hitl 1974 !OLI_CODE ! snow heights not necessarily in height order 1975 !OLI_CODE !---------------------------------------------------------------- 1976 !OLI_CODE 1977 !OLI_CODE DO ns = 1, m_index+1 1978 !OLI_CODE DO n = 0, m_index - ns + 1 1979 !OLI_CODE IF (hitl(n) > hitl(n+1)) THEN ! swap order 1980 !OLI_CODE tmp = hitl(n) 1981 !OLI_CODE hitl(n) = hitl(n+1) 1982 !OLI_CODE hitl(n+1) = tmp 1983 !OLI_CODE tmp = aicetl(n) 1984 !OLI_CODE aicetl(n) = aicetl(n+1) 1985 !OLI_CODE aicetl(n+1) = tmp 1986 !OLI_CODE END IF 1987 !OLI_CODE END DO 1988 !OLI_CODE END DO 1989 !OLI_CODE 1990 !OLI_CODE !---------------------------------------------------------------- 1991 !OLI_CODE ! divide semi-filled layer into set of sublayers each vertically homogenous 1992 !OLI_CODE !---------------------------------------------------------------- 1993 !OLI_CODE 1994 !OLI_CODE !hitl(3)---------------------------------------------------------------- 1995 !OLI_CODE ! | * * * * * * * * 1996 !OLI_CODE ! |* * * * * * * * * 1997 !OLI_CODE !hitl(2)---------------------------------------------------------------- 1998 !OLI_CODE ! | * * * * * * * * | * * * * * * * * 1999 !OLI_CODE ! |* * * * * * * * * |* * * * * * * * * 2000 !OLI_CODE !hitl(1)---------------------------------------------------------------- 2001 !OLI_CODE ! | * * * * * * * * | * * * * * * * * | * * * * * * * * 2002 !OLI_CODE ! |* * * * * * * * * |* * * * * * * * * |* * * * * * * * * 2003 !OLI_CODE !hitl(0)---------------------------------------------------------------- 2004 !OLI_CODE ! aicetl(0) aicetl(1) aicetl(2) aicetl(3) 2005 !OLI_CODE 2006 !OLI_CODE ! move up over layers incrementing volume 2007 !OLI_CODE DO n = 1, m_index+1 2008 !OLI_CODE 2009 !OLI_CODE area = sum(aicetl(:)) - & ! total area of sub-layer 2010 !OLI_CODE (rhos(n)/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow 2011 !OLI_CODE 2012 !OLI_CODE vol = (hitl(n) - hitl(n-1)) * area ! thickness of sub-layer times area 2013 !OLI_CODE 2014 !OLI_CODE IF (vol >= rem_vol) THEN ! have reached the sub-layer with the depth within 2015 !OLI_CODE hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - & 2016 !OLI_CODE alfan(1) 2017 !OLI_CODE exit 2018 !OLI_CODE ELSE ! still in sub-layer below the sub-layer with the depth 2019 !OLI_CODE rem_vol = rem_vol - vol 2020 !OLI_CODE END IF 2021 !OLI_CODE 2022 !OLI_CODE END DO 2023 !OLI_CODE 2024 !OLI_CODE END IF 2025 !OLI_CODE 2026 !OLI_CODE END SUBROUTINE calc_hpond 2027 !OLI_CODE 2028 !OLI_CODE 2029 !OLI_CODE SUBROUTINE permeability_phi(ticen, salin, vicen, perm) 2030 !OLI_CODE !!------------------------------------------------------------------- 2031 !OLI_CODE !! *** ROUTINE permeability_phi *** 2032 !OLI_CODE !! 2033 !OLI_CODE !! ** Purpose : Determine the liquid fraction of brine in the ice 2034 !OLI_CODE !! and its permeability 2035 !OLI_CODE !!------------------------------------------------------------------- 2036 !OLI_CODE REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: & 2037 !OLI_CODE ticen, & ! energy of melting for each ice layer (J/m2) 2038 !OLI_CODE salin 2039 !OLI_CODE 2040 !OLI_CODE REAL (wp), INTENT(IN) :: & 2041 !OLI_CODE vicen ! ice volume 2042 !OLI_CODE 2043 !OLI_CODE REAL (wp), INTENT(OUT) :: & 2044 !OLI_CODE perm ! permeability 2045 !OLI_CODE 2046 !OLI_CODE REAL (wp) :: & 2047 !OLI_CODE Sbr ! brine salinity 2048 !OLI_CODE 2049 !OLI_CODE REAL (wp), DIMENSION(nlay_i) :: & 2050 !OLI_CODE Tin, & ! ice temperature 2051 !OLI_CODE phi ! liquid fraction 2052 !OLI_CODE 2053 !OLI_CODE INTEGER :: k 2054 !OLI_CODE 2055 !OLI_CODE REAL (wp) :: & 2056 !OLI_CODE c2 = 2.0_wp 2057 !OLI_CODE 2058 !OLI_CODE !----------------------------------------------------------------- 2059 !OLI_CODE ! Compute ice temperatures from enthalpies using quadratic formula 2060 !OLI_CODE !----------------------------------------------------------------- 2061 !OLI_CODE 2062 !OLI_CODE DO k = 1,nlay_i 2063 !OLI_CODE Tin(k) = ticen(k) 2064 !OLI_CODE END DO 2065 !OLI_CODE 2066 !OLI_CODE !----------------------------------------------------------------- 2067 !OLI_CODE ! brine salinity and liquid fraction 2068 !OLI_CODE !----------------------------------------------------------------- 2069 !OLI_CODE 2070 !OLI_CODE IF (maxval(Tin-rtt) <= -c2) THEN 2071 !OLI_CODE 2072 !OLI_CODE DO k = 1,nlay_i 2073 !OLI_CODE Sbr = - 1.2_wp & 2074 !OLI_CODE -21.8_wp * (Tin(k)-rtt) & 2075 !OLI_CODE - 0.919_wp * (Tin(k)-rtt)**2 & 2076 !OLI_CODE - 0.01878_wp * (Tin(k)-rtt)**3 2077 !OLI_CODE phi(k) = salin(k)/Sbr ! liquid fraction 2078 !OLI_CODE END DO ! k 2079 !OLI_CODE 2080 !OLI_CODE ELSE 2081 !OLI_CODE 2082 !OLI_CODE DO k = 1,nlay_i 2083 !OLI_CODE Sbr = -17.6_wp * (Tin(k)-rtt) & 2084 !OLI_CODE - 0.389_wp * (Tin(k)-rtt)**2 & 2085 !OLI_CODE - 0.00362_wp* (Tin(k)-rtt)**3 2086 !OLI_CODE phi(k) = salin(k)/Sbr ! liquid fraction 2087 !OLI_CODE END DO 2088 !OLI_CODE 2089 !OLI_CODE END IF 2090 !OLI_CODE 2091 !OLI_CODE !----------------------------------------------------------------- 2092 !OLI_CODE ! permeability 2093 !OLI_CODE !----------------------------------------------------------------- 2094 !OLI_CODE 2095 !OLI_CODE perm = 3.0e-08_wp * (minval(phi))**3 2096 !OLI_CODE 2097 !OLI_CODE END SUBROUTINE permeability_phi 2098 !OLI_CODE 2099 !OLI_CODE #else 2100 !OLI_CODE !!---------------------------------------------------------------------- 2101 !OLI_CODE !! Default option Dummy Module No LIM-3 sea-ice model 2102 !OLI_CODE !!---------------------------------------------------------------------- 2103 !OLI_CODE CONTAINS 2104 !OLI_CODE SUBROUTINE lim_mp_init ! Empty routine 2105 !OLI_CODE END SUBROUTINE lim_mp_init 2106 !OLI_CODE SUBROUTINE lim_mp ! Empty routine 2107 !OLI_CODE END SUBROUTINE lim_mp 2108 !OLI_CODE SUBROUTINE compute_mp_topo ! Empty routine 2109 !OLI_CODE END SUBROUTINE compute_mp_topo 2110 !OLI_CODE SUBROUTINE pond_area ! Empty routine 2111 !OLI_CODE END SUBROUTINE pond_area 2112 !OLI_CODE SUBROUTINE calc_hpond ! Empty routine 2113 !OLI_CODE END SUBROUTINE calc_hpond 2114 !OLI_CODE SUBROUTINE permeability_phy ! Empty routine 2115 !OLI_CODE END SUBROUTINE permeability_phy 2116 !OLI_CODE #endif 2117 !OLI_CODE !!====================================================================== 2118 !OLI_CODE END MODULE limmp_topo 2119 !OLI_CODE 263 END SUBROUTINE lim_mp_H12 264 265 SUBROUTINE lim_mp_init 266 !!------------------------------------------------------------------- 267 !! *** ROUTINE lim_mp_init *** 268 !! 269 !! ** Purpose : Physical constants and parameters linked to melt ponds 270 !! over sea ice 271 !! 272 !! ** Method : Read the nammp namelist and check the melt pond 273 !! parameter values called at the first timestep (nit000) 274 !! 275 !! ** input : Namelist nammp 276 !!------------------------------------------------------------------- 277 INTEGER :: ios, ioptio ! Local integer output status for namelist read 278 NAMELIST/nammp/ ln_pnd_H12, ln_pnd_fwb, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb 279 !!------------------------------------------------------------------- 280 281 REWIND( numnam_ice_ref ) ! Namelist nammp in reference namelist : Melt Ponds 282 READ ( numnam_ice_ref, nammp, IOSTAT = ios, ERR = 901) 283 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammp in reference namelist', lwp ) 284 285 REWIND( numnam_ice_cfg ) ! Namelist nammp in configuration namelist : Melt Ponds 286 READ ( numnam_ice_cfg, nammp, IOSTAT = ios, ERR = 902 ) 287 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammp in configuration namelist', lwp ) 288 IF(lwm) WRITE ( numoni, nammp ) 289 290 IF(lwp) THEN ! control print 291 WRITE(numout,*) 292 WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds' 293 WRITE(numout,*) '~~~~~~~~~~~~~~~~' 294 WRITE(numout,*) ' Namelist namicethd_pnd:' 295 WRITE(numout,*) ' Evolutive melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12 296 WRITE(numout,*) ' Melt ponds store fresh water or not ln_pnd_fwb = ', ln_pnd_fwb 297 WRITE(numout,*) ' Prescribed melt pond fraction and depth ln_pnd_Cst = ', ln_pnd_CST 298 WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd 299 WRITE(numout,*) ' Prescribed pond depth rn_hpnd = ', rn_hpnd 300 WRITE(numout,*) ' Melt ponds affect albedo or not ln_pnd_alb = ', ln_pnd_alb 301 ENDIF 302 ! 303 ! !== set the choice of ice pond scheme ==! 304 ioptio = 0 305 nice_pnd = np_pndNO 306 IF( ln_pnd_CST ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndCST ; ENDIF 307 IF( ln_pnd_H12 ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndH12 ; ENDIF 308 IF( ioptio > 1 ) CALL ctl_stop( 'ice_thd_pnd_init: choose one and only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' ) 309 310 SELECT CASE( nice_pnd ) 311 CASE( np_pndNO ) 312 IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when no ponds' ) ; ENDIF 313 IF(ln_pnd_alb) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF 314 CASE( np_pndCST) 315 IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when ln_pnd_CST=true' ) ; ENDIF 316 END SELECT 317 ! 318 END SUBROUTINE lim_mp_init 319 2120 320 #else 2121 321 !!---------------------------------------------------------------------- 2122 !! Default option Empty module NO LIM sea-ice model 2123 !!---------------------------------------------------------------------- 2124 CONTAINS 2125 SUBROUTINE lim_mp_init ! Empty routine 2126 END SUBROUTINE lim_mp_init 2127 SUBROUTINE lim_mp ! Empty routine 2128 END SUBROUTINE lim_mp 2129 SUBROUTINE lim_mp_topo ! Empty routine 2130 END SUBROUTINE lim_mp_topo 2131 SUBROUTINE lim_mp_cesm ! Empty routine 2132 END SUBROUTINE lim_mp_cesm 2133 SUBROUTINE lim_mp_area ! Empty routine 2134 END SUBROUTINE lim_mp_area 2135 SUBROUTINE lim_mp_perm ! Empty routine 2136 END SUBROUTINE lim_mp_perm 322 !! Default option Empty module NO ESIM sea-ice model 323 !!---------------------------------------------------------------------- 2137 324 #endif 2138 325
Note: See TracChangeset
for help on using the changeset viewer.