Changeset 4672
- Timestamp:
- Sep 4, 2023, 4:35:35 PM (9 months ago)
- Location:
- LMDZ6/trunk/libf
- Files:
-
- 4 edited
-
phylmd/blowing_snow_ini_mod.F90 (modified) (4 diffs)
-
phylmd/blowing_snow_sublim_sedim.F90 (modified) (4 diffs)
-
phylmd/surf_landice_mod.F90 (modified) (6 diffs)
-
phylmdiso/surf_landice_mod.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/blowing_snow_ini_mod.F90
r4529 r4672 8 8 real coef_eva_bs, expo_eva_bs, fallv_bs, zeta_bs 9 9 real prt_bs, pbst_bs, qbst_bs 10 integer :: niter_bs,iflag_saltation_bs10 integer :: iflag_saltation_bs 11 11 12 12 !$OMP THREADPRIVATE(prt_level, lunout) … … 14 14 !$OMP THREADPRIVATE(coef_eva_bs, expo_eva_bs, fallv_bs, zeta_bs) 15 15 !$OMP THREADPRIVATE(pbst_bs, prt_bs, qbst_bs) 16 !$OMP THREADPRIVATE( niter_bs,iflag_saltation_bs)16 !$OMP THREADPRIVATE(iflag_saltation_bs) 17 17 18 18 contains … … 45 45 CALL getin_p('qbst_bs',qbst_bs) 46 46 47 pbst_bs= 0.0 147 pbst_bs= 0.0003 48 48 CALL getin_p('pbst_bs',pbst_bs) 49 49 … … 64 64 CALL getin_p('expo_eva_bs',expo_eva_bs) 65 65 66 niter_bs = 567 CALL getin_p('niter_bs',niter_bs)68 69 66 iflag_saltation_bs=0 70 67 CALL getin_p('iflag_saltation_bs',iflag_saltation_bs) -
LMDZ6/trunk/libf/phylmd/blowing_snow_sublim_sedim.F90
r4664 r4672 9 9 10 10 use blowing_snow_ini_mod, only : coef_eva_bs,RTT,RD,RG,expo_eva_bs, fallv_bs 11 use blowing_snow_ini_mod, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2 , niter_bs11 use blowing_snow_ini_mod, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2 12 12 USE lmdz_lscp_tools, only : calc_qsat_ecmwf 13 13 … … 125 125 ! vapor, temperature, precip fluxes update 126 126 zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime 127 zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime 127 zq(i) = max(0., zq(i)) 128 !zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime 129 !zqbs(i) = max(0., zqbs(i)) 128 130 ! melting 129 131 zt(i) = zt(i) & … … 151 153 152 154 153 ! vapor, temperature, precip fluxes update 155 ! vapor, temperature, precip fluxes update following sublimation 154 156 zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime 155 zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime 157 zq(i) = max(0., zq(i)) 158 !zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime 159 !zqbs(i) = max(0., zqbs(i)) 156 160 zt(i) = zt(i) & 157 161 + (sedimn(i)-sedim(i)) & … … 173 177 ENDDO ! loop on ngrid 174 178 175 ! Now sedimention scheme (alike that in lscp) 176 DO n = 1, niter_bs 177 DO i = 1, ngrid 178 zrho(i,k) = pplay(i,k) / zt(i) / RD 179 zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG) 180 velo(i,k) = fallv_bs 181 dqsedim = dtime/REAL(niter_bs)/zdz(i)*zqbs(i)*velo(i,k) ! dqice/dt=1/rho*d(rho*wice*qice)/dz 182 precbs = MIN(MAX(dqsedim,0.0),zqbs(i)) 183 zqbs(i) = MAX(zqbs(i)-1.*precbs , 0.0) 184 ENDDO !loop on ngrid 185 ENDDO ! loop on niter_bs 186 187 ! add to non sublimated precip 188 DO i=1,ngrid 189 sedim(i) = sedim(i)+max(0.,zqbsi(i)-zqbs(i))*(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 179 ! Now sedimention scheme 180 181 ! exact resolution of the conservation equation for qbs with the updated flux from the top (after evap) 182 ! valid only if the fall velocity is constant 183 184 DO i = 1, ngrid 185 186 zrho(i,k) = pplay(i,k) / zt(i) / RD 187 zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG) 188 velo(i,k) = fallv_bs 189 zqbs(i) = zqbsi(i)*exp(-velo(i,k)/zdz(i)*dtime)+sedim(i)/zrho(i,k)/velo(i,k) 190 zqbs(i) = max(zqbs(i),0.) 191 ! flux update 192 sedim(i) = sedim(i) + zrho(i,k)*zdz(i)/dtime*(zqbsi(i)-zqbs(i)) 193 sedim(i) = max(0.,sedim(i)) 194 190 195 ENDDO 196 197 ! old version with bugs 198 ! DO n = 1, niter_bs 199 ! DO i = 1, ngrid 200 ! zrho(i,k) = pplay(i,k) / zt(i) / RD 201 ! zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG) 202 ! velo(i,k) = fallv_bs 203 ! dqsedim = dtime/REAL(niter_bs)/zdz(i)*zqbs(i)*velo(i,k) ! dqice/dt=1/rho*d(rho*wice*qice)/dz 204 ! precbs = MIN(MAX(dqsedim,0.0),zqbs(i)) 205 ! zqbs(i) = MAX(zqbs(i)-1.*precbs , 0.0) 206 ! ENDDO !loop on ngrid 207 ! ENDDO ! loop on niter_bs 208 ! ! add to non sublimated precip 209 ! DO i=1,ngrid 210 ! sedim(i) = sedim(i)+max(0.,zqbsi(i)-zqbs(i))*(paprs(i,k)-paprs(i,k+1))/(RG*dtime) 211 ! ENDDO 212 ! 213 214 191 215 192 216 ! Outputs: -
LMDZ6/trunk/libf/phylmd/surf_landice_mod.F90
r4538 r4672 140 140 REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow 141 141 REAL, DIMENSION (klon,6) :: alb6 142 REAL :: rho0, rhoice, ustart0, hsalt,esalt, rhod142 REAL :: rho0, rhoice, ustart0,esalt, rhod 143 143 REAL :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc 144 144 REAL :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard 145 REAL, DIMENSION(klon) :: ws1, rhos, ustart, qsalt 145 REAL, DIMENSION(klon) :: ws1, rhos, ustart, qsalt, hsalt, snowini 146 REAL :: maxerosion ! in kg/s 147 146 148 ! End definition 147 149 !**************************************************************************************** … … 364 366 ustart0 = 0.211 365 367 rhoice = 920.0 366 rho0 = 200.0368 rho0 = 300.0 367 369 rhomax=450.0 368 rhohard=4 00.0370 rhohard=450.0 369 371 tau_dens0=86400.0*10. ! 10 days by default, in s 370 tau_densmin=86400.0 ! 1 days according to in situ obs by C. Amory 372 tau_densmin=21600.0 ! 1/4 days according to in situ obs by C. Amory during blowing snow + 373 ! Marshall et al. 1999 (snow densification during rain) 371 374 372 375 ! computation of threshold friction velocity 373 376 ! which depends on surface snow density 374 do i= 1, knon377 do j = 1, knon 375 378 ! estimation of snow density 376 379 ! snow density increases with snow age and 377 ! increases even faster in case of sedimentation of blowing snow 378 tau_dens=max(tau_densmin, tau_dens0*exp(-abs(precip_bs(i))/pbst_bs-abs(precip_rain(i))/prt_bs)) 379 rhos(i)=rho0+(rhohard-rho0)*(1.-exp(-agesno(i)*86400.0/tau_dens)) 380 ! increases even faster in case of sedimentation of blowing snow or rain 381 tau_dens=max(tau_densmin, & tau_dens0*exp(-abs(precip_bs(j))/pbst_bs-abs(precip_rain(j))/prt_bs) & 382 *exp(-max(tsurf(j)-272.0,0.))) 383 rhos(j)=rho0+(rhohard-rho0)*(1.-exp(-agesno(j)*86400.0/tau_dens)) 380 384 ! blowing snow flux formula used in MAR 381 ws1( i)=(u1(i)**2+v1(i)**2)**0.5382 ustar( i)=(cdragm(i)*(u1(i)**2+v1(i)**2))**0.5383 ustart( i)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(i),0.))*exp(max(0.,rhos(i)-rhomax))385 ws1(j)=(u1(j)**2+v1(j)**2)**0.5 386 ustar(j)=(cdragm(j)*(u1(j)**2+v1(j)**2))**0.5 387 ustart(j)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(j),0.))*exp(max(0.,rhos(j)-rhomax)) 384 388 ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till 385 389 ! rhohard<450) … … 397 401 cc=2.0 398 402 lambdasalt=0.45 399 do i=1, knon400 rhod=p1lay( i)/RD/temp_air(i)401 nunu=max(ustar( i)/ustart(i),1.e-3)402 fluxsalt=rhod/RG*(ustar( i)**3)*(1.-nunu**(-2)) * &403 do j =1, knon 404 rhod=p1lay(j)/RD/temp_air(j) 405 nunu=max(ustar(j)/ustart(j),1.e-3) 406 fluxsalt=rhod/RG*(ustar(j)**3)*(1.-nunu**(-2)) * & 403 407 (aa+bb*nunu**(-2)+cc*nunu**(-1)) 404 csalt=fluxsalt/(2.8*ustart( i))405 hsalt =0.08436*ustar(i)**1.27406 qsalt( i)=1./rhod*csalt*lambdasalt*RG/(max(ustar(i)**2,1E-6)) &407 * exp(-lambdasalt*RG*hsalt /max(ustar(i)**2,1E-6))408 qsalt( i)=max(qsalt(i),0.)408 csalt=fluxsalt/(2.8*ustart(j)) 409 hsalt(j)=0.08436*ustar(j)**1.27 410 qsalt(j)=1./rhod*csalt*lambdasalt*RG/(max(ustar(j)**2,1E-6)) & 411 * exp(-lambdasalt*RG*hsalt(j)/max(ustar(j)**2,1E-6)) 412 qsalt(j)=max(qsalt(j),0.) 409 413 enddo 410 414 … … 412 416 else 413 417 ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001) 414 do i=1, knon415 esalt=1./(3.25*max(ustar( i),0.001))416 hsalt =0.08436*ustar(i)**1.27417 qsalt( i)=(max(ustar(i)**2-ustart(i)**2,0.))/(RG*hsalt)*esalt418 !ep=qsalt*cdragm( i)*sqrt(u1(i)**2+v1(i)**2)418 do j=1, knon 419 esalt=1./(3.25*max(ustar(j),0.001)) 420 hsalt(j)=0.08436*ustar(j)**1.27 421 qsalt(j)=(max(ustar(j)**2-ustart(j)**2,0.))/(RG*hsalt(j))*esalt 422 !ep=qsalt*cdragm(j)*sqrt(u1(j)**2+v1(j)**2) 419 423 enddo 420 424 endif 421 425 422 ! calculation of erosion ( emission flux towards the first atmospheric level)426 ! calculation of erosion (flux positive towards the surface here) 423 427 ! consistent with implicit resolution of turbulent mixing equation 424 do i=1, knon 425 rhod=p1lay(i)/RD/temp_air(i) 426 fluxbs(i)=rhod*ws1(i)*cdragm(i)*zeta_bs*(AcoefQBS(i)-qsalt(i)) & 427 / (1.-rhod*ws1(i)*zeta_bs*cdragm(i)*BcoefQBS(i)*dtime) 428 !fluxbs(i)= zeta_bs*rhod*ws1(i)*cdragm(i)*(qbs1(i)-qsalt(i)) 428 do j=1, knon 429 i = knindex(j) 430 rhod=p1lay(j)/RD/temp_air(j) 431 ! Nemoto and Nishimura 2004 show that steady-state saltation is achieved within an interval of about 10s 432 ! we thus prevent snowerosion (snow particle transfer from the saltation layer to the first model level) 433 ! to exceed (in intensity) integrated over 10s to exceed the total mass of snow particle in the saltation layer 434 ! (rho*qsalt*hsalt) 435 maxerosion=hsalt(j)*qsalt(j)*rhod/10. 436 fluxbs(j)=rhod*ws1(j)*cdragm(j)*(AcoefQBS(j)-qsalt(j)) & 437 / (1.-rhod*ws1(j)*cdragm(j)*BcoefQBS(j)*dtime) 438 !fluxbs(j)= rhod*ws1(j)*cdragm(j)*(qbs1(j)-qsalt(j)) 439 fluxbs(j)=max(-maxerosion,fluxbs(j)) 440 441 ! for outputs 442 443 zxustartlic(i) = ustart(j) 444 zxrhoslic(i) = rhos(j) 429 445 enddo 430 446 431 ! for outputs432 do j = 1, knon433 i = knindex(j)434 zxustartlic(i) = ustart(j)435 zxrhoslic(i) = rhos(j)436 enddo437 438 447 endif 439 448 … … 441 450 442 451 !**************************************************************************************** 443 ! Calculate s urface snow amount452 ! Calculate snow amount 444 453 ! 445 454 !**************************************************************************************** … … 451 460 evap_totsnow(:)=evap(:) 452 461 ENDIF 453 462 463 454 464 CALL fonte_neige(knon, is_lic, knindex, dtime, & 455 465 tsurf, precip_rain, precip_totsnow, & 456 466 snow, qsol, tsurf_new, evap_totsnow) 457 467 468 458 469 WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0. 459 470 zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/(snow(1:knon)+10.0))) -
LMDZ6/trunk/libf/phylmdiso/surf_landice_mod.F90
r4538 r4672 185 185 REAL,DIMENSION(klon) :: precip_totsnow, evap_totsnow 186 186 REAL, DIMENSION (klon,6) :: alb6 187 REAL :: rho0, rhoice, ustart0, hsalt,esalt, rhod187 REAL :: rho0, rhoice, ustart0,esalt, rhod 188 188 REAL :: lambdasalt,fluxsalt, csalt, nunu, aa, bb, cc 189 189 REAL :: tau_dens, tau_dens0, tau_densmin, rhomax, rhohard 190 REAL, DIMENSION(klon) :: ws1, rhos, ustart, qsalt 190 REAL, DIMENSION(klon) :: ws1, rhos, ustart, qsalt, hsalt, snowini 191 REAL :: maxerosion ! in kg/s 192 191 193 ! End definition 192 194 !**************************************************************************************** … … 453 455 ustart0 = 0.211 454 456 rhoice = 920.0 455 rho0 = 200.0457 rho0 = 300.0 456 458 rhomax=450.0 457 rhohard=4 00.0459 rhohard=450.0 458 460 tau_dens0=86400.0*10. ! 10 days by default, in s 459 tau_densmin=86400.0 ! 1 days according to in situ obs by C. Amory 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) 460 463 461 464 ! computation of threshold friction velocity 462 465 ! which depends on surface snow density 463 do i= 1, knon466 do j = 1, knon 464 467 ! estimation of snow density 465 468 ! snow density increases with snow age and 466 ! increases even faster in case of sedimentation of blowing snow 467 tau_dens=max(tau_densmin, tau_dens0*exp(-abs(precip_bs(i))/pbst_bs-abs(precip_rain(i))/prt_bs)) 468 rhos(i)=rho0+(rhohard-rho0)*(1.-exp(-agesno(i)*86400.0/tau_dens)) 469 ! 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-abs(precip_rain(j))/prt_bs) & 471 *exp(-max(tsurf(j)-272.0,0.))) 472 rhos(j)=rho0+(rhohard-rho0)*(1.-exp(-agesno(j)*86400.0/tau_dens)) 469 473 ! blowing snow flux formula used in MAR 470 ws1( i)=(u1(i)**2+v1(i)**2)**0.5471 ustar( i)=(cdragm(i)*(u1(i)**2+v1(i)**2))**0.5472 ustart( i)=ustart0*exp(max(rhoice/rho0-rhoice/rhos(i),0.))*exp(max(0.,rhos(i)-rhomax))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)) 473 477 ! we have multiplied by exp to prevent erosion when rhos>rhomax (usefull till 474 478 ! rhohard<450) … … 486 490 cc=2.0 487 491 lambdasalt=0.45 488 do i=1, knon489 rhod=p1lay( i)/RD/temp_air(i)490 nunu=max(ustar( i)/ustart(i),1.e-3)491 fluxsalt=rhod/RG*(ustar( i)**3)*(1.-nunu**(-2)) * &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)) * & 492 496 (aa+bb*nunu**(-2)+cc*nunu**(-1)) 493 csalt=fluxsalt/(2.8*ustart( i))494 hsalt =0.08436*ustar(i)**1.27495 qsalt( i)=1./rhod*csalt*lambdasalt*RG/(max(ustar(i)**2,1E-6)) &496 * exp(-lambdasalt*RG*hsalt /max(ustar(i)**2,1E-6))497 qsalt( i)=max(qsalt(i),0.)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.) 498 502 enddo 499 503 … … 501 505 else 502 506 ! default formulation from MAR model (Amory et al. 2021, Gallee et al. 2001) 503 do i=1, knon504 esalt=1./(3.25*max(ustar( i),0.001))505 hsalt =0.08436*ustar(i)**1.27506 qsalt( i)=(max(ustar(i)**2-ustart(i)**2,0.))/(RG*hsalt)*esalt507 !ep=qsalt*cdragm( i)*sqrt(u1(i)**2+v1(i)**2)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) 508 512 enddo 509 513 endif 510 514 511 ! calculation of erosion ( emission flux towards the first atmospheric level)515 ! calculation of erosion (flux positive towards the surface here) 512 516 ! consistent with implicit resolution of turbulent mixing equation 513 do i=1, knon 514 rhod=p1lay(i)/RD/temp_air(i) 515 fluxbs(i)=rhod*ws1(i)*cdragm(i)*zeta_bs*(AcoefQBS(i)-qsalt(i)) & 516 / (1.-rhod*ws1(i)*zeta_bs*cdragm(i)*BcoefQBS(i)*dtime) 517 !fluxbs(i)= zeta_bs*rhod*ws1(i)*cdragm(i)*(qbs1(i)-qsalt(i)) 518 enddo 519 520 ! for outputs 521 do j = 1, knon 522 i = knindex(j) 523 zxustartlic(i) = ustart(j) 524 zxrhoslic(i) = rhos(j) 517 do j=1, knon 518 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 10s 521 ! 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 layer 523 ! (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 outputs 531 532 zxustartlic(i) = ustart(j) 533 zxrhoslic(i) = rhos(j) 525 534 enddo 526 535
Note: See TracChangeset
for help on using the changeset viewer.
