MODULE vdifc_mod IMPLICIT NONE CONTAINS SUBROUTINE vdifc(ngrid,nlay,nsoil,nq,nqsoil,ppopsk, $ ptimestep,pcapcal,lecrit, $ pplay,pplev,pzlay,pzlev,pz0, $ pu,pv,ph,pq,ptsrf,ptsoil,pemis,pqsurf,qsoil, $ pore_icefraction,pdufi,pdvfi,pdhfi, $ pdqfi,pfluxsrf, $ pdudif,pdvdif,pdhdif,pdtsrf,pq2, $ pdqdif,pdqsdif,wstar, $ hfmax,pcondicea_co2microp,sensibFlux, $ dustliftday,local_time,watercap, dwatercap_dif) use tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number, & igcm_dust_submicron, igcm_h2o_vap, & igcm_h2o_ice, alpha_lift, igcm_co2, & igcm_hdo_vap, igcm_hdo_ice, & igcm_stormdust_mass, igcm_stormdust_number use surfdat_h, only: watercaptag, frost_albedo_threshold, dryness, & old_wsublimation_scheme USE comcstfi_h, ONLY: cpp, r, rcp, g, pi use watersat_mod, only: watersat use turb_mod, only: turb_resolved, ustar, tstar use compute_dtau_mod, only: ti_injection_sol,tf_injection_sol use hdo_surfex_mod, only: hdo_surfex c use geometry_mod, only: longitude_deg,latitude_deg ! Joseph use dust_param_mod, only: doubleq, submicron, lifting use write_output_mod, only: write_output use comslope_mod, ONLY: nslope,def_slope,def_slope_mean, & subslope_dist,major_slope,iflat use microphys_h, only: To use paleoclimate_mod, only: d_coef,h2o_ice_depth,lag_layer use comsoil_h, only: layer, mlayer,adsorption_soil use vdif_cd_mod, only: vdif_cd use lmdz_call_atke, only: call_atke use dust_windstress_lift_mod, only: dust_windstress_lift IMPLICIT NONE c======================================================================= c c subject: c -------- c Turbulent diffusion (mixing) for potential T, U, V and tracer c c Shema implicite c On commence par rajouter au variables x la tendance physique c et on resoult en fait: c x(t+1) = x(t) + dt * (dx/dt)phys(t) + dt * (dx/dt)difv(t+1) c c author: c ------ c Hourdin/Forget/Fournier c======================================================================= c----------------------------------------------------------------------- c declarations: c ------------- include "callkeys.h" c c arguments: c ---------- INTEGER,INTENT(IN) :: ngrid,nlay,nsoil,nqsoil REAL,INTENT(IN) :: ptimestep REAL,INTENT(IN) :: pplay(ngrid,nlay),pplev(ngrid,nlay+1) REAL,INTENT(IN) :: pzlay(ngrid,nlay),pzlev(ngrid,nlay+1) REAL,INTENT(IN) :: pu(ngrid,nlay),pv(ngrid,nlay) REAL,INTENT(IN) :: ph(ngrid,nlay) REAL,INTENT(IN) :: ptsrf(ngrid,nslope),pemis(ngrid,nslope) REAL,INTENT(IN) :: pdufi(ngrid,nlay),pdvfi(ngrid,nlay) REAL,INTENT(IN) :: pdhfi(ngrid,nlay) REAL,INTENT(IN) :: pfluxsrf(ngrid,nslope) REAL,INTENT(OUT) :: pdudif(ngrid,nlay),pdvdif(ngrid,nlay) REAL,INTENT(OUT) :: pdtsrf(ngrid,nslope),pdhdif(ngrid,nlay) REAL,INTENT(IN) :: pcapcal(ngrid,nslope) REAL,INTENT(INOUT) :: pq2(ngrid,nlay+1) REAL,INTENT(IN) :: ptsoil(ngrid,nsoil,nslope) c Argument added for condensation: REAL,INTENT(IN) :: ppopsk(ngrid,nlay) logical,INTENT(IN) :: lecrit REAL,INTENT(IN) :: pcondicea_co2microp(ngrid,nlay)! tendency due to CO2 condensation (kg/kg.s-1) REAL,INTENT(IN) :: pz0(ngrid) ! surface roughness length (m) c Argument added to account for subgrid gustiness : REAL,INTENT(IN) :: wstar(ngrid), hfmax(ngrid)!, zi(ngrid) c Traceurs : integer,intent(in) :: nq REAL,INTENT(IN) :: pqsurf(ngrid,nq,nslope) REAL :: zqsurf(ngrid) ! temporary water tracer real,intent(in) :: pq(ngrid,nlay,nq), pdqfi(ngrid,nlay,nq) real,intent(out) :: pdqdif(ngrid,nlay,nq) real,intent(out) :: pdqsdif(ngrid,nq,nslope) REAL,INTENT(in) :: dustliftday(ngrid) REAL,INTENT(inout) :: qsoil(ngrid,nsoil,nqsoil,nslope) !subsurface tracers REAL,INTENT(in) :: local_time(ngrid) c local: c ------ REAL :: pt(ngrid,nlay) INTEGER ilev,ig,ilay,nlev,islope,ik,lice REAL z4st,zdplanck(ngrid) REAL zkv(ngrid,nlay+1),zkh(ngrid,nlay+1) REAL zkq(ngrid,nlay+1) REAL zcdv(ngrid,nslope),zcdh(ngrid,nslope) REAL :: zcdv_true(ngrid,nslope) REAL :: zcdh_true(ngrid,nslope) ! drag coeff are used by the LES to recompute u* and hfx REAL :: zcdv_tmp(ngrid),zcdh_tmp(ngrid) ! drag coeffs for the major sub-grid surface REAL :: zcdv_true_tmp(ngrid),zcdh_true_tmp(ngrid) ! drag coeffs (computed with wind gustiness for the major sub-grid surface REAL zu(ngrid,nlay),zv(ngrid,nlay) REAL zh(ngrid,nlay) REAL ztsrf2(ngrid) REAL z1(ngrid),z2(ngrid) REAL za(ngrid,nlay),zb(ngrid,nlay) REAL zb0(ngrid,nlay) REAL zc(ngrid,nlay),zd(ngrid,nlay) REAL zcst1 REAL zu2(ngrid) EXTERNAL SSUM,SCOPY REAL SSUM LOGICAL,SAVE :: firstcall=.true. !$OMP THREADPRIVATE(firstcall) c variable added for CO2 condensation: c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ REAL hh , zhcond(ngrid,nlay) REAL,PARAMETER :: latcond=5.9e5 REAL,PARAMETER :: tcond1mb=136.27 REAL,SAVE :: acond,bcond !$OMP THREADPRIVATE(acond,bcond) c Subtimestep & implicit treatment of water vapor c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ REAL zdqsdif_surf(ngrid) ! subtimestep pdqsdif for water ice REAL ztsrf(ngrid) ! temporary surface temperature in tsub REAL zdtsrf(ngrid,nslope) ! surface temperature tendancy in tsub REAL surf_h2o_lh(ngrid,nslope) ! Surface h2o latent heat flux REAL zsurf_h2o_lh(ngrid,nslope) ! Tsub surface h2o latent heat flux c For latent heat release from ground water ice sublimation c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ REAL tsrf_lh(ngrid) ! temporary surface temperature with lh effect REAL lh ! latent heat, formulation given in the Technical Document: ! "Modeling water ice sublimation under Phoenix-like conditions", Montmessin et al. 2004 c Tracers : c ~~~~~~~ INTEGER iq REAL zq(ngrid,nlay,nq) REAL zq1temp(ngrid) REAL rho(ngrid) ! near surface air density REAL qsat(ngrid) REAL hdoflux(ngrid,nslope) ! value of vapour flux of HDO REAL hdoflux_meshavg(ngrid) ! value of vapour flux of HDO ! REAL h2oflux(ngrid) ! value of vapour flux of H2O REAL old_h2o_vap(ngrid) ! traceur d'eau avant traitement REAL saved_h2o_vap(ngrid) ! traceur d'eau avant traitement REAL kmixmin c Argument added for surface water ice budget: REAL,INTENT(IN) :: watercap(ngrid,nslope) REAL,INTENT(OUT) :: dwatercap_dif(ngrid,nslope) c Subtimestep to compute h2o latent heat flux: REAL :: dtmax = 0.5 ! subtimestep temp criterion INTEGER tsub ! adaptative subtimestep (seconds) REAL subtimestep !ptimestep/nsubtimestep INTEGER nsubtimestep(ngrid) ! number of subtimestep (int) c Mass-variation scheme : c ~~~~~~~ INTEGER j,l REAL zcondicea(ngrid,nlay) REAL zt(ngrid,nlay),ztcond(ngrid,nlay+1) REAL betam(ngrid,nlay),dmice(ngrid,nlay) REAL pdtc(ngrid,nlay) REAL zhs(ngrid,nlay) REAL,SAVE :: ccond !$OMP THREADPRIVATE(ccond) c Theta_m formulation for mass-variation scheme : c ~~~~~~~ INTEGER,SAVE :: ico2 INTEGER llnt(ngrid) REAL,SAVE :: m_co2, m_noco2, A , B REAL vmr_co2(ngrid,nlay) REAL qco2,mmean(ngrid,nlay) !$OMP THREADPRIVATE(ico2,m_co2,m_noco2,A,B) REAL,INTENT(OUT) :: sensibFlux(ngrid) !!MARGAUX REAL DoH_vap(ngrid,nlay) !! Sub-grid scale slopes REAL :: pdqsdif_tmp(ngrid,nq) ! Temporary for dust lifting REAL :: watercap_tmp(ngrid) REAL :: zq_slope_vap(ngrid,nlay,nq,nslope) REAL :: zq_tmp_vap(ngrid,nlay,nq) REAL :: ptsrf_tmp(ngrid) REAL :: pqsurf_tmp(ngrid,nq) REAL :: pdqsdif_tmphdo(ngrid,nq) REAL :: qsat_tmp(ngrid) INTEGER :: indmax character*2 str2 !! Subsurface exchanges LOGICAL :: exchange ! boolean to check if exchange between the subsurface and the atmosphere can occurs REAL :: zdqsdif_regolith(ngrid,nslope) ! Flux from subsurface (positive pointing outwards) (kg/m^2/s) REAL zq1temp_regolith(ngrid) ! Temporary atmospheric mixing ratio after exchange with subsurface (kg / kg) REAL zdqsdif_tot(ngrid) ! subtimestep pdqsdif for water ice LOGICAL :: writeoutput ! boolean to say to soilexchange.F if we are at the last iteration and thus if he can write in the diagsoil REAL, INTENT(out) :: pore_icefraction(ngrid,nsoil,nslope) ! ice filling fraction in the pores !! Subsurface ice interactions REAL Tice(ngrid,nslope) ! subsurface temperature where ice is located (K) REAL qsat_ssi(ngrid,nslope) ! qsat(Tice) (kg/kg) REAL resist(ngrid,nslope) ! subsurface ice flux reduction coef (1) REAL zdqsdif_ssi_atm(ngrid,nslope) ! SSI - atmosphere flux (kg/m^2/s^-1) at a given subtimestep REAL zdqsdif_ssi_frost(ngrid,nslope) ! SSI - frost flux (kg/m^2/s^-1) at a given subtimestep REAL zdqsdif_ssi_atm_tot(ngrid,nslope) ! SSI - atmosphere flux (kg/m^2/s^-1) summed over all subtimestep REAL zdqsdif_ssi_frost_tot(ngrid,nslope) ! SSI - frost flux (kg/m^2/s^-1) summed over all subtimestep REAL zdqsdif_ssi_tot(ngrid,nslope) ! Total flux of the interactions with SSI kg/m^2/s^-1) REAL :: tol_frost = 2e-3 ! tolerence for frost thicnkess (kg/m^2) to avoid numerical noise effect c ** un petit test de coherence c -------------------------- ! AS: OK firstcall absolute IF (firstcall) THEN c To compute: Tcond= 1./(bcond-acond*log(.0095*p)) (p in pascal) bcond=1./tcond1mb acond=r/latcond ccond=cpp/(g*latcond) PRINT*,'In vdifc: 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 firstcall=.false. ENDIF DO ig = 1,ngrid ptsrf_tmp(ig) = ptsrf(ig,major_slope(ig)) pqsurf_tmp(ig,:) = pqsurf(ig,:,major_slope(ig)) ENDDO c----------------------------------------------------------------------- c 1. initialisation c ----------------- nlev=nlay+1 ! initialize output tendencies to zero: pdudif(1:ngrid,1:nlay)=0 pdvdif(1:ngrid,1:nlay)=0 pdhdif(1:ngrid,1:nlay)=0 pdtsrf(1:ngrid,1:nslope)=0 zdtsrf(1:ngrid,1:nslope)=0 surf_h2o_lh(1:ngrid,1:nslope)=0 zsurf_h2o_lh(1:ngrid,1:nslope)=0 pdqdif(1:ngrid,1:nlay,1:nq)=0 pdqsdif(1:ngrid,1:nq,1:nslope)=0 pdqsdif_tmp(1:ngrid,1:nq)=0 zdqsdif_surf(1:ngrid)=0 dwatercap_dif(1:ngrid,1:nslope)=0 zdqsdif_regolith(1:ngrid,1:nslope)=0 zq1temp_regolith(1:ngrid)=0 zdqsdif_tot(1:ngrid)=0 pore_icefraction(:,:,:) = 0. zdqsdif_ssi_atm(:,:) = 0. zdqsdif_ssi_frost(:,:) = 0. zdqsdif_ssi_tot(:,:) = 0. zdqsdif_ssi_atm_tot(:,:) = 0. zdqsdif_ssi_frost_tot(:,:) = 0. c ** calcul de rho*dz et dt*rho/dz=dt*rho**2 g/dp c avec rho=p/RT=p/ (R Theta) (p/ps)**kappa c ---------------------------------------- DO ilay=1,nlay DO ig=1,ngrid za(ig,ilay)=(pplev(ig,ilay)-pplev(ig,ilay+1))/g ! Mass variation scheme: betam(ig,ilay)=-za(ig,ilay)*latcond/(cpp*ppopsk(ig,ilay)) ENDDO ENDDO zcst1=4.*g*ptimestep/(r*r) DO ilev=2,nlev-1 DO ig=1,ngrid zb0(ig,ilev)=pplev(ig,ilev)* s (pplev(ig,1)/pplev(ig,ilev))**rcp / s (ph(ig,ilev-1)+ph(ig,ilev)) zb0(ig,ilev)=zcst1*zb0(ig,ilev)*zb0(ig,ilev)/ s (pplay(ig,ilev-1)-pplay(ig,ilev)) ENDDO ENDDO DO ig=1,ngrid zb0(ig,1)=ptimestep*pplev(ig,1)/(r*ptsrf_tmp(ig)) ENDDO c ** diagnostique pour l'initialisation c ---------------------------------- IF(lecrit) THEN ig=ngrid/2+1 PRINT*,'Pression (mbar) ,altitude (km),u,v,theta, rho dz' DO ilay=1,nlay WRITE(*,'(6f11.5)') s .01*pplay(ig,ilay),.001*pzlay(ig,ilay), s pu(ig,ilay),pv(ig,ilay),ph(ig,ilay),za(ig,ilay) ENDDO PRINT*,'Pression (mbar) ,altitude (km),zb' DO ilev=1,nlay WRITE(*,'(3f15.7)') s .01*pplev(ig,ilev),.001*pzlev(ig,ilev), s zb0(ig,ilev) ENDDO ENDIF c ----------------------------------- c Potential Condensation temperature: c ----------------------------------- c Compute CO2 Volume mixing ratio c ------------------------------- if (callcond.and.(ico2.ne.0)) then DO ilev=1,nlay DO ig=1,ngrid qco2=MAX(1.E-30 & ,pq(ig,ilev,ico2)+pdqfi(ig,ilev,ico2)*ptimestep) c Mean air molecular mass = 1/(q(ico2)/m_co2 + (1-q(ico2))/m_noco2) mmean(ig,ilev)=1/(A*qco2 +B) vmr_co2(ig,ilev) = qco2*mmean(ig,ilev)/m_co2 ENDDO ENDDO else DO ilev=1,nlay DO ig=1,ngrid vmr_co2(ig,ilev)=0.95 ENDDO ENDDO end if c forecast of atmospheric temperature zt and frost temperature ztcond c -------------------------------------------------------------------- if (callcond) then DO ilev=1,nlay DO ig=1,ngrid ztcond(ig,ilev)= & 1./(bcond-acond*log(.01*vmr_co2(ig,ilev)*pplay(ig,ilev))) if (pplay(ig,ilev).lt.1e-4) ztcond(ig,ilev)=0.0 !mars Monica ! zhcond(ig,ilev) = ! & (1./(bcond-acond*log(.0095*pplay(ig,ilev))))/ppopsk(ig,ilev) zhcond(ig,ilev) = ztcond(ig,ilev)/ppopsk(ig,ilev) END DO END DO ztcond(:,nlay+1)=ztcond(:,nlay) else zhcond(:,:) = 0 ztcond(:,:) = 0 end if c----------------------------------------------------------------------- c 2. ajout des tendances physiques c ----------------------------- DO ilev=1,nlay DO ig=1,ngrid zu(ig,ilev)=pu(ig,ilev)+pdufi(ig,ilev)*ptimestep zv(ig,ilev)=pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep zh(ig,ilev)=ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep ! zh(ig,ilev)=max(zh(ig,ilev),zhcond(ig,ilev)) ENDDO ENDDO zq(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq)+ & pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep c----------------------------------------------------------------------- c 3. schema de turbulence c -------------------- c ** source d'energie cinetique turbulente a la surface c (condition aux limites du schema de diffusion turbulente c dans la couche limite c --------------------- CALL vdif_cd(ngrid,nlay,nslope,pz0,g,pzlay,pplay,pu,pv,wstar, & ptsrf,ph,mmean(:,1),zq(:,:,igcm_h2o_vap), & pqsurf(:,igcm_h2o_ice,:), .false., & zcdv_true,zcdh_true) zu2(:)=pu(:,1)*pu(:,1)+pv(:,1)*pv(:,1) DO islope = 1,nslope IF (callrichsl) THEN zcdv(:,islope)=zcdv_true(:,islope)*sqrt(zu2(:)+ & (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2) zcdh(:,islope)=zcdh_true(:,islope)*sqrt(zu2(:)+ & (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2) ELSE zcdv(:,islope)=zcdv_true(:,islope)*sqrt(zu2(:)) ! 1 / bulk aerodynamic momentum conductance zcdh(:,islope)=zcdh_true(:,islope)*sqrt(zu2(:)) ! 1 / bulk aerodynamic heat conductance ENDIF ENDDO ustar(:) = 0 tstar(:) = 0 DO ig = 1,ngrid zcdv_tmp(ig) = zcdv(ig,major_slope(ig)) zcdh_tmp(ig) = zcdh(ig,major_slope(ig)) zcdv_true_tmp(ig) = zcdv_true(ig,major_slope(ig)) zcdh_true_tmp(ig) = zcdh_true(ig,major_slope(ig)) IF (callrichsl) THEN ustar(ig)=sqrt(zcdv_true(ig,major_slope(ig))) & *sqrt(zu2(ig)+(log(1.+0.7*wstar(ig) + & 2.3*wstar(ig)**2))**2) IF (zcdh_true(ig,major_slope(ig)) .ne. 0.) THEN ! When Cd=Ch=0, u*=t*=0 tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig)) & *zcdh_tmp(ig)/ustar(ig) ENDIF ELSE ustar(ig)=sqrt(zcdv_true(ig,major_slope(ig))) & *sqrt(zu2(ig)) tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig)) & *zcdh_true(ig,major_slope(ig)) & /sqrt(zcdv_true(ig,major_slope(ig))) ENDIF ENDDO c ** schema de diffusion turbulente dans la couche limite c ---------------------------------------------------- pt(:,:)=ph(:,:)*ppopsk(:,:) if (callyamada4) then call yamada4(ngrid,nlay,nq,ptimestep,g,r,pplev,pt s ,pzlev,pzlay,pu,pv,ph,pq,zcdv_true_tmp,pq2,zkv,zkh,zkq,ustar s ,9) elseif (callatke) then call call_atke(ptimestep,ngrid,nlay,zcdv_true_tmp, s zcdh_true_tmp,pu(:,1),pv(:,1),ptsrf_tmp, s pu,pv,pt,zq(:,1,igcm_h2o_vap),pplay,pplev, s pzlay,pzlev,pq2,zkv(:,1:nlay),zkh(:,1:nlay)) zkv(:,nlay+1) = zkv(:,nlay) zkh(:,nlay+1) = zkh(:,nlay) else call vdif_kc(ngrid,nlay,nq,ptimestep,g,pzlev,pzlay s ,pu,pv,ph,zcdv_true_tmp s ,pq2,zkv,zkh,zq) endif if ((doubleq).and.(ngrid.eq.1)) then kmixmin = 80. !80.! minimum eddy mix coeff in 1D do ilev=2,nlay do ig=1,ngrid zkh(ig,ilev) = max(kmixmin,zkh(ig,ilev)) zkv(ig,ilev) = max(kmixmin,zkv(ig,ilev)) end do end do end if c ** diagnostique pour le schema de turbulence c ----------------------------------------- IF(lecrit) THEN PRINT* PRINT*,'Diagnostic for the vertical turbulent mixing' PRINT*,'Cd for momentum and potential temperature' PRINT*,zcdv_tmp(ngrid/2+1),zcdh_tmp(ngrid/2+1) PRINT*,'Mixing coefficient for momentum and pot.temp.' DO ilev=1,nlay PRINT*,zkv(ngrid/2+1,ilev),zkh(ngrid/2+1,ilev) ENDDO ENDIF c----------------------------------------------------------------------- c 4. inversion pour l'implicite sur u c -------------------------------- c ** l'equation est c u(t+1) = u(t) + dt * {(du/dt)phys}(t) + dt * {(du/dt)difv}(t+1) c avec c /zu/ = u(t) + dt * {(du/dt)phys}(t) (voir paragraphe 2.) c et c dt * {(du/dt)difv}(t+1) = dt * {(d/dz)[ Ku (du/dz) ]}(t+1) c donc les entrees sont /zcdv/ pour la condition a la limite sol c et /zkv/ = Ku zb(1:ngrid,2:nlay)=zkv(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay) zb(1:ngrid,1)=zcdv_tmp(1:ngrid)*zb0(1:ngrid,1) DO ig=1,ngrid z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) zc(ig,nlay)=za(ig,nlay)*zu(ig,nlay)*z1(ig) zd(ig,nlay)=zb(ig,nlay)*z1(ig) ENDDO DO ilay=nlay-1,1,-1 DO ig=1,ngrid z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) zc(ig,ilay)=(za(ig,ilay)*zu(ig,ilay)+ $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) zd(ig,ilay)=zb(ig,ilay)*z1(ig) ENDDO ENDDO DO ig=1,ngrid zu(ig,1)=zc(ig,1) ENDDO DO ilay=2,nlay DO ig=1,ngrid zu(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zu(ig,ilay-1) ENDDO ENDDO c----------------------------------------------------------------------- c 5. inversion pour l'implicite sur v c -------------------------------- c ** l'equation est c v(t+1) = v(t) + dt * {(dv/dt)phys}(t) + dt * {(dv/dt)difv}(t+1) c avec c /zv/ = v(t) + dt * {(dv/dt)phys}(t) (voir paragraphe 2.) c et c dt * {(dv/dt)difv}(t+1) = dt * {(d/dz)[ Kv (dv/dz) ]}(t+1) c donc les entrees sont /zcdv/ pour la condition a la limite sol c et /zkv/ = Kv DO ig=1,ngrid z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) zc(ig,nlay)=za(ig,nlay)*zv(ig,nlay)*z1(ig) zd(ig,nlay)=zb(ig,nlay)*z1(ig) ENDDO DO ilay=nlay-1,1,-1 DO ig=1,ngrid z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) zc(ig,ilay)=(za(ig,ilay)*zv(ig,ilay)+ $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) zd(ig,ilay)=zb(ig,ilay)*z1(ig) ENDDO ENDDO DO ig=1,ngrid zv(ig,1)=zc(ig,1) ENDDO DO ilay=2,nlay DO ig=1,ngrid zv(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zv(ig,ilay-1) ENDDO ENDDO c----------------------------------------------------------------------- c Using the wind modified by friction for lifting and sublimation c ---------------------------------------------------------------- ! This is computed above and takes into account surface-atmosphere flux ! enhancement by subgrid gustiness and atmospheric-stability related ! variations of transfer coefficients. ! Calculate Cd again with wind slowed by friction c ------------------------------------------- CALL vdif_cd(ngrid,nlay,nslope,pz0,g,pzlay,pplay,zu,zv,wstar, & ptsrf,ph,mmean(:,1),zq(:,:,igcm_h2o_vap), & pqsurf(:,igcm_h2o_ice,:), .true., & zcdv_true,zcdh_true) zu2(:)=zu(:,1)*zu(:,1)+zv(:,1)*zv(:,1) DO islope = 1,nslope IF (callrichsl) THEN zcdv(:,islope)=zcdv_true(:,islope)*sqrt(zu2(:)+ & (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2) zcdh(:,islope)=zcdh_true(:,islope)*sqrt(zu2(:)+ & (log(1.+0.7*wstar(:) + 2.3*wstar(:)**2))**2) ELSE zcdv(:,islope)=zcdv_true(:,islope)*sqrt(zu2(:)) ! 1 / bulk aerodynamic momentum conductance zcdh(:,islope)=zcdh_true(:,islope)*sqrt(zu2(:)) ! 1 / bulk aerodynamic heat conductance ENDIF ENDDO ustar(:) = 0 tstar(:) = 0 DO ig = 1,ngrid zcdv_tmp(ig) = zcdv(ig,major_slope(ig)) zcdh_tmp(ig) = zcdh(ig,major_slope(ig)) zcdv_true_tmp(ig) = zcdv_true(ig,major_slope(ig)) zcdh_true_tmp(ig) = zcdh_true(ig,major_slope(ig)) IF (callrichsl) THEN ustar(ig)=sqrt(zcdv_true(ig,major_slope(ig))) & *sqrt(zu2(ig)+(log(1.+0.7*wstar(ig) + & 2.3*wstar(ig)**2))**2) IF (zcdh_true(ig,major_slope(ig)) .ne. 0.) THEN ! When Cd=Ch=0, u*=t*=0 tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig)) & *zcdh_tmp(ig)/ustar(ig) ENDIF ELSE ustar(ig)=sqrt(zcdv_true(ig,major_slope(ig))) & *sqrt(zu2(ig)) tstar(ig)=(ph(ig,1)-ptsrf_tmp(ig)) & *zcdh_true(ig,major_slope(ig)) & /sqrt(zcdv_true(ig,major_slope(ig))) ENDIF ENDDO c----------------------------------------------------------------------- c 6. inversion pour l'implicite sur h sans oublier le couplage c avec le sol (conduction) c ------------------------ c ** l'equation est c h(t+1) = h(t) + dt * {(dh/dt)phys}(t) + dt * {(dh/dt)difv}(t+1) c avec c /zh/ = h(t) + dt * {(dh/dt)phys}(t) (voir paragraphe 2.) c et c dt * {(dh/dt)difv}(t+1) = dt * {(d/dz)[ Kh (dh/dz) ]}(t+1) c donc les entrees sont /zcdh/ pour la condition de raccord au sol c et /zkh/ = Kh c ------------- c Mass variation scheme: zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay) zb(1:ngrid,1)=zcdh_tmp(1:ngrid)*zb0(1:ngrid,1) c on initialise dm c pdtc(:,:)=0. zt(:,:)=0. dmice(:,:)=0. c ** calcul de (d Planck / dT) a la temperature d'interface c ------------------------------------------------------ z4st=4.*5.67e-8*ptimestep IF (tke_heat_flux .eq. 0.) THEN DO ig=1,ngrid indmax = major_slope(ig) zdplanck(ig)=z4st*pemis(ig,indmax)*ptsrf(ig,indmax)* & ptsrf(ig,indmax)*ptsrf(ig,indmax) ENDDO ELSE zdplanck(:)=0. ENDIF ! calcul de zc et zd pour la couche top en prenant en compte le terme ! de variation de masse (on fait une boucle pour que \E7a converge) ! Identification des points de grilles qui ont besoin de la correction llnt(:)=1 IF (.not.turb_resolved) THEN IF (callcond) THEN DO ig=1,ngrid DO l=1,nlay if(zh(ig,l) .lt. zhcond(ig,l)) then llnt(ig)=300 ! 200 and 100 do not go beyond month 9 with normal dissipation goto 5 endif ENDDO 5 continue ENDDO ENDIF ENDIF DO ig=1,ngrid indmax = major_slope(ig) ! Initialization of z1 and zd, which do not depend on dmice z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) zd(ig,nlay)=zb(ig,nlay)*z1(ig) DO ilay=nlay-1,1,-1 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) zd(ig,ilay)=zb(ig,ilay)*z1(ig) ENDDO ! Convergence loop DO j=1,llnt(ig) z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) zc(ig,nlay)=za(ig,nlay)*zh(ig,nlay) & -betam(ig,nlay)*dmice(ig,nlay) zc(ig,nlay)=zc(ig,nlay)*z1(ig) ! zd(ig,nlay)=zb(ig,nlay)*z1(ig) ! calcul de zc et zd pour les couches du haut vers le bas DO ilay=nlay-1,1,-1 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) zc(ig,ilay)=(za(ig,ilay)*zh(ig,ilay)+ $ zb(ig,ilay+1)*zc(ig,ilay+1)- $ betam(ig,ilay)*dmice(ig,ilay))*z1(ig) ! zd(ig,ilay)=zb(ig,ilay)*z1(ig) ENDDO c ** calcul de la temperature_d'interface et de sa tendance. c on ecrit que la somme des flux est nulle a l'interface c a t + \delta t, c c'est a dire le flux radiatif a {t + \delta t} c + le flux turbulent a {t + \delta t} c qui s'ecrit K (T1-Tsurf) avec T1 = d1 Tsurf + c1 c (notation K dt = /cpp*b/) c + le flux dans le sol a t c + l'evolution du flux dans le sol lorsque la temperature d'interface c passe de sa valeur a t a sa valeur a {t + \delta t}. c ---------------------------------------------------- z1(ig)=pcapcal(ig,indmax)*ptsrf(ig,indmax) s + cpp*zb(ig,1)*zc(ig,1) s + zdplanck(ig)*ptsrf(ig,indmax) s + pfluxsrf(ig,indmax)*ptimestep z2(ig)= pcapcal(ig,indmax)+cpp*zb(ig,1)*(1.-zd(ig,1)) s +zdplanck(ig) ztsrf2(ig)=z1(ig)/z2(ig) ! pdtsrf(ig)=(ztsrf2(ig)-ptsrf(ig))/ptimestep !incremented outside loop zhs(ig,1)=zc(ig,1)+zd(ig,1)*ztsrf2(ig) c ** et a partir de la temperature au sol on remonte c ----------------------------------------------- DO ilay=2,nlay zhs(ig,ilay)=zc(ig,ilay)+zd(ig,ilay)*zhs(ig,ilay-1) ENDDO DO ilay=1,nlay zt(ig,ilay)=zhs(ig,ilay)*ppopsk(ig,ilay) ENDDO c Condensation/sublimation in the atmosphere c ------------------------------------------ c (computation of zcondicea and dmice) IF (.NOT. co2clouds) then DO l=nlay , 1, -1 IF(zt(ig,l).LT.ztcond(ig,l)) THEN pdtc(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1)) & *ccond*pdtc(ig,l) dmice(ig,l)= dmice(ig,l) + zcondicea(ig,l)*ptimestep END IF ENDDO ELSE DO l=nlay , 1, -1 zcondicea(ig,l)= 0.!pcondicea_co2microp(ig,l)* c & (pplev(ig,l) - pplev(ig,l+1))/g dmice(ig,l)= 0.!dmice(ig,l) + zcondicea(ig,l)*ptimestep pdtc(ig,l)=0. ENDDO ENDIF ENDDO!of Do j=1,XXX pdtsrf(ig,indmax)=(ztsrf2(ig)-ptsrf(ig,indmax))/ptimestep ENDDO !of Do ig=1,ngrid DO ig=1,ngrid ! computing sensible heat flux (atm => surface) sensibFlux(ig)=cpp*zb(ig,1)/ptimestep*(zhs(ig,1)-ztsrf2(ig)) ENDDO c Now implicit sheme on each sub-grid subslope: IF (nslope.ne.1) then DO islope=1,nslope DO ig=1,ngrid IF(islope.ne.major_slope(ig)) then IF (tke_heat_flux .eq. 0.) THEN zdplanck(ig)=z4st*pemis(ig,islope)*ptsrf(ig,islope)**3 ELSE zdplanck(ig) = 0. ENDIF zb(ig,1)=zcdh(ig,islope)*zb0(ig,1) z1(ig)=pcapcal(ig,islope)*ptsrf(ig,islope) s + cpp*zb(ig,1)*zc(ig,1) s + zdplanck(ig)*ptsrf(ig,islope) s + pfluxsrf(ig,islope)*ptimestep z2(ig)= pcapcal(ig,islope)+cpp*zb(ig,1)*(1.-zd(ig,1)) s +zdplanck(ig) ztsrf2(ig)=z1(ig)/z2(ig) pdtsrf(ig,islope)=(ztsrf2(ig)-ptsrf(ig,islope))/ptimestep ENDIF ! islope != indmax ENDDO ! ig ENDDO !islope ENDIF !nslope.ne.1 c----------------------------------------------------------------------- c TRACERS c ------- c Calcul du flux vertical au bas de la premiere couche (dust) : c ----------------------------------------------------------- do ig=1,ngrid rho(ig) = zb0(ig,1) /ptimestep c zb(ig,1) = 0. end do c Dust lifting: if (lifting) then #ifndef MESOSCALE if (doubleq.AND.submicron) then do ig=1,ngrid c if(qsurf(ig,igcm_co2).lt.1) then pdqsdif_tmp(ig,igcm_dust_mass) = & -alpha_lift(igcm_dust_mass) pdqsdif_tmp(ig,igcm_dust_number) = & -alpha_lift(igcm_dust_number) pdqsdif_tmp(ig,igcm_dust_submicron) = & -alpha_lift(igcm_dust_submicron) c end if end do else if (doubleq) then if (dustinjection.eq.0) then !injection scheme 0 (old) !or 2 (injection in CL) do ig=1,ngrid if(pqsurf_tmp(ig,igcm_co2).lt.1) then ! pas de soulevement si glace CO2 pdqsdif_tmp(ig,igcm_dust_mass) = & -alpha_lift(igcm_dust_mass) pdqsdif_tmp(ig,igcm_dust_number) = & -alpha_lift(igcm_dust_number) end if end do elseif(dustinjection.eq.1)then ! dust injection scheme = 1 injection from surface do ig=1,ngrid if(pqsurf_tmp(ig,igcm_co2).lt.1) then ! pas de soulevement si glace CO2 IF((ti_injection_sol.LE.local_time(ig)).and. & (local_time(ig).LE.tf_injection_sol)) THEN if (rdstorm) then !Rocket dust storm scheme pdqsdif_tmp(ig,igcm_stormdust_mass) = & -alpha_lift(igcm_stormdust_mass) & *dustliftday(ig) pdqsdif_tmp(ig,igcm_stormdust_number) = & -alpha_lift(igcm_stormdust_number) & *dustliftday(ig) pdqsdif_tmp(ig,igcm_dust_mass)= 0. pdqsdif_tmp(ig,igcm_dust_number)= 0. else pdqsdif_tmp(ig,igcm_dust_mass)= & -dustliftday(ig)* & alpha_lift(igcm_dust_mass) pdqsdif_tmp(ig,igcm_dust_number)= & -dustliftday(ig)* & alpha_lift(igcm_dust_number) endif if (submicron) then pdqsdif_tmp(ig,igcm_dust_submicron) = 0. endif ! if (submicron) ELSE ! outside dust injection time frame pdqsdif_tmp(ig,igcm_dust_mass)= 0. pdqsdif_tmp(ig,igcm_dust_number)= 0. if (rdstorm) then pdqsdif_tmp(ig,igcm_stormdust_mass)= 0. pdqsdif_tmp(ig,igcm_stormdust_number)= 0. end if ENDIF end if ! of if(qsurf(ig,igcm_co2).lt.1) end do endif ! end if dustinjection else if (submicron) then do ig=1,ngrid pdqsdif_tmp(ig,igcm_dust_submicron) = & -alpha_lift(igcm_dust_submicron) end do else #endif call dust_windstress_lift(ngrid,nlay,nq,rho, & zcdh_true_tmp,zcdh_tmp, & pqsurf_tmp(:,igcm_co2),pdqsdif_tmp) #ifndef MESOSCALE endif !doubleq.AND.submicron #endif else pdqsdif_tmp(1:ngrid,1:nq) = 0. end if c OU calcul de la valeur de q a la surface (water) : c ---------------------------------------- c Inversion pour l'implicite sur q c Cas des traceurs qui ne sont pas h2o_vap c h2o_vap est traite plus loin avec un sous pas de temps c hdo_vap est traite ensuite car dependant de h2o_vap c -------------------------------- do iq=1,nq !for all tracers except water vapor if ((.not. water).or.(.not. iq.eq.igcm_h2o_vap).or. & (.not. iq.eq.igcm_hdo_vap)) then zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay) zb(1:ngrid,1)=0 DO ig=1,ngrid z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig) zd(ig,nlay)=zb(ig,nlay)*z1(ig) ENDDO DO ilay=nlay-1,2,-1 DO ig=1,ngrid z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+ $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) zd(ig,ilay)=zb(ig,ilay)*z1(ig) ENDDO ENDDO if ((iq.eq.igcm_h2o_ice) $ .or. (hdo.and.(iq.eq.igcm_hdo_ice) )) then DO ig=1,ngrid z1(ig)=1./(za(ig,1)+zb(ig,1)+ $ zb(ig,2)*(1.-zd(ig,2))) zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ $ zb(ig,2)*zc(ig,2)) *z1(ig) !special case h2o_ice ENDDO else ! every other tracer DO ig=1,ngrid z1(ig)=1./(za(ig,1)+zb(ig,1)+ $ zb(ig,2)*(1.-zd(ig,2))) zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ $ zb(ig,2)*zc(ig,2) + $ (-pdqsdif_tmp(ig,iq)) *ptimestep) *z1(ig) !tracer flux from surface ENDDO endif !((iq.eq.igcm_h2o_ice) c Starting upward calculations for simple mixing of tracer (dust) DO ig=1,ngrid zq(ig,1,iq)=zc(ig,1) DO ilay=2,nlay zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq) ENDDO ENDDO DO islope = 1,nslope DO ig = 1,ngrid pdqsdif(ig,iq,islope) = pdqsdif_tmp(ig,iq) & * cos(pi*def_slope_mean(islope)/180.) ENDDO ENDDO endif! ((.not. water).or.(.not. iq.eq.igcm_h2o_vap)) then enddo ! of do iq=1,nq c --------- h2o_vap -------------------------------- c Treatement of the water frost c We use a subtimestep to take into account the release of latent heat do iq=1,nq if ((water).and.(iq.eq.igcm_h2o_vap)) then DO islope = 1,nslope DO ig=1,ngrid zqsurf(ig)=pqsurf(ig,igcm_h2o_ice,islope)/ & cos(pi*def_slope_mean(islope)/180.) watercap_tmp(ig) = watercap(ig,islope)/ & cos(pi*def_slope_mean(islope)/180.) ENDDO ! ig=1,ngrid c Computation of the subtimestep call make_tsub(ngrid,pdtsrf(:,islope),zqsurf, & ptimestep,dtmax,watercaptag, & nsubtimestep) c Calculation for turbulent exchange (rho Cd,h U (qatm - qsat(Tsurf)) with the surface (for ice) c initialization of ztsrf, which is surface temperature in c the subtimestep. saved_h2o_vap(:)= zq(:,1,igcm_h2o_vap) DO ig=1,ngrid c nsubtimestep(ig)=1 !for debug subtimestep = ptimestep/nsubtimestep(ig) ztsrf(ig)=ptsrf(ig,islope) ! +pdtsrf(ig)*subtimestep zq_tmp_vap(ig,:,:) =zq(ig,:,:) c Beginning of the subtimestep DO tsub=1,nsubtimestep(ig) if(tsub.eq.nsubtimestep(ig)) writeoutput = .true. zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay) & /float(nsubtimestep(ig)) if(old_wsublimation_scheme) then zb(1:ngrid,1)=zcdv(1:ngrid,islope)*zb0(1:ngrid,1) & /float(nsubtimestep(ig)) else zb(1:ngrid,1)=zcdh(1:ngrid,islope)*zb0(1:ngrid,1) & /float(nsubtimestep(ig)) endif zb(1:ngrid,1)=dryness(1:ngrid)*zb(1:ngrid,1) z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) zc(ig,nlay)=za(ig,nlay)*zq_tmp_vap(ig,nlay,iq)*z1(ig) zd(ig,nlay)=zb(ig,nlay)*z1(ig) DO ilay=nlay-1,2,-1 z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) zc(ig,ilay)=(za(ig,ilay)*zq_tmp_vap(ig,ilay,iq)+ $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) zd(ig,ilay)=zb(ig,ilay)*z1(ig) ENDDO z1(ig)=1./(za(ig,1)+zb(ig,1)+ $ zb(ig,2)*(1.-zd(ig,2))) zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+ $ zb(ig,2)*zc(ig,2)) * z1(ig) call watersat(1,ztsrf(ig),pplev(ig,1),qsat(ig)) old_h2o_vap(ig)=zq_tmp_vap(ig,1,igcm_h2o_vap) zd(ig,1)=zb(ig,1)*z1(ig) zq1temp(ig)=zc(ig,1)+ zd(ig,1)*qsat(ig) if(old_wsublimation_scheme) then zdqsdif_surf(ig)=rho(ig)*dryness(ig)*zcdv(ig,islope) & *(zq1temp(ig)-qsat(ig)) else zdqsdif_surf(ig)=rho(ig)*dryness(ig)*zcdh(ig,islope) & *(zq1temp(ig)-qsat(ig)) endif zdqsdif_tot(ig) = zdqsdif_surf(ig) ! -------------------------------------------------------------------------------------------------------------------------------- ! We consider here the possible interactions between the subsurface and the atmosphere in the case of the surface is free of frost ! when computing the complete adsorption/desorption model ! -------------------------------------------------------------------------------------------------------------------------------- if(.not.watercaptag(ig)) then if (((-(zdqsdif_surf(ig))* & subtimestep).gt.zqsurf(ig)) & .and.(pqsurf(ig,igcm_co2,islope).eq.0.)) then exchange = .true. else exchange = .false. endif else exchange = .false. endif ! -------------------------------------------------------------------------------------------------------------------------------- ! If one consider adsorption, all the fluxes to/from the surface/subsurface/atmosphere are computed here ! -------------------------------------------------------------------------------------------------------------------------------- if (adsorption_soil) then call soilwater(1,nlay,nq,nsoil, nqsoil, & ztsrf(ig),ptsoil(ig,:,islope),subtimestep, & exchange,qsat(ig),zq_tmp_vap(ig,:,:), & za(ig,:),zb(ig,:),zc(ig,:),zd(ig,:), & zdqsdif_surf(ig), zqsurf(ig), & qsoil(ig,:,:,islope), pplev(ig,1), rho(ig), & writeoutput,zdqsdif_regolith(ig,islope), & zq1temp_regolith(ig), & pore_icefraction(ig,:,islope)) if(.not.watercaptag(ig)) then if (exchange) then zq1temp(ig) = zq1temp_regolith(ig) zdqsdif_tot(ig)= & -zqsurf(ig)/subtimestep else zdqsdif_tot(ig) = zdqsdif_surf(ig) + & zdqsdif_regolith(ig,islope) ! boundary condition = qsat, but pdqsdif is calculated to update qsurf (including loss of surface ice to the subsurface) endif ! of "if exchange = true" endif ! of "if not.watercaptag" endif ! adsorption ! -------------------------------------------------------------------------------------------------------------------------------------- ! Here we do the same, but without computing the complete adsorpption/desorption. Note that it work only if one does not use adsorption ! If no subsurface ice, then the models computes the surface flux/water vapor flux as usual ! -------------------------------------------------------------------------------------------------------------------------------------- ! ------------------------------------------------------------------------------------------------------------------------------------------ ! First, We consider here the possible interactions between the subsurface and the atmosphere in the case of the surface is free of frost ! ------------------------------------------------------------------------------------------------------------------------------------------ if(.not.watercaptag(ig).and.(.not.adsorption_soil)) then if ((-zdqsdif_tot(ig)*subtimestep) & .ge.(zqsurf(ig))) then zdqsdif_tot(ig)=-zqsurf(ig)/subtimestep if((h2o_ice_depth(ig,islope).gt.0) & .and. (zqsurf(ig) .lt. 1e-9) & .and.lag_layer) then ! Atm <-> subsurface exchange, we need to update the exchange coefficient zb by a factor 1/1+R; R = zice*Cd,h*/D and add the flux from the subsurface zqsurf(ig) = 0 if(old_wsublimation_scheme) then resist(ig,islope)=h2o_ice_depth(ig,islope) & *zcdv(ig,islope)/d_coef(ig,islope) else resist(ig,islope)=h2o_ice_depth(ig,islope) & *zcdh(ig,islope)/d_coef(ig,islope) endif zb(ig,1)=zb(ig,1)*1/(1+resist(ig,islope)) ! change zb to account subsurface ice ! Now we add the flux from the subsurface ice : rho Cd,h U*(1/1+R) * (q1-qsat_ssi(Tice)) call compute_Tice(nsoil, ptsoil(ig,:,islope), & ztsrf(ig), & h2o_ice_depth(ig,islope), & Tice(ig,islope)) ! compute ice temperature call watersat(1,Tice(ig,islope),pplev(ig,1), & qsat_ssi(ig,islope)) qsat_ssi(ig,islope)=ztsrf(ig)/Tice(ig,islope) & *qsat_ssi(ig,islope) ! And now we solve correctly the system z1(ig)=1./(za(ig,1)+zb(ig,1)+ & zb(ig,2)*(1.-zd(ig,2))) zc(ig,1)=(za(ig,1)*zq_tmp_vap(ig,1,iq)+ & zb(ig,2)*zc(ig,2) + & (-zdqsdif_tot(ig)) *subtimestep) & * z1(ig) zd(ig,1)=zb(ig,1)*z1(ig) zq1temp(ig)=zc(ig,1)+ zd(ig,1) & *qsat_ssi(ig,islope) ! Flux from the subsurface if(old_wsublimation_scheme) then zdqsdif_ssi_atm(ig,islope)=rho(ig) & *dryness(ig)*zcdv(ig,islope) & *1/(1+resist(ig,islope)) & *(zq1temp(ig)-qsat_ssi(ig,islope)) else zdqsdif_ssi_atm(ig,islope)=rho(ig)* & *dryness(ig) *zcdh(ig,islope) & *1/(1+resist(ig,islope)) & *(zq1temp(ig)-qsat_ssi(ig,islope)) endif else ! No atm <-> subsurface exchange, we do it the usual way zdqsdif_tot(ig)=-zqsurf(ig)/subtimestep z1(ig)=1./(za(ig,1)+ zb(ig,2)*(1.-zd(ig,2))) zc(ig,1)=(za(ig,1)* & zq_tmp_vap(ig,1,igcm_h2o_vap)+ & zb(ig,2)*zc(ig,2) + & (-zdqsdif_tot(ig)) *subtimestep) *z1(ig) zq1temp(ig)=zc(ig,1) endif ! Of subsurface <-> atmosphere exchange endif ! sublimating more than available frost & surface - frost exchange endif !if .not.watercaptag(ig) ! -------------------------------------------------------------------------------------------------------------------------------- ! We check possible frost subsurface ice interaction: since there is no subsurface water ice mass reservoir represented (yet), ! we do not include their effect on the mass of surface frost. ! -------------------------------------------------------------------------------------------------------------------------------- if((h2o_ice_depth(ig,islope).gt.0).and.lag_layer & .and.(.not.adsorption_soil)) then ! First case: still frost at the surface but no watercaptag if(((watercaptag(ig))).or. & (((-zdqsdif_tot(ig)*subtimestep) & .lt.(zqsurf(ig))) & .and. (zqsurf(ig).gt.tol_frost))) then ! Still frost at the surface: we consider the possibility to have subsurface <-> frost exchange ! The flux between frost and ssi is D/zice *(qsat(Tsurf)-qsat_ssi(Tice)) call compute_Tice(nsoil, ptsoil(ig,:,islope), & ztsrf(ig), & h2o_ice_depth(ig,islope), & Tice(ig,islope)) ! compute ice temperature call watersat(1,Tice(ig,islope),pplev(ig,1), & qsat_ssi(ig,islope)) qsat_ssi(ig,islope)=ztsrf(ig)/Tice(ig,islope) & *qsat_ssi(ig,islope) zdqsdif_ssi_frost(ig,islope)= d_coef(ig,islope) & /h2o_ice_depth(ig,islope) & *rho(ig)*dryness(ig) & *(qsat(ig)-qsat_ssi(ig,islope)) ! Line to comment for now since we don't have a mass of subsurface frost ! in our computation (otherwise, we would not conserve the H2O mass in ! the system) if ((-zdqsdif_tot(ig)+zdqsdif_ssi_frost(ig,islope)) & *subtimestep.ge.(zqsurf(ig))) then zdqsdif_ssi_frost(ig,islope) = & zqsurf(ig)/subtimestep & + zdqsdif_tot(ig) endif zdqsdif_tot(ig) = zdqsdif_tot(ig) - & zdqsdif_ssi_frost(ig,islope) endif ! watercaptag or frost at the surface endif ! interaction frost <-> subsurface ice c Starting upward calculations for water : zq_tmp_vap(ig,1,igcm_h2o_vap)=zq1temp(ig) c Take into account the H2O latent heat impact on the surface temperature if (latentheat_surfwater) then lh=(2834.3-0.28*(ztsrf(ig)-To)- & 0.004*(ztsrf(ig)-To)*(ztsrf(ig)-To))*1.e+3 zdtsrf(ig,islope)= zdqsdif_tot(ig)*lh & /pcapcal(ig,islope) endif ! (latentheat_surfwater) then DO ilay=2,nlay zq_tmp_vap(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay) & *zq_tmp_vap(ig,ilay-1,iq) ENDDO c Subtimestep water budget : ztsrf(ig) = ztsrf(ig)+(pdtsrf(ig,islope) & + zdtsrf(ig,islope))*subtimestep zqsurf(ig)= zqsurf(ig)+( & zdqsdif_tot(ig))*subtimestep if (zqsurf(ig)<0 .and.(.not.watercaptag(ig))) then zqsurf(ig)=0 endif zdqsdif_ssi_atm_tot(ig,islope) = & zdqsdif_ssi_atm_tot(ig,islope) & + zdqsdif_ssi_atm(ig,islope)*subtimestep zdqsdif_ssi_frost_tot(ig,islope) = & zdqsdif_ssi_frost_tot(ig,islope) & +zdqsdif_ssi_frost(ig,islope)*subtimestep c Monitoring instantaneous latent heat flux in W.m-2 : zsurf_h2o_lh(ig,islope) = zsurf_h2o_lh(ig,islope)+ & (zdtsrf(ig,islope)*pcapcal(ig,islope)) & *subtimestep c We ensure that surface temperature can't rise above the solid domain if there c is still ice on the surface (oldschool) if(zqsurf(ig) & +zdqsdif_tot(ig)*subtimestep & .gt.frost_albedo_threshold) then ! if there is still ice, T cannot exceed To zdtsrf(ig,islope) = min(zdtsrf(ig,islope), & (To-ztsrf(ig))/subtimestep) ! ice melt case endif c End of the subtimestep ENDDO ! tsub=1,nsubtimestep c Integration of subtimestep temp and water budget : c (btw could also compute the post timestep temp and ice c by simply adding the subtimestep trend instead of this) surf_h2o_lh(ig,islope)= zsurf_h2o_lh(ig,islope)/ptimestep pdtsrf(ig,islope)= (ztsrf(ig) - & ptsrf(ig,islope))/ptimestep pdqsdif(ig,igcm_h2o_ice,islope)= & (zqsurf(ig)- pqsurf(ig,igcm_h2o_ice,islope)/ & cos(pi*def_slope_mean(islope)/180.)) & /ptimestep zdqsdif_ssi_atm_tot(ig,islope)= & zdqsdif_ssi_atm_tot(ig,islope)/ptimestep zdqsdif_ssi_frost_tot(ig,islope)= & zdqsdif_ssi_frost_tot(ig,islope)/ptimestep zdqsdif_ssi_tot(ig,islope) = & zdqsdif_ssi_atm_tot(ig,islope) & + zdqsdif_ssi_frost_tot(ig,islope) c if subliming more than qsurf(ice) and on watercaptag, water c sublimates from watercap reservoir (dwatercap_dif is <0) if(watercaptag(ig)) then if ((-pdqsdif(ig,igcm_h2o_ice,islope)*ptimestep) & .gt.(pqsurf(ig,igcm_h2o_ice,islope) & /cos(pi*def_slope_mean(islope)/180.))) then dwatercap_dif(ig,islope)= & pdqsdif(ig,igcm_h2o_ice,islope)+ & (pqsurf(ig,igcm_h2o_ice,islope) / & cos(pi*def_slope_mean(islope)/180.))/ptimestep pdqsdif(ig,igcm_h2o_ice,islope)= & - (pqsurf(ig,igcm_h2o_ice,islope)/ & cos(pi*def_slope_mean(islope)/180.))/ptimestep endif! ((-pdqsdif(ig)*ptimestep) endif !(watercaptag(ig)) then zq_slope_vap(ig,:,:,islope) = zq_tmp_vap(ig,:,:) ENDDO ! of DO ig=1,ngrid ENDDO ! islope c Some grid box averages: interface with the atmosphere DO ig = 1,ngrid DO ilay = 1,nlay zq(ig,ilay,iq) = 0. DO islope = 1,nslope zq(ig,ilay,iq) = zq(ig,ilay,iq) + $ zq_slope_vap(ig,ilay,iq,islope) * $ subslope_dist(ig,islope) ENDDO ENDDO ENDDO ! Recompute values in kg/m^2 slopped DO ig = 1,ngrid DO islope = 1,nslope pdqsdif(ig,igcm_h2o_ice,islope) = & pdqsdif(ig,igcm_h2o_ice,islope) & * cos(pi*def_slope_mean(islope)/180.) dwatercap_dif(ig,islope) = & dwatercap_dif(ig,islope) & * cos(pi*def_slope_mean(islope)/180.) ENDDO ENDDO END IF ! of IF ((water).and.(iq.eq.igcm_h2o_vap)) c --------- end of h2o_vap ---------------------------- c --------- hdo_vap ----------------------------------- c hdo_ice has already been with along h2o_ice c amongst "normal" tracers (ie not h2o_vap) if (hdo.and.(iq.eq.igcm_hdo_vap)) then zb(1:ngrid,2:nlay)=zkh(1:ngrid,2:nlay)*zb0(1:ngrid,2:nlay) zb(1:ngrid,1)=0 DO ig=1,ngrid z1(ig)=1./(za(ig,nlay)+zb(ig,nlay)) zc(ig,nlay)=za(ig,nlay)*zq(ig,nlay,iq)*z1(ig) zd(ig,nlay)=zb(ig,nlay)*z1(ig) ENDDO DO ilay=nlay-1,2,-1 DO ig=1,ngrid z1(ig)=1./(za(ig,ilay)+zb(ig,ilay)+ $ zb(ig,ilay+1)*(1.-zd(ig,ilay+1))) zc(ig,ilay)=(za(ig,ilay)*zq(ig,ilay,iq)+ $ zb(ig,ilay+1)*zc(ig,ilay+1))*z1(ig) zd(ig,ilay)=zb(ig,ilay)*z1(ig) ENDDO ENDDO hdoflux_meshavg(:) = 0. DO islope = 1,nslope pdqsdif_tmphdo(:,:) = pdqsdif(:,:,islope) & /cos(pi*def_slope_mean(islope)/180.) call watersat(ngrid,pdtsrf(:,islope)*ptimestep + & ptsrf(:,islope),pplev(:,1),qsat_tmp) CALL hdo_surfex(ngrid,nlay,nq,ptimestep, & zt,pplay,zq,pqsurf(:,:,islope), & saved_h2o_vap,qsat_tmp, & pdqsdif_tmphdo, & dwatercap_dif(:,islope)/cos(pi*def_slope_mean(islope)/180.), & hdoflux(:,islope)) pdqsdif(:,:,islope) = pdqsdif_tmphdo(:,:) * & cos(pi*def_slope_mean(islope)/180.) DO ig = 1,ngrid hdoflux_meshavg(ig) = hdoflux_meshavg(ig) + & hdoflux(ig,islope)*subslope_dist(ig,islope) ENDDO !ig = 1,ngrid ENDDO !islope = 1,nslope DO ig=1,ngrid z1(ig)=1./(za(ig,1)+zb(ig,1)+ $ zb(ig,2)*(1.-zd(ig,2))) zc(ig,1)=(za(ig,1)*zq(ig,1,iq)+ $ zb(ig,2)*zc(ig,2) + $ (-hdoflux_meshavg(ig)) *ptimestep) *z1(ig) !tracer flux from surface ENDDO DO ig=1,ngrid zq(ig,1,iq)=zc(ig,1) DO ilay=2,nlay zq(ig,ilay,iq)=zc(ig,ilay)+zd(ig,ilay)*zq(ig,ilay-1,iq) ENDDO ENDDO endif ! (hdo.and.(iq.eq.igcm_hdo_vap)) c --------- end of hdo ---------------------------- enddo ! of do iq=1,nq c --------- end of tracers ---------------------------- call write_output("surf_h2o_lh", & "Ground ice latent heat flux", & "W.m-2",surf_h2o_lh(:,iflat)) call write_output("nsubtimestep", & "Number of subtimesteps in vdifc", & "int",nsubtimestep(:)) call write_output('zdqsdif_ssi_frost_tot', & 'Flux between frost and subsurface ice','kg.m-2.s-1', & zdqsdif_ssi_frost_tot(:,iflat)) call write_output('zdqsdif_ssi_atm_tot', & 'Flux between atmosphere and subsurface ice','kg.m-2.s-1', & zdqsdif_ssi_atm_tot(:,iflat)) call write_output('zdqsdif_ssi_tot', & 'Total flux echange with subsurface ice','kg.m-2.s-1', & zdqsdif_ssi_tot(:,iflat)) C Diagnostic output for HDO ! if (hdo) then ! CALL write_output('hdoflux', ! & 'hdoflux', ! & ' ',hdoflux_meshavg(:)) ! CALL write_output('h2oflux', ! & 'h2oflux', ! & ' ',h2oflux(:)) ! endif c----------------------------------------------------------------------- c 8. calcul final des tendances de la diffusion verticale c ---------------------------------------------------- DO ilev = 1, nlay DO ig=1,ngrid pdudif(ig,ilev)=( zu(ig,ilev)- $ (pu(ig,ilev)+pdufi(ig,ilev)*ptimestep) )/ptimestep pdvdif(ig,ilev)=( zv(ig,ilev)- $ (pv(ig,ilev)+pdvfi(ig,ilev)*ptimestep) )/ptimestep hh = ph(ig,ilev)+pdhfi(ig,ilev)*ptimestep $ + (latcond*dmice(ig,ilev)/cpp)/ppopsk(ig,ilev) pdhdif(ig,ilev)=( zhs(ig,ilev)- hh )/ptimestep ENDDO ENDDO pdqdif(1:ngrid,1:nlay,1:nq)=(zq(1:ngrid,1:nlay,1:nq)- & (pq(1:ngrid,1:nlay,1:nq) & +pdqfi(1:ngrid,1:nlay,1:nq) & *ptimestep))/ptimestep c ** diagnostique final c ------------------ IF(lecrit) THEN PRINT*,'In vdif' PRINT*,'Ts (t) and Ts (t+st)' WRITE(*,'(a10,3a15)') s 'theta(t)','theta(t+dt)','u(t)','u(t+dt)' PRINT*,ptsrf(ngrid/2+1,:),ztsrf2(ngrid/2+1) DO ilev=1,nlay WRITE(*,'(4f15.7)') s ph(ngrid/2+1,ilev),zhs(ngrid/2+1,ilev), s pu(ngrid/2+1,ilev),zu(ngrid/2+1,ilev) ENDDO ENDIF END SUBROUTINE vdifc c==================================== SUBROUTINE make_tsub(naersize,dtsurf,qsurf,ptimestep, $ dtmax,watercaptag,ntsub) c Pas de temps adaptatif en estimant le taux de sublimation c et en adaptant avec un critere "dtmax" du chauffage a accomoder c dtmax est regle empiriquement (pour l'instant) a 0.5 K use surfdat_h, only: frost_albedo_threshold integer,intent(in) :: naersize real,intent(in) :: dtsurf(naersize) real,intent(in) :: qsurf(naersize) logical,intent(in) :: watercaptag(naersize) real,intent(in) :: ptimestep real,intent(in) :: dtmax real :: ztsub(naersize) integer :: i integer,intent(out) :: ntsub(naersize) do i=1,naersize c If not on permanent ice or thin frost layer : c do not use a subtimestep (ntsub=1) if ((qsurf(i).lt.frost_albedo_threshold).and. & (.not.watercaptag(i))) then ntsub(i) = 1 else ztsub(i) = ptimestep * dtsurf(i) / dtmax ntsub(i) = ceiling(abs(ztsub(i))) endif ! (qsurf(i).eq.0) then c c write(78,*), dtsurf*ptimestep, dtsurf, ntsub enddo! 1=1,ngrid END SUBROUTINE make_tsub c==================================== SUBROUTINE compute_Tice(nsoil, ptsoil, ptsurf, ice_depth, Tice) c Compute subsurface ice temperature by interpolating the temperatures between the two adjacent cells. use comsoil_h, only: layer, mlayer implicit none integer,intent(in) :: nsoil ! Number of soil layers real,intent(in) :: ptsoil(nsoil) ! Soil temperature (K) real,intent(in) :: ptsurf ! Soil temperature (K) real,intent(in) :: ice_depth ! Ice depth (m) real,intent(out) :: Tice ! Ice temperatures (K) c Local integer :: ik ! Loop variables integer :: indexice ! Index of the ice c Code: indexice = -1 if(ice_depth.lt.mlayer(0)) then indexice = 0. else do ik = 0,nsoil-2 ! go through all the layers to find the ice locations if((mlayer(ik).le.ice_depth).and. & (mlayer(ik+1).gt.ice_depth)) then indexice = ik+1 exit endif enddo endif if(indexice.lt.0) then call abort_physic("vdifc - compute Tice", & "subsurface ice is below the last soil layer",1) else if(indexice .ge. 1) then ! Linear inteprolation between soil temperature Tice = (ptsoil(indexice)-ptsoil(indexice+1)) & /(mlayer(indexice-1)-mlayer(indexice)) & *(ice_depth-mlayer(indexice)) + ptsoil(indexice+1) else ! Linear inteprolation between the 1st soil temperature and the surface temperature Tice = (ptsoil(1) - ptsurf)/mlayer(0) & *(ice_depth-mlayer(0)) + ptsoil(1) endif ! index ice >=0 endif !indexice <0 END SUBROUTINE compute_Tice END MODULE vdifc_mod