!======================================================================= ! CALLTHERM_INTERFACE !======================================================================= ! Main interface to the Martian thermal plume model ! This interface handles sub-timesteps for this model ! A call to this interface must be inserted in the main 'physics' routine ! NB: for information: ! In the Mars LMD-GCM, the thermal plume model is called after ! the vertical turbulent mixing scheme (Mellor and Yamada) ! and the surface layer scheme (Richardson-based surface layer + subgrid gustiness) ! Other routines called before the thermals model are ! radiative transfer and (orographic) gravity wave drag. ! ----------------------------------------------------------------------- ! Author : A. Colaitis 2011-01-05 (with updates 2011-2013) ! after C. Rio and F. Hourdin ! Institution : Laboratoire de Meteorologie Dynamique (LMD) Paris, France ! ----------------------------------------------------------------------- ! Corresponding author : A. Spiga aymeric.spiga_AT_upmc.fr ! ----------------------------------------------------------------------- ! ASSOCIATED FILES ! --> thermcell_main_mars.F90 ! --> thermcell_dqup.F90 ! --> comtherm_h.F90 ! ----------------------------------------------------------------------- ! Reference paper: ! A. Colaïtis, A. Spiga, F. Hourdin, C. Rio, F. Forget, and E. Millour. ! A thermal plume model for the Martian convective boundary layer. ! Journal of Geophysical Research (Planets), 118:1468-1487, July 2013. ! http://dx.doi.org/10.1002/jgre.20104 ! http://arxiv.org/abs/1306.6215 ! ----------------------------------------------------------------------- ! Reference paper for terrestrial plume model: ! C. Rio and F. Hourdin. ! A thermal plume model for the convective boundary layer : Representation of cumulus clouds. ! Journal of the Atmospheric Sciences, 65:407-425, 2008. ! ----------------------------------------------------------------------- SUBROUTINE calltherm_interface (ngrid,nlayer,nq, & & tracer,igcm_co2, & & zzlev,zzlay, & & ptimestep,pu,pv,pt,pq,pdu,pdv,pdt,pdq,q2, & & pplay,pplev,pphi,zpopsk, & & pdu_th,pdv_th,pdt_th,pdq_th,lmax,zmaxth,pbl_dtke, & & pdhdif,hfmax,wstar,sensibFlux) use comtherm_h use tracer_mod, only: nqmx,noms implicit none ! SHARED VARIABLES. This needs adaptations in another climate model. #include "comcstfi.h" !contains physical constant values such as ! "g" : gravitational acceleration (m.s-2) ! "r" : recuced gas constant (J.K-1.mol-1) ! "cpp" : specific heat of the atmosphere (J.kg-1.K-1) !-------------------------------------------------------- ! Input Variables !-------------------------------------------------------- INTEGER, INTENT(IN) :: ngrid ! number of horizontal grid points INTEGER, INTENT(IN) :: nlayer ! number of vertical grid points INTEGER, INTENT(IN) :: nq ! number of tracer species REAL, INTENT(IN) :: ptimestep !timestep (s) REAL, INTENT(IN) :: pplev(ngrid,nlayer+1) !intermediate pressure levels (Pa) REAL, INTENT(IN) :: pplay(ngrid,nlayer) !Pressure at the middle of the layers (Pa) REAL, INTENT(IN) :: pphi(ngrid,nlayer) !Geopotential at the middle of the layers (m2s-2) REAL, INTENT(IN) :: pu(ngrid,nlayer),pv(ngrid,nlayer) !u,v components of the wind (ms-1) REAL, INTENT(IN) :: pt(ngrid,nlayer),pq(ngrid,nlayer,nq)!temperature (K) and tracer concentration (kg/kg) REAL, INTENT(IN) :: zzlay(ngrid,nlayer) ! altitude at the middle of the layers REAL, INTENT(IN) :: zzlev(ngrid,nlayer+1) ! altitude at layer boundaries LOGICAL, INTENT(IN) :: tracer ! =.true. if tracers are present and to be transported INTEGER, INTENT(IN) :: igcm_co2 ! index of the CO2 tracer in mixing ratio array ! --> 0 if no tracer is CO2 (or no tracer at all) ! --> this prepares special treatment for polar night mixing ! (see thermcell_main_mars) REAL, INTENT(IN) :: pdu(ngrid,nlayer),pdv(ngrid,nlayer) ! wind velocity change from routines called ! before thermals du/dt (m/s/s) REAL, INTENT(IN) :: pdq(ngrid,nlayer,nq) ! tracer concentration change from routines called ! before thermals dq/dt (kg/kg/s) REAL, INTENT(IN) :: pdt(ngrid,nlayer) ! temperature change from routines called before thermals dT/dt (K/s) REAL, INTENT(IN) :: q2(ngrid,nlayer+1) ! turbulent kinetic energy REAL, INTENT(IN) :: zpopsk(ngrid,nlayer) ! ratio of pressure at middle of layer to surface pressure, ! to the power r/cp, i.e. zpopsk=(pplay(ig,l)/pplev(ig,1))**rcp REAL, INTENT(IN) :: pdhdif(ngrid,nlayer) ! potential temperature change from turbulent diffusion scheme dT/dt (K/s) REAL, INTENT(IN) :: sensibFlux(ngrid) ! sensible heat flux computed from surface layer scheme !-------------------------------------------------------- ! Output Variables !-------------------------------------------------------- REAL, INTENT(OUT) :: pdu_th(ngrid,nlayer) ! wind velocity change from thermals du/dt (m/s/s) REAL, INTENT(OUT) :: pdv_th(ngrid,nlayer) ! wind velocity change from thermals dv/dt (m/s/s) REAL, INTENT(OUT) :: pdt_th(ngrid,nlayer) ! temperature change from thermals dT/dt (K/s) REAL, INTENT(OUT) :: pdq_th(ngrid,nlayer,nq) ! tracer change from thermals dq/dt (kg/kg/s) INTEGER, INTENT(OUT) :: lmax(ngrid) ! layer number reacher by thermals in grid point REAL, INTENT(OUT) :: zmaxth(ngrid) ! equivalent to lmax, but in (m), interpolated REAL, INTENT(OUT) :: pbl_dtke(ngrid,nlayer+1) ! turbulent kinetic energy change from thermals dtke/dt REAL, INTENT(OUT) :: wstar(ngrid) ! free convection velocity (m/s) !-------------------------------------------------------- ! Thermals local variables !-------------------------------------------------------- REAL zu(ngrid,nlayer), zv(ngrid,nlayer) REAL zt(ngrid,nlayer) REAL d_t_ajs(ngrid,nlayer) REAL d_u_ajs(ngrid,nlayer), d_q_ajs(ngrid,nlayer,nq) REAL d_v_ajs(ngrid,nlayer) REAL fm_therm(ngrid,nlayer+1), entr_therm(ngrid,nlayer) REAL detr_therm(ngrid,nlayer),detrmod(ngrid,nlayer) REAL zw2(ngrid,nlayer+1) REAL fraca(ngrid,nlayer+1),zfraca(ngrid,nlayer+1) REAL q_therm(ngrid,nlayer), pq_therm(ngrid,nlayer,nq) REAL q2_therm(ngrid,nlayer), dq2_therm(ngrid,nlayer) REAL lmax_real(ngrid) REAL masse(ngrid,nlayer) INTEGER l,ig,iq,ii(1),k CHARACTER (LEN=20) modname !-------------------------------------------------------- ! Local variables for sub-timestep !-------------------------------------------------------- REAL d_t_the(ngrid,nlayer), d_q_the(ngrid,nlayer,nq) INTEGER isplit REAL fact REAL zfm_therm(ngrid,nlayer+1),zdt REAL zentr_therm(ngrid,nlayer),zdetr_therm(ngrid,nlayer) REAL zheatFlux(ngrid,nlayer) REAL zheatFlux_down(ngrid,nlayer) REAL zbuoyancyOut(ngrid,nlayer) REAL zbuoyancyEst(ngrid,nlayer) REAL zzw2(ngrid,nlayer+1) REAL zmax(ngrid) INTEGER ndt,zlmax !-------------------------------------------------------- ! Diagnostics !-------------------------------------------------------- REAL heatFlux(ngrid,nlayer) REAL heatFlux_down(ngrid,nlayer) REAL buoyancyOut(ngrid,nlayer) REAL buoyancyEst(ngrid,nlayer) REAL hfmax(ngrid),wmax(ngrid) REAL pbl_teta(ngrid),dteta(ngrid,nlayer) REAL rpdhd(ngrid,nlayer) REAL wtdif(ngrid,nlayer),rho(ngrid,nlayer) REAL wtth(ngrid,nlayer) ! ********************************************************************** ! Initializations ! ********************************************************************** lmax(:)=0 pdu_th(:,:)=0. pdv_th(:,:)=0. pdt_th(:,:)=0. entr_therm(:,:)=0. detr_therm(:,:)=0. q2_therm(:,:)=0. dq2_therm(:,:)=0. pbl_dtke(:,:)=0. fm_therm(:,:)=0. zw2(:,:)=0. fraca(:,:)=0. zfraca(:,:)=0. if (tracer) then pdq_th(:,:,:)=0. end if d_t_ajs(:,:)=0. d_u_ajs(:,:)=0. d_v_ajs(:,:)=0. d_q_ajs(:,:,:)=0. heatFlux(:,:)=0. heatFlux_down(:,:)=0. buoyancyOut(:,:)=0. buoyancyEst(:,:)=0. zmaxth(:)=0. lmax_real(:)=0. ! ********************************************************************** ! Preparing inputs for the thermals: increment tendancies ! from other subroutines called before the thermals model ! ********************************************************************** zu(:,:)=pu(:,:)+pdu(:,:)*ptimestep ! u-component of wind velocity zv(:,:)=pv(:,:)+pdv(:,:)*ptimestep ! v-component of wind velocity zt(:,:)=pt(:,:)+pdt(:,:)*ptimestep ! temperature pq_therm(:,:,:)=0. ! tracer concentration if(qtransport_thermals) then if(tracer) then ! tracer is a logical that is true if tracer must be transported in the GCM physics pq_therm(:,:,:)=pq(:,:,:)+pdq(:,:,:)*ptimestep ! tracer concentration endif endif IF(dtke_thermals) THEN DO l=1,nlayer q2_therm(:,l)=0.5*(q2(:,l)+q2(:,l+1)) ENDDO ENDIF ! ********************************************************************** ! --> CALLTHERM ! SUB-TIMESTEP LOOP ! ********************************************************************** zdt=ptimestep/REAL(nsplit_thermals) !subtimestep DO isplit=1,nsplit_thermals !temporal loop on the subtimestep ! Initialization of intermediary variables zzw2(:,:)=0. zmax(:)=0. lmax(:)=0 if (nq .ne. 0 .and. igcm_co2 .ne. 0) then !initialize co2 tracer tendancy d_q_the(:,:,igcm_co2)=0. endif ! CALL to main thermal routine CALL thermcell_main_mars(ngrid,nlayer,nq & & ,tracer,igcm_co2 & & ,zdt & & ,pplay,pplev,pphi,zzlev,zzlay & & ,zu,zv,zt,pq_therm,q2_therm & & ,d_t_the,d_q_the & & ,zfm_therm,zentr_therm,zdetr_therm,lmax,zmax & & ,zzw2,fraca,zpopsk & & ,zheatFlux,zheatFlux_down & & ,zbuoyancyOut,zbuoyancyEst) fact=1./REAL(nsplit_thermals) ! Update thermals tendancies d_t_the(:,:)=d_t_the(:,:)*ptimestep*fact !temperature if (igcm_co2 .ne. 0) then d_q_the(:,:,igcm_co2)=d_q_the(:,:,igcm_co2)*ptimestep*fact !co2 mass mixing ratio endif zmaxth(:)=zmaxth(:)+zmax(:)*fact !thermals height lmax_real(:)=lmax_real(:)+float(lmax(:))*fact !thermals height index fm_therm(:,:)=fm_therm(:,:) & !upward mass flux & +zfm_therm(:,:)*fact entr_therm(:,:)=entr_therm(:,:) & !entrainment mass flux & +zentr_therm(:,:)*fact detr_therm(:,:)=detr_therm(:,:) & !detrainment mass flux & +zdetr_therm(:,:)*fact zfraca(:,:)=zfraca(:,:) + fraca(:,:)*fact !updraft fractional coverage heatFlux(:,:)=heatFlux(:,:) & !upward heat flux & +zheatFlux(:,:)*fact heatFlux_down(:,:)=heatFlux_down(:,:) & !downward heat flux & +zheatFlux_down(:,:)*fact buoyancyOut(:,:)=buoyancyOut(:,:) & !plume final buoyancy & +zbuoyancyOut(:,:)*fact buoyancyEst(:,:)=buoyancyEst(:,:) & !plume estimated buoyancy used for vertical velocity computation & +zbuoyancyEst(:,:)*fact zw2(:,:)=zw2(:,:) + zzw2(:,:)*fact !vertical velocity ! Save tendancies d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_the(:,:) !temperature tendancy (delta T) if (igcm_co2 .ne. 0) then d_q_ajs(:,:,igcm_co2)=d_q_ajs(:,:,igcm_co2)+d_q_the(:,:,igcm_co2) !tracer tendancy (delta q) endif ! Increment temperature and co2 concentration for next pass in subtimestep loop zt(:,:) = zt(:,:) + d_t_the(:,:) !temperature if (igcm_co2 .ne. 0) then pq_therm(:,:,igcm_co2) = & & pq_therm(:,:,igcm_co2) + d_q_the(:,:,igcm_co2) !co2 tracer endif ENDDO ! isplit !**************************************************************** ! Do we have thermals that are too high ? lmax(:)=nint(lmax_real(:)) zlmax=MAXVAL(lmax(:))+2 if (zlmax .ge. nlayer) then print*,'thermals have reached last layer of the model' print*,'this is not good !' endif ! Now that we have computed total entrainment and detrainment, we can ! advect u, v, and q in thermals. (potential temperature and co2 MMR ! have already been advected in thermcell_main because they are coupled ! to the determination of the thermals caracteristics). This is done ! separatly because u,v, and q are not used in thermcell_main for ! any thermals-related computation : they are purely passive. ! mass of cells do l=1,nlayer masse(:,l)=(pplev(:,l)-pplev(:,l+1))/g enddo ! recompute detrainment mass flux from entrainment and upward mass flux ! this ensure mass flux conservation detrmod(:,:)=0. do l=1,zlmax do ig=1,ngrid detrmod(ig,l)=fm_therm(ig,l)-fm_therm(ig,l+1) & & +entr_therm(ig,l) if (detrmod(ig,l).lt.0.) then entr_therm(ig,l)=entr_therm(ig,l)-detrmod(ig,l) detrmod(ig,l)=0. endif enddo enddo ! u component of wind velocity advection in thermals ! result is a derivative (d_u_ajs in m/s/s) ndt=10 call thermcell_dqup(ngrid,nlayer,ptimestep & & ,fm_therm,entr_therm,detrmod, & & masse,zu,d_u_ajs,ndt,zlmax) ! v component of wind velocity advection in thermals ! result is a derivative (d_v_ajs in m/s/s) call thermcell_dqup(ngrid,nlayer,ptimestep & & ,fm_therm,entr_therm,detrmod, & & masse,zv,d_v_ajs,ndt,zlmax) ! non co2 tracers advection in thermals ! result is a derivative (d_q_ajs in kg/kg/s) if (nq .ne. 0.) then DO iq=1,nq if (iq .ne. igcm_co2) then call thermcell_dqup(ngrid,nlayer,ptimestep & & ,fm_therm,entr_therm,detrmod, & & masse,pq_therm(:,:,iq),d_q_ajs(:,:,iq),ndt,zlmax) endif ENDDO endif ! tke advection in thermals ! result is a tendancy (d_u_ajs in J) if (dtke_thermals) then call thermcell_dqup(ngrid,nlayer,ptimestep & & ,fm_therm,entr_therm,detrmod, & & masse,q2_therm,dq2_therm,ndt,zlmax) endif ! compute wmax for diagnostics DO ig=1,ngrid wmax(ig)=MAXVAL(zw2(ig,:)) ENDDO ! ********************************************************************** ! ********************************************************************** ! ********************************************************************** ! CALLTHERM END ! ********************************************************************** ! ********************************************************************** ! ********************************************************************** ! ********************************************************************** ! Preparing outputs ! ********************************************************************** do l=1,zlmax pdu_th(:,l)=d_u_ajs(:,l) pdv_th(:,l)=d_v_ajs(:,l) enddo ! if tracers are transported in thermals, update output variables, else these are 0. if(qtransport_thermals) then if(tracer) then do iq=1,nq if (iq .ne. igcm_co2) then do l=1,zlmax pdq_th(:,l,iq)=d_q_ajs(:,l,iq) !non-co2 tracers d_q_ajs are dq/dt (kg/kg/s) enddo else do l=1,zlmax pdq_th(:,l,iq)=d_q_ajs(:,l,iq)/ptimestep !co2 tracer d_q_ajs is dq (kg/kg) enddo endif enddo endif endif ! if tke is transported in thermals, update output variable, else this is 0. IF(dtke_thermals) THEN DO l=2,nlayer pbl_dtke(:,l)=0.5*(dq2_therm(:,l-1)+dq2_therm(:,l)) ENDDO pbl_dtke(:,1)=0.5*dq2_therm(:,1) pbl_dtke(:,nlayer+1)=0. ENDIF ! update output variable for temperature. d_t_ajs is delta T in (K), pdt_th is dT/dt in (K/s) do l=1,zlmax pdt_th(:,l)=d_t_ajs(:,l)/ptimestep enddo ! ********************************************************************** ! SURFACE LAYER INTERFACE ! Compute the free convection velocity w* scale for surface layer gustiness ! speed parameterization. The computed value of w* will be used at the next ! timestep to modify surface-atmosphere exchange fluxes, because the surface ! layer scheme and diffusion are called BEFORE the thermals. (outside of these ! routines) ! ********************************************************************** ! Potential temperature gradient dteta(:,nlayer)=0. DO l=1,nlayer-1 DO ig=1, ngrid dteta(ig,l) = ((zt(ig,l+1)-zt(ig,l))/zpopsk(ig,l)) & & /(zzlay(ig,l+1)-zzlay(ig,l)) ENDDO ENDDO ! Computation of the PBL mixed layer temperature DO ig=1, ngrid ii=MINLOC(abs(dteta(ig,1:lmax(ig)))) pbl_teta(ig) = zt(ig,ii(1))/zpopsk(ig,ii(1)) ENDDO ! In order to have an accurate w*, we must add the heat flux from the ! diffusion scheme to the computation of the maximum heat flux hfmax ! Here pdhdif is the potential temperature change from the diffusion ! scheme (Mellor and Yamada, see paper section 6, paragraph 57) ! compute rho as it is after the diffusion rho(:,:)=pplay(:,:) & & /(r*(pt(:,:)+pdhdif(:,:)*zpopsk(:,:)*ptimestep)) ! integrate -rho*pdhdif rpdhd(:,:)=0. DO ig=1,ngrid DO l=1,lmax(ig) rpdhd(ig,l)=0. DO k=1,l rpdhd(ig,l)=rpdhd(ig,l)-rho(ig,k)*pdhdif(ig,k)* & & (zzlev(ig,k+1)-zzlev(ig,k)) ENDDO rpdhd(ig,l)=rpdhd(ig,l)-sensibFlux(ig)/cpp ENDDO ENDDO ! compute w'theta' (vertical turbulent flux of temperature) from ! the diffusion scheme wtdif(:,:)=rpdhd(:,:)/rho(:,:) ! Now we compute the contribution of the thermals to w'theta': ! compute rho as it is after the thermals rho(:,:)=pplay(:,:)/(r*(zt(:,:))) ! integrate -rho*pdhdif DO ig=1,ngrid DO l=1,lmax(ig) rpdhd(ig,l)=0. DO k=1,l rpdhd(ig,l)=rpdhd(ig,l)-rho(ig,k)*(pdt_th(ig,k)/zpopsk(ig,k))* & & (zzlev(ig,k+1)-zzlev(ig,k)) ENDDO rpdhd(ig,l)=rpdhd(ig,l)+ & & rho(ig,1)*(heatFlux(ig,1)+heatFlux_down(ig,1)) ENDDO ENDDO rpdhd(:,nlayer)=0. ! compute w'teta' from thermals wtth(:,:)=rpdhd(:,:)/rho(:,:) ! Add vertical turbulent heat fluxes from the thermals and the diffusion scheme ! and compute the maximum DO ig=1,ngrid hfmax(ig)=MAXVAL(wtth(ig,:)+wtdif(ig,:)) ENDDO ! Finally we can compute the free convection velocity scale ! We follow Spiga et. al 2010 (QJRMS) ! ------------ DO ig=1, ngrid IF (zmax(ig) .gt. 0.) THEN wstar(ig)=(g*zmaxth(ig)*hfmax(ig)/pbl_teta(ig))**(1./3.) ELSE wstar(ig)=0. ENDIF ENDDO END