CESM_Driver: map_atmocn_mct.F90 Source File

00001 module map_atmocn_mct
00002 
00003 !---------------------------------------------------------------------
00004 !
00005 ! Purpose:
00006 !
00007 ! Collect coupling routines for sequential coupling of OCN-ATM.
00008 !       
00009 ! Author: R. Jacob, M. Vertenstein
00010 !
00011 !---------------------------------------------------------------------
00012 
00013   use shr_kind_mod      ,only: R8 => SHR_KIND_R8, IN=>SHR_KIND_IN
00014   use shr_sys_mod
00015   use shr_const_mod
00016   use shr_mct_mod, only: shr_mct_sMatPInitnc
00017   use mct_mod
00018 
00019   use seq_comm_mct, only : logunit, loglevel
00020   use seq_cdata_mod
00021   use seq_flds_indices
00022   use seq_infodata_mod
00023   use m_die
00024 
00025   implicit none
00026   save
00027   private  ! except
00028 
00029 !--------------------------------------------------------------------------
00030 ! Public interfaces
00031 !--------------------------------------------------------------------------
00032 
00033   public :: map_ocn2atm_init_mct
00034   public :: map_atm2ocn_init_mct
00035   public :: map_ocn2atm_mct
00036   public :: map_atm2ocn_mct
00037 
00038 !--------------------------------------------------------------------------
00039 ! Public data
00040 !--------------------------------------------------------------------------
00041 
00042 !--------------------------------------------------------------------------
00043 ! Private data
00044 !--------------------------------------------------------------------------
00045 
00046   type(mct_rearr), private :: Re_ocn2atm
00047   type(mct_rearr), private :: Re_atm2ocn
00048   type(mct_sMatp), private :: sMatp_Fa2o
00049   type(mct_sMatp), private :: sMatp_Sa2o
00050   type(mct_sMatp), private :: sMatp_Fo2a
00051   type(mct_sMatp), private :: sMatp_So2a
00052 
00053 #ifdef CPP_VECTOR
00054   logical :: usevector = .true.
00055 #else
00056   logical :: usevector = .false.
00057 #endif
00058   
00059 #ifdef SYSUNICOS
00060   logical :: usealltoall = .true.
00061 #else
00062   logical :: usealltoall = .false.
00063 #endif
00064   logical, private :: samegrid_mapa2o
00065   
00066   integer         :: ni_a   ! number of longitudes on input grid
00067   integer         :: nj_a   ! number of latitudes  on input grid
00068   integer         :: ni_o   ! number of longitudes on output grid
00069   integer         :: nj_o   ! number of latitudes  on output grid
00070    
00071 !=======================================================================
00072 contains
00073 !=======================================================================
00074 
00075   subroutine map_atm2ocn_init_mct( cdata_a, cdata_o)
00076 
00077     !-----------------------------------------------------
00078     ! 
00079     ! Arguments
00080     !
00081     type(seq_cdata),intent(in) :: cdata_a
00082     type(seq_cdata),intent(in) :: cdata_o
00083     !
00084     ! Local Variables
00085     !
00086     integer                  :: ka, km            ! indices
00087     integer                  :: lsize             ! size of attribute vector
00088     type(seq_infodata_type), pointer :: infodata
00089     type(mct_gsMap), pointer :: gsMap_a
00090     type(mct_gsMap), pointer :: gsMap_o
00091     type(mct_ggrid), pointer :: dom_a
00092     type(mct_ggrid), pointer :: dom_o
00093     type(mct_aVect)          :: areasrc           ! ocn areas from mapping file
00094     type(mct_aVect)          :: areadst           ! atm areas from mapping file
00095     integer                  :: mpicom
00096     character(*),parameter :: subName = '(map_atm2ocn_init_mct) '
00097     !-----------------------------------------------------
00098 
00099     call seq_cdata_setptrs(cdata_a, gsMap=gsMap_a, dom=dom_a)
00100     call seq_cdata_setptrs(cdata_o, gsMap=gsMap_o, dom=dom_o)
00101     call seq_cdata_setptrs(cdata_a, mpicom=mpicom, infodata=infodata)
00102 
00103     call seq_infodata_GetData( infodata, samegrid_ao=samegrid_mapa2o)
00104 
00105     if (samegrid_mapa2o) then
00106 
00107        call mct_rearr_init(gsMap_a, gsMap_o, mpicom, Re_atm2ocn)
00108 
00109        !--- want to copy atm area into atm and ocn aream ---
00110 
00111        lsize = mct_gsMap_lsize(gsMap_a, mpicom)
00112        call mct_aVect_init( areasrc, rList="aream", lsize=lsize )
00113        
00114        lsize = mct_gsMap_lsize(gsMap_o, mpicom)
00115        call mct_aVect_init( areadst, rList="aream", lsize=lsize )
00116        
00117        ka = mct_aVect_indexRa(dom_a%data, "area" )
00118        km = mct_aVect_indexRA(areasrc   , "aream")
00119        areasrc%rAttr(km,:) = dom_a%data%rAttr(ka,:)
00120 
00121        call mct_rearr_rearrange(areasrc, areadst, Re_atm2ocn, VECTOR=usevector, &
00122           ALLTOALL=usealltoall)
00123 
00124        ka = mct_aVect_indexRA(areasrc   ,"aream")
00125        km = mct_aVect_indexRA(dom_a%data,"aream")
00126        dom_a%data%rAttr(km,:) = areasrc%rAttr(ka,:)
00127 
00128        ka = mct_aVect_indexRA(areadst   ,"aream")
00129        km = mct_aVect_indexRA(dom_o%data,"aream")
00130        dom_o%data%rAttr(km,:) = areadst%rAttr(ka,:)
00131 
00132        call mct_aVect_clean(areasrc)
00133        call mct_aVect_clean(areadst)      
00134 
00135     else      
00136 
00137        call shr_mct_sMatPInitnc(sMatp_Fa2o,gsMap_a,gsMap_o,"seq_maps.rc", &
00138           "atm2ocnFmapname:","atm2ocnFmaptype:",mpicom,&
00139           ni_i=ni_a,nj_i=nj_a,ni_o=ni_o,nj_o=nj_o)
00140 
00141        call shr_mct_sMatPInitnc(sMatp_Sa2o,gsMap_a,gsMap_o,"seq_maps.rc", &
00142           "atm2ocnSmapname:","atm2ocnSmaptype:",mpicom)
00143 
00144     endif
00145 
00146   end subroutine map_atm2ocn_init_mct
00147 
00148 !=======================================================================
00149 
00150   subroutine map_ocn2atm_init_mct( cdata_o, cdata_a)
00151 
00152     !-----------------------------------------------------
00153     !
00154     ! Arguments
00155     !
00156     type(seq_cdata),intent(in) :: cdata_o
00157     type(seq_cdata),intent(in) :: cdata_a
00158     !
00159     ! Local Variables
00160     !
00161     type(seq_infodata_type), pointer :: infodata
00162     integer                  :: ka, km            ! indices
00163     type(mct_gsMap), pointer :: gsMap_o           ! ocn gsMap
00164     type(mct_gsMap), pointer :: gsMap_a           ! atm gsMap
00165     type(mct_ggrid), pointer :: dom_o             ! ocn domain
00166     type(mct_ggrid), pointer :: dom_a             ! atm domain
00167     integer                  :: mpicom            ! communicator spanning atm and ocn
00168     integer                  :: lsize             ! size of attribute vector
00169     type(mct_aVect)          :: areasrc           ! ocn areas from mapping file
00170     type(mct_aVect)          :: areadst           ! atm areas from mapping file
00171     character(*),parameter :: subName = '(map_ocn2atm_init_mct) '
00172     !-----------------------------------------------------
00173 
00174     call seq_cdata_setptrs(cdata_o, gsMap=gsMap_o, dom=dom_o)
00175     call seq_cdata_setptrs(cdata_a, gsMap=gsMap_a, dom=dom_a)
00176     call seq_cdata_setptrs(cdata_o, mpicom=mpicom, infodata=infodata)
00177 
00178     call seq_infodata_GetData( infodata, samegrid_ao=samegrid_mapa2o)
00179 
00180     ! Initialize ocn->atm mapping or rearranging
00181 
00182     if (samegrid_mapa2o) then
00183 
00184        call mct_rearr_init(gsMap_o, gsMap_a, mpicom, Re_ocn2atm)
00185 
00186     else
00187 
00188        lsize = mct_gsMap_lsize(gsMap_o, mpicom)
00189        call mct_aVect_init( areasrc, rList="aream", lsize=lsize )
00190        
00191        lsize = mct_gsMap_lsize(gsMap_a, mpicom)
00192        call mct_aVect_init( areadst, rList="aream", lsize=lsize )
00193        
00194        call shr_mct_sMatPInitnc(sMatp_Fo2a, gsMap_o, gsMap_a, "seq_maps.rc", &
00195             "ocn2atmFmapname:", "ocn2atmFmaptype:", mpicom, &
00196             areasrc=areasrc, areadst=areadst)
00197 
00198        call shr_mct_sMatPInitnc(sMatp_So2a, gsMap_o, gsMap_a, "seq_maps.rc", &
00199             "ocn2atmSmapname:", "ocn2atmSmaptype:", mpicom)
00200 
00201        !--- copy aream from mapping files
00202 
00203        km = mct_aVect_indexRA(dom_o%data,"aream")
00204        ka = mct_aVect_indexRA(areasrc   ,"aream")
00205        dom_o%data%rAttr(km,:) = areasrc%rAttr(ka,:)
00206 
00207        km = mct_aVect_indexRA(dom_a%data,"aream")
00208        ka = mct_aVect_indexRA(areadst   ,"aream")
00209        dom_a%data%rAttr(km,:) = areadst%rAttr(ka,:)
00210 
00211        call mct_aVect_clean(areasrc)
00212        call mct_aVect_clean(areadst)      
00213 
00214     endif
00215 
00216  end subroutine map_ocn2atm_init_mct
00217 
00218 !=======================================================================
00219 
00220  subroutine map_atm2ocn_mct( cdata_a, av_a, cdata_o, av_o, &
00221                              fluxlist, statelist)
00222 
00223     !-----------------------------------------------------
00224     !
00225     ! Arguments
00226     !
00227     type(seq_cdata) ,intent(in)          :: cdata_a
00228     type(mct_aVect) ,intent(in)          :: av_a
00229     type(seq_cdata) ,intent(in)          :: cdata_o
00230     type(mct_aVect) ,intent(out)         :: av_o
00231     character(len=*),intent(in),optional :: fluxlist
00232     character(len=*),intent(in),optional :: statelist
00233     !
00234     ! Local Variables
00235     !
00236     type(seq_infodata_type), pointer :: infodata
00237     type(mct_ggrid), pointer :: dom_o
00238     type(mct_ggrid), pointer :: dom_a
00239     type(mct_gsMap), pointer :: gsmap_a
00240     integer                  :: mpicom
00241     integer                  :: lsize
00242     integer                  :: ku,kv
00243     type(mct_aVect)          :: av_o_f     ! temporary flux attribute vector
00244     type(mct_aVect)          :: av_o_s     ! temporary state attribute vector
00245     logical                  :: do_npfix          ! npfix on or off
00246     character(*),parameter :: subName = '(map_atm2ocn_mct) '
00247     !-----------------------------------------------------
00248 
00249     call seq_cdata_setptrs(cdata_o, infodata=infodata)
00250     call seq_infodata_GetData( infodata, npfix=do_npfix)
00251 
00252     if (samegrid_mapa2o) then
00253 
00254        if (present(fluxlist) .or. present(statelist)) then
00255           if (present(fluxlist)) then
00256              call mct_rearr_rearrange_fldlist(av_a, av_o, Re_atm2ocn, VECTOR=usevector, &
00257                   ALLTOALL=usealltoall, fldlist=fluxlist)
00258           endif
00259           if (present(statelist)) then
00260              call mct_rearr_rearrange_fldlist(av_a, av_o, Re_atm2ocn, VECTOR=usevector, &
00261                   ALLTOALL=usealltoall, fldlist=statelist)
00262           endif
00263        else
00264           call mct_rearr_rearrange(av_a, av_o, Re_atm2ocn, VECTOR=usevector, ALLTOALL=usealltoall)
00265        endif
00266 
00267     else
00268        
00269       if (present(fluxlist) .or. present(statelist)) then
00270          if (present(fluxlist)) then
00271             lsize = mct_aVect_lsize(av_o)
00272             call mct_aVect_init (av_o_f, rlist=fluxlist , lsize=lsize)
00273             call mct_sMat_avMult(av_a, sMatp_Fa2o, av_o_f, VECTOR=usevector, rList=fluxlist)
00274             call mct_aVect_copy (aVin=av_o_f, aVout=av_o, vector=usevector)
00275             call mct_aVect_clean(av_o_f)
00276          end if
00277          if (present(statelist)) then
00278             lsize = mct_aVect_lsize(av_o)
00279             call mct_aVect_init (av_o_s, rlist=statelist, lsize=lsize)
00280             call mct_sMat_avMult(av_a, sMatp_Sa2o, av_o_s, VECTOR=usevector, rList=statelist)
00281             call mct_aVect_copy (aVin=av_o_s, aVout=av_o, vector=usevector)
00282             call mct_aVect_clean(av_o_s)
00283          end if
00284       else
00285          !--- default is flux mapping
00286          call mct_sMat_avMult(av_a, sMatp_Fa2o, av_o, VECTOR=usevector)
00287       endif
00288          
00289       ! Correct a->o vector mapping near NP if appropriate
00290       
00291       if (do_npfix) then
00292       if (present(statelist)) then 
00293          ku = mct_aVect_indexRA(av_a, 'Sa_u', perrwith='quiet')
00294          kv = mct_aVect_indexRA(av_a, 'Sa_v', perrwith='quiet')
00295          if (ku /= 0 .and. kv /= 0) then
00296             call seq_cdata_setptrs(cdata_o, dom=dom_o)
00297             call seq_cdata_setptrs(cdata_a, dom=dom_a, gsmap=gsmap_a)
00298             call seq_cdata_setptrs(cdata_a, mpicom=mpicom)
00299             call map_npfixNew4R(av_a, av_o, 'Sa_u', 'Sa_v', gsmap_a, dom_a, dom_o, ni_a, nj_a, mpicom)
00300          end if
00301       end if
00302       end if
00303 
00304     endif
00305 
00306   end subroutine map_atm2ocn_mct
00307 
00308 !=======================================================================
00309 
00310   subroutine map_ocn2atm_mct( cdata_o, av_o, cdata_a, av_a, &
00311                               fractions_o, fractions_a, &
00312                               fluxlist, statelist )
00313 
00314     !-----------------------------------------------------
00315     !
00316     ! Arguments
00317     !
00318     type(seq_cdata) , intent(in)           :: cdata_o
00319     type(mct_aVect) , intent(in)           :: av_o
00320     type(seq_cdata) , intent(in)           :: cdata_a
00321     type(mct_aVect) , intent(out)          :: av_a
00322     type(mct_aVect) , intent(in), optional :: fractions_o
00323     type(mct_aVect) , intent(in), optional :: fractions_a
00324     character(len=*), intent(in), optional :: fluxlist
00325     character(len=*), intent(in), optional :: statelist
00326     !
00327     ! Local variables
00328     !
00329     type(mct_aVect)          :: temp
00330     type(mct_aVect)          :: av_a_f, av_a_s
00331     integer                  :: ka, ko, kSo_t
00332     integer                  :: numats,i,j,ier
00333     integer                  :: lsize
00334     real(R8),allocatable     :: recip(:)
00335     integer                  :: rcnt
00336     real(R8)                 :: rmax,rsum,rval
00337     character(*),parameter :: subName = '(map_ocn2atm_mct) '
00338     !-----------------------------------------------------
00339 
00340     if (samegrid_mapa2o) then
00341 
00342        if (present(fluxlist) .or. present(statelist)) then
00343           if (present(fluxlist)) then
00344              call mct_rearr_rearrange_fldlist(av_o, av_a, Re_ocn2atm, VECTOR=usevector, &
00345                   ALLTOALL=usealltoall, fldlist=fluxlist)
00346           end if
00347           if (present(statelist)) then
00348              call mct_rearr_rearrange_fldlist(av_o, av_a, Re_ocn2atm, VECTOR=usevector, &
00349                   ALLTOALL=usealltoall, fldlist=statelist)
00350           end if
00351        else
00352           call mct_rearr_rearrange(av_o, av_a, Re_ocn2atm, VECTOR=usevector, ALLTOALL=usealltoall)
00353        endif
00354 
00355     else
00356 
00357        ! Normalize input data with fraction of atmosphere on ocean grid
00358           
00359        if (present(fractions_o) .and. present(fractions_a)) then
00360 
00361           lsize = mct_aVect_lsize(av_o)
00362           numats = mct_aVect_nRAttr(av_o)
00363           ko = mct_aVect_indexRA(fractions_o,"ofrac")
00364           call mct_aVect_init(temp, av_o, lsize=lsize)
00365           do j=1,lsize
00366              do i=1,numats
00367                 temp%rAttr(i,j) = av_o%rAttr(i,j)* fractions_o%rAttr(ko,j)
00368              end do
00369           end do
00370           
00371           ! Perform remapping
00372           
00373           if (present(fluxlist) .or. present(statelist)) then
00374              if (present (fluxlist)) then 
00375                 lsize = mct_aVect_lsize(av_a)
00376                 call mct_aVect_init (av_a_f, rlist=fluxlist , lsize=lsize)
00377                 call mct_sMat_avMult(temp, sMatp_Fo2a, av_a_f, VECTOR=usevector, rList=fluxlist )
00378                 call mct_aVect_copy (aVin=av_a_f, aVout=av_a, vector=usevector)
00379                 call mct_aVect_clean(av_a_f)
00380              end if
00381              if (present(statelist)) then 
00382                 lsize = mct_aVect_lsize(av_a)
00383                 call mct_aVect_init (av_a_s, rlist=statelist, lsize=lsize)
00384                 call mct_sMat_avMult(temp, sMatp_So2a, av_a_s, VECTOR=usevector, rList=statelist)
00385                 call mct_aVect_copy (aVin=av_a_s, aVout=av_a, vector=usevector)
00386                 call mct_aVect_clean(av_a_s)
00387              end if
00388           else
00389              ! --- default is flux mapping
00390              call mct_sMat_avMult(temp, sMatp_Fo2a, av_a_f, VECTOR=usevector)
00391           endif
00392 
00393           ! Clean up temporary vector
00394           call mct_aVect_clean(temp)
00395           
00396 
00397           ! Denormalize output data with fraction of ocean on atmosphere grid
00398 
00399           lsize = mct_aVect_lsize(av_a)
00400           numats = mct_aVect_nRAttr(av_a)
00401           allocate(recip(lsize),stat=ier)
00402           if(ier/=0) call die(subName,'allocate recip',ier)
00403           
00404           ko    = mct_aVect_indexRA(fractions_a, "ofrac")
00405           do j=1,lsize
00406              recip(j) = 0.0_R8
00407              if (fractions_a%rAttr(ko,j) /= 0.0_R8) then
00408                 recip(j)= 1.0_R8 / fractions_a%rAttr(ko,j)
00409              end if
00410              do i =1,numats
00411                    av_a%rAttr(i,j) = av_a%rAttr(i,j) * recip(j)
00412              end do
00413           end do
00414 
00415           deallocate(recip,stat=ier)
00416           if(ier/=0) call die(subName,'deallocate recip',ier)
00417 
00418        else
00419          
00420           if (present(fluxlist) .or. present(statelist)) then
00421              if (present (fluxlist)) then 
00422                 lsize = mct_aVect_lsize(av_a)
00423                 call mct_aVect_init (av_a_f, rlist=fluxlist , lsize=lsize)
00424                 call mct_sMat_avMult(av_o, sMatp_Fo2a, av_a_f, VECTOR=usevector, rList=fluxlist )
00425                 call mct_aVect_copy (aVin=av_a_f, aVout=av_a, vector=usevector)
00426                 call mct_aVect_clean(av_a_f)
00427              end if
00428              if (present(statelist)) then 
00429                 lsize = mct_aVect_lsize(av_a)
00430                 call mct_aVect_init (av_a_s, rlist=statelist, lsize=lsize)
00431                 call mct_sMat_avMult(av_o, sMatp_So2a, av_a_s, VECTOR=usevector, rList=statelist)
00432                 call mct_aVect_copy (aVin=av_a_s, aVout=av_a, vector=usevector)
00433                 call mct_aVect_clean(av_a_s)
00434              end if
00435           else
00436              ! --- default is flux mapping
00437              call mct_sMat_avMult(av_o, sMatp_Fo2a, av_a, VECTOR=usevector)
00438           endif
00439           
00440        end if
00441 
00442    end if
00443 
00444 
00445  end subroutine map_ocn2atm_mct
00446 
00447 !=======================================================================
00448 
00449  subroutine map_npFixNew4R(buni,buno,fld1,fld2,gsmapi,domi,domo,ni_i,nj_i,mpicom)
00450 
00451    !===============================================================================
00452    !    Correct the north pole mapping of velocity fields from the atm to ocn
00453    !    grids.  This assumes the input grid is a regular lat/lon with the north
00454    !    pole surrounded by the last latitude line in the input array.  The
00455    !    longitudes in the last latitude must be ordered and equally spaced.
00456    !
00457    !    4R is a low memory version of 3R.
00458    !    This version (New4R) is the same as 3R except it uses a lot less memory
00459    !    and is a bit faster.  Like 3R, it saves data between calls and so
00460    !    assumes the input grid remains constant for all calls.  This is bfb
00461    !    with 3R on bluevista as of 2/15/07.
00462    !
00463    !    !REVISION HISTORY:
00464    !    2007-Feb-12 - T. Craig -- modified New3R to reduce memory
00465    !    2007-Apr-27 - M. Vertenstein - implemented in sequential system
00466    !===============================================================================
00467    
00468 #include <mpif.h>  
00469    ! 
00470    ! Arguments
00471    !
00472    type(mct_Avect),intent(in)   :: buni    ! input  attribute vec
00473    type(mct_Avect),intent(out)  :: buno    ! output attribute vec
00474    character(*)   ,intent(in)   :: fld1    ! name of first input field
00475    character(*)   ,intent(in)   :: fld2    ! name of second input field
00476    integer        ,intent(in)   :: ni_i    ! number of longitudes in input grid
00477    integer        ,intent(in)   :: nj_i    ! number of latitudes in input grid
00478    type(mct_gsMap),pointer      :: gsmapi  ! input gsmap
00479    type(mct_gGrid),pointer      :: domi    ! input domain 
00480    type(mct_gGrid),pointer      :: domo    ! output domain 
00481    integer        ,intent(in)   :: mpicom  ! mpi communicator group 
00482    !
00483    ! Local Variables
00484    !
00485    integer(IN)  :: n,m                           ! generic indices
00486    integer(IN)  :: n1,n2,n3                      ! generic indices
00487    integer(IN)  :: kui,kvi                       ! field indices
00488    integer(IN)  :: kuo,kvo                       ! field indices
00489    integer(IN)  :: kin                           ! index index
00490    integer(IN)  :: nmin,nmax                     ! indices of highest latitude in input
00491    integer(IN)  :: npts                          ! local number of points in an aV
00492    integer(IN)  :: num                           ! number of points at highest latitude
00493    integer(IN)  :: kloni                         ! longitude index on input domain
00494    integer(IN)  :: klati                         ! latitude index on input domain
00495    integer(IN)  :: klono                         ! longitude index on output domain
00496    integer(IN)  :: klato                         ! latitude index on output domain
00497    integer(IN)  :: index                         ! index value
00498    real(R8)     :: rindex                        ! index value
00499    real(R8)     :: latmax                        ! value of highest latitude
00500    real(R8)     :: olon,olat                     ! output bundle lon/lat
00501    real(R8)     :: ilon,ilat                     ! input bundle lon/lat
00502    real(R8)     :: npu,npv                       ! np velocity fields relative to lon
00503    real(R8)     :: theta1,theta2                 ! angles for trig functions
00504    real(R8),allocatable,save :: ilon1(:)         ! lon of input grid at highest latitude
00505    real(R8),allocatable,save :: ilat1(:)         ! lat of input grid at highest latitude
00506    real(R8),allocatable,save :: ilon2(:)         ! lon of input grid at highest latitude
00507    real(R8),allocatable,save :: ilat2(:)         ! lat of input grid at highest latitude
00508    real(R8),allocatable      :: rarray(:)        ! temporary array
00509    real(R8),allocatable      :: rarray2(:,:)     ! temporary array
00510    real(R8)     :: w1,w2,w3,w4                   ! weights
00511    real(R8)     :: f1,f2,f3,f4                   ! function values
00512    real(R8)     :: alpha,beta                    ! used to generate weights
00513    real(R8)     :: rtmp                          ! real temporary
00514    real(R8),allocatable :: lData(:,:)            ! last lat local input bundle data
00515                                                  ! also compressed global data
00516    real(R8)   ,allocatable,save :: alphafound(:) ! list of found alphas
00517    real(R8)   ,allocatable,save :: betafound(:)  ! list of found betas
00518    integer(IN),allocatable,save :: mfound(:)     ! list of found ms
00519    integer(IN),allocatable,save :: nfound(:)     ! list of found ns
00520    real(R8)   ,allocatable      :: rfound(:)     ! temporary for copy
00521    integer(IN),allocatable      :: ifound(:)     ! temporary for copy
00522    integer(IN),save             :: cntfound      ! number found
00523    integer(IN),save             :: cntf_tot      ! cntfound total
00524    integer(IN)  :: cnt                           ! loop counter
00525    logical      :: found                         ! search for new interpolation
00526    integer(IN)  :: rcode                         ! error code
00527    integer(IN)  :: np1                           ! n+1 or tmp
00528    real(R8)     :: ilon1x                        ! tmp
00529    integer(IN),pointer,save  :: gindex(:)        ! global index 
00530    logical,save :: first_call = .true.           ! flags 1st invocation of routine
00531    integer(IN),save   :: tnpf1,tnpf2,tnpf3,tnpf4,tnpf5,tnpf6,tnpf7,tnpf8,tnpf9
00532    real(R8),parameter :: shr_const_deg2rad = shr_const_pi/180.0_R8  ! deg to rads
00533    integer(IN) :: mype
00534 
00535    !--- formats ---
00536    character(*),parameter :: subName = '(map_npFixNew4R) '
00537    character(*),parameter :: F00 = "('(map_npFixNew4R) ',8a)"
00538    character(*),parameter :: F01 = "('(map_npFixNew4R) ',a,i12)"
00539    !-------------------------------------------------------------------------------
00540 
00541    call MPI_COMM_RANK(mpicom,mype,rcode)
00542 
00543    kui   = mct_aVect_indexRA(buni,fld1,perrWith=subName)
00544    kvi   = mct_aVect_indexRA(buni,fld2,perrWith=subName)
00545    kuo   = mct_aVect_indexRA(buno,fld1,perrWith=subName)
00546    kvo   = mct_aVect_indexRA(buno,fld2,perrWith=subName)
00547 
00548 !  tcx 3/19/08, don't use GlobGridNum, it's not set properly in models
00549 !   kin   = mct_aVect_indexIA(domi%data,"GlobGridNum",perrWith=subName)
00550    klati = mct_aVect_indexRA(domi%data,"lat"        ,perrWith=subName)
00551    kloni = mct_aVect_indexRA(domi%data,"lon"        ,perrWith=subName)
00552    klato = mct_aVect_indexRA(domo%data,"lat"        ,perrWith=subName)
00553    klono = mct_aVect_indexRA(domo%data,"lon"        ,perrWith=subName)
00554 
00555    ! ni_i, nj_i should be read in from mapping file
00556 
00557    nmin = (ni_i)*(nj_i-1) + 1
00558    nmax = ni_i*nj_i
00559    num  = ni_i
00560 
00561    !---------------------------------------------------------------------------
00562    ! Initialization only
00563    !---------------------------------------------------------------------------
00564 
00565    if (first_call) then
00566 
00567      if (loglevel > 0) write(logunit,F00) " compute bilinear weights & indicies for NP region."
00568 
00569      allocate(ilon1(num))
00570      allocate(ilon2(num))
00571      allocate(ilat1(num))
00572      allocate(ilat2(num))
00573 
00574      ilon1 = 0._r8
00575      ilon2 = 0._r8
00576      ilat1 = 0._r8
00577      ilat2 = 0._r8
00578 
00579      npts = mct_aVect_lSize(domi%data)
00580      call mct_gsMap_orderedPoints(gsMapi, mype, gindex)
00581      do m=1,npts
00582        if (gindex(m).ge.nmin) then               ! are on highest latitude
00583          n = gindex(m) - nmin + 1                ! n is 1->ni_i lon index on highest latitude
00584          rtmp = domi%data%rAttr(kloni,m)      ! rtmp is longitude value on highest latitude
00585          ilon1(n) = mod(rtmp+360._R8,360._R8) ! ilon1(n) is longitude val mapped from 0->360
00586          ilat1(n) = domi%data%rAttr(klati,m)  ! ilat1(n) values should all be the same (i.e. highest lat)
00587        endif
00588      enddo
00589 
00590      !--- all gather local data, MPI_SUM is low memory and simple
00591      !--- but is a performance penalty compared to gatherv and copy
00592      !--- or a fancy send/recv 
00593 
00594      allocate(rarray(num))
00595      rarray = ilat1
00596      call MPI_ALLREDUCE(rarray,ilat1,num,MPI_REAL8,MPI_SUM,mpicom,rcode)
00597      if (rcode.ne.0) then
00598        write(logunit,*) trim(subName),' ilat1 rcode error ',rcode
00599        call shr_sys_abort()
00600      endif
00601      rarray = ilon1
00602      call MPI_ALLREDUCE(rarray,ilon1,num,MPI_REAL8,MPI_SUM,mpicom,rcode)
00603      if (rcode.ne.0) then
00604        write(logunit,*) trim(subName),' ilon1 rcode error ',rcode
00605        call shr_sys_abort()
00606      endif
00607 
00608      do n = 1,num
00609        np1 = mod(n,num)+1
00610        ilat2(n) = ilat1(np1)
00611        ilon2(n) = ilon1(np1)
00612        if (ilon2(n) < ilon1(n)) ilon2(n) = ilon2(n) + 360._R8
00613      enddo
00614 
00615      do n = 1,num
00616        if (ilat1(n) /= ilat2(n)) then
00617           write(logunit,*) trim(subname),' ERROR: ilat1 ne ilat2 ',n,ilat1(n),ilat2(n)
00618           call shr_sys_abort(trim(subname)//' ERROR: ilat1 ne ilat2')
00619        endif
00620        if (ilon2(n) < ilon1(n)) then
00621           write(logunit,*) trim(subname),' ERROR: ilon2 lt ilon1 ',n,ilon1(n),ilon2(n)
00622           call shr_sys_abort(trim(subname)//' ERROR: ilon2 ilon1 error')
00623        endif
00624        ! tcraig check that lon diffs are reasonable 4x average dlon seems like reasonable limit
00625        if (ilon2(n) - ilon1(n) > (360.0_R8/(num*1.0_R8))*4.0) then
00626           write(logunit,*) trim(subname),' ERROR: ilon1,ilon2 ',n,ilon1(n),ilon2(n)
00627           call shr_sys_abort(trim(subname)//' ERROR: ilon2 ilon1 size diff')
00628        endif
00629      enddo
00630 
00631      latmax = maxval(ilat1)
00632 
00633      !--- compute weights and save them ---
00634 
00635      npts = mct_aVect_lSize(buno)
00636      allocate(mfound(npts),nfound(npts),alphafound(npts),betafound(npts))
00637      cntfound = 0
00638      do m = 1,npts
00639        olat = domo%data%rAttr(klato,m)
00640        if (olat >= latmax) then
00641          rtmp = domo%data%rAttr(klono,m)
00642          olon = mod(rtmp,360._R8)
00643          n = 1
00644          found = .false.
00645          do while (n <= num .and. .not.found )
00646            if (    olon         >= ilon1(n) .and. olon < ilon2(n) .or.   &
00647                    olon+360._R8 >= ilon1(n) .and. olon < ilon2(n)) then
00648 
00649 
00650 !tcx ilat2==ilat1 so don't average
00651 !--->        ilat = (ilat1(n) + ilat2(n)) * 0.5_R8
00652              ilat = ilat1(n)
00653              if (ilon2(n) == ilon1(n)) then
00654                alpha = 0.5_R8
00655              else if (    olon >= ilon1(n) .and. olon < ilon2(n)) then
00656                alpha = (olon - ilon1(n)) / (ilon2(n) - ilon1(n))
00657              else if (olon+360._R8>= ilon1(n) .and. olon < ilon2(n)) then
00658                alpha = (olon+360._R8 - ilon1(n)) / (ilon2(n) - ilon1(n))
00659              else
00660                write(logunit,*) subName,' ERROR: olon ',olon,ilon1(n),ilon2(n)
00661              endif
00662              if (ilat >= 90._R8) then
00663                beta  = 1.0_R8
00664              else
00665                beta  = (olat - ilat) / (90._R8 - ilat)
00666              endif
00667 
00668              cntfound = cntfound + 1
00669              mfound(cntfound) = m
00670              nfound(cntfound) = n
00671              alphafound(cntfound) = alpha
00672              betafound(cntfound) = beta
00673              found = .true.
00674 
00675            endif
00676            n = n + 1     ! normal increment
00677          enddo
00678          if ( .not.found ) then
00679            write(logunit,*) subName,' ERROR: found = false ',found,m,olon,olat
00680          endif
00681        endif
00682      end do
00683 
00684      allocate(ifound(npts))
00685      ifound(1:cntfound) = mfound(1:cntfound)
00686      deallocate(mfound)
00687      if (cntfound > 0) then
00688         allocate(mfound(cntfound))
00689         mfound(1:cntfound) = ifound(1:cntfound)
00690      endif
00691 
00692      ifound(1:cntfound) = nfound(1:cntfound)
00693      deallocate(nfound)
00694      if (cntfound > 0) then
00695         allocate(nfound(cntfound))
00696         nfound(1:cntfound) = ifound(1:cntfound)
00697      endif
00698      deallocate(ifound)
00699 
00700      allocate(rfound(npts))
00701      rfound(1:cntfound) = alphafound(1:cntfound)
00702      deallocate(alphafound)
00703      if (cntfound > 0) then
00704         allocate(alphafound(cntfound))
00705         alphafound(1:cntfound) = rfound(1:cntfound)
00706      endif
00707 
00708      rfound(1:cntfound) = betafound(1:cntfound)
00709      deallocate(betafound)
00710      if (cntfound > 0) then
00711         allocate(betafound(cntfound))
00712         betafound(1:cntfound) = rfound(1:cntfound)
00713      endif
00714      deallocate(rfound)
00715 
00716      call MPI_ALLREDUCE(cntfound,cntf_tot,1,MPI_INTEGER,MPI_SUM,mpicom,rcode)
00717      if (mype == 0) then
00718         write(logunit,F01) ' total npfix points found = ',cntf_tot
00719      endif
00720 
00721      first_call = .false.
00722 
00723    endif
00724 
00725    !---------------------------------------------------------------------------
00726    ! Return if there is nothing to do; must be no points on any pes
00727    ! If there are any npfix points, all pes must continue to the np u,v calc
00728    !---------------------------------------------------------------------------
00729 
00730    if (cntf_tot < 1) then
00731       return
00732    endif
00733 
00734    !---------------------------------------------------------------------------
00735    ! Non-initialization, run-time fix
00736    !---------------------------------------------------------------------------
00737 
00738    !--- barrier not required but interesting for timing. ---
00739    !  call shr_mpi_barrier(mpicom,subName//" barrier")
00740 
00741    !--- extract index, u, v from buni ---
00742 
00743    allocate(lData(3,num))
00744    lData = 0._R8
00745    npts = mct_aVect_lSize(buni)
00746    do n=1,npts
00747      if (gindex(n).ge.nmin) then
00748        m = gindex(n) - nmin + 1
00749        lData(1,m) = gindex(n)
00750        lData(2,m) = buni%rAttr(kui,n)
00751        lData(3,m) = buni%rAttr(kvi,n)
00752      endif
00753    enddo
00754 
00755    !--- all gather local data, MPI_SUM is low memory and simple
00756    !--- but is a performance penalty compared to gatherv and copy
00757    !--- KLUDGE - this should be looked at when it becomes a performance/memory
00758    !--- penalty   
00759 
00760    allocate(rarray2(3,num))
00761    rarray2=lData
00762    call MPI_ALLREDUCE(rarray2,lData,3*num,MPI_REAL8,MPI_SUM,mpicom,rcode)
00763    deallocate(rarray2)
00764 
00765    if (rcode.ne.0) then
00766      write(logunit,*) trim(subName),' rcode error ',rcode
00767      call shr_sys_abort()
00768    endif
00769 
00770    do n2=1,num
00771      if (lData(1,n2).lt.0.1_R8) then
00772        write(logunit,*) trim(subName),' error allreduce ',n2
00773      endif
00774    enddo
00775 
00776    !--- compute npu, npv (pole data) and initialize ilon,ilat arrays ---
00777 
00778    npu = 0._R8
00779    npv = 0._R8
00780    do n = 1,num
00781      theta1 = ilon1(n)*shr_const_deg2rad
00782      npu = npu + cos(theta1)*lData(2,n) &
00783                - sin(theta1)*lData(3,n)
00784      npv = npv + sin(theta1)*lData(2,n) &
00785                + cos(theta1)*lData(3,n)
00786    enddo
00787    npu = npu / real(num,R8)
00788    npv = npv / real(num,R8)
00789 
00790    !--- compute updated pole vectors ---
00791 
00792 !DIR$ CONCURRENT
00793    do cnt = 1,cntfound
00794       m     = mfound(cnt)
00795       n     = nfound(cnt)
00796       np1   = mod(n,num)+1
00797       alpha = alphafound(cnt)
00798       beta  = betafound(cnt)
00799 
00800       w1 = (1.0_R8-alpha)*(1.0_R8-beta)
00801       w2 = (    alpha)*(1.0_R8-beta)
00802       w3 = (    alpha)*(    beta)
00803       w4 = (1.0_R8-alpha)*(    beta)
00804 
00805       theta1 = ilon1(n)*shr_const_deg2rad
00806       theta2 = ilon2(n)*shr_const_deg2rad
00807 
00808       f1 = lData(2,n)
00809       f2 = lData(2,np1)
00810       f3 =  cos(theta1)*npu + sin(theta1)*npv
00811       f4 =  cos(theta2)*npu + sin(theta2)*npv
00812       rtmp = w1*f1 + w2*f2 + w3*f3 + w4*f4
00813       buno%rAttr(kuo,m) = w1*f1 + w2*f2 + w3*f3 + w4*f4
00814 
00815       f1 = lData(3,n)
00816       f2 = lData(3,np1)
00817       f3 = -sin(theta1)*npu + cos(theta1)*npv
00818       f4 = -sin(theta2)*npu + cos(theta2)*npv
00819       rtmp = w1*f1 + w2*f2 + w3*f3 + w4*f4
00820       buno%rAttr(kvo,m) = w1*f1 + w2*f2 + w3*f3 + w4*f4
00821    enddo
00822 
00823    deallocate(lData)
00824    first_call = .false.
00825 
00826 end subroutine map_npFixNew4R
00827 
00828 !===============================================================================
00829 
00830 end module map_atmocn_mct