module mod_interp CONTAINS ! !*********************************************************************** ! SUBROUTINE INTP_INI( & &IM1,JM1,IM2,JM2,lon1,lat1,LON2,LAT2,mask1,MASK2, & &IG1,JG1,ierr, & &ia,ja,ib,jb,ic,jc,KM1,ibndr1,jbndr1, & &xdlon1,xdlat1,XDLON2,XDLAT2,COSANG,SINANG, & &nvrbs,iiunit) ! ! Calculate IG1( , ),JG1( , ). By definition, for each G2 sea ! gridpoint (i,j) (MASK2(i,j)=0), whose coordinates are ! (LON2(i,j),LAT2(i,j)), ! ( IG1(i,j),JG1(i,j) ) = (i1,j1) ! is the point on G1 grid for which (i,j) lies in the quadrangle ! whose vertices are *valid* (see definition of validity below) G1 ! gridpoints (i1,j1), (i1+1,j1), (i1+1,j1+1), (i1,j1+1) if these ! vertices are not all land points; if they are all land points, ! IG1(i,j) is set to 0 and no value is assigned to JG1(i,j). ! A G1 gridpoint (i,j) is considered *valid* if it is not flagged ! with mask1(i,j)=-1 value. (Such are gridpoints where no valid ! coordinate values are available.) ! If ia,ja,ib,jb,ic,jc are present, the 4 vertices are, ! instead of the above, (i1,j1), (ia(i1,j1),ja(i1,j1)), ! (ib(i1,j1),jb(i1,j1)), (ic(i1,j1),jc(i1,j1)). ! If there is no such a quadrangle, i.e. if a valid G2 point (i,j) is ! outside G1 domain, IG1(i,j) is set to -2; no value is assigned to ! JG1(i,j). If (i,j) is a G2 land point (MASK2(i,j).eq.1) lying inside ! G1 domain, IG1(i,j) is set to -1; no value is assigned to JG1(i,j). ! If (i,j) is not a valid G2 point (which is supposed to be flagged ! with MASK2(i,j)=-1 value) IG1(i,j) is set to -3; no value is ! assigned to JG1(i,j). ! ! Also, if present, calculate COSANG( , ), SINANG( , ). If ! for each G2 gridpoint (i,j) ANG is the angle between x direction ! (the direction of growing first index) on G1 grid and x direction ! at point (i,j) on G2 grid (the angle may have either sign and is ! measured from the latter direction to the former) then ! COSANG(i,j)=cos(ANG), SINANG(i,j)=sin(ANG) . ! implicit none integer,intent(in):: IM1,JM1,IM2,JM2 !<- dimensions for G1 and G2 grid arrays real,dimension(IM1,JM1),intent(in):: lon1,lat1 !<- arrays of longitudes and latitudes for G1 grid real,dimension(IM2,JM2),intent(in):: LON2,LAT2 !<- arrays of longitudes and latitudes for G2 grid integer,dimension(IM1,JM1),intent(in):: mask1 !<- sea/land mask for G1 grid, 0 - sea, 1 - land; ! -1 for not valid gridpoints (see above) integer,dimension(IM2,JM2),intent(in):: MASK2 !<- sea/land mask for G2 grid, 0 - sea, 1 - land; ! -1 for not valid gridpoints (see above) INTEGER,DIMENSION(IM2,JM2),INTENT(OUT):: IG1,JG1 !<- arrays that determine in which G1 gridcells G2 ! gridpoints lie, see above INTEGER,INTENT(INOUT):: ierr !<- out: 0 for normal return, nonzero otherwise; ! in: if data are read from disk (see iiu) and if the ! read input data fail to coincide with the ! corresponding arguments then: if ierr=0 on input the ! subroutine will return abnormally and if ierr!=0 on ! input the data will be calculated as if with iiu<0 integer,optional,dimension(IM1,JM1),intent(in):: ia,ja,ib,jb,ic,jc !<- determine the 4 vertices of each G1 gridcell, ! see above integer,optional,intent(in):: KM1 !<- length of G1 grid boundary (see ibndr1, jbndr1). ! If present (absent), both ibndr1,jbndr1 must also be ! present (absent) integer,optional,dimension(0:),intent(in):: ibndr1,jbndr1 !<- arrays representing i- and j-indices of those ! G1 gridpoints which belong to the G1 domain ! boundary; the representation is supposed to be ! continuous, i.e. subsequent array elements represent ! adjacent gridpoints; the boundary must be closed: ! ibndr1(0)=ibndr1(KM1), jbndr1(0)=jbndr1(KM1). ! If these arguments are absent, it is assumed ! that the G1 domain boundary consists of all those ! gridpoints (i,j) for which i=1 or i=IM1 or j=1 or ! j=JM1. real,optional,dimension(IM1,JM1),intent(in):: xdlon1,xdlat1 !<- latitudes and longitudes of x-direction vectors ! at each G1 gridpoints. The arguments must be either ! both present or both absent. If they are absent it ! is assumed that x-direction at each G1 gridpoint ! (i0,j0) is the direction of i increasing from i0 ! with j=j0=const. These are input data for ! computation of COSANG, SINANG (see above/below) real,optional,dimension(IM2,JM2),intent(in):: XDLON2,XDLAT2 !<- latitudes and longitudes of x-direction vectors ! at each G2 gridpoints. The arguments must be either ! both present or both absent. If they are absent it ! is assumed that x-direction at each G2 gridpoint ! (i0,j0) is the direction of i increasing from i0 ! with j=j0=const. These are input data for ! computation of COSANG, SINANG (see above/below) REAL,OPTIONAL,DIMENSION(IM2,JM2),INTENT(OUT):: COSANG,SINANG !<- cosine and sine of the angle between the ! x direction on G2 grid and x direction on G1 ! grid at each G2 gridpoint. The angle is measured ! from G2 x direction to G1 x direction. (For ! definition of x directions, see above.) integer,optional,intent(in):: nvrbs !<- verbosity level. With nvrbs=0 or if the argument ! is absent the subroutine is silent unless it ! returns abnormally. Positive value recommended for ! debugging/testing integer,optional,intent(in):: iiunit !<- if iiunit>0 the OUT arguments are not calculated ! but read from disk unit number iiunit, along with ! copies of IN arguments, in the order: ! IM1,JM1,IM2,JM2,lon1,lat1,LON2,LAT2,mask1,MASK2, ! IG1,JG1 (1 record to read) and, if present, ! COSANG, SINANG (next record to read). ! The read values for IN arguments must coincide with ! the corresponding values supplied through actual ! arguments; the reading is done without opening/ ! rewinding, from the current position; ! if iiunit=0 or is absent no disk exchange is ! undertaken; computations proceed as usual; ! if iiunit<0 the computations proceed as usual and ! the above arguments in the above order are written ! to disk unit number -iiunit, in the same way as they ! are read when iiunit>0. integer i,j,i1,j1,i2,j2,i2p,j2p,i0,j0,i1a,j1a,i1b,j1b,i1c,j1c, & &is,ie,js,je,isp,iep,jsp,jep,istep,nvrb, & &ios,ios1,ierr_in,ierr1,iiu,ivrbs,n,n1,k,k1,km,kmm,kb1,KM2,KM2M ! integer,dimension(:),allocatable:: ibnd1,jbnd1 integer,dimension(0:IM1*JM1):: ibnd1,jbnd1 integer,dimension(IM2*JM2):: IPEANO2,JPEANO2 integer,dimension(0:IM1,0:JM1):: msk1 integer,dimension(0:IM2,0:JM2):: MSK2 real x1,x1a,x1b,x1c,y1,y1a,y1b,y1c, & &r1r2,dlon2,dlat2,dlon1on2,dlat1on2 real,parameter:: pi=3.1415926535897932,rad=pi/180. real,dimension(IM1,JM1):: dlon1,dlat1 logical inG1,inG1cell,cross,stuck integer IM1r,JM1r,IM2r,JM2r real,dimension(IM1,JM1):: lon1r,lat1r real,dimension(IM2,JM2):: LON2r,LAT2r integer,dimension(IM1,JM1):: mask1r integer,dimension(IM2,JM2):: MASK2r ! ierr_in=ierr ierr=0 ! if (PRESENT(nvrbs)) then nvrb=nvrbs else nvrb=0 end if if (PRESENT(COSANG).neqv.PRESENT(SINANG)) then print*,'INTP_INI: **ERROR:** both or neither '// & & 'SINANG, COSANG must be present' ierr=1 RETURN end if if ( (PRESENT(ia).neqv.PRESENT(ja)) .or. & & (PRESENT(ia).neqv.PRESENT(ib)) .or. & & (PRESENT(ia).neqv.PRESENT(jb)) .or. & & (PRESENT(ia).neqv.PRESENT(ic)) .or. & & (PRESENT(ia).neqv.PRESENT(jc)) ) then print*,'INTP_INI: **ERROR** ia,ja,ib,jb,ic,jc must be all '// & & 'present or all not present' ierr=1 RETURN end if if ( (PRESENT(xdlon1).neqv.PRESENT(xdlat1)) ) then print*,'INTP_INI: **ERROR:** both or neither '// & & 'xdlon1, xdlat1 must be present' ierr=1 RETURN end if if ( (PRESENT(XDLON2).neqv.PRESENT(XDLAT2)) ) then print*,'INTP_INI: **ERROR:** both or neither '// & & 'XDLON2, XDLAT2 must be present' ierr=1 RETURN end if ! IF ( (PRESENT(KM1).neqv.PRESENT(ibndr1)) .or.& &(PRESENT(KM1).neqv.PRESENT(jbndr1)) ) THEN print*,'INTP_INI: **ERROR:** all or none of arguments '//& & 'KM1, ibndr1,jbndr1 must be present' ierr=1 RETURN END IF if (PRESENT(COSANG)) then ivrbs=10 if (PRESENT(xdlon1)) then dlon1=xdlon1 dlat1=xdlat1 else do j=1,JM1 do i=1,IM1 if (i.eq.IM1) then dlon1(i,j)=lon1(i,j)-lon1(i-1,j) dlat1(i,j)=lat1(i,j)-lat1(i-1,j) else if (i.eq.1) then dlon1(i,j)=lon1(i+1,j)-lon1(i,j) dlat1(i,j)=lat1(i+1,j)-lat1(i,j) else dlon1(i,j)=lon1(i+1,j)-lon1(i-1,j) dlat1(i,j)=lat1(i+1,j)-lat1(i-1,j) end if end do end do end if end if ! ! if (allocated(ibnd1)) then ! print*,'INTP_INI: ibnd1 previously allocated, deallocating...' ! deallocate(ibnd1) ! end if ! if (allocated(jbnd1)) then ! print*,'INTP_INI: jbnd1 previously allocated, deallocating...' ! deallocate(jbnd1) ! end if ! iiu=0 if (PRESENT(iiunit)) iiu=iiunit IF (iiu.gt.0) THEN read (iiu,iostat=ios) IM1r,JM1r,IM2r,JM2r,& & lon1r,lat1r,LON2r,LAT2r,mask1r,MASK2r,& & IG1,JG1 ios1=0 if (PRESENT(COSANG)) read (iiu,iostat=ios1) COSANG,SINANG if (ios.ne.0 .or. ios1.ne.0) then print*,'INTP_INI: ***ERROR*** reading INTP_INI data ',& & ios,ios1 if (ierr_in.eq.0) then ierr=2 RETURN else print*,'INTP_INI: output data will be '//& & 'computed, because of the above error' iiu=-iiunit end if ! ierr=2 end if if (iiu.gt.0) then if (IM1r.ne.IM1.or.JM1r.ne.JM1.or.IM2r.ne.IM2.or.JM2r.ne.JM2& & .or.ANY(lon1r.ne.lon1).or.ANY(lat1r.ne.lat1)& & .or.ANY(LON2r.ne.LON2).or.ANY(LAT2r.ne.LAT2)& & .or.ANY(mask1r.ne.mask1).or.ANY(MASK2r.ne.MASK2)) then print*,'INTP_INI: **ERROR: read INTP_INI input data'//& & ' differ from input arguments' if (ierr_in.eq.0) then ierr=3 else print*,'INTP_INI: output data will be '//& & 'computed, because of the above error' iiu=-iiunit end if end if end if if (iiu.gt.0) RETURN END IF ! 1. Find boundary of G1 domain as 1D arrays ibnd1, jbnd1 (G1 domain is ! where mask1(i,j).ge.0) if (PRESENT(KM1)) km=KM1 IF (PRESENT(KM1) .and. km.ne.0) THEN kmm=0 ! only for printout !! km=KM1 if (ibndr1(0).ne.ibndr1(km) .or. jbndr1(0).ne.jbndr1(km)) then print*,'INTP_INI: **ERROR** domain boundary (input arguments '//& & 'ibndr1,jbndr1) not closed ',& & ibndr1(0),jbndr1(0),ibndr1(km),jbndr1(km) ierr=5 RETURN end if ! allocate(ibnd1(0:km),jbnd1(0:km)) ! print*,'INTP_INI: done allocate(ibnd1(0:km),jbnd1(0:km))' ibnd1=ibndr1 jbnd1=jbndr1 ELSE ! kmm=2*(IM1+JM1-2) kmm=IM1*JM1 ! allocate(ibnd1(0:kmm),jbnd1(0:kmm)) ! 1.1. Find 2D array's boundary point (i1,j1) (where i=1 or i=IM1 or ! j=1 or j=JM1) which is also a G1 domain point call FBP(IM1,JM1,mask1,i1,j1,ierr) if (ierr.ne.0) then print*,& & 'INTP_INI: **ERROR** FAILED to find init. bound. domain 1 point' RETURN end if i=i1 j=j1 ! print*,i,j ! 1.2. Find domain boundary msk1(1:IM1,1:JM1)=mask1 msk1(:,0)=-1 msk1(IM1+1,:)=-1 msk1(:,JM1+1)=-1 msk1(0,:)=-1 ibnd1(0)=i jbnd1(0)=j km=0 do k=1,kmm ! if (msk1(i,j+1).ge.0) then if (msk1(i,j-1).eq.-1) then if (msk1(i+1,j).ge.0) then i=i+1 else j=j+1 end if ! else if (msk1(i-1,j).ge.0) then else if (msk1(i+1,j).eq.-1) then if (msk1(i,j+1).ge.0) then j=j+1 else i=i-1 end if ! else if (msk1(i,j-1).ge.0) then else if (msk1(i,j+1).eq.-1) then if (msk1(i-1,j).ge.0) then i=i-1 else j=j-1 end if ! else if (msk1(i+1,j).ge.0) then else if (msk1(i-1,j).eq.-1) then if (msk1(i,j-1).ge.0) then j=j-1 else i=i+1 end if ! else else if (msk1(i-1,j-1).eq.-1) then j=j-1 else if (msk1(i+1,j-1).eq.-1) then i=i+1 else if (msk1(i+1,j+1).eq.-1) then j=j+1 else if (msk1(i-1,j+1).eq.-1) then i=i-1 else print*,'INTP_INI: **ERROR** BAD DOMAIN 1 BOUNDARY ',i,j ierr=5 RETURN end if ibnd1(k)=i jbnd1(k)=j !! msk1(i,j)=-2 ! print*,k,kmm,km,i,j if (i.eq.ibnd1(0) .and. j.eq.jbnd1(0)) then km=k exit end if end do if (km.eq.0) then print*,'INTP_INI: **ERROR** domain boundary not closed ',& & ibnd1(0),jbnd1(0),ibnd1(km),jbnd1(km) ierr=5 RETURN end if END IF if (nvrb.gt.0) print*,'INTP_INI: km=',km,' kmm=',kmm ! do k=1,km ! print*,k,ibnd1(k),jbnd1(k),mask1(ibnd1(k),jbnd1(k)) ! end do ! 2. Find discrete analog of Peano curve for G2 domain (integer ! arrays IPEANO2, JPEANO2) ! (Now it's not very Peano, as it may jump across empty areas ! i.e. ones where MASK2<0. It is rather a trivial Hamiltonian ! path in grid 2 considered as a graph, with additional edges ! drawn over empty areas.) IG1=-3 is=1 ie=IM2 istep=1 k=0 do j=1,JM2 do i=is,ie,istep if (MASK2(i,j).ge.0) then k=k+1 IPEANO2(k)=i JPEANO2(k)=j end if end do istep=-istep i=is is=ie ie=i end do KM2=k if (nvrb.gt.0) print*,'INTP_INI: KM2=',KM2 ! 3. For each valid G2 point (i2,j2) [i.e. for each (i2,j2) such ! that MASK2(i2,j2)>=0, i.e. for each (i2,j2) lying on the Peano ! curve], find the G1 cell such that (i2,j2) lies inside it ! 3.1. Find if 1st of Peano curve G2 points lies inside G1 domain. ! Finish with this point (i2,j2), i.e. find IG1(i2,j2), JG1(i2,j2) ! if it does lie inside domain 1, otherwise assign IG1(i2,j2)=-2 inG1=.false. i1=ibnd1(0) j1=jbnd1(0) i2=IPEANO2(1) j2=JPEANO2(1) jloop: do j=1,JM1 do i=1,IM1 call G2Point_G1Cell(i2,j2,i,j,inG1,i1,j1) if (inG1) EXIT jloop end do end do jloop if (.not. inG1) IG1(i2,j2)=-2 ! 3.2. Moving along the Peano curve, for each G2 point (i2,j2) ! find the G1 cell such that (i2,j2) lies inside it. Such ! a cell must exist if (i2,j2) lies inside G2 domain; in this ! case, assign IG1(i2,j2), JG1(i2,j2), which define the cell; ! otherwise assign IG1(i2,j2)=-2 ! print*,'PC: 1st point:',inG1,i2,j2,i1,j1 kloop: do k=2,KM2 ! 3.2.1. Find if the current G2 point (i2,j2) on Peano curve ! is in the same position with respect to G1 domain (i.e. ! inside or outside it) as the previous G2 point on Peano ! curve. If it is outside G1 domain assign IG1(i2,j2)=-2 and ! go to the next G2 point i2p=i2 j2p=j2 i2=IPEANO2(k) j2=JPEANO2(k) call G2Point_G1Cell(i2,j2,i1,j1,inG1cell,i1,j1) ! if (present(ia).and.inG1cell) ! & print*,'PC: G1 cell same as previous ',k,inG1,i2,j2,i1,j1 if (inG1cell) then inG1=.true. CYCLE end if cross=.false. n1=0 ! for debugging do kb1=1,km cross=Intersects(& & LON2(i2p,j2p),LAT2(i2p,j2p),LON2(i2,j2),LAT2(i2,j2),& & lon1(ibnd1(kb1-1),jbnd1(kb1-1)),& & lat1(ibnd1(kb1-1),jbnd1(kb1-1)),& & lon1(ibnd1(kb1),jbnd1(kb1)),lat1(ibnd1(kb1),jbnd1(kb1))) !! cross=cross .neqv. Intersects( if (cross) then inG1=.not.inG1 n1=n1+1 ! for debugging k1=kb1 end if ! inG1=inG1 .neqv. cross !<- does same as above !! if (cross) EXIT <- deleted. One segment of G2 Peano curve !! may intersect G1 boundary more than once, !! so the criterion of the curve getting !! inside G1 domain from outside, or vice !! versa, is an odd number of crossings end do !! inG1=inG1.neqv.cross if (inG1 .and. cross) then ! this is for some acceleration but i1=ibnd1(k1) ! it is to be seen if it yields any j1=jbnd1(k1) end if if (.not. inG1) then IG1(i2,j2)=-2 ! this is better than in INTERP_INIT: now all external points get IG1=-2 . Check if it is OK in Coupler and change comments in mod_interp.f.new CYCLE end if ! if (present(ia)) print*,'PC: k=',k,inG1,i2,j2,i1,j1 ! 3.2.2. Starting from (i1,j1), on perimeters of expanding squares, ! find (i1,j1) such as G2Point_In_G1Cell(i2,j2,i1,j1), i.e. ! such as G2 point (i2,j2) lies inside grid 1 cell defined ! by G1 point (i1,j1) ! call G2Point_G1Cell(i2,j2,i1,j1,inG1cell,i1,j1) ! print*,'PC: G1 cell found at 1st guess ',k,inG1,i2,j2,i1,j1 ! if (inG1cell) CYCLE is=i1 ie=i1 js=j1 je=j1 do while (.true.) isp=is iep=ie jsp=js jep=je is=max(isp-1,1) ie=min(iep+1,IM1) js=max(jsp-1,1) je=min(jep+1,JM1) stuck=.true. if (is.lt.isp) then stuck=.false. do j=js,je call G2Point_G1Cell(i2,j2,is,j,inG1cell,i1,j1) if (inG1cell) CYCLE kloop end do end if if (ie.gt.iep) then stuck=.false. do j=js,je call G2Point_G1Cell(i2,j2,ie,j,inG1cell,i1,j1) if (inG1cell) CYCLE kloop end do end if if (js.lt.jsp) then stuck=.false. do i=isp,iep call G2Point_G1Cell(i2,j2,i,js,inG1cell,i1,j1) if (inG1cell) CYCLE kloop end do end if if (je.gt.jep) then stuck=.false. do i=isp,iep call G2Point_G1Cell(i2,j2,i,je,inG1cell,i1,j1) if (inG1cell) CYCLE kloop end do end if ! End 3.2.2 if (stuck) then ierr=4 inG1cell=.false. jloop1: do j=1,JM1 do i=1,IM1 call G2Point_G1Cell(i2,j2,i,j,inG1cell,i1,j1) if (inG1cell) EXIT jloop1 end do end do jloop1 if (inG1cell) then print*,'INTP_INI: the G2 point DOES lie in a G1 cell '//& & 'but it was not found by expanding squares ',& & i1,j1,lon1(i1,j1),lon2(i1,j1) else print*,'INTP_INI: the G2 point does NOT lie in any G1 '//& & 'cell' end if print*,'INTP_INI: **ERROR** STUCK, the sought G1 cell '//& & 'must exist but is not found ',& & k,i2,j2,LON2(i2,j2),LAT2(i2,j2),MASK2(i2,j2),n1,IG1(i2,j2),& & IG1(i2p,j2p) RETURN cycle kloop end if end do ! print*,'PC: G1 cell found ',k,inG1,i2,j2,i1,j1 ! end do kloop ! IF (iiu.lt.0) THEN write (-iiu,iostat=ios) IM1,JM1,IM2,JM2,& & lon1,lat1,LON2,LAT2,mask1,MASK2,& & IG1,JG1 ios1=0 if (present(COSANG)) write (-iiu,iostat=ios1) COSANG,SINANG if (ios.ne.0 .or. ios1.ne.0) then print*,'INTP_INI: **ERROR** writing output data ',& & ios,ios1 ierr=6 end if END IF ! if (nvrb.gt.0) then print*,'INTP_INI: done, IM1,JM1,IM2,JM2: ',IM1,JM1,IM2,JM2 k=0 n=0 k1=0 i1=0 j1=0 j2=0 do j=1,JM2 do i=1,IM2 if (MASK2(i,j).lt.0) then k1=k1+1 else if (MASK2(i,j).eq.1) then i1=i1+1 end if if (IG1(i,j).gt.0) then k=k+1 else if (IG1(i,j).eq.-2) then j1=j1+1 if (MASK2(i,j).eq.0) j2=j2+1 else if (IG1(i,j).eq.0) then n=n+1 end if end do end do print*,'INTP_INI: Total # of G2 points (2D array size) ',IM2*JM2 print*,'INTP_INI: Total # of valid G2 points ',& & IM2*JM2-k1 print*,'INTP_INI: # of land G2 points ',i1 print*,'INTP_INI: # of sea G2 points ',& & IM2*JM2-k1-i1 print*,'INTP_INI: # of sea G2 points inside G1 domain ',& & IM2*JM2-k1-i1-j2 print*,'INTP_INI: # of interpolatable G2 (sea) points ',k print*,'INTP_INI: # of G2 points outside G1 domain ',j1 print*,'INTP_INI: # of sea G2 points outside G1 domain ',j2 print*,'INTP_INI: # of sea G2 points in land G1 cells ',n end if ! RETURN !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! CONTAINS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine FBP(IM,JM,mask,i0,j0,ierr) implicit none integer IM,JM,mask(IM,JM),i0,j0,ierr integer i,j ! ierr=0 do i=1,IM if (mask(i,1).ge.0) then i0=i j0=1 return end if end do do i=IM,1,-1 if (mask(i,JM).ge.0) then i0=i j0=JM return end if end do do j=1,JM if (mask(IM,j).ge.0) then i0=IM j0=j return end if end do do j=JM,1,-1 if (mask(1,j).ge.0) then i0=1 j0=j return end if end do ierr=10 return end subroutine FBP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine G2Point_G1Cell(i,j,i1,j1,inside,i1out,j1out) implicit none integer,intent(in):: i,j,i1,j1 logical,intent(out):: inside integer,intent(out):: i1out,j1out logical cw ! inside=.false. if (MASK2(i,j).lt.0) then ! temporary, for debugging only print*,'INTP_INI: **ERROR** invalid G2 point in G2Point_G1Cell' STOP end if if (present(ia)) then i1a=ia(i1,j1) if (i1a.lt.1) return if (i1a.gt.IM1) return j1a=ja(i1,j1) if (j1a.lt.1) return if (j1a.gt.JM1) return i1b=ib(i1,j1) if (i1b.lt.1) return if (i1b.gt.IM1) return j1b=jb(i1,j1) if (j1b.lt.1) return if (j1b.gt.JM1) return i1c=ic(i1,j1) if (i1c.lt.1) return if (i1c.gt.IM1) return j1c=jc(i1,j1) if (j1c.lt.1) return if (j1c.gt.JM1) return else if (i1.ge.IM1) return if (j1.ge.JM1) return i1a=i1+1 j1a=j1 i1b=i1+1 j1b=j1+1 i1c=i1 j1c=j1+1 end if !B Check if (i,j) is inside quadrangle ! (i1,j1),(i1a,j1a),(i1b,j1b),(i1c,j1c) ! It is not assumed here that this list of vertices ! is counterclockwise if (mask1(i1,j1).eq.-1) return if (mask1(i1a,j1a).eq.-1) return if (mask1(i1b,j1b).eq.-1) return if (mask1(i1c,j1c).eq.-1) return x1=lon1(i1,j1)-LON2(i,j) y1=lat1(i1,j1)-LAT2(i,j) x1a=lon1(i1a,j1a)-LON2(i,j) y1a=lat1(i1a,j1a)-LAT2(i,j) cw=x1*y1a-x1a*y1.lt.0. x1b=lon1(i1b,j1b)-LON2(i,j) y1b=lat1(i1b,j1b)-LAT2(i,j) if (x1a*y1b-x1b*y1a.lt.0. .neqv. cw) return x1c=lon1(i1c,j1c)-LON2(i,j) y1c=lat1(i1c,j1c)-LAT2(i,j) if (x1b*y1c-x1c*y1b.lt.0. .neqv. cw) return if (x1c*y1-x1*y1c.lt.0. .neqv. cw) return !E inside=.true. i1out=i1 j1out=j1 if (MASK2(i,j).eq.0) then if (mask1(i1,j1).gt.0 .and. mask1(i1a,j1a).gt.0 .and.& & mask1(i1b,j1b).gt.0 .and. mask1(i1c,j1c).gt.0) then IG1(i,j)=0 return else IG1(i,j)=i1 JG1(i,j)=j1 end if else ! MASK(i,j) must be 1 (land), since (i,j) is supposed ! a valid G2 point, since it is supposed to lie inside ! a G1 cell IG1(i,j)=-1 return end if if (PRESENT(COSANG)) then if (present(XDLON2)) then dlon2=XDLON2(i,j) dlat2=XDLAT2(i,j) else if (i.eq.IM2) then dlon2=LON2(i,j)-LON2(i-1,j) dlat2=LAT2(i,j)-LAT2(i-1,j) else if (i.eq.1) then dlon2=LON2(i+1,j)-LON2(i,j) dlat2=LAT2(i+1,j)-LAT2(i,j) else dlon2=LON2(i+1,j)-LON2(i-1,j) dlat2=LAT2(i+1,j)-LAT2(i-1,j) end if end if dlon2=dlon2*cos(LAT2(i,j)*rad) ! print*,'INTP_INI to call BLIN',k,i,j CALL BLIN(lon1(i1,j1),lon1(i1a,j1a),& & lon1(i1b,j1b),lon1(i1c,j1c),& & lat1(i1,j1),lat1(i1a,j1a),& & lat1(i1b,j1b),lat1(i1c,j1c),& & dlon1(i1,j1),dlon1(i1a,j1a),& & dlon1(i1b,j1b),dlon1(i1c,j1c), & & LON2(i,j),LAT2(i,j),& & dlon1on2,& & ivrbs,ierr1) if (ierr1.ne.0) then ierr=100 print*,'INTP_INI: **ERROR** in BLIN ',ierr RETURN end if dlon1on2=dlon1on2*cos(LAT2(i,j)*rad) CALL BLIN(lon1(i1,j1),lon1(i1a,j1a),& & lon1(i1b,j1b),lon1(i1c,j1c),& & lat1(i1,j1),lat1(i1a,j1a),& & lat1(i1b,j1b),lat1(i1c,j1c),& & dlat1(i1,j1),dlat1(i1a,j1a),& & dlat1(i1b,j1b),dlat1(i1c,j1c), & & LON2(i,j),LAT2(i,j),& & dlat1on2,& & ivrbs,ierr1) r1r2=sqrt((dlon2**2+dlat2**2)*(dlon1on2**2+dlat1on2**2)) SINANG(i,j)=(dlon2*dlat1on2-dlon1on2*dlat2)/r1r2 COSANG(i,j)=(dlon2*dlon1on2+dlat2*dlat1on2)/r1r2 ! end if return end subroutine G2Point_G1Cell !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! logical function Intersects(x1s,y1s,x1e,y1e,x2s,y2s,x2e,y2e) implicit none real x1s,y1s,x1e,y1e,x2s,y2s,x2e,y2e real dx1,dy1,dx2,dy2,dxs,dys,D,Dt,Du ! Intersects=.false. dx1=x1e-x1s dy1=y1e-y1s dx2=x2e-x2s dy2=y2e-y2s D=dx2*dy1-dx1*dy2 dxs=x2s-x1s dys=y2s-y1s Dt=dys*dx2-dxs*dy2 Du=dys*dx1-dxs*dy1 if (D.lt.0.) then D=-D Dt=-Dt Du=-Du end if if (Dt.le.0.) return if (Dt.gt.D) return if (Du.le.0.) return Intersects=Du.le.D return end function Intersects !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! END SUBROUTINE INTP_INI ! !*********************************************************************** ! SUBROUTINE INTERP(IM1,JM1,IM2,JM2,lon1,lat1,LON2,LAT2,& &mask1,MASK2,IG1,JG1,fg1,& &FG2,ierr,& &ia,ja,ib,jb,ic,jc,NOSPVAL) !======================================================================= ! ! Interpolation from G1 grid to G2 grid ! --------------------------------------- ! ! lon1(IM1,JM1), lat1(IM1,JM1) (IN) is G1 grid coordinates ! LON2(IM2,JM2), LAT2(IM2,JM2) (IN) is G2 grid coordinates ! ! mask1 (IN) is sea/land mask on G1 grid (0 - sea, 1 - land) ! MASK2 (IN) is sea/land mask on G2 grid (0 - sea, 1 - land) ! IG1(IM2,JM2), JG1(IM2,JM2) (IN): for each G2 gridpoint (i,j), ! whose coordinates are (LON2(i,j),LAT2(i,j)), ! ( IG1(i,j),JG1(i,j) ) = (i1,j1) ! is the point on G1 grid for which (i,j) lies in the ! quadrangle whose angles are G1 gridpoints (i1,j1), (i1+1,j1), ! (i1+1,j1+1), (i1,j1+1). ! If ia,ja,ib,jb,ic,jc are present, the 4 vertices are, ! instead of the above, (i1,j1), (ia(i1,j1),ja(i1,j1)), ! (ib(i1,j1),jb(i1,j1)), (ic(i1,j1),jc(i1,j1)) ! ! fg1 (IN) is data on G1 grid to be interpolated from ! FG2 (OUT) is data on G2 grid to be obtained by interpolation ! ! !======================================================================= implicit none integer,intent(in):: IM1,JM1,IM2,JM2 !<- sea INTP_INI real,dimension(IM1,JM1),intent(in):: lon1,lat1 !<- sea INTP_INI real,dimension(IM2,JM2),intent(in):: LON2,LAT2 !<- sea INTP_INI integer,dimension(IM1,JM1),intent(in):: mask1 !<- sea INTP_INI integer,dimension(IM2,JM2),intent(in):: MASK2 !<- sea INTP_INI integer,dimension(IM2,JM2),intent(in):: IG1,JG1 !<- sea INTP_INI real,dimension(IM1,JM1),intent(in):: fg1(IM1,JM1) !<- sea above REAL,DIMENSION(IM1,JM1),INTENT(OUT):: FG2(IM2,JM2) !<- sea above INTEGER,INTENT(OUT):: ierr !<- 0 for normal and nonzero for abnormal return integer,optional,dimension(IM1,JM1),intent(in):: ia,ja,ib,jb,ic,jc !<- sea INTP_INI logical,optional,intent(in):: NOSPVAL !<- if present and =.true., no value is assigned to ! FG2(i,j) for those G2 gridpoints (i,j) which lie ! outside G1 domain (i.e. do not lie in any G1 grid ! cell) or are G2 landpoints. Otherwise, special ! value -1.E30 is assigned at such G2 gridpoints integer ierr1,ivrbs,k,i,j,i1,j1,i1a,j1a,i1b,j1b,i1c,j1c real s,x,y,f,x1,y1,f1,x1a,y1a,f1a,x1b,y1b,f1b,x1c,y1c,f1c,& &very_large_negative parameter (very_large_negative=-1.E30) logical SPVAL ! ! print*,'C: INTERP ENTERED, IM1,JM1,IM2,JM2: ',IM1,JM1,IM2,JM2 if ( (PRESENT(ia).neqv.PRESENT(ja)) .or.& & (PRESENT(ia).neqv.PRESENT(ib)) .or.& & (PRESENT(ia).neqv.PRESENT(jb)) .or.& & (PRESENT(ia).neqv.PRESENT(ic)) .or.& & (PRESENT(ia).neqv.PRESENT(jc)) ) then print*,'INTERP: **ERROR** ia,ja,ib,jb,ic,jc must be all '//& & 'present or all not present' ierr=1 RETURN end if if (PRESENT(NOSPVAL)) then SPVAL=.not.NOSPVAL else SPVAL=.true. end if ierr=0 ivrbs=10 ! DO j=1,JM2 DO i=1,IM2 ! i1=IG1(i,j) j1=JG1(i,j) IF (MASK2(i,j).ne.0 .or. i1.le.0) THEN if (SPVAL) FG2(i,j)=very_large_negative CYCLE END IF if (present(ia)) then i1a=ia(i1,j1) ! here, it is supposed that the previous j1a=ja(i1,j1) ! IF statement ensures that here i1>0, i1b=ib(i1,j1) ! which, in its turn, ensures that j1b=jb(i1,j1) ! i1,i1a,i1b,i1c stay within [1,IM1] and i1c=ic(i1,j1) ! j1,j1a,j1b,j1c stay within [1,JM1] j1c=jc(i1,j1) ! else i1a=i1+1 j1a=j1 i1b=i1+1 j1b=j1+1 i1c=i1 j1c=j1+1 end if ! !b Bilinear interpolation is used to obtain a G2 gridpoint !b value from the values at the four surrounding G1 gridpoints. !b To prevent interpolation from land, G1 land gridpoint !b values are substituted by the arithmetic mean of the k remaining !b G1 sea gridpoint values (0itmaxmax are not filled with any value, nor is it(,) ! computed for these points. ! If it(IM,JM) is not present, itmax, if present, is assigned the ! maximum of minimum number of gridsteps to a gridpoint with a ! previously defined value, which may be either an initially defined ! (input) value (where intpmask(,)=.true.) or an already extrapolated ! value. itmaxmax, if present, puts a limit on itmax so that no ! extrapolation is done if the nearest gridpoint with a previously ! defined value is farther than itmaxmax gridsteps. Thus, note ! that arguments itmax, itmaxmax have different meaning depending ! on the presence of argument it . ! Optional argument f_r serves those gripoints which failed to ! be filled by extrapolation and smooth averaging at extrapolation ! points is carried out; see also below. implicit none integer,intent(in):: IM,JM ! see real,intent(inout):: f(IM,JM) ! above logical,intent(in):: mask(IM,JM),intpmask(IM,JM) ! integer,intent(out):: ierr ! if there are lakes and lakeval ! is present, retuns a negative; ! otherwise if return is abnormal, ! including presence of a lake ! when lakeval is not present, ! returns a positive; otherwise ! returns 0 integer,optional,intent(in):: itmaxmax ! see integer,optional,intent(out):: itmax ! above real,optional,intent(in):: lakeval ! integer,optional,intent(in):: ivrbs ! prints lake info if >0 integer,optional,intent(out):: it(IM,JM) ! see above real,optional,intent(in):: f_r(IM,JM) ! reference values for f. ! If present, fills the ! unfilled gridpoints and ! extrapolated values f(i,j) ! are replaced by weighted ! average betw. themselves ! and f_r(i,j), with smooth ! transition. Must only be ! present together with ! itmaxmax and it . integer i,j,i1,j1,i2,j2,ii,jj,i1p,j1p,i2p,j2p,k,k1,n,m,ivr,itmm logical,dimension(IM,JM):: cc,imask logical nolake real s,w ! if (PRESENT(f_r) .and. .not.(PRESENT(itmaxmax).and. & &PRESENT(it))) then print*,'EXTRAP: f_r present but itmaxmax or it is not, aborted' ierr=3 RETURN end if ierr=0 imask=intpmask if (PRESENT(ivrbs)) then ivr=ivrbs else ivr=0 end if if (PRESENT(itmax)) itmax=0 DO j=1,JM iloop: DO i=1,IM m=0 IF (mask(i,j) .AND. .NOT.imask(i,j)) THEN i1=i i2=i j1=j j2=j cc(i,j)=.true. n=0 s=0. do while (n.eq.0) m=m+1 if (PRESENT(itmaxmax)) then if (m.gt.itmaxmax) then if (PRESENT(f_r)) f(i,j)=f_r(i,j) CYCLE iloop end if end if i1p=i1 i2p=i2 j1p=j1 j2p=j2 i1=max(i1-1,1) i2=min(i2+1,IM) j1=max(j1-1,1) j2=min(j2+1,JM) if (i1p.eq.1.and.i2p.eq.IM.and.j1p.eq.1.and.j2p.eq.JM) then print '("EXTRAP: FAILED/TERMINATED. After ",i5," steps ",& & "connected component for i=",i5,", j=",i5,& & " became the entire domain")',m,i,j call PRM ierr=1 RETURN end if if (j1.lt.j1p) then cc(i1p:i2p,j1)=cc(i1p:i2p,j1p) .and. mask(i1p:i2p,j1) nolake=ANY(cc(i1p:i2p,j1)) if (i1.lt.i1p) cc(i1,j1)=.false. end if if (i2.gt.i2p) then cc(i2,j1p:j2p)=cc(i2p,j1p:j2p) .and. mask(i2,j1p:j2p) nolake=nolake.or.ANY(cc(i2,j1p:j2p)) if (j1.lt.j1p) cc(i2,j1)=.false. end if if (j2.gt.j2p) then cc(i1p:i2p,j2)=cc(i1p:i2p,j2p) .and. mask(i1p:i2p,j2) nolake=nolake.or.ANY(cc(i1p:i2p,j2)) if (i2.gt.i2p) cc(i2,j2)=.false. end if if (i1.lt.i1p) then cc(i1,j1p:j2p)=cc(i1p,j1p:j2p) .and. mask(i1,j1p:j2p) nolake=nolake.or.ANY(cc(i1,j1p:j2p)) if (j2.gt.j2p) cc(i1,j2)=.false. end if if (.not.nolake) then if (ivr.gt.0 .or. & & .not. (PRESENT(lakeval).or.PRESENT(f_r)) ) then print '("EXTRAP: Domain not connected. On step",i5,& & " i=",i5," j=",i5," turned out a lake point")',m,i,j call PRM end if if (PRESENT(lakeval) .or. PRESENT(f_r)) then imask(i,j)=.true. if (PRESENT(lakeval)) then f(i,j)=lakeval if (ivr.gt.0) then print*,'EXTRAP: lake value lakeval=',lakeval,& & ' assigned' end if else f(i,j)=f_r(i,j) if (ivr.gt.0) then print*,'EXTRAP: at lake point, f_r value ',f_r(i,j),& & ' assigned' end if end if ierr=ierr-1 CYCLE iloop else print*,'EXTRAP: connectedness is required when '//& & 'neither lakeval nor f_r is present. TERMINATED' ierr=2 RETURN end if end if do k=1,2 do ii=i1+1,i2 if (cc(ii-1,j1)) cc(ii,j1)=mask(ii,j1) end do do jj=j1+1,j2 if (cc(i2,jj-1)) cc(i2,jj)=mask(i2,jj) end do do ii=i2-1,i1,-1 if (cc(ii+1,j2)) cc(ii,j2)=mask(ii,j2) end do do jj=j2-1,j1,-1 if (cc(i1,jj+1)) cc(i1,jj)=mask(i1,jj) end do do jj=j1+1,j2 if (cc(i1,jj-1)) cc(i1,jj)=mask(i1,jj) end do do ii=i1+1,i2 if (cc(ii-1,j2)) cc(ii,j2)=mask(ii,j2) end do do jj=j2-1,j1,-1 if (cc(i2,jj+1)) cc(i2,jj)=mask(i2,jj) end do do ii=i2-1,i1,-1 if (cc(ii+1,j1)) cc(ii,j1)=mask(ii,j1) end do end do IF (PRESENT(it)) THEN k1=IM+JM if (PRESENT(itmaxmax)) then itmm=itmaxmax else itmm=k1 end if do ii=i1+1,i2 if (cc(ii,j1).and.imask(ii,j1)) then k=it(ii,j1)+m k1=min(k,k1) if (k.le.itmm) then s=s+f(ii,j1) n=n+1 end if end if end do do jj=j1+1,j2 if (cc(i2,jj).and.imask(i2,jj)) then k=it(i2,jj)+m k1=min(k,k1) if (k.le.itmm) then s=s+f(i2,jj) n=n+1 end if end if end do do ii=i2-1,i1,-1 if (cc(ii,j2).and.imask(ii,j2)) then k=it(ii,j2)+m k1=min(k,k1) if (k.le.itmm) then s=s+f(ii,j2) n=n+1 end if end if end do do jj=j2-1,j1,-1 if (cc(i1,jj).and.imask(i1,jj)) then k=it(i1,jj)+m k1=min(k,k1) if (k.le.itmm) then s=s+f(i1,jj) n=n+1 end if end if end do ELSE do ii=i1+1,i2 if (cc(ii,j1).and.imask(ii,j1)) then s=s+f(ii,j1) n=n+1 end if end do do jj=j1+1,j2 if (cc(i2,jj).and.imask(i2,jj)) then s=s+f(i2,jj) n=n+1 end if end do do ii=i2-1,i1,-1 if (cc(ii,j2).and.imask(ii,j2)) then s=s+f(ii,j2) n=n+1 end if end do do jj=j2-1,j1,-1 if (cc(i1,jj).and.imask(i1,jj)) then s=s+f(i1,jj) n=n+1 end if end do END IF end do f(i,j)=s/n if (PRESENT(it)) then it(i,j)=k1 if (PRESENT(f_r)) then w=max(0.,(itmaxmax+1.-k1)/(itmaxmax+1.)) f(i,j)=w*f(i,j)+(1.-w)*f_r(i,j) end if else k1=m end if if (PRESENT(itmax)) itmax=max(itmax,k1) imask(i,j)=.true. END IF END DO iloop END DO RETURN CONTAINS subroutine PRM i1=max(1,i-10) i2=min(IM,i+10) j1=max(1,j-10) j2=min(JM,j+10) do jj=j1,j2 print*,'EXTRAP: j=',jj,' i from ',i1,' to ',i2,' mask: ',& & mask(i1:i2,jj) end do do jj=j1,j2 print*,'EXTRAP: j=',jj,' i from ',i1,' to ',i2,' imask: ',& & imask(i1:i2,jj) end do end subroutine PRM END SUBROUTINE EXTRAP ! !************************************************************ ! END module mod_interp