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