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: 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 : 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 : 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