diagfld Subroutine

public subroutine diagfld()

Uses

  • proc~~diagfld~~UsesGraph proc~diagfld diagfld module~modfields modfields proc~diagfld->module~modfields module~modglobal modglobal proc~diagfld->module~modglobal module~modmpi modmpi proc~diagfld->module~modmpi module~modsurfdata modsurfdata proc~diagfld->module~modsurfdata decomp_2d decomp_2d module~modfields->decomp_2d mpi mpi module~modmpi->mpi

Arguments

None

Calls

proc~~diagfld~~CallsGraph proc~diagfld diagfld proc~avexy_ibm avexy_ibm proc~diagfld->proc~avexy_ibm proc~fromztop fromztop proc~diagfld->proc~fromztop mpi_allreduce mpi_allreduce proc~avexy_ibm->mpi_allreduce

Called by

proc~~diagfld~~CalledByGraph proc~diagfld diagfld proc~thermodynamics thermodynamics proc~thermodynamics->proc~diagfld proc~readinitfiles readinitfiles proc~readinitfiles->proc~thermodynamics program~dalesurban DALESURBAN program~dalesurban->proc~thermodynamics program~dalesurban->proc~readinitfiles

Source Code

  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