var/obsproc/src/setup.F90

   1       SUBROUTINE SETUP (domain_check_h, IPROJ, PHIC, XLONC, &
   2                                         TRUELAT1, TRUELAT2, &
   3 #ifdef BKG
   4                         DIS, XCNTR, YCNTR, XN, POLE, PSI1,  C2)
   5 #else
   6                         MAXNES, NESTIX, NESTJX, DIS, NUMC, NESTI, NESTJ, &
   7                         IXC, JXC, XCNTR, YCNTR, XN, POLE, PSI1,  C2, &
   8                         XIM11, XJM11)
   9 #endif
  10 !------------------------------------------------------------------------------!
  11 !
  12 ! MAP PROJECTION FACTORS CALCULATION
  13 !
  14 ! INPUT:
  15 ! -----
  16 !  IPROJ:     PROJECTION TYPE
  17 !  PHIC:      CENTRAL LATITUDE
  18 !  XLONC:     CENTRAL LONGITUDE
  19 !  TRUELAT1:  TRUE LATITUDE 1
  20 !  TRUELAT2:  TRUE LATITUDE 2
  21 !  MAXNES:    NUMBER OF DOMAINS < 10
  22 !  NESTIX (): HORIZONTAL DOMAINS DIMENSIONS IN I DIRECTION
  23 !  NESTJX (): HORIZONTAL DOMAINS DIMENSIONS IN J DIRECTION
  24 !  DIS ():    HORIZONTAL DOMAINS GRID SPACES IN METERS
  25 !  NUMC ():   MOTHER DOMAIN ID (1= COARSE DOMAIN)
  26 !  NESTI ():  DOMAINS' I-INDECES OF LOW LEFT CORNER IN MOTHER DOMAIN
  27 !  NESTJ ():  DOMAINS' J-INDECES OF LOW LEFT CORNER IN MOTHER DOMAIN
  28 !
  29 ! OUTPUT:
  30 ! ------
  31 !  IXC:         COARSE DOMAIN GRID I DIMENSION
  32 !  JXC:         COARSE DOMAIN GRID J DIMENSION
  33 !  XCNTR:       X-COORDINATE OF COARSE DOMAIN CENTER
  34 !  YCNTR:       Y-COORDINATE OF COARSE DOMAIN CENTER
  35 !  XN:          CONE PROJECTION
  36 !  POLE:        POLE LATITUDE
  37 !  PSI1:
  38 !  C2:
  39 !  XIM11 (10):  DOMAINS' INDECES OF LOW LEFT CORNER OF DOMAIN IDD
  40 !  XJM11 (10):  DOMAINS' INDECES OF LOW LEFT CORNER OF DOMAIN IDD
  41 !
  42 !
  43 ! Y.-R. GUO,       September 2000
  44 ! F. VANDENBERGHE, March 2001
  45 !------------------------------------------------------------------------------!
  46
  47       IMPLICIT NONE
  48
  49       INTEGER, INTENT (inout) :: iproj
  50       REAL,    INTENT (inout) :: phic,  xlonc, truelat1, truelat2
  51 #ifdef BKG
  52       REAL,    INTENT (out)   :: pole, xn, psi1, c2, xcntr, ycntr
  53       REAL,    INTENT (inout), DIMENSION (3) :: dis
  54 #else
  55       INTEGER, INTENT (in)    :: maxnes
  56       INTEGER, INTENT (in), DIMENSION (maxnes) :: nestix, nestjx, NUMC, &
  57                                               nesti, nestj
  58       REAL,    INTENT (inout), DIMENSION (maxnes) :: dis
  59
  60       INTEGER, INTENT (out)   :: ixc, jxc
  61
  62       REAL,    INTENT (out)   :: pole, xn
  63       REAL,    INTENT (out)   :: xim11 (maxnes), xjm11 (maxnes)
  64       REAL,    INTENT (out)   :: psi1, c2, xcntr, ycntr
  65 #endif
  66 !------------------------------------------------------------------------------!
  67 #ifdef BKG
  68       REAL :: psx, cell, cell2, r, phictr, sign
  69
  70       INTEGER, DIMENSION (3) :: nproj
  71
  72       CHARACTER (LEN=80) :: project
  73 #else
  74 !     REAL :: conv, a
  75       REAL :: iexp, aexp
  76       REAL :: psx, cell, cell2, r, phictr, hdsize, sign
  77       REAL :: xjc, xdis, yjc, ydis
  78
  79       REAL :: cntrj0, cntri0
  80       INTEGER :: ixcmax, jxcmax
  81       REAL, DIMENSION (maxnes) :: xsouth,xwest,xnorth,xeast
  82
  83       INTEGER :: incr
  84       INTEGER :: m, nm, nmc, nd1, nd2, mismatch
  85       INTEGER :: ioffst, joffst, iimxn, jjmxn
  86       INTEGER :: iratio (maxnes), nratio (maxnes)
  87
  88       INTEGER, DIMENSION (3) :: nproj
  89
  90       CHARACTER (LEN=80) :: project
  91 !------------------------------------------------------------------------------!
  92 !     DATA conv  / 57.29578 /
  93 !     DATA a     / 6370.    /
  94 #endif
  95       include 'constants.inc'
  96
  97       logical  domain_check_h
  98       DATA nproj / 1, 2, 3  /
  99 !     DATA nproj / 'LAMCON', 'POLSTR', 'MERCAT' /
 100 !------------------------------------------------------------------------------!
 101       if (iproj == 0 ) then
 102           PROJECT='CE'
 103           XCNTR = 0
 104           YCNTR = 0
 105 !          DIS(1) = 2.0 * pi * a / (NESTJX (1)-1)
 106           DIS(1) = 110.0
 107           if (domain_check_h) STOP 'domain_check_h'
 108           goto 1100
 109       endif
 110
 111 !     WRITE (0,'(A)') " SET UP MAP PARAMETERS"
 112
 113 #ifdef BKG
 114 ! DEFINE THE PARAMETERS OF MAP BASED ON THE IPROJ:
 115
 116       XN = -1.0E08
 117
 118       IF (PHIC.LT.0) THEN
 119          SIGN=-1.      ! SOUTH HEMESPHERE
 120
 121       ELSE
 122
 123          SIGN=1.       ! NORTH HEMESPHERE
 124
 125       ENDIF
 126
 127       POLE = 90.
 128
 129       IF (IPROJ .EQ. NPROJ(1)) THEN
 130           if (abs(TRUELAT1 - TRUELAT2) .gt. 0.1) then
 131             XN = ALOG10 (COS (TRUELAT1 / CONV)) - &
 132                  ALOG10 (COS (TRUELAT2 / CONV))
 133
 134             XN = XN/(ALOG10 (TAN ((45.0 - SIGN*TRUELAT1/2.0) / CONV)) - &
 135                      ALOG10 (TAN ((45.0 - SIGN*TRUELAT2/2.0) / CONV)))
 136           else
 137             XN = sign*sin(truelat1 / conv)
 138           endif
 139
 140           PSI1 = 90.-SIGN*TRUELAT1
 141
 142           PROJECT='LC'
 143
 144       ELSE IF (IPROJ .EQ. NPROJ(2)) THEN
 145
 146           XN = 1.0
 147
 148
 149 ! PSI1 IS THE PSEUDO-LATITUDE
 150
 151           PSI1 = 90.-SIGN*TRUELAT1
 152
 153           PROJECT = 'ST'
 154
 155       ELSE IF (IPROJ .EQ. NPROJ(3)) THEN
 156
 157           XN = 0.
 158
 159           PSI1 = 0.
 160           PROJECT = 'ME'
 161
 162       END IF
 163
 164       PSI1 = PSI1 / CONV
 165
 166       IF (PHIC .LT. 0.) THEN
 167           PSI1 = -PSI1
 168           POLE = -POLE
 169
 170       ENDIF
 171
 172 !--------CALCULATE R
 173
 174       IF (IPROJ .NE. NPROJ(3)) THEN
 175
 176           PSX = (POLE-PHIC)/CONV
 177
 178           IF (IPROJ .EQ. NPROJ(1)) THEN
 179
 180               CELL  = A*SIN (PSI1)/XN
 181               CELL2 = (TAN (PSX/2.)) / (TAN (PSI1/2.))
 182          ENDIF
 183
 184          IF (IPROJ .EQ. NPROJ(2)) THEN
 185             CELL  = A*SIN (PSX)/XN
 186             CELL2 = (1. + COS (PSI1))/(1. + COS (PSX))
 187          ENDIF
 188
 189          R = CELL*(CELL2)**XN
 190
 191          XCNTR = 0.0
 192          YCNTR = -R
 193
 194       ENDIF
 195
 196 !-----FOR MERCATOR PROJECTION, THE PROJECTION IS TRUE AT LAT AT PHI1
 197
 198       IF (IPROJ .EQ. NPROJ(3)) THEN
 199
 200          C2     = A*COS (PSI1)
 201          XCNTR  = 0.0
 202          PHICTR = PHIC/CONV
 203          CELL   = COS (PHICTR)/(1.0+SIN (PHICTR))
 204          YCNTR  = - C2*ALOG (CELL)
 205
 206       ENDIF
 207 1100  continue
 208
 209 #else
 210
 211 ! REMOVE DOMAIN EXTENSION IN MM5 INPUT/OUTPUT FILES
 212
 213       IEXP = 0
 214       AEXP =  -1.001*DIS(1)
 215
 216 ! REMOVE INITIALIZATION OF TRUE LAT
 217
 218 !     TRUELAT1=91.0
 219 !     TRUELAT2=91.0
 220
 221 ! DEFINE THE PARAMETERS OF MAP BASED ON THE IPROJ:
 222
 223       XN = -1.0E08
 224
 225       IF (PHIC.LT.0) THEN
 226          SIGN=-1.      ! SOUTH HEMESPHERE
 227
 228       ELSE
 229
 230          SIGN=1.       ! NORTH HEMESPHERE
 231
 232       ENDIF
 233
 234       POLE = 90.
 235
 236       IF (IPROJ .EQ. NPROJ(1)) THEN
 237
 238           IF (ABS (TRUELAT1) .GT. 90.) THEN
 239               TRUELAT1 = 60.
 240               TRUELAT2 = 30.
 241               TRUELAT1 = SIGN*TRUELAT1
 242               TRUELAT2 = SIGN*TRUELAT2
 243           ENDIF
 244
 245           IF (ABS(TRUELAT1 - TRUELAT2) .GT. 0.1) then
 246             XN = ALOG10 (COS (TRUELAT1 / CONV)) - &
 247                  ALOG10 (COS (TRUELAT2 / CONV))
 248
 249             XN = XN/(ALOG10 (TAN ((45.0 - SIGN*TRUELAT1/2.0) / CONV)) - &
 250                      ALOG10 (TAN ((45.0 - SIGN*TRUELAT2/2.0) / CONV)))
 251           ELSE
 252             XN = SIGN*SIN(TRUELAT1 / CONV)
 253           ENDIF
 254
 255           PSI1 = 90.-SIGN*TRUELAT1
 256
 257           PROJECT='LC'
 258
 259       ELSE IF (IPROJ .EQ. NPROJ(2)) THEN
 260
 261           XN = 1.0
 262
 263           IF (ABS(TRUELAT1) .GT. 90.) THEN
 264
 265               TRUELAT1 = 60.
 266               TRUELAT2 = 0.
 267
 268               TRUELAT1 = SIGN*TRUELAT1
 269               TRUELAT2 = SIGN*TRUELAT2
 270
 271           ENDIF
 272
 273
 274 ! PSI1 IS THE PSEUDO-LATITUDE
 275
 276           PSI1 = 90.-SIGN*TRUELAT1
 277
 278           PROJECT = 'ST'
 279
 280       ELSE IF (IPROJ .EQ. NPROJ(3)) THEN
 281
 282           XN = 0.
 283
 284           IF (ABS (TRUELAT1) .GT. 90.) THEN
 285
 286               TRUELAT1 = 0.
 287               TRUELAT2 = 0.
 288           ENDIF
 289
 290           IF (TRUELAT1 .NE. 0.) THEN
 291               WRITE (0,'(/,A)') &
 292              "MERCATOR PROJECTION IS ONLY SUPPORTED AT 0 DEGREE TRUE LATITUDE."
 293               WRITE (0,'(A,/)') &
 294              "TRUELAT1 IS RESET TO 0"
 295
 296                TRUELAT1 = 0.
 297           ENDIF
 298
 299           PSI1 = 0.
 300           PROJECT = 'ME'
 301
 302       END IF
 303
 304       PSI1 = PSI1 / CONV
 305
 306       IF (PHIC .LT. 0.) THEN
 307           PSI1 = -PSI1
 308           POLE = -POLE
 309
 310       ENDIF
 311
 312 !--------CALCULATE R
 313
 314       IF (IPROJ .NE. NPROJ(3)) THEN
 315
 316           PSX = (POLE-PHIC)/CONV
 317
 318           IF (IPROJ .EQ. NPROJ(1)) THEN
 319
 320               CELL  = A*SIN (PSI1)/XN
 321               CELL2 = (TAN (PSX/2.)) / (TAN (PSI1/2.))
 322          ENDIF
 323
 324          IF (IPROJ .EQ. NPROJ(2)) THEN
 325             CELL  = A*SIN (PSX)/XN
 326             CELL2 = (1. + COS (PSI1))/(1. + COS (PSX))
 327          ENDIF
 328
 329          R = CELL*(CELL2)**XN
 330
 331          XCNTR = 0.0
 332          YCNTR = -R
 333
 334       ENDIF
 335
 336 !-----FOR MERCATOR PROJECTION, THE PROJECTION IS TRUE AT LAT AT PHI1
 337
 338       IF (IPROJ .EQ. NPROJ(3)) THEN
 339
 340          C2     = A*COS (PSI1)
 341          XCNTR  = 0.0
 342          PHICTR = PHIC/CONV
 343          CELL   = COS (PHICTR)/(1.0+SIN (PHICTR))
 344          YCNTR  = - C2*ALOG (CELL)
 345
 346       ENDIF
 347
 348 1100  continue
 349
 350       WRITE (0,'(2(A,F8.1),A,2X,A,f10.3)') &
 351       " COARSE GRID CENTER IS AT X = ",XCNTR," KM AND Y = ",YCNTR," KM.", &
 352       " DIS(1)=", DIS(1)
 353
 354
 355 !   CHECK THE COMPATIBILITY OF NEST DOMAINS WITH THE COARSE DOMAINS
 356 !     AND CALCULATE THE IRATIOS, INORTHS, ISOUTHS, IWESTS AND IEASTS
 357
 358 !     A) EXTENDING THE COARSE DOMAIN IF IEXP = 1
 359
 360         IXC = NESTIX (1)
 361         JXC = NESTJX (1)
 362
 363         IOFFST = 0
 364         JOFFST = 0
 365
 366       IF (IEXP .EQ. 1) THEN
 367
 368           INCR = INT (AEXP/DIS (1) + 1.001)
 369
 370           IXC = NESTIX (1) + INCR*2
 371           JXC = NESTJX (1) + INCR*2
 372
 373           IOFFST = INCR
 374           JOFFST = INCR
 375
 376           WRITE (0,'(A,I3)') " GRID IS EXPANDED BY ",INCR, &
 377                              " GRID POINTS ON EACH EDGE."
 378 !          WRITE (0,'(2(A,I3))') &
 379 !         "EXTENDED IXC IS ",IXC," EXTENDED JXC IS ",JXC
 380
 381       ENDIF
 382
 383 !-----CENTER OF GRID IN THE COARSE DOMAIN
 384
 385       CNTRJ0 = FLOAT (JXC+1)/2.
 386       CNTRI0 = FLOAT (IXC+1)/2.
 387
 388 !      WRITE (0,'(2(A,I5))') &
 389 !     "COARSE DOMAIN SIZE IX = ",NESTIX(1)," JX = ", NESTJX(1)
 390 !  MIX, MJX ARE USED IN SUB. TFUDGE:
 391
 392       IXCMAX = IXC
 393       JXCMAX = JXC
 394
 395       DO M = 1, MAXNES
 396          IXCMAX = MAX0 (NESTIX(M),IXCMAX)
 397          JXCMAX = MAX0 (NESTJX(M),JXCMAX)
 398       ENDDO
 399
 400 !      PRINT 24, IXCMAX,JXCMAX
 401 !24    FORMAT('   ++> THE MAXIMUM DIMENSION = (',2I5,') <++')
 402
 403 !  CHECK IF POLE IS INSIDE THE DOMAIN OR NOT FOR LAMBERT PROJECTION:
 404
 405       HDSIZE = (IXC-1)*DIS(1)/2.
 406
 407 !     IF (HDSIZE.GT.ABS(YCNTR) .AND. IPROJ.EQ.NPROJ(1))THEN
 408 !        PRINT *,'-------------------------------------------------'
 409 !        PRINT *,'HALF DOMAIN SIZE IN Y-DIRECTION = ',HDSIZE
 410 !        PRINT *,'    DISTANCE FROM CNTER TO POLE = ',ABS(YCNTR)
 411 !        PRINT *,'NOT MAKE SENSE WITH THE POLE IS INSIDE THE DOMAIN '
 412 !        PRINT *,'    FOR LAMBERT CONFORMAL PROJECTION!'
 413 !        PRINT *,'=== PLEASE RE-SPECIFY THE CENTER OR DOMAIN SIZE. ==='
 414 !        STOP
 415 !     ENDIF
 416
 417 !     B) CALCULATING THE IRATIOS, INORTHS AND IEASTS:
 418
 419       IRATIO (1) = 1
 420       NRATIO (1) = 1
 421       XSOUTH (1) = 1.
 422       XWEST  (1) = 1.
 423       XNORTH (1) = FLOAT (IXC)
 424       XEAST  (1) = FLOAT (JXC)
 425
 426       XJC = (XEAST(1) + 1.0)/2.
 427
 428 !      PRINT 27,XSOUTH(1),XWEST(1),XNORTH(1),XEAST(1),DIS(1),&
 429 !               IRATIO(1),NRATIO(1)
 430 ! 27   FORMAT(1X,'XSOUTH(1)= ',F6.1,2X,'XWEST(1)= ',F6.1,2X, &
 431 !      'XNORTH(1)= ',F6.1,2X,'XEAST(1)= ',F6.1,2X,'DIS(1)= ',F6.1,&
 432 !      2X,'IRATIO(1)= ',I3,2X,'NRATIO(1)= ',I3)
 433
 434       MISMATCH = 0
 435
 436       DO 30 NM = 2, MAXNES
 437
 438 !  DOMAINS' CONSISTENCY CHECK:
 439
 440       NMC = NUMC (NM)
 441
 442       IF (AMOD ((DIS (NMC)+0.09),DIS (NM)) .GT. 0.1) THEN
 443
 444          MISMATCH = MISMATCH + 1
 445
 446 !        PRINT 29, NM,NMC
 447 !        PRINT 31,NM,DIS(NM),NMC,DIS(NMC)
 448 !29      FORMAT(2X,'DOMAIN ',I2,' HAS INCORRECT GRID SIZE ', &
 449 !      '  IT HAS TO BE THE MULTIPLE OF DOMAIN ',I2)
 450 !31      FORMAT(2X,'DOMAIN ',I2,' GRID SIZE= ',F6.1,' KM', &
 451 !               '  DOMAIN ',I2,' GRID SIZE= ',F6.1,' KM')
 452
 453         GO TO 30
 454
 455       ENDIF
 456
 457       IRATIO (NM) = NINT (DIS (NMC)/ DIS (NM))
 458       NRATIO (NM) = NINT (DIS (1)  / DIS (NM))
 459
 460 !   MAKE SURE THE 4 CORNER POINTS OF NEST DOMAINS ARE ON THE
 461 !   PREVIOUS DOMAIN GRID-POINTS
 462
 463       IF (MOD((NESTIX(NM)-1),IRATIO(NM)).NE.0) THEN
 464
 465         MISMATCH = MISMATCH + 1
 466
 467         IIMXN = (INT(FLOAT(NESTIX(NM)-1)/IRATIO(NM))+1)*IRATIO(NM) + 1
 468
 469 !        PRINT 32,NM,NESTIX(NM),IRATIO(NM),IIMXN
 470 !32      FORMAT(2X,'NESTIX(',I2,')=',I4,' AND IRATIO=',I2,&
 471 !      ' DOES NOT MATCH, YOU MAY SET NESTIX TO ',I4)
 472
 473       ENDIF
 474
 475       IF (MOD ((NESTJX (NM)-1),IRATIO (NM)).NE.0) THEN
 476
 477          MISMATCH = MISMATCH + 1
 478
 479         JJMXN = (INT(FLOAT(NESTJX(NM)-1)/IRATIO(NM))+1)*IRATIO(NM) + 1
 480
 481 !        PRINT 33,NM,NESTJX(NM),IRATIO(NM),JJMXN
 482 !33      FORMAT(2X,'NESTJX(',I2,')=',I4,' AND IRATIO=',I2,&
 483 !      ' DOES NOT MATCH, YOU MAY SET NESTJX TO ',I4)
 484       ENDIF
 485 !
 486 !-----REDEFINE LOCATION OF LOWER LEFT CORNER OF FINE MESH (IN TERMS
 487 !     OF EXTENDED COARSE MESH - DOMAIN 1 INDICES) IF USING EXTENDED G
 488
 489 !      PRINT 34,NM,NESTIX(NM),NESTJX(NM),DIS(NM),NESTI(NM),NESTJ(NM),&
 490 !      NUMC(NM),IRATIO(NM),NRATIO(NM),IEXP
 491 !34    FORMAT(/1X,'DOMAIN ',I2,2X,'IX=',I3,2X,'JX=',I3,2X, &
 492 !      'DS= ',F6.1,2X,'ICNS=',I4,2X,'JCNS=',I4,2X,&
 493 !      'NUMC= ',I2,2X,'IRATIO= ',I2,2X,'NRATIO= ',I2,2X,&
 494 !      'IEXP= ',I2)
 495
 496       XDIS = 0.0
 497       YDIS = 0.0
 498       ND1 = NM
 499       ND2 = NMC
 500
 501 40    CONTINUE
 502
 503       XDIS = (NESTI(ND1)-1)*DIS(ND2) + XDIS
 504       YDIS = (NESTJ(ND1)-1)*DIS(ND2) + YDIS
 505
 506       IF (ND2 .GT. 1) THEN
 507         ND1 = ND2
 508         ND2 = NUMC(ND2)
 509         GO TO 40
 510       ENDIF
 511
 512       XSOUTH(NM) = XDIS/DIS(1) + FLOAT(IOFFST) + 1
 513       XWEST(NM)  = YDIS/DIS(1) + FLOAT(JOFFST) + 1
 514       XNORTH(NM) = XSOUTH(NM) + FLOAT(NESTIX(NM)-1)*DIS(NM)/DIS(1)
 515       XEAST(NM)  = XWEST(NM) + FLOAT(NESTJX(NM)-1)*DIS(NM)/DIS(1)
 516
 517 !      PRINT 35
 518 !      PRINT 36,XSOUTH(NM),XWEST(NM),XNORTH(NM),XEAST(NM)
 519 !35    FORMAT(2X,'COARSE MESH INDICES FOR THE 4 CORNER POINTS ARE')
 520 !36    FORMAT(2X,'SOUTH=',F6.1,3X,'WEST=',F6.1,3X,'NORTH =',F6.1,3X, &
 521 !            'EAST = ',F6.1)
 522
 523 30    CONTINUE
 524
 525       IF (MISMATCH.GT.0) THEN
 526
 527 !       PRINT *, &
 528 !      'TERRAIN STOP IN SUBROUTINE SETUP DUE TO INCORRECT NEST DOMAIN SET UP'
 529 !       STOP 1111
 530
 531        ENDIF
 532
 533 !      Copy output
 534
 535
 536        XIM11 = XSOUTH
 537        XJM11 = XWEST
 538
 539 #endif
 540        END SUBROUTINE SETUP