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 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" #include "gases.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