Index: trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90 (revision 485) +++ trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90 (revision 486) @@ -267,5 +267,5 @@ do l=1,nlayer zp=700./pplay(ig,l) - aerosol(ig,l,1)=(dusttau/700.)*(pplev(ig,l)-pplev(ig,l+1)) & + aerosol(ig,l,iaer)=(dusttau/700.)*(pplev(ig,l)-pplev(ig,l+1)) & *max( exp(.03*(1.-max(zp,1.))) , 1.E-3 ) & *QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer) Index: trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90 (revision 485) +++ trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90 (revision 486) @@ -211,25 +211,45 @@ ! Effective radius and variance of the aerosols -! CO2 ice: - DO l=1,nlayer - DO ig=1,ngrid - reffrad(ig,l,1) = 1.e-4 - nueffrad(ig,l,1) = 0.1 + do iaer=1,naerkind ! these values will change once the microphysics gets to work ! UNLESS tracer=.false., in which case we should be working with ! a fixed aerosol layer, and be able to define reffrad in a ! .def file. To be improved! - ENDDO - ENDDO - -! H2O ice: - if(naerkind.eq.2)then - DO l=1,nlayer - DO ig=1,ngrid - reffrad(ig,l,naerkind) = 1.e-5 - nueffrad(ig,l,naerkind) = 0.1 - ENDDO - ENDDO - endif + + if(iaer.eq.1)then ! CO2 ice + do l=1,nlayer + do ig=1,ngrid + reffrad(ig,l,iaer) = 1.e-4 + nueffrad(ig,l,iaer) = 0.1 + enddo + enddo + endif + + if(iaer.eq.2)then ! H2O ice + do l=1,nlayer + do ig=1,ngrid + reffrad(ig,l,iaer) = 1.e-5 + nueffrad(ig,l,iaer) = 0.1 + enddo + enddo + endif + + if(iaer.eq.3)then ! dust + do l=1,nlayer + do ig=1,ngrid + reffrad(ig,l,iaer) = 1.e-5 + nueffrad(ig,l,iaer) = 0.1 + enddo + enddo + endif + + if(iaer.gt.3)then + print*,'Error in callcorrk, naerkind is too high.' + print*,'The code still needs generalisation to arbitrary' + print*,'aerosol kinds and number.' + call abort + endif + + enddo print*, "callcorrk: Correlated-k data base folder:",trim(datadir) @@ -248,5 +268,4 @@ Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed - if((igcm_h2o_vap.eq.0) .and. varactive)then print*,'varactive in callcorrk but no h2o_vap tracer.' @@ -268,5 +287,4 @@ enddo enddo - if(kastprof)then @@ -426,5 +444,4 @@ end do - ! longwave DO nw=1,L_NSPECTI @@ -465,5 +482,4 @@ end do end do - ! test / correct for freaky s. s. albedo values @@ -515,5 +531,4 @@ !tauaero(L_LEVELS+1,iaer) = 0. end do - !print*,'Note changed tauaero BCs in callcorrk!' ! Albedo and emissivity @@ -533,5 +548,4 @@ acosz=mu0(ig) ! cosine of sun incident angle endif - !----------------------------------------------------------------------- @@ -591,10 +605,9 @@ end do - - !----------------------------------------------------------------------- ! kcm mode only if(kastprof)then + ! initial values equivalent to mugaz DO l=1,nlayer muvarrad(2*l) = mugaz @@ -602,8 +615,8 @@ END DO - do k=1,L_LEVELS - qvar(k) = 0.0 - end do - print*,'ASSUMING qH2O=0 EVERYWHERE IN CALLCORRK!' + !do k=1,L_LEVELS + ! qvar(k) = 0.0 + !end do + !print*,'ASSUMING qH2O=0 EVERYWHERE IN CALLCORRK!' endif @@ -633,7 +646,5 @@ qvar(2*nlayermx+1)=qsurf(ig,i_var)*muvar(ig,1)/mH2O - - endif - + endif ! Keep values inside limits for which we have radiative transfer coefficients @@ -678,5 +689,4 @@ tmid(L_LEVELS) = tlevrad(L_LEVELS) - ! test for out-of-bounds pressure if(plevrad(3).lt.pgasmin)then @@ -695,12 +705,11 @@ print*,'Minimum temperature is outside the radiative' print*,'transfer kmatrix bounds, exiting.' - print*,'WARNING, OVERRIDING FOR TEST' - !call abort + !print*,'WARNING, OVERRIDING FOR TEST' + call abort elseif(tlevrad(k).gt.tgasmax)then print*,'Maximum temperature is outside the radiative' print*,'transfer kmatrix bounds, exiting.' - print*,'WARNING, OVERRIDING FOR TEST' - !print*, 'T=',pt - !call abort + !print*,'WARNING, OVERRIDING FOR TEST' + call abort endif enddo @@ -735,9 +744,7 @@ tmid,pmid,taugsurf,qvar,muvarrad) - call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, & acosz,stel_fract,gweight, & nfluxtopv,nfluxoutv_nu,nfluxgndv_nu, & - !acosz,stel_fract,gweight,nfluxtopv,nfluxgndv_nu, & fmnetv,fluxupv,fluxdnv,fzerov,taugsurf) @@ -753,7 +760,4 @@ end if - - - !----------------------------------------------------------------------- ! Longwave @@ -846,6 +850,6 @@ if(specOLR)then if(ngrid.ne.1)then - !call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W m^-2",3,OLR_nu) - !call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W m^-2",3,OSR_nu) + call writediagspecIR(ngrid,"OLR3D","OLR(lon,lat,band)","W m^-2",3,OLR_nu) + call writediagspecVI(ngrid,"OSR3D","OSR(lon,lat,band)","W m^-2",3,OSR_nu) endif endif @@ -892,5 +896,5 @@ endif - !!! see physiq.F for explanations about CLFvarying. This is temporary. + ! see physiq.F for explanations about CLFvarying. This is temporary. if (lastcall .and. .not.CLFvarying) then IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi ) Index: trunk/LMDZ.GENERIC/libf/phystd/callkeys.h =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/callkeys.h (revision 485) +++ trunk/LMDZ.GENERIC/libf/phystd/callkeys.h (revision 486) @@ -38,4 +38,5 @@ & , hydrology & & , sourceevol & + & , icetstep & & , albedosnow & & , maxicethick & @@ -101,2 +102,4 @@ real tau_relax real cloudlvl + real icetstep + Index: trunk/LMDZ.GENERIC/libf/phystd/callsedim.F =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/callsedim.F (revision 485) +++ trunk/LMDZ.GENERIC/libf/phystd/callsedim.F (revision 486) @@ -53,5 +53,5 @@ INTEGER l,ig, iq - real zqi(ngridmx,nlayermx) ! to locally store tracers + real zqi(ngridmx,nlayermx,nqmx) ! to locally store tracers real masse (ngridmx,nlayermx) ! Layer mass (kg.m-2) real epaisseur (ngridmx,nlayermx) ! Layer thickness (m) @@ -79,5 +79,5 @@ firstcall=.false. ENDIF ! of IF (firstcall) - + !======================================================================= ! Preliminary calculation of the layer characteristics @@ -91,39 +91,54 @@ end do - iq=igcm_h2o_ice +!====================================================================== +! Calculate the transport due to sedimentation for each tracer +! [This has been rearranged by L. Kerber to allow the sedimentation +! of general tracers.] + + zqi(1:ngrid,1:nlay,1:nqmx) = 0. + do iq=1,nq + if((radius(iq).gt.1.e-9).and.(iq.ne.igcm_co2_ice)) then +! (no sedim for gases; co2_ice sedim is done elsewhere) -! The value of q is updated after the other parameterisations +! store locally updated tracers - do l=1,nlay + do l=1,nlay do ig=1,ngrid - ! store locally updated tracers - zqi(ig,l)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep + zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep enddo enddo ! of do l=1,nlay + +!====================================================================== +! Sedimentation +!====================================================================== +! Water + if (water.and.(iq.eq.igcm_h2o_ice)) then + call newsedim(ngrid,nlay,ngrid*nlay,ptimestep, + & pplev,masse,epaisseur,pt,rfall,rho_q(iq),zqi,wq) - - -!======================================================================= -! Calculate the transport due to sedimentation for each tracer - - call newsedim(ngrid,nlay,ngrid*nlay,ptimestep, - & pplev,masse,epaisseur,pt,rfall,rho_q(iq),zqi,wq) +! General Case + else + call newsedim(ngrid,nlay,1,ptimestep, + & pplev,masse,epaisseur,pt,radius(iq),rho_q(iq), + & zqi,wq) + endif !======================================================================= ! Calculate the tendencies +!====================================================================== do ig=1,ngrid ! Ehouarn: with new way of tracking tracers by name, this is simply - pdqs_sed(ig,iq)=wq(ig,1)/ptimestep + pdqs_sed(ig,iq) = wq(ig,1)/ptimestep end do DO l = 1, nlay DO ig=1,ngrid - pdqsed(ig,l,iq)=(zqi(ig,l)- - $ (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep + pdqsed(ig,l,iq)=(zqi(ig,l,iq)- + & (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep ENDDO ENDDO - + endif ! of no gases no co2_ice + enddo ! of do iq=1,nq return end - Index: trunk/LMDZ.GENERIC/libf/phystd/condens_co2cloud.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/condens_co2cloud.F90 (revision 485) +++ (revision ) @@ -1,535 +1,0 @@ - subroutine condens_co2cloud(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,reffrad,cpp3D) - - use radinc_h, only : naerkind - use gases_h - - implicit none - -!================================================================== -! Purpose -! ------- -! Condense and/or sublime CO2 ice on the ground and in the -! atmosphere, and sediment the ice. -! -! Inputs -! ------ -! ngrid Number of vertical columns -! nlayer Number of layers -! pplay(ngrid,nlayer) Pressure layers -! pplev(ngrid,nlayer+1) Pressure levels -! pt(ngrid,nlayer) Temperature (in K) -! ptsrf(ngrid) Surface temperature -! -! pdt(ngrid,nlayermx) Time derivative before condensation/sublimation of pt -! pdtsrf(ngrid) Time derivative before condensation/sublimation of ptsrf -! pqsurf(ngrid,nq) Sedimentation flux at the surface (kg.m-2.s-1) -! -! Outputs -! ------- -! pdpsrf(ngrid) \ Contribution of condensation/sublimation -! pdtc(ngrid,nlayermx) / to the time derivatives of Ps, pt, and ptsrf -! pdtsrfc(ngrid) / -! -! Both -! ---- -! piceco2(ngrid) CO2 ice at the surface (kg/m2) -! psolaralb(ngrid) Albedo at the surface -! pemisurf(ngrid) Emissivity of the surface -! -! Authors -! ------- -! Francois Forget (1996) -! Converted to Fortran 90 and slightly modified by R. Wordsworth (2009) -! -!================================================================== - -#include "dimensions.h" -#include "dimphys.h" -#include "comcstfi.h" -#include "surfdat.h" -#include "comgeomfi.h" -#include "comvert.h" -#include "callkeys.h" -#include "tracer.h" - -!----------------------------------------------------------------------- -! Arguments - - INTEGER ngrid, nlayer, nq - - REAL ptimestep - REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1) - REAL pcapcal(ngrid) - REAL pt(ngrid,nlayer) - REAL ptsrf(ngrid) - REAL pphi(ngrid,nlayer) - REAL pdt(ngrid,nlayer),pdtsrf(ngrid),pdtc(ngrid,nlayer) - REAL pdtsrfc(ngrid),pdpsrf(ngrid) -! REAL piceco2(ngrid),psolaralb(ngrid,2),pemisurf(ngrid) - REAL piceco2(ngrid),psolaralb(ngrid),pemisurf(ngrid) - - REAL pu(ngrid,nlayer) , pv(ngrid,nlayer) - REAL pdu(ngrid,nlayer) , pdv(ngrid,nlayer) - REAL pduc(ngrid,nlayer) , pdvc(ngrid,nlayer) - REAL pq(ngridmx,nlayer,nq),pdq(ngrid,nlayer,nq) - REAL pdqc(ngrid,nlayer,nq), pdqsc(ngrid,nq) - - REAL reffrad(ngrid,nlayer,naerkind) - -! Allow variations in cp with location - REAL cpp3D(ngridmx,nlayermx) ! specific heat capacity at const. pressure - - - -!----------------------------------------------------------------------- -! Local variables - - INTEGER l,ig,icap,ilay,it,iq - - REAL*8 zt(ngridmx,nlayermx) - REAL zq(ngridmx,nlayermx,nqmx) - REAL zcpi - REAL ztcond (ngridmx,nlayermx) - REAL ztcondsol(ngridmx) - REAL zdiceco2(ngridmx) - REAL zcondicea(ngridmx,nlayermx), zcondices(ngridmx) - REAL zfallice(ngridmx), Mfallice(ngridmx) - REAL zmflux(nlayermx+1) - REAL zu(nlayermx),zv(nlayermx) - REAL ztsrf(ngridmx) - REAL ztc(nlayermx), ztm(nlayermx+1) - REAL zum(nlayermx+1) , zvm(nlayermx+1) - LOGICAL condsub(ngridmx) - REAL subptimestep - Integer Ntime - real masse (ngridmx,nlayermx), w(ngridmx,nlayermx,nqmx) - real wq(ngridmx,nlayermx+1) - real vstokes,reff - -! Special diagnostic variables - real tconda1(ngridmx,nlayermx) - real tconda2(ngridmx,nlayermx) - real zdtsig (ngridmx,nlayermx) - real zdt (ngridmx,nlayermx) - -!----------------------------------------------------------------------- -! Saved local variables - - REAL emisref(ngridmx) - REAL latcond - REAL ccond - REAL cpice - SAVE emisref, cpice - SAVE latcond,ccond - - LOGICAL firstcall - SAVE firstcall - REAL SSUM - EXTERNAL SSUM - - DATA latcond /5.9e5/ - DATA cpice /1000./ - DATA firstcall/.true./ - - REAL CBRT - EXTERNAL CBRT - - INTEGER,SAVE :: i_co2ice=0 ! co2 ice - CHARACTER(LEN=20) :: tracername ! to temporarily store text - - integer igas - integer,save :: igasco2=0 - character(len=3) :: gasname - - real reffradmin, reffradmax - - real ppco2 - -!----------------------------------------------------------------------- -! Initializations - - call zerophys(ngrid*nlayer*nq,pdqc) - call zerophys(ngrid*nlayer*nq,pdtc) - call zerophys(ngridmx*nlayermx*nqmx,zq) - call zerophys(ngridmx*nlayermx,zt) - - !reffradmin=1.e-7 - !reffradmax=5.e-4 - !reffradmin=0.5e-7 - !reffradmax=0.1e-3 ! FF data - reffradmin=0.1e-5 - reffradmax=0.1e-3 ! JB data -! improve this in the future... - -! Initialisation - IF (firstcall) THEN - - ! find CO2 ice tracer - do iq=1,nqmx - tracername=noms(iq) - if (tracername.eq."co2_ice") then - i_co2ice=iq - endif - enddo - - ! find CO2 gas - do igas=1,ngasmx - gasname=gnom(igas) -! gasname=noms(igas) ! was a bug - if (gasname.eq."CO2") then - igasco2=igas - endif - enddo - - write(*,*) "condense_co2cloud: i_co2ice=",i_co2ice - - if((i_co2ice.lt.1))then - print*,'In condens_co2cloud but no CO2 ice tracer, exiting.' - stop - endif - - ccond=cpp/(g*latcond) - print*,'In condens_co2cloud: ccond=',ccond,' latcond=',latcond - -! Prepare special treatment if gas is not pure CO2 - !if (addn2) then - ! m_co2 = 44.01E-3 ! CO2 molecular mass (kg/mol) - ! m_noco2 = 28.02E-3 ! N2 molecular mass (kg/mol) -! Compute A and B coefficient use to compute -! mean molecular mass Mair defined by -! 1/Mair = q(ico2)/m_co2 + (1-q(ico2))/m_noco2 -! 1/Mair = A*q(ico2) + B - ! A = (1/m_co2 - 1/m_noco2) - ! B = 1/m_noco2 - !endif - -! Minimum CO2 mixing ratio below which mixing occurs with layer above: - !qco2min =0.75 - - firstcall=.false. - ENDIF - zcpi=1./cpp - -!----------------------------------------------------------------------- -! Calculation of CO2 condensation / sublimation -! -! Variables used: -! piceco2(ngrid) amount of co2 ice on the ground (kg/m2) -! zcondicea(ngrid,l) condensation rate in layer l (kg/m2/s) -! zcondices(ngrid) condensation rate on the ground (kg/m2/s) -! zfallice(ngrid) flux of ice falling on surface (kg/m2/s) -! pdtc(ngrid,nlayermx) dT/dt due to phase changes (K/s) - - -! Tendencies initially set to 0 (except pdtc) - DO l=1,nlayer - DO ig=1,ngrid - zcondicea(ig,l) = 0. - pduc(ig,l) = 0 - pdvc(ig,l) = 0 - pdqc(ig,l,i_co2ice) = 0 - END DO - END DO - - DO ig=1,ngrid - Mfallice(ig) = 0. - zfallice(ig) = 0. - zcondices(ig) = 0. - pdtsrfc(ig) = 0. - pdpsrf(ig) = 0. - condsub(ig) = .false. - zdiceco2(ig) = 0. - ENDDO - -!----------------------------------------------------------------------- -! Atmospheric condensation - - -! Compute CO2 Volume mixing ratio -! ------------------------------- -! if (addn2) then -! DO l=1,nlayer -! DO ig=1,ngrid -! qco2=pq(ig,l,ico2)+pdq(ig,l,ico2)*ptimestep -!! 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.5 -! ENDDO -! ENDDO -! end if - -! Forecast the atmospheric frost temperature 'ztcond' - DO l=1,nlayer - DO ig=1,ngrid - ppco2=gfrac(igasco2)*pplay(ig,l) - call get_tcond_co2(ppco2,ztcond(ig,l)) - ENDDO - ENDDO - -! Initialize zq and zt at the beginning of the sub-timestep loop - DO l=1,nlayer - DO ig=1,ngrid - zt(ig,l)=pt(ig,l) - zq(ig,l,i_co2ice)=pq(ig,l,i_co2ice) - IF( zq(ig,l,i_co2ice).lt.-1.e-6 ) THEN - print*,'Uh-oh, zq = ',zq(ig,l,i_co2ice),'at ig,l=',ig,l - if(l.eq.1)then - print*,'Perhaps the atmosphere is collapsing on surface...?' - endif - END IF - ENDDO - ENDDO - -! Calculate the mass of each atmospheric layer (kg.m-2) - do ilay=1,nlayer - do ig=1, ngrid - masse(ig,ilay)=(pplev(ig,ilay) - pplev(ig,ilay+1)) /g - end do - end do - -! ----------------------------------------------- -! START CONDENSATION/SEDIMENTATION SUB-TIME LOOP -! ----------------------------------------------- - Ntime = 20 ! number of sub-timestep - subptimestep = ptimestep/float(Ntime) - - DO it=1,Ntime - -! Add the tendencies from other physical processes at each subtimstep - DO l=1,nlayer - DO ig=1,ngrid - zt(ig,l) = zt(ig,l) + pdt(ig,l) * subptimestep - zq(ig,l,i_co2ice) = zq(ig,l,i_co2ice) + pdq(ig,l,i_co2ice) * subptimestep - END DO - END DO - - -! Gravitational sedimentation - -! sedimentation computed from radius computed from q -! assuming that the ice is splitted in Nmix particle - do ilay=1,nlayer - do ig=1, ngrid - - reff = CBRT( 3*zq(ig,ilay,i_co2ice)/( 4*Nmix_co2*pi*rho_co2 )) - - ! there should be a more elegant way of doing this... - if(reff.lt.1.e-16) reff=1.e-16 - - ! update reffrad for radiative transfer - if(reff.lt.reffradmin)then - reffrad(ig,ilay,1)=reffradmin - !print*,'reff below optical limit' - elseif(reff.gt.reffradmax)then - reffrad(ig,ilay,1)=reffradmax - !print*,'reff above optical limit' - else - reffrad(ig,ilay,1)=reff - endif - - call stokes & - (pplev(ig,ilay),pt(ig,ilay), & - reff,vstokes,rho_co2) - - !w(ig,ilay,i_co2ice) = 0.0 - w(ig,ilay,i_co2ice) = vstokes * subptimestep * & - pplev(ig,ilay)/(r*pt(ig,ilay)) - - end do - end do - -! Computing q after sedimentation - - call vlz_fi(ngrid,zq(1,1,i_co2ice),2.,masse,w(1,1,i_co2ice),wq) - - -! Progressively accumulating the flux to the ground -! Mfallice is the total amount of ice fallen to the ground - do ig=1,ngrid - Mfallice(ig) = Mfallice(ig) + wq(ig,i_co2ice) - end do - - -! Condensation / sublimation in the atmosphere -! -------------------------------------------- -! (calculation of zcondicea, zfallice and pdtc) -! (MODIFICATIONS FOR EARLY MARS: falling heat neglected, condensation -! of CO2 into tracer i_co2ice) - - DO l=nlayer , 1, -1 - DO ig=1,ngrid - pdtc(ig,l)=0. - - ! tweak ccond if the gas is non-ideal - if(nonideal)then - ccond=cpp3D(ig,l)/(g*latcond) - endif - - - IF ((zt(ig,l).LT.ztcond(ig,l)).or.(zq(ig,l,i_co2ice).gt.0)) THEN - condsub(ig)=.true. - pdtc(ig,l) = (ztcond(ig,l) - zt(ig,l))/subptimestep - pdqc(ig,l,i_co2ice) = pdtc(ig,l)*ccond*g - -! Case when the ice from above sublimes entirely - IF ((zq(ig,l,i_co2ice).lt.-pdqc(ig,l,i_co2ice)*subptimestep) & - .AND. (zq(ig,l,i_co2ice).gt.0)) THEN - - pdqc(ig,l,i_co2ice) = -zq(ig,l,i_co2ice)/subptimestep - pdtc(ig,l) =-zq(ig,l,i_co2ice)/(ccond*g*subptimestep) - - END IF - -! Temperature and q after condensation - zt(ig,l) = zt(ig,l) + pdtc(ig,l) * subptimestep - zq(ig,l,i_co2ice) = zq(ig,l,i_co2ice) + pdqc(ig,l,i_co2ice) * subptimestep - END IF - - ENDDO - ENDDO - ENDDO ! end of subtimestep loop - -! Computing global tendencies after the subtimestep - DO l=1,nlayer - DO ig=1,ngrid - pdtc(ig,l) = & - (zt(ig,l) - (pt(ig,l) + pdt(ig,l)*ptimestep))/ptimestep - pdqc(ig,l,i_co2ice) = & - (zq(ig,l,i_co2ice)-(pq(ig,l,i_co2ice)+pdq(ig,l,i_co2ice)*ptimestep))/ptimestep - END DO - END DO - DO ig=1,ngrid - zfallice(ig) = Mfallice(ig)/ptimestep - END DO - - -!----------------------------------------------------------------------- -! Condensation/sublimation on the ground -! (calculation of zcondices and pdtsrfc) - -! forecast of ground temperature ztsrf and frost temperature ztcondsol - DO ig=1,ngrid - ppco2=gfrac(igasco2)*pplay(ig,1) - call get_tcond_co2(ppco2,ztcondsol(ig)) - - ztsrf(ig) = ptsrf(ig) - - if((ztsrf(ig).le.ztcondsol(ig)+2.0).and.(ngrid.eq.1))then - print*,'CO2 is condensing on the surface in 1D. This atmosphere is doomed.' - print*,'T_surf = ',ztsrf,'K' - print*,'T_cond = ',ztcondsol,'K' - open(116,file='surf_vals.out') - write(116,*) 0.0, pplev(1,1), 0.0, 0.0 - close(116) - call abort - endif - - ztsrf(ig) = ptsrf(ig) + pdtsrf(ig)*ptimestep - - ENDDO - - DO ig=1,ngrid - IF(ig.GT.ngrid/2+1) THEN - icap=2 - ELSE - icap=1 - ENDIF - -! Loop over where we have condensation / sublimation - IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR. & ! ground condensation - (zfallice(ig).NE.0.) .OR. & ! falling snow - ((ztsrf(ig) .GT. ztcondsol(ig)) .AND. & ! ground sublimation - ((piceco2(ig)+zfallice(ig)*ptimestep) .NE. 0.))) THEN - condsub(ig) = .true. - -! Condensation or partial sublimation of CO2 ice - zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) & - /(latcond*ptimestep) - pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep - -! If the entire CO_2 ice layer sublimes -! (including what has just condensed in the atmosphere) - IF((piceco2(ig)/ptimestep+zfallice(ig)).LE. & - -zcondices(ig))THEN - zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig) - pdtsrfc(ig)=(latcond/pcapcal(ig))* & - (zcondices(ig)) - END IF - -! Changing CO2 ice amount and pressure - - zdiceco2(ig) = zcondices(ig) + zfallice(ig) - 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 - ENDDO - -! Surface albedo and emissivity of the ground below the snow (emisref) -! -------------------------------------------------------------------- - do ig =1,ngrid - IF(ig.GT.ngrid/2+1) THEN - icap=2 - ELSE - icap=1 - ENDIF - - if(.not.piceco2(ig).ge.0.) THEN - if(piceco2(ig).le.-1.e-8) print*, & - 'WARNING in condense_co2cloud: piceco2(',ig,')=', piceco2(ig) - piceco2(ig)=0. - endif - if (piceco2(ig).gt.0) then - psolaralb(ig) = albedice(icap) - emisref(ig) = emisice(icap) - else - psolaralb(ig) = albedodat(ig) - emisref(ig) = emissiv - pemisurf(ig) = emissiv - end if - end do - - return - end subroutine condens_co2cloud - -!------------------------------------------------------------------------- - subroutine get_tcond_co2(p,tcond) -! Calculates the condensation temperature for CO2 - - implicit none - -#include "callkeys.h" - - real p, peff, tcond - real, parameter :: ptriple=518000.0 - - peff=p!/co2supsat - - if(peff.lt.ptriple)then - tcond = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula - else - tcond = 684.2-92.3*log(peff)+4.32*log(peff)**2 - ! liquid-vapour transition (based on CRC handbook 2003 data) - endif - return - - end subroutine get_tcond_co2 Index: trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90 (revision 485) +++ trunk/LMDZ.GENERIC/libf/phystd/condense_cloud.F90 (revision 486) @@ -187,13 +187,14 @@ enddo - write(*,*) "condense_co2cloud: i_co2ice=",i_co2ice + write(*,*) "condense_cloud: i_co2ice=",i_co2ice if((i_co2ice.lt.1))then - print*,'In condens_co2cloud but no CO2 ice tracer, exiting.' + print*,'In condens_cloud but no CO2 ice tracer, exiting.' + print*,'Still need generalisation to arbitrary species!' stop endif ccond=cpp/(g*latcond) - print*,'In condens_co2cloud: ccond=',ccond,' latcond=',latcond + print*,'In condens_cloud: ccond=',ccond,' latcond=',latcond ! Prepare special treatment if gas is not pure CO2 Index: trunk/LMDZ.GENERIC/libf/phystd/inifis.F =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/inifis.F (revision 485) +++ trunk/LMDZ.GENERIC/libf/phystd/inifis.F (revision 486) @@ -472,4 +472,9 @@ write(*,*) " sourceevol = ",sourceevol + write(*,*) "Ice evolution timestep ?" + icetstep=100.0 ! default value + call getin("icetstep",icetstep) + write(*,*) " icetstep = ",icetstep + write(*,*) "Snow albedo ?" albedosnow=0.5 ! default value Index: trunk/LMDZ.GENERIC/libf/phystd/newsedim.F =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/newsedim.F (revision 485) +++ trunk/LMDZ.GENERIC/libf/phystd/newsedim.F (revision 486) @@ -79,5 +79,5 @@ PRINT*,'ngrid =',ngrid PRINT*,'ngridmx =',ngridmx - STOP + STOP ENDIF firstcall=.false. @@ -111,5 +111,5 @@ c (stokes law corrected for low pressure by the Cunningham c slip-flow correction according to Rossow (Icarus 36, 1-50, 1978) - + do l=1,nlay do ig=1, ngrid @@ -121,16 +121,4 @@ endif -!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! TEMPORARY MODIF BY RDW !!!! - !rfall=5.e-6 - if(rfall.lt.1.e-7)then - rfall=1.e-7 - endif - !if(rfall.gt.5.e-5)then - ! rfall=5.e-5 - !endif - -!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - vstokes(ig,l) = b * rho * rfall**2 * & (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rfall) @@ -147,5 +135,4 @@ c va traverser le niveau intercouche l : "dztop" est sa hauteur c au dessus de l (m), "ptop" est sa pression (Pa)) - do l=1,nlay do ig=1, ngrid @@ -154,14 +141,28 @@ Ep=0 k=0 - + w(ig,l) = 0. !! JF+AS ajout initialisation (LK MARS) c ************************************************************** c Simple Method - w(ig,l) = - & (1- exp(-dztop*g/(r*pt(ig,l))))*pplev(ig,l) / g -cc write(*,*) 'OK simple method l,w =', l, w(ig,l) -cc write(*,*) 'OK simple method dztop =', dztop -c ************************************************************** +cc w(ig,l) = +cc & (1- exp(-dztop*g/(r*pt(ig,l))))*pplev(ig,l) / g +cc write(*,*) 'OK simple method l,w =', l, w(ig,l) +cc write(*,*) 'OK simple method dztop =', dztop + w(ig,l) = 1. - exp(-dztop*g/(r*pt(ig,l))) + !!! Diagnostic: JF. Fix: AS. Date: 05/11 + !!! Probleme arrondi avec la quantite ci-dessus + !!! ---> vaut 0 pour -dztop*g/(r*pt(ig,l)) trop petit + !!! ---> dans ce cas on utilise le developpement limite ! + !!! ---> exp(-x) = 1 - x lorsque x --> 0 avec une erreur de x^2 / 2 + + IF ( w(ig,l) .eq. 0. ) THEN + w(ig,l) = ( dztop*g/(r*pt(ig,l)) ) * pplev(ig,l) / g + ELSE + w(ig,l) = w(ig,l) * pplev(ig,l) / g + ENDIF +! LK borrowed simple method from Mars model (AS/JF) + +!************************************************************** cccc Complex method : - if (dztop.gt.epaisseur(ig,l)) then + if (dztop.gt.epaisseur(ig,l)) then cccc Cas ou on "epuise" la couche l : On calcule le flux cccc Venant de dessus en tenant compte de la variation de Vstokes @@ -176,7 +177,16 @@ enddo Ep = Ep - epaisseur(ig,l+k) - end if - ptop=pplev(ig,l+k)*exp(-(dztop-Ep)*g/(r*pt(ig,l+k))) - w(ig,l) = (pplev(ig,l) -Ptop)/g +! ptop=pplev(ig,l+k)*exp(-(dztop-Ep)*g/(r*pt(ig,l+k))) + ptop=exp(-(dztop-Ep)*g/(r*pt(ig,l+k))) + IF ( ptop .eq. 1. ) THEN + !PRINT*, 'newsedim: exposant trop petit ', ig, l + ptop=pplev(ig,l+k) * ( 1. - (dztop-Ep)*g/(r*pt(ig,l+k))) + ELSE + ptop=pplev(ig,l+k) * ptop + ENDIF + + w(ig,l) = (pplev(ig,l) - ptop)/g + + endif !!! complex method c cc write(*,*) 'OK new method l,w =', l, w(ig,l) @@ -188,4 +198,5 @@ cc if(l.eq.8)write(*,*)'l=8,k, w',k, w(1,l) c ************************************************************** + end do end do @@ -195,6 +206,4 @@ c & wq(1,6),wq(1,7),pqi(1,6) - RETURN END - Index: trunk/LMDZ.GENERIC/libf/phystd/physiq.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/physiq.F90 (revision 485) +++ trunk/LMDZ.GENERIC/libf/phystd/physiq.F90 (revision 486) @@ -1917,5 +1917,5 @@ qsurf_hist(ig,igcm_h2o_ice) = & !qsurf_hist(ig,igcm_h2o_ice) + delta_ice*100.0 - qsurf_hist(ig,igcm_h2o_ice) + delta_ice*10.0 + qsurf_hist(ig,igcm_h2o_ice) + delta_ice*icetstep ! if ice has gone -ve, set to zero Index: trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90 =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90 (revision 485) +++ trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90 (revision 486) @@ -6,5 +6,5 @@ #include "bands.h" -!======================================================================C +!====================================================================== ! ! RADINC.H RADiation INCludes @@ -13,5 +13,5 @@ ! number of spectral intervals. . . ! -!======================================================================C +!====================================================================== ! RADIATION parameters @@ -60,9 +60,8 @@ integer, parameter :: L_NLEVRAD = llm+1 - !!!! THIS IS SET IN sugas_corrk - !!!! [use of allocatable arrays] -- AS 12/2011 + ! These are set in sugas_corrk + ! [uses allocatable arrays] -- AS 12/2011 integer :: L_NPREF, L_NTREF, L_REFVAR, L_PINT -!!! THIS ONE DOES NOT CHANGE MUCH... IT CAN BE REGARDED AS A CONSTANT. integer, parameter :: L_NGAUSS = 17 @@ -72,6 +71,4 @@ integer, parameter :: NAERKIND = 2 real, parameter :: L_TAUMAX = 35 - !integer, parameter :: L_TAUMAX = 35 - !integer, parameter :: L_TAUMAX = 40 ! For Planck function integration: Index: trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F (revision 485) +++ trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F (revision 486) @@ -112,7 +112,4 @@ REAL cloudfrac(ngridmx,nlayermx) REAL hice(ngridmx),totcloudfrac(ngridmx) - - -c======================================================================= c======================================================================= @@ -613,5 +610,5 @@ DO idt=1,ndt - IF (idt.eq.ndt-day_step-1) then !test + IF (idt.eq.ndt) then !test lastcall=.true. call stellarlong(day*1.0,zls) Index: trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F =================================================================== --- trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F (revision 485) +++ trunk/LMDZ.GENERIC/libf/phystd/sfluxi.F (revision 486) @@ -31,6 +31,4 @@ real*8 FZEROI(L_NSPECTI) real*8 taugsurf(L_NSPECTI,L_NGAUSS-1), fzero - -! real*8 BSURFtest ! by RW for test real*8 fup_tmp(L_NSPECTI),fdn_tmp(L_NSPECTI) @@ -73,6 +71,4 @@ ! NTS = TSURF*10.0D0-499 ! NTT = TTOP *10.0D0-499 - -! BSURFtest=0.0 DO 501 NW=1,L_NSPECTI