| 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 | integer id2,idm1 ! a trick to compile testphys1d with debug option |
|---|
| 28 | |
|---|
| 29 | id2 = 2 |
|---|
| 30 | idm1 = iim-1 |
|---|
| 31 | |
|---|
| 32 | ! rearrange topography on a 2d array |
|---|
| 33 | do j=2,jjm |
|---|
| 34 | ig0= 1+(j-2)*iim |
|---|
| 35 | do i=1,iim |
|---|
| 36 | topogrid(i,j)=geopot(ig0+i)/g |
|---|
| 37 | latigrid(i,j)=lati(ig0+i) |
|---|
| 38 | longgrid(i,j)=long(ig0+i) |
|---|
| 39 | enddo |
|---|
| 40 | enddo |
|---|
| 41 | !poles : |
|---|
| 42 | topogrid(:,1) = geopot(1)/g |
|---|
| 43 | latigrid(:,1) = lati(1) |
|---|
| 44 | longgrid(:,1) = long(1) |
|---|
| 45 | topogrid(:,jjm+1) = geopot(ngridmx)/g |
|---|
| 46 | latigrid(:,jjm+1) = lati(ngridmx) |
|---|
| 47 | longgrid(:,jjm+1) = long(ngridmx) |
|---|
| 48 | |
|---|
| 49 | |
|---|
| 50 | |
|---|
| 51 | ! compute topography gradient |
|---|
| 52 | ! topogrid and rad are both in meters |
|---|
| 53 | do j=2,jjm |
|---|
| 54 | do i=1,iim |
|---|
| 55 | gradx(i,j) = (topogrid(i,j+1) - topogrid(i,j-1)) / (latigrid(i,j+1)-latigrid(i,j-1)) |
|---|
| 56 | gradx(i,j) = gradx(i,j) / rad |
|---|
| 57 | enddo |
|---|
| 58 | grady(1,j) = (topogrid(id2,j) - topogrid(iim,j)) / (2*pi+longgrid(id2,j)-longgrid(iim,j)) |
|---|
| 59 | grady(1,j) = grady(1,j) / rad |
|---|
| 60 | grady(iim,j) = (topogrid(1,j) - topogrid(idm1,j)) / (2*pi+longgrid(1,j)-longgrid(idm1,j)) |
|---|
| 61 | grady(iim,j) = grady(iim,j) / rad |
|---|
| 62 | do i=2,iim-1 |
|---|
| 63 | grady(i,j) = (topogrid(i+1,j) - topogrid(i-1,j)) / (longgrid(i+1,j)-longgrid(i-1,j)) |
|---|
| 64 | grady(i,j) = grady(i,j) / rad |
|---|
| 65 | enddo |
|---|
| 66 | enddo |
|---|
| 67 | ! poles : |
|---|
| 68 | gradx(:,1) = 0. |
|---|
| 69 | grady(:,1) = 0. |
|---|
| 70 | gradx(:,jjm+1) = 0. |
|---|
| 71 | grady(:,jjm+1) = 0. |
|---|
| 72 | |
|---|
| 73 | |
|---|
| 74 | |
|---|
| 75 | ! compute slope inclination and orientation : |
|---|
| 76 | theta_sl(:) = 0. |
|---|
| 77 | psi_sl(:) = 0. |
|---|
| 78 | do j=2,jjm |
|---|
| 79 | do i=1,iim |
|---|
| 80 | |
|---|
| 81 | ig0= 1+(j-2)*iim |
|---|
| 82 | |
|---|
| 83 | theta_val=atan(sqrt( (gradx(i,j))**2 + (grady(i,j))**2 )) |
|---|
| 84 | |
|---|
| 85 | psi_val=0. |
|---|
| 86 | if (gradx(i,j) .ne. 0.) psi_val= -pi/2. - atan(grady(i,j)/gradx(i,j)) |
|---|
| 87 | if (gradx(i,j) .ge. 0.) psi_val= psi_val - pi |
|---|
| 88 | psi_val = 3*pi/2. - psi_val |
|---|
| 89 | psi_val = psi_val*180./pi |
|---|
| 90 | psi_val = MODULO(psi_val,360.) |
|---|
| 91 | |
|---|
| 92 | theta_sl(ig0+i) = theta_val |
|---|
| 93 | psi_sl(ig0+i) = psi_val |
|---|
| 94 | |
|---|
| 95 | enddo |
|---|
| 96 | enddo |
|---|
| 97 | |
|---|
| 98 | |
|---|
| 99 | |
|---|
| 100 | end subroutine getslopes |
|---|
