- Timestamp:
- Feb 29, 2024, 7:42:12 PM (3 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.F90
r4828 r4835 2 2 3 3 contains 4 subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,qv,qb s,pplay,paprs,dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)4 subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,qv,qb,pplay,paprs,dtemp_bs,dqv_bs,dqb_bs,bsfl,precip_bs) 5 5 6 6 !============================================================================== … … 11 11 12 12 13 use lmdz_blowing_snow_ini, only : coef_eva_bs,RTT,RD,RG,expo_eva_bs, fallv_bs, qbsmin14 use lmdz_blowing_snow_ini, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, tbsmelt, taumeltbs013 use lmdz_blowing_snow_ini, only : iflag_sublim_bs, iflag_sedim_bs, coef_sub_bs,RTT,RD,RG,fallv_bs 14 use lmdz_blowing_snow_ini, only : qbmin, RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, RV, RPI, tbsmelt, taumeltbs0, rhobs, r_bs 15 15 USE lmdz_lscp_tools, only : calc_qsat_ecmwf 16 16 … … 29 29 real, intent(in), dimension(ngrid,nlay) :: temp 30 30 real, intent(in), dimension(ngrid,nlay) :: qv 31 real, intent(in), dimension(ngrid,nlay) :: qb s31 real, intent(in), dimension(ngrid,nlay) :: qb 32 32 real, intent(in), dimension(ngrid,nlay) :: pplay 33 33 real, intent(in), dimension(ngrid,nlay+1) :: paprs … … 40 40 41 41 real, intent(out), dimension(ngrid,nlay) :: dtemp_bs 42 real, intent(out), dimension(ngrid,nlay) :: dq _bs43 real, intent(out), dimension(ngrid,nlay) :: dqb s_bs42 real, intent(out), dimension(ngrid,nlay) :: dqv_bs 43 real, intent(out), dimension(ngrid,nlay) :: dqb_bs 44 44 real, intent(out), dimension(ngrid,nlay+1) :: bsfl 45 45 real, intent(out), dimension(ngrid) :: precip_bs … … 51 51 52 52 integer :: k,i,n 53 real :: cpd, cpw, dqsub, maxdqsub, dqbsmelt 54 real :: dqsedim,precbs, dqmelt, zmelt, taumeltbs 55 real :: maxdqmelt, rhoair, dz 56 real, dimension(ngrid) :: zt,zq,zqbs,qsi,dqsi,qsl, dqsl,qzero,sedim,sedimn 57 real, dimension(ngrid) :: zqbsi, zqbs_up, zmair 58 real, dimension(ngrid) :: zvelo 53 real :: cpd, cpw, dqbsub, maxdqbsub, dqbmelt, zmair 54 real :: dqsedim,precbs, dqvmelt, zmelt, taumeltbs 55 real :: maxdqvmelt, rhoair, dz, dqbsedim 56 real :: delta_p, b_p, a_p, c_p, c_sub, qvsub 57 real :: p0, T0, Dv, Aprime, Bprime, Ka 58 real, dimension(ngrid) :: ztemp,zqv,zqb,zpres,qsi,dqsi,qsl,dqsl,qzero,sedim 59 real, dimension(ngrid) :: zpaprsup, zpaprsdn, ztemp_up, zqb_up, zvelo 60 real, dimension(ngrid,nlay) :: temp_seri, qb_seri, qv_seri 59 61 60 62 !++++++++++++++++++++++++++++++++++++++++++++++++++ … … 64 66 qzero(:)=0. 65 67 dtemp_bs(:,:)=0. 66 dq _bs(:,:)=0.67 dqb s_bs(:,:)=0.68 dqv_bs(:,:)=0. 69 dqb_bs(:,:)=0. 68 70 zvelo(:)=0. 69 zt(:)=0.70 zq(:)=0.71 zqbs(:)=0.72 71 sedim(:)=0. 73 sedimn(:)=0. 74 72 precip_bs(:)=0. 73 bsfl(:,:)=0. 74 75 76 Ka=2.4e-2 ! thermal conductivity of the air, SI 77 p0=101325.0 ! ref pressure 78 79 80 DO k=1,nlay 81 DO i=1,ngrid 82 temp_seri(i,k)=temp(i,k) 83 qv_seri(i,k)=qv(i,k) 84 qb_seri(i,k)=qb(i,k) 85 ENDDO 86 ENDDO 87 88 89 ! Sedimentation scheme 90 !---------------------- 91 92 IF (iflag_sedim_bs .GT. 0) THEN 75 93 ! begin of top-down loop 76 94 DO k = nlay, 1, -1 77 95 78 96 DO i=1,ngrid 79 zt(i)=temp(i,k) 80 zq(i)=qv(i,k) 81 zqbs(i)=qbs(i,k) 97 ztemp(i)=temp_seri(i,k) 98 zqv(i)=qv_seri(i,k) 99 zqb(i)=qb_seri(i,k) 100 zpres(i)=pplay(i,k) 101 zpaprsup(i)=paprs(i,k+1) 102 zpaprsdn(i)=paprs(i,k) 82 103 ENDDO 83 104 84 105 ! thermalization of blowing snow precip coming from above 85 106 IF (k.LE.nlay-1) THEN 86 107 87 108 DO i = 1, ngrid 88 zmair (i)=(paprs(i,k)-paprs(i,k+1))/RG109 zmair=(zpaprsdn(i)-zpaprsup(i))/RG 89 110 ! RVTMP2=rcpv/rcpd-1 90 cpd=RCPD*(1.0+RVTMP2*zq (i))111 cpd=RCPD*(1.0+RVTMP2*zqv(i)) 91 112 cpw=RCPD*RVTMP2 92 ! zqb s_up: blowing snow mass that has to be thermalized with113 ! zqb_up: blowing snow mass that has to be thermalized with 93 114 ! layer's air so that precipitation at the ground has the 94 115 ! same temperature as the lowermost layer 95 zqbs_up(i) = sedim(i)*dtime/zmair(i) 116 zqb_up(i) = sedim(i)*dtime/zmair 117 ztemp_up(i)=temp_seri(i,k+1) 118 96 119 ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer 97 zt (i) = ( (temp(i,k+1)+dtemp_bs(i,k+1))*zqbs_up(i)*cpw + cpd*zt(i) ) &98 / (cpd + zqbs_up(i)*cpw)120 ztemp(i) = ( ztemp_up(i)*zqb_up(i)*cpw + cpd*ztemp(i) ) & 121 / (cpd + zqb_up(i)*cpw) 99 122 ENDDO 100 ENDIF 101 102 103 104 ! calulation saturation specific humidity 105 CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:)) 106 CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:)) 107 108 109 ! 3 processes: melting, sublimation and precipitation of blowing snow 123 124 ENDIF 125 126 DO i = 1, ngrid 127 128 rhoair = zpres(i) / ztemp(i) / RD 129 dz = (zpaprsdn(i)-zpaprsup(i)) / (rhoair*RG) 130 ! BS fall velocity assumed to be constant for now 131 zvelo(i) = fallv_bs 132 ! dqb/dt_sedim=1/rho(sedim/dz - rho w qb/dz) 133 ! implicit resolution 134 dqbsedim = (sedim(i)/rhoair/dz*dtime+zqb(i))/(1.+zvelo(i)*dtime/dz) - zqb(i) 135 ! flux and dqb update 136 zqb(i)=zqb(i)+dqbsedim 137 !sedim(i) = sedim(i)-dqbsedim/dtime*rhoair*dz 138 sedim(i)=rhoair*zvelo(i)*zqb(i) 139 140 141 ! variables update: 142 bsfl(i,k)=sedim(i) 143 144 qb_seri(i,k) = zqb(i) 145 qv_seri(i,k) = zqv(i) 146 temp_seri(i,k) = ztemp(i) 147 148 ENDDO ! Loop on ngrid 149 150 151 ENDDO ! vertical loop 152 153 154 !surface bs flux 155 DO i = 1, ngrid 156 precip_bs(i) = sedim(i) 157 ENDDO 158 159 ENDIF 160 161 162 163 164 !+++++++++++++++++++++++++++++++++++++++++++++++ 165 ! Sublimation and melting 166 !++++++++++++++++++++++++++++++++++++++++++++++ 167 IF (iflag_sublim_bs .GT. 0) THEN 168 169 DO k = 1, nlay 170 171 DO i=1,ngrid 172 ztemp(i)=temp_seri(i,k) 173 zqv(i)=qv_seri(i,k) 174 zqb(i)=qb_seri(i,k) 175 zpres(i)=pplay(i,k) 176 zpaprsup(i)=paprs(i,k+1) 177 zpaprsdn(i)=paprs(i,k) 178 ENDDO 179 180 ! calulation saturation specific humidity 181 CALL CALC_QSAT_ECMWF(ngrid,ztemp(:),qzero(:),zpres(:),RTT,2,.false.,qsi(:),dqsi(:)) 182 183 110 184 DO i = 1, ngrid 111 185 112 rhoair = pplay(i,k) / zt(i) / RD113 dz = ( paprs(i,k)-paprs(i,k+1)) / (rhoair*RG)186 rhoair = zpres(i) / ztemp(i) / RD 187 dz = (zpaprsdn(i)-zpaprsup(i)) / (rhoair*RG) 114 188 ! BS fall velocity assumed to be constant for now 115 189 zvelo(i) = fallv_bs 116 190 117 191 118 IF (zt (i) .GT. RTT) THEN119 192 IF (ztemp(i) .GT. RTT) THEN 193 120 194 ! if temperature is positive, we assume that part of the blowing snow 121 195 ! already present melts and evaporates with a typical time 122 196 ! constant taumeltbs 123 124 taumeltbs=taumeltbs0*exp(-max(0.,(zt (i)-RTT)/(tbsmelt-RTT)))125 dq melt=min(zqbs(i),-1.*zqbs(i)*(exp(-dtime/taumeltbs)-1.))126 maxdq melt= max(RCPD*(1.0+RVTMP2*(zq(i)+zqbs(i)))*(zt(i)-RTT)/(RLMLT+RLVTT),0.)127 dq melt=min(dqmelt,maxdqmelt)128 ! q update, melting + evaporation129 zq (i) = zq(i) + dqmelt197 198 taumeltbs=taumeltbs0*exp(-max(0.,(ztemp(i)-RTT)/(tbsmelt-RTT))) 199 dqvmelt=min(zqb(i),-1.*zqb(i)*(exp(-dtime/taumeltbs)-1.)) 200 maxdqvmelt= max(RCPD*(1.0+RVTMP2*(zqv(i)+zqb(i)))*(ztemp(i)-RTT)/(RLMLT+RLVTT),0.) 201 dqvmelt=min(dqvmelt,maxdqvmelt) 202 ! qv update, melting + evaporation 203 zqv(i) = zqv(i) + dqvmelt 130 204 ! temp update melting + evaporation 131 zt (i) = zt(i) - dqmelt * (RLMLT+RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))132 ! qb supdate melting + evaporation133 zqb s(i)=zqbs(i)-dqmelt134 135 ELSE 205 ztemp(i) = ztemp(i) - dqvmelt * (RLMLT+RLVTT)/RCPD/(1.0+RVTMP2*(zqv(i)+zqb(i))) 206 ! qb update melting + evaporation 207 zqb(i)=zqb(i)-dqvmelt 208 209 ELSE 136 210 ! negative celcius temperature 137 211 ! Sublimation scheme 138 139 rhoair=pplay(i,k)/zt(i)/RD 140 dqsub = 1./rhoair*coef_eva_bs*(1.0-zq(i)/qsi(i))*((rhoair*zvelo(i)*zqbs(i))**expo_eva_bs)*dtime 141 dqsub = MAX(0.0,MIN(dqsub,zqbs(i))) 142 143 ! Sublimation limit: we ensure that the whole mesh does not exceed saturation wrt ice 144 maxdqsub = MAX(0.0, qsi(i)-zq(i)) 145 dqsub = MIN(dqsub,maxdqsub) 146 212 213 214 ! Sublimation formulation for ice crystals from Pruppacher & Klett, Rutledge & Hobbs 1983 215 ! assuming monodispered crystal distrib 216 ! dqb/dt_sub=-coef_sub_bs*(1-qv/qsi)*nc*8*r_bs/(Aprime+Bprime) 217 ! assuming Mi=rhobs*4/3*pi*r_bs**3 218 ! rhoair qb=nc*Mi -> nc=rhoair qb/Mi 219 ! dqb/dt_sub=-coef_sub_bs*(1-qv/qsi)*6*rhoair*qb/(rhobs*pi*r_bs**2)/(Aprime+Bprime) 220 ! dqb/dt_sub=-c_sub(1-qv/qsi)*qb 221 ! c_sub=coef_sub_bs*6*rhoair/(rhobs*pi*r_bs**2)/(Aprime+Bprime) 222 ! 223 ! Note the strong coupling between specific contents of water vapor and blowing snow during sublimation 224 ! equations dqv/dt_sub and dqb/dt_sub must be solved jointly to prevent irrealistic increase of water vapor 225 ! at typical physics time steps 226 ! we thus solve the differential equation using an implicit scheme for both qb and qv 227 228 ! we do not consider deposition, only sublimation 229 IF (zqv(i) .LT. qsi(i)) THEN 230 rhoair=zpres(i)/ztemp(i)/RD 231 !diffusivity of water vapor 232 Dv=0.0001*0.211*(p0/zpres(i))*((ztemp(i)/RTT)**1.94) ! water vapor diffusivity in air, SI 233 Ka=(5.69+0.017*(ztemp(i)-RTT))*1.e-5*100.*4.184 ! thermal conductivity of the air, SI 234 Aprime=RLSTT/Ka/ztemp(i)*(RLSTT/RV/ztemp(i) -1.) 235 Bprime=1./(rhoair*Dv*qsi(i)) 236 c_sub=coef_sub_bs*6.*rhoair/(rhobs*RPI*r_bs**2)/(Aprime+Bprime) 237 c_p=-zqb(i) 238 b_p=1.+c_sub*dtime-c_sub*dtime/qsi(i)*zqb(i)-c_sub*dtime/qsi(i)*zqv(i) 239 a_p=c_sub*dtime/qsi(i) 240 delta_p=(b_p**2)-4.*a_p*c_p 241 dqbsub=(-b_p+sqrt(delta_p))/(2.*a_p) - zqb(i) 242 dqbsub = MIN(0.0,MAX(dqbsub,-zqb(i))) 243 244 ! Sublimation limit: we ensure that the whole mesh does not exceed saturation wrt ice 245 maxdqbsub = MAX(0.0, qsi(i)-zqv(i)) 246 dqbsub = MAX(dqbsub,-maxdqbsub) 247 ELSE 248 dqbsub=0. 249 ENDIF 147 250 148 251 ! vapor, temperature, precip fluxes update following sublimation 149 zq (i) = zq(i) + dqsub150 zqb s(i)=zqbs(i)-dqsub151 zt (i) = zt(i) - dqsub*RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))252 zqv(i) = zqv(i) - dqbsub 253 zqb(i) = zqb(i) + dqbsub 254 ztemp(i) = ztemp(i) + dqbsub*RLSTT/RCPD/(1.0+RVTMP2*(zqv(i)+zqb(i))) 152 255 153 256 ENDIF 154 257 155 ! Sedimentation scheme 156 !---------------------- 157 158 sedimn(i) = rhoair*zqbs(i)*zvelo(i) 159 160 ! exact numerical resolution of sedimentation 161 ! assuming fall velocity is constant 162 163 zqbs(i) = zqbs(i)*exp(-zvelo(i)/dz*dtime)+sedim(i)/rhoair/zvelo(i) 164 165 ! flux update 166 sedim(i) = sedimn(i) 167 168 169 170 ! if qbs<qbsmin, sublimate or melt and evaporate qbs 171 ! see Gerber et al. 2023, JGR Atmos for the choice of qbsmin 172 IF (zqbs(i) .LT. qbsmin) THEN 173 zq(i) = zq(i)+zqbs(i) 174 IF (zt(i) .LT. RTT) THEN 175 zt(i) = zt(i) - zqbs(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i))) 258 ! if qb<qbmin, sublimate or melt and evaporate qb 259 ! see Gerber et al. 2023, JGR Atmos for the choice of qbmin 260 261 IF (zqb(i) .LT. qbmin) THEN 262 zqv(i) = zqv(i)+zqb(i) 263 IF (ztemp(i) .LT. RTT) THEN 264 ztemp(i) = ztemp(i) - zqb(i) * RLSTT/RCPD/(1.0+RVTMP2*(zqv(i))) 176 265 ELSE 177 zt (i) = zt(i) - zqbs(i) * (RLVTT+RLMLT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))266 ztemp(i) = ztemp(i) - zqb(i) * (RLVTT+RLMLT)/RCPD/(1.0+RVTMP2*(zqv(i))) 178 267 ENDIF 179 zqbs(i)=0. 180 ENDIF 181 182 183 ! Outputs: 184 bsfl(i,k)=sedim(i) 185 dqbs_bs(i,k) = zqbs(i)-qbs(i,k) 186 dq_bs(i,k) = zq(i) - qv(i,k) 187 dtemp_bs(i,k) = zt(i) - temp(i,k) 188 189 ENDDO ! Loop on ngrid 190 191 192 ENDDO ! vertical loop 193 194 195 !surface bs flux 196 DO i = 1, ngrid 197 precip_bs(i) = sedim(i) 198 ENDDO 268 zqb(i)=0. 269 ENDIF 270 271 ! variables update 272 temp_seri(i,k)=ztemp(i) 273 qv_seri(i,k)=zqv(i) 274 qb_seri(i,k)=zqb(i) 275 ENDDO 276 ENDDO 277 278 ENDIF 279 280 281 ! OUTPUTS 282 !++++++++++ 283 284 ! 3D variables 285 DO k=1,nlay 286 DO i=1,ngrid 287 dqb_bs(i,k) = qb_seri(i,k) - qb(i,k) 288 dqv_bs(i,k) = qv_seri(i,k) - qv(i,k) 289 dtemp_bs(i,k) = temp_seri(i,k) - temp(i,k) 290 ENDDO 291 ENDDO 292 199 293 200 294
Note: See TracChangeset
for help on using the changeset viewer.
