MODULE calcratqs_multi_mod !============================================= ! module containing subroutines that take ! into account the effect of convection, orography, ! surface heterogeneities and subgrid-scale ! turbulence on ratqs, i.e. on the width of the ! total water subgrid distribution. !============================================= IMPLICIT NONE ! Include !============================================= INCLUDE "YOETHF.h" INCLUDE "YOMCST.h" CONTAINS !======================================================================== SUBROUTINE calcratqs_inter(klon,klev,iflag_ratqs,pdtphys, & ratqsbas, wake_deltaq, wake_s, q_seri,qtc_cv, sigt_cv, & ratqs_inter) USE ioipsl_getin_p_mod, ONLY : getin_p implicit none !======================================================================== ! L. d'Alençon, 25/02/2021 ! Cette subroutine calcule une valeur de ratqsbas interactive dépendant de la présence de poches froides dans l'environnement. ! Elle est appelée par la subroutine calcratqs lorsque iflag_ratqs = 10. !======================================================================== ! Declarations LOGICAL, SAVE :: first = .TRUE. !$OMP THREADPRIVATE(first) ! Input integer,intent(in) :: klon,klev,iflag_ratqs real,intent(in) :: pdtphys,ratqsbas real, dimension(klon,klev),intent(in) :: wake_deltaq, q_seri,qtc_cv, sigt_cv real, dimension(klon),intent(in) :: wake_s ! Output real, dimension(klon,klev),intent(inout) :: ratqs_inter ! local integer i,k real, dimension(klon,klev) :: wake_dq REAL, SAVE :: a_ratqs_cv !$OMP THREADPRIVATE(a_ratqs_cv) REAL, SAVE :: tau_ratqs_wake !$OMP THREADPRIVATE(tau_ratqs_wake) REAL, SAVE :: a_ratqs_wake !$OMP THREADPRIVATE(a_ratqs_wake) real, dimension(klon) :: max_wake_dq, max_dqconv,max_sigt !------------------------------------------------------------------------- ! Caclul de ratqs_inter !------------------------------------------------------------------------- ! if (first) then tau_ratqs_wake = 3600. ! temps de relaxation de la variabilité a_ratqs_wake = 3. ! paramètre pilotant l'importance du terme dépendant des poches froides a_ratqs_cv = 0.5 CALL getin_p('tau_ratqs_wake', tau_ratqs_wake) CALL getin_p('a_ratqs_wake', a_ratqs_wake) CALL getin_p('a_ratqs_cv', a_ratqs_cv) first=.false. endif max_wake_dq(:) = 0. max_dqconv (:) = 0 max_sigt(:) = 0. if (iflag_ratqs.eq.10) then do k=1,klev do i=1,klon max_wake_dq(i) = max(abs(wake_deltaq(i,k)),max_wake_dq(i)) max_sigt(i) = max(abs(sigt_cv(i,k)),max_sigt(i)) max_dqconv(i) = max(abs(q_seri(i,k) - qtc_cv(i,k)),max_dqconv(i)) enddo enddo do k=1,klev do i=1,klon ratqs_inter(i,k)= ratqs_inter(i,k)*exp(-pdtphys/tau_ratqs_wake) + & a_ratqs_wake*(max_wake_dq(i)*(wake_s(i)**0.5/(1.-wake_s(i))))*(1.-exp(-pdtphys/tau_ratqs_wake))/q_seri(i,1) if (ratqs_inter(i,k)0) then ! ratqs_inter(i,k) = abs(q_seri(i,k) - qtc_cv(i,k)) ! endif enddo enddo endif return end !------------------------------------------------------------------ SUBROUTINE calcratqs_oro(klon,klev,qsat,temp,pplay,paprs,ratqs_oro) ! Etienne Vignon, November 2021: effect of subgrid orography on ratqs USE phys_state_var_mod, ONLY: zstd USE phys_state_var_mod, ONLY: pctsrf USE indice_sol_mod, only: nbsrf, is_lic, is_ter IMPLICIT NONE ! Declarations !-------------- ! INPUTS INTEGER, INTENT(IN) :: klon ! number of horizontal grid points INTEGER, INTENT(IN) :: klev ! number of vertical layers REAL, DIMENSION(klon,klev), INTENT(IN) :: qsat ! saturation specific humidity [kg/kg] REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! air temperature [K] REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! air pressure, layer's center [Pa] REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! air pressure, lower inteface [Pa] ! OUTPUTS REAL, DIMENSION(klon,klev), INTENT(out) :: ratqs_oro ! ratqs profile due to subgrid orography ! LOCAL INTEGER :: i,k REAL, DIMENSION(klon) :: orogradT,xsi0 REAL, DIMENSION (klon,klev) :: zlay REAL :: Lvs, temp0 ! Calculation of the near-surface temperature gradient along the topography !-------------------------------------------------------------------------- ! at the moment, we fix it at a constant value (moist adiab. lapse rate) orogradT(:)=-6.5/1000. ! K/m ! Calculation of near-surface surface ratqs !------------------------------------------- DO i=1,klon temp0=temp(i,1) IF (temp0 .LT. RTT) THEN Lvs=RLSTT ELSE Lvs=RLVTT ENDIF xsi0(i)=zstd(i)*ABS(orogradT(i))*Lvs/temp0/temp0/RV ratqs_oro(i,1)=xsi0(i) END DO ! Vertical profile of ratqs assuming an exponential decrease with height !------------------------------------------------------------------------ ! calculation of geop. height AGL zlay(:,1)= RD*temp(:,1)/(0.5*(paprs(:,1)+pplay(:,1))) & *(paprs(:,1)-pplay(:,1))/RG DO k=2,klev DO i = 1, klon zlay(i,k)= zlay(i,k-1)+RD*0.5*(temp(i,k-1)+temp(i,k)) & /paprs(i,k)*(pplay(i,k-1)-pplay(i,k))/RG ratqs_oro(i,k)=MAX(0.0,pctsrf(i,is_ter)*xsi0(i)*exp(-zlay(i,k)/MAX(zstd(i),1.))) END DO END DO END SUBROUTINE calcratqs_oro !============================================= SUBROUTINE calcratqs_hetero(klon,klev,t2m,q2m,temp,q,pplay,paprs,ratqs_hetero) ! Etienne Vignon, November 2021 ! Effect of subgrid surface heterogeneities on ratqs USE phys_local_var_mod, ONLY: s_pblh USE phys_state_var_mod, ONLY: pctsrf USE indice_sol_mod USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF IMPLICIT NONE include "YOMCST.h" ! INPUTS INTEGER, INTENT(IN) :: klon ! number of horizontal grid points INTEGER, INTENT(IN) :: klev ! number of vertical layers REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: t2m ! 2m temperature for each tile [K] REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: q2m ! 2m specific humidity for each tile [kg/kg] REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! air temperature [K] REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! specific humidity [kg/kg] REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! air pressure, layer's center [Pa] REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! air pressure, lower inteface [Pa] ! OUTPUTS REAL, DIMENSION(klon,klev), INTENT(out) :: ratqs_hetero ! ratsq profile due to surface heterogeneities INTEGER :: i,k,nsrf REAL, DIMENSION(klon) :: xsi0, ratiom, qsat2m, dqsatdT REAL, DIMENSION (klon,klev) :: zlay ! Calculation of near-surface surface ratqs !------------------------------------------- ratiom(:)=0. xsi0(:)=0. DO nsrf=1,nbsrf CALL CALC_QSAT_ECMWF(klon,t2m(:,nsrf),q2m(:,nsrf),paprs(:,1),RTT,0,.false.,qsat2m,dqsatdT) ratiom(:)=ratiom(:)+pctsrf(:,nsrf)*(q2m(:,nsrf)/qsat2m(:)) xsi0(:)=xsi0(:)+pctsrf(:,nsrf)*((q2m(:,nsrf)/qsat2m(:)-ratiom(:))**2) END DO xsi0(:)=sqrt(xsi0(:))/(ratiom(:)+1E-6) ! Vertical profile of ratqs assuming an exponential decrease with height !------------------------------------------------------------------------ ! calculation of geop. height AGL zlay(:,1)= RD*temp(:,1)/(0.5*(paprs(:,1)+pplay(:,1))) & *(paprs(:,1)-pplay(:,1))/RG ratqs_hetero(:,1)=xsi0(:) DO k=2,klev DO i = 1, klon zlay(i,k)= zlay(i,k-1)+RD*0.5*(temp(i,k-1)+temp(i,k)) & /paprs(i,k)*(pplay(i,k-1)-pplay(i,k))/RG ratqs_hetero(i,k)=MAX(xsi0(i)*exp(-zlay(i,k)/(s_pblh(i)+1.0)),0.0) END DO END DO END SUBROUTINE calcratqs_hetero !============================================= SUBROUTINE calcratqs_tke(klon,klev,pdtphys,temp,q,qsat,pplay,paprs,tke,tke_dissip,lmix,wprime,ratqs_tke) ! References: ! ! Etienne Vignon: effect of subgrid turbulence on ratqs ! ! Field, P.R., Hill, A., Furtado, K., Korolev, A., 2014b. Mixed-phase clouds in a turbulent environment. Part ! 2: analytic treatment. Q. J. R. Meteorol. Soc. 21, 2651–2663. https://doi.org/10.1002/qj.2175. ! ! Furtado, K., Field, P.R., Boutle, I.A., Morcrette, C.R., Wilkinson, J., 2016. A physically-based, subgrid ! parametrization for the production and maintenance of mixed-phase clouds in a general circulation ! model. J. Atmos. Sci. 73, 279–291. https://doi.org/10.1175/JAS-D-15-0021. USE phys_local_var_mod, ONLY: omega IMPLICIT NONE ! INPUTS INTEGER, INTENT(IN) :: klon ! number of horizontal grid points INTEGER, INTENT(IN) :: klev ! number of vertical layers REAL, INTENT(IN) :: pdtphys ! physics time step [s] REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! air temperature [K] REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! specific humidity [kg/kg] REAL, DIMENSION(klon,klev), INTENT(IN) :: qsat ! saturation specific humidity [kg/kg] REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! air pressure, layer's center [Pa] REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! air pressure, lower inteface [Pa] REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke ! Turbulent Kinetic Energy [m2/s2] REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke_dissip ! Turbulent Kinetic Energy Dissipation rate [m2/s3] REAL, DIMENSION(klon,klev+1), INTENT(IN) :: lmix ! Turbulent mixing length REAL, DIMENSION(klon,klev+1), INTENT(IN) :: wprime ! Turbulent vertical velocity scale [m/s] ! OUTPUTS REAL, DIMENSION(klon,klev), INTENT(out) :: ratqs_tke ! ratsq profile due to subgrid TKE ! LOCAL INTEGER :: i, k REAL :: AA, DD, NW, AAprime, VARLOG,rho,Lvs,taue,lhomo,dissmin,maxvarlog REAL, DIMENSION(klon,klev) :: sigmaw,w REAL, PARAMETER :: C0=10.0 REAL, PARAMETER :: lmin=0.001 REAL, PARAMETER :: ratqsmin=1E-6 REAL, PARAMETER :: ratqsmax=0.5 ! Calculation of large scale and turbulent vertical velocities !--------------------------------------------------------------- DO k=1,klev DO i=1,klon rho=pplay(i,k)/temp(i,k)/RD w(i,k)=-rho*RG*omega(i,k) sigmaw(i,k)=0.5*(wprime(i,k+1)+wprime(i,k)) ! turbulent vertical velocity at the middle of model layers. END DO END DO ! Calculation of ratqs !--------------------------------------------------------------- ratqs_tke(:,1)=ratqsmin ! set to a very low value to avoid division by 0 in order parts ! of the code DO k=2,klev ! we start from second model level since TKE is not defined at k=1 DO i=1,klon IF (temp(i,k) .LT. RTT) THEN Lvs=RLSTT ELSE Lvs=RLVTT ENDIF dissmin=0.01*(0.5*(tke(i,k)+tke(i,k+1))/pdtphys) maxvarlog=LOG(1.0+ratqsmax**2)! to prevent ratqs from exceeding an arbitrary threshold value AA=RG*(Lvs/(RCPD*temp(i,k)*temp(i,k)*RV) - 1./(RD*temp(i,k))) lhomo=MAX(0.5*(lmix(i,k)+lmix(i,k+1)),lmin) taue=(lhomo*lhomo/MAX(0.5*(tke_dissip(i,k)+tke_dissip(i,k+1)),dissmin))**(1./3) ! Fields et al. 2014 DD=1.0/taue NW=(sigmaw(i,k)**2)*SQRT(2./(C0*MAX(0.5*(tke_dissip(i,k)+tke_dissip(i,k+1)),dissmin))) AAprime=AA*NW VARLOG=AAprime/2./DD VARLOG=MIN(VARLOG,maxvarlog) ratqs_tke(i,k)=SQRT(MAX(EXP(VARLOG)-1.0,ratqsmin)) END DO END DO END SUBROUTINE calcratqs_tke END MODULE calcratqs_multi_mod