module sACOS_mod contains real(kind=SELECTED_REAL_KIND(13)) function sACOS(V) ! safer aCOS that handles near 1.0 and near -1.0 values implicit none real(kind=SELECTED_REAL_KIND(13)), intent(in) :: V real(kind=SELECTED_REAL_KIND(13)) :: av real(kind=SELECTED_REAL_KIND(13)), parameter :: near1=1.0001 av=V if(av>-near1 .and. av<-1.0) av=-1 if(av< near1 .and. av> 1.0) av=1 sACOS=acos(av) end function sACOS end module sACOS_mod SUBROUTINE EARTH_LATLON ( HLAT,HLON,VLAT,VLON, & !Earth lat,lon at H and V points DLMD1,DPHD1,WBD1,SBD1, & !input res,west & south boundaries, CENTRAL_LAT,CENTRAL_LON, & ! central lat,lon, all in degrees, IDS,IDE,JDS,JDE,KDS,KDE, & IMS,IME,JMS,JME,KMS,KME, & ITS,ITE,JTS,JTE,KTS,KTE ) !============================================================================ ! use sACOS_mod IMPLICIT NONE INTEGER, INTENT(IN ) :: IDS,IDE,JDS,JDE,KDS,KDE INTEGER, INTENT(IN ) :: IMS,IME,JMS,JME,KMS,KME INTEGER, INTENT(IN ) :: ITS,ITE,JTS,JTE,KTS,KTE REAL(4), INTENT(IN ) :: DLMD1,DPHD1,WBD1,SBD1 REAL(4), INTENT(IN ) :: CENTRAL_LAT,CENTRAL_LON REAL(4), DIMENSION(IMS:IME,JMS:JME), INTENT(OUT) :: HLAT,HLON,VLAT,VLON ! local INTEGER,PARAMETER :: KNUM=SELECTED_REAL_KIND(13) INTEGER :: I,J REAL(KIND=KNUM) :: WB,SB,DLM,DPH,TPH0,STPH0,CTPH0 REAL(KIND=KNUM) :: TDLM,TDPH,TLMH,TLMV,TLMH0,TLMV0,TPHH,TPHV,DTR REAL(KIND=KNUM) :: STPH,CTPH,STPV,CTPV,PI_2 REAL(KIND=KNUM) :: SPHH,CLMH,FACTH,SPHV,CLMV,FACTV REAL(KIND=KNUM), DIMENSION(IMS:IME,JMS:JME) :: GLATH,GLONH,GLATV,GLONV !------------------------------------------------------------------------- ! PI_2 = ACOS(0.) DTR = PI_2/90. WB = WBD1 * DTR ! WB: western boundary in radians SB = SBD1 * DTR ! SB: southern boundary in radians DLM = DLMD1 * DTR ! DLM: dlamda in radians DPH = DPHD1 * DTR ! DPH: dphi in radians TDLM = DLM + DLM ! TDLM: 2.0*dlamda TDPH = DPH + DPH ! TDPH: 2.0*DPH ! For earth lat lon only TPH0 = CENTRAL_LAT*DTR ! TPH0: central lat in radians STPH0 = SIN(TPH0) CTPH0 = COS(TPH0) WRITE(0,*) 'WB,SB,DLM,DPH,DTR: ',WBD1,SBD1,DLM,DPH,DTR WRITE(0,*) 'IMS,IME,JMS,JME,KMS,KME',IMS,IME,JMS,JME,KMS,KME WRITE(0,*) 'IDS,IDE,JDS,JDE,KDS,KDE',IDS,IDE,JDS,JDE,KDS,KDE WRITE(0,*) 'ITS,ITE,JTS,JTE,KTS,KTE',ITS,ITE,JTS,JTE,KTS,KTE ! .H DO J = JTS,MIN(JTE,JDE-1) ! H./ This loop takes care of zig-zag ! ! \.H starting points along j TLMH0 = WB - TDLM + MOD(J+1,2) * DLM ! ./ TLMH (rotated lats at H points) TLMV0 = WB - TDLM + MOD(J,2) * DLM ! H (//ly for V points) TPHH = SB + (J-1)*DPH ! TPHH (rotated lons at H points) are simple trans. TPHV = TPHH ! TPHV (rotated lons at V points) are simple trans. STPH = SIN(TPHH) CTPH = COS(TPHH) STPV = SIN(TPHV) CTPV = COS(TPHV) ! .H DO I = ITS,MIN(ITE,IDE-1) ! / TLMH = TLMH0 + I*TDLM ! \.H .U .H ! !H./ ----><---- SPHH = CTPH0 * STPH + STPH0 * CTPH * COS(TLMH) ! DLM + DLM GLATH(I,J)=ASIN(SPHH) ! GLATH: Earth Lat in radians CLMH = CTPH*COS(TLMH)/(COS(GLATH(I,J))*CTPH0) & - TAN(GLATH(I,J))*TAN(TPH0) IF(CLMH .GT. 1.) CLMH = 1.0 IF(CLMH .LT. -1.) CLMH = -1.0 FACTH = 1. IF(TLMH .GT. 0.) FACTH = -1. GLONH(I,J) = -CENTRAL_LON*DTR + FACTH*ACOS(CLMH) ENDDO DO I = ITS,MIN(ITE,IDE-1) TLMV = TLMV0 + I*TDLM SPHV = CTPH0 * STPV + STPH0 * CTPV * COS(TLMV) GLATV(I,J) = ASIN(SPHV) CLMV = CTPV*COS(TLMV)/(COS(GLATV(I,J))*CTPH0) & - TAN(GLATV(I,J))*TAN(TPH0) IF(CLMV .GT. 1.) CLMV = 1. IF(CLMV .LT. -1.) CLMV = -1. FACTV = 1. IF(TLMV .GT. 0.) FACTV = -1. GLONV(I,J) = -CENTRAL_LON*DTR + FACTV*ACOS(CLMV) ENDDO ENDDO ! Conversion to degrees (may not be required, eventually) DO J = JTS, MIN(JTE,JDE-1) DO I = ITS, MIN(ITE,IDE-1) HLAT(I,J) = GLATH(I,J) / DTR HLON(I,J)= -GLONH(I,J)/DTR IF(HLON(I,J) .GT. 180.) HLON(I,J) = HLON(I,J) - 360. IF(HLON(I,J) .LT. -180.) HLON(I,J) = HLON(I,J) + 360. VLAT(I,J) = GLATV(I,J) / DTR VLON(I,J) = -GLONV(I,J) / DTR IF(VLON(I,J) .GT. 180.) VLON(I,J) = VLON(I,J) - 360. IF(VLON(I,J) .LT. -180.) VLON(I,J) = VLON(I,J) + 360. ENDDO ENDDO END SUBROUTINE EARTH_LATLON SUBROUTINE G2T2H_egrid( IIH,JJH, & ! output grid index and weights HBWGT1,HBWGT2, & ! output weights in terms of parent grid HBWGT3,HBWGT4, & CLAT1, CLON1, & ! input source domain center CURRENT_DOMAIN_ID, & ! input source domain ID DX,DY, & ! input source domain resolution HLAT1,HLON1, & ! input source point regular lat and lon HLAT2,HLON2, & ! input target point regular lat and lon IDS,IDE,JDS,JDE, & ! input source outmost domain dimensions ISS,ISE,JSS,JSE, & ! source dimensions ITS,ITE,JTS,JTE ) ! target dimensions ! !*** XUEJIN ZHANG, Qingfu Liu --- Initial version (07/14/2011) !Function: Bilnear interpolation weight and indexing for E-grid ! from (HLAT1,HLON1) to (HLAT2,HLON2) ! !************************************************************* ! !*** FIGURE 1. ARRANGEMENT OF 4 VERTEXES FROM SOURCE DOMAIN !*** !*** 4 !*** !*** h !*** 1 2 !*** !*** !*** 3 ! !points 1,2,3,4 are source points !point h is target point !************************************************************* ! !Variables: !IIH---------I index in term of 4 source points (output) !JJH---------J index in term of 4 source points (output) !HBWGT1------weight in term of source point 1 (output) !HBWGT2------weight in term of source point 2 (output) !HBWGT3------weight in term of source point 3 (output) !HBWGT4------weight in term of source point 4 (output) !HLAT1-------source regular latitude (input) !HLON1-------source regular longitude (input) !HLAT2-------target regular latitude (input) !HLON2-------target regular longitude (input) !CLAT1-------source domain center regular latitude (input) !CLON1-------source domain center regular longitude (input) !CLAT2-------target domain center regular latitude (input) !CLON2-------target domain center regular longitude (input) !MAX_DOM-----source maximum domain number (input) !CURRENT_DOMAIN_ID--source domain ID (input) !I_PARENT_START-----source nest start I in parent domain (input) !J_PARENT_START-----source nest start I in parent domain (input) !DX,DY-------source outmost domain resolution (input) !RATIO-------source parent-nest domain ratio (always set = 3) !PARENT_ID---source PARENT ID (input) !GRID_ID-----source GRID ID (input) ! !IDS---------target outmost domain I-dimension start (input) !IDE---------target outmost domain I-dimension end (input) !JDS---------target outmost domain J-dimension start (input) !JDE---------target outmost domain J-dimension end (input) !ISS---------source array I-dimension start (input) !ISE---------source array I-dimension end (input) !JSS---------source array J-dimension start (input) !JSE---------source array J-dimension end (input) !ITS---------target array I-dimension start (input) !ITE---------target array I-dimension end (input) !JTS---------target array J-dimension start (input) !JTE---------target array J-dimension end (input) ! !************************************************************* ! DECLARE VARIABLES IMPLICIT NONE integer :: IDS,IDE,JDS,JDE integer :: ISS,ISE,JSS,JSE integer :: ITS,ITE,JTS,JTE integer, dimension(ITS:ITE,JTS:JTE) :: IIH,JJH real, dimension(ITS:ITE,JTS:JTE) :: HBWGT1,HBWGT2,HBWGT3,HBWGT4 real :: CLAT1, CLON1 real :: CLAT2, CLON2 integer :: CURRENT_DOMAIN_ID real :: DX,DY real, dimension(ISS:ISE,JSS:JSE) :: HLAT1,HLON1 real, dimension(ITS:ITE,JTS:JTE) :: HLAT2,HLON2 ! ! Local variables integer :: I,J,K,N,II,JJ real :: wb, sb, wb0, sb0 real, dimension(ISS:ISE,JSS:JSE) :: HLAT_R,HLON_R real :: slat, slon, rlat, rlon, olat, olon real, dimension(ISS:ISE,JSS:JSE) :: HLAT1_R,HLON1_R real, dimension(ITS:ITE,JTS:JTE) :: HLAT2_R,HLON2_R real, dimension(ITS:ITE,JTS:JTE) :: HLAT2_RT,HLON2_RT integer, dimension(ITS:ITE,JTS:JTE) :: IIH_RT,JJH_RT integer :: IHH_RT,JHH_RT real :: hlat2_tmp, hlon2_tmp integer :: iih_t,jjh_t,ihh_te,jhh_te integer :: icase real :: pi,aa,bb,alpha,beta real :: res,dis,dis1,dis2 real :: cint1,cint2,cint3,cint4,cint5 real :: check_sum !****************************************** ! I. Create source domain coordinator * !****************************************** ! !*************************************************************** ! I.1 Calculate west and south boundary of the source domain * !*************************************************************** print*,'inside g2t2h_egrid' print*,'dx,dy',dx,dy print *, 'ids,ide', ids,ide pi=atan(1.0)*4.0 res=sqrt(dx*dx+dy*dy) cint1=-pi/4. cint2=pi/4. cint3=3.*pi/4. cint4=5.*pi/4. cint5=7.*pi/4.+0.0001 !********************************************************************************* ! I.2 Define hlat and hlon of the source domain on rotated lat-lon coordinator * !********************************************************************************* olat = clat1; olon = clon1 call Convert_LatLon(1, (ISE-ISS+1), (JSE-JSS+1), hlon1, hlat1, hlon1_r, hlat1_r, olon, olat) ! Check hlat_r, hlon_r ! print *, 'hlon_r = ', hlon1_r(iss,jss), hlon1_r(ise,jse) ! print *, 'hlat_r = ', hlat1_r(iss,jss), hlat1_r(ise,jse) rlat = hlat1_r(iss,jss); rlon = hlon1_r(iss,jss) call Convert_LatLon(2, 1, 1, slon, slat, rlon, rlat, olon, olat) ! print *, 'Check the overlap between read-in (hlat1,hlon1) and (slat, slon) at SW point of source domain' ! print *, 'Defined by g2t2H-egrid, SW hlat1 = ', slat, 'SW hlon1 = ', slon ! print *, 'Read-in from wrfout_d*, SW hlat1 = ', hlat1(iss,jss), 'SW hlon1 = ', hlon1(iss,jss) rlat = hlat1_r(ise,jse); rlon = hlon1_r(ise,jse) call Convert_LatLon(2, 1, 1, slon, slat, rlon, rlat, olon, olat) ! print *, 'Check the overlap between read-in (hlat1,hlon1) and (slat, slon) at NE point of source domain' ! print *, 'Defined by g2t2H-egrid, NE hlat1 = ', slat, 'NE hlon1 = ', slon ! print *, 'Read-in from wrfout_d*, NE hlat1 = ', hlat1(ise,jse), 'NE hlon1 = ', hlon1(ise,jse) !********************************************************************************** ! II. Convert hlat and hlon of the target domain onto source domain * ! rotated lat-lon coordinator * !********************************************************************************** call Convert_LatLon(1, (ITE-ITS+1), (JTE-JTS+1), hlon2, hlat2, hlon2_r, hlat2_r, olon, olat) ! print*,'hlon2_r(1,1),hlat2_r(1,1)=',hlon2_r(1,1),hlat2_r(1,1) ! print*,'hlon1_r(1,1),hlat1_r(1,1)=',hlon1_r(1,1),hlat1_r(1,1) ! print*,'hlon2_r(ite,jte),hlat2_r(ite,jte)=',hlon2_r(ite,jte),hlat2_r(ite,jte) ! print*,'hlon1_r(ise,jse),hlat1_r(ise,jse)=',hlon1_r(ise,jse),hlat1_r(ise,jse) !********************************************************************************** ! III. Find index of (hlat2, hlon2) of the target domain on source domain * !********************************************************************************** ! !Function: Find the VERTEX 1 in FIGURE 1 ! ii=0 ! ! Counting the number of points within the source domain ! Purpose: eliminate loops outside of the source domain to reduce the computing time ! iih_rt = 0; jjh_rt = 0 if(current_domain_id.eq.1)then do j = jts,jte do i = its,ite hlat2_rt(i,j) = hlat2_r(i,j) hlon2_rt(i,j) = hlon2_r(i,j) iih_rt(i,j) = i jjh_rt(i,j) = j ii = ii + 1 end do end do else do j = jts,jte do i = its,ite if(hlat2_r(i,j).ge.hlat1_r(iss,jss) .and. & hlat2_r(i,j).lt.hlat1_r(ise,jse)) then if(hlon2_r(i,j).ge.hlon1_r(iss,jss) .and. & hlon2_r(i,j).lt.hlon1_r(ise,jse)) then hlat2_rt(i,j) = hlat2_r(i,j) hlon2_rt(i,j) = hlon2_r(i,j) iih_rt(i,j) = i jjh_rt(i,j) = j ii = ii + 1 endif endif enddo enddo end if print *,' total target points within source domain, ', ii ! ! Start inside loops ! ! write(81,'(4A5,4A12)') 'iih','jjh','iih_t','jjh_t','hlat2_tmp','hlon2_tmp','hlat1_r','hlon1_r' ! write(82,'(4A5,4A12)') 'iih','jjh','iih_t','jjh_t','hlat2_tmp','hlon2_tmp','hlat1_r','hlon1_r' ! write(83,'(4A5,4A12)') 'iih','jjh','iih_t','jjh_t','hlat2_tmp','hlon2_tmp','hlat1_r','hlon1_r' ! write(84,'(4A5,4A12)') 'iih','jjh','iih_t','jjh_t','hlat2_tmp','hlon2_tmp','hlat1_r','hlon1_r' ! write(85,'(4A5,4A12)') 'iih','jjh','iih_t','jjh_t','hlat2_tmp','hlon2_tmp','hlat1_r','hlon1_r' ! write(86,'(4A5,4A12)') 'iih','jjh','iih_t','jjh_t','hlat2_tmp','hlon2_tmp','hlat1_r','hlon1_r' do i = its, ite do j = jts, jte ! if(i.eq.75.and.j.eq.459)print *, 'qq test=',iih_rt(i,j),jjh_rt(i,j),hlat2_rt(i,j),hlon2_rt(i,j) if(iih_rt(i,j).gt.0.and.jjh_rt(i,j).gt.0) then hlat2_tmp = hlat2_rt(iih_rt(i,j),jjh_rt(i,j)) hlon2_tmp = hlon2_rt(iih_rt(i,j),jjh_rt(i,j)) iih_t = int((hlon2_tmp-hlon1_r(iss,jss))/(2*dx))+1 jjh_t = 2*int((hlat2_tmp-hlat1_r(iss,jss))/(2*dy))+1 if(iih_t.ge.iss.and.iih_t.lt.ise.and.jjh_t.gt.jss.and.jjh_t.lt.jse)then aa=hlon2_tmp-hlon1_r(iih_t,jjh_t+1) bb=hlat2_tmp-hlat1_r(iih_t,jjh_t+1) dis=sqrt(aa*aa+bb*bb)/res if(dis.lt.0.005)then iih(i,j)=iih_t jjh(i,j)=jjh_t+1 hbwgt1(i,j) = 1.0 hbwgt2(i,j) = 0.0 hbwgt3(i,j) = 0.0 hbwgt4(i,j) = 0.0 else alpha=atan2(bb,aa) if(alpha.lt.(-pi/4))alpha=alpha+2.*pi if(alpha.ge.cint1.and.alpha.lt.cint2)then iih(i,j)=iih_t jjh(i,j)=jjh_t+1 beta=cint2-alpha dis1=dis*sin(beta) dis2=dis*cos(beta) if(dis1.lt.(-0.001).or.dis1.gt.1.001.or. & dis2.lt.(-0.001).or.dis2.gt.1.001)then print*,'qliu test wrong1=',i,j,dis,dis1,dis2,beta stop end if hbwgt1(i,j) = (1.-dis1)*(1.-dis2) hbwgt2(i,j) = dis1*dis2 hbwgt3(i,j) = dis1*(1.-dis2) hbwgt4(i,j) = (1.-dis1)*dis2 else if(alpha.ge.cint2.and.alpha.lt.cint3)then iih(i,j)=iih_t jjh(i,j)=jjh_t+2 beta=cint3-alpha dis1=dis*sin(beta) dis2=dis*cos(beta) if(dis1.lt.(-0.001).or.dis1.gt.1.001.or. & dis2.lt.(-0.001).or.dis2.gt.1.001)then print*,'qliu test wrong2=',i,j,dis,dis1,dis2,beta stop end if hbwgt1(i,j) = (1.-dis1)*dis2 hbwgt2(i,j) = dis1*(1.-dis2) hbwgt3(i,j) = (1.-dis1)*(1.-dis2) hbwgt4(i,j) = dis1*dis2 else if(alpha.ge.cint3.and.alpha.lt.cint4)then iih(i,j)=iih_t-1 jjh(i,j)=jjh_t+1 beta=cint4-alpha dis1=dis*sin(beta) dis2=dis*cos(beta) if(dis1.lt.(-0.0015).or.dis1.gt.1.0015.or. & dis2.lt.(-0.0015).or.dis2.gt.1.0015)then print*,'qliu test wrong3=',i,j,dis,dis1,dis2,beta stop end if hbwgt1(i,j) = dis1*dis2 hbwgt2(i,j) = (1.-dis1)*(1.-dis2) hbwgt3(i,j) = (1.-dis1)*dis2 hbwgt4(i,j) = dis1*(1.-dis2) else if(alpha.ge.cint4.and.alpha.le.cint5)then iih(i,j)=iih_t jjh(i,j)=jjh_t beta=cint5-alpha dis1=dis*sin(beta) dis2=dis*cos(beta) if(dis1.lt.(-0.001).or.dis1.gt.1.001.or. & dis2.lt.(-0.001).or.dis2.gt.1.001)then print*,'qliu test wrong4=',i,j,dis,dis1,dis2,beta stop end if hbwgt1(i,j) = dis1*(1.-dis2) hbwgt2(i,j) = (1.-dis1)*dis2 hbwgt3(i,j) = dis1*dis2 hbwgt4(i,j) = (1.-dis1)*(1.-dis2) end if endif endif endif end do end do ! do i = its,ite ! do j = jts,jte ! check_sum=abs(hbwgt1(i,j))+abs(hbwgt2(i,j))+abs(hbwgt3(i,j))+abs(hbwgt4(i,j)) ! if(check_sum.gt.1.01)then ! print*, 'qingfu test test',icase ! print*,iih(i,j),jjh(i,j),hbwgt1(i,j),hbwgt2(i,j),hbwgt3(i,j),hbwgt4(i,j) ! print *, hlon2_tmp,hlon1_r(iih_t,jjh_t),hlon1_r(iih_t+1,jjh_t) ! print *, hlon2_tmp,hlon1_r(iih_t,jjh_t+2),hlon1_r(iih_t+1,jjh_t+2) ! print *, hlat2_tmp,hlat1_r(iih_t,jjh_t),hlat1_r(iih_t+1,jjh_t) ! print *, hlat2_tmp,hlat1_r(iih_t,jjh_t+2),hlat1_r(iih_t+1,jjh_t+2) ! stop ! end if ! write(91,'(4I6,4F10.6)') iih(i,j),jjh(i,j),i,j,hbwgt1(i,j),hbwgt2(i,j),hbwgt3(i,j),hbwgt4(i,j) ! end do ! end do !********************************************************************************** ! IV. Find weight of (hlat2, hlon2) of the target domain on source domain * !********************************************************************************** ! End of g2t2h_egrid end subroutine g2t2h_egrid SUBROUTINE G2T2V_egrid( IIV,JJV, & ! output grid index and weights VBWGT1,VBWGT2, & ! output weights in terms of parent grid VBWGT3,VBWGT4, & CLAT1, CLON1, & ! input source domain center(-6h) CURRENT_DOMAIN_ID, & ! input source domain ID DX,DY, & ! input source outmost domain resolution VLAT1,VLON1, & ! input source point regular lat and lon VLAT2,VLON2, & ! input target point regular lat and lon IDS,IDE,JDS,JDE, & ! input target outmost domain dimensions ISS,ISE,JSS,JSE, & ! source dimensions ITS,ITE,JTS,JTE ) ! target dimensions ! !*** XUEJIN ZHANG, Qingfu Liu --- Initial version (07/14/2011) !Function: Bilnear interpolation weight and indexing for E-grid ! from (VLAT1,HLON1) to (VLAT2,HLON2) ! !************************************************************* ! !*** FIGURE 1. ARRANGEMENT OF 4 VERTEXES FROM SOURCE DOMAIN !*** !*** 4 !*** !*** v !*** 1 2 !*** !*** !*** 3 ! !points 1,2,3,4 are source points !point h is target point !************************************************************* ! !Variables: !IIV---------I index in term of 4 source points (output) !JJV---------J index in term of 4 source points (output) !VBWGT1------weight in term of source point 1 (output) !VBWGT2------weight in term of source point 2 (output) !VBWGT3------weight in term of source point 3 (output) !VBWGT4------weight in term of source point 4 (output) !VLAT1-------source regular latitude (input) !VLON1-------source regular longitude (input) !VLAT2-------target regular latitude (input) !VLON2-------target regular longitude (input) !CLAT1-------source domain center regular latitude (input) !CLON1-------source domain center regular longitude (input) !CLAT2-------target domain center regular latitude (input) !CLON2-------target domain center regular longitude (input) !MAX_DOM-----source maximum domain number (input) !CURRENT_DOMAIN_ID--source domain ID (input) !I_PARENT_START-----source nest start I in parent domain (input) !J_PARENT_START-----source nest start I in parent domain (input) !DX,DY-------source outmost domain resolution (input) !RATIO-------source parent-nest domain ratio (always set = 3) !PARENT_ID---source PARENT ID (input) !GRID_ID-----source GRID ID (input) ! !IDS---------target outmost domain I-dimension start (input) !IDE---------target outmost domain I-dimension end (input) !JDS---------target outmost domain J-dimension start (input) !JDE---------target outmost domain J-dimension end (input) !ISS---------source array I-dimension start (input) !ISE---------source array I-dimension end (input) !JSS---------source array J-dimension start (input) !JSE---------source array J-dimension end (input) !ITS---------target array I-dimension start (input) !ITE---------target array I-dimension end (input) !JTS---------target array J-dimension start (input) !JTE---------target array J-dimension end (input) ! !************************************************************* ! DECLARE VARIABLES IMPLICIT NONE integer :: IDS,IDE,JDS,JDE integer :: ISS,ISE,JSS,JSE integer :: ITS,ITE,JTS,JTE integer, dimension(ITS:ITE,JTS:JTE) :: IIV,JJV real, dimension(ITS:ITE,JTS:JTE) :: VBWGT1,VBWGT2,VBWGT3,VBWGT4 real :: CLAT1, CLON1 real :: CLAT2, CLON2 integer :: CURRENT_DOMAIN_ID real :: DX,DY real, dimension(ISS:ISE,JSS:JSE) :: VLAT1,VLON1 real, dimension(ITS:ITE,JTS:JTE) :: VLAT2,VLON2 ! ! Local variables integer :: I,J,K,N,ii,jj real :: wb, sb, wb0, sb0 real, dimension(ISS:ISE,JSS:JSE) :: VLAT_R,VLON_R real :: slat, slon, rlat, rlon, olat, olon real, dimension(ISS:ISE,JSS:JSE) :: VLAT1_R,VLON1_R real, dimension(ITS:ITE,JTS:JTE) :: VLAT2_R,VLON2_R real, dimension(ITS:ITE,JTS:JTE) :: VLAT2_RT,VLON2_RT integer, dimension(ITS:ITE,JTS:JTE) :: IIV_RT,JJV_RT integer :: IVV_RT,JVV_RT real :: vlat2_tmp, vlon2_tmp integer :: iiv_t,jjv_t,ivv_te,jvv_te integer :: icase real :: pi,aa,bb,alpha,beta real :: res,dis,dis1,dis2 real :: cint1,cint2,cint3,cint4,cint5 real :: check_sum !****************************************** ! I. Create source domain coordinator * !****************************************** ! !*************************************************************** ! I.1 Calculate west and south boundary of the source domain * !*************************************************************** print*,'inside g2t2v_egrid' print*,'dx,dy',dx,dy print *, 'ids,ide', ids,ide pi=atan(1.0)*4.0 res=sqrt(dx*dx+dy*dy) cint1=-pi/4. cint2=pi/4. cint3=3.*pi/4. cint4=5.*pi/4. cint5=7.*pi/4.+0.0001 !********************************************************************************* ! I.2 Define vlat and vlon of the source domain on rotated lat-lon coordinator * !********************************************************************************* olat = clat1; olon = clon1 call Convert_LatLon(1, (ISE-ISS+1), (JSE-JSS+1), vlon1, vlat1, vlon1_r, vlat1_r, olon, olat) ! Check vlat_r, vlon_r ! print *, 'vlon_r = ', vlon1_r(iss,jss), vlon1_r(ise,jse) ! print *, 'vlat_r = ', vlat1_r(iss,jss), vlat1_r(ise,jse) rlat = vlat1_r(iss,jss); rlon = vlon1_r(iss,jss) call Convert_LatLon(2, 1, 1, slon, slat, rlon, rlat, olon, olat) ! print *, 'Check the overlap between read-in (vlat1,vlon1) and (slat, slon) at SW point of source domain' ! print *, 'Defined by g2t2V-egrid, SW vlat1 = ', slat, 'SW vlon1 = ', slon ! print *, 'Read-in from wrfout_d*, SW vlat1 = ', vlat1(iss,jss), 'SW vlon1 = ', vlon1(iss,jss) rlat = vlat1_r(ise,jse); rlon = vlon1_r(ise,jse) call Convert_LatLon(2, 1, 1, slon, slat, rlon, rlat, olon, olat) ! print *, 'Check the overlap between read-in (vlat1,vlon1) and (slat, slon) at NE point of source domain' ! print *, 'Defined by g2t2V-egrid, NE vlat1 = ', slat, 'NE vlon1 = ', slon ! print *, 'Read-in from wrfout_d*, NE vlat1 = ', vlat1(ise,jse), 'NE vlon1 = ', vlon1(ise,jse) !********************************************************************************** ! II. Convert vlat and vlon of the target domain onto source domain * ! rotated lat-lon coordinator * !********************************************************************************** call Convert_LatLon(1, (ITE-ITS+1), (JTE-JTS+1), vlon2, vlat2, vlon2_r, vlat2_r, olon, olat) ! print*,'vlon2_r(1,1),vlat2_r(1,1)=',vlon2_r(1,1),vlat2_r(1,1) ! print*,'vlon1_r(1,1),vlat1_r(1,1)=',vlon1_r(1,1),vlat1_r(1,1) ! print*,'vlon2_r(ite,jte),vlat2_r(ite,jte)=',vlon2_r(ite,jte),vlat2_r(ite,jte) ! print*,'vlon1_r(ise,jse),vlat1_r(ise,jse)=',vlon1_r(ise,jse),vlat1_r(ise,jse) !********************************************************************************** ! III. Find index of (hlat2, hlon2) of the target domain on source domain * !********************************************************************************** ! !Function: Find the VERTEX 1 in FIGURE 1 ! ii=0 ! ! Counting the number of points within the source domain ! Purpose: eliminate loops outside of the source domain to reduce the computing time ! iiv_rt = 0; jjv_rt = 0 if(current_domain_id.eq.1)then do j = jts,jte do i = its,ite vlat2_rt(i,j) = vlat2_r(i,j) vlon2_rt(i,j) = vlon2_r(i,j) iiv_rt(i,j) = i jjv_rt(i,j) = j ii = ii + 1 end do end do else ii=0 do j = jts,jte do i = its,ite if(vlat2_r(i,j).ge.vlat1_r(iss,jss) .and. & vlat2_r(i,j).lt.vlat1_r(ise,jse)) then if(vlon2_r(i,j).ge.vlon1_r(iss,jss) .and. & vlon2_r(i,j).lt.vlon1_r(ise,jse)) then vlat2_rt(i,j) = vlat2_r(i,j) vlon2_rt(i,j) = vlon2_r(i,j) iiv_rt(i,j) = i jjv_rt(i,j) = j ii = ii + 1 endif endif enddo enddo end if print *,' total target points within source domain, ', ii do i = its, ite do j = jts, jte if(iiv_rt(i,j).gt.0.and.jjv_rt(i,j).gt.0) then vlat2_tmp = vlat2_rt(iiv_rt(i,j),jjv_rt(i,j)) vlon2_tmp = vlon2_rt(iiv_rt(i,j),jjv_rt(i,j)) iiv_t = int((vlon2_tmp-vlon1_r(iss,jss))/(2*dx))+1 jjv_t = 2*int((vlat2_tmp-vlat1_r(iss,jss))/(2*dy))+1 if(iiv_t.ge.iss.and.iiv_t.lt.ise.and.jjv_t.gt.jss.and.jjv_t.lt.jse)then aa=vlon2_tmp-vlon1_r(iiv_t+1,jjv_t+1) bb=vlat2_tmp-vlat1_r(iiv_t+1,jjv_t+1) dis=sqrt(aa*aa+bb*bb)/res if(dis.lt.0.005)then iiv(i,j)=iiv_t+1 jjv(i,j)=jjv_t+1 vbwgt1(i,j) = 1.0 vbwgt2(i,j) = 0.0 vbwgt3(i,j) = 0.0 vbwgt4(i,j) = 0.0 else alpha=atan2(bb,aa) if(alpha.lt.(-pi/4))alpha=alpha+2.*pi if(alpha.ge.cint1.and.alpha.lt.cint2)then iiv(i,j)=iiv_t+1 jjv(i,j)=jjv_t+1 beta=cint2-alpha dis1=dis*sin(beta) dis2=dis*cos(beta) if(dis1.lt.(-0.001).or.dis1.gt.1.001.or. & dis2.lt.(-0.001).or.dis2.gt.1.001)then print*,'qliu test wrong1=',i,j,dis,dis1,dis2,beta stop end if vbwgt1(i,j) = (1.-dis1)*(1.-dis2) vbwgt2(i,j) = dis1*dis2 vbwgt3(i,j) = dis1*(1.-dis2) vbwgt4(i,j) = (1.-dis1)*dis2 else if(alpha.ge.cint2.and.alpha.lt.cint3)then iiv(i,j)=iiv_t jjv(i,j)=jjv_t+2 beta=cint3-alpha dis1=dis*sin(beta) dis2=dis*cos(beta) if(dis1.lt.(-0.001).or.dis1.gt.1.001.or. & dis2.lt.(-0.001).or.dis2.gt.1.001)then print*,'qliu test wrong2=',i,j,dis,dis1,dis2,beta stop end if vbwgt1(i,j) = (1.-dis1)*dis2 vbwgt2(i,j) = dis1*(1.-dis2) vbwgt3(i,j) = (1.-dis1)*(1.-dis2) vbwgt4(i,j) = dis1*dis2 else if(alpha.ge.cint3.and.alpha.lt.cint4)then iiv(i,j)=iiv_t jjv(i,j)=jjv_t+1 beta=cint4-alpha dis1=dis*sin(beta) dis2=dis*cos(beta) if(dis1.lt.(-0.0015).or.dis1.gt.1.0015.or. & dis2.lt.(-0.0015).or.dis2.gt.1.0015)then print*,'qliu test wrong3=',i,j,dis,dis1,dis2,beta stop end if vbwgt1(i,j) = dis1*dis2 vbwgt2(i,j) = (1.-dis1)*(1.-dis2) vbwgt3(i,j) = (1.-dis1)*dis2 vbwgt4(i,j) = dis1*(1.-dis2) else if(alpha.ge.cint4.and.alpha.le.cint5)then iiv(i,j)=iiv_t jjv(i,j)=jjv_t beta=cint5-alpha dis1=dis*sin(beta) dis2=dis*cos(beta) if(dis1.lt.(-0.001).or.dis1.gt.1.001.or. & dis2.lt.(-0.001).or.dis2.gt.1.001)then print*,'qliu test wrong4=',i,j,dis,dis1,dis2,beta stop end if vbwgt1(i,j) = dis1*(1.-dis2) vbwgt2(i,j) = (1.-dis1)*dis2 vbwgt3(i,j) = dis1*dis2 vbwgt4(i,j) = (1.-dis1)*(1.-dis2) end if endif endif endif end do end do ! do i = its,ite ! do j = jts,jte ! check_sum=abs(vbwgt1(i,j))+abs(vbwgt2(i,j))+abs(vbwgt3(i,j))+abs(vbwgt4(i,j)) ! if(check_sum.gt.1.01)then ! print*, 'qingfu test test' ! print*,iiv(i,j),jjv(i,j),vbwgt1(i,j),vbwgt2(i,j),vbwgt3(i,j),vbwgt4(i,j) ! print *, vlon2_tmp,vlon1_r(iiv_t,jjv_t),vlon1_r(iiv_t+1,jjv_t) ! print *, vlon2_tmp,vlon1_r(iiv_t,jjv_t+2),vlon1_r(iiv_t+1,jjv_t+2) ! print *, vlat2_tmp,vlat1_r(iiv_t,jjv_t),vlat1_r(iiv_t+1,jjv_t) ! print *, vlat2_tmp,vlat1_r(iiv_t,jjv_t+2),vlat1_r(iiv_t+1,jjv_t+2) ! stop ! end if ! write(93,'(4I6,4F10.6)') iiv(i,j),jjv(i,j),i,j,vbwgt1(i,j),vbwgt2(i,j),vbwgt3(i,j),vbwgt4(i,j) ! end do ! end do !********************************************************************************** ! IV. Find weight of (vlat2, vlon2) of the target domain on source domain * !********************************************************************************** ! End of g2t2v_egrid end subroutine g2t2v_egrid !============================================================================== SUBROUTINE Convert_LatLon(mode, im, jm, slon, slat, rlon, rlat, olon, olat) !****************************************************************************** ! ! Convert between standard & rotated latlon coordinates ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! This subroutine converts latlon coordinates between two coordinate systems, ! the standard geographic system {S} & the rotated spherical coord system {R} ! obtained by rotating the origin of {S} to (olat,olon) as the new origin: ! (0,0) -> (olat,olon) = (0,0)'. The poles are not allowed for conversion. ! ! mode = 1 : convert standard coord (alat,alon) to rotated coord (blat,blon) ! mode = 2 : convert rotated coord (blat,blon) to standard coord (alat,alon) ! Author: Kao-San Yeh (AOML/HRD) !______________________________________________________________________________ Implicit None ! -------------------------------------------------------------------- [input] Integer mode, im, jm !* mode= 1/2: Std -> Rot / Rot -> Std coord Real(4) olon, olat !* origin of rotated grid in standard coord ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - [in/out] Real(4) slon(im,jm), slat(im,jm) !* standard latlon coordinates (deg) Real(4) rlon(im,jm), rlat(im,jm) !* rotated latlon coordinates (deg) ! -------------------------------------------------------------------- [local] Real(8) alon(im,jm), alat(im,jm) !* standard latlon coordinates (deg) Real(8) blon(im,jm), blat(im,jm) !* rotated latlon coordinates (deg) Real(8) pi_cnst, deg2rad, eps1, sqrt2_, xcos, xsin, ycos, ysin Real(8) olamda, otheta, alamda, atheta, blamda, btheta Real(8) aaa, bbb, ccc, ddd, eee, xxx, yyy, zzz Integer i, j ! * * * * * * * * * * pi_cnst = DAcos(-1.D0) ; deg2rad = DAcos(-1.D0)/180. eps1 = 1. - 1.D-6 ; sqrt2_ = 1./sqrt(2.) !* Enhance precision for trigonometric functions !* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ alon(:,:) = slon(:,:) ; alat(:,:) = slat(:,:) blon(:,:) = rlon(:,:) ; blat(:,:) = rlat(:,:) olamda = olon*deg2rad ; otheta = olat*deg2rad aaa = Dcos(otheta) ; bbb = Dsin(otheta) !* Convert standard Earth latlon to rotated Model latlon !* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IF ( mode == 1 ) THEN !* - - - - - - - - - - - - - - - - - - - - DO j = 1, jm DO i = 1, im yyy = alat(i,j)*deg2rad xxx = alon(i,j)*deg2rad - olamda xcos = Dcos(xxx) ; xsin = Dsin(xxx) ycos = Dcos(yyy) ; ysin = Dsin(yyy) ! - - - - - - - - - - - - - - - - - - zzz = aaa*ysin - bbb*ycos*xcos zzz = max(-eps1, min(zzz, eps1)) !* Poles undefined btheta = DAsin(zzz) ; eee = 1./Dcos(btheta) ccc = ( bbb*ysin + aaa*ycos*xcos )*eee ddd = ( ycos*xsin )*eee IF ( ccc > sqrt2_ ) THEN blamda = DAsin(ddd) ELSEIF ( ccc < -sqrt2_ ) THEN blamda = pi_cnst*Sign(1.D+0,ddd) - DAsin(ddd) ELSE ! ( |ccc| <= sqrt2_ ) ! blamda = DAcos(ccc)*Sign(1.D+0,ddd) ENDIF blat(i,j) = btheta/deg2rad blon(i,j) = blamda/deg2rad ENDDO ENDDO !* Convert rotated Model latlon to standard Earth latlon !* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ELSE !* ( mode == 2 ) - - - - - - - - - - - - - - - - - - - - - - DO j = 1, jm DO i = 1, im yyy = blat(i,j)*deg2rad xxx = blon(i,j)*deg2rad xcos = Dcos(xxx) ; xsin = Dsin(xxx) ycos = Dcos(yyy) ; ysin = Dsin(yyy) ! - - - - - - - - - - - - - - - - - - zzz = aaa*ysin + bbb*ycos*xcos zzz = max(-eps1, min(zzz, eps1)) !* Poles undefined atheta = DAsin(zzz) ; eee = 1./Dcos(atheta) ccc = (-bbb*ysin + aaa*ycos*xcos )*eee ddd = ( ycos*xsin )*eee IF ( ccc > sqrt2_ ) THEN alamda = DAsin(ddd) + olamda ELSEIF ( ccc < -sqrt2_ ) THEN alamda = pi_cnst*Sign(1.D+0,ddd) - DAsin(ddd) + olamda ELSE ! ( |ccc| <= sqrt2_ ) ! alamda = DAcos(ccc)*Sign(1.D+0,ddd) + olamda ENDIF alat(i,j) = atheta/deg2rad alon(i,j) = alamda/deg2rad ENDDO ENDDO ENDIF !* - - - - - - - - - - - - - - - - - - - - - - - - - - - - !* Recover single precision for output !* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ slon(:,:) = alon(:,:) ; slat(:,:) = alat(:,:) rlon(:,:) = blon(:,:) ; rlat(:,:) = blat(:,:) END SUBROUTINE Convert_LatLon