source: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/subslope_mola.F90 @ 3094

Last change on this file since 3094 was 2910, checked in by emillour, 22 months ago

Mars PCM:
Make subslope_mola a module and turn (extremely) large static arrays into
allocated ones and assume that the "mola64.nc" file is in "datadir" rather
than a hard coded location.
EM
File size: 18.6 KB
Line 
1MODULE subslope_mola_mod
2
3IMPLICIT NONE
4
5CONTAINS
6
7SUBROUTINE subslope_mola(ngridmx,nslope,def_slope,subslope_dist)
8
9USE geometry_mod, ONLY: boundslon, boundslat ! boundaries of the cell (rad)
10use regular_lonlat_mod, ONLY : init_regular_lonlat, &
11                                 east, west, north, south, &
12                                 north_east, north_west, &
13                                 south_west, south_east
14USE datafile_mod, ONLY: datadir
15
16IMPLICIT NONE
17
18include "dimensions.h"
19
20      double precision :: resol
21      parameter(resol = 64)
22      integer :: jjm_mola, iim_mola
23      parameter(jjm_mola=180*resol, iim_mola=2*jjm_mola)
24      integer :: ierr
25! returned status code (==0 if OK)
26      real,allocatable :: theta_mola(:,:) !theta_mola(iim_mola,jjm_mola)
27!ED18 : slopes inclination (see getslope.F90)
28      real,allocatable :: psi_mola(:,:) !psi_mola(iim_mola,jjm_mola)
29!ED18 : slopes orientation (idem)     
30      REAL :: lon1, lon2, lat1, lat2 !bounds of the square
31      REAL :: lonlat_tmp !intermediate
32      INTEGER, INTENT(IN) :: nslope !nb of criteria intervals
33      REAL, INTENT(IN) :: def_slope(nslope+1) !list of values for criteria
34      REAL :: slopes_dist(nslope) !distribution of the slopes (with criterium criter) within the square
35      INTEGER :: k
36! Building of under-mesh statistics
37      INTEGER, INTENT(IN) :: ngridmx
38      INTEGER :: ig
39      real, INTENT(OUT) :: subslope_dist(ngridmx,nslope)
40      REAL, PARAMETER :: pi = acos(-1.)
41! diagnostics
42      REAL :: max_crit,max_lon,max_lat
43! special RFZ test
44      LOGICAL,PARAMETER :: rfz = .false.
45      INTEGER,PARAMETER :: igps = 1000 !first ig for RFZ location
46     
47      ! allocate big arrays
48      allocate(theta_mola(iim_mola,jjm_mola),stat=ierr)
49      if (ierr/=0) then
50        write(*,*)"subslope_mola: allocation of theta_mola(:,:) failed!"
51        stop
52      endif
53      allocate(psi_mola(iim_mola,jjm_mola),stat=ierr)
54      if (ierr/=0) then
55        write(*,*)"subslope_mola: allocation of psi_mola(:,:) failed!"
56        stop
57      endif
58     
59 
60!-------------Building of theta_mola and psi_mola
61      ! Assume that the mola file is to ben found in "datadir"
62      CALL mola(trim(datadir)//"/",&
63                ierr,theta_mola,psi_mola,resol,iim_mola,jjm_mola)
64
65!-------------Building of distribution
66      max_crit = 0.
67      DO ig = 1,ngridmx
68
69        lon1 = boundslon(ig,north_west) * 180./pi
70        lat1 = boundslat(ig,north_west) * 180./pi
71        lon2 = boundslon(ig,south_east) * 180./pi
72        lat2 = boundslat(ig,south_east) * 180./pi
73
74    !  IF(lon1.lt.-180.) lon1 = lon1+360.
75    !  IF(lon1.gt.180.) lon1 = lon1-360.
76    !  IF(lon2.gt.180.) lon2 = lon2-360.
77    !   IF(lon2.lt.-180.) lon2 = lon2+360.
78
79       IF (rfz.and.(ig.gt.igps)) THEN
80         slopes_dist(:) = 0.
81         slopes_dist(3) = 1.           
82       ELSE
83         CALL slopes_stat(lon1,lon2,lat1,lat2,slopes_dist, &
84             def_slope,nslope,theta_mola,psi_mola,iim_mola,jjm_mola, &
85             max_crit,max_lon,max_lat)
86       END IF
87         
88        DO k = 1, nslope
89          subslope_dist(ig,k) =slopes_dist(k)
90        ENDDO
91
92      ENDDO
93! correction try
94
95      if(nslope.eq.7) then
96       DO ig = 1,ngridmx
97        subslope_dist(ig,4) = 1 -  (subslope_dist(ig,1)+subslope_dist(ig,2) + &
98         subslope_dist(ig,3)+subslope_dist(ig,5)+subslope_dist(ig,6)+subslope_dist(ig,7))
99       ENDDO
100      endif
101   
102      PRINT*,'Diagnostics mola : max_crit, lon, lat = ', &
103              max_crit,max_lon,max_lat
104
105     ! deallocate big arrays (this routine is only called once)
106     deallocate(theta_mola,psi_mola)
107     
108END SUBROUTINE subslope_mola
109
110
111!==========================================================================================
112
113 SUBROUTINE mola(dset,ierr,theta_mola,psi_mola,resol,iim_mola,jjm_mola)
114
115
116!Give the MOLA altitude (alt), given lat and lon coordinates
117!Using bilinear interpolation from 32 pixels/degree MOLA file
118!
119!  12/2016 swiched some internal computations to double precision EM.
120!
121      implicit none
122
123      include "netcdf.inc"
124logical output_messages
125      parameter (output_messages=.true.)
126      double precision resol
127!c      parameter(resol=16) ! MOLA pixel/degree resolution
128
129
130      integer jjm_mola, iim_mola   ! # of longitude and latitude MOLA data values
131
132! Arguments
133! inputs
134      character*(*) dset ! Path to MCD datafiles
135
136! outputs
137      integer ierr ! returned status code (==0 if OK)
138      real theta_mola(iim_mola,jjm_mola) !ED18 : slopes inclination (see 1getslope.F90)
139      real psi_mola(iim_mola,jjm_mola) !ED18 : slopes orientation (idem)     
140!ED18 : rmq peut-être que ces tableaux seraient de taille iim_mola-1 *
141!jjm_mola-1
142!car les bords ne peuvent pas avoir de pente
143
144! Local variables
145
146      real latitude  ! north latitude (degrees)m
147      real longitude ! east longitude (degrees)
148      character*140 molafile ! MOLA datafile
149!      data molafile/'mola_32.2.nc'/
150!c      data molafile/'mola16.nc'/
151!c      data molafile/'mola32.nc'/
152!c      real invresol
153!c      parameter(invresol=1./resol)
154      integer*2  altmola(iim_mola,jjm_mola) ! MOLA altitude dataset
155!      save altmola !! ED18 or I get a bug (why?)
156
157      integer*2  mintopo_check,maxtopo_check ! known min and max of MOLAdataset
158!c      parameter(mintopo_check=-8156,maxtopo_check=21191) ! mola_32.2.nc
159!c      parameter(mintopo_check=-8177,maxtopo_check=21191) ! mola16.nc
160      parameter(mintopo_check=-8206,maxtopo_check=21191) ! mola32.nc
161      double precision dlat, dlon ! , lontmp
162      integer i,j ! ,count
163      double precision topo(4) ! neighboring MOLA points (for bilinear interpolation)
164      integer latsup,latinf,loninf,lonsup ! indexes of neighboring points
165      integer*2 mintopo, maxtopo ! min and max of read dataset
166      integer nid,nvarid ! NetCDF file and variable IDs
167      double precision colat ! colatitude
168      double precision lat,lon ! longitude and latitude, local values (in degrees)
169      real topogrid(iim_mola,jjm_mola) ! altmola in 'real' version
170
171
172!C 1.1. Open MOLA file
173         molafile=dset//'mola64.nc'
174         if (output_messages) then
175           write(*,*)"Loading MOLA topography from file ", &
176                       trim(molafile)
177         endif
178         ierr = NF_OPEN (molafile, NF_NOWRITE,nid)
179         if (ierr.NE.NF_NOERR) then
180            if (output_messages) then
181              write(*,*)"Error in mola: Could not open file ", &
182                       trim(molafile)
183            endif
184            stop
185         endif
186
187
188!C 1.2. Load data
189         ierr = NF_INQ_VARID (nid, "alt", nvarid)
190
191         ! note that MOLA "alt" are given as "short" (16 bits integers)
192         ierr = NF_GET_VAR(nid, nvarid, altmola)
193!                  ierr = NF_GET_VAR_INT2(nid, nvarid, altmola)
194         if (ierr.ne.NF_NOERR) then
195           if (output_messages) then
196            write(*,*)"Error in mola: <alt> not found"
197           endif
198           stop
199         endif
200
201
202
203
204!c 1.3. Close MOLA file
205         ierr=NF_CLOSE(nid)
206
207!c 1.4 Check that the MOLA dataset was correctly loaded
208
209!         mintopo=mintopo_check
210!         maxtopo=maxtopo_check
211!         do i=1,iim_mola
212!            do j=1,jjm_molaR
213!
214!               mintopo=min(mintopo,altmola(i,j))
215!               maxtopo=max(maxtopo,altmola(i,j))
216!            enddo
217!         enddo
218!         if ((mintopo.ne.mintopo_check).or. &
219!            (maxtopo.ne.maxtopo_check)) then
220!           if (output_messages) then
221!            write(*,*)"***ERROR Mola file ",molafile, &
222!                       " is not well read"
223!            write(*,*) "Minimum found: ", mintopo
224!            write(*,*) "Minimum should be:",mintopo_check
225!            write(*,*) "Maximum found: ", maxtopo
226!            write(*,*) "Maximum should be:",maxtopo_check
227!           endif
228!          ierr=16
229!          return
230!         endif
231         if (output_messages) then
232           write(*,*) "Done reading MOLA data"
233         endif
234
235!       topogrid = 1000*altmola
236           topogrid = float(altmola)
237
238      CALL getslopes_mola32(topogrid,theta_mola,psi_mola &
239     ,iim_mola,jjm_mola)
240
241
242END SUBROUTINE mola
243
244
245!-----------------------------------------------------------------------------------------
246
247SUBROUTINE  getslopes_mola32(topogrid,theta_mola,psi_mola,iim_mola,jjm_mola)
248
249
250implicit none
251
252
253! This routine computes slope inclination and orientation for the GCM
254! (callslope=.true. in callphys.def)
255! It works fine with a non-regular grid for zoomed simulations.
256! slope inclination angle (deg)  0 == horizontal, 90 == vertical
257! slope orientation angle (deg)  0 == Northward,  90 == Eastward, 180 ==
258! Southward, 270 == Westward
259! TN 04/1013
260integer :: res = 64
261integer, intent(in) :: iim_mola, jjm_mola
262real theta_mola(iim_mola,jjm_mola)
263real psi_mola(iim_mola,jjm_mola)
264real topogrid(iim_mola,jjm_mola) ! topography on lat/lon grid
265real latigrid(iim_mola,jjm_mola),longgrid(iim_mola,jjm_mola) ! meshgrid of latitude and longitude values (radians)
266real theta_val ! slope inclination
267real psi_val   ! slope orientation
268real gradx(iim_mola,jjm_mola) ! x: latitude-wise topography gradient, increasing northward
269real grady(iim_mola,jjm_mola) ! y: longitude-wise topography gradient, increasing westward (eastward?!)
270integer i,j,ig0
271!real :: theta_max
272real, parameter :: pi = acos(-1.)
273! WARNING: for Mars planet only
274real, parameter :: rad = 3396200.
275
276!theta_max = 0.
277
278do i=1, iim_mola
279  do j=1, jjm_mola
280    latigrid(i,j) = (90.01563 - j/float(res))*pi/180.
281    longgrid(i,j) = (-0.015625 + i/float(res))*pi/180.
282  enddo
283enddo
284
285! compute topography gradient
286! topogrid and rad are both in meters
287
288do j=1,jjm_mola
289
290!Special treatment for variable boundaries 
291   grady(1,j) = (topogrid(2,j) - topogrid(iim_mola,j)) / (2*pi+longgrid(2,j)-longgrid(iim_mola,j))
292   grady(1,j) = grady(1,j) / rad
293   grady(iim_mola,j) = (topogrid(1,j) - topogrid(iim_mola-1,j)) / &
294                          (2*pi+longgrid(1,j)-longgrid(iim_mola-1,j))
295   grady(iim_mola,j) = grady(iim_mola,j) / rad
296
297!Normal case
298   do i=2,iim_mola-1
299     grady(i,j) = (topogrid(i+1,j) - topogrid(i-1,j)) / (longgrid(i+1,j)-longgrid(i-1,j))
300     grady(i,j) = grady(i,j) / rad
301   enddo
302
303enddo
304
305
306do i=1, iim_mola
307
308!Special treatment for variable boundaries
309  if (i.ne.iim_mola/2) then
310       
311    gradx(i,1) = (topogrid(mod(i+iim_mola/2,iim_mola),1) - topogrid(i,2)) / &
312       (pi - latigrid(i,2)  - latigrid(mod(i+iim_mola/2,iim_mola),1))
313    gradx(i,1) = gradx(i,1) / rad
314    gradx(i,jjm_mola) = (topogrid(i,jjm_mola-1) - topogrid(mod(i+iim_mola/2,iim_mola),jjm_mola)) / &
315       (pi + latigrid(i,jjm_mola-1) + latigrid(mod(i+iim_mola/2,iim_mola),jjm_mola))
316    gradx(i,jjm_mola) = gradx(i,jjm_mola) / rad
317   else
318    gradx(i,1) =( topogrid(iim_mola,1) - topogrid(i,2) ) / &
319       (pi - latigrid(i,2)  - latigrid(iim_mola,1))
320    gradx(i,1) = gradx(i,1) / rad
321    gradx(i,jjm_mola) = (topogrid(i,jjm_mola-1) - topogrid(iim_mola,jjm_mola)) / &
322       (pi + latigrid(i,jjm_mola-1) + latigrid(iim_mola,jjm_mola))
323    gradx(i,jjm_mola) = gradx(i,jjm_mola) / rad
324   endif
325
326!Normal case
327   do j=2,jjm_mola-1
328               
329     gradx(i,j) = (topogrid(i,j+1) - topogrid(i,j-1)) / (latigrid(i,j+1)-latigrid(i,j-1))
330     gradx(i,j) = gradx(i,j) / rad
331
332   enddo
333
334enddo
335
336
337! compute slope inclination and orientation :
338do j=1,jjm_mola
339   do i=1,iim_mola
340     theta_val=atan(sqrt( (gradx(i,j))**2 + (grady(i,j))**2 ))*180./pi
341     psi_val=0.
342     if (gradx(i,j) .ne. 0.) psi_val= -pi/2. - atan(grady(i,j)/gradx(i,j))
343     if (gradx(i,j) .ge. 0.) psi_val= psi_val - pi
344     psi_val = 3*pi/2. - psi_val
345     psi_val = psi_val*180./pi
346     psi_val = MODULO(psi_val,360.)
347     theta_mola(i,j) = theta_val
348     psi_mola(i,j)   = psi_val
349!    if( (theta_val.lt.0).or.(psi_val.lt.0)) PRINT*,'valeurs abberrantes de &
350!      theta ou psi (= ',theta_val,psi_val, ')'
351!    if(theta_val.ge.theta_max) theta_max = theta_val
352   enddo
353enddo
354
355!print*,'theta max = ',theta_max
356
357
358END SUBROUTINE getslopes_mola32
359
360!-------------------------------------------------------------------------------------------
361
362
363SUBROUTINE slopes_stat(lon1_gcm,lon2_gcm,lat1,lat2,slopes_dist &
364          ,def_slope,nslope,theta_sl,psi_sl,nbp_lon,nbp_lat, &
365           max_crit,max_lon,max_lat)
366
367      IMPLICIT NONE
368
369      REAL,INTENT(IN) :: lat1, lat2 !bounds of the square
370      REAL,INTENT(IN) :: lon1_gcm, lon2_gcm ! bounds longitude in GCM ref  (-180E to 180E)
371      INTEGER,INTENT(IN) :: nslope !nb of criteria intervals
372      REAL,INTENT(IN) :: def_slope(nslope+1) !list of values for criteria
373      INTEGER,INTENT(IN) :: nbp_lon, nbp_lat
374      REAL,INTENT(IN) :: theta_sl(nbp_lon,nbp_lat)
375      REAL,INTENT(IN) :: psi_sl(nbp_lon,nbp_lat)
376      REAL slopes_dist(nslope) !distribution of the slopes (with criterium criter) within the square
377      INTEGER slopes_count(nslope) !intermediate for counting slopes
378      REAL, PARAMETER :: lat0 = 90.01563 !to be written differently with resol
379      REAL, PARAMETER :: lon0 = -0.015625
380      REAL :: lon1,lon2 ! bounds longitude in MOLA ref  (0 to 360E)
381
382      INTEGER :: i,j,k
383      INTEGER nb_slopes !number of slope points taken into account in the stat
384      INTEGER nb_slopes_tot  !number of slope points within the square           
385      INTEGER :: i_lon1, i_lon2
386      INTEGER :: j_lat1, j_lat2
387!      REAL :: crit !function
388      REAL :: val_crit !intermediate for evaluating the criterium
389      REAL :: max_crit !indicator for maximum slope criterium computed
390      REAL :: max_lon, max_lat
391          integer :: res = 64
392
393! Two cases must be distinguished
394
395   IF(lon1_gcm*lon2_gcm.gt.0) THEN !lon1_gcm & lon2_gcm on the same side
396      lon1 = lon1_gcm
397      IF (lon1_gcm.lt.0) lon1 = lon1_gcm + 360.
398      lon2 = lon2_gcm
399      IF (lon2_gcm.lt.0) lon2 = lon2_gcm + 360.
400
401      IF((lat1.lt.-90).or.(lat2.gt.90).or.(lon1.lt.0).or.(lon2.gt.360) &
402       .or.(lat1.lt.lat2).or.(lon1.gt.lon2)) THEN
403        PRINT*,'ERROR <subslope_mola>: latitudes or longitudes out of bounds'
404        PRINT*,'latitudes must be between -90 and 90 with lat1 > lat2'
405        PRINT*,'longitudes must be between 0 and 360 with lon1 < lon2'
406        PRINT*,'lat1, lat2 = ',lat1,' ',lat2
407        PRINT*,'lon1, lon2 = ',lon1,' ',lon2
408      STOP
409      ENDIF
410
411      i_lon1 = ceiling(res*(lon1-lon0))
412      i_lon2 = floor(res*(lon2-lon0))
413      j_lat1 = ceiling(res*(lat0-lat1))
414      j_lat2 = min(floor(res*(lat0-lat2)),nbp_lat)
415
416      nb_slopes = 0
417      nb_slopes_tot = 0
418      slopes_count(:) = 0
419
420!      PRINT*,'i_lon1 i_lon2 j_lat1 j_lat2'
421!      PRINT*,i_lon1,i_lon2,j_lat1,j_lat2
422   
423      IF(i_lon2.le.i_lon1) THEN
424        print*,"lon1, j_lon1 = ",lon1,i_lon1
425        print*,"lon2,j_lon2 = ",lon2,i_lon2
426      ENDIF 
427
428      DO i=i_lon1,i_lon2
429        DO j=j_lat1,j_lat2
430
431          DO k=1,nslope
432
433            val_crit = crit(theta_sl(i,j),psi_sl(i,j))
434            IF ((val_crit.ge.def_slope(k)).and. &
435                  (val_crit.lt.def_slope(k+1))) THEN
436             slopes_count(k) = slopes_count(k) + 1
437             nb_slopes = nb_slopes + 1
438            ENDIF
439            IF(abs(val_crit).gt.max_crit) THEN
440              max_crit = abs(val_crit)
441              max_lon = float(i)/float(res) + lon0
442              max_lat = lat0 - float(j)/float(res)
443            ENDIF
444          ENDDO
445          nb_slopes_tot = nb_slopes_tot + 1
446
447        ENDDO
448      ENDDO
449
450    ELSE !lon1_gcm<0 & lon2_gcm>0
451      !Special case, west and east side are treated separately
452      lon1 = lon1_gcm
453      IF (lon1_gcm.lt.0) lon1 = lon1_gcm + 360.
454      lon2 = lon2_gcm
455      IF (lon2_gcm.lt.0) lon2 = lon2_gcm + 360.
456
457      IF((lat1.lt.-90).or.(lat2.gt.90).or.(lon1.lt.0).or.(lon2.gt.360) &
458       .or.(lat1.lt.lat2).or.(lon1.lt.lon2)) THEN
459        PRINT*,'ERROR <subslope_mola>: latitudes or longitudes out of bounds'
460        PRINT*,'latitudes must be between -90 and 90 with lat1 > lat2'
461        PRINT*,'longitudes should be between 0 and 360 with &
462                  lon1 > lon2 (special case)'
463        PRINT*,'lat1, lat2 = ',lat1,' ',lat2
464        PRINT*,'lon1, lon2 = ',lon1,' ',lon2
465      STOP
466      ENDIF
467
468      i_lon1 = ceiling(res*(lon1-lon0))
469      i_lon2 = floor(res*(lon2-lon0))
470      j_lat1 = ceiling(res*(lat0-lat1))
471      j_lat2 = min(floor(res*(lat0-lat2)),nbp_lat)
472
473      nb_slopes = 0
474      nb_slopes_tot = 0
475      slopes_count(:) = 0
476
477!      PRINT*,'i_lon1 i_lon2 j_lat1 j_lat2'
478!      PRINT*,i_lon1,i_lon2,j_lat1,j_lat2
479
480!      IF(i_lon2.le.i_lon1) THEN
481!        print*,"lon1, j_lon1 = ",lon1,i_lon1
482!        print*,"lon2,j_lon2 = ",lon2,i_lon2
483!      ENDIF
484     
485      !computing west side
486      DO i=i_lon1,res*360
487        DO j=j_lat1,j_lat2
488            val_crit = crit(theta_sl(i,j),psi_sl(i,j))
489          DO k=1,nslope
490
491            IF ((val_crit.ge.def_slope(k)).and. &
492                  (val_crit.lt.def_slope(k+1))) THEN
493             slopes_count(k) = slopes_count(k) + 1
494             nb_slopes = nb_slopes + 1
495            ENDIF
496       
497
498            IF(abs(val_crit).gt.max_crit) THEN
499              max_crit = abs(val_crit)
500              max_lon = float(i)/float(res) + lon0
501              max_lat = lat0 - float(j)/float(res)
502            ENDIF
503          ENDDO
504          nb_slopes_tot = nb_slopes_tot + 1
505
506        ENDDO
507      ENDDO
508
509      !computing east side
510      DO i=1,i_lon2
511        DO j=j_lat1,j_lat2
512            val_crit = crit(theta_sl(i,j),psi_sl(i,j))
513          DO k=1,nslope
514
515            IF ((val_crit.ge.def_slope(k)).and. &
516                  (val_crit.lt.def_slope(k+1))) THEN
517             slopes_count(k) = slopes_count(k) + 1
518             nb_slopes = nb_slopes + 1
519            ENDIF
520 
521
522            IF(abs(val_crit).gt.max_crit) THEN
523              max_crit = abs(val_crit)
524              max_lon = float(i)/float(res) + lon0
525              max_lat = lat0 - float(j)/float(res)
526            ENDIF
527          ENDDO
528          nb_slopes_tot = nb_slopes_tot + 1
529
530        ENDDO
531      ENDDO
532
533     
534    ENDIF
535
536      IF(nb_slopes.ne.nb_slopes_tot) THEN
537        PRINT*,'WARNING: some slopes within the square are out of &
538       criteria'
539        PRINT*,'nb_slopes_tot = ',nb_slopes_tot
540        PRINT*,'nb_slopes = ',nb_slopes
541        PRINT*,float(nb_slopes_tot-nb_slopes)*100 / float(nb_slopes_tot) &
542        ,'% missing'
543      ENDIF
544
545!      PRINT*,"nb_slopes = ",nb_slopes
546! Here we truncate the distribution after 6decimals to have the sum of each  distribution = 1
547      DO k=1,nslope
548        slopes_dist(k)  = float(slopes_count(k)) /float(nb_slopes)
549        slopes_dist(k) = float(floor(1000000*slopes_dist(k)))/float(1000000)
550
551      ENDDO
552
553
554      CONTAINS
555
556!***********************************************************************
557!Defining the criterium used to compare slopes thanks to crit function
558!***********************************************************************
559
560      FUNCTION crit(theta,psi)
561
562      IMPLICIT NONE
563
564      REAL crit
565      REAL,INTENT(IN) :: theta, psi
566      real, parameter :: pi = acos(-1.)
567
568!      crit = theta
569! EX de critère différent
570      crit = theta*cos(psi*pi/180)
571
572      END FUNCTION crit
573
574
575END SUBROUTINE slopes_stat
576
577
578
579END MODULE subslope_mola_mod
Note: See TracBrowser for help on using the repository browser.