-
Notifications
You must be signed in to change notification settings - Fork 34
/
barotropic.f90
142 lines (120 loc) · 5.69 KB
/
barotropic.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
Module velocity
implicit none
contains
subroutine ubar(dat,outdat,z_w,z_wu,Nroms,II,JJ,xi_rho,eta_rho)
! ----------------------------------
! Program : ubar
!
! Trond Kristiansen, March 04 2009, bug fix 09.11.2009
! Rutgers University, NJ.
! -------------------------------------------------------------------------------------------------------
!
! USAGE: Compile this routine using Intel Fortran compiler and create
! a python module using the command:
! f2py --verbose --fcompiler=intelem -c -m barotropic barotropic.f90
! or
! f2py --verbose c -m barotropic barotropic.f90
! The resulting module is imported to python using:
! import barotropic
! To call the function from python use:
! barotropic.ubar(dat,outdat,zr,zw,Nroms,II,JJ)
!
! where: dat is the data such as temperature (3D structure (z,y,x))
!
! outdat is a 3D output array with the correct size (Nroms,JJ,II)
! zr is the depth matrix for the output grid (Nroms,JJ,II)
! zs is the 1D SODA depth-matrix (e.g. zs=[5,10,20,30])
! Nroms is the total depth levels in output grid
! JJ is the total grid points in eta direction
! II is the total grid points in xi direction
! -------------------------------------------------------------------------------------------------------
integer eta_rho, xi_rho, II, JJ, ic, jc, kc, Nroms
REAL(4), dimension(Nroms,JJ,II) :: dat
REAL(4), dimension(JJ,II) :: outdat
REAL(4), dimension(Nroms+1,eta_rho,xi_rho) :: z_w
REAL(4), dimension(Nroms+1,JJ,II) :: z_wu
!f2py intent(in,overwrite) dat,z_w,z_wu,Nroms, JJ, II, xi_rho, eta_rho
!f2py intent(in,out,overwrite) outdat
!f2py intent(hide) ic,jc,kc
print*,'---> Started ubar calculations'
! average z_w to Arakawa-C u,v-points (z_wu, z_wv)
do jc=1,JJ
do ic=2,II+1
do kc=1,Nroms+1
z_wu(kc,jc,ic-1) = 0.5*(z_w(kc,jc,ic-1)+z_w(kc,jc,ic))
end do
end do
end do
do jc=1,JJ
do ic=1,II
outdat(jc,ic)=0.0
do kc=1,Nroms
outdat(jc,ic) = outdat(jc,ic) + dat(kc,jc,ic)*abs(z_wu(kc+1,jc,ic) - z_wu(kc,jc,ic))
end do
if (abs(z_wu(1,jc,ic)) > 0.0) then
outdat(jc,ic) = outdat(jc,ic)/abs(z_wu(1,jc,ic))
else
outdat(jc,ic) = 0.0
end if
end do
end do
end subroutine ubar
subroutine vbar(dat,outdat,z_w,z_wv,Nroms,II,JJ,xi_rho,eta_rho)
! ----------------------------------
! Program : vbar
!
!
! Trond Kristiansen, March 04 2009, bug fix 09.11.2009
! Rutgers University, NJ.
! -------------------------------------------------------------------------------------------------------
!
! USAGE: Compile this routine using Intel Fortran compiler and create
! a python module using the command:
! f2py --verbose --fcompiler=intelem -c -m barotropic barotropic.f90
! or
! f2py --verbose --fcompiler=intelem -DF2PY_REPORT_ON_ARRAY_COPY=1 -c -m barotropic barotropic.f90
! The resulting module is imported to python using:
! import barotropic
! To call the function from python use:
! barotropic.ubar(dat,bathymetry,outdat,zr,zw,Nroms,II,JJ)
!
! where: dat is the data such as temperature (3D structure (z,y,x))
!
! outdat is a 3D output array with the correct size (Nroms,JJ,II)
! zr is the depth matrix for the output grid (Nroms,JJ,II)
! zs is the 1D SODA depth-matrix (e.g. zs=[5,10,20,30])
! Nroms is the total depth levels in output grid
! JJ is the total grid points in eta direction
! II is the total grid points in xi direction
! -------------------------------------------------------------------------------------------------------
integer eta_rho, xi_rho, II, JJ, ic, jc, kc, Nroms
REAL(4), dimension(Nroms,JJ,II) :: dat
REAL(4), dimension(JJ,II) :: outdat
REAL(4), dimension(Nroms+1,eta_rho,xi_rho) :: z_w
REAL(4), dimension(Nroms+1,JJ,II) :: z_wv
!f2py intent(in,overwrite) dat, z_w, z_wv, Nroms, JJ, II, xi_rho, eta_rho
!f2py intent(in,out,overwrite) outdat
!f2py intent(hide) ic,jc,kc
print*,'---> Started vbar calculations'
do jc=2,JJ+1
do ic=1,II
do kc=1,Nroms+1
z_wv(kc,jc-1,ic) = 0.5*(z_w(kc,jc-1,ic)+z_w(kc,jc,ic))
end do
end do
end do
do jc=1,JJ
do ic=1,II
outdat(jc,ic)=0.0
do kc=1,Nroms
outdat(jc,ic) = outdat(jc,ic) + dat(kc,jc,ic)*abs(z_wv(kc+1,jc,ic) - z_wv(kc,jc,ic))
end do
if (abs(z_wv(1,jc,ic)) > 0.0) then
outdat(jc,ic) = outdat(jc,ic)/abs(z_wv(1,jc,ic))
else
outdat(jc,ic) = 0.0
end if
end do
end do
end subroutine vbar
end module velocity