#include "MAPL_Exceptions.h" #include "MAPL_ErrLog.h" #include "MAPL_Generic.h" module MAPL_ExtDataMask use ESMF use MAPL_KeywordEnforcerMod use ESMFL_Mod use MAPL_BaseMod use MAPL_ExceptionHandling use gFTL_StringVector use MAPL_NewArthParserMod use MAPL_Constants implicit none private type, public :: ExtDataMask character(len=:), allocatable :: mask_type character(len=:), allocatable :: mask_arguments contains procedure :: get_mask_variables procedure :: evaluate_mask procedure :: evaluate_region_mask procedure :: evaluate_zone_mask procedure :: evaluate_box_mask end type ExtDataMask interface ExtDataMask module procedure new_ExtDataMask end interface ExtDataMask contains function new_ExtDataMask(mask_expression,rc) result(new_mask) type(ExtDataMask) :: new_mask character(len=*), intent(in) :: mask_expression integer, optional, intent(out) :: rc integer :: status character(len=:), allocatable :: function_name character(len=:), allocatable :: arguments integer :: i1,len i1 = index(mask_expression,"(") _ASSERT(i1 > 0,'Incorrect format for function expression: missing "("') function_name = adjustl(mask_expression(:i1-1)) function_name = ESMF_UtilStringLowerCase(function_name, _RC) if (index(function_name,"regionmask") /= 0) then new_mask%mask_type = "regionmask" else if (index(function_name,"zonemask") /= 0) then new_mask%mask_type = "zonemask" else if (index(function_name,"boxmask") /= 0) then new_mask%mask_type = "boxmask" else _FAIL("Invalid mask type") end if len = len_trim(mask_expression) arguments = adjustl(mask_expression(i1+1:len-1)) i1 = index(arguments,",") _ASSERT(i1 > 0,'Incorrect format for function expression: missing ","') new_mask%mask_arguments = arguments _RETURN(_SUCCESS) end function function get_mask_variables(this,rc) result(variables_in_mask) class(ExtDataMask), intent(inout) :: this type(StringVector) :: variables_in_mask integer, intent(out), optional :: rc integer :: i1,i2 logical :: twovar character(len=:), allocatable :: tmpstring1,tmpstring2 if (this%mask_type == "regionmask") twovar = .true. if (this%mask_type == "zonemask") twovar = .false. if (this%mask_type == "boxmask") twovar = .false. i1 = index(this%mask_arguments,",") i2 = index(this%mask_arguments,";") if (twovar) then tmpstring1 = this%mask_arguments(1:i1-1) variables_in_mask = parser_variables_in_expression(tmpstring1) tmpstring2 = this%mask_arguments(i1+1:i2-1) call variables_in_mask%push_back(trim(tmpstring2)) else tmpstring1 = this%mask_arguments(1:i1-1) variables_in_mask = parser_variables_in_expression(tmpstring1) end if _RETURN(_SUCCESS) end function subroutine evaluate_mask(this,state,var_name,rc) class(ExtDataMask), intent(inout) :: this type(ESMF_State), intent(inout) :: state character(len=*), intent(in) :: var_name integer, optional, intent(out) :: rc integer :: status select case(this%mask_type) case("regionmask") call this%evaluate_region_mask(state,var_name,_RC) case("zonemask") call this%evaluate_zone_mask(state,var_name,_RC) case("boxmask") call this%evaluate_box_mask(state,var_name,_RC) end select _RETURN(_SUCCESS) end subroutine evaluate_mask subroutine evaluate_region_mask(this,state,var_name,rc) class(ExtDataMask), intent(inout) :: this type(ESMF_State), intent(inout) :: state character(len=*), intent(in) :: var_name integer, optional, intent(out) :: rc integer :: status integer :: k,i character(len=:), allocatable :: maskString,maskname,vartomask integer, allocatable :: regionNumbers(:), flag(:) integer, allocatable :: mask(:,:) real, pointer :: rmask(:,:) real, pointer :: rvar2d(:,:) real, pointer :: rvar3d(:,:,:) real, pointer :: var2d(:,:) real, pointer :: var3d(:,:,:) integer :: rank,ib,ie type(ESMF_Field) :: field,temp_field call ESMF_StateGet(state,var_name,field,_RC) call ESMF_FieldGet(field,rank=rank,_RC) temp_field = create_field_from_Field(field,_RC) ! get mask string ib = index(this%mask_arguments,";") maskString = this%mask_arguments(ib+1:) ! get mask name ie = index(this%mask_arguments,";") ib = index(this%mask_arguments,",") vartomask = this%mask_arguments(:ib-1) maskname = this%mask_arguments(ib+1:ie-1) call MAPL_StateEval(state,vartomask,temp_field,_RC) call MAPL_GetPointer(state,rmask,maskName,_RC) if (rank == 2) then !call MAPL_GetPointer(state,rvar2d,vartomask,_RC) call ESMF_FieldGet(temp_field,0,farrayptr=rvar2d,_RC) call MAPL_GetPointer(state,var2d,var_name,_RC) else if (rank == 3) then !call MAPL_GetPointer(state,rvar3d,vartomask,_RC) call ESMF_FieldGet(temp_field,0,farrayptr=rvar3d,_RC) call MAPL_GetPointer(state,var3d,var_name,_RC) else _FAIL('Rank must be 2 or 3') end if k=32 allocate(regionNumbers(k), flag(k), stat=status) _VERIFY(STATUS) regionNumbers = 0 call ExtDataExtractIntegers(maskString,k,regionNumbers,rc=status) _VERIFY(STATUS) flag(:) = 1 WHERE(regionNumbers(:) == 0) flag(:) = 0 k = SUM(flag) deallocate(flag,stat=status) _VERIFY(STATUS) ! Set local mask to 1 where gridMask matches each integer (within precision!) ! --------------------------------------------------------------------------- allocate(mask(size(rmask,1),size(rmask,2)),stat=status) _VERIFY(STATUS) mask = 0 DO i=1,k WHERE(regionNumbers(i)-0.01 <= rmask .AND. & rmask <= regionNumbers(i)+0.01) mask = 1 END DO if (rank == 2) then var2d = rvar2d where(mask == 0) var2d = 0.0 else if (rank == 3) then var3d = rvar3d do i=1,size(var3d,3) where(mask == 0) var3d(:,:,i) = 0.0 enddo end if deallocate( mask) call ESMF_FieldDestroy(temp_field, noGarbage=.true., _RC) _RETURN(_SUCCESS) end subroutine evaluate_region_mask subroutine evaluate_zone_mask(this,state,var_name,rc) class(ExtDataMask), intent(inout) :: this type(ESMF_State), intent(inout) :: state character(len=*), intent(in) :: var_name integer, optional, intent(out) :: rc integer :: status integer :: i character(len=:), allocatable :: vartomask,clatS,clatN real, pointer :: rvar2d(:,:) real, pointer :: rvar3d(:,:,:) real, pointer :: var2d(:,:) real, pointer :: var3d(:,:,:) real(REAL64), pointer :: lats(:,:) real(REAL64) :: limitS, limitN type(ESMF_Field) :: field,temp_field type(ESMF_Grid) :: grid integer :: rank,ib,is type(ESMF_CoordSys_Flag) :: coordSys call ESMF_StateGet(state,var_name,field,_RC) call ESMF_FieldGet(field,rank=rank,grid=grid,_RC) temp_field = create_field_from_Field(field,_RC) ib = index(this%mask_arguments,",") vartomask = this%mask_arguments(:ib-1) is = index(this%mask_arguments,",",back=.true.) clatS = this%mask_arguments(ib+1:is-1) clatN = this%mask_arguments(is+1:) READ(clatS,*,IOSTAT=status) limitS _VERIFY(status) READ(clatN,*,IOSTAT=status) limitN _VERIFY(status) call ESMF_GridGetCoord(grid, coordDim=2, localDE=0, & staggerloc=ESMF_STAGGERLOC_CENTER, farrayPtr=lats, rc=status) _VERIFY(status) call ESMF_GridGet(grid,coordsys=coordsys,_RC) if (coordsys == ESMF_COORDSYS_SPH_RAD) then limitN=limitN*MAPL_PI_R8/180.0d0 limitS=limitS*MAPL_PI_R8/180.0d0 end if call MAPL_StateEval(state,vartomask,temp_field,_RC) if (rank == 2) then !call MAPL_GetPointer(state,rvar2d,vartomask,_RC) call ESMF_FieldGet(temp_field,0,farrayptr=rvar2d,_RC) call MAPL_GetPointer(state,var2d,var_name,_RC) else if (rank == 3) then !call MAPL_GetPointer(state,rvar3d,vartomask,_RC) call ESMF_FieldGet(temp_field,0,farrayptr=rvar3d,_RC) call MAPL_GetPointer(state,var3d,var_name,_RC) else _FAIL('Rank must be 2 or 3') end if if (rank == 2) then var2d = 0.0 where(limitS <= lats .and. lats <=limitN) var2d = rvar2d else if (rank == 3) then var3d = 0.0 do i=1,size(var3d,3) where(limitS <= lats .and. lats <=limitN) var3d(:,:,i) = rvar3d(:,:,i) enddo end if call ESMF_FieldDestroy(temp_field, noGarbage=.true., _RC) _RETURN(_SUCCESS) end subroutine evaluate_zone_mask subroutine evaluate_box_mask(this,state,var_name,rc) class(ExtDataMask), intent(inout) :: this type(ESMF_State), intent(inout) :: state character(len=*), intent(in) :: var_name integer, optional, intent(out) :: rc integer :: status integer :: i character(len=:), allocatable :: vartomask,strtmp real, pointer :: rvar2d(:,:) real, pointer :: rvar3d(:,:,:) real, pointer :: var2d(:,:) real, pointer :: var3d(:,:,:) real(REAL64), pointer :: lats(:,:) real(REAL64), pointer :: lons(:,:) real(REAL64) :: limitS, limitN, limitE, limitW real(REAL64) :: limitE1, limitW1 real(REAL64) :: limitE2, limitW2 type(ESMF_Field) :: field,temp_field type(ESMF_Grid) :: grid integer :: rank,is,nargs integer :: counts(3) logical :: isCube, twoBox real, allocatable :: temp2d(:,:) character(len=ESMF_MAXSTR) :: args(5) type(ESMF_CoordSys_Flag) :: coordSys call ESMF_StateGet(state,var_name,field,_RC) call ESMF_FieldGet(field,rank=rank,grid=grid,_RC) call ESMF_GridGet(grid,coordsys=coordsys,_RC) temp_field = create_field_from_Field(field,_RC) strtmp = this%mask_arguments do nargs=1,5 is = index(strtmp,',') if (is >0) then args(nargs) = strtmp(:is-1) else args(nargs) = strtmp end if strtmp = strtmp(is+1:) end do varToMask=args(1) READ(args(2),*,IOSTAT=status) limitS _VERIFY(status) READ(args(3),*,IOSTAT=status) limitN _VERIFY(status) READ(args(4),*,IOSTAT=status) limitW _VERIFY(status) READ(args(5),*,IOSTAT=status) limitE _VERIFY(status) _ASSERT(limitE > limitW,'LimitE must be greater than limitW') _ASSERT(limitE /= limitW,'LimitE cannot equal limitW') _ASSERT(limitN /= limitS,'LimitN cannot equal LimitS') _ASSERT((limitE-limitW)<=360.0d0,'(LimitE - LimitW) must be less than or equal to 360') call ESMF_GridGetCoord(grid, coordDim=1, localDE=0, & staggerloc=ESMF_STAGGERLOC_CENTER, farrayPtr=lons, rc=status) _VERIFY(status) call ESMF_GridGetCoord(grid, coordDim=2, localDE=0, & staggerloc=ESMF_STAGGERLOC_CENTER, farrayPtr=lats, rc=status) _VERIFY(status) ! do some tests if cube goes from 0 to 360, lat-lon -180 to 180 call MAPL_GridGet(grid, globalCellCountPerDim=COUNTS,rc=status) _VERIFY(STATUS) if (counts(2)==6*counts(1)) then isCube=.true. else isCube=.false. end if twoBox = .false. if (isCube) then if (limitW < 0.0d0 .and. limitE >=0.0d0) then ! need two boxes twoBox=.true. limitW1=0.0d0 limitE1=limitE limitW2=limitW+360.0d0 limitE2=360.0d0 else if (limitW <0.0d0 .and. limitE <0.0d0) then ! just shift limitW1=limitW+360.d0 limitE1=limitE+360.d0 else ! normal case limitW1=limitW limitE1=limitE end if else if (limitW <= 180.0d0 .and. limitE > 180.0d0) then ! need two boxes twoBox=.true. limitW1=limitW limitE1=180.0d0 limitW2=-180.d0 limitE2=limitE-360.0d0 else if (limitW > 180.0d0 .and. limitE > 180.0d0) then ! just shift limitW1=limitW-360.d0 limitE1=limitE-360.d0 else ! normal case limitW1=limitW limitE1=limitE end if end if if (coordSys == ESMF_COORDSYS_SPH_RAD) then limitE1=limitE1*MAPL_PI_R8/180.0d0 limitW1=limitW1*MAPL_PI_R8/180.0d0 if (twoBox) then limitE2=limitE2*MAPL_PI_R8/180.0d0 limitW2=limitW2*MAPL_PI_R8/180.0d0 end if limitN=limitN*MAPL_PI_R8/180.0d0 limitS=limitS*MAPL_PI_R8/180.0d0 end if call MAPL_StateEval(state,varToMask,temp_field,_RC) if (rank == 2) then call ESMF_FieldGet(temp_field,0,farrayptr=rvar2d,_RC) call MAPL_GetPointer(state,var2d,var_name,_RC) else if (rank == 3) then call ESMF_FieldGet(temp_field,0,farrayptr=rvar3d,_RC) call MAPL_GetPointer(state,var3d,var_name,_RC) else _FAIL('Rank must be 2 or 3') end if if (rank == 2) then var2d = 0.0 where(limitS <= lats .and. lats <=limitN .and. limitW1 <= lons .and. lons <= limitE1 ) var2d = rvar2d else if (rank == 3) then var3d = 0.0 do i=1,size(var3d,3) where(limitS <= lats .and. lats <=limitN .and. limitW1 <= lons .and. lons <= limitE1 ) var3d(:,:,i) = rvar3d(:,:,i) enddo end if if (twoBox) then allocate(temp2d(size(var2d,1),size(var2d,2)),stat=status) _VERIFY(STATUS) if (rank == 2) then temp2d = 0.0 where(limitS <= lats .and. lats <=limitN .and. limitW2 <= lons .and. lons <= limitE2 ) temp2d = rvar2d var2d=var2d+temp2d else if (rank == 3) then do i=1,size(var3d,3) temp2d = 0.0 where(limitS <= lats .and. lats <=limitN .and. limitW2 <= lons .and. lons <= limitE2 ) temp2d = rvar3d(:,:,i) var3d(:,:,i)=var3d(:,:,i)+temp2d enddo end if deallocate(temp2d) end if call ESMF_FieldDestroy(temp_field, noGarbage=.true., _RC) _RETURN(_SUCCESS) end subroutine evaluate_box_mask !------------------------------------------------------------------------------ !> ! Extract integers from a character-delimited string, for example, "-1,45,256,7,10". In the context ! of Chem_Util, this is provided for determining the numerically indexed regions over which an ! emission might be applied. ! ! In multiple passes, the string is parsed for the delimiter, and the characters up to, but not ! including the delimiter are taken as consecutive digits of an integer. A negative sign ("-") is ! allowed. After the first pass, each integer and its trailing delimiter are lopped of the head of ! the (local copy of the) string, and the process is started over. ! ! The default delimiter is a comma (","). ! ! "Unfilled" iValues are zero. ! ! Return codes: !1. Zero-length string. !2. iSize needs to be increased. ! !#### Assumptions/bugs: ! ! A non-zero return code does not stop execution. ! Allowed numerals are: 0,1,2,3,4,5,6,7,8,9. ! A delimiter must be separated from another delimiter by at least one numeral. ! The delimiter cannot be a numeral or a negative sign. ! The character following a negative sign must be an allowed numeral. ! The first character must be an allowed numeral or a negative sign. ! The last character must be an allowed numeral. ! The blank character (" ") cannot serve as a delimiter. ! !#### Examples of strings that will work: !``` ! "1" ! "-1" ! "-1,2004,-3" ! "1+-2+3" ! "-1A100A5" !``` ! !#### Examples of strings that will not work: !``` ! "1,--2,3" ! "1,,2,3" ! "1,A,3" ! "1,-,2" ! "1,2,3,4," ! "+1" ! "1 3 6" !``` ! SUBROUTINE ExtDataExtractIntegers(string,iSize,iValues,delimiter,verbose,rc) IMPLICIT NONE CHARACTER(LEN=*), INTENT(IN) :: string !! Character-delimited string of integers INTEGER, INTENT(IN) :: iSize INTEGER, INTENT(INOUT) :: iValues(iSize) !! Space allocated for extracted integers CHARACTER(LEN=*), OPTIONAL :: delimiter !! 1-character delimiter LOGICAL, OPTIONAL, INTENT(IN) :: verbose !! Let me know iValues as they are found. !! DEBUG directive turns on the message even !! if verbose is not present or if !! verbose = .FALSE. INTEGER, OPTIONAL, INTENT(OUT) :: rc !! Return code INTEGER :: base,count,i,iDash,last,lenStr INTEGER :: multiplier,pos,posDelim,sign CHARACTER(LEN=255) :: str CHARACTER(LEN=1) :: char,delimChar LOGICAL :: Done LOGICAL :: tellMe ! Initializations ! --------------- count = 1 Done = .FALSE. iValues(:) = 0 base = ICHAR("0") iDash = ICHAR("-") ! Determine verbosity, letting the DEBUG ! directive override local specification ! -------------------------------------- tellMe = .FALSE. IF(PRESENT(verbose)) THEN IF(verbose) tellMe = .TRUE. END IF #ifdef DEBUG tellMe = .TRUE. #endif ! Check for zero-length string ! ---------------------------- lenStr = LEN_TRIM(string) IF(lenStr == 0) THEN _FAIL("ERROR - Found zero-length string.") END IF ! Default delimiter is a comma ! ---------------------------- delimChar = "," IF(PRESENT(delimiter)) delimChar(1:1) = delimiter(1:1) ! Work on a local copy ! -------------------- str = TRIM(string) ! One pass for each delimited integer ! ----------------------------------- Parse: DO lenStr = LEN_TRIM(str) ! Parse the string for the delimiter ! ---------------------------------- posDelim = INDEX(TRIM(str),TRIM(delimChar)) ! If the delimiter does not exist, ! one integer remains to be extracted. ! ------------------------------------ IF(posDelim == 0) THEN Done = .TRUE. last = lenStr ELSE last = posDelim-1 END IF multiplier = 10**last ! Examine the characters of this integer ! -------------------------------------- Extract: DO pos=1,last char = str(pos:pos) i = ICHAR(char) ! Account for a leading "-" ! ------------------------- IF(pos == 1) THEN IF(i == iDash) THEN sign = -1 ELSE sign = 1 END IF END IF ! "Power" of 10 for this character ! -------------------------------- multiplier = multiplier/10 IF(pos == 1 .AND. sign == -1) CYCLE Extract ! Integer comes from remaining characters ! --------------------------------------- i = (i-base)*multiplier iValues(count) = iValues(count)+i IF(pos == last) THEN iValues(count) = iValues(count)*sign END IF END DO Extract IF(Done) EXIT ! Lop off the leading integer and try again ! ----------------------------------------- str(1:lenStr-posDelim) = str(posDelim+1:lenStr) str(lenStr-posDelim+1:255) = " " count = count+1 ! Check size ! ---------- IF(count > iSize) THEN _FAIL("ERROR - iValues does not have enough elements.") END IF END DO Parse _RETURN(ESMF_SUCCESS) END SUBROUTINE ExtDataExtractIntegers function create_field_from_field(input_field,rc) result(output_field) type(ESMF_Field) :: output_field type(ESMF_Field), intent(in) :: input_field integer, optional, intent(out) :: rc integer :: status type(ESMF_Grid) :: grid integer :: rank type(ESMF_TypeKind_Flag) :: typekind integer :: lb(1),ub(1) call ESMF_FieldGet(input_field,grid=grid,rank=rank,typekind=typekind,ungriddedLBound=lb,ungriddedUBound=ub,_RC) if (rank==2) then output_field = ESMF_FieldCreate(grid,typekind,_RC) else if (rank==3) then output_field = ESMF_FieldCreate(grid,typekind,ungriddedLBound=lb,ungriddedUBound=ub,name="temp_field",_RC) end if _RETURN(_SUCCESS) end function end module MAPL_ExtDataMask