<p><b>duda</b> 2011-01-06 13:36:31 -0700 (Thu, 06 Jan 2011)</p><p>BRANCH COMMIT<br>
<br>
Synchronize test case initialization for hydrostatic core with<br>
atmos_physics branch.<br>
<br>
M src/core_hyd_atmos/module_test_cases.F<br>
M src/core_hyd_atmos/Registry<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_nonhydrostatic/src/core_hyd_atmos/Registry
===================================================================
--- branches/atmos_nonhydrostatic/src/core_hyd_atmos/Registry 2011-01-06 01:35:10 UTC (rev 678)
+++ branches/atmos_nonhydrostatic/src/core_hyd_atmos/Registry 2011-01-06 20:36:31 UTC (rev 679)
@@ -140,6 +140,10 @@
var persistent real uReconstructZonal ( nVertLevels nCells Time ) 1 o uReconstructZonal diag - -
var persistent real uReconstructMeridional ( nVertLevels nCells Time ) 1 o uReconstructMeridional diag - -
+# For moist initialization
+var persistent real qsat ( nVertLevels nCells Time ) 1 o qsat diag_physics - -
+var persistent real relhum ( nVertLevels nCells Time ) 1 o relhum diag_physics - -
+
# Tendency variables
var persistent real tend_h ( nVertLevels nCells Time ) 1 - h tend - -
var persistent real tend_u ( nVertLevels nEdges Time ) 1 - u tend - -
Modified: branches/atmos_nonhydrostatic/src/core_hyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_nonhydrostatic/src/core_hyd_atmos/module_test_cases.F 2011-01-06 01:35:10 UTC (rev 678)
+++ branches/atmos_nonhydrostatic/src/core_hyd_atmos/module_test_cases.F 2011-01-06 20:36:31 UTC (rev 679)
@@ -35,7 +35,8 @@
if (config_test_case == 3) write(0,*) ' normal-mode perturbation included '
block_ptr => domain % blocklist
do while (associated(block_ptr))
- call hyd_test_case_1(block_ptr % mesh, block_ptr % state % time_levs(1) % state, config_test_case)
+ call hyd_test_case_1(block_ptr % mesh, block_ptr % state % time_levs(1) % state, &
+ block_ptr % diag_physics, config_test_case)
do i=2,nTimeLevs
call copy_state(block_ptr % state % time_levs(i) % state, block_ptr % state % time_levs(1) % state)
end do
@@ -52,7 +53,7 @@
!----------------------------------------------------------------------------------------------------------
- subroutine hyd_test_case_1(grid, state, test_case)
+ subroutine hyd_test_case_1(grid, state, diag_physics, test_case)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Setup baroclinic wave test case from Jablonowski and Williamson 2008 (QJRMS)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -61,6 +62,7 @@
type (mesh_type), intent(inout) :: grid
type (state_type), intent(inout) :: state
+ type (diag_physics_type), intent(inout) :: diag_physics
integer, intent(in) :: test_case
real (kind=RKIND), parameter :: u0 = 35.0
@@ -70,24 +72,31 @@
real (kind=RKIND), parameter :: u_perturbation = 1., pert_radius = 0.1, latitude_pert = 40., longitude_pert = 20.
real (kind=RKIND), parameter :: theta_c = pii/4.0
real (kind=RKIND), parameter :: lambda_c = 3.0*pii/2.0
- real (kind=RKIND), parameter :: rh_max = 0.4 ! Maximum relative humidity
real (kind=RKIND), parameter :: k_x = 9. ! Normal mode wave number
real (kind=RKIND), dimension(:), pointer :: rdnu, rdnw, fnm, fnp, dbn, dnu, dnw
real (kind=RKIND), dimension(:), pointer :: surface_pressure
real (kind=RKIND), dimension(:,:), pointer :: pressure, theta, alpha, geopotential, h
- real (kind=RKIND), dimension(:,:,:), pointer :: scalars
- integer :: iCell, iEdge, vtx1, vtx2, ivtx, k, nz, nz1, index_qv
+ integer :: iCell, iEdge, vtx1, vtx2, ivtx, k, nz, nz1
real (kind=RKIND) :: u, v, flux, fluxk, lat1, lat2, eta_v, r_pert, u_pert, lat_pert, lon_pert, r
real (kind=RKIND) :: ptop, p0, phi
real (kind=RKIND) :: lon_Edge
- real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: rel_hum, temperature
+ real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: temperature
real (kind=RKIND) :: ptmp, es, qvs
integer :: iter
+!.. initialization of moisture:
+ integer:: index_qv
+! real (kind=RKIND),parameter :: rh_max = 0.40 ! Maximum relative humidity
+ real (kind=RKIND),parameter :: rh_max = 0.70 ! Maximum relative humidity
+ real (kind=RKIND):: ttemp
+ real (kind=RKIND),dimension(:,:), pointer:: qsat,relhum
+ real (kind=RKIND),dimension(:,:,:),pointer:: scalars
+!.. end initialization of moisture.
+
! real (kind=RKIND), dimension(grid % nVertLevelsP1 ) :: hyai, hybi, znu, znw, znwc, znwv, hyam, hybm
! real (kind=RKIND), dimension(grid % nVertLevelsP1 ) :: znuc, znuv, bn, divh, dpn, teta, phi
real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: hyai, hybi, znu, znw, znwc, znwv, hyam, hybm
@@ -149,10 +158,16 @@
alpha => state % alpha % array
geopotential => state % geopotential % array
h => state % h % array
+
+!.. initialization of moisture:
scalars => state % scalars % array
+ qsat => diag_physics % qsat % array
+ relhum => diag_physics % relhum % array
+ scalars(:,:,:) = 0.0
+ qsat(:,:) = 0.0
+ relhum(:,:) = 0.0
+!.. end initialization of moisture.
- scalars(:,:,:) = 0.
-
p0 = 100000.
bn (1) = 1.
znw(1) = 1.
@@ -366,27 +381,24 @@
theta(k,1)*(pressure(k,1)/p0)**(rgas/cp)
end do
-!
-! initialization for moisture
-!
- if (config_mp_physics /= 0) then
+!.. calculation of relative humidity:
+! if(config_mp_physics /= 0) then
+! do iCell = 1, grid % nCells
+! do k = 1, nz1
+! ptmp = 0.5*(pressure(k,iCell) + pressure(k+1,iCell))
+! if(ptmp < 50000.) then
+! relhum(k,iCell) = 0.0
+! else
+! relhum(k,iCell) = (1.-((p0-ptmp)/50000.)**1.25)
+! end if
+! relhum(k,iCell) = min(rh_max,relhum(k,iCell))
+! enddo
+! enddo
+! else
+! relhum(:,:) = 0.
+! endif
+!.. end calculation of relative humidity.
- do iCell=1,grid % nCells
- do k=1,nz1
- ptmp = 0.5*(pressure(k,iCell) + pressure(k+1,iCell))
- if (ptmp < 50000.) then
- rel_hum(k,iCell) = 0.0
- else
- rel_hum(k,iCell) = (1.-((p0-ptmp)/50000.)**1.25)
- end if
- rel_hum(k,iCell) = min(rh_max,rel_hum(k,iCell))
- end do
- end do
- else
- rel_hum(:,:) = 0.
- end if
-
-
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! iteration
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -411,25 +423,38 @@
*2.*u0*cos(znuv(k))**1.5 &
+(1.6*cos(phi)**3 &
*(sin(phi)**2+2./3.)-pii/4.)*a*omega)
+
+ if(config_mp_physics /= 0) then
+ !.. ldf (11-15-2010): move calculation of the relative humidity inside iteration loop to
+ !.. calculate the saturation mixing ratio and initialize the water vapor mixing ratio.
+ !.. We assume the temperature computed above is the virtual temperature when moisture is
+ !.. present:
+ if(ptmp < 50000.) then
+ relhum(k,iCell) = 0.0
+ else
+ relhum(k,iCell) = (1.-((p0-ptmp)/50000.)**1.25)
+ endif
+ relhum(k,iCell) = min(rh_max,relhum(k,iCell))
+
+ !.. conversion from virtual temperature to temperature:
+ ttemp = temperature(k,iCell)/(1+0.608*scalars(index_qv,k,iCell))
+
+ !.. calculation of water vapor mixing ratio:
+ if (ttemp > 273.15) then
+ es = 1000.*0.6112*exp(17.67*(ttemp-273.15)/(ttemp-29.65))
+ else
+ es = 1000.*0.6112*exp(21.8745584*(ttemp-273.15)/(ttemp-7.66))
+ end if
+ qsat(k,iCell) = (287.04/461.6)*es/(ptmp-es)
+ if(relhum(k,iCell) .eq. 0.0) qsat(k,iCell) = 0.0
+ scalars(index_qv,k,iCell) = relhum(k,iCell)*qsat(k,iCell)
+ endif
+
+ !conversion of virtual temperature to potential temperature:
+ theta (k,iCell) = temperature(k,iCell)/(1+0.608*scalars(index_qv,k,iCell))*(ptmp/p0)**(-rgas/cp)
+ alpha (k,iCell) = (rgas/p0)*theta(k,iCell)*(1.+1.61*scalars(index_qv,k,iCell))*(ptmp/p0)**cvpm
+ enddo
- temperature(k,iCell) = temperature(k,iCell)/(1.+0.61*scalars(index_qv,k,iCell))
-
- theta (k,iCell) = temperature(k,iCell)* &
- (0.5*(pressure(k,iCell)+pressure(k+1,iCell))/p0)**(-rgas/cp)
- alpha (k,iCell) = (rgas/p0)*theta(k,iCell)*(1.+1.61*scalars(index_qv,k,iCell))* &
- (0.5*(pressure(k,iCell)+pressure(k+1,iCell))/p0)**cvpm
-
- if (temperature(k,iCell) > 273.15) then
- es = 1000.*0.6112*exp(17.67*(temperature(k,iCell)-273.15)/(temperature(k,iCell)-29.65))
- else
- es = 1000.*0.6112*exp(21.8745584*(temperature(k,iCell)-273.16)/(temperature(k,iCell)-7.66))
- end if
- qvs = (287.04/461.6)*es/(ptmp-es)
-! qvs = 380.*exp(17.27*(temperature(k,iCell)-273.)/(temperature(k,iCell)-36.))/ptmp
-
- scalars(index_qv,k,iCell) = rel_hum(k,iCell)*qvs
- end do
-
do k=nz1,1,-1
pressure(k,iCell) = pressure(k+1,iCell)-dnw(k)*h(k,iCell)*(1.+scalars(index_qv,k,iCell))
geopotential(k,iCell) = geopotential(k+1,iCell)+dnw(k)*h(k,iCell)*alpha(k,iCell)
@@ -438,13 +463,23 @@
end do
end do
- write(6,*) 'ptop = ',ptop,' zt = ',geopotential(nz,1)/gravity
+ 202 format(2i6,10(1x,e15.8))
+ write(0,*) '--- initialization of water vapor:'
+ do iCell = 1, grid % nCells
+ if(iCell == 1 .or. iCell == grid % nCells) then
+ write(0,*)
+ do k = nz1, 1, -1
+ write(0,202) iCell,k,temperature(k,iCell),relhum(k,iCell),qsat(k,iCell),scalars(index_qv,k,iCell)
+ enddo
+ endif
+ enddo
- write(6,*) ' full sounding with moisture'
- do k=1,nz1
- write(6,*) k, geopotential(k,1)/gravity, 0.01*pressure(k,1), theta(k,1), &
- theta(k,1)*(pressure(k,1)/p0)**(rgas/cp)
- end do
+! write(6,*) 'ptop = ',ptop,' zt = ',geopotential(nz,1)/gravity
+! write(6,*) ' full sounding with moisture'
+! do k=1,nz1
+! write(6,*) k, geopotential(k,1)/gravity, 0.01*pressure(k,1), theta(k,1), &
+! theta(k,1)*(pressure(k,1)/p0)**(rgas/cp)
+! end do
! When initializing a scalar, be sure not to put unreasonably large values
! into indices in the moist class
</font>
</pre>