subroutine condense_n2(ngrid,nlayer,nq,ptimestep, & pcapcal,pplay,pplev,ptsrf,pt, & pdt,pdtsrf,pq,pdq, & pqsurf,pdqsurfc,albedo,pemisurf, & albedo_bareground,albedo_n2_ice_SPECTV, & pdtc,pdtsrfc,pdpsrfc,pdqc) use radinc_h, only : L_NSPECTV use gases_h, only: gfrac, igas_n2 use radii_mod, only : n2_reffrad use aerosol_mod, only : iaero_n2 USE surfdat_h, only: emisice, emissiv USE geometry_mod, only: latitude ! in radians USE tracer_h, only: noms, rho_n2 use comcstfi_mod, only: g, r, cpp implicit none !================================================================== ! Purpose ! ------- ! Condense and/or sublime N2 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_n2_ice_SPECTV(L_NSPECTV) Spectral albedo of N2 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_n2_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 N2 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 picen2(ngrid) ! Updated Surface Tracer (kg/m2). REAL ztcond (ngrid,nlayer) ! Atmospheric Temperatures of condensation of N2. REAL ztnuc (ngrid,nlayer) ! Atmospheric Nucleation Temperatures. REAL ztcondsol(ngrid) ! Temperatures of condensation of N2 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 ppn2 ! !------------------------------------------ ! Saved local variables !------------------------------------------ REAL,SAVE :: latcond=2.5e5 REAL,SAVE :: ccond REAL,SAVE,ALLOCATABLE,DIMENSION(:) :: emisref !$OMP THREADPRIVATE(latcond,ccond,emisref) LOGICAL,SAVE :: firstcall=.true. !$OMP THREADPRIVATE(firstcall) INTEGER,SAVE :: i_n2ice=0 ! n2 ice !$OMP THREADPRIVATE(i_n2ice) CHARACTER(LEN=20) :: tracername ! to temporarily store text !------------------------------------------------ ! Initialization at the first call !------------------------------------------------ IF (firstcall) THEN ! ALLOCATE(emisref(ngrid)) ! ! Find N2 ice tracer. ! do iq=1,nq ! tracername=noms(iq) ! if (tracername.eq."n2_ice") then ! i_n2ice=iq ! endif ! enddo ! write(*,*) "condense_n2: i_n2ice=",i_n2ice ! if((i_n2ice.lt.1))then ! print*,'In condens_cloud but no N2 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 N2 ! if (addn2) then ! m_n2 = 44.01E-3 ! N2 molecular mass (kg/mol) ! m_non2 = 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(in2)/m_n2 + (1-q(in2))/m_non2 ! 1/Mair = A*q(in2) + B ! A = (1/m_n2 - 1/m_non2) ! B = 1/m_non2 ! endif ! Minimum N2 mixing ratio below which mixing occurs with layer above : qn2min =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 N2 Volume mixing ratio ! ------------------------------- ! if (addn2) then ! DO l=1,nlayer ! DO ig=1,ngrid ! qn2=pq(ig,l,in2)+pdq(ig,l,in2)*ptimestep ! Mean air molecular mass = 1/(q(in2)/m_n2 + (1-q(in2))/m_non2) ! mmean=1/(A*qn2 +B) ! vmr_n2(ig,l) = qn2*mmean/m_n2 ! ENDDO ! ENDDO ! else ! DO l=1,nlayer ! DO ig=1,ngrid ! vmr_n2(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 ! ppn2=gfrac(igas_N2)*pplay(ig,l) ! call get_tcond_n2(ppn2,ztcond(ig,l)) ! call get_tnuc_n2(ppn2,ztnuc(ig,l)) ! ENDDO ! ENDDO ! Initialize zq and zt at the beginning of the sub-timestep loop and qsurf. ! DO ig=1,ngrid ! picen2(ig)=pqsurf(ig,i_n2ice) ! DO l=1,nlayer ! zt(ig,l)=pt(ig,l) ! zq(ig,l,i_n2ice)=pq(ig,l,i_n2ice) ! IF( zq(ig,l,i_n2ice).lt.-1.e-6 ) THEN ! print*,'Uh-oh, zq = ',zq(ig,l,i_n2ice),'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_n2ice) = zq(ig,l,i_n2ice) + pdq(ig,l,i_n2ice) * subptimestep ! END DO ! END DO ! Gravitational sedimentation starts. ! Sedimentation computed from radius computed from q in module radii_mod. ! call n2_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_n2) ! !w(ig,ilay,i_n2ice) = 0.0 ! w(ig,ilay,i_n2ice) = 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_n2ice),2.,masse,w(1,1,i_n2ice),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_n2ice) ! END DO !---------------------------------------------------------- ! Condensation / sublimation in the atmosphere !---------------------------------------------------------- ! (MODIFICATIONS FOR EARLY MARS: falling heat neglected, condensation of N2 into tracer i_n2ice) ! 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_n2ice).gt.1.E-10)) THEN ! pdtc(ig,l) = (ztcond(ig,l) - zt(ig,l))/subptimestep ! pdqc(ig,l,i_n2ice) = pdtc(ig,l)*ccond*g ! Case when the ice from above sublimes entirely ! IF ((zq(ig,l,i_n2ice).lt.-pdqc(ig,l,i_n2ice)*subptimestep) & ! .AND. (zq(ig,l,i_n2ice).gt.0)) THEN ! pdqc(ig,l,i_n2ice) = -zq(ig,l,i_n2ice)/subptimestep ! pdtc(ig,l) =-zq(ig,l,i_n2ice)/(ccond*g*subptimestep) ! END IF ! Temperature and q after condensation ! zt(ig,l) = zt(ig,l) + pdtc(ig,l) * subptimestep ! zq(ig,l,i_n2ice) = zq(ig,l,i_n2ice) + pdqc(ig,l,i_n2ice) * 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_n2ice) = & ! (zq(ig,l,i_n2ice)-(pq(ig,l,i_n2ice)+pdq(ig,l,i_n2ice)*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 picen2(ig)=pqsurf(ig,i_n2ice) ppn2=gfrac(igas_N2)*pplay(ig,1) call get_tcond_n2(ppn2,ztcondsol(ig)) ztsrf(ig) = ptsrf(ig) if((ztsrf(ig).le.ztcondsol(ig)+2.0).and.(ngrid.eq.1))then print*,'N2 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 ((picen2(ig)+zfallice(ig)*ptimestep) .NE. 0.))) THEN ! Condensation or partial sublimation of N2 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((picen2(ig)/ptimestep+zfallice(ig)).LE. & -zcondices(ig))THEN zcondices(ig) = -picen2(ig)/ptimestep - zfallice(ig) pdtsrfc(ig)=(latcond/pcapcal(ig))* & (zcondices(ig)) END IF ! Changing N2 ice amount and pressure picen2(ig) = picen2(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_n2' 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.picen2(ig).ge.0.) THEN if(picen2(ig).le.-1.e-8) print*, & 'WARNING : in condense_n2cloud: picen2(',ig,')=', picen2(ig) picen2(ig)=0. endif ! if (picen2(ig) .gt. 1.) then ! N2 Albedo condition changed to ~1 mm coverage. Change by MT2015. ! DO nw=1,L_NSPECTV ! albedo(ig,nw) = albedo_n2_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_n2 !------------------------------------------------------------------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- subroutine get_tcond_n2(p,tcond) ! Calculates the condensation temperature for N2 implicit none real p, tcond IF (p.ge.0.529995) then ! tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB tcond_n2 = (1.)/(0.01583606505-1.211938776e-3*log(7.963685594e-5*p*vmr)) ELSE ! tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT BT tcond_n2 = (1.)/(1./35.6-1.089633028e-3*log(1.968275338*p*vmr)) ENDIF return end subroutine get_tcond_n2 !------------------------------------------------------------------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- !------------------------------------------------------------------------- subroutine get_tnuc_n2(p,tnuc) ! Calculates the nucleation temperature for N2, based on a simple super saturation criterion. JL 2011. use callkeys_mod, only: n2supsat implicit none real p, peff, tnuc real, parameter :: ptriple=518000.0 peff=p/n2supsat !! AF24: this code below is for CO2, needs to be udpated for N2 ! 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_n2