CESM_Driver: map_atmlnd_mct.F90 Source File

00001 module map_atmlnd_mct
00002 
00003 !---------------------------------------------------------------------
00004 !
00005 ! Purpose:
00006 !
00007 ! Collect coupling routines for sequential coupling of LND-ATM.
00008 !       
00009 !
00010 ! Author: R. Jacob, M. Vertenstein
00011 ! Revision History:
00012 ! 30Mar06 - P. Worley - added optional arguments to MCT_Rearrange call
00013 ! 13Apr06 - M. Vertenstein - cleaned up interfaces 
00014 !
00015 !---------------------------------------------------------------------
00016 
00017   use shr_kind_mod      ,only: R8 => SHR_KIND_R8, IN=>SHR_KIND_IN
00018   use shr_sys_mod
00019   use shr_const_mod
00020   use shr_mct_mod, only: shr_mct_sMatPInitnc
00021   use mct_mod
00022 
00023   use seq_comm_mct, only : logunit, loglevel
00024   use seq_cdata_mod
00025   use seq_flds_indices
00026   use seq_infodata_mod
00027   use m_die
00028 
00029   implicit none
00030   save
00031   private  ! except
00032 
00033 !--------------------------------------------------------------------------
00034 ! Public interfaces
00035 !--------------------------------------------------------------------------
00036 
00037   public :: map_lnd2atm_init_mct
00038   public :: map_atm2lnd_init_mct
00039   public :: map_lnd2atm_mct
00040   public :: map_atm2lnd_mct
00041 
00042 !--------------------------------------------------------------------------
00043 ! Public data
00044 !--------------------------------------------------------------------------
00045 
00046 !--------------------------------------------------------------------------
00047 ! Private data
00048 !--------------------------------------------------------------------------
00049 
00050   type(mct_rearr), private :: Re_lnd2atm
00051   type(mct_rearr), private :: Re_atm2lnd
00052   type(mct_sMatp), private :: sMatp_Fa2l
00053   type(mct_sMatp), private :: sMatp_Sa2l
00054   type(mct_sMatp), private :: sMatp_Fl2a
00055   type(mct_sMatp), private :: sMatp_Sl2a
00056 
00057 #ifdef CPP_VECTOR
00058     logical :: usevector = .true.
00059 #else
00060     logical :: usevector = .false.
00061 #endif
00062 
00063 #ifdef SYSUNICOS
00064     logical :: usealltoall = .true.
00065 #else
00066     logical :: usealltoall = .false.
00067 #endif
00068   logical, private :: samegrid_mapa2l
00069 
00070 !=======================================================================
00071    contains
00072 !=======================================================================
00073 
00074   subroutine map_atm2lnd_init_mct( cdata_a, cdata_l)
00075 
00076     !--------------------------------------------------
00077     ! 
00078     ! Arguments
00079     !
00080     type(seq_cdata),intent(in) :: cdata_a
00081     type(seq_cdata),intent(in) :: cdata_l
00082     ! 
00083     ! Local variables
00084     !
00085     type(seq_infodata_type), pointer :: infodata
00086     type(mct_gsMap), pointer :: gsMap_a           ! atm gsMap
00087     type(mct_gsMap), pointer :: gsMap_l           ! lnd gsMap
00088     type(mct_gGrid), pointer :: dom_a             ! atm domain
00089     type(mct_gGrid), pointer :: dom_l             ! lnd domain
00090     integer                  :: mpicom            ! communicator spanning atm and lnd
00091     integer                  :: ka, km            ! indices
00092     integer                  :: lsize             ! size of attribute vector
00093     type(mct_aVect)          :: areasrc           ! atm areas from mapping file
00094     type(mct_aVect)          :: areadst           ! lnd areas set to atm areas
00095     character(*),parameter :: subName = '(map_atm2lnd_init_mct) '
00096     !--------------------------------------------------
00097 
00098     call seq_cdata_setptrs(cdata_a, gsMap=gsMap_a, dom=dom_a)
00099     call seq_cdata_setptrs(cdata_l, gsMap=gsMap_l, dom=dom_l)
00100     call seq_cdata_setptrs(cdata_a, mpicom=mpicom, infodata=infodata)
00101 
00102     call seq_infodata_GetData( infodata, samegrid_al=samegrid_mapa2l)
00103 
00104     lsize = mct_gsMap_lsize(gsMap_a, mpicom)
00105     call mct_aVect_init( areasrc, rList="aream", lsize=lsize )
00106        
00107     lsize = mct_gsMap_lsize(gsMap_l, mpicom)
00108     call mct_aVect_init( areadst, rList="aream", lsize=lsize )
00109 
00110     if (samegrid_mapa2l) then
00111 
00112        call mct_rearr_init(gsMap_a, gsMap_l, mpicom, Re_atm2lnd)
00113 
00114        ! --- copy atm area into land aream
00115 
00116        ka = mct_aVect_indexRA(areasrc   , "aream")
00117        km = mct_aVect_indexRa(dom_a%data, "aream" )
00118        areasrc%rAttr(ka,:) = dom_a%data%rAttr(km,:)
00119        call mct_rearr_rearrange(areasrc, areadst, Re_atm2lnd, VECTOR=usevector, &
00120           ALLTOALL=usealltoall)
00121 
00122     else
00123 
00124        call shr_mct_sMatPInitnc(sMatp_Fa2l,gsMap_a,gsMap_l,"seq_maps.rc", &
00125           "atm2lndFmapname:","atm2lndFmaptype:",mpicom, &
00126           areasrc=areasrc, areadst=areadst)
00127 
00128        call shr_mct_sMatPInitnc(sMatp_Sa2l,gsMap_a,gsMap_l,"seq_maps.rc", &
00129           "atm2lndSmapname:","atm2lndSmaptype:",mpicom)
00130 
00131     endif
00132 
00133     ka = mct_aVect_indexRA(areadst   ,"aream")
00134     km = mct_aVect_indexRA(dom_l%data,"aream")
00135     dom_l%data%rAttr(km,:) = areadst%rAttr(ka,:)
00136 
00137     call mct_aVect_clean(areasrc)      
00138     call mct_aVect_clean(areadst)      
00139 
00140   end subroutine map_atm2lnd_init_mct
00141 
00142 !=======================================================================
00143 
00144   subroutine map_lnd2atm_init_mct( cdata_l, cdata_a )
00145 
00146     !--------------------------------------------------
00147     ! 
00148     ! Arguments
00149     !
00150     type(seq_cdata),intent(in) :: cdata_l
00151     type(seq_cdata),intent(in) :: cdata_a
00152     ! 
00153     ! Local variables
00154     !
00155     type(seq_infodata_type), pointer :: infodata
00156     type(mct_gsMap), pointer :: gsMap_a           ! atm gsMap
00157     type(mct_gsMap), pointer :: gsMap_l           ! lnd gsMap
00158     type(mct_gGrid), pointer :: dom_l             ! lnd domain
00159     integer                  :: kf,iam,ierr,lsize
00160     integer                  :: mpicom            ! communicator spanning atm and lnd
00161     character(*),parameter :: subName = '(map_lnd2atm_init_mct) '
00162     !--------------------------------------------------
00163 
00164     call seq_cdata_setptrs(cdata_l, gsMap=gsMap_l, dom=dom_l)
00165     call seq_cdata_setptrs(cdata_a, gsMap=gsMap_a)
00166     call seq_cdata_setptrs(cdata_a, mpicom=mpicom,infodata=infodata)
00167     call mpi_comm_rank(mpicom,iam,ierr)
00168     
00169     call seq_infodata_GetData( infodata, samegrid_al=samegrid_mapa2l)
00170 
00171     ! Initialize lnd -> atm mapping or rearranger
00172 
00173     if (samegrid_mapa2l) then
00174 
00175        call mct_rearr_init(gsMap_l, gsMap_a, mpicom, Re_lnd2atm)
00176 
00177     else
00178 
00179        call shr_mct_sMatPInitnc(sMatp_Fl2a, gsMap_l, gsMap_a, "seq_maps.rc", &
00180             "lnd2atmFmapname:", "lnd2atmFmaptype:", mpicom)
00181 
00182        call shr_mct_sMatPInitnc(sMatp_Sl2a, gsMap_l, gsMap_a, "seq_maps.rc", &
00183             "lnd2atmSmapname:", "lnd2atmSmaptype:", mpicom)
00184 
00185     endif
00186 
00187   end subroutine map_lnd2atm_init_mct
00188 
00189 !=======================================================================
00190 
00191   subroutine map_atm2lnd_mct( cdata_a, av_a, cdata_l, av_l, fluxlist, statelist )
00192 
00193     !--------------------------------------------------
00194     ! 
00195     ! Arguments
00196     !
00197     type(seq_cdata), intent(in)  :: cdata_a
00198     type(mct_aVect), intent(in)  :: av_a
00199     type(seq_cdata), intent(in)  :: cdata_l
00200     type(mct_aVect), intent(out) :: av_l
00201     character(len=*),intent(in), optional :: statelist
00202     character(len=*),intent(in), optional :: fluxlist
00203     !
00204     ! Local
00205     ! 
00206     integer                :: lsize
00207     type(mct_aVect)        :: av_l_f     ! temporary flux attribute vector
00208     type(mct_aVect)        :: av_l_s     ! temporary state attribute vector
00209     character(*),parameter :: subName = '(map_atm2lnd_mct) '
00210     !--------------------------------------------------
00211 
00212     if (samegrid_mapa2l) then
00213 
00214        if (present(fluxlist) .or. present(statelist)) then
00215           if (present(fluxlist)) then
00216              call mct_rearr_rearrange_fldlist(av_a, av_l, Re_atm2lnd, VECTOR=usevector, &
00217                   ALLTOALL=usealltoall, fldlist=fluxlist)
00218           endif
00219           if (present(statelist)) then
00220               call mct_rearr_rearrange_fldlist(av_a, av_l, Re_atm2lnd, VECTOR=usevector, &
00221                   ALLTOALL=usealltoall, fldlist=statelist)
00222           endif
00223        else
00224           call mct_rearr_rearrange(av_a, av_l, Re_atm2lnd, VECTOR=usevector, ALLTOALL=usealltoall)
00225        end if
00226 
00227     else
00228 
00229        if (present(fluxlist) .or. present(statelist)) then
00230           if (present(fluxlist)) then
00231              lsize = mct_aVect_lsize(av_l)
00232              call mct_aVect_init (av_l_f, rlist=fluxlist , lsize=lsize)
00233              call mct_sMat_avMult(av_a, sMatp_Fa2l, av_l_f, VECTOR=usevector, rList=fluxlist)
00234              call mct_aVect_copy (aVin=av_l_f, aVout=av_l, vector=usevector)
00235              call mct_aVect_clean(av_l_f)
00236           end if
00237           if (present(statelist)) then
00238              lsize = mct_aVect_lsize(av_l)
00239              call mct_aVect_init (av_l_s, rlist=statelist, lsize=lsize)
00240              call mct_sMat_avMult(av_a, sMatp_Sa2l, av_l_s, VECTOR=usevector, rList=statelist)
00241              call mct_aVect_copy (aVin=av_l_s, aVout=av_l, vector=usevector)
00242              call mct_aVect_clean(av_l_s)
00243           end if
00244        else
00245           !--- default is flux mapping
00246           call mct_sMat_avMult(av_a, sMatp_Fa2l, av_l, VECTOR=usevector)
00247        endif
00248     endif
00249        
00250   end subroutine map_atm2lnd_mct
00251 
00252 !=======================================================================
00253 
00254   subroutine map_lnd2atm_mct( cdata_l, av_l, cdata_a, av_a, &
00255                               fractions_l, fractions_a, &
00256                               fluxlist, statelist )
00257 
00258     !--------------------------------------------------
00259     ! 
00260     ! Arguments
00261     !
00262     type(seq_cdata) ,intent(in)  :: cdata_l
00263     type(mct_aVect) ,intent(in)  :: av_l
00264     type(seq_cdata) ,intent(in)  :: cdata_a
00265     type(mct_aVect) ,intent(out) :: av_a
00266     type(mct_aVect) ,intent(in), optional :: fractions_l
00267     type(mct_aVect) ,intent(in), optional :: fractions_a
00268     character(len=*),intent(in), optional :: fluxlist
00269     character(len=*),intent(in), optional :: statelist
00270     !
00271     ! Local
00272     !
00273     integer  :: i,j,kl,lsize,numats,ier
00274     real(r8) :: recip
00275     type(mct_aVect)          :: av_a_f     ! temporary flux attribute vector
00276     type(mct_aVect)          :: av_a_s     ! temporary state attribute vector
00277     type(mct_aVect)          :: temp       ! temporary attribute vector
00278     character(*),parameter :: subName = '(map_lnd2atm_mct) '
00279     !--------------------------------------------------
00280 
00281     if (samegrid_mapa2l) then
00282 
00283        if (present(fluxlist) .or. present(statelist)) then
00284           if (present(fluxlist)) then
00285              call mct_rearr_rearrange_fldlist(av_l, av_a, Re_lnd2atm, VECTOR=usevector, &
00286                   ALLTOALL=usealltoall, fldlist=fluxlist)
00287           endif
00288           if (present(statelist)) then
00289              call mct_rearr_rearrange_fldlist(av_l, av_a, Re_lnd2atm, VECTOR=usevector, &
00290                   ALLTOALL=usealltoall, fldlist=statelist)
00291           endif
00292        else
00293           call mct_rearr_rearrange(av_l, av_a, Re_lnd2atm, VECTOR=usevector, ALLTOALL=usealltoall)
00294        end if
00295 
00296     else
00297 
00298        ! Normalize input data with fraction of land on land grid
00299 
00300        if (present(fractions_l) .and. present(fractions_a)) then
00301 
00302           lsize = mct_aVect_lsize(av_l)
00303           numats = mct_aVect_nRAttr(av_l)
00304           kl = mct_aVect_indexRA(fractions_l,"lfrin")
00305           call mct_aVect_init(temp, av_l, lsize=lsize)
00306           do j=1,lsize
00307              do i=1,numats
00308                 temp%rAttr(i,j) = av_l%rAttr(i,j)* fractions_l%rAttr(kl,j)
00309              end do
00310           end do
00311 
00312           ! Perform remapping
00313 
00314           if (present(fluxlist) .or. present(statelist)) then
00315              if (present (fluxlist)) then
00316                 lsize = mct_aVect_lsize(av_a)
00317                 call mct_aVect_init (av_a_f, rlist=fluxlist , lsize=lsize)
00318                 call mct_sMat_avMult(temp, sMatp_Fl2a, av_a_f, VECTOR=usevector, rList=fluxlist )
00319                 call mct_aVect_copy (aVin=av_a_f, aVout=av_a, vector=usevector)
00320                 call mct_aVect_clean(av_a_f)
00321              end if
00322              if (present(statelist)) then
00323                 lsize = mct_aVect_lsize(av_a)
00324                 call mct_aVect_init (av_a_s, rlist=statelist, lsize=lsize)
00325                 call mct_sMat_avMult(temp, sMatp_Sl2a, av_a_s, VECTOR=usevector, rList=statelist)
00326                 call mct_aVect_copy (aVin=av_a_s, aVout=av_a, vector=usevector)
00327                 call mct_aVect_clean(av_a_s)
00328              end if
00329           else
00330              ! --- default is flux mapping
00331              call mct_sMat_avMult(temp, sMatp_Fl2a, av_a_f, VECTOR=usevector)
00332           endif
00333 
00334           ! Clean up temporary vector
00335           call mct_aVect_clean(temp)
00336 
00337           ! Denormalize output data with fraction of land on atmosphere grid 
00338           lsize = mct_aVect_lsize(av_a)
00339           numats = mct_aVect_nRAttr(av_a)
00340 
00341           kl    = mct_aVect_indexRA(fractions_a, "lfrin")
00342           do j=1,lsize
00343              if (fractions_a%rAttr(kl,j) /= 0.0_R8) then
00344                 recip = 1.0_R8 / fractions_a%rAttr(kl,j)
00345              else
00346                 recip = 0.0_R8
00347              end if
00348              do i =1,numats
00349                    av_a%rAttr(i,j) = av_a%rAttr(i,j) * recip
00350              end do
00351           end do
00352 
00353        else
00354 
00355           if (present(fluxlist) .or. present(statelist)) then
00356              if (present (fluxlist)) then
00357                 lsize = mct_aVect_lsize(av_a)
00358                 call mct_aVect_init (av_a_f, rlist=fluxlist , lsize=lsize)
00359                 call mct_sMat_avMult(av_l, sMatp_Fl2a, av_a_f, VECTOR=usevector, rList=fluxlist )
00360                 call mct_aVect_copy (aVin=av_a_f, aVout=av_a, vector=usevector)
00361                 call mct_aVect_clean(av_a_f)
00362              end if
00363              if (present(statelist)) then
00364                 lsize = mct_aVect_lsize(av_a)
00365                 call mct_aVect_init (av_a_s, rlist=statelist, lsize=lsize)
00366                 call mct_sMat_avMult(av_l, sMatp_Sl2a, av_a_s, VECTOR=usevector, rList=statelist)
00367                 call mct_aVect_copy (aVin=av_a_s, aVout=av_a, vector=usevector)
00368                 call mct_aVect_clean(av_a_s)
00369              end if
00370           else
00371              ! --- default is flux mapping
00372              call mct_sMat_avMult(av_l, sMatp_Fl2a, av_a, VECTOR=usevector)
00373           endif
00374 
00375        endif
00376 
00377     endif
00378 
00379   end subroutine map_lnd2atm_mct
00380 
00381 end module map_atmlnd_mct