module morphevolution implicit none save #ifdef HAVE_CONFIG_H #include "config.h" #endif contains subroutine transus(s,par) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Copyright (C) 2007 UNESCO-IHE, WL|Delft Hydraulics and Delft University ! ! Dano Roelvink, Ap van Dongeren, Ad Reniers, Jamie Lescinski, ! ! Jaap van Thiel de Vries, Robert McCall ! ! ! ! d.roelvink@unesco-ihe.org ! ! UNESCO-IHE Institute for Water Education ! ! P.O. Box 3015 ! ! 2601 DA Delft ! ! The Netherlands ! ! ! ! This library is free software; you can redistribute it and/or ! ! modify it under the terms of the GNU Lesser General Public ! ! License as published by the Free Software Foundation; either ! ! version 2.1 of the License, or (at your option) any later version. ! ! ! ! This library is distributed in the hope that it will be useful, ! ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! ! Lesser General Public License for more details. ! ! ! ! You should have received a copy of the GNU Lesser General Public ! ! License along with this library; if not, write to the Free Software ! ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 ! ! USA ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use params use spaceparams use xmpi_module use interp use paramsconst ! use vsmumod implicit none type(spacepars),target :: s type(parameters) :: par integer :: i,isig,numb integer :: j,jg real*8 :: exp_ero real*8,dimension(:),allocatable,save :: chain,cumchain real*8,dimension(:,:),allocatable,save :: vmag2,uau,uav,um,vm,ueu_sed,uev_sed,veu_sed,vev_sed real*8,dimension(:,:),allocatable,save :: ccvt,dcdz,dsigt,aref real*8,dimension(:,:),allocatable,save :: cc,ccb,cu,cv,Sus,Svs real*8,dimension(:,:),allocatable,save :: cub,cvb,Sub,Svb,pbbedu,pbbedv real*8,dimension(:,:),allocatable,save :: suq3d,svq3d,eswmax,eswbed,sigs,deltas real*8,dimension(:,:,:),allocatable,save :: dsig,ccv,sdif,cuq3d,cvq3d,fac real*8,dimension(:,:),allocatable,save :: sinthm,costhm real*8 :: delta,delta_x,shields,ftheta,psi_x,Sbtot ! Lodewijk: for direction of sediment transport (bed slope effect) real*8 :: Ssmtot, dzbds, Sbmtot ! Lodewijk: for magnitude of sediment transport (bed slope effect) !include 's.ind' !include 's.inp' if (.not. allocated(vmag2)) then allocate(vmag2 (s%nx+1,s%ny+1)) allocate(uau (s%nx+1,s%ny+1)) allocate(uav (s%nx+1,s%ny+1)) allocate(ueu_sed (s%nx+1,s%ny+1)) allocate(uev_sed (s%nx+1,s%ny+1)) allocate(veu_sed (s%nx+1,s%ny+1)) allocate(vev_sed (s%nx+1,s%ny+1)) allocate(cu (s%nx+1,s%ny+1)) allocate(cv (s%nx+1,s%ny+1)) allocate(cc (s%nx+1,s%ny+1)) allocate(ccb (s%nx+1,s%ny+1)) allocate(fac (s%nx+1,s%ny+1,par%ngd)) allocate(Sus (s%nx+1,s%ny+1)) allocate(Svs (s%nx+1,s%ny+1)) allocate(cub (s%nx+1,s%ny+1)) allocate(cvb (s%nx+1,s%ny+1)) allocate(Sub (s%nx+1,s%ny+1)) allocate(Svb (s%nx+1,s%ny+1)) allocate(pbbedu (s%nx+1,s%ny+1)) allocate(pbbedv (s%nx+1,s%ny+1)) allocate(ccvt (s%nx+1,s%ny+1)) allocate(dcdz (s%nx+1,s%ny+1)) allocate(dsigt (s%nx+1,s%ny+1)) allocate(dsig(s%nx+1,s%ny+1,par%kmax)) allocate(ccv(s%nx+1,s%ny+1,par%kmax)) allocate(sdif(s%nx+1,s%ny+1,par%kmax)) allocate(um (s%nx+1,s%ny+1)) allocate(vm (s%nx+1,s%ny+1)) allocate(deltas(s%nx+1,s%ny+1)) allocate(sigs(s%nx+1,s%ny+1)) allocate(eswmax(s%nx+1,s%ny+1)) allocate(eswbed(s%nx+1,s%ny+1)) allocate(suq3d(s%nx+1,s%ny+1)) allocate(svq3d(s%nx+1,s%ny+1)) allocate(cuq3d(s%nx+1,s%ny+1,par%kmax)) allocate(cvq3d(s%nx+1,s%ny+1,par%kmax)) allocate(aref(s%nx+1,s%ny+1)) allocate(chain(par%kmax)) allocate(cumchain(par%kmax)) allocate (sinthm(s%nx+1,s%ny+1)) allocate (costhm(s%nx+1,s%ny+1)) delta_x = 0.d0 ! Lodewijk shields = 0.d0 ! Lodewijk ftheta = 0.d0 ! Lodewijk psi_x = 0.d0 ! Lodewijk Sbtot = 0.d0 ! Lodewijk delta = 0.d0 ! Lodewijk uau = 0.d0 uav = 0.d0 um = 0.d0 vm = 0.d0 chain = 0.0d0 cumchain = 0.0d0 fac = 1.d0 exp_ero = 0.d0 ! generate sigma grid shape... do isig=2,par%kmax chain(isig) = par%sigfac**(isig-1) cumchain(isig) = cumchain(isig-1)+chain(isig) enddo endif !numb = par%num_grid_cells !allocation ! use eulerian velocities vmag2 = s%ue**2+s%ve**2 cu = 0.0d0 cv = 0.0d0 cub = 0.0d0 cvb = 0.0d0 s%dcsdx = 0.0d0 s%dcsdy = 0.0d0 sinthm = sin(s%thetamean-s%alfaz) costhm = cos(s%thetamean-s%alfaz) ! short wave runup if (par%swrunup==1 .and. par%struct==1) then call hybrid(s,par) endif ! compute turbulence due to wave breaking if (par%lwt==1 .or. par%turb == TURB_BORE_AVERAGED .or. par%turb == TURB_WAVE_AVERAGED) then call waveturb(s,par) endif if (par%swave==1) then ! include wave skewness and assymetry in sediment advection velocity if (par%waveform==WAVEFORM_RUESSINK_VANRIJN)then call RvR(s,par) elseif (par%waveform==WAVEFORM_VANTHIEL) then call vT(s,par) endif endif ! calculate equilibrium concentration/sediment source if ((par%form==FORM_SOULSBY_VANRIJN) .or. (par%form==FORM_VANTHIEL_VANRIJN))then ! Soulsby van Rijn ! Soulsby van Rijn and Van Thiel de Vries & Reniers 2008 formulations call sedtransform(s,par) end if ! compute reduction factor for sediment sources due to presence of hard layers do jg = 1,par%ngd do j=1,s%ny+1 do i=1,s%nx+1 exp_ero = par%morfac*par%dt/(1.d0-par%por)*s%hh(i,j)*(s%ceqsg(i,j,jg)*s%pbbed(i,j,1,jg)/s%Tsg(i,j,jg) & + s%ceqbg(i,j,jg)*s%pbbed(i,j,1,jg)/par%dt) ! limit erosion to available sediment on top of hard layer wwvv changed tiny into epsilon fac(i,j,jg) = min(1.d0,s%structdepth(i,j)*s%pbbed(i,j,1,jg)/max(epsilon(0.d0),exp_ero) ) !if (fac(i,j,jg)*exp_ero > dzbed(i,j,1)*pbbed(i,j,1,jg)) then ! limit erosion to available sand in top layer ! fac(i,j,jg) = min(fac(i,j,jg),dzbed(i,j,1)*pbbed(i,j,1,jg)/max(tiny(0.d0),exp_ero) ) ! write(*,*)'WARNING: expected erosion from top layer is larger than available sand in top layer' !endif enddo enddo enddo ! compute diffusion coefficient s%Dc = par%facDc*(par%nuh+par%nuhfac*s%hh*(s%DR/par%rho)**(1.d0/3.d0)) do jg = 1,par%ngd cc = s%ccg(:,:,jg) if (s%D50(jg)>0.002d0) then ! RJ: set ceqsg to zero for gravel. ! Dano: try without this fix cc = 0.d0 ! Can be used to test total transport mode endif ! ! X-direction ! ! Get ua in u points and split out in u and v direction uau(1:s%nx,:) = 0.5*(s%ua(1:s%nx,:)*costhm(1:s%nx,:)+s%ua(2:s%nx+1,:)*costhm(1:s%nx,:)) uav(1:s%nx,:) = 0.5*(s%ua(1:s%nx,:)*sinthm(1:s%nx,:)+s%ua(2:s%nx+1,:)*sinthm(1:s%nx,:)) if (par%nz>1) then ueu_sed(1:s%nx,:) = 0.5*(s%ue_sed(1:s%nx,:)+s%ue_sed(2:s%nx+1,:)) else ueu_sed=s%ueu endif veu_sed(1:s%nx,:) = 0.5*(s%ve_sed(1:s%nx,:)+s%ve_sed(2:s%nx+1,:)) if (xmpi_isbot) then veu_sed(s%nx+1,:) = veu_sed(s%nx,:) endif ! Compute vmagu including ua ! s%vmagu = sqrt((s%uu+uau)**2+(s%vu+uav)**2) s%vmagu = sqrt((ueu_sed+uau)**2+(veu_sed+uav)**2) ! sediment advection velocity for suspended load and bed load respectively ! REMARK: when vreps does not equal vv; no mass conservation ! s%ureps = s%ueu+uau ! s%urepb = s%ueu+uau ! RJ maybe reduce this velocity? s%ureps = ueu_sed+uau s%urepb = ueu_sed+uau ! do j=1,s%ny+1 do i=1,s%nx if(s%ureps(i,j)>0.d0) then ! test cu(i,j)=cc(i,j) cu(i,j)=par%thetanum*cc(i,j)+(1.d0-par%thetanum)*cc(min(i+1,s%nx),j) cub(i,j)=par%thetanum*s%pbbed(i,j,1,jg)*s%ceqbg(i,j,jg)+(1.d0-par%thetanum)& *s%pbbed(min(i+1,s%nx),j,1,jg)*s%ceqbg(min(i+1,s%nx),j,jg) !cub(i,j)=par%thetanum*ccb(i,j)+(1.d0-par%thetanum)*ccb(min(i+1,nx),j) elseif(s%ureps(i,j)<0.d0) then cu(i,j)=par%thetanum*cc(i+1,j)+(1.d0-par%thetanum)*cc(max(i,2),j) cub(i,j)=par%thetanum*s%pbbed(i+1,j,1,jg)*s%ceqbg(i+1,j,jg)+(1.d0-par%thetanum)& *s%pbbed(max(i,2),j,1,jg)*s%ceqbg(max(i,2),j,jg) !cub(i,j)=par%thetanum*ccb(i+1,j)+(1.d0-par%thetanum)*ccb(max(i,2),j) else cu(i,j)=0.5d0*(cc(i,j)+cc(i+1,j)) cub(i,j)=0.5d0*(s%pbbed(i,j,1,jg)*s%ceqbg(i,j,jg)+s%pbbed(i+1,j,1,jg)*s%ceqbg(i+1,j,jg)) !cub(i,j)=0.5d0*(ccb(i,j)+ccb(i+1,j)) endif s%dcsdx(i,j)=(cc(i+1,j)-cc(i,j))/s%dsu(i,j) enddo enddo ! wwvv dcdx(nx:1,:) is still untouched, correct this ofr the parallel case cu(s%nx+1,:) = cc(s%nx+1,:) !Robert ! wwvv fix this in parallel case ! wwvv in parallel version, there will be a discrepancy between the values ! of dzbdx(nx+1,:). !wwvv so fix that ! Sus = 0.d0 Sub = 0.d0 ! ! suspended load, Lodewijk: no bed slope effect (yet) Sus=par%sus*(cu*s%ureps*s%hu-s%Dc*s%hu*s%dcsdx)*s%wetu ! bed load, Lodewijk: no bed slope effect (yet) Sub=par%bed*(cub*s%urepb*s%hu)*s%wetu ! ! Originally bed slope effect of XBeach : par%bdslpeffmag = 1 ! Original one, but only on bed load : par%bdslpeffmag = 2 if (par%bdslpeffmag == BDSLPEFFMAG_ROELV_TOTAL) then if (par%bermslope>0) then where (s%H/s%hu>1.0d0.or.(s%hu<1.d0.and. par%instat==INSTAT_STAT_TABLE)) Sus = Sus-par%sus*(10.d0*par%facsl*cu*s%vmagu*s%hu*(s%dzbdx-par%bermslope))*s%wetu elsewhere Sus = Sus-par%sus*(par%facsl*cu*s%vmagu*s%hu*s%dzbdx)*s%wetu end where else Sus = Sus-par%sus*(par%facsl*cu*s%vmagu*s%hu*s%dzbdx)*s%wetu endif endif if (par%bdslpeffmag == BDSLPEFFMAG_ROELV_TOTAL .or. par%bdslpeffmag == BDSLPEFFMAG_ROELV_BED) then if (par%bermslope>0) then where (s%H/s%hu>1.0d0 .or. (s%hu<1.d0.and. par%instat==INSTAT_STAT_TABLE)) Sub = Sub-par%bed*(10.0d0*par%facsl*cub*s%vmagu*s%hu*(s%dzbdx-par%bermslope))*s%wetu elsewhere Sub = Sub-par%bed*(par%facsl*cub*s%vmagu*s%hu*s%dzbdx)*s%wetu end where else Sub = Sub-par%bed*(par%facsl*cub*s%vmagu*s%hu*s%dzbdx)*s%wetu endif endif ! ! ! Y-direction ! ! Jaap: get ua in v points and split out in u and v direction if (s%ny>0) then uau(:,1:s%ny) = 0.5*(s%ua(:,1:s%ny)*costhm(:,1:s%ny)+s%ua(:,2:s%ny+1)*costhm(:,1:s%ny)) uav(:,1:s%ny) = 0.5*(s%ua(:,1:s%ny)*sinthm(:,1:s%ny)+s%ua(:,2:s%ny+1)*sinthm(:,1:s%ny)) uau(:,s%ny+1) = uau(:,s%ny) ! Jaap uav(:,s%ny+1) = uav(:,s%ny) ! Jaap if (par%nz>1) then vev_sed(:,1:s%ny) = 0.5*(s%ve_sed(:,1:s%ny)+s%ve_sed(:,2:s%ny+1)) else vev_sed=s%vev endif uev_sed(:,1:s%ny) = 0.5*(s%ue_sed(:,1:s%ny)+s%ue_sed(:,2:s%ny+1)) if (xmpi_isright) then uev_sed(:,1:s%ny+1) = uev_sed(:,1:s%ny) endif else uau=s%ua*costhm uav=s%ua*sinthm uev_sed=s%ue_sed vev_sed=s%ve_sed endif ! Jaap: compute vmagv including ua ! s%vmagv = sqrt((s%uv+uau)**2+(s%vv+uav)**2) s%vmagv = sqrt((uev_sed+uau)**2+(vev_sed+uav)**2) ! sediment advection velocity for suspended load and bed load respectively ! REMARK: when vreps does not equal vv; no mass conservation ! s%vreps = s%vev+uav ! s%vrepb = s%vev+uav ! RJ maybe reduce this velocity? Should be s%vv instead of s%vev? s%vreps = vev_sed+uav s%vrepb = vev_sed+uav ! RJ maybe reduce this velocity? Should be s%vv instead of s%vev? ! if (s%ny>0) then do j=1,s%ny do i=1,s%nx+1 if(s%vreps(i,j)>0) then cv(i,j)=par%thetanum*cc(i,j)+(1.d0-par%thetanum)*cc(i,min(j+1,s%ny)) cvb(i,j)=par%thetanum*s%pbbed(i,j,1,jg)*s%ceqbg(i,j,jg)+(1.d0-par%thetanum)& *s%pbbed(i,min(j+1,s%ny),1,jg)*s%ceqbg(i,min(j+1,s%ny),jg) !cvb(i,j)=par%thetanum*ccb(i,j)+(1.d0-par%thetanum)*ccb(i,min(j+1,ny)) elseif(s%vreps(i,j)<0) then cv(i,j)=par%thetanum*cc(i,j+1)+(1.d0-par%thetanum)*cc(i,max(j,2)) cvb(i,j)=par%thetanum*s%pbbed(i,j+1,1,jg)*s%ceqbg(i,j+1,jg)+(1.d0-par%thetanum)& *s%pbbed(i,max(j,2),1,jg)*s%ceqbg(i,max(j,2),jg) !cvb(i,j)=par%thetanum*ccb(i,j+1)+(1.d0-par%thetanum)*ccb(i,max(j,2)) else cv(i,j)=0.5d0*(cc(i,j)+cc(i,j+1)) !Jaap: cc instead of cv cvb(i,j)=0.5d0*(s%pbbed(i,j,1,jg)*s%ceqbg(i,j,jg)+s%pbbed(i,j+1,1,jg)*s%ceqbg(i,j+1,jg)) !cvb(i,j)=0.5d0*(ccb(i,j)+ccb(i,j+1)) end if s%dcsdy(i,j)=(cc(i,j+1)-cc(i,j))/s%dnv(i,j) !Jaap end do end do ! wwvv dcdy(:,ny+1) is not filled in, so in parallel case: cv(:,s%ny+1) = cc(:,s%ny+1) !Robert ! wwvv in parallel version, there will be a discrepancy between the values ! of dzbdy(:,ny+1). ! wwvv so fix that else cv = cc cvb = s%ceqbg(:,:,jg) endif ! s%ny>0 ! ! Compute sedimnent transport in v-direction ! Svs = 0.d0 Svb = 0.d0 ! Suspended load Svs=par%sus*(cv*s%vreps*s%hv-s%Dc*s%hv*s%dcsdy)*s%wetv ! Bed load Svb=par%bed*(cvb*s%vrepb*s%hv)*s%wetv ! ! Originally bed slope magnitude effect of XBeach : par%bdslpeffmag = 1 ! Original one, but only on bed load : par%bdslpeffmag = 2 if (par%bdslpeffmag == BDSLPEFFMAG_ROELV_TOTAL) then Svs = Svs-par%sus*(par%facsl*cv*s%vmagv*s%hv*s%dzbdy)*s%wetv endif if (par%bdslpeffmag == BDSLPEFFMAG_ROELV_TOTAL .or. par%bdslpeffmag == BDSLPEFFMAG_ROELV_BED ) then Svb = Svb-par%bed*(par%facsl*cvb*s%vmagv*s%hv*s%dzbdy)*s%wetv endif ! ! ! Bed slope magnitude effect (as Souslby intended) and change direction transport (see Van Rijn 1993 (section 7.2.6)) ! ! if (par%bdslpeffmag == BDSLPEFFMAG_SOULS_TOTAL .or. par%bdslpeffmag == BDSLPEFFMAG_SOULS_BED) then do j=1,s%ny+1 do i=1,s%nx+1 if ((dabs(Sub(i,j)) > 0.000001d0) .or. (dabs(Svb(i,j)) > 0.000001d0)) then Sbmtot = dsqrt( Sub(i,j)**2.d0 + Svb(i,j)**2.d0 ) dzbds = s%dzbdx(i,j)*Sub(i,j)/Sbmtot + s%dzbdy(i,j)*Svb(i,j)/Sbmtot ! dzbdn = s%dzbdx*Svb/Sbtot + s%dzbdy*Sub/Sbtot Sub(i,j) = Sub(i,j)*(1.d0 - par%facsl*dzbds) Svb(i,j) = Svb(i,j)*(1.d0 - par%facsl*dzbds) ! endif if (((dabs(Sus(i,j)) > 0.000001d0) .or. (dabs(Svs(i,j)) > 0.000001d0)) & .and. par%bdslpeffmag == BDSLPEFFMAG_SOULS_TOTAL) then Ssmtot = dsqrt( Sus(i,j)**2.d0 + Svs(i,j)**2.d0 ) dzbds = s%dzbdx(i,j)*Sus(i,j)/Ssmtot + s%dzbdy(i,j)*Svs(i,j)/Ssmtot ! dzbdn = s%dzbdx*Svb/Sbtot + s%dzbdy*Sub/Sbtot Sus(i,j) = Sus(i,j)*(1.d0 - par%facsl*dzbds) Svs(i,j) = Svs(i,j)*(1.d0 - par%facsl*dzbds) endif enddo enddo endif ! ! Lodewijk: modify the direction of the bed load transport based on the bed slope, see Van Rijn 1993 (section 7.2.6) if (par%bdslpeffdir == BDSLPEFFDIR_TALMON) then do j=1,s%ny+1 do i=1,s%nx+1 if (((dabs(s%urepb(i,j)) > 0.0001d0) .or. (dabs(s%vrepb(i,j)) > 0.0001d0)) & .and. ((dabs(s%taubx(i,j)) > 0.0001d0) .or. (dabs(s%tauby(i,j)) > 0.0001d0))) then if (s%urepb(i,j) < 0.d0) then delta_x = datan(s%vrepb(i,j)/s%urepb(i,j))+par%px ! Angle between fluid velocity vector and the s%x-axis else delta_x = datan(s%vrepb(i,j)/s%urepb(i,j)) ! Angle between fluid velocity vector and the s%x-axis endif delta = (par%rhos-par%rho)/par%rho shields = sqrt(s%taubx(i,j)**2 + s%tauby(i,j)**2)/(delta*par%rho*par%g*s%D50(jg)) ! shields = (urepb(i,j)**2.d0+vrepb(i,j)**2.d0)*s%cf(i,j)/(par%g*D50(jg)*delta) ftheta = 1.d0/(9.d0*(s%D50(jg)/s%hh(i,j))**0.3d0*sqrt(shields)) ! Talmon psi_x = datan( (dsin(delta_x)-ftheta*s%dzbdy(i,j)) / (dcos(delta_x)-ftheta*s%dzbdx(i,j)) ) psi_x = par%bdslpeffdirfac*(psi_x - delta_x)+delta_x Sbtot = dsqrt( Sub(i,j)**2.d0 + Svb(i,j)**2.d0 ) ! Magnitude of sediment transport without direction modifcation ! Decompose the sediment transport again, know with the knowledge of the direction of the sediment transport vector Sub(i,j) = Sbtot * dcos(psi_x) Svb(i,j) = Sbtot * dsin(psi_x) else Sub(i,j) = 0.d0 Svb(i,j) = 0.d0 endif enddo enddo endif ! ! ! do j=1,s%ny+1 do i=1,s%nx if(Sub(i,j)>0.d0) then pbbedu(i,j) = s%pbbed(i,j,1,jg) elseif(Sub(i,j)<0.d0) then pbbedu(i,j)= s%pbbed(i+1,j,1,jg) else pbbedu(i,j)=0.5d0*(s%pbbed(i,j,1,jg)+s%pbbed(i+1,j,1,jg)) endif enddo enddo ! Sub = pbbedu*Sub ! do j=1,s%ny do i=1,s%nx+1 if(Svb(i,j)>0) then pbbedv(i,j)=s%pbbed(i,j,1,jg) else if(Svb(i,j)<0) then pbbedv(i,j)=s%pbbed(i,j+1,1,jg) else pbbedv(i,j)=0.5d0*(s%pbbed(i,j,1,jg)+s%pbbed(i,j+1,1,jg)) end if end do end do ! Svb = pbbedv*Svb ! ! BRJ: implicit concentration update (compute sources first, sink must be computed after updating actual sed.conc.) ! if (s%ny>0) then do j=2,s%ny do i=2,s%nx ! Changed to hh from hold by RJ (13072009) !**2/max(hh(i,j),par%hmin) s%ero(i,j,jg) = fac(i,j,jg)*s%hh(i,j)*s%ceqsg(i,j,jg)*s%pbbed(i,j,1,jg)/s%Tsg(i,j,jg) ! depo_ex(i,j,jg) = max(hold(i,j),0.01d0)*cc(i,j)/Tsg(i,j,jg) ! BRJ: the volume in the water column is updated and not the volume concentration. cc(i,j) = (par%dt*s%Tsg(i,j,jg))/(par%dt+s%Tsg(i,j,jg))* & (s%hold(i,j)*cc(i,j)/par%dt -((Sus(i,j)*s%dnu(i,j)-Sus(i-1,j)*s%dnu(i-1,j)+& Svs(i,j)*s%dsv(i,j)-Svs(i,j-1)*s%dsv(i,j-1))*s%dsdnzi(i,j)-& s%ero(i,j,jg))) cc(i,j)=max(cc(i,j),0.0d0) ! Jaap: negative cc's are possible... cc(i,j)=min(cc(i,j),par%cmax*s%hh(i,j)) s%depo_ex(i,j,jg) = cc(i,j)/s%Tsg(i,j,jg) enddo enddo else j=1 do i=2,s%nx ! Changed to hh from hold by RJ (13072009) !**2/max(hh(i,j),par%hmin) s%ero(i,j,jg) = fac(i,j,jg)*s%hh(i,j)*s%ceqsg(i,j,jg)*s%pbbed(i,j,1,jg)/s%Tsg(i,j,jg) ! depo_ex(i,j,jg) = max(hold(i,j),0.01d0)*cc(i,j)/Tsg(i,j,jg) ! BRJ: the volume in the water column is updated and not the volume concentration. cc(i,j) = (par%dt*s%Tsg(i,j,jg))/(par%dt+s%Tsg(i,j,jg))* & (s%hold(i,j)*cc(i,j)/par%dt -((Sus(i,j)*s%dnu(i,j)-Sus(i-1,j)*s%dnu(i-1,j))*s%dsdnzi(i,j)-& s%ero(i,j,jg))) cc(i,j)=max(cc(i,j),0.0d0) ! Jaap: negative cc's are possible... cc(i,j)=min(cc(i,j),par%cmax*s%hh(i,j)) s%depo_ex(i,j,jg) = cc(i,j)/s%Tsg(i,j,jg) enddo endif cc = cc/s%hh ! do lateral boundaries... if(xmpi_istop)then cc(1,:)=cc(2,:) s%ero(1,:,jg)=s%ero(2,:,jg) s%depo_ex(1,:,jg)=s%depo_ex(2,:,jg) endif if(xmpi_isleft .and. s%ny>0)then cc(:,1)=cc(:,2) s%ero(:,1,jg)=s%ero(:,2,jg) s%depo_ex(:,1,jg)=s%depo_ex(:,2,jg) endif if(xmpi_istop)then cc(s%nx+1,:)=cc(s%nx+1-1,:) s%ero(s%nx+1,:,jg)=s%ero(s%nx,:,jg) s%depo_ex(s%nx+1,:,jg)=s%depo_ex(s%nx,:,jg) endif if(xmpi_isright .and. s%ny>0)then cc(:,s%ny+1)=cc(:,s%ny+1-1) s%ero(:,s%ny+1,jg)=s%ero(:,s%ny,jg) s%depo_ex(:,s%ny+1,jg)=s%depo_ex(:,s%ny,jg) endif ! wwvv fix the first and last rows and columns of cc in parallel case #ifdef USEMPI call xmpi_shift_ee(cc) #endif ! ! wwvv border columns and rows of ccg Svg and Sug have to be communicated s%ccg(:,:,jg) = cc s%Svsg(:,:,jg) = Svs s%Susg(:,:,jg) = Sus s%Svbg(:,:,jg) = Svb s%Subg(:,:,jg) = Sub end do ! number of sediment fractions s%vmag=sqrt(max(vmag2,par%umin)) end subroutine transus !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine bed_update(s,par,zs0data) use params use interp use spaceparams use xmpi_module implicit none type(spacepars),target :: s type(parameters) :: par integer :: i,j,j1,jg,ii,ie,id,je,jd,jdz,ndz, hinterland,numb real*8 :: zs0data integer , dimension(:,:,:),allocatable,save :: indSus,indSub,indSvs,indSvb real*8 :: dzb,dzmax,dzt,dzleft,sdz,dzavt,fac,Savailable,dAfac real*8 , dimension(:,:),allocatable,save :: Sout,hav,zb_initial2,struct_initial2,tracker_vert_change,hsum,Hwavesum,taubxsum,zbsum,vmagsum,glmsum,eulsum,zssum,structdepth_in,structdepth_adj,structdepth_out real*8 , dimension(par%ngd) :: edg,edg1,edg2,dzg !real*8, dimension(365,1),save :: real*8 , dimension(:),pointer :: dz real*8 , dimension(:,:),pointer :: pb integer :: n_aval real*8 :: morf_adj !include 's.ind' !include 's.inp' if (.not. allocated(Sout)) then allocate(Sout(s%nx+1,s%ny+1)) allocate(hav(s%nx+1,s%ny+1)) allocate(indSus(s%nx+1,s%ny+1,par%ngd)) allocate(indSub(s%nx+1,s%ny+1,par%ngd)) allocate(indSvs(s%nx+1,s%ny+1,par%ngd)) allocate(indSvb(s%nx+1,s%ny+1,par%ngd)) endif ! Super fast 1D if (s%ny==0) then j1 = 1 else j1 = 2 endif s%dzbnow = 0.d0 dzb = 0.d0 if (par%t>=par%morstart .and. par%t < par%morstop .and. par%morfac > .999d0) then ! ! bed_predict ! ! reduce sediment transports when hard layer comes to surface ! this step is mainly necessary at the transition from hard layers to sand if (par%struct == 1) then do jg = 1,par%ngd indSus = 0 indSub = 0 indSvs = 0 indSvb = 0 Sout = 0.d0 do j=j1,s%ny+1 do i=2,s%nx+1 ! fluxes at i,j if (s%Subg(i,j,jg) > 0.d0) then ! bed load s%u-direction indSub(i,j,jg) = 1 Sout(i,j) = Sout(i,j) + s%Subg(i,j,jg)*s%dnu(i,j) endif ! fluxes at i-1,j if (s%Subg(i-1,j,jg) < 0.d0 ) then ! bed load s%u-direction Sout(i,j) = Sout(i,j) - s%Subg(i-1,j,jg)*s%dnu(i-1,j) endif if (par%sourcesink==0) then ! fluxes at i,j if (s%Susg(i,j,jg) > 0.d0 ) then ! suspended load s%u-direction indSus(i,j,jg) = 1 Sout(i,j) = Sout(i,j) + s%Susg(i,j,jg)*s%dnu(i,j) endif ! fluxes at i-1,j if (s%Susg(i-1,j,jg) < 0.d0 ) then ! suspended load s%u-direction Sout(i,j) = Sout(i,j) - s%Susg(i-1,j,jg)*s%dnu(i-1,j) endif endif enddo enddo if (s%ny>0) then do j=j1,s%ny+1 do i=2,s%nx+1 if (s%Svbg(i,j,jg) > 0.d0 ) then ! bed load s%v-direction indSvb(i,j,jg) = 1 Sout(i,j) = Sout(i,j) + s%Svbg(i,j,jg)*s%dsv(i,j) endif ! fluxes at i,j-1 if (s%Svbg(i,j-1,jg) < 0.d0 ) then ! bed load s%v-direction Sout(i,j) = Sout(i,j) - s%Svbg(i,j-1,jg)*s%dsv(i,j-1) endif if (par%sourcesink==0) then if (s%Svsg(i,j,jg) > 0.d0 ) then ! suspended load s%v-direction indSvs(i,j,jg) = 1 Sout(i,j) = Sout(i,j) + s%Svsg(i,j,jg)*s%dsv(i,j) endif ! fluxes at i,j-1 if (s%Svsg(i,j-1,jg) < 0.d0 ) then ! suspended load s%v-direction Sout(i,j) = Sout(i,j) - s%Svsg(i,j-1,jg)*s%dsv(i,j-1) endif endif ! sourcesink = 0 enddo !s%nx+1 enddo !s%ny+1 endif !s%ny>0 ! do j=j1,s%ny+1 do i=2,s%nx+1 Savailable = s%structdepth(i,j)*s%pbbed(i,j,1,jg)/par%morfac/par%dt*(1.d0-par%por)/s%dsdnzi(i,j) ! reduction factor for cell outgoing sediment transports wwvv changed tiny into epsilon fac = min(1.d0,Savailable/max(Sout(i,j),epsilon(0.d0)) ) ! fix sediment transports for the presence of a hard layer; remind indSus etc are 1 in cases of cell outgoing transports ! updated S oell outgoing transports cell incoming transports if (fac < 1.d0)then s%Subg(i,j,jg) = fac*indSub(i,j,jg)*s%Subg(i,j,jg) + (1-indSub(i,j,jg))*s%Subg(i,j,jg) s%Subg(i-1,j,jg) = fac*(1-indSub(i-1,j,jg))*s%Subg(i-1,j,jg) + indSub(i-1,j,jg)*s%Subg(i-1,j,jg) if (s%ny>0) then s%Svbg(i,j,jg) = fac*indSvb(i,j,jg)*s%Svbg(i,j,jg) + (1-indSvb(i,j,jg))*s%Svbg(i,j,jg) s%Svbg(i,j-1,jg) = fac*(1-indSvb(i,j-1,jg))*s%Svbg(i,j-1,jg) + indSvb(i,j-1,jg)*s%Svbg(i,j-1,jg) endif if (par%sourcesink==0) then s%Susg(i,j,jg) = fac*indSus(i,j,jg)*s%Susg(i,j,jg) + (1-indSus(i,j,jg))*s%Susg(i,j,jg) s%Susg(i-1,j,jg) = fac*(1-indSus(i-1,j,jg))*s%Susg(i-1,j,jg) + indSus(i-1,j,jg)*s%Susg(i-1,j,jg) if (s%ny>0) then s%Svsg(i,j,jg) = fac*indSvs(i,j,jg)*s%Svsg(i,j,jg) + (1-indSvs(i,j,jg))*s%Svsg(i,j,jg) s%Svsg(i,j-1,jg) = fac*(1-indSvs(i,j-1,jg))*s%Svsg(i,j-1,jg) + indSvs(i,j-1,jg)*s%Svsg(i,j-1,jg) endif !s%ny = 0 endif ! sourcesink = 0 endif !fac<1.d0 enddo ! s%nx+1 enddo !s%ny + 1 enddo !par%ngd endif !struct == 1 if (s%ny>0) then do j=2,s%ny do i=2,s%nx ! bed level changes per fraction in this morphological time step in meters sand including pores ! positive in case of erosion if (par%sourcesink==0) then dzg=par%morfac*par%dt/(1.d0-par%por)*( & ! Dano, dz from sus transport gradients ( s%Susg(i,j,:)*s%dnu(i,j)-s%Susg(i-1,j,:)*s%dnu(i-1,j) +& s%Svsg(i,j,:)*s%dsv(i,j)-s%Svsg(i,j-1,:)*s%dsv(i,j-1) +& ! dz from bed load transport gradients s%Subg(i,j,:)*s%dnu(i,j)-s%Subg(i-1,j,:)*s%dnu(i-1,j)+& s%Svbg(i,j,:)*s%dsv(i,j)-s%Svbg(i,j-1,:)*s%dsv(i,j-1) )*s%dsdnzi(i,j) ) elseif (par%sourcesink==1) then dzg=par%morfac*par%dt/(1.d0-par%por)*( & s%ero(i,j,:)-s%depo_ex(i,j,:) +& ( s%Subg(i,j,:)*s%dnu(i,j)-s%Subg(i-1,j,:)*s%dnu(i-1,j)+& s%Svbg(i,j,:)*s%dsv(i,j)-s%Svbg(i,j-1,:)*s%dsv(i,j-1) )*s%dsdnzi(i,j) ) endif if (par%ngd==1) then ! Simple bed update in case one fraction s%zb(i,j) = s%zb(i,j)-sum(dzg) s%dzbnow(i,j) = s%dzbnow(i,j)-sum(dzg) ! naamgeveing? s%dzbdt(i,j) = s%dzbnow(i,j)/par%dt s%sedero(i,j) = s%sedero(i,j)-sum(dzg) s%structdepth(i,j) = max(0.d0,s%structdepth(i,j)-sum(dzg)) else ! erosion/deposition rate of sand mass (m/s) ! positive in case of erosion edg = dzg*(1.d0-par%por)/par%dt dz=>s%dzbed(i,j,:) pb=>s%pbbed(i,j,:,:) call update_fractions(par,s,i,j,dz,pb,edg,sum(dzg)) endif enddo ! s%nx+1 enddo ! s%ny+1 else j=1 do i=2,s%nx ! bed level changes per fraction in this morphological time step in meters sand including pores ! positive in case of erosion if (par%sourcesink==0) then dzg=par%morfac*par%dt/(1.d0-par%por)*( & ! Dano, dz from sus transport gradients ( s%Susg(i,j,:)*s%dnu(i,j)-s%Susg(i-1,j,:)*s%dnu(i-1,j) +& s%Subg(i,j,:)*s%dnu(i,j)-s%Subg(i-1,j,:)*s%dnu(i-1,j) )*s%dsdnzi(i,j) +& (s%Svsg(i,j,:)+s%Svbg(i,j,:))*par%lsgrad) elseif (par%sourcesink==1) then dzg=par%morfac*par%dt/(1.d0-par%por)*( & s%ero(i,j,:)-s%depo_ex(i,j,:) +& ( s%Subg(i,j,:)*s%dnu(i,j)-s%Subg(i-1,j,:)*s%dnu(i-1,j) )*s%dsdnzi(i,j) +& (s%Svsg(i,j,:)+s%Svbg(i,j,:))*par%lsgrad) endif if (par%ngd==1) then ! Simple bed update in case one fraction s%zb(i,j) = s%zb(i,j)-sum(dzg) s%dzbnow(i,j) = s%dzbnow(i,j)-sum(dzg) s%dzbdt(i,j) = s%dzbnow(i,j)/par%dt s%sedero(i,j) = s%sedero(i,j)-sum(dzg) s%structdepth(i,j) = max(0.d0,s%structdepth(i,j)-sum(dzg)) else ! multiple fractions... ! erosion/deposition rate of sand mass (m/s) ! positive in case of erosion edg = dzg*(1.d0-par%por)/par%dt dz=>s%dzbed(i,j,:) pb=>s%pbbed(i,j,:,:) call update_fractions(par,s,i,j,dz,pb,edg,sum(dzg)) endif enddo ! s%nx+1 endif !s%ny>0 #ifdef USEMPI call xmpi_shift_ee(s%zb) #endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Avalanching ! ! ! if (par%avalanching==1) then do ii=1,nint(par%morfac) !aval=.false. n_aval=0 !Dano fix for MPI s%dzbdx=0.d0 s%dzbdy=0.d0 do j=1,s%ny+1 do i=1,s%nx if (s%structdepth(i,j) >= par%ALimit) then s%dzbdx(i,j)=(s%zb(i+1,j)-s%zb(i,j))/s%dsu(i,j) endif enddo enddo ! hav = s%hh ! Fix Hav for short wave runup: if (par%swrunup == 1) then do j = 1,s%ny+1 hinterland = 0 do i = 1, s%nx+1 if (s%structdepth(i,j) >= par%ALimit) then if (hinterland == 0 .and. s%runup(j)+s%zs(nint(s%iwl(j)),j) < s%zb(i,j) + par%eps) then hinterland = 1; endif if (s%wetz(i,j) == 1) then hav(i,j) = max(par%eps,(s%hh(i,j) + s%runup(j))) elseif (hinterland == 0) then hav(i,j) = max(par%eps,s%runup(j)+s%zs(nint(s%iwl(j)),j)-s%zb(i,j) ) else hav(i,j) = par%eps; endif endif enddo enddo endif ! do i=2,s%nx !-1 Jaap -1 gives issues for bed updating at mpi boundaries do j=1,s%ny+1 if (s%structdepth(i,j) >= par%ALimit) then !if (max( max(hh(i,j),par%delta*H(i,j)), max(hh(i+1,j),par%delta*H(i+1,j)) )>par%hswitch+par%eps) then if(max(hav(i,j),hav(i+1,j))>par%hswitch+par%eps) then ! Jaap instead of s%hh dzmax=par%wetslp; ! tricks: seaward of istruct (transition from sand to structure) wetslope is set to 0.03; if (i>nint(s%istruct(j))) then !dzmax = 0.03d0 dzmax = max(0.03d0,abs(s%dzbdx(i,j))*0.99d0) endif else dzmax=par%dryslp; end if if(abs(s%dzbdx(i,j))>dzmax+1.d-10 .and. s%structdepth(i+nint(max(0.d0,sign(1.d0,s%dzbdx(i,j)))),j)>par%eps) then n_aval=n_aval+1 dzb=sign(1.0d0,s%dzbdx(i,j))*(abs(s%dzbdx(i,j))-dzmax)*s%dsu(i,j) if (dzb >= 0.d0) then ie = i+1 ! index erosion point id = i ! index deposition point dAfac = s%dsdnzi(i,j)/s%dsdnzi(i+1,j) ! take into account varying grid sizes dzb=min(dzb,par%dzmax*par%dt/s%dsu(i,j)) ! make sure dzb is not in conflict with par%dzmax dzb=min(dzb,s%structdepth(i+1,j)) ! make sure dzb is not larger than sediment layer thickness else ie = i ! index erosion point id = i+1 ! index deposition point dAfac = s%dsdnzi(i+1,j)/s%dsdnzi(i,j) ! take into account varying grid sizes dzb=max(dzb,-par%dzmax*par%dt/s%dsu(i,j)) dzb=max(dzb,-s%structdepth(i,j)) endif if (par%ngd == 1) then ! Simple bed update in case one fraction dzleft = abs(dzb) s%zb(id,j) = s%zb(id,j)+dzleft*dAfac s%zb(ie,j) = s%zb(ie,j)-dzleft s%dzbnow(id,j) = s%dzbnow(id,j)+dzleft*dAfac ! naamgeveing? s%dzbnow(ie,j) = s%dzbnow(ie,j)-dzleft s%sedero(id,j) = s%sedero(id,j)+dzleft*dAfac s%sedero(ie,j) = s%sedero(ie,j)-dzleft s%structdepth(id,j) = max(0.d0,s%structdepth(id,j)+dzleft*dAfac) s%structdepth(ie,j) = max(0.d0,s%structdepth(ie,j)-dzleft) s%zs(id,j) = s%zs(id,j)+dzleft*dAfac s%zs(ie,j) = s%zs(ie,j)-dzleft s%dzav(id,j)= s%dzav(id,j)+dzleft*dAfac s%dzav(ie,j)= s%dzav(ie,j)-dzleft else ! multiple fractions... ! now fix fractions.... dz => s%dzbed(ie,j,:) pb => s%pbbed(ie,j,:,:) ! figure out how many depth layers (ndz) are eroded in point iii sdz = 0 ndz = 0 do while (sdz 0 ) edg1(jg) = -s%sedcal(jg)*edg2(jg)*dAfac ! deposition (dzt always < 0 ) enddo dzavt = dzavt + sum(edg2)*par%dt/(1.d0-par%por) call update_fractions(par,s,ie,j,s%dzbed(ie,j,:),s%pbbed(ie,j,:,:),edg2,dzavt) ! update bed in eroding point call update_fractions(par,s,id,j,s%dzbed(id,j,:),s%pbbed(id,j,:,:),edg1,-dzavt*dAfac) ! update bed in deposition point enddo ! update water levels and dzav s%zs(ie,j) = s%zs(ie,j)-dzavt s%dzav(ie,j)= s%dzav(ie,j)-dzavt s%zs(id,j) = s%zs(id,j)+dzavt*dAfac s%dzav(id,j)= s%dzav(id,j)+dzavt*dAfac end if ! yes/no multiple fractions end if endif end do end do !JJ: update y slopes after avalanching in X-direction seems more appropriate do j=1,s%ny do i=1,s%nx+1 if (s%structdepth(i,j) >= par%ALimit) then s%dzbdy(i,j)=(s%zb(i,j+1)-s%zb(i,j))/s%dnv(i,j) endif enddo enddo do j=2,s%ny !-1 Jaap -1 gives issues for bed updating at mpi boundaries do i=1,s%nx+1 if (s%structdepth(i,j) >= par%ALimit) then if(max(s%hh(i,j),s%hh(i,j+1))>par%hswitch+par%eps) then dzmax=par%wetslp else dzmax=par%dryslp end if if(abs(s%dzbdy(i,j))>dzmax .and. s%structdepth(i,j+nint(max(0.d0,sign(1.d0,s%dzbdy(i,j)))))>par%eps) then ! Jaap n_aval=n_aval+1 dzb=sign(1.0d0,s%dzbdy(i,j))*(abs(s%dzbdy(i,j))-dzmax)*s%dnv(i,j) ! if (dzb >= 0.d0) then je = j+1 ! index erosion point jd = j ! index deposition point dAfac = s%dsdnzi(i,j)/s%dsdnzi(i,j+1) ! take into account varying grid sizes dzb=min(dzb,par%dzmax*par%dt/s%dnv(i,j)) dzb=min(dzb,s%structdepth(i,j+1)) else je = j ! index erosion point jd = j+1 ! index deposition point dAfac = s%dsdnzi(i,j+1)/s%dsdnzi(i,j) ! take into account varying grid sizes dzb=max(dzb,-par%dzmax*par%dt/s%dnv(i,j)) dzb=max(dzb,-s%structdepth(i,j)) endif if (par%ngd == 1) then ! Simple bed update in case one fraction dzleft = abs(dzb) s%zb(i,jd) = s%zb(i,jd)+dzleft*dAfac s%zb(i,je) = s%zb(i,je)-dzleft s%dzbnow(i,jd) = s%dzbnow(i,jd)+dzleft*dAfac ! naamgeveing? s%dzbnow(i,je) = s%dzbnow(i,je)-dzleft s%sedero(i,jd) = s%sedero(i,jd)+dzleft*dAfac s%sedero(i,je) = s%sedero(i,je)-dzleft s%structdepth(i,jd) = max(0.d0,s%structdepth(i,jd)+dzleft*dAfac) s%structdepth(i,je) = max(0.d0,s%structdepth(i,je)-dzleft) s%zs(i,jd) = s%zs(i,jd)+dzleft*dAfac s%zs(i,je) = s%zs(i,je)-dzleft s%dzav(i,jd)= s%dzav(i,jd)+dzleft*dAfac s%dzav(i,je)= s%dzav(i,je)-dzleft else ! multiple fractions... dz => s%dzbed(i,je,:) pb => s%pbbed(i,je,:,:) ! figure out how many depth layers (ndz) are affected sdz = 0 ndz = 0 do while (sdz 0 ) edg1(jg) = -s%sedcal(jg)*edg2(jg)*dAfac ! deposition (dzt always < 0 ) enddo dzavt = dzavt + sum(edg2)*par%dt/(1.d0-par%por) call update_fractions(par,s,i,je,s%dzbed(i,je,:),s%pbbed(i,je,:,:),edg2,dzavt) ! upwind point call update_fractions(par,s,i,jd,s%dzbed(i,jd,:),s%pbbed(i,jd,:,:),edg1,-dzavt*dAfac) ! downwind point enddo ! update water levels and dzav s%zs(i,je) = s%zs(i,je)-dzavt s%dzav(i,je)= s%dzav(i,je)-dzavt s%zs(i,jd) = s%zs(i,jd)+dzavt*dAfac s%dzav(i,jd)= s%dzav(i,jd)+dzavt*dAfac endif !yes/no multiple fractions end if endif end do end do ! write(*,*)'t ',par%t,'ii ',ii,' n_aval ',n_aval,' rank ',xmpi_rank ! Dano: in parallel version bed update must take place AFTER EACH ITERATION #ifdef USEMPI call xmpi_allreduce(n_aval ,MPI_SUM) #endif if (n_aval>0) then #ifdef USEMPI call xmpi_shift_ee(s%zb) #endif else exit endif end do else s%dzbdx=0.d0 s%dzbdy=0.d0 do j=1,s%ny+1 do i=1,s%nx s%dzbdx(i,j)=(s%zb(i+1,j)-s%zb(i,j))/s%dsu(i,j) enddo enddo do j=1,s%ny do i=1,s%nx+1 s%dzbdy(i,j)=(s%zb(i,j+1)-s%zb(i,j))/s%dnv(i,j) enddo enddo end if ! ! bed boundary conditions ! ! Dano: the following Neumann boundaries introduce a mass error equal to ! the bed level change due to the b.c. times the cell area. if(xmpi_isleft .and. s%ny>0) then s%zb(:,1) = s%zb(:,2) s%dzbdt(:,1) = s%dzbdt(:,2) s%dzbnow(:,1) = s%dzbnow(:,2) s%sedero(:,1) = s%sedero(:,2) s%structdepth(:,1) = s%structdepth(:,2) s%pbbed(:,1,:,:)=s%pbbed(:,2,:,:) s%z0bed(:,1)=s%z0bed(:,2) s%dzbed(:,1,:)=s%dzbed(:,2,:) endif if(xmpi_isright .and. s%ny>0) then s%zb(:,s%ny+1) = s%zb(:,s%ny) s%dzbdt(:,s%ny+1) = s%dzbdt(:,s%ny) s%dzbnow(:,s%ny+1) = s%dzbnow(:,s%ny) s%sedero(:,s%ny+1) = s%sedero(:,s%ny) s%structdepth(:,s%ny+1) = s%structdepth(:,s%ny) s%pbbed(:,s%ny+1,:,:)=s%pbbed(:,s%ny,:,:) s%z0bed(:,s%ny+1)=s%z0bed(:,s%ny) s%dzbed(:,s%ny+1,:)=s%dzbed(:,s%ny,:) endif ! Robert: in parallel version bed update must take place on internal boundaries: #ifdef USEMPI call xmpi_shift_ee(s%zb) #endif ! Update representative sed.diameter at the bed for flow friction and output if (par%ngd>1) then do j=j1,max(s%ny,1) do i=2,s%nx s%D50top = sum(s%pbbed(i,j,1,:)*s%D50) s%D90top = sum(s%pbbed(i,j,1,:)*s%D90) enddo enddo endif endif ! if par%t>par%morstart ! !Add-on for Thermal Routine Implementation - Updated 07/12/2017 ! if (par%struct == 1) then if (par%counter_steps < 1) then allocate (zb_initial2(s%nx+1,s%ny+1)) !Bed location at first timestep allocate (struct_initial2(s%nx+1,s%ny+1)) !Structdepth at first timestep allocate (hsum(s%nx+1,s%ny+1)) allocate (zbsum(s%nx+1,s%ny+1)) allocate (zssum(s%nx+1,s%ny+1)) allocate (vmagsum(s%nx+1,s%ny+1)) allocate (glmsum(s%nx+1,s%ny+1)) allocate (eulsum(s%nx+1,s%ny+1)) allocate (Hwavesum(s%nx+1,s%ny+1)) allocate (taubxsum(s%nx+1,s%ny+1)) allocate (structdepth_adj(s%nx+1,s%ny+1)) allocate (structdepth_in(s%nx+1,s%ny+1)) allocate (structdepth_out(s%nx+1,s%ny+1)) call csv_zs0read(zs0data) ! Mean water level for each 30 minute period do j=1,s%ny+1 do i=1,s%nx+1 zbsum(i,j) = 0.0D+00 zssum(i,j) = 0.0D+00 vmagsum(i,j) = 0.0D+00 glmsum(i,j) = 0.0D+00 eulsum(i,j) = 0.0D+00 hsum(i,j) = 0.0D+00 Hwavesum(i,j) = 0.0D+00 taubxsum(i,j) = 0.0D+00 structdepth_adj(i,j) = 0.0D+00 structdepth_in(i,j) = 0.0D+00 structdepth_in(i,j) = s%structdepth(i,j) !State if structdepth at the beginning of the thermal routine end do end do do j=1,s%ny+1 do i=1,s%nx+1 zb_initial2(i,j) = 0.0D+00 struct_initial2(i,j) = 0.0D+00 open(unit=203,FILE="bed.dep", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="asis", ACTION="read") open(unit=204,FILE="permadepth.dep", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="asis", ACTION="read") read(203,300) zb_initial2(i,j) read(204,300) struct_initial2(i,j) 300 format (f10.5) end do end do endif ! do j=1,s%ny+1 do i=1,s%nx+1 zbsum(i,j) = zbsum(i,j) + s%zb(i,j) zssum(i,j) = zssum(i,j) + abs(s%zs(i,j)-zs0data) vmagsum(i,j) = vmagsum(i,j) + s%vmag(i,j) glmsum(i,j) = glmsum(i,j) + abs(s%u(i,j)) eulsum(i,j) = eulsum(i,j) + abs(s%ue(i,j)) hsum(i,j) = hsum(i,j)+s%hh(i,j) Hwavesum(i,j) = Hwavesum(i,j)+abs(s%H(i,j)) taubxsum(i,j) = taubxsum(i,j)+abs(s%taubx(i,j)) end do end do ! morf_adj=0.0D+00 if(par%morfacopt==1) then ! Setting the morfac parameter correctly based on morfacopt setting !!!!!!! morf_adj=1.0D+00 else if(par%morfacopt==0) then morf_adj=1.0D+00 end if par%counter_ub = par%counter_ub + 1 par%counter_steps = par%counter_steps + 1 if (par%timecalc >= (par%therm_p/morf_adj)) then ! if (par%timet < (par%therm_p/morf_adj)+5) then allocate(tracker_vert_change(s%nx+1,s%ny+1)) ! Arrayed used to keep track of change in thaw during storm event do j=1,s%ny+1 do i=1,s%nx+1 tracker_vert_change(i,j)=0.0D+00 end do end do end if ! if (par%sediment_removal == 1) then if(par%timet >= (par%tstop - (par%therm_p/morf_adj))) then !if (par%timet > (par%therm_p/morf_adj)+5) then do j=1,s%ny+1 do i=1,s%nx+1 structdepth_adj(i,j) = structdepth_in(i,j) - struct_initial2(i,j) - tracker_vert_change(i,j) !structdepth_adj(i,j) = structdepth_in(i,j) - structdepth_out(i,j) if (structdepth_adj(i,j) > 0) then !Deposition (structdepth_in(i,j) - structdepth_out(i,j) > 0) if (s%zb(i,j) > zb_initial2(i,j)) then s%structdepth(i,j) = struct_initial2(i,j) + (structdepth_adj(i,j)*(1-par%sed_removal)) !s%structdepth(i,j) = structdepth_out(i,j) + (structdepth_adj(i,j)*(1-par%nb)) s%zb(i,j) = zb_initial2(i,j) + (structdepth_adj(i,j)*(1-par%sed_removal)) endif end if enddo enddo endif endif call heat_transfer(s,par,zbsum,hsum,Hwavesum,taubxsum,vmagsum,glmsum,eulsum,zssum,structdepth_out,tracker_vert_change) par%timecalc = par%timecalc*0 par%counter_ub = par%counter_ub*0 do j=1,s%ny+1 do i=1,s%nx+1 hsum(i,j) = hsum(i,j)*0 zbsum(i,j) = zbsum(i,j)*0 zssum(i,j) = zssum(i,j)*0 vmagsum(i,j) = vmagsum(i,j)*0 glmsum(i,j) = glmsum(i,j)*0 eulsum(i,j) = eulsum(i,j)*0 Hwavesum(i,j) = Hwavesum(i,j)*0 taubxsum(i,j) = taubxsum(i,j)*0 enddo enddo call write_test(s,par,tracker_vert_change) !deallocate(tracker_vert_change) !deallocate(hsum) !deallocate(zbsum) !deallocate(zssum) !deallocate(vmagsum) !deallocate(glmsum) !deallocate(eulsum) !deallocate(Hwavesum) !deallocate(taubxsum) endif par%timecalc=par%timecalc+par%dt ! Logger for thawdepth. Used to determine the time period used for the thermal routine calculations par%timet=par%timet+par%dt endif ! ! ! end subroutine bed_update !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine update_fractions(par,s,i,j,dz,pb,edg,dzb) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Copyright (C) 2009 Technische Universiteit Delft ! Bram van Prooijen ! b.c.vanprooijen@tudelft.nl ! +31(0)15 2784070 ! Faculty of Civil Engineering and Geosciences ! department of Hydraulic Engineering ! PO Box 5048 ! 2600 GA Delft ! The Netherlands !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use params use spaceparams use xmpi_module implicit none type(spacepars),target :: s type(parameters) :: par integer :: i,j,jg,jd real*8 :: ED,zbold,dzbt,fac,dzb_loc real*8, intent(in) :: dzb real*8 , dimension(par%ngd),intent(in) :: edg real*8 , dimension(s%nd(i,j)),intent(inout) :: dz real*8 , dimension(s%nd(i,j),par%ngd),intent(inout) :: pb real*8 , dimension(:),allocatable,save :: Ap,b real*8 , dimension(:,:),allocatable,save :: Sm,A if (.not. allocated(Ap)) then allocate(Ap (s%nd(1,1))) allocate(b (s%nd(1,1))) allocate(Sm (s%nd(1,1),par%ngd)) allocate(A (s%nd(1,1),3)) endif !TODO, dzb_loc is not initialized can be Nan, leading to infinite loop dzb_loc = dzb !!!initialize Sm ED = sum(edg) ! do t_sub=1,nt_sub !loop over subtimesteps do while (abs(dzb_loc) .gt. 0.d0) dzbt = min(dzb_loc,dz(par%nd_var)) ! make sure erosion (dzg is positive) is limited to thickness of variable layer dzbt = max(dzbt,-par%frac_dz*dz(par%nd_var+1)) ! make sure deposition (dzg is negative) is limited to thickness of first layer below variable layer fac = dzbt/dzb_loc ! factor over mass change in cells to limit erosion and deposition dzb_loc = dzb_loc-dzbt ! update dzb do jg=1,par%ngd A=0. Ap=0. b=0. Sm(:,jg)=pb(:,jg)*dz*(1.d0-par%por) !!!build matrix A select case(par%nd_var) case(1) !in this case: A=0 case(2) A (1,1:3)= (/0.d0 , min(ED,0.d0) , max(ED,0.d0) /) A (2,1:3)= (/-min(ED,0.d0) , -max(ED,0.d0) , 0.d0 /) case(3:1000) A (1,1:3)= (/0.d0 , min(ED,0.d0) , max(ED,0.d0) /) A (2:par%nd_var-1,1)=-min(ED,0.d0) A (2:par%nd_var-1,2)=-abs(ED) A (2:par%nd_var-1,3)=max(ED,0.d0) A (par%nd_var,1:3)= (/-min(ED,0.d0) , -max(ED,0.d0) , 0.d0 /) end select !!!determine RHS ! Ap makes sure that layer nd_var varies in thickness in time ! Ap = 0 with single fraction (in case(1)) Ap(1) = sum(A(1,2:3)*pb(1:2,jg)) do jd = 2,par%nd_var Ap(jd) = sum(A (jd,:)*pb(jd-1:jd+1,jg)) enddo Ap(par%nd_var+1) = sum(A(par%nd_var+1,1:2)*pb(par%nd_var:par%nd_var+1,jg)) ! b represents the actual erosion and deposition in the top layer. ! However, the thickness of the top layer remains constant and instead layer nd_var will breath in thickness b(1) = -edg(jg) !!!update Sm ! Sm is the sediment mass per fraction per layer ! Sm(1:par%nd_var+1,jg) = Sm(1:par%nd_var+1,jg) + par%dt/nt_sub*(Ap(1:par%nd_var+1)+b(1:par%nd_var+1)) Sm(1:par%nd_var+1,jg) = Sm(1:par%nd_var+1,jg) + par%dt*fac*(Ap(1:par%nd_var+1)+b(1:par%nd_var+1)) enddo !fractions ! From Sm we can compute the new fraction ratios per layer and the layer thickness... do jd=1,par%nd_var+1 if (sum(Sm(jd,:))>0) then pb(jd,:) = Sm(jd,:)/sum(Sm(jd,:)) else pb(jd,:) = pb(jd,:) endif dz(jd) = sum(Sm(jd,:))/(1.d0-par%por) enddo !!! modify grid !merge two upper layers in case of erosion if (dz(par%nd_var) .lt. par%merge*dz(par%nd_var+1)) then forall (jg=1:par%ngd) pb(par%nd_var,jg) = (dz(par%nd_var)*pb(par%nd_var,jg) + dz(par%nd_var+1)* & pb(par%nd_var+1,jg))/(dz(par%nd_var)+dz(par%nd_var+1)) pb(par%nd_var+1:s%nd(i,j)-1,jg) = pb(par%nd_var+2:s%nd(i,j),jg) pb(s%nd(i,j),jg) = pb(s%nd(i,j),jg) endforall s%z0bed(i,j) = s%z0bed(i,j)-dz(par%nd_var+1) dz(par%nd_var) = dz(par%nd_var+1)+dz(par%nd_var) endif !split upper layer in case of sedimentation if (dz(par%nd_var)>par%split*dz(par%nd_var+1)) then pb(par%nd_var+1:s%nd(i,j),:) = pb(par%nd_var:s%nd(i,j)-1,:) s%z0bed(i,j) = s%z0bed(i,j)+dz(par%nd_var+1) dz(par%nd_var) = dz(par%nd_var)-dz(par%nd_var+1) endif enddo ! nt_sub pb = max(0.d0,min(pb,1.d0)) zbold = s%zb(i,j) s%zb(i,j) = s%z0bed(i,j)+sum(dz) s%dzbnow(i,j) = s%dzbnow(i,j)+(s%zb(i,j)-zbold) s%sedero(i,j) = s%sedero(i,j)+(s%zb(i,j)-zbold) s%structdepth(i,j) = max(0.d0,s%structdepth(i,j)+(s%zb(i,j)-zbold)) end subroutine update_fractions !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine sedtransform(s,par) use params use spaceparams use readkey_module use xmpi_module implicit none type(spacepars),target :: s type(parameters) :: par integer :: i,j,jg real*8 :: z0,dcf,dcfin,ML real*8 :: Te,Sster,cc1,cc2,wster,Ass real*8 :: kl,alpha,alpha1,alpha2,beta,psi real*8 , save :: delta,kvis,onethird,twothird,phi real*8 , dimension(:),allocatable ,save :: dster,ws0 real*8 , dimension(:,:),allocatable ,save :: vmg,Cd,Asb,dhdx,dhdy,Ts,hfac real*8 , dimension(:,:),allocatable ,save :: urms2,Ucr,Ucrc,Ucrw,term1,B2,srfTotal,srfRhee,vero,Ucrb,Ucrs real*8 , dimension(:,:),allocatable ,save :: uandv,b,fslope,hloc,ceqs,ceqb,fallvelredfac real*8 , dimension(:,:,:),allocatable,save :: w !include 's.ind' !include 's.inp' if (.not. allocated(vmg)) then allocate (vmg (s%nx+1,s%ny+1)) allocate (term1 (s%nx+1,s%ny+1)) allocate (B2 (s%nx+1,s%ny+1)) allocate (Cd (s%nx+1,s%ny+1)) allocate (Asb (s%nx+1,s%ny+1)) allocate (Ucr (s%nx+1,s%ny+1)) allocate (Ucrc (s%nx+1,s%ny+1)) allocate (Ucrw (s%nx+1,s%ny+1)) allocate (urms2 (s%nx+1,s%ny+1)) allocate (hloc (s%nx+1,s%ny+1)) allocate (Ts (s%nx+1,s%ny+1)) allocate (ceqs (s%nx+1,s%ny+1)) allocate (ceqb (s%nx+1,s%ny+1)) allocate (srfTotal(s%nx+1,s%ny+1)) ! Lodewijk allocate (Ucrb (s%nx+1,s%ny+1)) ! Lodewijk allocate (Ucrs (s%nx+1,s%ny+1)) ! Lodewijk allocate (srfRhee(s%nx+1,s%ny+1)) ! Lodewijk allocate (vero (s%nx+1,s%ny+1)) ! Lodewijk allocate (fallvelredfac(s%nx+1,s%ny+1)) ! Lodewijk allocate (w (s%nx+1,s%ny+1,par%ngd)) ! Lodewijk allocate (dster (par%ngd)) allocate (ws0 (par%ngd)) allocate (dhdx (s%nx+1,s%ny+1)) allocate (dhdy (s%nx+1,s%ny+1)) allocate (uandv (s%nx+1,s%ny+1)) allocate (b (s%nx+1,s%ny+1)) allocate (fslope(s%nx+1,s%ny+1)) allocate (hfac (s%nx+1,s%ny+1)) vmg = 0.d0 onethird=1.d0/3.d0 twothird=2.d0/3.d0 phi = par%reposeangle/180*par%px ! Angle of internal friction ! Robert: do only once, not necessary every time step do jg=1,par%ngd ! cjaap: compute fall velocity with simple expression from Ahrens (2000) Te = 20.d0 kvis = 4.d0/(20.d0+Te)*1d-5 ! Van rijn, 1993 Sster = s%D50(jg)/(4*kvis)*sqrt((par%rhos/par%rho-1)*par%g*s%D50(jg)) cc1 = 1.06d0*tanh(0.064d0*Sster*exp(-7.5d0/Sster**2)) cc2 = 0.22d0*tanh(2.34d0*Sster**(-1.18d0)*exp(-0.0064d0*Sster**2)) wster = cc1+cc2*Sster ws0(jg) = wster*sqrt((par%rhos/par%rho-1.d0)*par%g*s%D50(jg)) ! RJ: for modeling gravel delta = (par%rhos-par%rho)/par%rho dster(jg)=(delta*par%g/1.d-12)**onethird*s%D50(jg) if (par%fallvelred==0) then w(:,:,jg) = ws0(jg) endif enddo endif ! Lodewijk: calculate fall velocity reduction coefficient based on concentration of previous time step ! Rowe (1987) made an estimation of the exponent alpha by fitting a logarithmic function on a dataset of Richardson and Zaki (1954). if (par%fallvelred==1) then do jg=1,par%ngd alpha = 2.35d0*(2.d0+0.175d0*(ws0(jg)*s%D50(jg)/kvis)**0.75d0)/(1.d0+0.175d0*(ws0(jg)*s%D50(jg)/kvis)**0.75d0) fallvelredfac = (1.d0-s%ccg(:,:,jg))**(alpha) w(:,:,jg) = ws0(jg)*fallvelredfac enddo endif ! ! hloc = max(hh,0.01d0) ! Jaap hloc = max(s%hh,0.01) ! Compute mean fall velocity par%ws=0.d0 do jg=1,par%ngd par%ws=par%ws+w(1,1,jg) enddo par%ws=par%ws/par%ngd ! ! compute near bed turbulence ! ! due to short waves if (par%swave==1) then ! wave breaking induced turbulence due to short waves do j=1,s%ny+1 do i=1,s%nx+1 ! compute mixing length ! ML = 2*R(i,j)*par%Trep/(par%rho*c(i,j)*max(H(i,j),par%eps)) ML = dsqrt(2*s%R(i,j)*par%Trep/(par%rho*max(s%c(i,j),1d-10))) ! ML = 0.9d0*H(i,j) ML = min(ML,hloc(i,j)); ! exponential decay turbulence over depth dcfin = exp(min(100.d0,hloc(i,j)/max(ML,0.01d0))) dcf = min(1.d0,1.d0/(dcfin-1.d0)) ! Jaap: new approach: compute kb based on waveturb result s%kb(i,j) = s%kturb(i,j)*dcf ! if (par%turb == TURB_BORE_AVERAGED) then s%kb(i,j) = s%kb(i,j)*par%Trep/s%Tbore(i,j) endif enddo ! Jaap: rundown jet creating additional turbulence if (par%swrunup==1)then s%kb(nint(s%istruct(j)),j) = s%kb(nint(s%istruct(j)),j) + par%jetfac* & (s%E(nint(s%istruct(j)),j)*s%strucslope(j)*sqrt(par%g/s%hh(nint(s%istruct(j)),j)))**twothird endif enddo elseif (par%nonh==1) then do j=1,s%ny+1 do i=1,s%nx+1 !ML=max(s%rolthick(i,j),0.07d0*max(s%hh(i,j),par%hmin)) !s%kb(i,j) = s%kturb(i,j)*ML/max(s%hh(i,j),par%hmin) ! Simpler expression s%kb(i,j) = s%kturb(i,j) ! even simpler expression :-) enddo enddo endif !par%swave == 1 ! switch to include long wave stirring if (par%lws==1) then vmg = dsqrt(s%ue**2+s%ve**2) elseif (par%lws==0) then ! vmg lags on actual mean flow; but long wave contribution to mean flow is included... vmg = (1.d0-1.d0/par%cats/par%Trep*par%dt)*vmg + (1.d0/par%cats/par%Trep*par%dt)*dsqrt(s%ue**2+s%ve**2) endif urms2 = s%urms**2.d0+1.45d0*s%kb do jg = 1,par%ngd Ts = par%tsfac*hloc/w(:,:,jg) s%Tsg(:,:,jg) = max(Ts,par%Tsmin) ! ! calculate treshold velocity Ucr ! if (par%form==FORM_SOULSBY_VANRIJN) then ! Soulsby van Rijn if(s%D50(jg)<=0.0005d0) then Ucr=0.19d0*s%D50(jg)**0.1d0*log10(4.d0*hloc/s%D90(jg)) else !Dano see what happens with coarse material Ucr=8.5d0*s%D50(jg)**0.6d0*log10(4.d0*hloc/s%D90(jg)) end if elseif (par%form==FORM_VANTHIEL_VANRIJN) then ! Van Thiel de Vries & Reniers 2008 if(s%D50(jg)<=0.0005) then Ucrc=0.19d0*s%D50(jg)**0.1d0*log10(4.d0*hloc/s%D90(jg)) !Shields Ucrw=0.24d0*(delta*par%g)**0.66d0*s%D50(jg)**0.33d0*par%Trep**0.33d0 !Komar and Miller (1975) else if(s%D50(jg)<=0.002) then Ucrc=8.5d0*s%D50(jg)**0.6d0*log10(4.d0*hloc/s%D90(jg)) !Shields Ucrw=0.95d0*(delta*par%g)**0.57d0*s%D50(jg)**0.43*par%Trep**0.14 !Komar and Miller (1975) else if(s%D50(jg)>0.002) then Ucrc=1.3d0*sqrt(delta*par%g*s%D50(jg))*(hloc/s%D50(jg))**(0.5d0*onethird) !Maynord (1978) --> also Neill (1968) where 1.3d0 = 1.4d0 Ucrw=0.95d0*(delta*par%g)**0.57d0*s%D50(jg)**0.43*par%Trep**0.14 !Komar and Miller (1975) end if B2 = vmg/max(vmg+dsqrt(urms2),par%eps) Ucr = B2*Ucrc + (1-B2)*Ucrw !Van Rijn 2007 (Bed load transport paper) end if ! Lodewijk: implementation Rhee (2010), reduction sediment transport due dilatancy srfRhee(:,:) = 0.d0 srfTotal(:,:) = 1.d0 if (par%dilatancy == 1) then vero(:,:) = max(0.d0,-1.d0*s%dzbdt) ! Erosion velocity, for now asume it equal to -s%dzbdt of the previous time step kl = par%g/(160.d0*kvis)*(s%D15(jg)**2.d0)*((par%por**3.d0)/(1.d0-par%por)**2.d0) ! Permeability, Adel 1987, which is based on the Ergun equation ! bed porosity, estimated here as maximum porosity, to be added to user input ! A=3/4 for single particles and A=1/(1-n0) for a continuum ! Reduction factor on the critical Shields parameter by dilatancy (Van Rhee, 2010) srfRhee(:,:) = vero(:,:)/kl*(par%pormax-par%por)/(1.d0-par%pormax)*par%rheeA/delta endif ! Lodewijk: implementation bed slope reduction on critical flow velocity if (par%bdslpeffini == BDSLPEFFINI_NONE) then srfTotal(:,:) = 1.d0 + srfRhee(:,:) elseif (par%bdslpeffini == BDSLPEFFINI_TOTAL .or. par%bdslpeffini == BDSLPEFFINI_BED) then do j=1,s%ny+1 do i=1,s%nx+1 ! Prevent NaN values if too small values if ((dabs(s%ue(i,j)) > 0.000001d0 .or. dabs(s%ve(i,j)) > 0.000001d0) .and. & (dabs(s%dzbdx(i,j))>0.000001d0 .or. dabs(s%dzbdy(i,j))>0.000001d0)) then ! Angle between the x-axis and the flow velocity ! REMARK: also include waves in the velocity? if (s%ue(i,j) < 0.d0) then alpha1 = datan(s%ve(i,j)/s%ue(i,j)) + par%px else alpha1 = datan(s%ve(i,j)/s%ue(i,j)) endif ! Angle between the x-axis and the bed slope vector directed in down-slope direction, derived in thesis Lodewijk if (s%dzbdy(i,j) >= 0.d0) then alpha2 = -datan(s%dzbdx(i,j)/s%dzbdy(i,j))+1.5d0*par%px else alpha2 = -datan(s%dzbdx(i,j)/s%dzbdy(i,j))+0.5d0*par%px endif psi = alpha1-(alpha2-par%px) ! Angle between the flow direction and the up-slope directed vector if (dabs(s%dzbdx(i,j))<0.000001d0) then ! A smaller slope could result in a NaN for beta ! Beta purely based on dzbdy beta = datan(dabs(s%dzbdy(i,j))) else beta = datan(dabs(s%dzbdx(i,j)/dsin(datan(s%dzbdx(i,j)/s%dzbdy(i,j))))) ! Maximum absolute bed slope angle, derived in thesis Lodewijk endif beta = min(beta,phi) ! Take min to exclude NaN's if (par%dilatancy == 1) then srfTotal(i,j) = (dcos(psi)*dsin(beta)+dsqrt( & ((srfRhee(i,j))**2+2*srfRhee(i,j)*dcos(beta)+dcos(beta)**2 & )*dtan(phi)**2-dsin(psi)**2*dsin(beta)**2 & ))/dtan(phi) ! Soulsby (1997), modified by Lodewijk (see Thesis) else srfTotal(i,j) = (dcos(psi)*dsin(beta)+ & dsqrt(dcos(beta)**2*dtan(phi)**2-dsin(psi)**2*dsin(beta)**2))/dtan(phi) ! Soulsby (1997) endif endif end do end do endif ! Calculate the new critical velocity based on the modification factors on the Shields parameter, bed slope only on bed load Ucrb(:,:) = Ucr(:,:)*sqrt(srfTotal) ! Lodewijk if (par%bdslpeffini == BDSLPEFFINI_TOTAL) then Ucrs = Ucrb else Ucrs = Ucr*(1.d0+sqrt(srfRhee)) ! Lodewijk, no effect on suspended load by bed slope endif if (par%form==FORM_SOULSBY_VANRIJN) then ! Soulsby van Rijn ! ! drag coefficient z0 = par%z0 Cd=(0.40d0/(log(max(hloc,10.d0*z0)/z0)-1.0d0))**2 ! ! transport parameters Asb=0.005d0*hloc*(s%D50(jg)/hloc/(delta*par%g*s%D50(jg)))**1.2d0 ! bed load coefficent Ass=0.012d0*s%D50(jg)*dster(jg)**(-0.6d0)/(delta*par%g*s%D50(jg))**1.2d0 ! suspended load coeffient ! term1=(vmg**2+0.018/Cd*par%sws*urms2) ! Make 0.018/Cd is always smaller than the flow friction coefficient ! ! reduce sediment suspensions for (inundation) overwash conditions with critical flow velocitties ! vmag2=min(vmag2,par%smax*par%C**2*D50(jg)*delta) ! vmag2=min(vmag2,par%smax*par%g/par%cf*D50(jg)*delta) ! In terms of cf ! term1=sqrt(vmag2+0.018d0/Cd*urms2) ! nearbed-velocity ! the two above lines are comment out and replaced by a limit on total velocity u2+urms2, robert 1/9 and ap 28/11 ! term1=min(term1,par%smax*par%g/s%cf*s%D50(jg)*delta) term1=sqrt(term1) ! ceqb = 0.d0*term1 !initialize ceqb ceqs = 0.d0*term1 !initialize ceqs do j=1,s%ny+1 do i=1,s%nx ! Lodewijk: sepperate bed load from suspended load since Ucr not anymore the same if(term1(i,j)>Ucrb(i,j) .and. hloc(i,j)>par%eps) then ceqb(i,j)=Asb(i,j)*(term1(i,j)-Ucrb(i,j))**2.4d0 end if if(term1(i,j)>Ucrs(i,j) .and. hloc(i,j)>par%eps) then ceqs(i,j)=Ass*(term1(i,j)-Ucrs(i,j))**2.4d0 end if end do end do ! elseif (par%form==FORM_VANTHIEL_VANRIJN) then ! Van Thiel de Vries & Reniers 2008 ! ! transport parameters Asb=0.015d0*hloc*(s%D50(jg)/hloc)**1.2d0/(delta*par%g*s%D50(jg))**0.75d0 !bed load coefficent Ass=0.012d0*s%D50(jg)*dster(jg)**(-0.6d0)/(delta*par%g*s%D50(jg))**1.2d0 !suspended load coeffient ! ! Jaap: par%sws to set short wave stirring to zero ! Jaap: Van Rijn use Peak orbital flow velocity --> 0.64 corresponds to 0.4 coefficient regular waves Van Rijn (2007) term1=vmg**2+0.64d0*par%sws*urms2 ! reduce sediment suspensions for (inundation) overwash conditions with critical flow velocitties term1=min(term1,par%smax*par%g/max(s%cf,1d-10)*s%D50(jg)*delta) term1=sqrt(term1) ! ceqb = 0.d0*term1 !initialize ceqb ceqs = 0.d0*term1 !initialize ceqs do j=1,s%ny+1 do i=1,s%nx ! Lodewijk: sepperate bed load from suspended load since Ucr not anymore the same if(term1(i,j)>Ucrb(i,j) .and. hloc(i,j)>par%eps) then ceqb(i,j)=Asb(i,j)*(term1(i,j)-Ucrb(i,j))**1.5 end if if(term1(i,j)>Ucrs(i,j) .and. hloc(i,j)>par%eps) then ceqs(i,j)=Ass*(term1(i,j)-Ucrs(i,j))**2.4 end if end do end do ! end if ! ceqb = min(ceqb/hloc,par%cmax/2) ! maximum equilibrium bed concentration s%ceqbg(:,:,jg) = (1-par%bulk)*ceqb*s%sedcal(jg)*s%wetz ceqs = min(ceqs/hloc,par%cmax/2) ! maximum equilibrium suspended concentration s%ceqsg(:,:,jg) = (ceqs+par%bulk*ceqb)*s%sedcal(jg)*s%wetz ! Jaap: old brute method to prevent strong coastline erosion ! where (hloc<=par%hmin) ceqsg(:,:,jg) = 0.d0 enddo ! end of grain size classes ! end of grain size classes end subroutine sedtransform !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine waveturb(s,par) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Copyright (C) 2007 UNESCO-IHE, WL|Delft Hydraulics and Delft University ! ! Dano Roelvink, Ap van Dongeren, Ad Reniers, Jamie Lescinski, ! ! Jaap van Thiel de Vries, Robert McCall ! ! ! ! d.roelvink@unesco-ihe.org ! ! UNESCO-IHE Institute for Water Education ! ! P.O. Box 3015 ! ! 2601 DA Delft ! ! The Netherlands ! ! ! ! This library is free software; you can redistribute it and/or ! ! modify it under the terms of the GNU Lesser General Public ! ! License as published by the Free Software Foundation; either ! ! version 2.1 of the License, or (at your option) any later version. ! ! ! ! This library is distributed in the hope that it will be useful, ! ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! ! Lesser General Public License for more details. ! ! ! ! You should have received a copy of the GNU Lesser General Public ! ! License along with this library; if not, write to the Free Software ! ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 ! ! USA ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use params use spaceparams use xmpi_module implicit none type(spacepars),target :: s type(parameters) :: par integer :: i integer :: j real*8 :: ML, disturb real*8, save :: twothird real*8,dimension(:,:),allocatable,save :: ksource, kturbu,kturbv,Sturbu,Sturbv,dzsdt_cr !include 's.ind' !include 's.inp' if (.not. allocated(kturbu)) then allocate(ksource (s%nx+1,s%ny+1)) allocate(kturbu (s%nx+1,s%ny+1)) allocate(kturbv (s%nx+1,s%ny+1)) allocate(Sturbu (s%nx+1,s%ny+1)) allocate(Sturbv (s%nx+1,s%ny+1)) allocate(dzsdt_cr (s%nx+1,s%ny+1)) twothird=2.d0/3.d0 endif ! use lagrangian velocities kturbu = 0.0d0 !Jaap kturbv = 0.0d0 !Jaap ! dzsdt_cr=par%beta*s%c dzsdt_cr=par%beta*sqrt(par%g*s%hh) ! Update roller thickness !s%rolthick=s%rolthick+par%dt*(abs(s%dzsdt)-dzsdt_cr) Dano: not abs! s%rolthick=s%rolthick+par%dt*(s%dzsdt-dzsdt_cr) s%rolthick=max(s%rolthick,0.d0) s%rolthick=min(s%rolthick,s%hh) where (s%wetz==0) s%rolthick=0.d0 endwhere ! Jaap compute sources and sinks ! long wave source ksource = 0 if (par%lwt == 1) then ksource = par%g*s%rolthick*par%beta*sqrt(par%g*s%hh) !only important in shallow water, where s%c=sqrt(gh) endif if (par%turbadv == TURBADV_NONE) then s%kturb = (s%DR/par%rho)**twothird ! See Battjes, 1975 / 1985 else ksource = ksource + s%DR/par%rho ! Jaap do long wave turb approach for both short waves and long waves ! ! Turbulence in uu-points ! do j=1,s%ny+1 do i=1,s%nx if(s%uu(i,j)>0.d0) then kturbu(i,j)=par%thetanum*s%kturb(i,j)+(1.d0-par%thetanum)*s%kturb(min(i+1,s%nx),j) elseif(s%uu(i,j)<0.d0) then kturbu(i,j)=par%thetanum*s%kturb(i+1,j)+(1.d0-par%thetanum)*s%kturb(max(i,2),j) else kturbu(i,j)=0.5d0*(s%kturb(i,j)+s%kturb(i+1,j)) endif enddo enddo if (xmpi_isbot) kturbu(s%nx+1,:) = s%kturb(s%nx+1,:) ! ! Turbulence in vv-points ! do j=1,s%ny do i=1,s%nx+1 if(s%vv(i,j)>0) then kturbv(i,j)=par%thetanum*s%kturb(i,j)+(1.d0-par%thetanum)*s%kturb(i,min(j+1,s%ny)) else if(s%vv(i,j)<0) then kturbv(i,j)=par%thetanum*s%kturb(i,j+1)+(1.d0-par%thetanum)*s%kturb(i,max(j,2)) else kturbv(i,j)=0.5d0*(kturbv(i,j)+kturbv(i,j+1)) endif enddo enddo kturbv(:,s%ny+1) = s%kturb(:,s%ny+1) !Robert ! ! Turbulence advection in X and Y direction ! if (par%turbadv == TURBADV_LAGRANGIAN) then Sturbu=kturbu*s%uu*s%hu*s%wetu Sturbv=kturbv*s%vv*s%hv*s%wetv elseif (par%turbadv == TURBADV_EULERIAN) then Sturbu=kturbu*s%ueu*s%hu*s%wetu Sturbv=kturbv*s%vev*s%hv*s%wetv endif ! ! Update turbulence ! if (s%ny>0) then do j=2,s%ny+1 do i=2,s%nx+1 if (par%betad>0) then disturb=par%betad*s%kturb(i,j)**1.5d0 else ML=max(s%rolthick(i,j),0.07d0*max(s%hh(i,j),par%hmin)) disturb=0.08d0*max(s%hh(i,j),par%hmin)/ML*s%kturb(i,j)**1.5d0 endif s%kturb(i,j) = (s%hold(i,j)*s%kturb(i,j)-par%dt*( & (Sturbu(i,j)*s%dnu(i,j)-Sturbu(i-1,j)*s%dnu(i-1,j)+& Sturbv(i,j)*s%dsv(i,j)-Sturbv(i,j-1)*s%dsv(i,j-1))*s%dsdnzi(i,j)-& (ksource(i,j)-disturb)))/max(s%hh(i,j),par%hmin) s%kturb(i,j)=max(s%kturb(i,j),0.0d0) enddo enddo else j=1 do i=2,s%nx+1 s%kturb(i,j) = (s%hold(i,j)*s%kturb(i,j)-par%dt*( & (Sturbu(i,j)*s%dnu(i,j)-Sturbu(i-1,j)*s%dnu(i-1,j))*s%dsdnzi(i,j)-& (ksource(i,j)-par%betad*s%kturb(i,j)**1.5d0)))/max(s%hh(i,j),par%hmin) s%kturb(i,j)=max(s%kturb(i,j),0.0d0) enddo endif s%kturb = s%kturb/max(s%hh,0.01d0) ! Jaap only required for advection mode? if (xmpi_istop) s%kturb(1,:)=s%kturb(2,:) if (xmpi_isbot) s%kturb(s%nx+1,:)=s%kturb(s%nx+1-1,:) if (s%ny>0) then if (xmpi_isleft) s%kturb(:,1)=s%kturb(:,2) if (xmpi_isright) s%kturb(:,s%ny+1)=s%kturb(:,s%ny+1-1) endif #ifdef USEMPI call xmpi_shift_ee(s%kturb) #endif endif ! turbadv == 'none' end subroutine waveturb !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine RvR(s,par) use params use spaceparams use xmpi_module implicit none type(spacepars),target :: s type(parameters) :: par real*8 , save :: m1,m2,m3,m4,m5,m6,alpha,beta real*8 , dimension(:,:),allocatable,save :: Urs,Bm,B1 !include 's.ind' !include 's.inp' ! only in first timestep.. if (.not. allocated(Urs)) then allocate (Urs (s%nx+1,s%ny+1)) allocate (Bm (s%nx+1,s%ny+1)) allocate (B1 (s%nx+1,s%ny+1)) m1 = 0; ! a = 0 m2 = 0.7939; ! b = 0.79 +/- 0.023 m3 = -0.6065; ! c = -0.61 +/- 0.041 m4 = 0.3539; ! d = -0.35 +/- 0.032 m5 = 0.6373; ! e = 0.64 +/- 0.025 m6 = 0.5995; ! f = 0.60 +/- 0.043 alpha = -log10(exp(1.d0))/m4 beta = exp(m3/m4) endif Urs = 3.d0/8.d0*sqrt(2.d0)*s%H*s%k/(s%k*s%hh)**3 !Ursell number Urs = max(Urs,0.000000000001d0) Bm = m1 + (m2-m1)/(1.d0+beta*Urs**alpha) !Boltzmann sigmoid (eq 6) B1 = (-90.d0+90.d0*tanh(m5/Urs**m6))*par%px/180.d0 s%Sk = Bm*cos(B1) !Skewness (eq 8) s%As = Bm*sin(B1) !Asymmetry(eq 9) s%ua = par%sws*(par%facSk*s%Sk-par%facAs*s%As)*s%urms ! multiply Sk and As with wetz to get zeros at h = 0 for output s%Sk = s%Sk*s%wetz s%As = s%As*s%wetz end subroutine RvR !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine vT(s,par) use params use spaceparams use readkey_module use xmpi_module implicit none type(spacepars),target :: s type(parameters) :: par integer :: i,j integer , save :: nh,nt integer :: ih0,it0,ih1,it1 real*8 :: p,q,f0,f1,f2,f3 real*8 :: t0fac,siguref,duddtmax,dudtmax,duddtmean,dudtmean,detadxmean real*8 , save :: dh,dt real*8 , dimension(:,:),allocatable ,save :: h0,t0,detadxmax ! Robert: RF table now included in source code, rather than read from file ! Rienecker Fenton table with amongst others amplitudes non-linear components obtained with stream function theory include 'RF.inc' ! Robert: 'RF.inc' contains definition of RF as real*8(5,33,40) with "parameter" attribute ! so RF values may not be modified! To save memory, only rows 13,14,15,16 and 18 of the ! original matrix are stored. So new row 1 corresponds with old row 13, etc. !include 's.ind' !include 's.inp' ! only in first timestep.. if (.not. allocated(h0)) then allocate (h0 (s%nx+1,s%ny+1)) allocate (t0 (s%nx+1,s%ny+1)) allocate (detadxmax (s%nx+1,s%ny+1)) dh = 0.03d0 dt = 1.25d0 nh = floor(0.99d0/dh); nt = floor(50.d0/dt); endif ! non-linearity of short waves is listed in table as function of dimensionless wave height h0 and dimensionless wave period t0 ! compute dimensionless wave height and wave period in each grid point.. h0 = min(nh*dh,max(dh,min(s%H,s%hh)/s%hh)) t0 = min(nt*dt,max(dt,par%Trep*sqrt(par%g/s%hh))) ! estimate Sk, As and ua by interpolating table values do j=1,s%ny+1 do i=1,s%nx+1 ! interpolate table values.... ih0=floor(h0(i,j)/dh); it0=floor(T0(i,j)/dt); ih1=min(ih0+1,nh); it1=min(it0+1,nt); p=(h0(i,j)-ih0*dh)/dh; q=(T0(i,j)-it0*dt)/dt; f0=(1-p)*(1-q); f1=p*(1-q); f2=q*(1-p); f3=p*q; ! Skewness and assymetry s%Sk(i,j) = f0*RF(1,ih0,it0)+f1*RF(1,ih1,it0)+ f2*RF(1,ih0,it1)+f3*RF(1,ih1,it1) s%As(i,j) = f0*RF(2,ih0,it0)+f1*RF(2,ih1,it0)+ f2*RF(2,ih0,it1)+f3*RF(2,ih1,it1) ! Sediment advection velocity from Skewness and Assymetry ! ua(i,j) = par%sws*par%facua*(Sk(i,j)-As(i,j))*urms(i,j) s%ua(i,j) = par%sws*(par%facSk*s%Sk(i,j)-par%facAs*s%As(i,j))*s%urms(i,j) ! Estimate bore period Tbore and mean slope bore front to feeded back in roller energy balance ! correct slope in case 1.25>T0>50 if (t0(i,j)==50.d0) then t0fac = 50.d0/max((par%Trep*sqrt(par%g/s%hh(i,j))),50.d0) elseif (t0(i,j)==1.25)then t0fac = 1.25d0/min((par%Trep*sqrt(par%g/s%hh(i,j))),1.25d0) else t0fac = 1.d0 endif ! detadxmax for Tbore... ! dimnesionless maximum acceleration under bore front duddtmax = f0*RF(3,ih0,it0)+f1*RF(3,ih1,it0)+ f2*RF(3,ih0,it1)+f3*RF(3,ih1,it1) siguref = f0*RF(4,ih0,it0)+f1*RF(4,ih1,it0)+ f2*RF(4,ih0,it1)+f3*RF(4,ih1,it1) ! translate dimensionless duddtmax to real world dudtmax ! /scale with variance and go from [-] to [m/s^2] /tableb./dimensionless dudtmax dudtmax = s%urms(i,j)/max(par%eps,siguref)*sqrt(par%g/s%hh(i,j))*t0fac*duddtmax detadxmax(i,j) = dudtmax*sinh(s%k(i,j)*s%hh(i,j))/max(max(s%c(i,j),sqrt(s%H(i,j)*par%g)),1d-10)/s%sigm(i,j) ! detadxmean for roller energy balance dissipation... if (par%rfb==1) then duddtmean = f0*RF(5,ih0,it0)+f1*RF(5,ih1,it0)+ f2*RF(5,ih0,it1)+f3*RF(5,ih1,it1) dudtmean = s%urms(i,j)/max(par%eps,siguref)*sqrt(par%g/s%hh(i,j))*t0fac*duddtmean detadxmean = dudtmean*sinh(s%k(i,j)*s%hh(i,j))/max(s%c(i,j),sqrt(s%H(i,j)*par%g))/s%sigm(i,j) s%BR(i,j) = par%BRfac*sin(atan(detadxmean)) endif enddo enddo s%Tbore = max(par%Trep/25.d0,min(par%Trep/4.d0,s%H/(max(max(s%c,sqrt(s%H*par%g)),1d-10)*max(detadxmax,par%eps)))) s%Tbore = par%Tbfac*s%Tbore ! multiply Sk and As with wetz to get zeros at h = 0 for output s%Sk = s%Sk*s%wetz s%As = s%As*s%wetz end subroutine vT subroutine setbathy_update(s, par) use params use spaceparams use interp implicit none type(spacepars) :: s type(parameters) :: par integer :: i,j,dummy real*8,dimension(s%nx+1,s%ny+1) :: zbnew ! interpolate from file do j=1,s%ny+1 do i=1,s%nx+1 call LINEAR_INTERP(s%tsetbathy,s%setbathy(i,j,:),par%nsetbathy, & par%t,zbnew(i,j),dummy) enddo enddo ! update water level s%zs = s%zs+zbnew-s%zb ! update bed level s%zb = zbnew ! update wet and dry cells where (s%zs 0.15). ! If not present Hrunup will converge to H at the water line (where H = 0 per definition) do j1=indx,i-1 ! TODO: Strange condition. In case of no toe, or no steep slope, ! there will still be extra turbulence at L1 meter from the water line... ! cross shore location structure toe if (s%dzbdx(j1,j)<0.15d0 .or. s%structdepth(j1,j)>0.1d0) then indx = j1 endif enddo ! update Hrunup and runup x-location s%Hrunup(j) = s%H(indx,j) ! short wave height at revetment toe s%xHrunup(j) = s%xz(indx,j) ! cross-shore location revetment toe s%istruct(j) = indx*1.d0 ! cross-shore index revetment toe s%iwl(j) = (i-1)*1.d0 ! cross-shore location waterline (inlcuding lw-s%runup) ! now iteratively compute runup hav1d = s%hh(:,j) runup_old = huge(0.d0) s%runup(j) = 0; nIter = 0; maxIter = 50; do while (abs(s%runup(j)-runup_old)>0.01d0 .and. nIter < maxIter) nIter = nIter +1; runup_old = s%runup(j) slopeind = 0 where (hav1d>par%eps .and. s%dzbdx(:,j)>0.15d0) slopeind = 1 endwhere s%strucslope(j) = sum(s%dzbdx(indx:s%nx,j)*s%dsu(indx:s%nx,j)*slopeind(indx:s%nx))/ & max(par%eps,sum(s%dsu(indx:s%nx,j)*slopeind(indx:s%nx))) if (s%strucslope(j) > 0.d0) then irrb = s%strucslope(j)/sqrt(2*par%px*max(s%Hrunup(j),par%eps)/par%g/par%Trep**2) s%runup(j) = par%facrun*min(irrb,2.3d0)*s%Hrunup(j)*cos(2*par%px/par%Trep*par%t) else s%runup(j) = 0.d0; endif ! This triggers s%runup to be calculated for the complete hinterland (also behind a dune). ! Only values calculated for the first dune front are used in the avalanching algorithm (morphevolution) hav1d = s%wetz(:,j)*max(par%eps,(s%hh(:,j) + s%runup(j))) + & (1.d0-s%wetz(:,j))*max(par%eps,s%runup(j)+s%zs(i-1,j)-s%zb(:,j) ) enddo endif enddo enddo end subroutine hybrid !Introduction of Thermal Routines built by UAA Staff !Thermal Model attempts to add thermal capabilities to the Xbeach mechanical erosion model !First routine addresses heat transfer from air to bluff - Updated 05/30/2017 ! ! This routine is no longer used as improvements have been implemented to improve on heat transfer calculations!!!!! ! subroutine HT_air_bluff(s,par,hsum) ! ! use params ! use interp ! use spaceparams ! use xmpi_module ! use paramsconst ! ! implicit none ! ! type(spacepars),target :: s ! type(parameters) :: par ! ! real*8,dimension(320,1) :: hsum ! real*8,dimension(:,:),allocatable :: hhmean !,tdc_cond ! real*8 :: has,Qas,Lb ! integer :: i,j ! ! has = 30.d0 ! Lb=0.d0 ! Qas=0.d0 ! allocate(hhmean(s%nx+1,s%ny+1)) ! hhmean(i,j)=0.d0 ! ! do j=1,s%ny+1 ! do i=1,s%nx+1 ! hhmean(i,j)=hsum(i,j)/par%counter_ub ! if(hhmean(i,j)>par%water_cut_off) then ! s%structdepth(i,j)=s%structdepth(i,j)+0 ! else if(hhmean(i,j)<=par%water_cut_off) then ! Qas=has*(par%Ta-par%Ts_bluff)*par%morfac ! Lb = par%nb*par%qi*par%Li ! s%structdepth(i,j) = s%structdepth(i,j)+((Qas/Lb)*(par%timecalc*par%morfac)) ! endif ! enddo ! enddo ! deallocate(hhmean) ! deallocate(tdc_cond) ! ! end subroutine HT_air_bluff This routine is no longer used as improvements have been implemented into subroutine heat_transfer subroutine heat_transfer(s,par,zbsum,hsum,Hwavesum,taubxsum,vmagsum,glmsum,eulsum,zssum,structdepth_out,tracker_vert_change) use params use interp use spaceparams use xmpi_module use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par integer :: i,j,i3,numb real*8,dimension(:,:),allocatable :: tracker_vert_change real*8,dimension(:,:),allocatable :: fwmean,ustarmean,R,E,Ubmean,zbmean,vmagmean,glmmean,eulmean,zsmean real*8,dimension(:,:),allocatable :: hhmean,taubxmean,Hwavemean,thaw_depth_change_vert,hws,hsum,Hwavesum,taubxsum,zbsum,vmagsum,glmsum,eulsum,zssum,structdepth_out !real*8,dimension(365,1) :: real*8 :: P,Lb ! P = 0.d0 allocate(Ubmean(s%nx+1,s%ny+1)) allocate(fwmean(s%nx+1,s%ny+1)) allocate(ustarmean(s%nx+1,s%ny+1)) allocate(E(s%nx+1,s%ny+1)) allocate(hhmean(s%nx+1,s%ny+1)) allocate(zbmean(s%nx+1,s%ny+1)) allocate(zsmean(s%nx+1,s%ny+1)) allocate(Hwavemean(s%nx+1,s%ny+1)) allocate(taubxmean(s%nx+1,s%ny+1)) !allocate(Qwsed(s%nx+1,s%ny+1)) allocate(vmagmean(s%nx+1,s%ny+1)) allocate(glmmean(s%nx+1,s%ny+1)) allocate(eulmean(s%nx+1,s%ny+1)) allocate(hws(s%nx+1,s%ny+1)) !Qwsed(i,j)=0.0D+00 do j=1,s%ny+1 do i=1,s%nx+1 Ubmean(i,j)=0.0D+00 fwmean(i,j)=0.0D+00 ustarmean(i,j)=0.0D+00 E(i,j)=0.0D+00 hws(i,j)=0.0D+00 hhmean(i,j)=0.0D+00 zbmean(i,j)=0.0D+00 zsmean(i,j)=0.0D+00 Hwavemean(i,j)=0.0D+00 taubxmean(i,j)=0.0D+00 vmagmean(i,j)=0.0D+00 glmmean(i,j)=0.0D+00 eulmean(i,j)=0.0D+00 end do end do P=(par%v*par%Cw)/par%Kw ! if(par%nonh==0 .AND. par%swave==1) then ! do j=1,s%ny+1 do i=1,s%nx+1 hhmean(i,j)=hsum(i,j)/par%counter_ub Hwavemean(i,j)=Hwavesum(i,j)/par%counter_ub taubxmean(i,j)=taubxsum(i,j)/par%counter_ub vmagmean(i,j)=vmagsum(i,j)/par%counter_ub glmmean(i,j)=glmsum(i,j)/par%counter_ub eulmean(i,j)=eulsum(i,j)/par%counter_ub Ubmean(i,j)=0.5d0*Hwavemean(i,j)*(par%g/hhmean(i,j))**0.5d0 ustarmean(i,j)=(taubxmean(i,j)/par%qs)**0.5d0 E(i,j)=(0.52d0*(P**0.8d0)*((ustarmean(i,j)*par%ks)/par%v)**0.45d0) if(Ubmean(i,j)==0.0D+00 .OR. taubxmean(i,j)==0.0D+00) then fwmean(i,j)=0.0D+00 else fwmean(i,j) = (2.d0*taubxmean(i,j))/(par%qs*(Ubmean(i,j))**2.d0) end if if(hhmean(i,j)>par%water_cut_off) then hws(i,j) = (0.5d0*fwmean(i,j)*par%Cw*Ubmean(i,j))/(1+E(i,j)*(0.5d0*fwmean(i,j))**0.5d0) ! Lb = 0.0D+00 ! Lb =par%nb*par%Li ! Qwsed(i,j) = hws(i,j)*(par%Tw-par%Ts_bluff) ! s%structdepth(i,j) = s%structdepth(i,j)+((Qwsed(i,j)/Lb)*par%timecalc) else if (hhmean(i,j)<=par%water_cut_off) then hws(i,j) = 0.0D+00 ! s%structdepth(i,j)=s%structdepth(i,j)+0 endif enddo enddo hws(1,1)=hws(2,1) call csv_write(s,par,hws,Hwavemean,taubxmean,fwmean,Ubmean) !Run for swave=1 scenario !do j=1,s%ny+1 ! do i=1,s%nx+1 ! enddo ! enddo else if(par%nonh==1 .AND. par%swave==0) then ! Important one! !allocate(Qwsed(s%nx+1,s%ny+1)) !Qwsed(i,j)=0.d0 do j=1,s%ny+1 do i=1,s%nx+1 hhmean(i,j)=hsum(i,j)/par%counter_ub zbmean(i,j)=zbsum(i,j)/par%counter_ub zsmean(i,j)=zssum(i,j)/par%counter_ub vmagmean(i,j)=vmagsum(i,j)/par%counter_ub glmmean(i,j)=glmsum(i,j)/par%counter_ub eulmean(i,j)=eulsum(i,j)/par%counter_ub !!!!!!NEW!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if(hhmean(i,j)<=par%water_cut_off) then Hwavemean(i,j)=0.0D+00 else Hwavemean(i,j)=2*(zsmean(i,j)) endif Ubmean(i,j)=abs(0.5d0*Hwavemean(i,j)*(par%g/hhmean(i,j))**0.5d0) if(Ubmean(i,j)==0.0D+00) then fwmean(i,j)=0.0D+00 else if (Ubmean(i,j) > 0.0D+00) then fwmean(i,j) = 0.01d0 endif if (Ubmean(i,j) > par%ubset) then Ubmean(i,j) = Ubmean(i,j)*0 Ubmean(i,j) = par%ubset end if taubxmean(i,j) = (fwmean(i,j)*(par%qs*(Ubmean(i,j))**2.d0))/2.d0 ustarmean(i,j)=(taubxmean(i,j)/par%qs)**0.5d0 E(i,j)=(0.52d0*(P**0.8d0)*((ustarmean(i,j)*par%ks)/par%v)**0.45d0) if(hhmean(i,j)>par%water_cut_off) then hws(i,j) = (0.5d0*fwmean(i,j)*par%Cw*Ubmean(i,j))/(1+E(i,j)*(0.5d0*fwmean(i,j))**0.5d0) else if (hhmean(i,j)<=par%water_cut_off) then hws(i,j) = 0.0D+00 endif end do end do hws(1,1)=hws(2,1) call csv_write(s,par,hws,Hwavemean,taubxmean,fwmean,Ubmean) !call csv_read(hws) call Coast_Thermal(s,par,hws,zbmean,vmagmean,glmmean,eulmean,thaw_depth_change_vert) do j=1,s%ny+1 do i=1,s%nx+1 !hhmean(i,j)=hsum(i,j)/par%counter_ub s%structdepth(i,j) = s%structdepth(i,j)+thaw_depth_change_vert(i,j) tracker_vert_change(i,j) = tracker_vert_change(i,j) + thaw_depth_change_vert(i,j) end do end do else if(par%nonh==0 .AND. par%swave==0) then call csv_read(s,par,hws) call Coast_Thermal(s,par,hws,zbmean,vmagmean,glmmean,eulmean,thaw_depth_change_vert) do j=1,s%ny+1 do i=1,s%nx+1 if (hws(i,j)==0.0D+00) then s%structdepth(i,j)=s%structdepth(i,j)+thaw_depth_change_vert(i,j) end if end do end do end if deallocate(hws) deallocate(Ubmean) deallocate(fwmean) deallocate(ustarmean) deallocate(E) deallocate(hhmean) deallocate(zbmean) deallocate(vmagmean) deallocate(glmmean) deallocate(eulmean) deallocate(Hwavemean) deallocate(taubxmean) deallocate(zsmean) ! !deallocate(Qwsed) do j=1,s%ny+1 do i=1,s%nx+1 structdepth_out(i,j) = 0.0D+00 structdepth_out(i,j) = s%structdepth(i,j) !State of structural depth at end of thermal routine enddo enddo if ((par%nonh==1 .AND. par%swave==0) .OR. (par%nonh==0 .AND. par%swave==0)) then deallocate(thaw_depth_change_vert) end if ! if(s%hh(i,j)>par%water_cut_off) then ! fw(i,j) = (2*tb)/(q*Ub**2) !Ub(i,j)=sub(i,j)/counter_ub(i,j) !tb(i,j)=stb(i,j)/counter_ub(i,j) !R(i,j)=(ustar(i,j)*ks)/v !Thawing due to ambient air temperature ! material(i,j)=1 ! wet soil. Allowed to erode at this point. ! else if(s%hh(i,j)<=0) then ! material(i,j)=2 ! material(i,j)=1 ! wet soil. Allowed to erode at this point. ! material(i,j)=2 ! dry soil. Not allowed to erode at this point. ! endif ! do j=1,s%ny+1 ! do i=1,s%nx+1 ! if (material(i,j)==1) then !fw(i,j) = (2*tb)/(q*Ub**2) ! fw = 0.01 ! hws = (0.5*fw*Cw*Ub)/(1+E*(0.5*fw)**0.5) ! Qwsed = hws*(Tw-Ts_bed) ! s%structdepth(i,j) = (Qwsed/Lb)*par%timecalc ! else if (material(i,j)==2) then ! s%structdepth(i,j)=0 ! endif ! enddo ! enddo !s%structdepth(i,j) = s%thawdepth(i,j) end subroutine heat_transfer ! subroutine HT_water_soil(s,par,hsum) !Routine is not used. All current heat transfer components are available in subroutine heat_transfer use params use interp use spaceparams use xmpi_module use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par integer :: i,j,numb real*8 :: ynew,yold real*8,dimension(365,1) :: hsum real*8,dimension(:,:),allocatable :: hhmean allocate(hhmean(s%nx+1,s%ny+1)) hhmean(i,j)=0.d0 ynew=0.d0 yold=0.d0 do j=1,s%ny+1 do i=1,s%nx+1 hhmean(i,j)=hsum(i,j)/par%counter_ub if(hhmean(i,j)>par%water_cut_off) then ynew =-(2.5066D-32*(par%timet)**6)+(1.4379D-26*(par%timet)**5)-(3.2241D-21*(par%timet)**4)+(3.5801D-16*(par%timet)**3)-(2.0842D-11*(par%timet)**2)+(7.2871D-7*(par%timet)) yold = -(2.5066D-32*(par%timet-par%timecalc)**6)+(1.4379D-26*(par%timet-par%timecalc)**5)-(3.2241D-21*(par%timet-par%timecalc)**4)+(3.5801D-16*(par%timet-par%timecalc)**3)-(2.0842D-11*(par%timet-par%timecalc)**2)+(7.2871D-7*(par%timet-par%timecalc)) s%structdepth(i,j) = s%structdepth(i,j) + ynew - yold else if (hhmean(i,j)<=par%water_cut_off) then s%structdepth(i,j) = s%structdepth(i,j)+0 endif enddo enddo call csv_write2(s,par,ynew,yold) end subroutine HT_water_soil !Routine is not used. All current heat transfer components are available in subroutine heat_transfer subroutine csv_write2(s,par,ynew,yold) use params use interp use spaceparams use xmpi_module use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par real*8 :: ynew,yold logical :: advance, exist integer :: i,j inquire(file="Timet.txt", exist=exist) inquire(file="Timecalc.txt", exist=exist) inquire(file="ynew.txt", exist=exist) inquire(file="yold.txt", exist=exist) if (exist) then open(unit=50,FILE="Timet.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=51,FILE="Timecalc.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=52,FILE="ynew.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=53,FILE="yold.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") else open(unit=50,FILE="Timet.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=51,FILE="Timecalc.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=52,FILE="ynew.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=53,FILE="yold.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") end if write(50,20) par%timet write(51,20) par%timecalc write(52,21) ynew write(53,21) yold 20 format (f14.6) 21 format (f12.9) close(50) close(51) close(52) close(53) end subroutine csv_write2 subroutine csv_write(s,par,hws,Hwavemean,taubxmean,fwmean,Ubmean) use params use interp use spaceparams use xmpi_module use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par real*8,dimension(:,:),allocatable :: Hwavemean,taubxmean,fwmean,Ubmean,hws integer :: lun logical :: advance, exist integer :: i,j inquire(file="Heat_T_Coeff.txt", exist=exist) inquire(file="Wave_H_Mean", exist=exist) inquire(file="Bed_Shear_Mean.txt", exist=exist) inquire(file="Wave_Fric_Mean.txt", exist=exist) inquire(file="Orb_Vel_Mean.txt", exist=exist) if (exist) then open(unit=37,FILE="Heat_T_Coeff.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=38,FILE="Wave_H_Mean.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=39,FILE="Bed_Shear_Mean.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=40,FILE="Wave_Fric_Mean.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=41,FILE="Orb_Vel_Mean.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") else open(unit=37,FILE="Heat_T_Coeff.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=38,FILE="Wave_H_Mean.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=39,FILE="Bed_Shear_Mean.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=40,FILE="Wave_Fric_Mean.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=41,FILE="Orb_Vel_Mean.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") end if do i=1,s%nx+1 write(37,5) (hws(i,j),j=1,1) write(38,5) (Hwavemean(i,j),j=1,1) write(39,5) (taubxmean(i,j),j=1,1) write(40,5) (fwmean(i,j),j=1,1) write(41,5) (Ubmean(i,j),j=1,1) end do 5 format (f12.6) !write(37,*) "Next" close(37) close(38) close(39) close(40) close(41) end subroutine csv_write subroutine csv_read(s,par,hws) use params use interp use spaceparams use xmpi_module use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par real*8,dimension(:,:),allocatable :: hws integer :: lun logical :: advance, exist integer :: i,j inquire(file="Heat_T_Coeff.txt", exist=exist) if (exist) then open(unit=85,FILE="Heat_T_Coeff.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="asis", ACTION="read") else open(unit=85,FILE="Heat_T_Coeff.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="append", ACTION="write") end if do i=1,s%nx+1 read(85,110) (hws(i,j),j=1,1) end do 110 format (f16.6) !close(85) end subroutine csv_read subroutine csv_zs0read(zs0data) use params use interp use spaceparams use xmpi_module use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par real*8 :: zs0data integer :: lun logical :: advance, exist integer :: i,j ! inquire(file="zs0data.txt", exist=exist) if (exist) then open(unit=86,FILE="zs0data.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="asis", ACTION="read") else open(unit=86,FILE="zs0data.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="append", ACTION="write") end if read(86,127) zs0data 127 format (f16.6) !close(86) end subroutine csv_zs0read !Addition of Heat Transfer Routine Developed by Michelle and Guru !Heat transfer coefficients will be sent to this routine and thawdepth will be kicked out ! ! ! ! ! ! ! subroutine Coast_Thermal(s,par,hws,zbmean,vmagmean,glmmean,eulmean,thaw_depth_change_vert) !Developed by UAA !*****************************************************************************80 ! !! Coast_Thermal runs a thermal model of a coastal erosion problem. ! This routine is intended to interface with XBeach, where xbeach will provide the mechanical ! erosion modeling of a cross-shore beach profile, passing this model parameters associated ! with the coast line, this model will return a depth of thaw along the beach profile for the next ! erosion step in xbeach. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! 20 April 2016 by Michelle Wilber ! 26 January 2012 ! ! Original Author: ! ! original code was FD2D_HEAT_EXPLICIT_PRB - 26 January 2012 John Burkardt ! use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s !Xbeach structured arrays type(parameters) :: par !Xbeach structured arrays !external bc_test01 real ( kind = 8 ) cfl real ( kind = 8 ) dt !real ( kind = 8 ) Tw !integer ( kind = 4 ) nx !real ( kind = 8 ) hb(nx) !real*8,dimension(1,1) :: real*8 :: morf_adj,surfzone_Vol_total,v_factor,wet_length real*8 :: offshore_temp,heat_vel,adv_flux_total real*8 :: net_flux_total,rad_flux_total,twater_next,soil_flux_total real*8,save :: twater integer , dimension(:,:),allocatable :: temp_node_adj real*8,dimension(:,:),allocatable :: thaw_depth_change,thaw_depth,zbmean,vmagmean,glmmean,eulmean,thaw_depth_change_vert,hws real*8,dimension(:,:),allocatable,save :: h,hcond,hws_old !Enthalpy in array h is saved between each 30 minute session!!!!!!! real*8,dimension(:,:),allocatable :: h_new,h_newcond,h_test real*8,dimension(:,:),allocatable,save :: temp,thaw_depth_old,zb_old,tempcond,zb_initial !Temperature and thawdepth is saved in arrays temp and thaw_depth_old between each 30 minute session!!!!!!! real*8,dimension(:,:),allocatable :: temp_new,temp_newcond !Array reset between each 30 minute session real*8,dimension(:,:),allocatable :: tmat real*8,dimension(:),allocatable :: temp_top !Array reset between each 30 minute session real*8,dimension(:,:),allocatable :: temp_change,h_change,node_change,temp_change_cond,h_change_cond,c_heat_conv,c_heat_cond, cfl_conv, cfl_cond,node_adj,temp_final,thaw_tv,thaw_tc,thaw_hv,thaw_hc, delta_new real*8,dimension(1,1) :: surface_temp,surface_temp_storm real*8,dimension(1,1),save :: air_temp_storm real ( kind = 8 ) :: temp_f !freezeing temperature of water in soil integer ( kind = 4 ) :: i,i2,i3,do_min,do_max,i4 integer ( kind = 4 ) :: j,j2 real ( kind = 8 ), allocatable :: t(:) real ( kind = 8 ) :: t_min real ( kind = 8 ), allocatable :: y(:) real ( kind = 8 ) :: y_min ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FD1D_HEAT_EXPLICIT:' write ( *, '(a)' ) ' Compute an approximate solution to the time-dependent' write ( *, '(a)' ) ' one dimensional heat equation:' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' dH/dt - alpha * d2H/dx2 = f(x,t)' write ( *, '(a)' ) ' ' !write ( *, '(a)' ) ' Run a simple test case.' j=1 i2=0 i3=0 do_min=0 do_max=0 ! c = 2.0D+06 ! from the lit for frozen soil c, unfrozen more like 2.7D+06 ! heat capacity temp_f = 0.0D+00 ! freezing temp of soil water ! latent = 3.34D+08 * 2.0D-01 !latent heat of water by volume * volumetric water content J/m3 - assumes 20% water by volume ! morf_adj=0.0D+00 if(par%morfacopt==1) then ! Setting the morfac parameter correctly based on morfacopt setting !!!!!!! morf_adj=1.0D+00 else if(par%morfacopt==0) then morf_adj=1.0D+00 end if ! y_min = 0.0D+00 allocate (y(par%y_num)) ! Par%y_num is the number of equally spaced nodes to use between 0 and 1. do i2=1,par%y_num !Setting array y to zero at the start of the routine. y(i2)=0.0D+00 end do ! t_min = 0.0D+00 dt = (par%t_max - t_min ) / real (par%t_num - 1, kind = 8 ) allocate (t(par%t_num)) do j2=1,par%t_num !Setting array t to zero at the start of the routine. t(j2)=0.0D+00 end do call r8vec_linspace (s,par, par%t_num, t_min, par%t_max, t) !i3 is a pointer to the nodes used in Xbeach call r8vec_linspace (s,par, par%y_num, y_min, par%y_max, y) !i3 is a pointer to the nodes used in Xbeach ! ! Get the CFL coefficient. ! call fd1d_heat_explicit_cfl (s,par,t_min, y_min) if (par%timet > (par%therm_p/morf_adj)+5) then if(par%ht_calibration==1) then surface_temp(1,1) = 0.0D+00 call ht_calibration(s,par,air_temp_storm,surface_temp,temp_f,morf_adj) else if (par%ht_calibration==0) then surface_temp_storm(1,1) = 0.0D+00 call bluff_surface_temp(s,par,air_temp_storm,surface_temp_storm,temp_f,morf_adj) end if end if allocate (delta_new(s%nx+1,s%ny+1)) allocate (temp_final(par%y_num,s%nx+1)) !Temperature for thawdepth routine allocate (h_new(par%y_num,s%nx+1)) !Enthalpy next step allocate (h_test(par%y_num,s%nx+1)) !Enthalpy next step allocate (h_newcond(par%y_num,s%nx+1)) !Enthalpy next step allocate (temp_new(par%y_num,s%nx+1)) !Temperature next step allocate (temp_newcond(par%y_num,s%nx+1)) !Temperature next step allocate (temp_top(s%nx+1)) !Temperature in the top node in the vertical allocate (thaw_depth(s%nx+1,s%ny+1)) !Current thawdepth allocate (thaw_depth_change(s%nx+1,s%ny+1)) !Current thawdepth change in the vertical allocate (thaw_depth_change_vert(s%nx+1,s%ny+1)) !Current thawdepth change allocate (c_heat_conv(par%y_num,s%nx+1)) !Heat capacity for frozen and unfrozen soil in convective heat transfer routine allocate (c_heat_cond(par%y_num,s%nx+1)) !Heat capacity for frozen and unfrozen soil in conduction heat transfer routine allocate (cfl_conv(par%y_num,s%nx+1)) !Courant array for convection routine allocate (cfl_cond(par%y_num,s%nx+1)) !Courant array for conduction routine do i2=1,par%y_num do i3=1,s%nx+1 temp_final(i2,i3) = 0.0D+00 !Reset of arrays h_new(i2,i3) = 0.0D+00 !Reset of arrays h_newcond(i2,i3) = 0.0D+00 !Reset of arrays temp_new(i2,i3) = 0.0D+00 !Reset of arrays temp_newcond(i2,i3) = 0.0D+00 !Reset of arrays temp_top(i3) = 0.0D+00 !Reset of arrays c_heat_conv(i2,i3) = 0.0D+00 !Reset of arrays c_heat_cond(i2,i3) = 0.0D+00 !Reset of arrays cfl_conv(i2,i3) = 0.0D+00 !Reset of arrays cfl_cond(i2,i3) = 0.0D+00 !Reset of arrays thaw_depth(i3,j) = 0.0D+00 !Current thawdepth array thaw_depth_change(i3,j) = 0.0D+00 !Current thawdepth array thaw_depth_change_vert(i3,j) = 0.0D+00 !Current thawdepth array delta_new(i3,j) = 0.0D+00 end do end do ! allocate (node_change(s%nx+1,s%ny+1)) allocate (node_adj(s%nx+1,s%ny+1)) allocate (temp_node_adj(s%nx+1,s%ny+1)) do j=1,s%ny+1 do i3=1,s%nx+1 node_change(i3,j)=0.0D+00 node_adj(i3,j)=0.0D+00 temp_node_adj(i3,j)=0.0D+00 end do end do allocate (temp_change(par%y_num,s%nx+1)) allocate (temp_change_cond(par%y_num,s%nx+1)) allocate (h_change(par%y_num,s%nx+1)) allocate (h_change_cond(par%y_num,s%nx+1)) allocate (thaw_tv(par%y_num,s%nx+1)) allocate (thaw_tc(par%y_num,s%nx+1)) allocate (thaw_hv(par%y_num,s%nx+1)) allocate (thaw_hc(par%y_num,s%nx+1)) do i2=1,par%y_num do i3=1,s%nx+1 temp_change(i2,i3)=0.0D+00 temp_change_cond(i2,i3)=0.0D+00 h_change(i2,i3)=0.0D+00 h_change_cond(i2,i3)=0.0D+00 thaw_tv(i2,i3)=0.0D+00 thaw_tc(i2,i3)=0.0D+00 thaw_hv(i2,i3)=0.0D+00 thaw_hc(i2,i3)=0.0D+00 end do end do ! if (par%timet < (par%therm_p/morf_adj)+5) then !Allocation of arrays only occurs during the 1st 30 minute session as data needs to be preserved throughout simulation!!! allocate (h(par%y_num,s%nx+1)) !Enthalpy for convection allocate (hcond(par%y_num,s%nx+1)) !Enthalpy for conduction allocate (temp(par%y_num,s%nx+1)) !Temperature for convection allocate (tempcond(par%y_num,s%nx+1)) !Temperature for conduction do i2=1,par%y_num do i3=1,s%nx+1 h(i2,i3) = 0.0D+00 !Reset of arrays temp(i2,i3) = 0.0D+00 !Reset of arrays hcond(i2,i3) = 0.0D+00 !Reset of arrays tempcond(i2,i3) = 0.0D+00 !Reset of arrays end do end do endif ! if (par%timet < (par%therm_p/morf_adj)+5) then allocate (zb_old(s%nx+1,s%ny+1)) allocate (zb_initial(s%nx+1,s%ny+1)) allocate (hws_old(s%nx+1,s%ny+1)) allocate (thaw_depth_old(s%nx+1,s%ny+1)) !Thawdepth in previous step do j=1,s%ny+1 do i3=1,s%nx+1 hws_old(i3,j)=0.0D+00 hws_old(i3,j)=hws(i3,j) ! Old convective heat coefficient is updated so model can keep track of changes to wet cells. thaw_depth_old(i3,j) = 0.0D+00 !thawdepth array from previous 30 minute period zb_old(i3,j) = 0.0D+00 zb_initial(i3,j) = 0.0D+00 !zb_old(i3,j)=s%zb(i3,j) !Allocation and assignment of bed elevation from the 30 minute period. open(unit=200,FILE="thaw_depth_old_input.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="asis", ACTION="read") open(unit=201,FILE="bed.dep", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="asis", ACTION="read") open(unit=202,FILE="bed.dep", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="asis", ACTION="read") read(200,126) thaw_depth_old(i3,j) read(201,126) zb_old(i3,j) read(202,126) zb_initial(i3,j) 126 format (f10.5) end do end do end if ! call read_water_input(s,par,offshore_temp,rad_flux_total,hws,wet_length) !!!CALLING ON INPUT VALUES FOR WATER BUDGET i3=0 if(par%ht_calibration==1) then do_min=310 do_max=310 else if (par%ht_calibration==0) then do_min=par%start_node do_max=par%end_node end if !!!!!!!!!!!!!!!!!!!!!!!! if (par%timet < (par%therm_p/morf_adj)+5) then do i3=1,do_max if(par%calm_condition==1) then call ic_test01_cond (s,par,tempcond,i3) else call ic_test01_conv (s,par,temp,i3) !Initialization for NON-CALM CONDITIONS call ic_test01_cond (s,par,tempcond,i3) end if end do end if !!!!!!!!!!!!!!!!!!!!!!! do i3=do_min,do_max !Do loop for all the nodes in the cross-shore direction! ! if (par%timet < (par%therm_p/morf_adj)+5) then !Initial conditions and boundary conditions only run during 1st 30 minute session ! if(par%calm_condition==1) then !Initialization for CALM CONDITIONS !call ic_test01_cond (s,par,tempcond,i3) call bc_test01conduction (s,par,tempcond,h_newcond,c_heat_cond,temp_f,i3) ! do i2=1,par%y_num if(tempcond(i2,i3) >=temp_f) then cfl_cond(i2,i3)=((par%ht_ku/par%c_unfrozen_soil)*dt)/(par%y_max/(par%y_num-1))**2 c_heat_cond(i2,i3)=par%c_unfrozen_soil else cfl_cond(i2,i3)=((par%ht_kf/par%c_frozen_soil)*dt)/(par%y_max/(par%y_num-1))**2 c_heat_cond(i2,i3)=par%c_frozen_soil end if end do do i2=1,par%y_num if (tempcond(i2,i3) <= temp_f) then hcond(i2,i3) = tempcond(i2,i3) * c_heat_cond(i2,i3) !i2 is a pointer in the vertical and i3 in the horizontal else hcond(i2,i3) = tempcond(i2,i3) * c_heat_cond(i2,i3) + par%Li*par%nb !i2 is a pointer in the vertical and i3 in the horizontal endif end do !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! else !call ic_test01_conv (s,par,temp,i3) !Initialization for NON-CALM CONDITIONS !call ic_test01_cond (s,par,tempcond,i3) do i2=1,par%y_num if(temp(i2,i3) >=temp_f) then cfl_conv(i2,i3)=((par%ht_ku/par%c_unfrozen_soil)*dt)/(par%y_max/(par%y_num-1))**2 !!!!CFL AND HEAT CAPACITY ASSIGNMENTS c_heat_conv(i2,i3)=par%c_unfrozen_soil !CONVECTIVE BOUNDARY else cfl_conv(i2,i3)=((par%ht_kf/par%c_frozen_soil)*dt)/(par%y_max/(par%y_num-1))**2 c_heat_conv(i2,i3)=par%c_frozen_soil end if end do do i2=1,par%y_num if(tempcond(i2,i3) >=temp_f) then cfl_cond(i2,i3)=((par%ht_ku/par%c_unfrozen_soil)*dt)/(par%y_max/(par%y_num-1))**2 !!!!CFL AND HEAT CAPACITY ASSIGNMENTS c_heat_cond(i2,i3)=par%c_unfrozen_soil !CONDUCTIVE BOUNDARY else cfl_cond(i2,i3)=((par%ht_kf/par%c_frozen_soil)*dt)/(par%y_max/(par%y_num-1))**2 c_heat_cond(i2,i3)=par%c_frozen_soil end if end do call bc_test01(s,par,hws, cfl_conv, dt, y, t, temp,temp_top,h_new,c_heat_conv,temp_f,i3) !INITIALIZATION OF BOUNDARY CONDITIONS for NON-CALM CONDITIONS FOR 1ST call bc_test01conduction(s,par,tempcond,h_newcond,c_heat_cond,temp_f,i3) do i2=1,par%y_num if (temp(i2,i3) <= temp_f) then h(i2,i3) = temp(i2,i3) * c_heat_conv(i2,i3) !i2 is a pointer in the vertical and i3 in the horizontal else h(i2,i3) = temp(i2,i3) * c_heat_conv(i2,i3) + par%Li*par%nb !i2 is a pointer in the vertical and i3 in the horizontal endif if (tempcond(i2,i3) <= temp_f) then hcond(i2,i3) = tempcond(i2,i3) * c_heat_cond(i2,i3) !i2 is a pointer in the vertical and i3 in the horizontal else hcond(i2,i3) = tempcond(i2,i3) * c_heat_cond(i2,i3) + par%Li*par%nb !i2 is a pointer in the vertical and i3 in the horizontal endif end do endif end if if(par%calm_condition==0) then if(par%ht_bed_updating==1) then call erosion_sedimentation(s,par,zbmean,temp,tempcond,temp_change,temp_change_cond,h,hcond,h_change,h_change_cond,zb_old,node_change,i3,delta_new,thaw_depth_old,morf_adj,thaw_depth) !!!BED UPDATING FOR HEAT TRANSFER ROUTINE end if j=1 if(hws_old(i3,j) == 0.0D+00 .and. hws(i3,j) > 0.0D+00) then do i2=1,par%y_num h(i2,i3) = h(i2,i3)*0 temp(i2,i3) = temp(i2,i3)*0 h(i2,i3) = hcond(i2,i3) temp(i2,i3) = tempcond(i2,i3) !!!WETTING AND DRYING TRANSFER BETWEEN TEMPERATURE ARRAYS end do else if(hws_old(i3,j) > 0.0D+00 .and. hws(i3,j) == 0.0D+00) then do i2=1,par%y_num tempcond(i2,i3) = tempcond(i2,i3)*0 hcond(i2,i3) = hcond(i2,i3)*0 hcond(i2,i3) = h(i2,i3) tempcond(i2,i3) = temp(i2,i3) !!!WETTING AND DRYING TRANSFER BETWEEN TEMPERATURE ARRAYS end do end if else end if ! do j2 = 2, par%t_num if(par%calm_condition==1) then call fd1d_heat_explicit_conduction (s,par,y,t,dt,cfl_cond,temp_f,tempcond,temp_newcond,hcond,h_newcond,c_heat_cond,surface_temp,surface_temp_storm,i3,morf_adj) do i2=1,par%y_num hcond(i2,i3) = h_newcond(i2,i3) !Enthalpy array updated with new values tempcond(i2,i3) = temp_newcond(i2,i3) !Temperature array updated with new values end do else if (hws(i3,j) > 0.0D+00) then call fd1d_heat_explicit (s,par,hws,y,t,dt,cfl_conv,temp_f,temp_top,temp,temp_new,h,h_new,c_heat_conv,twater,h_test,i3,morf_adj) do i2=1,par%y_num h(i2,i3) = h_new(i2,i3) !Enthalpy array updated with new values temp(i2,i3) = temp_new(i2,i3) !Temperature array updated with new values end do elseif (hws(i3,j) == 0.0D+00) then call fd1d_heat_explicit_conduction (s,par,y,t,dt,cfl_cond,temp_f,tempcond,temp_newcond,hcond,h_newcond,c_heat_cond,surface_temp,surface_temp_storm,i3,morf_adj) do i2=1,par%y_num hcond(i2,i3) = h_newcond(i2,i3) !Enthalpy array updated with new values tempcond(i2,i3) = temp_newcond(i2,i3) !Temperature array updated with new values end do end if end if end do !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!HEAT TRANSFER ROUTINE IS RUN!! ! if(par%calm_condition==1) then ! Start of wetting and drying if statements!!!!! do i2=1,par%y_num temp_final(i2,i3)=tempcond(i2,i3) end do else j=1 if(hws_old(i3,j)==0.0D+00 .and. hws(i3,j)==0.0D+00) then do i2=1,par%y_num temp_final(i2,i3)=tempcond(i2,i3) end do else if(hws_old(i3,j) > 0.0D+00 .and. hws(i3,j) > 0.0D+00) then do i2=1,par%y_num temp_final(i2,i3)=temp(i2,i3) end do else if(hws_old(i3,j) == 0.0D+00 .and. hws(i3,j) > 0.0D+00) then do i2=1,par%y_num temp_final(i2,i3) = temp(i2,i3) end do else if(hws_old(i3,j) > 0.0D+00 .and. hws(i3,j) == 0.0D+00) then do i2=1,par%y_num temp_final(i2,i3) = tempcond(i2,i3) end do end if end if ! End of wetting and drying if statements!!!!! ! call find_thaw_depth(s,par,zbmean,temp_f,y,thaw_depth,i3,thaw_depth_change,thaw_depth_old,zb_old,node_adj,temp_final,delta_new,thaw_tv,thaw_tc,thaw_hv,thaw_hc,temp,tempcond,h,hcond,temp_node_adj,morf_adj,thaw_depth_change_vert) !Change in thawdepth after 30 minute session sent to HT_water_bluff subroutine to update structdepth array!!! ! j=1 zb_old(i3,j)=zb_old(i3,j)*0 ! Setting zb_old equal to zero zb_old(i3,j)=s%zb(i3,j) hws_old(i3,j)=hws_old(i3,j)*0 ! Updating zb_old with current bed elevation hws_old(i3,j)=hws(i3,j) end do !End of the do loop in the cross-shore direction! ! j=1 i4 = par%start_node do i3=1,i4-1 do i2=1,par%y_num !!!!!!!!!!!!!! Manual correction for when certain nodes are diverging !!!!!!!!!!!!!!! temp(i2,i3) = temp(i2,i3)*0 temp_final(i2,i3) = temp_final(i2,i3)*0 tempcond(i2,i3) = tempcond(i2,i3)*0 h(i2,i3) = h(i2,i3)*0 hcond(i2,i3) = hcond(i2,i3)*0 end do end do do i3=1,i4-1 do i2=1,par%y_num temp(i2,i3) = temp(i2,i4) tempcond(i2,i3) = tempcond(i2,i4) h(i2,i3) = h(i2,i4) hcond(i2,i3) = hcond(i2,i4) temp_final(i2,i3) = temp_final(i2,i4) end do end do do i3=1,i4-1 thaw_depth(i3,1) = thaw_depth(i3,1)*0 thaw_depth_old(i3,1) = thaw_depth_old(i3,1)*0 thaw_depth_change(i3,1) = thaw_depth_change(i3,1)*0 s%structdepth(i3,1) = s%structdepth(i3,1)*0 thaw_depth_change_vert(i3,1) = thaw_depth_change_vert(i3,1)*0 end do do i3=1,i4-1 thaw_depth(i3,1) = thaw_depth(i4,1) thaw_depth_old(i3,1) = thaw_depth_old(i4,1) thaw_depth_change(i3,1) = thaw_depth_change(i4,1) s%structdepth(i3,1) = s%structdepth(i4,1) thaw_depth_change_vert(i3,1) = thaw_depth_change_vert(i4,1) end do ! !!!!!!!!!!!!!! End of manual correction for when certain nodes are diverging !!!!!!!!!!!!!!! call water_heat_budget(s,par,temp_final,net_flux_total,hws,rad_flux_total,twater_next,twater,offshore_temp,soil_flux_total,vmagmean,glmmean,eulmean,heat_vel,adv_flux_total,morf_adj,surfzone_Vol_total,v_factor) !!!!COMPUTING A NEW WATER TEMPEREATURE ! call r8vec_write(s,par,thaw_depth,i3,thaw_depth_change,node_change,temp,h,hcond,tempcond,zbmean,net_flux_total,twater_next,rad_flux_total,vmagmean,glmmean,eulmean,soil_flux_total,heat_vel,adv_flux_total,zb_old,thaw_depth_old,delta_new,surfzone_Vol_total,wet_length,temp_f,h_test,surface_temp) !Subroutine now also creates file with change in thaw depth. ! ! if(par%timet >= (par%tstop - (par%therm_p/morf_adj))) then ! if (par%bed_adjustment==1) then call correction_bed_change(s,par,temp,tempcond,h,hcond,thaw_depth,delta_new,zb_initial,hws) !Routine to address the buildup of sediments in the nearshore zone else if (par%bed_adjustment==0) then endif call correction_thaw_struct(s,par,temp,tempcond,h,hcond,thaw_depth,delta_new) !Address any differences between thawdepth and structdepth call temp_enthrope_write(s,par,temp,tempcond,h,hcond,thaw_depth) end if deallocate(temp_top) deallocate(temp_new) deallocate(temp_newcond) deallocate(h_new) deallocate(h_newcond) deallocate(t) deallocate(y) deallocate(thaw_depth) deallocate(thaw_depth_change) deallocate(temp_final) deallocate(temp_change) deallocate(temp_change_cond) deallocate(h_change) deallocate(h_change_cond) deallocate(node_change) deallocate(thaw_tv) deallocate(thaw_tc) deallocate(thaw_hv) deallocate(thaw_hc) deallocate(c_heat_conv) deallocate(c_heat_cond) deallocate(cfl_conv) deallocate(cfl_cond) deallocate(delta_new) ! end subroutine Coast_Thermal subroutine correction_bed_change(s,par,temp,tempcond,h,hcond,thaw_depth,delta_new,zb_initial,hws) use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par real*8,dimension(:,:),allocatable :: temp,tempcond,h,hcond,corr,temp_change,temp_change_cond,h_change,h_change_cond,thaw_depth,delta_new,zb_initial,corr_bed,hws integer :: i,j,i2 ! allocate (corr_bed(s%nx+1,s%ny+1)) allocate (temp_change(par%y_num,s%nx+1)) allocate (temp_change_cond(par%y_num,s%nx+1)) allocate (h_change(par%y_num,s%nx+1)) allocate (h_change_cond(par%y_num,s%nx+1)) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do j=1,s%ny+1 do i=1,s%nx+1 if (hws(i,j) > 0) then corr_bed(i,j) = 0.0D+00 corr_bed(i,j) = s%zb(i,j) - zb_initial(i,j) if (corr_bed(i,j) < 0) then corr_bed(i,j) = 0.0D+00 endif corr_bed(i,j) = nint(corr_bed(i,j) / (par%y_max/(par%y_num-1))) s%structdepth(i,j) = s%structdepth(i,j) - (corr_bed(i,j) * (par%y_max/(par%y_num-1))) thaw_depth(i,j) = thaw_depth(i,j) - (corr_bed(i,j) * (par%y_max/(par%y_num-1)) * COSD(delta_new(i,j))) s%zb(i,j) = s%zb(i,j) - (corr_bed(i,j) * (par%y_max/(par%y_num-1))) end if end do end do ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! corr_bed(i,j) = abs(corr_bed(i,j)) do i2=1,par%y_num-(corr_bed(i,j)+1) temp_change(i2,i) = temp(i2+corr_bed(i,j),i) temp_change_cond(i2,i) = tempcond(i2+corr_bed(i,j),i) h_change(i2,i) = h(i2+corr_bed(i,j),i) h_change_cond(i2,i) = hcond(i2+corr_bed(i,j),i) end do do i2=par%y_num-corr_bed(i,j),par%y_num temp_change(i2,i) = temp(par%y_num,i) temp_change_cond(i2,i) = tempcond(par%y_num,i) h_change(i2,i) = h(i2-1,i) h_change_cond(i2,i) = hcond(i2-1,i) end do !temp_change(1,i3)=temp(1,i3) h_change(1,i)=h(1,i) h_change_cond(1,i)=hcond(1,i) h_change(par%y_num,i)=h(par%y_num,i) h_change_cond(par%y_num,i)=hcond(par%y_num,i) do i2=1,par%y_num temp(i2,i) = temp(i2,i)*0 ! Setting temp and h arrays equal to zero tempcond(i2,i) = tempcond(i2,i)*0 h(i2,i) = h(i2,i)*0 hcond(i2,i) = hcond(i2,i)*0 end do do i2=1,par%y_num temp(i2,i) = temp_change(i2,i) ! Updating temp and h arrays with new temperature profile tempcond(i2,i) = temp_change_cond(i2,i) h(i2,i) = h_change(i2,i) hcond(i2,i) = h_change_cond(i2,i) end do end subroutine correction_bed_change ! ! subroutine correction_thaw_struct(s,par,temp,tempcond,h,hcond,thaw_depth,delta_new) use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par real*8,dimension(:,:),allocatable :: temp,tempcond,h,hcond,corr,temp_change,temp_change_cond,h_change,h_change_cond,thaw_depth,delta_new integer :: i,j,i2 allocate (corr(s%nx+1,s%ny+1)) allocate (temp_change(par%y_num,s%nx+1)) allocate (temp_change_cond(par%y_num,s%nx+1)) allocate (h_change(par%y_num,s%nx+1)) allocate (h_change_cond(par%y_num,s%nx+1)) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do j=1,s%ny+1 do i=1,s%nx+1 corr(i,j) = 0.0D+00 corr(i,j) = (s%structdepth(i,j) * COSD(delta_new(i,j))) - thaw_depth(i,j) corr(i,j) = nint(corr(i,j) / (par%y_max/(par%y_num-1))) end do end do !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do j=1,s%ny+1 do i=1,s%nx+1 if (corr(i,j) > 0 ) then !Need to add positive nodes to thawdepth do i2=1,corr(i,j)+1 temp_change(i2,i) = temp(2,i) temp_change_cond(i2,i) = tempcond(2,i) h_change(i2,i) = h(2,i) h_change_cond(i2,i) = hcond(2,i) end do do i2=corr(i,j)+2,par%y_num temp_change(i2,i) = temp(i2-corr(i,j),i) temp_change_cond(i2,i) = tempcond(i2-corr(i,j),i) h_change(i2,i) = h(i2-corr(i,j),i) h_change_cond(i2,i) = hcond(i2-corr(i,j),i) end do temp_change(par%y_num,i) = temp(par%y_num,i) temp_change_cond(par%y_num,i) = tempcond(par%y_num,i) h_change(par%y_num,i) = h(par%y_num,i) h_change_cond(par%y_num,i) = hcond(par%y_num,i) h_change(1,i) = h(1,i) h_change_cond(1,i) = hcond(1,i) do i2=1,par%y_num temp(i2,i) = temp(i2,i)*0 ! Setting temp and h arrays equal to zero tempcond(i2,i) = tempcond(i2,i)*0 h(i2,i) = h(i2,i)*0 hcond(i2,i) = hcond(i2,i)*0 end do do i2=1,par%y_num temp(i2,i) = temp_change(i2,i) ! Updating temp and h arrays with new temperature profile tempcond(i2,i) = temp_change_cond(i2,i) h(i2,i) = h_change(i2,i) hcond(i2,i) = h_change_cond(i2,i) end do thaw_depth(i,j) = thaw_depth(i,j) + (corr(i,j)*(par%y_max/(par%y_num-1))) ! It's like a correction for sedimentation ! else if (corr(i,j) < 0 ) then !Need to reduce the number of positive nodes corr(i,j) = abs(corr(i,j)) do i2=1,par%y_num-(corr(i,j)+1) temp_change(i2,i) = temp(i2+corr(i,j),i) temp_change_cond(i2,i) = tempcond(i2+corr(i,j),i) h_change(i2,i) = h(i2+corr(i,j),i) h_change_cond(i2,i) = hcond(i2+corr(i,j),i) end do do i2=par%y_num-corr(i,j),par%y_num temp_change(i2,i) = temp(par%y_num,i) temp_change_cond(i2,i) = tempcond(par%y_num,i) h_change(i2,i) = h(i2-1,i) h_change_cond(i2,i) = hcond(i2-1,i) end do !temp_change(1,i3)=temp(1,i3) h_change(1,i)=h(1,i) h_change_cond(1,i)=hcond(1,i) h_change(par%y_num,i)=h(par%y_num,i) h_change_cond(par%y_num,i)=hcond(par%y_num,i) do i2=1,par%y_num temp(i2,i) = temp(i2,i)*0 ! Setting temp and h arrays equal to zero tempcond(i2,i) = tempcond(i2,i)*0 h(i2,i) = h(i2,i)*0 hcond(i2,i) = hcond(i2,i)*0 end do do i2=1,par%y_num temp(i2,i) = temp_change(i2,i) ! Updating temp and h arrays with new temperature profile tempcond(i2,i) = temp_change_cond(i2,i) h(i2,i) = h_change(i2,i) hcond(i2,i) = h_change_cond(i2,i) end do thaw_depth(i,j) = thaw_depth(i,j) - (corr(i,j)*(par%y_max/(par%y_num-1))) ! It's like a correction for erosion else end if end do end do end subroutine correction_thaw_struct !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! WATER HEAT BUDGET subroutine read_water_input(s,par,offshore_temp,rad_flux_total,hws,wet_length) use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par ! real*8 :: offshore_temp real*8 :: wet_length,rad_flux_total real*8 :: NCAR_Total logical :: exist integer :: i,j,i2 real*8,dimension(:,:),allocatable :: hws ! Calculation of wet length wet_length = 0 do j=1,s%ny+1 do i=1,s%nx if (hws(i,j) > 0.0) then wet_length = wet_length + (s%x(i+1,j)-s%x(i,j)) else end if end do end do offshore_temp=0.0D+00 rad_flux_total=0.0D+00 !NCAR_S=0.0D+00 !NCAR_L=0.0D+00 !NCAR_LA=0.0D+00 !NCAR_SE=0.0D+00 NCAR_Total=0.0D+00 !!!!!!!!!!!!!!!!!!!!!!!!!!NCAR!!!!!!!!!!!!!!!!!!!!!!!!!!!! !inquire(file="Short_Wave_Flux.txt", exist=exist) !inquire(file="Long_Wave_Flux.txt", exist=exist) !inquire(file="Latent_Heat_Flux.txt", exist=exist) !inquire(file="Sensitive_Heat_Flux.txt", exist=exist) inquire(file="Flux_Total.txt", exist=exist) if (exist) then !open(unit=100,FILE="Short_Wave_Flux.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="asis", ACTION="read") !open(unit=101,FILE="Long_Wave_Flux.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="asis", ACTION="read") !open(unit=102,FILE="Latent_Heat_Flux.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="asis", ACTION="read") !open(unit=103,FILE="Sensitive_Heat_Flux.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="asis", ACTION="read") open(unit=104,FILE="Flux_Total.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="asis", ACTION="read") else !open(unit=100,FILE="Short_Wave_Flux.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="append", ACTION="write") !open(unit=101,FILE="Long_Wave_Flux.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="append", ACTION="write") !open(unit=102,FILE="Latent_Heat_Flux.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="append", ACTION="write") !open(unit=103,FILE="Sensitive_Heat_Flux.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=104,FILE="Flux_Total.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="append", ACTION="write") end if ! read(100,120) (NCAR_S(i,j),j=1,1) ! read(101,120) (NCAR_L(i,j),j=1,1) ! read(102,120) (NCAR_LA(i,j),j=1,1) ! read(103,120) (NCAR_SE(i,j),j=1,1) !read(100,120) NCAR_S !read(101,120) NCAR_L !read(102,120) NCAR_LA !read(103,120) NCAR_SE read(104,120) NCAR_Total 120 format (f16.6) !rad_flux_total = (NCAR_S+NCAR_L+NCAR_LA+NCAR_SE)*wet_length rad_flux_total = NCAR_Total*wet_length !!!!!!!!!!!!!!!!!!!!!!!!!!NCAR!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! inquire(file="Offshore_Daily_Temp.txt", exist=exist) if (exist) then open(unit=105,FILE="Offshore_Daily_Temp.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="asis", ACTION="read") else open(unit=105,FILE="Offshore_Daily_Temp.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="append", ACTION="write") end if read(105,121) offshore_temp 121 format (f16.6) end subroutine read_water_input subroutine water_heat_budget(s,par,temp_final,net_flux_total,hws,rad_flux_total,twater_next,twater,offshore_temp,soil_flux_total,vmagmean,glmmean,eulmean,heat_vel,adv_flux_total,morf_adj,surfzone_Vol_total,v_factor) use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par real*8 :: morf_adj,ratio real*8 :: offshore_temp,heat_vel,v_factor real*8 :: soil_flux_total,twater,twater_next,rad_flux_total,adv_flux_total,net_flux_total,hnet,surfzone_Vol_total,offshore_depth real*8,dimension(:,:),allocatable :: soil_flux_rate,temp_final,vmagmean,glmmean,eulmean,hws,surfzone_Vol logical :: exist integer :: i,j,i2 j=1 v_factor = 0 offshore_depth = 0 do j=1,s%ny+1 do i=2,5 offshore_depth = offshore_depth + (s%zs(i,j)-s%zb(i,j)) !Approximate depth at offshore boundary end do end do v_factor = offshore_depth/4 ! Four nodes used to determine average depth allocate(soil_flux_rate(s%nx+1,s%ny+1)) allocate(surfzone_Vol(s%nx+1,s%ny+1)) soil_flux_total=0.0D+00 surfzone_Vol_total=0.0D+00 adv_flux_total=0.0D+00 net_flux_total=0.0D+00 twater_next=0.0D+00 hnet=0.0D+00 ! Calculate volume of nearshore zone do j=1,s%ny+1 do i=1,s%nx if (s%zs(i,j)-s%zb(i,j) > 0) then surfzone_Vol(i,j) = (s%zs(i,j)-s%zb(i,j))*(s%x(i+1,j)-s%x(i,j)) surfzone_Vol_total = surfzone_Vol_total + surfzone_Vol(i,j) else end if end do end do !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do j=1,s%ny+1 do i=1,s%nx+1 soil_flux_rate(i,j)=0.0D+00 end do end do if (par%timet < (par%therm_p/morf_adj)+5) then twater=par%Tw do j=1,s%ny+1 do i=1,s%nx soil_flux_rate(i,j) = hws(i,j) * (temp_final(1,i)-twater) end do end do else if (par%timet > (par%therm_p/morf_adj)+5) then ! twater=twater_next do j=1,s%ny+1 do i=1,s%nx soil_flux_rate(i,j) = hws(i,j) * (temp_final(1,i)-twater) end do end do end if do j=1,s%ny+1 do i=1,s%nx soil_flux_total = soil_flux_total + (((soil_flux_rate(i,j) + soil_flux_rate(i+1,j))/2) * (s%x(i+1,j)-s%x(i,j))) end do end do heat_vel=(sum(glmmean(1:16,1))/16) adv_flux_total = v_factor*par%qs*(par%Cw/1000.0D+00)*(heat_vel)*(offshore_temp-twater) ratio = adv_flux_total/(soil_flux_total + rad_flux_total) if (ratio > 4) then ! 4 is an arbitrary value used ratio = 4 elseif (ratio < -4) then ratio = -4 else endif if (ratio == 4 .or. ratio == -4) then adv_flux_total = ratio * (soil_flux_total + rad_flux_total) !Correction to avoid excessively high advective flux that breaks the model! else endif net_flux_total = adv_flux_total + rad_flux_total + soil_flux_total hnet=(par%qs*(par%Cw/1000.0D+00)*surfzone_Vol_total) !(par%net_heat_coeff) twater_next=(net_flux_total*par%therm_p + hnet*twater)/hnet ! Should be the correct calculation twater=twater*0 twater=twater_next deallocate(soil_flux_rate) deallocate(surfzone_Vol) end subroutine water_heat_budget !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!WATER HEAT BUDGET subroutine ht_calibration(s,par,air_temp_storm,surface_temp,temp_f,morf_adj) use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par real ( kind = 8 ) :: temp_f real*8,dimension(1,1) :: surface_temp,air_temp_storm logical :: exist integer :: i,j,i2,i3 real*8 :: morf_adj if ((par%timet > (par%therm_p/morf_adj)) .and. (par%timet < 3*(par%therm_p/morf_adj))) then !1 air_temp_storm(1,1)=0.0D+00 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if if ((par%timet > 86400) .and. (par%timet < 87000)) then !2 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if if ((par%timet > 172800) .and. (par%timet < 173400)) then !3 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if if ((par%timet > 259200) .and. (par%timet < 259800)) then !4 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if if ((par%timet > 345600) .and. (par%timet < 346200)) then !5 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if if ((par%timet > 432000) .and. (par%timet < 432600)) then !6 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if if ((par%timet > 518400) .and. (par%timet < 519000)) then !7 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if if ((par%timet > 604800) .and. (par%timet < 605400)) then !8 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if if ((par%timet > 691200) .and. (par%timet < 691800)) then !9 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if if ((par%timet > 777600) .and. (par%timet < 778200)) then !10 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if if ((par%timet > 864000) .and. (par%timet < 864600)) then !11 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if if ((par%timet > 950400) .and. (par%timet < 951000)) then !12 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if ! if (air_temp_storm(1,1) < temp_f) then surface_temp(1,1) = par%n_air_surf * air_temp_storm(1,1) else if (air_temp_storm(1,1) >= temp_f) then surface_temp(1,1) = par%n_air_surf * air_temp_storm(1,1) end if 116 format (f6.3) end subroutine ht_calibration subroutine bluff_surface_temp(s,par,air_temp_storm,surface_temp_storm,temp_f,morf_adj) use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par real ( kind = 8 ) :: temp_f real*8,dimension(1,1) :: air_temp_storm real*8,dimension(1,1) :: surface_temp_storm logical :: exist integer :: i,j,i2,i3 real*8 :: morf_adj ! inquire(file="Surface_Temp_Calibration.txt", exist=exist) if ((par%timet > (par%therm_p/morf_adj)) .and. (par%timet < 3*(par%therm_p/morf_adj))) then !1 air_temp_storm(1,1)=0.0D+00 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if if ((par%timet > 86400) .and. (par%timet < 87000)) then !2 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if if ((par%timet > 172800) .and. (par%timet < 173400)) then !3 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if if ((par%timet > 259200) .and. (par%timet < 259800)) then !4 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if if ((par%timet > 345600) .and. (par%timet < 346200)) then !5 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if if ((par%timet > 432000) .and. (par%timet < 432600)) then !6 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if if ((par%timet > 518400) .and. (par%timet < 519000)) then !7 open(unit=73,FILE="air_temp_storm.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="asis", ACTION="read") do i=1,1 read(73,116) (air_temp_storm(i,j),j=1,1) end do end if if (air_temp_storm(1,1) < temp_f) then surface_temp_storm(1,1) = air_temp_storm(1,1) elseif (air_temp_storm(1,1) >= temp_f) then surface_temp_storm(1,1) = par%n_air_surf * air_temp_storm(1,1) end if 116 format (f6.3) end subroutine bluff_surface_temp ! ! subroutine erosion_sedimentation(s,par,zbmean,temp,tempcond,temp_change,temp_change_cond,h,hcond,h_change,h_change_cond,zb_old,node_change,i3,delta_new,thaw_depth_old,morf_adj,thaw_depth) use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par real ( kind = 8 ),allocatable :: y(:) real ( kind = 8 ),allocatable :: t(:) real ( kind = 8 ),allocatable :: temp_change(:,:),temp(:,:),node_change(:,:),temp_change_cond(:,:),tempcond(:,:),zbmean(:,:),thaw_depth_old(:,:),thaw_depth(:,:),delta_new(:,:) real ( kind = 8 ),allocatable :: h_change(:,:),h(:,:),zb_old(:,:),h_change_cond(:,:),hcond(:,:) real ( kind = 8 ) :: erosion_factor,sedimentation_factor,morf_adj integer :: i,j,i2,i3 ! j=1 if(i3==1) then delta_new(i3,j)=ATAND(((s%zb(i3+1,j))-s%zb(i3,j))/(s%x(i3+1,j)-s%x(i3,j))) else if(i3==s%nx+1) then delta_new(i3,j)=ATAND(((s%zb(i3,j))-s%zb(i3-1,j))/(s%x(i3,j)-s%x(i3-1,j))) else delta_new(i3,j)=ATAND(((s%zb(i3+1,j))-s%zb(i3-1,j))/(s%x(i3+1,j)-s%x(i3-1,j))) endif if(delta_new(i3,j) >= par%slope_threshold) then erosion_factor=par%E_factor_high sedimentation_factor=par%S_factor_high else erosion_factor=par%E_factor_low sedimentation_factor=par%S_factor_low endif if (s%zb(i3,j) - zb_old(i3,j) < 0.0D+00) then node_change(i3,j)=abs(((s%zb(i3,j) - zb_old(i3,j))/(par%y_max/(par%y_num-1))) * COSD(delta_new(i3,j)) * erosion_factor) !Erosion factor to be calibrated for each model ! if((node_change(i3,j)*(par%y_max/(par%y_num-1))) > thaw_depth_old(i3,j)) then !If morfac=1 then erosion factor should be 1 node_change(i3,j) = floor(thaw_depth_old(i3,j)/(par%y_max/(par%y_num-1))) else if(par%morfac >=2.0D+00) then node_change(i3,j) = nint(node_change(i3,j)) else node_change(i3,j) = nint(node_change(i3,j)) endif endif ! else if (s%zb(i3,j) - zb_old(i3,j) > 0.0D+00) then node_change(i3,j)=abs(((s%zb(i3,j) - zb_old(i3,j))/(par%y_max/(par%y_num-1))) * COSD(delta_new(i3,j)) *sedimentation_factor) !Sedimentation factor to be calibrated for each model if(par%morfac >=2.0D+00) then !If morfac=1 then sedimentation factor should be 1 node_change(i3,j) = nint(node_change(i3,j)) else node_change(i3,j) = nint(node_change(i3,j)) endif else node_change(i3,j) = 0.0D+00 endif ! if (s%zb(i3,j) - zb_old(i3,j) > 0.0D+00 .AND. node_change(i3,j) >= 0.5D+00 ) then !Sedimentation do i2=1,node_change(i3,j)+1 ! temp_change(i2,i3) = temp(2,i3) temp_change_cond(i2,i3) = tempcond(2,i3) h_change(i2,i3) = h(2,i3) h_change_cond(i2,i3) = hcond(2,i3) end do do i2=node_change(i3,j)+2,par%y_num temp_change(i2,i3) = temp(i2-node_change(i3,j),i3) temp_change_cond(i2,i3) = tempcond(i2-node_change(i3,j),i3) h_change(i2,i3) = h(i2-node_change(i3,j),i3) h_change_cond(i2,i3) = hcond(i2-node_change(i3,j),i3) end do temp_change(par%y_num,i3) = temp(par%y_num,i3) temp_change_cond(par%y_num,i3) = tempcond(par%y_num,i3) h_change(par%y_num,i3) = h(par%y_num,i3) h_change_cond(par%y_num,i3) = hcond(par%y_num,i3) h_change(1,i3) = h(1,i3) h_change_cond(1,i3) = hcond(1,i3) do i2=1,par%y_num temp(i2,i3) = temp(i2,i3)*0 ! Setting temp and h arrays equal to zero tempcond(i2,i3) = tempcond(i2,i3)*0 h(i2,i3) = h(i2,i3)*0 hcond(i2,i3) = hcond(i2,i3)*0 end do do i2=1,par%y_num temp(i2,i3) = temp_change(i2,i3) ! Updating temp and h arrays with new temperature profile tempcond(i2,i3) = temp_change_cond(i2,i3) h(i2,i3) = h_change(i2,i3) hcond(i2,i3) = h_change_cond(i2,i3) end do ! else if (s%zb(i3,j) - zb_old(i3,j) < 0.0D+00 .AND. node_change(i3,j) >= 0.5D+00) then !Erosion do i2=1,par%y_num-(node_change(i3,j)+1) temp_change(i2,i3) = temp(i2+node_change(i3,j),i3) temp_change_cond(i2,i3) = tempcond(i2+node_change(i3,j),i3) h_change(i2,i3) = h(i2+node_change(i3,j),i3) h_change_cond(i2,i3) = hcond(i2+node_change(i3,j),i3) end do do i2=par%y_num-node_change(i3,j),par%y_num temp_change(i2,i3) = temp(par%y_num,i3) temp_change_cond(i2,i3) = tempcond(par%y_num,i3) h_change(i2,i3) = h(i2-1,i3) h_change_cond(i2,i3) = hcond(i2-1,i3) end do !temp_change(1,i3)=temp(1,i3) h_change(1,i3)=h(1,i3) h_change_cond(1,i3)=hcond(1,i3) h_change(par%y_num,i3)=h(par%y_num,i3) h_change_cond(par%y_num,i3)=hcond(par%y_num,i3) do i2=1,par%y_num temp(i2,i3) = temp(i2,i3)*0 ! Setting temp and h arrays equal to zero tempcond(i2,i3) = tempcond(i2,i3)*0 h(i2,i3) = h(i2,i3)*0 hcond(i2,i3) = hcond(i2,i3)*0 end do do i2=1,par%y_num temp(i2,i3) = temp_change(i2,i3) ! Updating temp and h arrays with new temperature profile tempcond(i2,i3) = temp_change_cond(i2,i3) h(i2,i3) = h_change(i2,i3) hcond(i2,i3) = h_change_cond(i2,i3) end do else return end if end subroutine erosion_sedimentation subroutine bc_test01 (s,par,hws,cfl_conv,dt,y,t,temp,temp_top,h_new,c_heat_conv,temp_f,i3) !*****************************************************************************80 ! !! BC_TEST01 evaluates the boundary conditions for problem 1. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 25 January 2012 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) X_NUM, the number of nodes. ! ! Input, real ( kind = 8 ) X(X_NUM), the node coordinates. ! ! Input, real ( kind = 8 ) T, the current time. ! ! Input, real ( kind = 8 ) temp(X_NUM), the current temperature values. ! ! Output, real ( kind = 8 ) temp(X_NUM), the current temperature values, after boundary ! conditions have been imposed. ! use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par !integer ( kind = 4 ) x_num real ( kind = 8 ) :: temp_f !real ( kind = 8 ) hb ! convection heat transfer coefficient !real ( kind = 8 ) Tw !temperature of water real ( kind = 8 ) :: cfl ! Courant Friedrichs Loewy coefficient real ( kind = 8 ) :: dt !time step !integer ( kind = 4 ) y_num real ( kind = 8 ),allocatable :: y(:) real ( kind = 8 ),allocatable :: t(:) real ( kind = 8 ),allocatable :: temp(:,:),c_heat_conv(:,:),cfl_conv(:,:),hws(:,:),h_new(:,:) !real ( kind = 8 ) x(x_num) real ( kind = 8 ) :: temp_top(:) ! new top bc real ( kind = 8 ) :: dy ! spatial step integer :: i,j,i2,i3 j=1 !convection BC at surface: dy = y(2) - y(1) temp_top(i3) = temp(1,i3) + 2 * dt / (c_heat_conv(1,i3) * dy) * (hws(i3,j) * (par%Tw - temp(1,i3))) + 2 * cfl_conv(i2,i3) * (temp(2,i3) - temp(1,i3)) !i3 is pointer for all nodes in the horizontal temp(1,i3) = temp_top(i3) !i3 is pointer for all nodes in the horizontal if (temp(1,i3) <= temp_f) then h_new(1,i3) = c_heat_conv(1,i3) * temp(1,i3) elseif (temp(1,i3) > temp_f) then h_new(1,i3) = c_heat_conv(1,i3) * temp(1,i3) + par%Li*par%nb end if temp(par%y_num,i3) = par%bottom_temp h_new(par%y_num,i3) = temp(par%y_num,i3) * c_heat_conv(par%y_num,i3) ! return end subroutine bc_test01 subroutine bc_test01conduction (s,par,tempcond,h_newcond,c_heat_cond,temp_f,i3) !***************************************************************************** ! use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par real ( kind = 8 ) :: temp_f real ( kind = 8 ),allocatable :: tempcond(:,:),h_newcond(:,:),c_heat_cond(:,:) integer :: i,j,i2,i3 tempcond(1,i3) = par%Ta*par%n_air_surf !i3 is pointer for all nodes in the horizontal if (tempcond(1,i3) <= temp_f) then h_newcond(1,i3) = c_heat_cond(1,i3) * tempcond(1,i3) elseif (tempcond(1,i3) > temp_f) then h_newcond(1,i3) = c_heat_cond(1,i3) * tempcond(1,i3) + par%Li*par%nb end if tempcond(par%y_num,i3) = par%bottom_temp h_newcond(par%y_num,i3) = tempcond(par%y_num,i3) * c_heat_cond(par%y_num,i3) ! return end subroutine bc_test01conduction subroutine bc_test02conduction (s,par,temp_newcond,surface_temp,surface_temp_storm,c_heat_cond,h_newcond,temp_f,i3,morf_adj) ! !***************************************************************************** ! use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par real ( kind = 8 ),allocatable :: temp_newcond(:,:) real ( kind = 8 ) :: temp_f real*8,dimension(1,1) :: surface_temp,surface_temp_storm real ( kind = 8 ),allocatable :: h_newcond(:,:) !The new enthalpy value at each node real ( kind = 8 ),allocatable :: c_heat_cond(:,:) integer :: i,j,i2,i3 logical :: advance, exist real*8 :: morf_adj ! set bottom to assumed ground temp of -2 degrees C if(par%ht_calibration==1) then !!!!Don't run calibration with Morfac greater than 1 !!!!!!!!! if((par%timet >= par%therm_p) .and. (par%timet< (2*par%therm_p))) then temp_newcond(1,i3) = par%Ta*par%n_air_surf !i3 is pointer for all nodes in the horizontal if (temp_newcond(1,i3) <= temp_f) then h_newcond(1,i3) = c_heat_cond(1,i3) * temp_newcond(1,i3) elseif (temp_newcond(1,i3) > temp_f) then h_newcond(1,i3) = c_heat_cond(1,i3) * temp_newcond(1,i3) + par%Li*par%nb end if ! temp_newcond(par%y_num,i3) = par%bottom_temp h_newcond(par%y_num,i3) = temp_newcond(par%y_num,i3) * c_heat_cond(par%y_num,i3) ! else if(par%timet >= (2*par%therm_p)) then temp_newcond(1,i3) = surface_temp(1,1) if (temp_newcond(1,i3) <= temp_f) then h_newcond(1,i3) = c_heat_cond(1,i3) * temp_newcond(1,i3) elseif (temp_newcond(1,i3) > temp_f) then h_newcond(1,i3) = c_heat_cond(1,i3) * temp_newcond(1,i3) + par%Li*par%nb end if temp_newcond(par%y_num,i3) = par%bottom_temp h_newcond(par%y_num,i3) = temp_newcond(par%y_num,i3) * c_heat_cond(par%y_num,i3) end if else if(par%ht_calibration==0) then if((par%timet >= (par%therm_p/morf_adj)) .and. (par%timet < (2*(par%therm_p/morf_adj)))) then temp_newcond(1,i3) = par%Ta*par%n_air_surf if (temp_newcond(1,i3) <= temp_f) then h_newcond(1,i3) = c_heat_cond(1,i3) * temp_newcond(1,i3) elseif (temp_newcond(1,i3) > temp_f) then h_newcond(1,i3) = c_heat_cond(1,i3) * temp_newcond(1,i3) + par%Li*par%nb end if temp_newcond(par%y_num,i3) = par%bottom_temp h_newcond(par%y_num,i3) = temp_newcond(par%y_num,i3) * c_heat_cond(par%y_num,i3) ! else if((par%timet/morf_adj) >= (2*(par%therm_p/morf_adj))) then temp_newcond(1,i3) = surface_temp_storm(1,1) if (temp_newcond(1,i3) <= temp_f) then h_newcond(1,i3) = c_heat_cond(1,i3) * temp_newcond(1,i3) elseif (temp_newcond(1,i3) > temp_f) then h_newcond(1,i3) = c_heat_cond(1,i3) * temp_newcond(1,i3) + par%Li*par%nb end if temp_newcond(par%y_num,i3) = par%bottom_temp h_newcond(par%y_num,i3) = temp_newcond(par%y_num,i3) * c_heat_cond(par%y_num,i3) end if end if return end subroutine bc_test02conduction ! subroutine bc_test02 (s,par,hws,cfl_conv,dt,y,t,temp,temp_new,temp_top,c_heat_conv,twater,h_new,temp_f,i3,morf_adj) !New routine to help with organization use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par real ( kind = 8 ) :: cfl,twater,watertemp ! Courant Friedrichs Loewy coefficient real ( kind = 8 ) :: dt !time step real ( kind = 8 ) :: c ! heat capacity of soil !integer ( kind = 4 ) y_num real ( kind = 8 ) :: temp_f real ( kind = 8 ),allocatable :: y(:) real ( kind = 8 ),allocatable :: t(:) real ( kind = 8 ),allocatable :: temp(:,:),temp_new(:,:),c_heat_conv(:,:),cfl_conv(:,:),hws(:,:) real ( kind = 8 ),allocatable :: h_new(:,:) !The new enthalpy value at each node !real ( kind = 8 ) x(x_num) real ( kind = 8 ) :: temp_top(:) ! new top bc real ( kind = 8 ) :: dy ! spatial step integer :: i,j,i2,i3 real*8 :: morf_adj j=1 dy = y(2) - y(1) watertemp=0.0D+00 ! if (par%timet < (par%therm_p/morf_adj)+5) then watertemp=par%Tw else if (par%timet > (par%therm_p/morf_adj)+5) then watertemp=twater end if temp_top(i3) = temp(1,i3) + 2 * dt / (c_heat_conv(1,i3) * dy) * (hws(i3,j) * (watertemp - temp(1,i3))) + 2 * cfl_conv(1,i3) * (temp(2,i3) - temp(1,i3)) temp_new(1,i3) = temp_top(i3) if (temp_new(1,i3) <= temp_f) then h_new(1,i3) = c_heat_conv(1,i3) * temp_new(1,i3) elseif (temp_new(1,i3) > temp_f) then h_new(1,i3) = c_heat_conv(1,i3) * temp_new(1,i3) + par%Li*par%nb end if temp_new(par%y_num,i3) = par%bottom_temp ! set bottom to assumed ground temp of -2 degrees C h_new(par%y_num,i3) = temp_new(par%y_num,i3) * c_heat_conv(par%y_num,i3) end subroutine bc_test02 subroutine ic_test01_conv (s,par,temp,i3) ! ! use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par !integer ( kind = 4 ) x_num real ( kind = 8 ), allocatable :: temp(:,:) real ( kind = 8 ), allocatable :: y(:) real ( kind = 8 ), allocatable :: t(:) !integer ( kind = 4 ) y_num !real ( kind = 8 ) h(y_num) !real ( kind = 8 ) x(x_num) integer :: i,j,i2,i3 logical :: exist j=1 if (par%first_run_flag == 0) then inquire(file="Final_Temp_Input.txt", exist=exist) !2nd Run if (exist) then !2nd Run open(unit=88,FILE="Final_Temp_Input.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="asis", ACTION="read") !2nd Run else !2nd Run open(unit=88,FILE="Final_Temp_Input.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="append", ACTION="write") !2nd Run end if !2nd Run do i2=1,par%y_num !2nd Run read(88,114) (temp(i2,i3)) !2nd Run end do !2nd Run 114 format (f10.6) !2nd Run !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! else if (par%first_run_flag == 1) then inquire(file="Temp_Profile_Conv.txt", exist=exist) !1st Run if (exist) then !1st Run open(unit=70,FILE="Temp_Profile_Conv.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="rewind", ACTION="read") !1st Run else !1st Run open(unit=70,FILE="Temp_Profile_Conv.txt", FORM="formatted", ACCESS="Sequential", STATUS="unknown", POSITION="rewind", ACTION="write") !1st Run end if !1st Run do i2=1,par%y_num !1st Run read(70,115) (temp(i2,i3)) !1st Run end do !1st Run 115 format (f10.4) !1st Run close(70) endif ! end subroutine ic_test01_conv subroutine ic_test01_cond (s,par,tempcond,i3) use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par !integer ( kind = 4 ) x_num real ( kind = 8 ), allocatable :: tempcond(:,:) real ( kind = 8 ), allocatable :: y(:) real ( kind = 8 ), allocatable :: t(:) !integer ( kind = 4 ) y_num !real ( kind = 8 ) h(y_num) !real ( kind = 8 ) x(x_num) integer :: i,j,i2,i3 logical :: exist j=1 if (par%first_run_flag == 0) then inquire(file="Final_Tempcond_Input.txt", exist=exist) !2nd Run if (exist) then !2nd Run open(unit=89,FILE="Final_Tempcond_Input.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="asis", ACTION="read") !2nd Run else !2nd Run open(unit=89,FILE="Final_Tempcond_Input.txt", FORM="formatted", ACCESS="sequential", STATUS="old", POSITION="append", ACTION="write") !2nd Run end if !2nd Run do i2=1,par%y_num !2nd Run read(89,115) (tempcond(i2,i3)) !2nd Run end do !2nd Run 115 format (f10.6) !2nd Run !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! else if (par%first_run_flag == 1) then inquire(file="Temp_Profile_Cond.txt", exist=exist) !1st Run if (exist) then !1st Run open(unit=71,FILE="Temp_Profile_Cond.txt", FORM="formatted", ACCESS="sequential", STATUS="unknown", POSITION="rewind", ACTION="read") !1st Run else !1st Run open(unit=71,FILE="Temp_Profile_Cond.txt", FORM="formatted", ACCESS="Sequential", STATUS="unknown", POSITION="rewind", ACTION="write") !1st Run end if !1st Run do i2=1,par%y_num !1st Run read(71,117) (tempcond(i2,i3)) !1st Run end do !1st Run 117 format (f10.4) !1st Run close(71) !1st Run endif end subroutine ic_test01_cond subroutine rhs_test01 (s,par, y, t, value,i3) !This routine is not used!!! !*****************************************************************************80 ! !! RHS_TEST01 evaluates the right hand side for problem 1. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 25 January 2012 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) X_NUM, the number of nodes. ! ! Input, real ( kind = 8 ) X(X_NUM), the node coordinates. ! ! Input, real ( kind = 8 ) T, the current time. ! ! Output, real ( kind = 8 ) VALUE(X_NUM), the source term. ! use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par integer ( kind = 4 ) i3 real ( kind = 8 ) t real ( kind = 8 ) value(par%y_num) real ( kind = 8 ) y(par%y_num) ! value(1:par%y_num) = 0.0D+00 end subroutine rhs_test01 !This routine is not used!!! subroutine fd1d_heat_explicit (s,par,hws,y,t,dt,cfl_conv,temp_f,temp_top,temp,temp_new,h,h_new,c_heat_conv,twater,h_test,i3,morf_adj) !*****************************************************************************80 ! !! FD1D_HEAT_EXPLICIT: Finite difference solution of 1D heat equation, with phase change ! ! Discussion: ! ! This program takes one time step to solve the 1D heat equation ! with an explicit method. ! ! This program solves ! ! dUdT - k * d2UdX2 = F(X,T) ! ! over the interval [A,B] with boundary conditions ! ! U(A,T) = UA(T), ! U(B,T) = UB(T), ! ! over the time interval [T0,T1] with initial conditions ! ! U(X,T0) = U0(X) ! ! The code uses the finite difference method to approximate the ! second derivative in space, and an explicit forward Euler approximation ! to the first derivative in time. ! ! The finite difference form can be written as ! ! U(X,T+dt) - U(X,T) ( U(X-dx,T) - 2 U(X,T) + U(X+dx,T) ) ! ------------------ = F(X,T) + alpha * ------------------------------------ ! dt dx * dx ! ! or, assuming we have solved for all values of U at time T, we have ! ! U(X,T+dt) = U(X,T) ! + cfl * ( U(X-dx,T) - 2 U(X,T) + U(X+dx,T) ) + dt * F(X,T) ! ! Here "cfl" is the Courant-Friedrichs-Loewy coefficient: ! ! cfl = alpha * dt / dx / dx ! ! In order for accurate results to be computed by this explicit method, ! the CFL coefficient must be less than 0.5. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 25 January 2012 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) X_NUM, the number of points to use in the ! spatial dimension. ! ! Input, real ( kind = 8 ) X(X_NUM), the coordinates of the nodes. ! ! Input, real ( kind = 8 ) T, the current time. ! ! Input, real ( kind = 8 ) DT, the size of the time step. ! ! Input, real ( kind = 8 ) CFL, the Courant-Friedrichs-Loewy coefficient, ! computed by FD1D_HEAT_EXPLICIT_CFL. ! ! Input, real ( kind = 8 ) H(X_NUM), the solution at the current time. ! ! Input, external RHS, the function which evaluates the right hand side. ! ! Input, external BC, the function which evaluates the boundary conditions. ! ! Output, real ( kind = 8 ) H_NEW(X_NUM), the solution at time T+DT. ! use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par real ( kind = 8 ) :: cfl,twater !real ( kind = 8 ) c real ( kind = 8 ) :: f(par%y_num) real ( kind = 8 ) :: dt real ( kind = 8 ), allocatable :: h(:,:) !The enthalpy value at each node real ( kind = 8 ), allocatable :: t(:) real ( kind = 8 ), allocatable :: y(:) real ( kind = 8 ), allocatable :: h_new(:,:),h_test(:,:) !The new enthalpy value at each node real ( kind = 8 ), allocatable :: temp_top(:),hws(:,:) real ( kind = 8 ) :: temp_f real ( kind = 8 ), allocatable :: temp(:,:) !Temperature array real ( kind = 8 ), allocatable :: temp_new(:,:),c_heat_conv(:,:),cfl_conv(:,:) !New temperature array integer :: i,j,i2,i3 real*8 :: morf_adj do i2=1,par%y_num if(temp(i2,i3) >=temp_f) then cfl_conv(i2,i3)=((par%ht_ku/par%c_unfrozen_soil)*dt)/(par%y_max/(par%y_num-1))**2 c_heat_conv(i2,i3)=par%c_unfrozen_soil else cfl_conv(i2,i3)=((par%ht_kf/par%c_frozen_soil)*dt)/(par%y_max/(par%y_num-1))**2 c_heat_conv(i2,i3)=par%c_frozen_soil end if end do ! !h_new(2:par%y_num-1,i3) = h(2:par%y_num-1,i3) + cfl_conv(2:par%y_num-1,i3) * c_heat_conv(2:par%y_num-1,i3) * (temp(1:par%y_num-2,i3) - 2.0D+00 * temp(2:par%y_num-1,i3) + temp(3:par%y_num,i3) ) ! do i2=2,par%y_num-1 h_new(i2,i3) = h(i2,i3) + (cfl_conv(i2-1,i3) * c_heat_conv(i2-1,i3) * temp(i2-1,i3)) + (cfl_conv(i2+1,i3) * c_heat_conv(i2+1,i3) * temp(i2+1,i3)) - 2.0D+00 * (cfl_conv(i2,i3) * c_heat_conv(i2,i3) * temp(i2,i3)) h_test(i2,i3) = h_new(i2,i3) end do ! call bc_test02 (s,par,hws,cfl_conv,dt,y,t,temp,temp_new,temp_top,c_heat_conv,twater,h_new,temp_f,i3,morf_adj) ! do i2=2,par%y_num-1 if (h_new(i2,i3) < c_heat_conv(i2,i3) * temp_f) then temp_new(i2,i3) = h_new(i2,i3)/c_heat_conv(i2,i3) else if ((h_new(i2,i3) >= temp_f * c_heat_conv(i2,i3)) .and. (h_new(i2,i3) < temp_f * c_heat_conv(i2,i3) + par%Li*par%nb)) then temp_new(i2,i3) = temp_f else temp_new(i2,i3) = (h_new(i2,i3) - (c_heat_conv(i2,i3) * temp_f) - par%Li*par%nb)/ c_heat_conv(i2,i3) ! really this is (Hn - (cs-cl)Tf - L) / cl end if end do !temp_new(1,i3) = temp(1,i3) !temp_new(par%y_num,i3) = temp(par%y_num,i3) end subroutine fd1d_heat_explicit subroutine fd1d_heat_explicit_conduction (s,par,y,t,dt,cfl_cond,temp_f,tempcond,temp_newcond,hcond,h_newcond,c_heat_cond,surface_temp,surface_temp_storm,i3,morf_adj) use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par real ( kind = 8 ) :: cfl !real ( kind = 8 ) c real ( kind = 8 ) :: f(par%y_num) real ( kind = 8 ) :: dt real*8,dimension(1,1) :: surface_temp, surface_temp_storm real ( kind = 8 ), allocatable :: hcond(:,:),h_newcond(:,:) !The enthalpy value at each node real ( kind = 8 ), allocatable :: t(:) real ( kind = 8 ), allocatable :: y(:) real ( kind = 8 ) :: temp_f real ( kind = 8 ), allocatable :: tempcond(:,:) !Temperature array real ( kind = 8 ), allocatable :: temp_newcond(:,:),c_heat_cond(:,:),cfl_cond(:,:) !New temperature array integer :: i,j,i2,i3 real*8 :: morf_adj do i2=1,par%y_num if(tempcond(i2,i3) >=temp_f) then cfl_cond(i2,i3)=((par%ht_ku/par%c_unfrozen_soil)*dt)/(par%y_max/(par%y_num-1))**2 c_heat_cond(i2,i3)=par%c_unfrozen_soil else cfl_cond(i2,i3)=((par%ht_kf/par%c_frozen_soil)*dt)/(par%y_max/(par%y_num-1))**2 c_heat_cond(i2,i3)=par%c_frozen_soil end if end do ! h_newcond(2:par%y_num-1,i3) = hcond(2:par%y_num-1,i3) + cfl_cond(2:par%y_num-1,i3) * c_heat_cond(2:par%y_num-1,i3) * (tempcond(1:par%y_num-2,i3) - 2.0D+00 * tempcond(2:par%y_num-1,i3) + tempcond(3:par%y_num,i3) ) do i2=2,par%y_num-1 h_newcond(i2,i3) = hcond(i2,i3) + (cfl_cond(i2-1,i3) * c_heat_cond(i2-1,i3) * tempcond(i2-1,i3)) + (cfl_cond(i2+1,i3) * c_heat_cond(i2+1,i3) * tempcond(i2+1,i3)) - 2.0D+00 * (cfl_cond(i2,i3) * c_heat_cond(i2,i3) * tempcond(i2,i3)) end do call bc_test02conduction (s,par,temp_newcond,surface_temp,surface_temp_storm,c_heat_cond,h_newcond,temp_f,i3,morf_adj) do i2=2,par%y_num-1 if (h_newcond(i2,i3) < c_heat_cond(i2,i3) * temp_f) then temp_newcond(i2,i3) = h_newcond(i2,i3)/c_heat_cond(i2,i3) else if ((h_newcond(i2,i3) >= temp_f * c_heat_cond(i2,i3)) .and. (h_newcond(i2,i3) < temp_f * c_heat_cond(i2,i3) + par%Li*par%nb)) then temp_newcond(i2,i3) = temp_f else temp_newcond(i2,i3) = (h_newcond(i2,i3) - (c_heat_cond(i2,i3) * temp_f) - par%Li*par%nb)/ c_heat_cond(i2,i3) ! really this is (Hn - (cs-cl)Tf - L) / cl end if end do !temp_newcond(1,i3) = tempcond(1,i3) !temp_newcond(par%y_num,i3) = tempcond(par%y_num,i3) end subroutine fd1d_heat_explicit_conduction subroutine fd1d_heat_explicit_cfl (s,par,t_min,y_min) !*****************************************************************************80 ! !! FD1D_HEAT_EXPLICIT_CFL: compute the Courant-Friedrichs-Loewy coefficient. ! ! Discussion: ! ! The equation to be solved has the form: ! ! dUdT - alpha * d2UdX2 = F(X,T) ! ! over the interval [X_MIN,X_MAX] with boundary conditions ! ! U(X_MIN,T) = U_X_MIN(T), ! U(X_MIN,T) = U_X_MAX(T), ! ! over the time interval [T_MIN,T_MAX] with initial conditions ! ! U(X,T_MIN) = U_T_MIN(X) ! ! The code uses the finite difference method to approximate the ! second derivative in space, and an explicit forward Euler approximation ! to the first derivative in time. ! ! The finite difference form can be written as ! ! U(X,T+dt) - U(X,T) ( U(X-dx,T) - 2 U(X,T) + U(X+dx,T) ) ! ------------------ = F(X,T) + alpha * ------------------------------------ ! dt dx * dx ! ! or, assuming we have solved for all values of U at time T, we have ! ! U(X,T+dt) = U(X,T) ! + cfl * ( U(X-dx,T) - 2 U(X,T) + U(X+dx,T) ) + dt * F(X,T) ! ! Here "cfl" is the Courant-Friedrichs-Loewy coefficient: ! ! cfl = alpha * dt / dx / dx ! ! In order for accurate results to be computed by this explicit method, ! the CFL coefficient must be less than 0.5! ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 24 January 2012 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! George Lindfield, John Penny, ! Numerical Methods Using MATLAB, ! Second Edition, ! Prentice Hall, 1999, ! ISBN: 0-13-012641-1, ! LC: QA297.P45. ! ! Parameters: ! ! Input, real ( kind = 8 ) alpha, the thermal conductivity coefficient divided by heat capacity. ! ! Input, integer ( kind = 4 ) T_NUM, the number of time values, including ! the initial value. ! ! Input, real ( kind = 8 ) T_MIN, T_MAX, the minimum and maximum times. ! ! Input, integer ( kind = 4 ) X_NUM, the number of nodes. ! ! Input, real ( kind = 8 ) X_MIN, X_MAX, the minimum and maximum spatial ! coordinates. ! ! Output, real ( kind = 8 ) CFL, the Courant-Friedrichs-Loewy coefficient. ! use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par integer ( kind = 8 ) i3 real ( kind = 8 ) cfl_frozen_soil,cfl_unfrozen_soil real ( kind = 8 ) dy real ( kind = 8 ) dt !real ( kind = 8 ) t_max real ( kind = 8 ) t_min !integer ( kind = 4 ) t_num !real ( kind = 8 ) y_max real ( kind = 8 ) y_min !integer ( kind = 4 ) y_num dy = (par%y_max - y_min ) / real (par%y_num - 1, kind = 8 ) dt = (par%t_max - t_min ) / real (par%t_num - 1, kind = 8 ) ! !( t_max - t_min ) / real ( t_num - 1, kind = 8 ) ! ! Check the CFL condition, print out its value, and quit if it is too large. ! cfl_frozen_soil = ((par%ht_kf/par%c_frozen_soil)*dt)/(par%y_max/(par%y_num-1))**2 cfl_unfrozen_soil = ((par%ht_ku/par%c_unfrozen_soil)*dt)/(par%y_max/(par%y_num-1))**2 write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) ' CFL stability criterion value for frozen soil = ', cfl_frozen_soil write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) ' CFL stability criterion value for unfrozen soil = ', cfl_unfrozen_soil write ( *, '(a)' ) ' ' if(0.5D+00 <= cfl_frozen_soil) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FD1D_HEAT_EXPLICIT_CFL - Fatal error!' write ( *, '(a)' ) ' CFL condition failed.' write ( *, '(a)' ) ' 0.5 <= alpha * dT / dX / dX = CFL.' stop else if(0.5D+00 <= cfl_unfrozen_soil) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FD1D_HEAT_EXPLICIT_CFL - Fatal error!' write ( *, '(a)' ) ' CFL condition failed.' write ( *, '(a)' ) ' 0.5 <= alpha * dT / dX / dX = CFL.' stop end if end subroutine fd1d_heat_explicit_cfl subroutine find_thaw_depth (s,par,zbmean,temp_f,y,thaw_depth,i3,thaw_depth_change,thaw_depth_old,zb_old,node_adj,temp_final,delta_new,thaw_tv,thaw_tc,thaw_hv,thaw_hc,temp,tempcond,h,hcond,temp_node_adj,morf_adj,thaw_depth_change_vert) ! !*****************************************************************************80 ! !! Find_thaw_depth finds the depth to thaw (the depth node at which the temperature ! switches from being above freezing to below freezing). ! use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par integer ( kind = 4 ) :: i,i2,j,i3 !integer ( kind = 4 ) y_num real ( kind = 8 ),allocatable :: temp_final(:,:),thaw_depth_old(:,:),thaw_depth(:,:),thaw_depth_change(:,:),zb_old(:,:),node_adj(:,:),zbmean(:,:),delta_new(:,:),thaw_depth_change_vert(:,:) real*8,dimension(:,:),allocatable :: thaw_tv,thaw_tc,thaw_hv,thaw_hc,temp,tempcond,h,hcond integer , dimension(:,:),allocatable :: temp_node_adj real ( kind = 8 ),allocatable :: y(:) real ( kind = 8 ) :: temp_f !the freezing temperature of the soil (moisture) real*8 :: morf_adj,erosion_factor,sedimentation_factor ! if(delta_new(i3,j) >= par%slope_threshold) then erosion_factor=par%E_factor_high sedimentation_factor=par%S_factor_high else erosion_factor=par%E_factor_low sedimentation_factor=par%S_factor_low endif j=1 if(s%zb(i3,j) - zb_old(i3,j) < 0.0D+00) then node_adj(i3,j)=abs(((s%zb(i3,j) - zb_old(i3,j))/(par%y_max/(par%y_num-1))) * COSD(delta_new(i3,j)) * erosion_factor) ! if((node_adj(i3,j)*(par%y_max/(par%y_num-1))) > thaw_depth_old(i3,j)) then node_adj(i3,j) = floor(thaw_depth_old(i3,j)/(par%y_max/(par%y_num-1))) else if(par%morfac >2.0D+00) then node_adj(i3,j) = nint(node_adj(i3,j)) else node_adj(i3,j) = nint(node_adj(i3,j)) endif endif ! else if (s%zb(i3,j) - zb_old(i3,j) > 0.0D+00) then node_adj(i3,j)=abs(((s%zb(i3,j) - zb_old(i3,j))/(par%y_max/(par%y_num-1))) * COSD(delta_new(i3,j)) * sedimentation_factor) if(par%morfac >2.0D+00) then node_adj(i3,j) = nint(node_adj(i3,j)) else node_adj(i3,j) = nint(node_adj(i3,j)) endif else node_adj(i3,j) = 0.0D+00 endif ! if (temp_final(1,i3) <= temp_f) then! the thaw depth is 0 (or x(1)) ! For CONVECTION CASE!!!!!! thaw_depth(i3,j) = y(1) else do i2 = 2, par%y_num if (temp_final(i2,i3) <= temp_f) then !the ground is frozen. We are assuming the thaw depth is between this and the previous node. Linearly extrapolate for the depth of the thaw thaw_depth(i3,j) = y(i2) + (temp_f - temp_final(i2,i3)) * (y(i2) - y(i2-1)) / (temp_final(i2,i3) - temp_final(i2-1,i3)) if(s%zb(i3,j) - zb_old(i3,j) > 0.0D+00) then thaw_depth_change(i3,j)=((thaw_depth(i3,j) - abs(node_adj(i3,j))*(par%y_max/(par%y_num-1))) - thaw_depth_old(i3,j))*par%tf else if(s%zb(i3,j) - zb_old(i3,j) < 0.0D+00) then thaw_depth_change(i3,j)=((thaw_depth(i3,j) + abs(node_adj(i3,j))*(par%y_max/(par%y_num-1))) - thaw_depth_old(i3,j))*par%tf else thaw_depth_change(i3,j) = thaw_depth(i3,j) - thaw_depth_old(i3,j) !Change in thawdepth due to heat transfer alone end if !thaw_depth_change(i3,j) = thaw_depth_change(i3,j)/COSD(delta_new) thaw_depth_old(i3,j)=thaw_depth_old(i3,j)*0 thaw_depth_old(i3,j)=thaw_depth(i3,j) thaw_depth_change_vert(i3,j) = thaw_depth_change(i3,j)/COSD(delta_new(i3,j)) !call thermal_acceleration(s,par,thaw_depth_change,temp,tempcond,h,hcond,thaw_tv,thaw_tc,thaw_hv,thaw_hc,temp_node_adj,i3,morf_adj) return !exit out of do loop at the top thaw boundary else !this is a bit wierd, but should work, every iteration that is thawed will set thaw depth to the !bottom, but this will be overridden if frozen soil is encountered in a later iteration thaw_depth(i3,j) = y(par%y_num) ! really, thaw depth is greater than the deepest x in this case... end if end do end if if(s%zb(i3,j) - zb_old(i3,j) >= 0.0D+00) then thaw_depth_change(i3,j)=((thaw_depth(i3,j) - abs(node_adj(i3,j))*(par%y_max/(par%y_num-1))) - thaw_depth_old(i3,j))*par%tf else if(s%zb(i3,j) - zb_old(i3,j) < 0.0D+00) then thaw_depth_change(i3,j)=((thaw_depth(i3,j) + abs(node_adj(i3,j))*(par%y_max/(par%y_num-1))) - thaw_depth_old(i3,j))*par%tf else thaw_depth_change(i3,j) = thaw_depth(i3,j) - thaw_depth_old(i3,j) !Change in thawdepth due to heat transfer alone end if !thaw_depth_change(i3,j) = thaw_depth_change(i3,j)/COSD(delta_new) thaw_depth_old(i3,j)=thaw_depth_old(i3,j)*0 thaw_depth_old(i3,j)=thaw_depth(i3,j) thaw_depth_change_vert(i3,j) = thaw_depth_change(i3,j)/COSD(delta_new(i3,j)) end subroutine find_thaw_depth subroutine thermal_acceleration(s,par,thaw_depth_change,temp,tempcond,h,hcond,thaw_tv,thaw_tc,thaw_hv,thaw_hc,temp_node_adj,i3,morf_adj) !Thermal acceleration is not used at this time!!!! use params !Thermal acceleration is not used at this time!!!! Difficult to make this concept work! use spaceparams !Thermal acceleration is not used at this time!!!! Difficult to make this concept work! use xmpi_module !Thermal acceleration is not used at this time!!!! Difficult to make this concept work! use interp !Thermal acceleration is not used at this time!!!! Difficult to make this concept work! use paramsconst !Thermal acceleration is not used at this time!!!! Difficult to make this concept work! implicit none type(spacepars),target :: s type(parameters) :: par ! integer ( kind = 4 ) :: i,i2,j,i3 integer , dimension(:,:),allocatable :: temp_node_adj real ( kind = 8 ),allocatable :: thaw_depth_change(:,:) real*8,dimension(:,:),allocatable :: thaw_tv,thaw_tc,thaw_hv,thaw_hc real*8,dimension(:,:),allocatable :: temp,tempcond,h,hcond real ( kind = 8 ),allocatable :: y(:) real ( kind = 8 ) :: temp_f real*8 :: morf_adj j=1 if (par%morfac >= 2.0D+00) then temp_node_adj(i3,j) = (thaw_depth_change(i3,j)/(par%y_max/(par%y_num-1))) j=1 do i2=1,temp_node_adj(i3,j)+1 thaw_tv(i2,i3) = temp(1,i3) thaw_hv(i2,i3) = h(2,i3) thaw_hc(i2,i3) = hcond(2,i3) thaw_tc(i2,i3) = tempcond(1,i3) end do j=1 do i2=temp_node_adj(i3,j)+2,par%y_num thaw_tv(i2,i3) = temp(i2-temp_node_adj(i3,j),i3) thaw_hv(i2,i3) = h(i2-temp_node_adj(i3,j),i3) thaw_hc(i2,i3) = hcond(i2-temp_node_adj(i3,j),i3) thaw_tc(i2,i3) = tempcond(i2-temp_node_adj(i3,j),i3) end do thaw_tv(par%y_num,i3) = temp(par%y_num,i3) thaw_tc(par%y_num,i3) = tempcond(par%y_num,i3) thaw_hv(par%y_num,i3) = h(par%y_num,i3) thaw_hc(par%y_num,i3) = hcond(par%y_num,i3) thaw_hv(1,i3) = h(1,i3) thaw_hc(1,i3) = hcond(1,i3) do i2=1,par%y_num temp(i2,i3) = temp(i2,i3)*0 ! Setting temp and h arrays equal to zero tempcond(i2,i3) = tempcond(i2,i3)*0 h(i2,i3) = h(i2,i3)*0 hcond(i2,i3) = hcond(i2,i3)*0 end do do i2=1,par%y_num temp(i2,i3) = thaw_tv(i2,i3) ! Updating temp and h arrays with new temperature profile tempcond(i2,i3) = thaw_tc(i2,i3) h(i2,i3) = thaw_hv(i2,i3) hcond(i2,i3) = thaw_hc(i2,i3) end do else if (par%morfac < 2.0D+00) then end if end subroutine thermal_acceleration subroutine get_unit (s,par,iunit,i3) !*****************************************************************************80 ! !! GET_UNIT returns a free FORTRAN unit number. ! ! Discussion: ! ! A "free" FORTRAN unit number is an integer between 1 and 99 which ! is not currently associated with an I/O device. A free FORTRAN unit ! number is needed in order to open a file with the OPEN command. ! ! If IUNIT = 0, then no free FORTRAN unit could be found, although ! all 99 units were checked (except for units 5, 6 and 9, which ! are commonly reserved for console I/O). ! ! Otherwise, IUNIT is an integer between 1 and 99, representing a ! free FORTRAN unit. Note that GET_UNIT assumes that units 5 and 6 ! are special, and will never return those values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 26 October 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer ( kind = 4 ) IUNIT, the free unit number. ! use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par integer ( kind = 4 ) i,i3 integer ( kind = 4 ) ios integer ( kind = 4 ) iunit logical lopen iunit = 0 do i = 1, 99 if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then inquire ( unit = i, opened = lopen, iostat = ios ) if ( ios == 0 ) then if ( .not. lopen ) then iunit = i return end if end if end if end do end subroutine get_unit subroutine r8mat_write (s,par, output_filename, m, n, table,i3) !*****************************************************************************80 ! !! R8MAT_WRITE writes an R8MAT file. ! ! Discussion: ! ! An R8MAT is an array of R8 values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 31 May 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) OUTPUT_FILENAME, the output file name. ! ! Input, integer ( kind = 4 ) M, the spatial dimension. ! ! Input, integer ( kind = 4 ) N, the number of points. ! ! Input, real ( kind = 8 ) TABLE(M,N), the table data. ! use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par integer ( kind = 4 ) m integer ( kind = 4 ) n integer ( kind = 4 ) j,i3 character ( len = * ) output_filename integer ( kind = 4 ) output_status integer ( kind = 4 ) output_unit character ( len = 30 ) string real ( kind = 8 ) table(m,n) ! ! Open the file. ! call get_unit (s,par,output_unit,i3) open ( unit = output_unit, file = output_filename, & status = 'replace', iostat = output_status ) if ( output_status /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_WRITE - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the output file "' // & trim ( output_filename ) // '" on unit ', output_unit output_unit = -1 stop end if ! ! Create a format string. ! ! For less precision in the output file, try: ! ! '(', m, 'g', 14, '.', 6, ')' ! if ( 0 < m .and. 0 < n ) then write ( string, '(a1,i8,a1,i8,a1,i8,a1)' ) '(', m, 'g', 24, '.', 16, ')' ! ! Write the data. ! do j = 1, n write ( output_unit, string ) table(1:m,j) end do end if ! ! Close the file. ! close ( unit = output_unit ) end subroutine r8mat_write subroutine r8vec_linspace (s,par, n, a_first, a_last, a) !*****************************************************************************80 ! !! R8VEC_LINSPACE creates a vector of linearly spaced values. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 14 March 2011 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of entries in the vector. ! ! Input, real ( kind = 8 ) A_FIRST, A_LAST, the first and last entries. ! ! Output, real ( kind = 8 ) A(N), a vector of linearly spaced data. ! use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par integer ( kind = 4 ) n,i3 real ( kind = 8 ) a(n) real ( kind = 8 ) a_first real ( kind = 8 ) a_last integer ( kind = 4 ) i if ( n == 1 ) then a(1) = ( a_first + a_last ) / 2.0D+00 else do i = 1, n a(i) = ( real ( n - i, kind = 8 ) * a_first & + real ( i - 1, kind = 8 ) * a_last ) & / real ( n - 1, kind = 8 ) end do end if end subroutine r8vec_linspace subroutine r8vec_write(s,par,thaw_depth,i3,thaw_depth_change,node_change,temp,h,hcond,tempcond,zbmean,net_flux_total,twater_next,rad_flux_total,vmagmean,glmmean,eulmean,soil_flux_total,heat_vel,adv_flux_total,zb_old,thaw_depth_old,delta_new,surfzone_Vol_total,wet_length,temp_f,h_test,surface_temp) !*****************************************************************************80 ! !! R8VEC_WRITE writes an R8VEC file. ! ! Discussion: ! ! An R8VEC is a vector of R8 values. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 10 July 2011 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) OUTPUT_FILENAME, the output file name. ! ! Input, integer ( kind = 4 ) N, the number of points. ! ! Input, real ( kind = 8 ) X(N), the data. ! use params use spaceparams use xmpi_module use interp use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par !integer ( kind = 4 ) n real*8,dimension(1,1) :: surface_temp real ( kind = 8 ) :: net_flux_total,twater_next,rad_flux_total,soil_flux_total,heat_vel,adv_flux_total,surfzone_Vol_total,wet_length,temp_f real ( kind = 8 ),allocatable :: thaw_depth(:,:),thaw_depth_change(:,:),node_change(:,:),temp(:,:),h(:,:),hcond(:,:),tempcond(:,:),zbmean(:,:),zb_old(:,:),thaw_depth_old(:,:),delta_new(:,:),h_test(:,:) real*8,dimension(:,:),allocatable :: vmagmean,glmmean,eulmean integer ( kind = 4 ) :: i,j,i3,i2 !character ( len = * ) output_filename !integer ( kind = 4 ) output_status !integer ( kind = 4 ) output_unit logical :: advance, exist ! inquire(file="Thaw_depth.txt", exist=exist) inquire(file="Thaw_depth_change.txt", exist=exist) inquire(file="Thaw_depth_old.txt", exist=exist) inquire(file="Node_change.txt", exist=exist) inquire(file="Temp.txt", exist=exist) inquire(file="Enthalpy.txt", exist=exist) inquire(file="Enthalpy_Conduction.txt", exist=exist) inquire(file="Temperature_Conduction.txt", exist=exist) inquire(file="zbmean.txt", exist=exist) inquire(file="net_flux_total.txt", exist=exist) inquire(file="Temp_Water_New.txt", exist=exist) inquire(file="rad_flux_total.txt", exist=exist) inquire(file="adv_flux_total.txt", exist=exist) inquire(file="soil_flux_total.txt", exist=exist) inquire(file="Bed_Level.txt", exist=exist) inquire(file="zb_old_data.txt", exist=exist) inquire(file="structdepth.txt", exist=exist) inquire(file="delta_new.txt", exist=exist) inquire(file="surfzone_volume.txt", exist=exist) inquire(file="wet_length.txt", exist=exist) inquire(file="temp_f.txt", exist=exist) inquire(file="surf_out.txt", exist=exist) if (exist) then open(unit=54,FILE="Thaw_depth.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=55,FILE="Thaw_depth_change.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=56,FILE="Node_change.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=57,FILE="Temp_Submerged.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=58,FILE="Enthalpy_Submerged.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=60,FILE="Enthalpy_Conduction.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=61,FILE="Temperature_Conduction.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=63,FILE="zbmean.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=64,FILE="net_flux_total.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=65,FILE="Temp_Water_New.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=66,FILE="rad_flux_total.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=67,FILE="adv_flux_total.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=68,FILE="soil_flux_total.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=69,FILE="Bed_Level.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=128,FILE="zb_old_data.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=129,FILE="Thaw_depth_old.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=130,FILE="structdepth.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=131,FILE="delta_new.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=132,FILE="surfzone_volume.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=133,FILE="wet_length.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=134,FILE="temp_f.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=135,FILE="surf_out.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") else open(unit=54,FILE="Thaw_depth.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=55,FILE="Thaw_depth_change.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=56,FILE="Node_change.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=57,FILE="Temp_Submerged.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=58,FILE="Enthalpy_Submerged.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=60,FILE="Enthalpy_Conduction.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=61,FILE="Temperature_Conduction.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=63,FILE="zbmean.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=64,FILE="net_flux_total.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=65,FILE="Temp_Water_New.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=66,FILE="rad_flux_total.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=67,FILE="adv_flux_total.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=68,FILE="soil_flux_total.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=69,FILE="Bed_Level.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=128,FILE="zb_old_data.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=129,FILE="Thaw_depth_old.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=130,FILE="structdepth.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=131,FILE="delta_new.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=132,FILE="surfzone_volume.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=133,FILE="wet_length.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=134,FILE="temp_f.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=135,FILE="surf_out.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") end if do i=1,s%nx+1 write(54,111) (thaw_depth(i,j),j=1,1) write(55,111) (thaw_depth_change(i,j),j=1,1) write(56,111) (node_change(i,j),j=1,1) write(69,111) (s%zb(i,j),j=1,1) write(128,111) (zb_old(i,j),j=1,1) write(129,111) (thaw_depth_old(i,j),j=1,1) write(130,111) (s%structdepth(i,j),j=1,1) write(131,111) (delta_new(i,j),j=1,1) end do 111 format (f14.8) do i2=1,par%y_num write (57,112) (temp(i2,i3),i3=par%tnode_subm,par%tnode_subm) ! 1000,1000 write (61,112) (tempcond(i2,i3),i3=par%tnode_bluff,par%tnode_bluff) ! 3300,3300 end do 112 format (f16.8) do i2=1,par%y_num write (58,113) (h(i2,i3),i3=par%tnode_subm,par%tnode_subm) !1000,1000 write (60,113) (hcond(i2,i3),i3=par%tnode_bluff,par%tnode_bluff) !3300,3300 end do 113 format (f20.8) do i=1,s%nx+1 write(63,118) (zbmean(i,j),j=1,1) end do 118 format (f10.6) ! write(64,119) net_flux_total write(66,119) rad_flux_total write(68,119) soil_flux_total write(65,119) twater_next write(67,119) adv_flux_total write(132,119) surfzone_Vol_total write(133,119) wet_length write(134,119) temp_f write(135,119) surface_temp(1,1) 119 format (f16.4) close(54) close(55) close(56) close(57) close(58) close(60) close(61) close(63) close(64) close(65) close(66) close(67) close(68) close(69) close(128) close(129) close(130) close(131) close(132) close(133) close(134) close(135) !if ( output_status /= 0 ) then !write ( *, '(a)' ) ' ' !write ( *, '(a)' ) 'R8VEC_WRITE - Fatal error!' !write ( *, '(a,i8)' ) ' Could not open the output file "' // & !trim ( output_filename ) // '" on unit ', output_unit !output_unit = -1 !stop !end if ! end subroutine r8vec_write ! subroutine temp_enthrope_write(s,par,temp,tempcond,h,hcond,thaw_depth) use params use interp use spaceparams use xmpi_module use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par real*8,dimension(:,:),allocatable :: temp,tempcond,h,hcond,thaw_depth !integer ( kind = 4 ) y_num logical :: advance, exist integer :: i,j,j2,i2,i3 do i3=1,s%nx+1 j=1 inquire(file="Final_Temp_Output.txt", exist=exist) inquire(file="Final_Tempcond_Output.txt", exist=exist) inquire(file="Final_h_Output.txt", exist=exist) inquire(file="Final_hcond_Output.txt", exist=exist) if (exist) then open(unit=74,FILE="Final_Temp_Output.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=75,FILE="Final_Tempcond_Output.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=76,FILE="Final_h_Output.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=77,FILE="Final_hcond_Output.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") else open(unit=74,FILE="Final_Temp_Output.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=75,FILE="Final_Tempcond_Output.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=76,FILE="Final_h_Output.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=77,FILE="Final_hcond_Output.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") end if do i2=1,par%y_num write(74,122) temp(i2,i3) write(75,122) tempcond(i2,i3) write(76,122) h(i2,i3) write(77,122) hcond(i2,i3) end do end do 122 format (f10.4) ! j=1 inquire(file="Structthaw_Output.txt", exist=exist) inquire(file="Bed_Level_Output.txt", exist=exist) inquire(file="Thaw_Depth_Output.txt", exist=exist) if (exist) then open(unit=78,FILE="Structthaw_Output.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=79,FILE="Bed_Level_Output.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") open(unit=80,FILE="Thaw_Depth_Output.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") else open(unit=78,FILE="Structthaw_Output.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=79,FILE="Bed_Level_Output.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") open(unit=80,FILE="Thaw_Depth_Output.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") end if do i3=1,s%nx+1 write(78,124) (s%structdepth(i3,j),j=1,1) write(79,124) (s%zb(i3,j),j=1,1) write(80,124) (thaw_depth(i3,j),j=1,1) end do 124 format (f10.4) close(74) close(75) close(76) close(77) close(78) close(79) close(80) end subroutine temp_enthrope_write subroutine write_test(s,par,tracker_vert_change) use params use interp use spaceparams use xmpi_module use paramsconst implicit none type(spacepars),target :: s type(parameters) :: par real*8,dimension(:,:),allocatable :: tracker_vert_change !integer ( kind = 4 ) y_num logical :: advance, exist integer :: i,j,j2,i2,i3 !real*8 :: v_factor,heat_vel inquire(file="tracker_vert_change.txt", exist=exist) if (exist) then open(unit=81,FILE="tracker_vert_change.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") !open(unit=82,FILE="write_test2.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") !open(unit=83,FILE="write_test3.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") !open(unit=84,FILE="write_test4.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") !open(unit=87,FILE="write_test5.txt", FORM="formatted", ACCESS="Sequential", STATUS="old", POSITION="append", ACTION="write") else open(unit=81,FILE="tracker_vert_change.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") !open(unit=82,FILE="write_test2.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") !open(unit=83,FILE="write_test3.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") !open(unit=84,FILE="write_test4.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") !open(unit=87,FILE="write_test5.txt", FORM="formatted", ACCESS="Sequential", STATUS="new", ACTION="write") end if do i=1,s%nx+1 write(81,123) (tracker_vert_change(i,j),j=1,1) end do !write(82,123) (hws_old(i3,j),j=1,1) 123 format (f16.4) !do i=1,s%nx+1 ! write(83,350) (temp_final(1,i)) !end do !350 format (f16.4) ! write(84,351) (v_factor) ! write(87,351) (heat_vel) !351 format (f16.4) !close(81) end subroutine write_test end module morphevolution