- Timestamp:
- Feb 27, 2024, 5:48:56 PM (9 months ago)
- File:
-
- 1 edited
-
LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90 (modified) (26 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90
r4832 r4833 63 63 !--Integer for interating over klon 64 64 INTEGER :: i 65 !--hum_to_flux: coef to convert a specific quantity to a flux 66 REAL, DIMENSION(klon) :: hum_to_flux 65 !--dhum_to_dflux: coef to convert a specific quantity to a flux 66 REAL, DIMENSION(klon) :: dhum_to_dflux 67 !-- 68 REAL, DIMENSION(klon) :: rho, dz 67 69 68 70 !--Saturation values … … 82 84 dqssubl = 0. 83 85 84 !-- hum_to_flux = rho * dz/dt = 1 / g * dP/dt 85 hum_to_flux(:) = ( paprsdn(:) - paprsup(:) ) / RG / dtime 86 !-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt 87 dhum_to_dflux(:) = ( paprsdn(:) - paprsup(:) ) / RG / dtime 88 rho(:) = pplay(:) / temp(:) / RD 89 dz(:) = ( paprsdn(:) - paprsup(:) ) / RG / rho(:) 86 90 87 91 !--Calculation of saturation specific humidity … … 117 121 !--same temperature as the lowermost layer 118 122 !--we convert the flux into a specific quantity qprecip 119 qprecip(i) = ( rain(i) + snow(i) ) / hum_to_flux(i)123 qprecip(i) = ( rain(i) + snow(i) ) / dhum_to_dflux(i) 120 124 !-- t(i,k+1) + d_t(i,k+1): new temperature of the overlying layer 121 125 temp(i) = ( tempupnew(i) * qprecip(i) * cpw + cpair * temp(i) ) & … … 133 137 !--Evaporation of liquid precipitation coming from above 134 138 !--in the clear sky only 135 !--dP/dz=beta*(1-q/qsat)*(P**expo_eva) (lines 1-2) 136 !--multiplying by dz = - dP / g / rho (line 3-4) 139 !--dprecip/dz = -beta*(1-qvap/qsat)*(precip**expo_eva) 137 140 !--formula from Sundqvist 1988, Klemp & Wilhemson 1978 138 draineva = - precipfracclr(i) * coef_eva * (1. - qvap(i) / qsatl(i)) & 139 * ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_eva & 140 * temp(i) * RD / pplay(i) & 141 * ( paprsdn(i) - paprsup(i) ) / RG 142 143 !--Evaporation is limited by 0 and by the total water amount in 144 !--the precipitation 145 draineva = MIN(0., MAX(draineva, -rainclr(i))) 141 !--Explicit formulation 142 !draineva = - precipfracclr(i) * coef_eva * (1. - qvap(i) / qsatl(i)) * dz(i) & 143 ! * ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_eva & 144 !draineva = MAX( - rainclr(i), draineva) 145 !--Exact explicit formulation (rainclr is resolved exactly, qvap explicitly) 146 !--which does not need a barrier on rainclr, because included in the formula 147 draineva = precipfracclr(i) * ( MAX(0., & 148 - coef_eva * ( 1. - expo_eva ) * (1. - qvap(i) / qsatl(i)) * dz(i) & 149 + ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) )**( 1. - expo_eva ) & 150 ) )**( 1. / ( 1. - expo_eva ) ) - rainclr(i) 151 152 !--Evaporation is limited by 0 153 draineva = MIN(0., draineva) 146 154 147 155 148 156 !--Sublimation of the solid precipitation coming from above 149 157 !--(same formula as for liquid precip) 150 dsnowsub = - precipfracclr(i) * coef_sub * (1. - qvap(i) / qsati(i)) & 151 * ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_sub & 152 * temp(i) * RD / pplay(i) & 153 * ( paprsdn(i) - paprsup(i) ) / RG 154 155 !--Sublimation is limited by 0 and by the total water amount in 156 !--the precipitation 158 !--Explicit formulation 159 !dsnowsub = - precipfracclr(i) * coef_sub * (1. - qvap(i) / qsatl(i)) * dz(i) & 160 ! * ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_sub & 161 !dsnowsub = MAX( - snowclr(i), dsnowsub) 162 !--Exact explicit formulation (snowclr is resolved exactly, qvap explicitly) 163 !--which does not need a barrier on snowclr, because included in the formula 164 dsnowsub = precipfracclr(i) * ( MAX(0., & 165 - coef_sub * ( 1. - expo_sub ) * (1. - qvap(i) / qsatl(i)) * dz(i) & 166 + ( snowclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) )**( 1. - expo_sub ) & 167 ) )**( 1. / ( 1. - expo_sub ) ) - snowclr(i) 168 169 !--Sublimation is limited by 0 157 170 ! TODO: change max when we will allow for vapor deposition in supersaturated regions 158 dsnowsub = MIN(0., MAX(dsnowsub, -snowclr(i)))171 dsnowsub = MIN(0., dsnowsub) 159 172 160 173 !--Evaporation limit: we ensure that the layer's fraction below … … 166 179 167 180 dprecip_evasub_max = MIN(0., ( qvap(i) - qsat(i) ) * precipfracclr(i)) & 168 * hum_to_flux(i)181 * dhum_to_dflux(i) 169 182 dprecip_evasub_tot = draineva + dsnowsub 170 183 … … 179 192 180 193 !--New solid and liquid precipitation fluxes after evap and sublimation 181 dqrevap = draineva / hum_to_flux(i)182 dqssubl = dsnowsub / hum_to_flux(i)194 dqrevap = draineva / dhum_to_dflux(i) 195 dqssubl = dsnowsub / dhum_to_dflux(i) 183 196 184 197 … … 196 209 197 210 !--Add tendencies 198 211 !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) 199 212 rainclr(i) = MAX(0., rainclr(i) + draineva) 200 213 snowclr(i) = MAX(0., snowclr(i) + dsnowsub) … … 306 319 307 320 INTEGER :: i 308 REAL, DIMENSION(klon) :: hum_to_flux321 REAL, DIMENSION(klon) :: dhum_to_dflux 309 322 REAL, DIMENSION(klon) :: qtot !--includes vap, liq, ice and precip 310 323 … … 371 384 dqifreez= 0. 372 385 373 !-- hum_to_flux: coef to convert a delta in specific quantity to a flux374 !-- hum_to_flux = rho * dz/dt = 1 / g * dP/dt375 hum_to_flux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime386 !--dhum_to_dflux: coef to convert a delta in specific quantity to a flux 387 !-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt 388 dhum_to_dflux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime 376 389 qtot(i) = qvap(i) + qliq(i) + qice(i) & 377 + ( raincld(i) + rainclr(i) + snowcld(i) + snowclr(i) ) / hum_to_flux(i)390 + ( raincld(i) + rainclr(i) + snowcld(i) + snowclr(i) ) / dhum_to_dflux(i) 378 391 379 392 !------------------------------------------------------------ … … 457 470 458 471 !--We add the tendencies 472 !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) 459 473 precipfraccld(i) = precipfraccld(i) + dprecipfraccld 460 474 precipfracclr(i) = precipfracclr(i) + dprecipfracclr 461 462 rainclr(i) = MAX( 0. , rainclr(i) + drainclr ) 463 snowclr(i) = MAX( 0. , snowclr(i) + dsnowclr ) 464 raincld(i) = MAX( 0. , raincld(i) - drainclr ) 465 snowcld(i) = MAX( 0. , snowcld(i) - dsnowclr ) 475 rainclr(i) = MAX(0., rainclr(i) + drainclr) 476 snowclr(i) = MAX(0., snowclr(i) + dsnowclr) 477 raincld(i) = MAX(0., raincld(i) - drainclr) 478 snowcld(i) = MAX(0., snowcld(i) - dsnowclr) 466 479 467 480 !--If vertical heterogeneity is taken into account, we use … … 501 514 !--then simplified. 502 515 503 !--The sticking efficacy is perfect.516 !--The collection efficiency is perfect. 504 517 Eff_rain_liq = 1. 505 518 coef_col = gamma_col * 3. / 4. / rho_rain / r_rain * Eff_rain_liq 506 IF ( (raincld(i) .GT. 0.) .AND. (coef_col .GT. 0.)) THEN519 IF ( raincld(i) .GT. 0. ) THEN 507 520 !-- ATTENTION Double implicit version 508 521 !-- BEWARE the formule below is FALSE (because raincld is a flux, not a delta flux) 509 !qrain_tmp = raincld(i) / hum_to_flux(i)522 !qrain_tmp = raincld(i) / dhum_to_dflux(i) 510 523 !coef_tmp = coef_col * dtime * ( qrain_tmp / precipfraccld(i) + qliq(i) / eff_cldfra ) 511 524 !dqlcol = qliq(i) * ( 1. / ( 1. + 0.5 * ( coef_tmp - 1. + SQRT( & … … 521 534 !--Add tendencies 522 535 qliq(i) = qliq(i) + dqlcol 523 raincld(i) = MAX( 0. , raincld(i) - dqlcol * hum_to_flux(i))536 raincld(i) = raincld(i) - dqlcol * dhum_to_dflux(i) 524 537 525 538 !--Diagnostic tendencies … … 532 545 Eff_snow_ice = 0.05 * EXP( 0.1 * ( temp(i) - RTT ) ) 533 546 coef_agg = gamma_agg * 3. / 4. / rho_snow / r_snow * Eff_snow_ice 534 IF ( (snowcld(i) .GT. 0.) .AND. (coef_agg .GT. 0.)) THEN547 IF ( snowcld(i) .GT. 0. ) THEN 535 548 !-- ATTENTION Double implicit version? 536 549 !--Barriers so that the processes do not consume more liquid/ice than … … 543 556 !--Add tendencies 544 557 qice(i) = qice(i) + dqiagg 545 snowcld(i) = MAX( 0. , snowcld(i) - dqiagg * hum_to_flux(i))558 snowcld(i) = snowcld(i) - dqiagg * dhum_to_dflux(i) 546 559 547 560 !--Diagnostic tendencies … … 607 620 qliq(i) = qliq(i) + dqlauto 608 621 qice(i) = qice(i) + dqiauto 609 raincld(i) = MAX( 0. , raincld(i) - dqlauto * hum_to_flux(i))610 snowcld(i) = MAX( 0. , snowcld(i) - dqiauto * hum_to_flux(i))622 raincld(i) = raincld(i) - dqlauto * dhum_to_dflux(i) 623 snowcld(i) = snowcld(i) - dqiauto * dhum_to_dflux(i) 611 624 612 625 !--Diagnostic tendencies … … 628 641 Eff_snow_liq = 0.2 629 642 coef_rim = gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq 630 IF ( (snowcld(i) .GT. 0.) .AND. (coef_rim .GT. 0.)) THEN643 IF ( snowcld(i) .GT. 0. ) THEN 631 644 !-- ATTENTION Double implicit version? 632 645 !--Barriers so that the processes do not consume more liquid than … … 638 651 639 652 !--Add tendencies 653 !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) 640 654 qliq(i) = qliq(i) + dqlrim 641 snowcld(i) = MAX( 0. , snowcld(i) - dqlrim * hum_to_flux(i))655 snowcld(i) = snowcld(i) - dqlrim * dhum_to_dflux(i) 642 656 643 657 !--Temperature adjustment with the release of latent … … 706 720 707 721 !--Barrier to limit the melting flux to the clr snow flux in the mesh 708 dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / hum_to_flux(i))722 dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / dhum_to_dflux(i)) 709 723 ENDIF 710 724 … … 720 734 721 735 !--Barrier to limit the melting flux to the cld snow flux in the mesh 722 dqscldmelt = MAX( dqscldmelt , -snowcld(i) / hum_to_flux(i))736 dqscldmelt = MAX(dqscldmelt , - snowcld(i) / dhum_to_dflux(i)) 723 737 ENDIF 724 738 … … 738 752 739 753 !--Add tendencies 740 741 rainclr(i) = MAX( 0. , rainclr(i) - dqsclrmelt * hum_to_flux(i))742 raincld(i) = MAX( 0. , raincld(i) - dqscldmelt * hum_to_flux(i))743 snowclr(i) = MAX( 0. , snowclr(i) + dqsclrmelt * hum_to_flux(i))744 snowcld(i) = MAX( 0. , snowcld(i) + dqscldmelt * hum_to_flux(i))754 !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) 755 rainclr(i) = MAX(0., rainclr(i) - dqsclrmelt * dhum_to_dflux(i)) 756 raincld(i) = MAX(0., raincld(i) - dqscldmelt * dhum_to_dflux(i)) 757 snowclr(i) = MAX(0., snowclr(i) + dqsclrmelt * dhum_to_dflux(i)) 758 snowcld(i) = MAX(0., snowcld(i) + dqscldmelt * dhum_to_dflux(i)) 745 759 746 760 … … 792 806 dqrtotfreez_step1 = 0. 793 807 794 IF ( (qice(i) .GT. 0.) .AND. (raincld(i) .GT. 0.)) THEN808 IF ( ( qice(i) .GT. 0. ) .AND. ( raincld(i) .GT. 0. ) ) THEN 795 809 dqrclrfreez = 0. 796 810 dqrcldfreez = 0. … … 815 829 816 830 !--We add the part of rain water that freezes, limited by a temperature barrier 817 !-- this quantity is calculated assuming that the number of drop that freeze correspond to the number 818 !-- of crystals collected (and assuming uniform distributions of ice crystals and rain drops) 819 dqrcldfreez = MAX(dqifreez*rho_rain*r_rain/(rho_ice*r_ice),-raincld(i)/hum_to_flux(i)) 831 !--this quantity is calculated assuming that the number of drop that freeze correspond to the number 832 !--of crystals collected (and assuming uniform distributions of ice crystals and rain drops) 833 !--The ice specific humidity that collide with rain is dqi = dNi 4/3 PI rho_ice r_ice**3 834 !--The rain that collide with ice is, similarly, dqr = dNr 4/3 PI rho_rain r_rain**3 835 !--The assumption above corresponds to dNi = dNr, i.e., 836 !-- dqr = dqi * (4/3 PI rho_rain * r_rain**3) / (4/3 PI rho_ice * r_ice**3) 837 dqrcldfreez = dqifreez * rho_rain * r_rain**3. / ( rho_ice * r_ice**3. ) 838 dqrcldfreez = MAX(dqrcldfreez, - raincld(i) / dhum_to_dflux(i)) 820 839 dqrcldfreez = MAX(dqrcldfreez, dqrfreez_max) 821 840 dqrtotfreez_step1 = dqrcldfreez 822 841 823 842 !--Add tendencies 843 !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) 824 844 qice(i) = qice(i) + dqifreez 825 raincld(i) = MAX( 0. , raincld(i) + dqrcldfreez * hum_to_flux(i))826 snowcld(i) = MAX( 0. , snowcld(i) - dqrcldfreez * hum_to_flux(i) - dqifreez * hum_to_flux(i))845 raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i)) 846 snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i) - dqifreez * dhum_to_dflux(i)) 827 847 temp(i) = temp(i) - dqrtotfreez_step1 * RLMLT / RCPD & 828 / ( 1. + RVTMP2 * qtot(i) )848 / ( 1. + RVTMP2 * qtot(i) ) 829 849 830 850 ENDIF … … 850 870 IF ( rainclr(i) .GT. 0. ) THEN 851 871 !--Exact solution of dqrain/dt = -qrain/tau_freez 852 dqrclrfreez = rainclr(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )872 dqrclrfreez = rainclr(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. ) 853 873 ENDIF 854 874 … … 856 876 IF ( raincld(i) .GT. 0. ) THEN 857 877 !--Exact solution of dqrain/dt = -qrain/tau_freez 858 dqrcldfreez = raincld(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )878 dqrcldfreez = raincld(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. ) 859 879 ENDIF 860 880 … … 874 894 875 895 !--Add tendencies 876 rainclr(i) = MAX( 0. , rainclr(i) + dqrclrfreez * hum_to_flux(i) ) 877 raincld(i) = MAX( 0. , raincld(i) + dqrcldfreez * hum_to_flux(i) ) 878 snowclr(i) = MAX( 0. , snowclr(i) - dqrclrfreez * hum_to_flux(i) ) 879 snowcld(i) = MAX( 0. , snowcld(i) - dqrcldfreez * hum_to_flux(i) ) 896 !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) 897 rainclr(i) = MAX(0., rainclr(i) + dqrclrfreez * dhum_to_dflux(i)) 898 raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i)) 899 snowclr(i) = MAX(0., snowclr(i) - dqrclrfreez * dhum_to_dflux(i)) 900 snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i)) 880 901 881 902 … … 885 906 / ( 1. + RVTMP2 * qtot(i) ) 886 907 887 dqrtotfreez=dqrtotfreez_step1+dqrtotfreez_step2888 908 !--Diagnostic tendencies 909 dqrtotfreez = dqrtotfreez_step1 + dqrtotfreez_step2 889 910 dqrfreez(i) = dqrtotfreez / dtime 890 911 dqsfreez(i) = -(dqrtotfreez + dqifreez) / dtime
Note: See TracChangeset
for help on using the changeset viewer.
