Changeset 4226 for LMDZ6/trunk
- Timestamp:
- Jul 23, 2022, 6:04:45 PM (2 years ago)
- File:
-
- 1 edited
-
LMDZ6/trunk/libf/phylmd/lscp_mod.F90 (modified) (65 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/lscp_mod.F90
r4114 r4226 23 23 !-------------------------------------------------------------------------------- 24 24 ! Aim: Large Scale Clouds and Precipitation (LSCP) 25 ! 25 ! 26 26 ! This code is a new version of the fisrtilp.F90 routine, which is itself a 27 27 ! merge of 'first' (superrsaturation physics, P. LeVan K. Laval) … … 46 46 ! Code structure: 47 47 ! 48 ! P0> Thermalization of the 48 ! P0> Thermalization of the precipitation coming from the overlying layer 49 49 ! P1> Evaporation of the precipitation (falling from the k+1 level) 50 50 ! P2> Cloud formation (at the k level) … … 86 86 ! (which is not trivial) 87 87 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 88 88 ! 89 89 USE dimphy 90 90 USE print_control_mod, ONLY: prt_level, lunout … … 112 112 include "nuage.h" 113 113 114 115 114 ! INPUT VARIABLES: 116 115 !----------------- … … 125 124 INTEGER, INTENT(IN) :: iflag_cld_th ! flag that determines the distribution of convective clouds 126 125 INTEGER, INTENT(IN) :: iflag_ice_thermo! flag to activate the ice thermodynamics 127 ! CR: if iflag_ice_thermo=2, only convection is active 126 ! CR: if iflag_ice_thermo=2, only convection is active 128 127 LOGICAL, INTENT(IN) :: ok_ice_sursat ! flag to determine if ice sursaturation is activated 129 128 130 129 LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active 131 130 132 131 133 132 !Inputs associated with thermal plumes 134 133 135 134 REAL, DIMENSION(klon,klev), INTENT(IN) :: ztv ! virtual potential temperature [K] 136 135 REAL, DIMENSION(klon,klev), INTENT(IN) :: zqta ! specific humidity within thermals [kg/kg] … … 144 143 145 144 REAL, DIMENSION(klon,klev), INTENT(INOUT):: ratqs ! function of pressure that sets the large-scale 146 ! cloud PDF (sigma=ratqs*qt)147 148 145 149 146 ! Input sursaturation en glace … … 166 163 REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: psfl ! large-scale snowfall flux in the column [kg/s/m2] 167 164 168 169 165 ! cofficients of scavenging fraction (for off-line) 170 166 … … 178 174 REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_nucl ! scavenging fraction due tu nucleation [-] 179 175 180 181 182 183 176 ! PROGRAM PARAMETERS: 184 177 !-------------------- 185 178 186 187 REAL, SAVE :: seuil_neb=0.001 ! cloud fraction threshold: a cloud really exists when exceeded 179 REAL, SAVE :: seuil_neb=0.001 ! cloud fraction threshold: a cloud really exists when exceeded 188 180 !$OMP THREADPRIVATE(seuil_neb) 189 181 190 INTEGER, SAVE :: ninter=5 182 INTEGER, SAVE :: ninter=5 ! number of iterations to calculate autoconversion to precipitation 191 183 !$OMP THREADPRIVATE(ninter) 192 184 193 INTEGER,SAVE :: iflag_evap_prec=1 ! precipitation evap. flag. 0: nothing, 1: "old way", 2: Max cloud fraction above to calculate the max of reevaporation 194 ! 4: LTP'method i.e. evaporation in the clear-sky fraction of the mesh only 185 INTEGER,SAVE :: iflag_evap_prec=1 ! precipitation evaporation flag. 0: nothing, 1: "old way", 186 ! 2: Max cloud fraction above to calculate the max of reevaporation 187 ! 4: LTP'method i.e. evaporation in the clear-sky fraction of the mesh only 195 188 !$OMP THREADPRIVATE(iflag_evap_prec) 196 197 198 REAL t_coup ! temperature threshold which determines the phase 199 PARAMETER (t_coup=234.0) ! for which the saturation vapor pressure is calculated 200 201 REAL DDT0 ! iteration parameter 189 190 REAL t_coup ! temperature threshold which determines the phase 191 PARAMETER (t_coup=234.0) ! for which the saturation vapor pressure is calculated 192 193 REAL DDT0 ! iteration parameter 202 194 PARAMETER (DDT0=.01) 203 195 204 REAL ztfondue 196 REAL ztfondue ! parameter to calculate melting fraction of precipitation 205 197 PARAMETER (ztfondue=278.15) 206 198 207 REAL, SAVE :: rain_int_min=0.001 199 REAL, SAVE :: rain_int_min=0.001 ! Minimum local rain intensity [mm/s] before the decrease in associated precipitation fraction 208 200 !$OMP THREADPRIVATE(rain_int_min) 209 210 201 211 202 … … 218 209 REAL ctot(klon,klev) 219 210 REAL ctot_vol(klon,klev) 220 REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 211 REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 221 212 REAL zdqsdT_raw(klon) 222 REAL gammasat(klon),dgammasatdt(klon) 213 REAL gammasat(klon),dgammasatdt(klon) ! coefficient to make cold condensation at the correct RH and derivative wrt T 223 214 REAL Tbef(klon),qlbef(klon),DT(klon),num,denom 224 215 REAL cste 225 216 REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon) 226 217 REAL Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon) 227 REAL erf 218 REAL erf 228 219 REAL qcloud(klon), icefrac_mpc(klon,klev), qincloud_mpc(klon) 229 220 REAL zrfl(klon), zrfln(klon), zqev, zqevt … … 248 239 ! temporary: alpha parameter for scavenging 249 240 ! 4 possible scavenging processes 250 REAL a_tr_sca(4) 251 save a_tr_sca 252 !$OMP THREADPRIVATE(a_tr_sca) 241 REAL,SAVE :: a_tr_sca(4) 242 !$OMP THREADPRIVATE(a_tr_sca) 253 243 REAL zalpha_tr 254 244 REAL zfrac_lessi … … 257 247 REAL zmair(klon), zcpair, zcpeau 258 248 REAL zlh_solid(klon), zm_solid ! for liquid -> solid conversion 259 REAL zrflclr(klon), zrflcld(klon) 260 REAL d_zrfl_clr_cld(klon), d_zifl_clr_cld(klon) 249 REAL zrflclr(klon), zrflcld(klon) 250 REAL d_zrfl_clr_cld(klon), d_zifl_clr_cld(klon) 261 251 REAL d_zrfl_cld_clr(klon), d_zifl_cld_clr(klon) 262 REAL ziflclr(klon), ziflcld(klon) 263 REAL znebprecipclr(klon), znebprecipcld(klon) 264 REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon) 252 REAL ziflclr(klon), ziflcld(klon) 253 REAL znebprecipclr(klon), znebprecipcld(klon) 254 REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon) 265 255 REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon) 266 256 REAL velo(klon,klev), vr … … 270 260 INTEGER i, k, n, kk 271 261 INTEGER n_i(klon), iter 272 INTEGER ncoreczq 262 INTEGER ncoreczq 273 263 INTEGER mpc_bl_points(klon,klev) 274 264 INTEGER,SAVE :: itap=0 … … 277 267 LOGICAL lognormale(klon) 278 268 LOGICAL convergence(klon) 279 LOGICAL appel1er 280 SAVE appel1er 269 LOGICAL,SAVE :: appel1er 281 270 !$OMP THREADPRIVATE(appel1er) 282 271 … … 290 279 !=============================================================================== 291 280 292 ! Few initial checks S281 ! Few initial checks 293 282 294 283 IF (iflag_t_glace.EQ.0) THEN … … 303 292 304 293 IF (.NOT.thermcep) THEN 305 306 294 abort_message = 'lscp cannot be used when thermcep=false' 295 CALL abort_physic(modname,abort_message,1) 307 296 ENDIF 308 297 309 298 IF (iflag_fisrtilp_qsat .LT. 0) THEN 310 311 299 abort_message = 'lscp cannot be used with iflag_fisrtilp<0' 300 CALL abort_physic(modname,abort_message,1) 312 301 ENDIF 313 302 314 315 303 ! Few initialisations 316 304 317 305 itap=itap+1 318 znebprecip(:)=0. 306 znebprecip(:)=0.0 319 307 ctot_vol(1:klon,1:klev)=0.0 320 308 rneblsvol(1:klon,1:klev)=0.0 321 smallestreal=1.e-9 322 znebprecipclr(:)=0. 323 znebprecipcld(:)=0. 309 smallestreal=1.e-9 310 znebprecipclr(:)=0.0 311 znebprecipcld(:)=0.0 324 312 mpc_bl_points(:,:)=0 325 313 … … 341 329 WRITE(lunout,*) 'I would prefer a 6 min sub-timestep' 342 330 ENDIF 343 331 344 332 !AA Temporary initialisation 345 333 a_tr_sca(1) = -0.5 … … 351 339 DO k = 1, klev 352 340 DO i = 1, klon 353 pfrac_nucl(i,k)=1. 354 pfrac_1nucl(i,k)=1. 355 pfrac_impa(i,k)=1. 356 beta(i,k)=0. 341 pfrac_nucl(i,k)=1.0 342 pfrac_1nucl(i,k)=1.0 343 pfrac_impa(i,k)=1.0 344 beta(i,k)=0.0 357 345 ENDDO 358 346 ENDDO … … 362 350 appel1er = .FALSE. 363 351 364 ENDIF !(appel1er) 352 ENDIF !(appel1er) 365 353 366 354 ! AA for 'safety' reasons … … 379 367 radliq(:,:) = 0.0 380 368 radicefrac(:,:) = 0.0 381 frac_nucl(:,:) = 1. 382 frac_impa(:,:) = 1. 369 frac_nucl(:,:) = 1.0 370 frac_impa(:,:) = 1.0 383 371 384 372 rain(:) = 0.0 … … 393 381 ziflprev(:)=0.0 394 382 zneb(:) = seuil_neb 395 zrflclr(:) = 0.0 396 ziflclr(:) = 0.0 397 zrflcld(:) = 0.0 398 ziflcld(:) = 0.0 399 tot_zneb(:) = 0.0 400 tot_znebn(:) = 0.0 401 d_tot_zneb(:) = 0.0 402 qzero(:) =0.0383 zrflclr(:) = 0.0 384 ziflclr(:) = 0.0 385 zrflcld(:) = 0.0 386 ziflcld(:) = 0.0 387 tot_zneb(:) = 0.0 388 tot_znebn(:) = 0.0 389 d_tot_zneb(:) = 0.0 390 qzero(:) = 0.0 403 391 404 392 !--ice sursaturation 405 gamma_ss(:,:) = 1. 406 qss(:,:) = 0. 407 rnebss(:,:) = 0. 393 gamma_ss(:,:) = 1.0 394 qss(:,:) = 0.0 395 rnebss(:,:) = 0.0 408 396 Tcontr(:,:) = missing_val 409 397 qcontr(:,:) = missing_val … … 419 407 420 408 DO k = klev, 1, -1 421 422 409 423 410 ! Initialisation temperature and specific humidity … … 447 434 ! RVTMP2=rcpv/rcpd-1 448 435 zcpair=RCPD*(1.0+RVTMP2*zq(i)) 449 zcpeau=RCPD*RVTMP2 450 436 zcpeau=RCPD*RVTMP2 437 451 438 ! zmqc: precipitation mass that has to be thermalized with 452 439 ! layer's air so that precipitation at the ground has the … … 456 443 zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) & 457 444 / (zcpair + zmqc(i)*zcpeau) 458 445 459 446 ENDDO 460 447 461 ELSE 448 ELSE 462 449 463 450 DO i = 1, klon … … 474 461 ! -------------------------------------------------------------------- 475 462 476 477 463 IF (iflag_evap_prec.GE.1) THEN 478 464 … … 486 472 487 473 DO i = 1, klon 488 474 489 475 ! if precipitation 490 476 IF (zrfl(i)+zifl(i).GT.0.) THEN 491 477 492 478 ! LTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec=4). 493 IF (iflag_evap_prec.EQ.4) THEN 494 zrfl(i) = zrflclr(i) 495 zifl(i) = ziflclr(i) 496 ENDIF 479 IF (iflag_evap_prec.EQ.4) THEN 480 zrfl(i) = zrflclr(i) 481 zifl(i) = ziflclr(i) 482 ENDIF 497 483 498 484 IF (iflag_evap_prec.EQ.1) THEN … … 501 487 znebprecip(i)=MAX(zneb(i),znebprecip(i)) 502 488 ENDIF 503 504 505 IF (iflag_evap_prec.EQ.4) THEN 506 ! Max evaporation not to saturate the whole mesh 507 zqev0 = MAX(0.0, zqs(i)-zq(i)) 508 ELSE 509 ! Max evap not to saturate the fraction below the cloud 510 zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i)) 511 ENDIF 512 489 490 IF (iflag_evap_prec.EQ.4) THEN 491 ! Max evaporation not to saturate the whole mesh 492 zqev0 = MAX(0.0, zqs(i)-zq(i)) 493 ELSE 494 ! Max evap not to saturate the fraction below the cloud 495 zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i)) 496 ENDIF 513 497 514 498 ! Evaporation of liquid precipitation coming from above … … 530 514 ENDIF 531 515 532 533 516 zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) & 534 *RG*dtime/(paprs(i,k)-paprs(i,k+1))535 517 *RG*dtime/(paprs(i,k)-paprs(i,k+1)) 518 536 519 ! sublimation of the solid precipitation coming from above 537 520 IF (iflag_evap_prec.EQ.3) THEN … … 547 530 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG 548 531 ENDIF 549 532 550 533 zqevti = MAX(0.0,MIN(zqevti,zifl(i))) & 551 *RG*dtime/(paprs(i,k)-paprs(i,k+1)) 534 *RG*dtime/(paprs(i,k)-paprs(i,k+1)) 552 535 553 536 ! A. JAM … … 556 539 ! does not reach saturation. In this case, we 557 540 ! redistribute zqev0 conserving the ratio liquid/ice 558 541 559 542 IF (zqevt+zqevti.GT.zqev0) THEN 560 543 zqev=zqev0*zqevt/(zqevt+zqevti) … … 589 572 zifl(i) = zifln(i) 590 573 591 IF (iflag_evap_prec.EQ.4) THEN 592 zrflclr(i) = zrfl(i) 593 ziflclr(i) = zifl(i) 594 IF(zrflclr(i) + ziflclr(i).LE.0) THEN 595 znebprecipclr(i) = 0. 596 ENDIF 597 zrfl(i) = zrflclr(i) + zrflcld(i) 598 zifl(i) = ziflclr(i) + ziflcld(i) 599 ENDIF 574 IF (iflag_evap_prec.EQ.4) THEN 575 zrflclr(i) = zrfl(i) 576 ziflclr(i) = zifl(i) 577 IF(zrflclr(i) + ziflclr(i).LE.0) THEN 578 znebprecipclr(i) = 0.0 579 ENDIF 580 zrfl(i) = zrflclr(i) + zrflcld(i) 581 zifl(i) = ziflclr(i) + ziflcld(i) 582 ENDIF 600 583 601 584 … … 605 588 zmelt = MIN(MAX(zmelt,0.),1.) 606 589 ! Ice melting 607 IF (iflag_evap_prec.EQ.4) THEN 608 zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i) 609 zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i) 610 zrfl(i)=zrflclr(i)+zrflcld(i) 611 ELSE 612 zrfl(i)=zrfl(i)+zmelt*zifl(i) 613 ENDIF 590 IF (iflag_evap_prec.EQ.4) THEN 591 zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i) 592 zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i) 593 zrfl(i)=zrflclr(i)+zrflcld(i) 594 ELSE 595 zrfl(i)=zrfl(i)+zmelt*zifl(i) 596 ENDIF 614 597 615 598 ! Latent heat of melting with precipitation thermalization … … 617 600 *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) 618 601 619 IF (iflag_evap_prec.EQ.4) THEN 620 ziflclr(i)=ziflclr(i)*(1.-zmelt) 621 ziflcld(i)=ziflcld(i)*(1.-zmelt) 622 zifl(i)=ziflclr(i)+ziflcld(i) 623 ELSE 624 zifl(i)=zifl(i)*(1.-zmelt) 602 IF (iflag_evap_prec.EQ.4) THEN 603 ziflclr(i)=ziflclr(i)*(1.-zmelt) 604 ziflcld(i)=ziflcld(i)*(1.-zmelt) 605 zifl(i)=ziflclr(i)+ziflcld(i) 606 ELSE 607 zifl(i)=zifl(i)*(1.-zmelt) 625 608 ENDIF 626 609 … … 632 615 633 616 ENDDO 634 617 635 618 ENDIF ! (iflag_evap_prec>=1) 636 619 … … 638 621 ! End precip evaporation 639 622 ! -------------------------------------------------------------------- 640 641 642 643 623 644 624 ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter 645 625 !------------------------------------------------------- 646 626 647 627 qtot(:)=zq(:)+zmqc(:) 648 628 CALL CALC_QSAT_ECMWF(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:)) … … 657 637 ENDDO 658 638 659 660 639 ! -------------------------------------------------------------------- 661 640 ! P2> Cloud formation … … 670 649 ! ------------------------------------------------------------------- 671 650 672 673 651 ! P2.1> With the PDFs (log-normal, bigaussian) 674 652 ! cloud propertues calculation with the initial values of t and q … … 697 675 ratqs,zqs,t) 698 676 699 700 677 !Jean Jouhaud's version, Decembre 2018 701 !------------------------------------- 702 678 !------------------------------------- 679 703 680 ELSEIF (iflag_cloudth_vert.EQ.6) THEN 704 681 … … 740 717 741 718 ELSE 742 ! When iflag_cld_th=5, we always assume 719 ! When iflag_cld_th=5, we always assume 743 720 ! bi-gaussian distribution 744 721 lognormale = .FALSE. … … 746 723 ENDIF 747 724 748 749 725 DT(:) = 0. 750 726 n_i(:)=0 … … 755 731 ! condensation with qsat(T) variation (adaptation) 756 732 ! Iterative Loop to converge towards qsat 757 758 733 759 734 DO iter=1,iflag_fisrtilp_qsat+1 760 735 761 736 DO i=1,klon 762 737 763 738 ! convergence = .true. until when convergence is not satisfied 764 739 convergence(i)=ABS(DT(i)).GT.DDT0 765 740 766 741 IF ((convergence(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN 767 742 768 743 ! if not convergence: 769 744 … … 781 756 ! Rneb, qzn and zcond for lognormal PDFs 782 757 qtot(i)=zq(i)+zmqc(i) 783 758 784 759 ENDIF 785 760 … … 787 762 788 763 ! Calculation of saturation specific humidity and ce fraction 789 CALL CALC_QSAT_ECMWF(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))790 CALL CALC_GAMMASAT(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:))764 CALL calc_qsat_ecmwf(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:)) 765 CALL calc_gammasat(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:)) 791 766 ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs 792 767 zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:) 793 CALL ICEFRAC_LSCP(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:)) 794 795 768 CALL icefrac_lscp(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:)) 796 769 797 770 DO i=1,klon … … 810 783 zpdf_e2(i)=sign(min(ABS(zpdf_e2(i)),5.),zpdf_e2(i)) 811 784 zpdf_e2(i)=1.-erf(zpdf_e2(i)) 812 785 813 786 !--ice sursaturation by Audran 814 787 IF ((.NOT.ok_ice_sursat).OR.(Tbef(i).GT.t_glace_min)) THEN … … 826 799 fcontrP(i,k)=0.0 !--idem 827 800 qss(i,k)=0.0 !--idem 828 801 829 802 ELSE 830 803 831 804 !------------------------------------ 832 805 ! ICE SUPERSATURATION … … 834 807 835 808 CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, t(i,k), zq(i), & 836 gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k), &837 rneb(i,k), zqn(i), rnebss(i,k), qss(i,k), &809 gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k), & 810 rneb(i,k), zqn(i), rnebss(i,k), qss(i,k), & 838 811 Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) ) 839 812 … … 862 835 qlbef(i)=max(0.,zqn(i)-zqs(i)) 863 836 num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT & 864 &+zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i)837 +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i) 865 838 denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) & 866 839 -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k) & 867 &*qlbef(i)*dzfice(i)840 *qlbef(i)*dzfice(i) 868 841 DT(i)=num/denom 869 842 n_i(i)=n_i(i)+1 … … 874 847 875 848 ENDDO ! iter=1,iflag_fisrtilp_qsat+1 876 877 878 849 879 850 ! P2.3> Final quantities calculation … … 886 857 ! ---------------------------------------------------------------- 887 858 888 889 890 859 ! Water vapor update, Phase determination and subsequent latent heat exchange 891 860 892 861 ! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs) 893 CALL ICEFRAC_LSCP(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:))862 CALL icefrac_lscp(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:)) 894 863 895 864 DO i=1, klon … … 923 892 ENDIF 924 893 925 926 894 ! water vapor update 927 zq(i) = zq(i) - zcond(i) 928 895 zq(i) = zq(i) - zcond(i) 896 929 897 ENDIF 930 898 … … 935 903 ENDDO 936 904 937 938 905 ! If vertical heterogeneity, change volume fraction 939 906 IF (iflag_cloudth_vert .GE. 3) THEN … … 945 912 ! P3> Precipitation formation 946 913 ! ---------------------------------------------------------------- 947 914 948 915 ! LTP: 949 916 IF (iflag_evap_prec==4) THEN … … 972 939 973 940 !2) Clear to cloudy air 974 d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) & 975 - d_tot_zneb(i) - zneb(i))) 941 d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) - d_tot_zneb(i) - zneb(i))) 976 942 IF (znebprecipclr(i) .GT. 0) THEN 977 943 d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i) … … 983 949 984 950 !Update variables 985 znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i) 951 znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i) 986 952 znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i) 987 953 zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i) … … 1021 987 ENDDO 1022 988 1023 CALL FALLICE_VELOCITY(klon,iwc(:),zt(:),zrho(:,k),pplay(:,k),ptconv(:,k),velo(:,k))989 CALL fallice_velocity(klon,iwc(:),zt(:),zrho(:,k),pplay(:,k),ptconv(:,k),velo(:,k)) 1024 990 1025 991 DO i = 1, klon … … 1068 1034 1069 1035 ENDIF 1070 1036 1071 1037 zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain , 0.0) 1072 1038 zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow , 0.0) … … 1076 1042 radliql(i,k) = radliql(i,k) + zoliql(i)/REAL(ninter+1) 1077 1043 radliqi(i,k) = radliqi(i,k) + zoliqi(i)/REAL(ninter+1) 1078 1044 1079 1045 ENDIF ! rneb >0 1080 1046 … … 1083 1049 ENDDO ! n = 1,niter 1084 1050 1085 1086 1087 1051 ! Precipitation flux calculation 1088 1052 1089 1053 DO i = 1, klon 1090 1054 1091 1092 1055 ziflprev(i)=zifl(i) 1093 1056 1094 1057 IF (rneb(i,k) .GT. 0.0) THEN 1095 1096 1058 1097 1059 ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux: … … 1126 1088 ! iced water budget 1127 1089 zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.) 1128 1129 IF (iflag_evap_prec.EQ.4) THEN 1130 zrflcld(i) = zrflcld(i)+zqprecl(i) & 1131 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 1132 ziflcld(i) = ziflcld(i)+ zqpreci(i) & 1133 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 1134 znebprecipcld(i) = rneb(i,k) 1135 zrfl(i) = zrflcld(i) + zrflclr(i) 1136 zifl(i) = ziflcld(i) + ziflclr(i) 1137 ELSE 1090 1091 IF (iflag_evap_prec.EQ.4) THEN 1092 zrflcld(i) = zrflcld(i)+zqprecl(i) & 1093 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 1094 ziflcld(i) = ziflcld(i)+ zqpreci(i) & 1095 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 1096 znebprecipcld(i) = rneb(i,k) 1097 zrfl(i) = zrflcld(i) + zrflclr(i) 1098 zifl(i) = ziflcld(i) + ziflclr(i) 1099 ELSE 1138 1100 zrfl(i) = zrfl(i)+ zqprecl(i) & 1139 1101 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 1140 1102 zifl(i) = zifl(i)+ zqpreci(i) & 1141 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 1103 *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 1142 1104 ENDIF 1143 1105 1144 1106 ENDIF ! rneb>0 1145 1107 1146 1147 1108 ENDDO 1148 1109 1149 1150 1151 1152 ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min 1110 ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min 1153 1111 ! if iflag_evap_pre=4 1154 IF (iflag_evap_prec.EQ.4) THEN 1155 1156 DO i=1,klon 1157 1158 IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN 1112 IF (iflag_evap_prec.EQ.4) THEN 1113 1114 DO i=1,klon 1115 1116 IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN 1159 1117 znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ & 1160 1118 (MAX(znebprecipclr(i),seuil_neb)*rain_int_min), ziflclr(i)/(MAX(znebprecipclr(i),seuil_neb)*rain_int_min))) 1161 1119 ELSE 1162 znebprecipclr(i)=0. 1163 ENDIF 1164 1165 IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN 1120 znebprecipclr(i)=0.0 1121 ENDIF 1122 1123 IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN 1166 1124 znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/ & 1167 1125 (MAX(znebprecipcld(i),seuil_neb)*rain_int_min), ziflcld(i)/(MAX(znebprecipcld(i),seuil_neb)*rain_int_min))) 1168 ELSE 1169 znebprecipcld(i)=0. 1170 ENDIF 1171 1172 ENDDO 1173 1174 ENDIF 1126 ELSE 1127 znebprecipcld(i)=0.0 1128 ENDIF 1129 1130 ENDDO 1131 1132 ENDIF 1175 1133 1176 1134 ! End of precipitation formation … … 1194 1152 ! we recaulate radliqi to account for contributions from the precipitation flux 1195 1153 1196 1197 1154 IF ((ok_radliq_snow) .AND. (k .LT. klev)) THEN 1198 1155 ! for the solid phase (crystals + snowflakes) … … 1205 1162 radliqi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)+ziflprev(i)/zrho(i,k+1)/velo(i,k+1) 1206 1163 ELSE 1207 radliqi(i,k)=zoliq(i)*zfice(i)+zqpreci(i) 1164 radliqi(i,k)=zoliq(i)*zfice(i)+zqpreci(i) 1208 1165 ENDIF 1209 1166 radliq(i,k)=radliql(i,k)+radliqi(i,k) … … 1218 1175 ENDDO 1219 1176 1220 1221 1222 1177 ! ---------------------------------------------------------------- 1223 1178 ! P4> Wet scavenging 1224 1179 ! ---------------------------------------------------------------- 1225 1226 1180 1227 1181 !Scavenging through nucleation in the layer … … 1248 1202 pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi) 1249 1203 frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi 1250 1204 1251 1205 ! Nucleation with a factor of -1 instead of -0.5 1252 1206 zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i)) … … 1254 1208 1255 1209 ENDIF 1256 1257 ENDDO 1258 1259 1210 1211 ENDDO 1212 1260 1213 ! Scavenging through impaction in the underlying layer 1261 1214 … … 1281 1234 1282 1235 ENDDO 1283 1236 1284 1237 !--save some variables for ice supersaturation 1285 1238 ! … … 1294 1247 qclr(i,k) = MAX(1.e-10,zq(i) - qvc(i,k) - qss(i,k)) !--added by OB because of pathologic cases with the lognormal 1295 1248 qcld(i,k) = qvc(i,k) + zcond(i) 1296 END 1249 ENDDO 1297 1250 !q_sat 1298 CALL CALC_QSAT_ECMWF(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,zqsatl(:,k),zdqs(:))1299 CALL CALC_QSAT_ECMWF(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,zqsats(:,k),zdqs(:))1300 1301 END 1251 CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,zqsatl(:,k),zdqs(:)) 1252 CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,zqsats(:,k),zdqs(:)) 1253 1254 ENDDO 1302 1255 1303 1256 !======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
