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

Last change on this file since 2909 was 2909, checked in by romain.vande, 2 years ago

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