subroutine condense_co2(ngrid,nlayer,nq,ptimestep, & pcapcal,pplay,pplev,ptsrf,pt, & pdt,pdtsrf,pq,pdq, & pqsurf,pdqsurfc,albedo,pemisurf, & albedo_bareground,albedo_co2_ice_SPECTV, & pdtc,pdtsrfc,pdpsrfc,pdqc) use radinc_h, only : L_NSPECTV, naerkind use gases_h, only: gfrac, igas_co2 use radii_mod, only : co2_reffrad use aerosol_mod, only : iaero_co2 USE surfdat_h, only: emisice, emissiv USE comgeomphy, only: latitude ! in radians USE tracer_h, only: noms, rho_co2 use comcstfi_mod, only: g, r, cpp 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 vertical layers. ! nq Number of tracers. ! ptimestep Duration of the physical timestep (s). ! pplay(ngrid,nlayer) Pressure layers (Pa). ! pplev(ngrid,nlayer+1) Pressure levels (Pa). ! pt(ngrid,nlayer) Atmospheric Temperatures (K). ! ptsrf(ngrid) Surface temperatures (K). ! pq(ngrid,nlayer,nq) Atmospheric tracers mixing ratios (kg/kg of air). ! pqsurf(ngrid,nq) Surface tracers (kg/m2). ! ! pdt(ngrid,nlayer) Time derivative before condensation/sublimation of pt. ! pdtsrf(ngrid) Time derivative before condensation/sublimation of ptsrf. ! pdq(ngrid,nlayer,nq) Time derivative before condensation/sublimation of ! ! albedo_bareground(ngrid) Albedo of the bare ground. ! albedo_co2_ice_SPECTV(L_NSPECTV) Spectral albedo of CO2 ice. ! ! Outputs ! ------- ! pdpsrfc(ngrid) \ Contribution of condensation/sublimation ! pdtc(ngrid,nlayer) \ to the time derivatives of ! pdtsrfc(ngrid) / Surface Pressure, Atmospheric Temperatures, ! pdqsurfc(ngrid) / Surface Temperatures, Surface Tracers, ! pdqc(ngrid,nlayer,nq) / and Atmospheric Tracers.* ! ! pemisurf(ngrid) Emissivity of the surface. ! ! Both ! ---- ! albedo(ngrid,L_NSPECTV) Spectral albedo of the surface. ! ! Authors ! ------- ! Francois Forget (1996) ! Converted to Fortran 90 and slightly modified by R. Wordsworth (2009) ! Includes simplifed nucleation by J. Leconte (2011) ! !================================================================== !-------------------------- ! Arguments !-------------------------- INTEGER,INTENT(IN) :: ngrid INTEGER,INTENT(IN) :: nlayer INTEGER,INTENT(IN) :: nq REAL,INTENT(IN) :: ptimestep REAL,INTENT(IN) :: pcapcal(ngrid) REAL,INTENT(IN) :: pplay(ngrid,nlayer) REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) REAL,INTENT(IN) :: ptsrf(ngrid) REAL,INTENT(IN) :: pt(ngrid,nlayer) REAL,INTENT(IN) :: pdt(ngrid,nlayer) REAL,INTENT(IN) :: pdtsrf(ngrid) REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) REAL,INTENT(IN) :: pqsurf(ngrid,nq) REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq) REAL,INTENT(IN) :: albedo_bareground(ngrid) REAL,INTENT(IN) :: albedo_co2_ice_SPECTV(L_NSPECTV) REAL,INTENT(INOUT) :: albedo(ngrid,L_NSPECTV) REAL,INTENT(OUT) :: pemisurf(ngrid) REAL,INTENT(OUT) :: pdtc(ngrid,nlayer) REAL,INTENT(OUT) :: pdtsrfc(ngrid) REAL,INTENT(OUT) :: pdpsrfc(ngrid) REAL,INTENT(OUT) :: pdqc(ngrid,nlayer,nq) REAL,INTENT(OUT) :: pdqsurfc(ngrid) !------------------------------ ! Local variables !------------------------------ INTEGER l,ig,icap,ilay,iq,nw,igas,it REAL reffrad(ngrid,nlayer) ! Radius (m) of the CO2 ice particles. REAL*8 zt(ngrid,nlayer) ! Updated Atmospheric Temperatures (K). REAL ztsrf(ngrid) ! Updated Surface Temperatures (K). REAL zq(ngrid,nlayer,nq) ! Updated Atmospheric tracers mixing ratios (kg/kg of air). REAL piceco2(ngrid) ! Updated Surface Tracer (kg/m2). REAL ztcond (ngrid,nlayer) ! Atmospheric Temperatures of condensation of CO2. REAL ztnuc (ngrid,nlayer) ! Atmospheric Nucleation Temperatures. REAL ztcondsol(ngrid) ! Temperatures of condensation of CO2 at the surface. REAL zcondices(ngrid) ! Condensation rate on the ground (kg/m2/s). REAL zfallice(ngrid) ! Flux of ice falling on the surface (kg/m2/s). REAL Mfallice(ngrid) ! Total amount of ice fallen to the ground during the timestep (kg/m2). REAL wq(ngrid,nlayer+1) ! Total amount of ice fallen to the ground during the timestep (kg/m2). REAL subptimestep ! Duration of the subtimestep (s) for the sedimentation. Integer Ntime ! Number of subtimesteps. REAL masse (ngrid,nlayer) ! Mass of atmospheric layers (kg/m2) REAL w(ngrid,nlayer,nq) ! REAL vstokes,reff ! REAL ppco2 ! !------------------------------------------ ! Saved local variables !------------------------------------------ REAL,SAVE :: latcond=5.9e5 REAL,SAVE :: ccond REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: emisref !$OMP THREADPRIVATE(latcond,ccond,emisref) LOGICAL,SAVE :: firstcall=.true. !$OMP THREADPRIVATE(firstcall) INTEGER,SAVE :: i_co2ice=0 ! co2 ice !$OMP THREADPRIVATE(i_co2ice) CHARACTER(LEN=20) :: tracername ! to temporarily store text !------------------------------------------------ ! Initialization at the first call !------------------------------------------------ IF (firstcall) THEN ALLOCATE(emisref(ngrid)) ! Find CO2 ice tracer. do iq=1,nq tracername=noms(iq) if (tracername.eq."co2_ice") then i_co2ice=iq endif enddo write(*,*) "condense_co2: i_co2ice=",i_co2ice if((i_co2ice.lt.1))then 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_cloud: 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 !------------------------------------------------ ! Tendencies initially set to 0 !------------------------------------------------ pdqc(1:ngrid,1:nlayer,1:nq) = 0. pdtc(1:ngrid,1:nlayer) = 0. zq(1:ngrid,1:nlayer,1:nq) = 0. zt(1:ngrid,1:nlayer) = 0. Mfallice(1:ngrid) = 0. zfallice(1:ngrid) = 0. zcondices(1:ngrid) = 0. pdtsrfc(1:ngrid) = 0. pdpsrfc(1:ngrid) = 0. pdqsurfc(1:ngrid) = 0. !---------------------------------- ! 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' and nucleation temperature 'ztnuc'. DO l=1,nlayer DO ig=1,ngrid ppco2=gfrac(igas_CO2)*pplay(ig,l) call get_tcond_co2(ppco2,ztcond(ig,l)) call get_tnuc_co2(ppco2,ztnuc(ig,l)) ENDDO ENDDO ! Initialize zq and zt at the beginning of the sub-timestep loop and qsurf. DO ig=1,ngrid piceco2(ig)=pqsurf(ig,i_co2ice) DO l=1,nlayer 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) ! Add the tendencies from other physical processes at each subtimstep. DO it=1,Ntime 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 starts. ! Sedimentation computed from radius computed from q in module radii_mod. call co2_reffrad(ngrid,nlayer,nq,zq,reffrad) DO ilay=1,nlayer DO ig=1,ngrid reff = reffrad(ig,ilay) 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,nlayer,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 !---------------------------------------------------------- ! (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. ! ztcond-> ztnuc in test beneath to nucleate only when super saturation occurs(JL 2011) IF ((zt(ig,l).LT.ztnuc(ig,l)).or.(zq(ig,l,i_co2ice).gt.1.E-10)) THEN 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 !----------------------------------------------------------------------- ! Forecast of ground temperature ztsrf and frost temperature ztcondsol. DO ig=1,ngrid ppco2=gfrac(igas_CO2)*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 ! 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 piceco2(ig) = piceco2(ig) + pdqsurfc(ig)*ptimestep pdqsurfc(ig) = zcondices(ig) + zfallice(ig) pdpsrfc(ig) = -pdqsurfc(ig)*g IF(ABS(pdpsrfc(ig)*ptimestep).GT.pplev(ig,1)) THEN PRINT*,'STOP in condens in condense_co2' PRINT*,'condensing more than total mass' PRINT*,'Grid point ',ig PRINT*,'Ps = ',pplev(ig,1) PRINT*,'d Ps = ',pdpsrfc(ig) STOP ENDIF END IF ENDDO ! end of ngrid loop. !--------------------------------------------------------------------------------------------- ! Surface albedo and emissivity of the ground below the snow (emisref) !--------------------------------------------------------------------------------------------- DO ig=1,ngrid IF(latitude(ig).LT.0.) THEN icap=2 ! Southern Hemisphere ELSE icap=1 ! Nortnern hemisphere 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. 1.) then ! CO2 Albedo condition changed to ~1 mm coverage. Change by MT2015. DO nw=1,L_NSPECTV albedo(ig,nw) = albedo_co2_ice_SPECTV(nw) ENDDO emisref(ig) = emisice(icap) else DO nw=1,L_NSPECTV albedo(ig,nw) = albedo_bareground(ig) ! Note : If you have some water, it will be taken into account in the "hydrol" routine. ENDDO emisref(ig) = emissiv pemisurf(ig) = emissiv end if END DO return end subroutine condense_co2 !------------------------------------------------------------------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- subroutine get_tcond_co2(p,tcond) ! Calculates the condensation temperature for CO2 implicit none real p, peff, tcond real, parameter :: ptriple=518000.0 peff=p 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 !------------------------------------------------------------------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- subroutine get_tnuc_co2(p,tnuc) ! Calculates the nucleation temperature for CO2, based on a simple super saturation criterion. JL 2011. use callkeys_mod, only: co2supsat implicit none real p, peff, tnuc real, parameter :: ptriple=518000.0 peff=p/co2supsat if(peff.lt.ptriple) then tnuc = (-3167.8)/(log(.01*peff)-23.23) ! Fanale's formula else tnuc = 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_tnuc_co2