1 | subroutine getslopes(ngrid,geopot) |
---|
2 | |
---|
3 | use comgeomfi_h, only: long, lati |
---|
4 | use slope_mod, only: theta_sl, psi_sl |
---|
5 | USE comcstfi_h |
---|
6 | implicit none |
---|
7 | |
---|
8 | #include "dimensions.h" |
---|
9 | |
---|
10 | ! This routine computes slope inclination and orientation for the GCM (callslope=.true. in callphys.def) |
---|
11 | ! It works fine with a non-regular grid for zoomed simulations. |
---|
12 | ! slope inclination angle (deg) 0 == horizontal, 90 == vertical |
---|
13 | ! slope orientation angle (deg) 0 == Northward, 90 == Eastward, 180 == Southward, 270 == Westward |
---|
14 | ! TN 04/1013 |
---|
15 | |
---|
16 | integer,intent(in) :: ngrid ! nnumber of atmospheric columns |
---|
17 | real,intent(in) :: geopot(ngrid) ! geopotential on phy grid |
---|
18 | real topogrid(iim,jjm+1) ! topography on lat/lon grid with poles and only one -180/180 point |
---|
19 | real latigrid(iim,jjm+1),longgrid(iim,jjm+1) ! meshgrid of latitude and longitude values (radians) |
---|
20 | real theta_val ! slope inclination |
---|
21 | real psi_val ! slope orientation |
---|
22 | real gradx(iim,jjm+1) ! x: latitude-wise topography gradient, increasing northward |
---|
23 | real grady(iim,jjm+1) ! y: longitude-wise topography gradient, increasing westward |
---|
24 | integer i,j,ig0 |
---|
25 | integer id2,idm1 ! a trick to compile testphys1d with debug option |
---|
26 | |
---|
27 | id2 = 2 |
---|
28 | idm1 = iim-1 |
---|
29 | |
---|
30 | ! rearrange topography on a 2d array |
---|
31 | do j=2,jjm |
---|
32 | ig0= 1+(j-2)*iim |
---|
33 | do i=1,iim |
---|
34 | topogrid(i,j)=geopot(ig0+i)/g |
---|
35 | latigrid(i,j)=lati(ig0+i) |
---|
36 | longgrid(i,j)=long(ig0+i) |
---|
37 | enddo |
---|
38 | enddo |
---|
39 | !poles : |
---|
40 | topogrid(:,1) = geopot(1)/g |
---|
41 | latigrid(:,1) = lati(1) |
---|
42 | longgrid(:,1) = long(1) |
---|
43 | topogrid(:,jjm+1) = geopot(ngrid)/g |
---|
44 | latigrid(:,jjm+1) = lati(ngrid) |
---|
45 | longgrid(:,jjm+1) = long(ngrid) |
---|
46 | |
---|
47 | |
---|
48 | |
---|
49 | ! compute topography gradient |
---|
50 | ! topogrid and rad are both in meters |
---|
51 | do j=2,jjm |
---|
52 | do i=1,iim |
---|
53 | gradx(i,j) = (topogrid(i,j+1) - topogrid(i,j-1)) / (latigrid(i,j+1)-latigrid(i,j-1)) |
---|
54 | gradx(i,j) = gradx(i,j) / rad |
---|
55 | enddo |
---|
56 | grady(1,j) = (topogrid(id2,j) - topogrid(iim,j)) / (2*pi+longgrid(id2,j)-longgrid(iim,j)) |
---|
57 | grady(1,j) = grady(1,j) / rad |
---|
58 | grady(iim,j) = (topogrid(1,j) - topogrid(idm1,j)) / (2*pi+longgrid(1,j)-longgrid(idm1,j)) |
---|
59 | grady(iim,j) = grady(iim,j) / rad |
---|
60 | do i=2,iim-1 |
---|
61 | grady(i,j) = (topogrid(i+1,j) - topogrid(i-1,j)) / (longgrid(i+1,j)-longgrid(i-1,j)) |
---|
62 | grady(i,j) = grady(i,j) / rad |
---|
63 | enddo |
---|
64 | enddo |
---|
65 | ! poles : |
---|
66 | gradx(:,1) = 0. |
---|
67 | grady(:,1) = 0. |
---|
68 | gradx(:,jjm+1) = 0. |
---|
69 | grady(:,jjm+1) = 0. |
---|
70 | |
---|
71 | |
---|
72 | |
---|
73 | ! compute slope inclination and orientation : |
---|
74 | theta_sl(:) = 0. |
---|
75 | psi_sl(:) = 0. |
---|
76 | do j=2,jjm |
---|
77 | do i=1,iim |
---|
78 | |
---|
79 | ig0= 1+(j-2)*iim |
---|
80 | |
---|
81 | theta_val=atan(sqrt( (gradx(i,j))**2 + (grady(i,j))**2 )) |
---|
82 | |
---|
83 | psi_val=0. |
---|
84 | if (gradx(i,j) .ne. 0.) psi_val= -pi/2. - atan(grady(i,j)/gradx(i,j)) |
---|
85 | if (gradx(i,j) .ge. 0.) psi_val= psi_val - pi |
---|
86 | psi_val = 3*pi/2. - psi_val |
---|
87 | psi_val = psi_val*180./pi |
---|
88 | psi_val = MODULO(psi_val,360.) |
---|
89 | |
---|
90 | theta_sl(ig0+i) = theta_val |
---|
91 | psi_sl(ig0+i) = psi_val |
---|
92 | |
---|
93 | enddo |
---|
94 | enddo |
---|
95 | |
---|
96 | |
---|
97 | |
---|
98 | end subroutine getslopes |
---|