<p><b>duda</b> 2011-01-25 17:17:05 -0700 (Tue, 25 Jan 2011)</p><p>BRANCH COMMIT<br>
<br>
Add code to extend the set of fields that are interpolated for<br>
nonhydrostatic real-data initialization to include all fields needed for<br>
a dynamics-only run, plus a few land-surface fields. <br>
<br>
Several parts of the initialization code are still hard-wired to use GFS<br>
data processed by the ungrib.exe program from the WPS.<br>
<br>
<br>
M src/core_init_nhyd_atmos/module_test_cases.F<br>
M src/core_init_nhyd_atmos/Registry<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_init_nhyd_atmos/Registry
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/Registry 2011-01-25 00:54:44 UTC (rev 704)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/Registry 2011-01-26 00:17:05 UTC (rev 705)
@@ -122,8 +122,12 @@
var persistent real cf2 ( ) 0 io cf2 mesh - -
var persistent real cf3 ( ) 0 io cf3 mesh - -
-# land use
+# static terrestrial fields
+var persistent real ter ( nCells ) 0 io ter mesh - -
+var persistent integer landmask ( nCells ) 0 io landmask mesh - -
var persistent integer lu_index ( nCells ) 0 io lu_index mesh - -
+var persistent integer soilcat_top ( nCells ) 0 io soilcat_top mesh - -
+var persistent integer soilcat_bot ( nCells ) 0 io soilcat_bot mesh - -
# description of the vertical grid structure
@@ -156,6 +160,22 @@
var persistent real psfc_fg ( nCells ) 0 o psfc fg - -
var persistent real pmsl_fg ( nCells ) 0 o pmsl fg - -
+# Horizontally interpolated from first-guess data
+# and should be read in by model
+var persistent real sm000010 ( nCells ) 0 o sm000010 fg - -
+var persistent real sm010040 ( nCells ) 0 o sm010040 fg - -
+var persistent real sm040100 ( nCells ) 0 o sm040100 fg - -
+var persistent real sm100200 ( nCells ) 0 o sm100200 fg - -
+var persistent real st000010 ( nCells ) 0 o st000010 fg - -
+var persistent real st010040 ( nCells ) 0 o st010040 fg - -
+var persistent real st040100 ( nCells ) 0 o st040100 fg - -
+var persistent real st100200 ( nCells ) 0 o st100200 fg - -
+var persistent real skintemp ( nCells ) 0 o skintemp fg - -
+var persistent real snow ( nCells ) 0 o snow fg - -
+var persistent real seaice ( nCells ) 0 o seaice fg - -
+var persistent real gfs_z ( nVertLevels nCells ) 0 o gfs_z fg - -
+var persistent real gfs_p ( nVertLevelsP1 nCells ) 0 o gfs_p fg - -
+
# Prognostic variables: read from input, saved in restart, and written to output
var persistent real u ( nVertLevels nEdges Time ) 2 io u state - -
var persistent real w ( nVertLevelsP1 nCells Time ) 2 io w state - -
Modified: branches/atmos_physics/src/core_init_nhyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/module_test_cases.F 2011-01-25 00:54:44 UTC (rev 704)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/module_test_cases.F 2011-01-26 00:17:05 UTC (rev 705)
@@ -5,6 +5,7 @@
use constants
use dmpar
use advection
+ use sort
contains
@@ -1992,7 +1993,7 @@
type (proj_info) :: proj
real (kind=RKIND), dimension(:), pointer :: rdzw, dzu, rdzu, fzm, fzp
- real (kind=RKIND), dimension(:), pointer :: vert_level, latPoints, lonPoints
+ real (kind=RKIND), dimension(:), pointer :: vert_level, latPoints, lonPoints, ter
real (kind=RKIND), dimension(:,:), pointer :: zgrid, zx, zz, hx
real (kind=RKIND), dimension(:,:), pointer :: pressure, ppb, pb, rho, rb, rr, tb, rtb, p, pp, dss, t, rt
real (kind=RKIND), dimension(:), pointer :: destField1d
@@ -2001,17 +2002,26 @@
real (kind=RKIND), dimension(:,:,:), pointer :: scalars
real (kind=RKIND), dimension(:,:,:), pointer :: deriv_two
+ real (kind=RKIND) :: target_z
integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, nz, nz1, itr, itrp, cell1, cell2, nCellsSolve
integer :: nInterpPoints, ndims
!This is temporary variable here. It just need when calculate tangential velocity v.
integer :: eoe, j
- integer, dimension(:), pointer :: nEdgesOnEdge
- integer, dimension(:,:), pointer :: edgesOnEdge, CellsOnEdge, cellsOnCell
- real (kind=RKIND), dimension(:), pointer :: dvEdge, AreaCell
+ integer, dimension(:), pointer :: nEdgesOnEdge, nEdgesOnCell
+ integer, dimension(:,:), pointer :: edgesOnEdge, cellsOnEdge, edgesOnCell, cellsOnCell
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, AreaCell
real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge
real (kind=RKIND), dimension(:,:), pointer :: v
+ real (kind=RKIND), dimension(:,:), pointer :: sorted_arr
+ real(kind=RKIND), dimension(:), pointer :: hs, hs1
+ real(kind=RKIND) :: hm, zh, dzmin, dzmina, dzminf, sm
+ integer :: nsmterrain, kz, sfc_k
+ logical :: hybrid, smooth
+
+ logical :: config_do_static
+
! For interpolating terrain and land use
integer :: nx, ny, nzz, iPoint, subx, suby
integer :: isigned, endian, wordsize, istatus
@@ -2033,7 +2043,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) :: ptmp, es, qvs, xnutr, znut, ptemp
+ real (kind=RKIND) :: ptmp, es, rs, rgas_moist, qvs, xnutr, znut, ptemp, rcv
integer :: iter
real (kind=RKIND), dimension(grid % nVertLevels + 1) :: hyai, hybi, znu, znw, znwc, znwv, hyam, hybm
@@ -2043,8 +2053,6 @@
real (kind=RKIND), dimension(grid % nVertLevels) :: zu, dzw, rdzwp, rdzwm
real (kind=RKIND), dimension(grid % nVertLevels) :: eta, etav, teta, ppi, tt
- real (kind=RKIND), dimension(grid % nCells) :: hs
-
real (kind=RKIND) :: d1, d2, d3, cf1, cf2, cf3, cof1, cof2, psurf
! storage for (lat,z) arrays for zonal velocity calculation
@@ -2053,14 +2061,18 @@
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(grid % nVertLevels + 1) :: fsum
real (kind=RKIND), dimension(nlat) :: lat_2d
real (kind=RKIND) :: dlat
real (kind=RKIND) :: z_edge, z_edge3, d2fdx2_cell1, d2fdx2_cell2
weightsOnEdge => grid % weightsOnEdge % array
nEdgesOnEdge => grid % nEdgesOnEdge % array
+ nEdgesOnCell => grid % nEdgesOnCell % array
edgesOnEdge => grid % edgesOnEdge % array
+ edgesOnCell => grid % edgesOnCell % array
dvEdge => grid % dvEdge % array
+ dcEdge => grid % dcEdge % array
AreaCell => grid % AreaCell % array
CellsOnEdge => grid % CellsOnEdge % array
cellsOnCell => grid % cellsOnCell % array
@@ -2080,6 +2092,7 @@
zx => grid % zx % array
zz => grid % zz % array
hx => grid % hx % array
+ ter => grid % ter % array
dss => grid % dss % array
pb => diag % exner_base % array
@@ -2102,10 +2115,24 @@
nz = nz1 + 1
nCellsSolve = grid % nCellsSolve
+ xnutr = 0.
+ zd = 12000.
+ znut = eta_t
+ etavs = (1.-0.252)*pii/2.
+ rcv = rgas/(cp-rgas)
+ r_earth = a
+ p0 = 1.e+05
+
+
+ config_do_static = .true.
+
!
! Scale all distances and areas from a unit sphere to one with radius a
!
+
+ if (config_do_static) then
+
grid % xCell % array = grid % xCell % array * a
grid % yCell % array = grid % yCell % array * a
grid % zCell % array = grid % zCell % array * a
@@ -2126,16 +2153,7 @@
call initialize_advection_rk(grid)
call initialize_deformation_weights(grid)
- xnutr = 0.
- zd = 12000.
- znut = eta_t
- etavs = (1.-0.252)*pii/2.
- r_earth = a
- p0 = 1.e+05
-
-
-#if 0
!
! Interpolate HGT
!
@@ -2151,7 +2169,7 @@
allocate(rarray(nx,ny,nzz))
allocate(nhs(grid % nCells))
nhs(:) = 0
- grid % hx % array(:,:) = 0.0
+ grid % ter % array(:) = 0.0
! do jTileStart=1,20401,ny-6
do jTileStart=1,961,ny-6
@@ -2184,7 +2202,7 @@
grid % nCells, grid % maxEdges, grid % nEdgesOnCell % array, grid % cellsOnCell % array, &
grid % latCell % array, grid % lonCell % array)
- grid % hx % array(1,iPoint) = grid % hx % array(1,iPoint) + rarray(i,j,1)
+ grid % ter % array(iPoint) = grid % ter % array(iPoint) + rarray(i,j,1)
nhs(iPoint) = nhs(iPoint) + 1
end do
@@ -2194,42 +2212,13 @@
end do
do iCell=1, grid % nCells
- grid % hx % array(1,iCell) = grid % hx % array(1,iCell) / real(nhs(iCell))
+ grid % ter % array(iCell) = grid % ter % array(iCell) / real(nhs(iCell))
end do
deallocate(rarray)
deallocate(nhs)
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! SMOOTH TOPOGRAPHY TO CREATE 3D HX FIELD
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- do k=2,grid%nVertLevels
- grid % hx % array(k,:) = grid % hx % array(k-1,:)
- do i=1,4
- do iCell=1,grid%nCells
- hs(iCell) = grid % hx % array(k,iCell)
- do j = 1,grid % nEdgesOnCell % array(iCell)
- hs(iCell) = hs(iCell) + grid % hx % array(k,cellsOnCell(j,iCell))
- end do
- hs(iCell) = hs(iCell) / real(grid%nEdgesOnCell%array(iCell)+1)
- end do
- do iCell=1,grid%nCells
- grid % hx % array(k,iCell) = hs(iCell)
- do j = 1,grid % nEdgesOnCell % array(iCell)
- grid % hx % array(k,iCell) = grid % hx % array(k,iCell) + hs(cellsOnCell(j,iCell))
- end do
- grid % hx % array(k,iCell) = grid % hx % array(k,iCell) / real(grid%nEdgesOnCell%array(iCell)+1)
- end do
- end do
- end do
-
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! CREATE VERTICAL GRID
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-
!
! Interpolate LU_INDEX
!
@@ -2291,8 +2280,438 @@
deallocate(rarray)
deallocate(ncat)
+ !
+ ! Derive LANDMASK
+ !
+ grid % landmask % array(:) = 0
+ do iCell=1, grid % nCells
+ if (grid % lu_index % array(iCell) /= 16) grid % landmask % array(iCell) = 1
+ end do
+
!
+ ! Interpolate SOILCAT_TOP
+ !
+ nx = 1200
+ ny = 1200
+ nzz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 1
+ scalefactor = 1.0
+ allocate(rarray(nx,ny,nzz))
+ allocate(ncat(16,grid % nCells))
+ ncat(:,:) = 0
+ grid % soilcat_top % array(:) = 0.0
+
+ do jTileStart=1,20401,ny
+ jTileEnd = jTileStart + ny - 1
+ do iTileStart=1,42001,nx
+ iTileEnd = iTileStart + nx - 1
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'/soiltype_top_30s/',iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
+write(0,*) trim(fname)
+
+ call read_geogrid(fname, len_trim(fname), &
+ rarray, &
+ nx, ny, nzz, &
+ isigned, endian, scalefactor, wordsize, istatus)
+write(0,*) istatus
+
+ iPoint = 1
+ do j=1,ny
+ do i=1,nx
+ lat_pt = -89.99583 + (jTileStart + j - 2) * 0.0083333333
+ lon_pt = -179.99583 + (iTileStart + i - 2) * 0.0083333333
+ lat_pt = lat_pt * pii / 180.0
+ lon_pt = lon_pt * pii / 180.0
+
+ iPoint = nearest_cell(lat_pt, lon_pt, &
+ iPoint, &
+ grid % nCells, grid % maxEdges, grid % nEdgesOnCell % array, grid % cellsOnCell % array, &
+ grid % latCell % array, grid % lonCell % array)
+
+ ncat(int(rarray(i,j,1)),iPoint) = ncat(int(rarray(i,j,1)),iPoint) + 1
+
+ end do
+ end do
+
+ end do
+ end do
+
+ do iCell=1, grid % nCells
+ grid % soilcat_top % array(iCell) = 1
+ do i=2,16
+ if (ncat(i,iCell) > ncat(grid % soilcat_top % array(iCell),iCell)) then
+ grid % soilcat_top % array(iCell) = i
+ end if
+ end do
+ end do
+
+ deallocate(rarray)
+ deallocate(ncat)
+
+
+ !
+ ! Interpolate SOILCAT_BOT
+ !
+ nx = 1200
+ ny = 1200
+ nzz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 1
+ scalefactor = 1.0
+ allocate(rarray(nx,ny,nzz))
+ allocate(ncat(16,grid % nCells))
+ ncat(:,:) = 0
+ grid % soilcat_bot % array(:) = 0.0
+
+ do jTileStart=1,20401,ny
+ jTileEnd = jTileStart + ny - 1
+ do iTileStart=1,42001,nx
+ iTileEnd = iTileStart + nx - 1
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'/soiltype_bot_30s/',iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
+write(0,*) trim(fname)
+
+ call read_geogrid(fname, len_trim(fname), &
+ rarray, &
+ nx, ny, nzz, &
+ isigned, endian, scalefactor, wordsize, istatus)
+write(0,*) istatus
+
+ iPoint = 1
+ do j=1,ny
+ do i=1,nx
+ lat_pt = -89.99583 + (jTileStart + j - 2) * 0.0083333333
+ lon_pt = -179.99583 + (iTileStart + i - 2) * 0.0083333333
+ lat_pt = lat_pt * pii / 180.0
+ lon_pt = lon_pt * pii / 180.0
+
+ iPoint = nearest_cell(lat_pt, lon_pt, &
+ iPoint, &
+ grid % nCells, grid % maxEdges, grid % nEdgesOnCell % array, grid % cellsOnCell % array, &
+ grid % latCell % array, grid % lonCell % array)
+
+ ncat(int(rarray(i,j,1)),iPoint) = ncat(int(rarray(i,j,1)),iPoint) + 1
+
+ end do
+ end do
+
+ end do
+ end do
+
+ do iCell=1, grid % nCells
+ grid % soilcat_bot % array(iCell) = 1
+ do i=2,16
+ if (ncat(i,iCell) > ncat(grid % soilcat_bot % array(iCell),iCell)) then
+ grid % soilcat_bot % array(iCell) = i
+ end if
+ end do
+ end do
+
+ deallocate(rarray)
+ deallocate(ncat)
+
+ end if ! config_do_static
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! BEGIN ADOPT GFS TERRAIN HEIGHT
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+ call read_met_init(trim(config_met_prefix), .false., trim(config_init_date), istatus)
+
+ if (istatus /= 0) then
+ write(0,*) 'Error reading initial met data'
+ return
+ end if
+
+ call read_next_met_field(field, istatus)
+ do while (istatus == 0)
+ if (index(field % field, 'SOILHGT') /= 0) then
+
+ !
+ ! Set up projection
+ !
+ call map_init(proj)
+
+ if (field % iproj == PROJ_LATLON) then
+ call map_set(PROJ_LATLON, proj, &
+ latinc = field % deltalat, &
+ loninc = field % deltalon, &
+ knowni = 1.0_4, &
+ knownj = 1.0_4, &
+ lat1 = field % startlat, &
+ lon1 = field % startlon)
+ end if
+
+
+ if (index(field % field, 'SOILHGT') /= 0) then
+ nInterpPoints = grid % nCells
+ latPoints => grid % latCell % array
+ lonPoints => grid % lonCell % array
+ destField1d => grid % ter % array
+ ndims = 1
+ end if
+
+ do i=1,nInterpPoints
+ lat = latPoints(i)*DEG_PER_RAD
+ lon = lonPoints(i)*DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ if (x < 0.0) then
+ lon = lon + 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ end if
+ if (ndims == 1) then
+ destField1d(i) = four_pt(field % nx, field % ny, field % slab, x, y)
+ else if (ndims == 2) then
+ destField2d(k,i) = four_pt(field % nx, field % ny, field % slab, x, y)
+ end if
+ end do
+ end if
+
+ deallocate(field % slab)
+ call read_next_met_field(field, istatus)
+ end do
+
+ call read_met_close()
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! END ADOPT GFS TERRAIN HEIGHT
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+
+ !
+ ! Vertical grid setup
+ !
+ allocate(hs (grid%nCells))
+ allocate(hs1(grid%nCells))
+
+! Fourth order smoother for terrain
+
+! nsmterrain = 2
+ nsmterrain = 0 ! NO SMOOTHING WHEN USING GFS TERRAIN HEIGHT DIRECTLY
+
+ do i=1,nsmterrain
+
+ do iCell=1,grid%nCells
+ hs(iCell) = 0.
+ do j = 1,nEdgesOnCell(iCell)
+ hs(iCell) = hs(iCell) + dvEdge(edgesOnCell(j,iCell)) &
+ / dcEdge(edgesOnCell(j,iCell)) &
+ * (ter(cellsOnCell(j,iCell))-ter(iCell))
+ end do
+ hs(iCell) = ter(iCell) + 0.125*hs(iCell)
+ end do
+
+ do iCell=1,grid %nCells
+ ter(iCell) = 0.
+ do j = 1,nEdgesOnCell(iCell)
+ ter(iCell) = ter(iCell) + dvEdge(edgesOnCell(j,iCell)) &
+ / dcEdge(edgesOnCell(j,iCell)) &
+ * (hs(cellsOnCell(j,iCell))-hs(iCell))
+ end do
+ ter(iCell) = hs(iCell) - 0.125*ter(iCell)
+ end do
+
+ end do
+
+ do iCell=1,grid % nCells
+ hx(:,iCell) = ter(iCell)
+ end do
+
+ hm = maxval(ter(:))
+ write(6,*) "max ter = ", hm
+
+! Metrics for hybrid coordinate and vertical stretching
+
+ str = 1.5
+! str = 1.
+ zt = 24000.
+ dz = zt/float(nz1)
+
+ do k=1,nz
+ zw(k) = (real(k-1)/real(nz1))**str*zt
+ if (k > 1) dzw(k-1) = zw(k)-zw(k-1)
+ end do
+
+! ah(k) governs the transition between terrain-following
+! and pure height coordinates
+! ah(k) = 1 is a smoothed terrain-following coordinate
+! ah(k) = 1.-zw(k)/zt is the basic terrain-following coordinate
+! ah(k) = 0 is a height coordinate
+
+! hybrid = .true.
+ hybrid = .false.
+
+ kz = nz
+ if (hybrid) then
+
+ zh = zt
+
+ do k=1,nz
+ if (zw(k) < zh) then
+ ah(k) = cos(.5*pii*zw(k)/zh)**6
+
+!!! ah(k) = ah(k)*(1.-zw(k)/zt)
+
+ else
+ ah(k) = 0.
+ kz = min(kz,k)
+ end if
+ end do
+
+ else
+
+ do k=1,nz
+ ah(k) = 1.-zw(k)/zt
+ end do
+
+ end if
+
+ do k=1,nz
+ write(6,*) k,zw(k), ah(k)
+ end do
+
+ do k=1,nz1
+ dzw (k) = zw(k+1)-zw(k)
+ rdzw(k) = 1./dzw(k)
+ zu(k ) = .5*(zw(k)+zw(k+1))
+ end do
+ do k=2,nz1
+ dzu (k) = .5*(dzw(k)+dzw(k-1))
+ rdzu(k) = 1./dzu(k)
+ fzp (k) = .5* dzw(k )/dzu(k)
+ fzm (k) = .5* dzw(k-1)/dzu(k)
+ rdzwp(k) = dzw(k-1)/(dzw(k )*(dzw(k)+dzw(k-1)))
+ rdzwm(k) = dzw(k )/(dzw(k-1)*(dzw(k)+dzw(k-1)))
+ end do
+
+!********** how are we storing cf1, cf2 and cf3?
+
+ 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))
+
+ write(0,*) ' cf1, cf2, cf3 = ',cf1,cf2,cf3
+
+ grid % cf1 % scalar = cf1
+ grid % cf2 % scalar = cf2
+ grid % cf3 % scalar = cf3
+
+! Smoothing algorithm for coordinate surfaces
+
+ smooth = .true.
+! smooth = .false.
+
+ if (smooth) then
+
+ dzmin = 0.2
+! dzmin = 0.4
+
+ do k=2,kz-1
+ hx(k,:) = hx(k-1,:)
+ dzminf = zw(k)-zw(k-1)
+
+! dzmin = max(.5,1.-.5*zw(k)/hm)
+
+ sm = .0625*min(.5*zw(k)/hm,1.)
+
+ do i=1,50
+ do iCell=1,grid %nCells
+ hs1(iCell) = 0.
+ do j = 1,nEdgesOnCell(iCell)
+
+ hs1(iCell) = hs1(iCell) + dvEdge(edgesOnCell(j,iCell)) &
+ / dcEdge(edgesOnCell(j,iCell)) &
+ * (hx(k,cellsOnCell(j,iCell))-hx(k,iCell))
+ end do
+ hs1(iCell) = hx(k,iCell) + sm*hs1(iCell)
+
+ hs(iCell) = 0.
+ do j = 1,nEdgesOnCell(iCell)
+ hs(iCell) = hs(iCell) + dvEdge(edgesOnCell(j,iCell)) &
+ / dcEdge(edgesOnCell(j,iCell)) &
+ * (hs1(cellsOnCell(j,iCell))-hs1(iCell))
+ end do
+ hs(iCell) = hs1(iCell) - 0.*hs(iCell)
+
+ end do
+ dzmina = minval(hs(:)-hx(k-1,:))
+ if (dzmina >= dzmin*(zw(k)-zw(k-1))) then
+ hx(k,:)=hs(:)
+ dzminf = dzmina
+ else
+ exit
+ end if
+ end do
+ write(6,*) k,i,sm,dzminf/(zw(k)-zw(k-1)),dzmina/(zw(k)-zw(k-1))
+ end do
+
+ do k=kz,nz
+ hx(k,:) = 0.
+ end do
+ else
+
+ do k=2,nz1
+ dzmina = minval(zw(k)+ah(k)*hx(k,:)-zw(k-1)-ah(k-1)*hx(k-1,:))
+ write(6,*) k,dzmina/(zw(k)-zw(k-1))
+ end do
+
+ end if
+
+ deallocate(hs )
+ deallocate(hs1)
+
+! Height of coordinate levels (calculation of zgrid)
+
+ do iCell=1,grid % nCells
+ do k=1,nz
+ zgrid(k,iCell) = zw(k) + ah(k)*hx(k,iCell)
+ end do
+ do k=1,nz1
+ zz (k,iCell) = (zw(k+1)-zw(k))/(zgrid(k+1,iCell)-zgrid(k,iCell))
+ end do
+ end do
+
+ do i=1, grid % nEdges
+ iCell1 = grid % CellsOnEdge % array(1,i)
+ iCell2 = grid % CellsOnEdge % array(2,i)
+ do k=1,nz
+ zx (k,i) = (zgrid(k,iCell2)-zgrid(k,iCell1)) / grid % dcEdge % array(i)
+ end do
+ end do
+ do i=1, grid % nCells
+ do k=1,nz1
+ ztemp = .5*(zgrid(k+1,i)+zgrid(k,i))
+ dss(k,i) = 0.
+ ztemp = zgrid(k,i)
+ if (ztemp.gt.zd+.1) then
+ dss(k,i) = dss(k,i)+xnutr*sin(.5*pii*(ztemp-zd)/(zt-zd))**2
+ end if
+ end do
+ enddo
+
+! do k=1,nz1
+! write(0,*) ' k, zgrid(k,1),hx(k,1) ',k,zgrid(k,1),hx(k,1)
+! enddo
+
+! do k=1,nz1
+! write(0,*) ' k, zx(k,1) ',k,zx(k,1)
+! enddo
+
+ write(0,*) ' grid metrics setup complete '
+
+
+ !
! Horizontally interpolate meteorological data
!
allocate(vert_level(config_nfglevels))
@@ -2315,11 +2734,34 @@
index(field % field, 'PMSL') /= 0 .or. &
index(field % field, 'PSFC') /= 0 .or. &
index(field % field, 'SOILHGT') /= 0 .or. &
- index(field % field, 'PRES') /= 0) then
+ index(field % field, 'SM000010') /= 0 .or. &
+ index(field % field, 'SM010040') /= 0 .or. &
+ index(field % field, 'SM040100') /= 0 .or. &
+ index(field % field, 'SM100200') /= 0 .or. &
+ index(field % field, 'ST000010') /= 0 .or. &
+ index(field % field, 'ST010040') /= 0 .or. &
+ index(field % field, 'ST040100') /= 0 .or. &
+ index(field % field, 'ST100200') /= 0 .or. &
+ index(field % field, 'PRES') /= 0 .or. &
+ index(field % field, 'SNOW') /= 0 .or. &
+ index(field % field, 'SEAICE') /= 0 .or. &
+ index(field % field, 'SKINTEMP') /= 0) then
- if (index(field % field, 'PMSL') == 0 .and. &
- index(field % field, 'PSFC') == 0 .and. &
- index(field % field, 'SOILHGT') == 0) then
+ if (index(field % field, 'SM000010') /= 0 .or. &
+ index(field % field, 'SM010040') /= 0 .or. &
+ index(field % field, 'SM040100') /= 0 .or. &
+ index(field % field, 'SM100200') /= 0 .or. &
+ index(field % field, 'ST000010') /= 0 .or. &
+ index(field % field, 'ST010040') /= 0 .or. &
+ index(field % field, 'ST040100') /= 0 .or. &
+ index(field % field, 'ST100200') /= 0 .or. &
+ index(field % field, 'SNOW') /= 0 .or. &
+ index(field % field, 'SEAICE') /= 0 .or. &
+ index(field % field, 'SKINTEMP') /= 0) then
+ k = 1
+ else if (index(field % field, 'PMSL') == 0 .and. &
+ index(field % field, 'PSFC') == 0 .and. &
+ index(field % field, 'SOILHGT') == 0) then
do k=1,config_nfglevels
if (vert_level(k) == field % xlvl .or. vert_level(k) == -1.0) exit
end do
@@ -2411,6 +2853,83 @@
lonPoints => grid % lonCell % array
destField1d => fg % soilz % array
ndims = 1
+ else if (index(field % field, 'SM000010') /= 0) then
+write(0,*) 'Interpolating SM000010'
+ nInterpPoints = grid % nCells
+ latPoints => grid % latCell % array
+ lonPoints => grid % lonCell % array
+ destField1d => fg % sm000010 % array
+ ndims = 1
+ else if (index(field % field, 'SM010040') /= 0) then
+write(0,*) 'Interpolating SM010040'
+ nInterpPoints = grid % nCells
+ latPoints => grid % latCell % array
+ lonPoints => grid % lonCell % array
+ destField1d => fg % sm010040 % array
+ ndims = 1
+ else if (index(field % field, 'SM040100') /= 0) then
+write(0,*) 'Interpolating SM040100'
+ nInterpPoints = grid % nCells
+ latPoints => grid % latCell % array
+ lonPoints => grid % lonCell % array
+ destField1d => fg % sm040100 % array
+ ndims = 1
+ else if (index(field % field, 'SM100200') /= 0) then
+write(0,*) 'Interpolating SM100200'
+ nInterpPoints = grid % nCells
+ latPoints => grid % latCell % array
+ lonPoints => grid % lonCell % array
+ destField1d => fg % sm100200 % array
+ ndims = 1
+ else if (index(field % field, 'ST000010') /= 0) then
+write(0,*) 'Interpolating ST000010'
+ nInterpPoints = grid % nCells
+ latPoints => grid % latCell % array
+ lonPoints => grid % lonCell % array
+ destField1d => fg % st000010 % array
+ ndims = 1
+ else if (index(field % field, 'ST010040') /= 0) then
+write(0,*) 'Interpolating ST010040'
+ nInterpPoints = grid % nCells
+ latPoints => grid % latCell % array
+ lonPoints => grid % lonCell % array
+ destField1d => fg % st010040 % array
+ ndims = 1
+ else if (index(field % field, 'ST040100') /= 0) then
+write(0,*) 'Interpolating ST040100'
+ nInterpPoints = grid % nCells
+ latPoints => grid % latCell % array
+ lonPoints => grid % lonCell % array
+ destField1d => fg % st040100 % array
+ ndims = 1
+ else if (index(field % field, 'ST100200') /= 0) then
+write(0,*) 'Interpolating ST100200'
+ nInterpPoints = grid % nCells
+ latPoints => grid % latCell % array
+ lonPoints => grid % lonCell % array
+ destField1d => fg % st100200 % array
+ ndims = 1
+ else if (index(field % field, 'SNOW') /= 0) then
+write(0,*) 'Interpolating SNOW'
+ nInterpPoints = grid % nCells
+ latPoints => grid % latCell % array
+ lonPoints => grid % lonCell % array
+ destField1d => fg % snow % array
+ ndims = 1
+ else if (index(field % field, 'SEAICE') /= 0) then
+write(0,*) 'Interpolating SEAICE'
+ nInterpPoints = grid % nCells
+ latPoints => grid % latCell % array
+ lonPoints => grid % lonCell % array
+ destField1d => fg % seaice % array
+ ndims = 1
+ else if (index(field % field, 'SKINTEMP') /= 0) then
+write(0,*) 'Interpolating SKINTEMP'
+ nInterpPoints = grid % nCells
+ latPoints => grid % latCell % array
+ lonPoints => grid % lonCell % array
+ destField1d => fg % skintemp % array
+ ndims = 1
end if
do i=1,nInterpPoints
@@ -2435,7 +2954,19 @@
call read_met_close()
+ ! Fix for isobaric data
+ if (minval(fg % p % array(:,:)) == 0.0) then
+ write(0,*) 'Setting pressure field for isobaric data'
+ do k=1,config_nfglevels
+ if (vert_level(k) /= 200100.0) then
+ fg % p % array(k,:) = vert_level(k)
+ else
+ fg % p % array(k,:) = fg % psfc % array(:)
+ end if
+ end do
+ end if
+
!
! Compute normal wind component and store in fg%u
!
@@ -2450,31 +2981,304 @@
!
! Vertically interpolate meteorological data
!
+ allocate(sorted_arr(2,config_nfglevels))
+
do iCell=1,grid%nCells
+
+ ! T
+ sorted_arr(:,:) = -999.0
+ do k=1,config_nfglevels
+ sorted_arr(1,k) = fg % z % array(k,iCell)
+!NOSFC if (vert_level(k) == 200100.0) sorted_arr(1,k) = fg % soilz % array(iCell)
+ if (vert_level(k) == 200100.0) sorted_arr(1,k) = 99999.0
+ sorted_arr(2,k) = fg % t % array(k,iCell)
+ end do
+ call quicksort(config_nfglevels, sorted_arr)
do k=1,grid%nVertLevels
- ! HGT
- ! PRES
- ! THETA
- ! RH
- ! T
+ target_z = 0.5 * (grid % zgrid % array(k,iCell) + grid % zgrid % array(k+1,iCell))
+ state % theta % array(k,iCell) = vertical_interp(target_z, config_nfglevels, sorted_arr, order=1, extrap=1)
end do
+
+
+ ! RH
+ sorted_arr(:,:) = -999.0
+ do k=1,config_nfglevels
+ sorted_arr(1,k) = fg % z % array(k,iCell)
+!NOSFC if (vert_level(k) == 200100.0) sorted_arr(1,k) = fg % soilz % array(iCell)
+ if (vert_level(k) == 200100.0) sorted_arr(1,k) = 99999.0
+ sorted_arr(2,k) = fg % rh % array(k,iCell)
+ end do
+ call quicksort(config_nfglevels, sorted_arr)
+ do k=1,grid%nVertLevels
+ target_z = 0.5 * (grid % zgrid % array(k,iCell) + grid % zgrid % array(k+1,iCell))
+ state % scalars % array(state % index_qv,k,iCell) = vertical_interp(target_z, config_nfglevels, sorted_arr, order=1, extrap=0)
+ end do
+
+
+ ! GHT
+ sorted_arr(:,:) = -999.0
+ do k=1,config_nfglevels
+ sorted_arr(1,k) = fg % z % array(k,iCell)
+!NOSFC if (vert_level(k) == 200100.0) sorted_arr(1,k) = fg % soilz % array(iCell)
+ if (vert_level(k) == 200100.0) sorted_arr(1,k) = 99999.0
+ sorted_arr(2,k) = fg % z % array(k,iCell)
+ end do
+ call quicksort(config_nfglevels, sorted_arr)
+ do k=1,grid%nVertLevels
+ target_z = 0.5 * (grid % zgrid % array(k,iCell) + grid % zgrid % array(k+1,iCell))
+ fg % gfs_z % array(k,iCell) = vertical_interp(target_z, config_nfglevels, sorted_arr, order=1, extrap=1)
+ end do
+
+
+ ! PRESSURE
+ sorted_arr(:,:) = -999.0
+ do k=1,config_nfglevels
+ sorted_arr(1,k) = fg % z % array(k,iCell)
+ if (vert_level(k) == 200100.0) then
+!NOSFC sorted_arr(1,k) = fg % soilz % array(iCell)
+ sorted_arr(1,k) = 99999.0
+ sfc_k = k
+ end if
+ sorted_arr(2,k) = log(fg % p % array(k,iCell))
+ end do
+ call quicksort(config_nfglevels, sorted_arr)
+ do k=1,grid%nVertLevels
+ target_z = 0.5 * (grid % zgrid % array(k,iCell) + grid % zgrid % array(k+1,iCell))
+ diag % pressure_base % array(k,iCell) = exp(vertical_interp(target_z, config_nfglevels, sorted_arr, order=1, extrap=1))
+ end do
+
+
+ ! PRESSURE
+ sorted_arr(:,:) = -999.0
+ do k=1,config_nfglevels
+ sorted_arr(1,k) = fg % z % array(k,iCell)
+ if (vert_level(k) == 200100.0) then
+!NOSFC sorted_arr(1,k) = fg % soilz % array(iCell)
+ sorted_arr(1,k) = 99999.0
+ sfc_k = k
+ end if
+ sorted_arr(2,k) = log(fg % p % array(k,iCell))
+ end do
+ call quicksort(config_nfglevels, sorted_arr)
+ do k=1,grid%nVertLevels+1
+ target_z = grid % zgrid % array(k,iCell)
+ fg % gfs_p % array(k,iCell) = exp(vertical_interp(target_z, config_nfglevels, sorted_arr, order=1, extrap=1))
+ end do
+
end do
+
do iEdge=1,grid%nEdges
+
+ ! U
+ sorted_arr(:,:) = -999.0
+ do k=1,config_nfglevels
+ sorted_arr(1,k) = 0.5 * (fg % z % array(k,cellsOnEdge(1,iEdge)) + fg % z % array(k,cellsOnEdge(2,iEdge)))
+!NOSFC if (vert_level(k) == 200100.0) sorted_arr(1,k) = 0.5 * (fg % soilz % array(cellsOnEdge(1,iEdge)) + fg % soilz % array(cellsOnEdge(2,iEdge)))
+ if (vert_level(k) == 200100.0) sorted_arr(1,k) = 99999.0
+ sorted_arr(2,k) = fg % u % array(k,iEdge)
+ end do
+ call quicksort(config_nfglevels, sorted_arr)
do k=1,grid%nVertLevels
- ! U
+ target_z = 0.25 * (grid % zgrid % array(k,cellsOnEdge(1,iEdge)) + grid % zgrid % array(k+1,cellsOnEdge(1,iEdge)) + grid % zgrid % array(k,cellsOnEdge(2,iEdge)) + grid % zgrid % array(k+1,cellsOnEdge(2,iEdge)))
+ state % u % array(k,iEdge) = vertical_interp(target_z, config_nfglevels, sorted_arr, order=1, extrap=0)
end do
+
end do
!
+ ! Adjust surface pressure for difference in topography
+ !
+ do sfc_k=1,config_nfglevels
+ if (vert_level(sfc_k) == 200100.) exit
+ end do
+ do iCell=1,grid%nCells
+
+ ! We need to extrapolate
+ sorted_arr(:,:) = -999.0
+ do k=1,config_nfglevels
+ sorted_arr(1,k) = fg % z % array(k,iCell)
+ if (vert_level(k) == 200100.0) then
+!NOSFC sorted_arr(1,k) = fg % soilz % array(iCell)
+ sorted_arr(1,k) = 99999.0
+ end if
+ sorted_arr(2,k) = log(fg % p % array(k,iCell))
+ end do
+ call quicksort(config_nfglevels, sorted_arr)
+ target_z = grid % zgrid % array(1,iCell)
+ fg % psfc % array(iCell) = exp(vertical_interp(target_z, config_nfglevels, sorted_arr, order=1, extrap=1))
+
+ end do
+
+ deallocate(sorted_arr)
+
+
+
+ !
! Diagnose fields needed in initial conditions file (u, w, rho, theta, scalars)
!
+ do iCell=1,grid%nCells
+ do k=1,grid%nVertLevels
+ es = 6.112 * exp((17.27*(state % theta % array(k,iCell) - 273.16))/(state % theta % array(k,iCell) - 35.86))
+ rs = 0.622 * es / (diag % pressure_base % array(k,iCell) - es)
+ state % scalars % array(state % index_qv,k,iCell) = rs * state % scalars % array(state % index_qv,k,iCell)
- deallocate(vert_level)
+ rgas_moist = (1.0 + 0.61 * state % scalars % array(state % index_qv,k,iCell)) * rgas
+ state % rho % array(k,iCell) = diag % pressure_base % array(k,iCell) / (rgas_moist * state % theta % array(k,iCell))
+ state % rho % array(k,iCell) = state % rho % array(k,iCell) / (1.0 + state % scalars % array(state % index_qv,k,iCell))
+ state % rho % array(k,iCell) = state % rho % array(k,iCell) / zz(k,iCell)
+
+ state % theta % array(k,iCell) = state % theta % array(k,iCell) * (p0 / diag % pressure_base % array(k,iCell)) ** (rgas / cp)
+ state % theta % array(k,iCell) = state % theta % array(k,iCell) * (1.0 + 1.61*state % scalars % array(state % index_qv,k,iCell))
+ end do
+ end do
+
+#if 0
+ do iCell=1,grid%nCells
+
+ fsum(:) = 0.0
+ do k=2,grid%nVertLevels+1
+ rgas_moist = (1.0 + 0.61 * state % scalars % array(state % index_qv,k-1,iCell)) * rgas
+!RGAS rgas_moist = rgas
+ ptmp = fg % gfs_p % array(k,iCell)
+ state % theta % array(k-1,iCell) = -gravity * (zgrid(k,iCell) - zgrid(k-1,iCell)) / (rgas_moist * (log(ptmp / fg % gfs_p % array(1,iCell)) + fsum(k-1)))
+ fsum(k) = fsum(k-1) + gravity * (zgrid(k,iCell) - zgrid(k-1,iCell)) / (rgas_moist * state % theta % array(k-1,iCell))
+
+ state % theta % array(k-1,iCell) = state % theta % array(k-1,iCell) * (p0 / diag % pressure_base % array(k-1,iCell)) ** (rgas / cp)
+ state % theta % array(k-1,iCell) = state % theta % array(k-1,iCell) * (1.0 + 1.61*state % scalars % array(state % index_qv,k-1,iCell))
+ end do
+ end do
+
#endif
+#if 0
+ do iCell=1,grid%nCells
+ do k=2,grid%nVertLevels+1
+ rgas_moist = (1.0 + 0.61 * state % scalars % array(state % index_qv,k-1,iCell)) * rgas
+ state % theta % array(k-1,iCell) = -gravity * (zgrid(k,iCell) - zgrid(k-1,iCell)) / (rgas_moist * log(fg % gfs_p % array(k,iCell) / fg % gfs_p % array(k-1,iCell)))
+ state % theta % array(k-1,iCell) = state % theta % array(k-1,iCell) * (p0 / diag % pressure_base % array(k-1,iCell)) ** (rgas / cp)
+ state % theta % array(k-1,iCell) = state % theta % array(k-1,iCell) * (1.0 + 1.61*state % scalars % array(state % index_qv,k-1,iCell))
+ end do
+ end do
+#endif
+
+
+
+ !
+ ! Reference state based on a dry isothermal atmosphere
+ !
+ do iCell=1,grid % nCells
+ do k=1,nz1
+ ztemp = 0.5*(zgrid(k+1,iCell)+zgrid(k,iCell))
+ ppb(k,iCell) = p0*exp(-gravity*ztemp/(rgas*t0b)) ! pressure_base
+ pb (k,iCell) = (ppb(k,iCell)/p0)**(rgas/cp) ! exner_base
+ rb (k,iCell) = ppb(k,iCell)/(rgas*t0b*zz(k,iCell)) ! rho_base
+ tb (k,iCell) = t0b/pb(k,iCell) ! theta_base
+ rtb(k,iCell) = rb(k,iCell)*tb(k,iCell) ! rtheta_base
+ p (k,iCell) = pb(k,iCell) ! exner
+ pp (k,iCell) = 0. ! pressure_p
+ rr (k,iCell) = 0. ! rho_p
+ end do
+ end do
+
+ do iEdge=1,grid%nEdges
+ do k=1,grid%nVertLevels
+ diag % ru % array(k,iEdge) = state % u % array(k,iEdge) * 0.5*(state % rho % array(k,cellsOnEdge(1,iEdge)) + state % rho % array(k,cellsOnEdge(2,iEdge)))
+ end do
+ end do
+
+
+ ! For z-metric term in omega equation
+ do iEdge = 1,grid % nEdges
+ cell1 = CellsOnEdge(1,iEdge)
+ cell2 = CellsOnEdge(2,iEdge)
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
+
+ do k = 1, grid%nVertLevels
+
+ if (config_theta_adv_order == 2) then
+
+ z_edge = (zgrid(k,cell1)+zgrid(k,cell2))/2.
+
+ else !theta_adv_order == 3 or 4
+
+ d2fdx2_cell1 = deriv_two(1,1,iEdge) * zgrid(k,cell1)
+ d2fdx2_cell2 = deriv_two(1,2,iEdge) * zgrid(k,cell2)
+ do i=1, grid % nEdgesOnCell % array (cell1)
+ if ( grid % CellsOnCell % array (i,cell1) > 0) &
+ d2fdx2_cell1 = d2fdx2_cell1 + deriv_two(i+1,1,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell1))
+ end do
+ do i=1, grid % nEdgesOnCell % array (cell2)
+ if ( grid % CellsOnCell % array (i,cell2) > 0) &
+ d2fdx2_cell2 = d2fdx2_cell2 + deriv_two(i+1,2,iEdge) * zgrid(k,grid % CellsOnCell % array (i,cell2))
+ end do
+
+ z_edge = 0.5*(zgrid(k,cell1) + zgrid(k,cell2)) &
+ - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 + d2fdx2_cell2) / 12.
+
+ if (config_theta_adv_order == 3) then
+ z_edge3 = - (grid % dcEdge % array(iEdge) **2) * (d2fdx2_cell1 - d2fdx2_cell2) / 12.
+ else
+ z_edge3 = 0.
+ end if
+
+ end if
+
+ zb(k,1,iEdge) = (z_edge-zgrid(k,cell1))*dvEdge(iEdge)/AreaCell(cell1)
+ zb(k,2,iEdge) = (z_edge-zgrid(k,cell2))*dvEdge(iEdge)/AreaCell(cell2)
+ zb3(k,1,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell1)
+ zb3(k,2,iEdge)= z_edge3*dvEdge(iEdge)/AreaCell(cell2)
+
+ if (k /= 1) then
+ zf(k,1,iEdge) = ( fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1) )*zb(k,1,iEdge)
+ zf(k,2,iEdge) = ( fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2) )*zb(k,2,iEdge)
+ zf3(k,1,iEdge)= ( fzm(k)*zz(k,cell1)+fzp(k)*zz(k-1,cell1) )*zb3(k,1,iEdge)
+ zf3(k,2,iEdge)= ( fzm(k)*zz(k,cell2)+fzp(k)*zz(k-1,cell2) )*zb3(k,2,iEdge)
+ end if
+
+ end do
+
+ end if
+ end do
+
+ diag % rw % array = 0.
+ state % w % array = 0.
+ do iEdge = 1,grid % nEdges
+
+ cell1 = CellsOnEdge(1,iEdge)
+ cell2 = CellsOnEdge(2,iEdge)
+
+ if (cell1 <= nCellsSolve .or. cell2 <= nCellsSolve ) then
+ do k = 2, grid%nVertLevels
+ flux = (fzm(k)*diag % ru % array(k,iEdge)+fzp(k)*diag % ru % array(k-1,iEdge))
+ diag % rw % array(k,cell2) = diag % rw % array(k,cell2) + zf(k,2,iEdge)*flux
+ diag % rw % array(k,cell1) = diag % rw % array(k,cell1) - zf(k,1,iEdge)*flux
+
+ if (config_theta_adv_order ==3) then
+ diag % rw % array(k,cell2) = diag % rw % array(k,cell2) &
+ - sign(1.,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,2,iEdge)*flux
+ diag % rw % array(k,cell1) = diag % rw % array(k,cell1) &
+ + sign(1.,diag % ru % array(k,iEdge))*config_coef_3rd_order*zf3(k,1,iEdge)*flux
+ end if
+
+ end do
+ end if
+
+ end do
+
+ ! Compute w from rho and rw
+ do iCell=1,grid%nCells
+ do k=2,grid%nVertLevels
+ state % w % array(k,iCell) = diag % rw % array(k,iCell) &
+ / (fzp(k) * state % rho % array(k-1,iCell) + fzm(k) * state % rho % array(k,iCell))
+ end do
+ end do
+
+ deallocate(vert_level)
+
+
end subroutine nhyd_test_case_gfs
@@ -2640,4 +3444,91 @@
end function nearest_edge
+
+ real function vertical_interp(target_z, nz, zf, order, extrap, surface_val, sealev_val)
+
+ implicit none
+
+ real (kind=RKIND), intent(in) :: target_z
+ integer, intent(in) :: nz
+ real (kind=RKIND), dimension(2,nz), intent(in) :: zf ! zf(1,:) is column of vertical coordinate values, zf(2,:) is column of field values
+ integer, intent(in), optional :: order
+ integer, intent(in), optional :: extrap
+ real (kind=RKIND), intent(in), optional :: surface_val
+ real (kind=RKIND), intent(in), optional :: sealev_val
+
+ integer :: k, lm, lp
+ real (kind=RKIND) :: wm, wp
+ real (kind=RKIND) :: slope
+
+ integer :: interp_order, extrap_type
+ real (kind=RKIND) :: surface, sealevel
+
+
+ if (present(order)) then
+ interp_order = order
+ else
+ interp_order = 2
+ end if
+
+ if (present(extrap)) then
+ extrap_type = extrap
+ else
+ interp_order = 1
+ end if
+
+ if (present(surface_val)) then
+ surface = surface_val
+ else
+ surface = 200100.0
+ end if
+
+ if (present(sealev_val)) then
+ sealevel = sealev_val
+ else
+ sealevel = 201300.0
+ end if
+
+ !
+ ! Extrapolation required
+ !
+ if (target_z < zf(1,1)) then
+ if (extrap_type == 0) then
+ vertical_interp = zf(2,1)
+ else if (extrap_type == 1) then
+ slope = (zf(2,2) - zf(2,1)) / (zf(1,2) - zf(1,1))
+ vertical_interp = zf(2,1) + slope * (target_z - zf(1,1))
+ end if
+ return
+ end if
+ if (target_z >= zf(1,nz)) then
+ if (extrap_type == 0) then
+ vertical_interp = zf(2,nz)
+ else if (extrap_type == 1) then
+ slope = (zf(2,nz) - zf(2,nz-1)) / (zf(1,nz) - zf(1,nz-1))
+ vertical_interp = zf(2,nz) + slope * (target_z - zf(1,nz))
+ end if
+ return
+ end if
+
+
+ !
+ ! No extrapolation required
+ !
+ do k=1,nz-1
+ if (target_z >= zf(1,k) .and. target_z < zf(1,k+1)) then
+ lm = k
+ lp = k+1
+ wm = (zf(1,k+1) - target_z) / (zf(1,k+1) - zf(1,k))
+ wp = (target_z - zf(1,k)) / (zf(1,k+1) - zf(1,k))
+ exit
+ end if
+ end do
+
+ vertical_interp = wm*zf(2,lm) + wp*zf(2,lp)
+
+ return
+
+ end function vertical_interp
+
end module test_cases
</font>
</pre>