<p><b>laura@ucar.edu</b> 2010-07-23 15:55:43 -0600 (Fri, 23 Jul 2010)</p><p>updated with calculation of initial relative humidity and water vapor; updated with parallel capabilities<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_nhyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_physics/src/core_nhyd_atmos/module_test_cases.F 2010-07-23 21:52:47 UTC (rev 427)
+++ branches/atmos_physics/src/core_nhyd_atmos/module_test_cases.F 2010-07-23 21:55:43 UTC (rev 428)
@@ -3,6 +3,7 @@
use grid_types
use configure
use constants
+ use dmpar
contains
@@ -51,7 +52,7 @@
block_ptr => domain % blocklist
do while (associated(block_ptr))
write(0,*) ' calling test case setup '
- call nhyd_test_case_squall_line(block_ptr % mesh, block_ptr % time_levs(1) % state, config_test_case)
+ call nhyd_test_case_squall_line(domain % dminfo, block_ptr % mesh, block_ptr % time_levs(1) % state, config_test_case)
write(0,*) ' returned from test case setup '
do i=2,nTimeLevs
call copy_state(block_ptr % time_levs(1) % state, block_ptr % time_levs(i) % state)
@@ -111,7 +112,7 @@
real (kind=RKIND) :: r_earth, etavs, ztemp, zd, zt, dz, gam, delt, str
- real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: rel_hum, temperature, qv
+ real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: rh, temperature, qv
real (kind=RKIND) :: ptmp, es, qvs, xnutr, znut, ptemp
integer :: iter
@@ -124,6 +125,17 @@
real (kind=RKIND) :: d1, d2, d3, cf1, cf2, cf3, cof1, cof2, psurf
+ ! storage for (lat,z) arrays for zonal velocity calculation
+
+ integer, parameter :: nlat=361
+ real (kind=RKIND), dimension(grid % nVertLevels + 1) :: zz_1d, zgrid_1d, hx_1d
+ real (kind=RKIND), dimension(grid % nVertLevels) :: flux_zonal
+ real (kind=RKIND), dimension(nlat, grid % nVertLevels) :: u_2d, etavs_2d
+ real (kind=RKIND), dimension(nlat) :: lat_2d
+ real (kind=RKIND) :: dlat
+
+ integer:: iq
+
!
! Scale all distances and areas from a unit sphere to one with radius a
!
@@ -174,7 +186,7 @@
t => state % theta % array
rt => grid % rtheta_p % array
-
+ scalars => state % scalars % array
scalars(:,:,:) = 0.
xnutr = 0.
@@ -202,7 +214,7 @@
enddo
enddo
- ! metrics for hybrid coordinate and vertical stretching
+ ! Metrics for hybrid coordinate and vertical stretching
str = 1.5
zt = 45000.
@@ -267,6 +279,10 @@
write(0,*) ' cf1, cf2, cf3 = ',cf1,cf2,cf3
+ grid % cf1 % scalar = cf1
+ grid % cf2 % scalar = cf2
+ grid % cf3 % scalar = cf3
+
do iCell=1,grid % nCells
do k=1,nz
zgrid(k,iCell) = (1.-ah(k))*(sh(k)*(zt-hx(k,iCell))+hx(k,iCell)) &
@@ -304,7 +320,124 @@
enddo
write(0,*) ' grid metrics setup complete '
+
+!************** section for 2d (lat,z) calc for zonal velocity
+
+ dlat = 0.5*pii/float(nlat-1)
+ do i = 1,nlat
+
+ lat_2d(i) = float(i-1)*dlat
+! write(0,*) ' zonal setup, latitude = ',lat_2d(i)*180./pii
+
+ do k=1,nz
+ phi = lat_2d(i)
+ hx_1d(k) = u0/gravity*cos(etavs)**1.5 &
+ *((-2.*sin(phi)**6 &
+ *(cos(phi)**2+1./3.)+10./63.) &
+ *(u0)*cos(etavs)**1.5 &
+ +(1.6*cos(phi)**3 &
+ *(sin(phi)**2+2./3.)-pii/4.)*r_earth*omega_e)
+ enddo
+
+ do k=1,nz
+ zgrid_1d(k) = (1.-ah(k))*(sh(k)*(zt-hx_1d(k))+hx_1d(k)) &
+ + ah(k) * sh(k)* zt
+ end do
+ do k=1,nz1
+ zz_1d (k) = (zw(k+1)-zw(k))/(zgrid_1d(k+1)-zgrid_1d(k))
+ end do
+
+ do k=1,nz1
+ ztemp = .5*(zgrid_1d(k+1)+zgrid_1d(k))
+ ppb(k,i) = p0*exp(-gravity*ztemp/(rgas*t0b))
+ pb (k,i) = (ppb(k,i)/p0)**(rgas/cp)
+ rb (k,i) = ppb(k,i)/(rgas*t0b*zz_1d(k))
+ tb (k,i) = t0b/pb(k,i)
+ rtb(k,i) = rb(k,i)*tb(k,i)
+ p (k,i) = pb(k,i)
+ pp (k,i) = 0.
+ rr (k,i) = 0.
+ end do
+
+
!
+! iterations to converge temperature as a function of pressure
+!
+
+ do itr = 1,10
+
+ do k=1,nz1
+ eta (k) = (ppb(k,i)+pp(k,i))/p0
+ etav(k) = (eta(k)-.252)*pii/2.
+ if(eta(k).ge.znut) then
+ teta(k) = t0*eta(k)**(rgas*dtdz/gravity)
+ else
+ teta(k) = t0*eta(k)**(rgas*dtdz/gravity) + delta_t*(znut-eta(k))**5
+ end if
+ end do
+ ! phi = grid % latCell % array (i)
+ phi = lat_2d (i)
+ do k=1,nz1
+ tt(k) = 0.
+ tt(k) = teta(k)+.75*eta(k)*pii*u0/rgas*sin(etav(k)) &
+ *sqrt(cos(etav(k)))* &
+ ((-2.*sin(phi)**6 &
+ *(cos(phi)**2+1./3.)+10./63.) &
+ *2.*u0*cos(etav(k))**1.5 &
+ +(1.6*cos(phi)**3 &
+ *(sin(phi)**2+2./3.)-pii/4.)*r_earth*omega_e)
+
+
+ ztemp = .5*(zgrid_1d(k)+zgrid_1d(k+1))
+ ptemp = ppb(k,i) + pp(k,i)
+
+ end do
+
+ do itrp = 1,25
+ do k=1,nz1
+ rr(k,i) = (pp(k,i)/(rgas*zz_1d(k)) &
+ -rb(k,i)*(tt(k)-t0b))/tt(k)
+ end do
+
+ ppi(1) = p0-.5*dzw(1)*gravity &
+ *(1.25*(rr(1,i)+rb(1,i))*(1.+qv(1,i)) &
+ -.25*(rr(2,i)+rb(2,i))*(1.+qv(2,i)))
+
+ ppi(1) = ppi(1)-ppb(1,i)
+ do k=1,nz1-1
+ ppi(k+1) = ppi(k)-.5*dzu(k+1)*gravity* &
+ (rr(k ,i)+(rr(k ,i)+rb(k ,i))*qv(k ,i) &
+ +rr(k+1,i)+(rr(k+1,i)+rb(k+1,i))*qv(k+1,i))
+ end do
+
+ do k=1,nz1
+ pp(k,i) = .2*ppi(k)+.8*pp(k,i)
+ end do
+
+ end do ! end inner iteration loop itrp
+
+ end do ! end outer iteration loop itr
+
+ do k=1,nz1
+ etavs_2d(i,k) = (0.5*(ppb(k,i)+ppb(k,i)+pp(k,i)+pp(k,i))/p0 - 0.252)*pii/2.
+! u_2d(i,k) = u0*(sin(2.*lat_2d(i))**2) *(cos(etavs_2d(i,k))**1.5)
+ u_2d(i,k) = u0*(sin(2.*lat_2d(i))**2) *(cos(etavs_2d(i,k))**1.5)*(rb(k,i)+rr(k,i))
+ end do
+
+ end do ! end loop over latitudes for 2D zonal wind field calc
+
+! do i=1,nlat
+! do k=1,nz1
+! u_2d(i,k) = u_2d(i,k) - u0*(sin(2.*lat_2d(i))**2) *(cos(etavs_2d(nlat/2,k))**1.5)
+! end do
+! end do
+!
+! write(22,*) nz1,nlat,u_2d
+
+!******************************************************************
+ write(6,*)
+ write(6,*) '--- RELATIVE HUMIDITY:', zt
+!
!---- baroclinc wave initialization ---------------------------------
!
! reference sounding based on dry isothermal atmosphere
@@ -323,14 +456,58 @@
rr (k,i) = 0.
end do
- if(i == 1) then
- do k=1,nz1
- write(0,*) ' k, ppb, pb, rb, tb (k,1) ',k,ppb(k,1),pb(k,1),rb(k,1)*zz(k,1),tb(k,1)
- enddo
- end if
+! if(i == 1) then
+! do k=1,nz1
+! write(0,*) ' k, ppb, pb, rb, tb (k,1) ',k,ppb(k,1),pb(k,1),rb(k,1)*zz(k,1),tb(k,1)
+! enddo
+! end if
+
+#ifdef DO_PHYSICS
!
+!... initialization for moisture (following initialization in hydrostatic core):
+! note that the initial profile of relative humidity is different from that
+! in the hydrostatic core since the grid is different.
+
+ 201 format(3i6,8(1x,e15.8))
+ 202 format(2i6,8(1x,e15.8))
+ if(config_microp_scheme .ne. 'off' .or. &
+ config_conv_shallow_scheme .ne. 'off' .or. &
+ config_conv_deep_scheme .ne. 'off' .or. &
+ config_pbl_scheme .ne. 'off' .or. &
+ config_eddy_scheme .ne. 'off' .or. &
+ config_radt_lw_scheme .ne. 'off' .or. &
+ config_radt_sw_scheme .ne. 'off') then
+
+!... relative humidity as function of pressure:
+! do k = 1, nz1
+! ptemp = ppb(k,i) + pp(k,i)
+! if(ptemp < 50000.) then
+! rh(k,i) = 0.0
+! elseif(p0-ptemp .lt. 0.) then
+! rh(k,i) = 1.0
+! else
+! rh(k,i) = (1.-((p0-ptemp)/50000.)**1.25)
+! endif
+! rh(k,i) = min(rh_max,rh(k,i))
+! enddo
+
+!... relative humidity as function of height:
+ do k = 1, nz1
+ ztemp = .5*(zgrid(k+1,i)+zgrid(k,i))
+ if(ztemp .gt. 5000.) then
+ rh(k,i) = 0.0
+ else
+ rh(k,i) = (1.-0.75*(ztemp/zt)**1.25)
+ endif
+ rh(k,i) = min(rh_max,rh(k,i))
+! write(6,202) i,k,ztemp,rh(k,i)
+ enddo
+ endif
+#endif
+!
! iterations to converge temperature as a function of pressure
!
+! write(6,*)
do itr = 1,10
do k=1,nz1
@@ -357,10 +534,24 @@
!write(0,*) ' k, tt(k) ',k,tt(k)
ztemp = .5*(zgrid(k,i)+zgrid(k+1,i))
ptemp = ppb(k,i) + pp(k,i)
-! qv(k,i) = env_qv( ztemp, tt(k), ptemp, 0 )
- qv(k,i) = 0.
+#ifdef DO_PHYSICS
+!
+!... computes the water vapor mixing ratio (following initialization in hydrostatic core):
+! note that the initial profile of water vapor is different from that in the hydrostatic
+! core since temperature and relative humidity are different.
+ if (tt(k) > 273.15) then
+ es = 1000.*0.6112*exp(17.67*(tt(k)-273.15)/(tt(k)-29.65))
+ else
+ es = 1000.*0.6112*exp(21.8745584*(tt(k)-273.16)/(tt(k)-7.66))
+ end if
+ qvs = (287.04/461.6)*es/(ptemp-es)
+ scalars(index_qv,k,i) = rh(k,i)*qvs
+! write(6,201) itr,i,k,ztemp,ptemp,es,qvs,rh(k,i),scalars(index_qv,k,i)
+#endif
end do
+! write(0,*)
+
! do k=2,nz1
! cqw(k,i) = 1./(1.+.5*(qv(k,i)+qv(k-1,i)))
! end do
@@ -397,11 +588,11 @@
rho (k,i) = rb(k,i) + rr(k,i)
end do
- if(i == 1) then
- do k=1,nz1
- write(0,*) ' k, p, t, rt ',k,p(k,1),t(k,1),rt(k,1)
- enddo
- end if
+! if(i == 1) then
+! do k=1,nz1
+! write(0,*) ' k, p, t, rt ',k,p(k,1),t(k,1),rt(k,1)
+! enddo
+! end if
end do ! end loop over cells
@@ -431,6 +622,7 @@
u_pert = 0.0
end if
+ call calc_flux_zonal(u_2d,etavs_2d,lat_2d,flux_zonal,lat1,lat2,grid % dvEdge % array(iEdge),a,u0,nz1,nlat)
do k=1,grid % nVertLevels
!! etavs = (0.5*(ppb(k,iCell1)+ppb(k,iCell2)+pp(k,iCell1)+pp(k,iCell2))/p0 - 0.252)*pii/2.
@@ -438,7 +630,13 @@
etavs = (0.5*(ppb(k,440)+ppb(k,440)+pp(k,440)+pp(k,440))/p0 - 0.252)*pii/2. ! 10262 mesh
! etavs = (0.5*(ppb(k,505)+ppb(k,505)+pp(k,505)+pp(k,505))/p0 - 0.252)*pii/2. ! 40962 mesh
- fluxk = u0*flux*(cos(etavs)**1.5)
+! fluxk = u0*flux*(cos(etavs)**1.5)
+
+ fluxk = u0*flux_zonal(k)/(0.5*(rb(k,iCell1)+rb(k,iCell2)+rr(k,iCell1)+rr(k,iCell2)))
+
+! if(k.eq.18) then
+! write(21,*) ' iEdge, u1, u2 ',iEdge,fluxk,u0*flux_zonal(k)
+! end if
!! fluxk = u0*flux*(cos(znuv(k))**(1.5))
!! fluxk = u0 * cos(grid % angleEdge % array(iEdge)) * (sin(lat1+lat2)**2) *(cos(etavs)**1.5)
state % u % array(k,iEdge) = fluxk + u_pert
@@ -487,17 +685,76 @@
enddo
! stop
+! write(6,*)
+! write(6,*) '--- END JW:'
+! do iCell = 1, grid%nCells
+! do k = 1,grid % nVertLevels
+! write(6,202) iCell,k,(state%scalars%array(iq,k,iCell),iq=moist_start,moist_end)
+! enddo
+! enddo
+
end subroutine nhyd_test_case_jw
+ subroutine calc_flux_zonal(u_2d,etavs_2d,lat_2d,flux_zonal,lat1_in,lat2_in,dvEdge,a,u0,nz1,nlat)
+
+ implicit none
+ integer, intent(in) :: nz1,nlat
+ real (kind=RKIND), dimension(nlat,nz1), intent(in) :: u_2d,etavs_2d
+ real (kind=RKIND), dimension(nlat), intent(in) :: lat_2d
+ real (kind=RKIND), dimension(nz1), intent(out) :: flux_zonal
+ real (kind=RKIND), intent(in) :: lat1_in, lat2_in, dvEdge, a, u0
+
+ integer :: k,i
+ real (kind=RKIND) :: lat1, lat2, w1, w2
+ real (kind=RKIND) :: dlat,da,db
+
+ lat1 = abs(lat1_in)
+ lat2 = abs(lat2_in)
+ if(lat2 <= lat1) then
+ lat1 = abs(lat2_in)
+ lat2 = abs(lat1_in)
+ end if
+
+ do k=1,nz1
+ flux_zonal(k) = 0.
+ end do
+
+ do i=1,nlat-1
+ if( (lat1 <= lat_2d(i+1)) .and. (lat2 >= lat_2d(i)) ) then
+
+ dlat = lat_2d(i+1)-lat_2d(i)
+ da = (max(lat1,lat_2d(i))-lat_2d(i))/dlat
+ db = (min(lat2,lat_2d(i+1))-lat_2d(i))/dlat
+ w1 = (db-da) -0.5*(db-da)**2
+ w2 = 0.5*(db-da)**2
+
+ do k=1,nz1
+ flux_zonal(k) = flux_zonal(k) + w1*u_2d(i,k) + w2*u_2d(i+1,k)
+ end do
+
+ end if
+
+ end do
+
+! renormalize for setting cell-face fluxes
+
+ do k=1,nz1
+ flux_zonal(k) = sign(1.,lat2_in-lat1_in)*flux_zonal(k)*dlat*a/dvEdge/u0
+ end do
+
+ end subroutine calc_flux_zonal
+
+
!----------------------------------------------------------------------------------------------------------
- subroutine nhyd_test_case_squall_line(grid, state, test_case)
+ subroutine nhyd_test_case_squall_line(dminfo, grid, state, test_case)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Setup squall line and supercell test case
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
implicit none
+ type (dm_info), intent(in) :: dminfo
type (grid_meta), intent(inout) :: grid
type (grid_state), intent(inout) :: state
integer, intent(in) :: test_case
@@ -543,7 +800,7 @@
real (kind=RKIND), dimension(grid % nVertLevels ) :: zu, dzw, rdzwp, rdzwm
real (kind=RKIND), dimension(grid % nVertLevels ) :: eta, etav, teta, ppi, tt
- real (kind=RKIND) :: d1, d2, d3, cf1, cf2, cf3
+ real (kind=RKIND) :: d1, d2, d3, cf1, cf2, cf3, cof1, cof2
real (kind=RKIND) :: ztr, thetar, ttr, thetas, um, us, zts, pitop, ptopb, rcp, rcv
real (kind=RKIND) :: radx, radz, zcent, xmid, delt, xloc, rad, temp, pres, yloc, ymid, a_scale
@@ -680,13 +937,23 @@
!********** how are we storing cf1, cf2 and cf3?
- d1 = .5*dzw(1)
- d2 = dzw(1)+.5*dzw(2)
- d3 = dzw(1)+dzw(2)+.5*dzw(3)
- cf1 = d2*d3*(d3-d2)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
- cf2 = d1*d3*(d1-d3)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
- cf3 = d1*d2*(d2-d1)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+ COF1 = (2.*DZU(2)+DZU(3))/(DZU(2)+DZU(3))*DZW(1)/DZU(2)
+ COF2 = DZU(2) /(DZU(2)+DZU(3))*DZW(1)/DZU(3)
+ CF1 = FZP(2) + COF1
+ CF2 = FZM(2) - COF1 - COF2
+ CF3 = COF2
+! d1 = .5*dzw(1)
+! d2 = dzw(1)+.5*dzw(2)
+! d3 = dzw(1)+dzw(2)+.5*dzw(3)
+! cf1 = d2*d3*(d3-d2)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+! cf2 = d1*d3*(d1-d3)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+! cf3 = d1*d2*(d2-d1)/(d2*d3*(d3-d2)+d1*d3*(d1-d3)+d1*d2*(d2-d1))
+
+ grid % cf1 % scalar = cf1
+ grid % cf2 % scalar = cf2
+ grid % cf3 % scalar = cf3
+
do iCell=1,grid % nCells
do k=1,nz
zgrid(k,iCell) = ah(k)*(zc(k)*(1.-hx(k,iCell)/zt)+hx(k,iCell)) &
@@ -784,6 +1051,9 @@
end do
end if
end do
+
+ call dmpar_bcast_reals(dminfo, nz1, grid % u_init % array)
+
!
! reference sounding based on dry atmosphere
!
@@ -796,6 +1066,8 @@
end do
pitop = pitop-.5*dzw(nz1)*gravity/(cp*tb(nz1,1)*zz(nz1,1))
+ call dmpar_bcast_real(dminfo, pitop)
+
ptopb = p0*pitop**(1./rcp)
write(6,*) 'ptopb = ',.01*ptopb
@@ -858,6 +1130,8 @@
ptop = p0*pitop**(1./rcp)
write(0,*) 'ptop = ',.01*ptop
+ call dmpar_bcast_real(dminfo, pitop)
+
do i = 1, grid % nCells
pp(nz1,i) = ptop-ptopb+.5*dzw(nz1)*gravity* &
@@ -906,6 +1180,9 @@
grid % qv_init % array(k) = scalars(index_qv,k,1)
end do
+
+ call dmpar_bcast_reals(dminfo, nz1, grid % t_init % array)
+ call dmpar_bcast_reals(dminfo, nz1, grid % qv_init % array)
!
do i=1,grid % ncells
</font>
</pre>