Changeset 5602
- Timestamp:
- Apr 3, 2025, 4:53:58 PM (2 months ago)
- Location:
- LMDZ6/branches/contrails/libf/phylmd
- Files:
-
- 5 edited
-
lmdz_aviation.f90 (modified) (7 diffs)
-
lmdz_lscp.f90 (modified) (5 diffs)
-
lmdz_lscp_condensation.f90 (modified) (25 diffs)
-
lmdz_lscp_ini.f90 (modified) (5 diffs)
-
physiq_mod.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/contrails/libf/phylmd/lmdz_aviation.f90
r5601 r5602 73 73 SUBROUTINE contrails_formation( & 74 74 klon, dtime, pplay, temp, qsat, qsatl, gamma_cond, flight_dist, & 75 cl dfra, pdf_loc, pdf_scale, pdf_alpha, keepgoing, &75 clrfra, pdf_loc, pdf_scale, pdf_alpha, keepgoing, & 76 76 Tcritcont, qcritcont, potcontfraP, potcontfraNP, & 77 77 dcfa_ini, dqia_ini, dqta_ini) … … 96 96 REAL, INTENT(IN) , DIMENSION(klon) :: gamma_cond ! condensation threshold w.r.t. qsat [-] 97 97 REAL, INTENT(IN) , DIMENSION(klon) :: flight_dist ! aviation distance flown concentration [m/s/m3] 98 REAL, INTENT(IN) , DIMENSION(klon) :: cl dfra ! cloudfraction [-]98 REAL, INTENT(IN) , DIMENSION(klon) :: clrfra ! clear sky fraction [-] 99 99 REAL, INTENT(IN) , DIMENSION(klon) :: pdf_loc ! location parameter of the clear sky PDF [%] 100 100 REAL, INTENT(IN) , DIMENSION(klon) :: pdf_scale ! scale parameter of the clear sky PDF [%] … … 166 166 167 167 168 IF ( ( ( 1. - cldfra(i)) .GT. eps ) .AND. ( temp(i) .LT. Tcritcont(i) ) ) THEN168 IF ( ( clrfra(i) .GT. eps ) .AND. ( temp(i) .LT. Tcritcont(i) ) ) THEN 169 169 170 170 pdf_x = qcritcont(i) / qsatl(i) * 100. … … 173 173 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 174 174 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2 175 pdf_fra_above_qcritcont = EXP( - pdf_y ) * ( 1. - cldfra(i))176 pdf_q_above_qcritcont = ( pdf_e3 * ( 1. - cldfra(i)) &175 pdf_fra_above_qcritcont = EXP( - pdf_y ) * clrfra(i) 176 pdf_q_above_qcritcont = ( pdf_e3 * clrfra(i) & 177 177 + pdf_loc(i) * pdf_fra_above_qcritcont ) * qsatl(i) / 100. 178 178 … … 182 182 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 183 183 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2 184 pdf_fra_above_qnuc = EXP( - pdf_y ) * ( 1. - cldfra(i))185 pdf_q_above_qnuc = ( pdf_e3 * ( 1. - cldfra(i)) &184 pdf_fra_above_qnuc = EXP( - pdf_y ) * clrfra(i) 185 pdf_q_above_qnuc = ( pdf_e3 * clrfra(i) & 186 186 + pdf_loc(i) * pdf_fra_above_qnuc ) * qsatl(i) / 100. 187 187 … … 191 191 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 192 192 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2 193 pdf_fra_above_qsat = EXP( - pdf_y ) * ( 1. - cldfra(i))194 pdf_q_above_qsat = ( pdf_e3 * ( 1. - cldfra(i)) &193 pdf_fra_above_qsat = EXP( - pdf_y ) * clrfra(i) 194 pdf_q_above_qsat = ( pdf_e3 * clrfra(i) & 195 195 + pdf_loc(i) * pdf_fra_above_qsat ) * qsatl(i) / 100. 196 196 … … 330 330 qvapincld_depsub_contrails = qsat + ( qvapincld - qsat ) & 331 331 * EXP( - dtime * qiceincld / kappa / rm_ice**2 ) 332 IF ( qvapincld_depsub_contrails .GE. ( qvapincld + qiceincld ) ) & 333 qvapincld_depsub_contrails = 0. 332 334 ENDIF ! qvapincld .GT. qsat 333 335 -
LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp.f90
r5601 r5602 8 8 SUBROUTINE lscp(klon,klev,dtime,missing_val, & 9 9 paprs, pplay, omega, temp, qt, ql_seri, qi_seri, & 10 ptconv, ratqs, sigma_qtherm,&10 ptconv, cldfracv, qcondcv, ratqs, sigma_qtherm, & 11 11 d_t, d_q, d_ql, d_qi, rneb, rneblsvol, & 12 12 pfraclr, pfracld, & … … 153 153 ! CR: if iflag_ice_thermo=2, only convection is active 154 154 LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active 155 REAL, DIMENSION(klon,klev), INTENT(IN) :: cldfracv ! cloud fraction from deep convection [-] 156 REAL, DIMENSION(klon,klev), INTENT(IN) :: qcondcv ! in-cloud condensed specific humidity from deep convection [kg/kg] 155 157 156 158 !Inputs associated with thermal plumes … … 741 743 742 744 CALL condensation_ice_supersat( & 743 klon, dtime, missing_val, &744 pplay(:,k), paprs(:,k), paprs(:,k+1), &745 klon, dtime, pplay(:,k), paprs(:,k), paprs(:,k+1), & 746 cldfracv(:,k), qcondcv(:,k), & 745 747 cf_seri(:,k), qvc_seri(:,k), qliq_in, qice_in, & 746 748 shear, tke_dissip(:,k), cell_area, stratomask(:,k), & … … 949 951 950 952 IF (ok_plane_contrail) THEN 951 !--Contrails do not precipitate. We remove then from the variables temporarily 953 !!--Contrails do not precipitate. We remove then from the variables temporarily 954 !DO i = 1, klon 955 ! rneb(i,k) = rneb(i,k) - contfra(i) 956 ! zoliqi(i) = zoliqi(i) - ( qcont(i) - zqs(i) * contfra(i) ) 957 !ENDDO 958 !--Contrails precipitate as natural clouds. We save the partition of ice 959 !--between natural clouds and contrails 960 !--NB. we use qcont as a temporary variable to save this partition 952 961 DO i = 1, klon 953 rneb(i,k) = rneb(i,k) - contfra(i) 954 zoliqi(i) = zoliqi(i) - ( qcont(i) - zqs(i) * contfra(i) ) 962 IF ( zoliqi(i) .GT. 0. ) THEN 963 qcont(i) = ( qcont(i) - zqs(i) * contfra(i) ) / zoliqi(i) 964 ELSE 965 qcont(i) = 0. 966 ENDIF 955 967 ENDDO 956 968 ENDIF … … 990 1002 991 1003 IF (ok_plane_contrail) THEN 992 !--Contrails are reintroduced in the variables 1004 !!--Contrails are reintroduced in the variables 1005 !DO i = 1, klon 1006 ! rneb(i,k) = rneb(i,k) + contfra(i) 1007 ! zoliqi(i) = zoliqi(i) + ( qcont(i) - zqs(i) * contfra(i) ) 1008 !ENDDO 1009 !--Contrails fraction is left unchanged, but contrails water has changed 993 1010 DO i = 1, klon 994 rneb(i,k) = rneb(i,k) + contfra(i) 995 zoliqi(i) = zoliqi(i) + ( qcont(i) - zqs(i) * contfra(i) ) 1011 IF ( zoliqi(i) .LE. 0. ) THEN 1012 contfra(i) = 0. 1013 qcont(i) = 0. 1014 ELSE 1015 qcont(i) = zqs(i) * contfra(i) + zoliqi(i) * qcont(i) 1016 ENDIF 996 1017 ENDDO 997 1018 ENDIF -
LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_condensation.f90
r5601 r5602 94 94 !********************************************************************************** 95 95 SUBROUTINE condensation_ice_supersat( & 96 klon, dtime, missing_val, pplay, paprsdn, paprsup, &96 klon, dtime, pplay, paprsdn, paprsup, cldfracv, qcondcv, & 97 97 cldfra_in, qvc_in, qliq_in, qice_in, shear, pbl_eps, cell_area, stratomask, & 98 temp, qtot , qsat, gamma_cond, ratqs, keepgoing, &98 temp, qtot_in, qsat, gamma_cond, ratqs, keepgoing, & 99 99 cldfra, qincld, qvc, issrfra, qissr, dcf_sub, dcf_con, dcf_mix, & 100 100 dqi_adj, dqi_sub, dqi_con, dqi_mix, dqvc_adj, dqvc_sub, dqvc_con, dqvc_mix, & … … 125 125 USE lmdz_lscp_ini, ONLY: beta_pdf_lscp, temp_thresh_pdf_lscp 126 126 USE lmdz_lscp_ini, ONLY: std100_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp 127 USE lmdz_lscp_ini, ONLY: coef_mixing_lscp, coef_shear_lscp , chi_mixing_lscp127 USE lmdz_lscp_ini, ONLY: coef_mixing_lscp, coef_shear_lscp 128 128 129 129 USE lmdz_lscp_ini, ONLY: ok_plane_contrail, aspect_ratio_contrails 130 130 USE lmdz_lscp_ini, ONLY: coef_mixing_contrails, coef_shear_contrails 131 USE lmdz_lscp_ini, ONLY: linear_contrails_lifetime131 USE lmdz_lscp_ini, ONLY: chi_mixing_contrails, linear_contrails_lifetime 132 132 USE lmdz_aviation, ONLY: contrails_formation 133 133 … … 139 139 INTEGER, INTENT(IN) :: klon ! number of horizontal grid points 140 140 REAL, INTENT(IN) :: dtime ! time step [s] 141 REAL, INTENT(IN) :: missing_val ! missing value for output142 141 ! 143 142 REAL, INTENT(IN) , DIMENSION(klon) :: pplay ! layer pressure [Pa] 144 143 REAL, INTENT(IN) , DIMENSION(klon) :: paprsdn ! pressure at the lower interface [Pa] 145 144 REAL, INTENT(IN) , DIMENSION(klon) :: paprsup ! pressure at the upper interface [Pa] 145 REAL, INTENT(IN) , DIMENSION(klon) :: cldfracv ! cloud fraction from deep convection [-] 146 REAL, INTENT(IN) , DIMENSION(klon) :: qcondcv ! in-cloud condensed specific humidity from deep convection [kg/kg] 146 147 REAL, INTENT(IN) , DIMENSION(klon) :: cldfra_in ! cloud fraction [-] 147 148 REAL, INTENT(IN) , DIMENSION(klon) :: qvc_in ! gridbox-mean water vapor in cloud [kg/kg] … … 153 154 REAL, INTENT(IN) , DIMENSION(klon) :: stratomask ! fraction of stratosphere in the mesh (1 or 0) 154 155 REAL, INTENT(IN) , DIMENSION(klon) :: temp ! temperature [K] 155 REAL, INTENT(IN) , DIMENSION(klon) :: qtot 156 REAL, INTENT(IN) , DIMENSION(klon) :: qtot_in ! total specific humidity (without precip) [kg/kg] 156 157 REAL, INTENT(IN) , DIMENSION(klon) :: qsat ! saturation specific humidity [kg/kg] 157 158 REAL, INTENT(IN) , DIMENSION(klon) :: gamma_cond ! condensation threshold w.r.t. qsat [-] … … 219 220 INTEGER :: i 220 221 LOGICAL :: ok_warm_cloud 221 REAL, DIMENSION(klon) :: q cld, qzero222 REAL, DIMENSION(klon) :: qtot, qcld, qzero 222 223 REAL, DIMENSION(klon) :: clrfra, qclr 223 REAL, DIMENSION(klon) :: clrfra_pdf, qclr_pdf224 224 ! 225 225 ! for lognormal … … 247 247 REAL :: L_shear, shear_fra 248 248 REAL :: qvapinclr_lim 249 REAL :: pdf_fra_above_nuc, pdf_q_above_nuc 249 250 REAL :: pdf_fra_above_lim, pdf_q_above_lim 250 REAL :: pdf_fra_below_lim , pdf_q_below_lim251 REAL :: pdf_fra_below_lim 251 252 ! 252 253 ! for cell properties … … 254 255 255 256 qzero(:) = 0. 257 qtot = qtot_in 256 258 257 259 !--Calculation of qsat w.r.t. liquid … … 317 319 ok_warm_cloud = ( temp(i) .GT. temp_nowater ) 318 320 321 !--We remove the convection water from the total available water 322 qtot(i) = qtot(i) - ( qcondcv(i) + qsat(i) ) * cldfracv(i) 323 319 324 !--The following barriers ensure that the traced cloud properties 320 325 !--are consistent. In some rare cases, i.e. the cloud water vapor 321 326 !--can be greater than the total water in the gridbox 322 cldfra(i) = MAX(0., MIN(1. , cldfra_in(i)))327 cldfra(i) = MAX(0., MIN(1. - cldfracv(i), cldfra_in(i))) 323 328 qcld(i) = MAX(0., MIN(qtot(i), qliq_in(i) + qice_in(i) + qvc_in(i))) 324 329 qvc(i) = MAX(0., MIN(qcld(i), qvc_in(i))) 325 330 326 !-- We initclear fraction properties327 clrfra(i) = 1.- cldfra(i)331 !--Initialise clear fraction properties 332 clrfra(i) = ( 1. - cldfracv(i) ) - cldfra(i) 328 333 qclr(i) = qtot(i) - qcld(i) 329 334 … … 384 389 contfra(i) = 0. 385 390 qcont(i) = 0. 391 clrfra(i) = clrfra(i) + dcfa_sub(i) 392 qclr(i) = qclr(i) + dqta_sub(i) 386 393 ELSE 387 394 !--We remove contrails from the main class … … 414 421 qcld(i) = 0. 415 422 qvc(i) = 0. 423 clrfra(i) = clrfra(i) - dcf_sub(i) 424 qclr(i) = qclr(i) - dqvc_sub(i) - dqi_sub(i) 416 425 417 426 !--If all the ice has been sublimated, we sublimate … … 443 452 qcld(i) = 0. 444 453 qvc(i) = 0. 454 clrfra(i) = clrfra(i) - dcf_sub(i) 455 qclr(i) = qclr(i) - dqvc_sub(i) - dqi_sub(i) 445 456 ENDIF 446 457 ELSE … … 482 493 483 494 !--Add tendencies 484 cldfra(i) = MAX(0., cldfra(i) + dcf_sub(i)) 485 qcld(i) = MAX(0., qcld(i) + dqt_sub) 486 qvc(i) = MAX(0., qvc(i) + dqvc_sub(i)) 487 495 cldfra(i) = cldfra(i) + dcf_sub(i) 496 qcld(i) = qcld(i) + dqt_sub 497 qvc(i) = qvc(i) + dqvc_sub(i) 498 clrfra(i) = clrfra(i) - dcf_sub(i) 499 qclr(i) = qclr(i) - dqt_sub 488 500 ENDIF ! qvapincld_new .GT. 0. 489 501 … … 497 509 !--This section relies on a distribution of water in the clear-sky region of 498 510 !--the mesh. 499 clrfra_pdf(i) = 1. - cldfra(i) - contfra(i)500 qclr_pdf(i) = qtot(i) - qcld(i) - qcont(i)501 511 502 512 !--If there is a clear-sky region 503 IF ( clrfra _pdf(i) .GT. eps ) THEN513 IF ( clrfra(i) .GT. eps ) THEN 504 514 505 515 !--New PDF 506 rhl_clr = qclr_pdf(i) / clrfra_pdf(i) / qsatl(i) * 100. 516 rhl_clr = qclr(i) / clrfra(i) / qsatl(i) * 100. 517 rhl_clr = MAX(0., MIN(1000., rhl_clr)) 507 518 508 519 !--Calculation of the properties of the PDF … … 516 527 ! * MAX( temp(i) - temp_thresh_pdf_lscp, 0. ) 517 528 !-- tuning coef * (cell area**0.25) * (function of temperature) 518 pdf_e1 = beta_pdf_lscp * ( clrfra _pdf(i) * cell_area(i) )**0.25 &529 pdf_e1 = beta_pdf_lscp * ( clrfra(i) * cell_area(i) )**0.25 & 519 530 * MAX( temp(i) - temp_thresh_pdf_lscp, 0. ) 520 531 IF ( rhl_clr .GT. 50. ) THEN … … 524 535 ENDIF 525 536 pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. ) 526 pdf_alpha(i) = EXP( MIN(1000., rhl_clr)/ 100. ) * pdf_e3537 pdf_alpha(i) = EXP( rhl_clr / 100. ) * pdf_e3 527 538 528 539 IF ( ok_warm_cloud ) THEN … … 550 561 IF ( cf_cond .GT. eps ) THEN 551 562 552 dcf_con(i) = clrfra_pdf(i) * cf_cond 553 dqt_con = clrfra_pdf(i) * qt_cond 563 dcf_con(i) = clrfra(i) * cf_cond 564 dqt_con = clrfra(i) * qt_cond 565 566 !--Barriers (should be useless 567 dcf_con(i) = MIN(dcf_con(i), clrfra(i)) 568 dqt_con = MIN(dqt_con, qclr(i)) 554 569 555 570 IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN … … 573 588 574 589 !--Add tendencies 575 cldfra(i) = MIN(1., cldfra(i) + dcf_con(i))576 qcld(i) = MIN(qtot(i), qcld(i) + dqt_con)577 qvc(i) = MIN(qcld(i), qvc(i) + dqvc_con(i))578 clrfra(i) = MAX(0., clrfra(i) - dcf_con(i))579 qclr(i) = MAX(0., qclr(i) - dqt_con)590 cldfra(i) = cldfra(i) + dcf_con(i) 591 qcld(i) = qcld(i) + dqt_con 592 qvc(i) = qvc(i) + dqvc_con(i) 593 clrfra(i) = clrfra(i) - dcf_con(i) 594 qclr(i) = qclr(i) - dqt_con 580 595 581 596 ENDIF ! ok_warm_cloud, cf_cond .GT. eps 582 583 !--We then calculate the part that is greater than qsat584 !--and lower than gamma_cond * qsat, which is the ice supersaturated region585 pdf_x = qsat(i) / qsatl(i) * 100.586 pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i)587 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y )588 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2589 issrfra(i) = EXP( - pdf_y ) * clrfra_pdf(i)590 qissr(i) = ( pdf_e3 * clrfra_pdf(i) + pdf_loc(i) * issrfra(i) ) * qsatl(i) / 100.591 592 issrfra(i) = issrfra(i) - dcf_con(i)593 qissr(i) = qissr(i) - dqvc_con(i) - dqi_con(i)594 595 597 ENDIF ! ( 1. - cldfra(i) ) .GT. eps 596 598 … … 702 704 !--because the clear sky fraction can only be reduced by condensation. 703 705 !--Thus the `pdf_xxx` variables are well defined. 706 pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100. 707 pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i) 708 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 709 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2 710 pdf_fra_above_nuc = EXP( - pdf_y ) * clrfra(i) 711 pdf_q_above_nuc = ( pdf_e3 * clrfra(i) & 712 + pdf_loc(i) * pdf_fra_above_nuc ) * qsatl(i) / 100. 713 704 714 pdf_x = qvapinclr_lim / qsatl(i) * 100. 705 715 pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i) 706 716 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 707 717 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2 708 pdf_fra_above_lim = EXP( - pdf_y ) * clrfra _pdf(i)709 pdf_q_above_lim = ( pdf_e3 * clrfra _pdf(i) &718 pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) 719 pdf_q_above_lim = ( pdf_e3 * clrfra(i) & 710 720 + pdf_loc(i) * pdf_fra_above_lim ) * qsatl(i) / 100. 711 721 712 pdf_fra_below_lim = clrfra_pdf(i) - pdf_fra_above_lim 713 pdf_q_below_lim = qclr_pdf(i) - pdf_q_above_lim 714 !--We remove the already condensated part (it is well defined) 715 pdf_fra_above_lim = pdf_fra_above_lim - dcf_con(i) 716 pdf_q_above_lim = pdf_q_above_lim - dqvc_con(i) - dqi_con(i) 717 718 !--sigma_mix quantifies 722 pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim 723 pdf_fra_above_lim = pdf_fra_above_lim - pdf_fra_above_nuc 724 pdf_q_above_lim = pdf_q_above_lim - pdf_q_above_nuc 725 726 !--sigma_mix is the ratio of the surroundings of the clouds where mixing 727 !--increases the size of the cloud, to the total surroundings of the clouds. 728 !--This implies that ( 1. - sigma_mix ) quantifies the ratio where mixing 729 !--decreases the size of the clouds 719 730 sigma_mix = pdf_fra_above_lim / ( pdf_fra_below_lim + pdf_fra_above_lim ) 720 731 … … 735 746 736 747 !--Add tendencies 737 cldfra(i) = MAX(0., MIN(1., cldfra(i) + dcf_mix(i)))738 qcld(i) = MAX(0., MIN(qtot(i), qcld(i) + dqvc_mix(i) + dqi_mix(i)))739 qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_mix(i)))740 clrfra(i) = MAX(0., MIN(1., clrfra(i) - dcf_mix(i)))741 qclr(i) = MAX(0., MIN(qtot(i), qclr(i) - dqvc_mix(i) - dqi_mix(i)))748 cldfra(i) = cldfra(i) + dcf_mix(i) 749 qcld(i) = qcld(i) + dqvc_mix(i) + dqi_mix(i) 750 qvc(i) = qvc(i) + dqvc_mix(i) 751 clrfra(i) = clrfra(i) - dcf_mix(i) 752 qclr(i) = qclr(i) - dqvc_mix(i) - dqi_mix(i) 742 753 743 754 ENDIF ! ( cldfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps ) … … 850 861 !--because the clear sky fraction can only be reduced by condensation. 851 862 !--Thus the `pdf_xxx` variables are well defined. 863 pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100. 864 pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i) 865 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 866 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2 867 pdf_fra_above_nuc = EXP( - pdf_y ) * clrfra(i) 868 pdf_q_above_nuc = ( pdf_e3 * clrfra(i) & 869 + pdf_loc(i) * pdf_fra_above_nuc ) * qsatl(i) / 100. 870 852 871 pdf_x = qvapinclr_lim / qsatl(i) * 100. 853 872 pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i) 854 873 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 855 874 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2 856 pdf_fra_above_lim = EXP( - pdf_y ) * clrfra _pdf(i)857 pdf_q_above_lim = ( pdf_e3 * clrfra _pdf(i) &875 pdf_fra_above_lim = EXP( - pdf_y ) * clrfra(i) 876 pdf_q_above_lim = ( pdf_e3 * clrfra(i) & 858 877 + pdf_loc(i) * pdf_fra_above_lim ) * qsatl(i) / 100. 859 878 860 pdf_fra_below_lim = clrfra_pdf(i) - pdf_fra_above_lim 861 pdf_q_below_lim = qclr_pdf(i) - pdf_q_above_lim 862 !--We remove the already condensated part (it is well defined) 863 pdf_fra_above_lim = pdf_fra_above_lim - dcf_con(i) 864 pdf_q_above_lim = pdf_q_above_lim - dqvc_con(i) - dqi_con(i) 865 866 !--sigma_mix quantifies 867 sigma_mix = pdf_fra_above_lim / ( pdf_fra_below_lim + pdf_fra_above_lim ) 879 pdf_fra_below_lim = clrfra(i) - pdf_fra_above_lim 880 pdf_fra_above_lim = pdf_fra_above_lim - pdf_fra_above_nuc 881 pdf_q_above_lim = pdf_q_above_lim - pdf_q_above_nuc 882 883 !--sigma_mix is the ratio of the surroundings of the clouds where mixing 884 !--increases the size of the cloud, to the total surroundings of the clouds. 885 !--This implies that ( 1. - sigma_mix ) quantifies the ratio where mixing 886 !--decreases the size of the clouds 887 !--For aviation, we increase the chance that the air surrounding contrails 888 !--is moist. This is quantified with chi_mixing_contrails 889 sigma_mix = chi_mixing_contrails * pdf_fra_above_lim & 890 / ( pdf_fra_below_lim + chi_mixing_contrails * pdf_fra_above_lim ) 868 891 869 892 IF ( pdf_fra_above_lim .GT. eps ) THEN … … 883 906 884 907 !--Add tendencies 885 contfra(i) = MAX(0., MIN(1., contfra(i) + dcfa_mix(i)))886 qcont(i) = MAX(0., MIN(qtot(i), qcont(i) + dqta_mix(i)))887 clrfra(i) = MAX(0., MIN(1., clrfra(i) - dcfa_mix(i)))888 qclr(i) = MAX(0., MIN(qtot(i), qclr(i) - dqta_mix(i)))908 contfra(i) = contfra(i) + dcfa_mix(i) 909 qcont(i) = qcont(i) + dqta_mix(i) 910 clrfra(i) = clrfra(i) - dcfa_mix(i) 911 qclr(i) = qclr(i) - dqta_mix(i) 889 912 890 913 ENDIF ! ( contfra(i) .GT. eps ) .AND. ( clrfra(i) .GT. eps ) … … 908 931 dqvc_con(i) = dqvc_con(i) / dtime 909 932 dqvc_mix(i) = dqvc_mix(i) / dtime 933 934 !--Diagnose ISSRs 935 IF ( clrfra(i) .GT. eps ) THEN 936 !--We then calculate the part that is greater than qsat 937 !--and lower than gamma_cond * qsat, which is the ice supersaturated region 938 pdf_x = gamma_cond(i) * qsat(i) / qsatl(i) * 100. 939 pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i) 940 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 941 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2 942 pdf_fra_above_nuc = EXP( - pdf_y ) * clrfra(i) 943 pdf_q_above_nuc = ( pdf_e3 * clrfra(i) & 944 + pdf_loc(i) * pdf_fra_above_nuc ) * qsatl(i) / 100. 945 946 pdf_x = qsat(i) / qsatl(i) * 100. 947 pdf_y = ( MAX( pdf_x - pdf_loc(i), 0. ) / pdf_scale(i) ) ** pdf_alpha(i) 948 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha(i) , pdf_y ) 949 pdf_e3 = pdf_scale(i) * ( 1. - pdf_e3 ) * pdf_e2 950 issrfra(i) = EXP( - pdf_y ) * clrfra(i) 951 qissr(i) = ( pdf_e3 * clrfra(i) + pdf_loc(i) * issrfra(i) ) * qsatl(i) / 100. 952 953 issrfra(i) = issrfra(i) - pdf_fra_above_nuc 954 qissr(i) = qissr(i) - pdf_q_above_nuc 955 ENDIF 910 956 911 957 !------------------------------------------- … … 937 983 CALL contrails_formation( & 938 984 klon, dtime, pplay, temp, qsat, qsatl, gamma_cond, & 939 flight_dist, cl dfra, pdf_loc, pdf_scale, pdf_alpha, keepgoing, &985 flight_dist, clrfra, pdf_loc, pdf_scale, pdf_alpha, keepgoing, & 940 986 Tcritcont, qcritcont, potcontfraP, potcontfraNP, & 941 987 dcfa_ini, dqia_ini, dqta_ini) … … 951 997 issrfra(i) = MAX(0., issrfra(i) - dcfa_ini(i)) 952 998 qissr(i) = MAX(0., qissr(i) - dqta_ini(i)) 953 cldfra(i) = MAX(0., MIN(1. , cldfra(i) + dcfa_ini(i)))999 cldfra(i) = MAX(0., MIN(1. - cldfracv(i), cldfra(i) + dcfa_ini(i))) 954 1000 qcld(i) = MAX(0., MIN(qtot(i), qcld(i) + dqta_ini(i))) 955 1001 qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dcfa_ini(i) * qsat(i))) … … 1135 1181 qvapincld_depsub = qsat + ( qvapincld - qsat ) & 1136 1182 * EXP( - dtime * qiceincld / kappa / rm_ice**2 ) 1183 IF ( qvapincld_depsub .GE. ( qvapincld + qiceincld ) ) & 1184 qvapincld_depsub = 0. 1137 1185 ENDIF ! qvapincld .GT. qsat 1138 1186 -
LMDZ6/branches/contrails/libf/phylmd/lmdz_lscp_ini.f90
r5589 r5602 200 200 REAL, SAVE, PROTECTED :: coef_shear_lscp=0.72 ! [-] additional coefficient for the shearing process (subprocess of the mixing process) 201 201 !$OMP THREADPRIVATE(coef_shear_lscp) 202 203 REAL, SAVE, PROTECTED :: chi_mixing_lscp=1. ! [-] factor for the macro distribution of ISSRs wrt clouds in a gridbox204 !$OMP THREADPRIVATE(chi_mixing_lscp)205 202 !--End of the parameters for condensation and ice supersaturation 206 203 … … 219 216 !$OMP THREADPRIVATE(coef_shear_contrails) 220 217 221 REAL, SAVE, PROTECTED :: chi_mixing_contrails ! [-] factor for the macro distribution of ISSRs wrt clouds in a gridbox218 REAL, SAVE, PROTECTED :: chi_mixing_contrails=1. ! [-] factor for increasing the chance that moist air is surrounding contrails 222 219 !$OMP THREADPRIVATE(chi_mixing_contrails) 223 220 … … 486 483 CALL getin_p('coef_mixing_lscp',coef_mixing_lscp) 487 484 CALL getin_p('coef_shear_lscp',coef_shear_lscp) 488 CALL getin_p('chi_mixing_lscp',chi_mixing_lscp)489 485 ! for aviation 490 486 CALL getin_p('aspect_ratio_contrails',aspect_ratio_contrails) … … 493 489 coef_shear_contrails=coef_shear_lscp 494 490 CALL getin_p('coef_shear_contrails',coef_shear_contrails) 495 chi_mixing_contrails=chi_mixing_lscp496 491 CALL getin_p('chi_mixing_contrails',chi_mixing_contrails) 497 492 CALL getin_p('rm_ice_crystals_contrails',rm_ice_crystals_contrails) … … 585 580 WRITE(lunout,*) 'lscp_ini, coef_mixing_lscp:', coef_mixing_lscp 586 581 WRITE(lunout,*) 'lscp_ini, coef_shear_lscp:', coef_shear_lscp 587 WRITE(lunout,*) 'lscp_ini, chi_mixing_lscp:', chi_mixing_lscp588 582 ! for aviation 589 583 WRITE(lunout,*) 'lscp_ini, aspect_ratio_contrails:', aspect_ratio_contrails -
LMDZ6/branches/contrails/libf/phylmd/physiq_mod.F90
r5601 r5602 3945 3945 3946 3946 CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay,omega, & 3947 t_seri, q_seri, ql_seri_lscp, qi_seri_lscp, ptconv, ratqs, sigma_qtherm, & 3947 t_seri, q_seri, ql_seri_lscp, qi_seri_lscp, & 3948 ptconv, rnebcon0, clwcon0, ratqs, sigma_qtherm, & 3948 3949 d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, & 3949 3950 pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb, &
Note: See TracChangeset
for help on using the changeset viewer.
