MODULE dustdevil_mod IMPLICIT NONE INTEGER,SAVE :: dust_devil_scheme ! dust devil scheme number ! 0 (default): old Renno et al. 1998 (JAS 55, 3244-3252) scheme !$OMP THREADPRIVATE(dust_devil_scheme) CONTAINS SUBROUTINE dustdevil(ngrid,nlay,nq, pplev,pu,pv,pt, ptsurf,pq2, & pdqdev,pdqs_dev) use tracer_mod, only: alpha_devil use surfdat_h, only: z0_default USE comcstfi_h, ONLY: g, cpp, r, rcp USE mod_phys_lmdz_para, ONLY: is_master, bcast USE ioipsl_getin_p_mod, ONLY : getin_p IMPLICIT NONE !======================================================================= ! ! Parameterization of dust devil activities ! to estimate dust lifting ! The dust devil activity is estimated based on ! Renno et al. 1998 (JAS 55, 3244-3252) ! ! It is proportional to (1-b)*Fs ! ! With b= [ps**(rcp+1) - ptop**(rcp+1)] / [(ps-ptop)*(rcp+1)* ps**rcp] ! with ptop pressure of the top of the boundary layer ! rcp = R/cp ! ! and Fs the surface sensible heat flux = Cd*|U|*(T(1) -Tsurf) ! !======================================================================= ! arguments: ! ---------- INTEGER,INTENT(IN) :: ngrid,nlay REAL,INTENT(IN) :: pplev(ngrid,nlay+1) REAL,INTENT(IN) :: pt(ngrid,nlay) REAL,INTENT(IN) :: pu(ngrid,nlay) REAL,INTENT(IN) :: pv(ngrid,nlay) REAL,INTENT(IN) :: pq2(ngrid,nlay+1) REAL,INTENT(IN) :: ptsurf(ngrid) ! Traceurs : INTEGER,INTENT(IN) :: nq real,intent(out) :: pdqdev(ngrid,nlay,nq) real,intent(out) :: pdqs_dev(ngrid,nq) ! local: ! ------ INTEGER :: ig,l,iq real,save :: Cd !$OMP THREADPRIVATE(Cd) real :: z1 LOGICAL,SAVE :: firstcall=.true. !$OMP THREADPRIVATE(firstcall) character(len=20),parameter :: rname="dustdevil" REAL :: devila(ngrid) integer :: ltop(ngrid) real :: b,rho,Fs,wind REAL,SAVE :: q2top=.5 ! value of q2 at the top of PBL REAL,SAVE :: seuil=.3 ! value of minimum dust devil activity ! for dust lifting !$OMP THREADPRIVATE(q2top,seuil) !----------------------------------------------------------------------- ! initialisation ! -------------- ! AS: OK firstcall absolute IF (firstcall) THEN ! get scheme number (default=0, old scheme) dust_devil_scheme=0 call getin_p("dust_devil_scheme",dust_devil_scheme) ! sanity check on available dust_devil_scheme values if (dust_devil_scheme.ne.0) then write(*,*) trim(rname)//" wrong value for dust_devil_scheme:",& dust_devil_scheme call abort_physic(rname,"bad dust_devil_scheme value",1) endif if(is_master) then write(*,*) 'In dustdevil :' write(*,*) ' q2top= ',q2top,' seuil= ', seuil ! A rough estimation of the horizontal drag coefficient Cd ! (the scale heigh is taken to be 13 km since we are typically ! dealing with daytime temperature around 250K. ! z1= -0.5*13.e3*log(pplev(1,2)/pplev(1,1)) Cd = (0.4/log(z1/z0_default))**2 ! Temporaire ! open(77,file='devil') endif !is_master ! share the info with all cores CALL bcast(z1) CALL bcast(Cd) firstcall=.false. ENDIF ! firstcall !----------------------------------------------------------------------- ! Initialisation of outputs do iq=1,nq do l=1,nlay do ig=1,ngrid pdqdev(ig,l,iq)= 0 end do end do end do pdqs_dev(1:ngrid,1:nq)=0 if (dust_devil_scheme==0) then !----------------------------------------------------------------------- ! Determining the top of the boundary layer ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ do ig=1,ngrid do l=2,nlay-1 if (pq2(ig,l).lt.q2top)then ltop(ig)=l goto 99 end if enddo 99 continue ! *************************************** !c if (ptsurf(ig).gt.255)then ! write(*,*) 'tsurf, ztop (km): ', ig, ! & ptsurf(ig), -12.*log(pplev(ig,ltop(ig))/pplev(ig,1)) ! if ((ptsurf(ig).gt.50.).and.( ! & (-12.*log(pplev(ig,ltop(ig))/pplev(ig,1))).gt.60.))then ! do l=1,nlay ! write(*,*) l, ! & -12.*log(pplev(ig,l)/pplev(ig,1)),pq2(ig,l) ! end do ! end if !c end if ! *************************************** enddo ! *************************************** ! do ig=100,148 ! write(*,*)'time,tsurf,ztop', localtime(ig),ptsurf(ig), ! & -12.*log(pplev(ig,ltop(ig))/pplev(ig,1)) ! end do ! *************************************** ! Calculation : dust devil intensity ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ do ig=1,ngrid ! -------------------------------------------------- ! 1) Version 1 : sensible heat flux using actual wind : ! Wind magnitude: ! wind = sqrt(pu(ig,1)**2+pv(ig,1)**2) ! ! -------------------------------------------------- ! 2) Version 2 : sensible heat flux using wind = 15 m/s wind = 15. ! ---------------------------------------------------- ! Density : rho=pplev(ig,1)/(R*pt(ig,1)) ! Sensible heat flux (W.m-2) (>0 if up) Fs= rho*cpp*Cd * wind * (ptsurf(ig) -pt(ig,1)) b= (pplev(ig,1)**(rcp+1) - pplev(ig,ltop(ig))**(rcp+1)) / & ( (pplev(ig,1)-pplev(ig,ltop(ig)))*(rcp+1)*pplev(ig,1)**rcp) ! b_diag(ig) = b ! TEMPORAIRE (pour diagnostique) ! Energy flux available to drive dust devil (W.m-2) : (1-b)*Fs ! Dust devil activity : devila(ig)= max( 0. , (1-b)*Fs - seuil ) enddo ! ! lifted dust (kg m-2 s-1) (<0 when lifting) ! ~~~~~~~~~~ do iq=1,nq do ig=1,ngrid pdqs_dev(ig,iq)= - alpha_devil(iq) * devila(ig) end do end do ! Injection of dust in the atmosphere (up to the top of pbl) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ do iq=1,nq do ig=1,ngrid if (devila(ig).ne.0.) then do l=1,ltop(ig) pdqdev(ig,l,iq)=-pdqs_dev(ig,iq)*g/ & (pplev(ig,1)-pplev(ig,ltop(ig))) end do end if end do end do endif ! of if (dust_devil_scheme==0) END SUBROUTINE dustdevil END MODULE dustdevil_mod