subroutine rain(ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,rneb) use watercommon_h, only: To, RLVTT, RCPD, RCPV, RV, RVTMP2 implicit none !================================================================== ! ! Purpose ! ------- ! Calculates H2O precipitation using simplified microphysics. ! ! Authors ! ------- ! Adapted from the LMDTERRE code by R. Wordsworth (2009) ! Original author Z. X. Li (1993) ! !================================================================== #include "dimensions.h" #include "dimphys.h" #include "tracer.h" #include "comcstfi.h" #include "callkeys.h" ! Pre-arguments (for universal model) real pq(ngridmx,nlayermx,nqmx) ! tracer (kg/kg) real qsurf(ngridmx,nqmx) ! tracer at the surface (kg.m-2) REAL pdt(ngridmx,nlayermx),pdq(ngridmx,nlayermx,nqmx) real dqrain(ngridmx,nlayermx,nqmx) ! tendency of H2O precipitation (kg/kg.s-1) real dqsrain(ngridmx) ! rain flux at the surface (kg.m-2.s-1) real dqssnow(ngridmx) ! snow flux at the surface (kg.m-2.s-1) REAL d_t(ngridmx,nlayermx) ! temperature increment ! Arguments REAL ptimestep ! time interval REAL pplev(ngridmx,nlayermx+1) ! inter-layer pressure REAL pplay(ngridmx,nlayermx) ! mid-layer pressure REAL t(ngridmx,nlayermx) ! input temperature (K) REAL zt(ngridmx,nlayermx) ! working temperature (K) REAL ql(ngridmx,nlayermx) ! liquid water (Kg/Kg) REAL q(ngridmx,nlayermx) ! specific humidity (Kg/Kg) REAL rneb(ngridmx,nlayermx) ! cloud fraction REAL d_q(ngridmx,nlayermx) ! water vapor increment REAL d_ql(ngridmx,nlayermx) ! liquid water / ice increment ! Subroutine options REAL seuil_neb ! Nebulosity threshold PARAMETER (seuil_neb=0.001) ! REAL ct ! Inverse of cloud precipitation time ! PARAMETER (ct=1./1800.) ! PARAMETER (ct=1./1849.479) REAL cl ! Precipitation threshold PARAMETER (cl=2.0e-4) INTEGER ninter PARAMETER (ninter=5) logical simple ! Use very simple Emanuel scheme? parameter(simple=.false.) logical evap_prec ! Does the rain evaporate? parameter(evap_prec=.true.) ! for simple scheme real t_crit PARAMETER (t_crit=218.0) real lconvert ! for precipitation evaporation (old scheme) real eeff1 real eeff2 ! parameter (eeff1=0.95) ! parameter (eeff2=0.98) parameter (eeff1=0.5) parameter (eeff2=0.8) ! Local variables INTEGER i, k, n REAL zqs(ngridmx,nlayermx), zdelta, zcor REAL zrfl(ngridmx), zrfln(ngridmx), zqev, zqevt REAL zoliq(ngridmx) REAL ztglace REAL zdz(ngridmx),zrho(ngridmx),ztot(ngridmx), zrhol(ngridmx) REAL zchau(ngridmx),zfroi(ngridmx),zfrac(ngridmx),zneb(ngridmx) real ttemp, ptemp real tnext(ngridmx,nlayermx) real l2c(ngridmx,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 firstcall SAVE firstcall ! Online functions REAL fallv, zzz ! falling speed of ice crystals fallv (zzz) = 3.29 * ((zzz)**0.16) DATA firstcall /.true./ 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 !print*,ptimestep !print*,1./ct !if(.not.simple)then ! IF (ABS(ptimestep-1./ct).GT.0.001) THEN ! PRINT*, 'Must talk to Laurent Li!!!' ! call abort ! ENDIF !endif firstcall = .false. ENDIF ! GCM -----> subroutine variables DO k = 1, nlayermx DO i = 1, ngridmx 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 ! Determine the cold clouds by their temperature ztglace = To - 15.0 ! Initialise the outputs DO k = 1, nlayermx DO i = 1, ngridmx d_t(i,k) = 0.0 d_q(i,k) = 0.0 d_ql(i,k) = 0.0 ENDDO ENDDO DO i = 1, ngridmx zrfl(i) = 0.0 zrfln(i) = 0.0 ENDDO ! calculate saturation mixing ratio DO k = 1, nlayermx DO i = 1, ngridmx ttemp = zt(i,k) ptemp = pplay(i,k) call watersat(ttemp,ptemp,zqs(i,k)) ENDDO ENDDO ! get column / layer conversion factor DO k = 1, nlayermx DO i = 1, ngridmx !l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/(g*ptimestep) 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, ngridmx IF (zrfl(i) .GT.0.) THEN zqev = MAX (0.0, (zqs(i,k)-q(i,k)))/ptimestep! BC modif here zqevt = 2.0e-5*(1.0-q(i,k)/zqs(i,k)) & *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) !zqev = MAX (0.0, (zqs(i,k)-q(i,k))*eeff1 ) !zqevt = (zrfl(i)/l2c(i,k))*eeff2 !zqev = MIN (zqev, zqevt) !zrfln(i) = zrfl(i) - zqev*l2c(i,k) !zrfln(i) = zrfl(i) - 1.5e-5*(1.0-q(i,k)/zqs(i,k))*sqrt(zrfl(i)) !zrfln(i) = min(zrfln(i),zrfl(i)) ! this is what is actually written in the manual 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) ENDIF ENDDO ENDIF DO i = 1, ngridmx zoliq(i) = 0.0 ENDDO if(simple)then DO i = 1, ngridmx ttemp = zt(i,k) IF (ttemp .ge. To) 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),0.0) zrfl(i) = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep ENDIF ENDIF ENDDO else DO i = 1, ngridmx 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)-ztglace) / (To-ztglace) zfrac(i) = MAX(zfrac(i), 0.0) zfrac(i) = MIN(zfrac(i), 1.0) zneb(i) = MAX(rneb(i,k), seuil_neb) ENDIF ENDDO DO n = 1, ninter DO i = 1, ngridmx IF (rneb(i,k).GT.0.0) THEN zchau(i) = (1./FLOAT(ninter)) * zoliq(i) & * (1.0-EXP(-(zoliq(i)/zneb(i)/cl)**2)) * zfrac(i) ! warning: this may give dodgy results for physics calls .ne. 48 per day... ! this is the ONLY place where zneb, ct and cl are used zrhol(i) = zrho(i) * zoliq(i) / zneb(i) zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) & *fallv(zrhol(i)) * (1.0-zfrac(i)) ! zfroi behaves oddly... ! * 0.1 * (1.0-zfrac(i)) 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 ! Change in cloud density and surface H2O values DO i = 1, ngridmx 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 simple 9999 continue ! Rain or snow on the ground DO i = 1, ngridmx if(zrfl(i).lt.0.0)then print*,'Droplets of negative rain are falling...' call abort endif IF (t(i,1) .LT. To) 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, ngridmx 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