diff --git a/src/FLEXPART.f90 b/src/FLEXPART.f90
index 78fbbf35a2726813a1d54b7a69c3831330de6427..70f2b155c6ba77ca348e7aa9fc569aed2ecd7ffb 100644
--- a/src/FLEXPART.f90
+++ b/src/FLEXPART.f90
@@ -248,7 +248,7 @@ subroutine read_options_and_initialise_flexpart
! Set the upper level for where the convection will be working
!*************************************************************
- call set_upperlevel_convect
+ call set_conv_top
if (numbnests.ge.1) then
! If nested wind fields are used, allocate arrays
@@ -295,7 +295,7 @@ subroutine read_options_and_initialise_flexpart
! Allocate dry deposition fields if necessary
!*********************************************
call alloc_drydepo
- call alloc_convection
+ call alloc_convect
call alloc_getfields
call alloc_interpol
diff --git a/src/advance_mod.f90 b/src/advance_mod.f90
index 34e1b3e42af7d6dbed8e334d067e4ee1e47e3291..2cb3a8e9e2c00df906967320c6af04b21d9c221a 100644
--- a/src/advance_mod.f90
+++ b/src/advance_mod.f90
@@ -27,11 +27,11 @@ module advance_mod
eps_eta=1.e-4
real :: &
eps
- private :: advance_abovePBL,advance_pbl,advance_petterson_corr,&
- advance_updateXY,advance_adjusttopheight
+ private :: adv_above_pbl, adv_in_pbl, petterson_corr, update_xy, pushpartdown
contains
subroutine advance(itime,ipart,thread)
+
!*****************************************************************************
! *
! Calculation of turbulent particle trajectories utilizing a *
@@ -128,9 +128,10 @@ subroutine advance(itime,ipart,thread)
weta_settling, & ! Settling velocity in eta coordinates
tropop, & ! height of troposphere
dxsave,dysave, & ! accumulated displacement in long and lat
- dawsave,dcwsave ! accumulated displacement in wind directions
+ dawsave,dcwsave ! accumulated displacement in wind directions
logical :: &
- abovePBL ! flag that will be set to 'true' if computation needs to be completed above PBL
+ abovePBL
+ ! flag will be set to 'true' if computation needs to be completed above PBL
eps=nxmax/3.e5
@@ -181,13 +182,14 @@ subroutine advance(itime,ipart,thread)
! Determine the lower left corner and its distance to the current position
! Calculate variables for time interpolation
!*******************************************
- call initialise_interpol_mod(itime,real(part(ipart)%xlon),real(part(ipart)%ylat),&
- real(part(ipart)%z),real(part(ipart)%zeta))
+ call init_interpol(itime, &
+ real(part(ipart)%xlon),real(part(ipart)%ylat),&
+ real(part(ipart)%z), real(part(ipart)%zeta))
! Compute maximum mixing height around particle position
!*******************************************************
- ! Compute the height of the troposphere and the PBL at the x-y location of the particle
+ ! Compute height of troposphere and PBL at x-y location of particle
call interpol_htropo_hmix(tropop,h)
zeta=real(part(ipart)%z)/h
@@ -195,13 +197,15 @@ subroutine advance(itime,ipart,thread)
! If particle is in the PBL, interpolate once and then make a
! time loop until end of interval is reached
!*************************************************************
- ! In the PBL we use meters instead of eta coordinates for the vertical transport
+ ! In the PBL we use meters instead of eta coordinates for vertical transport
+
abovePBL=.true.
if (zeta.le.1.) then
+
abovePBL=.false.
- call advance_pbl(itime,itimec,&
+ call adv_in_pbl(itime,itimec,&
dxsave,dysave,dawsave,dcwsave,abovePBL,nrand,ipart,thread)
- if ((wind_coord_type.eq.'ETA').and.(lsettling)) then
+ if (wind_coord_type.eq.'ETA' .and. lsettling) then
call w_to_weta(itime,real(part(ipart)%idt),part(ipart)%xlon, &
part(ipart)%ylat,part(ipart)%z,part(ipart)%zeta, &
part(ipart)%settling,weta_settling)
@@ -218,10 +222,10 @@ subroutine advance(itime,ipart,thread)
! Interpolate the wind
!*********************
- if (abovePBL) then
- call advance_abovePBL(itime,itimec,&
- dxsave,dysave,ux,vy,tropop,nrand,ipart)
- endif ! Above PBL computation
+
+ if (abovePBL) call adv_above_pbl(itime,itimec,dxsave,dysave, &
+ ux,vy,tropop,nrand,ipart)
+ ! Above PBL computation
!****************************************************************
! Add mesoscale random disturbances
@@ -229,18 +233,21 @@ subroutine advance(itime,ipart,thread)
! computation time
!****************************************************************
-
! Mesoscale wind velocity fluctuations are obtained by scaling
! with the standard deviation of the grid-scale winds surrounding
! the particle location, multiplied by a factor fturbmeso.
! The autocorrelation time constant is taken as half the
! time interval between wind fields
!****************************************************************
- if (.not. turboff) then ! mesoscale turbulence is found to give issues, so turned off
+
+ if (.not. turboff) then
+ ! mesoscale turbulence is found to give issues, so turned off
if (lmesoscale_turb) then
- call interpol_mesoscale(itime,real(part(ipart)%xlon),real(part(ipart)%ylat), &
+ call interpol_mesoscale(itime, &
+ real(part(ipart)%xlon),real(part(ipart)%ylat), &
real(part(ipart)%z),real(part(ipart)%zeta))
- call turb_mesoscale(nrand,dxsave,dysave,ipart,usig,vsig,wsig,wsigeta,eps_eta)
+ call turb_mesoscale(nrand,dxsave,dysave,ipart, &
+ usig,vsig,wsig,wsigeta,eps_eta)
endif
!*************************************************************
@@ -250,16 +257,17 @@ subroutine advance(itime,ipart,thread)
!*************************************************************
call windalign(dxsave,dysave,dawsave,dcwsave,ux,vy)
- dxsave=dxsave+ux ! comment by mc: comment this line to stop the particles horizontally for test reasons
+ dxsave=dxsave+ux
+ ! comment by MC: comment this line to stop particles horizontally for tests
dysave=dysave+vy
endif
- call advance_updateXY(dxsave,dysave,ipart)
- if (part(ipart)%nstop.eqv..true.) return
+ call update_xy(dxsave,dysave,ipart)
+ if (part(ipart)%nstop) return
! If particle above highest model level, set it back into the domain
!*******************************************************************
- call advance_adjusttopheight(ipart)
+ call pushpartdown(ipart)
!************************************************************************
! Now we could finish, as this was done in FLEXPART versions up to 4.0.
@@ -274,12 +282,12 @@ subroutine advance(itime,ipart,thread)
! is the "old" wind as required by the scheme); otherwise do nothing
!*************************************************************************
- if (part(ipart)%idt.ne.abs(lsynctime)) return
+ if (part(ipart)%idt .ne. abs(lsynctime)) return
- ! The Petterssen scheme can only be applied if the ending time of the time step
- ! (itime+ldt*ldirect) is still between the two wind fields held in memory;
- ! otherwise do nothing
- !******************************************************************************
+ ! The Petterssen scheme can only be applied if the ending time of the time step
+ ! (itime+ldt*ldirect) is still between the two wind fields held in memory;
+ ! otherwise do nothing
+ !******************************************************************************
if (abs(itime+part(ipart)%idt*ldirect).gt.abs(memtime(2))) return
@@ -289,16 +297,18 @@ subroutine advance(itime,ipart,thread)
! ngr = ngrid
! call find_ngrid(part(ipart)%xlon,part(ipart)%ylat)
- if (nglobal.and.(real(part(ipart)%ylat).gt.switchnorthg)) then
+ if (nglobal .and. real(part(ipart)%ylat).gt.switchnorthg) then
ngr=-1
- else if (sglobal.and.(real(part(ipart)%ylat).lt.switchsouthg)) then
+ else if (sglobal.and. real(part(ipart)%ylat).lt.switchsouthg) then
ngr=-2
else
ngr=0
! Temporary fix for nested layer edges: replaced eps with dxn and dyn (LB)
do j=numbnests,1,-1
- if ((real(part(ipart)%xlon).gt.xln(j)+dxn(j)).and.(real(part(ipart)%xlon).lt.xrn(j)-dxn(j)).and. &
- (real(part(ipart)%ylat).gt.yln(j)+dyn(j)).and.(real(part(ipart)%ylat).lt.yrn(j)-dyn(j))) then
+ if (real(part(ipart)%xlon).gt.xln(j)+dxn(j) .and. &
+ real(part(ipart)%xlon).lt.xrn(j)-dxn(j) .and. &
+ real(part(ipart)%ylat).gt.yln(j)+dyn(j) .and. &
+ real(part(ipart)%ylat).lt.yrn(j)-dyn(j)) then
ngr=j
exit
endif
@@ -307,11 +317,11 @@ subroutine advance(itime,ipart,thread)
if (ngr.ne.ngrid) return
- call advance_petterson_corr(itime,ipart)
+ call petterson_corr(itime,ipart)
+
end subroutine advance
-subroutine advance_abovePBL(itime,itimec,dxsave,dysave,&
- ux,vy,tropop,nrand,ipart)
+subroutine adv_above_pbl(itime,itimec,dxsave,dysave,ux,vy,tropop,nrand,ipart)
implicit none
integer, intent(in) :: &
@@ -399,40 +409,51 @@ subroutine advance_abovePBL(itime,itimec,dxsave,dysave,&
dysave=dysave+(v+vy)*dt
select case (wind_coord_type)
+
case ('ETA')
+
if (wp.ne.0.) then
call update_zeta_to_z(itime,ipart)
call update_z(ipart,wp*dt*real(ldirect))
- if (part(ipart)%z.lt.0.) call set_z(ipart,min(h-eps2,-1.*part(ipart)%z)) ! if particle below ground -> reflection
+ if (part(ipart)%z.lt.0.) call set_z(ipart,min(h-eps2,-1.*part(ipart)%z))
+ ! if particle below ground -> reflection
call update_z_to_zeta(itime,ipart)
endif
call update_zeta(ipart,weta*dt*real(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)
+
case ('METER')
+
call update_z(ipart,(w+wp)*dt*real(ldirect))
if (part(ipart)%z.lt.0.) call set_z(ipart,min(h-eps2,-1.*part(ipart)%z))
+
case default
+
call update_z(ipart,(w+wp)*dt*real(ldirect))
if (part(ipart)%z.lt.0.) call set_z(ipart,min(h-eps2,-1.*part(ipart)%z))
+
end select
-end subroutine advance_abovePBL
-subroutine advance_pbl(itime,itimec,&
- dxsave,dysave,dawsave,dcwsave,abovePBL,nrand,ipart,thread)
+end subroutine adv_above_pbl
+
+subroutine adv_in_pbl(itime,itimec, dxsave,dysave,dawsave,dcwsave, abovePBL, &
+ nrand,ipart,thread)
+
use drydepo_mod, only: drydepo_probability
implicit none
logical, intent(inout) :: &
- abovePBL ! flag that will be set to 'true' if computation needs to be completed above PBL
+ abovePBL
+ ! flag will be set to 'true' if computation needs to be completed above PBL
integer, intent(in) :: &
itime, & ! time index
ipart, & ! particle index
thread ! number of the omp thread
real, intent(inout) :: &
dxsave,dysave, & ! accumulated displacement in long and lat
- dawsave,dcwsave ! accumulated displacement in wind directions
+ dawsave,dcwsave ! accumulated displacement in wind directions
integer, intent(inout) :: &
itimec, & ! next timestep
nrand ! random number used for turbulence
@@ -440,7 +461,7 @@ subroutine advance_pbl(itime,itimec,&
dt, & ! real(ldt)
xts,yts,zts,ztseta, & ! local 'real' copy of the particle position
rhoa, & ! air density, used in CBL
- rhograd ! vertical gradient of the air density, used in CBL
+ rhograd ! vertical gradient of air density, used in CBL
integer :: &
loop, & ! loop variable for time in the PBL
nsp,insp ! loop variable for species
@@ -452,13 +473,15 @@ subroutine advance_pbl(itime,itimec,&
! BEGIN TIME LOOP
!================
! For wind_coord_type=ETA:
- ! Within this loop, only METER coordinates are used, and the new z value will be updated
- ! to ETA coordinates at the end
- !***************************************************************************************
+ ! Within this loop, only METER coordinates are used, and the new z value will
+ ! be updated to ETA coordinates at the end
+ !****************************************************************************
+
call update_zeta_to_z(itime,ipart)
loop=0
- pbl_loop : do
+ pbl_loop: do
+
loop=loop+1
if (method.eq.1) then
part(ipart)%idt=min(part(ipart)%idt,abs(lsynctime-itimec+itime))
@@ -473,7 +496,8 @@ subroutine advance_pbl(itime,itimec,&
zts=real(part(ipart)%z)
zeta=zts/h
- if (loop.eq.1) then ! Temporal interpolation only done for the first iteration
+ if (loop.eq.1) then ! Temporal interpolation only for the first iteration
+
if (ngrid.le.0) then
xts=real(part(ipart)%xlon)
yts=real(part(ipart)%ylat)
@@ -483,6 +507,7 @@ subroutine advance_pbl(itime,itimec,&
endif
else
+
! Determine the level below the current position for u,v,rho
!***********************************************************
call find_z_level_meters(zts)
@@ -491,8 +516,8 @@ subroutine advance_pbl(itime,itimec,&
! calculate it
!*****************************************************
call interpol_pbl_misslev()
- endif
+ endif
! Vertical interpolation of u,v,w,rho and drhodz
!***********************************************
@@ -500,25 +525,35 @@ subroutine advance_pbl(itime,itimec,&
! Vertical distance to the level below and above current position
! both in terms of (u,v) and (w) fields
!****************************************************************
+
call interpol_pbl_short(zts,rhoa,rhograd) ! Vertical interpolation
! Compute the turbulent disturbances
! Determine the sigmas and the timescales
!****************************************
+
if (.not.turboff) then
- call turb_pbl(ipart,nrand,dt,zts,rhoa,rhograd,thread) ! Note: zts and nrand get updated
+ call turb_pbl(ipart,nrand,dt,zts,rhoa,rhograd,thread)
+ ! Note: zts and nrand get updated
+
! Determine time step for next integration
!*****************************************
if (turbswitch) then
- part(ipart)%idt=int(min(tlw,h/max(2.*abs(part(ipart)%turbvel%w*sigw),1.e-5), &
- 0.5/abs(dsigwdz))*ctl)
+ part(ipart)%idt = int( &
+ min( tlw, &
+ h/max( 2.*abs(part(ipart)%turbvel%w*sigw), 1.e-5 ), &
+ 0.5/abs(dsigwdz) &
+ ) *ctl)
else
- part(ipart)%idt=int(min(tlw,h/max(2.*abs(part(ipart)%turbvel%w),1.e-5))*ctl)
+ part(ipart)%idt = int( &
+ min( tlw, &
+ h/max( 2.*abs(part(ipart)%turbvel%w), 1.e-5) &
+ ) *ctl)
endif
else
- part(ipart)%turbvel%u=0.0
- part(ipart)%turbvel%v=0.0
- part(ipart)%turbvel%w=0.0
+ part(ipart)%turbvel%u=0.
+ part(ipart)%turbvel%v=0.
+ part(ipart)%turbvel%w=0.
endif
part(ipart)%idt=max(part(ipart)%idt,mintime)
@@ -557,22 +592,31 @@ subroutine advance_pbl(itime,itimec,&
dawsave=dawsave+part(ipart)%turbvel%u*dt
dcwsave=dcwsave+part(ipart)%turbvel%v*dt
! How can I change the w to w(eta) efficiently?
+
select case (wind_coord_type)
+
case ('ETA')
+
call update_z(ipart,w*dt*real(ldirect))
zts=real(part(ipart)%z)
- ! HSO/AL: Particle managed to go over highest level -> interpolation error in goto 700
+ ! HSO/AL: Particle managed to go over highest level -> interpolation
+ ! error in goto 700
! alias interpol_wind (division by zero)
- if (zts.ge.height(nz)) call set_z(ipart,height(nz)-100.*eps) ! Manually for z instead
+ if (zts.ge.height(nz)) call set_z(ipart,height(nz)-100.*eps)
+ ! Manually for z instead
+
case ('METER')
+
call update_z(ipart,w*dt*real(ldirect))
- call advance_adjusttopheight(ipart)
+ call pushpartdown(ipart)
+
end select
zts=real(part(ipart)%z)
if (zts.gt.h) then
call update_z_to_zeta(itime,ipart)
- if (itimec.ne.itime+lsynctime) abovePBL=.true. ! complete the current interval above PBL
+ if (itimec.ne.itime+lsynctime) abovePBL=.true.
+ ! complete the current interval above PBL
return
endif
@@ -580,19 +624,22 @@ subroutine advance_pbl(itime,itimec,&
!************************************
call drydepo_probability(part(ipart)%prob,dt,zts,vdepo)
- if (zts.lt.0.) call set_z(ipart,min(h-eps2,-1.*part(ipart)%z)) ! if particle below ground -> reflection
-
+ if (zts.lt.0.) call set_z(ipart,min(h-eps2,-1.*part(ipart)%z))
+ ! if particle below ground -> reflection
if (itimec.eq.(itime+lsynctime)) then
- ! Converting the z position that changed through turbulence motions to eta coords
+ ! Convert z position that changed by turbulent motions to eta coords
call update_z_to_zeta(itime,ipart)
return ! finished
endif
+
end do pbl_loop
+
call update_z_to_zeta(itime,ipart)
-end subroutine advance_pbl
-subroutine advance_petterson_corr(itime,ipart)
+end subroutine adv_in_pbl
+
+subroutine petterson_corr(itime,ipart)
implicit none
@@ -648,22 +695,32 @@ subroutine advance_petterson_corr(itime,ipart)
select case (wind_coord_type)
case ('ETA')
+
call update_zeta_to_z(itime+part(ipart)%idt,ipart)
call update_z_to_zeta(itime+part(ipart)%idt,ipart)
zts=real(part(ipart)%z)
- call get_settling(itime+part(ipart)%idt,xts,yts,zts,nsp,part(ipart)%settling) !bugfix
- call w_to_weta(itime+part(ipart)%idt,real(part(ipart)%idt),part(ipart)%xlon, &
- part(ipart)%ylat,part(ipart)%z,part(ipart)%zeta, &
- part(ipart)%settling,weta_settling)
+ call get_settling( &
+ itime+part(ipart)%idt,xts,yts,zts,nsp,part(ipart)%settling) !bugfix
+ call w_to_weta( &
+ itime+part(ipart)%idt, real(part(ipart)%idt), part(ipart)%xlon, &
+ part(ipart)%ylat, part(ipart)%z, part(ipart)%zeta, &
+ part(ipart)%settling, weta_settling)
weta=weta+weta_settling
- !woldeta=real(part(ipart)%zeta-part(ipart)%zeta_prev)/real(part(ipart)%idt*ldirect)
+ !woldeta=
+ !real(part(ipart)%zeta-part(ipart)%zeta_prev)/real(part(ipart)%idt*ldirect)
+
case ('METER')
- call get_settling(itime+part(ipart)%idt,xts,yts,zts,nsp,part(ipart)%settling)
+
+ call get_settling( &
+ itime+part(ipart)%idt,xts,yts,zts,nsp,part(ipart)%settling)
w=w+part(ipart)%settling
case default
- call get_settling(itime+part(ipart)%idt,xts,yts,zts,nsp,part(ipart)%settling)
+
+ call get_settling( &
+ itime+part(ipart)%idt,xts,yts,zts,nsp,part(ipart)%settling)
w=w+part(ipart)%settling
+
end select
end if
endif
@@ -673,37 +730,45 @@ subroutine advance_petterson_corr(itime,ipart)
! (use half of it to correct position according to Petterssen)
!*************************************************************
- u=(u-uold)/2.
- v=(v-vold)/2.
+ u=(u-uold)*0.5
+ v=(v-vold)*0.5
select case (wind_coord_type)
+
case ('ETA')
+
weta=(weta-woldeta)/2.
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)
case ('METER')
+
w=(w-wold)/2.
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
+ if (part(ipart)%z.lt.0.) call set_z(ipart,min(h-eps2,-1.*part(ipart)%z)) ! if particle below ground -> reflection
case default
+
w=(w-wold)/2.
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))
+
end select
! Finally, correct the old position
!**********************************
- call advance_updateXY(u*part(ipart)%idt,v*part(ipart)%idt,ipart)
+
+ call update_xy(u*part(ipart)%idt,v*part(ipart)%idt,ipart)
! If particle above highest model level, set it back into the domain
!*******************************************************************
- call advance_adjusttopheight(ipart)
-end subroutine advance_petterson_corr
-subroutine advance_updateXY(xchange,ychange,ipart)
+ call pushpartdown(ipart)
+
+end subroutine petterson_corr
+
+subroutine update_xy(xchange,ychange,ipart)
implicit none
@@ -718,20 +783,30 @@ subroutine advance_updateXY(xchange,ychange,ipart)
eps=nxmax/3.e5
if (ngrid.ge.0) then
+
cosfact=dxconst/cos((real(part(ipart)%ylat)*dy+ylat0)*pi180)
call update_xlon(ipart,real(xchange*cosfact*real(ldirect),kind=dp))
call update_ylat(ipart,real(ychange*dyconst*real(ldirect),kind=dp))
+
else if (ngrid.eq.-1) then ! around north pole
- xlon=xlon0+real(part(ipart)%xlon)*dx !comment by mc: compute old particle position
+
+ xlon=xlon0+real(part(ipart)%xlon)*dx !comment by MC: compute old part pos.
ylat=ylat0+real(part(ipart)%ylat)*dy
- call cll2xy(northpolemap,ylat,xlon,xpol,ypol) !convert old particle position in polar stereographic
- gridsize=1000.*cgszll(northpolemap,ylat,xlon) !calculate size in m of grid element in polar stereographic coordinate
- xpol=xpol+xchange/gridsize*real(ldirect) !position in grid unit polar stereographic
+ call cll2xy(northpolemap,ylat,xlon,xpol,ypol)
+ !convert old particle position in polar stereographic
+ gridsize=1000.*cgszll(northpolemap,ylat,xlon)
+ !calculate size in m of grid element in polar stereographic coordinate
+ xpol=xpol+xchange/gridsize*real(ldirect)
+ !position in grid unit polar stereographic
ypol=ypol+ychange/gridsize*real(ldirect)
- call cxy2ll(northpolemap,xpol,ypol,ylat,xlon) !convert to lat long coordinate
- call set_xlon(ipart,real((xlon-xlon0)/dx,kind=dp))!convert to grid units in lat long coordinate, comment by mc
+ call cxy2ll(northpolemap,xpol,ypol,ylat,xlon)
+ !convert to lat long coordinate
+ call set_xlon(ipart,real((xlon-xlon0)/dx,kind=dp))
+ !convert to grid units in lat long coordinate, comment by mc
call set_ylat(ipart,real((ylat-ylat0)/dy,kind=dp))
+
else if (ngrid.eq.-2) then ! around south pole
+
xlon=xlon0+real(part(ipart)%xlon)*dx
ylat=ylat0+real(part(ipart)%ylat)*dy
call cll2xy(southpolemap,ylat,xlon,xpol,ypol)
@@ -741,44 +816,54 @@ subroutine advance_updateXY(xchange,ychange,ipart)
call cxy2ll(southpolemap,xpol,ypol,ylat,xlon)
call set_xlon(ipart,real((xlon-xlon0)/dx,kind=dp))
call set_ylat(ipart,real((ylat-ylat0)/dy,kind=dp))
+
endif
! If global data are available, use cyclic boundary condition
!************************************************************
+
if (xglobal) then
- if (part(ipart)%xlon.ge.real(nxmin1,kind=dp)) call update_xlon(ipart,-real(nxmin1,kind=dp))
- if (part(ipart)%xlon.lt.0.) call update_xlon(ipart,real(nxmin1,kind=dp))
- if (part(ipart)%xlon.le.real(eps,kind=dp)) call set_xlon(ipart,real(eps,kind=dp))
- if (abs(part(ipart)%xlon-real(nxmin1,kind=dp)).le.eps) call set_xlon(ipart,real(nxmin1-eps,kind=dp))
+ if (part(ipart)%xlon .ge. real(nxmin1, kind=dp)) &
+ call update_xlon(ipart,-real(nxmin1, kind=dp))
+ if (part(ipart)%xlon .lt. 0.) call update_xlon(ipart,real(nxmin1, kind=dp))
+ if (part(ipart)%xlon .le. real(eps, kind=dp)) &
+ call set_xlon(ipart,real(eps, kind=dp))
+ if (abs( part(ipart)%xlon - real(nxmin1, kind=dp)) .le. eps) &
+ call set_xlon(ipart,real(nxmin1-eps,kind=dp))
endif
! HSO/AL: Prevent particles from disappearing at the pole
!******************************************************************
if (sglobal .and. part(ipart)%ylat.lt.0. ) then
- call set_xlon(ipart,mod(part(ipart)%xlon+real(nxmin1/2.,kind=dp),real(nxmin1,kind=dp)))
+ call set_xlon(ipart, &
+ mod( part(ipart)%xlon + real(nxmin1*0.5, kind=dp), real(nxmin1, kind=dp)))
call set_ylat(ipart,-part(ipart)%ylat)
- ! In extremely rare cases, the ylat exceeds the bounds, so we set it back into the domain here
- if ( part(ipart)%ylat.gt.real(nymin1,kind=dp) ) then
- call set_ylat(ipart,real(nymin1,kind=dp)-mod(part(ipart)%ylat,real(nymin1,kind=dp)))
- endif
- else if (nglobal .and. part(ipart)%ylat.gt.real(nymin1,kind=dp) ) then
- call set_xlon(ipart,mod(part(ipart)%xlon+real(nxmin1/2.,kind=dp),real(nxmin1,kind=dp)))
- call set_ylat(ipart,2.*real(nymin1,kind=dp)-part(ipart)%ylat)
+ ! In extremely rare cases, the ylat exceeds the bounds,
+ ! so we set it back into the domain here
+ if ( part(ipart)%ylat.gt.real(nymin1,kind=dp) ) &
+ call set_ylat(ipart, &
+ real(nymin1, kind=dp) - mod( part(ipart)%ylat, real(nymin1, kind=dp)))
+ else if (nglobal .and. part(ipart)%ylat .gt. real(nymin1, kind=dp) ) then
+ call set_xlon(ipart, &
+ mod( part(ipart)%xlon + real(nxmin1*0.5, kind=dp), real(nxmin1, kind=dp)))
+ call set_ylat(ipart,2.*real(nymin1,kind=dp) - part(ipart)%ylat)
endif
! Check position: If trajectory outside model domain, terminate it
!*****************************************************************
- ! Not necessary to check when using global domain, but some problems in the meteo data could cause particles
- ! to go crazy.
+ ! Not necessary to check when using global domain, but some problems in the
+ ! meteo data could cause particles to go crazy.
! if (gdomainfill) return
- if ((part(ipart)%xlon.lt.0.).or.(part(ipart)%xlon.ge.real(nxmin1,kind=dp)).or.(part(ipart)%ylat.lt.0.).or. &
- (part(ipart)%ylat.gt.real(nymin1,kind=dp))) then
+
+ if (part(ipart)%xlon.lt.0. .or. part(ipart)%xlon.ge.real(nxmin1,kind=dp) &
+ .or. part(ipart)%ylat.lt.0. .or. part(ipart)%ylat.gt.real(nymin1,kind=dp)) then
part(ipart)%nstop=.true.
return
endif
-end subroutine advance_updateXY
+
+end subroutine update_xy
-subroutine advance_adjusttopheight(ipart)
+subroutine pushpartdown(ipart)
implicit none
@@ -788,15 +873,24 @@ subroutine advance_adjusttopheight(ipart)
eps=nxmax/3.e5
select case (wind_coord_type)
+
case ('ETA')
- if (part(ipart)%zeta.le.real(uvheight(nz),kind=dp)) then
+
+ if (part(ipart)%zeta.le.real(uvheight(nz),kind=dp)) &
call set_zeta(ipart,uvheight(nz)+eps_eta)
- endif
+
case ('METER')
- if (part(ipart)%z.ge.real(height(nz),kind=dp)) call set_z(ipart,height(nz)-100.*eps)
+
+ if (part(ipart)%z.ge.real(height(nz),kind=dp)) &
+ call set_z(ipart,height(nz)-100.*eps)
+
case default
- if (part(ipart)%z.ge.real(height(nz),kind=dp)) call set_z(ipart,height(nz)-100.*eps)
+
+ if (part(ipart)%z.ge.real(height(nz),kind=dp)) &
+ call set_z(ipart,height(nz)-100.*eps)
+
end select
-end subroutine advance_adjusttopheight
+
+end subroutine pushpartdown
end module advance_mod
diff --git a/src/binary_output_mod.f90 b/src/binary_output_mod.f90
index 57a7b933d62f2226c5ae98f97db07ed8c7ed82a8..0a194f56e9786d8f4dec1d2a3e8c5bb1bad66fa5 100644
--- a/src/binary_output_mod.f90
+++ b/src/binary_output_mod.f90
@@ -6,9 +6,9 @@
! This module contains routines that output gridded data to binary files. *
! *
! Not all routines that should have a netcdf equivalent, have one yet: *
- ! writeheader_bin_sfc_nest,writeheader_bin_sfc,concoutput_sfc, *
- ! concoutput_sfc_nest,initcond_output,initcond_output_inversion, *
- ! concoutput_inv_nest *
+ ! writeheader_bin_sfc_nest,writeheader_bin_sfc,concoutput_sfc, *
+ ! concoutput_sfc_nest,initcond_output,initcond_output_inv, *
+ ! concoutput_inv_nest *
! *
! L. Bakels 2022 *
! *
@@ -4388,7 +4388,7 @@ subroutine initcond_output(itime)
end do
end subroutine initcond_output
-subroutine initcond_output_inversion(itime)
+subroutine initcond_output_inv(itime)
! i
!*****************************************************************************
! *
@@ -4542,6 +4542,6 @@ subroutine initcond_output_inversion(itime)
if (listart) then
listart=.false.
endif
-end subroutine initcond_output_inversion
+end subroutine initcond_output_inv
-end module binary_output_mod
\ No newline at end of file
+end module binary_output_mod
diff --git a/src/cbl_mod.f90 b/src/cbl_mod.f90
index c1fa2887e26564c1712b682485b2f054cc33588c..a618437bcfc5c636382c82eb901458a945f43323 100644
--- a/src/cbl_mod.f90
+++ b/src/cbl_mod.f90
@@ -6,7 +6,7 @@ module cbl_mod
private :: cuberoot
- public :: cbl,re_initialize_particle,initialize_cbl_vel
+ public :: cbl,reinit_particle,init_cbl_vel
contains
@@ -235,7 +235,7 @@ function cuberoot(x) result(y)
y=sign((abs(x))**third,x)
end function cuberoot
-subroutine re_initialize_particle(zp,ust,wst,h,sigmaw,wp,nrand,ol)
+subroutine reinit_particle(zp,ust,wst,h,sigmaw,wp,nrand,ol)
! i i i i i io io i
!******************************************************************************
@@ -330,9 +330,9 @@ subroutine re_initialize_particle(zp,ust,wst,h,sigmaw,wp,nrand,ol)
end if
return
-end subroutine re_initialize_particle
+end subroutine reinit_particle
-subroutine initialize_cbl_vel(idum,zp,ust,wst,h,sigmaw,wp,ol,ithread)
+subroutine init_cbl_vel(idum,zp,ust,wst,h,sigmaw,wp,ol,ithread)
! i/o i i i i i o i
use par_mod, only:pi
@@ -414,6 +414,6 @@ subroutine initialize_cbl_vel(idum,zp,ust,wst,h,sigmaw,wp,ol,ithread)
end if
return
-end subroutine initialize_cbl_vel
+end subroutine init_cbl_vel
end module cbl_mod
diff --git a/src/com_mod.f90 b/src/com_mod.f90
index 27ed37d30951d1d18604fc4314a9a4bc69dd3855..29c51cddfb6742a64a1494f872030dada862334d 100644
--- a/src/com_mod.f90
+++ b/src/com_mod.f90
@@ -499,7 +499,7 @@ module com_mod
contains
- subroutine com_mod_alloc_part(nmpart)
+ subroutine mpi_alloc_part(nmpart)
!*******************************************************************************
! Dynamic allocation of arrays
!
@@ -521,6 +521,6 @@ contains
allocate(part_av_uu(nmpart),part_av_vv(nmpart),part_av_energy(nmpart))
end if
- end subroutine com_mod_alloc_part
+ end subroutine mpi_alloc_part
end module com_mod
diff --git a/src/conv_mod.f90 b/src/conv_mod.f90
index 47f12c9aeb5dcfba6ac9c945612a8268525fb16d..b295c2ebc9a2a2e78865b5de72935e3a911211f4 100644
--- a/src/conv_mod.f90
+++ b/src/conv_mod.f90
@@ -74,7 +74,7 @@ module conv_mod
contains
-subroutine alloc_convection
+subroutine alloc_convect
implicit none
if (.not.lconvection.eq.1) return
! ! nconvlevmax=nuvzmax-1
@@ -95,9 +95,9 @@ subroutine alloc_convection
! QP(NA),EP(NA),TH(NA),WT(NA),EVAP(NA),CLW(NA), &
! SIGP(NA),TP(NA),CPN(NA), &
! LV(NA),LVCP(NA),H(NA),HP(NA),GZ(NA),HM(NA))
-end subroutine alloc_convection
+end subroutine alloc_convect
-subroutine dealloc_convection
+subroutine dealloc_convect
implicit none
if (.not.lconvection.eq.1) return
! deallocate(pconv,phconv,dpr,pconv_hpa,phconv_hpa,ft,fq,sub, &
@@ -106,13 +106,12 @@ subroutine dealloc_convection
! deallocate(fup,fdown,ment,M,MP,QENT,ELIJ,SIJ,TVP,TV, &
! WATER,QP,EP,TH,WT,EVAP,CLW,SIGP,TP,CPN,LV,LVCP, &
! H,HP,GZ,HM)
-end subroutine dealloc_convection
+end subroutine dealloc_convect
-subroutine set_upperlevel_convect()
- ! Determine the uppermost level for which the convection scheme shall be applied
- ! by assuming that there is no convection above 50 hPa (for standard SLP)
- !*****************************************************************************
- implicit none
+subroutine set_conv_top()
+! Determine the uppermost level for which the convection scheme shall be applied
+! by assuming that there is no convection above 50 hPa (for standard SLP)
+!*****************************************************************************
integer :: i
real :: pint
@@ -124,10 +123,11 @@ subroutine set_upperlevel_convect()
nconvlev=i
if (nconvlev.gt.nconvlevmax-1) then
nconvlev=nconvlevmax-1
- write(*,*) 'Attention, convection only calculated up to ', &
+ write(*,*) 'INFORMATION: Convection only calculated up to ', &
akz(nconvlev)+bkz(nconvlev)*1013.25,' hPa'
endif
-end subroutine set_upperlevel_convect
+
+end subroutine set_conv_top
subroutine convmix(itime)
! i
@@ -199,10 +199,10 @@ subroutine convmix(itime)
! if no particles are present return after initialization
!********************************************************
- call get_alive_part_num(alivepart)
+ call get_alivepart_num(alivepart)
if (alivepart.le.0 ) return
- call get_total_part_num(totpart)
+ call get_totalpart_num(totpart)
allocate( igrid(totpart) )
allocate( ipoint(totpart) )
allocate( igridn(totpart,maxnests) )
diff --git a/src/coord_ec_mod.f90 b/src/coord_ec_mod.f90
index b051cfdfbeca98598a93ec4544346bb3e443b6d7..5bdfc0035f1c6714ccf572b9289f7bd2c6b134d9 100644
--- a/src/coord_ec_mod.f90
+++ b/src/coord_ec_mod.f90
@@ -95,9 +95,9 @@ subroutine z_to_zeta(itime,xt,yt,zold,zteta)
endif
call find_ngrid(xt,yt)
- call determine_grid_coordinates(real(xt),real(yt))
+ call find_grid_indices(real(xt),real(yt))
call find_grid_distances(real(xt),real(yt))
- call find_time_variables(itime)
+ call find_time_vars(itime)
! Integration method as used in the original verttransform_ecmwf.f90
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -213,9 +213,9 @@ subroutine zeta_to_z(itime,xt,yt,zteta,ztout)
! Convert eta z coordinate to meters
!***********************************
call find_ngrid(xt,yt)
- call determine_grid_coordinates(real(xt),real(yt))
+ call find_grid_indices(real(xt),real(yt))
call find_grid_distances(real(xt),real(yt))
- call find_time_variables(itime)
+ call find_time_vars(itime)
k=nz-1
frac=1.
@@ -343,9 +343,9 @@ subroutine z_to_zeta_lin(itime,xt,yt,zold,zteta)
prx,pr1,pr2 ! pressure of encompassing levels
call find_ngrid(xt,yt)
- call determine_grid_coordinates(real(xt),real(yt))
+ call find_grid_indices(real(xt),real(yt))
call find_grid_distances(real(xt),real(yt))
- call find_time_variables(itime)
+ call find_time_vars(itime)
! Integration method as used in the original verttransform_ecmwf.f90
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -439,9 +439,9 @@ subroutine zeta_to_z_lin(itime,xt,yt,zteta,ztout)
! Convert eta z coordinate to meters
!***********************************
call find_ngrid(xt,yt)
- call determine_grid_coordinates(real(xt),real(yt))
+ call find_grid_indices(real(xt),real(yt))
call find_grid_distances(real(xt),real(yt))
- call find_time_variables(itime)
+ call find_time_vars(itime)
k=nz-1
frac=1.
diff --git a/src/drydepo_mod.f90 b/src/drydepo_mod.f90
index 492bdb9ff99d504f21b8fc754f12e715f9a52e25..231348cfce9f134540664f32f20255692d56b19e 100644
--- a/src/drydepo_mod.f90
+++ b/src/drydepo_mod.f90
@@ -773,7 +773,7 @@ subroutine get_vdep_prob(itime,xt,yt,zt,prob)
! Determine nested grid coordinates
!**********************************
- call determine_grid_coordinates(xt,yt)
+ call find_grid_indices(xt,yt)
! Determine probability of deposition
!************************************
diff --git a/src/initialise_mod.f90 b/src/initialise_mod.f90
index 82713fc74997a32faa9490b1a42d84a5838204ca..e9fe50988cfd1d0a12d27848d8016529154af596 100644
--- a/src/initialise_mod.f90
+++ b/src/initialise_mod.f90
@@ -147,7 +147,7 @@ subroutine releaseparticle(itime)
if (totpart.gt.count%allocated) call alloc_particles(totpart-count%allocated)
end if
- call get_total_part_num(istart)
+ call get_totalpart_num(istart)
minpart=1
do i=1,numpoint
if ((itime.ge.ireleasestart(i)).and. &! are we within release interval?
@@ -215,7 +215,7 @@ subroutine releaseparticle(itime)
zaux=zpoint2(i)-zpoint1(i)
do j=1,numrel ! loop over particles to be released this time
- call get_new_part_index(ipart)
+ call get_newpart_index(ipart)
call spawn_particle(itime, ipart)
! Particle coordinates are determined by using a random position within the release volume
@@ -445,18 +445,18 @@ subroutine releaseparticle(itime)
end do
endif
- call get_total_part_num(numpart)
+ call get_totalpart_num(numpart)
end do ! numrel
endif ! releasepoint
end do ! numpoint
- call get_total_part_num(iend)
+ call get_totalpart_num(iend)
! NetCDF only: write initial positions of new particles
#ifdef USE_NCF
if ((iend-istart.gt.0).and.(ipout.ge.1)) then
- call wrt_particles_initialpos(itime,istart,iend)
+ call wrt_part_initialpos(itime,istart,iend)
call output_particle(itime,.true.)
endif
#endif
@@ -514,8 +514,8 @@ subroutine readpartpositions
!*****************************************
if (lnetcdfout.eq.1) then
#ifdef USE_NCF
- call readpartpositions_ncf(ibtime,ibdate)
- call get_total_part_num(numpart)
+ call read_partpos_ncf(ibtime,ibdate)
+ call get_totalpart_num(numpart)
numparticlecount=numpart
return
#endif
@@ -795,7 +795,7 @@ subroutine read_heightlevels(height_tmp,nmixz_tmp)
write(*,*) ' #### FROM VERTTRANSFORM_MOD. #### '
end subroutine read_heightlevels
-subroutine initialize_particle(itime,ipart)
+subroutine init_particle(itime,ipart)
! i i o o o
! o o o i i i o
!*****************************************************************************
@@ -879,7 +879,7 @@ subroutine initialize_particle(itime,ipart)
call find_ngrid(xt,yt)
! Compute maximum mixing height around particle position
!*******************************************************
- call determine_grid_coordinates(xt,yt)
+ call find_grid_indices(xt,yt)
h=max(hmix(ix ,jy,1,memind(1)), &
hmix(ixp,jy ,1,memind(1)), &
@@ -935,7 +935,7 @@ subroutine initialize_particle(itime,ipart)
if(-h/ol.gt.5) then
!if (ol.lt.0.) then
!if (ol.gt.0.) then !by mc : only for test correct is lt.0
- call initialize_cbl_vel( &
+ call init_cbl_vel( &
iseed1(thread),zt,ust,wst,h,sigw,part(ipart)%turbvel%w,ol,thread)
else
part(ipart)%turbvel%w=part(ipart)%turbvel%w*sigw
@@ -1027,7 +1027,7 @@ subroutine initialize_particle(itime,ipart)
part(ipart)%mesovel%w=rannumb(nrand+2)*wsig
end select
endif
-end subroutine initialize_particle
+end subroutine init_particle
subroutine init_domainfill
!
@@ -1705,7 +1705,7 @@ subroutine boundcond_domainfill(itime,loutend)
endif
do m=1,mmass
- call get_new_part_index(ipart)
+ call get_newpart_index(ipart)
call spawn_particle(itime, ipart)
! Assign particle positions
@@ -1920,7 +1920,7 @@ subroutine boundcond_domainfill(itime,loutend)
endif
do m=1,mmass
- call get_new_part_index(ipart)
+ call get_newpart_index(ipart)
call spawn_particle(itime, ipart)
! Assign particle positions
diff --git a/src/interpol_mod.f90 b/src/interpol_mod.f90
index 709f721c61a2b838702beb07906cabe321feb720..55e35ec2cd3df655e870b0602ffc995028c5d48d 100644
--- a/src/interpol_mod.f90
+++ b/src/interpol_mod.f90
@@ -33,7 +33,7 @@ module interpol_mod
logical :: lbounds(2),lbounds_w(2),lbounds_uv(2) ! marking particles below or above bounds
private :: interpol_wind_meter,interpol_wind_eta
- private :: standard_deviation_meter,standard_deviation_eta
+ private :: stdev_meter,stdev_eta
private :: interpol_partoutput_val_eta,interpol_partoutput_val_meter
interface hor_interpol
@@ -46,7 +46,7 @@ module interpol_mod
interface find_ngrid
- procedure find_ngrid_dp, find_ngrid_float
+ procedure find_ngrid_dp, find_ngrid_sp
end interface find_ngrid
!$OMP THREADPRIVATE(uprof,vprof,wprof,usigprof,vsigprof,wsigprof, &
!$OMP rhoprof,rhogradprof,u,v,w,usig,vsig,wsig, &
@@ -70,28 +70,27 @@ subroutine dealloc_interpol
! rhoprof,rhogradprof,indzindicator)
end subroutine dealloc_interpol
-subroutine initialise_interpol_mod(itime,xt,yt,zt,zteta)
+subroutine init_interpol(itime,xt,yt,zt,zteta)
+
! This routine initialises all important values used in the interpol module
! This includes:
! - The current grid number in which the particle is positioned
! - The interpolation fractions of the grid (x,y,z) and of time
- implicit none
-
integer, intent(in) :: itime ! time step
real, intent(in) :: xt,yt ! particle positions
real, intent(in) :: zt ! height in meters
real, intent(in) :: zteta ! height in eta coordinates
call find_ngrid(xt,yt)
- call determine_grid_coordinates(xt,yt)
+ call find_grid_indices(xt,yt)
call find_grid_distances(xt,yt)
- call find_time_variables(itime)
+ call find_time_vars(itime)
call find_z_level(zt,zteta)
-end subroutine initialise_interpol_mod
-subroutine determine_grid_coordinates(xt,yt)
- implicit none
+end subroutine init_interpol
+
+subroutine find_grid_indices(xt,yt)
real, intent(in) :: xt,yt ! particle positions
@@ -128,7 +127,8 @@ subroutine determine_grid_coordinates(xt,yt)
write(*,*) 'WARNING: interpol_mod.f90 ixp >= nxmax. xt,yt:',xt,yt
ixp=ixp-nxmax
end if
-end subroutine determine_grid_coordinates
+
+end subroutine find_grid_indices
subroutine find_grid_distances(xt,yt)
@@ -149,21 +149,21 @@ subroutine find_grid_distances(xt,yt)
p2=ddx*rddy
p3=rddx*ddy
p4=ddx*ddy
+
end subroutine find_grid_distances
-subroutine find_time_variables(itime)
-
- implicit none
+subroutine find_time_vars(itime)
integer, intent(in) :: itime ! time step
dt1=real(itime-memtime(1))
dt2=real(memtime(2)-itime)
dtt=1./(dt1+dt2)
-end subroutine find_time_variables
+
+end subroutine find_time_vars
subroutine find_z_level(zt,zteta)
- implicit none
+
real, intent(in) :: &
zt, & ! height in meters
zteta ! height in eta
@@ -177,10 +177,11 @@ subroutine find_z_level(zt,zteta)
case default
call find_z_level_meters(zt)
end select
+
end subroutine find_z_level
subroutine find_z_level_meters(zt)
- implicit none
+
real, intent(in) :: zt ! height in meters
integer :: i
@@ -205,20 +206,22 @@ subroutine find_z_level_meters(zt)
endif
end do
endif
+
end subroutine find_z_level_meters
subroutine find_z_level_eta(zteta)
- implicit none
+
real, intent(in) :: zteta ! height in eta coordinates
integer :: i ! loop variable
call find_z_level_eta_w(zteta)
call find_z_level_eta_uv(zteta)
+
end subroutine find_z_level_eta
subroutine find_z_level_eta_w(zteta)
- implicit none
+
real, intent(in) :: zteta ! height in eta coordinates
integer :: i ! loop variable
@@ -244,10 +247,11 @@ subroutine find_z_level_eta_w(zteta)
endif
end do
endif
+
end subroutine find_z_level_eta_w
subroutine find_z_level_eta_uv(zteta)
- implicit none
+
real, intent(in) :: zteta ! height in eta coordinates
integer :: i ! loop variable
@@ -272,32 +276,37 @@ subroutine find_z_level_eta_uv(zteta)
endif
end do
endif
+
end subroutine find_z_level_eta_uv
-subroutine find_vert_variables(vertlevels,zpos,zlevel,dz1,dz2,bounds,wlevel)
+subroutine find_vert_vars(vertlevels,zpos,zlevel,dz1,dz2,bounds,wlevel)
+
!*****************************************************************************
! *
! This subroutine computes the vertical interpolation variables *
- ! logarithmically, unless log_interpol=.false. in the par_mod *
+ ! logarithmically, unless log_interpol=.false. in the par_mod *
! *
! Author: L. Bakels *
!*****************************************************************************
- implicit none
- real, intent(in) :: vertlevels(:) ! vertical levels in coordinate system
- real, intent(in) :: zpos ! verticle particle position
- integer, intent(in) :: zlevel ! vertical level of interest
- logical, intent(in) :: bounds(2),wlevel ! flag marking if particles are outside bounds
- real, intent(inout) :: dz1,dz2 ! fractional distance to point 1 (closer to ground) and 2
+ real, intent(in) :: vertlevels(:) ! vertical levels in coordinate system
+ real, intent(in) :: zpos ! verticle particle position
+ integer, intent(in) :: zlevel ! vertical level of interest
+ logical, intent(in) :: bounds(2),wlevel ! flag marking if particles are
+ ! outside bounds
+ real, intent(inout) :: dz1,dz2 ! fractional distance to point 1
+ ! (closer to ground) and 2
real :: dz,dh1,dh,pfact
- real :: psint1(2),psint,pr1,pr2,pr_test ! pressure of encompassing levels
+ real :: psint1(2),psint,pr1,pr2,pr_test
+ ! pressure of encompassing levels
integer :: m
! Only do logarithmic interpolation when using ETA coordinates, since the
! levels are following pressure, while METER levels are linear.
!##############################################################
+
if (.not. log_interpol) then
- call find_vert_variables_lin(vertlevels,zpos,zlevel,dz1,dz2,bounds,wlevel)
+ call find_vert_vars_lin(vertlevels,zpos,zlevel,dz1,dz2,bounds,wlevel)
return
endif
@@ -360,17 +369,20 @@ subroutine find_vert_variables(vertlevels,zpos,zlevel,dz1,dz2,bounds,wlevel)
! dz1=(log(zpos)-log(vertlevels(zlevel)))*dz
! dz2=(log(vertlevels(zlevel+1))-log(zpos))*dz
! endif
-end subroutine find_vert_variables
-subroutine find_vert_variables_lin(vertlevels,zpos,zlevel,dz1,dz2,bounds,wlevel)
- implicit none
- real, intent(in) :: vertlevels(:) ! vertical levels in coordinate system
- real, intent(in) :: zpos ! verticle particle position
- integer, intent(in) :: zlevel ! vertical level of interest
- logical, intent(in) :: bounds(2),wlevel ! flag marking if particles are outside bounds
- real, intent(inout) :: dz1,dz2 ! fractional distance to point 1 (closer to ground) and 2
+end subroutine find_vert_vars
+
+subroutine find_vert_vars_lin(vertlevels,zpos,zlevel,dz1,dz2,bounds,wlevel)
+
+ real, intent(in) :: vertlevels(:) ! vertical levels in coordinate system
+ real, intent(in) :: zpos ! verticle particle position
+ integer, intent(in) :: zlevel ! vertical level of interest
+ logical, intent(in) :: bounds(2),wlevel! flag marking if particles are outside
+ ! bounds
+ real, intent(inout) :: dz1,dz2 ! fractional distance to point 1
+ ! (closer to ground) and 2
real :: dz,dh1,dh,pfact
- real :: psint1(2),psint,pr1,pr2,temp ! pressure of encompassing levels
+ real :: psint1(2),psint,pr1,pr2,temp ! pressure of encompassing levels
! If the particle is below bounds (bounds(1)==.true.):
if (bounds(1)) then
@@ -385,17 +397,14 @@ subroutine find_vert_variables_lin(vertlevels,zpos,zlevel,dz1,dz2,bounds,wlevel)
dz1=(zpos-vertlevels(zlevel))*dz
dz2=(vertlevels(zlevel+1)-zpos)*dz
endif
-end subroutine find_vert_variables_lin
+
+end subroutine find_vert_vars_lin
subroutine find_ngrid_dp(xt,yt)
- implicit none
- real :: &
- eps
- real(kind=dp), intent(in) :: &
- xt,yt ! particle positions on grid
- integer :: &
- j
+ real eps
+ real(kind=dp), intent(in) :: xt,yt ! particle positions on grid
+ integer :: j
eps=nxmax/3.e5
if (nglobal.and.(real(yt).gt.switchnorthg)) then
@@ -406,59 +415,57 @@ subroutine find_ngrid_dp(xt,yt)
ngrid=0
! Temporary fix for nested layer edges: replaced eps with dxn and dyn (LB)
do j=numbnests,1,-1
- if ((real(xt).gt.xln(j)+dxn(j)).and.(real(xt).lt.xrn(j)-dxn(j)).and. &
- (real(yt).gt.yln(j)+dyn(j)).and.(real(yt).lt.yrn(j)-dyn(j))) then
+ if (real(xt).gt.xln(j)+dxn(j) .and. real(xt).lt.xrn(j)-dxn(j) .and. &
+ real(yt).gt.yln(j)+dyn(j) .and. real(yt).lt.yrn(j)-dyn(j)) then
ngrid=j
exit
endif
end do
endif
+
end subroutine find_ngrid_dp
-subroutine find_ngrid_float(xt,yt)
+subroutine find_ngrid_sp(xt,yt)
- implicit none
- real :: &
- eps
- real, intent(in) :: &
- xt,yt ! particle positions on grid
- integer :: &
- j
+ real :: eps
+ real, intent(in) :: xt,yt ! particle positions on grid
+ integer :: j
eps=nxmax/3.e5
- if (nglobal.and.(yt.gt.switchnorthg)) then
+ if (nglobal .and. yt.gt.switchnorthg) then
ngrid=-1
- else if (sglobal.and.(yt.lt.switchsouthg)) then
+ else if (sglobal .and. yt.lt.switchsouthg) then
ngrid=-2
else
ngrid=0
! Temporary fix for nested layer edges: replaced eps with dxn and dyn (LB)
do j=numbnests,1,-1
- if ((xt.gt.xln(j)+dxn(j)).and.(xt.lt.xrn(j)-dxn(j)).and. &
- (yt.gt.yln(j)+dyn(j)).and.(yt.lt.yrn(j)-dyn(j))) then
+ if (xt.gt.xln(j)+dxn(j) .and. xt.lt.xrn(j)-dxn(j) .and. &
+ yt.gt.yln(j)+dyn(j) .and. yt.lt.yrn(j)-dyn(j)) then
ngrid=j
exit
endif
end do
endif
-end subroutine find_ngrid_float
-subroutine hor_interpol_4d(field,output,zlevel,indexh,ztot)
+end subroutine find_ngrid_sp
- implicit none
+subroutine hor_interpol_4d(field,output,zlevel,indexh,ztot)
- integer, intent(in) :: zlevel,ztot,indexh ! interpolation z level, z
- real, intent(in) :: field(0:nxmax-1,0:nymax-1,ztot,numwfmem) ! input field to interpolate over
- real, intent(inout) :: output ! interpolated values
+ integer, intent(in) :: zlevel,ztot,indexh ! interpolation z level, z
+ real, intent(in) :: field(0:nxmax-1,0:nymax-1,ztot,numwfmem)
+ ! input field to interpolate
+ real, intent(inout) :: output ! interpolated values
output=p1*field(ix ,jy ,zlevel,indexh) &
+ p2*field(ixp,jy ,zlevel,indexh) &
+ p3*field(ix ,jyp,zlevel,indexh) &
+ p4*field(ixp,jyp,zlevel,indexh)
+
end subroutine hor_interpol_4d
subroutine hor_interpol_2d(field,output)
- implicit none
+
real, intent(in) :: field(0:nxmax-1,0:nymax-1) ! 2D imput field
real, intent(inout) :: output ! Interpolated value
@@ -470,145 +477,159 @@ end subroutine hor_interpol_2d
subroutine hor_interpol_4d_nest(field,output,zlevel,indexh,ztot)
- implicit none
-
- integer, intent(in) :: zlevel,ztot,indexh ! interpolation z level, z
- real, intent(in) :: field(0:nxmaxn-1,0:nymaxn-1,ztot,numwfmem,numbnests) ! input field to interpolate over
- real, intent(inout) :: output ! interpolated values
+ integer, intent(in) :: zlevel,ztot,indexh ! interpolation z level, z
+ real, intent(in) :: field(0:nxmaxn-1,0:nymaxn-1,ztot,numwfmem,numbnests)
+ ! input field to interpolate
+ real, intent(inout) :: output ! interpolated values
output=p1*field(ix ,jy ,zlevel,indexh,ngrid) &
+ p2*field(ixp,jy ,zlevel,indexh,ngrid) &
+ p3*field(ix ,jyp,zlevel,indexh,ngrid) &
+ p4*field(ixp,jyp,zlevel,indexh,ngrid)
+
end subroutine hor_interpol_4d_nest
subroutine hor_interpol_2d_nest(field,output)
- implicit none
-
- real, intent(in) :: field(0:nxmaxn-1,0:nymaxn-1,numbnests) ! input field to interpolate over
- real, intent(inout) :: output ! interpolated values
+ real, intent(in) :: field(0:nxmaxn-1,0:nymaxn-1,numbnests)
+ ! input field to interpolate
+ real, intent(inout) :: output ! interpolated values
output=p1*field(ix ,jy ,ngrid) &
+ p2*field(ixp,jy ,ngrid) &
+ p3*field(ix ,jyp,ngrid) &
+ p4*field(ixp,jyp,ngrid)
+
end subroutine hor_interpol_2d_nest
subroutine temp_interpol(time1,time2,output)
- implicit none
-
real, intent(in) :: time1,time2 ! input data at two timesteps
real, intent(inout) :: output ! interpolated data
- output=(time1*dt2+time2*dt1)*dtt
+ output= (time1*dt2 + time2*dt1) * dtt
+
end subroutine temp_interpol
subroutine vert_interpol(input1,input2,dz1,dz2,output)
- implicit none
-
- real, intent(in) :: input1,input2 ! input data at two vertical levels, 1 being closer to ground
+ real, intent(in) :: input1,input2 ! input data at two vertical levels,
+ ! 1 being closer to ground
real, intent(in) :: dz1,dz2 ! logarithmic interpolation values
real, intent(inout) :: output ! interpolated data
- output = input1*dz2 + input2*dz1!input1**dz2 * input2**dz1
+ output = input1*dz2 + input2*dz1 ! input1**dz2 * input2**dz1
+
end subroutine vert_interpol
subroutine bilin_spatial_interpol(field,output,zlevel,dz1,dz2,ztot)
- implicit none
- integer, intent(in) :: zlevel,ztot ! interpolation z level
- real, intent(in) :: field(0:nxmax-1,0:nymax-1,ztot,numwfmem) ! input field to interpolate over
+
+ integer, intent(in) :: zlevel,ztot ! interpolation z level
+ real, intent(in) :: field(0:nxmax-1,0:nymax-1,ztot,numwfmem)
+ ! input field to interpolate
real, intent(in) :: dz1,dz2
- real, intent(inout) :: output(2) ! interpolated values
+ real, intent(inout) :: output(2) ! interpolated values
integer :: m,n,indzh
real :: output1(2)
do m=1,2
+
do n=1,2
indzh=zlevel+n-1
call hor_interpol_4d(field,output1(n),indzh,memind(m),ztot)
end do
+
!**********************************
! 2.) Linear vertical interpolation on logarithmic scale
!**********************************
+
call vert_interpol(output1(1),output1(2),dz1,dz2,output(m))
+
end do
+
end subroutine bilin_spatial_interpol
subroutine bilin_spatial_interpol_nest(field,output,zlevel,dz1,dz2,ztot)
- implicit none
- integer, intent(in) :: zlevel,ztot ! interpolation z level
- real, intent(in) :: field(0:nxmaxn-1,0:nymaxn-1,ztot,numwfmem,numbnests) ! input field to interpolate over
+
+ integer, intent(in) :: zlevel,ztot ! interpolation z level
+ real, intent(in) :: field(0:nxmaxn-1,0:nymaxn-1,ztot,numwfmem,numbnests)
+ ! input field to interpolate
real, intent(in) :: dz1,dz2
- real, intent(inout) :: output(2) ! interpolated values
+ real, intent(inout) :: output(2) ! interpolated values
integer :: m,n,indzh
real :: output1(2)
do m=1,2
+
do n=1,2
indzh=zlevel+n-1
call hor_interpol_4d_nest(field,output1(n),indzh,memind(m),ztot)
end do
+
!**********************************
! 2.) Linear vertical interpolation on logarithmic scale
!**********************************
+
call vert_interpol(output1(1),output1(2),dz1,dz2,output(m))
+
end do
+
end subroutine bilin_spatial_interpol_nest
subroutine compute_sl_sq(field,sl,sq,zlevel,indexh,ztot)
- implicit none
-
- integer, intent(in) :: zlevel,ztot,indexh ! interpolation z levels
- real, intent(in) :: field(0:nxmax-1,0:nymax-1,ztot,numwfmem) ! input field to interpolate over
- real, intent(inout) :: sl,sq ! standard deviation
+ integer, intent(in) :: zlevel,ztot,indexh ! interpolation z levels
+ real, intent(in) :: field(0:nxmax-1,0:nymax-1,ztot,numwfmem)
+ ! input field to interpolate
+ real, intent(inout) :: sl,sq ! standard deviation
sl=sl+field(ix ,jy ,zlevel,indexh)+field(ixp,jy ,zlevel,indexh) &
+field(ix ,jyp,zlevel,indexh)+field(ixp,jyp,zlevel,indexh)
+
sq=sq+field(ix ,jy ,zlevel,indexh)*field(ix ,jy ,zlevel,indexh)+ &
field(ixp,jy ,zlevel,indexh)*field(ixp,jy ,zlevel,indexh)+ &
field(ix ,jyp,zlevel,indexh)*field(ix ,jyp,zlevel,indexh)+ &
field(ixp,jyp,zlevel,indexh)*field(ixp,jyp,zlevel,indexh)
+
end subroutine compute_sl_sq
subroutine compute_sl_sq_nest(field,sl,sq,zlevel,indexh,ztot)
- implicit none
-
- integer, intent(in) :: zlevel,ztot,indexh ! interpolation z levels
- real, intent(in) :: field(0:nxmaxn-1,0:nymaxn-1,ztot,numwfmem,numbnests) ! input field to interpolate over
- real, intent(inout) :: sl,sq ! standard deviation
+ integer, intent(in) :: zlevel,ztot,indexh ! interpolation z levels
+ real, intent(in) :: field(0:nxmaxn-1,0:nymaxn-1,ztot,numwfmem,numbnests)
+ ! input field to interpolate
+ real, intent(inout) :: sl,sq ! standard deviation
sl=sl+field(ix ,jy ,zlevel,indexh,ngrid)+field(ixp,jy ,zlevel,indexh,ngrid) &
+field(ix ,jyp,zlevel,indexh,ngrid)+field(ixp,jyp,zlevel,indexh,ngrid)
+
sq=sq+field(ix ,jy ,zlevel,indexh,ngrid)*field(ix ,jy ,zlevel,indexh,ngrid)+ &
field(ixp,jy ,zlevel,indexh,ngrid)*field(ixp,jy ,zlevel,indexh,ngrid)+ &
field(ix ,jyp,zlevel,indexh,ngrid)*field(ix ,jyp,zlevel,indexh,ngrid)+ &
field(ixp,jyp,zlevel,indexh,ngrid)*field(ixp,jyp,zlevel,indexh,ngrid)
+
end subroutine compute_sl_sq_nest
-subroutine standard_deviation(sl,sq,ndivide,output)
- implicit none
+subroutine stdev(sl,sq,divisor,output)
- real, intent(in) :: sl,sq,ndivide
+ real, intent(in) :: sl,sq,divisor
real, intent(out) :: output
real :: xaux
real,parameter :: eps=1.0e-30
- xaux=sq-sl*sl/ndivide
+ xaux= sq - sl*sl/divisor
if (xaux.lt.eps) then
output=0.
else
- output=sqrt(xaux/(ndivide-1.))
+ output=sqrt(xaux/(divisor-1.))
endif
-end subroutine standard_deviation
+
+end subroutine stdev
! Interpolation functions
!************************
+
subroutine interpol_pbl(itime,xt,yt,zt,zteta)
! i i i i
!*****************************************************************************
@@ -639,8 +660,6 @@ subroutine interpol_pbl(itime,xt,yt,zt,zteta)
use turbulence_mod
- implicit none
-
integer, intent(in) :: itime
real, intent(in) :: xt,yt,zt,zteta
integer :: m,n,indexh
@@ -651,21 +670,24 @@ subroutine interpol_pbl(itime,xt,yt,zt,zteta)
! Auxiliary variables needed for interpolation
real :: ust1(2),wst1(2),oli1(2),oliaux
+
!********************************************
! Multilinear interpolation in time and space
!********************************************
! ngrid and grid coordinates have already been definded, and are included
! in the input (for nested: xtn,ytn; for not nested: xts,yts)
- !************************************************************************
+ !************************************************************************
! Determine the lower left corner and its distance to the current position
!*************************************************************************
+
call find_grid_distances(xt,yt)
! Calculate variables for time interpolation
!*******************************************
- call find_time_variables(itime)
+
+ call find_time_vars(itime)
!********************************************************
! 1. Interpolate u*, w* and Obukhov length for turbulence
@@ -687,6 +709,7 @@ subroutine interpol_pbl(itime,xt,yt,zt,zteta)
call hor_interpol_nest(olin,oli1(m),1,memind(m),1)
end do
endif
+
! b) Temporal interpolation
call temp_interpol(ust1(1),ust1(2),ust)
call temp_interpol(wst1(1),wst1(2),wst)
@@ -772,13 +795,16 @@ subroutine interpol_pbl(itime,xt,yt,zt,zteta)
! Only necessary for the Petterssen correction
if (wind_coord_type.eq.'ETA') then
- call find_vert_variables(wheight,zteta,indzeta,dz1weta,dz2weta,lbounds_w,.true.)
- call vert_interpol(wprofeta(indzeta),wprofeta(indzpeta),dz1weta,dz2weta,weta)
+ call find_vert_vars(wheight,zteta,indzeta, &
+ dz1weta,dz2weta,lbounds_w,.true.)
+ call vert_interpol(wprofeta(indzeta),wprofeta(indzpeta), &
+ dz1weta,dz2weta,weta)
endif
+
end subroutine interpol_pbl
-subroutine interpol_pbl_misslev()
- !
+subroutine interpol_pbl_misslev
+
!*****************************************************************************
! *
! This subroutine interpolates u,v,w, density and density gradients. *
@@ -801,12 +827,9 @@ subroutine interpol_pbl_misslev()
! Constants: *
! *
!*****************************************************************************
- implicit none
- integer :: n,iw(2)
real :: uh1(2),vh1(2),wh1(2),rho1(2),rhograd1(2)
- integer :: m
-
+ integer :: m,n,iw(2)
! Within the PBL, only METER coordinates are used
! with the exception of mesoscale turbulence,
@@ -816,6 +839,7 @@ subroutine interpol_pbl_misslev()
!********************************************
! Multilinear interpolation in time and space
!********************************************
+
iw(:)=(/ indz, indzp /)
do n=1,2
if (indzindicator(iw(n))) then
@@ -850,38 +874,39 @@ subroutine interpol_pbl_misslev()
indzindicator(iw(n))=.false.
endif
end do
+
end subroutine interpol_pbl_misslev
subroutine interpol_pbl_short(zt,rhoa,rhograd)
- implicit none
+
real, intent(in) :: zt
real, intent(inout) :: rhoa,rhograd
real :: dz1,dz2
- call find_vert_variables(height,zt,indz,dz1,dz2,lbounds,.false.)
+ call find_vert_vars(height,zt,indz,dz1,dz2,lbounds,.false.)
call vert_interpol(wprof(indz),wprof(indzp),dz1,dz2,w)
call vert_interpol(uprof(indz),uprof(indzp),dz1,dz2,u)
call vert_interpol(vprof(indz),vprof(indzp),dz1,dz2,v)
call vert_interpol(rhoprof(indz),rhoprof(indzp),dz1,dz2,rhoa)
call vert_interpol(rhogradprof(indz),rhogradprof(indzp),dz1,dz2,rhograd)
+
end subroutine interpol_pbl_short
subroutine interpol_mesoscale(itime,xt,yt,zt,zteta)
- use turbulence_mod
- implicit none
+ use turbulence_mod
integer, intent(in) :: itime
real, intent(in) :: xt,yt,zt,zteta
- integer :: m,indexh
integer :: iw(2),iuv(2),iweta(2)
+ integer :: m,indexh
! Where in the grid? Stereographic (ngrid<0) or nested (ngrid>0)
!***************************************************************
call find_ngrid(xt,yt)
- call determine_grid_coordinates(xt,yt)
+ call find_grid_indices(xt,yt)
! Determine the level below the current position
!***********************************************
@@ -893,19 +918,21 @@ subroutine interpol_mesoscale(itime,xt,yt,zt,zteta)
call find_z_level_eta(zteta)
iuv(:)=(/ induv, indpuv /)
iweta(:)=(/ indzeta, indzpeta /)
- call standard_deviation_eta(iw,iuv,iweta)
+ call stdev_eta(iw,iuv,iweta)
case ('METER')
iw(:)=(/ indz, indzp /)
- call standard_deviation_meter(iw)
+ call stdev_meter(iw)
case default
write(*,*) 'ERROR: wind_coord_type is not allowed ', wind_coord_type
write(*,*) 'Choose ETA or METER.'
stop
end select
+
end subroutine interpol_mesoscale
subroutine interpol_wind(itime,xt,yt,zt,zteta,pp)
! i i i i
+
!*****************************************************************************
! *
! This subroutine interpolates the wind data to current trajectory position.*
@@ -932,12 +959,8 @@ subroutine interpol_wind(itime,xt,yt,zt,zteta,pp)
! *
!*****************************************************************************
-
- implicit none
-
integer, intent(in) :: itime,pp
- real, intent(in) :: xt,yt,zt
- real, intent(in) :: zteta
+ real, intent(in) :: xt,yt,zt,zteta
integer :: iw(2),iuv(2),iweta(2)
@@ -945,7 +968,7 @@ subroutine interpol_wind(itime,xt,yt,zt,zteta,pp)
!***************************************************************
call find_ngrid(xt,yt)
- call determine_grid_coordinates(xt,yt)
+ call find_grid_indices(xt,yt)
! ! Multilinear interpolation in time and space
! !********************************************
@@ -955,7 +978,7 @@ subroutine interpol_wind(itime,xt,yt,zt,zteta,pp)
! Calculate variables for time interpolation
!*******************************************
- call find_time_variables(itime)
+ call find_time_vars(itime)
! Interpolate over the windfields depending on the prefered
! coordinate system
@@ -969,7 +992,7 @@ subroutine interpol_wind(itime,xt,yt,zt,zteta,pp)
iuv(:) = (/ induv, indpuv /)
iweta(:)= (/ indzeta, indzpeta /)
call interpol_wind_eta(zteta,iuv,iweta)
- !call standard_deviation_wind_eta(iw,iuv,iweta)
+ !call stdev_wind_eta(iw,iuv,iweta)
case ('METER')
! Determine the level below the current position for u,v
!*******************************************************
@@ -977,17 +1000,19 @@ subroutine interpol_wind(itime,xt,yt,zt,zteta,pp)
iw(:)=(/ indz, indzp /)
call interpol_wind_meter(zt,iw)
- !call standard_deviation_wind_meter(iw)
+ !call stdev_wind_meter(iw)
case default
write(*,*) 'ERROR: wind_coord_type is not allowed ', wind_coord_type
write(*,*) 'Choose ETA or METER.'
stop
end select
+
end subroutine interpol_wind
subroutine interpol_wind_short(itime,xt,yt,zt,zteta)
- ! i i i i
+! i i i i i
+
!*****************************************************************************
! *
! This subroutine interpolates the wind data to current trajectory position.*
@@ -1011,12 +1036,8 @@ subroutine interpol_wind_short(itime,xt,yt,zt,zteta)
! *
!*****************************************************************************
-
- implicit none
-
integer, intent(in) :: itime
- real, intent(in) :: xt,yt,zt
- real, intent(in) :: zteta
+ real, intent(in) :: xt,yt,zt,zteta
integer :: iw(2),iuv(2),iweta(2)
!********************************************
@@ -1026,12 +1047,12 @@ subroutine interpol_wind_short(itime,xt,yt,zt,zteta)
! Where in the grid? Stereographic (ngrid<0) or nested (ngrid>0)
!***************************************************************
call find_ngrid(xt,yt)
- call determine_grid_coordinates(xt,yt)
+ call find_grid_indices(xt,yt)
call find_grid_distances(xt,yt)
! Calculate variables for time interpolation
!*******************************************
- call find_time_variables(itime)
+ call find_time_vars(itime)
! Interpolate over the windfields depending on the prefered
! coordinate system
@@ -1060,13 +1081,15 @@ subroutine interpol_wind_short(itime,xt,yt,zt,zteta)
write(*,*) 'Choose ETA or METER.'
stop
end select
+
end subroutine interpol_wind_short
subroutine interpol_partoutput_val(fieldname,output,j)
- implicit none
+
integer, intent(in) :: j ! particle number
character(2), intent(in) :: fieldname ! input field to interpolate over
real, intent(inout) :: output
+
! Interpolate over the windfields depending on the prefered
! coordinate system
!**********************************************************
@@ -1078,18 +1101,16 @@ subroutine interpol_partoutput_val(fieldname,output,j)
case default
call interpol_partoutput_val_meter(fieldname,output,j)
end select
+
end subroutine interpol_partoutput_val
subroutine interpol_htropo_hmix(tropop,h)
- implicit none
- real, intent(inout) :: &
- tropop, & ! height of troposphere
- h ! mixing height
- real :: &
- h1(2) ! mixing height of 2 timesteps
- integer :: &
- mind, & ! windfield index
- i,j,k,m ! loop variables
+
+ real, intent(inout) :: tropop ! height of troposphere
+ real, intent(inout) :: h ! mixing height
+ real :: h1(2) ! mixing height of 2 timesteps
+ integer :: mind ! windfield index
+ integer :: i,j,k,m ! loop variables
h=0.
if (ngrid.le.0) then
@@ -1120,44 +1141,47 @@ subroutine interpol_htropo_hmix(tropop,h)
tropop=tropopausen(nix,njy,1,memind(1),ngrid)
endif
- if (interpolhmix) h=(h1(1)*dt2+h1(2)*dt1)*dtt
+ if (interpolhmix) h= (h1(1)*dt2 + h1(2)*dt1)*dtt
+
end subroutine interpol_htropo_hmix
subroutine interpol_density(itime,ipart,output)
- implicit none
-
integer, intent(in) :: itime,ipart ! time and particle index
real, intent(inout) :: output ! output density (rhoi)
integer :: ind
- real :: dz1,dz2
- real :: rhoprof(2)
+ real :: dz1,dz2,rhoprof(2)
! Where in the grid? Stereographic (ngrid<0) or nested (ngrid>0)
!***************************************************************
call find_ngrid(part(ipart)%xlon,part(ipart)%ylat)
- call determine_grid_coordinates(real(part(ipart)%xlon),real(part(ipart)%ylat))
+ call find_grid_indices(real(part(ipart)%xlon),real(part(ipart)%ylat))
call find_grid_distances(real(part(ipart)%xlon),real(part(ipart)%ylat))
- call find_time_variables(itime)
+ call find_time_vars(itime)
- ! Take density from 2nd wind field in memory (accurate enough, no time interpolation needed)
- !*****************************************************************************
+ ! Take density from 2nd wind field in memory
+ !(accurate enough, no time interpolation needed)
+ !***********************************************
+
select case (wind_coord_type)
case ('ETA')
call find_z_level_eta(real(part(ipart)%zeta))
- call find_vert_variables(uvheight,real(part(ipart)%zeta),induv,dz1,dz2,lbounds_uv,.false.)
+ call find_vert_vars(uvheight,real(part(ipart)%zeta),induv, &
+ dz1,dz2,lbounds_uv,.false.)
if (ngrid.le.0) then
do ind=induv,indpuv
call hor_interpol(rhoeta,rhoprof(ind-induv+1),ind,memind(2),nzmax)
end do
else
do ind=induv,indpuv
- call hor_interpol_nest(rhoetan,rhoprof(ind-induv+1),ind,memind(2),nzmax)
+ call hor_interpol_nest(rhoetan,rhoprof(ind-induv+1),ind,memind(2), &
+ nzmax)
end do
endif
case ('METER')
call find_z_level_meters(real(part(ipart)%z))
- call find_vert_variables(height,real(part(ipart)%z),indz,dz1,dz2,lbounds,.false.)
+ call find_vert_vars(height,real(part(ipart)%z),indz, &
+ dz1,dz2,lbounds,.false.)
if (ngrid.le.0) then
do ind=indz,indzp
call hor_interpol(rho,rhoprof(ind-indz+1),ind,memind(2),nzmax)
@@ -1170,7 +1194,9 @@ subroutine interpol_density(itime,ipart,output)
case default
stop 'wind_coord_type not defined in conccalc.f90'
end select
+
call vert_interpol(rhoprof(1),rhoprof(2),dz1,dz2,output)
+
end subroutine interpol_density
!*********************
@@ -1179,13 +1205,15 @@ end subroutine interpol_density
! Interpolation of wind fields
!*****************************
subroutine interpol_wind_eta(zteta,iuv,iweta)
- implicit none
+
+!* PRIVATE FUNCTION *
real, intent(in) :: zteta
integer,intent(in) :: iuv(2),iweta(2)
integer :: n,m
real :: uh(2),vh(2),wetah(2),uh1(2),vh1(2),wetah1(2)
real :: dz1uv,dz2uv,dz1weta,dz2weta
+
!**********************************************************************
! 1.) Bilinear horizontal interpolation
! This has to be done separately for 6 fields (Temporal(2)*Vertical(3))
@@ -1193,8 +1221,8 @@ subroutine interpol_wind_eta(zteta,iuv,iweta)
! Vertical distance to the level below and above current position
!****************************************************************
- call find_vert_variables(uvheight,zteta,induv,dz1uv,dz2uv,lbounds_uv,.false.)
- call find_vert_variables(wheight,zteta,indzeta,dz1weta,dz2weta,lbounds_w,.true.)
+ call find_vert_vars(uvheight,zteta,induv,dz1uv,dz2uv,lbounds_uv,.false.)
+ call find_vert_vars(wheight,zteta,indzeta,dz1weta,dz2weta,lbounds_w,.true.)
! Loop over 2 time steps and 2 levels
!************************************
@@ -1231,13 +1259,16 @@ subroutine interpol_wind_eta(zteta,iuv,iweta)
call vert_interpol(wetah1(1),wetah1(2),dz1weta,dz2weta,wetah(m))
end do
endif
+
call temp_interpol(uh(1),uh(2),u)
call temp_interpol(vh(1),vh(2),v)
call temp_interpol(wetah(1),wetah(2),weta)
+
end subroutine interpol_wind_eta
subroutine interpol_wind_meter(zt,iw)
- implicit none
+
+!* PRIVATE FUNCTION *
real, intent(in) :: zt
integer,intent(in) :: iw(2)
@@ -1252,7 +1283,7 @@ subroutine interpol_wind_meter(zt,iw)
! Vertical distance to the level below and above current position
!****************************************************************
- call find_vert_variables(height,zt,indz,dz1w,dz2w,lbounds,.false.)
+ call find_vert_vars(height,zt,indz,dz1w,dz2w,lbounds,.false.)
! Loop over 2 time steps and 2 levels
!************************************
@@ -1284,13 +1315,17 @@ subroutine interpol_wind_meter(zt,iw)
call vert_interpol(vh1(1),vh1(2),dz1w,dz2w,vh(m))
end do
endif
+
call temp_interpol(wh(1),wh(2),w)
call temp_interpol(uh(1),uh(2),u)
call temp_interpol(vh(1),vh(2),v)
+
end subroutine interpol_wind_meter
subroutine interpol_partoutput_val_eta(fieldname,output,j)
- implicit none
+
+!* PRIVATE FUNCTION *
+
integer, intent(in) :: j ! particle number
character(2), intent(in) :: fieldname ! input field to interpolate over
real, intent(inout) :: output
@@ -1298,7 +1333,8 @@ subroutine interpol_partoutput_val_eta(fieldname,output,j)
if (int(dz1out).eq.-1) then
call find_z_level_eta(real(part(j)%zeta))
- call find_vert_variables(uvheight,real(part(j)%zeta),induv,dz1out,dz2out,lbounds_uv,.false.)
+ call find_vert_vars(uvheight,real(part(j)%zeta),induv,dz1out,dz2out, &
+ lbounds_uv,.false.)
endif
select case(fieldname)
@@ -1346,7 +1382,7 @@ subroutine interpol_partoutput_val_eta(fieldname,output,j)
call temp_interpol(field1(1),field1(2),output)
case('WW','ww')
call find_z_level_meters(real(part(j)%z))
- call find_vert_variables(height,real(part(j)%z),indz,dz1out,dz2out,lbounds,.false.)
+ call find_vert_vars(height,real(part(j)%z),indz,dz1out,dz2out,lbounds,.false.)
if (ngrid.le.0) then
call bilin_spatial_interpol(ww,field1,induv,dz1out,dz2out,nzmax)
else
@@ -1362,10 +1398,13 @@ subroutine interpol_partoutput_val_eta(fieldname,output,j)
endif
call temp_interpol(field1(1),field1(2),output)
end select
+
end subroutine interpol_partoutput_val_eta
subroutine interpol_partoutput_val_meter(fieldname,output,j)
- implicit none
+
+!* PRIVATE FUNCTION *
+
integer, intent(in) :: j ! particle number
character(2), intent(in) :: fieldname ! input field to interpolate over
real, intent(inout) :: output
@@ -1373,7 +1412,7 @@ subroutine interpol_partoutput_val_meter(fieldname,output,j)
if (int(dz1out).eq.-1) then
call find_z_level_meters(real(part(j)%z))
- call find_vert_variables(height,real(part(j)%z),indz,dz1out,dz2out,lbounds,.false.)
+ call find_vert_vars(height,real(part(j)%z),indz,dz1out,dz2out,lbounds,.false.)
endif
select case(fieldname)
@@ -1434,17 +1473,18 @@ subroutine interpol_partoutput_val_meter(fieldname,output,j)
endif
call temp_interpol(field1(1),field1(2),output)
end select
+
end subroutine interpol_partoutput_val_meter
subroutine interpol_mixinglayer_eta(zt,zteta,rhoa,rhograd)
- implicit none
+
real, intent(in) :: zt,zteta
real, intent(inout) :: rhoa,rhograd
real :: dz1w,dz2w,dz1uv,dz2uv,dz1weta,dz2weta
- call find_vert_variables(height,zt,indz,dz1w,dz2w,lbounds,.false.)
- call find_vert_variables(uvheight,zteta,induv,dz1uv,dz2uv,lbounds_uv,.false.)
- call find_vert_variables(wheight,zteta,indzeta,dz1weta,dz2weta,lbounds_w,.true.)
+ call find_vert_vars(height,zt,indz,dz1w,dz2w,lbounds,.false.)
+ call find_vert_vars(uvheight,zteta,induv,dz1uv,dz2uv,lbounds_uv,.false.)
+ call find_vert_vars(wheight,zteta,indzeta,dz1weta,dz2weta,lbounds_w,.true.)
call vert_interpol(wprof(indz),wprof(indzp),dz1w,dz2w,w)
call vert_interpol(uprof(induv),uprof(indpuv),dz1uv,dz2uv,u)
@@ -1452,18 +1492,22 @@ subroutine interpol_mixinglayer_eta(zt,zteta,rhoa,rhograd)
call vert_interpol(rhoprof(induv),rhoprof(indpuv),dz1uv,dz2uv,rhoa)
call vert_interpol(rhogradprof(induv),rhogradprof(indpuv),dz1uv,dz2uv,rhograd)
call vert_interpol(wprofeta(indzeta),wprofeta(indzpeta),dz1weta,dz2weta,weta)
+
end subroutine interpol_mixinglayer_eta
-subroutine standard_deviation_eta(iw,iuv,iweta)
+subroutine stdev_eta(iw,iuv,iweta)
+
+!* PRIVATE FUNCTION *
+
! Standard deviation of surrounding grid points
! Only used in mesoscale turbulence calculations
!***********************************************
- implicit none
integer,intent(in) :: iw(2),iuv(2),iweta(2)
- real :: wsl,wsq,wxaux,usl,usq,uxaux,vsl,vsq,vxaux,wetasl,wetasq,wetaxaux
+ real :: wsl,wsq,wxaux,usl,usq,uxaux,vsl,vsq,vxaux,wetasl,wetasq,wetaxaux
integer :: n,m
real,parameter :: eps=1.0e-30
+
! Standard deviations
!********************
wsl=0.
@@ -1500,17 +1544,20 @@ subroutine standard_deviation_eta(iw,iuv,iweta)
end do
endif
- call standard_deviation(wsl,wsq,16.,wsig)
- call standard_deviation(usl,usq,16.,usig)
- call standard_deviation(vsl,vsq,16.,vsig)
- call standard_deviation(wetasl,wetasq,16.,wsigeta)
-end subroutine standard_deviation_eta
+ call stdev(wsl,wsq,16.,wsig)
+ call stdev(usl,usq,16.,usig)
+ call stdev(vsl,vsq,16.,vsig)
+ call stdev(wetasl,wetasq,16.,wsigeta)
+
+end subroutine stdev_eta
+
+subroutine stdev_meter(iw)
+
+!* PRIVATE FUNCTION *
-subroutine standard_deviation_meter(iw)
! Standard deviation of surrounding grid points
! Only used in mesoscale turbulence calculations
!***********************************************
- implicit none
integer,intent(in) :: iw(2)
real :: wsl,wsq,wxaux,usl,usq,uxaux,vsl,vsq,vxaux
@@ -1549,9 +1596,10 @@ subroutine standard_deviation_meter(iw)
end do
endif
- call standard_deviation(wsl,wsq,16.,wsig)
- call standard_deviation(usl,usq,16.,usig)
- call standard_deviation(vsl,vsq,16.,vsig)
-end subroutine standard_deviation_meter
+ call stdev(wsl,wsq,16.,wsig)
+ call stdev(usl,usq,16.,usig)
+ call stdev(vsl,vsq,16.,vsig)
+
+end subroutine stdev_meter
end module interpol_mod
diff --git a/src/netcdf_output_mod.f90 b/src/netcdf_output_mod.f90
index 9eb6f88b2a91457151754294e82b8f74a2b6b043..f96e2f3eb728a0f948750cb2c0e543582bff2215 100644
--- a/src/netcdf_output_mod.f90
+++ b/src/netcdf_output_mod.f90
@@ -119,8 +119,8 @@ module netcdf_output_mod
public :: writeheader_ncf, concoutput_sfc_ncf_nest, concoutput_ncf, &
concoutput_ncf_nest, concoutput_sfc_ncf, writeheader_partoutput, partoutput_ncf, &
- open_partoutput_file, close_partoutput_file, readpartpositions_ncf, particle_initialpos, &
- wrt_particles_initialpos, topo_written, mass_written, partinit_ncf, open_partinit_file, &
+ open_partoutput_file, close_partoutput_file, read_partpos_ncf, particle_initialpos, &
+ wrt_part_initialpos, topo_written, mass_written, partinit_ncf, open_partinit_file, &
read_init_cond_nc, partinitpointer1, tpointer
contains
@@ -341,7 +341,7 @@ subroutine writeheader_ncf(lnest)
! If starting from a restart file, new data will be added to the existing grid file
if ((ipin.eq.1).or.(ipin.eq.4)) then
- call read_gridIDs(lnest)
+ call read_grid_id(lnest)
return
endif
@@ -747,7 +747,7 @@ subroutine writeheader_ncf(lnest)
return
end subroutine writeheader_ncf
-subroutine read_gridIDs(lnest)
+subroutine read_grid_id(lnest)
implicit none
logical, intent(in) :: lnest
@@ -805,7 +805,7 @@ subroutine read_gridIDs(lnest)
call nf90_err(nf90_close(ncid))
-end subroutine read_gridIDs
+end subroutine read_grid_id
subroutine concoutput_ncf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
@@ -1708,7 +1708,7 @@ subroutine particle_initialpos(itime,idate,itime_start,idate_start)
call nf90_err(nf90_close(ncid))
end subroutine particle_initialpos
-subroutine wrt_particles_initialpos(itime,istart,iend)
+subroutine wrt_part_initialpos(itime,istart,iend)
!*****************************************************************************
! *
@@ -1758,7 +1758,7 @@ subroutine wrt_particles_initialpos(itime,istart,iend)
call nf90_err(nf90_close(ncid))
partinitpointer = partinitpointer+newpart
-end subroutine wrt_particles_initialpos
+end subroutine wrt_part_initialpos
subroutine partinit_ncf(itime,field,fieldname,imass,ncid)
@@ -2353,7 +2353,7 @@ subroutine partoutput_ncf(itime,field,fieldname,imass,ncid)
! call nf90_err(nf90_close(ncid))
end subroutine partoutput_ncf
-subroutine readpartpositions_ncf(ibtime,ibdate)
+subroutine read_partpos_ncf(ibtime,ibdate)
!*****************************************************************************
! *
@@ -2466,7 +2466,7 @@ subroutine readpartpositions_ncf(ibtime,ibdate)
end do
call nf90_err(nf90_close(ncidend))
-end subroutine readpartpositions_ncf
+end subroutine read_partpos_ncf
subroutine read_init_cond_nc()
@@ -2652,7 +2652,7 @@ subroutine read_init_cond_nc()
endif
end do
write(*,FMT='(A,ES14.7)') ' Total mass to be released:', sum(xmass(1:numpoint,1:nspec))
- call get_total_part_num(numpart)
+ call get_totalpart_num(numpart)
numparticlecount=numpart
call nf90_err(nf90_close(ncidend))
diff --git a/src/output_mod.f90 b/src/output_mod.f90
index 60dd1daec8c0379a1e2777ae64b9958a3ae4cd58..3abcdd29ac7b2765ee796b7afcd916c6898c1603 100644
--- a/src/output_mod.f90
+++ b/src/output_mod.f90
@@ -23,7 +23,7 @@ module output_mod
character(len=256) :: restart_filename1,restart_filename2,restart_filename3
contains
-subroutine initialise_output(itime,filesize)
+subroutine init_output(itime,filesize)
implicit none
@@ -124,7 +124,7 @@ subroutine initialise_output(itime,filesize)
end do
#endif
endif
-end subroutine initialise_output
+end subroutine init_output
subroutine finalise_output(itime)
! Complete the calculation of initial conditions for particles not yet terminated
@@ -150,7 +150,7 @@ subroutine finalise_output(itime)
if (linit_cond.ge.1) then
if(linversionout.eq.1) then
- call initcond_output_inversion(itime) ! dump initial cond. field
+ call initcond_output_inv(itime) ! dump initial cond. field
else
call initcond_output(itime) ! dump initial cond. fielf
endif
@@ -260,7 +260,7 @@ subroutine output_restart(itime,loutnext,outnum)
if(stat == 0) close(1234, status='delete')
end subroutine output_restart
-subroutine output_heightlevels(height_tmp,nmixz_tmp)
+subroutine output_heightlevs(height_tmp,nmixz_tmp)
implicit none
real,intent(in) :: height_tmp(nzmax)
@@ -281,7 +281,7 @@ subroutine output_heightlevels(height_tmp,nmixz_tmp)
write(unitheightlevels) height_tmp(kz)
end do
close(unitheightlevels)
-end subroutine output_heightlevels
+end subroutine output_heightlevs
subroutine output_particle(itime,initial_output)
! i
@@ -353,7 +353,7 @@ subroutine output_particle(itime,initial_output)
!$OMP PARALLEL PRIVATE(i,j,m,tmp,ns,i_av,cartxyz_comp,cartxyz,np)
! Some variables needed for temporal interpolation
!*************************************************
- call find_time_variables(itime)
+ call find_time_vars(itime)
!$OMP DO
do i=1,numpart
@@ -372,7 +372,7 @@ subroutine output_particle(itime,initial_output)
! Where in the grid? Stereographic (ngrid<0) or nested (ngrid>0)
!***************************************************************
call find_ngrid(real(part(i)%xlon),real(part(i)%ylat))
- call determine_grid_coordinates(real(part(i)%xlon),real(part(i)%ylat))
+ call find_grid_indices(real(part(i)%xlon),real(part(i)%ylat))
call find_grid_distances(real(part(i)%xlon),real(part(i)%ylat))
! First set dz1out from interpol_mod to -1 so it only is calculated once per particle
!************************************************************************************
@@ -595,7 +595,7 @@ subroutine output_particle(itime,initial_output)
#endif
end subroutine output_particle
-subroutine output_concentrations(itime,loutstart,loutend,loutnext,outnum)
+subroutine output_conc(itime,loutstart,loutend,loutnext,outnum)
use unc_mod
use prepoutgrid_mod
use par_mod
@@ -719,7 +719,7 @@ subroutine output_concentrations(itime,loutstart,loutend,loutnext,outnum)
outnum=outnum+weight
call conccalc(itime,weight)
endif
-end subroutine output_concentrations
+end subroutine output_conc
subroutine conccalc(itime,weight)
! i i
@@ -1345,7 +1345,7 @@ subroutine conccalc(itime,weight)
end do
end subroutine conccalc
-subroutine partpos_average(itime,j)
+subroutine partpos_avg(itime,j)
!**********************************************************************
! This subroutine averages particle quantities, to be used for particle
@@ -1382,7 +1382,7 @@ subroutine partpos_average(itime,j)
! Some variables needed for temporal interpolation
!*************************************************
- call find_time_variables(itime)
+ call find_time_vars(itime)
xlon=xlon0+real(part(j)%xlon)*dx
ylat=ylat0+real(part(j)%ylat)*dy
@@ -1393,7 +1393,7 @@ subroutine partpos_average(itime,j)
! Where in the grid? Stereographic (ngrid<0) or nested (ngrid>0)
!***************************************************************
call find_ngrid(real(part(j)%xlon),real(part(j)%ylat))
- call determine_grid_coordinates(real(part(j)%xlon),real(part(j)%ylat))
+ call find_grid_indices(real(part(j)%xlon),real(part(j)%ylat))
call find_grid_distances(real(part(j)%xlon),real(part(j)%ylat))
! First set dz1out from interpol_mod to -1 so it only is calculated once per particle
@@ -1490,6 +1490,6 @@ subroutine partpos_average(itime,j)
cart_comp=.false.
return
-end subroutine partpos_average
+end subroutine partpos_avg
end module output_mod
diff --git a/src/particle_mod.f90 b/src/particle_mod.f90
index 048caecd54d0ff357358dfc6640f2b7d623c5e7d..0696063a3e9965d488b5214ef1757f99aa4099e0 100644
--- a/src/particle_mod.f90
+++ b/src/particle_mod.f90
@@ -115,9 +115,9 @@ module particle_mod
terminate_particle, &
spawn_particle, &
spawn_particles, &
- get_total_part_num, &
- get_alive_part_num, &
- get_new_part_index, &
+ get_totalpart_num, &
+ get_alivepart_num, &
+ get_newpart_index, &
particle_allocated, &
update_xlon, &
update_ylat, &
@@ -125,35 +125,35 @@ module particle_mod
count
interface update_xlon
- procedure update_xlon_dp, update_xlon_float, update_xlon_int
+ procedure update_xlon_dp, update_xlon_sp, update_xlon_int
end interface update_xlon
interface set_xlon
- procedure set_xlon_dp, set_xlon_float, set_xlon_int
+ procedure set_xlon_dp, set_xlon_sp, set_xlon_int
end interface set_xlon
interface update_ylat
- procedure update_ylat_dp, update_ylat_float, update_ylat_int
+ procedure update_ylat_dp, update_ylat_sp, update_ylat_int
end interface update_ylat
interface set_ylat
- procedure set_ylat_dp, set_ylat_float, set_ylat_int
+ procedure set_ylat_dp, set_ylat_sp, set_ylat_int
end interface set_ylat
interface update_z
- procedure update_z_dp,update_z_float
+ procedure update_z_dp,update_z_sp
end interface update_z
interface set_z
- procedure set_z_dp,set_z_float
+ procedure set_z_dp,set_z_sp
end interface set_z
interface update_zeta
- procedure update_zeta_dp,update_zeta_float
+ procedure update_zeta_dp,update_zeta_sp
end interface update_zeta
interface set_zeta
- procedure set_zeta_dp,set_zeta_float
+ procedure set_zeta_dp,set_zeta_sp
end interface set_zeta
contains
@@ -176,7 +176,7 @@ contains
end function particle_allocated
- subroutine get_new_part_index(ipart)
+ subroutine get_newpart_index(ipart)
!**************************************************
! Returns the first free spot to put a new particle
@@ -186,9 +186,9 @@ contains
ipart = count%spawned + 1
- end subroutine get_new_part_index
+ end subroutine get_newpart_index
- subroutine get_total_part_num(npart)
+ subroutine get_totalpart_num(npart)
!********************************************
! Returns total number of particles spawned *
@@ -198,9 +198,9 @@ contains
npart = count%spawned
- end subroutine get_total_part_num
+ end subroutine get_totalpart_num
- subroutine get_alive_part_num(npart)
+ subroutine get_alivepart_num(npart)
!**********************************************
! Returns number of particles currently alive *
@@ -210,7 +210,7 @@ contains
npart = count%alive
- end subroutine get_alive_part_num
+ end subroutine get_alivepart_num
subroutine spawn_particles(itime, nmpart)
@@ -423,7 +423,7 @@ contains
end subroutine update_xlon_dp
- subroutine update_xlon_float(ipart,xchange)
+ subroutine update_xlon_sp(ipart,xchange)
!**************************************
! Updates the longitude of the particle
@@ -434,7 +434,7 @@ contains
part(ipart)%xlon = part(ipart)%xlon + real(xchange,kind=dp)
- end subroutine update_xlon_float
+ end subroutine update_xlon_sp
subroutine update_xlon_int(ipart,xchange)
@@ -464,7 +464,7 @@ contains
end subroutine set_xlon_dp
- subroutine set_xlon_float(ipart,xvalue)
+ subroutine set_xlon_sp(ipart,xvalue)
!**************************************
! Sets the longitude of the particle
@@ -475,7 +475,7 @@ contains
part(ipart)%xlon = real(xvalue,kind=dp)
- end subroutine set_xlon_float
+ end subroutine set_xlon_sp
subroutine set_xlon_int(ipart,xvalue)
@@ -506,7 +506,7 @@ contains
end subroutine update_ylat_dp
- subroutine update_ylat_float(ipart,ychange)
+ subroutine update_ylat_sp(ipart,ychange)
!**************************************
! Updates the latitude of the particle
@@ -517,7 +517,7 @@ contains
part(ipart)%ylat = part(ipart)%ylat + real(ychange,kind=dp)
- end subroutine update_ylat_float
+ end subroutine update_ylat_sp
subroutine update_ylat_int(ipart,ychange)
@@ -547,7 +547,7 @@ contains
end subroutine set_ylat_dp
- subroutine set_ylat_float(ipart,yvalue)
+ subroutine set_ylat_sp(ipart,yvalue)
!**************************************
! Sets the latitude of the particle
@@ -558,7 +558,7 @@ contains
part(ipart)%ylat = real(yvalue,kind=dp)
- end subroutine set_ylat_float
+ end subroutine set_ylat_sp
subroutine set_ylat_int(ipart,yvalue)
@@ -590,7 +590,7 @@ contains
end subroutine update_z_dp
- subroutine update_z_float(ipart,zchange)
+ subroutine update_z_sp(ipart,zchange)
!**************************************
! Updates the height of the particle
@@ -603,7 +603,7 @@ contains
part(ipart)%meterupdate=.false.
part(ipart)%etaupdate=.true.
- end subroutine update_z_float
+ end subroutine update_z_sp
subroutine update_zeta_dp(ipart,zchange)
@@ -620,7 +620,7 @@ contains
end subroutine update_zeta_dp
- subroutine update_zeta_float(ipart,zchange)
+ subroutine update_zeta_sp(ipart,zchange)
!**************************************
! Updates the height of the particle
@@ -633,7 +633,7 @@ contains
part(ipart)%etaupdate=.false.
part(ipart)%meterupdate=.true.
- end subroutine update_zeta_float
+ end subroutine update_zeta_sp
! End update z positions
! Update z positions
@@ -652,7 +652,7 @@ contains
end subroutine set_z_dp
- subroutine set_z_float(ipart,zvalue)
+ subroutine set_z_sp(ipart,zvalue)
!**************************************
! Updates the height of the particle
@@ -665,7 +665,7 @@ contains
part(ipart)%meterupdate=.false.
part(ipart)%etaupdate=.true.
- end subroutine set_z_float
+ end subroutine set_z_sp
subroutine set_zeta_dp(ipart,zvalue)
@@ -682,7 +682,7 @@ contains
end subroutine set_zeta_dp
- subroutine set_zeta_float(ipart,zvalue)
+ subroutine set_zeta_sp(ipart,zvalue)
!**************************************
! Updates the height of the particle
@@ -695,7 +695,7 @@ contains
part(ipart)%etaupdate=.false.
part(ipart)%meterupdate=.true.
- end subroutine set_zeta_float
+ end subroutine set_zeta_sp
! End update z positions
end module particle_mod
diff --git a/src/timemanager_mod.f90 b/src/timemanager_mod.f90
index ead957b4bdcd21138590eb4bedb43543584d2efc..1c0b2fb8c50c70f72f7c1ccc09cd683f62fa01ed 100644
--- a/src/timemanager_mod.f90
+++ b/src/timemanager_mod.f90
@@ -202,7 +202,7 @@ subroutine timemanager
if ((itime.ne.itime_init).and.(mod(itime,loutrestart).eq.0)) call output_restart(itime,loutnext,outnum)
if (itime.ne.0) write(*,*) part(1)%xlon,part(1)%ylat,part(1)%z,part(1)%zeta
- call initialise_output(itime,filesize)
+ call init_output(itime,filesize)
! Get necessary wind fields if not available
!*******************************************
@@ -286,7 +286,7 @@ subroutine timemanager
!$OMP END PARALLEL
count%alive=alive_tmp
count%spawned=spawned_tmp
- call get_total_part_num(numpart)
+ call get_totalpart_num(numpart)
else
call releaseparticle(itime)
endif
@@ -301,7 +301,7 @@ subroutine timemanager
! If middle of averaging period of output fields is reached, accumulated
! deposited mass radioactively decays
!***********************************************************************
- if (DEP.and.(itime.eq.loutnext).and.(ldirect.gt.0)) call radioactive_decay() !OMP, unc_mod.f90 (needs test)
+ if (DEP.and.(itime.eq.loutnext).and.(ldirect.gt.0)) call deposit_decay() !OMP, unc_mod.f90 (needs test)
! Is the time within the computation interval, if not, skip
@@ -324,7 +324,7 @@ subroutine timemanager
endif
! Check whether concentrations are to be calculated and outputted
!****************************************************************
- call output_concentrations(itime,loutstart,loutend,loutnext,outnum)
+ call output_conc(itime,loutstart,loutend,loutnext,outnum)
call SYSTEM_CLOCK(count_clock, count_rate, count_max)
s_writepart = s_writepart + ((count_clock - count_clock0)/real(count_rate)-s_temp)
endif
@@ -385,7 +385,7 @@ subroutine timemanager
!***********************************
if ((part(j)%tstart.eq.itime).or.(itime.eq.0)) then
call update_zeta_to_z(itime, j)
- call initialize_particle(itime,j)
+ call init_particle(itime,j)
endif
! Memorize particle positions
@@ -423,7 +423,7 @@ subroutine timemanager
call advance(itime,j,thread)
if (part(j)%nstop.eqv..true.) cycle
- if (n_average.gt.0) call partpos_average(itime,j)
+ if (n_average.gt.0) call partpos_avg(itime,j)
! Calculate the gross fluxes across layer interfaces
!***************************************************
@@ -479,34 +479,43 @@ subroutine timemanager
terminated_tmp=terminated_tmp+1
else
- ! Dry deposition and radioactive decay for each species
- ! Also check maximum (of all species) of initial mass remaining on the particle;
- ! if it is below a threshold value, terminate particle
- !*****************************************************************************
+! Dry deposition and radioactive decay for each species
+! Also check maximum (of all species) of initial mass remaining on the particle;
+! if it is below a threshold value, terminate particle
+!*****************************************************************************
xmassfract=0.
do ks=1,nspec
- if (DRYDEPSPEC(ks)) then ! dry deposition (and radioactive decay)
+
+ if (DRYDEPSPEC(ks)) then ! dry deposition (and radioactive decay)
+
call drydepo_massloss(j,ks,ldeltat,drydeposit(ks))
- else if (decay(ks).gt.0.) then ! no dry deposition, but radioactive decay
- part(j)%mass(ks)=part(j)%mass(ks)*exp(-real(abs(lsynctime))*decay(ks))
+
+ else if (decay(ks).gt.0.) then ! no dry depo, but radioactive decay
+
+ part(j)%mass(ks) = part(j)%mass(ks) * &
+ exp( -real(abs(lsynctime)) * decay(ks) )
+
endif
+
! Skip check on mass fraction when npoint represents particle number
if (mdomainfill.eq.0.and.mquasilag.eq.0) then
- if ((ipin.eq.3).or.(ipin.eq.4)) then
- if (part(j)%mass_init(ks).gt.0) then
- xmassfract=max(xmassfract,part(j)%mass(ks)/part(j)%mass_init(ks))
- endif
+ if (ipin.eq.3 .or. ipin.eq.4) then
+ if (part(j)%mass_init(ks).gt.0) &
+ xmassfract = max( xmassfract, &
+ part(j)%mass(ks) / part(j)%mass_init(ks) )
else if (xmass(part(j)%npoint,ks).gt.0.) then
- xmassfract=max(xmassfract,real(npart(part(j)%npoint))* &
- part(j)%mass(ks)/xmass(part(j)%npoint,ks))
+ xmassfract = max( xmassfract, real( npart(part(j)%npoint) ) * &
+ part(j)%mass(ks) / xmass(part(j)%npoint,ks) )
endif
else
xmassfract=1.0
end if
+
end do
- if (xmassfract.le.minmassfrac) then ! terminate all particles carrying less mass
+ if (xmassfract.le.minmassfrac) then
+ ! terminate all particles carrying less mass
call terminate_particle(j,itime)
alive_tmp=alive_tmp-1
terminated_tmp=terminated_tmp+1
@@ -615,7 +624,7 @@ subroutine timemanager
call dealloc_windf
call dealloc_domainfill
call dealloc_drydepo
- call dealloc_convection
+ call dealloc_convect
call dealloc_getfields
call dealloc_interpol
call dealloc_random
diff --git a/src/turbulence_mod.f90 b/src/turbulence_mod.f90
index d63115ecd34ac685d1763c8287f1562ee5e897c3..d347dec2a53bcf43160a206639934484a7851838 100644
--- a/src/turbulence_mod.f90
+++ b/src/turbulence_mod.f90
@@ -109,7 +109,7 @@ subroutine turb_pbl(ipart,nrand,dt,zts,rhoa,rhograd,thread)
bth*rannumb(nrand)*sqrt(dtf))*icbt_r
delz=wp*dtf
if ((flagrein.eq.1).or.(wp.ne.wp).or.((wp-1.).eq.wp)) then
- call re_initialize_particle(zts,ust,wst,h,sigw,old_wp_buf,nrand,ol)
+ call reinit_particle(zts,ust,wst,h,sigw,old_wp_buf,nrand,ol)
wp=old_wp_buf
delz=wp*dtf
nan_count(thread+1)=nan_count(thread+1)+1
diff --git a/src/unc_mod.f90 b/src/unc_mod.f90
index 312117c99a37d57ccc7fcf1bf8a4b342777fbc64..e536ef7f9127f6807343bc32121b7e98edcb1deb 100644
--- a/src/unc_mod.f90
+++ b/src/unc_mod.f90
@@ -42,7 +42,7 @@ module unc_mod
contains
-subroutine radioactive_decay()
+subroutine deposit_decay()
! Accumulated deposited mass radioactively decays
use com_mod
@@ -94,6 +94,6 @@ subroutine radioactive_decay()
end do
!$OMP END DO
!$OMP END PARALLEL
-end subroutine radioactive_decay
+end subroutine deposit_decay
end module unc_mod
diff --git a/src/verttransform_mod.f90 b/src/verttransform_mod.f90
index bd8de09961aa11ff196fbc72112e1f14021b0910..9bd5b70a8e04b48c1010f970d9f65b3f916ac3e3 100644
--- a/src/verttransform_mod.f90
+++ b/src/verttransform_mod.f90
@@ -110,7 +110,7 @@ subroutine verttransform_ecmwf(n,uuh,vvh,wwh,pvh)
! Search for a point with high surface pressure (i.e. not above significant topography)
! Then, use this point to construct a reference z profile, to be used at all times
!*****************************************************************************
- call verttransform_init(n)
+ call vertransf_init(n)
! Do not repeat initialization of the Cartesian z grid
!*****************************************************
@@ -240,7 +240,7 @@ subroutine verttransform_nest(n,uuhn,vvhn,wwhn,pvhn)
end subroutine verttransform_nest
-subroutine verttransform_init(n)
+subroutine vertransf_init(n)
use par_mod
use com_mod
@@ -300,8 +300,8 @@ subroutine verttransform_init(n)
endif
end do
- call output_heightlevels(height,nmixz)
-end subroutine verttransform_init
+ call output_heightlevs(height,nmixz)
+end subroutine vertransf_init
subroutine verttransf_ec_windf(n,uuh,vvh,wwh,pvh,rhoh,prsh,pinmconv)
use par_mod
@@ -1125,7 +1125,7 @@ subroutine verttransform_gfs(n,uuh,vvh,wwh,pvh)
! Search for a point with high surface pressure (i.e. not above significant topography)
! Then, use this point to construct a reference z profile, to be used at all times
!*****************************************************************************
- call verttransform_init(n)
+ call vertransf_init(n)
! Do not repeat initialization of the Cartesian z grid
!*****************************************************
diff --git a/src/wetdepo_mod.f90 b/src/wetdepo_mod.f90
index a74abc4b68fda3b891f5b55d4d2440564bbb4a9a..ab1e0fa168a27bd2b61c49ed578888958b1abcbe 100644
--- a/src/wetdepo_mod.f90
+++ b/src/wetdepo_mod.f90
@@ -307,7 +307,7 @@ subroutine get_wetscav(itime,ltsample,loutnext,jpart,ks,grfraction,inc_count,blc
if (yts.ge.real(ny-1)) yts=real(ny-1)-0.00001
endif
- call determine_grid_coordinates(xts,yts)
+ call find_grid_indices(xts,yts)
call find_grid_distances(xts,yts)
if (ngrid.le.0) then