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