-
Notifications
You must be signed in to change notification settings - Fork 156
/
Copy pathhdiff.f90
145 lines (126 loc) · 3.9 KB
/
hdiff.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
133
134
135
136
137
138
139
140
141
142
143
144
145
subroutine hdiff(u_x,u_y,v_x,v_y,t_x,t_y,u_t,v_t,t_t,mype)
!$$$ subprogram documentation block
! . . . .
! subprogram: hdiff calculate horiziontal diffusion
! prgmmr: rancic
!
! abstract: compute horizontal diffusion for wind, and virtual
! temperature
!
! program history log:
! 2010-02-25 rancic
!
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use constants, only: zero,half,one,two
use gridmod, only: lat2,lon2,nsig
use tends4pertmod, only: time_step
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: mype
real(r_kind),dimension(lat2,lon2,nsig),intent(in ) :: u_x,u_y,v_x,v_y,t_x,t_y
real(r_kind),dimension(lat2,lon2,nsig),intent(inout) :: u_t,v_t,t_t
! Declare local variables
real(r_kind),dimension(lat2,lon2,nsig):: uc_x,uc_y,vc_x,vc_y,tc_x,tc_y
real(r_kind),dimension(lat2,lon2,nsig):: u_xx,u_yy,v_xx,v_yy,t_xx,t_yy
real(r_kind) uave,vave,up,vp
real(r_kind) khdff,kht,facdec
integer(i_kind) i,j,k
khdff=1.73e04_r_kind
kht=khdff*time_step
t_xx=zero; t_yy=zero
u_xx=zero; u_yy=zero
v_xx=zero; v_yy=zero
do k=1,nsig
if(k<nsig-2) then
facdec=exp(-3.*(k-1)/(nsig-1))
else
facdec=one
end if
facdec=facdec*kht
do j=1,lon2
do i=1,lat2
uc_x(i,j,k)=u_x(i,j,k)*facdec
vc_x(i,j,k)=v_x(i,j,k)*facdec
tc_x(i,j,k)=t_x(i,j,k)*facdec
uc_y(i,j,k)=u_y(i,j,k)*facdec
vc_y(i,j,k)=v_y(i,j,k)*facdec
tc_y(i,j,k)=t_y(i,j,k)*facdec
end do
end do
end do
call mp_compact_dlon2(uc_x,u_xx,.false.,nsig,mype)
call mp_compact_dlon2(vc_x,v_xx,.false.,nsig,mype)
call mp_compact_dlon2(tc_x,t_xx,.false.,nsig,mype)
call mp_compact_dlat2(uc_y,u_yy,.false.,nsig,mype)
call mp_compact_dlat2(vc_y,v_yy,.false.,nsig,mype)
call mp_compact_dlat2(tc_y,t_yy,.false.,nsig,mype)
do k=1,nsig
do j=1,lon2
do i=1,lat2
u_t(i,j,k)=u_t(i,j,k)+u_xx(i,j,k)+u_yy(i,j,k)
v_t(i,j,k)=v_t(i,j,k)+v_xx(i,j,k)+v_yy(i,j,k)
t_t(i,j,k)=t_t(i,j,k)+t_xx(i,j,k)+t_yy(i,j,k)
end do
end do
end do
end subroutine hdiff
subroutine hdiff_ad(u_x,u_y,v_x,v_y,t_x,t_y,u_t,v_t,t_t,mype)
use kinds,only: r_kind,i_kind
use constants, only: zero,half,one,two
use gridmod, only: lat2,lon2,nsig
use tends4pertmod, only: time_step
implicit none
! Declare passed variables
integer(i_kind) ,intent(in ) :: mype
real(r_kind),dimension(lat2,lon2,nsig),intent( out) :: u_x,u_y,v_x,v_y,t_x,t_y
real(r_kind),dimension(lat2,lon2,nsig),intent(inout) :: u_t,v_t,t_t
! Declare local variables
integer(i_kind) i,j,k
real(r_kind),dimension(lat2,lon2,nsig):: uc_x,uc_y,vc_x,vc_y,tc_x,tc_y
real(r_kind),dimension(lat2,lon2,nsig):: u_xx,u_yy,v_xx,v_yy,t_xx,t_yy
real(r_kind) khdff,kht,facdec
!
! Preliminaries
!
khdff=1.73e04_r_kind
kht=khdff*time_step
!
! Start adjoint
!
t_xx=t_t ; t_yy=t_t
u_xx=u_t ; u_yy=u_t
v_xx=v_t ; v_yy=v_t
tc_x=zero ; tc_y=zero
uc_x=zero ; uc_y=zero
vc_x=zero ; vc_y=zero
call mp_compact_dlon2_ad(tc_x,t_xx,.false.,nsig,mype)
call mp_compact_dlon2_ad(vc_x,v_xx,.false.,nsig,mype)
call mp_compact_dlon2_ad(uc_x,u_xx,.false.,nsig,mype)
call mp_compact_dlat2_ad(tc_y,t_yy,.false.,nsig,mype)
call mp_compact_dlat2_ad(vc_y,v_yy,.false.,nsig,mype)
call mp_compact_dlat2_ad(uc_y,u_yy,.false.,nsig,mype)
do k=nsig,1,-1
if(k<nsig-2) then
facdec=exp(-3.*(k-1)/(nsig-1))
else
facdec=one
end if
facdec=facdec*kht
do j=1,lon2
do i=1,lat2
u_x(i,j,k)=uc_x(i,j,k)*facdec
v_x(i,j,k)=vc_x(i,j,k)*facdec
t_x(i,j,k)=tc_x(i,j,k)*facdec
u_y(i,j,k)=uc_y(i,j,k)*facdec
v_y(i,j,k)=vc_y(i,j,k)*facdec
t_y(i,j,k)=tc_y(i,j,k)*facdec
end do
end do
end do
end subroutine hdiff_ad