MODULE MODULE_MAP
!------------------------------------------------------------------------------!
! MM5 Map projection parameters needed to compute MM5 coordinates (xj, yi)
! from latitude and longitude.
! get_map_params: Read the parameters from an MM5 V2/V3 input/output file
! set_map_param: Set the map parameters from the input read in MM5 file
!
! F. VANDENBERGHE, March 2001
!------------------------------------------------------------------------------!
! INTEGER :: IDD, IPROJ, IXC, JXC
! REAL :: PHIC, XLONC, XN, TRUELAT1, TRUELAT2, POLE, &
! DIS (10), XIM11 (10), XJM11 (10), &
! PSI1, C2, XCNTR, YCNTR
USE MODULE_MM5
!------------------------------------------------------------------------------!
CONTAINS
SUBROUTINE get_map_params (iunit)
!------------------------------------------------------------------------------!
! PURPOSE:
! -------
! READ THE MAP PROJECTION PARAMETERS: IPROJ, IXC, JXC, PHIC, XLONC, XN,
! TRUELAT1, TRUELAT2, POLE, DSC, DSM, XIM11, XJM11
! FROM A MM5 V2/V3 INPUT/OUTPUT FILE.
!
! Input:
! iunit: Logical file unit (file must be preably opened
!
! Output:
! IXC: COARSE DOMAIN GRID DIMENSION IN I (N-S) DIRECTION
! JXC: COARSE DOMAIN GRID DIMENSION IN J (E-W) DIRECTION
! IPROJ: MAP PROJECTION (1 = LAMBER, 2 = POLAR, 3 = MERCATOR)
! DSC: COARSE DOMAIN GRID DISTANCE IN KM
! PHIC: COARSE DOMAIN CENTER LATITUDE
! XLONC: COARSE DOMAIN CENTER LONGITUDE
! XN: CONE FACTOR
! TRUELAT1: TRUE LATITUDE 1 (DEGREE)
! TRUELAT2: TRUE LATITUDE 2 (DEGREE)
! POLE: POLE POSITION IN DEGREE LATIITUDE
! DSM : DOMAIN GRID DISTANCE IN KM
! XIM11: I LOCATION IN THE COARSE DOMAIN OF THE CURRENT DOMAIN POINT (1,1)
! XJM11: J LOCATION IN THE COARSE DOMAIN OF THE CURRENT DOMAIN POINT (1,1)
!
!------------------------------------------------------------------------------!
IMPLICIT NONE
INTEGER, INTENT (in) :: iunit
! INTEGER, INTENT (out) :: iproj, ixc, jxc
! REAL, INTENT (out) :: phic, pole, truelat1, xn, &
! xlonc, truelat2, dsc, dsm, xim11, xjm11
INTEGER :: io_error
CHARACTER (LEN= 8) :: message
INTEGER(kind=4) :: flag
integer(kind=4), dimension(50,20) :: bhi
real(kind=4), dimension(20,20) :: bhr
character(80), dimension(50,20) :: bhic
character(80), dimension(20,20) :: bhrc
!------------------------------------------------------------------------------!
REWIND (iunit)
! 1. READ V3 HEADER
! ==================
! 1.1 V3 very first record must be an integer
! --------------------------------------
flag = 2
READ (iunit, IOSTAT = io_error) flag
! 1.2 Then come headers
! -----------------
IF ((io_error .EQ. 0) .AND. (flag .EQ. 0)) THEN
READ (iunit, IOSTAT = io_error) bhi, bhr, bhic, bhrc
! 1.3 Integer map parameters
! ----------------------
IXC = bhi ( 5,1)
JXC = bhi ( 6,1)
IPROJ = bhi ( 7,1)
#ifdef BKG
NUMC = 1
NESTIX= bhi ( 5,1)
NESTJX= bhi ( 6,1)
NESTI = 1
NESTJ = 1
IDD = bhi(13,1)
NUMC (IDD) = bhi(14,1)
NESTIX(IDD) = bhi(16,1)
NESTJX(IDD) = bhi(17,1)
NESTI (IDD) = bhi(18,1)
NESTJ (IDD) = bhi(19,1)
! 1.4 Real map parameters
! -------------------
ptop = bhr(2,2)
ps0 = bhr(2,5)
ts0 = bhr(3,5)
tlp = bhr(4,5)
DIS (1) = bhr (1, 1) * 0.001
DIS (IDD) = bhr (9, 1) * 0.001
#else
IDD = bhi (13,1)
! 1.4 Real map parameters
! -------------------
DIS (1) = bhr (1, 1)
DIS (IDD) = bhr (9, 1)
#endif
PHIC = bhr (2, 1)
XLONC = bhr (3, 1)
XN = bhr (4, 1)
TRUELAT1 = bhr (5, 1)
TRUELAT2 = bhr (6, 1)
POLE = bhr ( 7, 1)
XIM11 (IDD) = bhr (10, 1)
XJM11 (IDD) = bhr (11, 1)
ELSE
write(unit=*, fmt='(a,i4)') &
'Error reading big header of unit:', iunit
call abort()
ENDIF
! 3. OTHER FILE FORMATS ARE UNKNOWN
! ==================================
IF (io_error .NE. 0) THEN
WRITE (message,'(I3)') iunit
CALL error_handler ('get_mm5_params.F90',&
' Cannot read file unit ', TRIM (message), .TRUE.)
ENDIF
! 4. PRINT OUT
! ============
WRITE (0,'(4(A,A,F7.2,/))',ADVANCE='no') &
"Map parameters: "," Phic = ",PHIC, &
" "," Xlonc= ",XLONC, &
" "," True1= ",TRUELAT1,&
" "," True2= ",TRUELAT2
WRITE (0,'(A,A,I4)',ADVANCE='no') &
" "," Iproj= ",IPROJ
SELECT CASE (IPROJ)
CASE (1); WRITE (0,'(A)') " (LAMBERT CONFORMAL)"
CASE (2); WRITE (0,'(A)') " (POLAR STEREO.)"
CASE (3); WRITE (0,'(A)') " (MERCATOR)"
CASE DEFAULT ; WRITE (0,'(A)') " (UNKNOWN)"
END SELECT
END SUBROUTINE GET_MAP_PARAMS
! SUBROUTINE SET_MAP_PARAMS (iproj, ixc, jxc, &
! phic, xlonc, xn, truelat1, truelat2, pole,&
! dsc, dsm, xim11, xjm11, &
! c2, psi1, xcntr, ycntr)
SUBROUTINE SET_MAP_PARAMS
!------------------------------------------------------------------------------!
! Purpose:
! -------
! Set the map projection parameters: c2, psi1, xcntr, ycntr
!
! Input:
! -----
! iproj: Map projection (1 = lamber, 2 = polar, 3 = mercator)
! phic: Coarse domain center latitude (in degree)
! pole: Pole position latitude (in degree)
! truelat1: True latitude for domain 1 (in degree)
!
! Output:
! ------
! c2:
! psi1:
! xcntr:
! ycntr:
!
!------------------------------------------------------------------------------!
IMPLICIT NONE
! INTEGER, INTENT (in) :: iproj, ixc, jxc
! REAL, INTENT (in) :: phic, pole, truelat1, xn, &
! xlonc, truelat2, dsc, dsm, xim11, xjm11
! REAL, INTENT (out) :: c2, psi1, xcntr, ycntr
REAL :: psx, cell, cell2, r, phictr
! Pi and earth radius in km
! REAL, PARAMETER :: PI = 3.1415926535897932346
! REAL, PARAMETER :: CONV = 180. / PI
! REAL, PARAMETER :: A = 6378.15
INCLUDE 'constants.inc'
!------------------------------------------------------------------------------C
! WRITE (0,'(A)') ' SET MAP PARAMETERS'
! DEFINE PSI1:
IF (IPROJ .EQ. 1 .OR. IPROJ.EQ.2) THEN
IF(PHIC .LT. 0) THEN
PSI1 = 90.+TRUELAT1
PSI1 = -PSI1
ELSE
PSI1 = 90.-TRUELAT1
ENDIF
ELSE
PSI1 = 0.
ENDIF
PSI1 = PSI1/CONV
! --------CALCULATE R
IF (IPROJ .NE. 3) THEN
PSX = (POLE-PHIC)/CONV
IF (IPROJ.EQ.1) THEN
CELL = A*SIN (PSI1)/XN
CELL2 = (TAN (PSX/2.))/(TAN (PSI1/2.))
WRITE (0,*) ' CELL = ' ,A,PSI1,XN,CELL
ENDIF
IF (IPROJ.EQ.2) THEN
CELL = A*SIN (PSX)/XN
CELL2 = (1. + COS (PSI1))/(1. + COS (PSX))
ENDIF
R = CELL*(CELL2)**XN
XCNTR = 0.0
YCNTR = -R
WRITE (0,*) ' YCNTR = ' ,CELL,CELL2,XN
ENDIF
!
! -----FOR MERCATOR PROJECTION, THE PROJECTION IS TRUE AT LAT AT PHI1
IF (IPROJ .EQ. 3) THEN
C2 = A*COS (PSI1)
XCNTR = 0.0
PHICTR = PHIC/CONV
CELL = COS (PHICTR)/(1.0+SIN (PHICTR))
YCNTR = - C2*ALOG (CELL)
ENDIF
WRITE (0,'(3(A,A,F7.2,/))',ADVANCE='no') &
" "," C2 = ",C2, &
" "," Psi1 = ",Psi1, &
" "," Xcntr= ",Xcntr
WRITE (0,'(1(A,A,F9.2))') &
" "," Ycntr= ",Ycntr
! 4. END
! =======
RETURN
END SUBROUTINE SET_MAP_PARAMS
SUBROUTINE LLXY (XLATI,XLONI,X,Y)
!------------------------------------------------------------------------------!
!
! ROUTINE LLXY
! **************
!
!
! PURPOSE: CALCULATES THE (X,Y) LOCATION (DOT) IN THE MESOSCALE GRIDS
! ------- FROM LATITUDES AND LONGITUDES
!
!
! INPUT:
! -----
! XLAT: LATITUDES
! XLON: LONGITUDES
!
! OUTPUT:
! -----
! X: THE COORDINATE IN X (J)-DIRECTION.
! Y: THE COORDINATE IN Y (I)-DIRECTION.
!
!------------------------------------------------------------------------------!
IMPLICIT NONE
REAL, INTENT(IN) :: XLATI, XLONI
! REAL, INTENT(IN) :: ID
REAL, INTENT(OUT) :: X, Y
REAL :: DXLON
REAL :: XLAT, XLON
REAL :: XX, YY, XC, YC
REAL :: CELL, PSI0, PSX, R, FLP
REAL :: CENTRI, CENTRJ
INTEGER :: NRATIO
real :: bb
INCLUDE 'constants.inc'
!------------------------------------------------------------------------------!
!
XLON = XLONI
XLAT = XLATI
!
XLAT = MAX (XLAT, -89.95)
XLAT = MIN (XLAT, +89.95)
!
DXLON = XLON - XLONC
IF (DXLON.GT. 180) DXLON = DXLON - 360.
IF (DXLON.LT.-180) DXLON = DXLON + 360.
IF (IPROJ.EQ.3) THEN
!
XC = XCNTR
YC = YCNTR
!
! For Mercator projection, high latitude has no definition:
IF (abs(XLAT) .gt.89.95) THEN
PRINT '(A,F8.3,A,F8.3)', &
"**WARNING** MERCATOR PROJECTION, BUT abs(XLAT)= |", &
XLAT,"| > 89.95 ?? XLON=",xlon
Y = -9.9999e8
X = -9.9999e8
RETURN
ENDIF
CELL = COS(XLAT/CONV)/(1.0+SIN(XLAT/CONV))
YY = -C2*ALOG(CELL)
XX = C2*DXLON/CONV
ELSE
!
PSI0 = ( POLE - PHIC )/CONV
XC = 0.0
!
!------ CALCULATE X,Y COORDS. RELATIVE TO POLE
!
FLP = XN*DXLON/CONV
PSX = ( POLE - XLAT )/CONV
if (IPROJ == 2) then
bb = 2.0*(COS(PSI1/2.0)**2)
YC = -A*bb*TAN(PSI0/2.0)
R = -A*bb*TAN(PSX/2.0)
else
bb = -A/XN*SIN(PSI1)
YC = bb*(TAN(PSI0/2.0)/TAN(PSI1/2.0))**XN
R = bb*(TAN(PSX /2.0)/TAN(PSI1/2.0))**XN
endif
IF (PHIC .LT. 0.0) THEN
XX = R*SIN(FLP)
YY = R*COS(FLP)
ELSE
XX = -R*SIN(FLP)
YY = R*COS(FLP)
END IF
ENDIF
!
! TRANSFORM (1,1) TO THE ORIGIN
!
CENTRI = FLOAT (IXC + 1)/2.0 ! the location of the center in the coarse
CENTRJ = FLOAT (JXC + 1)/2.0 ! domain
!
X = ( XX - XC )/DIS (1) + CENTRJ ! the (X,Y) coordinates in the coarse
Y = ( YY - YC )/DIS (1) + CENTRI ! domain
!
NRATIO = NINT(DIS (1) / DIS (IDD))
!
X = (X - XJM11 (IDD))*NRATIO + 1.0
Y = (Y - XIM11 (IDD))*NRATIO + 1.0
!
END SUBROUTINE LLXY
SUBROUTINE XYLL (XX,YY,XLAT,XLON)
!--------------------------------------------------------------------------CCC
!
! ROUTINE XYLL
! **************
!
!
! PURPOSE: CALCULATES THE LATITUDES AND LONGITUDES FROM THE
! ------- (X,Y) LOCATION (DOT) IN THE MESOSCALE GRIDS.
!
! INPUT:
! -----
! X: THE COORDINATE IN X (J)-DIRECTION.
! Y: THE COORDINATE IN Y (I)-DIRECTION.
!
! OUTPUT:
! -----
!
! XLAT: LATITUDES
! XLON: LONGITUDES
!----------------------------------------------------------------------------CC
!
IMPLICIT NONE
!
! ARGUMENT
!
REAL, INTENT(IN) :: XX, YY
REAL, INTENT(OUT) :: XLAT,XLON
! INTEGER, INTENT(IN) :: ID
!
! LOCAL
!
REAL :: CNTRI, CNTRJ
REAL :: X, Y
REAL :: FLP, FLPP, R, RXN, CELL, CEL1, CEL2, PSX
INCLUDE 'constants.inc'
!------------------------------------------------------------------------------C
CNTRI = float(IXC+1)/2.
CNTRJ = float(JXC+1)/2.
!
!-----CALCULATE X AND Y POSITIONS OF GRID AT DOT POINT
!
X = XCNTR+(XX-1.0)*DIS(IDD)+(XJM11(IDD)-CNTRJ)*DIS(1)
Y = YCNTR+(YY-1.0)*DIS(IDD)+(XIM11(IDD)-CNTRI)*DIS(1)
!-----NOW CALCULATE LAT AND LON OF THIS POINT
!
IF (IPROJ.NE.3) THEN
!
IF (Y.EQ.0.) THEN
!
IF (X.GE.0.0) FLP = 90.0/CONV
IF (X.LT.0.0) FLP = -90.0/CONV
!
ELSE
!
IF (PHIC.LT.0.0)THEN
FLP = ATAN2(X,Y)
ELSE
FLP = ATAN2(X,-Y)
ENDIF
!
ENDIF
!
FLPP = (FLP/XN)*CONV+XLONC
!
IF (FLPP.LT.-180.) FLPP = FLPP + 360
IF (FLPP.GT.180.) FLPP = FLPP - 360.
!
XLON = FLPP
!
!--------NOW SOLVE FOR LATITUDE
!
R = SQRT(X*X+Y*Y)
!
IF (PHIC.LT.0.0) R = -R
!
IF (IPROJ.EQ.1) THEN
CELL = (R*XN)/(A*SIN(PSI1))
RXN = 1.0/XN
CEL1 = TAN(PSI1/2.)*(CELL)**RXN
ENDIF
!
IF (IPROJ.EQ.2) THEN
CELL = R/A
CEL1 = CELL/(1.0+COS(PSI1))
ENDIF
!
CEL2 = ATAN(CEL1)
PSX = 2.*CEL2*CONV
XLAT = POLE-PSX
!
ENDIF
!
!-----CALCULATIONS FOR MERCATOR LAT,LON
!
IF (IPROJ.EQ.3) THEN
!
XLON = XLONC + ((X-XCNTR)/C2)*CONV
!
IF (XLON.LT.-180.) XLON = XLON + 360
IF (XLON.GT.180.) XLON = XLON - 360.
CELL = EXP(Y/C2)
XLAT = 2.*(CONV*ATAN(CELL))-90.0
!
ENDIF
!
END SUBROUTINE XYLL
!==============================================================================!
END MODULE MODULE_MAP