diff --git a/src/bulkmicro_kk.f90 b/src/bulkmicro_kk.f90 index 9a35dad4..4e6e2f2d 100644 --- a/src/bulkmicro_kk.f90 +++ b/src/bulkmicro_kk.f90 @@ -62,7 +62,7 @@ subroutine do_bulkmicro_kk thlpmcr, qtpmcr, qrp, Nrp) call bulkmicrotend #ifdef DALES_GPU - call sedimentation_rain_gpu(qr, Nr, rhof, dzf, qrbaes, qrroof, qrmask, delt, & + call sedimentation_rain_gpu(qr, Nr, rhof, dzf, qrbase, qrroof, qrmask, delt, & Dvr, xr, qrp, Nrp, precep) #else call sedimentation_rain(qr, Nr, rhof, dzf, qrbase, qrroof, qrmask, delt, & @@ -393,8 +393,8 @@ subroutine sedimentation_rain(qr, Nr, rhof, dzf, qrbase, qrroof, qrmask, & do j = 2, j1 do i = 2, i1 if (qrmask(i,j,k)) then - sed_qr = max(0.0_field_r, 0.006_field_r * 1E6 * Dvr(i,j,k) - 0.2_field_r) * qr_spl(i,j,k) * rhof(k) - sed_Nr = max(0.0_field_r, 0.0035_field_r * 1E6 * Dvr(i,j,k) - 0.1_field_r) * Nr_spl(i,j,k) + sed_qr = max(0.0_field_r, 0.006_field_r * 1E6_field_r * Dvr(i,j,k) - 0.2_field_r) * qr_spl(i,j,k) * rhof(k) + sed_Nr = max(0.0_field_r, 0.0035_field_r * 1E6_field_r * Dvr(i,j,k) - 0.1_field_r) * Nr_spl(i,j,k) qr_spl(i,j,k) = qr_spl(i,j,k) - sed_qr * dt_spl / (dzf(k) * rhof(k)) Nr_spl(i,j,k) = Nr_spl(i,j,k) - sed_Nr * dt_spl / dzf(k) @@ -539,7 +539,7 @@ subroutine sedimentation_rain_gpu(qr, Nr, rhof, dzf, qrbase, qrroof, qrmask, & do j = 2, j1 do i = 2, i1 if (qrmask(i,j,k)) then - precep(i,j,k) = max(0.0_field_r, 0.006_field_r * 1E6 * Dvr(i,j,k) - 0.2_field_r) * qr_spl(i,j,k) + precep(i,j,k) = max(0.0_field_r, 0.006_field_r * 1E6_field_r * Dvr(i,j,k) - 0.2_field_r) * qr_spl(i,j,k) endif enddo enddo @@ -557,8 +557,8 @@ subroutine sedimentation_rain_gpu(qr, Nr, rhof, dzf, qrbase, qrroof, qrmask, & do j = 2, j1 do i = 2, i1 if (qrmask(i,j,k)) then - sed_qr = max(0.0_field_r, 0.006_field_r *1E6 * Dvr(i,j,k) - 0.2_field_r) * qr_spl(i,j,k) * rhof(k) - sed_Nr = max(0.0_field_r, 0.0035_field_r *1E6 * Dvr(i,j,k) - 0.1_field_r) * Nr_spl(i,j,k) + sed_qr = max(0.0_field_r, 0.006_field_r *1E6_field_r * Dvr(i,j,k) - 0.2_field_r) * qr_spl(i,j,k) * rhof(k) + sed_Nr = max(0.0_field_r, 0.0035_field_r *1E6_field_r * Dvr(i,j,k) - 0.1_field_r) * Nr_spl(i,j,k) qr_tmp(i,j,k) = qr_tmp(i,j,k) - sed_qr * dt_spl / (dzf(k) * rhof(k)) Nr_tmp(i,j,k) = Nr_tmp(i,j,k) - sed_Nr * dt_spl / dzf(k) @@ -572,8 +572,8 @@ subroutine sedimentation_rain_gpu(qr, Nr, rhof, dzf, qrbase, qrroof, qrmask, & do j = 2, j1 do i = 2, i1 if (qrmask(i,j,k)) then - sed_qr = max(0.0_field_r, 0.006_field_r *1E6 * Dvr(i,j,k) - 0.2_field_r) * qr_spl(i,j,k) * rhof(k) - sed_Nr = max(0.0_field_r, 0.0035_field_r *1E6 * Dvr(i,j,k) - 0.1_field_r) * Nr_spl(i,j,k) + sed_qr = max(0.0_field_r, 0.006_field_r *1E6_field_r * Dvr(i,j,k) - 0.2_field_r) * qr_spl(i,j,k) * rhof(k) + sed_Nr = max(0.0_field_r, 0.0035_field_r *1E6_field_r * Dvr(i,j,k) - 0.1_field_r) * Nr_spl(i,j,k) !$acc atomic update qr_tmp(i,j,k) = qr_tmp(i,j,k) - sed_qr*dt_spl/(dzf(k)*rhof(k)) diff --git a/src/bulkmicro_sb.f90 b/src/bulkmicro_sb.f90 index 378d079e..5f25b942 100644 --- a/src/bulkmicro_sb.f90 +++ b/src/bulkmicro_sb.f90 @@ -87,9 +87,9 @@ subroutine do_bulkmicro_sb mur, xr, qrp, Nrp, delt, qtpmcr, thlpmcr) call bulkmicrotend #ifdef DALES_GPU - call sedimentation_rain_gpu(qr, Nr, rhof, dzf, qrbase, qrroof, Dvr, lbdr, & - mur, delt, xr, l_lognormal, l_mur_cst, & - mur_cst, qrp, Nrp, qrmask, precep) + call sedimentation_rain_gpu(qr, Nr, rhof, dzf, qrbase, qrroof, qrmask, & + l_lognormal, l_mur_cst, mur_cst, delt, Dvr, lbdr, & + mur, xr, qrp, Nrp, precep) #else call sedimentation_rain(qr, Nr, rhof, dzf, qrbase, qrroof, qrmask, & l_lognormal, l_mur_cst, mur_cst, delt, Dvr, lbdr, & @@ -583,8 +583,8 @@ subroutine sedimentation_rain(qr, Nr, rhof, dzf, qrbase, qrroof, qrmask, & do i = 2, i1 if (qrmask(i,j,k)) then - wfall_qr = max(0., (a_tvsb - b_tvsb * (1 + c_tvsb / lbdr(i,j,k))**(-1 * (mur(i,j,k)+4)))) - wfall_Nr = max(0., (a_tvsb - b_tvsb * (1 + c_tvsb / lbdr(i,j,k))**(-1 * (mur(i,j,k)+1)))) + wfall_qr = max(0._field_r, (a_tvsb - b_tvsb * (1 + c_tvsb / lbdr(i,j,k))**(-1 * (mur(i,j,k)+4)))) + wfall_Nr = max(0._field_r, (a_tvsb - b_tvsb * (1 + c_tvsb / lbdr(i,j,k))**(-1 * (mur(i,j,k)+1)))) sed_qr = wfall_qr * qr_spl(i,j,k) * rhof(k) ! m/s * kg/m3 sed_Nr = wfall_Nr * Nr_spl(i,j,k) @@ -1049,4 +1049,4 @@ real function erfint(beta, D, D_min, D_max, sig2,nnn ) erfint = max(0., D**nn*exp(0.5*nn**2*sig2)*0.5*(erfymax-erfymin)) end function erfint -end module bulkmicro_sb \ No newline at end of file +end module bulkmicro_sb diff --git a/src/modgpu.f90 b/src/modgpu.f90 index ee31fc3c..535c3f9f 100644 --- a/src/modgpu.f90 +++ b/src/modgpu.f90 @@ -49,7 +49,7 @@ subroutine update_gpu dudxls, dudyls, dudtls, dvdxls, dvdyls, & dvdtls, dthvdz, qvsl, qvsi, esl, qsat use modglobal, only: dzf, dzh, zh, zf, delta, deltai, & - dsv, rd, rv, esatmtab, esatitab, esatltab + rd, rv, esatmtab, esatitab, esatltab use modsurfdata, only: z0m, z0h, obl, tskin, qskin, Cm, Cs, & ustar, dudz, dvdz, thlflux, qtflux, & dqtdz, dthldz, svflux, svs, horv, ra, rs, wsvsurf @@ -78,7 +78,7 @@ subroutine update_gpu !$acc& dthldtls, dqtdxls, dqtdyls, dqtdtls, & !$acc& dudxls, dudyls, dudtls, dvdxls, dvdyls, & !$acc& dvdtls, dthvdz, qvsl, qvsi, esl, qsat, & - !$acc& dzf, dzh, zh, zf, delta, deltai, dsv, & + !$acc& dzf, dzh, zh, zf, delta, deltai, & !$acc& z0m, z0h, obl, tskin, qskin, Cm, Cs, & !$acc& ustar, dudz, dvdz, thlflux, qtflux, & !$acc& dqtdz, dthldz, svflux, svs, horv, ra, rs, wsvsurf, & @@ -122,7 +122,7 @@ subroutine update_host dudxls, dudyls, dudtls, dvdxls, dvdyls, & dvdtls, dthvdz, qvsl, qvsi, esl, qsat use modglobal, only: dzf, dzh, zh, zf, delta, deltai, & - dsv, rd, rv, esatmtab, esatitab, esatltab + rd, rv, esatmtab, esatitab, esatltab use modsurfdata, only: z0m, z0h, obl, tskin, qskin, Cm, Cs, & ustar, dudz, dvdz, thlflux, qtflux, & dqtdz, dthldz, svflux, svs, horv, ra, rs, wsvsurf @@ -153,7 +153,7 @@ subroutine update_host !$acc& dthldtls, dqtdxls, dqtdyls, dqtdtls, & !$acc& dudxls, dudyls, dudtls, dvdxls, dvdyls, & !$acc& dvdtls, dthvdz, qvsl, qvsi, esl, qsat, & - !$acc& dzf, dzh, zh, zf, delta, deltai, dsv, & + !$acc& dzf, dzh, zh, zf, delta, deltai, & !$acc& z0m, z0h, obl, tskin, qskin, Cm, Cs, & !$acc& ustar, dudz, dvdz, thlflux, qtflux, & !$acc& dqtdz, dthldz, svflux, svs, horv, ra, rs, wsvsurf, & diff --git a/src/modopenboundary.f90 b/src/modopenboundary.f90 index cc36e1a4..598eaed1 100644 --- a/src/modopenboundary.f90 +++ b/src/modopenboundary.f90 @@ -787,7 +787,7 @@ subroutine applyboundaryf(a,sx,ex,sy,ey,sz,ez,ih,jh,ib,val,nx1,nx2,lmax0,turb,pr valtarget = (fp*val(j,k,itp)+fm*val(j,k,itm)+turb(j,k))*coefdir a(2-ih:sx-1,j+1,k) = ( 2.*dx*valtarget - & a(sx,j+1,k)*(coefdir*dx+2.*coefneu) ) / (coefdir*dx-2.*coefneu) - if(lmax0==1) a(sx-1,j+1,k) = max(0.,a(sx-1,j+1,k)) + if(lmax0==1) a(sx-1,j+1,k) = max(0._field_r,a(sx-1,j+1,k)) endif end do end do @@ -804,7 +804,7 @@ subroutine applyboundaryf(a,sx,ex,sy,ey,sz,ez,ih,jh,ib,val,nx1,nx2,lmax0,turb,pr valtarget = (fp*val(j,k,itp)+fm*val(j,k,itm)+turb(j,k))*coefdir a(ex+1:i1+ih,j+1,k) = ( 2.*dx*valtarget - & a(ex,j+1,k)*(coefdir*dx-2.*coefneu) ) / (coefdir*dx+2.*coefneu) - if(lmax0==1) a(ex+1,j+1,k) = max(a(ex+1,j+1,k),0.) + if(lmax0==1) a(ex+1,j+1,k) = max(a(ex+1,j+1,k),0._field_r) endif end do end do @@ -821,7 +821,7 @@ subroutine applyboundaryf(a,sx,ex,sy,ey,sz,ez,ih,jh,ib,val,nx1,nx2,lmax0,turb,pr valtarget = (fp*val(i,k,itp)+fm*val(i,k,itm)+turb(i,k))*coefdir a(i+1,2-jh:sy-1,k) = ( 2.*dy*valtarget - & a(i+1,sy,k)*(coefdir*dy+2.*coefneu) ) / (coefdir*dy-2.*coefneu) - if(lmax0==1) a(i+1,sy-1,k) = max(a(i+1,sy-1,k),0.) + if(lmax0==1) a(i+1,sy-1,k) = max(a(i+1,sy-1,k),0._field_r) endif end do end do @@ -838,7 +838,7 @@ subroutine applyboundaryf(a,sx,ex,sy,ey,sz,ez,ih,jh,ib,val,nx1,nx2,lmax0,turb,pr valtarget = (fp*val(i,k,itp)+fm*val(i,k,itm)+turb(i,k))*coefdir a(i+1,ey+1:j1+jh,k) = ( 2.*dy*valtarget - & a(i+1,ey,k)*(coefdir*dy-2.*coefneu) ) / (coefdir*dy+2.*coefneu) - if(lmax0==1) a(i+1,ey+1,k) = max(a(i+1,ey+1,k),0.) + if(lmax0==1) a(i+1,ey+1,k) = max(a(i+1,ey+1,k),0._field_r) endif end do end do @@ -862,7 +862,7 @@ subroutine applyboundaryf(a,sx,ex,sy,ey,sz,ez,ih,jh,ib,val,nx1,nx2,lmax0,turb,pr valtarget = (fp*val(i,j,itp)+fm*val(i,j,itm)+turb(i,j))*coefdir+ddz*coefneu a(i+1,j+1,ez) = ( 2.*dzh(ez)*valtarget - & a(i+1,j+1,ez-1)*(coefdir*dzh(ez)-2.*coefneu) ) / (coefdir*dzh(ez)+2.*coefneu) - if(lmax0==1) a(i+1,j+1,ez) = max(a(i+1,j+1,ez),0.) + if(lmax0==1) a(i+1,j+1,ez) = max(a(i+1,j+1,ez),0._field_r) endif end do end do @@ -881,7 +881,9 @@ subroutine applyboundaryh(ib,nx1,nx2,turb) integer, intent(in) :: nx1,nx2,ib real(field_r), intent(in), dimension(nx1,nx2) :: turb integer :: i,j,k,itmc,itmn,itpc,itpn,ipatch,jpatch,kpatch - real :: tm,tp,fpc,fmc,fpn,fmn,unext,uwallcurrent,ipos,jpos,tau + real :: tm,tp,fpc,fmc,fpn,fmn,unext,ipos,jpos,tau + real(field_r) :: uwallcurrent + itmc=1 itmn=1 if(ntboundary>1) then diff --git a/src/modsimpleice.f90 b/src/modsimpleice.f90 index 4885dfd8..e09018e9 100644 --- a/src/modsimpleice.f90 +++ b/src/modsimpleice.f90 @@ -302,7 +302,7 @@ subroutine autoconvert tc=tmp0(i,j,k)-tmelt ! Temperature wrt melting point times=min(1.e3,(3.56*tc+106.7)*tc+1.e3) ! Time scale for ice autoconversion auti=qli/times - aut = min(autl + auti,ql0(i,j,k)/delt) + aut = min(autl + auti,ql0(i,j,k)/delt*1.0) ! convert RHS to double for nvidia compiler qrp(i,j,k) = qrp(i,j,k)+aut qtpmcr(i,j,k) = qtpmcr(i,j,k)-aut thlpmcr(i,j,k) = thlpmcr(i,j,k)+(rlv/(cp*exnf(k)))*aut @@ -321,7 +321,7 @@ subroutine autoconvert autl=max(0.,timekessl*(qll-qll0)) tc=tmp0(i,j,k)-tmelt auti=max(0.,betakessi*exp(0.025*tc)*(qli-qli0)) - aut = min(autl + auti,ql0(i,j,k)/delt) + aut = min(autl + auti,ql0(i,j,k)/delt*1.0) ! convert RHS to double for nvidia compiler qrp(i,j,k) = qrp(i,j,k)+aut qtpmcr(i,j,k) = qtpmcr(i,j,k)-aut thlpmcr(i,j,k) = thlpmcr(i,j,k)+(rlv/(cp*exnf(k)))*aut @@ -368,7 +368,7 @@ subroutine accrete accr=(gaccrl+gaccri)*qrr/(qrr+1.e-9) accs=(gaccsl+gaccsi)*qrs/(qrs+1.e-9) accg=(gaccgl+gaccgi)*qrg/(qrg+1.e-9) - acc= min(accr+accs+accg,ql0(i,j,k)/delt) ! total growth by accretion + acc= min(accr+accs+accg,ql0(i,j,k)/delt*1.0) ! total growth by accretion ! convert RHS to double for nvidia compiler qrp(i,j,k) = qrp(i,j,k)+acc qtpmcr(i,j,k) = qtpmcr(i,j,k)-acc thlpmcr(i,j,k) = thlpmcr(i,j,k)+(rlv/(cp*exnf(k)))*acc @@ -411,7 +411,7 @@ subroutine evapdep evapdepg=(4.*pi/(betag*rhof(k)))*(ssi-1.)*ventg*thfun ! total growth by deposition and evaporation ! limit with qr and ql after accretion and autoconversion - devap= max(min(evapfactor*(evapdepr+evapdeps+evapdepg),ql0(i,j,k)/delt+qrp(i,j,k)),-qr(i,j,k)/delt-qrp(i,j,k)) + devap= max(min(evapfactor*(evapdepr+evapdeps+evapdepg),ql0(i,j,k)/delt+qrp(i,j,k)*1.0),-qr(i,j,k)/delt-qrp(i,j,k)*1.0) ! qrp(i,j,k) = qrp(i,j,k)+devap qtpmcr(i,j,k) = qtpmcr(i,j,k)-devap thlpmcr(i,j,k) = thlpmcr(i,j,k)+(rlv/(cp*exnf(k)))*devap diff --git a/src/modsimpleice2.f90 b/src/modsimpleice2.f90 index 97d6959c..e95da8d3 100644 --- a/src/modsimpleice2.f90 +++ b/src/modsimpleice2.f90 @@ -166,16 +166,16 @@ subroutine simpleice2 implicit none integer:: i,j,k - real:: qrsmall, qrsum - real :: qll,qli,ddisp,lwc,autl,tc,times,auti,aut ! autoconvert - real :: qrr,qrs,qrg, gaccrl,gaccsl,gaccgl,gaccri,gaccsi,gaccgi,accr,accs,accg,acc !accrete - real :: ssl,ssi,ventr,vents,ventg,thfun,evapdepr,evapdeps,evapdepg,devap !evapdep - real :: dt_spl,wfallmax,vtr,vts,vtg,vtf ! precipitation - real :: tmp_lambdar, tmp_lambdas, tmp_lambdag + real(field_r) :: qrsmall, qrsum + real(field_r) :: qll,qli,ddisp,lwc,autl,tc,times,auti,aut ! autoconvert + real(field_r) :: qrr,qrs,qrg, gaccrl,gaccsl,gaccgl,gaccri,gaccsi,gaccgi,accr,accs,accg,acc !accrete + real(field_r) :: ssl,ssi,ventr,vents,ventg,thfun,evapdepr,evapdeps,evapdepg,devap !evapdep + real(field_r) :: dt_spl,wfallmax,vtr,vts,vtg,vtf ! precipitation + real(field_r) :: tmp_lambdar, tmp_lambdas, tmp_lambdag integer :: jn integer :: n_spl !< sedimentation time splitting loop - real :: ilratio_,lambdar_,lambdas_,lambdag_, rsgratio_, sgratio_ ! local values instead of global arrays + real(field_r) :: ilratio_,lambdar_,lambdas_,lambdag_, rsgratio_, sgratio_ ! local values instead of global arrays logical :: qrmask_, qcmask_ logical :: rain_present, snow_present, graupel_present ! logicals for presence of different forms of water in the current cell logical :: any_qr, any_snow_graupel ! logicals for precense of any precipitation, and for presense of snow/graupel in the whole system @@ -356,9 +356,9 @@ subroutine simpleice2 ! qll=ql0(i,j,k)*ilratio(i,j,k) ! qli=ql0(i,j,k)-qll - autl=max(0.,timekessl*(qll-qll0)) + autl=max(0._field_r,timekessl*(qll-qll0)) tc=tmp0(i,j,k)-tmelt - auti=max(0.,betakessi*exp(0.025*tc)*(qli-qli0)) + auti=max(0._field_r,betakessi*exp(0.025_field_r*tc)*(qli-qli0)) aut = min(autl + auti,ql0(i,j,k)/delt) qrp(i,j,k) = qrp(i,j,k) + aut qtpmcr(i,j,k) = qtpmcr(i,j,k) - aut