diff --git a/src/advance_mod.f90 b/src/advance_mod.f90
index e9691ab957d1a353c8978eb7d17b802f32564d69..36cb4fd110663a59052cc1c81873595700394be1 100644
--- a/src/advance_mod.f90
+++ b/src/advance_mod.f90
@@ -752,12 +752,12 @@ subroutine petterssen_corr(itime,ipart)
v=(v-vold)*0.5
#ifdef ETA
- weta=(weta-woldeta)/2.
+ weta=(weta-woldeta)*0.5
call update_zeta(ipart,weta*real(part(ipart)%idt*ldirect))
if (part(ipart)%zeta.ge.1.) call set_zeta(ipart,1.-(part(ipart)%zeta-1.))
if (part(ipart)%zeta.eq.1.) call update_zeta(ipart,-eps_eta)
#else
- w=(w-wold)/2.
+ w=(w-wold)*0.5
call update_z(ipart,w*real(part(ipart)%idt*ldirect))
if (part(ipart)%z.lt.0.) call set_z(ipart,min(h-eps2,-1.*part(ipart)%z)) ! if particle below ground -> reflection
#endif
diff --git a/src/binary_output_mod.f90 b/src/binary_output_mod.f90
index df48ce7fe338b4ae8c96279e568fc74d5d844b68..39087de7cb252016dca4b4ebbdc843faefd19d95 100644
--- a/src/binary_output_mod.f90
+++ b/src/binary_output_mod.f90
@@ -806,9 +806,9 @@ subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
mind=memind(2)
do kz=1,numzgrid
if (kz.eq.1) then
- halfheight=outheight(1)/2.
+ halfheight=outheight(1)*0.5
else
- halfheight=(outheight(kz)+outheight(kz-1))/2.
+ halfheight=(outheight(kz)+outheight(kz-1))*0.5
endif
do kzz=2,nz
if ((height(kzz-1).lt.halfheight).and. &
@@ -1426,9 +1426,9 @@ subroutine concoutput_nest(itime,outnum)
mind=memind(2)
do kz=1,numzgrid
if (kz.eq.1) then
- halfheight=outheight(1)/2.
+ halfheight=outheight(1)*0.5
else
- halfheight=(outheight(kz)+outheight(kz-1))/2.
+ halfheight=(outheight(kz)+outheight(kz-1))*0.5
endif
do kzz=2,nz
if ((height(kzz-1).lt.halfheight).and. &
@@ -2065,9 +2065,9 @@ subroutine concoutput_inversion(itime,outnum,gridtotalunc,wetgridtotalunc, &
do kz=1,numzgrid
if (kz.eq.1) then
- halfheight=outheight(1)/2.
+ halfheight=outheight(1)*0.5
else
- halfheight=(outheight(kz)+outheight(kz-1))/2.
+ halfheight=(outheight(kz)+outheight(kz-1))*0.5
endif
do kzz=2,nz
if ((height(kzz-1).lt.halfheight).and. &
@@ -2612,9 +2612,9 @@ subroutine concoutput_inversion_nest(itime,outnum)
do kz=1,numzgrid
if (kz.eq.1) then
- halfheight=outheight(1)/2.
+ halfheight=outheight(1)*0.5
else
- halfheight=(outheight(kz)+outheight(kz-1))/2.
+ halfheight=(outheight(kz)+outheight(kz-1))*0.5
endif
do kzz=2,nz
if ((height(kzz-1).lt.halfheight).and. &
@@ -3086,9 +3086,9 @@ subroutine concoutput_sfc(itime,outnum,gridtotalunc,wetgridtotalunc, &
do kz=1,numzgrid
if (kz.eq.1) then
- halfheight=outheight(1)/2.
+ halfheight=outheight(1)*0.5
else
- halfheight=(outheight(kz)+outheight(kz-1))/2.
+ halfheight=(outheight(kz)+outheight(kz-1))*0.5
endif
do kzz=2,nz
if ((height(kzz-1).lt.halfheight).and. &
@@ -3729,9 +3729,9 @@ subroutine concoutput_sfc_nest(itime,outnum)
do kz=1,numzgrid
if (kz.eq.1) then
- halfheight=outheight(1)/2.
+ halfheight=outheight(1)*0.5
else
- halfheight=(outheight(kz)+outheight(kz-1))/2.
+ halfheight=(outheight(kz)+outheight(kz-1))*0.5
endif
do kzz=2,nz
if ((height(kzz-1).lt.halfheight).and. &
diff --git a/src/cbl_mod.f90 b/src/cbl_mod.f90
index 32592ac560895cb0fb9eae41a5a07fb4486fcd84..9a9331a637cb26da33dd7d3bdf3cf8c94decbe52 100644
--- a/src/cbl_mod.f90
+++ b/src/cbl_mod.f90
@@ -88,7 +88,7 @@ subroutine cbl(wp,zp,wst,h,rhoa,rhograd,sigmaw,dsigmawdz,tlw,ptot,Q,phi,ath,bth,
! stability transition function see Cassiani et al(2015) BLM
transition=1.
!if (ol.lt.-50) transition=((sin(((ol+100.)/100.)*pi))-1.)/2.
- if (-h/ol.lt.15) transition=((sin((((-h/ol)+10.)/10.)*pi)))/2.+0.5
+ if (-h/ol.lt.15) transition=((sin((((-h/ol)+10.)/10.)*pi)))*0.5+0.5
! momento secondo
@@ -284,7 +284,7 @@ subroutine reinit_particle(zp,wst,h,sigmaw,wp,nrand,ol)
z=zp/h
transition=1.
- if (-h/ol.lt.15) transition=((sin((((-h/ol)+10.)/10.)*pi)))/2.+0.5
+ if (-h/ol.lt.15) transition=((sin((((-h/ol)+10.)/10.)*pi)))*0.5+0.5
w2=sigmaw*sigmaw
w3=(((1.2*z*((1.-z)**(3./2.)))+eps)*wst**3)*transition
@@ -374,7 +374,7 @@ subroutine init_cbl_vel(idum,zp,wst,h,sigmaw,wp,ol,ithread)
transition=1.
- if (-h/ol.lt.15) transition=((sin((((-h/ol)+10.)/10.)*pi)))/2.+0.5 !see also in cbl.f90
+ if (-h/ol.lt.15) transition=((sin((((-h/ol)+10.)/10.)*pi)))*0.5+0.5 !see also in cbl.f90
w2=sigmaw*sigmaw
w3=(((1.2*z*((1.-z)**(3./2.)))+eps)*wst**3) *transition
diff --git a/src/cmapf_mod.f90 b/src/cmapf_mod.f90
index f5aa61e8eaa2935b8ca1d8a9dc084ff07feda00a..cf2e216c1385b7bdea912094c886e883030984d7 100644
--- a/src/cmapf_mod.f90
+++ b/src/cmapf_mod.f90
@@ -228,7 +228,7 @@ real function cgszll (strcmp, xlat)
endif
else
slat = sin(radpdg * xlat)
- ymerc = log((1. + slat) / (1. - slat))/2.
+ ymerc = log((1. + slat) / (1. - slat))*0.5
!efact = exp(ymerc)
!cgszll = 2. * strcmp(7) * exp (strcmp(1) * ymerc)
!c / (efact + 1./efact)
@@ -270,7 +270,7 @@ real function cgszxy (strcmp, x,y)
(1./5. + temp * &
(1./7. ))))
else
- ymerc = - log( 1. - efact ) /2. /strcmp(1)
+ ymerc = - log( 1. - efact ) *0.5 /strcmp(1)
endif
if (ymerc .gt. 6.) then
if (strcmp(1) .gt. almst1) then
@@ -398,7 +398,7 @@ subroutine cnxyll (strcmp, xi,eta, xlat,xlong)
(1./7. ))))
else
! Code for moderate values of gamma
- ymerc = - log ( 1. - arg1 ) /2. / gamma
+ ymerc = - log ( 1. - arg1 ) *0.5 / gamma
endif
! Convert ymerc to latitude
temp = exp( - abs(ymerc) )
diff --git a/src/drydepo_mod.f90 b/src/drydepo_mod.f90
index 9ce54e021c134677a55c7032f41f5ec47943472d..0c70351ee7bdbc982c971a3f9fa3e699905f1174 100644
--- a/src/drydepo_mod.f90
+++ b/src/drydepo_mod.f90
@@ -906,7 +906,7 @@ subroutine getvdep(n,ix,jy,ust,temp,pa,L,gr,rh,rr,snow,vdepo)
ylat=jy*dy+ylat0
if (ylat.lt.0) then
- jul=jul+365/2
+ jul=jul+365.*0.5
endif
@@ -1093,7 +1093,7 @@ subroutine getvdep_nest(n,ix,jy,ust,temp,pa, &
ylat=jy*dy+ylat0
if (ylat.lt.0) then
- jul=jul+365/2
+ jul=jul+365.*0.5
endif
@@ -1373,8 +1373,8 @@ subroutine partdep(nc,density,fract,schmi,vset,ra,ustar,nyl,rhoa,vdep_tmp)
alpha1=0.45+10.0/(exp(2.5*log10(dfdr))+30.0)
beta1=1.-37.0/(exp(3.0*log10(dfdr))+100.0)
kn1=10.**(alpha1*(-log10(Fn(ic)))**beta1)
- ks=(ks1(ic)+ks2(ic))/2.
- kn=(kn1+kn2(ic))/2.
+ ks=(ks1(ic)+ks2(ic))*0.5
+ kn=(kn1+kn2(ic))*0.5
endif
do i=1,20
diff --git a/src/flux_mod.f90 b/src/flux_mod.f90
index f2c0af93ef445bd8b17f1e83985882291c65e5f8..acaad622f7241f59419fcc6cad5315c265fabec2 100644
--- a/src/flux_mod.f90
+++ b/src/flux_mod.f90
@@ -78,8 +78,8 @@ subroutine calcfluxes(itime,nage,jpart,xold,yold,zold,thread)
#ifdef ETA
call update_zeta_to_z(itime,jpart)
#endif
- xmean=(xold+real(part(jpart)%xlon))/2.
- ymean=(yold+real(part(jpart)%ylat))/2.
+ xmean=(xold+real(part(jpart)%xlon))*0.5
+ ymean=(yold+real(part(jpart)%ylat))*0.5
ixave=int((xmean*dx+xoutshift)/dxout)
jyave=int((ymean*dy+youtshift)/dyout)
@@ -138,7 +138,7 @@ subroutine calcfluxes(itime,nage,jpart,xold,yold,zold,thread)
! 1) Particle does not cross domain boundary
- if (abs(xold-part(jpart)%xlon).lt.real(nx)/2.) then
+ if (abs(xold-part(jpart)%xlon).lt.real(nx)*0.5) then
ix1=int((xold*dx+xoutshift)/dxout+0.5)
ix2=int((part(jpart)%xlon*dx+xoutshift)/dxout+0.5)
do k=1,nspec
diff --git a/src/getfields_mod.f90 b/src/getfields_mod.f90
index fe60124a5437d8428674edbb86261300ebfe0ce0..2e256233c3acceb73c971be7b170691f5f3d0063 100644
--- a/src/getfields_mod.f90
+++ b/src/getfields_mod.f90
@@ -1440,8 +1440,8 @@ real function obukhov(ps,tsfc,tdsfc,tlev,ustar,hf,akm,bkm,plev)
tv=tsfc*(1.+0.378*e/ps) ! virtual temperature
rhoa=ps/(r_air*tv) ! air density
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
- ak1=(akm(1)+akm(2))/2.
- bk1=(bkm(1)+bkm(2))/2.
+ ak1=(akm(1)+akm(2))*0.5
+ bk1=(bkm(1)+bkm(2))*0.5
plev=ak1+bk1*ps ! Pressure level 1
end if
theta=tlev*(100000./plev)**(r_air/cpa) ! potential temperature
diff --git a/src/initialise_mod.f90 b/src/initialise_mod.f90
index 798f2acf07830c72462c25207347950759eb19c1..e2a8fb79c11782d996aadbe486fa7d579cea8216 100644
--- a/src/initialise_mod.f90
+++ b/src/initialise_mod.f90
@@ -154,7 +154,7 @@ subroutine releaseparticles(itime)
! Determine the local day and time
!*********************************
- xlonav=xlon0+(xpoint2(i)+xpoint1(i))/2.*dx ! longitude needed to determine local time
+ xlonav=xlon0+(xpoint2(i)+xpoint1(i))*0.5*dx ! longitude needed to determine local time
if (xlonav.lt.-180.) xlonav=xlonav+360.
if (xlonav.gt.180.) xlonav=xlonav-360.
jullocal=jul+real(xlonav,kind=dp)/360._dp ! correct approximately for time zone to obtain local time
@@ -194,7 +194,7 @@ subroutine releaseparticles(itime)
rfraction=abs(real(npart(i))*real(lsynctime)/ &
real(ireleaseend(i)-ireleasestart(i)))
if ((itime.eq.ireleasestart(i)).or. &
- (itime.eq.ireleaseend(i))) rfraction=rfraction/2.
+ (itime.eq.ireleaseend(i))) rfraction=rfraction*0.5
! Take the species-average time correction factor in order to scale the
! number of particles released this time
@@ -369,7 +369,7 @@ subroutine kindz_to_z(ipart)
if (press.lt.presspart) then
if (kz.eq.1) then
- call set_z(ipart,height(1)/2.)
+ call set_z(ipart,height(1)*0.5)
else
dp1=pressold-presspart
dp2=presspart-press
@@ -1078,7 +1078,7 @@ subroutine init_domainfill
deltacol=(pp(1)-pp(nz))/real(ncolumn)
- pnew=pp(1)+deltacol/2.
+ pnew=pp(1)+deltacol*0.5
jj=0
do j=1,ncolumn ! looping over the number of particles within the column
@@ -1253,7 +1253,7 @@ subroutine init_domainfill
write(*,*) 'Total number of particles at model start: ',numpart
write(*,*) 'Maximum number of particles per column: ',numcolumn
write(*,*) 'If ',fractus,' <1, better use more particles'
- fractus=sqrt(max(fractus,1.))/2.
+ fractus=sqrt(max(fractus,1.))*0.5
do ljy=ny_sn(1),ny_sn(2) ! loop about latitudes
do lix=nx_we(1),nx_we(2) ! loop about longitudes
@@ -1285,7 +1285,7 @@ subroutine init_domainfill
!*******************************************
deltacol=(pp(1)-pp(nz))/real(ncolumn)
- pnew=pp(1)+deltacol/2.
+ pnew=pp(1)+deltacol*0.5
do j=1,ncolumn
pnew=pnew-deltacol
do kz=1,nz-1
@@ -1472,13 +1472,13 @@ subroutine boundcond_domainfill(itime,loutend)
!*****************************************************************************
if (j.eq.1) then
- deltaz=(zcolumn_we(k,jy,2)+zcolumn_we(k,jy,1))/2.
+ deltaz=(zcolumn_we(k,jy,2)+zcolumn_we(k,jy,1))*0.5
else if (j.eq.numcolumn_we(k,jy)) then
! In order to avoid taking a very high column for very many particles,
! use the deltaz from one particle below instead
- deltaz=(zcolumn_we(k,jy,j)-zcolumn_we(k,jy,j-2))/2.
+ deltaz=(zcolumn_we(k,jy,j)-zcolumn_we(k,jy,j-2))*0.5
else
- deltaz=(zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1))/2.
+ deltaz=(zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1))*0.5
endif
if ((jy.eq.ny_sn(1)).or.(jy.eq.ny_sn(2))) then
boundarea=deltaz*111198.5/2.*dy
@@ -1557,8 +1557,8 @@ subroutine boundcond_domainfill(itime,loutend)
! reduced by the mass of this (these) particle(s)
!******************************************************************************
- if (acc_mass_we(k,jy,j).ge.xmassperparticle/2.) then
- mmass=int((acc_mass_we(k,jy,j)+xmassperparticle/2.)/ &
+ if (acc_mass_we(k,jy,j).ge.xmassperparticle*0.5) then
+ mmass=int((acc_mass_we(k,jy,j)+xmassperparticle*0.5)/ &
xmassperparticle)
acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)- &
real(mmass)*xmassperparticle
@@ -1687,18 +1687,18 @@ subroutine boundcond_domainfill(itime,loutend)
!*****************************************************************************
if (j.eq.1) then
- deltaz=(zcolumn_sn(k,ix,2)+zcolumn_sn(k,ix,1))/2.
+ deltaz=(zcolumn_sn(k,ix,2)+zcolumn_sn(k,ix,1))*0.5
else if (j.eq.numcolumn_sn(k,ix)) then
! deltaz=height(nz)-(zcolumn_sn(k,ix,j-1)+
! + zcolumn_sn(k,ix,j))/2.
! In order to avoid taking a very high column for very many particles,
! use the deltaz from one particle below instead
- deltaz=(zcolumn_sn(k,ix,j)-zcolumn_sn(k,ix,j-2))/2.
+ deltaz=(zcolumn_sn(k,ix,j)-zcolumn_sn(k,ix,j-2))*0.5
else
- deltaz=(zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1))/2.
+ deltaz=(zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1))*0.5
endif
if ((ix.eq.nx_we(1)).or.(ix.eq.nx_we(2))) then
- boundarea=deltaz*111198.5/2.*cosfact*dx
+ boundarea=deltaz*111198.5*0.5*cosfact*dx
else
boundarea=deltaz*111198.5*cosfact*dx
endif
@@ -1773,8 +1773,8 @@ subroutine boundcond_domainfill(itime,loutend)
! reduced by the mass of this (these) particle(s)
!******************************************************************************
- if (acc_mass_sn(k,ix,j).ge.xmassperparticle/2.) then
- mmass=int((acc_mass_sn(k,ix,j)+xmassperparticle/2.)/ &
+ if (acc_mass_sn(k,ix,j).ge.xmassperparticle*0.5) then
+ mmass=int((acc_mass_sn(k,ix,j)+xmassperparticle*0.5)/ &
xmassperparticle)
acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)- &
real(mmass)*xmassperparticle
diff --git a/src/netcdf_output_mod.f90 b/src/netcdf_output_mod.f90
index e464d0e5383f5ec7ec64dcda1b925118e2007a87..97b465fc53a3c3b4755069415c0627c73f586898 100644
--- a/src/netcdf_output_mod.f90
+++ b/src/netcdf_output_mod.f90
@@ -964,9 +964,9 @@ subroutine concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridto
!$OMP DO
do kz=1,numzgrid
if (kz.eq.1) then
- halfheight=outheight(1)/2.
+ halfheight=outheight(1)*0.5
else
- halfheight=(outheight(kz)+outheight(kz-1))/2.
+ halfheight=(outheight(kz)+outheight(kz-1))*0.5
endif
do kzz=2,nz
if ((height(kzz-1).lt.halfheight).and. &
@@ -1351,9 +1351,9 @@ subroutine concoutput_nest_netcdf(itime,outnum)
!$OMP DO
do kz=1,numzgrid
if (kz.eq.1) then
- halfheight=outheight(1)/2.
+ halfheight=outheight(1)*0.5
else
- halfheight=(outheight(kz)+outheight(kz-1))/2.
+ halfheight=(outheight(kz)+outheight(kz-1))*0.5
endif
do kzz=2,nz
if ((height(kzz-1).lt.halfheight).and. &
diff --git a/src/pbl_profile_mod.f90 b/src/pbl_profile_mod.f90
index 4f4584aeae7a2af5b86b432b38e1e1d29a3da140..bc30358256621e4419bb0fd17bb9a873d6fab19b 100644
--- a/src/pbl_profile_mod.f90
+++ b/src/pbl_profile_mod.f90
@@ -54,7 +54,7 @@ function psih (z,l)
- b*c/d + 1.
else
x=(1.-16.*zeta)**(.25)
- psih=2.*log((1.+x*x)/2.)
+ psih=2.*log((1.+x*x)*0.5)
end if
end if
@@ -78,9 +78,9 @@ real function psim(z,al)
if(zeta.le.0.) then
! UNSTABLE CASE
x=(1.-15.*zeta)**0.25
- a1=((1.+x)/2.)**2
- a2=(1.+x**2)/2.
- psim=log(a1*a2)-2.*atan(x)+pi/2.
+ a1=((1.+x)*0.5)**2
+ a2=(1.+x**2)*0.5
+ psim=log(a1*a2)-2.*atan(x)+pi*0.5
else
! STABLE CASE
psim=-4.7*zeta
@@ -171,9 +171,9 @@ subroutine pbl_profile(ps,td2m,zml1,t2m,tml1,u10m,uml1,stress,hf)
!! Successive approximation will
al=50. !! not converge
ustar=(vonkarman*deltau)/ &
- (log(zml1/10.)-psim(zml1,al)+psim(10.,al))
+ (log(zml1*0.1)-psim(zml1,al)+psim(10.,al))
thetastar=(vonkarman*deltat/r1)/ &
- (log(zml1/2.)-psih(zml1,al)+psih(2.,al))
+ (log(zml1*0.5)-psih(zml1,al)+psih(2.,al))
hf=rhoa*cpa*ustar*thetastar
stress=ustar*ustar*rhoa
return
@@ -183,9 +183,9 @@ subroutine pbl_profile(ps,td2m,zml1,t2m,tml1,u10m,uml1,stress,hf)
do iter=1,maxiter
alold=al
ustar=(vonkarman*deltau)/ &
- (log(zml1/10.)-psim(zml1,al)+psim(10.,al))
+ (log(zml1*0.1)-psim(zml1,al)+psim(10.,al))
thetastar=(vonkarman*deltat/r1)/ &
- (log(zml1/2.)-psih(zml1,al)+psih(2.,al))
+ (log(zml1*0.5)-psih(zml1,al)+psih(2.,al))
al=(tmean*ustar**2)/(ga*vonkarman*thetastar)
aldiff=abs((al-alold)/alold)
if(aldiff.lt.0.01) exit !! Successive approximation successful
diff --git a/src/readoptions_mod.f90 b/src/readoptions_mod.f90
index fbc1f90f313188be17a1f7413e9b5230d3c26ca1..eddb689710d0b4312e901140518a647fd25119a8 100644
--- a/src/readoptions_mod.f90
+++ b/src/readoptions_mod.f90
@@ -1617,9 +1617,9 @@ subroutine readoutgrid
! Determine the half levels, i.e. middle levels of the output grid
!*****************************************************************
- outheighthalf(1)=outheight(1)/2.
+ outheighthalf(1)=outheight(1)*0.5
do j=2,numzgrid
- outheighthalf(j)=(outheight(j-1)+outheight(j))/2.
+ outheighthalf(j)=(outheight(j-1)+outheight(j))*0.5
end do
xoutshift=xlon0-outlon0
@@ -2475,7 +2475,7 @@ subroutine readreleases
jul1=juldate(id1,it1)
jul2=juldate(id2,it2)
- julm=(jul1+jul2)/2.
+ julm=(jul1+jul2)*0.5
if (jul1.gt.jul2) then
write(*,*) 'FLEXPART MODEL ERROR'
write(*,*) 'Release stops before it begins.'
@@ -2975,7 +2975,7 @@ subroutine readspecies(id_spec,pos_spec)
Fs(pos_spec)=f*e**(1.3)*(dquer(pos_spec)**3/(psa*pia*pla)) ! Stokes' regime
endif
! Pre-compute ks and kn values needed for horizontal and average orientation (B&B Figure 12 k_(s,max))
- ks1(pos_spec)=(Fs(pos_spec)**(1./3.) + Fs(pos_spec)**(-1./3.))/2.
+ ks1(pos_spec)=(Fs(pos_spec)**(1./3.) + Fs(pos_spec)**(-1./3.))*0.5
ks2(pos_spec)=0.5*((Fs(pos_spec)**0.05)+(Fs(pos_spec)**(-0.36)))
kn2(pos_spec)=10.**(alpha2*(-log10(Fn(pos_spec)))**beta2)
diff --git a/src/settling_mod.f90 b/src/settling_mod.f90
index 3496dac46da62ab48fc92e91c21cdbe3f034c2cd..1b8cbafd0e2367e43648c70486e31a0dfaec1af6 100644
--- a/src/settling_mod.f90
+++ b/src/settling_mod.f90
@@ -268,8 +268,8 @@ subroutine get_settling(xt,yt,zt,nsp,settling)
alpha1=0.45+10.0/(exp(2.5*log10(dfdr))+30.0)
beta1=1.-37.0/(exp(3.0*log10(dfdr))+100.0)
kn1=10.**(alpha1*(-log10(Fn(nsp)))**beta1)
- ks=(ks1(nsp)+ks2(nsp))/2.
- kn=(kn1+kn2(nsp))/2.
+ ks=(ks1(nsp)+ks2(nsp))*0.5
+ kn=(kn1+kn2(nsp))*0.5
endif
do i=1,20
diff --git a/src/verttransform_mod.f90 b/src/verttransform_mod.f90
index 9810af594ae1a630872ebe7559d91659a83edbc6..9952cf746bd9546eb8f0fc6d43875456069e4f02 100644
--- a/src/verttransform_mod.f90
+++ b/src/verttransform_mod.f90
@@ -642,8 +642,8 @@ subroutine verttransform_ecmwf_windfields(n,nxlim,nylim,uuh,vvh,wwh,pvh,rhoh,prs
dz2=etauvheight(ix,jy,kz,n)-height(iz)
dz=dz1+dz2
- dzdx1=(etauvheight(ixp,jy,kz-1,n)-etauvheight(ix1,jy,kz-1,n))/2.
- dzdx2=(etauvheight(ixp,jy,kz,n)-etauvheight(ix1,jy,kz,n))/2.
+ dzdx1=(etauvheight(ixp,jy,kz-1,n)-etauvheight(ix1,jy,kz-1,n))*0.5
+ dzdx2=(etauvheight(ixp,jy,kz,n)-etauvheight(ix1,jy,kz,n))*0.5
dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
@@ -724,7 +724,7 @@ subroutine verttransform_ecmwf_polar(n)
ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
vv(nx/2-1,nymin1,iz,n))-xlonr
else
- ddpol=pi/2.-xlonr
+ ddpol=pi*0.5-xlonr
endif
if(ddpol.lt.0.) ddpol=2.*pi+ddpol
if(ddpol.gt.2.*pi) ddpol=ddpol-2.*pi
@@ -766,7 +766,7 @@ subroutine verttransform_ecmwf_polar(n)
ddpol=pi+atan(uueta(nx/2-1,nymin1,iz,n)/ &
vveta(nx/2-1,nymin1,iz,n))-xlonr
else
- ddpol=pi/2-xlonr
+ ddpol=pi*0.5-xlonr
endif
if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
@@ -879,7 +879,7 @@ subroutine verttransform_ecmwf_polar(n)
ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
vv(nx/2-1,0,iz,n))+xlonr
else
- ddpol=pi/2.-xlonr
+ ddpol=pi*0.5-xlonr
endif
if(ddpol.lt.0.) ddpol=2.*pi+ddpol
if(ddpol.gt.2.*pi) ddpol=ddpol-2.*pi
@@ -924,7 +924,7 @@ subroutine verttransform_ecmwf_polar(n)
ddpol=pi+atan(uueta(nx/2-1,0,iz,n)/ &
vveta(nx/2-1,0,iz,n))+xlonr
else
- ddpol=pi/2-xlonr
+ ddpol=pi*0.5-xlonr
endif
if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
@@ -1686,7 +1686,7 @@ subroutine verttransform_gfs(n,uuh,vvh,wwh,pvh)
elseif (vv(nx/2-1,nymin1,iz,n).gt.0.) then
ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
else
- ddpol=pi/2.-xlonr
+ ddpol=pi*0.5-xlonr
endif
if(ddpol.lt.0.) ddpol=2.*pi+ddpol
if(ddpol.gt.2.*pi) ddpol=ddpol-2.*pi
@@ -1754,7 +1754,7 @@ subroutine verttransform_gfs(n,uuh,vvh,wwh,pvh)
elseif (vv(nx/2-1,0,iz,n).gt.0.) then
ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))-xlonr
else
- ddpol=pi/2.-xlonr
+ ddpol=pi*0.5-xlonr
endif
if(ddpol.lt.0.) ddpol=2.*pi+ddpol
if(ddpol.gt.2.*pi) ddpol=ddpol-2.*pi
@@ -2023,7 +2023,7 @@ subroutine verttransform_ecmwf_heights(nxlim,nylim, &
end do
do kz=2,nwz-1
- wzlev(ix,jy,kz)=(uvzlev(ix,jy,kz+1)+uvzlev(ix,jy,kz))/2.
+ wzlev(ix,jy,kz)=(uvzlev(ix,jy,kz+1)+uvzlev(ix,jy,kz))*0.5
end do
wzlev(ix,jy,nwz)=wzlev(ix,jy,nwz-1)+ &
uvzlev(ix,jy,nuvz)-uvzlev(ix,jy,nuvz-1)