MODULE co2condens_mod IMPLICIT NONE logical, save :: scavco2cond = .false. ! flag for using scavenging_by_co2 !$OMP THREADPRIVATE(scavco2cond) real, save :: CO2cond_ps = 1. ! Coefficient to control the surface pressure change CONTAINS SUBROUTINE co2condens(ngrid,nlayer,nq,nslope,ptimestep, $ pcapcal,pplay,pplev,ptsrf,pt, $ pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, $ piceco2,perennial_co2ice, $ psolaralb,pemisurf,rdust, $ pdtc,pdtsrfc,pdpsrf,pduc,pdvc,pdqc, $ fluxsurf_sw,zls, $ zdqssed_co2,pcondicea_co2microp, $ pdqsc) use tracer_mod, only: noms, igcm_h2o_ice, igcm_h2o_vap, & igcm_dust_mass, igcm_dust_number, & igcm_ccn_mass, igcm_ccn_number, & igcm_hdo_ice, igcm_hdo_vap, & nqperes,nqfils, ! MVals: variables isotopes & qperemin,masseqmin, & igcm_co2 use surfdat_h, only: emissiv use geometry_mod, only: latitude, ! grid point latitudes (rad) & longitude_deg, latitude_deg use planete_h, only: obliquit use comcstfi_h, only: cpp, g, r, pi use dust_param_mod, only: freedust use write_output_mod, only: write_output #ifndef MESOSCALE USE vertical_layers_mod, ONLY: ap, bp #endif use comslope_mod, ONLY: subslope_dist,def_slope_mean USE paleoclimate_mod, ONLY: paleoclimate 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 ; updated 1996 -- 2018 c ------ c adapted to external CO2 ice clouds scheme by Deborah Bardet (2018) ' 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 INTEGER,INTENT(IN) :: nslope ! number of subslope REAL,INTENT(IN) :: ptimestep ! physics timestep (s) REAL,INTENT(IN) :: pcapcal(ngrid,nslope) 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,nslope) ! 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,nslope) ! 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(IN) :: zdqssed_co2(ngrid) ! CO2 flux at the surface (kg.m-2.s-1) REAL,INTENT(IN) :: pcondicea_co2microp(ngrid,nlayer)! tendency due to CO2 condensation (kg/kg.s-1) REAL,INTENT(INOUT) :: piceco2(ngrid,nslope) ! CO2 ice on the surface (kg.m-2) REAL,INTENT(INOUT) :: perennial_co2ice(ngrid,nslope) ! Perennial CO2 ice on the surface (kg.m-2) REAL,INTENT(INOUT) :: psolaralb(ngrid,2,nslope) ! albedo of the surface REAL,INTENT(INOUT) :: pemisurf(ngrid,nslope) ! emissivity of the surface REAL,INTENT(IN) :: rdust(ngrid,nlayer) ! dust effective radius ! tendencies due to CO2 condensation/sublimation: REAL,INTENT(OUT) :: pdtc(ngrid,nlayer) ! tendency on temperature (K/s) REAL,INTENT(OUT) :: pdtsrfc(ngrid,nslope) ! 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 REAL,INTENT(OUT) :: pdqsc(ngrid,nq,nslope) ! tendency on surface 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 ----------------- INTEGER i,j INTEGER l,ig,iq,icap 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,nslope) REAL zdiceco2_mesh_avg(ngrid) REAL zcondicea(ngrid,nlayer) ! condensation rate in layer l (kg/m2/s) REAL zcondices(ngrid,nslope) ! condensation rate on the ground (kg/m2/s) REAL zcondices_mesh_avg(ngrid) ! condensation rate on the ground (kg/m2/s) REAL zfallice(ngrid,nlayer+1) ! amount of ice falling from layer l (kg/m2/s) REAL condens_layer(ngrid,nlayer) ! co2clouds: condensation rate in layer l (kg/m2/s) REAL condens_column(ngrid) ! co2clouds: sum(condens_layer(ig,:)) (kg/m2/s) REAL zfallheat REAL zmflux(nlayer+1) REAL zu(nlayer),zv(nlayer) REAL zqc(nlayer,nq),zq1(nlayer) REAL ztsrf(ngrid,nslope) 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 REAL availco2 LOGICAL condsub(ngrid,nslope) real :: emisref(ngrid,nslope) REAL zdq_scav(ngrid,nlayer,nq) ! tendency due to scavenging by co2 REAL zq(ngrid,nlayer,nq) 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. c local saved variables integer,save :: ico2 ! index of CO2 tracer real,save :: 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 :: m_co2, m_noco2, A , B LOGICAL,SAVE :: firstcall = .true. !,firstcall2=.true. !$OMP THREADPRIVATE(ico2,qco2,mmean,acond,bcond,ccond,m_co2,m_noco2) !$OMP THREADPRIVATE(A,B,firstcall) c D.BARDET: for debug real ztc3D(ngrid,nlayer) REAL ztm3D(ngrid,nlayer) REAL zmflux3D(ngrid,nlayer) c MVals: variables isotopes REAL Ratio1(nlayer),Ratiom1(nlayer+1) REAL masseq(nlayer),wq(nlayer+1) INTEGER ifils,iq2 c Subslope: REAL :: emisref_tmp(ngrid) ! temporary emissiviy of a sub-grid surface before computing emissivity with snow [1] REAL :: pemisurf_tmp(ngrid) ! temporary emissiviy of a sub-grid surface after computing emissivity with snow [1] REAL :: alb_tmp(ngrid,2) ! temporary albedo before computing albedo with TES/fixed value [1] REAL :: zcondices_tmp(ngrid) ! temporary condensation rate [kg/m^2/s] REAL :: piceco2_tmp(ngrid) ! temporary amount of CO2 frost on the surface [kg/m^2] REAL :: perennial_co2ice_tmp(ngrid) ! temporary amount of perennial CO2 frost on the surface [kg/m^2] LOGICAL :: condsub_tmp(ngrid) ! Boolean to check if CO2 ice is condensing / sublimating on the sub grid surface [1] REAL :: zfallice_tmp(ngrid,nlayer+1) ! temporary amount of ice falling from layer l for a specific sub-grid surface [kg/m^2/s] REAL :: condens_layer_tmp(ngrid,nlayer) ! temporary condensation rate in layer l (co2 cloud) for a specific sub-grid surface [kg/m^2/s] INTEGER :: islope ! index for loop variables REAL :: pdqsc_tmp(ngrid,nq) ! tendency on surface tracers (grid-mesh average) after scavenging c---------------------------------------------------------------------- c Initialisation c -------------- c ! AS: firstcall OK absolute IF (firstcall) THEN bcond=1./tcond1mb ccond=cpp/(g*latcond) acond=r/latcond firstcall=.false. write(*,*) 'CO2condens: improved_ztcond=',improved_ztcond PRINT*,'In co2condens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond PRINT*,'acond,bcond,ccond',acond,bcond,ccond ico2=0 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 ENDIF ! of IF (firstcall) zcpi=1./cpp if (paleoclimate) piceco2 = piceco2 + perennial_co2ice 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 c ------------------------------------- zcondicea(1:ngrid,1:nlayer) = 0. zfallice(1:ngrid,1:nlayer+1) = 0. pduc(1:ngrid,1:nlayer) = 0 pdvc(1:ngrid,1:nlayer) = 0 pdtc(1:ngrid,1:nlayer) = 0. pdqsc(1:ngrid,1:nq,1:nslope) = 0 pdqsc_tmp(1:ngrid,1:nq) = 0 pdqc(1:ngrid,1:nlayer,1:nq) = 0 zcondices(1:ngrid,1:nslope) = 0. zcondices_mesh_avg(1:ngrid)=0. pdtsrfc(1:ngrid,1:nslope) = 0. pdpsrf(1:ngrid) = 0. condsub(1:ngrid,1:nslope) = .false. zdiceco2(1:ngrid,1:nslope) = 0. zdiceco2_mesh_avg(1:ngrid)=0. zfallheat=0 zdq_scav(:,:,:)=0. c Update tendencies from previous processes c ------------------------------------- DO l=1,nlayer DO ig=1,ngrid zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep do iq=1,nq zq(ig,l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep enddo ENDDO ENDDO 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 endif IF (.NOT. co2clouds) then c forecast of atmospheric temperature zt and frost temperature ztcond c -------------------------------------------------------------------- DO l=1,nlayer DO ig=1,ngrid ! 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 condens_layer(:,:) = zcondicea(:,:) condens_column(:) = 0. if (scavco2cond) then call scavenging_by_co2(ngrid,nlayer,nq,ptimestep,pplev,zq, & rdust,zcondicea,zfallice,zdq_scav, & pdqsc_tmp) DO islope = 1,nslope pdqsc(:,:,islope)= pdqsc_tmp(:,:)* & cos(pi*def_slope_mean(islope)/180.) ENDDO endif call write_output("co2condens_zfallice", & "CO2 ice tendency on the surface", & "kg.m-2.s-1",zfallice(:,1)) ELSE ! if co2 clouds condens_layer(:,:) = 0. condens_column(:) = 0. DO l=nlayer , 1, -1 DO ig=1,ngrid condens_layer(ig,l) = pcondicea_co2microp(ig,l)* & (pplev(ig,l) - pplev(ig,l+1))/g ENDDO ENDDO DO ig=1,ngrid condens_column(ig) = sum(condens_layer(ig,:)) zfallice(ig,1) = zdqssed_co2(ig) DO islope = 1,nslope piceco2(ig,islope) = piceco2(ig,islope) + & zdqssed_co2(ig)*ptimestep * & cos(pi*def_slope_mean(islope)/180.) ENDDO ENDDO call write_output("co2condens_zfallice", & "CO2 ice tendency on the surface", & "kg.m-2.s-1",zdqssed_co2(:)) ! otherwise we have not 1 day(1proc) = 1 day (24procs) test ENDIF ! end of if co2clouds call write_output("co2condens_pdtc", & "Temperature tendency due to CO2 condensation", & "K.s-1",pdtc(:,:)) ! call write_output("condens_layer", ! & "", ! & " ",condens_layer(:,:)) 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))) DO islope = 1,nslope ztsrf(ig,islope) = ptsrf(ig,islope) + & pdtsrf(ig,islope)*ptimestep ENDDO ENDDO c c Condensation/sublimation on the ground c -------------------------------------- c (compute zcondices and pdtsrfc) c c No microphysics of CO2 clouds DO ig=1,ngrid IF(latitude(ig).lt.0) THEN ! Southern hemisphere icap=2 ELSE ! Northern hemisphere icap=1 ENDIF DO islope = 1,nslope c Need first to put piceco2_slope(ig,islope) in kg/m^2 flat piceco2(ig,islope) = piceco2(ig,islope) & /cos(pi*def_slope_mean(islope)/180.) c c Loop on where we have condensation/ sublimation IF ((ztsrf(ig,islope) .LT. ztcondsol(ig)) .OR. ! ground cond $ (zfallice(ig,1).NE.0.) .OR. ! falling snow $ ((ztsrf(ig,islope) .GT. ztcondsol(ig)) .AND. ! ground sublim. $ ((piceco2(ig,islope)+zfallice(ig,1)*ptimestep) $ .NE. 0.))) THEN condsub(ig,islope) = .true. IF (zfallice(ig,1).gt.0 .and. .not. co2clouds) then ! potential eneregy release due to the impact of the snowfall ! appendix of forget et al. 1999 zfallheat = zfallice(ig,1) * (pphi(ig,1) + & cpice*(ztcond(ig,1)-ztcondsol(ig)))/latcond ELSE zfallheat = 0. ENDIF c condensation or partial sublimation of CO2 ice c """"""""""""""""""""""""""""""""""""""""""""""" if(ztsrf(ig,islope).LT. ztcondsol(ig)) then c condensation zcondices(ig,islope)=pcapcal(ig,islope) & *(ztcondsol(ig)-ztsrf(ig,islope)) & /((latcond+cpp*(zt(ig,1)-ztcondsol(ig)))*ptimestep) & - zfallheat else c sublimation zcondices(ig,islope)=pcapcal(ig,islope) & *(ztcondsol(ig)-ztsrf(ig,islope)) & /(latcond*ptimestep) & - zfallheat endif pdtsrfc(ig,islope) = (ztcondsol(ig) - ztsrf(ig,islope)) & /ptimestep #ifdef MESOSCALE print*, "not enough CO2 tracer in 1st layer to condense" print*, ">>> to be implemented in the mesoscale case" print*, "because this uses ap levels..." #else c If there is not enough CO2 tracer in 1st layer to condense c """""""""""""""""""""""""""""""""""""""""""""""""""""" IF(ico2.ne.0) then c Available CO2 tracer in layer 1 at end of timestep (kg/m2) #ifndef MESOSCALE availco2 = pq(ig,1,ico2)*((ap(1)-ap(2))+ & (bp(1)-bp(2))*(pplev(ig,1)/g - & (zcondices(ig,islope) + zfallice(ig,1)) & *ptimestep)) if ((zcondices(ig,islope) + condens_layer(ig,1)) & *ptimestep & .gt.availco2) then zcondices(ig,islope) = availco2/ptimestep - & condens_layer(ig,1) pdtsrfc(ig,islope) = (latcond/pcapcal(ig,islope))* & (zcondices(ig,islope)+zfallheat) end if #else availco2 = pq(ig,1,igcm_co2) PRINT*, "MESOSCALE: CO2 tracer AT FIRST LEVEL IS NOT' & ' CORRECTED FROM SIGMA LEVELS" #endif ENDIF #endif c If the entire CO2 ice layer sublimes on the slope c """""""""""""""""""""""""""""""""""""""""""""""""""" c (including what has just condensed in the atmosphere) if (co2clouds) then IF((piceco2(ig,islope)/ptimestep).LE. & -zcondices(ig,islope))THEN zcondices(ig,islope) = -piceco2(ig,islope)/ptimestep pdtsrfc(ig,islope)=(latcond/pcapcal(ig,islope)) * & (zcondices(ig,islope)+zfallheat) END IF else IF((piceco2(ig,islope)/ptimestep+zfallice(ig,1)).LE. & -zcondices(ig,islope))THEN zcondices(ig,islope) = -piceco2(ig,islope) & /ptimestep - zfallice(ig,1) pdtsrfc(ig,islope)=(latcond/pcapcal(ig,islope)) * & (zcondices(ig,islope)+zfallheat) END IF end if c Changing CO2 ice amount and pressure per slope: c """""""""""""""""""""""""""""""""""" zdiceco2(ig,islope) = zcondices(ig,islope)+zfallice(ig,1) & + condens_column(ig) if (co2clouds) then ! add here only direct condensation/sublimation piceco2(ig,islope) = piceco2(ig,islope) + & zcondices(ig,islope)*ptimestep else ! add here also CO2 ice in the atmosphere piceco2(ig,islope) = piceco2(ig,islope) + & zdiceco2(ig,islope)*ptimestep end if zcondices_mesh_avg(ig) = zcondices_mesh_avg(ig) + & zcondices(ig,islope)* subslope_dist(ig,islope) zdiceco2_mesh_avg(ig) = zdiceco2_mesh_avg(ig) + & zdiceco2(ig,islope)* subslope_dist(ig,islope) END IF ! if there is condensation/sublimation piceco2(ig,islope) = piceco2(ig,islope) & * cos(pi*def_slope_mean(islope)/180.) ENDDO !islope pdpsrf(ig) = -zdiceco2_mesh_avg(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*,'Longitude(degrees): ',longitude_deg(ig) PRINT*,'Latitude(degrees): ',latitude_deg(ig) PRINT*,'Ps = ',pplev(ig,1) PRINT*,'d Ps = ',pdpsrf(ig) CALL abort_physic('co2condens', & 'condensing more than total mass', 1) ENDIF ENDDO ! of DO ig=1,ngrid c ******************************************************************** c Surface albedo and emissivity of the surface below the snow (emisref) c ******************************************************************** ! Redistribute piceco2 into piceco2 (frost) and perennial_co2ice ! -------------------------------------------------------------- if (paleoclimate) then where (piceco2 > perennial_co2ice) ! Perennial co2 ice has not been affected ! It means: ! - In case of sublimation, only frost is lost ! - In case of condensation, only frost accumulates new ice piceco2 = piceco2 - perennial_co2ice else where ! Perennial co2 ice has been affected ! It means that frost disappeared with sublimation and perennial ice is being lost perennial_co2ice = piceco2 piceco2 = 0. end where endif ! Check that amont of CO2 ice is not problematic ! ---------------------------------------------- DO ig = 1,ngrid DO islope = 1,nslope if (piceco2(ig,islope) < 0.) then if (piceco2(ig,islope) <= -5.e-8) print*, $ 'WARNING co2condens piceco2(',ig,',',islope, $ ') =',piceco2(ig,islope) piceco2(ig,islope) = 0. endif ENDDO ENDDO if (paleoclimate) then DO ig = 1,ngrid DO islope = 1,nslope if (perennial_co2ice(ig,islope) < 0.) then if (perennial_co2ice(ig,islope) <= -5.e-8) print*, $ 'WARNING co2condens perennial_co2ice(',ig,',',islope, $ ') =',perennial_co2ice(ig,islope) perennial_co2ice(ig,islope) = 0. endif ENDDO ENDDO endif ! Set albedo and emissivity of the surface ! ---------------------------------------- DO islope = 1,nslope piceco2_tmp = piceco2(:,islope) alb_tmp = psolaralb(:,:,islope) emisref_tmp = 0. perennial_co2ice_tmp = perennial_co2ice(:,islope) CALL albedocaps(zls,ngrid,piceco2_tmp,perennial_co2ice_tmp, & alb_tmp,emisref_tmp) perennial_co2ice(:,islope) = perennial_co2ice_tmp psolaralb(:,1,islope) = alb_tmp(:,1) psolaralb(:,2,islope) = alb_tmp(:,2) emisref(:,islope) = emisref_tmp ENDDO ! Set pemisurf() to emissiv when there is bare surface (needed for co2snow) ! ------------------------------------------------------------------------- where (piceco2 == 0.) pemisurf = emissiv ! 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 (any(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 zqc(l,iq)=zq(ig,l,iq)+zdq_scav(ig,l,iq)*ptimestep ! zdq_scav=0 if co2clouds=true enddo enddo c Mass fluxes through the sigma levels (kg.m-2.s-1) (>0 when up) c """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" zmflux(1) = -zcondices_mesh_avg(ig) - zdqssed_co2(ig) DO l=1,nlayer zmflux(l+1) = zmflux(l) - condens_layer(ig,l) #ifndef MESOSCALE & + (bp(l)-bp(l+1))*(-pdpsrf(ig)/g) 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.) then zmflux(l+1)=0. end if #else zmflux(l+1) = zmflux(l) - condens_layer(ig,l) if (abs(zmflux(l+1)).lt.1E-13) zmflux(l+1)=0. PRINT*, "MESOSCALE: FLUX THROUGH SIGMA LEVELS from"// & "dPS HAVE TO BE IMPLEMENTED" #endif END DO #ifdef MESOSCALE print*,"absurd mass set because bp not available" print*,"TO BE FIXED" #else c Mass of each layer at the end of timestep c ----------------------------------------- DO l=1,nlayer masse(l)=( pplev(ig,l) - pplev(ig,l+1) + & (bp(l)-bp(l+1))*pdpsrf(ig)*ptimestep)/g END DO #endif 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: if (zmflux(1).lt.0) then c Surface condensation ztm(1) = ztc(1) else c Surface sublimation: ztm(1) = ztcondsol(ig) endif 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) ! MVals: loop over the fathers ("peres") do iq=1,nqperes do l=1,nlayer zq1(l)=zqc(l,iq) enddo zqm1(1)=zqm(1,iq) call vl1d(nlayer,zq1,2.,masse,w,zqm1) do l=2,nlayer zqc(l,iq)=zq1(l) zqm(l,iq)=zqm1(l) enddo ! MVals: loop over the sons ("fils") if (nqfils(iq).gt.0) then if (iq.eq.igcm_h2o_ice) then iq2=igcm_hdo_ice else if (iq.eq.igcm_h2o_vap) then iq2=igcm_hdo_vap else call abort_physic("co2condens_mod","invalid isotope",1) endif do l=1,nlayer if (zqc(l,iq).gt.qperemin) then Ratio1(l)=zqc(l,iq2)/zqc(l,iq) else Ratio1(l)=0. endif masseq(l)=max(masse(l)*zqc(l,iq),masseqmin) wq(l)=w(l)*zqm(l,iq) enddo Ratiom1(1)=zqm(1,iq2) call vl1d(nlayer,Ratio1,2.,masseq,wq,Ratiom1) zqm(1,iq2)=Ratiom1(1)*zqc(1,iq) do l=2,nlayer zqm(l,iq2)=Ratiom1(l)*zqm(l,iq) enddo endif !if (nqfils(iq).gt.0) then enddo !iq=1,nqperes 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)= zqc(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 IF(.not. co2clouds) THEN c Tendencies on T zdtsig(ig,l) = (1/masse(l)) * & ( zmflux(l)*(ztm(l) - ztc(l)) & - zmflux(l+1)*(ztm(l+1) - ztc(l)) & + condens_layer(ig,l)*(ztcond(ig,l)-ztc(l)) ) ELSE zdtsig(ig,l) = (1/masse(l)) * & ( zmflux(l)*(ztm(l) - ztc(l)) & - zmflux(l+1)*(ztm(l+1) - ztc(l))) ENDIF c D.BARDET: for diagnotics zmflux3D(ig,l)=zmflux(l) ztm3D(ig,l)=ztm(l) ztc3D(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 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) - zqc(l,iq)) & - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq)) & + condens_layer(ig,l)*(zqc(l,iq)-1.) ) END DO else DO l=1,nlayer pdqc(ig,l,iq)= (1/masse(l)) * & ( zmflux(l)*(zqm(l,iq) - zqc(l,iq)) & - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq)) & + condens_layer(ig,l)*zqc(l,iq) ) pdqc(ig,l,iq)=pdqc(ig,l,iq)+zdq_scav(ig,l,iq) ! zdq_scav=0 if co2clouds=true END DO end if enddo end if ! if (condsub) END DO ! loop on ig c *************************************************************** c CO2 snow / clouds scheme c *************************************************************** DO islope = 1,nslope emisref_tmp(:) = emisref(:,islope) condsub_tmp(:) = condsub(:,islope) condens_layer_tmp(:,:) = condens_layer(:,:)* & cos(pi*def_slope_mean(islope)/180.) zcondices_tmp(:) = zcondices(:,islope)* & cos(pi*def_slope_mean(islope)/180.) zfallice_tmp(:,:) = zfallice(:,:)* & cos(pi*def_slope_mean(islope)/180.) pemisurf_tmp(:) = pemisurf(:,islope) call co2snow(ngrid,nlayer,ptimestep,emisref_tmp,condsub_tmp, & pplev,condens_layer_tmp,zcondices_tmp,zfallice_tmp, & pemisurf_tmp) pemisurf(:,islope) = pemisurf_tmp(:) ENDDO 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 write_output('tconda1', c &'Taux de condensation CO2 atmospherique /Pa', c & 'kg.m-2.Pa-1.s-1',tconda1) c call write_output('tconda2', c &'Taux de condensation CO2 atmospherique /m', c & 'kg.m-3.s-1',tconda2) ! output falling co2 ice in 1st layer: ! call write_output('fallice', ! &'Precipitation of co2 ice', ! & 'kg.m-2.s-1',zfallice(1,1)) #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.).and.(.not.(paleoclimate))) 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 ztcondsol(ngrid)= & 1./(bcond-acond*log(.01*vmr_co2(ngrid,1)* & (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep))) DO islope = 1,nslope pdtsrfc(ngrid,islope)=(ztcondsol(ngrid)- & ztsrf(ngrid,islope))/ptimestep ENDDO ! islope = 1,nslope endif endif #endif END SUBROUTINE co2condens 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 CALL abort_physic('co2condens', & 'top layer is disapearing !', 1) 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 co2condens !!) 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 END SUBROUTINE vl1d c ***************************************************************** SUBROUTINE scavenging_by_co2(ngrid,nlayer,nq,ptimestep,pplev,pq, & rdust,pcondicea,pfallice,pdq_scav,pdqsc) c c c Calcul de la quantite de poussiere lessivee par les nuages de CO2 c c -------------------------------------------------------------------- use tracer_mod, only: nqmx, igcm_h2o_vap, igcm_h2o_ice, & igcm_dust_mass, igcm_dust_number, & igcm_ccn_mass, igcm_ccn_number, & rho_dust, nuice_sed, nuice_ref,r3n_q use comcstfi_h, only: g use dust_param_mod, only: freedust IMPLICIT NONE include "callkeys.h" ! for the flags water and microphys c c c Arguments: 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) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa) REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) REAL,INTENT(IN) :: rdust(ngrid,nlayer) ! dust effective radius REAL,INTENT(IN) :: pcondicea(ngrid,nlayer) ! condensation rate in layer l (kg/m2/s) REAL,INTENT(IN) :: pfallice(ngrid,nlayer+1) ! amount of ice falling from layer l (kg/m2/s) REAL,INTENT(OUT) :: pdq_scav(ngrid,nlayer,nq) ! tendency due to scavenging by co2 REAL,INTENT(OUT) :: pdqsc(ngrid,nq) ! tendency on surface tracers c Locals: INTEGER l,ig REAL scav_ratio_dust, scav_ratio_wice ! ratio of the dust/water ice mass mixing ratios in condensing CO2 ice and in air REAL scav_dust_mass(nlayer+1) ! dust flux (mass) scavenged towards the lower layer (kg/m2/s) (POSITIVE WHEN DOWNWARD) REAL scav_dust_number(nlayer+1) ! dust flux (number) scavenged towards the lower layer (kg/m2/s) (POSITIVE WHEN DOWNWARD) REAL scav_ccn_mass(nlayer+1) ! ccn flux (mass) scavenged towards the lower layer REAL scav_ccn_number(nlayer+1) ! ccn flux (number) scavenged towards the lower layer REAL scav_h2o_ice(nlayer+1) ! water ice flux (mass) scavenged towards the lower layer REAL massl ! mass of the layer l at point ig (kg/m2) c Initialization: scav_ratio_dust = 20 !1 !10 !100 !1000 ! the scavenging ratio value of 20 is a good compromise to remove the dust in the polar night scav_ratio_wice = scav_ratio_dust ! while not drying up the water cycle (which occurs at scav_ratio_wice values above 50 at least) pdq_scav(:,:,:)=0. pdqsc(:,:)=0. DO ig=1,ngrid scav_dust_mass(nlayer+1)=0. scav_dust_number(nlayer+1)=0. scav_ccn_mass(nlayer+1)=0. scav_ccn_number(nlayer+1)=0. scav_h2o_ice(nlayer+1)=0. DO l=nlayer , 1, -1 massl=(pplev(ig,l)-pplev(ig,l+1))/g IF(pcondicea(ig,l).GT.0.)THEN ! if CO2 condenses and traps dust/water ice ! Calculation of the tendencies if (freedust) then pdq_scav(ig,l,igcm_dust_mass)=-pq(ig,l,igcm_dust_mass) & /ptimestep*(1-exp( & -scav_ratio_dust*pcondicea(ig,l)*ptimestep/massl)) pdq_scav(ig,l,igcm_dust_number)=pdq_scav(ig,l,igcm_dust_mass) & *r3n_q/rdust(ig,l) endif if (freedust.AND.microphys) then pdq_scav(ig,l,igcm_ccn_mass)=-pq(ig,l,igcm_ccn_mass) & /ptimestep*(1-exp( & -scav_ratio_wice*pcondicea(ig,l)*ptimestep/massl)) pdq_scav(ig,l,igcm_ccn_number)=pdq_scav(ig,l,igcm_ccn_mass) & *r3n_q/rdust(ig,l) endif if (water) then pdq_scav(ig,l,igcm_h2o_ice)=-pq(ig,l,igcm_h2o_ice) & /ptimestep*(1-exp( & -scav_ratio_wice*pcondicea(ig,l)*ptimestep/massl)) endif ELSE IF(pcondicea(ig,l).LT.0.)THEN ! if CO2 sublimates and releases dust/water ice ! Calculation of the tendencies if (freedust) then pdq_scav(ig,l,igcm_dust_mass)=-pcondicea(ig,l)/massl* & scav_dust_mass(l+1)/pfallice(ig,l+1) pdq_scav(ig,l,igcm_dust_number)=-pcondicea(ig,l)/massl* & scav_dust_number(l+1)/pfallice(ig,l+1) endif if (freedust.AND.microphys) then pdq_scav(ig,l,igcm_ccn_mass)=-pcondicea(ig,l)/massl* & scav_ccn_mass(l+1)/pfallice(ig,l+1) pdq_scav(ig,l,igcm_ccn_number)=-pcondicea(ig,l)/massl* & scav_ccn_number(l+1)/pfallice(ig,l+1) endif if (water) then pdq_scav(ig,l,igcm_h2o_ice)=-pcondicea(ig,l)/massl* & scav_h2o_ice(l+1)/pfallice(ig,l+1) endif END IF ! Calculation of the scavenged dust/wice flux towards the lower layers if (freedust) then scav_dust_mass(l)=-pdq_scav(ig,l,igcm_dust_mass)*massl & +scav_dust_mass(l+1) scav_dust_number(l)=-pdq_scav(ig,l,igcm_dust_number)*massl & +scav_dust_number(l+1) endif if (freedust.AND.microphys) then scav_ccn_mass(l)=-pdq_scav(ig,l,igcm_ccn_mass)*massl & +scav_ccn_mass(l+1) scav_ccn_number(l)=-pdq_scav(ig,l,igcm_ccn_number)*massl & +scav_dust_number(l+1) endif if (water) then scav_h2o_ice(l)=-pdq_scav(ig,l,igcm_h2o_ice)*massl & +scav_h2o_ice(l+1) endif ENDDO ! Calculation of the surface tendencies if (freedust) then pdqsc(ig,igcm_dust_mass)=pdqsc(ig,igcm_dust_mass) & +scav_dust_mass(1) pdqsc(ig,igcm_dust_number)=pdqsc(ig,igcm_dust_number) & +scav_dust_number(1) endif if (freedust.AND.microphys) then pdqsc(ig,igcm_dust_mass)=pdqsc(ig,igcm_dust_mass) & +scav_ccn_mass(1) pdqsc(ig,igcm_dust_number)=pdqsc(ig,igcm_dust_number) & +scav_ccn_number(1) endif if (water) then pdqsc(ig,igcm_h2o_ice)=scav_h2o_ice(1) endif ENDDO ! loop on ngrid END SUBROUTINE scavenging_by_co2 END MODULE co2condens_mod