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