<p><b>duda</b> 2012-07-26 17:46:47 -0600 (Thu, 26 Jul 2012)</p><p>BRANCH COMMIT<br>
<br>
Correct coordinate surface smoothing for parallel execution, and control <br>
this smoothing by the existing namelist parameter config_smooth_surfaces.<br>
<br>
<br>
M src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F<br>
</p><hr noshade><pre><font color="gray">Modified: branches/dcmip/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F
===================================================================
--- branches/dcmip/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-07-26 23:39:51 UTC (rev 2067)
+++ branches/dcmip/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-07-26 23:46:47 UTC (rev 2068)
@@ -4642,22 +4642,30 @@
!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
- real (kind=RKIND), dimension(:), pointer :: dvEdge, AreaCell, xCell, yCell
+ integer, dimension(:), pointer :: nEdgesOnEdge, nEdgesOnCell
+ integer, dimension(:,:), pointer :: edgesOnEdge, CellsOnEdge, edgesOnCell
+ real (kind=RKIND), dimension(:), pointer :: dvEdge, AreaCell, xCell, yCell, dcEdge
real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge
- integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, nz, nz1, itr, itrp, cell1, cell2, nCellsSolve
+ integer :: iCell, iCell1, iCell2 , iEdge, vtx1, vtx2, ivtx, i, k, kz, nz, nz1, itr, itrp, cell1, cell2, nCellsSolve
integer :: index_qv
real (kind=RKIND) :: ptop, p0, pis, flux, d2fdx2_cell1, d2fdx2_cell2
real (kind=RKIND) :: ztemp, zd, zt, dz, str
+ real(kind=RKIND), dimension(:), pointer :: hs, hs1
real (kind=RKIND), dimension(grid % nVertLevels, grid % nCells) :: rh
real (kind=RKIND) :: es, qvs, xnutr, ptemp
- integer :: iter
+ integer :: iter, nsm
+ integer, dimension(:,:), pointer :: cellsOnCell
+ type (field1DReal):: tempField
+
+ type (block_type), pointer :: block
+ type (parallel_info), pointer :: parinfo
+ type (dm_info), pointer :: dminfo
+
real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: zc, zw, ah
real (kind=RKIND), dimension(grid % nVertLevels ) :: zu, dzw, rdzwp, rdzwm
real (kind=RKIND), allocatable, dimension(:) :: psiVertex
@@ -4666,12 +4674,18 @@
real (kind=RKIND) :: um, us, rcp, rcv
real (kind=RKIND) :: xmid, temp, pres, a_scale, xac, xlac, shear, tsurf, usurf
- real (kind=RKIND) :: xi, yi, ri, xa, xc, yc, xla, zinv, xn2, xn2m, xn2l, sm, dzh, dzht, dzmin, z_edge, z_edge3
+ real (kind=RKIND) :: xi, yi, ri, xa, xc, yc, xla, zinv, xn2, xn2m, xn2l, sm, dzh, dzht, dzmin, dzmina, dzminf, &
+ dzmina_global, z_edge, z_edge3, sm0
integer, dimension(grid % nCells, 2) :: next_cell
real (kind=RKIND), dimension(grid % nCells) :: hxzt, pitop, ptopb
logical, parameter :: terrain_smooth = .false.
+ block => grid % block
+ parinfo => block % parinfo
+ dminfo => block % domain % dminfo
+
+
!
! Scale all distances
!
@@ -4695,9 +4709,13 @@
weightsOnEdge => grid % weightsOnEdge % array
nEdgesOnEdge => grid % nEdgesOnEdge % 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
+ nEdgesOnCell => grid % nEdgesOnCell % array
deriv_two => grid % deriv_two % array
nz1 = grid % nVertLevels
@@ -4861,72 +4879,105 @@
hx(nz,iCell) = zt
-!***** SHP -> get the temporary point information for the neighbor cell ->> should be changed!!!!!
+ enddo
+ write(0,*) ' hx computation complete '
+
!!! MGD WE NEED TO REPLACE THIS TERRAIN SMOOTHING WITH TC9
-if (terrain_smooth) then
- do i=1,grid % nCells
- !option 1
- !IF(yCell(i).eq.yCell(iCell).and.xCell(i).eq.xCell(iCell)-sqrt(3.)*grid % dcEdge % array(1)) next_cell(iCell,1) = i
- !IF(yCell(i).eq.yCell(iCell).and.xCell(i).eq.xCell(iCell)+sqrt(3.)*grid % dcEdge % array(1)) next_cell(iCell,2) = i
- !option 2
- next_cell(iCell,1) = iCell - 8 ! note ny=4
- next_cell(iCell,2) = iCell + 8 ! note ny=4
+ kz = nz
- if (xCell(iCell).le. 3.*grid % dcEdge % array(1)) then
- next_cell(iCell,1) = 1
- else if (xCell(iCell).ge. maxval(xCell(:))-3.*grid % dcEdge % array(1)) then
- next_cell(iCell,2) = 1
- end if
+ if (config_smooth_surfaces) then
- end do
-end if
- enddo
- write(0,*) ' hx computation complete '
+ write(0,*) ' '
+ write(0,*) ' Smoothing vertical coordinate surfaces'
+ write(0,*) ' '
+ allocate(hs (grid % nCells+1))
+ allocate(hs1(grid % nCells+1))
-! smoothing grid for the upper level >> but not propoer for parallel programing
-if (terrain_smooth) then
- dzmin=.7
- do k=2,nz1
- sm = .25*min((zc(k)-zc(k-1))/dz,1.0_RKIND)
- do i=1,grid % nCells
- hx(k,i) = hx(k-1,i)
- end do
+ dzmin = 0.5
+ sm0 = 0.5
+ nsm = 30
- do iter = 1,20 !iteration for smoothing
+ write(6,*) 'dzmin = ',dzmin,' sm0 = ',sm0,' nsm = ',nsm
- do i=1,grid % nCells
- hxzt(i) = hx(k,i) + sm*(hx(k,next_cell(i,2))-2.*hx(k,i)+hx(k,next_cell(i,1)))
- end do
- dzh = zc(k) - zc(k-1)
- do i=1,grid % nCells
- dzht = zc(k)+hxzt(i) - zc(k-1)-hx(k-1,i)
- if(dzht.lt.dzh) dzh = dzht
- end do
+ do k=2,kz-1
+ hx(k,:) = hx(k-1,:)
+ dzminf = zw(k)-zw(k-1)
- if(dzh.gt.dzmin*(zc(k)-zc(k-1))) then
- do i=1,grid % nCells
- hx(k,i) = hxzt(i)
+! dzmin = max(0.5_RKIND,1.-.5*zw(k)/hm)
+
+ sm = sm0*max( min(.5*zw(k)/hm,1.0_RKIND), .05 )
+
+ do i=1,nsm
+ 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
- else
- goto 99 !SHP - this algorithm should be changed
- end if
- end do !end of iteration for smoothing
-99 write(0,*) "PASS-SHP"
- end do
-end if
+ tempField % block => block
+ tempField % dimSizes(1) = grid % nCells
+ tempField % sendList => parinfo % cellsToSend
+ tempField % recvList => parinfo % cellsToRecv
+ tempField % array => hs
+ call mpas_dmpar_exch_halo_field(tempField)
+
+ ! dzmina = minval(hs(:)-hx(k-1,:))
+ dzmina = minval(zw(k)+ah(k)*hs(1:grid%nCellsSolve)-zw(k-1)-ah(k-1)*hx(k-1,1:grid%nCellsSolve))
+ call mpas_dmpar_min_real(dminfo, dzmina, dzmina_global)
+ ! write(0,*) ' k,i, dzmina, dzmin, zw(k)-zw(k-1) ', k,i, dzmina, dzmin, zw(k)-zw(k-1)
+ if (dzmina_global >= dzmin*(zw(k)-zw(k-1))) then
+ hx(k,:)=hs(:)
+ dzminf = dzmina_global
+ else
+ exit
+ end if
+ end do
+ write(0,*) 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
+
+ deallocate(hs )
+ deallocate(hs1)
+
+ else
+
+ do k=2,nz1
+ dzmina = minval(zw(k)+ah(k)*hx(k,:)-zw(k-1)-ah(k-1)*hx(k-1,:))
+ write(0,*) k,dzmina/(zw(k)-zw(k-1))
+ end do
+
+ end if
+
+
do iCell=1,grid % nCells
do k=1,nz
- if (terrain_smooth) then
- zgrid(k,iCell) = ah(k)*(zc(k)*(1.-hx(k,iCell)/zt)+hx(k,iCell)) &
- + (1.-ah(k)) * zc(k)
+ if (config_smooth_surfaces) then
+ zgrid(k,iCell) = ah(k)*(zc(k)*(1.-hx(k,iCell)/zt)+hx(k,iCell)) &
+ + (1.-ah(k)) * zc(k)
else
- zgrid(k,iCell) = ah(k)*(zc(k)*(1.-hx(1,iCell)/zt)+hx(1,iCell)) &
- + (1.-ah(k)) * zc(k)
+ zgrid(k,iCell) = ah(k)*(zc(k)*(1.-hx(1,iCell)/zt)+hx(1,iCell)) &
+ + (1.-ah(k)) * zc(k)
end if
end do
do k=1,nz1
@@ -5363,7 +5414,7 @@
integer :: index_qv
real (kind=RKIND) :: ptop, p0, pis, flux, d2fdx2_cell1, d2fdx2_cell2
- real(kind=RKIND), allocatable, dimension(:) :: hs, hs1
+ real(kind=RKIND), dimension(:), pointer :: hs, hs1
real (kind=RKIND) :: ztemp, zd, zt, dz, str, zh, hmax
@@ -5371,21 +5422,32 @@
real (kind=RKIND) :: es, qvs, xnutr, ptemp
integer :: iter, nsm, kz
+ type (field1DReal):: tempField
+
+ type (block_type), pointer :: block
+ type (parallel_info), pointer :: parinfo
+ type (dm_info), pointer :: dminfo
+
real (kind=RKIND), dimension(grid % nVertLevels + 1 ) :: zc, zw, ah
real (kind=RKIND), dimension(grid % nVertLevels ) :: zu, dzw, rdzwp, rdzwm
real (kind=RKIND), allocatable, dimension(:) :: psiVertex
real (kind=RKIND) :: d1, d2, d3, cof1, cof2, cf1, cf2, cf3
real (kind=RKIND) :: um, us, rcp, rcv, gamma, xa, zinb, zint, tinv, th_inb, th_int
- real (kind=RKIND) :: xmid, temp, pres, a_scale, rad, shear, tsurf, usurf, sm0, dzmina, dzminf
+ real (kind=RKIND) :: xmid, temp, pres, a_scale, rad, shear, tsurf, usurf, sm0, dzmina, dzmina_global, dzminf
real (kind=RKIND) :: xi, yi, r1m, r2m, xc, yc, xla, zinv, xn2, xn2m, xn2l, sm, dzh, dzht, dzmin, z_edge, z_edge3
integer, dimension(grid % nCells, 2) :: next_cell
real (kind=RKIND), dimension(grid % nCells) :: pitop, ptopb
- logical, parameter :: terrain_smooth = .false. , hybrid = .false.
-! logical, parameter :: terrain_smooth = .true. , hybrid = .true.
+ logical, parameter :: hybrid = .false.
+! logical, parameter :: hybrid = .true.
+ block => grid % block
+ parinfo => block % parinfo
+ dminfo => block % domain % dminfo
+
+
!
! Scale all distances
!
@@ -5599,29 +5661,31 @@
hmax = maxval(hx(1,:))
write(6,*) "max terrain height = ",hmax
- if (terrain_smooth) then
+ if (config_smooth_surfaces) then
+ write(0,*) ' '
+ write(0,*) ' Smoothing vertical coordinate surfaces'
+ write(0,*) ' '
+
allocate(hs (grid % nCells+1))
allocate(hs1(grid % nCells+1))
- dzmin = .5
-
- sm0 = .05
-
+ dzmin = 0.5
+ sm0 = 0.5
nsm = 30
write(6,*) 'dzmin = ',dzmin,' sm0 = ',sm0,' nsm = ',nsm
do k=2,kz-1
-
hx(k,:) = hx(k-1,:)
dzminf = zw(k)-zw(k-1)
- sm = sm0*max( min(.5*zw(k)/hm,1.), .05 )
+! dzmin = max(0.5_RKIND,1.-.5*zw(k)/hm)
+ sm = sm0*max( min(.5*zw(k)/hm,1.0_RKIND), .05 )
+
do i=1,nsm
-
- do iCell=1,grid %nCells
+ do iCell=1,grid % nCells
hs1(iCell) = 0.
do j = 1,nEdgesOnCell(iCell)
@@ -5632,26 +5696,35 @@
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
+ ! 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(zw(k)+ah(k)*hs(:)-zw(k-1)-ah(k-1)*hx(k-1,:))
- if(dzmina.ge.dzmin*(zw(k)-zw(k-1))) then
+
+ tempField % block => block
+ tempField % dimSizes(1) = grid % nCells
+ tempField % sendList => parinfo % cellsToSend
+ tempField % recvList => parinfo % cellsToRecv
+ tempField % array => hs
+
+ call mpas_dmpar_exch_halo_field(tempField)
+
+ ! dzmina = minval(hs(:)-hx(k-1,:))
+ dzmina = minval(zw(k)+ah(k)*hs(1:grid%nCellsSolve)-zw(k-1)-ah(k-1)*hx(k-1,1:grid%nCellsSolve))
+ call mpas_dmpar_min_real(dminfo, dzmina, dzmina_global)
+ ! write(0,*) ' k,i, dzmina, dzmin, zw(k)-zw(k-1) ', k,i, dzmina, dzmin, zw(k)-zw(k-1)
+ if (dzmina_global >= dzmin*(zw(k)-zw(k-1))) then
hx(k,:)=hs(:)
- dzminf = dzmina
+ dzminf = dzmina_global
else
- go to 95
+ exit
end if
end do
- 95 continue
-! write(6,*) k,i-1,sm,dzminf/(zw(k)-zw(k-1)),dzmina/(zw(k)-zw(k-1))
- write(6,30) k,zw(k),i-1,sm,dzminf/(zw(k)-zw(k-1)),dzmina/(zw(k)-zw(k-1))
- 30 format(I5,F10.2,I5,3F9.5)
+ write(0,*) k,i,sm,dzminf/(zw(k)-zw(k-1)),dzmina/(zw(k)-zw(k-1))
end do
do k=kz,nz
@@ -5665,11 +5738,12 @@
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))
+ write(0,*) k,dzmina/(zw(k)-zw(k-1))
end do
end if
+
do iCell=1,grid % nCells
do k=1,nz
zgrid(k,iCell) = zw(k) + ah(k)*hx(k,iCell)
</font>
</pre>