! ! NAME ! grid_interpol - from sh/gg/ll to gg/ll ! ! HISTORY ! v1.? ! Original ! v1.1 ! Bug fixed: southpole in gg -> ll ! module Grid_Interpol implicit none ! --- in/out ------------------------------ private public :: Interpol public :: NewInterpol public :: InterpolMask public :: Aver public :: IntArea public :: IntLat, IntLon !public :: ShTruncation, ShRefinement public :: Tgg2llFracs, Tll2ggFracs, Init, Done, FracSum ! --- const ---------------------------------- character(len=*), parameter :: mname = 'grid_interpol' ! --- types ------------------------------- type Tgg2llFracs integer :: nlon, nlat, np integer, pointer :: ncov(:,:) integer, pointer :: indx(:,:,:) real, pointer :: frac(:,:,:) real, pointer :: A_gg(:) ! m2 real, pointer :: A_ll(:,:) ! m2 end type Tgg2llFracs type Tll2ggFracs integer :: nlon, nlat, np integer, pointer :: ncov(:) ! (1:np) integer, pointer :: ii(:,:) ! (1:np,1:max(ncov(k))) integer, pointer :: jj(:,:) ! (1:np,1:max(ncov(k))) real, pointer :: ff(:,:) ! (1:np,1:max(ncov(k))) real, pointer :: A_gg(:) ! m2 real, pointer :: A_ll(:,:) ! m2 end type Tll2ggFracs ! --- interfaces -------------------------- interface Init module procedure gg2ll_Init module procedure ll2gg_Init end interface interface Done module procedure gg2ll_Done module procedure ll2gg_Done end interface interface FracSum module procedure gg2ll_FracSum module procedure ll2gg_FracSum end interface interface Interpol module procedure Interpol_sh_gg module procedure Interpol_shi_gg module procedure Interpol_sh_ll module procedure Interpol_shi_ll module procedure Interpol_gg_ll module procedure Interpol_ll_gg end interface interface NewInterpol module procedure NewInterpol_shi_gg end interface interface Aver module procedure Aver_gg_ll module procedure Aver_sh_ll end interface interface IntArea module procedure IntArea_sh_ll_f module procedure IntArea_shi_ll_f module procedure IntArea_sh_ll_fgh module procedure IntArea_shi_ll_fgh module procedure IntArea_sh_ll_fh module procedure IntArea_shi_ll_fh end interface interface IntLat module procedure IntLat_sh_ll module procedure IntLat_shi_ll end interface interface IntLon module procedure IntLon_sh_ll module procedure IntLon_shi_ll end interface contains ! ========================================================= ! === ! === evaluate spectral fields ! === ! ========================================================= ! from spectral to reduced gausian grid subroutine Interpol_sh_gg( sh, ggi, gg, status ) use Grid_Type_sh, only : TshGrid, Eval_Lons, Check use Grid_Type_gg, only : TggGridInfo, Check ! --- in/out ---------------------------------- type(TshGrid), intent(in) :: sh type(TggGridInfo), intent(in) :: ggi real, intent(out) :: gg(ggi%np) integer, intent(out) :: status ! --- const -------------------------------- character(len=*), parameter :: rname = mname//'/Interpol_sh_gg' ! --- local ----------------------------------- !real, allocatable :: llgrid(:,:) real, allocatable :: llgrid(:) integer :: nlon integer :: jn !, js ! --- begin ----------------------------------- call Check( sh ) call Check( ggi, gg ) !allocate( llgrid(maxval(ggi%nlon),2) ) allocate( llgrid(maxval(ggi%nlon)) ) ! northern rows: !do jn = 1, ggi%nlat/2 do jn = 1, ggi%nlat ! southern row: !js = ggi%nlat + 1 - jn ! only if one of the rows is marked: !if ( ggi%latflag(jn) .or. ggi%latflag(js) ) then if ( ggi%latflag(jn) ) then nlon = ggi%nlon(jn) !call Eval_Lons( llgrid(1:nlon,1:2), sh, ggi%lat(jn), nlon, 0.0, nlon ) !gg(ggi%i1(jn):ggi%im(jn)) = llgrid(1:nlon,1) !gg(ggi%i1(js):ggi%im(js)) = llgrid(1:nlon,2) call Eval_Lons( llgrid(1:nlon), sh, ggi%lat(jn), nlon, 0.0, nlon, status ) gg(ggi%i1(jn):ggi%im(jn)) = llgrid(1:nlon) end if end do deallocate( llgrid ) end subroutine Interpol_sh_gg ! * subroutine NewInterpol_shi_gg( shi, shc, ggi, gg, status ) use Grid_Type_sh, only : TshGridInfo, Check, Eval_Lons use Grid_Type_sh, only : TshGrid, Init, Done, Set use Grid_Type_gg, only : TggGridInfo, Check ! --- in/out ---------------------------------- type(TshGridInfo), intent(in) :: shi complex, intent(in) :: shc(shi%np) type(TggGridInfo), intent(in) :: ggi real, intent(out) :: gg(ggi%np) integer, intent(out) :: status ! --- const -------------------------------- character(len=*), parameter :: rname = mname//'/NewInterpol_shi_gg' ! --- local ----------------------------------- real, pointer :: llgrid(:) integer :: nlon integer :: jn !, js ! --- begin ----------------------------------- allocate( llgrid(maxval(ggi%nlon)) ) ! northern rows: !do jn = 1, ggi%nlat/2 do jn = 1, ggi%nlat ! southern row: !js = ggi%nlat + 1 - jn ! only if one of the rows is marked: !if ( ggi%latflag(jn) .or. ggi%latflag(js) ) then if ( ggi%latflag(jn) ) then nlon = ggi%nlon(jn) !call Eval_Lons( llgrid(1:nlon,1:2), sh, ggi%lat(jn), nlon, 0.0, nlon ) !gg(ggi%i1(jn):ggi%im(jn)) = llgrid(1:nlon,1) !gg(ggi%i1(js):ggi%im(js)) = llgrid(1:nlon,2) call Eval_Lons( shi, shc, ggi%lat(jn), nlon, 0.0, nlon, llgrid(1:nlon), status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if gg(ggi%i1(jn):ggi%im(jn)) = llgrid(1:nlon) end if end do ! done deallocate( llgrid ) ! ok status = 0 end subroutine NewInterpol_shi_gg ! * subroutine Interpol_shi_gg( shi, shc, ggi, gg, status ) use Grid_Type_sh, only : TshGridInfo, Check, Eval_Lons use Grid_Type_sh, only : TshGrid, Init, Done, Set use Grid_Type_gg, only : TggGridInfo, Check ! --- in/out ---------------------------------- type(TshGridInfo), intent(in) :: shi complex, intent(in) :: shc(shi%np) type(TggGridInfo), intent(in) :: ggi real, intent(out) :: gg(ggi%np) integer, intent(out) :: status ! --- const -------------------------------- character(len=*), parameter :: rname = mname//'/Interpol_shi_gg' ! --- local ----------------------------------- type(TshGrid) :: sh !real, allocatable :: llgrid(:,:) real, allocatable :: llgrid(:) integer :: nlon integer :: jn !, js ! --- begin ----------------------------------- ! store input in old type grid: call Init( sh ) call Set( sh, shi%T, shc ) !allocate( llgrid(maxval(ggi%nlon),2) ) allocate( llgrid(maxval(ggi%nlon)) ) ! northern rows: !do jn = 1, ggi%nlat/2 do jn = 1, ggi%nlat ! southern row: !js = ggi%nlat + 1 - jn ! only if one of the rows is marked: !if ( ggi%latflag(jn) .or. ggi%latflag(js) ) then if ( ggi%latflag(jn) ) then nlon = ggi%nlon(jn) !call Eval_Lons( llgrid(1:nlon,1:2), sh, ggi%lat(jn), nlon, 0.0, nlon ) !gg(ggi%i1(jn):ggi%im(jn)) = llgrid(1:nlon,1) !gg(ggi%i1(js):ggi%im(js)) = llgrid(1:nlon,2) call Eval_Lons( llgrid(1:nlon), sh, ggi%lat(jn), nlon, 0.0, nlon, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if gg(ggi%i1(jn):ggi%im(jn)) = llgrid(1:nlon) end if end do ! done deallocate( llgrid ) call Done( sh ) ! ok status = 0 end subroutine Interpol_shi_gg ! === ! from sh to ll subroutine Interpol_sh_ll( sh, lli, ll, status ) use Grid_Type_sh, only : TshGrid, Eval_Lons, Check use Grid_Type_ll, only : TllGridInfo, Check ! --- in/out ------------------------------- type(TshGrid), intent(in) :: sh type(TllGridInfo), intent(in) :: lli real, intent(out) :: ll(lli%im,lli%jm) integer, intent(out) :: status ! --- const -------------------------------- character(len=*), parameter :: rname = mname//'/Interpol_sh_ll' ! --- local -------------------------------- integer :: j ! --- begin -------------------------------- call Check( sh ) call Check( lli, 'n', ll, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! loop over latitudes do j = 1, lli%jm ! evaluate at oposite latidutes: call Eval_Lons( ll(:,j), sh, lli%lat(j), int(360.0/lli%dlon_deg), & lli%lon(1), lli%im, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do ! ok status = 0 end subroutine Interpol_sh_ll ! * ! from shi/sh to ll subroutine Interpol_shi_ll( shi, shc, lli, ll, status ) use Grid_Type_sh, only : TshGridInfo, Check, Eval_Lons use Grid_Type_sh, only : TshGrid, Init, Done, Set use Grid_Type_ll, only : TllGridInfo, Check ! --- in/out ------------------------------- type(TshGridInfo), intent(in) :: shi complex, intent(in) :: shc(shi%np) type(TllGridInfo), intent(in) :: lli real, intent(out) :: ll(lli%im,lli%jm) integer, intent(out) :: status ! --- const -------------------------------- character(len=*), parameter :: rname = mname//'/Interpol_shi_ll' ! --- local -------------------------------- type(TshGrid) :: sh integer :: j ! --- begin -------------------------------- ! store input in old type grid: call Init( sh ) call Set( sh, shi%T, shc ) ! check lat/lon arguments: call Check( lli, 'n', ll, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! loop over latitudes do j = 1, lli%jm ! evaluate spectral field: call Eval_Lons( ll(:,j), sh, lli%lat(j), int(360.0/lli%dlon_deg), & lli%lon(1), lli%im, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do ! done call Done( sh ) ! ok status = 0 end subroutine Interpol_shi_ll ! ========================================================= ! === ! === interpolate from gg ! === ! ========================================================= ! Return a logical array with size gg; ! true if a point is required for interpolation from gg to ll ! ! Include 'depth' extra rows. ! ! Also set flags in ggi for each row. ! subroutine InterpolMask( ggi, gg, lli, depth ) use Grid_Type_gg, only : TggGridInfo, Check use Grid_Type_ll, only : TllGridInfo, Check ! --- in/out ------------------------------- type(TggGridInfo), intent(inout) :: ggi logical, intent(out) :: gg(ggi%np) type(TllGridInfo), intent(in) :: lli integer, intent(in) :: depth ! --- local -------------------------------- integer :: j, jn, js ! --- begin ------------------------------ ! note: gg lats from north to south! !print *, 'llilat:',lli%blat_deg(0), lli%blat_deg(lli%nlat) ! search south most gg lat outside lli js = 1 do if ( js == ggi%nlat ) exit if ( ggi%lat(js) < lli%blat(0) ) exit js = js + 1 end do js = min( ggi%nlat, js+depth ) ! search north-most gg lat row outside lli jn = ggi%nlat do if ( jn == 1 ) exit if ( ggi%lat(jn) > lli%blat(lli%nlat) ) exit jn = jn - 1 end do jn = max( 1, jn-depth ) ! set all rows between js and jn to true: gg = .false. ggi%latflag = .false. do j = jn, js ggi%latflag(j) = .true. gg(ggi%i1(j):ggi%im(j)) = .true. end do end subroutine InterpolMask ! ************************************************************* ! ! Determine fraction of gg cell that covers a ll cell. ! ! Returns three arrays: ! ! integer :: ncov(im,jm) ! integer :: indx(im,jm,max_nconv) ! real :: frac(im,jm,max_nconv) ! ! For an ll cell (i,j), 'ncov' gives the number of gg cells covering it. ! Their gg indices are stored in : indx(i,j,1:ncov(i,j)), ! corresponding fractions are in : frac(i,j,1:ncov(i,j)) . ! subroutine gg2ll_Init( gg2ll, ggi, lli, status ) use num, only : Interval use grid_tools, only : pi, ll_area_frac use grid_type_gg, only : TggGridInfo, AreaOper use grid_type_ll, only : TllGridInfo, AreaOper ! --- in/out --------------------------- type(Tgg2llFracs), intent(out) :: gg2ll type(TggGridInfo), intent(in) :: ggi type(TllGridInfo), intent(in) :: lli integer, intent(out) :: status ! --- const --------------------------------- character(len=*), parameter :: rname = mname//'/gg2ll_Init' ! --- local ---------------------------- integer :: max_ncov_lon, max_ncov_lat, max_ncov integer :: i, j integer :: iflag integer :: gi, gi1, gi2 integer :: gj, gj1, gj2 integer :: nlon real :: gblon(0:ggi%nlon_reg) real :: west1, east1, south1, north1 real :: west2, east2, south2, north2 real :: frac integer :: ncov ! --- begin ---------------------------- ! save dimension gg2ll%nlon = lli%nlon gg2ll%nlat = lli%nlat gg2ll%np = ggi%np ! estimate maximum number of gg cells covering an ll cell: max_ncov_lon = ceiling(lli%dlon/minval(ggi%dlon)) max_ncov_lat = ceiling(lli%dlat/minval(ggi%dlat)) max_ncov = (max_ncov_lon+1) * (max_ncov_lat+1) ! allocate arrays allocate( gg2ll%ncov(lli%nlon,lli%nlat) ) allocate( gg2ll%indx(lli%nlon,lli%nlat,max_ncov) ) allocate( gg2ll%frac(lli%nlon,lli%nlat,max_ncov) ) ! zero by default: gg2ll%ncov = 0 gg2ll%indx = 0 gg2ll%frac = 0.0 gj1 = 0 gj2 = 0 ! loop over ll cells from north to south (gg direction): do j = lli%nlat, 1, -1 do i = 1, lli%nlon ! extract boundaries of ll cell: west2 = lli%blon(i-1) east2 = lli%blon(i) if ( east2 < 0.0 ) then west2 = west2 + 2*pi ! (0,2pi) east2 = east2 + 2*pi ! (0,2pi) end if south2 = lli%blat(j-1) north2 = lli%blat(j) ! search gg rows including north and south ll lat; ! negative lats to get increasing values ... call Interval( -ggi%blat, -north2 , gj1, iflag ) if ( iflag /= 0 ) stop 'BUG IN gg2ll_Init : wrong iflag gj1' call Interval( -ggi%blat, -south2, gj2, iflag ) if ( iflag /= 0 ) stop 'BUG IN gg2ll_Init : wrong iflag gj2' gi1 = 0 gi2 = 0 ! loop over gg lat rows: do gj = gj1, gj2 ! boundary lons nlon = ggi%nlon(gj) do gi = 0, nlon gblon(gi) = (gi-0.5)*ggi%dlon(gj) end do ! search cells including west and east bound of ll cell if ( west2 < gblon(0) ) then call Interval( gblon(0:nlon), west2+2*pi, gi1, iflag ) else if ( west2 > gblon(nlon) ) then call Interval( gblon(0:nlon), west2-2*pi, gi1, iflag ) else call Interval( gblon(0:nlon), west2 , gi1, iflag ) end if if ( iflag /= 0 ) then print *, 'gblon=', gblon(0:nlon) print *, 'west2=',west2 print *, 'iflag=',iflag stop 'BUG IN gg2ll_Init : wrong iflag gi1' end if if ( east2 < gblon(0) ) then call Interval( gblon(0:nlon), east2+2*pi, gi2, iflag ) else if ( east2 > gblon(nlon) ) then call Interval( gblon(0:nlon), east2-2*pi, gi2, iflag ) else call Interval( gblon(0:nlon), east2 , gi2, iflag ) end if if ( iflag /= 0 ) then print *, 'gblon=', gblon(0:nlon) print *, 'east2=',east2 stop 'BUG IN gg2ll_Init : wrong iflag gi2' end if ! loop over all gg cells in current row: gi = gi1 do ! extract boundaries of gg cell: west1 = (gi-1.5)*ggi%dlon(gj) ! (0,2pi) east1 = (gi-0.5)*ggi%dlon(gj) ! (0,2pi) south1 = ggi%blat(gj) north1 = ggi%blat(gj-1) ! shift if gg cell is right from [0,2pi] if ( west1 > east2 ) then west1 = west1 - 2*pi ! (0,2pi) east1 = east1 - 2*pi ! (0,2pi) end if ! determine covarage fraction: frac = ll_area_frac( west1, east1, south1, north1, & west2, east2, south2, north2 ) ! fill fraction: if ( frac > 0.0 .and. frac <= 1.0 ) then ncov = gg2ll%ncov(i,j) + 1 gg2ll%ncov(i,j) = ncov gg2ll%indx(i,j,ncov) = ggi%i1(gj)-1 + gi gg2ll%frac(i,j,ncov) = frac else if ( abs(frac) < 1.0e-4 ) then ! almost no coverage ... else if ( abs(frac-1.0) < 1.0e-4 ) then !print *, 'WARNING in module grid_interpol, gg2ll_Init' !print *, 'frac=',frac !print *, ' 1:', west1, east1, south1, north1 !print *, ' 2:', west2, east2, south2, north2 ncov = gg2ll%ncov(i,j) + 1 gg2ll%ncov(i,j) = ncov gg2ll%indx(i,j,ncov) = ggi%i1(gj)-1 + gi gg2ll%frac(i,j,ncov) = 1.0 else print *, 'ERROR in module grid_interpol, gg2ll_Init' print *, 'frac=',frac print *, ' 1:', west1, east1, south1, north1 print *, ' 2:', west2, east2, south2, north2 stop end if if ( gi == gi2 ) exit gi = gi + 1 if ( gi == nlon+1 ) gi = 1 end do end do end do end do ! store cell area's allocate( gg2ll%A_gg(ggi%np) ) call AreaOper( ggi, gg2ll%A_gg, '=', 'm2' ) allocate( gg2ll%A_ll(lli%nlon,lli%nlat) ) call AreaOper( lli, gg2ll%A_ll, '=', 'm2', status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! ok status = 0 end subroutine gg2ll_Init ! *** subroutine gg2ll_Done( gg2ll ) ! --- in/out --------------------------- type(Tgg2llFracs), intent(inout) :: gg2ll ! --- begin ---------------------------- ! deallocate arrays deallocate( gg2ll%ncov ) deallocate( gg2ll%indx ) deallocate( gg2ll%frac ) deallocate( gg2ll%A_gg ) deallocate( gg2ll%A_ll ) end subroutine gg2ll_Done ! *** ! ! ncov(i,j) ! ll(i,j) = sum gg(indx(i,j,k)) * frac(i,j,k) ! k=1 ! subroutine gg2ll_FracSum( gg2ll, gg, ll, status, key ) ! --- in/out --------------------------- type(Tgg2llFracs), intent(in) :: gg2ll real, intent(in) :: gg(:) real, intent(out) :: ll(:,:) integer, intent(out) :: status character(len=*), intent(in), optional :: key ! --- const ---------------------------------- character(len=*), parameter :: rname = mname//'/gg2ll_FracSum' ! --- local ---------------------------- character(len=10) :: the_key integer :: i, j, k integer :: ncov, indx real :: fac ! --- begin ---------------------------- the_key = 'none' if ( present(key) ) the_key = key if ( any(shape(gg2ll%ncov)/=shape(ll)) ) then write (*,'("ERROR - shapes of ll and gg2ll do not match:")') write (*,'("ERROR - ll : ",2i4)') shape(ll) write (*,'("ERROR - gg2ll : ",2i4)') shape(gg2ll%ncov) write (*,'("ERROR in ",a)') rname; status=1; return end if if ( minval(gg2ll%indx)<0 .or. maxval(gg2ll%indx)>size(gg) ) then write (*,'("ERROR - indices of gg array in gg2ll out of range:")') write (*,'("ERROR - gg : 1 - ",i4)') size(gg) write (*,'("ERROR - gg2ll : ",i4," - ",i4)') minval(gg2ll%indx), maxval(gg2ll%indx) write (*,'("ERROR in ",a)') rname; status=1; return end if if ( minval(gg2ll%frac)<0.0 .or. maxval(gg2ll%frac)>1.0 ) then write (*,'("ERROR - fraction in gg2ll out of range:")') write (*,'("ERROR - gg2ll: ",i4," - ",i4)') minval(gg2ll%frac), maxval(gg2ll%frac) write (*,'("ERROR in ",a)') rname; status=1; return end if ll = 0.0 do j = 1, gg2ll%nlat do i = 1, gg2ll%nlon ncov = gg2ll%ncov(i,j) if ( ncov > 0 ) then do k = 1, ncov indx = gg2ll%indx(i,j,k) select case ( the_key ) case ( 'none' ) fac = 1.0 case ( 'area-aver', 'area-sum' ) fac = gg2ll%A_gg(indx) case default write (*,'("ERROR - key `",a,"` not supported")') trim(the_key) write (*,'("ERROR in ",a)') rname; status=1; return end select ll(i,j) = ll(i,j) + gg(indx) * fac * gg2ll%frac(i,j,k) end do end if select case ( the_key ) case ( 'none', 'area-sum' ) ! nothing to be done case ( 'area-aver' ) ll(i,j) = ll(i,j) / gg2ll%A_ll(i,j) case default write (*,'("ERROR - key `",a,"` not supported")') trim(the_key) write (*,'("ERROR in ",a)') rname; status=1; return end select end do end do ! ok status = 0 end subroutine gg2ll_FracSum ! ************************************************************* ! ! Determine fraction of gg cell that covers a ll cell. ! ! Returns three arrays: ! ! integer :: ncov(im,jm) ! integer :: indx(im,jm,max_nconv) ! real :: frac(im,jm,max_nconv) ! ! For a gg cell 'k', 'ncov(k)' gives the number of ll cells covering it. ! Their ll indices are stored in : ii(k,1:ncov(k)), jj(k,1:ncov(k)) ; ! corresponding fractions are in : frac(k,1:ncov(k)) . ! subroutine ll2gg_Init( ll2gg, lli, ggi, status ) use num, only : Interval use grid_tools, only : pi, ll_area_frac use grid_type_gg, only : TggGridInfo, AreaOper use grid_type_ll, only : TllGridInfo, AreaOper ! --- in/out --------------------------- type(Tll2ggFracs), intent(out) :: ll2gg type(TllGridInfo), intent(in) :: lli type(TggGridInfo), intent(in) :: ggi integer, intent(out) :: status ! --- const ------------------------------------ character(len=*), parameter :: rname = mname//'/ll2gg_Init' ! --- local ---------------------------- integer :: max_ncov_lon, max_ncov_lat, max_ncov integer :: gi, gi1, gi2 integer :: gj, gj1, gj2 integer :: li, li1, li2 integer :: lj, lj1, lj2 real :: west2, east2, south2, north2 real :: west1, east1, south1, north1 integer :: iflag integer :: gp real :: frac integer :: ncov ! --- begin ---------------------------- ! save dimension ll2gg%nlon = lli%nlon ll2gg%nlat = lli%nlat ll2gg%np = ggi%np ! estimate maximum number of ll cells covering a gg cell: max_ncov_lon = ceiling(maxval(ggi%dlon)/lli%dlon) max_ncov_lat = ceiling(maxval(ggi%dlat)/lli%dlat) max_ncov = (max_ncov_lon+1) * (max_ncov_lat+1) ! allocate arrays allocate( ll2gg%ncov(ggi%np) ) allocate( ll2gg%ii(ggi%np,max_ncov) ) allocate( ll2gg%jj(ggi%np,max_ncov) ) allocate( ll2gg%ff(ggi%np,max_ncov) ) ! zero by default: ll2gg%ncov = 0 ll2gg%ii = 0 ll2gg%jj = 0 ll2gg%ff = 0.0 ! index in gg arrays: gp = 0 ! row range in ll grid: lj1 = 0 lj2 = 0 ! loop over gg rows: do gj = 1, ggi%nlat ! extract north/south boundaries of gg cell: north2 = ggi%blat(gj-1) ! [-pi,pi] south2 = ggi%blat(gj) ! [-pi,pi] ! search ll rows including north and south gg lat; call Interval( lli%blat, south2, lj1, iflag ) if ( iflag /= 0 ) stop 'BUG IN ll2gg_Init : wrong iflag lj1' call Interval( lli%blat, north2, lj2, iflag ) if ( iflag /= 0 ) stop 'BUG IN ll2gg_Init : wrong iflag lj2' ! loop over cells in row: do gi = 1, ggi%nlon(gj) ! next index in 1d row ... gp = gp + 1 ! set east/west boundaries of gg cell: west2 = (gi-1.5) * ggi%dlon(gj) ! [0,2pi] east2 = (gi-0.5) * ggi%dlon(gj) ! [0,2pi] ! loop over ll lat rows: do lj = lj1, lj2 ! search ll cells including west and east bound of gg cell if ( west2 > lli%blon(lli%nlon) ) then call Interval( lli%blon, west2-2*pi, li1, iflag ) else call Interval( lli%blon, west2 , li1, iflag ) end if if ( iflag /= 0 ) then print *, 'lli%blon=', lli%blon print *, 'west2=',west2 print *, 'iflag=',iflag stop 'BUG IN ll2gg_Init : wrong iflag li1' end if if ( east2 > lli%blon(lli%nlon) ) then call Interval( lli%blon, east2-2*pi, li2, iflag ) else call Interval( lli%blon, east2 , li2, iflag ) end if if ( iflag /= 0 ) then print *, 'lli%blon=', lli%blon print *, 'east2=',east2 print *, 'east2-2pi=',east2-2*pi stop 'BUG IN ll2gg_Init : wrong iflag li2' end if ! loop over ll lon cels: li = li1 do ! extract boundaries of ll cell: west1 = lli%blon(li-1) ! [-pi,pi] east1 = lli%blon(li ) ! [-pi,pi] south1 = lli%blat(lj-1) ! [-pi,pi] north1 = lli%blat(lj ) ! [-pi,pi] ! shift if completely left from gg cell: if ( east1 < west2 ) then west1 = west1 + 2*pi east1 = east1 + 2*pi end if ! determine covarage fraction: frac = ll_area_frac( west1, east1, south1, north1, & west2, east2, south2, north2 ) ! fill fraction: if ( frac > 0.0 .and. frac <= 1.0 ) then ncov = ll2gg%ncov(gp) + 1 ll2gg%ncov(gp) = ncov ll2gg%ii(gp,ncov) = li ll2gg%jj(gp,ncov) = lj ll2gg%ff(gp,ncov) = frac else if ( abs(frac) < 1.0e-4 ) then ! almost no coverage ... else if ( abs(frac-1.0) < 1.0e-4 ) then !print *, 'WARNING in module grid_interpol, ll2gg_Init' !print *, 'frac=',frac !print *, ' 1:', west1, east1, south1, north1 !print *, ' 2:', west2, east2, south2, north2 ncov = ll2gg%ncov(gp) + 1 ll2gg%ncov(gp) = ncov ll2gg%ii(gp,ncov) = li ll2gg%jj(gp,ncov) = lj ll2gg%ff(gp,ncov) = 1.0 else print *, 'ERROR in module grid_interpol, ll2gg_Init' print *, 'frac=',frac print *, ' 1:', west1, east1, south1, north1 print *, ' 2:', west2, east2, south2, north2 stop end if if ( li == li2 ) exit li = li + 1 if ( li == lli%nlon+1 ) li = 1 end do ! ll i end do ! ll j end do ! gg i end do ! gg j ! store cell area's allocate( ll2gg%A_gg(ggi%np) ) call AreaOper( ggi, ll2gg%A_gg, '=', 'm2' ) ! allocate( ll2gg%A_ll(lli%nlon,lli%nlat) ) call AreaOper( lli, ll2gg%A_ll, '=', 'm2', status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! ok status = 0 end subroutine ll2gg_Init ! *** subroutine ll2gg_Done( ll2gg ) ! --- in/out --------------------------- type(Tll2ggFracs), intent(inout) :: ll2gg ! --- begin ---------------------------- ! deallocate arrays deallocate( ll2gg%ncov ) deallocate( ll2gg%ii ) deallocate( ll2gg%jj ) deallocate( ll2gg%ff ) deallocate( ll2gg%A_gg ) deallocate( ll2gg%A_ll ) end subroutine ll2gg_Done ! *** ! ! ncov(p) ! gg(p) = sum ll(ii(p,k),jj(p,k)) * ff(p,k) ! k=1 ! subroutine ll2gg_FracSum( ll2gg, ll, gg, key ) ! --- in/out --------------------------- type(Tll2ggFracs), intent(in) :: ll2gg real, intent(in) :: ll(:,:) real, intent(out) :: gg(:) character(len=*), intent(in), optional :: key ! --- local ---------------------------- character(len=10) :: the_key integer :: gp, i, j, k integer :: ncov real :: fac ! --- begin ---------------------------- the_key = 'none' if ( present(key) ) the_key = key if ( size(ll2gg%ncov) /= size(gg) ) then print *, 'shapes of gg and ll2gg do not match:' print *, ' ll : ', size(gg) print *, ' ll2gg: ', size(ll2gg%ncov) stop 'FATAL ERROR IN ll2gg_FracSum' end if if ( minval(ll2gg%ii)<0 .or. maxval(ll2gg%ii)>size(ll,1) ) then print *, 'indices of ll array in ll2gg out of range:' print *, ' ll i : 1 - ', size(ll,1) print *, ' ll2gg: ', minval(ll2gg%ii),'-',maxval(ll2gg%ii) stop 'FATAL ERROR IN ll2gg_FracSum' end if if ( minval(ll2gg%jj)<0 .or. maxval(ll2gg%jj)>size(ll,2) ) then print *, 'indices of ll array in ll2gg out of range:' print *, ' ll j : 1 - ', size(ll,2) print *, ' ll2gg: ', minval(ll2gg%jj),'-',maxval(ll2gg%jj) stop 'FATAL ERROR IN ll2gg_FracSum' end if if ( minval(ll2gg%ff)<0.0 .or. maxval(ll2gg%ff)>1.0 ) then print *, 'fraction in ll2gg out of range:' print *, ' ll2gg: ', minval(ll2gg%ff),'-',maxval(ll2gg%ff) stop 'FATAL ERROR IN ll2gg_FracSum' end if ! init to zero: gg = 0.0 ! loop over gg do gp = 1, ll2gg%np ncov = ll2gg%ncov(gp) if ( ncov > 0 ) then do k = 1, ncov i = ll2gg%ii(gp,k) j = ll2gg%jj(gp,k) select case ( the_key ) case ( 'none' ) fac = 1.0 case ( 'area-aver', 'area-sum' ) fac = ll2gg%A_ll(i,j) case default print *, 'Sorry, key "'//trim(the_key)//'" not supported.' stop 'BUG IN ll2gg_FracSum' end select !if (abs(ll(i,j))>0.0) then !print *, gp,':',gg(gp) ,'+', ll(i,j) ,'*', fac ,'*', ll2gg%ff(gp,k),'=',gg(gp) + ll(i,j) * fac * ll2gg%ff(gp,k) !endif gg(gp) = gg(gp) + ll(i,j) * fac * ll2gg%ff(gp,k) end do end if select case ( the_key ) case ( 'none', 'area-sum' ) ! nothing to be done !if (abs(gg(gp))>0.0) then !print *, ' ' !endif case ( 'area-aver' ) gg(gp) = gg(gp) / ll2gg%A_gg(gp) !if (abs(gg(gp))>0.0) then !print *, gp,':',gg(gp),'/',ll2gg%A_gg(gp),'=',gg(gp) / ll2gg%A_gg(gp) !print *, ' ' !endif case default print *, 'Sorry, key "'//trim(the_key)//'" not supported.' stop 'BUG IN ll2gg_FracSum' end select end do ! gg cells end subroutine ll2gg_FracSum ! ************************************************************* ! gg to ll subroutine Interpol_gg_ll( ggi, gg, lli, ll, status ) use Binas, only : deg2rad use Num, only : Interp_Lin, CircInterp_Lin use Grid_Type_gg, only : TggGridInfo, Check use Grid_Type_ll, only : TllGridInfo, Check ! --- in/out ------------------------------- type(TggGridInfo), intent(in) :: ggi real, intent(in) :: gg(ggi%np) type(TllGridInfo), intent(in) :: lli real, intent(out) :: ll(lli%im,lli%jm) integer, intent(out) :: status ! --- const ---------------------------------- character(len=*), parameter :: rname = mname//'/Interpol_gg_ll' ! --- local -------------------------------- integer :: nlon, nlon_max integer :: nlat integer :: i1, im real, allocatable :: lons(:) real, allocatable :: row(:) real, allocatable :: gl(:,:) real, allocatable :: gl_lat(:) integer :: i, j integer :: j1, jm integer :: ilast ! --- begin -------------------------------- call Check( ggi, gg ) call Check( lli, 'n', ll, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if nlat = ggi%nlat nlon_max = maxval(ggi%nlon) ! ll in lon, gg in lat allocate( gl(size(ll,1),0:nlat+1) ); gl = 0.0 ! latitudes from south->north allocate( gl_lat(0:nlat+1) ) gl_lat(0) = -90.0 * deg2rad ! south pole (rad) do j = 1, nlat gl_lat(nlat+1-j) = ggi%lat(j) ! rad end do gl_lat(nlat+1) = 90.0 * deg2rad ! north pole (rad) ! row in gg grid; doubled from -360.0 to 360.0 allocate( lons(nlon_max) ); lons = 0.0 allocate( row(nlon_max) ) ; row = 0.0 ! select first and last Gaussian lat: j1 = nlat do if ( (j1 == 1) .or. (ggi%lat(j1) > maxval(lli%lat)) ) exit j1 = j1 - 1 end do jm = 1 do if ( (jm == nlat) .or. (ggi%lat(jm) < minval(lli%lat)) ) exit jm = jm + 1 end do ! loop over Gaussian latitudes, from north to south! do j = j1, jm ! number of lon points at this latitude: nlon = ggi%nlon(j) ! start and end indices in grid point array i1 = ggi%i1(j) im = ggi%im(j) ! lons in [0,2pi) do i = 1, nlon lons(i) = (i-1)*360.0/nlon end do lons(1:nlon) = lons(1:nlon) * deg2rad ! grid values (doubled) row(1:nlon) = gg(i1:im) ! set north pole (j=1 in gg, j=nlat+1 in gl) if ( j == 1 ) then gl(:,nlat+1) = sum(row(1:nlon))/nlon end if ! set south pole if ( j == nlat ) then gl(:,0) = sum(row(1:nlon))/nlon end if ! Interpolate over lon (circular arrays); ! swap lats to ensure south -> north: ilast = 1 do i = 1, lli%im call CircInterp_Lin( lons(1:nlon), 360.0*deg2rad, row(1:nlon), & lli%lon(i), gl(i,nlat+1-j), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do end do ! Linear interpolation over lat: ilast = 1 do j = 1, lli%jm do i = 1, lli%im call Interp_Lin( gl_lat, gl(i,:), lli%lat(j), ll(i,j), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do end do ! free memory deallocate( gl ) deallocate( gl_lat ) deallocate( lons ) deallocate( row ) ! ok status = 0 end subroutine Interpol_gg_ll ! ========================================================= ! === ! === interpolate from ll ! === ! ========================================================= ! ll to gg ! ! First interpol in lat direction to Gaussian lats, ! then interpolate in lon direction. ! NOTE: MIPSpro compiler gives errors if argument names are the ! same as used for Interpol_gg_ll ... subroutine Interpol_ll_gg( lli, ll, xggi, xgg, status ) use Binas, only : deg2rad use Num, only : Interp_Lin, CircInterp_Lin use Grid_Type_gg, only : TggGridInfo, Check use Grid_Type_ll, only : TllGridInfo, Check ! --- in/out ------------------------------- type(TllGridInfo), intent(in) :: lli real, intent(in) :: ll(lli%im,lli%jm) type(TggGridInfo), intent(in) :: xggi real, intent(out) :: xgg(1:xggi%np) integer, intent(out) :: status ! --- const ---------------------------------- character(len=*), parameter :: rname = mname//'/Interpol_ll_gg' ! --- local -------------------------------- integer :: nlon integer :: nlat integer :: i, j integer :: i0 integer :: ilast real :: dlon_deg real :: glat, glon real :: period real, allocatable :: row(:) ! --- begin -------------------------------- call Check( xggi, xgg ) call Check( lli, 'n', ll, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if nlat = xggi%nlat ! ll grid interpolated to gaussian lat: allocate( row(1:lli%im) ) ! loop over Gaussian latitudes: do j = 1, nlat ! number of lon points at this latitude: nlon = xggi%nlon(j) ! current Gaussian lat glat = xggi%lat(j) !! error if ll grid does not cover this gaussian lat !if ( glat < minval(lli%lat) .or. glat > maxval(lli%lat) ) then ! write (*,'("ERROR - ll grid does not cover current gg latitude:")') ! write (*,'("ERROR - ll lat range : ",2f12.4)') minval(lli%lat), maxval(lli%lat) ! write (*,'("ERROR - gg latitude : ",f12.4)') glat ! write (*,'("ERROR in ",a)') rname; status=1 !end if ! * interpolate ll grid to this gaussian lat: ! check for direction of ll lats: if ( lli%lat(1) < lli%lat(lli%jm) ) then ! 'normal' : grid stored from south to north ilast = 1 do i = 1, lli%im call Interp_Lin( lli%lat, ll(i,:), glat, row(i), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do else ! grid stored from north to south; fake with negative lats: ilast = 1 do i = 1, lli%im call Interp_Lin( -lli%lat, ll(i,:), -glat, row(i), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do end if ! base index of current row in xgg array: i0 = xggi%i1(j)-1 ! interpolate from ll lons to lons in this xgg row: period = 360.0*deg2rad ilast = 1 do i = 1, nlon ! curren lon in xgg grid: glon = (i-1)*xggi%dlon(j) ! periodic interpolation: call CircInterp_Lin( lli%lon, period, row, glon, xgg(i0+i), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do end do ! free memory deallocate( row ) ! ok status = 0 end subroutine Interpol_ll_gg ! ========================================================= ! === ! === area average gg ! === ! ========================================================= subroutine Aver_gg_ll( ggi, gg, lli, ll ) use Binas, only : deg2rad use Num, only : IntervalQuad_Lin, IntervalQuad_Cos_Lin use Grid_Type_gg, only : TggGridInfo, Check use Grid_Type_ll, only : TllGridInfo, Check ! --- in/out ------------------------------- type(TggGridInfo), intent(in) :: ggi real, intent(in) :: gg(ggi%np) type(TllGridInfo), intent(in) :: lli real, intent(out) :: ll(lli%im,lli%jm) integer :: status ! --- const ---------------------------------- character(len=*), parameter :: rname = mname//'/ Aver_gg_ll' ! --- local -------------------------------- integer :: nlon, nlon_max integer :: nlat integer :: i1, im real, allocatable :: lons(:) real, allocatable :: row(:) real, allocatable :: gl(:,:) real, allocatable :: gl_lat(:) real, allocatable :: gl_dim2(:) integer :: i, j integer :: ilast ! --- begin -------------------------------- call Check( ggi, gg ) call Check( lli, 'n', ll, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; stop; end if nlat = ggi%nlat nlon_max = maxval(ggi%nlon) ! ll in lon, gg in lat allocate( gl(size(ll,1),0:nlat+1) ) allocate( gl_dim2(0:nlat+1) ) ! latitudes from south->north allocate( gl_lat(0:nlat+1) ) ! row in gg grid; doubled from -360.0 to 360.0 allocate( lons(2*nlon_max) ) allocate( row(2*nlon_max) ) ! loop over Gaussian latitudes, from north to south! do j = 1, nlat ! number of lon points at this latitude: nlon = ggi%nlon(j) ! start and end indices in grid point array i1 = ggi%i1(j) im = ggi%im(j) ! lons from -360 to 360 do i = 1, nlon lons(i) = -360.0 + (i-1)*360.0/nlon lons(nlon+i) = (i-1)*360.0/nlon end do lons(1:2*nlon) = lons(1:2*nlon) * deg2rad ! grid values (doubled) row(1:2*nlon) = (/ gg(i1:im), gg(i1:im) /) ! integrate over dlon assuming linear interpolation; ! swap lats to ensure south -> north; ! result in [g] rad ilast = 1 do i = 1, lli%im call IntervalQuad_Lin( lons(1:2*nlon), row(1:2*nlon), & lli%blon(i-1), lli%blon(i), gl(i,nlat+1-j), & ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; stop; end if end do gl_lat(nlat+1-j) = ggi%lat(j) ! rad ! take care of poles (reverse lats); ! set polar value to average of surrounding points: if ( j == 1 ) then ! north pole gl_lat(nlat+1) = 90.0 * deg2rad gl(:,nlat+1) = lli%dlon * sum(gg(i1:im))/(im-i1+1) end if if ( j == nlat ) then ! south pole gl_lat(0) = -90.0 * deg2rad gl(:,0) = lli%dlon * sum(gg(i1:im))/(im-i1+1) end if end do ! integrate over dlat assuming linear interpolation; ! weight with cos(lat) to account for smaller cells near the poles, ! result in [g] rad^2 ilast = 1 do j = 1, lli%jm do i = 1, lli%im !>>> bsf15k bug fix >>>>>>>>>>>>>>>>>>>>>>>>>>>>>> !call IntervalQuad_Cos_Lin( gl_lat, gl(i,:), lli%blat(j-1), lli%blat(j), ll(i,j), ilast ) gl_dim2 = gl(i,:) call IntervalQuad_Cos_Lin( gl_lat, gl_dim2, lli%blat(j-1), lli%blat(j), & ll(i,j), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; stop; end if !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< end do end do ! area average; lli%area(j) for cells in row j in rad^2 do j = 1, lli%jm ll(:,j) = ll(:,j) / lli%area(j) end do ! free memory deallocate( gl ) deallocate( gl_dim2 ) deallocate( gl_lat ) deallocate( lons ) deallocate( row ) end subroutine Aver_gg_ll ! ========================================================= ! === ! === area average sh ! === ! ========================================================= ! Deterimine spectral truncation for area integration. integer function ShTruncation( T, dlon, dlat ) ! --- in/out -------------------------- integer, intent(in) :: T real, intent(in) :: dlon, dlat ! deg ! --- begin -------------------------- !print *, 'WARNING - spectral fields are not truncated ...' ! o no truncation: ShTruncation = T ! o choose minium T based on grid resolution: !ShTruncation = min( T, ceiling( (360.0/min(dlon,dlat))/2 - 1 ) ) !print *, 'WARNING - adhoc truncation of spectral fields for testing ...' !ShTruncation = 21 end function ShTruncation ! Deterimine refinement for averaging spectral fields ! over distance 'cellspacing' in degrees. ! Cell is divided in number of intervals returned by this function. integer function ShRefinement( T, cellspacing ) ! --- in/out ---------------------------------------- integer, intent(in) :: T real, intent(in) :: cellspacing ! --- local ----------------------------------------- real :: shres integer :: nstep ! --- beging ----------------------------------------- ! o fixec number of intervals per cell !ShRefinement = 20 !ShRefinement = 1 !ShRefinement = nint( cellspacing * 2 ) !write (*,'(" WARNING: ShRefinement = ",i4)') ShRefinement ! o resultion specified by spectral truncation ! truncation T <--> resolution 360.0/(2(T+1)) ! nstep points within resolution implied by T shres = 360.0 / (2*(T+1)) nstep = 2 ShRefinement = max( 1, ceiling(nstep*cellspacing/shres) ) end function ShRefinement ! *********************************** subroutine Aver_sh_ll( sh, lli, ll, status ) use grid_tools, only : deg2rad, pi use Num, only : IntervalQuad_Lin, IntervalQuad_Cos_Lin use Grid_Type_sh, only : TshGrid, Check, Eval_Lons use Grid_Type_ll, only : TllGridInfo, Check ! --- in/out ------------------------------- type(TshGrid), intent(in) :: sh type(TllGridInfo), intent(in) :: lli real, intent(out) :: ll(lli%im,lli%jm) integer, intent(out) :: status ! --- const ---------------------------------- character(len=*), parameter :: rname = mname//'/Aver_sh_ll' ! --- local -------------------------------- integer :: i, j, jf integer :: ilast integer :: nlon_fine integer :: T integer :: refinement_i, refinement_j real, allocatable :: llgrid(:) real, allocatable :: lons_fine(:), row_fine(:) real, allocatable :: llf(:,:), lat(:) ! --- begin -------------------------------- call Check( sh ) T = sh%T call Check( lli, 'n', ll, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; stop; end if ! determine refinement (5 points per sh resolution) refinement_i = ShRefinement( T, lli%dlon_deg ) refinement_j = ShRefinement( T, lli%dlon_deg ) ! number of lons in fine grid on complete circle: nlon_fine = 360.0/lli%dlon_deg * refinement_i ! store evaluation of spectral field: allocate( llgrid(nlon_fine) ) ! lons on complete circle from westb+[0,2pi) allocate( row_fine(0:nlon_fine) ) allocate( lons_fine(0:nlon_fine) ) do i = 0, nlon_fine lons_fine(i) = i*2*pi/nlon_fine end do lons_fine = lli%blon(0) + lons_fine ! ll in lon, fine in lat allocate( llf(lli%im,0:refinement_j) ) allocate( lat(0:refinement_j) ) ! loop over latitudes in ll grid: do j = 1, lli%jm ! loop over latitudes in fine grid: do jf = 0, refinement_j ! latitude in fine grid: lat(jf) = lli%blat(j-1) + jf*lli%dlat/refinement_j ! evaluate row: ! (oposite latitudes, but use only one) call Eval_Lons( llgrid, sh, lat(jf), & nlon_fine, lons_fine(0), nlon_fine, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! copy result, cyclic: row_fine = (/ llgrid(:), llgrid(1) /) ! integral in lon direction assuming linear interpolation, ! result in [sh] rad ilast = 1 do i = 1, lli%im call IntervalQuad_Lin( lons_fine, row_fine, lli%blon(i-1), lli%blon(i), llf(i,jf), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do end do ! integral in lat direction; ! weight with cos(lat) to account for smaller cells near the poles, ! result in [sh] rad^2 do i = 1, lli%im ilast = 1 call IntervalQuad_Cos_Lin( lat, llf(i,:), lli%blat(j-1), lli%blat(j), ll(i,j), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do ! area average; lli%area(j) for cells in row j in rad^2 ll(:,j) = ll(:,j) / lli%area(j) end do ! free memory deallocate( llgrid ) deallocate( lons_fine ) deallocate( row_fine ) deallocate( llf ) deallocate( lat ) end subroutine Aver_sh_ll ! =========================================================== ! ! key == 'aver' ! ! int int F(k) dA / A ! ! key == 'exp,aver' ! ! int int exp(F(k)) dA / A ! subroutine IntArea_shi_ll_f( key, shi, shc, lli, ll, status ) use grid_tools, only : deg2rad, pi use Num, only : IntervalQuad_Lin, IntervalQuad_Cos_Lin use Grid_Type_sh, only : TshGrid, TshGridInfo use Grid_Type_sh, only : Init, Done, Check, Set, Truncate, SpN use Grid_Type_sh, only : sh_Pnm, Eval_Lons use Grid_Type_ll, only : TllGridInfo, Check ! --- in/out ----------------------------------------- character(len=*), intent(in) :: key type(TshGridInfo), intent(in) :: shi complex, intent(in) :: shc(:) type(TllGridInfo), intent(in) :: lli real, intent(out) :: ll(lli%im,lli%jm) integer, intent(out) :: status ! --- const ---------------------------------- character(len=*), parameter :: rname = mname//'/IntArea_shi_ll_f' ! --- local ------------------------------------------ integer :: i, j, jf integer :: ilast type(TshGrid) :: sh integer :: Ttr real, allocatable :: Pnm(:) integer :: refinement_i, refinement_j integer :: nlon_fine_360, nlon_fine real, allocatable :: ff(:) real, allocatable :: lons_fine(:), row_fine(:) real, allocatable :: llf(:,:), lat(:) real, allocatable :: llf_dim2(:) !logical :: aver_to_prev ! --- begin ------------------------------------------ ! store input in old type grid: call Init( sh ) call Set( sh, shi%T, shc ) ! use truncation up to grid resolution: Ttr = ShTruncation( shi%T, lli%dlon_deg, lli%dlat_deg ) call Truncate( sh, Ttr ) ! allocate space for legendre coeff: allocate( Pnm(SpN(Ttr)) ) ! determine refinement (5 points per sh resolution) refinement_i = ShRefinement( Ttr, lli%dlon_deg ) refinement_j = ShRefinement( Ttr, lli%dlat_deg ) ! number of lons in fine grid on complete circle: nlon_fine_360 = 360.0/lli%dlon_deg * refinement_i nlon_fine = lli%im * refinement_i ! store evaluation of spectral field: allocate( ff(0:nlon_fine) ) ! lons on arc westb+[0,..) allocate( row_fine(0:nlon_fine) ) allocate( lons_fine(0:nlon_fine) ) do i = 0, nlon_fine lons_fine(i) = i*2*pi/nlon_fine_360 end do lons_fine = lli%blon(0) + lons_fine ! ll in lon, fine in lat allocate( llf(lli%im,0:refinement_j) ) allocate( llf_dim2(0:refinement_j) ) allocate( lat(0:refinement_j) ) ! loop over latitudes in ll grid: !aver_to_prev = .false. do j = 1, lli%jm ! loop over latitudes in fine grid: do jf = 0, refinement_j ! latitude in fine grid: lat(jf) = lli%blat(j-1) + jf*lli%dlat/refinement_j !! southpole ? !if ( abs(lat(jf) - (-pi/2)) < 1.0e-4 ) then ! ! fill with average of next row ! aver_to_prev = .true. ! cycle !end if ! !! northpole ? !if ( abs(lat(jf) - (pi/2)) < 1.0e-4 ) then ! ! fill with average of previous row: ! llf(:,jf) = sum(llf(:,jf-1)) / size(llf,1) ! exit !end if ! evaluate Legendre functions: call sh_Pnm( Pnm, Ttr, lat(jf), status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! evaluate rows: call Eval_Lons( ff, sh, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! combine fields: select case ( key ) ! ! int int F(k) dA / A ! case ( 'aver' ) row_fine = ff / lli%area(j) ! ! int int exp(F(k)) dA / A ! case ( 'exp,aver' ) row_fine = exp(ff) / lli%area(j) ! ! error ... ! case default write (*,'("ERROR - unknown key `",a,"`")') trim(key) write (*,'("ERROR in ",a)') rname; status=1; return end select ! integral in lon direction assuming linear interpolation, ! result in [sh] rad ilast = 1 do i = 1, lli%im call IntervalQuad_Lin( lons_fine, row_fine, lli%blon(i-1), lli%blon(i), llf(i,jf), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do !! copy to southpole .. !if ( (jf==1) .and. aver_to_prev ) then ! llf(:,0) = sum(llf(:,1)) / size(llf,1) ! aver_to_prev = .false. !end if end do ! loop over rows in fine grid ! integral in lat direction; ! weight with cos(lat) to account for smaller grid cells: do i = 1, lli%im ilast = 1 !call IntervalQuad_Cos_Lin( lat, llf(i,:,l), lli%blat(j-1), lli%blat(j), ll(i,j,l), ilast ) !>>> bsf15k bug fix >>>>>>>>>>>>>>>>>>>>>>>>>>>>>> llf_dim2 = llf(i,:) call IntervalQuad_Cos_Lin( lat, llf_dim2, lli%blat(j-1), lli%blat(j), ll(i,j), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< end do end do ! loop over rows ! free memory deallocate( Pnm ) deallocate( ff ) deallocate( lons_fine ) deallocate( row_fine ) deallocate( llf ) deallocate( llf_dim2 ) deallocate( lat ) call Done( sh ) end subroutine IntArea_shi_ll_f ! *** subroutine IntArea_sh_ll_f( key, F, lli, ll, status ) use grid_tools, only : deg2rad, pi use Num, only : IntervalQuad_Lin, IntervalQuad_Cos_Lin use Grid_Type_sh, only : TshGrid, Init, Done, SpN, sh_Pnm, Eval_Lons, Set, Check use Grid_Type_ll, only : TllGridInfo, Check ! --- in/out ----------------------------------------- character(len=*), intent(in) :: key type(TshGrid), intent(in) :: F(:) type(TllGridInfo), intent(in) :: lli real, intent(out) :: ll(lli%im,lli%jm,size(F)) integer, intent(out) :: status ! --- const -------------------------------- character(len=*), parameter :: rname = mname//'/IntArea_sh_ll_f' ! --- local ------------------------------------------ integer :: i, j, jf integer :: l, lm integer :: ilast integer :: T type(TshGrid) :: sh real, allocatable :: Pnm(:) integer :: refinement_i, refinement_j integer :: nlon_fine_360, nlon_fine real, allocatable :: ff(:) real, allocatable :: lons_fine(:), row_fine(:) real, allocatable :: llf(:,:,:), lat(:) real, allocatable :: llf_dim2(:) ! --- begin ------------------------------------------ ! number of levels: lm = size(F) ! check arguments: T = F(1)%T do l = 2, lm call Check( F(l), T ) end do ! use truncation up to grid resolution: T = ShTruncation( T, lli%dlon_deg, lli%dlat_deg ) call Init( sh, T ) ! allocate space for legendre coeff: allocate( Pnm(SpN(T)) ) ! determine refinement (5 points per sh resolution) refinement_i = ShRefinement( T, lli%dlon_deg ) refinement_j = ShRefinement( T, lli%dlat_deg ) ! number of lons in fine grid on complete circle: nlon_fine_360 = 360.0/lli%dlon_deg * refinement_i nlon_fine = lli%im * refinement_i ! store evaluation of spectral field: allocate( ff(0:nlon_fine) ) ! lons on arc westb+[0,..) allocate( row_fine(0:nlon_fine) ) allocate( lons_fine(0:nlon_fine) ) do i = 0, nlon_fine lons_fine(i) = i*2*pi/nlon_fine_360 end do lons_fine = lli%blon(0) + lons_fine ! ll in lon, fine in lat allocate( llf(lli%im,0:refinement_j,lm) ) allocate( llf_dim2(0:refinement_j) ) allocate( lat(0:refinement_j) ) ! loop over latitudes in ll grid: !aver_to_prev = .false. do j = 1, lli%jm ! loop over latitudes in fine grid: do jf = 0, refinement_j ! latitude in fine grid: lat(jf) = lli%blat(j-1) + jf*lli%dlat/refinement_j !! southpole ? !if ( abs(lat(jf) - (-pi/2)) < 1.0e-4 ) then ! ! fill with average of next row ! aver_to_prev = .true. ! cycle !end if ! !! northpole ? !if ( abs(lat(jf) - (pi/2)) < 1.0e-4 ) then ! ! fill with average of previous row: ! do l = 1, lm ! llf(:,jf,l) = sum(llf(:,jf-1,l)) / size(llf,1) ! end do ! exit !end if ! evaluate Legendre functions: call sh_Pnm( Pnm, T, lat(jf), status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! evaluate rows: do l = 1, lm call Set( sh, T, F(l) ) call Eval_Lons( ff, sh, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! combine fields: select case ( key ) ! ! int int F(k) dA / A ! case ( 'aver' ) row_fine = ff / lli%area(j) ! ! int int exp(F(k)) dA / A ! case ( 'exp,aver' ) row_fine = exp(ff) / lli%area(j) ! ! error ... ! case default write (*,'("ERROR - unsupported integrand key : ",a)') trim(key) write (*,'("ERROR in ",a)') rname; status=1; return end select ! integral in lon direction assuming linear interpolation, ! result in [sh] rad ilast = 1 do i = 1, lli%im call IntervalQuad_Lin( lons_fine, row_fine, lli%blon(i-1), lli%blon(i), llf(i,jf,l), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do end do ! loop over levels !! copy to southpole .. !if ( (jf==1) .and. aver_to_prev ) then ! do l = 1, lm ! llf(:,0,l) = sum(llf(:,1,l)) / size(llf,1) ! end do ! aver_to_prev = .false. !end if end do ! loop over rows in fine grid ! integral in lat direction; ! weight with cos(lat) to account for smaller grid cells: do l = 1, lm do i = 1, lli%im ilast = 1 !call IntervalQuad_Cos_Lin( lat, llf(i,:,l), lli%blat(j-1), lli%blat(j), ll(i,j,l), ilast ) !>>> bsf15k bug fix >>>>>>>>>>>>>>>>>>>>>>>>>>>>>> llf_dim2 = llf(i,:,l) call IntervalQuad_Cos_Lin( lat, llf_dim2, lli%blat(j-1), lli%blat(j), ll(i,j,l), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< end do end do end do ! loop over rows ! free memory deallocate( Pnm ) deallocate( ff ) deallocate( lons_fine ) deallocate( row_fine ) deallocate( llf ) deallocate( llf_dim2 ) deallocate( lat ) call Done( sh ) ! ok status = 0 end subroutine IntArea_sh_ll_f ! ! key == 'F*(da+db*exp(H))*cos' ! ! int int F(k) (da+db*exp(H)) cos(lat) dA ! ! key == 'F*G*(db*exp(H))/cos' ! ! int int F(k) G db*exp(H) / cos(lat) dA ! ! (for these keys, results are added to ll) ! ! subroutine IntArea_sh_ll_fgh( key, F, G, H, da, db, lli, ll, status ) use grid_tools, only : deg2rad, pi use Num, only : IntervalQuad_Lin use Grid_Type_sh, only : TshGrid, Init, Done, SpN, sh_Pnm, Eval_Lons, Set, Check use Grid_Type_ll, only : TllGridInfo, Check ! --- in/out ----------------------------------------- character(len=*), intent(in) :: key type(TshGrid), intent(in) :: F(:) type(TshGrid), intent(in) :: G, H real, intent(in) :: da(size(F)), db(size(F)) type(TllGridInfo), intent(in) :: lli real, intent(inout) :: ll(lli%im,lli%jm,size(F)) integer, intent(out) :: status ! --- const -------------------------------------- character(len=*), parameter :: rname = mname//'/IntArea_sh_ll_fgh' ! --- local ------------------------------------------ integer :: i, j, jf integer :: l, lm integer :: ilast integer :: T type(TshGrid) :: sh real, allocatable :: Pnm(:) integer :: refinement_i, refinement_j integer :: nlon_fine_360, nlon_fine real, allocatable :: ff(:) real, allocatable :: gg(:) real, allocatable :: exp_hh(:) real, allocatable :: lons_fine(:), row_fine(:) real, allocatable :: llf(:,:,:), lat(:) real, allocatable :: llf_dim2(:) logical :: aver_to_prev real :: res ! --- begin ------------------------------------------ ! number of levels: lm = size(F) ! check arguments: call Check( H ) T = H%T do l = 1, lm call Check( F(l), T ) end do call Check( G, T ) ! use truncation up to grid resolution: T = ShTruncation( T, lli%dlon_deg, lli%dlat_deg ) call Init( sh, T ) ! allocate space for legendre coeff: allocate( Pnm(SpN(T)) ) ! determine refinement (5 points per sh resolution) refinement_i = ShRefinement( T, lli%dlon_deg ) refinement_j = ShRefinement( T, lli%dlat_deg ) ! number of lons in fine grid on complete circle: nlon_fine_360 = 360.0/lli%dlon_deg * refinement_i nlon_fine = lli%im * refinement_i ! store evaluation of spectral field: allocate( ff(0:nlon_fine) ) allocate( gg(0:nlon_fine) ) allocate( exp_hh(0:nlon_fine) ) ! lons on arc westb+[0,..) allocate( row_fine(0:nlon_fine) ) allocate( lons_fine(0:nlon_fine) ) do i = 0, nlon_fine lons_fine(i) = i*2*pi/nlon_fine_360 end do lons_fine = lli%blon(0) + lons_fine ! ll in lon, fine in lat allocate( llf(lli%im,0:refinement_j,lm) ) allocate( lat(0:refinement_j) ) allocate( llf_dim2(0:refinement_j) ) ! *** integrals in lon direction ! loop over latitudes in ll grid: aver_to_prev = .false. do j = 1, lli%jm ! loop over latitudes in fine grid: do jf = 0, refinement_j ! latitude in fine grid: lat(jf) = lli%blat(j-1) + jf*lli%dlat/refinement_j ! southpole ? if ( abs(lat(jf) - (-pi/2)) < 1.0e-4 ) then ! fill with average of next row aver_to_prev = .true. cycle end if ! northpole ? if ( abs(lat(jf) - (pi/2)) < 1.0e-4 ) then ! fill with average of previous row: do l = 1, lm llf(:,jf,l) = sum(llf(:,jf-1,l)) / size(llf,1) end do exit end if ! evaluate Legendre functions: call sh_Pnm( Pnm, T, lat(jf), status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! evaluate rows: call Set( sh, T, G ) call Eval_Lons( gg, sh, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if call Set( sh, T, H ) call Eval_Lons( exp_hh, sh, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if exp_hh = exp(exp_hh) do l = 1, lm ! combine fields: select case ( key ) ! ! int int F(k) (da+db*exp(H)) cos(lat) dA ! case ( 'F*(da+db*exp(H))*cos' ) call Set( sh, T, F(l) ) call Eval_Lons( ff, sh, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if row_fine = ff * ( da(l) + exp_hh*db(l) ) * cos(lat(jf)) ! ! int int F(k) G db*exp(H) / cos(lat) dA ! case ( 'F*G*(db*exp(H))/cos' ) if ( db(l) == 0.0 ) then row_fine = 0.0 else call Set( sh, T, F(l) ) call Eval_Lons( ff, sh, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if row_fine = ff * gg * ( exp_hh*db(l) ) / cos(lat(jf)) end if ! ! error ... ! case default print *, 'IntArea_sh_ll_fgh - unknown key "'//trim(key)//'"' stop end select ! integral in lon direction assuming linear interpolation, ! result in [sh] rad ilast = 1 do i = 1, lli%im call IntervalQuad_Lin( lons_fine, row_fine, lli%blon(i-1), lli%blon(i), llf(i,jf,l), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do end do ! loop over levels ! copy to southpole .. if ( (jf==1) .and. aver_to_prev ) then do l = 1, lm llf(:,0,l) = sum(llf(:,1,l)) / size(llf,1) end do aver_to_prev = .false. end if end do ! loop over rows in fine grid ! *** integral in lat direction ! 3D field: do l = 1, lm ilast = 1 do i = 1, lli%im !call lquad( lat, llf(i,:,l), lli%blat(j-1), lli%blat(j), res, ilast ) !>>> bsf15k bug fix >>>>>>>>>>>>>>>>>>>>>>>>>>>>>> llf_dim2 = llf(i,:,l) call IntervalQuad_Lin( lat, llf_dim2, lli%blat(j-1), lli%blat(j), res, ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! add result ll(i,j,l) = ll(i,j,l) + res end do end do end do ! loop over rows ! free memory deallocate( Pnm ) deallocate( ff ) deallocate( gg ) deallocate( exp_hh ) deallocate( lons_fine ) deallocate( row_fine ) deallocate( llf ) deallocate( llf_dim2 ) deallocate( lat ) call Done( sh ) ! ok status = 0 end subroutine IntArea_sh_ll_fgh ! *** subroutine IntArea_shi_ll_fgh( key, shi_in, nlev, F, G, H, da, db, lli, ll, status ) use grid_tools, only : deg2rad, pi use Num, only : IntervalQuad_Lin use Grid_Type_sh, only : TshGridInfo, Init, Done, sh_Pnm, Eval_Lons, Set, Check use Grid_Type_ll, only : TllGridInfo, Check ! --- in/out ----------------------------------------- character(len=*), intent(in) :: key type(TshGridInfo), intent(in) :: shi_in integer, intent(in) :: nlev complex, intent(in) :: F(shi_in%np,nlev) complex, intent(in) :: G(shi_in%np) complex, intent(in) :: H(shi_in%np) real, intent(in) :: da(nlev), db(nlev) type(TllGridInfo), intent(in) :: lli real, intent(inout) :: ll(lli%im,lli%jm,nlev) integer, intent(out) :: status ! --- const -------------------------------------- character(len=*), parameter :: rname = mname//'/IntArea_sh_ll_fgh' ! --- local ------------------------------------------ integer :: i, j, jf integer :: l integer :: ilast type(TshGridInfo) :: shi integer :: T real, pointer :: Pnm(:) integer :: refinement_i, refinement_j integer :: nlon_fine_360, nlon_fine real, pointer :: ff(:) real, pointer :: gg(:) real, pointer :: exp_hh(:) real, pointer :: lons_fine(:), row_fine(:) real, pointer :: llf(:,:,:), lat(:) real, pointer :: llf_dim2(:) logical :: aver_to_prev real :: res ! --- begin ------------------------------------------ ! use truncation up to grid resolution: T = ShTruncation( shi_in%T, lli%dlon_deg, lli%dlat_deg ) !T = ShTruncation( 159, lli%dlon_deg, lli%dlat_deg ) ! spectral info: call Init( shi, T, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! allocate space for legendre coeff: allocate( Pnm(shi%np) ) ! determine refinement (5 points per sh resolution) refinement_i = ShRefinement( shi%T, lli%dlon_deg ) refinement_j = ShRefinement( shi%T, lli%dlat_deg ) ! number of lons in fine grid on complete circle: nlon_fine_360 = 360.0/lli%dlon_deg * refinement_i nlon_fine = lli%im * refinement_i ! store evaluation of spectral field: allocate( ff(0:nlon_fine) ) allocate( gg(0:nlon_fine) ) allocate( exp_hh(0:nlon_fine) ) ! lons on arc westb+[0,..) allocate( row_fine(0:nlon_fine) ) allocate( lons_fine(0:nlon_fine) ) do i = 0, nlon_fine lons_fine(i) = i*2*pi/nlon_fine_360 end do lons_fine = lli%blon(0) + lons_fine ! ll in lon, fine in lat allocate( llf(lli%im,0:refinement_j,nlev) ) allocate( lat(0:refinement_j) ) allocate( llf_dim2(0:refinement_j) ) ! *** integrals in lon direction ! loop over latitudes in ll grid: aver_to_prev = .false. do j = 1, lli%jm ! loop over latitudes in fine grid: do jf = 0, refinement_j ! latitude in fine grid: lat(jf) = lli%blat(j-1) + jf*lli%dlat/refinement_j ! southpole ? if ( abs(lat(jf) - (-pi/2)) < 1.0e-4 ) then ! fill with average of next row aver_to_prev = .true. cycle end if ! northpole ? if ( abs(lat(jf) - (pi/2)) < 1.0e-4 ) then ! fill with average of previous row: do l = 1, nlev llf(:,jf,l) = sum(llf(:,jf-1,l)) / size(llf,1) end do exit end if ! evaluate Legendre functions: call sh_Pnm( Pnm, shi%T, lat(jf), status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! evaluate rows: !print *, 'aaa10b shi_in ', shi_in%T, shi_in%np !print *, 'aaa10b G size ', size(G) !print *, 'aaa10b G val ', G(1:10) !print *, 'aaa10b shi ', shi%T, shi%np !print *, 'aaa10b Pnm size ', size(Pnm) !print *, 'aaa10b Pnm val ', Pnm(1:10) call Eval_Lons( shi_in, G, shi, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, gg, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if call Eval_Lons( shi_in, H, shi, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, exp_hh, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if exp_hh = exp(exp_hh) ! loop over levels: !xOMP PARALLEL & !xOMP default ( none ) & !xOMP shared ( nlev, key, T, F, Pnm, nlon_fine_360, nlon_fine, lons_fine ) & !xOMP shared ( ff, gg, da, exp_hh, db, lat, jf, lli, llf ) & !xOMP private ( l, sh, row_fine, ilast, status ) !xOMP DO do l = 1, nlev ! combine fields: select case ( key ) ! ! int int F(k) (da+db*exp(H)) cos(lat) dA ! case ( 'F*(da+db*exp(H))*cos' ) call Eval_Lons( shi_in, F(:,l), shi, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, ff, status ) !if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if row_fine = ff * ( da(l) + exp_hh*db(l) ) * cos(lat(jf)) ! ! int int F(k) G db*exp(H) / cos(lat) dA ! case ( 'F*G*(db*exp(H))/cos' ) if ( db(l) == 0.0 ) then row_fine = 0.0 else call Eval_Lons( shi_in, F(:,l), shi, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, ff, status ) !if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if row_fine = ff * gg * ( exp_hh*db(l) ) / cos(lat(jf)) end if ! ! error ... ! case default print *, 'IntArea_sh_ll_fgh - unknown key "'//trim(key)//'"' !stop end select ! integral in lon direction assuming linear interpolation, ! result in [sh] rad ilast = 1 do i = 1, lli%im call IntervalQuad_Lin( lons_fine, row_fine, lli%blon(i-1), lli%blon(i), llf(i,jf,l), ilast, status ) !if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do end do ! levels !xOMP END DO !xOMP END PARALLEL ! copy to southpole .. if ( (jf==1) .and. aver_to_prev ) then do l = 1, nlev llf(:,0,l) = sum(llf(:,1,l)) / size(llf,1) end do aver_to_prev = .false. end if end do ! loop over rows in fine grid ! *** integral in lat direction ! loop over levels: !xOMP PARALLEL DO do l = 1, nlev ilast = 1 do i = 1, lli%im !call lquad( lat, llf(i,:,l), lli%blat(j-1), lli%blat(j), res, ilast ) !>>> bsf15k bug fix >>>>>>>>>>>>>>>>>>>>>>>>>>>>>> llf_dim2 = llf(i,:,l) call IntervalQuad_Lin( lat, llf_dim2, lli%blat(j-1), lli%blat(j), res, ilast, status ) !if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! add result ll(i,j,l) = ll(i,j,l) + res end do end do !xOMP END PARALLEL DO end do ! loop over rows ! free memory deallocate( Pnm ) deallocate( ff ) deallocate( gg ) deallocate( exp_hh ) deallocate( lons_fine ) deallocate( row_fine ) deallocate( llf ) deallocate( llf_dim2 ) deallocate( lat ) call Done( shi )!, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! ok status = 0 end subroutine IntArea_shi_ll_fgh ! *** ! ! key == '[F*(da+db*exp(H))*cos]/[()*cos]' ! ! int int F(k) (da+db*exp(H)) cos(lat) dA / int int (da+db*exp(H)) cos(lat) dA ! ! Result is put in ll, not added. ! Uses cos integration in lat direction. ! subroutine IntArea_sh_ll_fh( key, F, H, da, db, lli, ll, status ) use grid_tools, only : deg2rad, pi use Num, only : IntervalQuad_Lin, IntervalQuad_Cos_Lin use Grid_Type_sh, only : TshGrid, Init, Done, SpN, sh_Pnm, Eval_Lons, Set, Check use Grid_Type_ll, only : TllGridInfo, Check ! --- in/out ----------------------------------------- character(len=*), intent(in) :: key type(TshGrid), intent(in) :: F(:) type(TshGrid), intent(in) :: H real, intent(in) :: da(size(F)), db(size(F)) type(TllGridInfo), intent(in) :: lli real, intent(inout) :: ll(lli%im,lli%jm,size(F)) integer, intent(out) :: status ! --- const -------------------------------------- character(len=*), parameter :: rname = mname//'/IntArea_sh_ll_fh' ! --- local ------------------------------------------ integer :: i, j, jf integer :: l, lm integer :: ilast integer :: T type(TshGrid) :: sh real, allocatable :: Pnm(:) integer :: refinement_i, refinement_j integer :: nlon_fine_360, nlon_fine real, allocatable :: ff(:) real, allocatable :: exp_hh(:) real, allocatable :: lons_fine(:), row_fine(:) real, allocatable :: llf(:,:,:), lat(:) real, allocatable :: llf_dim2(:) real :: res real, allocatable :: llf_expH(:,:) ! --- begin ------------------------------------------ ! number of levels: lm = size(F) ! check arguments: call Check( H ) T = H%T do l = 1, lm call Check( F(L), T ) end do T = ShTruncation( T, lli%dlon_deg, lli%dlat_deg ) call Init( sh, T ) ! allocate space for legendre coeff: allocate( Pnm(SpN(T)) ) ! determine refinement (5 points per sh resolution) refinement_i = ShRefinement( T, lli%dlon_deg ) refinement_j = ShRefinement( T, lli%dlat_deg ) ! number of lons in fine grid on complete circle: nlon_fine_360 = 360.0/lli%dlon_deg * refinement_i nlon_fine = lli%im * refinement_i ! store evaluation of spectral field: allocate( ff(0:nlon_fine) ) allocate( exp_hh(0:nlon_fine) ) ! lons on arc westb+[0,..) allocate( row_fine(0:nlon_fine) ) allocate( lons_fine(0:nlon_fine) ) do i = 0, nlon_fine lons_fine(i) = i*2*pi/nlon_fine_360 end do lons_fine = lli%blon(0) + lons_fine ! ll in lon, fine in lat allocate( llf(lli%im,0:refinement_j,lm) ) allocate( llf_expH(lli%im,0:refinement_j) ) allocate( lat(0:refinement_j) ) allocate( llf_dim2(0:refinement_j) ) ! *** integrals in lon direction ! loop over latitudes in ll grid: do j = 1, lli%jm ! loop over latitudes in fine grid: do jf = 0, refinement_j ! latitude in fine grid: lat(jf) = lli%blat(j-1) + jf*lli%dlat/refinement_j ! evaluate Legendre functions: call sh_Pnm( Pnm, T, lat(jf), status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! evaluate rows: call Set( sh, T, H ) call Eval_Lons( exp_hh, sh, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! exponent: exp_hh = exp(exp_hh) ! compute lon integral over exp(H), as a part of integral over second term ! (use linear interpolation): ilast = 1 do i = 1, lli%im call IntervalQuad_Lin( lons_fine, exp_hh, lli%blon(i-1), lli%blon(i), llf_expH(i,jf), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do do l = 1, lm call Set( sh, T, F(l) ) call Eval_Lons( ff, sh, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! combine fields: select case ( key ) ! ! int int F(k) (da+db*exp(H)) cos(lat) dA / int int (da+db*exp(H)) cos(lat) dA ! case ( '[F*(da+db*exp(H))*cos]/[()*cos]' ) ! lon integral of first term: row_fine = ff * ( da(l) + exp_hh*db(l) ) ! ! error ... ! case default print *, 'IntArea_sh_ll_fh - unknown key "'//trim(key)//'"' stop end select ! integral in lon direction assuming linear interpolation, ! result in [sh] rad ilast = 1 do i = 1, lli%im call IntervalQuad_Lin( lons_fine, row_fine, lli%blon(i-1), lli%blon(i), llf(i,jf,l), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do end do ! loop over levels end do ! loop over rows in fine grid ! *** integral in lat direction ! 3D field: do l = 1, lm ilast = 1 do i = 1, lli%im !call IntervalQuad_Cos_Lin( lat, llf(i,:,l), lli%blat(j-1), lli%blat(j), res, ilast ) !>>> bsf15k bug fix >>>>>>>>>>>>>>>>>>>>>>>>>>>>>> llf_dim2 = llf(i,:,l) call IntervalQuad_Cos_Lin( lat, llf_dim2, lli%blat(j-1), lli%blat(j), res, ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! set result ll(i,j,l) = res end do end do ! weight with integral over exp_hh : ilast = 1 do i = 1, lli%im !>>> bsf15k bug fix >>>>>>>>>>>>>>>>>>>>>>>>>>>>>> !call IntervalQuad_Cos_Lin( lat, llf_expH(i,:), lli%blat(j-1), lli%blat(j), res, ilast ) llf_dim2 = llf_expH(i,:) call IntervalQuad_Cos_Lin( lat, llf_dim2, lli%blat(j-1), lli%blat(j), res, ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! weight with int int (da+db*exp(H))*cos(lat) dA do l = 1, lm ll(i,j,l) = ll(i,j,l) / (da(l)*lli%area(j) + db(l)*res) end do end do end do ! loop over rows ! free memory deallocate( Pnm ) deallocate( ff ) deallocate( exp_hh ) deallocate( lons_fine ) deallocate( row_fine ) deallocate( llf ) deallocate( llf_dim2 ) deallocate( llf_expH ) deallocate( lat ) call Done( sh ) ! ok status = 0 end subroutine IntArea_sh_ll_fh ! *** subroutine IntArea_shi_ll_fh( key, shi_in, nlev, F, H, da, db, lli, ll, status ) use grid_tools, only : deg2rad, pi use Num, only : IntervalQuad_Lin, IntervalQuad_Cos_Lin use Grid_Type_sh, only : TshGrid, Init, Done use Grid_Type_sh, only : SpN, sh_Pnm, Eval_Lons, Set, Check, SpN, Truncate use Grid_Type_sh, only : TshGridInfo use Grid_Type_ll, only : TllGridInfo, Check ! --- in/out ----------------------------------------- character(len=*), intent(in) :: key type(TshGridInfo), intent(in) :: shi_in integer, intent(in) :: nlev complex, intent(in) :: F(shi_in%np,nlev) complex, intent(in) :: H(shi_in%np) real, intent(in) :: da(nlev), db(nlev) type(TllGridInfo), intent(in) :: lli real, intent(inout) :: ll(lli%im,lli%jm,nlev) integer, intent(out) :: status ! --- const ------------------------------------------ character(len=*), parameter :: rname = mname//'/IntArea_shi_ll_fh' ! --- local ------------------------------------------ integer :: i, j, jf integer :: l !, lm integer :: ilast type(TshGridInfo) :: shi ! type(TshGrid) :: sh integer :: Ttr real, allocatable :: Pnm(:) integer :: refinement_i, refinement_j integer :: nlon_fine_360, nlon_fine real, allocatable :: ff(:) real, allocatable :: exp_hh(:) real, allocatable :: lons_fine(:), row_fine(:) real, allocatable :: llf(:,:,:), lat(:) real, allocatable :: llf_dim2(:) real :: res real, allocatable :: llf_expH(:,:) ! --- begin ------------------------------------------ ! ! number of levels: ! lm = size(F,2) ! ! ! check arguments: ! if ( (size(F,1) /= shi%np) .or. (size(H) /= shi%np) ) then ! write (*,'("ERROR - number of complex coeff does not match with sh grid definition:")') ! write (*,'("ERROR - shi%np : ",i6)') shi%np ! write (*,'("ERROR - size(F,1) : ",i6)') size(F,1) ! write (*,'("ERROR - size(H) : ",i6)') size(H) ! write (*,'("ERROR in ",a)') rname; status=1; return ! end if ! ! ! input temporary stored in old type grid: ! call Init( sh, shi%T ) ! use truncation up to grid resolution: Ttr = ShTruncation( shi_in%T, lli%dlon_deg, lli%dlat_deg ) ! call Truncate( sh, Ttr ) call Init( shi, Ttr, status ) ! allocate space for legendre coeff: ! allocate( Pnm(SpN(Ttr)) ) allocate( Pnm(shi%np) ) ! determine refinement (5 points per sh resolution) refinement_i = ShRefinement( shi%T, lli%dlon_deg ) refinement_j = ShRefinement( shi%T, lli%dlat_deg ) ! number of lons in fine grid on complete circle: nlon_fine_360 = 360.0/lli%dlon_deg * refinement_i nlon_fine = lli%im * refinement_i ! store evaluation of spectral field: allocate( ff(0:nlon_fine) ) allocate( exp_hh(0:nlon_fine) ) ! lons on arc westb+[0,..) allocate( row_fine(0:nlon_fine) ) allocate( lons_fine(0:nlon_fine) ) do i = 0, nlon_fine lons_fine(i) = i*2*pi/nlon_fine_360 end do lons_fine = lli%blon(0) + lons_fine ! ll in lon, fine in lat !allocate( llf(lli%im,0:refinement_j,lm) ) allocate( llf(lli%im,0:refinement_j,nlev) ) allocate( llf_expH(lli%im,0:refinement_j) ) allocate( lat(0:refinement_j) ) allocate( llf_dim2(0:refinement_j) ) ! *** integrals in lon direction ! loop over latitudes in ll grid: do j = 1, lli%jm ! loop over latitudes in fine grid: do jf = 0, refinement_j ! latitude in fine grid: lat(jf) = lli%blat(j-1) + jf*lli%dlat/refinement_j ! evaluate Legendre functions: !call sh_Pnm( Pnm, Ttr, lat(jf), status ) call sh_Pnm( Pnm, shi%T, lat(jf), status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! evaluate rows: !call Set( sh, shi%T, H ) !call Truncate( sh, Ttr ) !call Eval_Lons( exp_hh, sh, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, status ) !if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if call Eval_Lons( shi_in, H, shi, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, exp_hh, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if exp_hh = exp(exp_hh) ! compute lon integral over exp(H), as a part of integral over second term ! (use linear interpolation): ilast = 1 do i = 1, lli%im call IntervalQuad_Lin( lons_fine, exp_hh, lli%blon(i-1), lli%blon(i), llf_expH(i,jf), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do ! combine fields: select case ( key ) ! ! int int F(k) (da+db*exp(H)) cos(lat) dA / int int (da+db*exp(H)) cos(lat) dA ! case ( '[F*(da+db*exp(H))*cos]/[()*cos]' ) !$OMP PARALLEL & !$OMP default ( none ) & !$OMP shared ( shi_in, F, shi, Pnm, nlon_fine_360, lons_fine, nlon_fine, lli, llf ) & !$OMP shared ( exp_hh, da, db, jf ) & !$OMP shared ( nlev ) & !$OMP private ( l, ff, status, row_fine, ilast ) !$OMP DO do l = 1, nlev !do l = 1, lm !call Set( sh, shi%T, F(:,l) ) !call Truncate( sh, Ttr ) call Eval_Lons( shi_in, F(:,l), shi, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, ff, status ) !if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! lon integral of first term: row_fine = ff * ( da(l) + exp_hh*db(l) ) ! integral in lon direction assuming linear interpolation, ! result in [sh] rad ilast = 1 do i = 1, lli%im call IntervalQuad_Lin( lons_fine, row_fine, lli%blon(i-1), lli%blon(i), llf(i,jf,l), ilast, status ) !if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do end do ! loop over levels !$OMP END DO !$OMP END PARALLEL ! ! error ... ! case default print *, 'IntArea_shi_ll_fh - unknown key "'//trim(key)//'"' stop end select end do ! loop over rows in fine grid ! *** integral in lat direction ! 3D field: do l = 1, nlev !do l = 1, lm ilast = 1 do i = 1, lli%im !call IntervalQuad_Cos_Lin( lat, llf(i,:,l), lli%blat(j-1), lli%blat(j), res, ilast ) !>>> bsf15k bug fix >>>>>>>>>>>>>>>>>>>>>>>>>>>>>> llf_dim2 = llf(i,:,l) call IntervalQuad_Cos_Lin( lat, llf_dim2, lli%blat(j-1), lli%blat(j), res, ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! set result ll(i,j,l) = res end do end do ! weight with integral over exp_hh : ilast = 1 do i = 1, lli%im !>>> bsf15k bug fix >>>>>>>>>>>>>>>>>>>>>>>>>>>>>> !call IntervalQuad_Cos_Lin( lat, llf_expH(i,:), lli%blat(j-1), lli%blat(j), res, ilast ) llf_dim2 = llf_expH(i,:) call IntervalQuad_Cos_Lin( lat, llf_dim2, lli%blat(j-1), lli%blat(j), res, ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< ! weight with int int (da+db*exp(H))*cos(lat) dA do l = 1, nlev !do l = 1, lm ll(i,j,l) = ll(i,j,l) / (da(l)*lli%area(j) + db(l)*res) end do end do end do ! loop over rows ! free memory deallocate( Pnm ) deallocate( ff ) deallocate( exp_hh ) deallocate( lons_fine ) deallocate( row_fine ) deallocate( llf ) deallocate( llf_dim2 ) deallocate( llf_expH ) deallocate( lat ) !call Done( sh ) ! ok status = 0 end subroutine IntArea_shi_ll_fh ! *** subroutine IntLat_sh_ll( key, F, H, da, db, lli, ll_u, status ) use grid_tools, only : deg2rad, pi use Num, only : IntervalQuad_Lin use Grid_Type_sh, only : TshGrid, Init, Done, Set, Check use Grid_Type_sh, only : sh_Pnm, Eval_Lons, SpN use Grid_Type_ll, only : TllGridInfo, Check ! --- in/out ----------------------------------------- character(len=*), intent(in) :: key type(TshGrid), intent(in) :: F(:) type(TshGrid), intent(in) :: H real, intent(in) :: da(size(F)) real, intent(in) :: db(size(F)) type(TllGridInfo), intent(in) :: lli real, intent(inout) :: ll_u(0:lli%im,lli%jm,size(F)) integer, intent(out) :: status ! --- const -------------------------------------- character(len=*), parameter :: rname = mname//'/IntLat_sh_ll' ! --- local ------------------------------------------ integer :: i, j, jf integer :: l, lm integer :: ilast integer :: T type(TshGrid) :: sh real, allocatable :: Pnm(:) integer :: refinement_j integer :: nlon_fine_360, nlon_fine real, allocatable :: ff(:) real, allocatable :: exp_hh(:) real, allocatable :: llf(:,:,:), lat(:), llf_row(:) logical :: aver_to_prev ! --- begin ------------------------------------------ ! number of levels: lm = size(F) ! check arguments: call Check( H ) T = H%T do l = 1, lm call Check( F(l), T ) end do T = ShTruncation( T, lli%dlon_deg, lli%dlat_deg ) call Init( sh, T ) ! allocate space for legendre coeff: allocate( Pnm(SpN(T)) ) ! determine refinement based on spectral resolution ! and length of required integration interval: refinement_j = ShRefinement( T, lli%dlat_deg ) ! number of lons in fine grid on complete circle: nlon_fine_360 = 360.0/lli%dlon_deg nlon_fine = lli%im ! store evaluation of spectral field: allocate( ff(0:nlon_fine) ) allocate( exp_hh(0:nlon_fine) ) ! ll in lon, fine in lat allocate( llf(0:lli%im,0:refinement_j,lm) ) allocate( llf_row(0:refinement_j) ) allocate( lat(0:refinement_j) ) ! loop over latitudes in ll grid: aver_to_prev = .false. do j = 1, lli%jm ! loop over latitudes in fine grid: do jf = 0, refinement_j ! latitude in fine grid: lat(jf) = lli%blat(j-1) + jf*lli%dlat/refinement_j ! southpole ? if ( abs(lat(jf) - (-pi/2)) < 1.0e-4 ) then ! fill with average of next row aver_to_prev = .true. cycle end if ! northpole ? if ( abs(lat(jf) - (pi/2)) < 1.0e-4 ) then ! fill with average of previous row: do l = 1, lm llf(:,jf,l) = sum(llf(:,jf-1,l)) / size(llf,1) end do exit end if ! evaluate Legendre functions: call sh_Pnm( Pnm, T, lat(jf), status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! evaluate rows: call Set( sh, T, H ) call Eval_Lons( exp_hh, sh, Pnm, nlon_fine_360, lli%blon(0), nlon_fine+1, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if exp_hh = exp(exp_hh) do l = 1, lm call Set( sh, T, F(l) ) call Eval_Lons( ff, sh, Pnm, nlon_fine_360, lli%blon(0), nlon_fine+1, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! combine fields: select case ( key ) ! ! [ int exp(H) dlat ] / deltalat ! case ( 'exp(H),aver' ) llf(:,jf,l) = exp_hh / lli%dlat ! ! int F (da+db*exp(H)) / cos(lat) dlat ! case ( '(da+exp*db)/cos' ) llf(:,jf,l) = ff * ( da(l) + exp_hh*db(l) ) / cos(lat(jf)) ! ! error ... ! case default print *, 'IntLat_sh_ll - unknown key "'//trim(key)//'"' stop end select end do ! loop over levels ! copy to southpole .. if ( (jf==1) .and. aver_to_prev ) then do l = 1, lm llf(:,0,l) = sum(llf(:,1,l)) / size(llf,1) end do aver_to_prev = .false. end if end do ! loop over rows in fine grid ! integral in lat direction: do l = 1, lm do i = 0, lli%im ilast = 1 llf_row = llf(i,:,l) call IntervalQuad_Lin( lat, llf_row, lli%blat(j-1), lli%blat(j), ll_u(i,j,l), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do end do end do ! loop over rows ! free memory call Done( sh ) deallocate( Pnm ) deallocate( ff ) deallocate( exp_hh ) deallocate( llf ) deallocate( llf_row ) deallocate( lat ) ! ok status = 0 end subroutine IntLat_sh_ll ! *** ! *** subroutine IntLat_shi_ll( key,shi_in,nlev, F, H, da, db, lli, ll_u, status ) use grid_tools, only : deg2rad, pi use Num, only : IntervalQuad_Lin use Grid_Type_sh, only : TshGridInfo, TshGrid, Init, Done, Set, Check use Grid_Type_sh, only : sh_Pnm, Eval_Lons, SpN use Grid_Type_ll, only : TllGridInfo, Check ! --- in/out ----------------------------------------- character(len=*), intent(in) :: key type(TshGridInfo), intent(in) :: shi_in integer, intent(in) :: nlev complex, intent(in) :: F(shi_in%np,nlev) complex, intent(in) :: H(shi_in%np) real, intent(in) :: da(nlev) real, intent(in) :: db(nlev) type(TllGridInfo), intent(in) :: lli real, intent(inout) :: ll_u(0:lli%im,lli%jm,nlev) integer, intent(out) :: status ! --- const -------------------------------------- character(len=*), parameter :: rname = mname//'/IntLat_shi_ll' ! --- local ------------------------------------------ integer :: i, j, jf integer :: l, lm integer :: ilast type(TshGridInfo) :: shi integer :: T ! type(TshGrid) :: sh real, allocatable :: Pnm(:) integer :: refinement_j integer :: nlon_fine_360, nlon_fine real, allocatable :: ff(:) real, allocatable :: exp_hh(:) real, allocatable :: llf(:,:,:), lat(:), llf_row(:) logical :: aver_to_prev ! --- begin ------------------------------------------ T = ShTruncation( shi_in%T, lli%dlon_deg, lli%dlat_deg ) ! T = ShTruncation( T, lli%dlon_deg, lli%dlat_deg ) call Init( shi, T,status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! allocate space for legendre coeff: allocate( Pnm(shi%np) ) ! determine refinement based on spectral resolution ! and length of required integration interval: refinement_j = ShRefinement( shi%T, lli%dlat_deg ) ! number of lons in fine grid on complete circle: nlon_fine_360 = 360.0/lli%dlon_deg nlon_fine = lli%im ! store evaluation of spectral field: allocate( ff(0:nlon_fine) ) allocate( exp_hh(0:nlon_fine) ) ! ll in lon, fine in lat allocate( llf(0:lli%im,0:refinement_j,nlev) ) allocate( llf_row(0:refinement_j) ) allocate( lat(0:refinement_j) ) ! loop over latitudes in ll grid: aver_to_prev = .false. do j = 1, lli%jm ! loop over latitudes in fine grid: do jf = 0, refinement_j ! latitude in fine grid: lat(jf) = lli%blat(j-1) + jf*lli%dlat/refinement_j ! southpole ? if ( abs(lat(jf) - (-pi/2)) < 1.0e-4 ) then ! fill with average of next row aver_to_prev = .true. cycle end if ! northpole ? if ( abs(lat(jf) - (pi/2)) < 1.0e-4 ) then ! fill with average of previous row: do l = 1, nlev llf(:,jf,l) = sum(llf(:,jf-1,l)) / size(llf,1) end do exit end if ! evaluate Legendre functions: call sh_Pnm( Pnm, shi%T, lat(jf), status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! evaluate rows: call Eval_Lons( shi_in, H, shi, Pnm, nlon_fine_360, lli%blon(0), nlon_fine+1,exp_hh, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if exp_hh = exp(exp_hh) do l = 1, nlev call Eval_Lons( shi_in, (/F(:,l)/), shi, Pnm, nlon_fine_360, lli%blon(0), nlon_fine+1,ff, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! combine fields: select case ( key ) ! ! [ int exp(H) dlat ] / deltalat ! case ( 'exp(H),aver' ) llf(:,jf,l) = exp_hh / lli%dlat ! ! int F (da+db*exp(H)) / cos(lat) dlat ! case ( '(da+exp*db)/cos' ) llf(:,jf,l) = ff * ( da(l) + exp_hh*db(l) ) / cos(lat(jf)) ! ! error ... ! case default print *, 'IntLat_sh_ll - unknown key "'//trim(key)//'"' stop end select end do ! loop over levels ! copy to southpole .. if ( (jf==1) .and. aver_to_prev ) then do l = 1, nlev llf(:,0,l) = sum(llf(:,1,l)) / size(llf,1) end do aver_to_prev = .false. end if end do ! loop over rows in fine grid ! integral in lat direction: do l = 1, nlev do i = 0, lli%im ilast = 1 llf_row = llf(i,:,l) call IntervalQuad_Lin( lat, llf_row, lli%blat(j-1), lli%blat(j), ll_u(i,j,l), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do end do end do ! loop over rows ! free memory call Done( shi ) deallocate( Pnm ) deallocate( ff ) deallocate( exp_hh ) deallocate( llf ) deallocate( llf_row ) deallocate( lat ) ! ok status = 0 end subroutine IntLat_shi_ll ! *** subroutine IntLon_sh_ll( key, F, H, da, db, lli, ll_v, status ) use grid_tools, only : deg2rad, pi use Num, only : IntervalQuad_Lin use Grid_Type_sh, only : TshGrid, Init, Done, SpN, Set, Check use Grid_Type_sh, only : sh_Pnm, Eval_Lons use Grid_Type_ll, only : TllGridInfo, Check ! --- in/out ----------------------------------------- character(len=*), intent(in) :: key type(TshGrid), intent(in) :: F(:) type(TshGrid), intent(in) :: H real, intent(in) :: da(size(F)) real, intent(in) :: db(size(F)) type(TllGridInfo), intent(in) :: lli real, intent(inout) :: ll_v(lli%im,0:lli%jm,size(F)) integer, intent(out) :: status ! --- const -------------------------------------- character(len=*), parameter :: rname = mname//'/IntLon_sh_ll' ! --- local ------------------------------------------ integer :: i, j integer :: l, lm integer :: ilast integer :: T type(TshGrid) :: sh real, allocatable :: Pnm(:) integer :: refinement_i integer :: nlon_fine_360, nlon_fine real, allocatable :: ff(:) real, allocatable :: exp_hh(:) real, allocatable :: lons_fine(:), row_fine(:) logical :: pole ! --- begin ------------------------------------------ ! number of levels: lm = size(F) ! check arguments: call Check( H ) T = H%T do l = 1, lm call Check( F(l), T ) end do ! use truncated harmonics ? T = ShTruncation( T, lli%dlon_deg, lli%dlat_deg ) call Init( sh, T ) ! allocate space for legendre coeff: allocate( Pnm(SpN(T)) ) ! determine refinement based on spectral resolution ! and length of required integration interval: refinement_i = ShRefinement( T, lli%dlon_deg ) ! number of lons in fine grid on complete circle: nlon_fine_360 = 360.0/lli%dlon_deg * refinement_i nlon_fine = lli%im * refinement_i ! store evaluation of spectral field: allocate( ff(0:nlon_fine) ) allocate( exp_hh(0:nlon_fine) ) ! lons on arc westb+[0,..) allocate( row_fine(0:nlon_fine) ) allocate( lons_fine(0:nlon_fine) ) do i = 0, nlon_fine lons_fine(i) = i*2*pi/nlon_fine_360 end do lons_fine = lli%blon(0) + lons_fine ! loop over boundary latitudes in ll grid: do j = 0, lli%jm ! pole ? then int f(x) dx = f pole = ( abs(lli%blat(j) - (-pi/2)) < 1.0e-4 ) .or. & ( abs(lli%blat(j) - ( pi/2)) < 1.0e-4 ) ! evaluate Legendre functions: call sh_Pnm( Pnm, T, lli%blat(j), status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! evaluate rows: call Set( sh, T, H ) call Eval_Lons( exp_hh, sh, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if exp_hh = exp(exp_hh) ! loop over levels do l = 1, lm call Set( sh, T, F(l) ) call Eval_Lons( ff, sh, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! combine fields: select case ( key ) ! ! [ int exp(H) dlon ] / deltalon ! case ( 'exp(H),aver' ) if ( pole ) then row_fine = exp_hh else row_fine = exp_hh / lli%dlon end if ! ! int F (da+db*exp(H)) dlon ! case ( '(da+exp*db)' ) if ( pole ) then row_fine = 0.0 else row_fine = ff * ( da(l) + exp_hh*db(l) ) end if ! ! error ... ! case default print *, 'IntLon_sh_ll - unknown key "'//trim(key)//'"' stop end select ! special treatment of poles: if ( pole ) then ! same value at al 'longitudes', thus for example the average ... ll_v(:,j,l) = sum(row_fine)/size(row_fine) else ! integral in lon direction assuming linear interpolation, ! result in [sh] rad ilast = 1 do i = 1, lli%im call IntervalQuad_Lin( lons_fine, row_fine, lli%blon(i-1), lli%blon(i), ll_v(i,j,l), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do end if end do ! loop over levels end do ! loop over rows ! free memory call Done( sh ) deallocate( Pnm ) deallocate( ff ) deallocate( exp_hh ) deallocate( lons_fine ) deallocate( row_fine ) ! ok status = 0 end subroutine IntLon_sh_ll ! *** subroutine IntLon_shi_ll( key, shi_in,nlev, F, H, da, db, lli, ll_v, status ) use grid_tools, only : deg2rad, pi use Num, only : IntervalQuad_Lin use Grid_Type_sh, only : TshGridInfo, TshGrid, Init, Done, SpN, Set, Check use Grid_Type_sh, only : sh_Pnm, Eval_Lons use Grid_Type_ll, only : TllGridInfo, Check ! --- in/out ----------------------------------------- character(len=*), intent(in) :: key type(TshGridInfo), intent(in) :: shi_in integer, intent(in) :: nlev complex, intent(in) :: F(shi_in%np,nlev) complex, intent(in) :: H(shi_in%np) real, intent(in) :: da(nlev) real, intent(in) :: db(nlev) type(TllGridInfo), intent(in) :: lli real, intent(inout) :: ll_v(lli%im,0:lli%jm,nlev) integer, intent(out) :: status ! --- const -------------------------------------- character(len=*), parameter :: rname = mname//'/IntLon_shi_ll' ! --- local ------------------------------------------ integer :: i, j integer :: l integer :: ilast type(TshGridInfo) :: shi integer :: T ! type(TshGrid) :: sh real, allocatable :: Pnm(:) integer :: refinement_i integer :: nlon_fine_360, nlon_fine real, allocatable :: ff(:) real, allocatable :: exp_hh(:) real, allocatable :: lons_fine(:), row_fine(:) logical :: pole ! --- begin ------------------------------------------ ! use truncated harmonics ? T = ShTruncation( shi_in%T, lli%dlon_deg, lli%dlat_deg ) ! T = ShTruncation( T, lli%dlon_deg, lli%dlat_deg ) call Init( shi, T ,status) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! allocate space for legendre coeff: allocate( Pnm(shi%np) ) ! determine refinement based on spectral resolution ! and length of required integration interval: refinement_i = ShRefinement( shi%T, lli%dlon_deg ) ! number of lons in fine grid on complete circle: nlon_fine_360 = 360.0/lli%dlon_deg * refinement_i nlon_fine = lli%im * refinement_i ! store evaluation of spectral field: allocate( ff(0:nlon_fine) ) allocate( exp_hh(0:nlon_fine) ) ! lons on arc westb+[0,..) allocate( row_fine(0:nlon_fine) ) allocate( lons_fine(0:nlon_fine) ) do i = 0, nlon_fine lons_fine(i) = i*2*pi/nlon_fine_360 end do lons_fine = lli%blon(0) + lons_fine ! loop over boundary latitudes in ll grid: do j = 0, lli%jm ! pole ? then int f(x) dx = f pole = ( abs(lli%blat(j) - (-pi/2)) < 1.0e-4 ) .or. & ( abs(lli%blat(j) - ( pi/2)) < 1.0e-4 ) ! evaluate Legendre functions: call sh_Pnm( Pnm, shi%T, lli%blat(j), status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! evaluate rows: ! call Set( sh, T, H ) ! call Eval_Lons( exp_hh, sh, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, status ) call Eval_Lons( shi_in,H,shi,Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1,exp_hh, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if exp_hh = exp(exp_hh) ! loop over levels do l = 1, nlev ! call Set( sh, T, F(l) ) ! call Eval_Lons( ff, sh, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, status ) call Eval_Lons( shi_in,F(:,l), shi, Pnm, nlon_fine_360, lons_fine(0), nlon_fine+1, ff,status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if ! combine fields: select case ( key ) ! ! [ int exp(H) dlon ] / deltalon ! case ( 'exp(H),aver' ) if ( pole ) then row_fine = exp_hh else row_fine = exp_hh / lli%dlon end if ! ! int F (da+db*exp(H)) dlon ! case ( '(da+exp*db)' ) if ( pole ) then row_fine = 0.0 else row_fine = ff * ( da(l) + exp_hh*db(l) ) end if ! ! error ... ! case default print *, 'IntLon_sh_ll - unknown key "'//trim(key)//'"' stop end select ! special treatment of poles: if ( pole ) then ! same value at al 'longitudes', thus for example the average ... ll_v(:,j,l) = sum(row_fine)/size(row_fine) else ! integral in lon direction assuming linear interpolation, ! result in [sh] rad ilast = 1 do i = 1, lli%im call IntervalQuad_Lin( lons_fine, row_fine, lli%blon(i-1), lli%blon(i), ll_v(i,j,l), ilast, status ) if (status/=0) then; write (*,'("ERROR in ",a)') rname; status=1; return; end if end do end if end do ! loop over levels end do ! loop over rows ! free memory call Done( shi ) deallocate( Pnm ) deallocate( ff ) deallocate( exp_hh ) deallocate( lons_fine ) deallocate( row_fine ) ! ok status = 0 end subroutine IntLon_shi_ll end module Grid_Interpol