!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine idea_phys(im,ix,levs,prsi,prsl, &
& adu,adv,adt,adr,ntrac,dtp,lat, &
& solhr,slag,sdec,cdec,sinlat,coslat, &
& xlon,xlat,oro,cozen,swh,hlw,dt6dt, &
& thermodyn_id,sfcpress_id,gen_coord_hybrid,me,&
& mpi_ior,mpi_comm)
!-----------------------------------------------------------------------
! add temp, wind changes due to viscosity and thermal conductivity
! also solar heating
! Apr 06 2012 Henry Juang, initial implement for NEMS
! Jul 26 2012 Jun Wang, add mpi info
! Sep 06 2012 Jun Wang, add changing pressure to cb
! Dec 2012 Jun Wang, change to new rad_merge (from Rashid and Fei)
! May 2013 Jun Wang, tmp updated after rad_merge
! Jun 2013 S. Moorthi Some optimization and cosmetic changes
! Oct 2013 Henry Juang, correct the sequence to get prsi from model top
!-----------------------------------------------------------------------
use physcons, amo2=>con_amo2,avgd => con_avgd
use idea_composition
!
implicit none
! Argument
integer, intent(in) :: im ! number of data points in adt (first dim)
integer, intent(in) :: ix ! max data points in adt (first dim)
integer, intent(in) :: levs ! number of pressure levels
integer, intent(in) :: lat ! latitude index
integer, intent(in) :: ntrac ! number of tracer
integer, intent(in) :: me ! my pe
integer, intent(in) :: mpi_ior ! mpi real for io
integer, intent(in) :: mpi_comm ! mpi communicator
!
real, intent(in) :: dtp ! time step in second
real, intent(inout) :: prsi(ix,levs+1) ! pressure
real, intent(inout) :: prsl(ix,levs) ! pressure
real, intent(in) :: hlw(ix,levs) ! long wave rad (K/s)
real, intent(in) :: swh(ix,levs) ! short wave rad (K/s)
real, intent(in) :: cozen(im) ! time avg(1 hour) cos zenith angle
real, intent(in) :: oro(im) ! surface height (m)
real, intent(in) :: solhr,slag,sdec,cdec ! for solar zenith angle
real, intent(in) :: xlon(im),xlat(im),coslat(im),sinlat(im)
real, intent(inout) :: adr(ix,levs,ntrac) ! tracer
real, intent(inout) :: adt(ix,levs) ! temperature
real, intent(inout) :: adu(ix,levs) ! real u
real, intent(inout) :: adv(ix,levs) ! real v
real, intent(inout) :: dt6dt(ix,levs,6)! diagnostic array
integer,intent(in) :: thermodyn_id, sfcpress_id
logical,intent(in) :: gen_coord_hybrid
! Local variables
real,parameter :: pa2cb=0.001,cb2pa=1000.
!
real cp(ix,levs),cospass(im),dt(ix,levs),rtime1,hold1,n(ix,levs)
real o_n(ix,levs),o2_n(ix,levs),n2_n(ix,levs),o3_n(ix,levs), &
& am(ix,levs),dudt(ix,levs),dvdt(ix,levs),dtdt(ix,levs),xmu(im), &
& dtco2c(ix,levs),dtco2h(ix,levs) &
&,dth2oh(ix,levs),dth2oc(ix,levs),dto3(ix,levs),rho(ix,levs) &
&,wtot(ix,levs),zg(ix,levs) &
&,amin,amax,grav(ix,levs) &
&,prslk(ix,levs),prsik(ix,levs+1),phil(ix,levs),phii(ix,levs+1)
! solar
real utsec,sda,maglat(im),maglon(im),btot(im), &
& dipang(im),essa(im),dlat,dlon
integer i,k,dayno,j1,j2
! change to real windl !hmhj already real wind
!hmhj do i=1,im
!hmhj adu(i,1:levs)=adu(i,1:levs)/coslat(i)
!hmhj adv(i,1:levs)=adv(i,1:levs)/coslat(i)
!hmhj enddo ! i
! get phil geopotential from temperature
! change prsi and prsl to centibar from pascal
do k=1,levs
do i=1,im
prsi(i,k) = prsi(i,k)*pa2cb
prsl(i,k) = prsl(i,k)*pa2cb
enddo
enddo
do i=1,im
prsi(i,levs+1) = prsi(i,levs+1)*pa2cb
enddo
!hmhj call GET_PHI_gc_h(im,ix,levs,ntrac,adt,adr,prsi,phii,phil)
call get_phi(im,ix,levs,ntrac,adt,adr, &
& thermodyn_id, sfcpress_id, &
& gen_coord_hybrid, &
& prsi,prsik,prsl,prslk,phii,phil)
! get height
call phi2z(im,ix,levs,phil,oro,zg,grav)
! print*,'wwwz',zg(1,1:150)
! print*,'wwwg',grav(1,1:150)
! print*,'wwwp',phil(1,1:150),oro(1)
!
! get composition at layers (/cm3) and rho (kg/m3)
call idea_tracer(im,ix,levs,ntrac,2,grav,prsi,prsl,adt,adr, &
& dtp,o_n,o2_n,n2_n,n,rho,am)
! calculate cp
call getcp_idea(im,ix,levs,ntrac,adr,cp, &
& thermodyn_id,gen_coord_hybrid)
! dissipation
call idea_phys_dissipation(im,ix,levs,grav,prsi,prsl, &
& adu,adv,adt,o_n,o2_n,n2_n,dtp,cp,dt6dt)
!
! get cos solar zenith angle (instant)
call presolar(im,ix,solhr,slag, &
& sinlat,coslat,sdec,cdec,xlon,xlat &
& ,cospass,dayno,utsec,sda &
& ,maglat,maglon,btot,dipang,essa)
! get solar heating and NO cooling then get temp adjustment
call idea_sheat(im,ix,levs,adt,dt,cospass,o_n,o2_n,n2_n,rho, &
& cp,lat,dayno,prsl,zg,grav,am,maglat,dt6dt)
! rtime1=3600.*6.
do k=1,levs
do i=1,im
adt(i,k) = adt(i,k) + dt(i,k)*dtp
! dt3dt(i,k,1) = dt(i,k)*rtime1
!
! ion_drag - change to /m3
o_n(i,k) = o_n(i,k) * 1.e6
o2_n(i,k) = o2_n(i,k) * 1.e6
n2_n(i,k) = n2_n(i,k) * 1.e6
n(i,k) = n(i,k) * 1.e6
enddo
enddo
call idea_ion(solhr,cospass,zg,o_n,o2_n,n2_n,cp, &
& adu,adv,adt,dudt,dvdt,dtdt,rho,xlat,xlon,ix,im,levs,&
& dayno,utsec,sda,maglon,maglat,btot,dipang,essa)
! do i=1,im
! dlat=xlat(i)*180./3.14159
! dlon=xlon(i)*180./3.14159
! if(abs(dlat-60.).le.1..and.abs(dlon-270.).le.1.) then
! print*,'www0',solhr,dudt(i,140)*dtp,dvdt(i,140)*dtp, &
! &dtdt(i,140)*dtp,adu(i,140),adv(i,140)
! endif
do k=1,levs
do i=1,im
adu(i,k) = adu(i,k) + dtp*dudt(i,k)
adv(i,k) = adv(i,k) + dtp*dvdt(i,k)
adt(i,k) = adt(i,k) + dtp*dtdt(i,k)
enddo
enddo
! change u,V back !hmhj no need to change back, they are real wind
!hmhj do i=1,im
!hmhj adu(i,1:levs)=adu(i,1:levs)*coslat(i)
!hmhj adv(i,1:levs)=adv(i,1:levs)*coslat(i)
!hmhj enddo ! i
! radiation
! co2 cooling, heating
call idea_co2(im,ix,levs,nlev_co2,ntrac,grav,cp,adr,adt, &
& dtco2c,cospass,dtco2h)
!hmhj&'/mtb/save/wx20fw/fcst07rd',dtco2c,cospass,dtco2h)
! h2o cooling heating 110-41 down ward
call idea_h2o(im,ix,levs,nlev_h2o,nlevc_h2o,ntrac,grav,cp, &
& adr,adt,dth2oh,cospass,dth2oc)
dt6dt(1:im,1:levs,4) = dtco2c
! dt6dt(1:im,1:levs,5) = dth2oc
! dt6dt(1:im,1:levs,6) = dth2oh
! dth2oc=0.
! dth2oh=0.
! o2 o3 heating
call o3pro(im,ix,levs,ntrac,adr,am,n,o3_n)
call idea_o2_o3(im,ix,levs,cospass,adt,o2_n,o3_n,rho,cp, &
& zg,grav,dto3)
! get xmu
do i=1,im
if(cospass(i) > 0.0001 .and. cozen(i) > 0.0001) then
xmu(i) = cospass(i)/cozen(i)
else
xmu(i) = 0.
endif
enddo
! dt6dt
! do i=1,im
! do k=1,levs
! dt6dt(i,k,2)=dtco2c(i,k)
! dt6dt(i,k,3)=dtco2h(i,k)
! dt6dt(i,k,4)=dth2oc(i,k)
! dt6dt(i,k,5)=dth2oh(i,k)
! dt6dt(i,k,6)=dto3(i,k)
! enddo
! enddo
! merge
call rad_merge(im,ix,levs,hlw,swh,prsi,prsl,wtot,
& xmu,dtco2c,dtco2h,dth2oh,dth2oc,dto3,dt6dt)
do k=1,levs
do i=1,im
adt(i,k) = adt(i,k) + dtp*wtot(i,k)
! dt6dt
dt6dt(i,k,2) = wtot(i,k)
! change prsi and prsl back to pascal
prsi(i,k) = prsi(i,k)*cb2pa
prsl(i,k) = prsl(i,k)*cb2pa
enddo
enddo
do i=1,im
prsi(i,levs+1) = prsi(i,levs+1)*cb2pa
enddo
return
end subroutine
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
subroutine getcp_idea(im,ix,levs,ntrac,adr,xcp, &
& thermodyn_id,gen_coord_hybrid)
!
use tracer_const
!hmhj use resol_def , only: thermodyn_id
!hmhj use namelist_def , only: gen_coord_hybrid
!
implicit none
integer, intent(in) :: im ! number of data points in adr (first dim)
integer, intent(in) :: ix ! max data points in adr (first dim)
integer, intent(in) :: levs ! number of pressure levels
integer, intent(in) :: ntrac ! number of tracer
real, intent(in) :: adr(ix,levs,ntrac) ! tracer kg/kg
real, intent(out) :: xcp(ix,levs) !CP (J/kg/k)
integer thermodyn_id
logical gen_coord_hybrid
!
! local
real sumq(ix,levs),work1
integer i,j,k,ntb
sumq = 0.0
xcp = 0.0
if( gen_coord_hybrid .and. thermodyn_id == 3 ) then
ntb = 1
elseif (ntrac >= 4) then
ntb = 4
else
return
endif
do i=ntb,ntrac
if( cpi(i) /= 0.0 ) then
do k=1,levs
do j=1,im
work1 = adr(j,k,i)
sumq(j,k) = sumq(j,k) + work1
xcp(j,k) = xcp(j,k) + work1*cpi(i)
enddo
enddo
endif
enddo
do k=1,levs
do j=1,im
xcp(j,k) = xcp(j,k) + (1.-sumq(j,k))*cpi(0)
enddo
enddo
return
end
subroutine rad_merge(im,ix,levs,hlw,swh,prsi,prsl,wtot, &
& xmu,dtco2c,dtco2h,dth2oh,dth2oc,dto3,dt6dt)
!
implicit none
integer, intent(in) :: im ! number of data points in hlw,dt..(first dim)
integer, intent(in) :: ix ! max data points in hlw,... (first dim)
integer, intent(in) :: levs ! number of pressure levels
real, parameter :: xb=7.5, xt=8.5, rdx=1./(xt-xb)
real, intent(in) :: hlw(ix,levs) ! GFS lw rad (K/s)
real, intent(in) :: swh(ix,levs) ! GFS sw rad (K/s)
real, intent(in) :: prsi(ix,levs+1) ! pressure
real, intent(in) :: prsl(ix,levs) ! pressure
real, intent(in) :: xmu(im) ! mormalized cos zenith angle
real, intent(in) :: dtco2c(ix,levs) ! idea co2 cooling(K/s)
real, intent(in) :: dtco2h(ix,levs) ! idea co2 heating(K/s)
real, intent(in) :: dth2oc(ix,levs) ! idea h2o cooling(K/s)
real, intent(in) :: dth2oh(ix,levs) ! idea h2o heating(K/s)
real, intent(in) :: dto3(ix,levs) ! idea o3 heating(K/s)
real, intent(out) :: wtot(ix,levs) ! GFS idea combined rad
real, intent(inout) :: dt6dt(ix,levs,6)
! local
real xk,wl,wh
integer i,k,j
!
do k=1,levs
do i=1,im
xk = log(prsi(i,1)/prsl(i,k))
wh = dtco2c(i,k)+dth2oc(i,k)+dtco2h(i,k)+dth2oh(i,k)+dto3(i,k)
wl = hlw(i,k)+swh(i,k)*xmu(i)
if(xk < xb) then
wtot(i,k) = wl
elseif(xk >= xb .and. xk <= xt) then
wtot(i,k) = (wl*(xt-xk) + wh*(xk-xb))*rdx
else
wtot(i,k) = wh
endif
enddo
enddo
return
end
!
subroutine getmax(ain,n1,n,m,rmin,j1,rmax,j2)
real ain(n1,m)
rmin = 1.e36
rmax = -1.e36
i1 = 500
j1 = 500
i2 = 500
j2 = 500
do j=1,m
do i=1,n
if(rmin > ain(i,j)) then
rmin = ain(i,j)
i1 = i
j1 = j
endif
if(rmax < ain(i,j)) then
rmax = ain(i,j)
i2 = i
j2 = j
endif
enddo
enddo
return
end
subroutine getmax2(ain,ain1,n1,n,m,rmax,j2)
real ain(n1,m),ain1(n1,m)
rmax = -1.e36
i1 = 500
j1 = 500
i2 = 500
j2 = 500
do j=1,m
do i=1,n
sq = sqrt(ain(i,j)*ain(i,j) + ain1(i,j)*ain1(i,j))
if(rmax < sq) then
rmax = sq
i2 = i
j2 = j
endif
enddo
enddo
return
end
subroutine phi2z(im,ix,levs,phi,soro,z,grav)
! Subroutine to calculate geometric height and gravity from geopotential
! in a hydrostatic atmosphere, assuming a spherically symmetric planet
! and Newton's gravity.
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! File history
! Feb 26, 2010: Rashid Akmaev
! Loosely based on Hojun Wang's phi2hgt but generalized to rid of
! recursive calculations, include surface orography, and calculate
! gravity.
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Define constants
! - Earth radius (m) and
! - gravity at sea level (m/s**2)
! If used with GFS/WAM codes "use" this module
use physcons, only: re => con_rerth, g0 => con_g
implicit none
! If the module is not available, comment out the "use" line above and
! uncomment this line
! real, parameter:: re = 6.3712e+6, g0 = 9.80665e+0
real, parameter:: g0re = g0*re, g0re2 = g0*re*re
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Subroutine parameters
! INPUTS
! - array dimensions (following GFS conventios): first actual, first
! maximum, number of levels
integer, intent(in):: im,ix,levs
! - geopotential (m**2/s**2)
! - surface orography (m)
real, intent(in):: phi(ix,levs),soro(im)
! OUTPUTS
! - height (m)
real, intent(out):: z(ix,levs)
! Optional output
! - gravity (m/s**2)
! real, intent(out), optional:: grav(ix,levs)
real, intent(out):: grav(ix,levs)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Local variables
integer:: i,l
real:: phis(im)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Calculate surface geopotential
do i = 1,im
phis(i) = g0re*soro(i)/(re+soro(i))
enddo
! Calculate height
z(:,:) = 0.
do l = 1,levs
do i = 1,im
z(i,l) = re*(phis(i)+phi(i,l))/(g0re-(phis(i)+phi(i,l)))
enddo
enddo
! ***Optionally*** calculate gravity
! if(present(grav)) then
grav(:,:) = 0.
do l = 1,levs
do i = 1,im
grav(i,l) = g0re2/((re+z(i,l))*(re+z(i,l)))
enddo
enddo
! endif
end subroutine phi2z
!----------------------------------------------------------------------------
subroutine gravco2(levs,phi,soro,gg)
! Subroutine is modified from phi2z above to compute gravity for co2cin,
! the first gaussian point is chosen to represent the whole data domain
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! File history
! Dec 26, 2012: Jun Wang modified from phi2z from Rashid
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Define constants
! - Earth radius (m) and
! - gravity at sea level (m/s**2)
! If used with GFS/WAM codes "use" this module
use physcons, only: re => con_rerth, g0 => con_g
implicit none
! If the module is not available, comment out the "use" line above and
real, parameter:: g0re = g0*re, g0re2 = g0*re**2
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Subroutine parameters
! INPUTS
integer, intent(in):: levs
! - geopotential (m**2/s**2)
! - surface orography (m)
real, intent(in):: phi(levs),soro
! OUTPUTS
! Optional output
! - gravity (m/s**2)
real, intent(out):: gg(levs)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Local variables
integer:: i,l
real:: phis
real:: z(levs)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Calculate surface geopotential
phis = g0re*soro/(re+soro)
print *,'in grevco2 phis=',phis,'phi=',phi(1:100:10),'soro=',soro,
& 're=',re
! Calculate height
z(:) = 0.
do l = 1,levs
z(l) = re*(phis+phi(l))/(g0re-(phis+phi(l)))
enddo
! ***Optionally*** calculate gravity
gg(:) = 0.
do l = 1,levs
gg(l) = g0re2/((re+z(l))*(re+z(l)))
enddo
print *,'in grevco2 gg=',gg(1:100:10)
!
end subroutine gravco2
!----------------------------------------------------------------------------
subroutine getphilvl(levs,ntrac,ps,t,q,dp,gen_coord_hybrid,
& thermodyn_id,phil,prsi)
! Subroutine computes phi on a single point on model levels from p,tmp,
! and trcers for general
! hybrid for enthalpy
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! File history
! Dec 26, 2010: Jun Wang
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
use tracer_const, only : ri
implicit none
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Subroutine parameters
! INPUTS
integer, intent(in):: levs,ntrac
logical, intent(in):: gen_coord_hybrid
integer, intent(in):: thermodyn_id
!
! Local variables
real,parameter :: pa2cb=0.001,zero=0.0
! - sfc pressure (pascal)
! - pressure thickness (pascal)
! - tmp (k)
! - tracers
real, intent(in):: ps,t(levs),dp(levs),q(levs,ntrac)
! OUTPUTS
! Optional output
! - model layer enthalpy
real, intent(out):: phil(levs),prsi(levs+1)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Local variables
real:: tem,dphi,phii,sumq(levs),xr(levs)
integer :: k,n
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! init
phii = zero
! Calculate enthalpy
!
! print *,'in getphilvl,thermodyn_id=',thermodyn_id,
! & thermodyn_id.eq.3
if( gen_coord_hybrid ) then
if( thermodyn_id == 3 ) then ! Enthalpy case
!get r
sumq = zero
xr = zero
do n=1,ntrac
if( ri(n) > 0.0 ) then
do k=1,levs
xr(k) = xr(k) + q(k,n) * ri(n)
sumq(k) = sumq(k) + q(k,n)
enddo
endif
enddo
do k=1,levs
xr(k) = (1.-sumq(k))*ri(0) + xr(k)
enddo
!
!hmhj prsi(1) = ps*pa2cb
!hmhj do k=1,levs
!hmhj prsi(k+1) = prsi(k)-dp(k)*pa2cb
!hmhj enddo
!hmhj if compute prsi, we from ptop=0 with dp down to psfc
prsi(levs+1) = 0.
do k=levs,1,-1
prsi(k) = prsi(k+1) + dp(k)*pa2cb
enddo
! print *,'in getphilvl,prsi=',prsi(1:100:10)
!
do k = 1,levs
tem = xr(k) * T(k)
dphi = (prsi(k) - prsi(k+1)) * TEM
& /(prsi(k) + prsi(k+1))
phil(k) = phii + dphi
phii = phil(k) + dphi
enddo
!
else
print *,'ERROR: No phil is compute, this routine is ',
& 'for gen-hybrid with enthalpy'
!
endif
endif
!
end subroutine getphilvl
!------------------------------------------------------------------