Changeset 4821
- Timestamp:
- Feb 14, 2024, 9:26:23 PM (11 months ago)
- Location:
- LMDZ6/trunk/libf
- Files:
-
- 4 edited
-
phylmd/lmdz_blowing_snow_ini.F90 (modified) (4 diffs)
-
phylmd/lmdz_blowing_snow_sublim_sedim.F90 (modified) (5 diffs)
-
phylmd/physiq_mod.F90 (modified) (1 diff)
-
phylmdiso/physiq_mod.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_ini.F90
r4724 r4821 3 3 implicit none 4 4 5 save 6 integer :: prt_level=0,lunout 7 real RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG 8 real coef_eva_bs, expo_eva_bs, fallv_bs, zeta_bs 9 real prt_bs, pbst_bs, qbst_bs 10 integer :: iflag_saltation_bs 5 real, save, protected :: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG 6 real, save, protected :: coef_eva_bs, expo_eva_bs, fallv_bs, zeta_bs 7 real, save, protected :: prt_bs, pbst_bs, qbst_bs 8 9 integer, save, protected :: iflag_saltation_bs 11 10 12 !$OMP THREADPRIVATE(prt_level, lunout)13 11 !$OMP THREADPRIVATE(RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG) 14 12 !$OMP THREADPRIVATE(coef_eva_bs, expo_eva_bs, fallv_bs, zeta_bs) … … 16 14 !$OMP THREADPRIVATE(iflag_saltation_bs) 17 15 18 real tbsmelt ! parameter to calculate melting fraction of BS sedimentation 19 parameter (tbsmelt=278.15) 16 real, save, protected :: tbsmelt=278.15 ! parameter to calculate melting fraction of BS sedimentation 17 real, save, protected :: taumeltbs0=1800.0 ! Melting time scale of blowing snow at 273.15K 18 real, save, protected :: qbsmin=1.E-10 ! Minimum blowing snow specific content 20 19 21 real taumeltbs0 ! Melting time scale of blowing snow at 273.15K 22 parameter (taumeltbs0=1800.0) 23 24 real qbsmin ! Minimum blowing snow specific content 25 parameter (qbsmin=1.E-10) 20 !$OMP THREADPRIVATE(tbsmelt, taumeltbs0, qbsmin) 26 21 27 22 28 23 contains 29 24 30 subroutine blowing_snow_ini(prt_level_in,lunout_in, & 31 RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in,& 25 subroutine blowing_snow_ini(RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in,& 32 26 RVTMP2_in, RTT_in,RD_in,RG_in) 33 27 34 28 USE ioipsl_getin_p_mod, ONLY : getin_p 35 29 36 integer, intent(in) :: prt_level_in,lunout_in37 30 real, intent(in) :: RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in 38 31 real, intent(in) :: RVTMP2_in, RTT_in, RD_in, RG_in 39 32 40 33 41 prt_level=prt_level_in42 34 RG=RG_in 43 35 RD=RD_in … … 49 41 RG=RG_in 50 42 RVTMP2=RVTMP2_in 51 lunout=lunout_in52 43 53 44 … … 60 51 prt_bs= 0.0003 61 52 CALL getin_p('prt_bs',prt_bs) 62 63 53 64 54 zeta_bs= 3. -
LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.F90
r4724 r4821 51 51 52 52 integer :: k,i,n 53 real :: zqev0, zqevi, zqevti, zcpair, zcpeau, dqbsmelt54 real :: dqsedim,precbs, d eltaqchaud, zmelt, taumeltbs55 real :: maxd eltaqchaud53 real :: cpd, cpw, dqsub, maxdqsub, dqbsmelt 54 real :: dqsedim,precbs, dqmelt, zmelt, taumeltbs 55 real :: maxdqmelt, rhoair, dz 56 56 real, dimension(ngrid) :: zt,zq,zqbs,qsi,dqsi,qsl, dqsl,qzero,sedim,sedimn 57 real, dimension(ngrid) :: zqbsi, zmqc, zmair, zdz58 real, dimension(ngrid ,nlay) :: velo, zrho57 real, dimension(ngrid) :: zqbsi, zqbs_up, zmair 58 real, dimension(ngrid) :: zvelo 59 59 60 60 !++++++++++++++++++++++++++++++++++++++++++++++++++ … … 66 66 dq_bs(:,:)=0. 67 67 dqbs_bs(:,:)=0. 68 velo(:,:)=0.68 zvelo(:)=0. 69 69 zt(:)=0. 70 70 zq(:)=0. 71 71 zqbs(:)=0. 72 72 sedim(:)=0. 73 sedimn(:)=0. 73 74 74 75 ! begin of top-down loop … … 81 82 ENDDO 82 83 83 84 ! thermalization of blowing snow precip coming from above 84 85 IF (k.LE.nlay-1) THEN 85 86 86 ! thermalization of blowing snow precip coming from above87 87 DO i = 1, ngrid 88 88 zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG 89 89 ! RVTMP2=rcpv/rcpd-1 90 zcpair=RCPD*(1.0+RVTMP2*zq(i))91 zcpeau=RCPD*RVTMP292 ! z mqc: precipitationmass that has to be thermalized with90 cpd=RCPD*(1.0+RVTMP2*zq(i)) 91 cpw=RCPD*RVTMP2 92 ! zqbs_up: blowing snow mass that has to be thermalized with 93 93 ! layer's air so that precipitation at the ground has the 94 94 ! same temperature as the lowermost layer 95 z mqc(i) = (sedim(i))*dtime/zmair(i)95 zqbs_up(i) = sedim(i)*dtime/zmair(i) 96 96 ! 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))*z mqc(i)*zcpeau + zcpair*zt(i) ) &98 / ( zcpair + zmqc(i)*zcpeau)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) 99 99 ENDDO 100 ELSE101 DO i = 1, ngrid102 zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG103 zmqc(i) = 0.104 ENDDO105 106 100 ENDIF 107 101 … … 111 105 CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:)) 112 106 CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:)) 113 ! sublimation calculation114 ! SUndqvist formula dP/dz=beta*(1-q/qsat)*sqrt(P)115 107 108 109 ! 3 processes: melting, sublimation and precipitation of blowing snow 116 110 DO i = 1, ngrid 117 118 zrho(i,k) = pplay(i,k) / zt(i) / RD119 zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)120 ! BS fall velocity121 velo(i,k) = fallv_bs122 123 111 124 112 IF (zt(i) .GT. RTT) THEN 125 113 126 ! if positive celcius temperature, we assume 127 ! that part of the the blowing snow flux melts and evaporates 128 129 ! vapor, bs, temperature, precip fluxes update 130 zmelt = ((zt(i)-RTT)/(tbsmelt-RTT)) 131 zmelt = MIN(MAX(zmelt,0.),1.) 132 sedimn(i)=sedim(i)*zmelt 133 deltaqchaud=-(sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime 134 ! max evap such as celcius temperature cannot become negative 135 maxdeltaqchaud= max(RCPD*(1.0+RVTMP2*(zq(i)+zqbs(i)))*(zt(i)-RTT)/(RLMLT+RLVTT),0.) 136 137 deltaqchaud=min(max(deltaqchaud,0.),maxdeltaqchaud) 138 zq(i) = zq(i) + deltaqchaud 139 140 ! melting + evaporation 141 zt(i) = zt(i) - deltaqchaud * (RLMLT+RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i))) 142 143 sedim(i)=sedimn(i) 144 145 ! if temperature still positive, we assume that part of the blowing snow 146 ! already present in the mesh melts and evaporates with a typical time 147 ! constant between taumeltbs0 and 0 148 149 IF ( zt(i) .GT. RTT ) THEN 150 taumeltbs=taumeltbs0*exp(-max(0.,(zt(i)-RTT)/(tbsmelt-RTT))) 151 deltaqchaud=min(zqbs(i),zqbs(i)/taumeltbs*dtime) 152 maxdeltaqchaud= max(RCPD*(1.0+RVTMP2*(zq(i)+zqbs(i)))*(zt(i)-RTT)/(RLMLT+RLVTT),0.) 153 deltaqchaud=min(max(deltaqchaud,0.),maxdeltaqchaud) 154 zq(i) = zq(i) + deltaqchaud 155 156 ! melting + evaporation 157 zt(i) = zt(i) - deltaqchaud * (RLMLT+RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i))) 158 ! qbs update 159 zqbs(i)=max(zqbs(i)-deltaqchaud,0.) 160 ENDIF 161 162 ! now sedimentation scheme with an exact numerical resolution 163 ! (correct if fall velocity is constant) 164 zqbsi(i)=zqbs(i) 165 zqbs(i) = zqbsi(i)*exp(-velo(i,k)/zdz(i)*dtime)+sedim(i)/zrho(i,k)/velo(i,k) 166 zqbs(i) = max(zqbs(i),0.) 167 168 ! flux update 169 sedim(i) = sedim(i) + zrho(i,k)*zdz(i)/dtime*(zqbsi(i)-zqbs(i)) 170 sedim(i) = max(0.,sedim(i)) 171 114 ! if temperature is positive, we assume that part of the blowing snow 115 ! already present melts and evaporates with a typical time 116 ! constant taumeltbs 117 118 taumeltbs=taumeltbs0*exp(-max(0.,(zt(i)-RTT)/(tbsmelt-RTT))) 119 dqmelt=min(zqbs(i),zqbs(i)/taumeltbs*dtime) 120 maxdqmelt= max(RCPD*(1.0+RVTMP2*(zq(i)+zqbs(i)))*(zt(i)-RTT)/(RLMLT+RLVTT),0.) 121 dqmelt=min(max(dqmelt,0.),maxdqmelt) 122 ! q update, melting + evaporation 123 zq(i) = zq(i) + dqmelt 124 ! temp update melting + evaporation 125 zt(i) = zt(i) - dqmelt * (RLMLT+RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i))) 126 ! qbs update melting + evaporation 127 zqbs(i)=zqbs(i)-dqmelt 128 172 129 ELSE 173 130 ! negative celcius temperature 174 131 ! Sublimation scheme 175 132 176 zqevti = coef_eva_bs*(1.0-zq(i)/qsi(i))*(sedim(i)**expo_eva_bs) &177 *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG178 zqevti = MAX(0.0,MIN(zqevti,sedim(i)))*RG*dtime/(paprs(i,k)-paprs(i,k+1))133 rhoair=pplay(i,k)/zt(i)/RD 134 dqsub = 1./rhoair*coef_eva_bs*(1.0-zq(i)/qsi(i))*((rhoair*zvelo(i)*zqbs(i))**expo_eva_bs)*dtime 135 dqsub = MAX(0.0,MIN(dqsub,zqbs(i))) 179 136 180 137 ! Sublimation limit: we ensure that the whole mesh does not exceed saturation wrt ice 181 zqev0 = MAX(0.0, qsi(i)-zq(i)) 182 zqevi = MIN(zqev0,zqevti) 183 184 ! New solid precipitation fluxes 185 sedimn(i) = Max(0.,sedim(i) - zqevi*(paprs(i,k)-paprs(i,k+1))/RG/dtime) 138 maxdqsub = MAX(0.0, qsi(i)-zq(i)) 139 dqsub = MIN(dqsub,maxdqsub) 186 140 187 141 188 142 ! vapor, temperature, precip fluxes update following sublimation 189 zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime 190 zq(i) = max(0., zq(i)) 191 zt(i) = zt(i) & 192 + (sedimn(i)-sedim(i)) & 193 * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & 194 * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i))) 195 196 sedim(i)=sedimn(i) 197 zqbsi(i)=zqbs(i) 198 199 ! now sedimentation scheme with an exact numerical resolution 200 ! (correct if fall velocity is constant) 201 202 zqbs(i) = zqbsi(i)*exp(-velo(i,k)/zdz(i)*dtime)+sedim(i)/zrho(i,k)/velo(i,k) 203 zqbs(i) = max(zqbs(i),0.) 204 205 ! flux update 206 sedim(i) = sedim(i) + zrho(i,k)*zdz(i)/dtime*(zqbsi(i)-zqbs(i)) 207 sedim(i) = max(0.,sedim(i)) 208 209 210 ENDIF 211 212 213 ! if qbs<qbsmin, sublimate or melt and evaporate qbs 214 ! see Gerber et al. 2023, JGR Atmos for the choice of qbsmin 215 IF (zqbs(i) .LT. qbsmin) THEN 143 zq(i) = zq(i) + dqsub 144 zqbs(i)=zqbs(i)-dqsub 145 zt(i) = zt(i) - dqsub*RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i))) 146 147 ENDIF 148 149 ! Sedimentation scheme 150 151 rhoair = pplay(i,k) / zt(i) / RD 152 dz = (paprs(i,k)-paprs(i,k+1)) / (rhoair*RG) 153 ! BS fall velocity assumed to be constant for now 154 zvelo(i) = fallv_bs 155 156 sedimn(i) = rhoair*zqbs(i)*zvelo(i) 157 158 ! exact numerical resolution of sedimentation 159 ! assuming fall velocity is constant 160 161 zqbs(i) = zqbs(i)*exp(-zvelo(i)/dz*dtime)+sedim(i)/rhoair/zvelo(i) 162 163 ! flux update 164 sedim(i) = sedimn(i) 165 166 167 168 ! if qbs<qbsmin, sublimate or melt and evaporate qbs 169 ! see Gerber et al. 2023, JGR Atmos for the choice of qbsmin 170 IF (zqbs(i) .LT. qbsmin) THEN 216 171 zq(i) = zq(i)+zqbs(i) 217 172 IF (zt(i) .LT. RTT) THEN … … 221 176 ENDIF 222 177 zqbs(i)=0. 223 ENDIF 224 225 226 ENDDO ! loop on ngrid 227 228 178 ENDIF 229 179 230 180 231 181 ! Outputs: 232 DO i = 1, ngrid233 182 bsfl(i,k)=sedim(i) 234 183 dqbs_bs(i,k) = zqbs(i)-qbs(i,k) 235 184 dq_bs(i,k) = zq(i) - qv(i,k) 236 185 dtemp_bs(i,k) = zt(i) - temp(i,k) 237 ENDDO 186 187 ENDDO ! Loop on ngrid 238 188 239 189 -
LMDZ6/trunk/libf/phylmd/physiq_mod.F90
r4819 r4821 1859 1859 CALL ratqs_ini(klon,klev,iflag_thermals,lunout,nbsrf,is_lic,is_ter,RG,RV,RD,RCPD,RLSTT,RLVTT,RTT) 1860 1860 CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_sursat,iflag_ratqs,fl_cor_ebil,RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RG,RPI) 1861 CALL blowing_snow_ini(prt_level,lunout, & 1862 RCPD, RLSTT, RLVTT, RLMLT, & 1861 CALL blowing_snow_ini(RCPD, RLSTT, RLVTT, RLMLT, & 1863 1862 RVTMP2, RTT,RD,RG) 1864 1863 ! Test de coherence sur oc_cdnc utilisé uniquement par cloud_optics_prop -
LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
r4819 r4821 1951 1951 CALL ratqs_ini(klon,klev,iflag_thermals,lunout,nbsrf,is_lic,is_ter,RG,RV,RD,RCPD,RLSTT,RLVTT,RTT) 1952 1952 CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_sursat,iflag_ratqs,fl_cor_ebil,RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RG,RPI) 1953 CALL blowing_snow_ini(prt_level,lunout, & 1954 RCPD, RLSTT, RLVTT, RLMLT, & 1953 CALL blowing_snow_ini(RCPD, RLSTT, RLVTT, RLMLT, & 1955 1954 RVTMP2, RTT,RD,RG) 1956 1955 ! Test de coherence sur oc_cdnc utilisé uniquement par cloud_optics_prop
Note: See TracChangeset
for help on using the changeset viewer.
