subroutine rain(ngrid,nq,ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,rneb) ! to use 'getin' use ioipsl_getincom use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds, RLVTT, RCPD, RCPV, RV, RVTMP2,Psat_water,Tsat_water,rhowater use radii_mod, only: h2o_cloudrad USE tracer_h implicit none !================================================================== ! ! Purpose ! ------- ! Calculates H2O precipitation using simplified microphysics. ! ! Authors ! ------- ! Adapted from the LMDTERRE code by R. Wordsworth (2009) ! Added rain vaporization in case of T>Tsat ! Original author Z. X. Li (1993) ! !================================================================== #include "dimensions.h" #include "dimphys.h" #include "comcstfi.h" #include "callkeys.h" ! Arguments integer,intent(in) :: ngrid ! number of atmospherci columns integer,intent(in) :: nq ! number of tracers real,intent(in) :: ptimestep ! time interval real,intent(in) :: pplev(ngrid,nlayermx+1) ! inter-layer pressure (Pa) real,intent(in) :: pplay(ngrid,nlayermx) ! mid-layer pressure (Pa) real,intent(in) :: t(ngrid,nlayermx) ! input temperature (K) real,intent(in) :: pdt(ngrid,nlayermx) ! input tendency on temperature (K/s) real,intent(in) :: pq(ngrid,nlayermx,nq) ! tracers (kg/kg) real,intent(in) :: pdq(ngrid,nlayermx,nq) ! input tendency on tracers real,intent(out) :: d_t(ngrid,nlayermx) ! temperature tendency (K/s) real,intent(out) :: dqrain(ngrid,nlayermx,nq) ! tendency of H2O precipitation (kg/kg.s-1) real,intent(out) :: dqsrain(ngrid) ! rain flux at the surface (kg.m-2.s-1) real,intent(out) :: dqssnow(ngrid) ! snow flux at the surface (kg.m-2.s-1) real,intent(in) :: rneb(ngrid,nlayermx) ! cloud fraction REAL zt(ngrid,nlayermx) ! working temperature (K) REAL ql(ngrid,nlayermx) ! liquid water (Kg/Kg) REAL q(ngrid,nlayermx) ! specific humidity (Kg/Kg) REAL d_q(ngrid,nlayermx) ! water vapor increment REAL d_ql(ngrid,nlayermx) ! liquid water / ice increment ! Subroutine options REAL,PARAMETER :: seuil_neb=0.001 ! Nebulosity threshold INTEGER,save :: precip_scheme ! id number for precipitaion scheme ! for simple scheme (precip_scheme=1) REAL,SAVE :: rainthreshold ! Precipitation threshold in simple scheme ! for sundquist scheme (precip_scheme=2-3) REAL,SAVE :: cloud_sat ! Precipitation threshold in non simple scheme REAL,SAVE :: precip_timescale ! Precipitation timescale ! for Boucher scheme (precip_scheme=4) REAL,SAVE :: Cboucher ! Precipitation constant in Boucher 95 scheme REAL,PARAMETER :: Kboucher=1.19E8 REAL,SAVE :: c1 INTEGER,PARAMETER :: ninter=5 logical,parameter :: evap_prec=.true. ! Does the rain evaporate? ! for simple scheme real,parameter :: t_crit=218.0 real lconvert ! Local variables INTEGER i, k, n REAL zqs(ngrid,nlayermx),Tsat(ngrid,nlayermx), zdelta, zcor REAL zrfl(ngrid), zrfln(ngrid), zqev, zqevt REAL zoliq(ngrid) REAL zdz(ngrid),zrho(ngrid),ztot(ngrid), zrhol(ngrid) REAL zchau(ngrid),zfroi(ngrid),zfrac(ngrid),zneb(ngrid) real reffh2oliq(ngrid,nlayermx),reffh2oice(ngrid,nlayermx) real ttemp, ptemp, psat_tmp real tnext(ngrid,nlayermx) real l2c(ngrid,nlayermx) real dWtot ! Indices of water vapour and water ice tracers INTEGER, SAVE :: i_vap=0 ! water vapour INTEGER, SAVE :: i_ice=0 ! water ice LOGICAL,SAVE :: firstcall=.true. ! Online functions REAL fallv, fall2v, zzz ! falling speed of ice crystals fallv (zzz) = 3.29 * ((zzz)**0.16) fall2v (zzz) =10.6 * ((zzz)**0.31) !for use with radii IF (firstcall) THEN i_vap=igcm_h2o_vap i_ice=igcm_h2o_ice write(*,*) "rain: i_ice=",i_ice write(*,*) " i_vap=",i_vap PRINT*, 'in rain.F, ninter=', ninter PRINT*, 'in rain.F, evap_prec=', evap_prec write(*,*) "Precipitation scheme to use?" precip_scheme=1 ! default value call getin("precip_scheme",precip_scheme) write(*,*) " precip_scheme = ",precip_scheme if (precip_scheme.eq.1) then write(*,*) "rainthreshold in simple scheme?" rainthreshold=0. ! default value call getin("rainthreshold",rainthreshold) write(*,*) " rainthreshold = ",rainthreshold else if (precip_scheme.eq.2.or.precip_scheme.eq.3) then write(*,*) "cloud water saturation level in non simple scheme?" cloud_sat=2.6e-4 ! default value call getin("cloud_sat",cloud_sat) write(*,*) " cloud_sat = ",cloud_sat write(*,*) "precipitation timescale in non simple scheme?" precip_timescale=3600. ! default value call getin("precip_timescale",precip_timescale) write(*,*) " precip_timescale = ",precip_timescale else if (precip_scheme.eq.4) then write(*,*) "multiplicative constant in Boucher 95 precip scheme" Cboucher=1. ! default value call getin("Cboucher",Cboucher) write(*,*) " Cboucher = ",Cboucher c1=1.00*1.097/rhowater*Cboucher*Kboucher endif firstcall = .false. ENDIF ! GCM -----> subroutine variables DO k = 1, nlayermx DO i = 1, ngrid zt(i,k) = t(i,k)+pdt(i,k)*ptimestep ! a big fat bug was here q(i,k) = pq(i,k,i_vap)+pdq(i,k,i_vap)*ptimestep ql(i,k) = pq(i,k,i_ice)+pdq(i,k,i_ice)*ptimestep !q(i,k) = pq(i,k,i_vap)!+pdq(i,k,i_vap) !ql(i,k) = pq(i,k,i_ice)!+pdq(i,k,i_ice) if(q(i,k).lt.0.)then ! if this is not done, we don't conserve water q(i,k)=0. endif if(ql(i,k).lt.0.)then ql(i,k)=0. endif ENDDO ENDDO ! Initialise the outputs DO k = 1, nlayermx DO i = 1, ngrid d_t(i,k) = 0.0 d_q(i,k) = 0.0 d_ql(i,k) = 0.0 ENDDO ENDDO DO i = 1, ngrid zrfl(i) = 0.0 zrfln(i) = 0.0 ENDDO ! calculate saturation mixing ratio DO k = 1, nlayermx DO i = 1, ngrid ttemp = zt(i,k) ptemp = pplay(i,k) ! call watersat(ttemp,ptemp,zqs(i,k)) call Psat_water(ttemp,ptemp,psat_tmp,zqs(i,k)) call Tsat_water(ptemp,Tsat(i,k)) ENDDO ENDDO ! get column / layer conversion factor DO k = 1, nlayermx DO i = 1, ngrid l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/g ENDDO ENDDO ! Vertical loop (from top to bottom) ! We carry the rain with us and calculate that added by warm/cold precipitation ! processes and that subtracted by evaporation at each level. DO 9999 k = nlayermx, 1, -1 IF (evap_prec) THEN ! note no rneb dependence! DO i = 1, ngrid IF (zrfl(i) .GT.0.) THEN if(zt(i,k).gt.Tsat(i,k))then !! treat the case where all liquid water should boil zqev=MIN((zt(i,k)-Tsat(i,k))*RCPD*l2c(i,k)/RLVTT/ptimestep,zrfl(i)) zrfl(i)=MAX(zrfl(i)-zqev,0.) d_q(i,k)=zqev/l2c(i,k)*ptimestep d_t(i,k) = - d_q(i,k) * RLVTT/RCPD else zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*l2c(i,k)/ptimestep !there was a bug here zqevt= 2.0e-5*(1.0-q(i,k)/zqs(i,k)) & !default was 2.e-5 *sqrt(zrfl(i))*l2c(i,k)/pplay(i,k)*zt(i,k)*R ! BC modif here zqevt = MAX (zqevt, 0.0) zqev = MIN (zqev, zqevt) zqev = MAX (zqev, 0.0) zrfln(i)= zrfl(i) - zqev zrfln(i)= max(zrfln(i),0.0) d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep !d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged! zrfl(i) = zrfln(i) end if ENDIF ENDDO ENDIF DO i = 1, ngrid zoliq(i) = 0.0 ENDDO if(precip_scheme.eq.1)then DO i = 1, ngrid ttemp = zt(i,k) IF (ttemp .ge. T_h2O_ice_liq) THEN lconvert=rainthreshold ELSEIF (ttemp .gt. t_crit) THEN lconvert=rainthreshold*(1.- t_crit/ttemp) lconvert=MAX(0.0,lconvert) ELSE lconvert=0. ENDIF IF (ql(i,k).gt.1.e-9) then zneb(i) = MAX(rneb(i,k), seuil_neb) IF ((ql(i,k)/zneb(i)).gt.lconvert)THEN ! precipitate! d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0) zrfl(i) = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep ENDIF ENDIF ENDDO elseif (precip_scheme.ge.2) then DO i = 1, ngrid IF (rneb(i,k).GT.0.0) THEN zoliq(i) = ql(i,k) zrho(i) = pplay(i,k) / ( zt(i,k) * R ) zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g) zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds) zfrac(i) = MAX(zfrac(i), 0.0) zfrac(i) = MIN(zfrac(i), 1.0) zneb(i) = MAX(rneb(i,k), seuil_neb) ENDIF ENDDO !recalculate liquid water particle radii call h2o_cloudrad(ngrid,ql,reffh2oliq,reffh2oice) SELECT CASE(precip_scheme) !precip scheme from Sundquist 78 CASE(2) DO n = 1, ninter DO i = 1, ngrid IF (rneb(i,k).GT.0.0) THEN ! this is the ONLY place where zneb, precip_timescale and cloud_sat are used zchau(i) = (ptimestep/(FLOAT(ninter)*precip_timescale)) * zoliq(i) & * (1.0-EXP(-(zoliq(i)/zneb(i)/cloud_sat)**2)) * zfrac(i) zrhol(i) = zrho(i) * zoliq(i) / zneb(i) zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) & *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly... ztot(i) = zchau(i) + zfroi(i) IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0 ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i)) zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0) ENDIF ENDDO ENDDO !precip scheme modified from Sundquist 78 (in q**3) CASE(3) DO n = 1, ninter DO i = 1, ngrid IF (rneb(i,k).GT.0.0) THEN ! this is the ONLY place where zneb, precip_timescale and cloud_sat are used zchau(i) = (ptimestep/(FLOAT(ninter)*precip_timescale*cloud_sat**2)) * (zoliq(i)/zneb(i))**3 zrhol(i) = zrho(i) * zoliq(i) / zneb(i) zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) & *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly... ztot(i) = zchau(i) + zfroi(i) IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0 ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i)) zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0) ENDIF ENDDO ENDDO !precip scheme modified from Boucher 95 CASE(4) DO n = 1, ninter DO i = 1, ngrid IF (rneb(i,k).GT.0.0) THEN ! this is the ONLY place where zneb and c1 are used zchau(i) = ptimestep/FLOAT(ninter) *c1* zrho(i) & *(zoliq(i)/zneb(i))**2*reffh2oliq(i,k)*zneb(i)* zfrac(i) zrhol(i) = zrho(i) * zoliq(i) / zneb(i) zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) & *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly... ztot(i) = zchau(i) + zfroi(i) IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0 ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i)) zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0) ENDIF ENDDO ENDDO END SELECT ! precip_scheme ! Change in cloud density and surface H2O values DO i = 1, ngrid IF (rneb(i,k).GT.0.0) THEN d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep zrfl(i) = zrfl(i)+ MAX(ql(i,k)-zoliq(i),0.0)*l2c(i,k)/ptimestep ENDIF ENDDO endif ! if precip_scheme=1 9999 continue ! Rain or snow on the ground DO i = 1, ngrid if(zrfl(i).lt.0.0)then print*,'Droplets of negative rain are falling...' call abort endif IF (t(i,1) .LT. T_h2O_ice_liq) THEN dqssnow(i) = zrfl(i) dqsrain(i) = 0.0 ELSE dqssnow(i) = 0.0 dqsrain(i) = zrfl(i) ! liquid water = ice for now ENDIF ENDDO ! now subroutine -----> GCM variables DO k = 1, nlayermx DO i = 1, ngrid if(evap_prec)then dqrain(i,k,i_vap) = d_q(i,k)/ptimestep d_t(i,k) = d_t(i,k)/ptimestep else dqrain(i,k,i_vap) = 0.0 d_t(i,k) = 0.0 endif dqrain(i,k,i_ice) = d_ql(i,k)/ptimestep ENDDO ENDDO RETURN end subroutine rain