From 2dd042a5c6c19ce41fc0884782e9ada733495a03 Mon Sep 17 00:00:00 2001 From: Niki Zadeh Date: Mon, 15 Oct 2018 17:30:54 -0400 Subject: [PATCH] Diag decimation prototype, first attemp at general algorithm - According to Alistair, the decimation method could be solely deduced from the axes%x_cell_method, axes%y_cell_method and probably the area_cell_method at the time of send_data - This is the summary of the algoritm f(Id,Jd) = \sum_{i,j} f(Id+i,Jd+j) * weight(Id+i,Jd+j) / [ \sum_{i,j} weight(Id+i,Jd+j)] Id,Jd are the decimated (coarse grid) indices run over the coarsened compute grid, i and j run from 0 to dl-1 (dl being the decimation level) if and jf weight(if,jf) run over the original fine computre grid x_cell_method y_cell_method area_cell_method weight(if,jf) example --------------------------------------------------------------------- ------------- mean mean mean A(if,jf)*h(if,jf) theta point mean mean dy(if,jf)*h(if,jf) u mean point mean dx(if,jf)*h(if,jf) v mean mean sum A(if,jf) h*theta sum sum sum 1 volcello point sum sum 1 umo --- src/framework/MOM_diag_mediator.F90 | 318 ++++++++++++++++++---------- 1 file changed, 201 insertions(+), 117 deletions(-) diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index 02033958f0..bed00b2355 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -3434,8 +3434,6 @@ subroutine decimate_diag_indices_get(f1,f2, dl, diag_cs,isv,iev,jsv,jev) endif end subroutine decimate_diag_indices_get - - subroutine decimate_diag_field_3d(locfield, locfield_decim, dl, diag_cs, diag,isv,iev,jsv,jev, mask) real, dimension(:,:,:), pointer :: locfield @@ -3446,20 +3444,46 @@ subroutine decimate_diag_field_3d(locfield, locfield_decim, dl, diag_cs, diag,is integer, intent(out):: isv,iev,jsv,jev real, optional,target, intent(in) :: mask(:,:,:) !< If present, use this real array as the data mask. !locals - real, dimension(:,:,:), pointer :: locmask => NULL() - integer :: isl,iel,jsl,jel + real, dimension(:,:,:), pointer :: locmask + integer :: isl,iel,jsl,jel, xy_method + + locmask => NULL() + !Get the correct indices corresponding to input field + !Shape of the input diag field isl=1; iel=size(locfield,1)/dl - jsl=1; jel=size(locfield,2)/dl + jsl=1; jel=size(locfield,2)/dl + !Get the shape of the decimated field isv,iev,jsv,jev call decimate_diag_indices_get(iel,jel, dl, diag_cs,isv,iev,jsv,jev) + !Set the non-decimated mask, it must be associated and initialized if (present(mask)) then locmask => mask - call decimate_field(locfield, locfield_decim, dl, method='pave', mask=locmask, area=diag_cs%G%areaT) elseif (associated(diag%axes%mask3d)) then - call decimate_field(locfield, locfield_decim, dl, method='pave', mask=diag%axes%mask3d, area=diag_cs%G%areaT) + locmask => diag%axes%mask3d else call MOM_error(FATAL, "decimate_diag_field_3d: Cannot decimate without a mask!!! ") endif + !Determine the decimation method + !Make a two digit integer case id = x_case*10+y_case where x_case and/or y_case is 1 for sum and 2 for mean + xy_method = 0 !subsample and pick the SW corner value + if (trim(diag%axes%y_cell_method)=='sum') then + xy_method = xy_method + 1 + elseif (trim(diag%axes%y_cell_method)=='mean') then + xy_method = xy_method + 2 + endif + if (trim(diag%axes%x_cell_method)=='sum') then + xy_method = xy_method + 10 + elseif (trim(diag%axes%x_cell_method)=='mean') then + xy_method = xy_method + 20 + endif +! if (trim(diag%axes%area_cell_method)=='sum') then +! xy_method = xy_method + 100 +! elseif (trim(diag%axes%area_cell_method)=='mean') then +! xy_method = xy_method + 200 +! endif + + call decimate_field(locfield, locfield_decim, dl, xy_method, locmask, diag_cs) + end subroutine decimate_diag_field_3d subroutine decimate_diag_field_2d(locfield, locfield_decim, dl, diag_cs, diag,isv,iev,jsv,jev, mask) @@ -3471,157 +3495,217 @@ subroutine decimate_diag_field_2d(locfield, locfield_decim, dl, diag_cs, diag,is integer, intent(out):: isv,iev,jsv,jev real, optional,target, intent(in) :: mask(:,:) !< If present, use this real array as the data mask. !locals - real, dimension(:,:), pointer :: locmask => NULL() - integer :: isl,iel,jsl,jel + real, dimension(:,:), pointer :: locmask + integer :: isl,iel,jsl,jel, xy_method + + locmask => NULL() + !Get the correct indices corresponding to input field + !Shape of the input diag field isl=1; iel=size(locfield,1)/dl - jsl=1; jel=size(locfield,2)/dl + jsl=1; jel=size(locfield,2)/dl + !Get the shape of the decimated field isv,iev,jsv,jev call decimate_diag_indices_get(iel,jel, dl, diag_cs,isv,iev,jsv,jev) + !Set the non-decimated mask, it must be associated and initialized if (present(mask)) then locmask => mask - call decimate_field(locfield, locfield_decim, dl, method='pave', mask=locmask) elseif (associated(diag%axes%mask2d)) then - call decimate_field(locfield, locfield_decim, dl, method='pave', mask=diag%axes%mask2d) + locmask => diag%axes%mask2d else call MOM_error(FATAL, "decimate_diag_field_2d: Cannot decimate without a mask!!! ") endif - + !Determine the decimation method + !Make a two digit integer case id = x_case*10+y_case where x_case and/or y_case is 1 for sum and 2 for mean + xy_method = 0 !subsample and pick the SW corner value + if (trim(diag%axes%y_cell_method)=='sum') then + xy_method = xy_method + 1 + elseif (trim(diag%axes%y_cell_method)=='mean') then + xy_method = xy_method + 2 + endif + if (trim(diag%axes%x_cell_method)=='sum') then + xy_method = xy_method + 10 + elseif (trim(diag%axes%x_cell_method)=='mean') then + xy_method = xy_method + 20 + endif +! if (trim(diag%axes%area_cell_method)=='sum') then +! xy_method = xy_method + 100 +! elseif (trim(diag%axes%area_cell_method)=='mean') then +! xy_method = xy_method + 200 +! endif + + call decimate_field(locfield, locfield_decim, dl, xy_method, locmask, diag_cs) + end subroutine decimate_diag_field_2d -subroutine decimate_field_3d(field_in, field_out, level, method, mask, area) +subroutine decimate_field_3d(field_in, field_out, dl, method, mask, diag_cs) real, dimension(:,:,:) , pointer :: field_in real, dimension(:,:,:) , allocatable :: field_out - integer , intent(in) :: level - character(len=4), optional, intent(in) :: method !< sampling method, one of samp,pave,aave,vave - real, dimension(:,:,:), optional , pointer :: mask - real, dimension(:,:), optional , intent(in) :: area - !locals + integer , intent(in) :: dl + integer, optional, intent(in) :: method !< sampling method, one of 00,01,02,10,20,11,22 + real, dimension(:,:,:), pointer :: mask + type(diag_ctrl), intent(in) :: diag_CS !< Structure used to regulate diagnostic output + !locals + character(len=240) :: mesg integer :: i,j,ii,jj,is,js,i0,j0 integer :: isl,iel,jsl,jel integer :: k,ks,ke - real :: ave,tot_non_zero,a1 - character(len=4) :: samplemethod - samplemethod = 'samp' - if(present(method)) samplemethod = method - + real :: ave,total_weight,weight + real :: epsilon = 1.0e-20 !Always start from the first element is=1 js=1 ks = lbound(field_in,3) ; ke = ubound(field_in,3) - isl=1; iel=size(field_in,1)/level - jsl=1; jel=size(field_in,2)/level - allocate(field_out(isl:iel,jsl:jel,ks:ke)) + isl=1; iel=size(field_in,1)/dl + jsl=1; jel=size(field_in,2)/dl + allocate(field_out(isl:iel,jsl:jel,ks:ke)) - - select case (samplemethod) - case ('samp') !subsample the SW corner cell + if(method .eq. 0) then !point average the fields in dl^2 cells do k= ks,ke ; do j=jsl,jel ; do i=isl,iel - i0 = is+level*(i-isl) - j0 = js+level*(j-jsl) - field_out(i,j,k) = field_in(i0,j0,k) + i0 = is+dl*(i-isl) + j0 = js+dl*(j-jsl) + ave = 0.0 + total_weight = 0.0 + do ii=i0,i0+dl-1 ; do jj=j0,j0+dl-1 + total_weight = total_weight + mask(ii,jj,k) + ave=ave+field_in(ii,jj,k)*mask(ii,jj,k) + enddo; enddo + field_out(i,j,k) = ave/max(1.0,total_weight) !Avoid zero mask at all aggregating cells where ave=0.0 enddo; enddo; enddo - case ('pave') !point average of the cells - if(present(mask)) then - do k= ks,ke ; do j=jsl,jel ; do i=isl,iel - i0 = is+level*(i-isl) - j0 = js+level*(j-jsl) - ave = 0.0 - tot_non_zero = 0.0 - do ii=i0,i0+level-1 ; do jj=j0,j0+level-1 - tot_non_zero = tot_non_zero + mask(ii,jj,k) - ave=ave+field_in(ii,jj,k)*mask(ii,jj,k) - enddo; enddo - field_out(i,j,k) = ave/max(1.0,tot_non_zero) !Avoid zero mask at all aggregating cells where ave=0.0 - enddo; enddo; enddo - else - do k= ks,ke ; do j=jsl,jel ; do i=isl,iel - i0 = is+level*(i-isl) - j0 = js+level*(j-jsl) - ave = 0.0 - tot_non_zero = 0.0 - do ii=i0,i0+level-1 ; do jj=j0,j0+level-1 - tot_non_zero = tot_non_zero + 1 - ave=ave+field_in(ii,jj,k) - enddo; enddo - field_out(i,j,k) = ave/tot_non_zero - enddo; enddo; enddo - endif - case ('aave') !area average of the cells + elseif(method .eq. 1) then !point in x, sum in y + do k= ks,ke ; do j=jsl,jel ; do i=isl,iel + i0 = is+dl*(i-isl) + j0 = js+dl*(j-jsl) + ave = 0.0 + total_weight = 0.0 + ii=i0 + do jj=j0,j0+dl-1 + total_weight = total_weight + mask(ii,jj,k) + ave=ave+field_in(ii,jj,k)*mask(ii,jj,k) + enddo + field_out(i,j,k) = ave/(total_weight+epsilon) !Avoid zero mask at all aggregating cells where ave=0.0 + enddo; enddo; enddo + elseif(method .eq. 10) then !sum in x, point in y + do k= ks,ke ; do j=jsl,jel ; do i=isl,iel + i0 = is+dl*(i-isl) + j0 = js+dl*(j-jsl) + ave = 0.0 + total_weight = 0.0 + jj=j0 + do ii=i0,i0+dl-1 + total_weight = total_weight + mask(ii,jj,k) + ave=ave+field_in(ii,jj,k)*mask(ii,jj,k) + enddo + field_out(i,j,k) = ave/(total_weight+epsilon) !Avoid zero mask at all aggregating cells where ave=0.0 + enddo; enddo; enddo + elseif(method .eq. 2) then !point in x, normal area average in y do k= ks,ke ; do j=jsl,jel ; do i=isl,iel - field_out(i,j,k) = 0.0 - i0 = is+level*(i-isl) - j0 = js+level*(j-jsl) + i0 = is+dl*(i-isl) + j0 = js+dl*(j-jsl) ave = 0.0 - tot_non_zero = 0.0 - do ii=i0,i0+level-1 ; do jj=j0,j0+level-1 - a1 = area(ii,jj)*mask(ii,jj,k) - tot_non_zero = tot_non_zero + a1 - ave=ave + field_in(ii,jj,k) * a1 + total_weight = 0.0 + ii=i0 + do jj=j0,j0+dl-1 + weight =mask(ii,jj,k)*diag_cs%G%dyCu(ii,jj)*diag_cs%h(ii,jj,k) + total_weight = total_weight +weight + ave=ave+field_in(ii,jj,k)*weight + enddo + field_out(i,j,k) = ave/(total_weight+epsilon) !Avoid zero mask at all aggregating cells where ave=0.0 + enddo; enddo; enddo + elseif(method .eq. 20) then !normal area average in x, point in y + do k= ks,ke ; do j=jsl,jel ; do i=isl,iel + i0 = is+dl*(i-isl) + j0 = js+dl*(j-jsl) + ave = 0.0 + total_weight = 0.0 + jj=j0 + do ii=i0,i0+dl-1 + weight = mask(ii,jj,k)*diag_cs%G%dxCv(ii,jj)*diag_cs%h(ii,jj,k) + total_weight = total_weight + weight + ave=ave+field_in(ii,jj,k)*weight + enddo + field_out(i,j,k) = ave/(total_weight+epsilon) !Avoid zero mask at all aggregating cells where ave=0.0 + enddo; enddo; enddo + elseif(method .eq. 22) then !Volume average the fields in x and y + do k= ks,ke ; do j=jsl,jel ; do i=isl,iel + i0 = is+dl*(i-isl) + j0 = js+dl*(j-jsl) + ave = 0.0 + total_weight = 0.0 + do ii=i0,i0+dl-1 ; do jj=j0,j0+dl-1 + weight = mask(ii,jj,k)*diag_cs%G%areaT(ii,jj)*diag_cs%h(ii,jj,k) + total_weight = total_weight + weight + ave=ave+field_in(ii,jj,k)*weight enddo; enddo - if(tot_non_zero .gt. 0.0) field_out(i,j,k) = ave/tot_non_zero + field_out(i,j,k) = ave/(total_weight+epsilon) !Avoid zero mask at all aggregating cells where ave=0.0 enddo; enddo; enddo - case default - call MOM_error(FATAL, "decimate_field_3d: unknown sampling method "//trim(samplemethod)) - end select - + elseif(method .eq. 11) then !sum the fields in x and y + do k= ks,ke ; do j=jsl,jel ; do i=isl,iel + i0 = is+dl*(i-isl) + j0 = js+dl*(j-jsl) + ave = 0.0 + total_weight = 0.0 + do ii=i0,i0+dl-1 ; do jj=j0,j0+dl-1 + total_weight = total_weight + mask(ii,jj,k) + ave=ave+field_in(ii,jj,k)*mask(ii,jj,k) + enddo; enddo + field_out(i,j,k) = ave/(total_weight+epsilon) !Avoid zero mask at all aggregating cells where ave=0.0 + enddo; enddo; enddo + else + write (mesg,*) " unknown sampling method: ",method + call MOM_error(FATAL, "decimate_field_3d: "//trim(mesg)) + endif + end subroutine decimate_field_3d -subroutine decimate_field_2d(field_in, field_out, level, method, mask) +subroutine decimate_field_2d(field_in, field_out, dl, method, mask, diag_cs) real, dimension(:,:) , pointer :: field_in real, dimension(:,:) , allocatable :: field_out - integer , intent(in) :: level - character(len=4), optional, intent(in) :: method !< sampling method, one of samp,pave,aave,vave - real, dimension(:,:), optional , pointer :: mask + integer , intent(in) :: dl + integer, optional, intent(in) :: method !< sampling method, one of 00,01,02,10,20,11,22 + real, dimension(:,:), pointer :: mask + type(diag_ctrl), intent(in) :: diag_CS !< Structure used to regulate diagnostic output !locals + character(len=240) :: mesg integer :: i,j,ii,jj,is,js,i0,j0 integer :: isl,iel,jsl,jel - real :: ave,tot_non_zero - character(len=4) :: samplemethod - samplemethod = 'samp' - if(present(method)) samplemethod = method - + real :: ave,total_weight,weight + real :: epsilon = 1.0e-20 !Always start from the first element is=1 js=1 - isl=1; iel=size(field_in,1)/level - jsl=1; jel=size(field_in,2)/level + isl=1; iel=size(field_in,1)/dl + jsl=1; jel=size(field_in,2)/dl allocate(field_out(isl:iel,jsl:jel)) - select case (samplemethod) - case ('samp') !subsample the SW corner cell + if(method .eq. 0) then !point average the fields in dl^2 cells do j=jsl,jel ; do i=isl,iel - i0 = is+level*(i-isl) - j0 = js+level*(j-jsl) - field_out(i,j) = field_in(i0,j0) - enddo; enddo - case ('pave') !point average of the cells - if(present(mask)) then - do j=jsl,jel ; do i=isl,iel - i0 = is+level*(i-isl) - j0 = js+level*(j-jsl) - ave = 0.0 - tot_non_zero = 0.0 - do ii=i0,i0+level-1 ; do jj=j0,j0+level-1 - tot_non_zero = tot_non_zero + mask(ii,jj) - ave=ave+field_in(ii,jj)*mask(ii,jj) - enddo; enddo - field_out(i,j) = ave/max(1.0,tot_non_zero) !Avoid zero mask at all aggregating cells where ave=0.0 + i0 = is+dl*(i-isl) + j0 = js+dl*(j-jsl) + ave = 0.0 + total_weight = 0.0 + do ii=i0,i0+dl-1 ; do jj=j0,j0+dl-1 + total_weight = total_weight + mask(ii,jj) + ave=ave+field_in(ii,jj)*mask(ii,jj) enddo; enddo - else !Niki: How are we supposed to decimate/average without a mask? What if field_in is on land at one or more aggregating cells? - do j=jsl,jel ; do i=isl,iel - i0 = is+level*(i-isl) - j0 = js+level*(j-jsl) - ave = 0.0 - tot_non_zero = 0.0 - do ii=i0,i0+level-1 ; do jj=j0,j0+level-1 - tot_non_zero = tot_non_zero + 1 - ave=ave+field_in(ii,jj) - enddo; enddo - field_out(i,j) = ave/tot_non_zero + field_out(i,j) = ave/max(1.0,total_weight) !Avoid zero mask at all aggregating cells where ave=0.0 + enddo; enddo + elseif(method .eq. 22) then !Volume average the fields in x and y + do j=jsl,jel ; do i=isl,iel + i0 = is+dl*(i-isl) + j0 = js+dl*(j-jsl) + ave = 0.0 + total_weight = 0.0 + do ii=i0,i0+dl-1 ; do jj=j0,j0+dl-1 + weight = mask(ii,jj)*diag_cs%G%areaT(ii,jj)*diag_cs%h(ii,jj,1) + total_weight = total_weight + weight + ave=ave+field_in(ii,jj)*weight enddo; enddo - endif - case default - call MOM_error(FATAL, "decimate_field_2d: unknown sampling method "//trim(samplemethod)) - end select + field_out(i,j) = ave/(total_weight+epsilon) !Avoid zero mask at all aggregating cells where ave=0.0 + enddo; enddo + + else + write (mesg,*) " unknown sampling method: ",method + call MOM_error(FATAL, "decimate_field_2d: "//trim(mesg)) + endif end subroutine decimate_field_2d