SUBROUTINE newcondens(ngrid,nlayer,nq,ptimestep, $ pcapcal,pplay,pplev,ptsrf,pt, $ pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, $ piceco2,psolaralb,pemisurf, $ pdtc,pdtsrfc,pdpsrf,pduc,pdvc,pdqc, $ fluxsurf_sw,zls) use tracer_mod, only: noms use surfdat_h, only: emissiv, phisfi use geometry_mod, only: latitude ! grid point latitudes (rad) use planete_h USE comcstfi_h #ifndef MESOSCALE USE comvert_mod, ONLY: bp #endif IMPLICIT NONE c======================================================================= c subject: c -------- c Condensation/sublimation of CO2 ice on the ground and in the c atmosphere c (Scheme described in Forget et al., Icarus, 1998) c c author: Francois Forget 1994-1996 c ------ c c input: c ----- c ngrid nombre de points de verticales c (toutes les boucles de la physique sont au c moins vectorisees sur ngrid) c nlayer nombre de couches c pplay(ngrid,nlayer) Pressure levels c pplev(ngrid,nlayer+1) Pressure levels c pt(ngrid,nlayer) temperature (en K) c ptsrf(ngrid) temperature de surface c c \ c pdt(ngrid,nlayer)\ derivee temporelle physique avant condensation c / ou sublimation pour pt,ptsrf c pdtsrf(ngrid) / c c output: c ------- c c pdpsrf(ngrid) \ derivee temporelle physique (contribution de c pdtc(ngrid,nlayer) / la condensation ou sublimation) pour Ps,pt,ptsrf c pdtsrfc(ngrid) / c c Entree/sortie c ------------- c c piceco2(ngrid) : quantite de glace co2 au sol (kg/m2) c psolaralb(ngrid,2) : albedo au sol c pemisurf(ngrid) : emissivite du sol c c======================================================================= c c 0. Declarations : c ------------------ c include "callkeys.h" c----------------------------------------------------------------------- c Arguments : c --------- INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns INTEGER,INTENT(IN) :: nlayer ! number of vertical layers INTEGER,INTENT(IN) :: nq ! number of tracers REAL,INTENT(IN) :: ptimestep ! physics timestep (s) REAL,INTENT(IN) :: pcapcal(ngrid) REAL,INTENT(IN) :: pplay(ngrid,nlayer) !mid-layer pressure (Pa) REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa) REAL,INTENT(IN) :: ptsrf(ngrid) ! surface temperature (K) REAL,INTENT(IN) :: pt(ngrid,nlayer) ! atmospheric temperature (K) REAL,INTENT(IN) :: pphi(ngrid,nlayer) ! geopotential (m2.s-2) REAL,INTENT(IN) :: pdt(ngrid,nlayer) ! tendency on temperature from ! previous physical processes (K/s) REAL,INTENT(IN) :: pdu(ngrid,nlayer) ! tendency on zonal wind (m/s2) ! from previous physical processes REAL,INTENT(IN) :: pdv(ngrid,nlayer) ! tendency on meridional wind (m/s2) ! from previous physical processes REAL,INTENT(IN) :: pdtsrf(ngrid) ! tendency on surface temperature from ! previous physical processes (K/s) REAL,INTENT(IN) :: pu(ngrid,nlayer) ! zonal wind (m/s) REAL,INTENT(IN) :: pv(ngrid,nlayer) ! meridional wind (m/s) REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (../kg_air) REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq) ! tendency on tracers from ! previous physical processes REAL,INTENT(INOUT) :: piceco2(ngrid) ! CO2 ice on the surface (kg.m-2) REAL,INTENT(INOUT) :: psolaralb(ngrid,2) ! albedo of the surface REAL,INTENT(INOUT) :: pemisurf(ngrid) ! emissivity of the surface ! tendencies due to CO2 condensation/sublimation: REAL,INTENT(OUT) :: pdtc(ngrid,nlayer) ! tendency on temperature (K/s) REAL,INTENT(OUT) :: pdtsrfc(ngrid) ! tendency on surface temperature (K/s) REAL,INTENT(OUT) :: pdpsrf(ngrid) ! tendency on surface pressure (Pa/s) REAL,INTENT(OUT) :: pduc(ngrid,nlayer) ! tendency on zonal wind (m.s-2) REAL,INTENT(OUT) :: pdvc(ngrid,nlayer) ! tendency on meridional wind (m.s-2) REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq) ! tendency on tracers ! added to calculate flux dependent albedo: REAL,intent(in) :: fluxsurf_sw(ngrid,2) real,intent(in) :: zls ! solar longitude (rad) c c Local variables : c ----------------- c variables used for albedo parametrization c -------------------------------------------- INTEGER i,j c REAL Fluxmean(jjp1) INTEGER l,ig,iq,icap,nmix LOGICAL transparency, fluxdependent c flag transparency if you want to make the co2ice semi-transparent PARAMETER(transparency=.true.) c flag fluxdependent if you want the co2ice albedo to be dependent on c the incident solar flux PARAMETER(fluxdependent=.false.) REAL slopy,alpha,constA,constB,constT,albediceF_new(ngrid) REAL zt(ngrid,nlayer) REAL zcpi REAL ztcond (ngrid,nlayer+1) ! CO2 condensation temperature (atm) REAL ztcondsol(ngrid) ! CO2 condensation temperature (surface) REAL zdiceco2(ngrid) REAL zcondicea(ngrid,nlayer) REAL zcondices(ngrid) REAL zfallice(ngrid,nlayer+1) , zfallheat REAL zmflux(nlayer+1) REAL zu(nlayer),zv(nlayer) REAL zq(nlayer,nq),zq1(nlayer) REAL ztsrf(ngrid) REAL ztc(nlayer), ztm(nlayer+1) REAL zum(nlayer+1) , zvm(nlayer+1) REAL zqm(nlayer+1,nq),zqm1(nlayer+1) REAL masse(nlayer),w(nlayer+1) REAL Sm(nlayer),Smq(nlayer,nq),mixmas,qmix LOGICAL condsub(ngrid) real :: emisref(ngrid) c variable speciale diagnostique real tconda1(ngrid,nlayer) real tconda2(ngrid,nlayer) c REAL zdiceco2a(ngrid) ! for diagnostic only real zdtsig (ngrid,nlayer) real zdt (ngrid,nlayer) real vmr_co2(ngrid,nlayer) ! co2 volume mixing ratio ! improved_ztcond flag: If set to .true. (AND running with a 'co2' tracer) ! then condensation temperature is computed using partial pressure of CO2 logical,parameter :: improved_ztcond=.true. ! Bound co2 (tracer) values... logical,parameter :: bound_qco2=.false. real,parameter :: qco2max=1.1 real,parameter :: qco2mini=0.1 real :: zqco2 c local saved variables integer,save :: ico2 ! index of CO2 tracer real,save :: qco2min,qco2,mmean real,parameter :: latcond=5.9e5 ! (J/kg) Latent heat of solid CO2 ice real,parameter :: tcond1mb=136.27 ! condensation temperature (K) at 1 mbar real,parameter :: cpice=1000. ! (J.kg-1.K-1) specific heat of CO2 ice REAL,SAVE :: acond,bcond,ccond ! REAL,SAVE :: albediceF(ngrid) real,save :: m_co2, m_noco2, A , B LOGICAL,SAVE :: firstcall = .true. !,firstcall2=.true. integer flag c---------------------------------------------------------------------- c Initialisation c -------------- c IF (firstcall) THEN bcond=1./tcond1mb ccond=cpp/(g*latcond) acond=r/latcond firstcall=.false. write(*,*) 'Newcondens: improved_ztcond=',improved_ztcond write(*,*) 'Newcondens: bound_qco2=',bound_qco2 PRINT*,'In newcondens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond PRINT*,'acond,bcond,ccond',acond,bcond,ccond ico2=0 if (tracer) then c Prepare Special treatment if one of the tracer is CO2 gas do iq=1,nq if (noms(iq).eq."co2") then ico2=iq m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol) m_noco2 = 33.37E-3 ! Non condensible mol mass (kg/mol) c Compute A and B coefficient use to compute c mean molecular mass Mair defined by c 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2 c 1/Mair = A*q(ico2) + B A =(1/m_co2 - 1/m_noco2) B=1/m_noco2 endif enddo c minimum CO2 mix. ratio below which mixing occur with layer above: qco2min =0.75 end if ENDIF ! of IF (firstcall) zcpi=1./cpp c c====================================================================== c Calcul of CO2 condensation sublimation c ============================================================ c c Used variable : c piceco2(ngrid) : amount of co2 ice on the ground (kg/m2) c zcondicea(ngrid,l): condensation rate in layer l (kg/m2/s) c zcondices(ngrid): condensation rate on the ground (kg/m2/s) c zfallice(ngrid,l):amount of ice falling from layer l (kg/m2/s) c c pdtc(ngrid,nlayer) : dT/dt due to cond/sub c c c Tendencies set to 0 (except pdtc) c ------------------------------------- DO l=1,nlayer DO ig=1,ngrid zcondicea(ig,l) = 0. zfallice(ig,l) = 0. pduc(ig,l) = 0 pdvc(ig,l) = 0 END DO END DO DO iq=1,nq DO l=1,nlayer DO ig=1,ngrid pdqc(ig,l,iq) = 0 END DO END DO END DO DO ig=1,ngrid zfallice(ig,nlayer+1) = 0. zcondices(ig) = 0. pdtsrfc(ig) = 0. pdpsrf(ig) = 0. condsub(ig) = .false. zdiceco2(ig) = 0. ENDDO zfallheat=0 c ************************* c ATMOSPHERIC CONDENSATION c ************************* c Compute CO2 Volume mixing ratio c ------------------------------- if (improved_ztcond.and.(ico2.ne.0)) then DO l=1,nlayer DO ig=1,ngrid qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2) mmean=1/(A*qco2 +B) vmr_co2(ig,l) = qco2*mmean/m_co2 ENDDO ENDDO else DO l=1,nlayer DO ig=1,ngrid vmr_co2(ig,l)=0.95 ENDDO ENDDO end if c forecast of atmospheric temperature zt and frost temperature ztcond c -------------------------------------------------------------------- DO l=1,nlayer DO ig=1,ngrid zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep ! ztcond(ig,l)=1./(bcond-acond*log(.0095*pplay(ig,l))) if (pplay(ig,l).ge.1e-4) then ztcond(ig,l)= & 1./(bcond-acond*log(.01*vmr_co2(ig,l)*pplay(ig,l))) else ztcond(ig,l)=0.0 !mars Monica endif ENDDO ENDDO ztcond(:,nlayer+1)=ztcond(:,nlayer) c Condensation/sublimation in the atmosphere c ------------------------------------------ c (calcul of zcondicea , zfallice and pdtc) c DO l=nlayer , 1, -1 DO ig=1,ngrid pdtc(ig,l)=0. IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN condsub(ig)=.true. IF (zfallice(ig,l+1).gt.0) then zfallheat=zfallice(ig,l+1)* & (pphi(ig,l+1)-pphi(ig,l) + & cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/latcond ELSE zfallheat=0. ENDIF pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1)) & *ccond*pdtc(ig,l)- zfallheat c Case when the ice from above sublimes entirely c """"""""""""""""""""""""""""""""""""""""""""""" IF (zfallice(ig,l+1).lt.- zcondicea(ig,l)) then pdtc(ig,l)=(-zfallice(ig,l+1)+zfallheat)/ & (ccond*(pplev(ig,l)-pplev(ig,l+1))) zcondicea(ig,l)= -zfallice(ig,l+1) END IF zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1) END IF ENDDO ENDDO c ************************* c SURFACE CONDENSATION c ************************* c forecast of ground temperature ztsrf and frost temperature ztcondsol c -------------------------------------------------------------------- DO ig=1,ngrid ztcondsol(ig)= & 1./(bcond-acond*log(.01*vmr_co2(ig,1)*pplev(ig,1))) ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep ENDDO c c Condensation/sublimation on the ground c -------------------------------------- c (compute zcondices and pdtsrfc) c DO ig=1,ngrid IF(latitude(ig).lt.0) THEN ! Southern hemisphere icap=2 ELSE ! Northern hemisphere icap=1 ENDIF c c Loop on where we have condensation/ sublimation IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. ! ground cond $ (zfallice(ig,1).NE.0.) .OR. ! falling snow $ ((ztsrf(ig) .GT. ztcondsol(ig)) .AND. ! ground sublim. $ ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN condsub(ig) = .true. IF (zfallice(ig,1).gt.0) then zfallheat=zfallice(ig,1)* & (pphi(ig,1)- phisfi(ig) + & cpice*(ztcond(ig,1)-ztcondsol(ig)))/latcond ELSE zfallheat=0. ENDIF c condensation or partial sublimation of CO2 ice c """"""""""""""""""""""""""""""""""""""""""""""" zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) & /(latcond*ptimestep) - zfallheat pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep c If the entire CO_2 ice layer sublimes c """""""""""""""""""""""""""""""""""""""""""""""""""" c (including what has just condensed in the atmosphere) IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE. & -zcondices(ig))THEN zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1) pdtsrfc(ig)=(latcond/pcapcal(ig))* & (zcondices(ig)+zfallheat) END IF c Changing CO2 ice amount and pressure : c """""""""""""""""""""""""""""""""""" zdiceco2(ig) = zcondices(ig) + zfallice(ig,1) piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep pdpsrf(ig) = -zdiceco2(ig)*g IF(ABS(pdpsrf(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 = ',pdpsrf(ig) STOP ENDIF END IF ! if there is condensation/sublimmation ENDDO ! of DO ig=1,ngrid c ******************************************************************** c Surface albedo and emissivity of the surface below the snow (emisref) c ******************************************************************** c Prepare the case where albedo varies with insolation: c ---------------------------------------------------- ! if (fluxdependent) then ! c Calcul du flux moyen (zonal mean) ! do j=1,jjp1 ! Fluxmean(j)=0 ! do i=1,iim ! ig=1+(j-2)*iim +i ! if(j.eq.1) ig=1 ! if(j.eq.jjp1) ig=ngrid ! Fluxmean(j)=Fluxmean(j)+fluxsurf_sw(ig,1) ! $ +fluxsurf_sw(ig,2) ! enddo ! Fluxmean(j)=Fluxmean(j)/float(iim) ! enddo ! c const A and B used to calculate the albedo which depends on solar flux c albedice=constA+constB*Flux c constT = time step to calculate the solar flux when flux decreases ! constA=0.26 c constA=0.33 c constA=0.186 ! constB=0.00187 ! constT=10 ! endif ! of if (fluxdependent) ! Check that amont of CO2 ice is not problematic DO ig=1,ngrid if(.not.piceco2(ig).ge.0.) THEN if(piceco2(ig).le.-5.e-8) print*, $ 'WARNING newcondens piceco2(',ig,')=', piceco2(ig) piceco2(ig)=0. endif ENDDO ! Set albedo and emissivity of the surface ! ---------------------------------------- CALL albedocaps(zls,ngrid,piceco2,psolaralb,emisref) c Calcul de l'albedo c ------------------ ! do ig =1,ngrid ! IF(ig.GT.ngrid/2+1) THEN ! icap=2 ! ELSE ! icap=1 ! ENDIF ! IF(firstcall2) THEN ! albediceF(ig)=albedice(icap) ! ENDIF c if there is still co2ice ccccccccccccccccccccccc ! if (piceco2(ig).gt.0) then ! emisref(ig) = emisice(icap) c if flux dependent albedo is used c -------------------------------- ! if (fluxdependent) then ! j=INT((ig-2)/iim)+2 ! if(ig.eq.1) j=1 ! if(ig.eq.ngrid) j=jjp1 c albediceF_new(ig)=MIN(constA+constB*Fluxmean(j), c $ constA+constB*250) ! albediceF_new(ig)=constA+constB*Fluxmean(j) ! if (albediceF(ig).gt.albediceF_new(ig)) then ! albediceF(ig)=albediceF(ig)+ ptimestep/(daysec* ! $ constT)*(albediceF_new(ig)-albediceF(ig)) ! else ! albediceF(ig)=albediceF_new(ig) ! endif c if part of the ice is transparent c slopy = pente de la droite: alpha = y*co2ice/1620 c pour une valeur superieur a une epaisseur de glace donnee c ici, epaisseur limite = 10cm ! if (transparency) then ! slopy=1/(1620*0.10) ! alpha=MIN(slopy*piceco2(ig),1.) ! psolaralb(ig,1) = alpha*albediceF(ig) ! $ +(1-alpha)*albedodat(ig) ! psolaralb(ig,2) = psolaralb(ig,1) ! else ! psolaralb(ig,1) = albediceF(ig) ! psolaralb(ig,2) = psolaralb(ig,1) ! endif ! else c transparency set to true and fluxdependent set to false ! if (transparency) then ! slopy=1/(1620*0.10) ! alpha=MIN(slopy*piceco2(ig),1.) ! psolaralb(ig,1) = alpha*albedice(icap) ! $ +(1-alpha)*albedodat(ig) ! psolaralb(ig,2) = psolaralb(ig,1) ! else c simplest case: transparency and flux dependent set to false ! psolaralb(ig,1) = albedice(icap) ! psolaralb(ig,2) = albedice(icap) ! endif ! endif c no more co2ice, albedo = ground albedo ! else ! psolaralb(ig,1) = albedodat(ig) ! psolaralb(ig,2) = albedodat(ig) ! emisref(ig) = emissiv ! pemisurf(ig) = emissiv ! endif ! end do ! end of the ig loop ! set pemisurf() to emissiv when there is bare surface (needed for co2snow) DO ig=1,ngrid if (piceco2(ig).eq.0) then pemisurf(ig)=emissiv endif ENDDO ! firstcall2=.false. c *************************************************************** c Correction to account for redistribution between sigma or hybrid c layers when changing surface pressure (and warming/cooling c of the CO2 currently changing phase). c ************************************************************* DO ig=1,ngrid if (condsub(ig)) then do l=1,nlayer ztc(l) =zt(ig,l) +pdtc(ig,l) *ptimestep zu(l) =pu(ig,l) +pdu( ig,l) *ptimestep zv(l) =pv(ig,l) +pdv( ig,l) *ptimestep do iq=1,nq zq(l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep enddo end do c Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up) c """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" zmflux(1) = -zcondices(ig) DO l=1,nlayer zmflux(l+1) = zmflux(l) -zcondicea(ig,l) #ifndef MESOSCALE & + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1)) c 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. #else if (abs(zmflux(l+1)).lt.1E-13) zmflux(l+1)=0. #endif END DO c Mass of each layer c ------------------ DO l=1,nlayer masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g END DO c Corresponding fluxes for T,U,V,Q c """""""""""""""""""""""""""""""" c averaging operator for TRANSPORT c """""""""""""""""""""""""""""""" c Value transfert at the surface interface when condensation c sublimation: ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep zum(1) = 0 zvm(1) = 0 do iq=1,nq zqm(1,iq)=0. ! most tracer do not condense ! enddo c Special case if one of the tracer is CO2 gas if (ico2.ne.0) zqm(1,ico2)=1. ! flux is 100% CO2 c Van Leer scheme: DO l=1,nlayer+1 w(l)=-zmflux(l)*ptimestep END DO call vl1d(nlayer,ztc,2.,masse,w,ztm) call vl1d(nlayer,zu ,2.,masse,w,zum) call vl1d(nlayer,zv ,2.,masse,w,zvm) do iq=1,nq do l=1,nlayer zq1(l)=zq(l,iq) enddo zqm1(1)=zqm(1,iq) call vl1d(nlayer,zq1,2.,masse,w,zqm1) do l=2,nlayer zq( l,iq)=zq1(l) zqm(l,iq)=zqm1(l) enddo enddo c Surface condensation affects low winds if (zmflux(1).lt.0) then zum(1)= zu(1) * (w(1)/masse(1)) zvm(1)= zv(1) * (w(1)/masse(1)) if (w(1).gt.masse(1)) then ! ensure numerical stability zum(1)= (zu(1)-zum(2))*masse(1)/w(1) + zum(2) zvm(1)= (zv(1)-zvm(2))*masse(1)/w(1) + zvm(2) end if end if ztm(nlayer+1)= ztc(nlayer) ! should not be used, but... zum(nlayer+1)= zu(nlayer) ! should not be used, but... zvm(nlayer+1)= zv(nlayer) ! should not be used, but... do iq=1,nq zqm(nlayer+1,iq)= zq(nlayer,iq) enddo #ifdef MESOSCALE !!!! AS: This part must be commented in the mesoscale model !!!! AS: ... to avoid instabilities. !!!! AS: you have to compile with -DMESOSCALE to do so #else c Tendencies on T, U, V, Q c """""""""""""""""""""""" DO l=1,nlayer c Tendencies on T zdtsig(ig,l) = (1/masse(l)) * & ( zmflux(l)*(ztm(l) - ztc(l)) & - zmflux(l+1)*(ztm(l+1) - ztc(l)) & + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l)) ) pdtc(ig,l) = pdtc(ig,l) + zdtsig(ig,l) c Tendencies on U pduc(ig,l) = (1/masse(l)) * & ( zmflux(l)*(zum(l) - zu(l)) & - zmflux(l+1)*(zum(l+1) - zu(l)) ) c Tendencies on V pdvc(ig,l) = (1/masse(l)) * & ( zmflux(l)*(zvm(l) - zv(l)) & - zmflux(l+1)*(zvm(l+1) - zv(l)) ) END DO #endif c Tendencies on Q do iq=1,nq ! if (noms(iq).eq.'co2') then if (iq.eq.ico2) then c SPECIAL Case when the tracer IS CO2 : DO l=1,nlayer pdqc(ig,l,iq)= (1/masse(l)) * & ( zmflux(l)*(zqm(l,iq) - zq(l,iq)) & - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) & + zcondicea(ig,l)*(zq(l,iq)-1.) ) END DO else DO l=1,nlayer pdqc(ig,l,iq)= (1/masse(l)) * & ( zmflux(l)*(zqm(l,iq) - zq(l,iq)) & - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) & + zcondicea(ig,l)*zq(l,iq) ) END DO end if enddo c -------------------------------------------------------- c Roughly Simulate Molecular mixing when CO2 is too depleted by c Surface condensation (mixing starts if qco2 < qco2min ) c FF 06/2004 c WARNING : this is now done in convadj, better (FF 02/2005) c -------------------------------------------------------- flag=0 ! now done in convadj : must be =0 if (flag.eq.1) then if(ico2.gt.0) then ! relevant only if one tracer is CO2 if(pq(ig,1,ico2)+(pdq(ig,1,ico2)+pdqc(ig,1,ico2))*ptimestep & .lt.qco2min) then do iq=1,nq zq(1,iq)=pq(ig,1,iq) & +(pdq(ig,1,iq)+pdqc(ig,1,iq))*ptimestep Smq(1,iq) = masse(1)*zq(1,iq) end do Sm(1) = masse(1) do l =2,nlayer do iq=1,nq zq(l,iq)=pq(ig,l,iq) & +(pdq(ig,l,iq)+pdqc(ig,l,iq))*ptimestep smq(l,iq) = smq(l-1,iq) + masse(l)*zq(l,iq) end do sm(l) = sm(l-1) + masse(l) if(zq(l,ico2).gt.qco2min) then c mixmas: mass of atmosphere that must be mixed to reach qco2min mixmas = (sm(l-1)*zq(l,ico2)-Smq(l-1,ico2)) & /(zq(l,ico2)-qco2min) if((mixmas.le.sm(l)))then c OK if mixed mass less than mass of layers affected nmix=l ! number of layer affected by mixing goto 99 end if end if end do 99 continue do iq=1,nq qmix=zq(nmix,iq) & +(Smq(nmix-1,iq)-zq(nmix,iq)*Sm(nmix-1))/mixmas do l=1,nmix-1 pdqc(ig,l,iq)= & (qmix-pq(ig,l,iq))/ptimestep - pdq(ig,l,iq) end do c layer only partly mixed : pdqc(ig,nmix,iq)=( & qmix+(Sm(nmix)-mixmas)*(zq(nmix,iq)-qmix)/masse(nmix) & -pq(ig,nmix,iq))/ptimestep -pdq(ig,nmix,iq) end do end if end if endif ! (flag.eq.1) end if ! if (condsub) END DO ! loop on ig c *************************************************************** c CO2 snow / clouds scheme c *************************************************************** call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev, & zcondicea,zcondices,zfallice,pemisurf) c *************************************************************** c Ecriture des diagnostiques c *************************************************************** c DO l=1,nlayer c DO ig=1,ngrid c Taux de cond en kg.m-2.pa-1.s-1 c tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1)) c Taux de cond en kg.m-3.s-1 c tconda2(ig,l)=tconda1(ig,l)*pplay(ig,l)*g/(r*pt(ig,l)) c END DO c END DO c call WRITEDIAGFI(ngrid,'tconda1', c &'Taux de condensation CO2 atmospherique /Pa', c & 'kg.m-2.Pa-1.s-1',3,tconda1) c call WRITEDIAGFI(ngrid,'tconda2', c &'Taux de condensation CO2 atmospherique /m', c & 'kg.m-3.s-1',3,tconda2) ! output falling co2 ice in 1st layer: ! call WRITEDIAGFI(ngrid,'fallice', ! &'Precipitation of co2 ice', ! & 'kg.m-2.s-1',2,zfallice(1,1)) !! Specific stuff to bound co2 tracer values .... if (bound_qco2.and.(ico2.ne.0)) then do ig=1,ngrid do l=1,nlayer zqco2=pq(ig,l,ico2) & +(pdq(ig,l,ico2)+pdqc(ig,l,ico2))*ptimestep if (zqco2.gt.qco2max) then ! correct pdqc: pdqc(ig,l,ico2)=((qco2max-pq(ig,l,ico2))/ptimestep) & -pdq(ig,l,ico2) write(*,*) "newcondens: adapting pdqc(ig,l,ico2)", & " so that co2 conc. does not exceed",qco2max write(*,*) " ig:",ig," l:",l endif ! of if (zqco2.gt.qco2max) if (zqco2.lt.qco2mini) then ! correct pdqc: pdqc(ig,l,ico2)=((qco2mini-pq(ig,l,ico2))/ptimestep) & -pdq(ig,l,ico2) write(*,*) "newcondens: adapting pdqc(ig,l,ico2)", & " so that co2 conc. is not less than",qco2mini write(*,*) " ig:",ig," l:",l endif ! of if (zqco2.lt.qco2mini) end do enddo endif ! of if (bound_qco2.and.(ico2.ne.0)) then #ifndef MESOSCALE ! Extra special case for surface temperature tendency pdtsrfc: ! we want to fix the south pole temperature to CO2 condensation temperature if (caps.and.(obliquit.lt.27.)) then ! check if last grid point is the south pole if (abs(latitude(ngrid)-(-pi/2.)).lt.1.e-5) then ! NB: Updated surface pressure, at grid point 'ngrid', is ! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep ! write(*,*) "newcondens: South pole: latitude(ngrid)=", ! & latitude(ngrid) ztcondsol(ngrid)= & 1./(bcond-acond*log(.01*vmr_co2(ngrid,1)* & (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep))) pdtsrfc(ngrid)=(ztcondsol(ngrid)-ztsrf(ngrid))/ptimestep endif endif #endif return end c ***************************************************************** SUBROUTINE vl1d(nlayer,q,pente_max,masse,w,qm) c c c Operateur de moyenne inter-couche pour calcul de transport type c Van-Leer " pseudo amont " dans la verticale c q,w sont des arguments d'entree pour le s-pg .... c masse : masse de la couche Dp/g c w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2) c pente_max = 2 conseillee c c c -------------------------------------------------------------------- IMPLICIT NONE c c c c Arguments: c ---------- integer nlayer real masse(nlayer),pente_max REAL q(nlayer),qm(nlayer+1) REAL w(nlayer+1) c c Local c --------- c INTEGER l c real dzq(nlayer),dzqw(nlayer),adzqw(nlayer),dzqmax real sigw, Mtot, MQtot integer m c integer ismax,ismin c On oriente tout dans le sens de la pression c W > 0 WHEN DOWN !!!!!!!!!!!!! do l=2,nlayer dzqw(l)=q(l-1)-q(l) adzqw(l)=abs(dzqw(l)) enddo do l=2,nlayer-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(nlayer)=0. do l = 1,nlayer-1 c Regular scheme (transfered mass < layer mass) c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then sigw=w(l+1)/masse(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.masse(l)) then sigw=w(l+1)/masse(l) qm(l+1)=(q(l)-0.5*(1.+sigw)*dzq(l)) c Extended scheme (transfered mass > layer mass) c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else if(w(l+1).gt.0.) then m=l+1 Mtot = masse(m) MQtot = masse(m)*q(m) do while ((m.lt.nlayer).and.(w(l+1).gt.(Mtot+masse(m+1)))) m=m+1 Mtot = Mtot + masse(m) MQtot = MQtot + masse(m)*q(m) end do if (m.lt.nlayer) then sigw=(w(l+1)-Mtot)/masse(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 disapearing !' stop end if else ! if(w(l+1).lt.0) m = l-1 Mtot = masse(m+1) MQtot = masse(m+1)*q(m+1) if (m.gt.0) then ! because some compilers will have problems ! evaluating masse(0) do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m)))) m=m-1 Mtot = Mtot + masse(m+1) MQtot = MQtot + masse(m+1)*q(m+1) if (m.eq.0) exit end do endif if (m.gt.0) then sigw=(w(l+1)+Mtot)/masse(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 c boundary conditions (not used in newcondens !!) c qm(nlayer+1)=0. c if(w(1).gt.0.) then c qm(1)=q(1) c else c qm(1)=0. c end if return end