SUBROUTINE mass_redistribution(ngrid,nlayer,nq,ptimestep, & rnat,pcapcal,pplay,pplev,pt,ptsrf,pq,pqs, & pu,pv,pdt,pdtsrf,pdq,pdu,pdv,pdmassmr, & pdtmr,pdtsrfmr,pdpsrfmr,pdumr,pdvmr,pdqmr,pdqsmr) USE watercommon_h, Only: Tsat_water,RLVTT IMPLICIT NONE !======================================================================= ! subject: ! -------- ! Mass and momentum fluxes through sigma levels as the surface pressure is modified are also taken into account ! ! author: Jeremy Leconte 2012 (from F.Forget 1998) ! ------ ! ! input: ! ----- ! ngrid nombre de points de verticales ! (toutes les boucles de la physique sont au ! moins vectorisees sur ngrid) ! nlayer nombre de couches ! pplay(ngrid,nlayer) Pressure levels ! pplev(ngrid,nlayer+1) Pressure levels ! nq Number of tracers ! ! pt(ngrid,nlayer) temperature (en K) ! pq(ngrid,nlayer,nq) tracer specific concentration (kg/kg of air) ! pu,pv (ngrid,nlayer) wind velocity (m/s) ! ! ! pdX physical tendency of X before mass redistribution ! ! pdmassmr air Mass added to the atmosphere in each layer (kg/m2/s) ! ! output: ! ------- ! ! pdXmr(ngrid) physical tendency of X after mass redistribution ! ! ! !======================================================================= ! ! 0. Declarations : ! ------------------ ! #include "dimensions.h" #include "dimphys.h" #include "comcstfi.h" #include "surfdat.h" #include "comgeomfi.h" #include "comvert.h" #include "paramet.h" #include "callkeys.h" #include "tracer.h" !----------------------------------------------------------------------- ! Arguments : ! --------- INTEGER ngrid, nlayer, nq REAL ptimestep REAL pcapcal(ngridmx) INTEGER rnat(ngridmx) REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1) REAL pt(ngrid,nlayer),pdt(ngrid,nlayer) REAL ptsrf(ngrid),pdtsrf(ngrid) REAL pdtmr(ngrid,nlayer) REAL pu(ngrid,nlayer) , pv(ngrid,nlayer) REAL pdu(ngrid,nlayer) , pdv(ngrid,nlayer) REAL pdmassmr(ngrid,nlayer) REAL pdumr(ngrid,nlayer) , pdvmr(ngrid,nlayer) REAL pq(ngridmx,nlayer,nq),pdq(ngrid,nlayer,nq) REAL pqs(ngridmx,nq) REAL pdqmr(ngrid,nlayer,nq),pdqsmr(ngrid,nq) REAL pdpsrfmr(ngrid),pdtsrfmr(ngrid) ! ! Local variables : ! ----------------- ! Boiling/sublimation REAL Tsat(ngridmx),zmassboil(ngridmx) ! vertical reorganization of sigma levels REAL zzu(nlayermx),zzv(nlayermx) REAL zzq(nlayermx,nq),zzt(nlayermx) ! Dummy variables INTEGER n,l,ig,iq REAL zdtsig(ngridmx,nlayermx) REAL zmass(ngridmx,nlayermx),zzmass(nlayermx),w(nlayermx+1) REAL zdmass_sum(ngridmx,nlayermx+1) REAL zmflux(nlayermx+1) REAL zq1(nlayermx) REAL ztsrf(ngridmx) REAL ztm(nlayermx+1) REAL zum(nlayermx+1) , zvm(nlayermx+1) REAL zqm(nlayermx+1,nqmx),zqm1(nlayermx+1) ! local saved variables LOGICAL, SAVE :: firstcall=.true. !---------------------------------------------------------------------- ! Initialisation ! -------------- ! IF (firstcall) THEN firstcall=.false. ENDIF ! !====================================================================== ! Calcul of h2o condensation ! ============================================================ ! ! Used variable : ! pdmassmr : air Mass added to the atmosphere in each layer per unit time (kg/m2/s) ! zdmass_sum(ngrid,l) : total air mass added to the atm above layer l per unit time (kg/m2/s) ! ! ! Surface tracer Tendencies set to 0 ! ------------------------------------- pdqsmr(1:ngridmx,1:nqmx)=0. ztsrf(1:ngridmx) = ptsrf(1:ngridmx) + pdtsrf(1:ngridmx)*ptimestep DO ig=1,ngrid zdmass_sum(ig,nlayermx+1)=0. DO l = nlayermx, 1, -1 zmass(ig,l) = (pplev(ig,l)-pplev(ig,l+1))/g zdmass_sum(ig,l)= zdmass_sum(ig,l+1)+pdmassmr(ig,l) END DO END DO if (water) then do ig=1,ngridmx call Tsat_water(pplev(ig,1)+zdmass_sum(ig,1)*g*ptimestep,Tsat(ig)) enddo call writediagfi(ngrid,'Tsat','saturation temperature at surface','',2,Tsat) do ig=1,ngridmx if (ztsrf(ig).gt.Tsat(ig)) then zmassboil(ig)=(ptsrf(ig)-Tsat(ig))*pcapcal(ig)/RLVTT/ptimestep if ((zmassboil(ig)*ptimestep.gt.pqs(ig,igcm_h2o_vap)).and.(rnat(ig).eq.1)) then zmassboil(ig)=pqs(ig,igcm_h2o_vap)/ptimestep endif zmassboil(ig)=zmassboil(ig)*0.3 !momentary, should be 1. JL12 pdqsmr(ig,igcm_h2o_vap)=-zmassboil(ig) pdtsrfmr(ig)=-zmassboil(ig)*RLVTT/pcapcal(ig) ztsrf(ig)=ptsrf(ig)+pdtsrfmr(ig)*ptimestep else zmassboil(ig)=0. pdtsrfmr(ig)=0. endif enddo endif ! ************************* ! UPDATE SURFACE ! ************************* ! Changing pressure at the surface: ! """""""""""""""""""""""""""""""""""" pdpsrfmr(1:ngridmx) = (zdmass_sum(1:ngridmx,1)+zmassboil(1:ngridmx))*g do ig = 1, ngridmx IF(ABS(pdpsrfmr(ig)*ptimestep).GT.pplev(ig,1)) THEN PRINT*,'STOP in condens' PRINT*,'condensing more than total mass' PRINT*,'Grid point ',ig PRINT*,'Ps = ',pplev(ig,1) PRINT*,'d Ps = ',pdpsrfmr(ig)*ptimestep STOP ENDIF enddo ! of DO ig=1,ngrid ! *************************************************************** ! Correction to account for redistribution between sigma or hybrid ! layers when changing surface pressure ! zzx quantities have dimension (nlayermx) to speed up calculation ! ************************************************************* DO ig=1,ngrid zzt(1:nlayermx) = pt(ig,1:nlayermx) + pdt(ig,1:nlayermx) * ptimestep zzu(1:nlayermx) = pu(ig,1:nlayermx) + pdu(ig,1:nlayermx) * ptimestep zzv(1:nlayermx) = pv(ig,1:nlayermx) + pdv(ig,1:nlayermx) * ptimestep zzq(1:nlayermx,1:nqmx)=pq(ig,1:nlayermx,1:nqmx)+pdq(ig,1:nlayermx,1:nqmx)*ptimestep ! must add the water that has fallen??? ! Mass fluxes of air through the sigma levels (kg.m-2.s-1) (>0 when up) ! """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" zmflux(1) = zmassboil(ig) DO l=1,nlayer zmflux(l+1) = zmflux(l) + pdmassmr(ig,l) - (bp(l)-bp(l+1))*(zdmass_sum(ig,1)+zmflux(1)) ! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) zmflux(l+1)=0. END DO ! Mass of each layer ! ------------------ zzmass(1:nlayermx)=zmass(ig,1:nlayermx)*(1.+pdpsrfmr(ig)*ptimestep/pplev(ig,1)) ! Corresponding fluxes for T,U,V,Q ! """""""""""""""""""""""""""""""" ! averaging operator for TRANSPORT ! """""""""""""""""""""""""""""""" ! Value transfert at the surface interface when condensation ! sublimation: ztm(1) = ztsrf(ig) zum(1) = 0. zvm(1) = 0. zqm(1,1:nqmx)=0. ! most tracer do not condense ! if (water) zqm(1,igcm_h2o_vap)=1. ! flux is 100% h2o at surface ! Van Leer scheme: w(1:nlayermx+1)=-zmflux(1:nlayermx+1)*ptimestep call vl1d(zzt,2.,zzmass,w,ztm) call vl1d(zzu ,2.,zzmass,w,zum) call vl1d(zzv ,2.,zzmass,w,zvm) do iq=1,nqmx zq1(1:nlayermx)=zzq(1:nlayermx,iq) zqm1(1)=zqm(1,iq) ! print*,iq ! print*,zq1 call vl1d(zq1,2.,zzmass,w,zqm1) do l=2,nlayer zzq(l,iq)=zq1(l) zqm(l,iq)=zqm1(l) enddo enddo ! Surface condensation affects low winds if (zmflux(1).lt.0) then zum(1)= zzu(1) * (w(1)/zzmass(1)) zvm(1)= zzv(1) * (w(1)/zzmass(1)) if (w(1).gt.zzmass(1)) then ! ensure numerical stability zum(1)= (zzu(1)-zum(2))*zzmass(1)/w(1) + zum(2) zvm(1)= (zzv(1)-zvm(2))*zzmass(1)/w(1) + zvm(2) end if end if ztm(nlayer+1)= zzt(nlayer) ! should not be used, but... zum(nlayer+1)= zzu(nlayer) ! should not be used, but... zvm(nlayer+1)= zzv(nlayer) ! should not be used, but... zqm(nlayer+1,1:nqmx)= zzq(nlayer,1:nqmx) ! Tendencies on T, U, V, Q ! """""""""""""""""""""""" DO l=1,nlayer ! Tendencies on T pdtmr(ig,l) = (1/zzmass(l)) * & (zmflux(l)*(ztm(l) - zzt(l))-zmflux(l+1)*(ztm(l+1)-zzt(l))) !JL12 the last term in Newcondens has been set to zero because we are only dealing with redistribution here ! Tendencies on U pdumr(ig,l) = (1/zzmass(l)) *( zmflux(l)*(zum(l) - zzu(l)) - zmflux(l+1)*(zum(l+1) - zzu(l)) ) ! Tendencies on V pdvmr(ig,l) = (1/zzmass(l)) *( zmflux(l)*(zvm(l) - zzv(l)) - zmflux(l+1)*(zvm(l+1) - zzv(l)) ) END DO ! Tendencies on Q do iq=1,nqmx DO l=1,nlayer pdqmr(ig,l,iq)= (1/zzmass(l)) * & (zmflux(l)*(zqm(l,iq)-zzq(l,iq))- zmflux(l+1)*(zqm(l+1,iq)-zzq(l,iq)) - pdmassmr(ig,l)*zzq(l,iq)) END DO enddo END DO ! loop on ig return end ! ***************************************************************** SUBROUTINE vl1d(q,pente_max,zzmass,w,qm) ! ! ! Operateur de moyenne inter-couche pour calcul de transport type ! Van-Leer " pseudo amont " dans la verticale ! q,w sont des arguments d'entree pour le s-pg .... ! masse : masse de la couche Dp/g ! w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2) ! pente_max = 2 conseillee ! ! ! -------------------------------------------------------------------- IMPLICIT NONE #include "dimensions.h" ! Arguments: ! ---------- real zzmass(llm),pente_max REAL q(llm),qm(llm+1) REAL w(llm+1) ! ! Local ! --------- ! INTEGER l ! real dzq(llm),dzqw(llm),adzqw(llm),dzqmax real sigw, Mtot, MQtot integer m ! integer ismax,ismin ! On oriente tout dans le sens de la pression ! W > 0 WHEN DOWN !!!!!!!!!!!!! do l=2,llm dzqw(l)=q(l-1)-q(l) adzqw(l)=abs(dzqw(l)) enddo do l=2,llm-1 if(dzqw(l)*dzqw(l+1).gt.0.) then dzq(l)=0.5*(dzqw(l)+dzqw(l+1)) else dzq(l)=0. endif dzqmax=pente_max*min(adzqw(l),adzqw(l+1)) dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l)) enddo dzq(1)=0. dzq(llm)=0. do l = 1,llm-1 ! Regular scheme (transfered mass < layer mass) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if(w(l+1).gt.0. .and. w(l+1).le.zzmass(l+1)) then sigw=w(l+1)/zzmass(l+1) qm(l+1)=(q(l+1)+0.5*(1.-sigw)*dzq(l+1)) else if(w(l+1).le.0. .and. -w(l+1).le.zzmass(l)) then sigw=w(l+1)/zzmass(l) qm(l+1)=(q(l)-0.5*(1.+sigw)*dzq(l)) ! Extended scheme (transfered mass > layer mass) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else if(w(l+1).gt.0.) then m=l+1 Mtot = zzmass(m) MQtot = zzmass(m)*q(m) do while ((m.lt.llm).and.(w(l+1).gt.(Mtot+zzmass(m+1)))) m=m+1 Mtot = Mtot + zzmass(m) MQtot = MQtot + zzmass(m)*q(m) end do if (m.lt.llm) then sigw=(w(l+1)-Mtot)/zzmass(m+1) qm(l+1)= (1/w(l+1))*(MQtot + (w(l+1)-Mtot)*(q(m+1)+0.5*(1.-sigw)*dzq(m+1)) ) else ! w(l+1) = Mtot ! qm(l+1) = Mqtot / Mtot write(*,*) 'top layer is disappearing !',l,Mtot,w(l+1),qm(l+1) print*,zzmass print*,w print*,q print*,qm stop end if else ! if(w(l+1).lt.0) m = l-1 Mtot = zzmass(m+1) MQtot = zzmass(m+1)*q(m+1) if (m.gt.0) then ! because some compilers will have problems ! evaluating zzmass(0) do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+zzmass(m)))) m=m-1 Mtot = Mtot + zzmass(m+1) MQtot = MQtot + zzmass(m+1)*q(m+1) if (m.eq.0) exit end do endif if (m.gt.0) then sigw=(w(l+1)+Mtot)/zzmass(m) qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*(q(m)-0.5*(1.+sigw)*dzq(m)) ) else qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*qm(1)) end if end if enddo ! boundary conditions (not used in newcondens !!) ! qm(llm+1)=0. ! if(w(1).gt.0.) then ! qm(1)=q(1) ! else ! qm(1)=0. ! end if return end