- Timestamp:
- Jun 27, 2025, 7:30:06 PM (5 weeks ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90
r5717 r5728 97 97 cldfra_in, qvc_in, qliq_in, qice_in, shear, pbl_eps, cell_area, & 98 98 temp, qtot_in, qsat, gamma_cond, ratqs, keepgoing, pt_pron_clds, & 99 dzsed_abv, flsed_abv, cfsed_abv _in, &99 dzsed_abv, flsed_abv, cfsed_abv, & 100 100 dzsed_lincont_abv, flsed_lincont_abv, cfsed_lincont_abv, & 101 101 dzsed_circont_abv, flsed_circont_abv, cfsed_circont_abv, & … … 110 110 dcfl_ini, dqil_ini, dqtl_ini, dcfl_sub, dqil_sub, dqtl_sub, & 111 111 dcfl_cir, dqtl_cir, dcfl_mix, dqil_mix, dqtl_mix, & 112 dcfl_sed, dqtl_sed, dqil_sed, dcfc_sed, dqtc_sed, dqic_sed, & 112 113 dcfc_sub, dqic_sub, dqtc_sub, dcfc_mix, dqic_mix, dqtc_mix) 113 114 … … 172 173 REAL, INTENT(IN) , DIMENSION(klon) :: dzsed_abv ! sedimentated cloud height IN THE LAYER ABOVE [m] 173 174 REAL, INTENT(IN) , DIMENSION(klon) :: flsed_abv ! sedimentated ice flux FROM THE LAYER ABOVE [kg/s/m2] 174 REAL, INTENT(IN) , DIMENSION(klon) :: cfsed_abv _in! cloud fraction IN THE LAYER ABOVE [-]175 REAL, INTENT(IN) , DIMENSION(klon) :: cfsed_abv ! cloud fraction IN THE LAYER ABOVE [-] 175 176 REAL, INTENT(IN) , DIMENSION(klon) :: dzsed_lincont_abv ! sedimentated linear contrails height IN THE LAYER ABOVE [m] 176 177 REAL, INTENT(IN) , DIMENSION(klon) :: flsed_lincont_abv ! sedimentated ice flux in linear contrails FROM THE LAYER ABOVE [kg/s/m2] … … 251 252 REAL, INTENT(INOUT), DIMENSION(klon) :: dqil_mix ! linear contrails ice specific humidity tendency because of mixing [kg/kg/s] 252 253 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtl_mix ! linear contrails total specific humidity tendency because of mixing [kg/kg/s] 254 REAL, INTENT(INOUT), DIMENSION(klon) :: dcfl_sed ! linear contrails cloud fraction tendency because of ice sedimentation [s-1] 255 REAL, INTENT(INOUT), DIMENSION(klon) :: dqil_sed ! linear contrails ice specific humidity tendency because of ice sedimentation [kg/kg/s] 256 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtl_sed ! linear contrails total specific humidity tendency because of ice sedimentation [kg/kg/s] 253 257 REAL, INTENT(INOUT), DIMENSION(klon) :: dcfc_sub ! contrail cirrus cloud fraction tendency because of sublimation [s-1] 254 258 REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_sub ! contrail cirrus ice specific humidity tendency because of sublimation [kg/kg/s] … … 257 261 REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_mix ! contrail cirrus ice specific humidity tendency because of mixing [kg/kg/s] 258 262 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtc_mix ! contrail cirrus total specific humidity tendency because of mixing [kg/kg/s] 263 REAL, INTENT(INOUT), DIMENSION(klon) :: dcfc_sed ! contrail cirrus cloud fraction tendency because of ice sedimentation [s-1] 264 REAL, INTENT(INOUT), DIMENSION(klon) :: dqic_sed ! contrail cirrus ice specific humidity tendency because of ice sedimentation [kg/kg/s] 265 REAL, INTENT(INOUT), DIMENSION(klon) :: dqtc_sed ! contrail cirrus total specific humidity tendency because of ice sedimentation [kg/kg/s] 259 266 ! 260 267 ! Local … … 289 296 ! 290 297 ! for sedimentation 291 REAL, DIMENSION(klon) :: cfsed_abv,qised_abv298 REAL, DIMENSION(klon) :: qised_abv 292 299 REAL :: iwc, icefall_velo, qice_sedim 293 REAL :: sedfra_abv, sedfra1, sedfra2, sedfra3 300 REAL :: sedfra_abv, sedfra1, sedfra2, sedfra3, sedfra_tmp 301 REAL :: dcf_sed1, dcf_sed2, dqvc_sed1, dqvc_sed2, dqi_sed1, dqi_sed2 294 302 ! 295 303 ! for mixing … … 309 317 REAL :: contrails_conversion_factor 310 318 REAL, DIMENSION(klon) :: qised_lincont_abv, qised_circont_abv 319 REAL :: dcfl_sed1, dcfl_sed2, dqtl_sed1, dqtl_sed2, dqil_sed1, dqil_sed2 320 REAL :: dcfc_sed1, dcfc_sed2, dqtc_sed1, dqtc_sed2, dqic_sed1, dqic_sed2 311 321 312 322 qzero(:) = 0. … … 401 411 402 412 qised_abv(i) = 0. 403 cfsed_abv(i) = cfsed_abv_in(i)404 413 IF ( dzsed_abv(i) .GT. eps ) THEN 405 414 !--If ice sedimentation is activated, the quantity of sedimentated ice was added … … 514 523 ENDIF 515 524 525 qised_lincont_abv(i) = 0. 516 526 IF ( dzsed_lincont_abv(i) .GT. eps ) THEN 527 !--If ice sedimentation is activated, the quantity of sedimentated ice was added 528 !--to the total water vapor in the precipitation routine. Here we remove it 529 !--(it will be reincluded later) 517 530 qised_lincont_abv(i) = flsed_lincont_abv(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime 518 qised_abv(i) = qised_abv(i) - qised_lincont_abv(i) 519 cfsed_abv(i) = cfsed_abv(i) - cfsed_lincont_abv(i) 531 qclr(i) = qclr(i) - qised_lincont_abv(i) 520 532 ENDIF 521 533 534 qised_circont_abv(i) = 0. 522 535 IF ( dzsed_circont_abv(i) .GT. eps ) THEN 536 !--If ice sedimentation is activated, the quantity of sedimentated ice was added 537 !--to the total water vapor in the precipitation routine. Here we remove it 538 !--(it will be reincluded later) 523 539 qised_circont_abv(i) = flsed_circont_abv(i) / ( paprsdn(i) - paprsup(i) ) * RG * dtime 524 qised_abv(i) = qised_abv(i) - qised_circont_abv(i) 525 cfsed_abv(i) = cfsed_abv(i) - cfsed_circont_abv(i) 540 qclr(i) = qclr(i) - qised_circont_abv(i) 526 541 ENDIF 527 542 … … 537 552 dqil_mix(i) = 0. 538 553 dqtl_mix(i) = 0. 554 dcfl_sed(i) = 0. 555 dqil_sed(i) = 0. 556 dqtl_sed(i) = 0. 539 557 dcfc_sub(i) = 0. 540 558 dqic_sub(i) = 0. … … 543 561 dqic_mix(i) = 0. 544 562 dqtc_mix(i) = 0. 563 dcfc_sed(i) = 0. 564 dqic_sed(i) = 0. 565 dqtc_sed(i) = 0. 545 566 ELSE 546 567 lincontfra(i) = 0. … … 572 593 cldfra(i) = MAX(0., cldfra(i) - lincontfra(i)) 573 594 qcld(i) = MAX(0., qcld(i) - qlincont(i)) 574 qvc(i) = MAX(0., qvc(i) - qsat(i) * lincontfra(i))595 qvc(i) = MAX(0., MIN(qcld(i), qvc(i) - qsat(i) * lincontfra(i))) 575 596 ENDIF ! lincontfra(i) .GT. eps 576 597 … … 592 613 cldfra(i) = MAX(0., cldfra(i) - circontfra(i)) 593 614 qcld(i) = MAX(0., qcld(i) - qcircont(i)) 594 qvc(i) = MAX(0., qvc(i) - qsat(i) * circontfra(i))615 qvc(i) = MAX(0., MIN(qcld(i), qvc(i) - qsat(i) * circontfra(i))) 595 616 ENDIF ! circontfra(i) .GT. eps 596 617 … … 1377 1398 1378 1399 IF ( ok_ice_sedim ) THEN 1400 !--Initialisation 1401 dcf_sed1 = 0. 1402 dqvc_sed1 = 0. 1403 dqi_sed1 = 0. 1404 dcf_sed2 = 0. 1405 dqvc_sed2 = 0. 1406 dqi_sed2 = 0. 1407 dcfl_sed1 = 0. 1408 dqtl_sed1 = 0. 1409 dqil_sed1 = 0. 1410 dcfl_sed2 = 0. 1411 dqtl_sed2 = 0. 1412 dqil_sed2 = 0. 1413 dcfc_sed1 = 0. 1414 dqtc_sed1 = 0. 1415 dqic_sed1 = 0. 1416 dcfc_sed2 = 0. 1417 dqtc_sed2 = 0. 1418 dqic_sed2 = 0. 1379 1419 !--------------------------------------- 1380 1420 !-- ICE SEDIMENTATION -- … … 1393 1433 qice_sedim = ( qcld(i) - qvc(i) ) * dzsed(i) / dz 1394 1434 !--Tendencies 1395 dcf_sed (i)= - cldfra(i) * dzsed(i) / dz1396 dqi_sed (i)= - qice_sedim1397 dqvc_sed (i)= - qvc(i) * dzsed(i) / dz1435 dcf_sed1 = - cldfra(i) * dzsed(i) / dz 1436 dqi_sed1 = - qice_sedim 1437 dqvc_sed1 = - qvc(i) * dzsed(i) / dz 1398 1438 !--Convert to flux 1399 1439 flsed(i) = qice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime 1440 !--Save the vapor in the cloud (will be needed later) 1441 qvapincld = qvc(i) / cldfra(i) 1442 !--Add tendencies 1443 cldfra(i) = cldfra(i) + dcf_sed1 1444 qcld(i) = qcld(i) + dqvc_sed1 + dqi_sed1 1445 qvc(i) = qvc(i) + dqvc_sed1 1446 clrfra(i) = clrfra(i) - dcf_sed1 1447 qclr(i) = qclr(i) - dqvc_sed1 - dqi_sed1 1448 ENDIF 1449 ! 1450 IF ( lincontfra(i) .GT. eps ) THEN 1451 icefall_velo = fallice_linear_contrails 1452 dzsed_lincont(i) = icefall_velo * dtime * MAX(0., MIN(1., (2. - icefall_velo * dtime / dz))) 1453 cfsed_lincont(i) = lincontfra(i) 1454 qice_sedim = ( qlincont(i) - qsat(i) ) * dzsed_lincont(i) / dz 1455 !--Tendencies 1456 dcfl_sed1 = - lincontfra(i) * dzsed_lincont(i) / dz 1457 dqil_sed1 = - qice_sedim 1458 dqtl_sed1 = - qlincont(i) * dzsed_lincont(i) / dz 1459 !--Convert to flux 1460 flsed_lincont(i) = qice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime 1461 !--Add tendencies 1462 lincontfra(i) = lincontfra(i) + dcfl_sed1 1463 qlincont(i) = qlincont(i) + dqtl_sed1 1464 clrfra(i) = clrfra(i) - dcfl_sed1 1465 qclr(i) = qclr(i) - dqtl_sed1 1466 ENDIF 1467 ! 1468 IF ( circontfra(i) .GT. eps ) THEN 1469 icefall_velo = fallice_cirrus_contrails 1470 dzsed_circont(i) = icefall_velo * dtime * MAX(0., MIN(1., (2. - icefall_velo * dtime / dz))) 1471 cfsed_circont(i) = circontfra(i) 1472 qice_sedim = ( qcircont(i) - qsat(i) ) * dzsed_circont(i) / dz 1473 !--Tendencies 1474 dcfc_sed1 = - circontfra(i) * dzsed_circont(i) / dz 1475 dqic_sed1 = - qice_sedim 1476 dqtc_sed1 = - qcircont(i) * dzsed_circont(i) / dz 1477 !--Convert to flux 1478 flsed_circont(i) = qice_sedim * ( paprsdn(i) - paprsup(i) ) / RG / dtime 1479 !--Add tendencies 1480 circontfra(i) = circontfra(i) + dcfc_sed1 1481 qcircont(i) = qcircont(i) + dqtc_sed1 1482 clrfra(i) = clrfra(i) - dcfc_sed1 1483 qclr(i) = qclr(i) - dqtc_sed1 1400 1484 ENDIF 1401 1485 ! … … 1417 1501 !--ice sedimentation 1418 1502 !--(1) we use the clear-sky PDF to determine the humidity and increase the 1419 !--cloud fraction / sublimate falling ice 1503 !--cloud fraction / sublimate falling ice. NB it can also fall into 1504 !--another cloud class 1420 1505 !--(2) we do not increase cloud fraction and add the falling ice to the cloud 1421 1506 !--(3) we increase the cloud fraction but use in-cloud water vapor as … … 1430 1515 1431 1516 !--For region (1) 1432 IF ( ( sedfra1 .GT. eps ) .AND. ( clrfra(i) .GT. eps ) ) THEN 1433 1434 IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN 1435 tauinv_depsub = qiceincld**(1./3.) * kappa_depsub 1436 qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime * tauinv_depsub ) ) 1437 ELSE 1438 qvapinclr_lim = qsat(i) - qiceincld 1517 IF ( sedfra1 .GT. eps ) THEN 1518 IF ( clrfra(i) .GT. eps ) THEN 1519 IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN 1520 tauinv_depsub = qiceincld**(1./3.) * kappa_depsub 1521 qvapinclr_lim = qsat(i) - qiceincld / ( 1. - EXP( - dtime * tauinv_depsub ) ) 1522 ELSE 1523 qvapinclr_lim = qsat(i) - qiceincld 1524 ENDIF 1525 1526 rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. 1527 pdf_x = qvapinclr_lim / qsatl(i) * 100. 1528 pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) 1529 pdf_x = qsat(i) / qsatl(i) * 100. 1530 pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i) 1531 IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows 1532 pdf_fra_above_lim = 0. 1533 pdf_q_above_lim = 0. 1534 ELSEIF ( pdf_y .LT. -10. ) THEN 1535 pdf_fra_above_lim = clrfra(i) 1536 pdf_q_above_lim = qclr(i) 1537 ELSE 1538 pdf_y = EXP( pdf_y ) 1539 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 1540 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 1541 pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) 1542 pdf_q_above_lim = ( pdf_e3 * clrfra(i) & 1543 + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. 1544 ENDIF 1545 1546 pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim 1547 1548 !--The first three lines allow to favor ISSR rather than subsaturated sky for 1549 !--the growth. The fourth line indicates that the ice is falling only in 1550 !--the clear-sky area. The rest falls into other cloud classes 1551 sedfra_tmp = MAX(sedfra1 - pdf_fra_below_lim, MIN(pdf_fra_above_lim, & 1552 sedfra1 * chi_sedim * pdf_fra_above_lim & 1553 / ( pdf_fra_below_lim + chi_sedim * pdf_fra_above_lim ))) & 1554 * clrfra(i) / (totfra_in(i) - cldfra(i)) 1555 1556 IF ( pdf_fra_above_lim .GT. eps ) THEN 1557 IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN 1558 dcf_sed2 = dcf_sed2 + sedfra_tmp 1559 dqvc_sed2 = dqvc_sed2 + sedfra_tmp * pdf_q_above_lim / pdf_fra_above_lim 1560 dqi_sed2 = dqi_sed2 + sedfra_tmp * qiceincld 1561 ELSE 1562 dcf_sed2 = dcf_sed2 + sedfra_tmp 1563 dqvc_sed2 = dqvc_sed2 + sedfra_tmp * qsat(i) 1564 dqi_sed2 = dqi_sed2 + sedfra_tmp & 1565 * (qiceincld + pdf_q_above_lim / pdf_fra_above_lim - qsat(i)) 1566 ENDIF 1567 ENDIF 1568 ENDIF ! clrfra(i) .GT. eps 1569 1570 !--If the sedimented ice falls into a linear contrail, it increases its content 1571 IF ( lincontfra(i) .GT. eps ) THEN 1572 sedfra_tmp = sedfra1 * lincontfra(i) / (totfra_in(i) - cldfra(i)) 1573 dqil_sed2 = dqil_sed2 + sedfra_tmp * qiceincld 1439 1574 ENDIF 1440 1575 1441 rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. 1442 pdf_x = qvapinclr_lim / qsatl(i) * 100. 1443 pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) 1444 pdf_x = qsat(i) / qsatl(i) * 100. 1445 pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i) 1446 IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows 1447 pdf_fra_above_lim = 0. 1448 pdf_q_above_lim = 0. 1449 ELSEIF ( pdf_y .LT. -10. ) THEN 1450 pdf_fra_above_lim = clrfra(i) 1451 pdf_q_above_lim = qclr(i) 1452 ELSE 1453 pdf_y = EXP( pdf_y ) 1454 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 1455 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 1456 pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) 1457 pdf_q_above_lim = ( pdf_e3 * clrfra(i) & 1458 + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. 1576 !--If the sedimented ice falls into a contrail cirrus, it increases its content 1577 IF ( circontfra(i) .GT. eps ) THEN 1578 sedfra_tmp = sedfra1 * circontfra(i) / (totfra_in(i) - cldfra(i)) 1579 dqic_sed2 = dqic_sed2 + sedfra_tmp * qiceincld 1459 1580 ENDIF 1460 1461 sedfra1 = MAX(sedfra1 - pdf_fra_below_lim, MIN(pdf_fra_above_lim, & 1462 sedfra1 * chi_sedim * pdf_fra_above_lim & 1463 / ( pdf_fra_below_lim + chi_sedim * pdf_fra_above_lim ))) 1464 1465 IF ( pdf_fra_above_lim .GT. eps ) THEN 1466 IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN 1467 dcf_sed(i) = dcf_sed(i) + sedfra1 1468 dqvc_sed(i) = dqvc_sed(i) + sedfra1 * pdf_q_above_lim / pdf_fra_above_lim 1469 dqi_sed(i) = dqi_sed(i) + sedfra1 * qiceincld 1470 ELSE 1471 dcf_sed(i) = dcf_sed(i) + sedfra1 1472 dqvc_sed(i) = dqvc_sed(i) + sedfra1 * qsat(i) 1473 dqi_sed(i) = dqi_sed(i) + sedfra1 & 1474 * (qiceincld + pdf_q_above_lim / pdf_fra_above_lim - qsat(i)) 1475 ENDIF 1476 ENDIF 1477 ENDIF ! ( sedfra1 .GT. eps .AND. ( clrfra(i) .GT. eps ) 1581 ENDIF ! sedfra1 .GT. eps 1478 1582 1479 1583 !--For region (2) 1480 dqi_sed (i) = dqi_sed(i)+ sedfra2 * qiceincld1584 dqi_sed2 = dqi_sed2 + sedfra2 * qiceincld 1481 1585 1482 1586 !--For region (3) 1483 1587 IF ( cldfra(i) .GT. eps ) THEN 1484 dcf_sed (i) = dcf_sed(i)+ sedfra31485 dqvc_sed (i) = dqvc_sed(i) + sedfra3 * qvc(i) / cldfra(i)1486 dqi_sed (i) = dqi_sed(i)+ sedfra3 * qiceincld1588 dcf_sed2 = dcf_sed2 + sedfra3 1589 dqvc_sed2 = dqvc_sed2 + sedfra3 * qvapincld 1590 dqi_sed2 = dqi_sed2 + sedfra3 * qiceincld 1487 1591 ENDIF 1488 ENDIF ! qised_abv(i) .GT. 0. 1592 ENDIF ! sedfra_abv .GT. eps 1593 1594 IF ( ok_plane_contrail ) THEN 1595 sedfra_abv = MIN(dzsed_lincont_abv(i), dz) / dz * cfsed_lincont_abv(i) 1596 IF ( sedfra_abv .GT. eps ) THEN 1597 1598 !--Three regions to be defined: (1) clear-sky, (2) cloudy, and 1599 !--(3) region where it was previously cloudy but now clear-sky because of 1600 !--ice sedimentation 1601 !--(1) we use the clear-sky PDF to determine the humidity and increase the 1602 !--cloud fraction / sublimate falling ice. NB it can also fall into 1603 !--another cloud class 1604 !--(2) we do not increase cloud fraction and add the falling ice to the cloud 1605 !--(3) we increase the cloud fraction but use in-cloud water vapor as 1606 !--background vapor 1607 sedfra2 = MIN(cfsed_lincont(i), cfsed_lincont_abv(i)) & 1608 * MAX(0., MIN(dz, dzsed_lincont_abv(i)) - dzsed_lincont(i)) / dz 1609 sedfra3 = MIN(cfsed_lincont(i), cfsed_lincont_abv(i)) & 1610 * MIN(MIN(dz, dzsed_lincont_abv(i)), dzsed_lincont(i)) / dz 1611 sedfra1 = sedfra_abv - sedfra3 - sedfra2 1612 1613 qiceincld = qised_lincont_abv(i) / sedfra_abv 1614 1615 !--For region (1) 1616 IF ( sedfra1 .GT. eps ) THEN 1617 IF ( clrfra(i) .GT. eps ) THEN 1618 qvapinclr_lim = qsat(i) - qiceincld 1619 1620 rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. 1621 pdf_x = qvapinclr_lim / qsatl(i) * 100. 1622 pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) 1623 pdf_x = qsat(i) / qsatl(i) * 100. 1624 pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i) 1625 IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows 1626 pdf_fra_above_lim = 0. 1627 pdf_q_above_lim = 0. 1628 ELSEIF ( pdf_y .LT. -10. ) THEN 1629 pdf_fra_above_lim = clrfra(i) 1630 pdf_q_above_lim = qclr(i) 1631 ELSE 1632 pdf_y = EXP( pdf_y ) 1633 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 1634 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 1635 pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) 1636 pdf_q_above_lim = ( pdf_e3 * clrfra(i) & 1637 + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. 1638 ENDIF 1639 1640 pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim 1641 1642 !--The first three lines allow to favor ISSR rather than subsaturated sky for 1643 !--the growth. The fourth line indicates that the ice is falling only in 1644 !--the clear-sky area. The rest falls into other cloud classes 1645 sedfra_tmp = MAX(sedfra1 - pdf_fra_below_lim, MIN(pdf_fra_above_lim, & 1646 sedfra1 * chi_sedim * pdf_fra_above_lim & 1647 / ( pdf_fra_below_lim + chi_sedim * pdf_fra_above_lim ))) & 1648 * clrfra(i) / (totfra_in(i) - lincontfra(i)) 1649 1650 IF ( pdf_fra_above_lim .GT. eps ) THEN 1651 dcfl_sed2 = dcfl_sed2 + sedfra_tmp 1652 dqtl_sed2 = dqtl_sed2 + sedfra_tmp * qsat(i) 1653 dqil_sed2 = dqil_sed2 + sedfra_tmp & 1654 * (qiceincld + pdf_q_above_lim / pdf_fra_above_lim - qsat(i)) 1655 ENDIF 1656 ENDIF ! clrfra(i) .GT. eps 1657 1658 !--If the sedimented ice falls into a natural cirrus, it increases its content 1659 IF ( cldfra(i) .GT. eps ) THEN 1660 sedfra_tmp = sedfra1 * cldfra(i) / (totfra_in(i) - lincontfra(i)) 1661 dqi_sed2 = dqi_sed2 + sedfra_tmp * qiceincld 1662 ENDIF 1663 1664 !--If the sedimented ice falls into a contrail cirrus, it increases its content 1665 IF ( circontfra(i) .GT. eps ) THEN 1666 sedfra_tmp = sedfra1 * circontfra(i) / (totfra_in(i) - lincontfra(i)) 1667 dqic_sed2 = dqic_sed2 + sedfra_tmp * qiceincld 1668 ENDIF 1669 ENDIF ! sedfra1 .GT. eps 1670 1671 !--For region (2) 1672 dqil_sed2 = dqil_sed2 + sedfra2 * qiceincld 1673 1674 !--For region (3) 1675 IF ( lincontfra(i) .GT. eps ) THEN 1676 dcfl_sed2 = dcfl_sed2 + sedfra3 1677 dqtl_sed2 = dqtl_sed2 + sedfra3 * (qsat(i) + qiceincld) 1678 dqil_sed2 = dqil_sed2 + sedfra3 * qiceincld 1679 ENDIF 1680 ENDIF ! sedfra_abv .GT. eps 1681 1682 ! 1683 sedfra_abv = MIN(dzsed_circont_abv(i), dz) / dz * cfsed_circont_abv(i) 1684 IF ( sedfra_abv .GT. eps ) THEN 1685 1686 !--Three regions to be defined: (1) clear-sky, (2) cloudy, and 1687 !--(3) region where it was previously cloudy but now clear-sky because of 1688 !--ice sedimentation 1689 !--(1) we use the clear-sky PDF to determine the humidity and increase the 1690 !--cloud fraction / sublimate falling ice. NB it can also fall into 1691 !--another cloud class 1692 !--(2) we do not increase cloud fraction and add the falling ice to the cloud 1693 !--(3) we increase the cloud fraction but use in-cloud water vapor as 1694 !--background vapor 1695 sedfra2 = MIN(cfsed_circont(i), cfsed_circont_abv(i)) & 1696 * MAX(0., MIN(dz, dzsed_circont_abv(i)) - dzsed_circont(i)) / dz 1697 sedfra3 = MIN(cfsed_circont(i), cfsed_circont_abv(i)) & 1698 * MIN(MIN(dz, dzsed_circont_abv(i)), dzsed_circont(i)) / dz 1699 sedfra1 = sedfra_abv - sedfra3 - sedfra2 1700 1701 qiceincld = qised_circont_abv(i) / sedfra_abv 1702 1703 !--For region (1) 1704 IF ( sedfra1 .GT. eps ) THEN 1705 IF ( clrfra(i) .GT. eps ) THEN 1706 qvapinclr_lim = qsat(i) - qiceincld 1707 1708 rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. 1709 pdf_x = qvapinclr_lim / qsatl(i) * 100. 1710 pdf_loc = rhl_clr - pdf_scale(i) * pdf_gamma(i) 1711 pdf_x = qsat(i) / qsatl(i) * 100. 1712 pdf_y = LOG( MAX( ( pdf_x - pdf_loc ) / pdf_scale(i), eps) ) * pdf_alpha(i) 1713 IF ( pdf_y .GT. 10. ) THEN !--Avoid overflows 1714 pdf_fra_above_lim = 0. 1715 pdf_q_above_lim = 0. 1716 ELSEIF ( pdf_y .LT. -10. ) THEN 1717 pdf_fra_above_lim = clrfra(i) 1718 pdf_q_above_lim = qclr(i) 1719 ELSE 1720 pdf_y = EXP( pdf_y ) 1721 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 1722 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_gamma(i) 1723 pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) 1724 pdf_q_above_lim = ( pdf_e3 * clrfra(i) & 1725 + pdf_loc * pdf_fra_above_lim ) * qsatl(i) / 100. 1726 ENDIF 1727 1728 pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim 1729 1730 !--The first three lines allow to favor ISSR rather than subsaturated sky for 1731 !--the growth. The fourth line indicates that the ice is falling only in 1732 !--the clear-sky area. The rest falls into other cloud classes 1733 sedfra_tmp = MAX(sedfra1 - pdf_fra_below_lim, MIN(pdf_fra_above_lim, & 1734 sedfra1 * chi_sedim * pdf_fra_above_lim & 1735 / ( pdf_fra_below_lim + chi_sedim * pdf_fra_above_lim ))) & 1736 * clrfra(i) / (totfra_in(i) - circontfra(i)) 1737 1738 IF ( pdf_fra_above_lim .GT. eps ) THEN 1739 dcfc_sed2 = dcfc_sed2 + sedfra_tmp 1740 dqtc_sed2 = dqtc_sed2 + sedfra_tmp * qsat(i) 1741 dqic_sed2 = dqic_sed2 + sedfra_tmp & 1742 * (qiceincld + pdf_q_above_lim / pdf_fra_above_lim - qsat(i)) 1743 ENDIF 1744 ENDIF ! clrfra(i) .GT. eps 1745 1746 !--If the sedimented ice falls into a natural cirrus, it increases its content 1747 IF ( cldfra(i) .GT. eps ) THEN 1748 sedfra_tmp = sedfra1 * cldfra(i) / (totfra_in(i) - circontfra(i)) 1749 dqi_sed2 = dqi_sed2 + sedfra_tmp * qiceincld 1750 ENDIF 1751 1752 !--If the sedimented ice falls into a contrail cirrus, it increases its content 1753 IF ( lincontfra(i) .GT. eps ) THEN 1754 sedfra_tmp = sedfra1 * lincontfra(i) / (totfra_in(i) - circontfra(i)) 1755 dqil_sed2 = dqil_sed2 + sedfra_tmp * qiceincld 1756 ENDIF 1757 ENDIF ! sedfra1 .GT. eps 1758 1759 !--For region (2) 1760 dqic_sed2 = dqic_sed2 + sedfra2 * qiceincld 1761 1762 !--For region (3) 1763 IF ( circontfra(i) .GT. eps ) THEN 1764 dcfc_sed2 = dcfc_sed2 + sedfra3 1765 dqtc_sed2 = dqtc_sed2 + sedfra3 * (qsat(i) + qiceincld) 1766 dqic_sed2 = dqic_sed2 + sedfra3 * qiceincld 1767 ENDIF 1768 ENDIF ! sedfra_abv .GT. eps 1769 ENDIF ! ok_plane_contrails 1489 1770 1490 1771 !--Add tendencies 1491 cldfra(i) = MIN(totfra_in(i), cldfra(i) + dcf_sed(i)) 1492 qcld(i) = qcld(i) + dqvc_sed(i) + dqi_sed(i) 1493 qvc(i) = qvc(i) + dqvc_sed(i) 1494 clrfra(i) = MAX(0., clrfra(i) - dcf_sed(i)) 1772 cldfra(i) = cldfra(i) + dcf_sed2 1773 qcld(i) = qcld(i) + dqvc_sed2 + dqi_sed2 1774 qvc(i) = qvc(i) + dqvc_sed2 1775 clrfra(i) = clrfra(i) - dcf_sed2 1776 qclr(i) = qclr(i) - dqvc_sed2 - dqi_sed2 1495 1777 !--We re-include sublimated sedimentated ice into clear sky water vapor 1496 qclr(i) = qclr(i) - dqvc_sed(i) + qised_abv(i) - dqi_sed(i) 1778 qclr(i) = qclr(i) + qised_abv(i) 1779 !--Diagnostics 1780 dcf_sed(i) = dcf_sed1 + dcf_sed2 1781 dqvc_sed(i) = dqvc_sed1 + dqvc_sed2 1782 dqi_sed(i) = dqi_sed1 + dqi_sed2 1783 1784 if ( ok_plane_contrail ) THEN 1785 !--Add tendencies 1786 lincontfra(i) = lincontfra(i) + dcfl_sed2 1787 qlincont(i) = qlincont(i) + dqtl_sed2 1788 clrfra(i) = clrfra(i) - dcfl_sed2 1789 qclr(i) = qclr(i) - dqtl_sed2 1790 !--We re-include sublimated sedimentated ice into clear sky water vapor 1791 qclr(i) = qclr(i) + qised_lincont_abv(i) 1792 !--Diagnostics 1793 dcfl_sed(i) = dcfl_sed1 + dcfl_sed2 1794 dqtl_sed(i) = dqtl_sed1 + dqtl_sed2 1795 dqil_sed(i) = dqil_sed1 + dqil_sed2 1796 1797 !--Add tendencies 1798 circontfra(i) = circontfra(i) + dcfc_sed2 1799 qcircont(i) = qcircont(i) + dqtc_sed2 1800 clrfra(i) = clrfra(i) - dcfc_sed2 1801 qclr(i) = qclr(i) - dqtc_sed2 1802 !--We re-include sublimated sedimentated ice into clear sky water vapor 1803 qclr(i) = qclr(i) + qised_circont_abv(i) 1804 !--Diagnostics 1805 dcfc_sed(i) = dcfc_sed1 + dcfc_sed2 1806 dqtc_sed(i) = dqtc_sed1 + dqtc_sed2 1807 dqic_sed(i) = dqic_sed1 + dqic_sed2 1808 ENDIF 1497 1809 1498 1810 ENDIF ! ok_ice_sedim
Note: See TracChangeset
for help on using the changeset viewer.
