var/obsproc/src/module_map.F90

   1 MODULE MODULE_MAP
   2
   3 !------------------------------------------------------------------------------!
   4 ! MM5 Map projection parameters needed to compute MM5 coordinates (xj, yi)
   5 ! from latitude and longitude.
   6 ! get_map_params: Read the parameters from an MM5 V2/V3 input/output file
   7 ! set_map_param:  Set the map parameters from the input read in MM5 file
   8 !
   9 ! F. VANDENBERGHE, March 2001
  10 !------------------------------------------------------------------------------!
  11
  12
  13 ! INTEGER  :: IDD, IPROJ, IXC,   JXC
  14 ! REAL     :: PHIC,  XLONC, XN,    TRUELAT1, TRUELAT2, POLE, &
  15 !             DIS (10),   XIM11 (10), XJM11 (10), &
  16 !             PSI1,  C2,    XCNTR, YCNTR
  17
  18   USE MODULE_MM5
  19
  20 !------------------------------------------------------------------------------!
  21 CONTAINS
  22
  23       SUBROUTINE get_map_params (iunit)
  24 !------------------------------------------------------------------------------!
  25 ! PURPOSE:
  26 ! -------
  27 !   READ THE MAP PROJECTION PARAMETERS: IPROJ, IXC, JXC, PHIC, XLONC, XN,
  28 !   TRUELAT1, TRUELAT2, POLE, DSC, DSM, XIM11, XJM11
  29 !   FROM A MM5 V2/V3 INPUT/OUTPUT FILE.
  30 !
  31 ! Input:
  32 !   iunit:  Logical file unit (file must be preably opened
  33 !
  34 ! Output:
  35 !   IXC:      COARSE DOMAIN GRID DIMENSION IN I (N-S) DIRECTION
  36 !   JXC:      COARSE DOMAIN GRID DIMENSION IN J (E-W) DIRECTION
  37 !   IPROJ:    MAP PROJECTION (1 = LAMBER, 2 = POLAR, 3 = MERCATOR)
  38 !   DSC:      COARSE DOMAIN GRID DISTANCE IN KM
  39 !   PHIC:     COARSE DOMAIN CENTER LATITUDE
  40 !   XLONC:    COARSE DOMAIN CENTER LONGITUDE
  41 !   XN:       CONE FACTOR
  42 !   TRUELAT1: TRUE LATITUDE 1 (DEGREE)
  43 !   TRUELAT2: TRUE LATITUDE 2 (DEGREE)
  44 !   POLE:     POLE POSITION IN DEGREE LATIITUDE
  45 !   DSM :     DOMAIN GRID DISTANCE IN KM
  46 !   XIM11:    I LOCATION IN THE COARSE DOMAIN OF THE CURRENT DOMAIN POINT (1,1)
  47 !   XJM11:    J LOCATION IN THE COARSE DOMAIN OF THE CURRENT DOMAIN POINT (1,1)
  48 !
  49 !------------------------------------------------------------------------------!
  50       IMPLICIT NONE
  51
  52       INTEGER, INTENT (in)  :: iunit
  53 !     INTEGER, INTENT (out) :: iproj, ixc, jxc
  54 !     REAL,    INTENT (out) :: phic, pole, truelat1, xn, &
  55 !                              xlonc, truelat2, dsc, dsm, xim11, xjm11
  56
  57       INTEGER               :: io_error
  58       CHARACTER (LEN= 8)    :: message
  59
  60       INTEGER(kind=4)       :: flag
  61       integer(kind=4), dimension(50,20) :: bhi
  62       real(kind=4),    dimension(20,20) :: bhr
  63       character(80),   dimension(50,20) :: bhic
  64       character(80),   dimension(20,20) :: bhrc
  65
  66 !------------------------------------------------------------------------------!
  67
  68       REWIND (iunit)
  69
  70 ! 1.  READ V3 HEADER
  71 ! ==================
  72
  73 ! 1.1 V3 very first record must be an integer
  74 !     --------------------------------------
  75
  76       flag = 2
  77
  78       READ (iunit, IOSTAT = io_error) flag
  79
  80 ! 1.2 Then come headers
  81 !     -----------------
  82
  83       IF ((io_error .EQ. 0) .AND. (flag .EQ. 0)) THEN
  84
  85            READ (iunit, IOSTAT = io_error) bhi, bhr, bhic, bhrc
  86
  87 ! 1.3 Integer map parameters
  88 !     ----------------------
  89
  90            IXC   = bhi ( 5,1)
  91            JXC   = bhi ( 6,1)
  92            IPROJ = bhi ( 7,1)
  93 #ifdef BKG
  94
  95            NUMC  = 1
  96
  97            NESTIX= bhi ( 5,1)
  98            NESTJX= bhi ( 6,1)
  99
 100            NESTI = 1
 101            NESTJ = 1
 102
 103            IDD         = bhi(13,1)
 104            NUMC  (IDD) = bhi(14,1)
 105
 106            NESTIX(IDD) = bhi(16,1)
 107            NESTJX(IDD) = bhi(17,1)
 108
 109            NESTI (IDD) = bhi(18,1)
 110            NESTJ (IDD) = bhi(19,1)
 111
 112 ! 1.4 Real map parameters
 113 !     -------------------
 114
 115            ptop  = bhr(2,2)
 116            ps0   = bhr(2,5)
 117            ts0   = bhr(3,5)
 118            tlp   = bhr(4,5)
 119
 120            DIS (1)     = bhr (1, 1) * 0.001
 121            DIS (IDD)   = bhr (9, 1) * 0.001
 122 #else
 123            IDD   = bhi (13,1)
 124
 125 ! 1.4 Real map parameters
 126 !     -------------------
 127
 128            DIS (1)     = bhr (1, 1)
 129            DIS (IDD)   = bhr (9, 1)
 130 #endif
 131            PHIC        = bhr (2, 1)
 132            XLONC       = bhr (3, 1)
 133            XN          = bhr (4, 1)
 134            TRUELAT1    = bhr (5, 1)
 135            TRUELAT2    = bhr (6, 1)
 136            POLE        = bhr ( 7, 1)
 137            XIM11 (IDD) = bhr (10, 1)
 138            XJM11 (IDD) = bhr (11, 1)
 139
 140
 141       ELSE
 142
 143          write(unit=*, fmt='(a,i4)') &
 144               'Error reading big header of unit:', iunit
 145          call abort()
 146
 147       ENDIF
 148
 149
 150 ! 3.  OTHER FILE FORMATS ARE UNKNOWN
 151 ! ==================================
 152
 153       IF (io_error .NE. 0) THEN
 154           WRITE (message,'(I3)') iunit
 155           CALL error_handler ('get_mm5_params.F90',&
 156              ' Cannot read file unit ', TRIM  (message), .TRUE.)
 157       ENDIF
 158
 159
 160 ! 4. PRINT OUT
 161 ! ============
 162
 163       WRITE (0,'(4(A,A,F7.2,/))',ADVANCE='no')  &
 164       "Map parameters: "," Phic = ",PHIC,    &
 165       "                "," Xlonc= ",XLONC,   &
 166       "                "," True1= ",TRUELAT1,&
 167       "                "," True2= ",TRUELAT2
 168
 169       WRITE (0,'(A,A,I4)',ADVANCE='no') &
 170       "                "," Iproj= ",IPROJ
 171
 172       SELECT CASE (IPROJ)
 173              CASE (1);      WRITE (0,'(A)') " (LAMBERT CONFORMAL)"
 174              CASE (2);      WRITE (0,'(A)') " (POLAR STEREO.)"
 175              CASE (3);      WRITE (0,'(A)') " (MERCATOR)"
 176              CASE DEFAULT ; WRITE (0,'(A)') " (UNKNOWN)"
 177       END SELECT
 178
 179       END SUBROUTINE GET_MAP_PARAMS
 180
 181 !     SUBROUTINE SET_MAP_PARAMS (iproj, ixc,   jxc, &
 182 !                                phic,  xlonc, xn,    truelat1, truelat2, pole,&
 183 !                                dsc,   dsm,   xim11, xjm11, &
 184 !                                c2, psi1, xcntr, ycntr)
 185       SUBROUTINE SET_MAP_PARAMS
 186 !------------------------------------------------------------------------------!
 187 ! Purpose:
 188 ! -------
 189 !   Set the map projection parameters: c2, psi1, xcntr, ycntr
 190 !
 191 ! Input:
 192 ! -----
 193 !  iproj:    Map projection (1 = lamber, 2 = polar, 3 = mercator)
 194 !  phic:     Coarse domain center latitude (in degree)
 195 !  pole:     Pole position latitude (in degree)
 196 !  truelat1: True latitude for domain 1 (in degree)
 197 !
 198 ! Output:
 199 ! ------
 200 !  c2:
 201 !  psi1:
 202 !  xcntr:
 203 !  ycntr:
 204 !
 205 !------------------------------------------------------------------------------!
 206       IMPLICIT NONE
 207
 208 !     INTEGER, INTENT (in)  :: iproj, ixc, jxc
 209 !     REAL,    INTENT (in)  :: phic, pole, truelat1, xn, &
 210 !                              xlonc, truelat2, dsc, dsm, xim11, xjm11
 211 !     REAL,    INTENT (out) :: c2, psi1, xcntr, ycntr
 212
 213       REAL :: psx, cell, cell2, r, phictr
 214
 215       ! Pi and earth radius in km
 216
 217 !     REAL, PARAMETER :: PI   = 3.1415926535897932346
 218 !     REAL, PARAMETER :: CONV = 180. / PI
 219 !     REAL, PARAMETER :: A    = 6378.15
 220
 221       INCLUDE 'constants.inc'
 222 !------------------------------------------------------------------------------C
 223
 224 !     WRITE (0,'(A)') ' SET MAP PARAMETERS'
 225
 226 !   DEFINE PSI1:
 227
 228         IF (IPROJ .EQ. 1 .OR. IPROJ.EQ.2) THEN
 229
 230             IF(PHIC .LT. 0) THEN
 231                PSI1 = 90.+TRUELAT1
 232                PSI1 = -PSI1
 233             ELSE
 234                PSI1 = 90.-TRUELAT1
 235             ENDIF
 236
 237         ELSE
 238
 239             PSI1 = 0.
 240
 241         ENDIF
 242
 243         PSI1 = PSI1/CONV
 244
 245 ! --------CALCULATE R
 246
 247         IF (IPROJ .NE. 3) THEN
 248
 249             PSX = (POLE-PHIC)/CONV
 250
 251             IF (IPROJ.EQ.1) THEN
 252
 253                 CELL  = A*SIN (PSI1)/XN
 254                 CELL2 = (TAN (PSX/2.))/(TAN (PSI1/2.))
 255                 WRITE (0,*) ' CELL = ' ,A,PSI1,XN,CELL
 256
 257            ENDIF
 258
 259            IF (IPROJ.EQ.2) THEN
 260
 261                CELL  = A*SIN (PSX)/XN
 262                CELL2 = (1. + COS (PSI1))/(1. + COS (PSX))
 263
 264            ENDIF
 265
 266            R = CELL*(CELL2)**XN
 267
 268            XCNTR = 0.0
 269            YCNTR = -R
 270             WRITE (0,*) ' YCNTR = ' ,CELL,CELL2,XN
 271
 272         ENDIF
 273 !
 274 ! -----FOR MERCATOR PROJECTION, THE PROJECTION IS TRUE AT LAT AT PHI1
 275
 276         IF (IPROJ .EQ. 3) THEN
 277
 278             C2     = A*COS (PSI1)
 279             XCNTR  = 0.0
 280             PHICTR = PHIC/CONV
 281             CELL   = COS (PHICTR)/(1.0+SIN (PHICTR))
 282             YCNTR  = - C2*ALOG (CELL)
 283
 284          ENDIF
 285
 286       WRITE (0,'(3(A,A,F7.2,/))',ADVANCE='no') &
 287       "                "," C2   = ",C2,   &
 288       "                "," Psi1 = ",Psi1, &
 289       "                "," Xcntr= ",Xcntr
 290       WRITE (0,'(1(A,A,F9.2))') &
 291       "                "," Ycntr= ",Ycntr
 292
 293 ! 4.  END
 294 ! =======
 295       RETURN
 296
 297       END SUBROUTINE SET_MAP_PARAMS
 298
 299       SUBROUTINE LLXY (XLATI,XLONI,X,Y)
 300 !------------------------------------------------------------------------------!
 301 !
 302 !                 ROUTINE LLXY
 303 !                **************
 304 !
 305 !
 306 ! PURPOSE:  CALCULATES THE (X,Y) LOCATION (DOT) IN THE MESOSCALE GRIDS
 307 ! -------   FROM LATITUDES AND LONGITUDES
 308 !
 309 !
 310 !  INPUT:
 311 !  -----
 312 !   XLAT:    LATITUDES
 313 !   XLON:    LONGITUDES
 314 !
 315 ! OUTPUT:
 316 ! -----
 317 !   X:        THE COORDINATE IN X (J)-DIRECTION.
 318 !   Y:        THE COORDINATE IN Y (I)-DIRECTION.
 319 !
 320 !------------------------------------------------------------------------------!
 321
 322        IMPLICIT NONE
 323
 324        REAL, INTENT(IN)  :: XLATI, XLONI
 325 !      REAL, INTENT(IN)  :: ID
 326        REAL, INTENT(OUT) :: X, Y
 327        REAL              :: DXLON
 328        REAL              :: XLAT, XLON
 329        REAL              :: XX, YY, XC, YC
 330        REAL              :: CELL, PSI0, PSX, R, FLP
 331        REAL              :: CENTRI, CENTRJ
 332        INTEGER           :: NRATIO
 333
 334        real              :: bb
 335
 336        INCLUDE 'constants.inc'
 337
 338 !------------------------------------------------------------------------------!
 339 !
 340        XLON = XLONI
 341        XLAT = XLATI
 342 !
 343        XLAT = MAX (XLAT, -89.95)
 344        XLAT = MIN (XLAT, +89.95)
 345 !
 346        DXLON = XLON - XLONC
 347        IF (DXLON.GT. 180) DXLON = DXLON - 360.
 348        IF (DXLON.LT.-180) DXLON = DXLON + 360.
 349
 350        IF (IPROJ.EQ.3) THEN
 351 !
 352            XC = XCNTR
 353            YC = YCNTR
 354 !
 355 ! For Mercator projection, high latitude has no definition:
 356            IF (abs(XLAT) .gt.89.95) THEN
 357            PRINT '(A,F8.3,A,F8.3)', &
 358              "**WARNING** MERCATOR PROJECTION, BUT abs(XLAT)= |", &
 359              XLAT,"| > 89.95 ??   XLON=",xlon
 360              Y = -9.9999e8
 361              X = -9.9999e8
 362              RETURN
 363            ENDIF
 364
 365            CELL = COS(XLAT/CONV)/(1.0+SIN(XLAT/CONV))
 366            YY = -C2*ALOG(CELL)
 367            XX = C2*DXLON/CONV
 368
 369        ELSE
 370 !
 371               PSI0 = ( POLE - PHIC )/CONV
 372               XC = 0.0
 373 !
 374 !------ CALCULATE X,Y COORDS. RELATIVE TO POLE
 375 !
 376               FLP = XN*DXLON/CONV
 377
 378               PSX = ( POLE - XLAT )/CONV
 379
 380               if (IPROJ == 2) then
 381                  bb = 2.0*(COS(PSI1/2.0)**2)
 382                  YC = -A*bb*TAN(PSI0/2.0)
 383                   R = -A*bb*TAN(PSX/2.0)
 384               else
 385                  bb = -A/XN*SIN(PSI1)
 386                  YC = bb*(TAN(PSI0/2.0)/TAN(PSI1/2.0))**XN
 387                   R = bb*(TAN(PSX /2.0)/TAN(PSI1/2.0))**XN
 388               endif
 389
 390               IF (PHIC .LT. 0.0) THEN
 391                   XX = R*SIN(FLP)
 392                   YY = R*COS(FLP)
 393               ELSE
 394                   XX = -R*SIN(FLP)
 395                   YY = R*COS(FLP)
 396               END IF
 397
 398        ENDIF
 399 !
 400 !  TRANSFORM (1,1) TO THE ORIGIN
 401 !
 402        CENTRI = FLOAT (IXC + 1)/2.0  ! the location of the center in the coarse
 403        CENTRJ = FLOAT (JXC + 1)/2.0  ! domain
 404 !
 405        X = ( XX - XC )/DIS (1) + CENTRJ  ! the (X,Y) coordinates in the coarse
 406        Y = ( YY - YC )/DIS (1) + CENTRI  ! domain
 407 !
 408        NRATIO = NINT(DIS (1) / DIS (IDD))
 409 !
 410        X = (X - XJM11 (IDD))*NRATIO + 1.0
 411        Y = (Y - XIM11 (IDD))*NRATIO + 1.0
 412 !
 413       END SUBROUTINE LLXY
 414
 415       SUBROUTINE XYLL (XX,YY,XLAT,XLON)
 416 !--------------------------------------------------------------------------CCC
 417 !
 418 !                 ROUTINE XYLL
 419 !                **************
 420 !
 421 !
 422 ! PURPOSE:  CALCULATES THE LATITUDES AND LONGITUDES FROM THE
 423 ! -------   (X,Y) LOCATION (DOT) IN THE MESOSCALE GRIDS.
 424 !
 425 !  INPUT:
 426 !  -----
 427 !  X:         THE COORDINATE IN X (J)-DIRECTION.
 428 !  Y:         THE COORDINATE IN Y (I)-DIRECTION.
 429 !
 430 ! OUTPUT:
 431 ! -----
 432 !
 433 !   XLAT:     LATITUDES
 434 !   XLON:     LONGITUDES
 435 !----------------------------------------------------------------------------CC
 436 !
 437         IMPLICIT NONE
 438 !
 439 ! ARGUMENT
 440 !
 441         REAL, INTENT(IN)  :: XX, YY
 442         REAL, INTENT(OUT) :: XLAT,XLON
 443 !       INTEGER, INTENT(IN) :: ID
 444 !
 445 ! LOCAL
 446 !
 447         REAL :: CNTRI, CNTRJ
 448         REAL :: X, Y
 449         REAL :: FLP, FLPP, R, RXN, CELL, CEL1, CEL2, PSX
 450
 451         INCLUDE 'constants.inc'
 452 !------------------------------------------------------------------------------C
 453
 454         CNTRI = float(IXC+1)/2.
 455         CNTRJ = float(JXC+1)/2.
 456 !
 457 !-----CALCULATE X AND Y POSITIONS OF GRID AT DOT POINT
 458 !
 459        X = XCNTR+(XX-1.0)*DIS(IDD)+(XJM11(IDD)-CNTRJ)*DIS(1)
 460        Y = YCNTR+(YY-1.0)*DIS(IDD)+(XIM11(IDD)-CNTRI)*DIS(1)
 461
 462 !-----NOW CALCULATE LAT AND LON OF THIS POINT
 463 !
 464        IF (IPROJ.NE.3) THEN
 465 !
 466            IF (Y.EQ.0.) THEN
 467 !
 468               IF (X.GE.0.0) FLP =  90.0/CONV
 469               IF (X.LT.0.0) FLP = -90.0/CONV
 470 !
 471            ELSE
 472 !
 473               IF (PHIC.LT.0.0)THEN
 474                   FLP = ATAN2(X,Y)
 475               ELSE
 476                  FLP = ATAN2(X,-Y)
 477               ENDIF
 478 !
 479            ENDIF
 480 !
 481            FLPP = (FLP/XN)*CONV+XLONC
 482 !
 483            IF (FLPP.LT.-180.) FLPP = FLPP + 360
 484            IF (FLPP.GT.180.)  FLPP = FLPP - 360.
 485 !
 486            XLON = FLPP
 487 !
 488 !--------NOW SOLVE FOR LATITUDE
 489 !
 490          R = SQRT(X*X+Y*Y)
 491 !
 492          IF (PHIC.LT.0.0) R = -R
 493 !
 494          IF (IPROJ.EQ.1) THEN
 495              CELL = (R*XN)/(A*SIN(PSI1))
 496              RXN  = 1.0/XN
 497              CEL1 = TAN(PSI1/2.)*(CELL)**RXN
 498          ENDIF
 499 !
 500          IF (IPROJ.EQ.2) THEN
 501              CELL = R/A
 502              CEL1 = CELL/(1.0+COS(PSI1))
 503          ENDIF
 504 !
 505          CEL2 = ATAN(CEL1)
 506          PSX  = 2.*CEL2*CONV
 507          XLAT = POLE-PSX
 508 !
 509        ENDIF
 510 !
 511 !-----CALCULATIONS FOR MERCATOR LAT,LON
 512 !
 513        IF (IPROJ.EQ.3) THEN
 514 !
 515            XLON = XLONC + ((X-XCNTR)/C2)*CONV
 516 !
 517            IF (XLON.LT.-180.) XLON = XLON + 360
 518            IF (XLON.GT.180.)  XLON = XLON - 360.
 519
 520            CELL = EXP(Y/C2)
 521            XLAT = 2.*(CONV*ATAN(CELL))-90.0
 522 !
 523        ENDIF
 524 !
 525        END SUBROUTINE XYLL
 526 !==============================================================================!
 527 END MODULE MODULE_MAP