SUBROUTINE SETUP (domain_check_h, IPROJ, PHIC, XLONC, & TRUELAT1, TRUELAT2, & #ifdef BKG DIS, XCNTR, YCNTR, XN, POLE, PSI1, C2) #else MAXNES, NESTIX, NESTJX, DIS, NUMC, NESTI, NESTJ, & IXC, JXC, XCNTR, YCNTR, XN, POLE, PSI1, C2, & XIM11, XJM11) #endif !------------------------------------------------------------------------------! ! ! MAP PROJECTION FACTORS CALCULATION ! ! INPUT: ! ----- ! IPROJ: PROJECTION TYPE ! PHIC: CENTRAL LATITUDE ! XLONC: CENTRAL LONGITUDE ! TRUELAT1: TRUE LATITUDE 1 ! TRUELAT2: TRUE LATITUDE 2 ! MAXNES: NUMBER OF DOMAINS < 10 ! NESTIX (): HORIZONTAL DOMAINS DIMENSIONS IN I DIRECTION ! NESTJX (): HORIZONTAL DOMAINS DIMENSIONS IN J DIRECTION ! DIS (): HORIZONTAL DOMAINS GRID SPACES IN METERS ! NUMC (): MOTHER DOMAIN ID (1= COARSE DOMAIN) ! NESTI (): DOMAINS' I-INDECES OF LOW LEFT CORNER IN MOTHER DOMAIN ! NESTJ (): DOMAINS' J-INDECES OF LOW LEFT CORNER IN MOTHER DOMAIN ! ! OUTPUT: ! ------ ! IXC: COARSE DOMAIN GRID I DIMENSION ! JXC: COARSE DOMAIN GRID J DIMENSION ! XCNTR: X-COORDINATE OF COARSE DOMAIN CENTER ! YCNTR: Y-COORDINATE OF COARSE DOMAIN CENTER ! XN: CONE PROJECTION ! POLE: POLE LATITUDE ! PSI1: ! C2: ! XIM11 (10): DOMAINS' INDECES OF LOW LEFT CORNER OF DOMAIN IDD ! XJM11 (10): DOMAINS' INDECES OF LOW LEFT CORNER OF DOMAIN IDD ! ! ! Y.-R. GUO, September 2000 ! F. VANDENBERGHE, March 2001 !------------------------------------------------------------------------------! IMPLICIT NONE INTEGER, INTENT (inout) :: iproj REAL, INTENT (inout) :: phic, xlonc, truelat1, truelat2 #ifdef BKG REAL, INTENT (out) :: pole, xn, psi1, c2, xcntr, ycntr REAL, INTENT (inout), DIMENSION (3) :: dis #else INTEGER, INTENT (in) :: maxnes INTEGER, INTENT (in), DIMENSION (maxnes) :: nestix, nestjx, NUMC, & nesti, nestj REAL, INTENT (inout), DIMENSION (maxnes) :: dis INTEGER, INTENT (out) :: ixc, jxc REAL, INTENT (out) :: pole, xn REAL, INTENT (out) :: xim11 (maxnes), xjm11 (maxnes) REAL, INTENT (out) :: psi1, c2, xcntr, ycntr #endif !------------------------------------------------------------------------------! #ifdef BKG REAL :: psx, cell, cell2, r, phictr, sign INTEGER, DIMENSION (3) :: nproj CHARACTER (LEN=80) :: project #else ! REAL :: conv, a REAL :: iexp, aexp REAL :: psx, cell, cell2, r, phictr, hdsize, sign REAL :: xjc, xdis, yjc, ydis REAL :: cntrj0, cntri0 INTEGER :: ixcmax, jxcmax REAL, DIMENSION (maxnes) :: xsouth,xwest,xnorth,xeast INTEGER :: incr INTEGER :: m, nm, nmc, nd1, nd2, mismatch INTEGER :: ioffst, joffst, iimxn, jjmxn INTEGER :: iratio (maxnes), nratio (maxnes) INTEGER, DIMENSION (3) :: nproj CHARACTER (LEN=80) :: project !------------------------------------------------------------------------------! ! DATA conv / 57.29578 / ! DATA a / 6370. / #endif include 'constants.inc' logical domain_check_h DATA nproj / 1, 2, 3 / ! DATA nproj / 'LAMCON', 'POLSTR', 'MERCAT' / !------------------------------------------------------------------------------! if (iproj == 0 ) then PROJECT='CE' XCNTR = 0 YCNTR = 0 ! DIS(1) = 2.0 * pi * a / (NESTJX (1)-1) DIS(1) = 110.0 if (domain_check_h) STOP 'domain_check_h' goto 1100 endif ! WRITE (0,'(A)') " SET UP MAP PARAMETERS" #ifdef BKG ! DEFINE THE PARAMETERS OF MAP BASED ON THE IPROJ: XN = -1.0E08 IF (PHIC.LT.0) THEN SIGN=-1. ! SOUTH HEMESPHERE ELSE SIGN=1. ! NORTH HEMESPHERE ENDIF POLE = 90. IF (IPROJ .EQ. NPROJ(1)) THEN if (abs(TRUELAT1 - TRUELAT2) .gt. 0.1) then XN = ALOG10 (COS (TRUELAT1 / CONV)) - & ALOG10 (COS (TRUELAT2 / CONV)) XN = XN/(ALOG10 (TAN ((45.0 - SIGN*TRUELAT1/2.0) / CONV)) - & ALOG10 (TAN ((45.0 - SIGN*TRUELAT2/2.0) / CONV))) else XN = sign*sin(truelat1 / conv) endif PSI1 = 90.-SIGN*TRUELAT1 PROJECT='LC' ELSE IF (IPROJ .EQ. NPROJ(2)) THEN XN = 1.0 ! PSI1 IS THE PSEUDO-LATITUDE PSI1 = 90.-SIGN*TRUELAT1 PROJECT = 'ST' ELSE IF (IPROJ .EQ. NPROJ(3)) THEN XN = 0. PSI1 = 0. PROJECT = 'ME' END IF PSI1 = PSI1 / CONV IF (PHIC .LT. 0.) THEN PSI1 = -PSI1 POLE = -POLE ENDIF !--------CALCULATE R IF (IPROJ .NE. NPROJ(3)) THEN PSX = (POLE-PHIC)/CONV IF (IPROJ .EQ. NPROJ(1)) THEN CELL = A*SIN (PSI1)/XN CELL2 = (TAN (PSX/2.)) / (TAN (PSI1/2.)) ENDIF IF (IPROJ .EQ. NPROJ(2)) THEN CELL = A*SIN (PSX)/XN CELL2 = (1. + COS (PSI1))/(1. + COS (PSX)) ENDIF R = CELL*(CELL2)**XN XCNTR = 0.0 YCNTR = -R ENDIF !-----FOR MERCATOR PROJECTION, THE PROJECTION IS TRUE AT LAT AT PHI1 IF (IPROJ .EQ. NPROJ(3)) THEN C2 = A*COS (PSI1) XCNTR = 0.0 PHICTR = PHIC/CONV CELL = COS (PHICTR)/(1.0+SIN (PHICTR)) YCNTR = - C2*ALOG (CELL) ENDIF 1100 continue #else ! REMOVE DOMAIN EXTENSION IN MM5 INPUT/OUTPUT FILES IEXP = 0 AEXP = -1.001*DIS(1) ! REMOVE INITIALIZATION OF TRUE LAT ! TRUELAT1=91.0 ! TRUELAT2=91.0 ! DEFINE THE PARAMETERS OF MAP BASED ON THE IPROJ: XN = -1.0E08 IF (PHIC.LT.0) THEN SIGN=-1. ! SOUTH HEMESPHERE ELSE SIGN=1. ! NORTH HEMESPHERE ENDIF POLE = 90. IF (IPROJ .EQ. NPROJ(1)) THEN IF (ABS (TRUELAT1) .GT. 90.) THEN TRUELAT1 = 60. TRUELAT2 = 30. TRUELAT1 = SIGN*TRUELAT1 TRUELAT2 = SIGN*TRUELAT2 ENDIF IF (ABS(TRUELAT1 - TRUELAT2) .GT. 0.1) then XN = ALOG10 (COS (TRUELAT1 / CONV)) - & ALOG10 (COS (TRUELAT2 / CONV)) XN = XN/(ALOG10 (TAN ((45.0 - SIGN*TRUELAT1/2.0) / CONV)) - & ALOG10 (TAN ((45.0 - SIGN*TRUELAT2/2.0) / CONV))) ELSE XN = SIGN*SIN(TRUELAT1 / CONV) ENDIF PSI1 = 90.-SIGN*TRUELAT1 PROJECT='LC' ELSE IF (IPROJ .EQ. NPROJ(2)) THEN XN = 1.0 IF (ABS(TRUELAT1) .GT. 90.) THEN TRUELAT1 = 60. TRUELAT2 = 0. TRUELAT1 = SIGN*TRUELAT1 TRUELAT2 = SIGN*TRUELAT2 ENDIF ! PSI1 IS THE PSEUDO-LATITUDE PSI1 = 90.-SIGN*TRUELAT1 PROJECT = 'ST' ELSE IF (IPROJ .EQ. NPROJ(3)) THEN XN = 0. IF (ABS (TRUELAT1) .GT. 90.) THEN TRUELAT1 = 0. TRUELAT2 = 0. ENDIF IF (TRUELAT1 .NE. 0.) THEN WRITE (0,'(/,A)') & "MERCATOR PROJECTION IS ONLY SUPPORTED AT 0 DEGREE TRUE LATITUDE." WRITE (0,'(A,/)') & "TRUELAT1 IS RESET TO 0" TRUELAT1 = 0. ENDIF PSI1 = 0. PROJECT = 'ME' END IF PSI1 = PSI1 / CONV IF (PHIC .LT. 0.) THEN PSI1 = -PSI1 POLE = -POLE ENDIF !--------CALCULATE R IF (IPROJ .NE. NPROJ(3)) THEN PSX = (POLE-PHIC)/CONV IF (IPROJ .EQ. NPROJ(1)) THEN CELL = A*SIN (PSI1)/XN CELL2 = (TAN (PSX/2.)) / (TAN (PSI1/2.)) ENDIF IF (IPROJ .EQ. NPROJ(2)) THEN CELL = A*SIN (PSX)/XN CELL2 = (1. + COS (PSI1))/(1. + COS (PSX)) ENDIF R = CELL*(CELL2)**XN XCNTR = 0.0 YCNTR = -R ENDIF !-----FOR MERCATOR PROJECTION, THE PROJECTION IS TRUE AT LAT AT PHI1 IF (IPROJ .EQ. NPROJ(3)) THEN C2 = A*COS (PSI1) XCNTR = 0.0 PHICTR = PHIC/CONV CELL = COS (PHICTR)/(1.0+SIN (PHICTR)) YCNTR = - C2*ALOG (CELL) ENDIF 1100 continue WRITE (0,'(2(A,F8.1),A,2X,A,f10.3)') & " COARSE GRID CENTER IS AT X = ",XCNTR," KM AND Y = ",YCNTR," KM.", & " DIS(1)=", DIS(1) ! CHECK THE COMPATIBILITY OF NEST DOMAINS WITH THE COARSE DOMAINS ! AND CALCULATE THE IRATIOS, INORTHS, ISOUTHS, IWESTS AND IEASTS ! A) EXTENDING THE COARSE DOMAIN IF IEXP = 1 IXC = NESTIX (1) JXC = NESTJX (1) IOFFST = 0 JOFFST = 0 IF (IEXP .EQ. 1) THEN INCR = INT (AEXP/DIS (1) + 1.001) IXC = NESTIX (1) + INCR*2 JXC = NESTJX (1) + INCR*2 IOFFST = INCR JOFFST = INCR WRITE (0,'(A,I3)') " GRID IS EXPANDED BY ",INCR, & " GRID POINTS ON EACH EDGE." ! WRITE (0,'(2(A,I3))') & ! "EXTENDED IXC IS ",IXC," EXTENDED JXC IS ",JXC ENDIF !-----CENTER OF GRID IN THE COARSE DOMAIN CNTRJ0 = FLOAT (JXC+1)/2. CNTRI0 = FLOAT (IXC+1)/2. ! WRITE (0,'(2(A,I5))') & ! "COARSE DOMAIN SIZE IX = ",NESTIX(1)," JX = ", NESTJX(1) ! MIX, MJX ARE USED IN SUB. TFUDGE: IXCMAX = IXC JXCMAX = JXC DO M = 1, MAXNES IXCMAX = MAX0 (NESTIX(M),IXCMAX) JXCMAX = MAX0 (NESTJX(M),JXCMAX) ENDDO ! PRINT 24, IXCMAX,JXCMAX !24 FORMAT(' ++> THE MAXIMUM DIMENSION = (',2I5,') <++') ! CHECK IF POLE IS INSIDE THE DOMAIN OR NOT FOR LAMBERT PROJECTION: HDSIZE = (IXC-1)*DIS(1)/2. ! IF (HDSIZE.GT.ABS(YCNTR) .AND. IPROJ.EQ.NPROJ(1))THEN ! PRINT *,'-------------------------------------------------' ! PRINT *,'HALF DOMAIN SIZE IN Y-DIRECTION = ',HDSIZE ! PRINT *,' DISTANCE FROM CNTER TO POLE = ',ABS(YCNTR) ! PRINT *,'NOT MAKE SENSE WITH THE POLE IS INSIDE THE DOMAIN ' ! PRINT *,' FOR LAMBERT CONFORMAL PROJECTION!' ! PRINT *,'=== PLEASE RE-SPECIFY THE CENTER OR DOMAIN SIZE. ===' ! STOP ! ENDIF ! B) CALCULATING THE IRATIOS, INORTHS AND IEASTS: IRATIO (1) = 1 NRATIO (1) = 1 XSOUTH (1) = 1. XWEST (1) = 1. XNORTH (1) = FLOAT (IXC) XEAST (1) = FLOAT (JXC) XJC = (XEAST(1) + 1.0)/2. ! PRINT 27,XSOUTH(1),XWEST(1),XNORTH(1),XEAST(1),DIS(1),& ! IRATIO(1),NRATIO(1) ! 27 FORMAT(1X,'XSOUTH(1)= ',F6.1,2X,'XWEST(1)= ',F6.1,2X, & ! 'XNORTH(1)= ',F6.1,2X,'XEAST(1)= ',F6.1,2X,'DIS(1)= ',F6.1,& ! 2X,'IRATIO(1)= ',I3,2X,'NRATIO(1)= ',I3) MISMATCH = 0 DO 30 NM = 2, MAXNES ! DOMAINS' CONSISTENCY CHECK: NMC = NUMC (NM) IF (AMOD ((DIS (NMC)+0.09),DIS (NM)) .GT. 0.1) THEN MISMATCH = MISMATCH + 1 ! PRINT 29, NM,NMC ! PRINT 31,NM,DIS(NM),NMC,DIS(NMC) !29 FORMAT(2X,'DOMAIN ',I2,' HAS INCORRECT GRID SIZE ', & ! ' IT HAS TO BE THE MULTIPLE OF DOMAIN ',I2) !31 FORMAT(2X,'DOMAIN ',I2,' GRID SIZE= ',F6.1,' KM', & ! ' DOMAIN ',I2,' GRID SIZE= ',F6.1,' KM') GO TO 30 ENDIF IRATIO (NM) = NINT (DIS (NMC)/ DIS (NM)) NRATIO (NM) = NINT (DIS (1) / DIS (NM)) ! MAKE SURE THE 4 CORNER POINTS OF NEST DOMAINS ARE ON THE ! PREVIOUS DOMAIN GRID-POINTS IF (MOD((NESTIX(NM)-1),IRATIO(NM)).NE.0) THEN MISMATCH = MISMATCH + 1 IIMXN = (INT(FLOAT(NESTIX(NM)-1)/IRATIO(NM))+1)*IRATIO(NM) + 1 ! PRINT 32,NM,NESTIX(NM),IRATIO(NM),IIMXN !32 FORMAT(2X,'NESTIX(',I2,')=',I4,' AND IRATIO=',I2,& ! ' DOES NOT MATCH, YOU MAY SET NESTIX TO ',I4) ENDIF IF (MOD ((NESTJX (NM)-1),IRATIO (NM)).NE.0) THEN MISMATCH = MISMATCH + 1 JJMXN = (INT(FLOAT(NESTJX(NM)-1)/IRATIO(NM))+1)*IRATIO(NM) + 1 ! PRINT 33,NM,NESTJX(NM),IRATIO(NM),JJMXN !33 FORMAT(2X,'NESTJX(',I2,')=',I4,' AND IRATIO=',I2,& ! ' DOES NOT MATCH, YOU MAY SET NESTJX TO ',I4) ENDIF ! !-----REDEFINE LOCATION OF LOWER LEFT CORNER OF FINE MESH (IN TERMS ! OF EXTENDED COARSE MESH - DOMAIN 1 INDICES) IF USING EXTENDED G ! PRINT 34,NM,NESTIX(NM),NESTJX(NM),DIS(NM),NESTI(NM),NESTJ(NM),& ! NUMC(NM),IRATIO(NM),NRATIO(NM),IEXP !34 FORMAT(/1X,'DOMAIN ',I2,2X,'IX=',I3,2X,'JX=',I3,2X, & ! 'DS= ',F6.1,2X,'ICNS=',I4,2X,'JCNS=',I4,2X,& ! 'NUMC= ',I2,2X,'IRATIO= ',I2,2X,'NRATIO= ',I2,2X,& ! 'IEXP= ',I2) XDIS = 0.0 YDIS = 0.0 ND1 = NM ND2 = NMC 40 CONTINUE XDIS = (NESTI(ND1)-1)*DIS(ND2) + XDIS YDIS = (NESTJ(ND1)-1)*DIS(ND2) + YDIS IF (ND2 .GT. 1) THEN ND1 = ND2 ND2 = NUMC(ND2) GO TO 40 ENDIF XSOUTH(NM) = XDIS/DIS(1) + FLOAT(IOFFST) + 1 XWEST(NM) = YDIS/DIS(1) + FLOAT(JOFFST) + 1 XNORTH(NM) = XSOUTH(NM) + FLOAT(NESTIX(NM)-1)*DIS(NM)/DIS(1) XEAST(NM) = XWEST(NM) + FLOAT(NESTJX(NM)-1)*DIS(NM)/DIS(1) ! PRINT 35 ! PRINT 36,XSOUTH(NM),XWEST(NM),XNORTH(NM),XEAST(NM) !35 FORMAT(2X,'COARSE MESH INDICES FOR THE 4 CORNER POINTS ARE') !36 FORMAT(2X,'SOUTH=',F6.1,3X,'WEST=',F6.1,3X,'NORTH =',F6.1,3X, & ! 'EAST = ',F6.1) 30 CONTINUE IF (MISMATCH.GT.0) THEN ! PRINT *, & ! 'TERRAIN STOP IN SUBROUTINE SETUP DUE TO INCORRECT NEST DOMAIN SET UP' ! STOP 1111 ENDIF ! Copy output XIM11 = XSOUTH XJM11 = XWEST #endif END SUBROUTINE SETUP