MODULE hurricane_wind
IMPLICIT none
REAL, PARAMETER :: rearth = 6371.e3
REAL, PARAMETER :: rho_0 = 1024.e0
REAL, PARAMETER :: pi = 3.1415927
REAL, PARAMETER :: roa = 1.28
REAL, PARAMETER :: grav = 9.806e0
PRIVATE :: rearth, rho_0, pi, roa, grav
PRIVATE :: date2day, verinterp, interp1d, expwnd
CONTAINS
SUBROUTINE windmsg(ifplane,time_start,time,msgfile,im,jm,cor, east_u, &
north_u,east_v,north_v,wndux,wnduy,wndvx,wndvy)
USE ieee_arithmetic
INTEGER, PARAMETER :: pln = 100
INTEGER hour, lat, lon, mx, rmw
INTEGER readunit
INTEGER*4 startdate, date
INTEGER ni, i, j, n1, ipmax, nm1, np1, ii
INTEGER day,month,year,end_flag
INTEGER garb(5),rd1(4),rd2(4)
CHARACTER*19 name
CHARACTER*8 model_date
REAL,DIMENSION(pln) :: x,y,tm,pres,pres0,rmaxa,wspmax,radcls
REAL,DIMENSION(pln,4) :: r18v,r26v
REAL,DIMENSION(5) :: rref18v,rref26v,alphv
REAL,DIMENSION(14):: rad,ws,radm,wsm,angl
REAL :: timev
REAL cmp,t1,t2,f0,f1,l0,l1,r,a7,b,e,delp,x0,y0
REAL xcen, ycen, pres1, pres2, rcls, wsmax
REAL cl0, cnorth_e, x1, x2, y1, y2, utx, uty
REAL prsmin, wndmax, ws18, ws26
REAL wtan,wrad,rangl,alphaw,rmax
REAL wx, wy, wm, wind_scale
REAL rref18, rref26, wnd, utxa, utya
REAL deltax,deltax1,deltay,deltay1,dxdy,julday
CHARACTER(LEN=7),INTENT(IN) :: msgfile
CHARACTER(LEN=26), INTENT(IN) :: time_start
INTEGER, INTENT(IN) :: im, jm, ifplane
REAL, INTENT(IN) :: time
REAL, DIMENSION(im,jm), INTENT(IN) :: cor, east_u,north_u,east_v, &
north_v
REAL, DIMENSION(im,jm), INTENT(INOUT) :: wndux,wnduy,wndvx,wndvy
DATA rad/0.,.4,.7,.8,.95,1.,1.35, &
2.7,4.05,5.4,6.75,8.1,10.8,13.5/
DATA angl/0.,2.,4.,6.,7.,7.,14.,23.,24.,22.,21.,21.,21.,21./
rmax=50.e3
e=exp(1.)
timev=time
wind_scale = 1.0
model_date=time_start(3:4)//time_start(6:7)// &
time_start(9:10)//time_start(12:13)
read(model_date,'(i8)') startdate
readunit = 15
!----------------------- reading message file --------------------
17 FORMAT(a19,i6,1x,i4,1x,i3,1x,i4,1x,i3,1x,i3,3i5, &
1x,i2,1x,i3,1x,i4,1x,i4,1x,i4,1x,i4,1x,i4,1x,i4,1x,i4,1x,i4)
open(readunit,file=msgfile,status='old')
end_flag=0.
i=0
do while(end_flag.eq.0)
read(readunit,17) name,date,hour,lat,lon,garb,mx,rmw,rd1,rd2
!----------- falk 09-12-05 change output
! if(mod(iint,24).eq.0) then
! print*,'reading file ',filename
! write(6,17) name,date,hour,lat,lon,garb,mx,rmw,rd1,rd2
! end if
!rmy: uncomment the following code block for debug only
! if(my_task.eq.master_task) then
! write(*,*) 'name = ',name
! write(*,*) 'date = ',date
! write(*,*) 'hour = ',hour
! write(*,*) 'lat = ',lat
! write(*,*) 'lon = ',lon
! write(*,*) 'garb = ',garb
! write(*,*) 'mx = ',mx
! write(*,*) 'rmw = ',rmw
! write(*,*) 'rd1 = ',rd1
! write(*,*) 'rd2 = ',rd2
! end if
if(date.eq.0) then
end_flag = 1
close(readunit)
else
i=i+1
date=date*100+hour/100
CALL date2day(year,julday,date)
tm(i)=julday
x(i)=lon/10. !rmy: eastern hemsiphere
y(i)=lat/10. !rmy: northern hemisphere
if (east_u(1,1).lt.0.) x(i)=-x(i) !rmy: western hemisphere
if (north_u(1,1).lt.0.) y(i)=-y(i) !rmy: southern hemisphere
pres(i)=float(garb(3))
pres0(i)=float(garb(4))
radcls(i)=float(garb(5))*1000
wspmax(i)=float(mx)
rmaxa(i)=float(rmw)
do ni=1,4
n1=ni+(1-mod(ni,2))*sign(2,3-ni)
r18v(i,ni)=rd1(n1)
r26v(i,ni)=rd2(n1)
if(wspmax(i).le.26.or.r26v(i,ni).le.rmaxa(i)) &
r26v(i,ni)=-999
if(wspmax(i).le.18.or.r18v(i,ni).le.rmaxa(i)) &
r18v(i,ni)=-999
if(r26v(i,ni).gt.r18v(i,ni)) r26v(i,ni)=-999
end do
end if
end do
ipmax=i
close(readunit)
!--------------------- calculating starting day -----------------
CALL date2day(year,julday,startdate)
do i=1,ipmax
tm(i)=tm(i)-julday
end do
!++++++++++++++++++++++++++++++++++++++++++++++++++++++
! interpolation of hurricane path to determine the current position
cmp=tm(ipmax)
if(cmp.lt.time.or.time.lt.tm(1)) return
CALL verinterp(ipmax,1,tm,x,timev,f0)
CALL verinterp(ipmax,1,tm,y,timev,l0)
xcen=f0
ycen=l0
CALL verinterp(ipmax,1,tm,pres,timev,pres1)
CALL verinterp(ipmax,1,tm,pres0,timev,pres2)
CALL verinterp(ipmax,1,tm,radcls,timev,rcls)
delp=(pres2-pres1)*100.
if ( delp .lt. 100.0 ) delp = 100.0
CALL verinterp(ipmax,1,tm,wspmax,timev,wsmax)
prsmin=pres1
wndmax=wsmax
wsmax=wsmax*wind_scale
ws18=18*wind_scale
ws26=26*wind_scale
CALL verinterp(ipmax,1,tm,rmaxa,timev,rmax)
rmax=rmax*1.e3
do ni=1,4
CALL interp1d(-999.0,ipmax,1,tm,r18v(1,ni),timev,rref18v(ni))
if(rref18v(ni).ne.-999) rref18v(ni) = rref18v(ni)*1.e3
CALL interp1d(-999.0,ipmax,1,tm,r26v(1,ni),timev,rref26v(ni))
if(rref26v(ni).ne.-999) rref26v(ni) = rref26v(ni)*1.e3
alphv(ni) = (ni-1)*pi/2
end do
do ni=2,6
n1=mod(ni-1,4)+1
nm1=mod(ni-2,4)+1
np1=mod(ni,4)+1
if(rref18v(n1).eq.-999) then
if(rref18v(nm1).ne.-999) then
if(rref18v(np1).ne.-999) then
rref18v(n1)=0.5*(rref18v(nm1)+rref18v(np1))
else
rref18v(n1)=rref18v(nm1)
end if
else
if(rref18v(np1).ne.-999) then
rref18v(n1)=rref18v(np1)
else
rref18v(n1)=-999
end if
end if
end if
if(rref26v(n1).eq.-999) then
if(rref26v(nm1).ne.-999) then
if(rref26v(np1).ne.-999) then
rref26v(n1)=0.5*(rref26v(nm1)+rref26v(np1))
else
rref26v(n1)=rref26v(nm1)
end if
else
if(rref26v(np1).ne.-999) then
rref26v(n1)=rref26v(np1)
else
rref26v(n1)=-999
end if
end if
end if
end do
rref18v(5) = rref18v(1)
rref26v(5) = rref26v(1)
alphv(5) = alphv(4)+pi/2
!----------- falk 09-12-05 change output
!if(mod(iint,24).eq.0) then
! print*,'time=',time
! print*,'current hurricane position (x,y): ',f0,l0
! print*,'wsmax=',wsmax,'; delp=',delp,'; rmax=',rmax
! print*,'rref18v=',rref18v
! print*,'rref26v=',rref26v
! end if
x0=f0
y0=l0
f0=f0*2.*pi/360
l0=l0*2.*pi/360
cl0=cnorth_e*2.*pi/360
!--- f0,l0 are the carrent (lonitude,lattitude) coordinates of the hurricane
!--- calculating utx and uty (hurricane speed)
cmp=tm(ipmax)
do i=1,ipmax
if(abs(tm(i)-time).le.cmp) then
cmp=abs(tm(i)-time)
ii=i
end if
end do
!----------- falk 01-05-04 not to go out bnd
if((tm(ii)-time).le.0.and.ii.ne.ipmax) then
t1=tm(ii)
t2=tm(ii+1)
x1=x(ii)
x2=x(ii+1)
y1=y(ii)
y2=y(ii+1)
else
t2=tm(ii)
t1=tm(ii-1)
x2=x(ii)
x1=x(ii-1)
y2=y(ii)
y1=y(ii-1)
end if
if(ifplane.eq.1) then
deltax1=rearth*cos(cl0)*(x2-x1)*2.*pi/360
else
deltax1=rearth*cos(l0)*(x2-x1)*2.*pi/360
end if
deltay1=rearth*(y2-y1)*2.*pi/360
utx=deltax1/((t2-t1)*24.*3600.)
uty=deltay1/((t2-t1)*24.*3600.)
!--- calculating PARAMETERs for wind profile formula
b=wsmax**2*e*roa/delp
if ( b .gt. 3.0 ) b = 3.0
a7=rmax**b
do i=1,14
radm(i)=rmax*rad(i)
end do
do i=1,im
do j=1,jm
if(ieee_is_nan(east_u(i,j)).or. &
ieee_is_nan(north_u(i,j)).or. &
ieee_is_nan(east_v(i,j)) .or. &
ieee_is_nan(north_v(i,j)).or. &
f0.eq.0.or.l0.eq.0) then
print*,'Skipping wind generation in windmsg: &
f0 = , l0 = ',f0,l0
else
! calculating u-wind stress for i,j point of u-velocity grid
f1=east_u(i,j)*2.*pi/360
l1=north_u(i,j)*2.*pi/360
if(ifplane.eq.1) then
deltax=rearth*cos(cl0)*(f1-f0)
else
deltax=rearth*cos(l0)*(f1-f0)
end if
if(abs(deltax).lt.1.e-8) deltax=sign(1.e-8,deltax)
deltay=rearth*(l1-l0)
if(abs(deltay).lt.1.e-8) deltay=sign(1.e-8,deltay)
dxdy=deltax*deltay
r=sqrt(deltax**2+deltay**2)
alphaw=atan(abs(deltay/deltax))*sign(1.,dxdy) + &
(1-sign(1.,dxdy))*pi/2+(1-sign(1.,deltay))*pi/2
if(alphaw.ge.pi/4) then
alphaw=alphaw-pi/4
else
alphaw=alphaw-pi/4+2*pi
end if
CALL verinterp(5,1,alphv,rref18v,alphaw,rref18)
CALL verinterp(5,1,alphv,rref26v,alphaw,rref26)
CALL verinterp(14,1,radm,angl,r,rangl)
! calculating wind speed
if(r.ge.rcls*1.5) then
wnd=0
utxa=0
utya=0
else
if(rref18.le.0.and.rref26.le.0) then
wnd=sqrt(a7*b*delp*exp(-a7/r**b)/(roa*r**b)+r**2* &
cor(i,j)**2/4.)-r*cor(i,j)/2.
utxa=utx/2.
utya=uty/2.
else
wnd=expwnd(r,rmax,rref18,rref26,wsmax,ws18,ws26)
utxa=0.
utya=0.
end if
!rmy elect to make wind in f-plane runs perfectly axisymmetric
if(ifplane.eq.1) then
utxa=0.
utya=0.
end if
rangl=rangl*pi/180.
wtan=wnd*cos(rangl)
wrad=-wnd*sin(rangl)
wx = wrad*(deltax)/r-wtan*(deltay)/r
wy = wtan*(deltax)/r+wrad*(deltay)/r
wm = SQRT(wx*wx + wy*wy)
wind_scale = MAX(1.-MAX(wm-15.,0.)/25.*0.2,0.75)
wndux(i,j)=wndux(i,j) + wx*wind_scale + MIN(1.0, (rmax/r))*utx/2.0
wnduy(i,j)=wnduy(i,j) + wy*wind_scale + MIN(1.0, (rmax/r))*uty/2.0
endif
! calculating v-wind stress for i,j point of v-velocity grid
f1=east_v(i,j)*2.*pi/360
l1=north_v(i,j)*2.*pi/360
if(ifplane.eq.1) then
deltax=rearth*cos(cl0)*(f1-f0)
else
deltax=rearth*cos(l0)*(f1-f0)
end if
if(abs(deltax).lt.1.e-8) deltax=sign(1.e-8,deltax)
deltay=rearth*(l1-l0)
if(abs(deltay).lt.1.e-8) deltay=sign(1.e-8,deltay)
dxdy=deltax*deltay
r=sqrt(deltax**2+deltay**2)
alphaw=atan(abs(deltay/deltax))*sign(1.,dxdy) + &
(1-sign(1.,dxdy))*pi/2+(1-sign(1.,deltay))*pi/2
if(alphaw.ge.pi/4) then
alphaw=alphaw-pi/4
else
alphaw=alphaw-pi/4+2*pi
end if
CALL verinterp(5,1,alphv,rref18v,alphaw,rref18)
CALL verinterp(5,1,alphv,rref26v,alphaw,rref26)
CALL verinterp(14,1,radm,angl,r,rangl)
! calculating wind speed
if(r.ge.rcls*1.5) then
wnd=0
utxa=0
utya=0
else
if(rref18.le.0.and.rref26.le.0) then
wnd=sqrt(a7*b*delp*exp(-a7/r**b)/(roa*r**b)+r**2* &
cor(i,j)**2/4.)-r*cor(i,j)/2.
utxa=utx/2.
utya=uty/2.
else
wnd=expwnd(r,rmax,rref18,rref26,wsmax,ws18,ws26)
utxa=0.
utya=0.
end if
!rmy elect to make wind in f-plane runs perfectly axisymmetric
if(ifplane.eq.1) then
utxa=0.
utya=0.
end if
rangl=rangl*pi/180.
wtan=wnd*cos(rangl)
wrad=-wnd*sin(rangl)
wx = wrad*(deltax)/r-wtan*(deltay)/r
wy = wtan*(deltax)/r+wrad*(deltay)/r
wm = SQRT(wx*wx + wy*wy)
wind_scale = MAX(1.-MAX(wm-15.,0.)/25.*0.2,0.75)
wndvx(i,j)=wndvx(i,j) + wx*wind_scale + MIN(1.0, (rmax/r))*utx/2.0
wndvy(i,j)=wndvy(i,j) + wy*wind_scale + MIN(1.0, (rmax/r))*uty/2.0
end if
end if
end do
end do
return
END SUBROUTINE windmsg
!_______________________________________________________________________
SUBROUTINE windstr(im, jm, wndx, wndy, xcmp, ycmp, strcmp)
INTEGER :: i, j
LOGICAL, INTENT(IN) :: xcmp, ycmp
INTEGER, INTENT(IN) :: im, jm
REAL :: wm, z0w, cd, cd1, dcdl
REAL, DIMENSION(im,jm), INTENT(IN) :: wndx, wndy
REAL, DIMENSION(im,jm), INTENT(OUT) :: strcmp
IF ( xcmp .AND. ycmp ) THEN
PRINT*,' xcmp & ycmp can not have same value'
PRINT*,' Exiting windstr in hurricane_wind.f90'
RETURN
END IF
IF ( xcmp ) THEN
DO j = 1, jm
DO i = 1, im
wm = SQRT(wndx(i,j)*wndx(i,j) + wndy(i,j)*wndy(i,j))
IF (wm < 12.5) THEN
z0w = 0.0185/grav*(0.001*wm**2.+0.028*wm)**2.
cd = 0.4**2./(log(10./z0w))**2.
ELSE
cd1 = 3.58e-9*wm**3.-9.88e-7*wm**2.+7.81e-5*wm+7.9107e-4
dcdl = 1.12e-7*wm**2.+2.49e-5*wm-3.32e-4
cd = cd1-dcdl
END IF
strcmp(i,j) = -cd*roa*wm*wndx(i,j)/rho_0
END DO
END DO
ELSE IF( ycmp ) THEN
DO j = 1, jm
DO i = 1, im
wm = SQRT(wndx(i,j)*wndx(i,j) + wndy(i,j)*wndy(i,j))
IF (wm < 12.5) THEN
z0w = 0.0185/grav*(0.001*wm**2.+0.028*wm)**2.
cd = 0.4**2./(log(10./z0w))**2.
ELSE
cd1 = 3.58e-9*wm**3.-9.88e-7*wm**2.+7.81e-5*wm+7.9107e-4
dcdl = 1.12e-7*wm**2.+2.49e-5*wm-3.32e-4
cd = cd1-dcdl
END IF
strcmp(i,j) = -cd*roa*wm*wndy(i,j)/rho_0
END DO
END DO
END IF
END SUBROUTINE windstr
!_______________________________________________________________________
SUBROUTINE date2day(year,julday,date)
INTEGER*4 date
INTEGER ni
INTEGER dat2day(12),dat2dayl(12),day,month,year,hour
REAL julday
REAL*8 tmp
DATA dat2day/31,28,31,30,31,30,31,31,30,31,30,31/
DATA dat2dayl/31,29,31,30,31,30,31,31,30,31,30,31/
year=int(date/1000000.)
month=nint(100*(date/1000000.-int(date/1000000.)))
julday=0
if(mod(year,4).eq.0) then
do ni=1,month-1
julday=julday+dat2dayl(ni)
end do
else
do ni=1,month-1
julday=julday+dat2day(ni)
end do
end if
julday=julday+nint(100*(date/10000.-int(date/10000.)))
hour=date-nint(date/100.)*100
julday=julday+float(hour)/24.
return
END SUBROUTINE date2day
!_______________________________________________________________________
SUBROUTINE verinterp(nv,nvi,xv,yv,xvi,yvi)
!--------------------------------------------------------------
! this SUBROUTINE determines ni values yi at the points xi
! interpolating between the n values y at the points x
!--------------------------------------------------------------
INTEGER, INTENT(IN) :: nv, nvi
REAL, DIMENSION(nv), INTENT(IN) :: xv,yv
REAL, INTENT(IN) :: xvi
REAL, INTENT(OUT) :: yvi
REAL, PARAMETER :: maskv=-999.0
REAL tmpv(500)
INTEGER n,ivi
INTEGER iv,jv
if((xvi.gt.xv(nv).and.xvi.gt.xv(1)).or. &
(xvi.lt.xv(1).and.xvi.lt.xv(nv))) then
if(xvi.gt.xv(nv).and.xvi.gt.xv(1)) then
if(xv(nv).gt.xv(1)) then
yvi=yv(nv)
else
yvi=yv(1)
end if
else
if(xv(nv).gt.xv(1)) then
yvi=yv(1)
else
yvi=yv(nv)
end if
end if
else
do jv=1,nv-1
tmpv(jv)=(xvi-xv(jv))*(xvi-xv(jv+1))
end do
do jv=1,nv-1
if(tmpv(jv).le.0) then
ivi=jv
end if
end do
yvi=(yv(ivi)*abs(xv(ivi+1)-xvi)+yv(ivi+1)* &
abs(xvi-xv(ivi)))/abs(xv(ivi+1)-xv(ivi))
end if
return
END SUBROUTINE verinterp
!_______________________________________________________________________
SUBROUTINE interp1d(maskv,nv,nvi,xv,yv,xvi,yvi)
!--------------------------------------------------------------
! this SUBROUTINE determines ni values yi at the points xi
! interpolating between the n values y at the points x
! values equal to mask are ignored
!--------------------------------------------------------------
INTEGER, INTENT(IN) :: nv,nvi
REAL, DIMENSION(nv), INTENT(IN) :: xv, yv
REAL, INTENT(IN) :: maskv, xvi
REAL, INTENT(OUT) :: yvi
REAL tmpv(500)
INTEGER ivi
INTEGER iv,jv
do iv=1,nvi
if((xvi.gt.xv(nv).and.xvi.gt.xv(1)).or. &
(xvi.lt.xv(1).and.xvi.lt.xv(nv))) then
if(xvi.gt.xv(nv).and.xvi.gt.xv(1)) then
if(xv(nv).gt.xv(1)) then
yvi=yv(nv)
else
yvi=yv(1)
end if
else
if(xv(nv).gt.xv(1)) then
yvi=yv(1)
else
yvi=yv(nv)
end if
end if
else
do jv=1,nv-1
tmpv(jv)=(xvi-xv(jv))*(xvi-xv(jv+1))
end do
do jv=1,nv-1
if(tmpv(jv).le.0) ivi=jv
end do
if(yv(ivi).eq.maskv.or.yv(ivi+1).eq.maskv) then
yvi=maskv
else
yvi=(yv(ivi)*abs(xv(ivi+1)-xvi)+yv(ivi+1)* &
abs(xvi-xv(ivi)))/abs(xv(ivi+1)-xv(ivi))
end if
end if
end do
return
END SUBROUTINE interp1d
!_______________________________________________________________________
FUNCTION expwnd(r,rmax,rref18,rref26,wsmax,ws18,ws26)
REAL rref18,rref26,r,rmax,wsmax,ws18,ws26
REAL b,expwnd, r1, ws
r1=0.5*(rref18+rref26)
ws=0.5*(ws18+ws26)
if(rref18.le.0.) then
r1=rref26
ws=ws26
end if
if(rref26.le.0.) then
r1=rref18
ws=ws18
end if
if(r.ge.rmax) then
b=(rmax-r1)/log(ws/wsmax)
expwnd=wsmax*exp((rmax-r)/b)
else
expwnd=r*wsmax/rmax
end if
return
END FUNCTION expwnd
END MODULE