subroutine diagfld
use modglobal, only : ib,ie,ih,jb,je,jh,kb,ke,kh,khc,nsv,zh,zf,rslabs,grav,rlv,cp,rd,rv,pref0
use modfields, only : u0,v0,thl0,qt0,ql0,sv0,u0av,v0av,thl0av,qt0av,ql0av,sv0av, &
presf,presh,exnf,exnh,rhof,thvf,IIc,IIcs,IIu,IIus,IIv,IIvs
use modsurfdata,only : thls,ps
use modmpi, only : slabsum,myid,avexy_ibm
implicit none
real tv
integer k, n
!!*********************************************************
!! 1.0 calculate average profiles of u,v,thl,qt and ql *
!! assuming hydrostatic equilibrium *
!!*********************************************************
! initialise local MPI arrays
u0av = 0.0
v0av = 0.0
thl0av = 0.0
th0av = 0.0
qt0av = 0.0
ql0av = 0.0
sv0av = 0.
!CvH changed momentum array dimensions to same value as scalars!
! call slabsum(u0av ,kb,ke+kh,u0 ,ib-ih,ie+ih,jb-jh,je+jh,kb,ke+kh,ib,ie,jb,je,kb,ke+kh)
! call slabsum(u0av ,kb,ke+kh,u0(:,:,kb:ke+kh) ,ib-ih,ie+ih,jb-jh,je+jh,kb,ke+kh,ib,ie,jb,je,kb,ke+kh)
call avexy_ibm(u0av(kb:ke+kh),u0(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIu(ib:ie,jb:je,kb:ke+kh),IIus(kb:ke+kh),.false.)
! call slabsum(v0av ,kb,ke+kh,v0 ,ib-ih,ie+ih,jb-jh,je+jh,kb,ke+kh,ib,ie,jb,je,kb,ke+kh)
! call slabsum(v0av ,kb,ke+kh,v0(:,:,kb:ke+kh) ,ib-ih,ie+ih,jb-jh,je+jh,kb,ke+kh,ib,ie,jb,je,kb,ke+kh)
call avexy_ibm(v0av(kb:ke+kh),v0(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIv(ib:ie,jb:je,kb:ke+kh),IIvs(kb:ke+kh),.false.)
! call slabsum(thl0av,kb,ke+kh,thl0,ib-ih,ie+ih,jb-jh,je+jh,kb,ke+kh,ib,ie,jb,je,kb,ke+kh)
! call slabsum(thl0av,kb,ke+kh,thl0(:,:,kb:ke+kh),ib-ih,ie+ih,jb-jh,je+jh,kb,ke+kh,ib,ie,jb,je,kb,ke+kh)
call avexy_ibm(thl0av(kb:ke+kh),thl0(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc(ib:ie,jb:je,kb:ke+kh),IIcs(kb:ke+kh),.false.)
!write(*,*) 'thl0av(kb), thl0av(kb+1)', thl0av(kb), thl0av(kb+1)
!if (IIbl == 0) then ! as lEB applies blocks to kb and masking matrices average this to zero
! thl0av(kb) = thl0av(kb+1)
!end if
! call slabsum(qt0av ,kb,ke+kh,qt0 ,ib-ih,ie+ih,jb-jh,je+jh,kb,ke+kh,ib,ie,jb,je,kb,ke+kh)
! call slabsum(qt0av ,kb,ke+kh,qt0(:,:,kb:ke+kh) ,ib-ih,ie+ih,jb-jh,je+jh,kb,ke+kh,ib,ie,jb,je,kb,ke+kh)
call avexy_ibm(qt0av(kb:ke+kh),qt0(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc(ib:ie,jb:je,kb:ke+kh),IIcs(kb:ke+kh),.false.)
! call slabsum(ql0av ,kb,ke+kh,ql0 ,ib-ih,ie+ih,jb-jh,je+jh,kb,ke+kh,ib,ie,jb,je,kb,ke+kh)
call avexy_ibm(ql0av(kb:ke+kh),ql0(ib:ie,jb:je,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc(ib:ie,jb:je,kb:ke+kh),IIcs(kb:ke+kh),.false.)
exnf = 1-grav*zf/(cp*thls)
exnh = 1-grav*zh/(cp*thls)
th0av = thl0av + (rlv/cp)*ql0av/exnf
!write(*,*) 'thl0av',thl0av
!write(*,*) 'th0av',th0av
do n=1,nsv
! call slabsum(sv0av(kb,n),kb,ke+kh,sv0(ib-ih,jb-jh,kb,n),ib-ih,ie+ih,jb-jh,je+jh,kb,ke+kh,ib,ie,jb,je,kb,ke+kh)
call avexy_ibm(sv0av(kb:ke+khc,n),sv0(ib:ie,jb:je,kb:ke+khc,n),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc(ib:ie,jb:je,kb:ke+khc),IIcs(kb:ke+khc),.false.)
end do
! sv0av = sv0av/rslabs
!***********************************************************
! 2.0 calculate average profile of pressure at full and *
! half levels, assuming hydrostatic equilibrium. *
!***********************************************************
! 2.1 Use first guess of theta, then recalculate theta
call fromztop
exnf = (presf/pref0)**(rd/cp)
th0av = thl0av + (rlv/cp)*ql0av/exnf
! 2.2 Use new updated value of theta for determination of pressure
call fromztop
!***********************************************************
! 3.0 Construct density profiles and exner function *
! for further use in the program *
!***********************************************************
! 3.1 determine exner
exnh(kb) = (ps/pref0)**(rd/cp)
exnf(kb) = (presf(kb)/pref0)**(rd/cp)
do k=kb+1,ke+kh
exnf(k) = (presf(k)/pref0)**(rd/cp)
exnh(k) = (presh(k)/pref0)**(rd/cp)
end do
thvf(kb) = th0av(kb)*exnf(kb)*(1+(rv/rd-1)*qt0av(kb)-rv/rd*ql0av(kb))
rhof(kb) = presf(kb)/(rd*thvf(kb))
! 3.2 determine rho
do k=kb+1,ke !+kh ?
! write(*,*) "exnf(k)",exnf(k)
! write(*,*) "th0av(k)",th0av(k)
! write(*,*) "qt0av(k)",qt0av(k)
! write(*,*) "ql0av(k)",ql0av(k)
thvf(k) = th0av(k)*exnf(k)*(1.+(rv/rd-1)*qt0av(k)-rv/rd*ql0av(k))
rhof(k) = presf(k)/(rd*thvf(k))
end do
return
end subroutine diagfld