source: trunk/LMDZ.MARS/libf/phymars/slope_mod.F90

Last change on this file was 3203, checked in by jbclement, 11 months ago

Mars PCM:

  • Addition of the possibility to use subslopes parametization in 1D. To do so, put 'nslope=x' in "run.def" where x (1, 3, 5 or 7) is the number of slopes you want to, or modify it in an appropriate "startfi.nc". Then, default subslopes definition and distribution are used by the model. The already existing case of using one slope whose inclination and orientation are defined in "run.def" is still possible.
  • Some cleanings throughout the code, in particular related to unused variables.

JBC

File size: 9.1 KB
RevLine 
[3183]1MODULE slope_mod
[1047]2
3implicit none
4
[3183]5real, save, dimension(:), allocatable :: theta_sl ! slope angle versus horizontal (deg)
6real, save, dimension(:), allocatable :: psi_sl   ! slope orientation (deg)
[1047]7
[2578]8!$OMP THREADPRIVATE(theta_sl,psi_sl)
9
[3183]10!=======================================================================
[1047]11contains
[3183]12!=======================================================================
[1047]13
[3183]14SUBROUTINE getslopes(ngrid,geopot)
[1770]15
[3183]16use geometry_mod,       only: longitude, latitude ! in radians
17use comcstfi_h,         only: g, rad, pi
18use mod_phys_lmdz_para, only: is_parallel
19use mod_grid_phy_lmdz,  only: nbp_lon, nbp_lat
[1770]20
[3183]21implicit none
[1770]22
[3183]23! This routine computes slope inclination and orientation for the GCM (callslope=.true. in callphys.def)
24! It works fine with a non-regular grid for zoomed simulations.
25! slope inclination angle (deg)  0 == horizontal, 90 == vertical
26! slope orientation angle (deg)  0 == Northward,  90 == Eastward, 180 == Southward, 270 == Westward
27! TN 04/1013
[1770]28
[3183]29! Input arguments
30integer,                intent(in) :: ngrid  ! nnumber of atmospheric columns
31real, dimension(ngrid), intent(in) :: geopot ! geopotential on phy grid
[1770]32
[3183]33! Local variables
34real, dimension(nbp_lon,nbp_lat) :: topogrid           ! topography on lat/lon grid with poles and only one -180/180 point
35real, dimension(nbp_lon,nbp_lat) :: latigrid, longgrid ! meshgrid of latitude and longitude values (radians)
36real, dimension(nbp_lon,nbp_lat) :: gradx              ! x: latitude-wise topography gradient,  increasing northward
37real, dimension(nbp_lon,nbp_lat) :: grady              ! y: longitude-wise topography gradient, increasing westward
38real                             :: theta_val          ! slope inclination
39real                             :: psi_val            ! slope orientation
40integer                          :: i, j, ig0
41integer                          :: id2, idm1          ! a trick to compile testphys1d with debug option
42
43if (is_parallel) then
44    ! This routine only works in serial mode so stop now.
45    write(*,*) "getslopes Error: this routine is not designed to run in parallel"
46    call abort_physic("getslopes",'cannot be run in parallel',1)
47endif
48
49id2  = 2
50idm1 = nbp_lon-1
51
52! rearrange topography on a 2d array
53do j = 2,nbp_lat-1
54    ig0 = 1 + (j - 2)*nbp_lon
55    do i = 1,nbp_lon
56        topogrid(i,j) = geopot(ig0 + i)/g
57        latigrid(i,j) = latitude(ig0 + i)
58        longgrid(i,j) = longitude(ig0 + i)
59    enddo
60enddo
61
62! poles:
63topogrid(:,1) = geopot(1)/g
64latigrid(:,1) = latitude(1)
65longgrid(:,1) = longitude(1)
66topogrid(:,nbp_lat) = geopot(ngrid)/g
67latigrid(:,nbp_lat) = latitude(ngrid)
68longgrid(:,nbp_lat) = longitude(ngrid)
69
70! compute topography gradient
71! topogrid and rad are both in meters
72do j = 2,nbp_lat - 1
73    do i=1,nbp_lon
74        gradx(i,j) = (topogrid(i,j + 1) - topogrid(i,j - 1))/(latigrid(i,j + 1)-latigrid(i,j - 1))
75        gradx(i,j) = gradx(i,j)/rad
76    enddo
77    grady(1,j) = (topogrid(id2,j) - topogrid(nbp_lon,j))/(2*pi + longgrid(id2,j) - longgrid(nbp_lon,j))
78    grady(1,j) = grady(1,j) / rad
79    grady(nbp_lon,j) = (topogrid(1,j) - topogrid(idm1,j))/(2*pi + longgrid(1,j) - longgrid(idm1,j))
80    grady(nbp_lon,j) = grady(nbp_lon,j)/rad
81    do i = 2,nbp_lon - 1
82        grady(i,j) = (topogrid(i + 1,j) - topogrid(i-1,j))/(longgrid(i + 1,j) - longgrid(i - 1,j))
83        grady(i,j) = grady(i,j)/rad
84    enddo
85enddo
86
87! poles:
88gradx(:,1) = 0.
89grady(:,1) = 0.
90gradx(:,nbp_lat) = 0.
91grady(:,nbp_lat) = 0.
92
93! compute slope inclination and orientation:
94theta_sl = 0.
95psi_sl  = 0.
96do j = 2,nbp_lat - 1
97    do i = 1,nbp_lon
98        ig0 = 1 + (j - 2)*nbp_lon
99
100        theta_val = atan(sqrt((gradx(i,j))**2 + (grady(i,j))**2))
101
102        psi_val = 0.
103        if (gradx(i,j) /= 0.) psi_val = -pi/2. - atan(grady(i,j)/gradx(i,j))
104        if (gradx(i,j) >= 0.) psi_val = psi_val - pi
105        psi_val = 3*pi/2. - psi_val
106        psi_val = psi_val*180./pi
107        psi_val = modulo(psi_val,360.)
108
109        theta_sl(ig0 + i) = theta_val
110        psi_sl(ig0 + i)   = psi_val
111    enddo
112enddo
113
114end subroutine getslopes
115
116!=======================================================================
117
118SUBROUTINE param_slope(csza,declin,rho,latitude,taudust,albedo,theta_s,psi_s,fdir_0,ftot_0,ftot)
119!***********************************************************************
120!
121! SUBROUTINE:
122! param_slope
123!
124!
125! PURPOSE:
126! computes total solar irradiance on a given Martian slope
127!
128!
129! INPUTS:
130! csza          cosine solar zenith angle
131! declin        sun declination (rad)
132! rho           sun right ascension (rad)
133! latitude      latitude (deg)
134! taudust       dust optical depth at reference wavelength 0.67 mic.
135! albedo        spectrally integrated surface Lambertian reflection albedo
136! theta_s       slope inclination angle (deg)
137!                       0 is horizontal, 90 is vertical
138! phi_s         slope azimuth (deg)
139!                       0 >> Northward
140!                       90 >> Eastward
141!                       180 >> Southward
142!                       270 >> Westward
143! ftot_0        spectrally integrated total irradiance on an horizontal surface (W/m2)
144! fdir_0        spectrally integrated direct irradiance on an horizontal surface (W/m2)
145!
146!
147! OUTPUTS:
148! ftot          spectrally integrated total irradiance on the slope (W/m2)
149!
150! REFERENCE:
151! "Fast and accurate estimation of irradiance on Martian slopes"
152! A. Spiga & F. Forget
153! .....
154!
155! AUTHOR:
156! A. Spiga (spiga@lmd.jussieu.fr)
157! March 2008
158!
159!***********************************************************************
160
161use comcstfi_h, only: pi
162
163implicit none
164
165! Input arguments
166real, intent(in) :: csza, declin, rho, latitude, taudust, theta_s, psi_s, albedo, ftot_0 , fdir_0
167
168! Output arguments
169real, intent(out) :: ftot
170
171! Local variables
172real                 :: deg2rad, a, mu_s, sigma_s, fdir, fscat, fscat_0, fref, ratio
173real, dimension(4,2) :: mat_M, mat_N, mat_T
174real, dimension(2)   :: g_vector
175real, dimension(4)   :: s_vector
176!***********************************************************************
177! Prerequisite
178deg2rad = pi/180.
179if ((theta_s > 90.) .or. (theta_s < 0.)) then
180    write(*,*) 'please set theta_s between 0 and 90', theta_s
[3203]181    call abort_physic("param_slopes","invalid theta_s",1)
[3183]182endif
183
184! Solar Zenith angle (radian)
185if (csza < 0.01) then
186    !print *, 'sun below horizon'
187    !fdir_0=0.
188    fdir    = 0.
189    fscat_0 = 0.
190    fscat   = 0.   
191    fref    = 0.
192else
193! Low incidence fix
194!    if (csza < 0.15) csza = 0.15
195
196! 'Slope vs Sun' azimuth (radian)
197    if (cos(declin)*sin(rho) == 0. .and. &
198        sin(deg2rad*latitude)*cos(declin)*cos(rho) - cos(deg2rad*latitude)*sin(declin) == 0.) then
199        a = deg2rad*psi_s ! some compilator need specfying value for atan2(0,0) 
200    else
201        a = deg2rad*psi_s + atan2(cos(declin)*sin(rho),sin(deg2rad*latitude)*cos(declin)*cos(rho)-cos(deg2rad*latitude)*sin(declin))
202    endif
203
204! Cosine of slope-sun phase angle
205    mu_s = csza*cos(deg2rad*theta_s) - cos(a)*sin(deg2rad*theta_s)*sqrt(1-csza**2)
206    if (mu_s <= 0.) mu_s = 0.
207
208! Sky-view factor
209    sigma_s=0.5*(1. + cos(deg2rad*theta_s))
210
211! Direct flux on the slope
212    fdir = fdir_0*mu_s/csza
213
214! Reflected flux on the slope
215    fref = albedo*(1 - sigma_s)*ftot_0
216
217! Scattered flux on a flat surface
218    fscat_0 = ftot_0 - fdir_0
219
220! Scattering vector (slope vs sky)
221    s_vector = (/ 1., exp(-taudust), sin(deg2rad*theta_s), sin(deg2rad*theta_s)*exp(-taudust) /)
222
223! Geometry vector (slope vs sun)
224    g_vector = (/ mu_s/csza, 1. /)
225
226! Coupling matrix
227    if (csza >= 0.5) then
228        mat_M(:,1) = (/ -0.264,  1.309,  0.208, -0.828 /)
229        mat_M(:,2) = (/  1.291*sigma_s, -1.371*sigma_s, -0.581,  1.641 /)
230        mat_N(:,1) = (/  0.911, -0.777, -0.223,  0.623 /)
231        mat_N(:,2) = (/ -0.933*sigma_s,  0.822*sigma_s,  0.514, -1.195 /)
232    else
233        mat_M(:,1) = (/ -0.373,  0.792, -0.095,  0.398 /)
234        mat_M(:,2) = (/  1.389*sigma_s, -0.794*sigma_s, -0.325,  0.183 /)
235        mat_N(:,1) = (/  1.079,  0.275,  0.419, -1.855 /)
236        mat_N(:,2) = (/ -1.076*sigma_s, -0.357*sigma_s, -0.075,  1.844 /)
237    endif
238    mat_T = mat_M + csza*mat_N
239
240! Scattered flux slope ratio
241    if (deg2rad*theta_s <= 0.0872664626) then
242    ! low angles
243        s_vector = (/ 1., exp(-taudust) , sin(0.0872664626), sin(0.0872664626)*exp(-taudust) /)
244        ratio = dot_product(matmul(s_vector, mat_T),g_vector)
245        ratio = 1. + (ratio - 1.)*deg2rad*theta_s/0.0872664626
246    else
247    ! general case
248    ratio = dot_product(matmul(s_vector,mat_T),g_vector) 
249    ! NB: ratio = dot_product(s_vector,matmul(mat_T,g_vector)) is equivalent
250    endif
251
252! Scattered flux on the slope
253    fscat = ratio*fscat_0
254endif ! if (csza < 0.01)
255
256! Total flux on the slope
257ftot = fdir + fref + fscat
258
259! Display results
260!  print *, 'sca component 0 ', ftot_0-fdir_0
261!  print *, 'dir component 0 ', fdir_0
262!  print *, 'scattered component ', fscat
263!  print *, 'direct component ', fdir
264!  print *, 'reflected component ', fref
265
266END SUBROUTINE param_slope
267
268!=======================================================================
269
270SUBROUTINE ini_slope_mod(ngrid)
271
272implicit none
273
274integer, intent(in) :: ngrid ! number of atmospheric columns
275
276allocate(theta_sl(ngrid))
277allocate(psi_sl(ngrid))
278
279END SUBROUTINE ini_slope_mod
280
281!=======================================================================
282
283SUBROUTINE end_slope_mod
284
285implicit none
286
287if (allocated(theta_sl)) deallocate(theta_sl)
288if (allocated(psi_sl)) deallocate(psi_sl)
289
290END SUBROUTINE end_slope_mod
291
292END MODULE slope_mod
Note: See TracBrowser for help on using the repository browser.