|
XBeach
|
00001 module wetcells_module 00002 implicit none 00003 save 00004 private 00005 public compute_wetcells 00006 contains 00007 subroutine compute_wetcells(s,par) 00008 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00009 ! Copyright (C) 2007 UNESCO-IHE, WL|Delft Hydraulics and Delft University ! 00010 ! Dano Roelvink, Ap van Dongeren, Ad Reniers, Jamie Lescinski, ! 00011 ! Jaap van Thiel de Vries, Robert McCall ! 00012 ! ! 00013 ! d.roelvink@unesco-ihe.org ! 00014 ! UNESCO-IHE Institute for Water Education ! 00015 ! P.O. Box 3015 ! 00016 ! 2601 DA Delft ! 00017 ! The Netherlands ! 00018 ! ! 00019 ! This library is free software; you can redistribute it and/or ! 00020 ! modify it under the terms of the GNU Lesser General Public ! 00021 ! License as published by the Free Software Foundation; either ! 00022 ! version 2.1 of the License, or (at your option) any later version. ! 00023 ! ! 00024 ! This library is distributed in the hope that it will be useful, ! 00025 ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! 00026 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! 00027 ! Lesser General Public License for more details. ! 00028 ! ! 00029 ! You should have received a copy of the GNU Lesser General Public ! 00030 ! License along with this library; if not, write to the Free Software ! 00031 ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 ! 00032 ! USA ! 00033 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00034 use params, only: parameters 00035 use spaceparamsdef, only: spacepars 00036 use paramsconst 00037 00038 IMPLICIT NONE 00039 00040 type(spacepars),target :: s 00041 type(parameters),intent(in) :: par 00042 00043 integer :: i 00044 integer :: j 00045 integer :: upwinddist 00046 integer,dimension(:,:),allocatable,save :: weteb 00047 00048 if(.not.allocated(weteb)) allocate(weteb(s%nx+1,s%ny+1)) 00049 00050 ! wwvv in the next lines 00051 ! hu(i,j) is more or less a function of hh(i+1,j) 00052 ! In the parallel case, this needs some action because 00053 ! for processes not at the bottom, the last row of 00054 ! hu (hu(nx+1,:)) has to be taken from the neighbour below 00055 ! hu is used later on in this subroutine, so we have to insert 00056 ! an mpi call. 00057 ! The same for hum 00058 ! Of course, no action is necessary if hu(nx+1:) is never used... 00059 do j=1,s%ny+1 00060 do i=1,s%nx+1 !Ap 00061 ! Water depth in u-points do momentum equation: mean 00062 !! 00063 !! ARBJ: mean water depth or weighted water depth? How to deal with this in curvi-linear? 00064 !! 00065 s%hum(i,j)=max(.5d0*(s%hh(i,j)+s%hh(min(s%nx,i)+1,j)),par%eps) 00066 end do 00067 end do 00068 ! wwvv here the mpi code to communicate a row of hu 00069 ! we send to the neighbour above and receive from the neighbour 00070 ! below: 00071 ! Wetting and drying criterion (only do momentum balance) 00072 do j=1,s%ny+1 00073 do i=1,s%nx+1 00074 if(s%hu(i,j)>par%eps .and. s%hum(i,j)>par%eps) then ! Jaap and Pieter: If you want to compute correct advection term 00075 s%wetu(i,j)=1 ! then both s%hu and s%hum should be larger than par%eps. It is not 00076 else ! necessarily true that if s%hu>par%eps also s%hum>par%eps. 00077 s%wetu(i,j)=0 00078 end if 00079 end do 00080 end do 00081 ! wwvv about the same for hv, only in the left-right direction 00082 ! hv(i,j) is more or less a function of hh(i,j+1) 00083 ! so in the parallel case, hv(:,ny+1) has to be collected 00084 ! from the right neighbour 00085 ! the same for hvm 00086 do j=1,s%ny+1 00087 do i=1,s%nx+1 00088 ! Water depth in v-points do momentum equation: mean 00089 s%hvm(i,j)=max(.5d0*(s%hh(i,j)+s%hh(i,min(s%ny,j)+1)),par%eps) 00090 end do 00091 end do 00092 ! Wetting and drying criterion (only do momentum balance) 00093 do j=1,s%ny+1 00094 do i=1,s%nx+1 00095 if(s%hv(i,j)>par%eps .and. s%hvm(i,j)>par%eps) then 00096 s%wetv(i,j)=1 00097 else 00098 s%wetv(i,j)=0 00099 end if 00100 end do 00101 end do 00102 ! Jaap Wetting and drying criterion eta points 00103 do j=1,s%ny+1 00104 do i=1,s%nx+1 00105 !A eta point is wet if any of the surrounding velocity points is wet... 00106 ! wetz(i,j) = min(1,wetu(max(i,2)-1,j)+wetu(i,j)+wetv(i,j)+wetv(i,max(j,2)-1)) 00107 if(s%hh(i,j)>par%eps) then 00108 s%wetz(i,j)=1 00109 else 00110 s%wetz(i,j)=0 00111 end if 00112 end do 00113 end do 00114 00115 ! Wetting and drying for wave module 00116 where(s%hh+par%delta*s%H>par%eps .or. s%wetz==1) 00117 s%wete = 1 00118 elsewhere 00119 s%wete = 0 00120 endwhere 00121 ! also need to surround wet cells 00122 weteb = 0 00123 ! upwind distance in wave propagation routine 00124 if (par%scheme==SCHEME_UPWIND_1) then 00125 upwinddist = 1 00126 else 00127 upwinddist = 2 00128 endif 00129 forall(i=1:s%nx+1,j=1:s%ny+1,s%wete(i,j)==0) 00130 weteb(i,j) = max( maxval( s%wetz( max(i-upwinddist,1):min(i+upwinddist,s%nx+1) ,j ) ), & 00131 maxval( s%wetz( i, max(j-upwinddist,1):min(j+upwinddist,s%ny+1) ) ), & 00132 maxval( s%wetu( max(i-1,1):min(i,s%nx+1) ,j ) ), & 00133 maxval( s%wetv( i, max(j-1,1):min(j,s%ny+1) ) ) & 00134 ) 00135 !weteb(i,j) = maxval(s%wete(max(i-2,1):min(i+2,s%nx+1),max(j-2,1):min(j+2,s%ny+1)) 00136 !weteb(i,j) = max( maxval( s%wete( max(i-2,1):min(i+2,s%nx+1) ,j ) ), & 00137 ! maxval( s%wete( i ,max(j-2,1):min(j+2,s%ny+1) ) ) & 00138 ! ) 00139 endforall 00140 where(s%wete==0) 00141 s%wete=weteb 00142 endwhere 00143 00144 end subroutine compute_wetcells 00145 end module wetcells_module
1.7.6.1
