MODULE rocketduststorm_mod IMPLICIT NONE REAL, SAVE, ALLOCATABLE :: dustliftday(:) ! dust lifting rate (s-1) REAL, SAVE, ALLOCATABLE :: alpha_hmons(:) ! slope winds lifting mesh fraction CONTAINS !======================================================================= ! ROCKET DUST STORM - vertical transport and detrainment !======================================================================= ! calculating the vertical flux ! calling vl_storm : transport scheme of stormdust tracers ! detrainement of stormdust into nomal background dust ! ----------------------------------------------------------------------- ! Authors: C. Wang; F. Forget; T. Bertrand ! Institution: Laboratoire de Meteorologie Dynamique (LMD) Paris, France ! ----------------------------------------------------------------------- SUBROUTINE rocketduststorm(ngrid,nlayer,nq,ptime,ptimestep, & pq,pdqfi,pt,pdtfi,pplev,pplay,pzlev, & pzlay,pdtsw,pdtlw, & ! input for radiative transfer clearatm,icount,zday,zls, & tsurf,igout,totstormfract, & ! input sub-grid scale cloud clearsky,totcloudfrac, & ! output pdqrds,wspeed,dsodust,dsords) use tracer_mod, only: igcm_stormdust_mass,igcm_stormdust_number & ,igcm_dust_mass,igcm_dust_number & ,rho_dust USE comcstfi_h, only: r,g,cpp,rcp use dimradmars_mod, only: albedo,naerkind use comsaison_h, only: dist_sol,mu0,fract use surfdat_h, only: emis,co2ice,zmea, zstd, zsig, hmons use planetwide_mod, only: planetwide_maxval,planetwide_minval ! use rocketstorm_h, only: rdsinject use callradite_mod implicit none !-------------------------------------------------------- ! Input Variables !-------------------------------------------------------- INTEGER, INTENT(IN) :: ngrid ! number of horizontal grid points INTEGER, INTENT(IN) :: nlayer ! number of vertical grid points INTEGER, INTENT(IN) :: nq ! number of tracer species REAL, INTENT(IN) :: ptime REAL, INTENT(IN) :: ptimestep REAL, INTENT(IN) :: pq(ngrid,nlayer,nq) ! advected field nq REAL, INTENT(IN) :: pdqfi(ngrid,nlayer,nq)! tendancy field pq REAL, INTENT(IN) :: pt(ngrid,nlayer) ! temperature at mid-layer (K) REAL, INTENT(IN) :: pdtfi(ngrid,nlayer) ! tendancy temperature at mid-layer (K) REAL, INTENT(IN) :: pplay(ngrid,nlayer) ! pressure at middle of the layers REAL, INTENT(IN) :: pplev(ngrid,nlayer+1) ! pressure at intermediate levels REAL, INTENT(IN) :: pzlay(ngrid,nlayer) ! altitude at the middle of the layers REAL, INTENT(IN) :: pzlev(ngrid,nlayer+1) ! altitude at layer boundaries REAL, INTENT(IN) :: pdtsw(ngrid,nlayer) ! (K/s) env REAL, INTENT(IN) :: pdtlw(ngrid,nlayer) ! (K/s) env ! input for second radiative transfer logical, INTENT(IN) :: clearatm INTEGER, INTENT(INOUT) :: icount real, intent(in) :: zday real, intent(in) :: zls real, intent(in) :: tsurf(ngrid) integer, intent(in) :: igout real, intent(in) :: totstormfract(ngrid) ! sbgrid scale water ice clouds logical, intent(in) :: clearsky real, intent(in) :: totcloudfrac(ngrid) !-------------------------------------------------------- ! Output Variables !-------------------------------------------------------- REAL, INTENT(OUT) :: pdqrds(ngrid,nlayer,nq) ! tendancy field for dust when detraining REAL, INTENT(OUT) :: wspeed(ngrid,nlayer+1) ! vertical speed within the rocket dust storm REAL, INTENT(OUT) :: dsodust(ngrid,nlayer) ! density scaled opacity of env. dust REAL, INTENT(OUT) :: dsords(ngrid,nlayer) ! density scaled opacity of storm dust !-------------------------------------------------------- ! Local variables !-------------------------------------------------------- INTEGER l,ig,tsub,iq,ll ! chao local variables from callradite.F REAL zdtlw1(ngrid,nlayer) ! (K/s) storm REAL zdtsw1(ngrid,nlayer) ! (K/s) storm REAL zt(ngrid,nlayer) ! actual temperature at mid-layer (K) REAL zdtvert(nlayer) ! dT/dz , lapse rate REAL ztlev(nlayer) ! temperature at intermediate levels l+1/2 without last level REAL zdtlw1_lev(nlayer),zdtsw1_lev(nlayer) ! rad. heating rate at intermediate levels l+1/2 REAL zdtlw_lev(nlayer),zdtsw_lev(nlayer) ! rad. heating rate at intermediate levels l+1/2 REAL zqstorm_mass(ngrid,nlayer) ! tracer pq mass intermediate REAL zqstorm_mass_col(nlayer) REAL zqstorm_number(ngrid,nlayer) ! tracer field pq number intermediate REAL zqstorm_number_col(nlayer) REAL zqi_mass(ngrid,nlayer) ! tracer pq mass intermediate REAL zqi_number(ngrid,nlayer) ! tracer pq mass intermediate REAL zdqvlstorm_mass(ngrid,nlayer) ! tendancy pdq mass after vertical transport REAL zdqvlstorm_number(ngrid,nlayer) ! tendancy pdq number after vertical transport REAL zdqdetstorm_mass(ngrid,nlayer) ! tendancy field pq mass after detrainment only REAL zdqdetstorm_number(ngrid,nlayer) ! tendancy field pq number after detrainment only REAL zdqenv_mass(ngrid,nlayer) ! tendancy pdq mass from dust-> ! stormdust in slp REAL zdqenv_number(ngrid,nlayer) ! tendancy pdq number from dust-> ! stormdust in slp REAL masse(nlayer) REAL zq(ngrid,nlayer,nq) ! updated tracers REAL wrds(nlayer) ! vertical flux within the rocket dust storm REAL wqrdsmass(nlayer+1) ! flux mass from vl_storm REAL wqrdsnumber(nlayer+1) ! flux number from vl_storm INTEGER nsubtimestep !number of subtimestep when calling vl_storm REAL subtimestep !ptimestep/nsubtimestep REAL dtmax !considered time needed for dust to cross one layer !minimal value over a column logical storm(ngrid) !logical : true if you have some storm dust in the column ! real slope(ngrid) !logical : true if you don't have storm and have !a slope ! real wslplev(ngrid,nlayer) ! real wslp(ngrid,nlayer) REAL coefdetrain !coefficient for detrainment : % of stormdust detrained REAL,PARAMETER:: coefmin =0.025 !C 0 mu0lim and with storm=false) !! case 3 rocket dust storm (storm=true) !! 3. Vertical transport !! 4. Detrainment ! ********************************************************************** ! ********************************************************************** !! ********************************************************************** !! Firstcall: Evaluate slope wind mesh fraction IF (firstcall) then call planetwide_maxval(hmons,hmax ) call planetwide_minval(hmons,hmin ) do ig=1,ngrid ! It's hard to know the exact the scale of upward flow, ! we assume that the maximun is 10% of the mesh area. alpha_hmons(ig)= 0.1*(hmons(ig)-hmin)/(hmax-hmin) enddo firstcall = .false. ENDIF !firstcall ! ********************************************************************** ! 0. Initializations ! ********************************************************************** storm(:)=.false. pdqrds(:,:,:) = 0. zdqdetstorm_mass(:,:)=0. zdqdetstorm_number(:,:)=0. wspeed(:,:)=0. detrainment(:,:)=0. zqstorm_mass_col(:)=0. zqstorm_number_col(:)=0. lapserate(:,:)=0. deltahr(:,:)=0. rdsdustqvl0(:,:)=0. rdsdustqvl1(:,:)=0. zqstormdustslp(:,:)=0. znstormdustslp(:,:)=0. zqdustslp(:,:)=0. zndustslp(:,:)=0. zq(:,:,:) = 0. w0(:)=0. ! w1(:)=0. ! ztb1(:)=0. newzt(:,:)=0 wtemp(:,:)=0. wadiabatic(:,:)=0. slpbg(:)=0. ! buoyt(:)=0. tnew(:)=0. envtemp(:)=0. envt(:,:)=0. scheme(:)=0 alpha(:)=0. stormdust_m0(:,:)=0. stormdust_m1(:,:)=0. stormdust_m2(:,:)=0. ! totdust0(:)=0. ! totdust1(:)=0. !! no update for the stormdust tracer and temperature fields !! because previous callradite was for background dust zq(1:ngrid,1:nlayer,1:nq)=pq(1:ngrid,1:nlayer,1:nq) zt(1:ngrid,1:nlayer)=pt(1:ngrid,1:nlayer) !! get actual q for stormdust and dust tracers do l=1,nlayer do ig=1, ngrid zqi_mass(ig,l)=zq(ig,l,igcm_dust_mass) zqi_number(ig,l)=zq(ig,l,igcm_dust_number) zqstorm_mass(ig,l)=zq(ig,l,igcm_stormdust_mass) zqstorm_number(ig,l)=zq(ig,l,igcm_stormdust_number) !for diagnostics: stormdust_m0(ig,l)=zqstorm_mass(ig,l) enddo enddo ! of do l=1,nlayer !! Check if there is a rocket dust storm do ig=1,ngrid storm(ig)=.false. do l=1,nlayer if (zqstorm_mass(ig,l)/zqi_mass(ig,l) .gt. 1.E-4) then storm(ig)=.true. exit endif enddo enddo ! ********************************************************************* ! 1. Call the second radiative transfer for stormdust, obtain the extra heating ! ********************************************************************* CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo, & emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, & zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1, & fluxtop_sw1,tauref,tau,aerosol,dsodust,tauscaling, & taucloudtes,rdust,rice,nuice,co2ice,rstormdust, & totstormfract,clearatm,dsords, & clearsky,totcloudfrac) ! ********************************************************************** ! 2. Compute vertical velocity for storm dust ! ********************************************************************** DO ig=1,ngrid !! ********************************************************************** !! 2.1 case 1: Nothing to do when no storm and no slope ! IF ((mu0(ig) .LE. mu0lim) .AND. .NOT.(storm(ig)) ) then ! scheme(ig)=1 ! cycle ! endif IF ((alpha_hmons(ig) .EQ. 0.) .AND. .NOT.(storm(ig))) then scheme(ig)=1 cycle !!no slope endif ! whatever the situation is, we need the vertical velocity computed by ! the rds scheme zdtvert(1)=0. !This is the env. lapse rate DO l=1,nlayer-1 zdtvert(l+1)=(zt(ig,l+1)-zt(ig,l))/(pzlay(ig,l+1)-pzlay(ig,l)) lapserate(ig,l+1)=zdtvert(l+1) !for diagnosing ENDDO ! computing heating rates gradient at boundraies of each layer ! start from surface zdtlw1_lev(1)=0. zdtsw1_lev(1)=0. zdtlw_lev(1)=0. zdtsw_lev(1)=0. ztlev(1)=zt(ig,1) DO l=1,nlayer-1 ! Calculation for the dust storm fraction zdtlw1_lev(l+1)=(zdtlw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ & zdtlw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / & (pzlay(ig,l+1)-pzlay(ig,l)) zdtsw1_lev(l+1)=(zdtsw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ & zdtsw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / & (pzlay(ig,l+1)-pzlay(ig,l)) !MV18: calculation for the background dust fraction zdtlw_lev(l+1)=(pdtlw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ & pdtlw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / & (pzlay(ig,l+1)-pzlay(ig,l)) zdtsw_lev(l+1)=(pdtsw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ & pdtsw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / & (pzlay(ig,l+1)-pzlay(ig,l)) ztlev(l+1)=(zt(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ & zt(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l))) / & (pzlay(ig,l+1)-pzlay(ig,l)) ENDDO DO l=1,nlayer deltahr(ig,l)=(zdtlw1_lev(l)+zdtsw1_lev(l)) & -(zdtlw_lev(l)+zdtsw_lev(l)) wadiabatic(ig,l)=-deltahr(ig,l)/(g/cpp+ & max(zdtvert(l),-0.99*g/cpp)) !limit vertical wind in case of lapse rate close to adiabat wadiabatic(ig,l)=max(wadiabatic(ig,l),-wlim) wadiabatic(ig,l)=min(wadiabatic(ig,l),wlim) ENDDO !! ********************************************************************** !! 2.2 case 2: daytime slope wind scheme IF ((mu0(ig) .gt. mu0lim) .and. .not. storm(ig)) then scheme(ig)=2 alpha(ig) = alpha_hmons(ig) ! interpolate the env. temperature above the mountain top call intep_vtemp(nlayer,hmons(ig),zt(ig,:),pzlay(ig,:), & envtemp,slpbg(ig)) envt(ig,:)=envtemp(:) !for diagnosis ! second: estimate the vertical velocity at boundraies of each layer !wspeed(ig,1)=0. ! at surface, already initialized !!!!!!!the first layer of atmosphere!!!!!!!!!! IF (slpbg(ig) .gt. 0.) THEN !only postive buoyancy ! if slpbg(ig) lt 0, means the slope flow is colder than env. (night or ! early morning ?) !!!!!!!!!!! ! ideal method !w1(ig)=-sqrt(g*slpbg(ig)/zt(ig,1)*hmons(ig))*sin(zsig(ig)) ! new scheme, simply proportional to temperature and ! mountain height w0(ig)=-1.5e-4*g*slpbg(ig)/zt(ig,1)*hmons(ig) ! otherwise, w0(ig) =0. wspeed(ig,2)=w0(ig) ELSE wspeed(ig,2)=wadiabatic(ig,2) !! for early morning ? ENDIF ! prepare the integration, NOTE: if w is too small, may have artificials IF (abs(wspeed(ig,2)) .lt. 0.01 ) & wspeed(ig,2)=sign(0.01,wspeed(ig,2)) newzt(ig,1)= zt(ig,1) !temperature of the first layer atmosphere ! above the mountain top (radiative ! equilibrium on Mars) ! estimate the vertical velocities ! if w0(ig) >= 0, means downward motion, no upslope winds, we switch to ! assume that the extra heating integrally convert to ! vertical motion. if ( w0(ig) .ge. 0 ) then !! normal, it is impossible, !! because mu(ig)>0.1 here do l=3,nlayer wspeed(ig,l)=wadiabatic(ig,l) enddo else ! estimate the velocities by taking into account the heating due ! to storm dust, the cooling due to vertical motion ... !!!!!!!!!!!the simple scheme!!!!!!!!! do l=2,nlayer-1 !if superadiabatic layer if ( zdtvert(l) .lt. -g/cpp) then !case 1 ! test, also decrease adiabatically ? !newzt(ig,l)= & ! zt(ig,l-1)-g/cpp*(pzlay(ig,l)-pzlay(ig,l-1)) newzt(ig,l)=zt(ig,l) !wspeed(ig,l+1)=wspeed(ig,l) else !not superadiabatic newzt(ig,l)=newzt(ig,l-1)+(deltahr(ig,l)/ & (-wspeed(ig,l))-g/cpp)* & (pzlay(ig,l)-pzlay(ig,l-1)) endif ! end of if superadiabatic or not !wtemp(ig,l+1)=wspeed(ig,l)**2+2.*g*(pzlev(ig,l+1) & ! -pzlev(ig,l))*(k1*(newzt(ig,l) & ! -envtemp(l))/envtemp(l)) wtemp(ig,l+1)=(1.-2.*k1*(pzlev(ig,l+1)-pzlev(ig,l)))*& wspeed(ig,l)**2+2.*k2*g*(pzlev(ig,l+1) & -pzlev(ig,l))*((newzt(ig,l) & -envtemp(l))/envtemp(l)) if (wtemp(ig,l+1) .gt. 0.) then !case 2 wspeed(ig,l+1)=-sqrt(wtemp(ig,l+1)) ! if |wspeed| < |wadiabatic| then go to wadiabatic if (wspeed(ig,l+1) .gt. wadiabatic(ig,l+1)) then do ll=l,nlayer-1 newzt(ig,ll)=envtemp(ll) wspeed(ig,ll+1)=wadiabatic(ig,ll+1) enddo exit endif ! avoid artificials if (abs(wspeed(ig,l+1)) .lt. 0.01 ) & wspeed(ig,l+1)=sign(0.01,wspeed(ig,l+1)) else if (l .lt. nlayer) then !case 3 do ll=l,nlayer-1 newzt(ig,ll)=envtemp(ll) wspeed(ig,ll+1)=wadiabatic(ig,ll+1) enddo !overshot exit else wspeed(ig,l+1)=0. endif enddo endif !w0 ELSE !! ********************************************************************** !! 2.3 case 3: storm=true if (storm(ig)) then scheme(ig)=3 alpha(ig) = totstormfract(ig) do l=1,nlayer wspeed(ig,l)=wadiabatic(ig,l) enddo endif ! storm=1 ENDIF ! rds or slp !!!!!!!! estimate the amount of dust for diagnostics DO l=1,nlayer ! transfer background dust + storm dust (concentrated) zqstormdustslp(ig,l) =zqi_mass(ig,l)+ & zqstorm_mass(ig,l)/alpha(ig) znstormdustslp(ig,l) =zqi_number(ig,l)+ & zqstorm_number(ig,l)/alpha(ig) zqdustslp(ig,l) = zqi_mass(ig,l) zndustslp(ig,l) = zqi_number(ig,l) rdsdustqvl0(ig,l)=zqstormdustslp(ig,l) !for diagnosis ENDDO ! ********************************************************************** ! 3. Vertical transport ! ********************************************************************** do l=1,nlayer masse(l)=(pplev(ig,l)-pplev(ig,l+1))/g enddo !Estimation of "dtmax" (s) to be used for vertical transport dtmax=ptimestep !secu is a margin on subtimstep to avoid dust crossing many layers do l=2,nlayer if (wspeed(ig,l).lt.0.) then ! case up dtmax=min(dtmax,(pzlev(ig,l)-pzlev(ig,l-1))/ & (secu*abs(wspeed(ig,l)))) else if (wspeed(ig,l).gt.0.) then dtmax=min(dtmax,(pzlev(ig,l+1)-pzlev(ig,l))/ & (secu*abs(wspeed(ig,l)))) endif enddo nsubtimestep= int(ptimestep/dtmax) subtimestep=ptimestep/float(nsubtimestep) do l=1,nlayer wrds(l)=wspeed(ig,l)*pplev(ig,l)/(r*ztlev(l)) & *subtimestep enddo do l=1,nlayer zqstorm_mass_col(l)= zqstormdustslp(ig,l) !zqstorm_mass(ig,l) zqstorm_number_col(l)=znstormdustslp(ig,l) ! zqstorm_number(ig,l) enddo do tsub=1,nsubtimestep wqrdsmass(:)=0. wqrdsnumber(:)=0. CALL vl_storm(nlayer,zqstorm_mass_col,2., & masse,wrds ,wqrdsmass) CALL vl_storm(nlayer,zqstorm_number_col,2., & masse,wrds ,wqrdsnumber) enddo !!!!!generate the "extra" dust do l=1,nlayer rdsdustqvl1(ig,l)=zqstorm_mass_col(l) ! for diagnosis ! extra dust = storm dust !zqdustslp(ig,l)=zqi_mass(ig,l) !(1.-alpha(ig))*zqi_mass(ig,l) !zndustslp(ig,l)=zqi_number(ig,l) !(1.-alpha(ig))*zqi_number(ig,l) !zqstorm_mass_col(l)=alpha(ig)*zqstorm_mass_col(l) !zqstorm_number_col(l)=alpha(ig)*zqstorm_number_col(l) !with compensation if (zqstorm_mass_col(l) .lt. zqi_mass(ig,l) ) then ! the following two equations are easier to understand zqdustslp(ig,l)=(1.-alpha(ig))*zqi_mass(ig,l)+alpha(ig)* & zqstorm_mass_col(l) zndustslp(ig,l)=(1.-alpha(ig))*zqi_number(ig,l)+alpha(ig)& *zqstorm_number_col(l) !with a bug, should be zqi+alpha**** !zqdustslp(ig,l)=zqi_mass(ig,l)-alpha(ig)* & ! (zqstorm_mass_col(l)-zqi_mass(ig,l)) !zndustslp(ig,l)=zqi_number(ig,l)-alpha(ig)* & ! (zqstorm_number_col(l)-zqi_number(ig,l)) zqstorm_mass_col(l)=0. zqstorm_number_col(l)=0. else zqstorm_mass_col(l)=alpha(ig)* & (zqstorm_mass_col(l)-zqi_mass(ig,l)) zqstorm_number_col(l)=alpha(ig)* & (zqstorm_number_col(l)-zqi_number(ig,l)) ! the mass mixing ratio of environmental dust doesn't change. endif !diagnose stormdust_m1(ig,l)=zqstorm_mass_col(l) enddo !======================================================================= ! calculate the tendencies due to vertical transport do l=1,nlayer ! tendencies due to vertical transport zdqvlstorm_mass(ig,l)= (zqstorm_mass_col(l)- & zqstorm_mass(ig,l)) /ptimestep zdqvlstorm_number(ig,l)= (zqstorm_number_col(l)- & zqstorm_number(ig,l)) /ptimestep zdqenv_mass(ig,l)=(zqdustslp(ig,l)-zqi_mass(ig,l))/ptimestep zdqenv_number(ig,l)=(zndustslp(ig,l)-zqi_number(ig,l)) & /ptimestep ! chao for output only !qstormdustvl1(ig,l)=zqstorm_mass_col(l) !nstormdustvl1(ig,l)=zqstorm_number_col(l) !stormdust_m_col1(ig)=stormdust_m_col1(ig)+zqstorm_mass_col(l) & ! *(pplev(ig,l)-pplev(ig,l+1))/g !rdsdustqvl1(ig,l)=zqstorm_mass_col(l) enddo ! ********************************************************************** ! 4. Detrainment: convert dust storm to background dust ! ********************************************************************** do l=1,nlayer ! compute the coefficient of detrainment if ((max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1))) .lt. & detrainlim) .or. (zqdustslp(ig,l) .gt. & 10000.*zqstorm_mass_col(l))) then coefdetrain=1. else if (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1))) & .le. wlim) then ! case where detrainment depends on vertical wind ! coefdetrain=0.5*(((1-coefmin)/(detrainlim-3.)**2)* & ! (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1)))-3.)**2 & ! +coefmin) coefdetrain=1.*(((1-coefmin)/(detrainlim-wlim)**2)* & (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1)))-wlim)**2 & +coefmin) !coefdetrain=0.5 else if (max(abs(wspeed(ig,l)),abs(wspeed(ig,l+1))).gt. 10. )& then coefdetrain=0.025 else coefdetrain=coefmin endif detrainment(ig,l)=coefdetrain !for diagnose ! Calculate tendancies corresponding to pdq after detrainement ! pdqdet = tendancy corresponding to detrainment only zdqdetstorm_mass(ig,l)=-coefdetrain*zqstorm_mass_col(l) & /ptimestep zdqdetstorm_number(ig,l)=-coefdetrain*zqstorm_number_col(l) & /ptimestep ! pdqrds ( tendancy corresponding to vertical transport and ! detrainment) = zdqvlstorm + pdqdet pdqrds(ig,l,igcm_stormdust_mass)=zdqdetstorm_mass(ig,l) & +zdqvlstorm_mass(ig,l) pdqrds(ig,l,igcm_stormdust_number)=zdqdetstorm_number(ig,l) & +zdqvlstorm_number(ig,l) pdqrds(ig,l,igcm_dust_mass)= zdqenv_mass(ig,l) & -zdqdetstorm_mass(ig,l) pdqrds(ig,l,igcm_dust_number)= zdqenv_number(ig,l) & -zdqdetstorm_number(ig,l) !diagnose stormdust_m2(ig,l)=zqstorm_mass_col(l)-coefdetrain*zqstorm_mass_col(l) enddo ! nlayer ! endif !======================================================================= enddo ! end do ig=1,ngrid ! !chao check conservation here ! do l=1,nlayer ! do ig=1,ngrid ! totdust0(ig)=totdust0(ig)+ & ! zq(ig,l,igcm_stormdust_mass)* & ! ((pplev(ig,l) - pplev(ig,l+1)) / g) & ! + zq(ig,l,igcm_dust_mass)* & ! ((pplev(ig,l) - pplev(ig,l+1)) / g) ! totdust1(ig)=totdust1(ig)+ & ! (zq(ig,l,igcm_stormdust_mass) + & ! pdqrds(ig,l,igcm_stormdust_mass)*ptimestep)* & ! ((pplev(ig,l) - pplev(ig,l+1)) / g) & ! + ( zq(ig,l,igcm_dust_mass)+ & ! pdqrds(ig,l,igcm_dust_mass)*ptimestep)* & ! ((pplev(ig,l) - pplev(ig,l+1)) / g) ! enddo ! enddo ! call writediagfi(ngrid,'totdust0','total dust before rds', & ! ' ',2,totdust0) ! call writediagfi(ngrid,'totdust1','total dust after rds', & ! ' ',2,totdust1) !output for diagnosis call WRITEDIAGFI(ngrid,'detrainment', & 'coefficient of detrainment',' ',3,detrainment) !call WRITEDIAGFI(ngrid,'qstormvl1','mmr of stormdust after rds_vl', & ! & 'kg/kg',3,qstormdustvl1) call WRITEDIAGFI(ngrid,'lapserate','lapse rate in the storm', & & 'k/m',3,lapserate) call WRITEDIAGFI(ngrid,'deltahr','extra heating rates', & & 'K/s',3,deltahr) call WRITEDIAGFI(ngrid,'wold', & 'wind generated from adiabatic colling ', & & 'm/s',3,wadiabatic) call WRITEDIAGFI(ngrid,'newzt','perturbated temperature', & & 'K/s',3,newzt) call WRITEDIAGFI(ngrid,'zt','unperturbated temperature', & & 'K/s',3,zt) call WRITEDIAGFI(ngrid,'wtemp','under square root', & & 'K/s',3,wtemp) !call WRITEDIAGFI(ngrid,'stormdust_m_col1','mass of stormdust after rds_vl', & ! & 'kg/kg',2,stormdust_m_col1) !call WRITEDIAGFI(ngrid,'temprds','temp for calculating zdtvert', & ! & 'k',3,temprds) call WRITEDIAGFI(ngrid,'stormdust_m0','mass col of stormdust before rds_vl', & & 'kg/kg',3,stormdust_m0) call WRITEDIAGFI(ngrid,'stormdust_m1','mass col of stormdust after rds_vl', & & 'kg/kg',3,stormdust_m1) call WRITEDIAGFI(ngrid,'stormdust_m2','mass col of stormdust after rds_vl', & & 'kg/kg',3,stormdust_m2) ! call writediagfi(ngrid,'wslp','estimated slope winds','m/s',3,wslp) ! call writediagfi(ngrid,'wslp2','estimated slope winds2','m/s',3,wslp2) ! call writediagfi(ngrid,'zhb','estimated slope winds2','m/s',3,zhb) ! call writediagfi(ngrid,'bruntf','bouyancy frequency',' ',3,bruntf) ! call writediagfi(ngrid,'slpdepth','slope depth','m',2,slpdepth) ! call writediagfi(ngrid,'slpu','perbulation u','m/s',3,slpu) ! call writediagfi(ngrid,'slpzh','perbulation zh',' ',3,slpzh) ! call writediagfi(ngrid,'zqslp','zq in rocketduststorm','ikg/kg',3, & ! zq(:,:,igcm_dust_mass)) ! call writediagfi(ngrid,'zrdsqslp','zq in rocketduststorm','ikg/kg',3, & ! zq(:,:,igcm_stormdust_mass)) ! call writediagfi(ngrid,'wslplev','estimated slope winds','m/s',3,wslplev) ! call writediagfi(ngrid,'slope','identified slope wind effect',' ',2,slope) call writediagfi(ngrid,'w0','max of slope wind',' ',2,w0) ! call writediagfi(ngrid,'w1','max of slope wind',' ',2,w1) call writediagfi(ngrid,'mu0','cosine of solar incidence angle',& ' ',2,mu0) ! call writediagfi(ngrid,'storm','identified rocket dust storm',& ! ' ',2,real(storm)) call writediagfi(ngrid,'scheme','which scheme',& ' ',2,real(scheme)) call writediagfi(ngrid,'alpha','coefficient alpha',' ',2,alpha) ! call writediagfi(ngrid,'q2rds','alpha zq',' ', & ! 3,q2rds) call writediagfi(ngrid,'rdsdustqvl0','not vl storm slp', & 'kg/kg',3,zqstormdustslp) call writediagfi(ngrid,'rdsdustqvl1','vled storm slp', & 'kg/kg',3,rdsdustqvl1) call writediagfi(ngrid,'dustqvl0','not vl slp', & 'kg/kg',3,zqi_mass) call writediagfi(ngrid,'dustqvl1','vled slp', & 'kg/kg',3,zqdustslp) ! call WRITEDIAGFI(ngrid,'lmax_th2', & ! 'hauteur du thermique','point', & ! 2,real(lmax_th(:))) ! call WRITEDIAGFI(ngrid,'zmax_th2', & ! 'hauteur du thermique','m', & ! 2,zmax_th) ! call writediagfi(ngrid,'lslope','lenght of slope',' ',2,lslope) ! call writediagfi(ngrid,'hmon','identified slope wind effect',' ',2,hmon) call writediagfi(ngrid,'envt','interpolated env. temp.', & 'K',3,envt) call writediagfi(ngrid,'hmons','identified slope wind effect', & ' ',2,hmons) ! call writediagfi(ngrid,'slpbg','temp. diff along slope', & ! ' ',2,slpbg) ! call writediagfi(ngrid,'zmea','identified slope wind effect', & ! ' ',2,zmea) ! call writediagfi(ngrid,'zsig','identified slope wind effect', & ! ' ',2,zsig) ! call writediagfi(ngrid,'zhslpenv','difference of zh above mons',' ',2,zhslpenv) ! call writediagfi(ngrid,'lref','identified slope wind effect',' ',2,real(lref)) END SUBROUTINE rocketduststorm !******************************************************************************** !******************************************************************************** SUBROUTINE vl_storm(nlay,q,pente_max,masse,w,wq) ! ! Auteurs: P.Le Van, F.Hourdin, F.Forget ! ! ******************************************************************** ! Shema d'advection " pseudo amont " dans la verticale ! pour appel dans la physique (sedimentation) ! ******************************************************************** ! q rapport de melange (kg/kg)... ! masse : masse de la couche Dp/g ! w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2) ! pente_max = 2 conseillee ! ! ! -------------------------------------------------------------------- IMPLICIT NONE ! ! Arguments: ! ---------- integer,intent(in) :: nlay ! number of atmospheric layers real masse(nlay),pente_max REAL q(nlay) REAL w(nlay) REAL wq(nlay+1) ! ! Local ! --------- ! INTEGER i,l,j,ii ! real dzq(nlay),dzqw(nlay),adzqw(nlay),dzqmax real newmasse real sigw, Mtot, MQtot integer m REAL SSUM,CVMGP,CVMGT integer ismax,ismin ! On oriente tout dans le sens de la pression c'est a dire dans le ! sens de W do l=2,nlay dzqw(l)=q(l-1)-q(l) adzqw(l)=abs(dzqw(l)) enddo do l=2,nlay-1 #ifdef CRAY dzq(l)=0.5* , cvmgp(dzqw(l)+dzqw(l+1),0.,dzqw(l)*dzqw(l+1)) #else if(dzqw(l)*dzqw(l+1).gt.0.) then dzq(l)=0.5*(dzqw(l)+dzqw(l+1)) else dzq(l)=0. endif #endif dzqmax=pente_max*min(adzqw(l),adzqw(l+1)) dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l)) enddo dzq(1)=0. dzq(nlay)=0. ! write(*,*),'TB14 wq before up',wq(1,:) ! write(*,*),'TB14 q before up',q(1,:) ! --------------------------------------------------------------- ! .... calcul des termes d'advection verticale ....... ! --------------------------------------------------------------- ! calcul de - d( q * w )/ d(sigma) qu'on ajoute a dq pour calculer dq ! ! No flux at the model top: wq(nlay+1)=0. ! 1) Compute wq where w < 0 (up) (NOT USEFUL FOR SEDIMENTATION) ! =============================== ! Surface flux up: if(w(1).lt.0.) wq(1)=0. ! warning : not always valid do l = 1,nlay-1 ! loop different than when w>0 if(w(l+1).le.0)then ! Regular scheme (transfered mass < 1 layer) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if(-w(l+1).le.masse(l))then sigw=w(l+1)/masse(l) wq(l+1)=w(l+1)*(q(l)-0.5*(1.+sigw)*dzq(l)) ! Extended scheme (transfered mass > 1 layer) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else m = l-1 Mtot = masse(m+1) MQtot = masse(m+1)*q(m+1) if (m.le.0)goto 77 do while(-w(l+1).gt.(Mtot+masse(m))) ! do while(-w(l+1).gt.Mtot) m=m-1 Mtot = Mtot + masse(m+1) MQtot = MQtot + masse(m+1)*q(m+1) if (m.le.0)goto 77 end do 77 continue if (m.gt.0) then sigw=(w(l+1)+Mtot)/masse(m) wq(l+1)= (MQtot + (-w(l+1)-Mtot)* & (q(m)-0.5*(1.+sigw)*dzq(m)) ) else ! new w(l+1) = -Mtot wq(l+1) = -MQtot ! write(*,*) 'TB14 MQtot = ',MQtot,l ! old ! wq(l+1)= (MQtot + (-w(l+1)-Mtot)*qm(1)) ! write(*,*) 'a rather weird situation in vlz_fi !' ! stop end if endif endif ! w<0 (up) enddo do l = 1,nlay-1 ! loop different than when w>0 q(l)=q(l) + (wq(l+1)-wq(l))/masse(l) enddo ! write(*,*),'TB14 masse',masse(1,:) ! write(*,*),'TB14 wq before down after up',wq(1,:) ! write(*,*),'TB14 q before down',q(1,:) ! 2) Compute wq where w > 0 (down) (ALWAYS FOR SEDIMENTATION) ! =============================== ! Initialisation wq = 0 to consider now only downward flux wq(:)=0. ! do l = 1,nlay ! loop different than when w<0 if(w(l).gt.0.)then ! Regular scheme (transfered mass < 1 layer) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if(w(l).le.masse(l))then sigw=w(l)/masse(l) wq(l)=w(l)*(q(l)+0.5*(1.-sigw)*dzq(l)) ! write(*,*),'TB14 wq after up',wq(1,:) ! Extended scheme (transfered mass > 1 layer) ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ else m=l Mtot = masse(m) MQtot = masse(m)*q(m) if(m.ge.nlay)goto 88 do while(w(l).gt.(Mtot+masse(m+1))) m=m+1 Mtot = Mtot + masse(m) MQtot = MQtot + masse(m)*q(m) if(m.ge.nlay)goto 88 end do 88 continue if (m.lt.nlay) then sigw=(w(l)-Mtot)/masse(m+1) wq(l)=(MQtot + (w(l)-Mtot)* & (q(m+1)+0.5*(1.-sigw)*dzq(m+1)) ) else w(l) = Mtot wq(l) = MQtot end if end if end if ! w>0 (down) enddo do l = 1,nlay ! loop different than when w<0 q(l)=q(l) + (wq(l+1)-wq(l))/masse(l) enddo end subroutine vl_storm !******************************************************************************** SUBROUTINE intep_vtemp(nlayer,hm,temp,zlay,envtemp,slpb) USE comcstfi_h, only: g,cpp implicit none ! this subroutine is using for obtaining the environmental ! temperature profile when a subgrid mountain exists. ! input: integer,intent(in) :: nlayer real,intent(in) :: hm ! the height of mountain real,intent(in) :: temp(nlayer) !large scale temp. profile of one mesh real,intent(in) :: zlay(nlayer) ! height at the middle of each layer ! output: real,intent(out) :: envtemp(nlayer) ! environment temperature real,intent(out) :: slpb !the temperature difference between slope and ! env. at the half height of mountain ! local variables integer l,il,tmpl integer lnew !the layer of atmosphere above the mountain ! corresponding to the env. (for buoyancy ! calc. ) real newh(nlayer) !height at the middle of each layer ! account for the exist of mountain ! real g,cpp real halfh ! half the height of a mountain !initilize the array lnew=0 newh(:)=0. envtemp(:)=0. tmpl=1 do l=1,nlayer newh(l)=hm+zlay(l) do il=tmpl,nlayer-1 !MV18: added the -1 if (newh(l) .ge. zlay(il) .and. newh(l) .lt. zlay(il+1))then ! find the corresponding layer lnew=il ! interpolate envtemp(l) = temp(il)+(newh(l)-zlay(lnew))*& (temp(il+1)-temp(il))/(zlay(il+1)-zlay(il)) exit !go to the next layer else if (newh(l) .ge. zlay(nlayer)) then ! higher than the highest layer lnew=nlayer envtemp(l)=temp(nlayer) endif enddo ! this can accelerate the loop tmpl=lnew enddo halfh=0.5*hm if (halfh .le. zlay(1) ) then slpb=0. else if (halfh .gt. zlay(nlayer)) then !normally, impossible for halfh gt zlay(l), anyway... tmpl=nlayer !difference between surface and atmosphere (env.) slpb=temp(1)-(temp(nlayer-1)+((halfh-zlay(nlayer-1))* & (temp(tmpl)-temp(tmpl-1))/(zlay(tmpl)-zlay(tmpl-1)))) else do l=1,nlayer-1 if ((halfh .gt. zlay(l)) .and. (halfh .le. zlay(l+1)))then tmpl= l slpb=temp(1)-(temp(tmpl)+(halfh-zlay(tmpl))* & (temp(tmpl+1)-temp(tmpl))/(zlay(tmpl+1)-zlay(tmpl))) endif enddo endif end subroutine intep_vtemp ! initialization module variables subroutine ini_rocketduststorm_mod(ngrid) implicit none integer, intent(in) :: ngrid allocate(dustliftday(ngrid)) allocate(alpha_hmons(ngrid)) end subroutine ini_rocketduststorm_mod subroutine end_rocketduststorm_mod implicit none if (allocated(dustliftday)) deallocate(dustliftday) if (allocated(alpha_hmons)) deallocate(alpha_hmons) end subroutine end_rocketduststorm_mod END MODULE rocketduststorm_mod