subroutine initglobal
use modmpi, only:nprocs, myid, comm3d, my_real, mpierr
implicit none
integer :: advarr(4)
real phi, colat, silat, omega, omega_gs
integer :: i, k, n
character(80) chmess
!timestepping
if (courant < 0) then
select case (iadv_mom)
case (iadv_cd2)
courant = 1.5
case default
courant = 1.4
end select
if (any(iadv_sv(1:nsv) == iadv_kappa) .or. any((/iadv_thl, iadv_qt, iadv_tke/) == iadv_kappa)) then
courant = min(courant, 1.1)
elseif (any(iadv_sv(1:nsv) == iadv_upw) .or. any((/iadv_thl, iadv_qt, iadv_tke/) == iadv_upw)) then
courant = min(courant, 1.1)
elseif (any(iadv_sv(1:nsv) == iadv_cd2) .or. any((/iadv_thl, iadv_qt, iadv_tke/) == iadv_cd2)) then
courant = min(courant, 1.5)
end if
end if
! phsgrid
jmax = jtot/nprocs
isen = imax/nprocs
jsen = jmax
!set the number of ghost cells. NB: This switch has to run in order of required ghost cells
advarr = (/iadv_mom, iadv_tke, iadv_thl, iadv_qt/)
if (any(advarr == iadv_kappa)) then
ih = 2
jh = 2
! ih = 1
! jh = 1
kh = 1
elseif (any(advarr == iadv_cd2) .or. any(iadv_sv(1:nsv) == iadv_cd2)) then
ih = 1
jh = 1
kh = 1
ihc = 1
jhc = 1
khc = 1
end if
! J. Tomas added this for using only kappa scheme for sv(:)
if (any(iadv_sv(1:nsv) == iadv_kappa)) then
ihc = 2
jhc = 2
khc = 2
end if
offset = 1
ib = 2 - offset
ie = imax + 1 - offset
jb = 2 - offset
je = jmax + 1 - offset
jgb = jb ! global j range (starting at the same as j as the processor j range)
jge = jtot + 1 - offset ! global j range
kb = 1 - offset
ke = kmax - offset
! Global constants
! Select advection scheme for scalars. If not set in the options file, the momentum scheme is used
if (iadv_tke < 0) iadv_tke = iadv_mom
if (iadv_thl < 0) iadv_thl = iadv_mom
if (iadv_qt < 0) iadv_qt = iadv_mom
!CvH remove where
!where (iadv_sv<0) iadv_sv = iadv_mom
!tg3315 added - only uses kappa advection scheme...
do n = 1, nsv
iadv_sv(n) = iadv_kappa
end do
!ends here
phi = xlat*pi/180.
colat = cos(phi)
silat = sin(phi)
omega = 7.292e-5
omega_gs = 7.292e-5
om22 = 2.*omega*colat
om23 = 2.*omega*silat
om22_gs = 2.*omega_gs*colat
om23_gs = 2.*omega_gs*silat
! Variables
allocate (dsv(nsv))
write (cexpnr, '(i3.3)') iexpnr
! Create the physical grid variables
allocate (dzf(kb - kh:ke + kh))
allocate (dzf2(kb - kh:ke + kh))
allocate (dzfi(kb - kh:ke + kh))
allocate (dzfiq(kb - kh:ke + kh))
allocate (dzfi5(kb - kh:ke + kh))
allocate (dzh(kb:ke + kh))
allocate (dzhi(kb:ke + kh))
allocate (dzhiq(kb:ke + kh))
allocate (dzh2i(kb:ke + kh))
allocate (zh(kb:ke + kh))
allocate (zf(kb:ke + kh))
allocate (dxf(ib - ih:ie + ih))
allocate (dxf2(ib - ih:ie + ih))
allocate (dxfi(ib - ih:ie + ih))
allocate (dxfiq(ib - ih:ie + ih))
allocate (dxfi5(ib - ih:ie + ih))
allocate (dxh(ib:ie + ih))
allocate (dxhi(ib:ie + ih))
allocate (dxhiq(ib:ie + ih))
allocate (dxh2i(ib:ie + ih))
allocate (xh(ib:ie + ih))
allocate (xf(ib:ie + ih))
allocate (delta(ib - ih:ie + ih, kb:ke + kh))
rslabs = real(imax*jtot)
dy = ysize/float(jtot)
! MPI
! Note, that the loop for reading zf and calculating zh
! has been split so that reading is only done on PE 1
if (myid == 0) then
open (ifinput, file='prof.inp.'//cexpnr)
read (ifinput, '(a72)') chmess
read (ifinput, '(a72)') chmess
do k = kb, ke
read (ifinput, *) zf(k)
end do
close (ifinput)
! J. Tomas: Read the x-coordinates of the cell centers from xgrid.inp.XXX
open (ifinput, file='xgrid.inp.'//cexpnr)
read (ifinput, '(a72)') chmess
read (ifinput, '(a72)') chmess
do i = ib, ie
read (ifinput, *) xf(i)
end do
close (ifinput)
end if ! end if myid==0
! MPI broadcast kmax elements from zf
call MPI_BCAST(zf, kmax, MY_REAL, 0, comm3d, mpierr)
! MPI broadcast imax elements from xf
call MPI_BCAST(xf, imax, MY_REAL, 0, comm3d, mpierr)
zh(kb) = 0.0
do k = kb, ke
zh(k + 1) = zh(k) + 2.0*(zf(k) - zh(k))
end do
zf(ke + kh) = zf(ke) + 2.0*(zh(ke + kh) - zf(ke))
do k = kb, ke
dzf(k) = zh(k + 1) - zh(k)
end do
dzf(ke + 1) = dzf(ke)
dzf(kb - 1) = dzf(kb)
dzh(kb) = 2*zf(kb)
do k = kb + 1, ke + kh
dzh(k) = zf(k) - zf(k - 1)
end do
! j. tomas: same trick for x-direction...
xh(ib) = 0.0
do i = ib, ie
xh(i + 1) = xh(i) + 2.0*(xf(i) - xh(i))
end do
xf(ie + ih) = xf(ie) + 2.0*(xh(ie + ih) - xf(ie))
do i = ib, ie
dxf(i) = xh(i + 1) - xh(i)
end do
dxf(ie + 1) = dxf(ie)
dxf(ib - 1) = dxf(ib)
dxh(ib) = 2*xf(ib)
do i = ib + 1, ie + ih
dxh(i) = xf(i) - xf(i - 1)
end do
do k = kb, ke + kh
do i = ib - ih, ie + ih
delta(i, k) = (dxf(i)*dy*dzf(k))**(1./3.)
end do
end do
!--------------------------------------------------
! *** Check whether the grid is equidistant *****
!--------------------------------------------------
!if (myid == 0) then
!do k=kb,ke+kh
!if (.not.(dzf(k).eq.dzf(1)))
! write (6, *) &
! 'WARNING, You are working with a non-equidistant grid!!!!'
!end if
!end do
!end if ! end if myid==0
dzhi = 1./dzh
dzfi = 1./dzf
dzf2 = dzf*dzf
dxhi = 1./dxh
dxfi = 1./dxf
dxf2 = dxf*dxf
dyi = 1./dy
dy2 = dy*dy
dzhiq = 0.25*dzhi
dzfiq = 0.25*dzfi
dxhiq = 0.25*dxhi
dxfiq = 0.25*dxfi
dyiq = 0.25*dyi
dzh2i = dzhi*dzhi
dxh2i = dxhi*dxhi
dy2i = dyi*dyi
dzfi5 = 0.5*dzfi
dxfi5 = 0.5*dxfi
dyi5 = 0.5*dyi
! Grid used in kappa scheme advection (extra ghost nodes)
if (any(iadv_sv(1:nsv) == iadv_kappa)) then
allocate (dzfc(kb - khc:ke + khc))
allocate (dxfc(ib - ihc:ie + ihc))
allocate (dzfci(kb - khc:ke + khc))
allocate (dxfci(ib - ihc:ie + ihc))
allocate (dzhci(kb - 1:ke + khc))
allocate (dxhci(ib - 1:ie + ihc))
dzfc(kb - kh:ke + kh) = dzf(kb - kh:ke + kh)
dzfc(kb - khc) = dzfc(kb - kh)
dzfc(ke + khc) = dzfc(ke + kh)
dxfc(ib - ih:ie + ih) = dxf(ib - ih:ie + ih)
dxfc(ib - ihc) = dxfc(ib - ih)
dxfc(ie + ihc) = dxfc(ie + ih)
dzhci(kb:ke + kh) = dzhi(kb:ke + kh)
dzhci(kb - 1) = dzhci(kb)
dzhci(ke + khc) = dzhci(ke + kh)
dxhci(ib:ie + ih) = dxhi(ib:ie + ih)
dxhci(ib - 1) = dxhci(ib)
dxhci(ie + ihc) = dxhci(ie + ih)
dzfci = 1./dzfc
dxfci = 1./dxfc
end if
if (myid == 0) then
write (6, *) 'lev dz zf zh dzh delta(ib,k)'
do k = ke + 1, kb, -1
write (6, '(i4,5f8.5)') k, dzf(k), zf(k), zh(k), dzh(k), delta(ib, k)
end do
! same for x:
write (6, *) 'lev dxf xf xh dxh delta(i,kb)'
do i = ie + 1, ib, -1
write (6, '(i4,5f9.5)') i, dxf(i), xf(i), xh(i), dxh(i), delta(i, kb)
end do
end if
tnextrestart = trestart
tnextfielddump = tfielddump
! tnextstatsdump = tstatsdump
timeleft = runtime ! tg3315 previously btime + runtime
end subroutine initglobal