Index: trunk/LMDZ.MARS/libf/phymars/vdifc.F =================================================================== --- trunk/LMDZ.MARS/libf/phymars/vdifc.F (revision 472) +++ trunk/LMDZ.MARS/libf/phymars/vdifc.F (revision 473) @@ -72,4 +72,6 @@ c ------ + REAL ust(ngridmx),tst(ngridmx) + INTEGER ilev,ig,ilay,nlev @@ -112,4 +114,28 @@ REAL kmixmin + +c Mass-variation scheme : +c ~~~~~~~ + + INTEGER j,l + REAL zcondicea(ngrid,nlayermx) + REAL zt(ngridmx,nlayermx),ztcond(ngridmx,nlayermx+1) + REAL betam(ngridmx,nlayermx),dmice(ngridmx,nlayermx) + REAL pdtc(ngrid,nlayermx) + REAL zhs(ngridmx,nlayermx) + REAL cpice,ccond + SAVE ccond,cpice + DATA cpice /1000./ + +c Theta_m formulation for mass-variation scheme : +c ~~~~~~~ + + INTEGER ico2,llnt(ngridmx) + SAVE ico2 + REAL m_co2, m_noco2, A , B + SAVE A, B, m_co2, m_noco2 + REAL vmr_co2(ngridmx,nlayermx) + REAL qco2,mmean + c ** un petit test de coherence c -------------------------- @@ -126,6 +152,27 @@ bcond=1./tcond1mb acond=r/latcond + ccond=cpp/(g*latcond) PRINT*,'In vdifc: Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond - PRINT*,' acond,bcond',acond,bcond + 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,nqmx + 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 + end if firstcall=.false. @@ -149,4 +196,6 @@ DO ig=1,ngrid za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g +! Mass variation scheme: + betam(ig,ilay)=-za(ig,ilay)*latcond/(cpp*ppopsk(ig,ilay)) ENDDO ENDDO @@ -163,5 +212,5 @@ ENDDO DO ig=1,ngrid - zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig)) + zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf(ig)) ENDDO @@ -185,17 +234,51 @@ ENDIF +c ----------------------------------- c Potential Condensation temperature: c ----------------------------------- -c if (callcond) then -c DO ilev=1,nlay -c DO ig=1,ngrid -c zhcond(ig,ilev) = -c & (1./(bcond-acond*log(.0095*pplay(ig,ilev))))/ppopsk(ig,ilev) -c END DO -c END DO -c else - call zerophys(ngrid*nlay,zhcond) -c end if +c Compute CO2 Volume mixing ratio +c ------------------------------- + if (callcond.and.(ico2.ne.0)) then + DO ilev=1,nlay + DO ig=1,ngrid + qco2=pq(ig,ilev,ico2)+pdqfi(ig,ilev,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,ilev) = qco2*mmean/m_co2 + ENDDO + ENDDO + else + DO ilev=1,nlay + DO ig=1,ngrid + vmr_co2(ig,ilev)=0.95 + ENDDO + ENDDO + end if + +c forecast of atmospheric temperature zt and frost temperature ztcond +c -------------------------------------------------------------------- + + DO ilev=1,nlay + DO ig=1,ngrid + ztcond(ig,ilev)= + & 1./(bcond-acond*log(.01*vmr_co2(ig,ilev)*pplay(ig,ilev))) + if (pplay(ig,ilev).lt.1e-4) ztcond(ig,ilev)=0.0 !mars Monica + ENDDO + ENDDO + + ztcond(:,nlay+1)=ztcond(:,nlay) + + if (callcond) then + DO ilev=1,nlay + DO ig=1,ngrid +! zhcond(ig,ilev) = +! & (1./(bcond-acond*log(.0095*pplay(ig,ilev))))/ppopsk(ig,ilev) + zhcond(ig,ilev) = ztcond(ig,ilev)/ppopsk(ig,ilev) + END DO + END DO + else + call zerophys(ngrid*nlay,zhcond) + end if @@ -209,5 +292,5 @@ zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep - zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev)) +! zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev)) ENDDO ENDDO @@ -239,7 +322,13 @@ zcdv(:)=zcdv_true(:)*sqrt(zu2(:)+(0.7*wmax(:))**2) ! subgrid gustiness added by quadratic interpolation with a coeff beta zcdh(:)=zcdh_true(:)*sqrt(zu2(:)+(1.2*wmax(:))**2) ! LES comparisons. This parameter is linked to the thermals settings) +! ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)+(0.7*wmax(:))**2) +! tst(:)=(ph(:,1)-ptsrf(:))*zcdh(:)/ust(:) +! ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)) +! tst(:)=(ph(:,1)-ptsrf(:))*zcdh_true(:)/sqrt(zcdv_true(:)) ELSE zcdv(:)=zcdv_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic momentum conductance zcdh(:)=zcdh_true(:)*sqrt(zu2(:)) ! 1 / bulk aerodynamic heat conductance +! ust(:)=sqrt(zcdv_true(:))*sqrt(zu2(:)) +! tst(:)=(ph(:,1)-ptsrf(:))*zcdh(:)/ust(:) ENDIF @@ -254,11 +343,7 @@ ! & 2,zcdh_true) ! call WRITEDIAGFI(ngridmx,'u*', -! & 'friction velocity','m/s', -! & 2,sqrt(zcdv_true(:))*sqrt(zu2(:)+(0.7*wmax(:))**2)) -! call WRITEDIAGFI(ngridmx,'T*', -! & 'friction temperature','K', -! & 2,zcdh_true(:)*sqrt(zu2(:)+(1.2*wmax(:))**2)* -! & -(ptsrf(:)-ph(:,1)) -! & /(sqrt(zcdv_true(:))*sqrt(zu2(:)+(0.7*wmax(:))**2))) +! & 'friction velocity','m/s',2,ust) +! call WRITEDIAGFI(ngridmx,'T*', +! & 'friction temperature','K',2,tst) ! call WRITEDIAGFI(ngridmx,'rm-1', ! & 'aerodyn momentum conductance','m/s', @@ -405,22 +490,13 @@ c ------------- +c Mass variation scheme: CALL multipl((nlay-1)*ngrid,zkh(1,2),zb0(1,2),zb(1,2)) CALL multipl(ngrid,zcdh,zb0,zb) - DO ig=1,ngrid - z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) - zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay)*z1(ig) - zd(ig,nlay)=zb(ig,nlay)*z1(ig) - ENDDO - - DO ilay=nlay-1,1,-1 - DO ig=1,ngrid - z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ - $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) - zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+ - $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) - zd(ig,ilay)=zb(ig,ilay)*z1(ig) - ENDDO - ENDDO +c on initialise dm c + + pdtc(:,:)=0. + zt(:,:)=0. + dmice(:,:)=0. c ** calcul de (d Planck / dT) a la temperature d'interface @@ -431,4 +507,57 @@ zdplanck(ig)=z4st*pemis(ig)*ptsrf(ig)*ptsrf(ig)*ptsrf(ig) ENDDO + + +! calcul de zc et zd pour la couche top en prenant en compte le terme +! de variation de masse (on fait une boucle pour que ça converge) + +! Identification des points de grilles qui ont besoin de la correction + + llnt(:)=1 + + DO ig=1,ngrid + DO l=1,nlay + if(zh(ig,l) .lt. zhcond(ig,l)) then + llnt(ig)=300 +! 200 and 100 do not go beyond month 9 with normal dissipation + goto 5 + endif + ENDDO +5 continue + ENDDO + + DO ig=1,ngrid + +! Initialization of z1 and zd, which do not depend on dmice + + z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) + zd(ig,nlay)=zb(ig,nlay)*z1(ig) + + DO ilay=nlay-1,1,-1 + z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ + $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) + zd(ig,ilay)=zb(ig,ilay)*z1(ig) + ENDDO + +! Convergence loop + + DO j=1,llnt(ig) + + z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) + zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay) + & -betam(ig,nlay)*dmice(ig,nlay) + zc(ig,nlay)=zc(ig,nlay)*z1(ig) +! zd(ig,nlay)=zb(ig,nlay)*z1(ig) + +! calcul de zc et zd pour les couches du haut vers le bas + + DO ilay=nlay-1,1,-1 + z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ + $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) + zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+ + $ zb(ig,ilay+1)*zc(ig,ilay+1)- + $ betam(ig,ilay)*dmice(ig,ilay))*z1(ig) +! zd(ig,ilay)=zb(ig,ilay)*z1(ig) + ENDDO c ** calcul de la temperature_d'interface et de sa tendance. @@ -444,30 +573,43 @@ c ---------------------------------------------------- - DO ig=1,ngrid z1(ig)=pcapcal(ig)*ptsrf(ig)+cpp*zb(ig,1)*zc(ig,1) s +zdplanck(ig)*ptsrf(ig)+ pfluxsrf(ig)*ptimestep z2(ig)= pcapcal(ig)+cpp*zb(ig,1)*(1.-zd(ig,1))+zdplanck(ig) ztsrf2(ig)=z1(ig)/z2(ig) - pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep - -c Modif speciale CO2 condensation: -c tconds = 1./(bcond-acond*log(.0095*pplev(ig,1))) -c if ((callcond).and. -c & ((co2ice(ig).ne.0).or.(ztsrf2(ig).lt.tconds)))then -c zh(ig,1)=zc(ig,1)+zd(ig,1)*tconds -c else - zh(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig) -c end if - ENDDO +! pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep !incremented outside loop + + zhs(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig) c ** et a partir de la temperature au sol on remonte c ----------------------------------------------- - DO ilay=2,nlay - DO ig=1,ngrid - hh = max( zh(ig,ilay-1) , zhcond(ig,ilay-1) ) ! modif co2cond - zh(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*hh - ENDDO - ENDDO + DO ilay=2,nlay + zhs(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zhs(ig,ilay-1) + ENDDO + DO ilay=1,nlay + zt(ig,ilay)=zhs(ig,ilay)*ppopsk(ig,ilay) + ENDDO + +c Condensation/sublimation in the atmosphere +c ------------------------------------------ +c (computation of zcondicea and dmice) + + zcondicea(ig,:)=0. + pdtc(ig,:)=0. + + DO l=nlay , 1, -1 + IF(zt(ig,l).LT.ztcond(ig,l)) THEN + pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep + zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1)) + & *ccond*pdtc(ig,l) + dmice(ig,l)= dmice(ig,l) + zcondicea(ig,l)*ptimestep + END IF + ENDDO + + ENDDO !of Do j=1,XXX + + ENDDO !of Do ig=1,nlay + + pdtsrf(:)=(ztsrf2(:)-ptsrf(:))/ptimestep @@ -570,5 +712,5 @@ ! special case for water ice tracer: do not include ! h2o ice tracer from surface (which is set when handling - ! h2o vapour case (see further down). + ! h2o vapour case (see further down). DO ig=1,ngrid z1(ig)=1./(za(ig,1)+zb(ig,1)+ @@ -641,11 +783,10 @@ pdvdif(ig,ilev)=( zv(ig,ilev)- $ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep - hh = max(ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep , - $ zhcond(ig,ilev)) ! modif co2cond - pdhdif(ig,ilev)=( zh(ig,ilev)- hh )/ptimestep - ENDDO - ENDDO - - + hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep + $ + (latcond*dmice(ig,ilev)/cpp)/ppopsk(ig,ilev) + pdhdif(ig,ilev)=( zhs(ig,ilev)- hh )/ptimestep + ENDDO + ENDDO + if (tracer) then DO iq = 1, nq @@ -670,18 +811,9 @@ DO ilev=1,nlay WRITE(*,'(4f15.7)') - s ph(ngrid/2+1,ilev),zh(ngrid/2+1,ilev), + s ph(ngrid/2+1,ilev),zhs(ngrid/2+1,ilev), s pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev) ENDDO ENDIF - - if (calltherm .and. outptherm) then - if (ngrid .eq. 1) then - call WRITEDIAGFI(ngrid,'zh','zh inside vdifc', - & 'SI',1,ph(:,:)+pdhfi(:,:)*ptimestep) - call WRITEDIAGFI(ngrid,'zkh','zkh', - & 'SI',1,zkh) - endif - endif RETURN