! Copyright 2014 College of William and Mary ! ! Licensed under the Apache License, Version 2.0 (the "License"); ! you may not use this file except in compliance with the License. ! You may obtain a copy of the License at ! ! http://www.apache.org/licenses/LICENSE-2.0 ! ! Unless required by applicable law or agreed to in writing, software ! distributed under the License is distributed on an "AS IS" BASIS, ! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. ! See the License for the specific language governing permissions and ! limitations under the License. !===================================================================== !===================================================================== ! MORSELFE INITIALIZATION SUBROUTINES ! ! subroutine sed_alloc ! subroutine sed_init ! !===================================================================== !===================================================================== SUBROUTINE sed_alloc() !--------------------------------------------------------------------! ! This subroutine allocates and pre-initialize sediment model arrays ! ! ! ! Adapted from former subroutines initialize_scalars and ! ! initialize_ocean (former init_sed.F90), other allocation formerly ! ! done within schism_init.F90 were also moved within this subroutine. ! ! ! ! Author: florian ganthy (fganthy@lnec.pt ; florian.ganthy@gmail.com)! ! Date: 2013/06/07 ! ! ! ! History: ! ! *** Previous history (former subroutines) *** ! ! ! ! 2012/12 - F.Ganthy : form homogenisation of sediments routines ! ! 2012/12 - F.Ganthy : modifications for Bottom Composition ! ! Generation (BCG) purpose (added output ! ! files - bed median grain size, bed fraction)! ! 2013/01 - F.Ganthy : Implementation of roughness predictor ! ! 2013/03 - F.Ganthy : Implementation of wave-induced bedload ! ! transport ! ! 2013/04 - F.Ganthy : Implementation of wave-current bottom stress! ! 2013/05 - F.Ganthy : Updates related to ripple predictor ! ! 2013/05 - F.Ganthy : Added node-centered volume control area ! ! 2013/05 - F.Ganthy : Updates to the ripple predictor: ! ! - Changes on the total bedform ! ! roughness computation ! ! - Add wave-ripple computation from ! ! Nielsen (1992) ! ! ! ! *** Current history (former subroutines) *** ! ! ! 2020/02 - B.Mengual : > Allocation of erosion/deposition fluxes, ! ! porosity, morphological ramp value, ! ! Elfrink variables, wave parameters, ! ! bedload flux caused by wave acceleration ! ! > Implementation of a morphological ramp ! ! value read in imorphogrid.gr3 if sed_morph=1! ! > Introduction of sedlay_ini_opt option for ! ! the initialization of layer thickness and ! ! porosity ! ! > Different porosity initialization ! ! depending on poro_option ! ! 2020/07 - M.Pezerat,B.Mengual : introduction of tau_m for ! ! bottom shear stress computations ! ! ! !--------------------------------------------------------------------! USE sed_mod USE schism_glbl, ONLY: rkind,dav,dave,nea,npa,nvrt,mnei_p, & &ntrs,irange_tr,ipgl,ielg,i34,elnode,np_global, & & ifile_char,ifile_len,area,np,nne,indel,iself,nnp,indnd,nxq, & & isbnd,rough_p,errmsg,xnd,ynd,xcj,ycj,xctr,yctr,elside,ics,xel,yel, & &in_dir,out_dir,len_in_dir,len_out_dir USE schism_msgp, ONLY: myrank,parallel_abort,exchange_p2d IMPLICIT NONE !- Local variables --------------------------------------------------! LOGICAL :: lexist REAL(rkind), PARAMETER :: IniVal = 0.0d0 CHARACTER(len=48) :: inputfile INTEGER :: i,j,k,ie,jj,id,id1,id2,id3,nd,indx,m,nwild(3),nwild2(3) INTEGER :: ised,ic,itmp,istat REAL(rkind) :: aux1,aux2,xtmp,ytmp,tmp1,cff1,cff2,cff3,cff4,cff5,cff6 REAL(rkind) :: bed_frac_sum,ar1,ar2 real(rkind) :: signa !swild98 used for exchange (deallocate immediately afterwards) REAL(rkind),ALLOCATABLE :: swild98(:,:,:) !- Start Statement --------------------------------------------------! IF(myrank==0) write(16,*)'Entering sed_alloc' !--------------------------------------------------------------------! !* ARRAYS ALLOCATION (a few have been done in read_sed_input) !--------------------------------------------------------------------! !--------------------------------------------------------------! !* 1D arrays defined on elements !--------------------------------------------------------------! ALLOCATE(Zob(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: Zob allocation failure') ALLOCATE(dave(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: dave allocation failure') ALLOCATE(FX_r(nea), stat=i ) IF(i/=0) CALL parallel_abort('Sed: FX_r allocation failure') ALLOCATE(FY_r(nea), stat=i ) IF(i/=0) CALL parallel_abort('Sed: FY_r allocation failure') ALLOCATE(bustr(nea), stat=i ) IF(i/=0) CALL parallel_abort('Sed: bustr allocation failure') ALLOCATE(bvstr(nea), stat=i ) IF(i/=0) CALL parallel_abort('Sed: bvstr allocation failure') ALLOCATE(bed_thick(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: bed_thick allocation failure') ALLOCATE(hs(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: hs allocation failure') ALLOCATE(tp(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: tp allocation failure') ALLOCATE(wlpeak(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: wlpeak allocation failure') ALLOCATE(uorb(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: uorb allocation failure') ALLOCATE(uorbp(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: uorbp allocation failure') ALLOCATE(dirpeak(nea),stat=i) !Anouk IF(i/=0) CALL parallel_abort('Sed: dirpeak allocation failure') ALLOCATE(tau_m(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: tau_m allocation failure') ALLOCATE(tau_c(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: tau_c allocation failure') ALLOCATE(tau_w(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: tau_w allocation failure') ALLOCATE(tau_wc(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: tau_wc allocation failure') !BM ALLOCATE(wdir(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: wdir allocation failure') ALLOCATE(ubott(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: ubott allocation failure') ALLOCATE(vbott(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: vbott allocation failure') ALLOCATE(U_crest(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: U_crest allocation failure') ALLOCATE(U_trough(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: U_trough allocation failure') ALLOCATE(T_crest(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: T_crest allocation failure') ALLOCATE(T_trough(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: T_trough allocation failure') ALLOCATE(Qaccu(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: acceleration-induced bedload') ALLOCATE(Qaccv(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: acceleration-induced bedload') ALLOCATE(r_accu(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: acceleration-induced bedload') ALLOCATE(r_accv(nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: acceleration-induced bedload') !--------------------------------------------------------------! !* 2D arrays defined on elements !--------------------------------------------------------------! ALLOCATE(mcoefd(0:mnei_p,np),stat=i) IF (i/=0) CALL parallel_abort('Sed: mcoefd allocation failure') ALLOCATE(Hz(nvrt,nea),stat=i) IF(i/=0) CALL parallel_abort('Sed: Hz allocation failure') ALLOCATE(bottom(nea,MBOTP),stat=i) IF(i/=0) CALL parallel_abort('Sed: bottom allocation failure') ALLOCATE(sedcaty(nea,ntr_l),stat=i) !Tsinghua group:pick up sediment flux IF(i/=0) CALL parallel_abort('Sed: sedcaty allocation failure') ! BM ALLOCATE(eroflxel(nea,ntr_l),stat=i) IF(i/=0) CALL parallel_abort('Sed: erosion flux at elem') ALLOCATE(depflxel(nea,ntr_l),stat=i) IF(i/=0) CALL parallel_abort('Sed: deposition flux at elem') ALLOCATE(Uorbi_elfrink(nea,ech_uorb+1),stat=i) IF(i/=0) CALL parallel_abort('Sed: Uorbi_elfrink alloc failure') !--------------------------------------------------------------! !* 3D arrays defined on elements !--------------------------------------------------------------! ALLOCATE(bed(Nbed,nea,MBEDP),stat=i) IF(i/=0) CALL parallel_abort('Sed: bed allocation failure') ALLOCATE(bed_frac(Nbed,nea,ntr_l),stat=i) IF(i/=0) CALL parallel_abort('Sed: bed_frac allocation failure') !--------------------------------------------------------------! !* 4D arrays defined on elements !--------------------------------------------------------------! ALLOCATE(bed_mass(Nbed,nea,2,ntr_l),stat=i) IF(i/=0) CALL parallel_abort('Sed: bed_mass allocation failure') !--------------------------------------------------------------! !* 1D arrays defined on nodes !--------------------------------------------------------------! ALLOCATE(vc_area(npa),stat=i) IF(i/=0) CALL parallel_abort('sed_alloc: vc_area allocation failure') ALLOCATE(bedthick_overall(npa),bed_d50n(npa),stat=i) IF(i/=0) CALL parallel_abort('sed_alloc: bed_d50n allocation failure') ALLOCATE(bed_taun(npa),stat=i) IF(i/=0) CALL parallel_abort('sed_alloc: bed_taun allocation failure') ALLOCATE(bed_rough(npa),stat=i) IF(i/=0) CALL parallel_abort('sed_alloc: bed_rough allocation failure') ALLOCATE(lbc_sed(npa),stat=i) IF(i/=0) CALL parallel_abort('sed_alloc: lbc_sed allocation failure') ALLOCATE(bc_sed(npa),stat=i) IF(i/=0) CALL parallel_abort('sed_alloc: bc_sed allocation failure') ! BM ALLOCATE(poron(npa),stat=i) IF(i/=0) CALL parallel_abort('sed_alloc: porosity at node') ALLOCATE(eroflxn(npa),stat=i) IF(i/=0) CALL parallel_abort('sed_alloc:: erosion flux at node') ALLOCATE(depflxn(npa),stat=i) IF(i/=0) CALL parallel_abort('sed_alloc: deposition flux at node') ALLOCATE(Qaccun(npa),stat=i) IF(i/=0) CALL parallel_abort('sed_alloc: Qacc at node') ALLOCATE(Qaccvn(npa),stat=i) IF(i/=0) CALL parallel_abort('sed_alloc: Qacc at node') !morphological ramp value(-) ALLOCATE(imnp(npa),stat=i) IF(i/=0) CALL parallel_abort('sed_alloc: morphological ramp value') !--------------------------------------------------------------! !* 2D arrays defined on nodes !--------------------------------------------------------------! ALLOCATE(bedldu(npa,ntr_l),stat=i) IF(i/=0) CALL parallel_abort('sed_alloc: bedldu alloc failure') ALLOCATE(bedldv(npa,ntr_l),stat=i) IF(i/=0) CALL parallel_abort('sed_alloc: bedldv alloc failure') ALLOCATE(bed_fracn(npa,ntr_l),stat=i) IF(i/=0) CALL parallel_abort('sed_alloc: bed_fracn alloc failure') !--------------------------------------------------------------! !* Other 1D arrays !--------------------------------------------------------------! ALLOCATE ( isand(ntr_l), stat=i ) IF(i/=0) CALL parallel_abort('Main: allocation failure') !--------------------------------------------------------------------! !* Pre-initialization !--------------------------------------------------------------------! !--------------------------------------------------------------! !* 1D arrays defined on elements !--------------------------------------------------------------! Zob(:) = IniVal dave(:) = IniVal FX_r(:) = IniVal FY_r(:) = IniVal bustr(:) = IniVal bvstr(:) = IniVal bed_thick(:) = IniVal hs(:) = IniVal tp(:) = IniVal wlpeak(:) = IniVal uorb(:) = IniVal uorbp(:) = IniVal tau_m(:) = IniVal tau_c(:) = IniVal tau_w(:) = IniVal tau_wc(:) = IniVal !BM Qaccu(:)=IniVal Qaccv(:)=IniVal r_accu(:)=IniVal r_accv(:)=IniVal !--------------------------------------------------------------! !* 2D arrays defined on elements !--------------------------------------------------------------! mcoefd(:,:) = IniVal Hz(:,:) = IniVal bottom(:,:) = IniVal sedcaty(:,:)= IniVal !Tsinghua group !BM eroflxel(:,:)=IniVal depflxel(:,:)=IniVal !--------------------------------------------------------------! !* 3D arrays defined on elements !--------------------------------------------------------------! bed(:,:,:) = IniVal bed_frac(:,:,:) = IniVal !--------------------------------------------------------------! !* 3D arrays defined on elements !--------------------------------------------------------------! bed_mass(:,:,:,:) = IniVal !--------------------------------------------------------------! !* 1D arrays defined on nodes !--------------------------------------------------------------! vc_area(:) = IniVal bed_d50n(:) = IniVal bed_taun(:) = IniVal bed_rough(:) = IniVal ! BM poron(:) = IniVal eroflxn(:)=IniVal depflxn(:)=IniVal Qaccun(:)=IniVal Qaccvn(:)=IniVal !morphological ramp value(-) imnp(:)=1.d0 !Special initializations bc_sed(:) = -9999 lbc_sed(:) = .FALSE. !--------------------------------------------------------------! !* 2D arrays defined on nodes !--------------------------------------------------------------! bedldu(:,:) = IniVal bedldv(:,:) = IniVal bed_fracn(:,:) = IniVal !--------------------------------------------------------------! !* Other 1D arrays !--------------------------------------------------------------! isand(:) = 1 ! init. !--------------------------------------------------------------------! ! - Set tracer indices into 1:ntracers ! In SCHISM T and S are in different arrays from the tracers array... !--------------------------------------------------------------------! ic = irange_tr(1,5) DO i = 1,ntr_l isand(i) = ic ic = ic+1 END DO !--------------------------------------------------------------------! ! - Computes matrix coefficients for the JCG solver ! Used for the computation of depth variation induced by bedload !--------------------------------------------------------------------! mcoefd = 0.d0 aux1=22.0d0/108.0d0 aux2=7.0d0/108.0d0 DO i=1,np !residents DO j=1,nne(i) ie=indel(j,i) id=iself(j,i) if(i34(ie)==3) then mcoefd(0,i) = mcoefd(0,i)+area(ie)*aux1 !diagonal !Other 2 nodes do jj=1,2 !other 2 nodes nd=elnode(nxq(jj,id,i34(ie)),ie) indx=0 do m=1,nnp(i) if(indnd(m,i)==nd) then indx=m; exit endif enddo !m if(indx==0) call parallel_abort('SED_INIT: failed to find') mcoefd(indx,i)=mcoefd(indx,i)+area(ie)*aux2 enddo !jj ! IF(isbnd(1,i)==0.and.j==nne(i)) THEN !internal ball ! mcoefd(1,i) = mcoefd(1,i)+area(ie)*aux2 ! ELSE ! mcoefd(j+1,i) = mcoefd(j+1,i)+area(ie)*aux2 ! ENDIF else !quad id1=nxq(1,id,i34(ie)) id2=nxq(2,id,i34(ie)) id3=nxq(3,id,i34(ie)) if(ics==1) then ar1=signa(xnd(i),xcj(elside(id3,ie)),xctr(ie),ynd(i),ycj(elside(id3,ie)),yctr(ie)) ar2=signa(xnd(i),xctr(ie),xcj(elside(id2,ie)),ynd(i),yctr(ie),ycj(elside(id2,ie))) else !ll cff1=(xel(id,ie)+xel(id1,ie))/2.d0 !xcj(elside(id3,ie)) cff2=(yel(id,ie)+yel(id1,ie))/2.d0 !ycj(elside(id3,ie)) cff3=sum(xel(1:4,ie))/4.d0 !xctr cff4=sum(yel(1:4,ie))/4.d0 !yctr cff5=(xel(id,ie)+xel(id3,ie))/2.d0 !xcj(elside(id2,ie)) cff6=(yel(id,ie)+yel(id3,ie))/2.d0 !ycj(elside(id2,ie)) ar1=signa(xel(id,ie),cff1,cff3,yel(id,ie),cff2,cff4) ar2=signa(xel(id,ie),cff3,cff5,yel(id,ie),cff4,cff6) endif !ics if(ar1<=0.d0.or.ar2<=0.d0) call parallel_abort('SED_INIT:area<=0') mcoefd(0,i)=mcoefd(0,i)+(ar1+ar2)*7./12 !diagonal !Find indices do jj=1,3 nd=elnode(nxq(jj,id,i34(ie)),ie) indx=0 do m=1,nnp(i) if(indnd(m,i)==nd) then indx=m; exit endif enddo !m if(indx==0) call parallel_abort('SED_INIT: faile to find2') nwild(jj)=indx enddo !jj mcoefd(nwild(1),i)=mcoefd(nwild(1),i)+ar1/4+ar2/12 mcoefd(nwild(3),i)=mcoefd(nwild(3),i)+ar1/12+ar2/4 mcoefd(nwild(2),i)=mcoefd(nwild(2),i)+ar1/12+ar2/12 endif !i34 ENDDO ! END loop nne ENDDO ! END loop np !--------------------------------------------------------------------! ! - Control volume at each node used in filter ! Split quads into 2 tri's !--------------------------------------------------------------------! vc_area = 0.0d0 nwild2(1:3)=(/1,3,4/) !prep. indices for 2nd tri of quad DO i=1,nea if(i34(i)==4) then !2 areas !nwild(1:3)=elnode(1:3,i) !ar1=signa(xnd(nwild(1)),xnd(nwild(2)),xnd(nwild(3)),ynd(nwild(1)),ynd(nwild(2)),ynd(nwild(3))) ar1=signa(xel(1,i),xel(2,i),xel(3,i),yel(1,i),yel(2,i),yel(3,i)) !nwild(2)=elnode(3,i) !nwild(3)=elnode(4,i) !ar2=signa(xnd(nwild(1)),xnd(nwild(2)),xnd(nwild(3)),ynd(nwild(1)),ynd(nwild(2)),ynd(nwild(3))) ar2=signa(xel(1,i),xel(3,i),xel(4,i),yel(1,i),yel(3,i),yel(4,i)) if(ar1<=0.or.ar2<=0) call parallel_abort('SED_INIT:area2<=0') endif DO j=1,3 if(i34(i)==3) then nd=elnode(j,i) vc_area(nd)=vc_area(nd)+area(i)/3 else !quad !1st tri nd=elnode(j,i) vc_area(nd)=vc_area(nd)+ar1/3 !2nd tri nd=elnode(nwild2(j),i) vc_area(nd)=vc_area(nd)+ar2/3 endif !i34(i) ENDDO !j ENDDO ! i CALL exchange_p2d(vc_area) ! Read in total bed thickness at nodes OPEN(10,FILE=in_dir(1:len_in_dir)//'bedthick.ic',STATUS='OLD') READ(10,*); READ(10,*) DO i = 1,np_global READ(10,*) itmp,xtmp,ytmp,tmp1 IF(tmp1<0) CALL parallel_abort('Sed_init: bed_thick<0!') IF(ipgl(i)%rank==myrank) bedthick_overall(ipgl(i)%id)=tmp1 ENDDO !i=1,np_global CLOSE(10) ! BM: read morphological ramp value(-) in imorphogrid.gr3 IF (sed_morph==1) THEN INQUIRE(FILE='imorphogrid.gr3',EXIST=lexist) IF (lexist) THEN OPEN(10,FILE='imorphogrid.gr3',STATUS='OLD') READ(10,*) READ(10,*) j,itmp IF(itmp/=np_global) CALL parallel_abort('sed3d: Check np_global in imorphogrid.gr3') DO i = 1,np_global READ(10,*) itmp,xtmp,ytmp,tmp1 IF(tmp1<0.or.tmp1>1) CALL parallel_abort('sed3d: local imorpho <0 or >1!') IF(ipgl(i)%rank == myrank) imnp(ipgl(i)%id) = tmp1 ENDDO !np_global CLOSE(10) ELSE CALL parallel_abort('sed_morph=1 requires imorphogrid.gr3') ENDIF !lexist ENDIF !sed_morph ! For cold start only !if(ihot==0) then !======================================================== !--------------------------------------------------------------------! ! - Reading bed_frac files and convert bed_fraction from nodes to ! elements. ! For instance, bed_fraction is applied to all bed layers, i.e. no vertical variation !--------------------------------------------------------------------! ALLOCATE(swild98(Nbed,npa,ntr_l),stat=istat) ! * READING bed_frac_x.ic DO ised = 1,ntr_l WRITE(ifile_char,'(i03)')ised ifile_char=ADJUSTL(ifile_char); ifile_len=LEN_TRIM(ifile_char) inputfile='bed_frac_'//ifile_char(1:ifile_len)//'.ic' OPEN(10,FILE=in_dir(1:len_in_dir)//inputfile,STATUS='OLD') READ(10,*) !read in first line, no need to store it READ(10,*) !read in second line, no need to store it DO i = 1,np_global READ(10,*) itmp,xtmp,ytmp,tmp1 IF(tmp1<0.or.tmp1>1) CALL parallel_abort('Sed: bed_frac wrong!') IF(ipgl(i)%rank==myrank) swild98(:,ipgl(i)%id,ised) = tmp1 ENDDO !i=1,np_global CLOSE(10) ENDDO !ised=1,ntr_l !'-------------------------------------------------------------------! ! - Mapping bed fraction (conversion from node to elements) !--------------------------------------------------------------------! DO ised = 1,ntr_l DO k = 1,Nbed DO i = 1,nea bed_frac(k,i,ised) = sum(swild98(k,elnode(1:i34(i),i),ised))/i34(i) ENDDO ! END loop nea ENDDO ! END loop Nbed ENDDO ! END loop ntr_l DEALLOCATE(swild98,stat=istat) !--------------------------------------------------------------------! ! - Check the sum of bed fractions ! Sum should be equal to 1 but not exactly possible due to rounding ! - Sum must not exceed 1.01 ! - Sum must not be lower than 0.99 !--------------------------------------------------------------------! DO k=1,Nbed DO i=1,nea bed_frac_sum = 0.0d0 DO ised=1,ntr_l bed_frac_sum = bed_frac_sum+bed_frac(k,i,ised) ENDDO !End loop ntr_l IF (bed_frac_sum>1.01d0.or.bed_frac_sum<0.99d0) THEN WRITE(errmsg,*)'SED: sum of bed_frac/=1 at elem ',ielg(i),bed_frac_sum CALL parallel_abort(errmsg) ENDIF ENDDO !i ENDDO !k !--------------------------------------------------------------------! ! - Initialize the bed model (layer thickness, age and porosity) !--------------------------------------------------------------------! ! DO i=1,Nbed ! DO j=1,nea ! bed(i,j,ithck) = sum(bedthick_overall(elnode(1:i34(j),j)))/i34(j)/real(Nbed) !>=0 ! bed(i,j,iaged) = 0.0d0 ! IF (poro_option .EQ. 2) THEN ! porosity=f(geom std dev d50) ! CALL sed_comp_poro_noncoh(bed_frac(i,j,:),porosity) ! END IF ! bed(i,j,iporo) = porosity ! ENDDO ! End loop Nbed ! ENDDO !End loop nea ! > BM rewrite and introduce sedlay_ini_opt (see sediment.in) DO j=1,nea DO i=1,Nbed IF (sedlay_ini_opt==1 .AND. Nbed > 1) THEN IF (i==1) THEN bed(i,j,ithck) = toplay_inithick ELSE bed(i,j,ithck) = sum(bedthick_overall(elnode(1:i34(j),j))-toplay_inithick)/i34(j)/real(Nbed-1) !>=0 END IF ELSE bed(i,j,ithck) = sum(bedthick_overall(elnode(1:i34(j),j)))/i34(j)/real(Nbed) !>=0 END IF bed(i,j,iaged) = 0.0d0 IF (poro_option .EQ. 2) THEN ! porosity=f(geom std dev d50) CALL sed_comp_poro_noncoh(bed_frac(i,j,:),porosity) END IF bed(i,j,iporo) = porosity ENDDO ! End loop Nbed ENDDO !End loop nea !======================================================== ! endif !ireset IF(myrank==0) write(16,*)'Leaving sed_alloc' !--------------------------------------------------------------------! END SUBROUTINE sed_alloc !===================================================================== !===================================================================== SUBROUTINE sed_init !--------------------------------------------------------------------! ! This subroutine initialize variables and arrays for the sediment ! ! model: ! ! - computes coefficients for the jcg solver ! ! - computes control volume area at nodes ! ! - read bed fraction file (bed_frac_x.ic) ! ! - mapping of bed fraction ! ! - checking sum of bed fraction ! ! - initialize the bed model (layer thickness, age, porosity) ! ! - initialize the bed mass ! ! - initialize exposed sediment layer properties ! ! - initialize sediment roughness length ! ! - initialize total bed thickness ! ! - initialization of arrays defined at nodes ! ! - if debuging: write sediment initialization ! ! ! ! Adapted from former subroutines initialize_scalars, ! ! initialize_ocean (former init_sed.F90), and sed_init (former ! ! sed_init.F90) and other initializations formerly done within ! !schism_init.F90 were also moved within this subroutine. ! ! ! ! The former subroutine sed_init.F90 was adapted from ROMS routine ! ! ana_sediment.h ! ! Copyright (c) 2002-2007 The ROMS/TOMS Group ! ! Licensed under a MIT/X style license ! ! See License_ROMS.txt ! ! ! ! Author: florian ganthy (fganthy@lnec.pt ; florian.ganthy@gmail.com)! ! Date: 2013/06/07 ! ! ! ! History: ! ! *** Previous history (former subroutines) *** ! ! ! ! 2012/12 - F.Ganthy : form homogenisation of sediments routines ! ! 2012/12 - F.Ganthy : modifications for Bottom Composition ! ! Generation (BCG) purpose (added output ! ! files - bed median grain size, bed fraction)! ! 2013/01 - F.Ganthy : Implementation of roughness predictor ! ! 2013/03 - F.Ganthy : Implementation of wave-induced bedload ! ! transport ! ! 2013/04 - F.Ganthy : Implementation of wave-current bottom stress! ! 2013/05 - F.Ganthy : Updates related to ripple predictor ! ! 2013/05 - F.Ganthy : Added node-centered volume control area ! ! 2013/05 - F.Ganthy : Updates to the ripple predictor: ! ! - Changes on the total bedform ! ! roughness computation ! ! - Add wave-ripple computation from ! ! Nielsen (1992) ! ! ! ! *** Current history (former subroutines) *** ! ! 2013/06 - F.Ganthy : Also removed old initializations currently ! ! which can be found in former svn revision ! ! ! ! 2020/02 - B.Mengual : Active layer thickness initialized as the ! ! the minimum between the top layer one and ! ! actv_max (max thickness input parameter) ! ! ! !--------------------------------------------------------------------! USE sed_mod USE schism_glbl, ONLY: nea,npa,mnei_p,ntrs,irange_tr,ipgl,ielg,i34,elnode,np_global, & & ifile_char,ifile_len,area,np,nne,indel,iself,nnp,indnd,nxq, & & isbnd,rough_p,errmsg,xnd,ynd,xcj,ycj,xctr,yctr,elside,ics,xel,yel, & &in_dir,out_dir,len_in_dir,len_out_dir USE schism_msgp, ONLY: myrank,parallel_abort,exchange_p2d IMPLICIT NONE real(rkind) :: signa !- Local variables --------------------------------------------------! INTEGER :: i,j,k,ie,jj,id,id1,id2,id3,nd,indx,m,nwild(3),nwild2(3) INTEGER :: ised,ic,itmp,istat REAL(rkind) :: aux1,aux2,xtmp,ytmp,tmp1,cff1,cff2,cff3,cff4,cff5,cff6 REAL(rkind) :: bed_frac_sum,ar1,ar2 REAL(rkind),DIMENSION(npa) :: bdfc CHARACTER(len=48) :: inputfile !- Start Statement --------------------------------------------------! IF(myrank==0) WRITE(16,*)'Entering sed_init' !--------------------------------------------------------------------! ! - Initialize the bed mass ([kg/m2]=[m]*[kg.m-3]*[-]) !--------------------------------------------------------------------! DO k=1,Nbed DO j=1,nea DO ised=1,ntr_l bed_mass(k,j,:,ised) = bed(k,j,ithck)* & & Srho(ised)* & & (1.0d0-bed(k,j,iporo))* & & bed_frac(k,j,ised) ENDDO ! End loop ntr_l ENDDO ! End loop nea END DO ! End loop Nbed !--------------------------------------------------------------------! ! - Initialize exposed sediment layer properties: ! - total D50 (m) ! - total density (kg.m-3) ! - total settling velocity (m.s-1) ! - total critical shear stress (m2.s-2) !--------------------------------------------------------------------! DO j=1,nea cff1 = 1.0d0 cff2 = 1.0d0 cff3 = 1.0d0 cff4 = 1.0d0 cff5 = 0.0d0 ! Weighted geometric mean DO ised=1,ntr_l cff1 = cff1*Sd50(ised)**bed_frac(1,j,ised) cff2 = cff2*Srho(ised)**bed_frac(1,j,ised) cff3 = cff3*Wsed(ised)**bed_frac(1,j,ised) cff4 = cff4*tau_ce(ised)**bed_frac(1,j,ised) cff5 = cff5+bed_frac(1,j,ised) ENDDO ! End loop ntr_l if(cff5<0) then call parallel_abort('SED_INIT: cff5<0 (1)') else if(cff5==0) then !care-takers for all-eroded case IF (sed_debug .EQ. 1) THEN WRITE(12,*)'SED_INIT: all eroded at elem. ',ielg(j) END IF bottom(j,isd50) = Sd50(1) bottom(j,idens) = Srho(1) bottom(j,iwsed) = Wsed(1) bottom(j,itauc) = tau_ce(1) else bottom(j,isd50) = cff1**(1.0d0/cff5) bottom(j,idens) = cff2**(1.0d0/cff5) bottom(j,iwsed) = cff3**(1.0d0/cff5) bottom(j,itauc) = cff4**(1.0d0/cff5) endif !cff5 ENDDO ! j !Init. active layer thickness !bottom(:,iactv)=bed(1,:,ithck) bottom(:,iactv)=MIN(actv_max,bed(1,:,ithck)) ! BM !--------------------------------------------------------------------! ! - Initialize sediment main roughness length and bedform predictor ! related roughness length ! Here, the initialization is the same whether or not the bedfrom ! predictor is used !--------------------------------------------------------------------! ! Current ripple roughness length (Soulsby, 1997) bottom(:,izcr) = 0.0d0 ! Sand waves roughness length (Van Rijn, 1984) bottom(:,izsw) = 0.0d0 ! Wave ripples roughness length (Grant and Madsen, 1982 OR Nielsen, 1992) bottom(:,izwr) = 0.0d0 ! Bed load bottom roughness length (Grant and Madsen, 1982 OR Nielsen, 1992) bottom(:,izbld) = 0.0d0 DO i=1,nea ! Nikurasde roughness length bottom(i,izNik) = bottom(i,isd50)/12.0d0 ! Default roughness bottom(i,izdef) = sum(rough_p(elnode(1:i34(i),i)))/i34(i) ! Apparent initial roughness bottom(i,izapp) = bottom(i,izdef) ! Roughness length effectively used (even if bedform predictor is not used) Zob(i) = bottom(i,izdef) IF(Zob(i).LE.0.0d0)THEN WRITE(errmsg,*) 'SED_INIT: Zob <=0 at elem, rank:', Zob(i),i,myrank CALL parallel_abort(errmsg) ENDIF ENDDO ! i !--------------------------------------------------------------------! ! - Initialize total bed thickness by summing over all layers ! This array is not actively updated after this !--------------------------------------------------------------------! bed_thick(:)=0.0d0 DO i=1,Nbed DO j=1,nea bed_thick(j) = bed_thick(j)+bed(i,j,ithck) ENDDO !j ENDDO !i !--------------------------------------------------------------------! ! - Initialization of variables defined at nodes !--------------------------------------------------------------------! bed_d50n(:) = 0.0d0 bed_taun(:) = 0.0d0 bed_fracn(:,:) = 0.0d0 bdfc(:) = 0.0d0 DO i=1,nea DO j=1,i34(i) DO ised=1,ntr_l bed_fracn(elnode(j,i),ised) = bed_fracn(elnode(j,i),ised)+ & & bed_frac(1,i,ised) ENDDO !ised bed_d50n(elnode(j,i)) = bed_d50n(elnode(j,i))+bottom(i,isd50) bdfc(elnode(j,i)) = bdfc(elnode(j,i))+1.0d0 ENDDO !j ENDDO !i DO i=1,np if(bdfc(i)==0) call parallel_abort('SED_INIT: bdfc(i)==0') DO ised=1,ntr_l bed_fracn(i,ised) = bed_fracn(i,ised)/bdfc(i) ENDDO !ised bed_d50n(i) = bed_d50n(i)/bdfc(i) bed_rough(i) = rough_p(i) !roughness length ENDDO !i DO ised=1,ntr_l CALL exchange_p2d(bed_fracn(:,ised)) ENDDO !END loop ntr_l CALL exchange_p2d(bed_d50n) CALL exchange_p2d(bed_rough) !--------------------------------------------------------------------! ! - Sediment debug outputs !--------------------------------------------------------------------! ! IF(sed_debug==1) THEN DO i=1,Nbed DO j=1,nea DO k=1,ntr_l WRITE(12,*)'Nbed:',i,' nea:',j,' ntr_l:',k, & & ' bed_frac.ic:',real(bed_frac(i,j,k)) ENDDO !k ENDDO !j ENDDO !i DO i=1,Nbed DO j=1,nea WRITE(12,*)'Nbed:',i,' nea:',j,' bed_thick:',real(bed_thick(j)), & & ' bed(ithck):',real(bed(i,j,ithck)),' bed(iaged):', & & real(bed(i,j,iaged)),' bed(iporo):',real(bed(i,j,iporo)) ENDDO !j ENDDO !i DO i=1,nea WRITE(12,*)'nea:',i,' bottom(itauc):',real(bottom(i,itauc)), & & ' bottom(isd50):',real(bottom(i,isd50)), & & ' bottom(iwsed):',real(bottom(i,iwsed)), & & ' bottom(idens):',real(bottom(i,idens)) ENDDO !i ! ENDIF !End sed_debug !--------------------------------------------------------------------! IF(myrank==0) write(16,*)'Leaving sed_init' !--------------------------------------------------------------------! END SUBROUTINE sed_init