Changeset 4835
- Timestamp:
- Feb 29, 2024, 7:42:12 PM (2 months ago)
- Location:
- LMDZ6/trunk
- Files:
-
- 12 edited
-
DefLists/field_def_lmdz.xml (modified) (1 diff)
-
libf/phylmd/lmdz_blowing_snow_ini.F90 (modified) (4 diffs)
-
libf/phylmd/lmdz_blowing_snow_sublim_sedim.F90 (modified) (6 diffs)
-
libf/phylmd/phys_local_var_mod.F90 (modified) (3 diffs)
-
libf/phylmd/phys_output_ctrlout_mod.F90 (modified) (1 diff)
-
libf/phylmd/phys_output_write_mod.F90 (modified) (3 diffs)
-
libf/phylmd/physiq_mod.F90 (modified) (1 diff)
-
libf/phylmd/surf_landice_mod.F90 (modified) (5 diffs)
-
libf/phylmdiso/phys_local_var_mod.F90 (modified) (2 diffs)
-
libf/phylmdiso/phys_output_ctrlout_mod.F90 (modified) (1 diff)
-
libf/phylmdiso/physiq_mod.F90 (modified) (1 diff)
-
libf/phylmdiso/surf_landice_mod.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/DefLists/field_def_lmdz.xml
r4820 r4835 117 117 <field id="rhosnow_lic" long_name="snow density lic" unit="kg/m3" /> 118 118 <field id="ustart_lic" long_name="ustar threshold for snow erosion" unit="m/s" /> 119 <field id="qsalt_lic" long_name="blowing snow concentration in saltation layer" unit="kg/kg" /> 119 120 <field id="evap_ter" long_name="evaporation at surface ter" unit="kg/(s*m2)" /> 120 121 <field id="evap_lic" long_name="evaporation at surface lic" unit="kg/(s*m2)" /> -
LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_ini.F90
r4821 r4835 3 3 implicit none 4 4 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 5 real, save, protected :: RCPD, RV, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG, RPI 6 real, save, protected :: coef_sub_bs, fallv_bs, zeta_bs, c_esalt_bs 7 real, save, protected :: prt_bs, pbst_bs, qbst_bs, r_bs 8 integer, save, protected :: iflag_saltation_bs, iflag_sedim_bs, iflag_sublim_bs 10 9 11 !$OMP THREADPRIVATE(RCPD, R LSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG)12 !$OMP THREADPRIVATE(coef_ eva_bs, expo_eva_bs, fallv_bs, zeta_bs)10 !$OMP THREADPRIVATE(RCPD, RV, RLSTT, RLVTT, RLMLT, RVTMP2, RTT,RD,RG, RPI) 11 !$OMP THREADPRIVATE(coef_sub_bs, fallv_bs, r_bs, zeta_bs, c_esalt_bs) 13 12 !$OMP THREADPRIVATE(pbst_bs, prt_bs, qbst_bs) 14 !$OMP THREADPRIVATE(iflag_s altation_bs)13 !$OMP THREADPRIVATE(iflag_sedim_bs, iflag_sublim_bs) 15 14 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 15 real, save, protected :: tbsmelt=278.15 ! parameter to calculate melting fraction of BS sedimentation 16 real, save, protected :: taumeltbs0=600.0 ! Melting time scale of blowing snow at 273.15K 17 real, save, protected :: qbmin=1.E-10 ! Minimum blowing snow specific content 18 !$OMP THREADPRIVATE(tbsmelt, taumeltbs0, qbmin) 19 19 20 !$OMP THREADPRIVATE(tbsmelt, taumeltbs0, qbsmin) 20 real, save, protected :: tau_dens0_bs=864000. ! 10 days by default, in s 21 real, save, protected :: tau_densmin_bs= 21600. ! 1/4 days according to in situ obs by C. Amory during blowing snow + 22 ! Marshall et al. 1999 (snow densification during rain) 23 real, save, protected :: tau_eqsalt_bs= 10. ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt of about 10s 24 real, save, protected :: rhofresh_bs = 300.0 ! fresh snow density kg/m3 25 real, save, protected :: rhohard_bs = 450.0 ! hard snow density kg/m3 26 real, save, protected :: rhoice_bs = 920.0 ! ice density kg/m3 27 real, save, protected :: rhobs=900.0 ! blowing snow density (kg/m3) following Bintanja et al. 2001 part I 28 !$OMP THREADPRIVATE(rhoice_bs, rhofresh_bs, rhohard_bs, tau_dens0_bs, tau_densmin_bs, tau_eqsalt_bs, rhobs) 21 29 22 30 … … 24 32 25 33 subroutine blowing_snow_ini(RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in,& 26 RVTMP2_in, RTT_in,RD_in,RG_in )34 RVTMP2_in, RTT_in,RD_in,RG_in, RV_in, RPI_in) 27 35 28 36 USE ioipsl_getin_p_mod, ONLY : getin_p 29 37 30 real, intent(in) :: RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in 31 real, intent(in) :: RVTMP2_in, RTT_in, RD_in, RG_in 38 real, intent(in) :: RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in, RPI_in 39 real, intent(in) :: RVTMP2_in, RTT_in, RD_in, RG_in, RV_in 32 40 33 41 34 42 RG=RG_in 35 43 RD=RD_in 44 RV=RV_in 36 45 RCPD=RCPD_in 37 46 RLVTT=RLVTT_in … … 41 50 RG=RG_in 42 51 RVTMP2=RVTMP2_in 52 RPI=RPI_in 43 53 54 c_esalt_bs= 3.25 55 CALL getin_p('c_esalt_bs',c_esalt_bs) 44 56 45 57 qbst_bs= 0.001 … … 58 70 CALL getin_p('fallv_bs',fallv_bs) 59 71 60 coef_ eva_bs = 2.e-561 CALL getin_p('coef_ eva_bs',coef_eva_bs)72 coef_sub_bs = 0.1 73 CALL getin_p('coef_sub_bs',coef_sub_bs) 62 74 63 expo_eva_bs = 0.564 CALL getin_p(' expo_eva_bs',expo_eva_bs)75 iflag_sublim_bs=1 76 CALL getin_p('iflag_sublim_bs',iflag_sublim_bs) 65 77 66 iflag_s altation_bs=067 CALL getin_p('iflag_s altation_bs',iflag_saltation_bs)78 iflag_sedim_bs=1 79 CALL getin_p('iflag_sedim_bs',iflag_sedim_bs) 68 80 81 r_bs=150.0e-6 82 CALL getin_p('r_bs',r_bs) 69 83 70 84 end subroutine blowing_snow_ini -
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 -
LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90
r4830 r4835 335 335 REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: tpot, tpote, ue, uq, uwat, ve, vq, vwat, zxffonte 336 336 !$OMP THREADPRIVATE(tpot, tpote, ue, uq, uwat, ve, vq, vwat, zxffonte) 337 REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zxustartlic, zxrhoslic 338 !$OMP THREADPRIVATE(zxustartlic, zxrhoslic )337 REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zxustartlic, zxrhoslic, zxqsaltlic 338 !$OMP THREADPRIVATE(zxustartlic, zxrhoslic, zxqsaltlic) 339 339 REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zxfqcalving 340 340 !$OMP THREADPRIVATE(zxfqcalving) … … 842 842 ALLOCATE(zxtsol(klon), snow_lsc(klon), zxfqfonte(klon), zxqsurf(klon)) 843 843 ALLOCATE(zxrunofflic(klon)) 844 ALLOCATE(zxustartlic(klon), zxrhoslic(klon) )845 zxustartlic(:)=0. ; zxrhoslic(:)=0. 844 ALLOCATE(zxustartlic(klon), zxrhoslic(klon), zxqsaltlic(klon)) 845 zxustartlic(:)=0. ; zxrhoslic(:)=0. ; zxqsaltlic(:)=0. 846 846 ALLOCATE(rain_lsc(klon)) 847 847 ALLOCATE(rain_num(klon)) … … 1174 1174 DEALLOCATE(zxfqcalving, zxfluxlat) 1175 1175 DEALLOCATE(zxrunofflic) 1176 DEALLOCATE(zxustartlic, zxrhoslic )1176 DEALLOCATE(zxustartlic, zxrhoslic, zxqsaltlic) 1177 1177 DEALLOCATE(zxtsol, snow_lsc, zxfqfonte, zxqsurf) 1178 1178 DEALLOCATE(rain_lsc) -
LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90
r4819 r4835 388 388 TYPE(ctrl_out), SAVE :: o_rhosnow_lic = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), & 389 389 'rhosnow_lic', 'snow density lic', 'kg/m3', (/ ('', i=1, 10) /)) 390 TYPE(ctrl_out), SAVE :: o_qsalt_lic = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), & 391 'qsalt_lic', 'qb in saltation layer lic', 'kg/kg', (/ ('', i=1, 10) /)) 390 392 TYPE(ctrl_out), SAVE :: o_sens_prec_liq_oce = ctrl_out((/ 5, 5, 10, 10, 5, 10, 11, 11, 11, 11/), & 391 393 'sens_rain_oce', 'Sensible heat flux of liquid prec. over ocean', 'W/m2', (/ ('', i=1, 10) /)) -
LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90
r4830 r4835 48 48 o_psol, o_mass, o_qsurf, o_qsol, & 49 49 o_precip, o_rain_fall, o_rain_con, o_ndayrain, o_plul, o_pluc, o_plun, & 50 o_snow, o_msnow, o_fsnow, o_evap, o_snowerosion, o_ustart_lic, o_ rhosnow_lic, o_bsfall, &50 o_snow, o_msnow, o_fsnow, o_evap, o_snowerosion, o_ustart_lic, o_qsalt_lic, o_rhosnow_lic, o_bsfall, & 51 51 o_ep,o_epmax_diag, & ! epmax_cape 52 52 o_tops, o_tops0, o_topl, o_topl0, & … … 310 310 USE phys_local_var_mod, ONLY: zxfluxlat, slp, ptstar, pt0, zxtsol, zt2m, & 311 311 zn2mout, t2m_min_mon, t2m_max_mon, evap, & 312 snowerosion, zxustartlic, zxrhoslic, &312 snowerosion, zxustartlic, zxrhoslic, zxqsaltlic, & 313 313 l_mixmin,l_mix, tke_dissip, & 314 314 zu10m, zv10m, zq2m, zustar, zxqsurf, & … … 897 897 CALL histwrite_phy(o_ustart_lic, zxustartlic) 898 898 CALL histwrite_phy(o_rhosnow_lic, zxrhoslic) 899 CALL histwrite_phy(o_qsalt_lic, zxqsaltlic) 899 900 ENDIF 900 901 -
LMDZ6/trunk/libf/phylmd/physiq_mod.F90
r4830 r4835 1864 1864 CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_sursat,iflag_ratqs,fl_cor_ebil,RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RG,RPI) 1865 1865 CALL blowing_snow_ini(RCPD, RLSTT, RLVTT, RLMLT, & 1866 RVTMP2, RTT,RD,RG )1866 RVTMP2, RTT,RD,RG, RV, RPI) 1867 1867 ! Test de coherence sur oc_cdnc utilisé uniquement par cloud_optics_prop 1868 1868 IF (ok_newmicro) then -
LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90
r4828 r4835 31 31 USE cpl_mod, ONLY : cpl_send_landice_fields 32 32 USE calcul_fluxs_mod 33 USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic 33 USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic, zxqsaltlic 34 34 USE phys_output_var_mod, ONLY : snow_o,zfra_o 35 35 !FC 36 36 USE ioipsl_getin_p_mod, ONLY : getin_p 37 USE lmdz_blowing_snow_ini, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs38 37 USE lmdz_blowing_snow_ini, ONLY : c_esalt_bs, zeta_bs, pbst_bs, prt_bs, rhoice_bs, rhohard_bs 38 USE lmdz_blowing_snow_ini, ONLY : rhofresh_bs, tau_eqsalt_bs, tau_dens0_bs, tau_densmin_bs 39 39 #ifdef CPP_INLANDSIS 40 40 USE surf_inlandsis_mod, ONLY : surf_inlandsis … … 140 140 REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow 141 141 REAL, DIMENSION (klon,6) :: alb6 142 REAL :: rho0, rhoice, ustart0,esalt, rhod142 REAL :: esalt 143 143 REAL :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc 144 REAL :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard 145 REAL, DIMENSION(klon) :: ws1, rhos, ustart, qsalt, hsalt, snowini 146 REAL :: ustarmin, maxerosion, tau_eq_salt 144 REAL :: tau_dens, maxerosion, fluxbs_2 145 REAL, DIMENSION(klon) :: ws1, rhod, rhos, ustart0, ustart, qsalt, hsalt 147 146 148 147 ! End definition … … 331 330 ! 332 331 !**************************************************************************************** 333 334 CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)) 335 336 337 ! EV: following lines are obsolete since we set alb1 and alb2 to constant values 338 ! I therefore comment them 339 ! alb1(1:knon) = alb_neig(1:knon)*zfra(1:knon) + & 340 ! 0.6 * (1.0-zfra(1:knon)) 332 341 333 ! 342 334 !IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux" … … 361 353 362 354 363 364 ! Simple blowing snow param 365 if (ok_bs) then 366 ustarmin=1.e-3 367 ustart0 = 0.211 368 rhoice = 920.0 369 rho0 = 300.0 370 rhomax=450.0 371 rhohard=450.0 372 tau_dens0=86400.0*10. ! 10 days by default, in s 373 tau_densmin=21600.0 ! 1/4 days according to in situ obs by C. Amory during blowing snow + 374 ! Marshall et al. 1999 (snow densification during rain) 375 tau_eq_salt=10. 376 355 !**************************************************************************************** 356 ! Simple blowing snow param 357 !**************************************************************************************** 358 ! we proceed in 2 steps: 359 ! first we erode - if possible -the accumulated snow during the time step 360 ! then we update the density of the underlying layer and see if we can also erode 361 ! this layer 362 363 364 if (ok_bs) then 365 fluxbs(:)=0. 366 do j=1,klon 367 ws1(j)=(u1(j)**2+v1(j)**2)**0.5 368 ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5 369 rhod(j)=p1lay(j)/RD/temp_air(j) 370 ustart0(j) =(log(2.868)-log(1.625))/0.085*sqrt(cdragm(j)) 371 enddo 372 373 ! 1st step: erosion of fresh snow accumulated during the time step 374 do j=1, knon 375 376 rhos(j)=rhofresh_bs 377 ! blowing snow flux formula used in MAR 378 ustart(j)=ustart0(j)*exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))*exp(max(0.,rhos(j)-rhohard_bs)) 379 ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs 380 ! computation of qbs at the top of the saltation layer 381 ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001) 382 esalt=1./(c_esalt_bs*max(1.e-6,ustar(j))) 383 hsalt(j)=0.08436*(max(1.e-6,ustar(j))**1.27) 384 qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt 385 ! calculation of erosion (flux positive towards the surface here) 386 ! consistent with implicit resolution of turbulent mixing equation 387 ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s 388 ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level) 389 ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer 390 ! (rho*qsalt*hsalt) 391 ! during this first step we also lower bound the erosion to the amount of fresh snow accumulated during the time step 392 maxerosion=min(precip_snow(j),hsalt(j)*qsalt(j)*rhod(j)/tau_eqsalt_bs) 393 394 fluxbs(j)=rhod(j)*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) & 395 / (1.-rhod(j)*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime) 396 !fluxbs(j)= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j)) 397 fluxbs(j)=max(-maxerosion,fluxbs(j)) 398 enddo 399 400 401 ! we now compute the snow age of the overlying layer (snow surface after erosion of the fresh snow accumulated during the time step) 402 ! this is done through the routine albsno 403 CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)+fluxbs(:)) 404 405 ! 2nd step: 377 406 ! computation of threshold friction velocity 378 407 ! which depends on surface snow density … … 381 410 ! snow density increases with snow age and 382 411 ! increases even faster in case of sedimentation of blowing snow or rain 383 tau_dens=max(tau_densmin , tau_dens0*exp(-abs(precip_bs(j))/pbst_bs - &384 abs(precip_rain(j))/prt_bs)*exp(-max(tsurf(j)- 272.0,0.)))385 rhos(j)=rho 0+(rhohard-rho0)*(1.-exp(-agesno(j)*86400.0/tau_dens))412 tau_dens=max(tau_densmin_bs, tau_dens0_bs*exp(-abs(precip_bs(j))/pbst_bs - & 413 abs(precip_rain(j))/prt_bs)*exp(-max(tsurf(j)-RTT,0.))) 414 rhos(j)=rhofresh_bs+(rhohard_bs-rhofresh_bs)*(1.-exp(-agesno(j)*86400.0/tau_dens)) 386 415 ! blowing snow flux formula used in MAR 387 ws1(j)=(u1(j)**2+v1(j)**2)**0.5 388 ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5 389 ustart(j)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(j),0.))*exp(max(0.,rhos(j)-rhomax)) 390 ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till 391 ! rhohard<450) 416 ustart(j)=ustart0(j)*exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))*exp(max(0.,rhos(j)-rhohard_bs)) 417 ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs 418 ! computation of qbs at the top of the saltation layer 419 ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001) 420 esalt=c_esalt_bs*ustar(j) 421 hsalt(j)=0.08436*(max(1.e-6,ustar(j))**1.27) 422 qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt 423 ! calculation of erosion (flux positive towards the surface here) 424 ! consistent with implicit resolution of turbulent mixing equation 425 ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s 426 ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level) 427 ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer 428 ! (rho*qsalt*hsalt) 429 maxerosion=hsalt(j)*qsalt(j)*rhod(j)/tau_eqsalt_bs 430 fluxbs_2=rhod(j)*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) & 431 / (1.-rhod(j)*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime) 432 !fluxbs_2= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j)) 433 fluxbs_2=max(-maxerosion,fluxbs_2) 434 fluxbs(j)=fluxbs(j)+fluxbs_2 392 435 enddo 393 394 ! computation of qbs at the top of the saltation layer 395 ! two formulations possible 396 ! pay attention that qbs is a mixing ratio and has to be converted 397 ! to specific content 398 399 if (iflag_saltation_bs .eq. 1) then 400 ! expression from CRYOWRF (Sharma et al. 2022) 401 aa=2.6 402 bb=2.5 403 cc=2.0 404 lambdasalt=0.45 405 do j =1, knon 406 rhod=p1lay(j)/RD/temp_air(j) 407 nunu=max(ustar(j)/ustart(j),1.e-3) 408 fluxsalt=rhod/RG*(ustar(j)**3)*(1.-nunu**(-2)) * & 409 (aa+bb*nunu**(-2)+cc*nunu**(-1)) 410 csalt=fluxsalt/(2.8*ustart(j)) 411 hsalt(j)=0.08436*ustar(j)**1.27 412 qsalt(j)=1./rhod*csalt*lambdasalt*RG/(max(ustar(j)**2,ustarmin)) & 413 * exp(-lambdasalt*RG*hsalt(j)/max(ustar(j)**2,ustarmin)) 414 qsalt(j)=max(qsalt(j),0.) 415 enddo 416 417 418 else 419 ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001) 420 do j=1, knon 421 esalt=1./(3.25*max(ustarmin,ustar(j))) 422 hsalt(j)=0.08436*(max(ustarmin,ustar(j))**1.27) 423 qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt 424 !ep=qsalt*cdragm(j)*sqrt(u1(j)**2+v1(j)**2) 425 enddo 426 endif 427 428 ! calculation of erosion (flux positive towards the surface here) 429 ! consistent with implicit resolution of turbulent mixing equation 430 do j=1, knon 436 437 ! for outputs 438 do j=1, knon 431 439 i = knindex(j) 432 rhod=p1lay(j)/RD/temp_air(j)433 ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eq_salt of about 10s434 ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)435 ! integrated over tau_eq_salt to exceed the total mass of snow particle in the saltation layer436 ! (rho*qsalt*hsalt)437 maxerosion=hsalt(j)*qsalt(j)*rhod/tau_eq_salt438 fluxbs(j)=rhod*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) &439 / (1.-rhod*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime)440 !fluxbs(j)= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j))441 fluxbs(j)=min(0.,max(-maxerosion,fluxbs(j)))442 ! for outputs443 444 440 zxustartlic(i) = ustart(j) 445 441 zxrhoslic(i) = rhos(j) 446 enddo 447 448 endif 442 zxqsaltlic(i)=qsalt(j) 443 enddo 444 445 446 else 447 ! those lines are useful to calculate the snow age 448 CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)) 449 450 endif ! if ok_bs 449 451 450 452 -
LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90
r4830 r4835 363 363 REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: tpot, tpote, ue, uq, uwat, ve, vq, vwat, zxffonte 364 364 !$OMP THREADPRIVATE(tpot, tpote, ue, uq, uwat, ve, vq, vwat, zxffonte) 365 REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zxustartlic, zxrhoslic 366 !$OMP THREADPRIVATE(zxustartlic, zxrhoslic )365 REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zxustartlic, zxrhoslic, zxqsaltlic 366 !$OMP THREADPRIVATE(zxustartlic, zxrhoslic, zxqsaltlic) 367 367 REAL,ALLOCATABLE,SAVE,DIMENSION(:) :: zxfqcalving 368 368 !$OMP THREADPRIVATE(zxfqcalving) … … 958 958 ALLOCATE(zxtsol(klon), snow_lsc(klon), zxfqfonte(klon), zxqsurf(klon)) 959 959 ALLOCATE(zxrunofflic(klon)) 960 ALLOCATE(zxustartlic(klon), zxrhoslic(klon), zxqsaltlic(klon)) 961 zxustartlic(:)=0. ; zxrhoslic(:)=0. ; zxqsaltlic(:)=0. 960 962 ALLOCATE(rain_lsc(klon)) 961 963 ALLOCATE(rain_num(klon)) -
LMDZ6/trunk/libf/phylmdiso/phys_output_ctrlout_mod.F90
r4819 r4835 388 388 TYPE(ctrl_out), SAVE :: o_rhosnow_lic = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), & 389 389 'rhosnow_lic', 'snow density lic', 'kg/m3', (/ ('', i=1, 10) /)) 390 TYPE(ctrl_out), SAVE :: o_qsalt_lic = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), & 391 'qsalt_lic', 'qb in saltation layer lic', 'kg/kg', (/ ('', i=1, 10) /)) 390 392 TYPE(ctrl_out), SAVE :: o_sens_prec_liq_oce = ctrl_out((/ 5, 5, 10, 10, 5, 10, 11, 11, 11, 11/), & 391 393 'sens_rain_oce', 'Sensible heat flux of liquid prec. over ocean', 'W/m2', (/ ('', i=1, 10) /)) -
LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90
r4830 r4835 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 1953 CALL blowing_snow_ini(RCPD, RLSTT, RLVTT, RLMLT, & 1954 RVTMP2, RTT,RD,RG )1954 RVTMP2, RTT,RD,RG, RV, RPI) 1955 1955 ! Test de coherence sur oc_cdnc utilisé uniquement par cloud_optics_prop 1956 1956 IF (ok_newmicro) then -
LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90
r4724 r4835 35 35 USE cpl_mod, ONLY : cpl_send_landice_fields 36 36 USE calcul_fluxs_mod 37 USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic 37 USE phys_local_var_mod, ONLY : zxrhoslic, zxustartlic, zxqsaltlic 38 38 USE phys_output_var_mod, ONLY : snow_o,zfra_o 39 39 #ifdef ISO … … 50 50 !FC 51 51 USE ioipsl_getin_p_mod, ONLY : getin_p 52 USE lmdz_blowing_snow_ini, ONLY : zeta_bs, pbst_bs, prt_bs, iflag_saltation_bs53 52 USE lmdz_blowing_snow_ini, ONLY : c_esalt_bs, zeta_bs, pbst_bs, prt_bs, rhoice_bs, rhohard_bs 53 USE lmdz_blowing_snow_ini, ONLY : rhofresh_bs, tau_eqsalt_bs, tau_dens0_bs, tau_densmin_bs 54 54 #ifdef CPP_INLANDSIS 55 55 USE surf_inlandsis_mod, ONLY : surf_inlandsis … … 57 57 58 58 USE indice_sol_mod 59 60 59 61 60 ! INCLUDE "indicesol.h" … … 185 184 REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow 186 185 REAL, DIMENSION (klon,6) :: alb6 187 REAL :: rho0, rhoice, ustart0,esalt, rhod186 REAL :: esalt 188 187 REAL :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc 189 REAL :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard 190 REAL, DIMENSION(klon) :: ws1, rhos, ustart, qsalt, hsalt, snowini 191 REAL :: maxerosion ! in kg/s 188 REAL :: tau_dens, maxerosion, fluxbs_2 189 REAL, DIMENSION(klon) :: ws1, rhod, rhos, ustart0, ustart, qsalt, hsalt 192 190 193 191 ! End definition … … 420 418 ! 421 419 !**************************************************************************************** 422 423 CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)) 424 425 426 ! EV: following lines are obsolete since we set alb1 and alb2 to constant values 427 ! I therefore comment them 428 ! alb1(1:knon) = alb_neig(1:knon)*zfra(1:knon) + & 429 ! 0.6 * (1.0-zfra(1:knon)) 420 430 421 ! 431 422 !IM: plusieurs choix/tests sur l'albedo des "glaciers continentaux" … … 450 441 451 442 452 453 ! Simple blowing snow param 454 if (ok_bs) then 455 ustart0 = 0.211 456 rhoice = 920.0 457 rho0 = 300.0 458 rhomax=450.0 459 rhohard=450.0 460 tau_dens0=86400.0*10. ! 10 days by default, in s 461 tau_densmin=21600.0 ! 1/4 days according to in situ obs by C. Amory during blowing snow + 462 ! Marshall et al. 1999 (snow densification during rain) 463 443 !**************************************************************************************** 444 ! Simple blowing snow param 445 !**************************************************************************************** 446 ! we proceed in 2 steps: 447 ! first we erode - if possible -the accumulated snow during the time step 448 ! then we update the density of the underlying layer and see if we can also erode 449 ! this layer 450 451 452 if (ok_bs) then 453 fluxbs(:)=0. 454 do j=1,klon 455 ws1(j)=(u1(j)**2+v1(j)**2)**0.5 456 ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5 457 rhod(j)=p1lay(j)/RD/temp_air(j) 458 ustart0(j) =(log(2.868)-log(1.625))/0.085*sqrt(cdragm(j)) 459 enddo 460 461 ! 1st step: erosion of fresh snow accumulated during the time step 462 do j=1, knon 463 464 rhos(j)=rhofresh_bs 465 ! blowing snow flux formula used in MAR 466 ustart(j)=ustart0(j)*exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))*exp(max(0.,rhos(j)-rhohard_bs)) 467 ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs 468 ! computation of qbs at the top of the saltation layer 469 ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001) 470 esalt=1./(c_esalt_bs*max(1.e-6,ustar(j))) 471 hsalt(j)=0.08436*(max(1.e-6,ustar(j))**1.27) 472 qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt 473 ! calculation of erosion (flux positive towards the surface here) 474 ! consistent with implicit resolution of turbulent mixing equation 475 ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s 476 ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level) 477 ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer 478 ! (rho*qsalt*hsalt) 479 ! during this first step we also lower bound the erosion to the amount of fresh snow accumulated during the time step 480 maxerosion=min(precip_snow(j),hsalt(j)*qsalt(j)*rhod(j)/tau_eqsalt_bs) 481 482 fluxbs(j)=rhod(j)*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) & 483 / (1.-rhod(j)*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime) 484 !fluxbs(j)= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j)) 485 fluxbs(j)=max(-maxerosion,fluxbs(j)) 486 enddo 487 488 489 ! we now compute the snow age of the overlying layer (snow surface after erosion of the fresh snow accumulated during the time step) 490 ! this is done through the routine albsno 491 CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)+fluxbs(:)) 492 493 ! 2nd step: 464 494 ! computation of threshold friction velocity 465 495 ! which depends on surface snow density … … 468 498 ! snow density increases with snow age and 469 499 ! increases even faster in case of sedimentation of blowing snow or rain 470 tau_dens=max(tau_densmin , tau_dens0*exp(-abs(precip_bs(j))/pbst_bs &471 -abs(precip_rain(j))/prt_bs) *exp(-max(tsurf(j)-272.0,0.)))472 rhos(j)=rho 0+(rhohard-rho0)*(1.-exp(-agesno(j)*86400.0/tau_dens))500 tau_dens=max(tau_densmin_bs, tau_dens0_bs*exp(-abs(precip_bs(j))/pbst_bs - & 501 abs(precip_rain(j))/prt_bs)*exp(-max(tsurf(j)-RTT,0.))) 502 rhos(j)=rhofresh_bs+(rhohard_bs-rhofresh_bs)*(1.-exp(-agesno(j)*86400.0/tau_dens)) 473 503 ! blowing snow flux formula used in MAR 474 ws1(j)=(u1(j)**2+v1(j)**2)**0.5 475 ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5 476 ustart(j)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(j),0.))*exp(max(0.,rhos(j)-rhomax)) 477 ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till 478 ! rhohard<450) 504 ustart(j)=ustart0(j)*exp(max(rhoice_bs/rhofresh_bs-rhoice_bs/rhos(j),0.))*exp(max(0.,rhos(j)-rhohard_bs)) 505 ! we have multiplied by exp to prevent erosion when rhos>rhohard_bs 506 ! computation of qbs at the top of the saltation layer 507 ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001) 508 esalt=c_esalt_bs*ustar(j) 509 hsalt(j)=0.08436*(max(1.e-6,ustar(j))**1.27) 510 qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt 511 ! calculation of erosion (flux positive towards the surface here) 512 ! consistent with implicit resolution of turbulent mixing equation 513 ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within a time tau_eqsalt_bs of about 10s 514 ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level) 515 ! integrated over tau_eqsalt_bs to exceed the total mass of snow particle in the saltation layer 516 ! (rho*qsalt*hsalt) 517 maxerosion=hsalt(j)*qsalt(j)*rhod(j)/tau_eqsalt_bs 518 fluxbs_2=rhod(j)*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) & 519 / (1.-rhod(j)*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime) 520 !fluxbs_2= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j)) 521 fluxbs_2=max(-maxerosion,fluxbs_2) 522 fluxbs(j)=fluxbs(j)+fluxbs_2 479 523 enddo 480 481 ! computation of qbs at the top of the saltation layer 482 ! two formulations possible 483 ! pay attention that qbs is a mixing ratio and has to be converted 484 ! to specific content 485 486 if (iflag_saltation_bs .eq. 1) then 487 ! expression from CRYOWRF (Sharma et al. 2022) 488 aa=2.6 489 bb=2.5 490 cc=2.0 491 lambdasalt=0.45 492 do j =1, knon 493 rhod=p1lay(j)/RD/temp_air(j) 494 nunu=max(ustar(j)/ustart(j),1.e-3) 495 fluxsalt=rhod/RG*(ustar(j)**3)*(1.-nunu**(-2)) * & 496 (aa+bb*nunu**(-2)+cc*nunu**(-1)) 497 csalt=fluxsalt/(2.8*ustart(j)) 498 hsalt(j)=0.08436*ustar(j)**1.27 499 qsalt(j)=1./rhod*csalt*lambdasalt*RG/(max(ustar(j)**2,1E-6)) & 500 * exp(-lambdasalt*RG*hsalt(j)/max(ustar(j)**2,1E-6)) 501 qsalt(j)=max(qsalt(j),0.) 502 enddo 503 504 505 else 506 ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001) 507 do j=1, knon 508 esalt=1./(3.25*max(ustar(j),0.001)) 509 hsalt(j)=0.08436*ustar(j)**1.27 510 qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt 511 !ep=qsalt*cdragm(j)*sqrt(u1(j)**2+v1(j)**2) 512 enddo 513 endif 514 515 ! calculation of erosion (flux positive towards the surface here) 516 ! consistent with implicit resolution of turbulent mixing equation 517 do j=1, knon 524 525 ! for outputs 526 do j=1, knon 518 527 i = knindex(j) 519 rhod=p1lay(j)/RD/temp_air(j)520 ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within an interval of about 10s521 ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level)522 ! to exceed (in intensity) integrated over 10s to exceed the total mass of snow particle in the saltation layer523 ! (rho*qsalt*hsalt)524 maxerosion=hsalt(j)*qsalt(j)*rhod/10.525 fluxbs(j)=rhod*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) &526 / (1.-rhod*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime)527 !fluxbs(j)= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j))528 fluxbs(j)=max(-maxerosion,fluxbs(j))529 530 ! for outputs531 532 528 zxustartlic(i) = ustart(j) 533 529 zxrhoslic(i) = rhos(j) 534 enddo 535 536 endif 537 538 539 540 !**************************************************************************************** 541 ! Calculate surface snow amount 530 zxqsaltlic(i)=qsalt(j) 531 enddo 532 533 534 else 535 ! those lines are useful to calculate the snow age 536 CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)) 537 538 endif ! if ok_bs 539 540 541 542 543 544 545 !**************************************************************************************** 546 ! Calculate snow amount 542 547 ! 543 548 !****************************************************************************************
Note: See TracChangeset
for help on using the changeset viewer.
