! Biju Thomas(GSO/URI) on 2016/05/15 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