C ********************************************************************** C * * C * SOFTWARE LICENSING * C * * C * This program is free software; you can redistribute * C * it and/or modify it under the terms of the GNU * C * General Public License as published by the Free * C * Software Foundation, either Version 2 of the * C * license, or (at your option) any later version. * C * * C * This program is distributed in the hope that it * C * will be useful, but without any warranty; without * C * even the implied warranty of merchantability or * C * fitness for a particular purpose. See the GNU * C * General Public License for more details. * C * * C * A copy of the GNU General Public License is * C * available at http://www.gnu.org/copyleft/gpl.html * C * or by writing to the Free Software Foundation, Inc.,* C * 59 Temple Place - Suite 330, Boston, MA 02111, USA. * C * * C ********************************************************************** C----------------------------------------------------------------- C BEGINNING OF WAVE CODE !WAVE C----------------------------------------------------------------- subroutine wavemdo(dtw,dth,nwavemel) C------------------------------------------------------------ C Subroutine to calculate wave energy density as a function of C x, y and theta. Grid is orthogonal curvilinear. For rectilinear C grid, simply set dx(i,j) or dy(i,j) = constant C C Input: C z(k),dz(k) = vertical coordinate and spacing (non-d.) C zz(k),dzz(k) = middle cell coordinate and spacing C dx(i,j), dy(i,j) = spatial grid cell increments located ati C cell center (m) C h(i,j) = bottom depth C d(i,j) = water column depth (m) C fsm(i,j) = 1 for water cell, = 0 for land cell C u10x(i,j),u10y(i,j) = wind vector at 10 m height (m s-1) C u(i,j,k),v(i,j,k) = current vector (m s-1) c nitera = # of iterations in s.r. stepmpd C Output wave variables: C enb(i,j,m), en(i,j,m), enf(i,j,m) = past, present,future C energy densities (m+3 s-2) C ent(i,j) en averaged over all wave directions (m+3 s-2) C Hs(i,j) = significant wave height = 4*sqrt(ent/grav) (m) C thtav(i,j) = average wave propagation direction (rad) C Sxx,Sxy,Syy = wave radiation stress for momentum equation C related to wave velocity moments (m+3 s-2) C (tpx0,tpy0) = wave pressure stress at surface (m+2 s-2) C (tpx,tpy)=(tpx0,tpy0)*tpzdist = wave pressure below the C surface (m+2 s-2) C Internal wave variables (some could also be output) C theta(m) = wave directional coordinate located at cell edge (rad) C cs(m),sn(m) = cosine and sine of theta C thetam(m) = wave directional coordinate located at cell center C csm(m),snm(m) = cosine and sine of thetam C sigth(i,j,m)=theta dependent frequency, swell and/or wind C waves (s-1) C sigthav(i,j)=frequency, averaged on en(i,j,m) (s-1) C cth(i,j,m)=theta dependent phase speed (m s-1) C kth(i,j,m)= theta dependent wave number (m-1) C kthav(i,j)=wave number, averaged on en(i,j,m) (s-1) C kthavD(i,j)= average wave number *d(i,j) C cg(i,j,m) = spectrally averaged group speeds located at C cell center (m s-1) C sigp(i,j)=peak frequency of wind-driven waves (s-1) C cp(i,j) = peak frequency phase speed from sigp (m s-1) C kp(i,j) = peak frequency wave number from sigp (m-1) C kpD(i,j) = kp(i,j)*d(i,j) (non-d) C cgx(i,j,m), cgy(i,j,m) = propagation vector in x, y C directions located at cell edge (m s-1) C ud(i,j),vd(i,j) = Dopler velocities (m s-1) C cthth(i,j,m) = propagation speed in the theta direction C located at cell edge (rad s-1) C u10(i,j) = wind speed C ageinv(i,j)=u10*sigp/grav (= 0.83 for fully developed C waves) (non-d) C fspr(i,j)=spreading function for swin (non-d) C beta = constant in fspr; higher beta = narrower spread C but recommended value = 2,2 (non-d) C swin(i,j,m) = wave energy source term (m3 s-3) C sdis(i,j,m) = surface wave energy dissipation (m3 s-3) C bdis1(i,j,m) = bottom friction dissipation (m3 s-3) C bdis2(i,j,m) = bottom depth-induced dissipation (m3 s-3) C gama = empirical constant in formulation of bdis2 (non-d) C Sradx(i,j),Srady(i,j),Sradxy(i,j),Srad(i,j) = wave radiation C stress for wave energy equation related to wave velocity C moments (s-1) C Fcg(i,jm), Fct(i,j,m) = correction for cg and cth to C account for spectral averaging (non-d) C Cin, swcon = constants in swin curve fit (non-d) C adis, bdis = empirically determined constants in sdis term C (non-d) C swfct = factor to reduce dissipation for swells (non-d) C dif = horizontal advective diffusion coeficient (see note C in x-y step) (non-d) C------------------------------------------------------------ C CHANGES MADE BY RMarsooli_AUGUST2015 (COUPLING WITH sECOM) C 1- theta is named thetaa to avoid conflict with variables used in sECOM C 2- u10 is named u10wind to avoid conflict with variables used in sECOM C 3- u10x and u10y have been replaced by WSX2D and WSY2D & C which are x- and y-component of wind speed used in sECOM C 4- cs&sn are named css&snn to avoid conflict with variables used in sECOM C 5- mm is named mmm to avoid conflict with variables used in sECOM C 6- sdis is named sdiss to avoid conflict with variables used in sECOM C 7- cd is named cdd to avoid conflict with variables used in sECOM C 8- cp&kp are named cpp&kpp to avoid conflict with variables used in sECOM C 9- subroutine stress is named stress_mel to avoid conflict with sECOM C 10- subroutine speeds is named speeds_mel to avoid conflict with sECOM C 11- fsm is named fsmwave C------------------------------------------------------------- CRM implicit none C CRM include 'pom08.c' INCLUDE 'comdeck' !RMarsooli SAVE !RMarsooli C integer & nitera,mbig,mploc,mmloc,loop, & imax,jmax,init real & dxinv,dyinv,dx2inv,dy2inv,Hsmax, & diff,summ,fsprm,tht_ops,fb,Tp,bb,Hsbc, & cgxm,cgym,sorc,adv,sor,cgnd,fcoshinv,fsinhinv, & snmav,csmav,Tsig,Cf,Dcoeff,cdmax real ! 2-d variables & usdif(im,jm),vsdif(im,jm), & Sradx(im,jm),Sradxy(im,jm),Srady(im,jm),Srad(im,jm), & Srad0(im,jm) real ! 3-d variables & xflux(im,jm,mmm),yflux(im,jm,mmm),tflux(im,jm,mmm), & uak(im,jm,mmm), & Fcg(im,jm,mmm),Fct(im,jm,mmm),tpsw(im,jm,mmm), & tpsw2(im,jm,mmm),tpsw3(im,jm,mmm),tpsw4(im,jm,mmm), & tpsw5(im,jm,mmm), & westang,eastang,southang,northang real smothh,beta,Cin,swcon,adis,bdis,swfct,dif,gama REAL u10wind(im,jm),eninw(mmm),sigthinw(mmm) !RMarsooli real ucouple(im,jm,kb),vcouple(im,jm,kb) !RMarsooli DIMENSION COM(80) !RMarsooli, AUGUST2015 save dtw2 data smothh/.20/ !.20/ data beta/2.2/ ! spreading const. data Cin/370./ ! Sin const. data swcon/0.33/ ! Sin const. data adis/0.8/ !0.85/ !0.925/ ! Sdis const. data bdis/0.36e-4/ !0.36e-4/ !/0.18e-4/ ! Sdis const. data swfct/0.2/ ! factor to decrease swell diss. data nitera/3/ ! # of iterations in MPDATA data dif/0.5/ ! diffusion coeficient ! in x-y step section data gama/0.7/ ! depth-induced breaking const. data cdmax/.0030/ ! max. drag coefficient data initwave/0/ iint=INTX !RMarsooli C WAVE CELLS, RMarsooli_AUGUST2015 DO I=1,IM DO J=1,JM FSMW(I,J)=FSM(I,J)*DHMWAD(I,J) cc . *DUM(I,J)*DUM(I+1,J)*DVM(I,J)*DVM(I,J+1) ENDDO ENDDO DO N=1,NUMQBC IC=IQC(N) JC=JQC(N) FSMW(IC,JC)=0.0 ENDDO cc DO N=1,NUMEBC cc FSMW(IETA(N),JETA(N))=0.0 cc ENDDO DO N=1,NUMECR FSMW(ICHARD(N),JCHARD(N))=0.0 ENDDO C C------------------------------------------------------------------ if(iint.eq.1) then C------------------------------------------------------------------ rewind 10 open(11, file='specavs',form='formatted') dtw2=2.*dtw dthinv=1./dth ! inverse angle increment crm write(6,'(''wave parameters'')') crm write(6,'('' Cin ='',f8.1,'' adis ='',f8.3,'' bdis ='',e12.4)') crm & Cin,adis,bdis crm write(6,'('' beta ='',f10.3)') beta crm write(6,'('' nitera ='',i5 )') nitera crm write(6,'('' dif ='',f8.3)') dif crm write(6,'('' cdmax ='',f8.3)') cdmax do j=1,jm do i=1,im-1 CRM aren(i,j)=dx(i,j)*dy(i,j) aren(i,j)=H1(i,j)*H2(i,j) !RMarsooli areninv(i,j)=1./aren(i,j) enddo enddo write(6,'(5x,'' m theta cs sn thetam'', &'' csm snm'')') do m=1,mmm thetam(m)=dth*float(m-mmm/2) thetaa(m)=dth*float(m-mmm/2)-0.5*dth css(m)=cos(thetaa(m)) snn(m)=sin(thetaa(m)) csm(m)=cos(thetam(m)) snm(m)=sin(thetam(m)) write (6,'(5x,i5,6f10.3)') m, thetaa(m),css(m),snn(m), & thetam(m),csm(m),snm(m) enddo C C Open boundary conditions. (For a closed basin, not used) C Add north conditions as needed. crm Hsbc=0.2 !Sig. wave height at boundaries crm Tp=10 !wave period at boundaries crm westang=0. !western wave propagation angle; incoming crmC westang=pi !western wave propagation angle; outgoing crm eastang=0. !eastern wave propagation angle; outgoing crmC eastang=pi !eastern wave propagation angle; incoming crm southang=pi/2. !southern wave propagation angle; incoming crm do m=1,mm crm do j=1,jm crm enw(j,m)=grav*(Hsbc/4.)**2*fsprm(beta,m,mm,westang,thetam(m)) crm enw(j,m)=max(.0001,enw(j,m))*fsm(1,j) crm ene(j,m)=grav*(Hsbc/4.)**2*fsprm(beta,m,mm,eastang,thetam(m)) crm ene(j,m)=max(.0001,ene(j,m))*fsm(im,j) crm sigthe(j,m)=2.*pi/Tp crm sigthw(j,m)=2.*pi/Tp crm enddo crm enddo crm do m=1,mm crm do i=1,im crm ens(i,m)=grav*(1.00/4.)**2*fsprm(beta,m,mm,southang,thetam(m)) crm ens(i,m)=max(.0001,ens(i,m))*fsm(i,1) crm sigths(i,m)=2.*pi/Tp crm enddo crm enddo endif C-------------------------Initialize----------------------------- if(initwave.eq.0) then C---------------------------------------------------------------- do k=1,kb F1(k)=0. F2(k)=0. F3(k)=0. F4(k)=0. enddo do j=1,jm do i=1,im ud(i,j)=0. vd(i,j)=0. sdistot(i,j)=0. bdistot(i,j)=0. swintot(i,j)=0. thtav(i,j)=0. qdiss(i,j)=0.5 do m=1,mmm cth(i,j,m)=10. kth(i,j,m)=.1 kthD(i,j,m)=1. Fcg(i,j,m)=0.5 Fct(i,j,m)=1. bdis1(i,j,m)=0. bdis2(i,j,m)=0. enddo ageinv(i,j)=5. crm sigp(i,j)=1.5 sigp(i,j)=2.*pi/WPINITIAL !RMarsooli ustp(i,j)=0. cpp(i,j)=grav/sigp(i,j) kpp(i,j)=sigp(i,j)/cpp(i,j) kpD(i,j)=kpp(i,j)*d(i,j) kthavD(i,j)=kthD(i,j,1) Sradx(i,j)=0. Sradxy(i,j)=0. Srady(i,j)=0. Srad(i,j)=0. ud(i,j)=0. vd(i,j)=0. tht_wnd(i,j)=0. crm do k=1,kb crm Sxx(i,j,k)=0. crm Sxy(i,j,k)=0. crm Syy(i,j,k)=0. crm enddo enddo enddo do m=1,mmm do j=1,jm do i=1,im sigthb(i,j,m)=sigp(i,j) sigth(i,j,m)=sigthb(i,j,m) sigthf(i,j,m)=sigth(i,j,m) cth(i,j,m)=grav/sigth(i,j,m) kth(i,j,m)=sigth(i,j,m)**2/grav kthD(i,j,m)=kth(i,j,m)*d(i,j) C Set up intitial enb, en, enf crm Hs(i,j)=0.05 Hs(i,j)=WHINITIAL !RMarsooli enb(i,j,m)=grav*(Hs(i,j)/4.)**2 & *fsprm(beta,m,mmm,tht_wnd(i,j),thetam(m)) enb(i,j,m)=max(.000001,enb(i,j,m))*FSMW(i,j) en(i,j,m)=enb(i,j,m) enf(i,j,m)=en(i,j,m) C xflux(i,j,m)=0. yflux(i,j,m)=0. tflux(i,j,m)=0. fspr(i,j,m)=fsprm(beta,m,mmm,tht_wnd(i,j),thetam(m)) enddo enddo enddo do j=1,jm do i=1,im ent(i,j)=summ(enb,dth,i,j,im,jm,mmm) enddo enddo C C CRM call prxy('tht_wnd',0.,tht_wnd,im,iskp,jm,jskp,0.) CRM call prxy('Hs ',0.,Hs ,im,iskp,jm,jskp,0.) CRM call prxz('fspr',0.,fspr,im,iskp,jm,mm,2,jm/2,jm-2,0.,d,thetam) c call prxz('enf ',0.,enf ,im,iskp,jm,mm,2,100,199,0.,d,thetam) c call prxy('ent',0.,ent,im,iskp,jm,jskp,0.) C---------------------------------------------------------------- initwave=1 endif C---------------------------------------------------------------- call wave2cur (dth) !WAVE radiation: C dtw2=2.*dtw do j=2,jm-1 do i=2,im-1 CRM u10(i,j)=sqrt(u10x(i,j)**2+u10y(i,j)**2) CRM usdif(i,j)=(u10x(i,j)-u(i,j,1))*fsm(i,j) ! Used in s.r. stress CRM vsdif(i,j)=(u10y(i,j)-v(i,j,1))*fsm(i,j) c c RMarsooli IF(TOR.EQ.'BAROTROPIC') THEN UCURRENT0=0.5*(UA(I,J)+UA(I+1,J)) VCURRENT0=0.5*(VA(I,J)+VA(I,J+1)) ELSE UCURRENT0=0.5*(U(I,J,1)+U(I+1,J,1)) VCURRENT0=0.5*(V(I,J,1)+V(I,J+1,1)) ENDIF c Current velocity w.r.t true west-east and south-north axis UCURRENT=UCURRENT0*COSANG(I,J)-VCURRENT0*SINANG(I,J) VCURRENT=UCURRENT0*SINANG(I,J)+VCURRENT0*COSANG(I,J) u10wind(i,j)=SQRT(WSX2D(I,J)**2+WSY2D(I,J)**2) usdif(i,j)=(WSX2D(I,J)-UCURRENT)*FSMW(i,j) vsdif(i,j)=(WSY2D(I,J)-VCURRENT)*FSMW(i,j) ccc UCURRENT=UCURRENT0 !VELOCITY COMPONENT IN ZETA1 DIRECTION ccc VCURRENT=VCURRENT0 !VELOCITY COMPONENT IN ZETA2 DIRECTION ccc uwinds=WSX2D(i,j)*COSANG(i,j)+WSY2D(i,j)*SINANG(i,j) !WIND w.r.t ZETA1,TOWARD ccc vwinds=-WSX2D(i,j)*SINANG(i,j)+WSY2D(i,j)*COSANG(i,j) !WIND w.r.t ZETA2,TOWARD ccc u10wind(i,j)=SQRT(uwinds**2+vwinds**2) ccc usdif(i,j)=(uwinds-UCURRENT)*FSMW(i,j) ccc vsdif(i,j)=(vwinds-VCURRENT)*FSMW(i,j) cc if(abs(uwinds).lt.0.01) usdif(i,j)=0.0 cc if(abs(vwinds).lt.0.01) vsdif(i,j)=0.0 CRM if(u10(i,j).gt.0.0001) CRM & tht_wnd(i,j)=acos(u10x(i,j)/(u10(i,j)))*sign(1.,u10y(i,j)) if(u10wind(i,j).gt.0.0001) then C Wind direction w.r.t. true east (x direction) in rad, counterclockwise is positive tht_wnd(i,j)=acos(WSX2D(I,J)/(u10wind(i,j))) . *sign(1.,WSY2D(I,J)) c wind direction w.r.t. ZETA1, toward, c counterclockwise:+, clockwise:- ccc tht_wnd(i,j)=acos(uwinds/(u10wind(i,j))) ccc . *sign(1.,vwinds) endif CRM if(u10 (i,j).gt.4.00 ) then if(u10wind(i,j).gt.4.00 ) then do m=1,mmm fspr(i,j,m)=fsprm(beta,m,mmm,tht_wnd(i,j),thetam(m)) & *FSMW(i,j) enddo else do m=1,mmm fspr(i,j,m)=0. enddo endif enddo enddo call stress_mel(iint,nwavemel,Hs,ageinv,usdif,vsdif,cdd,tpx0,tpy0, & cdmax,ttx0,tty0,ustw,ustp,im,jm . ,FSMW,SIG,WN) !FSMW, SIG, and WN are added by RMarsooli C Calculate sigp for wind driven portion of enf do j=1,jm do i=1,im c ageinv(i,j)=1. if (FSMW(i,j).gt.0.) then CRM if(u10(i,j).gt.4.0) then call wvfcns(sigp(i,j),d(i,j),kpp(i,j),cpp(i,j)) CRM ageinv(i,j)=u10(i,j)*sigp(i,j)/grav ageinv(i,j)=u10wind(i,j)*sigp(i,j)/grav !RMarsooli if(u10wind(i,j).gt.4.0) then enwnd(i,j)= 0.000001 do m=2,mmm-1 if(fspr(i,j,m).gt.0.10) & enwnd(i,j)=enwnd(i,j)+en(i,j,m)*dth enddo Dcoeff=.0022*(0.7+.12*sqrt(enwnd(i,j))) CRM sigp(i,j) =(grav/u10(i,j)) CRM & *(Dcoeff*u10(i,j)**4/(grav*enwnd(i,j)))**.303 sigp(i,j) =(grav/u10wind(i,j)) !RMarsooli & *(Dcoeff*u10wind(i,j)**4/(grav*enwnd(i,j)))**.303 call wvfcns(sigp(i,j),d(i,j),kpp(i,j),cpp(i,j)) CRM ageinv(i,j)=u10(i,j)*sigp(i,j)/grav ageinv(i,j)=u10wind(i,j)*sigp(i,j)/grav !RMarsooli endif endif enddo enddo c call prxz('enf ',time,enf ,im,iskp,jm,mmm,2,111,jm-2, 0.,d,thetam) c call prxy('ent',time,ent,im,iskp,jm,jskp ,0. ) c call prxy('sigp',time,sigp,im,iskp,jm,jskp ,0. ) c call prxy('kp ',time,kp ,im,iskp,jm,jskp ,0. ) c call prxy('cpp ',time,cpp ,im,iskp,jm,jskp ,0. ) c call prxy('cg ',time,cg ,im,iskp,jm,jskp,0. ) C C---------------- Calculate propagation speeds----------------------- do j=1,jm do i=1,im if(FSMW(i,j).gt.0.) then ent(i,j)=summ(en,dth,i,j,im,jm,mmm) if(ent(i,j).lt.0.) then write(6,'(''entlt0'',4i5,e12.3)') iint,nwavemel,i,j,ent(i,j) stop endif Hs(i,j)=0.5*(4.*sqrt(ent(i,j)/grav)+Hs(i,j)) C Find kth and cth do m=2,mmm-1 call wvfcns(sigthb(i,j,m),d(i,j),kth(i,j,m),cth(i,j,m)) kthD(i,j,m)=kth(i,j,m)*d(i,j) C Find spectrally averaged corrections to cg and cthth C cg used in x-y step; cthth used in theta step call speeds_mel(kthD(i,j,m),Fcg(i,j,m),Fct(i,j,m)) cg(i,j,m)=Fcg(i,j,m)*cth(i,j,m) enddo endif enddo enddo C C------------------------------------------------------------- do j=2,jm-1 do i=2,im-1 if(FSMW(i,j).eq.0.) then tps(i,j)=(FSMW(i+1,j)+FSMW(i-1,j)+FSMW(i,j+1) . +FSMW(i,j-1)) do m=1,mmm cg(i,j,m)= & FSMW(i+1,j)*cg(i+1,j,m)+FSMW(i-1,j)*cg(i-1,j,m) & +FSMW(i,j+1)*cg(i,j+1,m)+FSMW(i,j-1)*cg(i,j-1,m) enddo if(tps(i,j).gt.0) then do m=1,mmm cg(i,j,m)=cg(i,j,m)/tps(i,j) enddo endif endif enddo enddo C--------------------------------------------------------------- C The solution is split between x-y steps, theta steps and C source steps. C --------------------------------------------------------------- C x-y step for wave energy C --------------------------------------------------------------- C If dif=0.2, then diff = 0.2 everywhere, but, on open water C cells imediately next to land cells, diff = 0.5 for local upwind C B.C. If dif=0.5, then diff=0.5 everywhere. C N.B. not a great difference if diff = 0.5 everywhere. do m=2,mmm-1 do j=2,jm do i=2,im cc cgx(i,j,m)=0.5*(cg(i,j,m)+cg(i-1,j,m))*csm(m) cc & +0.5*(ud(i,j)+ud(i-1,j)) cc cgy(i,j,m)=0.5*(cg(i,j,m)+cg(i,j-1,m))*snm(m) cc & +0.5*(vd(i,j)+vd(i,j-1)) c !RMarsooli cgx_temp=0.5*(cg(i,j,m)+cg(i-1,j,m))*csm(m) cgy_temp=0.5*(cg(i,j,m)+cg(i,j-1,m))*snm(m) cgx(i,j,m)=cgx_temp*COSANG(i,j)+cgy_temp*SINANG(i,j) !cgx w.r.t ZETA1 . +0.5*(ud(i,j)+ud(i-1,j)) cgy(i,j,m)=-cgx_temp*SINANG(i,j)+cgy_temp*COSANG(i,j) !cgy w.r.t ZETA2 . +0.5*(vd(i,j)+vd(i,j-1)) diff=0.5*(1.-dum(i,j))+dif*dum(i,j) xflux(i,j,m)= 0.5*cgx(i,j,m)*(enb(i,j,m)+enb(i-1,j,m)) & -diff*abs(cgx(i,j,m))*(enb(i,j,m)-enb(i-1,j,m)) diff=0.5*(1.-dvm(i,j))+dif*dvm(i,j) yflux(i,j,m)= 0.5*cgy(i,j,m)*(enb(i,j,m)+enb(i,j-1,m)) & -diff*abs(cgy(i,j,m))*(enb(i,j,m)-enb(i,j-1,m)) CRM xflux(i,j,m)=0.5*(dy(i,j)+dy(i-1,j))*xflux(i,j,m) !*dum(i,j) CRM yflux(i,j,m)=0.5*(dx(i,j)+dx(i,j-1))*yflux(i,j,m) !*dvm(i,j) xflux(i,j,m)=0.5*(H2(i,j)+H2(i-1,j))*xflux(i,j,m) !*dum(i,j) !RMarsooli yflux(i,j,m)=0.5*(H1(i,j)+H1(i,j-1))*yflux(i,j,m) !*dvm(i,j) !RMarsooli enddo enddo enddo do m=2,mmm-1 do j=2,jm-1 do i=2,im-1 enf(i,j,m)=enb(i,j,m) & -(xflux(i+1,j,m)-xflux(i,j,m) & +yflux(i,j+1,m)-yflux(i,j,m))*areninv(i,j)*dtw2 enf(i,j,m)=enf(i,j,m)*FSMW(i,j) enddo enddo enddo C-------------------Essential Cyclic B.C.s on m --------------- do j=1,jm do i=1,im sigth(i,j,1)=sigth(i,j,mmm-1) sigth(i,j,mmm)=sigth(i,j,2) enf(i,j,1)=enf(i,j,mmm-1) enf(i,j,mmm)=enf(i,j,2) enddo enddo C -- ------------------------------------------------------------- C theta step (MDATA) for wave energy C --------------------------------------------------------------- do m=1,mmm do j=1,jm do i=1,im uak(i,j,m)=css(m)*ud(i,j)+snn(m)*vd(i,j) enddo enddo enddo do j=2,jm-1 do i=2,im-1 if(FSMW(i,j).eq.1.) then CRM dx2inv=1./(dx(i+1,j)+dx(i-1,j)) CRM dy2inv=1./(dy(i,j+1)+dy(i,j-1)) dx2inv=1./(H1(i+1,j)+H1(i-1,j)) !RMarsooli dy2inv=1./(H2(i,j+1)+H2(i,j-1)) !RMarsooli do m=1,mmm cthth(i,j,m)= ! speed in theta coordinate & Fct(i,j,m)*grav/(2.*cth(i,j,m))*fcoshinv(kthD(i,j,m))**2 & *( snn(m)*(d(i+1,j)-d(i-1,j))*dx2inv . *FSMW(i+1,j)*FSMW(i-1,j) & -css(m)*(d(i,j+1)-d(i,j-1))*dy2inv . *FSMW(i,j+1)*FSMW(i,j-1)) & +snn(m)*(uak(i+1,j,m)-uak(i-1,j,m))*dx2inv & *FSMW(i+1,j)*FSMW(i-1,j) & -css(m)*(uak(i,j+1,m)-uak(i,j-1,m))*dy2inv & *FSMW(i,j+1)*FSMW(i,j-1) c Protect against refraction CFL violation, but, if invoked, may cause C errors in shoaling waters cthth(i,j,m)=min(cthth(i,j,m),0.99*dth/dtw2) enddo endif enddo enddo do m=1,mmm do j=1,jm cthth(im,j,m)=cthth(im-1,j,m) cthth(1,j,m)=cthth(2,j,m) enddo do i=1,im cthth(i,jm,m)=cthth(i,jm-1,m) cthth(i,1,m)=cthth(i,2,m) enddo enddo C call stepmpd(enb,en,enf,cthth,im,jm,mmm,dtw2,dthinv,nitera) C --------------------------------------------------------------- C Source step for wave energy C --------------------------------------------------------------- do i=2,im-1 do j=2,jm-1 if(FSMW(i,j).eq.1.) then C Calculate Hm and qdiss for later use in depth-induced breaking c cc RMarsooli, gama based on Ruessink et al. (2003) cc gama= 0.76*kthav(i,j)*d(i,j)+0.29 cc gama=max(0.6,min(gama,0.8)) cc Hm(i,j)=gama*d(i,j) fb=8.0*ent(i,j)/grav/Hm(i,j)**2 C fb=min(fb,1.0) ! necessary ? qdiss(i,j)=exp((qdiss(i,j)-1)/fb) C do m=2,mmm-1 C Create neg. "tail" to oppose wave travelling opposit to wind. C Const. 0.4 from Donelan 1999. tht_ops=tht_wnd(i,j)+pi tpsw(i,j,m)=fsprm(beta,m,mmm,tht_wnd(i,j),thetam(m)) & -0.4*fsprm(beta,m,mmm,tht_ops ,thetam(m)) swin(i,j,m)=Cin *exp(-swcon *ageinv(i,j))*FSMW(i,j) swin(i,j,m)=max(swin(i,j,m),10. )*ustp(i,j)**3*tpsw(i,j,m) C Surface dissipation. bdis is used for wind driven portion of en. C Reduced swfct*bdis is used for swell portion. bb=bdis if(abs(swin(i,j,m)).lt.1.e-6) bb=swfct*bdis sdiss(i,j,m)=adis*swin(i,j,m)+bb*sigth(i,j,m)*en(i,j,m) enf(i,j,m)=enf(i,j,m)+(swin(i,j,m)-sdiss(i,j,m))*dtw2 C Bottom dissipation, first due to friction, second to depth- C induced breaking a la Battjes and Janssen 1978. bdis1(i,j,m)=.003*(sigthb(i,j,m) & *sqrt(2.*enb(i,j,m)/grav)*fsinhinv(kthD(i,j,m)))**3 bdis2(i,j,m)=(en(i,j,m)/ent(i,j))*grav*sigth(i,j,m) & /(8.*pi)*qdiss(i,j)*(Hm(i,j))**2 enf(i,j,m)=enf(i,j,m)-(bdis1(i,j,m)+bdis2(i,j,m))*dtw2 enf(i,j,m)=max(.000001,enf(i,j,m))*FSMW(i,j) enddo C endif enddo enddo C RMarsooli Jan2016, energy dissipation due to vegetation drag C Based on Mendez and Losada (2004) IF(WETLAND.EQ.'INCLUDE') THEN DO II=1,NUMVEG I=IDXVEG(II) J=IDYVEG(II) if(FSMW(I,J).eq.1.and.D(I,J).GT.WETMIN.and. . WN(I,J).LT.62.8)then heightveg=AMIN1(D(I,J),HVEGA(I,J)) !----------Kobayashi Et Al. Formula------------ IF(VEGCDWAVE.EQ.'KOB_ETAL') THEN wavevel=0.000001 . +(3.1416*WVHT(I,J)/(WVPD(I,J)+SMALLVAL)) . *COSH(WN(I,J)*heightveg)/(SINH(WN(I,J)*D(I,J))+SMALLVAL) vegQ(I,J)=wavevel*BVEG(I,J)/0.000001 CDVEGW(I,J)=0.08+(2200./(vegQ(I,J)+SMALLVAL))**2.4 !----------Mendez and Losada Formula----------- ELSEIF(VEGCDWAVE.EQ.'MEND&LOS') THEN wavehrms=WVHT(I,J)/1.41 wavevel=0.000001 . +(wavehrms*sigp(I,J)/2.) . *COSH(WN(I,J)*heightveg)/(SINH(WN(I,J)*D(I,J)+SMALLVAL)) relheight=max(SMALLVAL , heightveg/D(I,J)) vegQ(I,J)=wavevel*(6.2832/(sigp(I,J)+SMALLVAL)) . /BVEG(I,J) vegQ(I,J)=vegQ(I,J)/(relheight**0.76) CDVEGW(I,J)=exp(-0.0138*vegQ(I,J))/(vegQ(I,J)**0.3+SMALLVAL) !--------Anderson and Smith (Re) Formula------- ELSEIF(VEGCDWAVE.EQ.'AN&SM_Re') THEN wavehrms=WVHT(I,J)/1.41 wavevel=0.000001 . +(wavehrms*sigp(I,J)/2.) . *COSH(WN(I,J)*heightveg)/(SINH(WN(I,J)*D(I,J)+SMALLVAL)) relheight=max(SMALLVAL , heightveg/D(I,J)) vegQ(I,J)=wavevel*BVEG(I,J)/0.000001 vegQ(I,J)=vegQ(I,J)/relheight**1.5 CDVEGW(I,J)=0.11+(2067.7/(vegQ(I,J)+SMALLVAL))**0.64 !--------Anderson and Smith (KC) Formula------- ELSEIF(VEGCDWAVE.EQ.'AN&SM_KC') THEN wavehrms=WVHT(I,J)/1.41 wavevel=0.000001 . +(wavehrms*sigp(I,J)/2.) . *COSH(WN(I,J)*heightveg)/(SINH(WN(I,J)*D(I,J)+SMALLVAL)) relheight=max(SMALLVAL , heightveg/D(I,J)) vegQ(I,J)=wavevel*(6.2832/(sigp(I,J)+SMALLVAL)) . /BVEG(I,J) vegQ(I,J)=vegQ(I,J)/relheight**1.5 CDVEGW(I,J)=0.97+(33.5/(vegQ(I,J)+SMALLVAL))**1.69 ELSEIF(VEGCDWAVE.EQ.'CONSTANT') THEN ELSE write(*,*) 'Error' stop ENDIF CDVEGW(I,J)=min(CDVEGW(I,J),10.) temp_var=kthav(I,J)*heightveg tempnum=(sinh(temp_var))**3 + 3.*(sinh(temp_var)) tempden=3.*kthav(i,j)*( cosh(kthav(I,J)*D(I,J)) )**3 bdisveg0=sqrt(2./3.1416/grav)*grav**2 . *CDVEGW(I,J)*NVEG(I,J)*BVEG(I,J) . *(kthav(I,J)/(sigthav(I,J)+SMALLVAL))**3 . *(tempnum/(tempden+SMALLVAL))*sqrt(ent(I,J)) do m=2,mmm-1 bdisveg(I,J,m)=bdisveg0*en(I,J,m) enf(I,J,m)=enf(I,J,m)-bdisveg(I,J,m)*dtw2 enf(I,J,m)=max(e_min,enf(I,J,m))*FSMW(I,J) enddo c c if(mod(INTX,4800).eq.0) c . write(136,'(f10.5,2i5,5f15.7)') THOUR,i,j,wavevel, c . WVHT(I,J),WVPD(I,J),WN(I,J),sigp(I,J) endif c ENDDO ENDIF CRM if(iint.eq.iprint) then CRM write(6,'(''enb ='',10e12.3)') (enb (12,6,m),m=1,26) CRM write(6,'(''sigthb='',10e12.3)') (sigthb(12,6,m),m=1,26) CRM write(6,'(''kthD ='',10e12.3)') (kthD (12,6,m),m=1,26) CRM write(6,'(''bdis1='',10e12.3)') (bdis1(12,6,m),m=1,26) CRM write(6,'(''bdis2='',10e12.3)') (bdis2(12,6,m),m=1,26) CRM write(6,'(''qdiss ='',10e12.3)') (qdiss (i,6),i=1,24) CRM write(6,'(''Hm ='',10e12.3)') (Hm (i,6),i=1,24) CRM endif C --------------------------------------------------------------- C Efficiency could be improved (slightly) for mode=0 or mode=2 C calculations by adding code which would bypass using u and v. CRM if(mode.eq.2) then IF (TOR.EQ.'BAROTROPIC') THEN !RMarsooli do k=1,kb do j=1,jm do i=1,im u(i,j,k)=ua(i,j) v(i,j,k)=va(i,j) enddo enddo enddo endif C Include interactive radiation stress/current terms CRM call cur2wave(im,jm,kb,dx,dy,fsm,z,dz,u,v,kthavD, CRM & Sradx,Sradxy,Srady,Srad,ud,vd,iint) call cur2wave(im,jm,kb,H1,H2,FSMW,z,dz,u,v,kthavD, !RMarsooli & Sradx,Sradxy,Srady,Srad,ud,vd,iint) do i=3,im-2 do j=3,jm-2 do m=2,mmm-1 C Eq.(14) has D instead of kth*D sorc= kth(i,j,m)*dt(i,j)*en(i,j,m) & * ( csm(m)*csm(m)*Sradx(i,j) & +csm(m)*snm(m)*Sradxy(i,j) & +snm(m)*snm(m)*Srady(i,j) & + Srad(i,j) ) sorc=sorc*FSMW(i-1,j)*FSMW(i+1,j)*FSMW(i,j-1) . *FSMW(i,j+1) enf(i,j,m)=enf(i,j,m)-sorc*dtw2 enddo enddo enddo C --------------------------------------------------------------- C x-y step for frequency C --------------------------------------------------------------- do m=1,mmm-1 do j=2,jm-1 do i=2,im-1 cc cgxm=cg(i,j,m)*csm(m)+ud(i,j) cc cgym=cg(i,j,m)*snm(m)+vd(i,j) c !RMarsooli cgx_temp=cg(i,j,m)*csm(m) !+ud(i,j) cgy_temp=cg(i,j,m)*snm(m) !+vd(i,j) cgxm=cgx_temp*COSANG(i,j)+cgy_temp*SINANG(i,j) !cgx w.r.t ZETA1 . +ud(i,j) cgym=-cgx_temp*SINANG(i,j)+cgy_temp*COSANG(i,j) !cgy w.r.t ZETA2 . +vd(i,j) adv=(0.5*cgxm*(sigthb(i+1,j,m)-sigthb(i-1,j,m)) & -0.5*abs(cgxm) & *(sigthb(i+1,j,m)-2.*sigthb(i,j,m)+sigthb(i-1,j,m)) ) CRM & /dx(i,j) & /H1(i,j) !RMarsooli & +(0.5*cgym*(sigthb(i,j+1,m)-sigthb(i,j-1,m)) & -0.5*abs(cgym) & *(sigthb(i,j+1,m)-2.*sigthb(i,j,m)+sigthb(i,j-1,m)) ) CRM & /dy(i,j) & /H2(i,j) !RMarsooli sigthf(i,j,m)=sigthb(i,j,m)-adv*FSMW(i,j)*dtw2 !---- Nudge in wind driven sigp only when fspr > 0. -------------- Cf= sqrt(fspr(i,j,m))*dtw2*sigthb(i,j,m) sigthf(i,j,m)=sigthf(i,j,m)/(1.+Cf)+sigp(i,j)*Cf/(1.+Cf) enddo enddo enddo C --------------------------------------------------------------- C Source step for frequency C --------------------------------------------------------------- C N.B. Effect of time derivative of dt neglected for now. do i=2,im-1 do j=2,jm-1 CRM dx2inv=1./(dx(i+1,j)+dx(i-1,j)) CRM dy2inv=1./(dy(i,j+1)+dy(i,j-1)) dx2inv=1./(H1(i+1,j)+H1(i-1,j)) !RMarsooli dy2inv=1./(H2(i,j+1)+H2(i,j-1)) !RMarsooli do m=2,mmm-1 sor=-cg(i,j,m)*kth(i,j,m)*0.5*( & (csm(m)*csm(m))*(ud(i+1,j)-ud(i-1,j))*dx2inv & +(csm(m)*snm(m))*(vd(i+1,j)-vd(i-1,j))*dx2inv & +(csm(m)*snm(m))*(ud(i,j+1)-ud(i,j-1))*dy2inv & +(snm(m)*snm(m))*(vd(i,j+1)-vd(i,j-1))*dy2inv ) & +(sigth(i,j,m)/d(i,j))*(Fcg(i,j,m)-0.5) & *( ud(i,j)*(d(i+1,j)-d(i-1,j))*dx2inv & +vd(i,j)*(d(i,j+1)-d(i,j-1))*dy2inv ) sor=sor*FSMW(i-1,j)*FSMW(i+1,j)*FSMW(i,j-1) . *FSMW(i,j+1) sigthf(i,j,m)=sigthf(i,j,m)+sor*dtw2 sigthf(i,j,m)=max(sigthf(i,j,m),0.1) enddo enddo enddo do j=2,jm-1 do m=2,mmm-1 sigthf(1,j,m)=sigthf(2,j,m) sigthf(im,j,m)=sigthf(im-1,j,m) enddo enddo do i=2,im-1 do m=2,mmm-1 sigthf(i,1,m)=sigthf(i,2,m) sigthf(i,jm,m)=sigthf(i,jm-1,m) enddo enddo C-------------------------------------------------------------- C Boundary Condtions C Make sure that fsm=1 on open boundaries and C fsm=0 for closed boundaries. C For solid boundaries waves should pass through the boundaries C assuming they are thus dissipated. Otherwise, one of the C following 2 B.C.s can be implemeneted (uncommented) depending C on the desired physics. C C---------------------Reflection at side walls ----------- c do i=1,im c do m=1,12 c enf(i,2,m)=enf(i,3,26-m) c enf(i,jm-1,26-m)=enf(i,jm-2,m) c enddo c enddo C--- For cyclic B.C.s on the north-south boundaries ---------- C--- to create j-independent problem ------------------------- c do m=2,mm c do i=1,im c enf(i,2,m)=enf(i,jm-2,m) !cyclic B.C.s c enf(i,jm-1,m)=enf(i,3,m) c enddo c enddo C--- For open B.C.s on the eastern boundary ------------------ CRM do j=1,jm CRM do m=1,mmm CRM cgnd=cgx(im,j,m)/(dx(im,j)+dx(im-1,j))*dtw CRM cgnd=cgx(im-1,j,m)/(H1(im,j)+H1(im-1,j))*dtw !RMarsooli CRM enf(im,j,m)=en(im,j,m) CRM & -(cgnd+abs(cgnd))*(en(im,j,m)-en(im-1,j,m)) CRM & -(cgnd-abs(cgnd))*(ene(j,m)-en(im,j,m)) CRM enddo CRM enddo C--- For open B.C.s on the western boundary ------------------- CRM do j=1,jm CRM do m=1,mmm CRM cgnd=cgx(2,j,m)/(dx(2,j)+dx(1,j))*dtw CRM cgnd=cgx(2,j,m)/(H1(2,j)+H1(1,j))*dtw !RMarsooli CRM enf(1,j,m)=en(1,j,m) CRM & -(cgnd+abs(cgnd))*(en(1,j,m)-enw(j,m)) CRM & -(cgnd-abs(cgnd))*(en(2,j,m)-en(1,j,m)) CRM enddo CRM enddo C--- For open B.C.s on the southern boundary ------------------- CRM do i=1,im CRM do m=1,mmm CCRM cgnd=cgy(i,2,m)/(dy(i,2)+dy(i,1))*dtw CRM cgnd=cgy(i,2,m)/(H2(i,2)+H2(i,1))*dtw !RMarsooli CRM enf(i,1,m)=en(i,1,m) CRM & -(cgnd+abs(cgnd))*(en(i,1,m)-ens(i,m)) CRM & -(cgnd-abs(cgnd))*(en(i,2,m)-en(i,1,m)) CRM enddo CRM enddo C BC for SWELLS, RMarsooli_2015 IF(NUMWMBC.EQ.0) THEN !NO WAVE DATA IS AVAILABLE NWAVEBC=NUMEBC ELSE NWAVEBC=NUMWMBC ENDIF DO N=1,NWAVEBC IF(NUMWMBC.EQ.0) THEN I=IETA(N) J=JETA(N) IC=ICON(N) JC=JCON(N) theta_swell=thtav(IC,JC) WHBFLX(N)=Hs(IC,JC) ELSE I=IWM(N) J=JWM(N) IC=IWMC(N) JC=JWMC(N) cc theta_temp=WDBFLX(N)*3.141593/180.+ANG(I,J) cc if(theta_temp.gt.6.2832) theta_temp=theta_temp-6.2832 cc theta_swell=1.5708-theta_temp cc if(theta_temp.gt.4.7124) theta_swell=7.854-theta_temp c WDBFLX is w.r.t. true North, toward, clockwise + c Convert it to Mellor convention: w.r.t true east, toward, clockwise:-, counterclockwise:+ theta_temp=WDBFLX(N)*3.141593/180.-1.5708 if(theta_temp.lt.0.) . theta_temp=WDBFLX(N)*3.141593/180.+4.7121 theta_swell=-theta_temp if(theta_swell.lt.-3.141593) . theta_swell=6.2832-theta_temp cc write(*,*) theta_swell c WDBFLX is w.r.t. true North, toward, clockwise + c Determine swell direction w.r.t. zeta1, (clockwise:-, counterclockwise:+) cc ZETAN=6.2832-ANG(I,J)+1.5708 cc IF(ANG(I,J).LT.1.5708) ZETAN=1.5708-ANG(I,J) cc theta_swell=-(WDBFLX(N)*3.141593/180.-ZETAN) cc if(theta_swell.LT.-3.141593) theta_swell=6.2832+theta_swell c WDBFLX in wave boundary condition file is w.r.t. true North, toward, clockwise + cccccc theta_swell=WDBFLX(N)*3.141593/180.+ANG(I,J) cccccc if(theta_swell.gt.6.283185) theta_swell=theta_swell-6.283185 cc diffmin=9999999. cc if(theta_swell.ge.0.) then cc do ii=1,mmm cc if(thetam(ii).ge.0.) then cc diffmin0=abs(theta_swell-thetam(ii)) cc if(diffmin0.lt.diffmin) then cc m=i cc diffmin=diffmin0 cc endif cc endif cc enddo cc else cc do ii=1,mmm cc if(thetam(ii).lt.0.) then cc diffmin0=abs(theta_swell-thetam(ii)) cc if(diffmin0.lt.diffmin) then cc m=i cc diffmin=diffmin0 cc endif cc endif cc enddo cc endif ENDIF do ii=1,mmm enf(I,J,ii)=0. enddo DO m=1,mmm enboundary=9.81*(WHBFLX(N)/4.)**2 & *fsprm(beta,m,mmm,theta_swell,thetam(m)) enboundary=max(.000001,enboundary*FSMW(I,J)) IF(J.EQ.JC) THEN !WEST AND EAST BOUNDARIES cgnd=cgx(IC,JC,m)/(H1(I,J)+H1(IC,JC))*dtw IF(IC.GT.I) THEN !WEST BOUNDARY enf(I,J,m)=en(I,J,m) & -(cgnd+abs(cgnd))*(en(I,J,m)-enboundary) & -(cgnd-abs(cgnd))*(en(IC,JC,m)-en(I,J,m)) ELSE !EAST BOUNDARY enf(I,J,m)=en(I,J,m) & -(cgnd+abs(cgnd))*(en(I,J,m)-en(IC,JC,m)) & -(cgnd-abs(cgnd))*(enboundary-en(I,J,m)) ENDIF ELSE !SOUTH AND NORTH BOUNDARIES cgnd=cgx(IC,JC,m)/(H2(I,J)+H2(IC,JC))*dtw c cgnd=cgx(IC,JC,m)/(DELY(I,J)+DELY(IC,JC))*dtw IF(JC.GT.J) THEN !SOUTH BOUNDARY enf(I,J,m)=en(I,J,m) & -(cgnd+abs(cgnd))*(en(I,J,m)-enboundary) & -(cgnd-abs(cgnd))*(en(IC,JC,m)-en(I,J,m)) ELSE !NORTH BOUNDARY enf(I,J,m)=en(I,J,m) & -(cgnd+abs(cgnd))*(en(I,J,m)-en(IC,JC,m)) & -(cgnd-abs(cgnd))*(enboundary-en(I,J,m)) ENDIF ENDIF ENDDO ENDDO C RMarsooli !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c do J=2,JM-1 c DO m=1,mmm c enf(2,J,m)=enf(3,J,m)*FSMW(2,J) c enddo c enddo cc DO i=2,IM-1 cc DO m=1,mmm cc enf(i,2,m)=enf(i,jm-2,m)*FSMW(i,2) !cyclic B.C.s cc enf(i,jm-1,m)=enf(i,3,m)*FSMW(i,jm-1) cc enddo cc enddo c do i=1,im-2 !im !reflection side walls c do m=1,6 !12 c enf(i,2,m)=enf(i,3,14-m) !26-m) c! enf(i,jm-1,26-m)=enf(i,jm-2,m) c enf(i,jm-1,14-m)=enf(i,jm-2,m) c enddo c enddo c C------------------------------------------------------------ C Mask out land values and apply Asselin filter C------------------------------------------------------------ do m=1,mmm do j=1,jm do i=1,im enf(i,j,m)=enf(i,j,m)*FSMW(i,j) en(i,j,m)=en(i,j,m) & +0.5*smothh*(enb(i,j,m)-2.*en(i,j,m)+enf(i,j,m)) sigth(i,j,m)=sigth(i,j,m) & +0.5*smothh*(sigthb(i,j,m)-2.*sigth(i,j,m)+sigthf(i,j,m)) en(i,j,m)=max(.000001,en(i,j,m))*FSMW(i,j) enf(i,j,m)=max(.000001,enf(i,j,m))*FSMW(i,j) enddo enddo enddo C------ Reset time sequence ----------------------------------- do m=1,mmm do j=1,jm do i=1,im enb(i,j,m)=en(i,j,m) en(i,j,m)=enf(i,j,m) sigthb(i,j,m)=sigth(i,j,m) sigth(i,j,m)=sigthf(i,j,m) enddo enddo enddo c ----Calculate total dissipation for use in the tke equation --- do j=1,jm do i=1,im sdistot(i,j)=summ(sdiss,dth,i,j,im,jm,mmm)*FSMW(i,j) do m=1,mmm tpsw(i,j,m)=bdis1(i,j,m)+bdis2(i,j,m) enddo bdistot(i,j)=summ(tpsw,dth,i,j,im,jm,mmm)*FSMW(i,j) enddo enddo C------------------------------------------------------------ C Determine some diagnostic wave properties C------------------------------------------------------------ do j=1,jm do i=1,im if(FSMW(i,j).eq.1.) then C Determine average wave propagation angle do m=1,mmm tpsw(i,j,m)=en(i,j,m)*csm(m) enddo csmav=summ(tpsw,dth,i,j,im,jm,mmm) do m=1,mmm tpsw(i,j,m)=en(i,j,m)*snm(m) enddo snmav=summ(tpsw,dth,i,j,im,jm,mmm) if(csmav.eq.0.) csmav=csmav+1.e-7 thtav(i,j)=atan(snmav/csmav) if (snmav.gt.0..and.csmav.lt.0.) thtav(i,j)=thtav(i,j)+pi if (snmav.lt.0..and.csmav.lt.0.) thtav(i,j)=thtav(i,j)-pi swintot(i,j)=summ(swin,dth,i,j,im,jm,mmm) C Determine average frequency and wave number do m=1,mmm tpsw(i,j,m)=en(i,j,m)*sigth(i,j,m) enddo ent(i,j)=summ(en,dth,i,j,im,jm,mmm) sigthav(i,j)=summ(tpsw ,dth,i,j,im,jm,mmm)/ent(i,j) do m=1,mmm tpsw(i,j,m)=en(i,j,m)*kth(i,j,m) enddo kthav(i,j)=summ(tpsw ,dth,i,j,im,jm,mmm)/ent(i,j) kthavD(i,j)=kthav(i,j)*dt(i,j) endif enddo enddo return end C-------------------------------------------------------------------- subroutine cur2wave(im,jm,kb,dx,dy,FSMW,z,dz,u,v,kthavD, & Sradx,Sradxy,Srady,Srad,ud,vd,iint) C-------------------------------------------------------------------- C Current to wave interaction terms C Calculates vertical averages using spectrally averaged functions C provided by subroutine specavs. C Terms described by Eq. (14) in Mellor et al (2008) C-------------------------------------------------------------------- implicit none C Input integer i,j,k,im,jm,kb,iint real kthavD(im,jm),FSMW(im,jm) real dx(im,jm),dy(im,jm),z(kb),dz(kb),u(im,jm,kb),v(im,jm,kb) C Output real Sradx(im,jm),Sradxy(im,jm),Srady(im,jm),Srad(im,jm), & Srad0(im,jm),ud(im,jm),vd(im,jm) C Internal real sumk real F1(kb),F2(kb),F3(kb),F4(kb), & tp1(kb),tp2(kb),tp3(kb), & gradxu(kb),gradxv(kb),gradyu(kb),gradyv(kb), & dxiv,dyiv do j=3,jm-2 do i=3,im-2 call specavs (kthavD(i,j),F1,F2,F3,F4,z,kb,i,j,iint) dxiv=1./dx(i,j) dyiv=1./dy(i,j) do k=1,kb-1 tp1(k)=0.25*(u(i,j,k)+u(i+1,j,k))*(F3(k)+F3(k+1))*dz(k) tp2(k)=0.25*(v(i,j,k)+v(i,j+1,k))*(F3(k)+F3(k+1))*dz(k) gradxu(k)=(u(i+1,j,k)-u(i,j,k))*dxiv*FSMW(i+1,j) . *FSMW(i,j) gradxv(k)=0.25*(v(i+1,j+1,k)-v(i-1,j,k) & +v(i+1,j,k)-v(i-1,j+1,k))*dxiv ! corrected 07/07/09 & *FSMW(i+1,j+1)*FSMW(i-1,j)*FSMW(i+1,j) . *FSMW(i-1,j+1) gradyu(k)=0.25*(u(i+1,j+1,k)-u(i,j-1,k) & +u(i,j+1,k)-u(i+1,j-1,k))*dyiv ! corrected 07/07/09 & *FSMW(i+1,j+1)*FSMW(i,j-1)*FSMW(i,j+1) . *FSMW(i+1,j-1) gradyv(k)=(v(i,j+1,k)-v(i,j,k))*dyiv*FSMW(i,j+1) . *FSMW(i,j) enddo ud(i,j)=sumk(tp1,1.,kb) ! Dopler velocity vd(i,j)=sumk(tp2,1.,kb) ! '' '' C do k=1,kb-1 tp1(k)=0.5*(F1(k)+F1(k+1))*gradxu(k)*dz(k) tp2(k)=0.5*(F1(k)+F1(k+1))*(gradxv(k)+gradyu(k))*dz(k) tp3(k)=0.5*(F1(k)+F1(k+1))*gradyv(k)*dz(k) enddo Sradx(i,j)=sumk(tp1,1.,kb) Sradxy(i,j)=sumk(tp2,1.,kb) Srady(i,j)=sumk(tp3,1.,kb) do k=1,kb-1 tp1(k)=0.5*(F2(k)+F2(k+1))*gradxu(k)*dz(k) tp2(k)=0.5*(F2(k)+F2(k+1))*gradyv(k)*dz(k) enddo Srad(i,j)=-(sumk(tp1,1.,kb)+sumk(tp2,1.,kb)) enddo enddo do j=1,jm ud(im-1,j)=ud(im-2,j) ud(im,j)=ud(im-1,j) ud(2,j)=ud(3,j) ud(1,j)=ud(2,j) vd(im-1,j)=vd(im-2,j) vd(im,j)=vd(im-1,j) vd(2,j)=vd(3,j) vd(1,j)=vd(2,j) enddo do i=1,im ud(i,jm-1)=ud(i,jm-2) ud(i,jm)=ud(i,jm-1) ud(i,2)=ud(i,3) ud(i,1)=ud(i,2) vd(i,jm-1)=vd(i,jm-2) vd(i,jm)=vd(i,jm-1) vd(i,2)=vd(i,3) vd(i,1)=vd(i,2) enddo return end C---------------------------------------------------------------- subroutine specavs (kpD,F1,F2,F3,F4,z,kb,i,j,iint) C This is a subroutine to deliver spectrally averaged F1 - F4. C averaged. C F1-F4 are spectrally averaged components of wave radiation C stresses. C N.B. This subroutine uses the file, specavs. implicit none integer me,ne parameter (me=17,ne=21) integer i,j,k,m,n,kb,init,iint real kpD,kpDd(me) real F1d(me,ne),F2d(me,ne),F3d(me,ne),F4d(me,ne) real F1i(ne),F2i(ne),F3i(ne),F4i(ne) real F1(kb),F2(kb),F3(kb),F4(kb) real zeta(ne),zetf(ne),z(kb) save init,kpDd,F1d,F2d,F3d,F4d,zeta,zetf data init/0/ C---------------------------------------------------------------- if(init.eq.0) then C---------------------------------------------------------------- do m=1,16 kpDd(m)=0.2*float(m-1) enddo crm write(6,'(''Read spectrally averaged functions from specavs'')') read(11,'(2x)') crm write(6,'(/,'' F1'')') do n=1,ne read(11,'(f6.2,15(f6.3),2f8.3)') & zeta(n),(F1d(m,n),m=2,me),zetf(n) F1d(1,n)=1.0/0.1 ! avoid infinity at kpD=0 crm write(6,'(f6.2,f6.2,15(f6.3),2f8.3)') crm & zeta(n),(F1d(m,n),m=1,me),zetf(n) enddo read(11,'(2x)') crm write(6,'(/,'' F2'')') do n=1,ne read(11,'(f6.2,15(f6.3),2f8.3)') & zeta(n),(F2d(m,n),m=2,me),zetf(n) F2d(1,n)=0. ! asymptote for kpD=0 c write(6,'(f6.2,16(f6.3),2f8.3)') c & zeta(n),(F2d(m,n),m=1,me),zetf(n) enddo read(11,'(2x)') crm write(6,'(/,'' F3'')') do n=1,ne read(11,'(f6.2,15(f6.3),2f8.3)') & zeta(n),(F3d(m,n),m=2,me),zetf(n) F3d(1,n)=1.5+zeta(n) c write(6,'(f6.2,16(f6.3),2f8.3)') c & zeta(n),(F3d(m,n),m=1,me),zetf(n) enddo read(11,'(2x)') crm write(6,'(/,'' F4'')') do n=1,ne read(11,'(f6.2,15(f6.3),2f8.3)') & zeta(n),(F4d(m,n),m=2,me),zetf(n) F4d(1,n)=1.+zeta(n) c write(6,'(f6.2,16(f6.3),2f8.3)') c & zeta(n),(F4d(m,n),m=1,me),zetf(n) enddo C---------------------------------------------------------------- init=1 endif C---------------------------------------------------------------- C Solve for F1, F2, F3, F4 C---------------------------------------------------------------- m=kpD/0.2+1 if (m.lt.16) then C Interpolate w.r.t. kpD =< 3. (0 < m< 17) do n=1,ne F1i(n)=F1d(m,n)+(F1d(m+1,n)-F1d(m,n))*(kpD-kpDd(m))/0.2 F2i(n)=F2d(m,n)+(F2d(m+1,n)-F2d(m,n))*(kpD-kpDd(m))/0.2 F3i(n)=F3d(m,n)+(F3d(m+1,n)-F3d(m,n))*(kpD-kpDd(m))/0.2 F4i(n)=F4d(m,n)+(F4d(m+1,n)-F4d(m,n))*(kpD-kpDd(m))/0.2 enddo C Interpolate w.r.t. z do k=1,kb-1 n=-z(k)/.05+1 F1(k)=F1i(n)+(F1i(n+1)-F1i(n))*(-z(k)+zeta(n))/.05 F2(k)=F2i(n)+(F2i(n+1)-F2i(n))*(-z(k)+zeta(n))/.05 F3(k)=F3i(n)+(F3i(n+1)-F3i(n))*(-z(k)+zeta(n))/.05 F4(k)=F4i(n)+(F4i(n+1)-F4i(n))*(-z(k)+zeta(n))/.05 enddo F1(kb)=F1i(ne) F2(kb)=F2i(ne) F3(kb)=F3i(ne) F4(kb)=F4i(ne) else C For m=17, kpD=100 and F1 is scaled different from m<17 C Use similarity for large kpD > 3. (m = 17) C e.g., F1(kpD, z) = (kpD/100)*F1(100, kpD*z/100) do k=1,kb F1(k)=0. F2(k)=0. F3(k)=0. F4(k)=0. n=1-z(k)*kpD/.1 if(n.lt.ne) then !asumptotes for kpD > 3 F1(k)=F1d(17,n)+(F1d(17,n+1)-F1d(17,n))/.001 & *(-z(k)*kpD/100.+zetf(n)) F2(k)=F1(k) F3(k)=F3d(17,n)+(F3d(17,n+1)-F3d(17,n))/.001 & *(-z(k)*kpD/100.+zetf(n)) F3(k)=F3(k)*kpD*0.01 F4(k)=F4d(17,n)+(F4d(17,n+1)-F4d(17,n))/.001 & *(-z(k)*kpD/100.+zetf(n)) endif enddo endif return end C---------------------------------------------------------------- subroutine speeds_mel(kpD,Fcg,Fct) C This is a subroutine to deliver spectrally averaged Fcg, Fct C Fcg = ratio of cg to cp where cp is spectral peak phase speed. C Fct = is factor in cthth (theta speed) to make it a spectrally C averaged. Fcg and Fct are functions of kpD. implicit none integer k real kpD,kpDd(17),Fcg,Fct real Fcgd(17),Fctd(17),c_cp(17),dum(17) data kpDd/ & 0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, & 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 100.0/ data c_cp/ & 1.000, 1.005, 1.001, 0.968, 0.951, 0.936, 0.924, 0.915, 0.910, & 0.906, 0.905, 0.904, 0.904, 0.904, 0.905, 0.905, 0.910/ data Fcgd/ & 1.000, 0.984, 0.928, 0.834, 0.756, 0.689, 0.635, 0.592, 0.560, & 0.535, 0.516, 0.502, 0.491, 0.483, 0.477, 0.472, 0.472/ data dum/ & 1.000, 0.987, 0.950, 0.897, 0.837, 0.776, 0.720, 0.671, 0.631, & 0.598, 0.573, 0.554, 0.539, 0.529, 0.521, 0.515, 0.500/ data Fctd/ & 1.000, 0.984, 0.954, 0.931, 0.892, 0.858, 0.830, 0.809, 0.795, & 0.787, 0.784, 0.786, 0.791, 0.801, 0.814, 0.829, 0.829/ k=kpD/0.2+1 CRM if(kpD.lt.1) call printall if (k.lt.16) then Fcg=Fcgd(k)+(Fcgd(k+1)-Fcgd(k))*(kpD-kpDd(k))/0.2 Fct=Fctd(k)+(Fctd(k+1)-Fctd(k))*(kpD-kpDd(k))/0.2 else k=16 Fcg=Fcgd(k)+(Fcgd(k+1)-Fcgd(k))*(kpD-kpDd(k))/97. Fct=Fctd(k)+(Fctd(k+1)-Fctd(k))*(kpD-kpDd(k))/97. endif return end C------------------------------------------------------------------ subroutine stress_mel (iint,nwavemel,Hs,ageinv,usdif,vsdif,cdd, & tpx0,tpy0,cdmax,ttx0,tty0,ustw,ustp,im,jm, . FSMW,SIG,WN) C----------------------------------------------------------------- C Input: usdif,vsdif = difference between 10m wind velocity and C surface velocity C Hs = significant wave height C ageinv = inverse wave age C Output: tpx0,tpy0 = wind pressure wave slope correlation = C wind induced surface stress C ttx0,tty0 = turbulence surface stress (non-zero only C for low winds (< about 4 m s-1) C Surface roughness due to Donelan (1990) C----------------------------------------------------------------- C implicit none integer i,j,im,jm real u1,u2,z0w,z0t,cdt,cdw,grav,pi,r,rhalf real Hs(im,jm),ageinv(im,jm),usdif(im,jm),vsdif(im,jm) real cdd(im,jm),tpx0(im,jm),tpy0(im,jm),ustw(im,jm),ustp(im,jm) real ttx0(im,jm),tty0(im,jm),FSMW(im,jm) real kappa,nu real SIG(im,jm),WN(im,jm) !RMarsooli data nu/1.8e-6/,r/0.0011630/,rhalf/.0034100/,kappa/0.41/ data z0w/1.e-6/,z0t/1.e-6/,grav/9.807/,pi/3.1415926/ do j=2,jm-1 do i=2,im-1 cdd(i,j)=0.0 !RMarsooli if(FSMW(i,j).gt.0.) then u2=usdif(i,j)**2+vsdif(i,j)**2 u1=sqrt(u2) z0w=1.38e-4*Hs(i,j)*(ageinv(i,j))**2.66+1.e-5 !Donelan z0t=0.18*nu/(ustw(i,j) +.000001) cdw=(kappa/log(10./z0w))**2 cdw=min(cdw,cdmax) cdt=(kappa/log(10./z0t))**2 C The transition from smooth surface turbulent flow and friction C drag to a wave surface and form drag is abrupt; this probably C needs improvement but nevertheless gives a continuos Cd vs U10 C result that agrees with data. if(cdw.gt.cdt) then tpx0(i,j)=r*cdw*u1*usdif(i,j) tpy0(i,j)=r*cdw*u1*vsdif(i,j) ustp(i,j)=sqrt(sqrt(tpx0(i,j)**2+tpy0(i,j)**2)) ttx0(i,j)=0. tty0(i,j)=0. cdt=0. else ttx0(i,j)=r*cdt*u1*usdif(i,j) tty0(i,j)=r*cdt*u1*vsdif(i,j) tpx0(i,j)=0. tpy0(i,j)=0. cdw=0. endif cdd(i,j)=cdw+cdt ustw(i,j)=sqrt(r*(cdw+cdt)*u2) endif enddo enddo return end C---------------------------------------------------------------- subroutine stepmpd(fb,f,ff,c,im,jm,mmm,dtw2,dthinv,nitera) C********************************************************************** C * * C * FUNCTION : Integrates conservative scalar equations. * C * * C * This is a first-order upstream scheme, which * C * reduces implicit diffusion using the Smolarkiewicz * C * iterative upstream scheme with an antidiffusive * C * velocity. * C * * C * It is based on the subroutines of Gianmaria Sannino * C * (Inter-university Computing Consortium, Rome, Italy)* C * and Vincenzo Artale (Italian National Agency for * C * New Technology and Environment, Rome, Italy), * C * downloaded from the POM FTP site on 1 Nov. 2001. * C * * C * * C * Reference: * C * * C * Smolarkiewicz, P.K.; A fully multidimensional * C * positive definite advection transport algorithm * C * with small implicit diffusion, Journal of * C * Computational Physics, 54, 325-362, 1984. * C * * C ********************************************************************** C implicit none integer itera,nitera,i,j,m,im,jm,mmm real sw,dtw2,dthinv,value_min,epsilon real fbmem(im,jm,mmm),ffmem(im,jm,mmm),c(im,jm,mmm) real flux(im,jm,mmm),speed(im,jm,mmm) real fb(im,jm,mmm),f(im,jm,mmm),ff(im,jm,mmm) real mol,abs_1,abs_2 real udx,u2dt,vdy,v2dt,wdz,w2dt value_min=1.e-9 epsilon=1.e-14 sw=1.0 C do m=1,mmm do j=1,jm do i=1,im speed(i,j,m)=c(i,j,m) fbmem(i,j,m)=fb(i,j,m) ffmem(i,j,m)=ff(i,j,m) end do end do end do C do itera=1,nitera C do m=2,mmm-1 do j=1,jm do i=1,im flux(i,j,m)=0.5e0 $ *((speed(i,j,m)+abs(speed(i,j,m))) $ *fbmem(i,j,m-1)+ $ (speed(i,j,m)-abs(speed(i,j,m))) $ *fbmem(i,j,m)) end do end do end do do j=1,jm do i=1,im flux(i,j,1)=flux(i,j,mmm-1) flux(i,j,mmm)=flux(i,j,2) end do end do C C Add net advective fluxes and step forward in time: C do m=1,mmm-1 do j=1,jm do i=1,im ff(i,j,m)=ffmem(i,j,m) $ +(flux(i,j,m)-flux(i,j,m+1))*dthinv*dtw2 end do end do end do C * Calculate the antidiffusive speed used to * C * reduce the numerical diffusion associated with the * C * upstream differencing scheme. * C * * do m=2,mmm-1 do j=1,jm do i=1,im if(ff(i,j,m).lt.value_min.or. $ ff(i,j,m-1).lt.value_min) then speed(i,j,m)=0.e0 else wdz=abs(speed(i,j,m)) w2dt=speed(i,j,m)*speed(i,j,m)*dtw2*dthinv mol=(ff(i,j,m)-ff(i,j,m-1)) $ /(ff(i,j,m)+ff(i,j,m-1)+epsilon) speed(i,j,m)=(wdz-w2dt)*mol*sw abs_1=abs(wdz) abs_2=abs(w2dt) if(abs_1.lt.abs_2)speed(i,j,m)=0.e0 end if end do end do end do C do m=1,mmm do j=1,jm do i=1,im ffmem(i,j,m)=ff(i,j,m) fbmem(i,j,m)=ff(i,j,m) end do end do end do C C End of Smolarkiewicz scheme end do return end C----------------------------------------------------------------- subroutine wave2cur (dth) C This routine supplies Sxx,Sxy.Syy for insertion into the C into the momentum equation when waves are coupled with C the circulation model. c implicit none CRM include 'pom08.c' INCLUDE 'comdeck' !RMarsooli SAVE !RMarsooli real tpxx(mmm),tpxy(mmm),tpyy(mmm),tp(mmm),summ1 real stpxx,stpxy,stpyy,stp,dxinv,dyinv real F1m(kb),F2m(kb),us(im,jm,kb) C if(iint.eq.3.and.nwavemel.eq.2) stop do j=2,jm-1 do i=2,im-1 if(FSMW(i,j).gt.0.) then do m=1,mmm tpxx(m)=kth(i,j,m)*en(i,j,m)*csm(m)*csm(m) tpxy(m)=kth(i,j,m)*en(i,j,m)*csm(m)*snm(m) tpyy(m)=kth(i,j,m)*en(i,j,m)*snm(m)*snm(m) tp(m)=kth(i,j,m)*en(i,j,m) enddo stpxx=summ1(tpxx,dth,mmm) stpxy=summ1(tpxy,dth,mmm) stpyy=summ1(tpyy,dth,mmm) stp =summ1(tp ,dth,mmm) call specavs (kthavD(i,j),F1,F2,F3,F4,z,kb,i,j,iint) do k=1,kb-1 F1m(k)=0.5*(F1(k)+F1(k+1)) F2m(k)=0.5*(F2(k)+F2(k+1)) c Sxx(i,j,k)=stpxx*F1m(k)+stp*F2m(k) c Sxx(i,j,k)=max(0.0,Sxx(i,j,k)) !RM tempsxx=stpxx*F1m(k)+stp*F2m(k) Sxx(i,j,k)=radlimit(tempsxx,Sxx(i,j,k)) c Sxy(i,j,k)=stpxy*F1m(k) c Sxy(i,j,k)=max(0.0,Sxy(i,j,k)) !RM tempsxy=stpxy*F1m(k) Sxy(i,j,k)=radlimit(tempsxy,Sxy(i,j,k)) c Syy(i,j,k)=stpyy*F1m(k)+stp*F2m(k) c Syy(i,j,k)=max(0.0,Syy(i,j,k)) !RM tempsyy=stpyy*F1m(k)+stp*F2m(k) Syy(i,j,k)=radlimit(tempsyy,Syy(i,j,k)) C tpzdist(i,j,k)=F4(k) ! subsurface momentum ! tranfer after multiplication ! by (tpx0,tpy0) after do 9000 C Calculate Stokes drift us(i,j,k)=ent(i,j)*kthav(i,j)/(sigthav(i,j)*dt(i,j)) & *(F4(k)-F4(k+1))/dz(k) ust(i,j,k)=us(i,j,k)*cos(thtav(i,j)) vst(i,j,k)=us(i,j,k)*sin(thtav(i,j)) enddo endif do k=1,kb-1 !!!!!!!!!!!!!!!!!!!RM Sxx(i,j,k)=Sxx(i,j,k)*FSMW(i,j) Sxy(i,j,k)=Sxy(i,j,k)*FSMW(i,j) Syy(i,j,k)=Syy(i,j,k)*FSMW(i,j) enddo enddo enddo write(111,'(4f15.7)') THOUR,Sxx(61,177,1),Sxy(61,177,1), . Syy(61,177,1) C RMarsooli C BOUNDARY CONDITIONS DO N=1,NUMEBC IE=IETA(N) JE=JETA(N) IC=ICON(N) JC=JCON(N) IF(FSM(IE+1,JE).EQ.0.0.AND.JE.EQ.JC) THEN DO K=1,KB Sxx(IE,JE,K)=Sxx(IE-1,JE,K) Sxy(IE,JE,K)=Sxy(IE-1,JE,K) Syy(IE,JE,K)=Syy(IE-1,JE,K) ENDDO ELSE IF(FSM(IE-1,JE).EQ.0.0.AND.JE.EQ.JC) THEN DO K=1,KB Sxx(IE,JE,K)=Sxx(IE+1,JE,K) Sxy(IE,JE,K)=Sxy(IE+1,JE,K) Syy(IE,JE,K)=Syy(IE+1,JE,K) ENDDO ELSE IF(FSM(IE,JE+1).EQ.0.0.AND.IE.EQ.IC) THEN DO K=1,KB Sxx(IE,JE,K)=Sxx(IE,JE-1,K) Sxy(IE,JE,K)=Sxy(IE,JE-1,K) Syy(IE,JE,K)=Syy(IE,JE-1,K) ENDDO ELSE IF(FSM(IE,JE-1).EQ.0.0.AND.IE.EQ.IC) THEN DO K=1,KB Sxx(IE,JE,K)=Sxx(IE,JE+1,K) Sxy(IE,JE,K)=Sxy(IE,JE+1,K) Syy(IE,JE,K)=Syy(IE,JE+1,K) ENDDO END IF ENDDO DO J=2,JM-1 DO K=1,KB Sxx(1,J,K)=Sxx(2,J,K) Sxy(1,J,K)=Sxy(2,J,K) Syy(1,J,K)=Syy(2,J,K) ust(1,J,K)=ust(2,J,K) vst(1,J,K)=vst(2,J,K) Sxx(IM,J,K)=Sxx(IM-1,J,K) Sxy(IM,J,K)=Sxy(IM-1,J,K) Syy(IM,J,K)=Syy(IM-1,J,K) ust(IM,J,K)=ust(IM-1,J,K) vst(IM,J,K)=vst(IM-1,J,K) ENDDO ENDDO DO I=2,IM-1 DO K=1,KB Sxx(I,1,K)=Sxx(I,2,K) Sxy(I,1,K)=Sxy(I,2,K) Syy(I,1,K)=Syy(I,2,K) ust(I,1,K)=ust(I,2,K) vst(I,1,K)=vst(I,2,K) Sxx(I,JM,K)=Sxx(I,JM-1,K) Sxy(I,JM,K)=Sxy(I,JM-1,K) Syy(I,JM,K)=Syy(I,JM-1,K) ust(I,JM,K)=ust(I,JM-1,K) vst(I,JM,K)=vst(I,JM-1,K) ENDDO ENDDO return end C -------------------------------------------------------------- subroutine wvfcns(sig,d,wvnu,c) implicit none real geeinv parameter(geeinv=1./9.807) C C This lookup table inputs sig and d and outputs wvnu, c and cg C using the dispersion relation C C sig = frequency (s-1) C d = water column depth (m) C wvnu = wave number (m-1) C c = phase speed (m s-1) C cg = group speed (m s-1) C f1 = sig**2*d/gee C f2 = wvnu*d C f3 = cg/c C C ------------------------------------------------------------------ integer k real sig,d,wvnu,c,cg,dknd,cgnd,fqnd real f1(81),f2(81),f3(81) data f1/ & 0.0000, 0.0500, 0.1000, 0.1500, 0.2000, 0.2500, 0.3000, & 0.3500, 0.4000, 0.4500, 0.5000, 0.5500, 0.6000, 0.6500, & 0.7000, 0.7500, 0.8000, 0.8500, 0.9000, 0.9500, 1.0000, & 1.0500, 1.1000, 1.1500, 1.2000, 1.2500, 1.3000, 1.3500, & 1.4000, 1.4500, 1.5000, 1.5500, 1.6000, 1.6500, 1.7000, & 1.7500, 1.8000, 1.8500, 1.9000, 1.9500, 2.0000, 2.0500, & 2.1000, 2.1500, 2.2000, 2.2500, 2.3000, 2.3500, 2.4000, & 2.4500, 2.5000, 2.5500, 2.6000, 2.6500, 2.7000, 2.7500, & 2.8000, 2.8500, 2.9000, 2.9500, 3.0000, 3.0500, 3.1000, & 3.1500, 3.2000, 3.2500, 3.3000, 3.3500, 3.4000, 3.4500, & 3.5000, 3.5500, 3.6000, 3.6500, 3.7000, 3.7500, 3.8000, & 3.8500, 3.9000, 3.9500, 4.0000/ data f2/ & 0.0000, 0.2255, 0.3216, 0.3973, 0.4627, 0.5218, 0.5767, & 0.6284, 0.6778, 0.7255, 0.7717, 0.8168, 0.8611, 0.9046, & 0.9476, 0.9902, 1.0324, 1.0744, 1.1163, 1.1580, 1.1997, & 1.2414, 1.2831, 1.3249, 1.3668, 1.4088, 1.4511, 1.4934, & 1.5360, 1.5788, 1.6218, 1.6651, 1.7085, 1.7523, 1.7962, & 1.8405, 1.8849, 1.9297, 1.9746, 2.0199, 2.0653, 2.1110, & 2.1569, 2.2031, 2.2495, 2.2960, 2.3428, 2.3898, 2.4369, & 2.4843, 2.5318, 2.5795, 2.6273, 2.6752, 2.7233, 2.7716, & 2.8199, 2.8684, 2.9170, 2.9657, 3.0144, 3.0633, 3.1123, & 3.1613, 3.2104, 3.2596, 3.3088, 3.3581, 3.4074, 3.4568, & 3.5063, 3.5558, 3.6053, 3.6548, 3.7044, 3.7541, 3.8037, & 3.8500, 3.9000, 3.9500, 4.0000/ data f3/ & 1.0000, 0.9834, 0.9671, 0.9510, 0.9352, 0.9196, 0.9042, & 0.8891, 0.8743, 0.8598, 0.8455, 0.8315, 0.8179, 0.8045, & 0.7914, 0.7786, 0.7662, 0.7541, 0.7422, 0.7308, 0.7196, & 0.7088, 0.6983, 0.6882, 0.6784, 0.6689, 0.6598, 0.6511, & 0.6426, 0.6345, 0.6268, 0.6193, 0.6122, 0.6054, 0.5990, & 0.5928, 0.5870, 0.5814, 0.5761, 0.5711, 0.5664, 0.5619, & 0.5577, 0.5538, 0.5500, 0.5465, 0.5432, 0.5401, 0.5373, & 0.5345, 0.5320, 0.5297, 0.5274, 0.5254, 0.5235, 0.5217, & 0.5200, 0.5185, 0.5171, 0.5157, 0.5145, 0.5134, 0.5123, & 0.5114, 0.5104, 0.5096, 0.5088, 0.5081, 0.5075, 0.5069, & 0.5063, 0.5058, 0.5053, 0.5049, 0.5045, 0.5041, 0.5038, & 0.5000, 0.5000, 0.5000, 0.5000/ fqnd=sig*sig*d*geeinv k=min(1+fqnd*20.,80.) dknd=f2(k)+(f2(k+1)-f2(k))*(fqnd-f1(k))*20. cgnd=f3(k)+(f3(k+1)-f3(k))*(fqnd-f1(k))*20. if(fqnd.lt.0.2) dknd=sqrt(fqnd) wvnu=dknd/d c=sig/wvnu cg=cgnd*c return end C-------------------------------------------------------------------- function fsprm(beta,m,mmm,tht_wnd,thetam) c real*8 beta,tht_wnd,thetam,fsprm !RMarsooli save tht,fdat,dtht,ifirst dimension tht(35),fdat(35) data pi/3.141592654/ data ifirst/0/ if(ifirst.eq.0) then crm write(6,'(''spreading function'',/, crm & '' m tht fdat'')') dtht=pi/float(32) do k=1,35 fdat(k)=0. tht(k)=dtht*(k-1) if(k.lt.17) fdat(k)=0.5*beta/cosh(beta*tht(k))**2 crm write(6,'(i5,2f10.4)') k,tht(k),fdat(k) enddo endif ifirst=1 C Create spreading function around the wind direction fsprm=0. thtrel=abs(thetam-tht_wnd) if(thtrel.gt.pi) thtrel=abs(thtrel-2.*pi) kth=thtrel/dtht+1 c if(kth.le.14) fsprm=fdat(kth) & +(fdat(kth+1)-fdat(kth))*(thtrel-tht(kth))/dtht return end C------------------------------------------------------------ function fcoshinv(x) dimension xd(15),fdat(15) data xd/0.,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7./ data fdat/1.0,0.887,0.648,0.425,0.265,0.163,0.099,0.060, & 0.037,0.022,0.013,0.008,0.005,0.003,0.002/ fcoshinv=0.0 if(x.lt.7.0) then ixd=int(2.0*x)+1 fcoshinv=fdat(ixd)+(fdat(ixd+1)-fdat(ixd))*(x-xd(ixd))*2.0 endif return end C---------------------------------------------------------------- function fsinhinv(x) dimension xd(15),fdat(15) data xd/0.,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7./ data fdat/1.919,1.919,0.851,0.470,0.276,0.165,0.100,0.060, & 0.037,0.022,0.013,0.008,0.005,0.003,0.002/ fsinhinv=0.0 if(x.lt.7.0) then ixd=int(2.0*x)+1 fsinhinv=fdat(ixd)+(fdat(ixd+1)-fdat(ixd))*(x-xd(ixd))*2.0 endif return end C---------------------------------------------------------------- function mploc(a,b,i,j,im,jm,mmm) dimension a(im,jm,mmm),b(im,jm,mmm) amx=0. mmx=1 do m=2,mmm if(b(i,j,m).gt.0.10) then if(a(i,j,m).gt.amx) then amx=a(i,j,m) mmx=m endif endif enddo mploc=mmx return end C---------------------------------------------------------------- function mmloc(a,i,j,im,jm,mmm) dimension a(im,jm,mmm) amx=0. mmx=1 do m=2,mmm if(a(i,j,m).gt.amx) then amx=a(i,j,m) mmx=m endif c if(i.eq.16.and.j.eq.31)then c write(6,'(''mmloc'',4i5,2e12.3)') i,j,m,mmx,amx,a(i,j,m) c endif enddo mmloc=mmx return end C---------------------------------------------------------------- function summ(a,dpi,i,j,im,jm,mmm) dimension a(im,jm,mmm) summ=0. do m=2,mmm-1 summ=summ+a(i,j,m)*dpi c if(i.eq.3.and.j.eq.23) write(6,'(i5,e12.3)')m,summ enddo return end C---------------------------------------------------------------- function summ1(a,dpi,mmm) dimension a(mmm) summ1=0. do m=2,mmm-1 summ1=summ1+a(m)*dpi enddo return end C---------------------------------------------------------------- function sumk(a,dpi,kb) dimension a(kb) sumk=0. do k=2,kb-1 sumk=sumk+a(k)*dpi enddo return end C------------------------RMarsooli function radlimit (tempval,radval) c The purpos of this function is to increase the stability of sECOM c This function limits the change of radiation stresses between two c consecetive time steps not more than 20 percent real tempval,radval tempval=max(0.00000010,tempval) if(tempval.gt.1.2*radval) then radlimit=1.2*radval+0.00000010 elseif(tempval.lt.0.8*radval) then radlimit=0.8*radval+0.00000010 else radlimit=tempval endif end C---------------------------------------------------------------------- C END of WAVE CODE !WAVE C----------------------------------------------------------------------