!------------------------------------------------------------------------------ ! Global Modeling and Assimilation Office (GMAO) ! ! Goddard Earth Observing System (GEOS) ! ! MAPL Component ! !------------------------------------------------------------------------------ ! #include "MAPL_Generic.h" #include "unused_dummy.H" ! !> !### MODULE: `MAPL_OrbGridCompMod` ! ! Author: GMAO SI-Team ! ! The module `MAPL_OrbGridCompMod` is an ESMF gridded component implementing ! the MAPL-5 Nominal Orbit Calculator. ! ! It was developed for MAPL-5 release Fortuna 2.3 and later. ! !#### History !- 30Nov2010 da Silva Initial version. ! MODULE MAPL_OrbGridCompMod ! ! !USES: ! Use ESMF Use ESMFL_Mod Use MAPL_BaseMod Use MAPL_GenericMod Use MAPL_Constants Use MAPL_CommsMod, only: MAPL_AM_I_ROOT Use MAPL_ErrorHandlingMod IMPLICIT NONE PRIVATE ! ! !PUBLIC MEMBER FUNCTIONS: PUBLIC SetServices ! Legacy state ! ------------ TYPE Orb_State PRIVATE type (ESMF_Config) :: CF ! Private Config integer :: no ! Number of masks character(len=ESMF_MAXSTR), pointer :: Instrument(:) character(len=ESMF_MAXSTR), pointer :: Satellite(:) real, pointer :: Swath(:) integer, pointer :: halo(:) logical :: verbose=.FALSE. real, pointer :: x(:,:), y(:,:) integer :: face END TYPE Orb_State ! Hook for the ESMF ! ----------------- TYPE Orb_Wrap TYPE (Orb_State), pointer :: PTR => null() END TYPE Orb_WRAP ! !------------------------------------------------------------------------------ CONTAINS !------------------------------------------------------------------------------ !> ! Sets IRF services for the Orb Grid Component. ! Sets Initialize, Run and Finalize services. ! SUBROUTINE SetServices ( GC, RC ) type(ESMF_GridComp), intent(INOUT) :: GC ! gridded component integer, intent(out), optional :: RC ! return code ! Local derived type aliases ! -------------------------- type (Orb_State), pointer :: self ! internal, that is type (Orb_wrap) :: wrap character(len=ESMF_MAXSTR) :: comp_name integer :: i, nCols _Iam_('SetServices') integer :: status logical :: found ! ------------ ! Get my name and set-up traceback handle ! --------------------------------------- call ESMF_GridCompGet( GC, name=comp_name, _RC ) Iam = TRIM(comp_name) // '::' // TRIM(Iam) ! Greetings ! --------- IF(MAPL_AM_I_ROOT()) THEN PRINT *, TRIM(Iam)//': ACTIVE' END IF ! Wrap internal state for storing in GC; rename legacyState ! ------------------------------------- allocate ( self, stat=STATUS ) _VERIFY(STATUS) wrap%ptr => self ! Load private Config Attributes ! ------------------------------ self%CF = ESMF_ConfigCreate(_RC) inquire(file="MAPL_OrbGridComp.rc", exist=found) if (found) then call ESMF_ConfigLoadFile ( self%CF,'MAPL_OrbGridComp.rc',rc=status) _VERIFY(STATUS) call ESMF_ConfigGetAttribute(self%CF, self%verbose, Label='verbose:', default=.false. , _RC ) ! ------------------------ ! Get Mask Definitions ! ------------------------ call ESMF_ConfigGetDim(self%CF, self%no, nCols, LABEL='Nominal_Orbits::',_RC) _ASSERT(self%no>0,'needs informative message') allocate(self%Instrument(self%no), self%Satellite(self%no), & self%Swath(self%no), self%halo(self%no), __STAT__) if ( self%verbose .AND. MAPL_AM_I_ROOT() ) then write(*,*)" Swath" write(*,*)"Instrument Satellite (km) Halo Width" write(*,*)"--------------- ----------- --------- -------------" end if call ESMF_ConfigFindLabel(self%CF, 'Nominal_Orbits::',_RC) do i = 1, self%no call ESMF_ConfigNextLine(self%CF,_RC) call ESMF_ConfigGetAttribute(self%CF,self%Instrument(i),_RC) call ESMF_ConfigGetAttribute(self%CF,self%Satellite(i),_RC) call ESMF_ConfigGetAttribute(self%CF,self%Swath(i),_RC) call ESMF_ConfigGetAttribute(self%CF,self%halo(i),_RC) if ( self%verbose .AND. MAPL_AM_I_ROOT() ) then write(*,'(1x,a15,4x,a11,4x,f9.1,4x,i3)') self%Instrument(i), self%Satellite(i), self%Swath(i), self%halo(i) end if end do else self%no = 0 endif ! ------------------------ ! ESMF Functional Services ! ------------------------ ! Set the Initialize, Run, Finalize entry points ! ---------------------------------------------- call MAPL_GridCompSetEntryPoint ( GC, ESMF_METHOD_INITIALIZE, Initialize_, RC=STATUS) _VERIFY(STATUS) call MAPL_GridCompSetEntryPoint ( GC, ESMF_METHOD_RUN, Run_, RC=STATUS) _VERIFY(STATUS) ! Store internal state in GC ! -------------------------- call ESMF_UserCompSetInternalState ( GC, 'Orb_state', wrap, STATUS ) _VERIFY(STATUS) ! ------------------ ! MAPL Data Services ! ! in addition to bundle add each instrument as export in case we want to write out in history do i=1,self%no call MAPL_AddExportSpec(GC, & SHORT_NAME = trim(self%Instrument(i)) , & UNITS = 'days' , & DIMS = MAPL_DimsHorzOnly , & RC = STATUS ) enddo call MAPL_AddExportSpec(GC, & SHORT_NAME = 'SATORB', & LONG_NAME = 'Satellite_orbits', & UNITS = 'days' , & DIMS = MAPL_DimsHorzOnly , & DATATYPE = MAPL_BundleItem , & RC = STATUS ) _VERIFY(STATUS) call MAPL_TimerAdd (gc,name="Run" ,rc=status) ! Generic Set Services ! -------------------- call MAPL_GenericSetServices ( GC, _RC ) ! All done ! -------- _RETURN(ESMF_SUCCESS) END SUBROUTINE SetServices !------------------------------------------------------------------------------ !> ! The Initialize method of the orbital component, sets up bundle. ! subroutine Initialize_( GC, IMPORT, EXPORT, CLOCK, RC ) type(ESMF_GridComp), intent(inout) :: GC !! Gridded component type(ESMF_State), intent(inout) :: IMPORT !! Import state type(ESMF_State), intent(inout) :: EXPORT !! Export state type(ESMF_Clock), intent(inout) :: CLOCK !! The clock integer, optional, intent( out) :: RC !! Error code ! ErrLog Variables character(len=ESMF_MAXSTR) :: IAm integer :: STATUS character(len=ESMF_MAXSTR) :: COMP_NAME ! Local type(ESMF_VM) :: VM type(ESMF_FIELD) :: FIELD type(ESMF_GRID) :: GRID type(ESMF_FIELDBUNDLE) :: BUNDLE type(Orb_state), pointer :: self ! Legacy state type(Orb_Wrap) :: wrap integer :: KND, HW, DIMS, LOCATION integer :: i ! New stuff for lat-lon grid needed if doing cube-sphere type (MAPL_MetaComp), pointer :: MAPL_OBJ character(len=ESMF_MAXSTR) :: gridtype_default character(len=ESMF_MAXSTR) :: gridtype ! extra things for cubed sphere integer :: IM, JM, face real(ESMF_KIND_R8), pointer :: EdgeLons(:,:), EdgeLats(:,:) type(ESMF_Info) :: infoh ! Begin... ! Get the target components name and set-up traceback handle. ! ----------------------------------------------------------- call ESMF_GridCompGet ( GC, name=COMP_NAME, VM=VM, RC=STATUS ) _VERIFY(STATUS) Iam = trim(COMP_NAME)//'::Initialize_' call MAPL_GenericInitialize ( GC, IMPORT, EXPORT, CLOCK, RC=STATUS) _VERIFY(STATUS) call ESMF_UserCompGetInternalState(gc, 'Orb_state', WRAP, STATUS) _VERIFY(STATUS) self => wrap%ptr if (self%no == 0) then _RETURN(ESMF_SUCCESS) endif call ESMF_StateGet(EXPORT,'SATORB',BUNDLE,RC=STATUS) _VERIFY(STATUS) call ESMF_GridCompGet ( GC, grid=GRID, RC=STATUS) _VERIFY(STATUS) ! set some info about the fields we will be adding . . . HW=0 DIMS=MAPL_DimsHorzOnly LOCATION=MAPL_VLocationCenter KND=MAPL_R4 do i = 1, self%no field=mapl_FieldCreateEmpty(trim(self%Instrument(i)),Grid,RC=STATUS) _VERIFY(STATUS) call MAPL_FieldAllocCommit(field, dims=dims, location=location, typekind=knd, hw=hw, RC=STATUS) _VERIFY(STATUS) call MAPL_FieldBundleAdd(Bundle,Field,RC=STATUS) _VERIFY(STATUS) enddo ! find out what type of grid we are on, if so gridtype_default='Lat-Lon' call ESMF_InfoGetFromHost(Grid,infoh,rc=status) _VERIFY(STATUS) call ESMF_InfoGet(infoh,key='GridType',value=gridtype,default=gridtype_default,rc=status) _VERIFY(STATUS) if (gridtype=='Cubed-Sphere') then call MAPL_GetObjectFromGC(GC,MAPL_OBJ,rc=status) _VERIFY(STATUS) call MAPL_Get(MAPL_OBJ, im=im, jm=jm, rc=status) _VERIFY(STATUS) allocate(EdgeLons(IM+1,JM+1),stat=status) _VERIFY(STATUS) allocate(EdgeLats(IM+1,JM+1),stat=status) _VERIFY(STATUS) call MAPL_GridGetCorners(Grid,EdgeLons,EdgeLats,rc=status) _VERIFY(STATUS) call check_face(IM+1,JM+1,EdgeLons,EdgeLats,FACE) self%face=face allocate(self%x(IM+1,JM+1),self%y(IM+1,JM+1)) call cube_xy(IM+1,JM+1,self%x,self%y,EdgeLons,EdgeLats,face) deallocate(EdgeLons) deallocate(EdgeLats) endif _RETURN(ESMF_SUCCESS) end subroutine Initialize_ ! !----------------------------------------------------------------------------------------------------- !> ! `Run_` --- Runs Orb ! SUBROUTINE Run_ ( gc, IMPORT, EXPORT, CLOCK, rc ) ! !USES: implicit NONE ! !INPUT PARAMETERS: type(ESMF_Clock), intent(inout) :: CLOCK !! The clock ! !OUTPUT PARAMETERS: type(ESMF_GridComp), intent(inout) :: GC !! Grid Component type(ESMF_State), intent(inout) :: IMPORT !! Import State type(ESMF_State), intent(inout) :: EXPORT !! Export State integer, optional :: rc !! Error return code: !! 0 - all is well !! 1 - ! local type (ESMF_VM) :: VM type (MAPL_MetaComp), pointer :: MAPL_OBJ integer :: IM,JM,LM real, pointer, dimension(:,:) :: LONS real, pointer, dimension(:,:) :: LATS real, pointer, dimension(:,:) :: PTR_TMP, PTR_TMP_EX integer :: iyr,imm,idd,ihr,imn,isc type(ESMF_TimeInterval) :: timeinterval type(ESMF_Time) :: IntervalTime integer, dimension(2) :: interval_nymd, interval_nhms, ihalo, jhalo ! dates and halos for masking routine integer :: deltat ! sub interval for satellite propagation real, dimension(3) :: swath ! satellite swath for masking routine real :: undef ! undefined for masking routine character(len=ESMF_MAXSTR) :: sat_name ! Satellite name type(Orb_state), pointer :: self ! Legacy state type(ESMF_Grid) :: Grid ! Grid type(ESMF_Time) :: Time ! Current time type(ESMF_Config) :: CF ! Universal Config integer :: k, nymd, nhms ! date, time character(len=ESMF_MAXSTR) :: comp_name character(len=ESMF_MAXSTR) :: gridtype_default character(len=ESMF_MAXSTR) :: gridtype type(ESMF_FieldBundle) :: BUNDLE type(ESMF_Info) :: infoh integer :: NORB integer :: IM_world,JM_world,counts(5),imsize integer :: status _Iam_('Run_') _UNUSED_DUMMY(IMPORT) ! Get my name and set-up traceback handle ! --------------------------------------- call ESMF_GridCompGet( GC, name=comp_name, VM=VM, _RC ) Iam = trim(comp_name) // '::Run' ! Extract relevant runtime information ! ------------------------------------ call extract_ ( GC, CLOCK, self, GRID, CF, time, nymd, nhms, timeinterval, _RC) if (self%no == 0) then _RETURN(ESMF_SUCCESS) endif call MAPL_GetObjectFromGC ( GC, MAPL_OBJ, RC=STATUS) _VERIFY(STATUS) call MAPL_TimerOn(MAPL_OBJ,"Run") call MAPL_Get(MAPL_OBJ, & IM = IM, & JM = JM, & LM = LM, & LONS = LONS, & LATS = LATS, & RC=STATUS ) _VERIFY(STATUS) ! Figure out what type of grid we are on gridtype_default='Lat-Lon' call ESMF_InfoGetFromHost(Grid,infoh,rc=status) _VERIFY(STATUS) call ESMF_InfoGet(infoh,key='GridType',value=gridtype,default=gridtype_default,rc=status) _VERIFY(STATUS) ! Get the time interval, and start and end time ! timeinterval=timeinterval/2 ! IntervalTime=time-timeinterval IntervalTime=time call ESMF_TimeGet(IntervalTime,yy=iyr,mm=imm,dd=idd,h=ihr,m=imn,s=isc,_RC) call MAPL_PackTime(nymd,iyr,imm,idd) call MAPL_PackTime(nhms,ihr,imn,isc) interval_nymd(1)=nymd interval_nhms(1)=nhms IntervalTime=time+timeinterval call ESMF_TimeGet(IntervalTime,yy=iyr,mm=imm,dd=idd,h=ihr,m=imn,s=isc,_RC) call MAPL_PackTime(nymd,iyr,imm,idd) call MAPL_PackTime(nhms,ihr,imn,isc) interval_nymd(2)=nymd interval_nhms(2)=nhms ! set swath to zero for now swath=0. call MAPL_GridGet(GRID, globalCellCountPerDim=COUNTS, RC=STATUS) _VERIFY(STATUS) IM_world = counts(1) JM_world = counts(2) if (JM_world == 6*IM_world) then imsize = 4*IM_World else imsize = IM_World endif ! swatch(3) is actually the size of the length to use for interpolation if( imsize.le.200 ) call ESMF_ConfigGetAttribute(self%CF,swath(3),LABEL='INTERPOLATION_WIDTH:', DEFAULT= 10.0 ,RC=STATUS) if( imsize.gt.200 .and. & imsize.le.400 ) call ESMF_ConfigGetAttribute(self%CF,swath(3),LABEL='INTERPOLATION_WIDTH:', DEFAULT= 5.0 ,RC=STATUS) if( imsize.gt.400 .and. & imsize.le.800 ) call ESMF_ConfigGetAttribute(self%CF,swath(3),LABEL='INTERPOLATION_WIDTH:', DEFAULT= 2.0 ,RC=STATUS) if( imsize.gt.800 .and. & imsize.le.1600 ) call ESMF_ConfigGetAttribute(self%CF,swath(3),LABEL='INTERPOLATION_WIDTH:', DEFAULT= 1.0 ,RC=STATUS) if( imsize.gt.1600 ) call ESMF_ConfigGetAttribute(self%CF,swath(3),LABEL='INTERPOLATION_WIDTH:', DEFAULT= 0.5 ,RC=STATUS) call ESMF_ConfigDestroy( self%CF, rc=status) _VERIFY(STATUS) ! define undef undef=MAPL_UNDEF ! set deltat in seconds deltat=10 ! get orb bundle call ESMF_StateGet(EXPORT,'SATORB',BUNDLE,RC=STATUS) _VERIFY(STATUS) call ESMF_FieldBundleGet(BUNDLE,fieldCOUNT=NORB,RC=STATUS) _VERIFY(STATUS) _ASSERT(NORB==self%no,'needs informative message') ! loop over each satellite and get it's mask do k=1,NORB call ESMFL_BundleGetPointerToData(BUNDLE,trim(self%instrument(k)),PTR_TMP,RC=STATUS) _VERIFY(STATUS) call MAPL_GetPointer(EXPORT,PTR_TMP_EX,trim(self%instrument(k)),RC=STATUS) _VERIFY(STATUS) sat_name=trim(self%Satellite(k)) swath(1)=self%swath(k) swath(2)=self%swath(k) ihalo = self%halo(k) jhalo = self%halo(k) if (gridtype == 'Lat-Lon') then if (associated(PTR_TMP)) call DoMasking_ (PTR_TMP, im, jm, lons, lats, undef, & sat_name, interval_nymd, interval_nhms, deltat, swath, & ihalo, jhalo, rc=status ) else if (gridtype == 'Cubed-Sphere') then if (associated(PTR_TMP)) call DoMasking_CS (PTR_TMP, im, jm, self%x, self%y, undef, & sat_name, interval_nymd, interval_nhms, deltat, swath, & ihalo, jhalo, self%face, rc=status ) endif ! if HISTORY is asking for mask to write this will be allocated if (associated(PTR_TMP_EX)) PTR_TMP_EX=PTR_TMP if (associated(PTR_TMP)) nullify(PTR_TMP) if (associated(PTR_TMP_EX)) nullify(PTR_TMP_EX) enddo ! All done ! -------- call MAPL_TimerOff(MAPL_OBJ,"Run") _RETURN(ESMF_SUCCESS) END SUBROUTINE Run_ !....................................................................... subroutine extract_ ( GC, CLOCK, self, GRID, CF, time, nymd, nhms, timeinterval, rc) type(ESMF_GridComp), intent(INout) :: GC ! Grid Comp object type(ESMF_Clock), intent(in) :: CLOCK ! Clock type(Orb_state), pointer :: self ! Legacy state type(ESMF_Grid), intent(out) :: GRID ! Grid type(ESMF_Config), intent(out) :: CF ! Universal Config type(ESMF_TIME), intent(out) :: Time ! Time type(ESMF_TimeInterval), intent(out) :: TimeInterval ! Time Intervale integer, intent(out) :: nymd, nhms ! date, time integer, intent(out), optional :: rc ! --- character(len=ESMF_MAXSTR) :: comp_name _Iam_('extract_') type(Orb_Wrap) :: wrap integer :: iyr, imm, idd, ihr, imn, isc integer :: status ! Get my name and set-up traceback handle ! --------------------------------------- call ESMF_GridCompGet( GC, NAME=comp_name, _RC ) Iam = trim(COMP_NAME) // '::extract_' rc = 0 ! Get my internal state ! --------------------- call ESMF_UserCompGetInternalState(gc, 'Orb_state', WRAP, STATUS) _VERIFY(STATUS) self => wrap%ptr ! Get the configuration ! --------------------- call ESMF_GridCompGet ( GC, config=CF, _RC ) ! Extract time as simple integers from clock ! ------------------------------------------ call ESMF_ClockGet(CLOCK,currTIME=TIME,timeStep=TimeInterval,_RC) call ESMF_TimeGet(TIME ,yy=iyr, mm=imm, dd=idd, h=ihr, m=imn, s=isc, _RC) call MAPL_PackTime(nymd,iyr,imm,idd) call MAPL_PackTime(nhms,ihr,imn,isc) ! Extract the ESMF Grid ! --------------------- call ESMF_GridCompGet ( GC, grid=GRID, _RC) _RETURN(ESMF_SUCCESS) end subroutine extract_ !------------------------------------------------------------------------ !> ! The routine `DoMasking_` ... ! subroutine DoMasking_ (field, im, jm, lons, lats, undef, & sat_name, nymd, nhms, dt, swath, & ihalo, jhalo, rc ) use MAPL_NominalOrbitsMod implicit NONE ! !INPUT PARAMETERS: integer, intent(in) :: im, jm integer, intent(in) :: ihalo(2), jhalo(2) real, intent(in) :: lons(im,jm) real, intent(in) :: lats(im,jm) real, intent(in) :: swath(3) real, intent(in) :: undef character(len=*), intent(in) :: sat_name !! Satellite name integer, intent(in) :: dt !! delta t in secs integer, intent(in) :: nymd(2) !! Beginning/ending date: YYYYMMDD integer, intent(in) :: nhms(2) !! Beginning/ending time: HHMMSS ! !OUTPUT PARAMETERS: real, intent(inout) :: field(im,jm) integer, intent(out), optional :: rc !! Error return code !! = 0 all is well !! = 3 memory allocation error ! ---- INTEGER, PARAMETER :: dp=SELECTED_REAL_KIND(15,307) ! Local workspace ! --------------- integer :: isegs, jsegs, mask(im,jm) real(dp):: V_mean real(dp) :: SwathWidth(2) = (/ 0.0, 0.0 /) integer nobs real(dp), pointer :: tlons(:) => null() real(dp), pointer :: tlats(:) => null() real(dp), pointer :: slons(:,:) => null() real(dp), pointer :: slats(:,:) => null() SwathWidth(1:2) = swath(1:2) ! type conversion ! Interpoaltion parameters; swath(3) is the swath resolution, in km ! ----------------------------------------------------------------- V_mean = 7.5 ! Mean sat speed in km/s; assumes a 90 minute period ! this should be a function of satellite. isegs = nint((swath(1)+swath(2)) / swath(3) ) ! segments across-track jsegs = V_mean * dt / swath(3) ! segments along-track ! Simple masking without swath ! ---------------------------- if ( swath(1) .eq. 0 .AND. swath(2) .eq. 0 ) then ! Compute tracks ! -------------- call orbits_track(tlons,tlats, sat_name, nymd, nhms, dt, rc) ! Simple masking ! -------------- mask=0 call orb_mask_lonlat(mask,im,jm,lons,lats,tlons,tlats,size(tlons),jsegs,-180.,180.) deallocate(tlons,tlats) else ! Compute tracks ! -------------- call orbits_swath(slons,slats, sat_name, nymd, nhms, dt,SwathWidth, rc, wrap=.true.) ! Swath masking ! ------------- nobs = size(slons)/3 mask = 0 call orb_swath_mask_lonlat(mask,im,jm,lons,lats,slons,slats,nobs,isegs,jsegs,-180.,180.) deallocate(slons,slats) end if ! Add optional halo ! ----------------- if ( (sum(ihalo) + sum(jhalo)) > 0) then call orb_halo(im,jm,mask,ihalo,jhalo) end if ! Apply mask to field ! ------------------- field = 1.0 where ( mask /= 1 ) field = undef ! All done ! -------- end subroutine DoMasking_ !-------------------------------------------------------------------------------- !> ! The routine `DoMasking_CS` ... ! subroutine DoMasking_CS (field, im, jm, x, y, undef, & sat_name, nymd, nhms, dt, swath, & ihalo, jhalo, face, rc ) use MAPL_NominalOrbitsMod implicit NONE ! !INPUT PARAMETERS: integer, intent(in) :: im, jm integer, intent(in) :: ihalo(2), jhalo(2) real, intent(in) :: x(im+1,jm+1) real, intent(in) :: y(im+1,jm+1) real, intent(in) :: swath(3) real, intent(in) :: undef character(len=*), intent(in) :: sat_name !! Satellite name integer, intent(in) :: dt !! delta t in secs integer, intent(in) :: nymd(2) !! Beginning/ending date: YYYYMMDD integer, intent(in) :: nhms(2) !! Beginning/ending time: HHMMSS integer, intent(in) :: face ! !OUTPUT PARAMETERS: real, intent(inout) :: field(im,jm) integer, optional :: rc !! Error return code !! = 0 all is well !! = 3 memory allocation error ! ---- INTEGER, PARAMETER :: dp=SELECTED_REAL_KIND(15,307) ! Local workspace ! --------------- integer :: isegs, jsegs, mask(im,jm) real(dp):: V_mean real(dp) :: SwathWidth(2) = (/ 0.0, 0.0 /) integer nobs, status real(dp), pointer :: tlons(:) => null() real(dp), pointer :: tlats(:) => null() real(dp), pointer :: slons(:,:) => null() real(dp), pointer :: slats(:,:) => null() ! character(len=12) :: Iam="DoMasking_CS" SwathWidth(1:2) = swath(1:2) ! type conversion ! Interpoaltion parameters; swath(3) is the swath resolution, in km ! ----------------------------------------------------------------- V_mean = 7.5 ! Mean sat speed in km/s; assumes a 90 minute period ! this should be a function of satellite. isegs = nint((swath(1)+swath(2)) / swath(3) ) ! segments across-track jsegs = V_mean * dt / swath(3) ! segments along-track ! Simple masking without swath ! ---------------------------- if ( swath(1) .eq. 0 .AND. swath(2) .eq. 0 ) then mask = 0 ! Compute tracks ! -------------- call orbits_track(tlons,tlats, sat_name, nymd, nhms, dt, rc) ! Simple masking ! -------------- call orb_mask_xy(mask,im,jm,x,y,tlons,tlats,size(tlons),jsegs,0.,360.,face,rc=status) ! _VERIFY(STATUS) deallocate(tlons,tlats) else mask = 0 ! Compute tracks ! -------------- call orbits_swath(slons,slats, sat_name, nymd, nhms, dt,SwathWidth, rc, wrap=.true.) ! Swath masking ! ------------- nobs = size(slons)/3 call orb_swath_mask_xy(mask,im,jm,x,y,slons,slats,nobs,isegs,jsegs,0.,360.,face) deallocate(slons,slats) end if ! Add optional halo ! ----------------- if ( (sum(ihalo) + sum(jhalo)) > 0) then call orb_halo(im,jm,mask,ihalo,jhalo) end if ! Apply mask to field ! ------------------- field = 1.0 where ( mask /= 1 ) field = undef ! All done ! -------- end subroutine DoMasking_CS subroutine flatten_latlon(lats,lons,lats_1d,lons_1d,im,jm) implicit none integer, intent(in) :: im,jm real, intent(in) :: lats(im,jm), lons(im,jm) real, intent(inout) :: lats_1d(jm), lons_1d(im) integer :: i do i =1,im lons_1d(i)=lons(i,1)*180./MAPL_PI enddo do i =1,jm lats_1d(i)=lats(1,i)*180./MAPL_PI enddo return end subroutine flatten_latlon subroutine flatten_xy(x,y,x_1d,y_1d,im,jm,im_1d,jm_1d,switch) implicit none integer, intent(in) :: im,jm, im_1d, jm_1d real, intent(in) :: x(im,jm), y(im,jm) real, intent(inout) :: y_1d(jm_1d), x_1d(im_1d) logical, intent(in) :: switch integer :: i if (.not.switch) then do i =1,im_1d x_1d(i)=x(i,1) enddo do i =1,jm_1d y_1d(i)=y(1,i) enddo else if (switch) then do i =1,im_1d x_1d(i)=x(1,i) enddo do i =1,jm_1d y_1d(i)=y(i,1) enddo endif return end subroutine flatten_xy subroutine orb_edges_1d(ecoords,coords,idim) implicit NONE integer, intent(in) :: idim real, intent(in) :: coords(idim) real, intent(out) :: ecoords(idim+1) ecoords(1) = coords(1) - 0.5 * ( coords(2) - coords(1) ) ecoords(2:idim) = 0.5 * ( coords(1:idim-1)+coords(2:idim) ) ecoords(idim+1) = coords(idim) + 0.5 * (coords(idim) - coords(idim-1)) return end subroutine orb_edges_1d subroutine orb_mask_xy(mask,im,jm,x,y,tlons,tlats,nobs,jsegs,lb,ub,face,rc) implicit NONE INTEGER, PARAMETER :: dp=SELECTED_REAL_KIND(15,307) integer, intent(in) :: im, jm, nobs, jsegs real, intent(in) :: x(im+1,jm+1) real, intent(in) :: y(im+1,jm+1) real, intent(in) :: lb,ub integer, intent(in) :: face integer, optional, intent(out) :: rc real, pointer :: ex(:), ey(:) real(dp) :: tlons(nobs), tlats(nobs) real(dp) :: beta integer, intent(out) :: mask(im,jm) real :: x_loc, y_loc integer i, j, m, n, nfail, inbox, imp1, jmp1, face_pnt integer im_1d, jm_1d, itmp logical switch real :: wcorner_x(4),wcorner_y(4) ! corners of world for this proc real :: lat,lon ! _Iam_('orb_mask_xy') if (present(rc)) then rc = ESMF_SUCCESS end if switch = .false. if ( abs(x(1,1)-x(2,1)) < abs(x(1,1)-x(1,2)) ) switch = .true. if (.not.switch) then allocate(ex(im+1),ey(jm+1)) im_1d = im+1 jm_1d = jm+1 imp1 = im+1 jmp1 = jm+1 endif if (switch) then allocate(ex(jm+1),ey(im+1)) im_1d = jm+1 jm_1d = im+1 imp1 = im+1 jmp1 = jm+1 end if call flatten_xy(x,y,ex,ey,imp1,jmp1,im_1d,jm_1d,switch) nfail = 0 ! loop over each point ! first define corners of the world wcorner_x(1)=minval(ex) wcorner_x(4)=minval(ex) wcorner_x(2)=maxval(ex) wcorner_x(3)=maxval(ex) wcorner_y(1)=minval(ey) wcorner_y(4)=maxval(ey) wcorner_y(2)=minval(ey) wcorner_y(3)=maxval(ey) do n = 2, nobs do m = 1, jsegs ! interpolate track along orbit for better resoluation between points beta = (m - 1.0 ) / ( jsegs - 1.0 ) if ( abs(tlons(n)-tlons(n-1))<180. ) then lon = (1.0-beta) * tlons(n-1) + beta * tlons(n) else if ( tlons(n) > tlons(n-1) ) then lon = (1.0-beta) * tlons(n-1) + beta * (tlons(n)-360.) else lon = (1.0-beta) * (tlons(n-1)-360.) + beta * tlons(n) end if if (lon < lb) lon=lon+360. if (lon > ub) lon=lon-360. lat = (1.0-beta) * tlats(n-1) + beta * tlats(n) ! check if lat and lon is in the piece of world on this processor lat = lat * MAPL_PI/180.0 lon = lon * MAPL_PI/180.0 call check_face_pnt(LON,LAT,face_pnt) if (face_pnt == face) then call cube_xy_point(x_loc,y_loc,LAT,LON,face_pnt) inbox = pnt_in_rect(x_loc,y_loc,wcorner_x,wcorner_y) ! if we found the point in the box if (inbox == 1) then i = ijsearch(ex,im_1d,x_loc,.false.) j = ijsearch(ey,jm_1d,y_loc,.false.) ! fill in mask for i,j if (switch) then itmp = i i = j j = itmp endif if ( i>0 .and. j>0 .and. i<=im .and. j<=jm) then mask(i,j)=1 end if endif endif enddo enddo deallocate(ex,ey) end subroutine orb_mask_xy subroutine orb_mask_lonlat(mask,im,jm,lons,lats,tlons,tlats,nobs,jsegs,lb,ub) implicit NONE INTEGER, PARAMETER :: dp=SELECTED_REAL_KIND(15,307) integer, intent(in) :: im, jm, nobs, jsegs real, intent(in) :: lons(im,jm) real, intent(in) :: lats(im,jm) real, intent(in) :: lb,ub real :: lons_1d(im), lats_1d(jm) real :: elons(im+1), elats(jm+1) real(dp) :: tlons(nobs), tlats(nobs) real(dp) :: beta integer, intent(out) :: mask(im,jm) integer i, j, m, n, inbox, imp1, jmp1 real :: wcorner_lat(4),wcorner_lon(4) ! corners of world for this proc real :: lat,lon ! Build edge coords ! ----------------- call flatten_latlon(lats,lons,lats_1d,lons_1d,im,jm) call orb_edges_1d(elons,lons_1d,im) call orb_edges_1d(elats,lats_1d,jm) ! since we will need these imp1=im+1 jmp1=jm+1 ! first define corners of the world wcorner_lat(1)=elats(1) wcorner_lat(4)=elats(jmp1) wcorner_lat(2)=elats(1) wcorner_lat(3)=elats(jmp1) wcorner_lon(1)=elons(1) wcorner_lon(4)=elons(1) wcorner_lon(2)=elons(imp1) wcorner_lon(3)=elons(imp1) do n = 2, nobs do m = 1, jsegs ! interpolate track along orbit for better resoluation between points beta = (m - 1.0 ) / ( jsegs - 1.0 ) if ( abs(tlons(n)-tlons(n-1))<180. ) then lon = (1.0-beta) * tlons(n-1) + beta * tlons(n) else if ( tlons(n) > tlons(n-1) ) then lon = (1.0-beta) * tlons(n-1) + beta * (tlons(n)-360.) else lon = (1.0-beta) * (tlons(n-1)-360.) + beta * tlons(n) end if if (lon < lb) lon=lon+360. if (lon > ub) lon=lon-360. lat = (1.0-beta) * tlats(n-1) + beta * tlats(n) ! check if lat and lon is in the piece of world on this processor inbox = pnt_in_rect(lat,lon,wcorner_lat,wcorner_lon) ! if we found the point in the box if (inbox == 1) then i = ijsearch(elons,im+1,lon,.false.) j = ijsearch(elats,jm+1,lat,.false.) ! fill in mask for i,j if (i >0 .and. j > 0) then mask(i,j)=1 endif endif enddo enddo end subroutine orb_mask_lonlat subroutine orb_swath_mask_xy(mask,im,jm,x,y,slons,slats,nobs,isegs,jsegs,lb,ub,face) implicit NONE INTEGER, PARAMETER :: dp=SELECTED_REAL_KIND(15,307) integer, intent(in) :: im, jm, nobs, isegs, jsegs real, intent(in) :: x(im+1,jm+1) real, intent(in) :: y(im+1,jm+1) real, intent(in) :: lb,ub integer, intent(in) :: face real, pointer :: ex(:), ey(:) real(dp) :: slons(3,nobs), slats(3,nobs) real(dp) :: alpha, beta, lon1, lon2, lat1, lat2 integer, intent(out) :: mask(im,jm) real :: x_loc, y_loc integer :: i, j, k, m, n, imp1, jmp1, inbox, itmp integer :: im_1d, jm_1d real :: wcorner_x(4),wcorner_y(4) real :: lat,lon real :: sdnom1,sdnom2,eplonl1,eplonl2,eplonr1,eplonr2 real :: eplatl1,eplatl2,eplatr1,eplatr2 real :: sp1,sp2 real :: sdnom,eplon1,eplon2,eplat1,eplat2 real :: latf integer :: face_pnt logical :: switch ! find indices have constant values of coordinate switch = .false. if ( abs(x(1,1)-x(2,1)) < abs(x(1,1)-x(1,2)) ) switch = .true. if (.not.switch) then allocate(ex(im+1),ey(jm+1)) im_1d = im+1 jm_1d = jm+1 imp1=im+1 jmp1=jm+1 endif if (switch) then allocate(ex(jm+1),ey(im+1)) im_1d = jm+1 jm_1d = im+1 imp1=im+1 jmp1=jm+1 end if call flatten_xy(x,y,ex,ey,imp1,jmp1,im_1d,jm_1d,switch) ! first define corners of the world wcorner_x(1)=minval(ex) wcorner_x(4)=minval(ex) wcorner_x(2)=maxval(ex) wcorner_x(3)=maxval(ex) wcorner_y(1)=minval(ey) wcorner_y(4)=maxval(ey) wcorner_y(2)=minval(ey) wcorner_y(3)=maxval(ey) do k = 1, isegs ! cros-track alpha = (k - 1.0 ) / ( isegs - 1.0 ) do n = 2, nobs if (abs(slons(1,n-1)-slons(3,n-1)) < 180.) then lon1 = (1.0-alpha) * slons(1,n-1) + alpha * slons(3,n-1) eplonl1 = slons(1,n-1) eplonr1 = slons(3,n-1) else if (slons(1,n-1) > slons(3,n-1)) then lon1 = (1.0-alpha) * slons(1,n-1) + alpha * (slons(3,n-1)+360.) eplonl1 = slons(1,n-1) eplonr1 = slons(3,n-1)+360. else lon1 = (1.0-alpha) * (slons(1,n-1)+360.) + alpha * slons(3,n-1) eplonl1 = slons(1,n-1)+360. eplonr1 = slons(3,n-1) endif if (abs(slons(1,n)-slons(3,n)) < 180.) then lon2 = (1.0-alpha) * slons(1,n) + alpha * slons(3,n) eplonl2 = slons(1,n) eplonr2 = slons(3,n) else if (slons(1,n) > slons(3,n)) then lon2 = (1.0-alpha) * slons(1,n) + alpha * (slons(3,n)+360.) eplonl2 = slons(1,n) eplonr2 = slons(3,n)+360. else lon2 = (1.0-alpha) * (slons(1,n)+360.) + alpha * slons(3,n) eplonl2 = slons(1,n)+360. eplonr2 = slons(3,n) endif ! interpolate along great circle unless endpoints of interpolation have same lon associate(d2r => MAPL_DEGREES_TO_RADIANS_R8, r2d => MAPL_RADIANS_TO_DEGREES) eplatl1 = slats(1,n-1) eplatr1 = slats(3,n-1) eplatl2 = slats(1,n) eplatr2 = slats(3,n) sdnom1 = sin((eplonl1-eplonr1)*d2r) sdnom2 = sin((eplonl2-eplonr2)*d2r) if (abs(sdnom1) /= 0.) then sp1 = sin((lon1-eplonr1)*d2r)/sdnom1 sp2 = sin((lon1-eplonl1)*d2r)/sdnom1 lat1 = atan(tan(eplatl1*d2r)*sp1 - tan(eplatr1*d2r)*sp2) lat1 = lat1*r2d else lat1 = (1.0-alpha) * slats(1,n-1) + alpha * slats(3,n-1) endif if (abs(sdnom2) /= 0.) then sp1 = sin((lon2-eplonr2)*d2r)/sdnom2 sp2 = sin((lon2-eplonl2)*d2r)/sdnom2 lat2 = atan(tan(eplatl2*d2r)*sp1 - tan(eplatr2*d2r)*sp2) lat2 = lat2*r2d else lat2 = (1.0-alpha) * slats(1,n) + alpha * slats(3,n) endif do m = 1, jsegs ! along track refinement beta = (m - 1.0 ) / ( jsegs - 1.0 ) if (abs(lon2-lon1) < 180.) then lon = (1.0-beta) * lon1 + beta * lon2 eplon1=lon1 eplon2=lon2 else if (lon2 > lon1) then lon = (1.0-beta) * (lon1+360.) + beta * lon2 eplon1=lon1+360. eplon2=lon2 else lon = (1.0-beta) * lon1 + beta * (lon2+360.) eplon1=lon1 eplon2=lon2+360. endif eplat1=lat1 eplat2=lat2 sdnom=sin((eplon1-eplon2)*d2r) if (abs(sdnom) /= 0. ) then sp1=sin((lon-eplon2)*d2r)/sdnom sp2=sin((lon-eplon1)*d2r)/sdnom latf = atan(tan(eplat1*d2r)*sp1-tan(eplat2*d2r)*sp2) latf = latf*r2d lat = latf else lat = (1.0-beta) * lat1 + beta * lat2 endif if (lon < lb) lon=lon+360. if (lon > ub) lon=lon-360. lat = lat * MAPL_PI/180.0 lon = lon * MAPL_PI/180.0 call check_face_pnt(LON,LAT,face_pnt) if (face_pnt == face) then call cube_xy_point(x_loc,y_loc,LAT,LON,face) inbox = pnt_in_rect(x_loc,y_loc,wcorner_x,wcorner_y) if (inbox == 1) then i = ijsearch(ex,im_1d,x_loc,.false.) j = ijsearch(ey,jm_1d,y_loc,.false.) if (switch) then itmp = i i = j j = itmp endif if ( i>0 .and. j>0 .and. i<=im .and. j<=jm) then mask(i,j)=1 end if end if endif end do ! msegs end associate end do ! nobs end do ! ksegs deallocate(ex,ey) end subroutine orb_swath_mask_xy !........................................................................................... subroutine orb_swath_mask_lonlat(mask,im,jm,lons,lats,slons,slats,nobs,isegs,jsegs,lb,ub) implicit NONE INTEGER, PARAMETER :: dp=SELECTED_REAL_KIND(15,307) integer, intent(in) :: im, jm, nobs, isegs, jsegs real, intent(in) :: lons(im,jm) real, intent(in) :: lats(im,jm) real, intent(in) :: lb,ub real :: lons_1d(im),lats_1d(jm) real :: elons(im+1), elats(jm+1) real(dp) :: slons(3,nobs), slats(3,nobs) real(dp) :: alpha, beta, lon1, lon2, lat1, lat2 integer, intent(out) :: mask(im,jm) integer i, j, k, m, n, nfail, imp1, jmp1, inbox real :: wcorner_lat(4),wcorner_lon(4) real :: lat,lon real :: sdnom1,sdnom2,eplonl1,eplonl2,eplonr1,eplonr2 real :: eplatl1,eplatl2,eplatr1,eplatr2 real :: sp1,sp2 real :: sdnom,eplon1,eplon2,eplat1,eplat2 real :: latf ! Build edge coords ! ----------------- call flatten_latlon(lats,lons,lats_1d,lons_1d,im,jm) call orb_edges_1d(elons,lons_1d,im) call orb_edges_1d(elats,lats_1d,jm) imp1=im+1 jmp1=jm+1 nfail = 0 ! first define corners of the world wcorner_lat(1)=elats(1) wcorner_lat(4)=elats(jmp1) wcorner_lat(2)=elats(1) wcorner_lat(3)=elats(jmp1) wcorner_lon(1)=elons(1) wcorner_lon(4)=elons(1) wcorner_lon(2)=elons(imp1) wcorner_lon(3)=elons(imp1) do k = 1, isegs ! cros-track alpha = (k - 1.0 ) / ( isegs - 1.0 ) do n = 2, nobs if (abs(slons(1,n-1)-slons(3,n-1)) < 180.) then lon1 = (1.0-alpha) * slons(1,n-1) + alpha * slons(3,n-1) eplonl1 = slons(1,n-1) eplonr1 = slons(3,n-1) else if (slons(1,n-1) > slons(3,n-1)) then lon1 = (1.0-alpha) * slons(1,n-1) + alpha * (slons(3,n-1)+360.) eplonl1 = slons(1,n-1) eplonr1 = slons(3,n-1)+360. else lon1 = (1.0-alpha) * (slons(1,n-1)+360.) + alpha * slons(3,n-1) eplonl1 = slons(1,n-1)+360. eplonr1 = slons(3,n-1) endif if (abs(slons(1,n)-slons(3,n)) < 180.) then lon2 = (1.0-alpha) * slons(1,n) + alpha * slons(3,n) eplonl2 = slons(1,n) eplonr2 = slons(3,n) else if (slons(1,n) > slons(3,n)) then lon2 = (1.0-alpha) * slons(1,n) + alpha * (slons(3,n)+360.) eplonl2 = slons(1,n) eplonr2 = slons(3,n)+360. else lon2 = (1.0-alpha) * (slons(1,n)+360.) + alpha * slons(3,n) eplonl2 = slons(1,n)+360. eplonr2 = slons(3,n) endif ! interpolate along great circle unless endpoints of interpolation have same lon associate(d2r => MAPL_DEGREES_TO_RADIANS_R8, r2d => MAPL_RADIANS_TO_DEGREES) eplatl1 = slats(1,n-1) eplatr1 = slats(3,n-1) eplatl2 = slats(1,n) eplatr2 = slats(3,n) sdnom1 = sin((eplonl1-eplonr1)*d2r) sdnom2 = sin((eplonl2-eplonr2)*d2r) if (abs(sdnom1) /= 0.) then sp1 = sin((lon1-eplonr1)*d2r)/sdnom1 sp2 = sin((lon1-eplonl1)*d2r)/sdnom1 lat1 = atan(tan(eplatl1*d2r)*sp1 - tan(eplatr1*d2r)*sp2) lat1 = lat1*r2d else lat1 = (1.0-alpha) * slats(1,n-1) + alpha * slats(3,n-1) endif if (abs(sdnom2) /= 0.) then sp1 = sin((lon2-eplonr2)*d2r)/sdnom2 sp2 = sin((lon2-eplonl2)*d2r)/sdnom2 lat2 = atan(tan(eplatl2*d2r)*sp1 - tan(eplatr2*d2r)*sp2) lat2 = lat2*r2d else lat2 = (1.0-alpha) * slats(1,n) + alpha * slats(3,n) endif do m = 1, jsegs ! along track refinement beta = (m - 1.0 ) / ( jsegs - 1.0 ) if (abs(lon2-lon1) < 180.) then lon = (1.0-beta) * lon1 + beta * lon2 eplon1=lon1 eplon2=lon2 else if (lon2 > lon1) then lon = (1.0-beta) * (lon1+360.) + beta * lon2 eplon1=lon1+360. eplon2=lon2 else lon = (1.0-beta) * lon1 + beta * (lon2+360.) eplon1=lon1 eplon2=lon2+360. endif eplat1=lat1 eplat2=lat2 sdnom=sin((eplon1-eplon2)*d2r) if (abs(sdnom) /= 0. ) then sp1=sin((lon-eplon2)*d2r)/sdnom sp2=sin((lon-eplon1)*d2r)/sdnom latf = atan(tan(eplat1*d2r)*sp1-tan(eplat2*d2r)*sp2) latf = latf*r2d lat = latf else lat = (1.0-beta) * lat1 + beta * lat2 endif if (lon < lb) lon=lon+360. if (lon > ub) lon=lon-360. inbox = pnt_in_rect(lat,lon,wcorner_lat,wcorner_lon) if (inbox == 1) then i = ijsearch(elons,im+1,lon,.false.) j = ijsearch(elats,jm+1,lat,.false.) if ( i>0 .and. i<=im .and. j>0 .and. j<=jm ) mask(i,j) = 1 end if end do ! msegs end associate end do ! nobs end do ! ksegs end subroutine orb_swath_mask_lonlat integer function pnt_in_rect(x0,y0,x,y) implicit none real, intent(in) :: x0,y0 real, intent(in), dimension(4) :: x,y ! local variable integer :: i, j, i1,i2, nintersect real :: s,t,x1,y1, denom, tol real, dimension(2) :: p0,p1,q0,q1,w,v,u tol=0. ! for now set x1,y1 x1=x0+0.2 y1=y0+0.2 q0(1)=x0 q0(2)=y0 q1(1)=x1 q1(2)=y1 nintersect = 0 ! loop over four sides do i = 1, 4 i1 = i i2 = i+1-(4*(i/4)) ! define edge p0(1)=x(i1) p0(2)=y(i1) p1(1)=x(i2) p1(2)=y(i2) ! define vectors w,v,u do j=1,2 w(j)=p0(j)-q0(j) v(j)=q1(j)-q0(j) u(j)=p1(j)-p0(j) enddo denom = v(1)*u(2)-v(2)*u(1) ! if the ray and edge are parallel they probably don't cross, should really ! check if they are on top of each other but probably won't happen if (abs(denom).gt.tol) then s = (v(2)*w(1)-v(1)*w(2))/denom t = -(u(1)*w(2)-u(2)*w(1))/denom ! now test if ray intersects segment t > 0 and 0 < s < 1 if (0. <= s .and. s <= 1. .and. t >= 0.) nintersect = nintersect + 1 ! in future put in test if t = 0, i.e. the point is on one of the edges endif enddo if (mod(nintersect,2) == 1) then ! if (nintersect == 1) then pnt_in_rect = 1 else pnt_in_rect = 0 endif end function pnt_in_rect integer function ijsearch(coords,idim,value,periodic) ! fast bisection version implicit NONE integer, intent(in) :: idim real, intent(in) :: coords(idim) real, intent(inout) :: value logical, intent(in) :: periodic integer i, i1, i2, k if ( periodic ) then if ( value>coords(idim) ) value = value - 360. endif if ( value .eq. coords(idim) ) then ijsearch = idim return endif ijsearch = -1 i1 = 1 i2 = idim if (coords(idim) > coords(1)) then do k = 1, idim ! it should never take take long i = (i1 + i2) / 2 if ( (value .ge. coords(i)) ) then if (value .lt. coords(i+1) ) then ijsearch = i exit else i1 = i end if else i2 = i endif end do else do k = 1, idim ! it should never take take long i = (i1 + i2) / 2 if ( (value .lt. coords(i)) ) then if (value .ge. coords(i+1) ) then ijsearch = i exit else i1 = i end if else i2 = i endif end do endif end function subroutine check_face(IM,JM,LONS,LATS,face) integer, intent(in) :: im,jm real(ESMF_KIND_R8), intent(in) :: LONS(IM,JM),LATS(IM,JM) integer, intent(inout) :: face real :: lon,lat integer :: i,j,k real :: s(6),smin, xyz(3), rsq3 integer :: fmin,ifmin(6),imax,imaxt rsq3=1.0/sqrt(3.) do i=1,IM do j=1,JM smin = 30.0 lon=LONS(i,j)+MAPL_PI/18.0 lat=LATS(i,j) xyz(1)=cos(lon)*cos(lat) xyz(2)=sin(lon)*cos(lat) xyz(3)=sin(lat) if (xyz(1) /= 0.0) then s(1)=rsq3/xyz(1) s(2)=-rsq3/xyz(1) else s(1)=1000.0 s(2)=1000.0 endif if (xyz(2) /= 0.0) then s(3)=rsq3/xyz(2) s(4)=-rsq3/xyz(2) else s(3)=1000.0 s(4)=1000.0 endif if (xyz(3) /= 0.0) then s(5)=rsq3/xyz(3) s(6)=-rsq3/xyz(3) else s(5)=1000.0 s(6)=1000.0 endif do k=1,6 if (s(k) > 0) then if (s(k) < smin) then smin = s(k) fmin = k endif endif enddo ifmin(fmin) = ifmin(fmin)+1 enddo enddo imax = 0 do k=1,6 imaxt=ifmin(k) if (imaxt > imax) then imax = imaxt face = k endif enddo end subroutine check_face subroutine cube_xy(IM,JM,x,y,LONS,LATS,face) integer, intent(in) :: IM,JM real, intent(inout) :: x(IM,JM),y(IM,JM) real(ESMF_KIND_R8), intent(in) :: LATS(IM,JM),LONS(IM,JM) integer, intent(in) :: face real :: rsq3,LAT,LON integer :: i,j rsq3 = 1.0/sqrt(3.0) do i=1,IM do j=1,JM LAT=LATS(I,J) LON=LONS(I,J)+MAPL_PI/18.0 select case(face) case (1) x(I,J) = rsq3*tan(LON) y(I,J) = rsq3*tan(LAT)/cos(LON) case (2) x(I,J) = rsq3*tan(LON) y(I,J) = -rsq3*tan(LAT)/cos(LON) case (3) x(I,J) = -rsq3*cos(LON)/sin(LON) y(I,J) = rsq3*tan(LAT)/sin(LON) case (4) x(I,J) = -rsq3*cos(LON)/sin(LON) y(I,J) = -rsq3*tan(LAT)/sin(LON) case (5) x(I,J) = rsq3*sin(LON)*cos(LAT)/sin(LAT) y(I,J) = -rsq3*cos(LON)*cos(LAT)/sin(LAT) case (6) x(I,J) = -rsq3*sin(LON)*cos(LAT)/sin(LAT) y(I,J) = -rsq3*cos(LON)*cos(LAT)/sin(LAT) end select enddo enddo end subroutine cube_xy subroutine cube_xy_point(x,y,LAT,LON,face) real, intent(inout) :: x,y real, intent(in) :: LAT,LON integer, intent(in) :: face real :: rsq3,llon,llat rsq3 = 1.0/sqrt(3.0) LLAT=LAT LLON=LON+MAPL_PI/18.0 select case(face) case (1) x = rsq3*tan(LLON) y = rsq3*tan(LLAT)/cos(LLON) case (2) x = rsq3*tan(LLON) y = -rsq3*tan(LLAT)/cos(LLON) case (3) x = -rsq3*cos(LLON)/sin(LLON) y = rsq3*tan(LLAT)/sin(LLON) case (4) x = -rsq3*cos(LLON)/sin(LLON) y = -rsq3*tan(LLAT)/sin(LLON) case (5) x = rsq3*sin(LLON)*cos(LLAT)/sin(LLAT) y = -rsq3*cos(LLON)*cos(LLAT)/sin(LLAT) case (6) x = -rsq3*sin(LLON)*cos(LLAT)/sin(LLAT) y = -rsq3*cos(LLON)*cos(LLAT)/sin(LLAT) end select end subroutine cube_xy_point subroutine check_face_pnt(LON,LAT,face) real, intent(in) :: LON,LAT integer, intent(inout) :: face real :: llon,llat integer :: k real :: s(6),smin, xyz(3), rsq3 integer :: fmin character(len=ESMF_MAXSTR) :: Iam Iam = 'check_face_pnt' rsq3=1.0/sqrt(3.) smin = 30.0 llon=lon+MAPL_PI/18.0 llat=lat xyz(1)=cos(llon)*cos(llat) xyz(2)=sin(llon)*cos(llat) xyz(3)=sin(llat) fmin = 7 if (xyz(1) /= 0.) then s(1)=rsq3/xyz(1) s(2)=-rsq3/xyz(1) else s(1)=1000. s(2)=1000. endif if (xyz(2) /= 0.) then s(3)=rsq3/xyz(2) s(4)=-rsq3/xyz(2) else s(3)=1000. s(4)=1000. endif if (xyz(3) /= 0.) then s(5)=rsq3/xyz(3) s(6)=-rsq3/xyz(3) else s(5)=1000. s(6)=1000. endif do k=1,6 if (s(k) > 0) then if (s(k) < smin) then smin = s(k) fmin = k endif endif enddo if (fmin /= 7) then face = fmin endif end subroutine check_face_pnt subroutine orb_halo(im,jm,mask,ihalo,jhalo,rc) integer, intent(in ) :: im integer, intent(in ) :: jm integer, intent(inout) :: mask(im,jm) integer, intent(in ) :: ihalo(2) integer, intent(in ) :: jhalo(2) integer, optional, intent(out ) :: rc character(len=ESMF_MAXSTR) :: Iam integer :: i, j, is, js integer :: tmask(im,jm) Iam = "orb_halo" tmask = 0 do i = 1, im do j = 1, jm if (mask(i,j) == 1) then do is = i-ihalo(1),i+ihalo(2) if (is >= 1 .and. is <= im) then do js = j-jhalo(1),j+jhalo(2) if (js >= 1 .and. js <=jm) tmask(is,js) = 1 end do ! js loop end if ! is in range end do ! is loop end if ! (i,j) has mask = 1 end do ! j loop end do ! i loop mask = tmask _RETURN(ESMF_SUCCESS) end subroutine orb_halo end module MAPL_OrbGridCompMod