! Data module for 2-D shallow water model on a staggered grid ! ! | | | ! h - u - h - u - h - ! | | | ! ↑ v q v q v ! y | | | ! h - u - h - u - h - ! x --> ! module sw2_data real, parameter :: pi = 3.14159, xmax = 1.015625, ymax = 1.015625 integer, parameter :: imax=65, jmax=65 integer, parameter :: mq=1, mu=2, mv=3, mp=4 real, parameter :: dgx=xmax/imax, dgy=ymax/jmax real, dimension(imax, jmax) :: h,u,v,f real, dimension(imax) :: xq = (/ (dgx*i+dgx/2., i=0, imax-1 ) /), & xh = (/ (dgx*i, i=0, imax-1 ) /) real, dimension(jmax) :: yq = (/ (dgy*j+dgy/2., j=0, jmax-1 ) /), & yh = (/ (dgy*j, j=0, jmax-1 ) /) logical, dimension(imax,jmax,4) :: mask real :: p_time ! present time. real :: f0, beta, amp, nu, delta contains !------------------- set up -------------------------- subroutine sw_setup(file_name) character(len=*) :: file_name character(len=16) :: cformat open(1,file=file_name) write(cformat, '(A,I4,A)') '(', imax, 'L1)' do j=1, jmax read (1,cformat) mask(:,j,mq) write(*,cformat) mask(:,j,mq) end do close(1) mask(:,:,mu) = mask(:,:,mq) .or. cshift(mask(:,:,mq),-1,2) mask(:,:,mv) = mask(:,:,mq) .or. cshift(mask(:,:,mq),-1,1) mask(:,:,mp) = mask(:,:,mu) .or. cshift(mask(:,:,mu),-1,1) do j=1, jmax f(:,j) = f0 + beta*(yq(j) - ymax/2.) end do end subroutine sw_setup !------------------ average_x ------------------------ function AX (m, var) real, dimension(imax,jmax) :: var, AX if(m==mp .or. m==mv) then do i=1, imax i2 = i+1 if(i==imax) i2=1 do j=1, jmax if(mask(i,j,m-2)) then AX(i,j) = (var(i2,j) + var(i,j))/2. else AX(i,j) = 0. end if end do end do else do i=1, imax i2 = i-1 if(i==1) i2=imax do j=1, jmax if(mask(i,j,m).and.mask(i2,j,m)) then AX(i,j) = (var(i,j) + var(i2,j))/2. else if(mask(i,j,m)) then AX(i,j) = var(i,j) else if(mask(i2,j,m)) then AX(i,j) = var(i2,j) else AX(i,j) = 0. end if end do end do end if end function AX !------------------ average_y ------------------------ function AY (m, var) real, dimension(imax,jmax) :: var, AY if(m==mp .or. m==mu) then do j=1, jmax j2 = j+1 if(j==jmax) j2=1 do i=1, imax if (mask(i,j,m-1)) then AY(i,j) = (var(i,j2) + var(i,j))/2. else AY(i,j) = 0. end if end do end do else do j=1, jmax j2 = j-1 if(j==1) j2=jmax do i=1, imax if (mask(i,j,m).and.mask(i,j2,m)) then AY(i,j) = (var(i,j) + var(i,j2))/2. else if(mask(i,j,m)) then AY(i,j) = var(i,j) else if(mask(i,j2,m)) then AY(i,j) = var(i,j2) else AY(i,j) = 0. end if end do end do end if end function AY !------------------ derivative_x ------------------------ function DX (m, var) real, dimension(imax,jmax) :: var, DX if(m==mp .or. m==mv) then do i=1, imax i2 = i+1 if(i==imax) i2=1 do j=1, jmax if(mask(i,j,m-2)) then DX(i,j) = (var(i2,j) - var(i,j))/dgx else DX(i,j) = 0. end if end do end do else do i=1, imax i2 = i-1 if(i==1) i2=imax do j=1, jmax if(mask(i,j,m).and.mask(i2,j,m)) then DX(i,j) = (var(i,j) - var(i2,j))/dgx else if(mask(i,j,m)) then DX(i,j) = 2*var(i,j)/dgx else if(mask(i2,j,m)) then DX(i,j) = -2*var(i2,j)/dgx else DX(i,j) = 0. end if end do end do end if end function DX !------------------ derivative_y ------------------------ function DY (m, var) real, dimension(imax,jmax) :: var, DY if(m==mp .or. m==mu) then do j=1, jmax j2 = j+1 if(j==jmax) j2=1 do i=1, imax if (mask(i,j,m-1)) then DY(i,j) = (var(i,j2) - var(i,j))/dgy else DY(i,j) = 0. end if end do end do else do j=1, jmax j2 = j-1 if(j==1) j2=jmax do i=1, imax if (mask(i,j,m).and.mask(i,j2,m)) then DY(i,j) = (var(i,j) - var(i,j2))/dgy else if(mask(i,j,m)) then DY(i,j) = 2*var(i,j)/dgy else if(mask(i,j2,m)) then DY(i,j) = -2*var(i,j2)/dgy else DY(i,j) = 0. end if end do end do end if end function DY !------------------ divergence ------------------------ function DIV(u,v) real, dimension(imax,jmax) :: DIV, u, v do j=1, jmax j2 = j-1 if(j==1) j2=jmax do i=1, imax i2 = i-1 if(i==1) i2=imax dd = 0. area = 0. if(mask(i,j,mq)) then dd = dd + u(i,j)/dgx + v(i,j)/dgy area = area + 0.5 end if if(mask(i,j2,mq)) then dd = dd + u(i,j)/dgx - v(i,j2)/dgy area = area + 0.5 end if if(mask(i2,j,mq)) then dd = dd - u(i2,j)/dgx + v(i,j)/dgy area = area + 0.5 end if if(mask(i2,j2,mq)) then dd = dd - u(i2,j)/dgx - v(i,j2)/dgy area = area + 0.5 end if if(area/=0) then DIV(i,j) = dd/area else DIV(i,j) = 0. end if end do end do end function DIV end module !----------------------------------------------------------------- ! Shallow water equation module. ! Potential-enstrophy-conserving scheme by Sadourny (1975) JAS !----------------------------------------------------------------- module sw2_eq use sw2_data real, dimension(imax,jmax) :: dh1, dh2, dh3, dh4, & & UF, UL, VF, VL, du1, du2, du3, du4, & & Q, B, HH, dv1, dv2, dv3, dv4 real :: dt integer :: int_count = 0 contains !------------- Runge-Kutta ------------ subroutine runge_kutta call time_derivative(u , v , h , du1, dv1, dh1) call time_derivative(u+du1/2., v+dv1/2., h+dh1/2., du2, dv2, dh2) call time_derivative(u+du2/2., v+dv2/2., h+dh2/2., du3, dv3, dh3) call time_derivative(u+du3 , v+dv3 , h+dh3 , du4, dv4, dh4) u = u + (du1/2. + du2 + du3 + du4/2.) / 3. v = v + (dv1/2. + dv2 + dv3 + dv4/2.) / 3. h = h + (dh1/2. + dh2 + dh3 + dh4/2.) / 3. int_count = int_count + 1 p_time = dt*int_count end subroutine !------------- Shallow Water Eq. ------------ subroutine time_derivative(u, v, h, du, dv, dh) real, dimension(imax,jmax) :: u, v, h, du, dv, dh UF = (1.+ AX(mp,h))*u ! Volume Transport (x) VF = (1.+ AY(mp,h))*v ! Volume Transport (y) B = h + (AX(mu,u*u) + AY(mv,v*v))/2. ! Bernoulli Function Q = (DX(mv,v)-DY(mu,u)+f) / (1.+AX(mv,AY(mp,h))) ! Potential Vorticity UL = nu*(DX(mp,DX(mu,u)) + DY(mq,DY(mu,u))) VL = nu*(DX(mq,DX(mv,v)) + DY(mp,DY(mv,v))) du = -dt * (-AY(mq,Q)*AY(mq,AX(mv,VF)) + DX(mp,B) - UL) dv = -dt * ( AX(mq,Q)*AX(mq,AY(mu,UF)) + DY(mp,B) - VL) dh = -dt * DIV(UF,VF) ! for initial problem ! dh = -dt *( DIV(UF,VF) - HH) ! for steady forcing end subroutine end module !----------------------------------------------------------------- ! Shallow water model (Fortran90 style) !----------------------------------------------------------------- program main use sw2_eq use dcl character(LEN=32) :: cfile !----------parameter SET --------- f0=20. beta=20. delta = 0. nu = 1.e-4 dt = 0.01 amp = 1.e-3 fac = 0.006 ! for beta=1 ! fac = 0.03 ! for beta=1 ! fac = 0.06 ! for beta=20 ! fac = 0.3 ! for beta=20 thermal call sw_setup('c1.dat') call InitH ! for wind driven ! call InitT ! for thermal forcing u = -AY(mq,AX(mv,DY(mp,h))/f) v = AX(mq,AY(mu,DX(mp,h))/f) ! call DclSetParm('IFONT' , 2) ! call DclSetParm('WINDOW_WIDTH' , 800) ! call DclSetParm('WINDOW_HEIGHT' , 600) call DclSetParm('WAIT' , .false.) call DclSetParm('WAIT_OPENING' , .false.) call DclSetParm('ALTERNATE' , .true.) call DclSetParm('USE_FULL_WINDOW' , .true.) call DclSetParm('DUMP' , .true.) call DclOpenGraphics(1) call DclSetAspectRatio(1.0, 0.75) call DclSetAxisFactor(1.0) call get_time(time0) call DrawGraph do i=1, 200 ! 200 for beta=20 1000 for beta=1 do ii = 1, 100 call runge_kutta if(mod(ii-1,50).eq.0) write(*,'(/,i3,\)') i if(mod(ii-1,1).eq.0) write(*,'(a,\)') '.' end do write(*,'(/,a,\)') ' ' call DrawGraph call get_time(time1) time = time1-time0 time0 = time1 tmass = sum(ax(mv,ay(mp,h)))*dgx*dgy write(*,*) 'time', time, ' mass=', tmass end do call DclCloseGraphics contains !----------------------------------------------------------------- subroutine get_time(time) character*10 ctime call date_and_time(time=ctime) read(ctime, '(2i2, f6.3)') ihour, minu, sec time = ihour*3600. + minu*60. + sec end subroutine !----------------------------------------------------------------- subroutine InitT HH=0. do j=0, 10 HH(:,jmax-j) = amp*exp(-dgy*f(1,jmax)*j) end do !--------------------- ! h = HH ! for initial problem !--------------------- hav= SUM(AX(mu,AY(mp,HH)))/count(mask(:,:,1)) HH = HH-hav HH = HH*fac ! for steady problem end subroutine !----------------------------------------------------------------- subroutine InitH i1 = jmax/4 i2 = jmax/4*3+1 j1 = jmax/4 j2 = jmax/4*3+1 iw = jmax/16 HH = 0. do j=1, jmax do i=1, imax if(i1+iw.lt.i .and. i.lt.i2-iw) then cx = 1. else if(i.lt.i1-iw .or. i2+iw.lt.i) then cx = 0. else if(i.le.i1+iw) then x = (i-(i1-iw))/iw*pi/2 cx = (1.-cos(x))/2. else if(i.ge.i2-iw) then x = ((i2+iw)-i)/iw*pi/2 cx = (1.-cos(x))/2. end if if(j1+iw.lt.j .and. j.lt.j2-iw) then cy = 1. else if(j.lt.j1-iw .or. j2+iw.lt.j) then cy = 0. else if(j.le.j1+iw) then y = (j-(j1-iw))/iw*pi/2 cy = (1.-cos(y))/2. else if(j.ge.j2-iw) then y = ((j2+iw)-j)/iw*pi/2 cy = (1.-cos(y))/2. end if HH(i,j) = amp*min(cx,cy) end do end do !--------------------- h = HH ! for initial problem !--------------------- ! hav= SUM(AX(mu,AY(mp,HH)))/count(mask(:,:,1)) ! HH = HH-hav ! HH = HH*fac ! for steady problem end subroutine !----------------------------------------------------------------- subroutine DrawGraph character(len=32) :: ctime real, dimension(imax,jmax):: var call DclNewFrame CALL UWSGXZ(.FALSE.) ! 将来的にNewFigに含まれる予定。 CALL UWSGYZ(.FALSE.) ! 将来的にNewFigに含まれる予定。 call DclSetWindow ( 0., 1., 0., 1.) call DclSetViewPort ( 0.1, .75, 0.05, 0.7) call DclSetTransFunction ! vmin = -0.6*amp ! vmax = 0.6*amp ! vmin = -0.4*amp ! vmax = 1.*amp ! vmin = -0.6*amp ! vmax = 0.8*amp vmin = 0.*amp vmax = 1.*amp var = h call DclSetParm ('INTERPRET_MISSING_VALUE', .true.) call DclGetParm ('MISSING_REAL', RMISS) where(.not.mask(:,:,mp)) var = RMISS call DclSetColorFile('rainbow.cmap') call DclSetColorRange (vmin, vmax) call DclSetParm ('CYCLE_COLOR', .false.) call DclSetMaskRange (0., amp/10.) call DclSetParm ('MASK_COLOR' , .true.) call DclPaintData (var) call DclSetParm ('DRAW_RIGHT_LABEL', .false.) call DclDrawScaledAxis write(ctime,'(A,F6.2)') 'TIME=', p_time call DclDrawTextNormalized(0.1, 0.72, ctime, height=0.025, centering=-1, index=3) CALL UWSGXZ(.FALSE.) ! 将来的にNewFigに含まれる予定。 CALL UWSGYZ(.FALSE.) ! 将来的にNewFigに含まれる予定。 call DclNewFig call DclSetParm ('DRAW_RIGHT_LABEL', .true.) call DclPaintColorBarY (0.80, 0.83, 0.1, 0.65, vmin, vmax, 'r') end subroutine end