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

MODULE subslope_mola_mod

IMPLICIT NONE

CONTAINS

SUBROUTINE subslope_mola(ngridmx,nslope,def_slope,subslope_dist)

USE geometry_mod, ONLY: boundslon, boundslat ! boundaries of the cell (rad) 
use regular_lonlat_mod, ONLY : init_regular_lonlat, &
                                 east, west, north, south, &
                                 north_east, north_west, &
                                 south_west, south_east
USE datafile_mod, ONLY: datadir

IMPLICIT NONE

include "dimensions.h"

      double precision :: resol
      parameter(resol = 64)
      integer :: jjm_mola, iim_mola 
      parameter(jjm_mola=180*resol, iim_mola=2*jjm_mola)
      integer :: ierr 
! returned status code (==0 if OK)
      real,allocatable :: theta_mola(:,:) !theta_mola(iim_mola,jjm_mola) 
!ED18 : slopes inclination (see getslope.F90)
      real,allocatable :: psi_mola(:,:) !psi_mola(iim_mola,jjm_mola) 
!ED18 : slopes orientation (idem)     
      REAL :: lon1, lon2, lat1, lat2 !bounds of the square
      REAL :: lonlat_tmp !intermediate
      INTEGER, INTENT(IN) :: nslope !nb of criteria intervals
      REAL, INTENT(IN) :: def_slope(nslope+1) !list of values for criteria
      REAL :: slopes_dist(nslope) !distribution of the slopes (with criterium criter) within the square
      INTEGER :: k
! Building of under-mesh statistics
      INTEGER, INTENT(IN) :: ngridmx
      INTEGER :: ig
      real, INTENT(OUT) :: subslope_dist(ngridmx,nslope)
      REAL, PARAMETER :: pi = acos(-1.)
! diagnostics
      REAL :: max_crit,max_lon,max_lat
! special RFZ test
      LOGICAL,PARAMETER :: rfz = .false.
      INTEGER,PARAMETER :: igps = 1000 !first ig for RFZ location
      
      ! allocate big arrays
      allocate(theta_mola(iim_mola,jjm_mola),stat=ierr)
      if (ierr/=0) then
        write(*,*)"subslope_mola: allocation of theta_mola(:,:) failed!"
        stop
      endif
      allocate(psi_mola(iim_mola,jjm_mola),stat=ierr)
      if (ierr/=0) then
        write(*,*)"subslope_mola: allocation of psi_mola(:,:) failed!"
        stop
      endif
      
  
!-------------Building of theta_mola and psi_mola 
      ! Assume that the mola file is to ben found in "datadir"
      CALL mola(trim(datadir)//"/",&
                ierr,theta_mola,psi_mola,resol,iim_mola,jjm_mola)

!-------------Building of distribution
      max_crit = 0.
      DO ig = 1,ngridmx

        lon1 = boundslon(ig,north_west) * 180./pi 
        lat1 = boundslat(ig,north_west) * 180./pi
        lon2 = boundslon(ig,south_east) * 180./pi
        lat2 = boundslat(ig,south_east) * 180./pi 

    !  IF(lon1.lt.-180.) lon1 = lon1+360.
    !  IF(lon1.gt.180.) lon1 = lon1-360.
    !  IF(lon2.gt.180.) lon2 = lon2-360.
    !   IF(lon2.lt.-180.) lon2 = lon2+360.

       IF (rfz.and.(ig.gt.igps)) THEN
         slopes_dist(:) = 0.
         slopes_dist(3) = 1.            
       ELSE
         CALL slopes_stat(lon1,lon2,lat1,lat2,slopes_dist, &
             def_slope,nslope,theta_mola,psi_mola,iim_mola,jjm_mola, &
             max_crit,max_lon,max_lat)
       END IF
          
        DO k = 1, nslope
          subslope_dist(ig,k) =slopes_dist(k)
        ENDDO

      ENDDO
! correction try

      if(nslope.eq.7) then
       DO ig = 1,ngridmx
        subslope_dist(ig,4) = 1 -  (subslope_dist(ig,1)+subslope_dist(ig,2) + &
         subslope_dist(ig,3)+subslope_dist(ig,5)+subslope_dist(ig,6)+subslope_dist(ig,7))
       ENDDO
      endif
    
      PRINT*,'Diagnostics mola : max_crit, lon, lat = ', &
              max_crit,max_lon,max_lat

     ! deallocate big arrays (this routine is only called once)
     deallocate(theta_mola,psi_mola)
     
END SUBROUTINE subslope_mola


!==========================================================================================

 SUBROUTINE mola(dset,ierr,theta_mola,psi_mola,resol,iim_mola,jjm_mola)


!Give the MOLA altitude (alt), given lat and lon coordinates
!Using bilinear interpolation from 32 pixels/degree MOLA file
!
!  12/2016 swiched some internal computations to double precision EM.
!
      implicit none

      include "netcdf.inc"
logical output_messages
      parameter (output_messages=.true.)
      double precision resol
!c      parameter(resol=16) ! MOLA pixel/degree resolution


      integer jjm_mola, iim_mola   ! # of longitude and latitude MOLA data values 

! Arguments
! inputs
      character*(*) dset ! Path to MCD datafiles 

! outputs
      integer ierr ! returned status code (==0 if OK)
      real theta_mola(iim_mola,jjm_mola) !ED18 : slopes inclination (see 1getslope.F90)
      real psi_mola(iim_mola,jjm_mola) !ED18 : slopes orientation (idem)     
!ED18 : rmq peut-être que ces tableaux seraient de taille iim_mola-1 *
!jjm_mola-1 
!car les bords ne peuvent pas avoir de pente

! Local variables

      real latitude  ! north latitude (degrees)m
      real longitude ! east longitude (degrees)
      character*140 molafile ! MOLA datafile 
!      data molafile/'mola_32.2.nc'/
!c      data molafile/'mola16.nc'/
!c      data molafile/'mola32.nc'/
!c      real invresol
!c      parameter(invresol=1./resol)
      integer*2  altmola(iim_mola,jjm_mola) ! MOLA altitude dataset
!      save altmola !! ED18 or I get a bug (why?)

      integer*2  mintopo_check,maxtopo_check ! known min and max of MOLAdataset
!c      parameter(mintopo_check=-8156,maxtopo_check=21191) ! mola_32.2.nc
!c      parameter(mintopo_check=-8177,maxtopo_check=21191) ! mola16.nc
      parameter(mintopo_check=-8206,maxtopo_check=21191) ! mola32.nc
      double precision dlat, dlon ! , lontmp
      integer i,j ! ,count
      double precision topo(4) ! neighboring MOLA points (for bilinear interpolation)
      integer latsup,latinf,loninf,lonsup ! indexes of neighboring points
      integer*2 mintopo, maxtopo ! min and max of read dataset
      integer nid,nvarid ! NetCDF file and variable IDs
      double precision colat ! colatitude
      double precision lat,lon ! longitude and latitude, local values (in degrees)
      real topogrid(iim_mola,jjm_mola) ! altmola in 'real' version


!C 1.1. Open MOLA file
         molafile=dset//'mola64.nc'
         if (output_messages) then
           write(*,*)"Loading MOLA topography from file ", &
                       trim(molafile)
         endif
         ierr = NF_OPEN (molafile, NF_NOWRITE,nid)
         if (ierr.NE.NF_NOERR) then
            if (output_messages) then
              write(*,*)"Error in mola: Could not open file ", &
                       trim(molafile)
            endif
            stop
         endif


!C 1.2. Load data
         ierr = NF_INQ_VARID (nid, "alt", nvarid)

         ! note that MOLA "alt" are given as "short" (16 bits integers)
         ierr = NF_GET_VAR(nid, nvarid, altmola)
!                  ierr = NF_GET_VAR_INT2(nid, nvarid, altmola)
         if (ierr.ne.NF_NOERR) then
           if (output_messages) then
            write(*,*)"Error in mola: <alt> not found"
           endif
           stop
         endif




!c 1.3. Close MOLA file
         ierr=NF_CLOSE(nid)

!c 1.4 Check that the MOLA dataset was correctly loaded

!         mintopo=mintopo_check
!         maxtopo=maxtopo_check
!         do i=1,iim_mola
!            do j=1,jjm_molaR
!
!               mintopo=min(mintopo,altmola(i,j))
!               maxtopo=max(maxtopo,altmola(i,j))
!            enddo
!         enddo
!         if ((mintopo.ne.mintopo_check).or. &
!            (maxtopo.ne.maxtopo_check)) then
!           if (output_messages) then
!            write(*,*)"***ERROR Mola file ",molafile, &
!                       " is not well read"
!            write(*,*) "Minimum found: ", mintopo
!            write(*,*) "Minimum should be:",mintopo_check
!            write(*,*) "Maximum found: ", maxtopo
!            write(*,*) "Maximum should be:",maxtopo_check
!           endif
!          ierr=16
!          return
!         endif
         if (output_messages) then
           write(*,*) "Done reading MOLA data"
         endif

!       topogrid = 1000*altmola
           topogrid = float(altmola)

      CALL getslopes_mola32(topogrid,theta_mola,psi_mola &
     ,iim_mola,jjm_mola)


END SUBROUTINE mola


!-----------------------------------------------------------------------------------------

SUBROUTINE  getslopes_mola32(topogrid,theta_mola,psi_mola,iim_mola,jjm_mola)


implicit none


! This routine computes slope inclination and orientation for the GCM
! (callslope=.true. in callphys.def)
! It works fine with a non-regular grid for zoomed simulations.
! slope inclination angle (deg)  0 == horizontal, 90 == vertical
! slope orientation angle (deg)  0 == Northward,  90 == Eastward, 180 ==
! Southward, 270 == Westward
! TN 04/1013
integer :: res = 64
integer, intent(in) :: iim_mola, jjm_mola
real theta_mola(iim_mola,jjm_mola)
real psi_mola(iim_mola,jjm_mola)
real topogrid(iim_mola,jjm_mola) ! topography on lat/lon grid 
real latigrid(iim_mola,jjm_mola),longgrid(iim_mola,jjm_mola) ! meshgrid of latitude and longitude values (radians)
real theta_val ! slope inclination
real psi_val   ! slope orientation
real gradx(iim_mola,jjm_mola) ! x: latitude-wise topography gradient, increasing northward
real grady(iim_mola,jjm_mola) ! y: longitude-wise topography gradient, increasing westward (eastward?!)
integer i,j,ig0
!real :: theta_max
real, parameter :: pi = acos(-1.)
! WARNING: for Mars planet only
real, parameter :: rad = 3396200.

!theta_max = 0.

do i=1, iim_mola
  do j=1, jjm_mola
    latigrid(i,j) = (90.01563 - j/float(res))*pi/180.
    longgrid(i,j) = (-0.015625 + i/float(res))*pi/180.
  enddo
enddo

! compute topography gradient
! topogrid and rad are both in meters

do j=1,jjm_mola

!Special treatment for variable boundaries  
   grady(1,j) = (topogrid(2,j) - topogrid(iim_mola,j)) / (2*pi+longgrid(2,j)-longgrid(iim_mola,j))
   grady(1,j) = grady(1,j) / rad
   grady(iim_mola,j) = (topogrid(1,j) - topogrid(iim_mola-1,j)) / &
                          (2*pi+longgrid(1,j)-longgrid(iim_mola-1,j))
   grady(iim_mola,j) = grady(iim_mola,j) / rad

!Normal case
   do i=2,iim_mola-1
     grady(i,j) = (topogrid(i+1,j) - topogrid(i-1,j)) / (longgrid(i+1,j)-longgrid(i-1,j))
     grady(i,j) = grady(i,j) / rad
   enddo

enddo


do i=1, iim_mola

!Special treatment for variable boundaries 
  if (i.ne.iim_mola/2) then
        
    gradx(i,1) = (topogrid(mod(i+iim_mola/2,iim_mola),1) - topogrid(i,2)) / &
       (pi - latigrid(i,2)  - latigrid(mod(i+iim_mola/2,iim_mola),1))
    gradx(i,1) = gradx(i,1) / rad
    gradx(i,jjm_mola) = (topogrid(i,jjm_mola-1) - topogrid(mod(i+iim_mola/2,iim_mola),jjm_mola)) / &
       (pi + latigrid(i,jjm_mola-1) + latigrid(mod(i+iim_mola/2,iim_mola),jjm_mola))
    gradx(i,jjm_mola) = gradx(i,jjm_mola) / rad
   else
    gradx(i,1) =( topogrid(iim_mola,1) - topogrid(i,2) ) / &
       (pi - latigrid(i,2)  - latigrid(iim_mola,1))
    gradx(i,1) = gradx(i,1) / rad
    gradx(i,jjm_mola) = (topogrid(i,jjm_mola-1) - topogrid(iim_mola,jjm_mola)) / &
       (pi + latigrid(i,jjm_mola-1) + latigrid(iim_mola,jjm_mola))
    gradx(i,jjm_mola) = gradx(i,jjm_mola) / rad
   endif

!Normal case
   do j=2,jjm_mola-1
                
     gradx(i,j) = (topogrid(i,j+1) - topogrid(i,j-1)) / (latigrid(i,j+1)-latigrid(i,j-1))
     gradx(i,j) = gradx(i,j) / rad

   enddo

enddo


! compute slope inclination and orientation :
do j=1,jjm_mola
   do i=1,iim_mola
     theta_val=atan(sqrt( (gradx(i,j))**2 + (grady(i,j))**2 ))*180./pi
     psi_val=0.
     if (gradx(i,j) .ne. 0.) psi_val= -pi/2. - atan(grady(i,j)/gradx(i,j))
     if (gradx(i,j) .ge. 0.) psi_val= psi_val - pi
     psi_val = 3*pi/2. - psi_val
     psi_val = psi_val*180./pi
     psi_val = MODULO(psi_val,360.)
     theta_mola(i,j) = theta_val
     psi_mola(i,j)   = psi_val
!    if( (theta_val.lt.0).or.(psi_val.lt.0)) PRINT*,'valeurs abberrantes de &
!      theta ou psi (= ',theta_val,psi_val, ')'
!    if(theta_val.ge.theta_max) theta_max = theta_val
   enddo
enddo

!print*,'theta max = ',theta_max


END SUBROUTINE getslopes_mola32

!-------------------------------------------------------------------------------------------


SUBROUTINE slopes_stat(lon1_gcm,lon2_gcm,lat1,lat2,slopes_dist &
          ,def_slope,nslope,theta_sl,psi_sl,nbp_lon,nbp_lat, &
           max_crit,max_lon,max_lat)

      IMPLICIT NONE

      REAL,INTENT(IN) :: lat1, lat2 !bounds of the square 
      REAL,INTENT(IN) :: lon1_gcm, lon2_gcm ! bounds longitude in GCM ref  (-180E to 180E)
      INTEGER,INTENT(IN) :: nslope !nb of criteria intervals
      REAL,INTENT(IN) :: def_slope(nslope+1) !list of values for criteria
      INTEGER,INTENT(IN) :: nbp_lon, nbp_lat
      REAL,INTENT(IN) :: theta_sl(nbp_lon,nbp_lat)
      REAL,INTENT(IN) :: psi_sl(nbp_lon,nbp_lat)
      REAL slopes_dist(nslope) !distribution of the slopes (with criterium criter) within the square 
      INTEGER slopes_count(nslope) !intermediate for counting slopes
      REAL, PARAMETER :: lat0 = 90.01563 !to be written differently with resol
      REAL, PARAMETER :: lon0 = -0.015625
      REAL :: lon1,lon2 ! bounds longitude in MOLA ref  (0 to 360E)

      INTEGER :: i,j,k
      INTEGER nb_slopes !number of slope points taken into account in the stat
      INTEGER nb_slopes_tot  !number of slope points within the square           
      INTEGER :: i_lon1, i_lon2
      INTEGER :: j_lat1, j_lat2
!      REAL :: crit !function
      REAL :: val_crit !intermediate for evaluating the criterium
      REAL :: max_crit !indicator for maximum slope criterium computed
      REAL :: max_lon, max_lat
          integer :: res = 64

! Two cases must be distinguished

   IF(lon1_gcm*lon2_gcm.gt.0) THEN !lon1_gcm & lon2_gcm on the same side 
      lon1 = lon1_gcm
      IF (lon1_gcm.lt.0) lon1 = lon1_gcm + 360.
      lon2 = lon2_gcm
      IF (lon2_gcm.lt.0) lon2 = lon2_gcm + 360.

      IF((lat1.lt.-90).or.(lat2.gt.90).or.(lon1.lt.0).or.(lon2.gt.360) &
       .or.(lat1.lt.lat2).or.(lon1.gt.lon2)) THEN
        PRINT*,'ERROR <subslope_mola>: latitudes or longitudes out of bounds'
        PRINT*,'latitudes must be between -90 and 90 with lat1 > lat2'
        PRINT*,'longitudes must be between 0 and 360 with lon1 < lon2'
        PRINT*,'lat1, lat2 = ',lat1,' ',lat2
        PRINT*,'lon1, lon2 = ',lon1,' ',lon2
      STOP
      ENDIF

      i_lon1 = ceiling(res*(lon1-lon0))
      i_lon2 = floor(res*(lon2-lon0))
      j_lat1 = ceiling(res*(lat0-lat1))
      j_lat2 = min(floor(res*(lat0-lat2)),nbp_lat)

      nb_slopes = 0
      nb_slopes_tot = 0
      slopes_count(:) = 0

!      PRINT*,'i_lon1 i_lon2 j_lat1 j_lat2'
!      PRINT*,i_lon1,i_lon2,j_lat1,j_lat2
    
      IF(i_lon2.le.i_lon1) THEN
        print*,"lon1, j_lon1 = ",lon1,i_lon1
        print*,"lon2,j_lon2 = ",lon2,i_lon2
      ENDIF  

      DO i=i_lon1,i_lon2
        DO j=j_lat1,j_lat2

          DO k=1,nslope

            val_crit = crit(theta_sl(i,j),psi_sl(i,j))
            IF ((val_crit.ge.def_slope(k)).and. &
                  (val_crit.lt.def_slope(k+1))) THEN
             slopes_count(k) = slopes_count(k) + 1
             nb_slopes = nb_slopes + 1
            ENDIF
            IF(abs(val_crit).gt.max_crit) THEN
              max_crit = abs(val_crit)
              max_lon = float(i)/float(res) + lon0
              max_lat = lat0 - float(j)/float(res)
            ENDIF
          ENDDO
          nb_slopes_tot = nb_slopes_tot + 1

        ENDDO
      ENDDO

    ELSE !lon1_gcm<0 & lon2_gcm>0
      !Special case, west and east side are treated separately
      lon1 = lon1_gcm
      IF (lon1_gcm.lt.0) lon1 = lon1_gcm + 360.
      lon2 = lon2_gcm
      IF (lon2_gcm.lt.0) lon2 = lon2_gcm + 360.

      IF((lat1.lt.-90).or.(lat2.gt.90).or.(lon1.lt.0).or.(lon2.gt.360) &
       .or.(lat1.lt.lat2).or.(lon1.lt.lon2)) THEN
        PRINT*,'ERROR <subslope_mola>: latitudes or longitudes out of bounds'
        PRINT*,'latitudes must be between -90 and 90 with lat1 > lat2'
        PRINT*,'longitudes should be between 0 and 360 with &
                  lon1 > lon2 (special case)'
        PRINT*,'lat1, lat2 = ',lat1,' ',lat2
        PRINT*,'lon1, lon2 = ',lon1,' ',lon2
      STOP
      ENDIF

      i_lon1 = ceiling(res*(lon1-lon0))
      i_lon2 = floor(res*(lon2-lon0))
      j_lat1 = ceiling(res*(lat0-lat1))
      j_lat2 = min(floor(res*(lat0-lat2)),nbp_lat)

      nb_slopes = 0
      nb_slopes_tot = 0
      slopes_count(:) = 0

!      PRINT*,'i_lon1 i_lon2 j_lat1 j_lat2'
!      PRINT*,i_lon1,i_lon2,j_lat1,j_lat2

!      IF(i_lon2.le.i_lon1) THEN
!        print*,"lon1, j_lon1 = ",lon1,i_lon1
!        print*,"lon2,j_lon2 = ",lon2,i_lon2
!      ENDIF
      
      !computing west side
      DO i=i_lon1,res*360
        DO j=j_lat1,j_lat2
            val_crit = crit(theta_sl(i,j),psi_sl(i,j))
          DO k=1,nslope

            IF ((val_crit.ge.def_slope(k)).and. &
                  (val_crit.lt.def_slope(k+1))) THEN
             slopes_count(k) = slopes_count(k) + 1
             nb_slopes = nb_slopes + 1
            ENDIF
       

            IF(abs(val_crit).gt.max_crit) THEN
              max_crit = abs(val_crit)
              max_lon = float(i)/float(res) + lon0
              max_lat = lat0 - float(j)/float(res)
            ENDIF
          ENDDO
          nb_slopes_tot = nb_slopes_tot + 1

        ENDDO
      ENDDO

      !computing east side
      DO i=1,i_lon2
        DO j=j_lat1,j_lat2
            val_crit = crit(theta_sl(i,j),psi_sl(i,j))
          DO k=1,nslope

            IF ((val_crit.ge.def_slope(k)).and. &
                  (val_crit.lt.def_slope(k+1))) THEN
             slopes_count(k) = slopes_count(k) + 1
             nb_slopes = nb_slopes + 1
            ENDIF
 

            IF(abs(val_crit).gt.max_crit) THEN
              max_crit = abs(val_crit)
              max_lon = float(i)/float(res) + lon0
              max_lat = lat0 - float(j)/float(res)
            ENDIF
          ENDDO
          nb_slopes_tot = nb_slopes_tot + 1

        ENDDO
      ENDDO

      
    ENDIF

      IF(nb_slopes.ne.nb_slopes_tot) THEN
        PRINT*,'WARNING: some slopes within the square are out of &
       criteria'
        PRINT*,'nb_slopes_tot = ',nb_slopes_tot
        PRINT*,'nb_slopes = ',nb_slopes
        PRINT*,float(nb_slopes_tot-nb_slopes)*100 / float(nb_slopes_tot) &
        ,'% missing'
      ENDIF

!      PRINT*,"nb_slopes = ",nb_slopes
! Here we truncate the distribution after 6decimals to have the sum of each  distribution = 1
      DO k=1,nslope
        slopes_dist(k)  = float(slopes_count(k)) /float(nb_slopes)
        slopes_dist(k) = float(floor(1000000*slopes_dist(k)))/float(1000000)

      ENDDO


      CONTAINS

!***********************************************************************
!Defining the criterium used to compare slopes thanks to crit function
!***********************************************************************

      FUNCTION crit(theta,psi)

      IMPLICIT NONE

      REAL crit
      REAL,INTENT(IN) :: theta, psi
      real, parameter :: pi = acos(-1.)

!      crit = theta
! EX de critère différent
      crit = theta*cos(psi*pi/180)

      END FUNCTION crit


END SUBROUTINE slopes_stat



END MODULE subslope_mola_mod