-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy patheulerian_module.f90
132 lines (109 loc) · 3.45 KB
/
eulerian_module.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
module eulerian_module
use kind_module, only: i4b, dp
use planet_module, only: a=>planet_radius
use grid_module, only: &
nlon, nlat, ntrunc, &
gu, gv, gphi, sphi, sphi_old, lon, lat, coslatr, wind
use time_module, only: nstep, hstep, deltat, ifile, hfile
use legendre_transform_module, only: &
legendre_analysis, legendre_synthesis, &
legendre_synthesis_dlon, legendre_synthesis_dlat
use io_module, only: io_save
implicit none
private
real(kind=dp), dimension(:,:), allocatable, private :: dgphi
integer(kind=i4b), private :: nsave = 0
complex(kind=dp), dimension(:,:), allocatable, private :: sphi1
private :: update
public :: eulerian_init, eulerian_timeint, eulerian_clean
contains
subroutine eulerian_init()
implicit none
integer(kind=i4b) :: j
allocate(dgphi(nlon,nlat), sphi1(0:ntrunc,0:ntrunc))
print *, "Saving initial value"
call io_save(ifile, 1, gphi, "replace")
call io_save(ifile, 2, gu*a, "old")
call io_save(ifile, 3, gv*a, "old")
call legendre_synthesis(sphi,gphi)
print *, "step=0 t=0"
print *, "Saving step=0"
call io_save(hfile, 1, gphi, "replace")
call io_save(hfile, 2, gu*a, "old")
call io_save(hfile, 3, gv*a, "old")
nsave = 1
print *, "step=1/2", " t=", real(0.5*deltat)
call update(0.0d0*deltat,0.5d0*deltat)
print *, "step=1", " t=", real(deltat)
call update(0.5d0*deltat,deltat)
if (hstep==1) then
call io_save(hfile, 3*nsave+1, gphi, "old")
call io_save(hfile, 3*nsave+2, gu*a, "old")
call io_save(hfile, 3*nsave+3, gv*a, "old")
nsave = nsave + 1
end if
end subroutine eulerian_init
subroutine eulerian_clean()
deallocate(dgphi, sphi1)
end subroutine eulerian_clean
subroutine eulerian_timeint()
implicit none
integer(kind=i4b) :: i, j
do i=2, nstep
print *, "step=", i, " t=", real(i*deltat)
call update((i-1)*deltat,2.0d0*deltat)
if (mod(i,hstep)==0) then
print *, "Saving step=", i
call io_save(hfile, 3*nsave+1, gphi, "old")
call io_save(hfile, 3*nsave+2, gu*a, "old")
call io_save(hfile, 3*nsave+3, gv*a, "old")
nsave = nsave + 1
end if
end do
end subroutine eulerian_timeint
subroutine update(t,dt)
use time_module, only: etf, kappa
use uv_module, only: uv_nodiv, uv_div
implicit none
real(kind=dp), intent(in) :: t, dt
real(kind=dp) :: knt
integer(kind=i4b) :: i, j, m, n
select case(wind)
case("nodiv ")
call uv_nodiv(t,lon,lat,gu,gv)
case("div ")
call uv_div(t,lon,lat,gu,gv)
end select
call legendre_synthesis(sphi_old, gphi)
! dF/dlon
call legendre_synthesis_dlon(sphi, dgphi)
do j=1, nlat
dgphi(:,j) = dgphi(:,j)*coslatr(j)
end do
gphi = gphi - dt*gu*dgphi
! cos(lat)dF/dlat
call legendre_synthesis_dlat(sphi, dgphi)
do j=1, nlat
dgphi(:,j) = dgphi(:,j)*coslatr(j)
end do
gphi = gphi - dt*gv*dgphi
call legendre_analysis(gphi, sphi1)
! dissipation
do n=1, ntrunc
knt = kappa*dt*(n*(n+1.0_dp))**2
do m=0, n
sphi1(n,m) = (1.0_dp-knt)*sphi1(n,m)
end do
end do
! update
do m=0, ntrunc
sphi_old(m:ntrunc,m) = sphi(m:ntrunc,m) + &
etf * (sphi_old(m:ntrunc,m)-2.0_dp*sphi(m:ntrunc,m)+sphi1(m:ntrunc,m))
end do
do m=0, ntrunc
do n=m, ntrunc
sphi(n,m) = sphi1(n,m)
end do
end do
end subroutine update
end module eulerian_module