Changeset 4072 for LMDZ6/trunk/libf/phylmd/lscp_mod.F90
- Timestamp:
- Feb 1, 2022, 6:34:34 PM (2 years ago)
- File:
-
- 1 edited
-
LMDZ6/trunk/libf/phylmd/lscp_mod.F90 (modified) (35 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/lscp_mod.F90
r4062 r4072 6 6 7 7 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 8 SUBROUTINE LSCP(dtime,missing_val, & 8 SUBROUTINE LSCP(dtime,missing_val, & 9 9 paprs,pplay,t,q,ptconv,ratqs, & 10 d_t, d_q, d_ql, d_qi, rneb, rneb_seri, & 10 d_t, d_q, d_ql, d_qi, rneb, rneb_seri, & 11 11 radliq, radicefrac, rain, snow, & 12 12 pfrac_impa, pfrac_nucl, pfrac_1nucl, & … … 25 25 ! 26 26 ! This code is a new version of the fisrtilp.F90 routine, which is itself a 27 ! fusionof 'first' (superrsaturation physics, P. LeVan K. Laval)27 ! merge of 'first' (superrsaturation physics, P. LeVan K. Laval) 28 28 ! and 'ilp' (il pleut, L. Li) 29 29 ! … … 126 126 INTEGER, INTENT(IN) :: iflag_ice_thermo! flag to activate the ice thermodynamics 127 127 ! CR: if iflag_ice_thermo=2, only convection is active 128 LOGICAL, INTENT(IN) :: ok_ice_sursat ! flag to determine if ice sursaturation is activated 128 LOGICAL, INTENT(IN) :: ok_ice_sursat ! flag to determine if ice sursaturation is activated 129 129 130 130 LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active … … 138 138 REAL, DIMENSION(klon,klev), INTENT(IN) :: zpspsk ! exner potential (p/100000)**(R/cp) 139 139 REAL, DIMENSION(klon,klev), INTENT(IN) :: ztla ! liquid temperature within thermals [K] 140 REAL, DIMENSION(klon,klev), INTENT(INOUT) :: zthl 140 REAL, DIMENSION(klon,klev), INTENT(INOUT) :: zthl ! liquid potential temperature [K] 141 141 142 142 ! INPUT/OUTPUT variables … … 145 145 REAL, DIMENSION(klon,klev), INTENT(INOUT):: ratqs ! function of pressure that sets the large-scale 146 146 ! cloud PDF (sigma=ratqs*qt) 147 147 148 148 149 ! Input sursaturation en glace … … 204 205 PARAMETER (ztfondue=278.15) 205 206 206 REAL, SAVE :: rain_int_min=0.001 ! Minimum local rain intensity [mm/s] before the decrease in associatesprecipitation fraction207 REAL, SAVE :: rain_int_min=0.001 ! Minimum local rain intensity [mm/s] before the decrease in associated precipitation fraction 207 208 !$OMP THREADPRIVATE(rain_int_min) 208 209 209 210 210 LOGICAL, SAVE :: ok_radliq_precip=.false. ! Inclusion of all precipitation in radliq (cloud water seen by radiation) 211 !$OMP THREADPRIVATE(ok_radliq_precip) 212 213 214 215 211 216 ! LOCAL VARIABLES: 212 217 !---------------- 213 218 214 219 215 REAL qsl , qsi220 REAL qsl(klon), qsi(klon) 216 221 REAL zct, zcl 217 222 REAL ctot(klon,klev) 218 223 REAL ctot_vol(klon,klev) 219 INTEGER mpc_bl_points(klon,klev)220 224 REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 221 225 REAL zdqsdT_raw(klon) … … 230 234 REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti 231 235 REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon) 232 REAL zoliq p(klon), zoliqi(klon)236 REAL zoliql(klon), zoliqi(klon) 233 237 REAL zt(klon) 234 REAL zdz(klon),zrho(klon), ztot, zrhol(klon)238 REAL zdz(klon),zrho(klon),iwc(klon) 235 239 REAL zchau,zfroi,zfice(klon),zneb(klon),znebprecip(klon) 236 REAL zmelt,z pluie,zice240 REAL zmelt,zrain,zsnow,zprecip 237 241 REAL dzfice(klon) 238 REAL zsolid,qtot,dqsl,dqsi 242 REAL zsolid 243 REAL qtot(klon), qzero(klon) 244 REAL dqsl(klon),dqsi(klon) 239 245 REAL smallestreal 240 246 ! Variables for Bergeron process … … 262 268 REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon) 263 269 REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon) 264 REAL velo(klon )270 REAL velo(klon,klev), vr 265 271 REAL qlmpc, qimpc, rnebmpc 266 REAL radliqi(klon,klev) 272 REAL radliqi(klon,klev), radliql(klon,klev) 267 273 268 274 INTEGER i, k, n, kk 269 275 INTEGER n_i(klon), iter 270 276 INTEGER ncoreczq 277 INTEGER mpc_bl_points(klon,klev) 271 278 INTEGER,SAVE :: itap=0 272 279 !$OMP THREADPRIVATE(itap) … … 334 341 CALL getin_p('seuil_neb',seuil_neb) 335 342 CALL getin_p('rain_int_min',rain_int_min) 343 CALL getin_p('ok_radliq_precip',ok_radliq_precip) 336 344 WRITE(lunout,*) 'lscp, ninter:', ninter 337 345 WRITE(lunout,*) 'lscp, iflag_evap_prec:', iflag_evap_prec 338 346 WRITE(lunout,*) 'lscp, seuil_neb:', seuil_neb 339 347 WRITE(lunout,*) 'lscp, rain_int_min:', rain_int_min 348 WRITE(lunout,*) 'lscp, ok_radliq_precip:', ok_radliq_precip 340 349 341 350 ! check for precipitation sub-time steps … … 397 406 tot_znebn(:) = 0.0 398 407 d_tot_zneb(:) = 0.0 408 qzero(:)=0.0 399 409 400 410 !--ice sursaturation … … 468 478 ! -------------------------------------------------------------------- 469 479 ! A part of the precipitation coming from above is evaporated/sublimated. 470 !471 480 ! -------------------------------------------------------------------- 472 481 … … 474 483 IF (iflag_evap_prec.GE.1) THEN 475 484 485 ! Calculation of saturation specific humidity 486 ! depending on temperature: 487 CALL CALC_QSAT_ECMWF(klon,zt(:),qzero(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:)) 488 ! wrt liquid water 489 CALL CALC_QSAT_ECMWF(klon,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:)) 490 ! wrt ice 491 CALL CALC_QSAT_ECMWF(klon,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:)) 476 492 477 493 DO i = 1, klon … … 480 496 IF (zrfl(i)+zifl(i).GT.0.) THEN 481 497 482 CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))483 484 498 ! LTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec=4). 485 499 IF (iflag_evap_prec.EQ.4) THEN … … 503 517 ENDIF 504 518 505 ! A. JAM 506 ! We consider separately qsatice and qstal 507 ! qsat wrt liquid phase according to thermcep 508 CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,1,.false.,qsl,dqsl) 509 510 519 511 520 ! Evaporation of liquid precipitation coming from above 512 521 ! dP/dz=beta*(1-q/qsat)*sqrt(P) … … 515 524 516 525 IF (iflag_evap_prec.EQ.3) THEN 517 zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl ) &526 zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl(i)) & 518 527 *SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) & 519 528 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG 520 529 ELSE IF (iflag_evap_prec.EQ.4) THEN 521 zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl ) &530 zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl(i)) & 522 531 *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) & 523 532 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG 524 533 ELSE 525 zqevt = 1.*coef_eva*(1.0-zq(i)/qsl )*SQRT(zrfl(i)) &534 zqevt = 1.*coef_eva*(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) & 526 535 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG 527 536 ENDIF … … 531 540 *RG*dtime/(paprs(i,k)-paprs(i,k+1)) 532 541 533 ! qsat wrt ice according to thermcep534 CALL CALC_QSAT_ECMWF(zt(i),0.,pplay(i,k),RTT,2,.false.,qsi,dqsi)535 536 542 ! sublimation of the solid precipitation coming from above 537 543 IF (iflag_evap_prec.EQ.3) THEN 538 zqevti = znebprecip(i)*coef_eva *(1.0-zq(i)/qsi) &544 zqevti = znebprecip(i)*coef_eva_i*(1.0-zq(i)/qsi(i)) & 539 545 *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) & 540 546 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG 541 547 ELSE IF (iflag_evap_prec.EQ.4) THEN 542 zqevti = znebprecipclr(i)*coef_eva *(1.0-zq(i)/qsi) &548 zqevti = znebprecipclr(i)*coef_eva_i*(1.0-zq(i)/qsi(i)) & 543 549 *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) & 544 550 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG 545 551 ELSE 546 zqevti = 1.*coef_eva *(1.0-zq(i)/qsi)*SQRT(zifl(i)) &552 zqevti = 1.*coef_eva_i*(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) & 547 553 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG 548 554 ENDIF … … 566 572 567 573 568 ! EV: agrees with JL's comments below so correct and comment569 ! JLD: I do not understand the lines below. We distribute the liquid570 ! and solid parts of the precipitation even if the layer does not get571 ! saturated. To my opinion, we should all replace with:572 ! zqev=zqevt573 ! zqevi=zqevti574 ! IF (zqevt+zqevti.GT.0.) THEN575 ! zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)576 ! zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)577 ! ELSE578 ! zqev=0.579 ! zqevi=0.580 ! ENDIF581 ! ENDIF582 583 574 ! New solid and liquid precipitation fluxes 584 575 zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) & … … 659 650 ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter 660 651 !------------------------------------------------------- 661 662 DO i = 1, klon 652 653 qtot(:)=zq(:)+zmqc(:) 654 CALL CALC_QSAT_ECMWF(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:)) 655 DO i = 1, klon 663 656 zdelta = MAX(0.,SIGN(1.,RTT-zt(i))) 664 qtot=zq(i)+zmqc(i)665 CALL CALC_QSAT_ECMWF(zt(i),qtot,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i))666 657 zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta) 667 668 658 IF (zq(i) .LT. 1.e-15) THEN 669 659 ncoreczq=ncoreczq+1 … … 671 661 ENDIF 672 662 673 663 ENDDO 674 664 675 665 … … 777 767 DO i=1,klon 778 768 779 ! convergence = .true. sinceconvergence is not satisfied769 ! convergence = .true. until when convergence is not satisfied 780 770 convergence(i)=ABS(DT(i)).GT.DDT0 781 771 … … 796 786 797 787 ! Rneb, qzn and zcond for lognormal PDFs 798 qtot=zq(i)+zmqc(i) 799 CALL CALC_QSAT_ECMWF(Tbef(i),qtot,pplay(i,k),RTT,0,.false.,zqs(i),zdqs(i)) 800 CALL CALC_GAMMASAT(Tbef(i),qtot,pplay(i,k),gammasat(i),dgammasatdt(i)) 801 802 ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs 803 zdqs(i) = gammasat(i)*zdqs(i)+zqs(i)*dgammasatdt(i) 788 qtot(i)=zq(i)+zmqc(i) 789 790 ENDIF 791 792 ENDDO 793 794 ! Calculation of saturation specific humidity and ce fraction 795 CALL CALC_QSAT_ECMWF(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:)) 796 CALL CALC_GAMMASAT(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:)) 797 ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs 798 zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:) 799 CALL ICEFRAC_LSCP(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:)) 800 801 802 803 DO i=1,klon 804 805 IF ((convergence(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN 804 806 805 807 zpdf_sig(i)=ratqs(i,k)*zq(i) … … 826 828 ENDIF 827 829 828 rnebss(i,k)=0.0 !--a jout OB (necessaire car boucle de convergence sur le temps)830 rnebss(i,k)=0.0 !--added by OB (needed because of the convergence loop on time) 829 831 fcontrN(i,k)=0.0 !--idem 830 832 fcontrP(i,k)=0.0 !--idem 831 833 qss(i,k)=0.0 !--idem 832 834 833 835 ELSE 836 834 837 !------------------------------------ 835 ! SURSATURATION EN GLACE838 ! ICE SUPERSATURATION 836 839 !------------------------------------ 837 840 838 841 CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, t(i,k), zq(i), & 839 gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k), & 840 rneb(i,k), zqn(i), rnebss(i,k), qss(i,k), & 842 gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k), & 843 rneb(i,k), zqn(i), rnebss(i,k), qss(i,k), & 841 844 Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) ) 842 845 … … 850 853 851 854 852 ! P2.2.2> Approximative calculation of temperature variation DT 853 ! due to condensation. 854 ! Calculated variables: 855 ! dT : temperature change due to condensation 856 ! --------------------------------------------------------------- 857 858 ! EV: calculation of icefrac in one sole function 859 CALL icefrac_lscp(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:)) 855 ! P2.2.2> Approximative calculation of temperature variation DT 856 ! due to condensation. 857 ! Calculated variables: 858 ! dT : temperature change due to condensation 859 !--------------------------------------------------------------- 860 860 861 861 862 IF (zfice(i).LT.1) THEN … … 896 897 897 898 ! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs) 898 CALL icefrac_lscp(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:))899 CALL ICEFRAC_LSCP(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:)) 899 900 900 901 DO i=1, klon … … 1002 1003 1003 1004 DO i = 1, klon 1004 IF (rneb(i,k).GT.0.0) THEN1005 1005 zoliq(i) = zcond(i) 1006 zrho(i) = pplay(i,k) / zt(i) / RD 1007 zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG) 1008 1009 zneb(i) = MAX(rneb(i,k), seuil_neb) 1006 zoliqi(i)= zoliq(i)*zfice(i) 1007 zoliql(i)= zoliq(i)*(1.-zfice(i)) 1008 zrho(i) = pplay(i,k) / zt(i) / RD 1009 zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG) 1010 iwc(i) = 0. 1011 zneb(i) = MAX(rneb(i,k),seuil_neb) 1010 1012 radliq(i,k) = zoliq(i)/REAL(ninter+1) 1011 1013 radicefrac(i,k) = zfice(i) 1012 1014 radliqi(i,k)=zoliq(i)*zfice(i)/REAL(ninter+1) 1013 ENDIF1015 radliql(i,k)=zoliq(i)*(1.-zfice(i))/REAL(ninter+1) 1014 1016 ENDDO 1015 1017 … … 1024 1026 DO i=1,klon 1025 1027 IF (rneb(i,k).GT.0.0) THEN 1026 zrhol(i) = zrho(i) * zoliq(i) / zneb(i)1028 iwc(i) = zrho(i) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content 1027 1029 ENDIF 1028 1030 ENDDO 1029 1031 1030 CALL FALLICE_VELOCITY(klon, zrhol(:),zt(:),zrho(:),pplay(:,k),ptconv(:,k),velo(:))1032 CALL FALLICE_VELOCITY(klon,iwc(:),zt(:),zrho(:),pplay(:,k),ptconv(:,k),velo(:,k)) 1031 1033 1032 1034 DO i = 1, klon … … 1034 1036 IF (rneb(i,k).GT.0.0) THEN 1035 1037 1036 ! Initialization of z pluie, zice and ztot:1037 z pluie=0.1038 z ice=0.1039 z tot=0.1038 ! Initialization of zrain, zsnow and zprecip: 1039 zrain=0. 1040 zsnow=0. 1041 zprecip=0. 1040 1042 1041 1043 IF (zneb(i).EQ.seuil_neb) THEN 1042 z tot= 0.01043 z ice= 0.01044 z pluie= 0.01044 zprecip = 0.0 1045 zsnow = 0.0 1046 zrain= 0.0 1045 1047 ELSE 1046 1048 ! water quantity to remove: zchau (Sundqvist, 1978) … … 1049 1051 zcl=cld_lc_con 1050 1052 zct=1./cld_tau_con 1051 zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i)*velo(i )*zfice(i)1053 zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i)*velo(i,k)*zfice(i) 1052 1054 ELSE 1053 1055 zcl=cld_lc_lsc 1054 1056 zct=1./cld_tau_lsc 1055 1057 zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) & ! dqice/dt=1/rho*d(rho*wice*qice)/dz 1056 *velo(i ) * zfice(i)1058 *velo(i,k) * zfice(i) 1057 1059 ENDIF 1058 1060 … … 1061 1063 ! surface fraction (which is larger and artificially 1062 1064 ! reduces the in-cloud water). 1065 1063 1066 IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN 1064 1067 zchau = zct *dtime/REAL(ninter) * zoliq(i) & … … 1069 1072 ENDIF 1070 1073 1071 z pluie= MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))1072 z ice= MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))1073 z tot = MAX(zpluie + zice,0.0)1074 zrain = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i))) 1075 zsnow = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i)) 1076 zprecip = MAX(zrain + zsnow,0.0) 1074 1077 1075 1078 ENDIF 1076 1079 1077 ztot = MIN(ztot,zoliq(i)) 1078 zice = MIN(zice,ztot) 1079 zpluie = MIN(zpluie,ztot) 1080 1081 zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie , 0.0) 1082 zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice , 0.0) 1083 zoliq(i) = MAX(zoliq(i)-ztot , 0.0) 1080 zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain , 0.0) 1081 zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow , 0.0) 1082 zoliq(i) = MAX(zoliq(i)-zprecip , 0.0) 1084 1083 1085 1084 radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1) 1085 radliql(i,k) = radliql(i,k) + zoliql(i)/REAL(ninter+1) 1086 1086 radliqi(i,k) = radliqi(i,k) + zoliqi(i)/REAL(ninter+1) 1087 1087 1088 ENDIF 1088 ENDIF ! rneb >0 1089 1089 1090 1090 ENDDO ! i = 1,klon … … 1093 1093 1094 1094 1095 ! Caculate the percentage of ice in "radliq" so cloud+precip seen by the radiation scheme 1095 ! Include the contribution to non-evaporated precipitation to radliq 1096 IF ((ok_radliq_precip) .AND. (k .LT. klev)) THEN 1097 ! for rain drops, we assume a fallspeed of 5m/s 1098 vr=5.0 1099 radliql(:,k)=radliql(:,k)+zrfl(:)/vr/zrho(:) 1100 ! for ice crystals, we take the fallspeed from level above 1101 radliqi(:,k)=radliqi(:,k)+zifl(:)/zrho(:)/velo(:,k+1) 1102 radliq(:,k)=radliql(:,k)+radliqi(:,k) 1103 ENDIF 1104 1105 ! caculate the percentage of ice in "radliq" so cloud+precip seen by the radiation scheme 1096 1106 DO i=1,klon 1097 IF (radliq(i,k) .GT. 0 ) THEN1107 IF (radliq(i,k) .GT. 0.) THEN 1098 1108 radicefrac(i,k)=MIN(MAX(radliqi(i,k)/radliq(i,k),0.),1.) 1099 1109 ENDIF … … 1101 1111 1102 1112 1103 ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux: 1104 ! If T<0C, liquid precip are converted into ice, which leads to an increase in 1105 ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly 1106 ! taken into account through a linearization of the equations and by approximating 1107 ! the liquid precipitation process with a "threshold" process. We assume that 1108 ! condensates are not modified during this operation. Liquid precipitation is 1109 ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases. 1110 ! Water vapor increases as well 1111 ! Note that compared to fisrtilp, we always assume iflag_bergeron=2 1112 1113 1113 ! Precipitation flux calculation 1114 1114 1115 1115 DO i = 1, klon 1116 1116 1117 1117 IF (rneb(i,k) .GT. 0.0) THEN 1118 1119 zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i) 1120 zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i)) 1121 zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) 1122 coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i) 1123 ! Computation of DT if all the liquid precip freezes 1124 DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1)) 1125 ! T should not exceed the freezing point 1126 ! that is Delta > RTT-zt(i) 1127 DeltaT = max( min( RTT-zt(i), DeltaT) , 0. ) 1128 zt(i) = zt(i) + DeltaT 1129 ! water vaporization due to temp. increase 1130 Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT 1131 ! we add this vaporized water to the total vapor and we remove it from the precip 1132 zq(i) = zq(i) + Deltaq 1133 ! The three "max" lines herebelow protect from rounding errors 1134 zcond(i) = max( zcond(i)- Deltaq, 0. ) 1135 ! liquid precipitation freezes oe evaporates 1136 Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT 1137 zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. ) 1138 ! iced water budget 1139 zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.) 1118 1119 1120 ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux: 1121 ! If T<0C, liquid precip are converted into ice, which leads to an increase in 1122 ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly 1123 ! taken into account through a linearization of the equations and by approximating 1124 ! the liquid precipitation process with a "threshold" process. We assume that 1125 ! condensates are not modified during this operation. Liquid precipitation is 1126 ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases. 1127 ! Water vapor increases as well 1128 ! Note that compared to fisrtilp, we always assume iflag_bergeron=2 1129 1130 zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i) 1131 zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i)) 1132 zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) 1133 coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i) 1134 ! Computation of DT if all the liquid precip freezes 1135 DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1)) 1136 ! T should not exceed the freezing point 1137 ! that is Delta > RTT-zt(i) 1138 DeltaT = max( min( RTT-zt(i), DeltaT) , 0. ) 1139 zt(i) = zt(i) + DeltaT 1140 ! water vaporization due to temp. increase 1141 Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT 1142 ! we add this vaporized water to the total vapor and we remove it from the precip 1143 zq(i) = zq(i) + Deltaq 1144 ! The three "max" lines herebelow protect from rounding errors 1145 zcond(i) = max( zcond(i)- Deltaq, 0. ) 1146 ! liquid precipitation converted to ice precip 1147 Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT 1148 zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. ) 1149 ! iced water budget 1150 zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.) 1151 1140 1152 1141 d_ql(i,k) = (1-zfice(i))*zoliq(i)1142 d_qi(i,k) = zfice(i)*zoliq(i)1143 1144 1153 IF (iflag_evap_prec.EQ.4) THEN 1145 1154 zrflcld(i) = zrflcld(i)+zqprecl(i) & … … 1165 1174 IF (iflag_evap_prec.EQ.4) THEN 1166 1175 1167 DO i=1, 1176 DO i=1,klon 1168 1177 1169 1178 IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN … … 1192 1201 ! Outputs: 1193 1202 ! Precipitation fluxes at layer interfaces 1194 ! and temperature and vaportendencies1203 ! and temperature and water species tendencies 1195 1204 DO i = 1, klon 1196 1205 psfl(i,k)=zifl(i) 1197 1206 prfl(i,k)=zrfl(i) 1207 d_ql(i,k) = (1-zfice(i))*zoliq(i) 1208 d_qi(i,k) = zfice(i)*zoliq(i) 1198 1209 d_q(i,k) = zq(i) - q(i,k) 1199 1210 d_t(i,k) = zt(i) - t(i,k) … … 1264 1275 ENDDO 1265 1276 1266 !--save some variables for ice su rsaturation1277 !--save some variables for ice supersaturation 1267 1278 ! 1268 1279 DO i = 1, klon 1269 ! pour la mémoire1280 ! for memory 1270 1281 rneb_seri(i,k) = rneb(i,k) 1271 1282 1272 ! pour lesdiagnostics1283 ! for diagnostics 1273 1284 rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k) 1274 1285 1275 1286 qvc(i,k) = zqs(i) * rneb(i,k) 1276 qclr(i,k) = MAX(1.e-10,zq(i) - qvc(i,k) - qss(i,k)) !--a jout OB a cause de cas pathologiques avec lognormale=F1287 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 1277 1288 qcld(i,k) = qvc(i,k) + zcond(i) 1278 1279 !q_sat 1280 CALL CALC_QSAT_ECMWF(Tbef(i),0.,pplay(i,k),RTT,1,.false.,zqsatl(i,k),zdqs(i)) 1281 CALL CALC_QSAT_ECMWF(Tbef(i),0.,pplay(i,k),RTT,2,.false.,zqsats(i,k),zdqs(i)) 1282 1283 ENDDO 1284 1285 ENDDO 1289 END DO 1290 !q_sat 1291 CALL CALC_QSAT_ECMWF(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,zqsatl(:,k),zdqs(:)) 1292 CALL CALC_QSAT_ECMWF(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,zqsats(:,k),zdqs(:)) 1293 1294 END DO 1286 1295 1287 1296 !======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
