- Timestamp:
- Feb 23, 2024, 8:07:46 PM (9 months ago)
- Location:
- LMDZ6/trunk/libf/phylmd
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90
r4830 r4832 161 161 !$OMP THREADPRIVATE(rho_rain) 162 162 163 REAL, SAVE, PROTECTED :: rho_ice=920. ! A COMMENTER TODO [kg/m3] 164 !$OMP THREADPRIVATE(rho_ice) 165 163 166 REAL, SAVE, PROTECTED :: rho_snow ! A COMMENTER TODO [kg/m3] 164 167 !$OMP THREADPRIVATE(rho_snow) … … 169 172 REAL, SAVE, PROTECTED :: r_snow=1.E-3 ! A COMMENTER TODO [m] 170 173 !$OMP THREADPRIVATE(r_snow) 174 175 REAL, SAVE, PROTECTED :: r_ice=50.E-6 ! A COMMENTER TODO [m] 176 !$OMP THREADPRIVATE(r_ice) 171 177 172 178 REAL, SAVE, PROTECTED :: tau_auto_snow_min=100. ! A COMMENTER TODO [s] -
LMDZ6/trunk/libf/phylmd/lmdz_lscp_poprecip.F90
r4830 r4832 136 136 !--multiplying by dz = - dP / g / rho (line 3-4) 137 137 !--formula from Sundqvist 1988, Klemp & Wilhemson 1978 138 139 138 draineva = - precipfracclr(i) * coef_eva * (1. - qvap(i) / qsatl(i)) & 140 139 * ( rainclr(i) / MAX(thresh_precip_frac, precipfracclr(i)) ) ** expo_eva & … … 197 196 198 197 !--Add tendencies 199 rainclr(i) = rainclr(i) + draineva 200 snowclr(i) = snowclr(i) + dsnowsub 198 199 rainclr(i) = MAX(0., rainclr(i) + draineva) 200 snowclr(i) = MAX(0., snowclr(i) + dsnowsub) 201 201 202 202 !--If there is no more precip fluxes, the precipitation fraction in clear … … 248 248 cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, rain_int_min, & 249 249 thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, & 250 rho_rain, rho_snow, r_rain, r_snow, 250 rho_rain, rho_snow, r_rain, r_snow, rho_ice, r_ice, & 251 251 tau_auto_snow_min, tau_auto_snow_max, & 252 252 thresh_precip_frac, eps, air_thermal_conduct, & … … 340 340 REAL :: dqrfreez_max 341 341 REAL :: tau_freez 342 REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez 342 REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez, dqrtotfreez_step1, dqrtotfreez_step2 343 343 REAL :: coef_freez 344 344 REAL :: dqifreez !--loss of ice cloud content due to collection of ice from rain [kg/kg/s] … … 371 371 dqifreez= 0. 372 372 373 !--hum_to_flux: coef to convert a specific quantity to a flux373 !--hum_to_flux: coef to convert a delta in specific quantity to a flux 374 374 !-- hum_to_flux = rho * dz/dt = 1 / g * dP/dt 375 375 hum_to_flux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime … … 398 398 !--evaporated (see barrier at the end of poprecip_precld) 399 399 !--NB. precipfraccld(i) is here the cloud fraction of the layer above 400 precipfractot = 1. - ( 1. - precipfractot ) * & 400 !precipfractot = 1. - ( 1. - precipfractot ) * & 401 ! ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) & 402 ! / ( 1. - MIN( precipfraccld(i), 1. - eps ) ) 403 404 405 IF ( precipfraccld(i) .GT. ( 1. - eps ) ) THEN 406 precipfractot = 1. 407 ELSE 408 precipfractot = 1. - ( 1. - precipfractot ) * & 401 409 ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) & 402 / ( 1. - MIN( precipfraccld(i), 1. - eps) )403 410 / ( 1. - precipfraccld(i) ) 411 ENDIF 404 412 405 413 !--precipfraccld(i) is here the cloud fraction of the layer above … … 451 459 precipfraccld(i) = precipfraccld(i) + dprecipfraccld 452 460 precipfracclr(i) = precipfracclr(i) + dprecipfracclr 453 rainclr(i) = rainclr(i) + drainclr 454 snowclr(i) = snowclr(i) + dsnowclr 455 raincld(i) = raincld(i) - drainclr 456 snowcld(i) = snowcld(i) - dsnowclr 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 ) 457 466 458 459 467 !--If vertical heterogeneity is taken into account, we use 460 468 !--the "true" volume fraction instead of a modified … … 513 521 !--Add tendencies 514 522 qliq(i) = qliq(i) + dqlcol 515 raincld(i) = raincld(i) - dqlcol * hum_to_flux(i)523 raincld(i) = MAX( 0. , raincld(i) - dqlcol * hum_to_flux(i) ) 516 524 517 525 !--Diagnostic tendencies … … 535 543 !--Add tendencies 536 544 qice(i) = qice(i) + dqiagg 537 snowcld(i) = snowcld(i) - dqiagg * hum_to_flux(i)545 snowcld(i) = MAX( 0. , snowcld(i) - dqiagg * hum_to_flux(i) ) 538 546 539 547 !--Diagnostic tendencies … … 599 607 qliq(i) = qliq(i) + dqlauto 600 608 qice(i) = qice(i) + dqiauto 601 raincld(i) = raincld(i) - dqlauto * hum_to_flux(i)602 snowcld(i) = snowcld(i) - dqiauto * hum_to_flux(i)609 raincld(i) = MAX( 0. , raincld(i) - dqlauto * hum_to_flux(i) ) 610 snowcld(i) = MAX( 0. , snowcld(i) - dqiauto * hum_to_flux(i) ) 603 611 604 612 !--Diagnostic tendencies … … 631 639 !--Add tendencies 632 640 qliq(i) = qliq(i) + dqlrim 633 snowcld(i) = snowcld(i) - dqlrim * hum_to_flux(i)641 snowcld(i) = MAX( 0. , snowcld(i) - dqlrim * hum_to_flux(i) ) 634 642 635 643 !--Temperature adjustment with the release of latent … … 699 707 !--Barrier to limit the melting flux to the clr snow flux in the mesh 700 708 dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / hum_to_flux(i)) 701 ENDIF709 ENDIF 702 710 703 711 !--In cloudy air … … 726 734 dqscldmelt = dqsmelt_max * dqscldmelt / dqstotmelt 727 735 dqstotmelt = dqsmelt_max 736 728 737 ENDIF 729 738 730 739 !--Add tendencies 731 rainclr(i) = rainclr(i) - dqsclrmelt * hum_to_flux(i) 732 raincld(i) = raincld(i) - dqscldmelt * hum_to_flux(i) 733 snowclr(i) = snowclr(i) + dqsclrmelt * hum_to_flux(i) 734 snowcld(i) = snowcld(i) + dqscldmelt * hum_to_flux(i) 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) ) 745 735 746 736 747 !--Temperature adjustment with the release of latent … … 749 760 !-- FREEZING 750 761 !--------------------------------------------------------- 751 !--Process through which rain freezes into snow. This is 752 !--parameterized as an exponential decrease of the rain 753 !--water content. 754 !--The formula is homemade. 762 !--Process through which rain freezes into snow. 763 !-- We parameterize it as a 2 step process: 764 !-- first: freezing following collision with ice crystals 765 !-- second: immersion freezing following (inspired by Bigg 1953) 766 !-- the latter is parameterized as an exponential decrease of the rain 767 !-- water content with a homemade formulya 755 768 !--This is based on a caracteritic time of freezing, which 756 769 !--exponentially depends on temperature so that it is … … 758 771 !--0 for RTT (=273.15 K). 759 772 !--NB.: this process needs a temperature adjustment 760 761 773 !--dqrfreez_max: maximum rain freezing so that temperature 762 774 !-- stays lower than 273 K [kg/kg] … … 768 780 IF ( ( rainclr(i) + raincld(i) ) .GT. 0. ) THEN 769 781 770 !--Computed according to 771 !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q 772 dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD & 782 783 !--1st step: freezing following collision with ice crystals 784 !--Sub-process of freezing which quantifies the collision between 785 !--ice crystals in suspension and falling rain droplets. 786 !--The rain droplets freeze, becoming graupel, and carrying 787 !--the ice crystal (which acted as an ice nucleating particle). 788 !--The formula is adapted from the riming formula. 789 !-- it works only in the cloudy part 790 791 dqifreez = 0. 792 dqrtotfreez_step1 = 0. 793 794 IF ((qice(i) .GT. 0.) .AND. (raincld(i) .GT. 0.)) THEN 795 dqrclrfreez = 0. 796 dqrcldfreez = 0. 797 798 !--Computed according to 799 !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q 800 dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD & 773 801 * ( 1. + RVTMP2 * qtot(i) )) 774 802 775 tau_freez = 1. / ( beta_freez &776 * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) )777 778 !--Initialisation779 dqrclrfreez = 0.780 dqrcldfreez = 0.781 782 !--In clear air783 IF ( rainclr(i) .GT. 0. ) THEN784 !--Exact solution of dqrain/dt = -qrain/tau_freez785 dqrclrfreez = rainclr(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )786 ENDIF787 788 !--In cloudy air789 IF ( raincld(i) .GT. 0. ) THEN790 !--Exact solution of dqrain/dt = -qrain/tau_freez791 dqrcldfreez = raincld(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. )792 793 !---------------------------------------------------------794 !-- FREEZING OF RAIN WITH CONTACT OF ICE CRYSTALS795 !---------------------------------------------------------796 803 !--Sub-process of freezing which quantifies the collision between 797 804 !--ice crystals in suspension and falling rain droplets. … … 800 807 !--The formula is adapted from the riming formula. 801 808 802 !--The sticking efficacy is perfect.809 !--The collision efficiency is assumed unity 803 810 Eff_rain_ice = 1. 804 811 coef_freez = gamma_freez * 3. / 4. / rho_rain / r_rain * Eff_rain_ice 805 !-- ATTENTION Double implicit version? + barriers if needed806 812 !--Exact version, which does not need a barrier because of 807 !--the exponential decrease 813 !--the exponential decrease. 808 814 dqifreez = qice(i) * ( EXP( - dtime * coef_freez * raincld(i) / precipfraccld(i) ) - 1. ) 809 815 810 !--We add the two freezing processes in cloud 811 dqrcldfreez = dqrcldfreez - dqifreez 812 ENDIF 813 814 !--Barriers 816 !--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)) 820 dqrcldfreez = MAX(dqrcldfreez, dqrfreez_max) 821 dqrtotfreez_step1 = dqrcldfreez 822 823 !--Add tendencies 824 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) ) 827 temp(i) = temp(i) - dqrtotfreez_step1 * RLMLT / RCPD & 828 / ( 1. + RVTMP2 * qtot(i) ) 829 830 ENDIF 831 832 !-- Second step immersion freezing of rain drops 833 !-- with a homemade timeconstant depending on temperature 834 835 dqrclrfreez = 0. 836 dqrcldfreez = 0. 837 dqrtotfreez_step2 = 0. 838 !--Computed according to 839 !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q 840 841 dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD & 842 * ( 1. + RVTMP2 * qtot(i) )) 843 844 845 tau_freez = 1. / ( beta_freez & 846 * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) ) 847 848 849 !--In clear air 850 IF ( rainclr(i) .GT. 0. ) THEN 851 !--Exact solution of dqrain/dt = -qrain/tau_freez 852 dqrclrfreez = rainclr(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. ) 853 ENDIF 854 855 !--In cloudy air 856 IF ( raincld(i) .GT. 0. ) THEN 857 !--Exact solution of dqrain/dt = -qrain/tau_freez 858 dqrcldfreez = raincld(i) / hum_to_flux(i) * ( EXP( - dtime / tau_freez ) - 1. ) 859 ENDIF 860 861 !--temperature barrie step 2 815 862 !--It is activated if the total is LOWER than the max 816 863 !--because everything is negative 817 dqrtotfreez = dqrclrfreez + dqrcldfreez 818 IF ( dqrtotfreez .LT. dqrfreez_max ) THEN 864 dqrtotfreez_step2 = dqrclrfreez + dqrcldfreez 865 866 IF ( dqrtotfreez_step2 .LT. dqrfreez_max ) THEN 819 867 !--We redistribute the max freezed rain keeping 820 868 !--the clear/cloud partition of the freezing rain 821 dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez 822 dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez 823 dqrtotfreez = dqrfreez_max869 dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez_step2 870 dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez_step2 871 dqrtotfreez_step2 = dqrfreez_max 824 872 ENDIF 825 873 826 874 827 875 !--Add tendencies 828 rainclr(i) = rainclr(i) + dqrclrfreez * hum_to_flux(i) 829 raincld(i) = raincld(i) + dqrcldfreez * hum_to_flux(i) 830 snowclr(i) = snowclr(i) - dqrclrfreez * hum_to_flux(i) 831 snowcld(i) = snowcld(i) - dqrcldfreez * hum_to_flux(i) 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) ) 880 832 881 833 882 !--Temperature adjustment with the uptake of latent 834 883 !--heat because of freezing 835 temp(i) = temp(i) - dqrtotfreez * RLMLT / RCPD &884 temp(i) = temp(i) - dqrtotfreez_step2 * RLMLT / RCPD & 836 885 / ( 1. + RVTMP2 * qtot(i) ) 837 886 887 dqrtotfreez=dqrtotfreez_step1+dqrtotfreez_step2 838 888 !--Diagnostic tendencies 839 889 dqrfreez(i) = dqrtotfreez / dtime 840 dqsfreez(i) = - dqrtotfreez/ dtime890 dqsfreez(i) = -(dqrtotfreez + dqifreez) / dtime 841 891 842 892 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.