CESM_Driver: seq_flux_mct.F90 Source File

00001 module seq_flux_mct
00002   
00003   use shr_kind_mod,      only: r8 => shr_kind_r8, in=>shr_kind_in
00004   use shr_sys_mod,       only: shr_sys_abort
00005   use shr_flux_mod,      only: shr_flux_atmocn
00006   use shr_orb_mod,       only: shr_orb_params, shr_orb_cosz, shr_orb_decl
00007 
00008   use mct_mod
00009   use seq_flds_indices
00010   use seq_comm_mct
00011   use seq_cdata_mod
00012   use seq_infodata_mod
00013 
00014   implicit none
00015   private       
00016   save
00017 
00018 !--------------------------------------------------------------------------
00019 ! Public interfaces
00020 !--------------------------------------------------------------------------
00021 
00022   public seq_flux_init_mct
00023   public seq_flux_atmocn_mct
00024 
00025 !--------------------------------------------------------------------------
00026 ! Private data
00027 !--------------------------------------------------------------------------
00028 
00029   real(r8), pointer     :: lats(:)    ! latitudes  (degrees)
00030   real(r8), pointer     :: lons(:)    ! longitudes (degrees)
00031   integer(in) , allocatable ::  mask (:)  ! ocn domain mask: 0 <=> inactive cell
00032 
00033   real(r8), allocatable ::  uocn (:)  ! ocn velocity, zonal
00034   real(r8), allocatable ::  vocn (:)  ! ocn velocity, meridional
00035   real(r8), allocatable ::  tocn (:)  ! ocean temperature
00036   real(r8), allocatable ::  zbot (:)  ! atm level height
00037   real(r8), allocatable ::  ubot (:)  ! atm velocity, zonal     
00038   real(r8), allocatable ::  vbot (:)  ! atm velocity, meridional
00039   real(r8), allocatable ::  thbot(:)  ! atm potential T
00040   real(r8), allocatable ::  shum (:)  ! atm specific humidity
00041   real(r8), allocatable ::  dens (:)  ! atm density
00042   real(r8), allocatable ::  tbot (:)  ! atm bottom surface T
00043   real(r8), allocatable ::  sen  (:)  ! heat flux: sensible 
00044   real(r8), allocatable ::  lat  (:)  ! heat flux: latent   
00045   real(r8), allocatable ::  lwup (:)  ! lwup over ocean
00046   real(r8), allocatable ::  evap (:)  ! water flux: evaporation
00047   real(r8), allocatable ::  taux (:)  ! wind stress, zonal
00048   real(r8), allocatable ::  tauy (:)  ! wind stress, meridional
00049   real(r8), allocatable ::  tref (:)  ! diagnostic:  2m ref T
00050   real(r8), allocatable ::  qref (:)  ! diagnostic:  2m ref Q
00051   real(r8), allocatable :: duu10n(:)  ! diagnostic: 10m wind speed squared
00052 
00053   real(r8), allocatable ::  ustar(:)  ! saved ustar
00054   real(r8), allocatable ::  re   (:)  ! saved re
00055   real(r8), allocatable ::  ssq  (:)  ! saved sq
00056 
00057   ! Conversion from degrees to radians
00058 
00059   real(r8),parameter :: const_pi      = SHR_CONST_PI       ! pi
00060   real(r8),parameter :: const_deg2rad = const_pi/180.0_R8  ! deg to rads
00061 
00062 !===============================================================================
00063 contains
00064 !===============================================================================
00065 
00066   subroutine seq_flux_init_mct( cdata, xao)
00067 
00068     !-----------------------------------------------------------------------
00069     !
00070     ! Arguments
00071     !
00072     type(seq_cdata), intent(in)  :: cdata
00073     type(mct_aVect), intent(out) :: xao  
00074     !
00075     ! Local variables
00076     !
00077 
00078     integer(in)                     :: nloc
00079     type(mct_gGrid), pointer        :: dom
00080     type(mct_gsMap), pointer        :: gsMap
00081     integer                         :: lsize
00082     integer                         :: mpicom
00083     integer                         :: ier
00084     real(r8), pointer     ::  rmask(:)  ! ocn domain mask
00085     character(*),parameter :: subName =   '(seq_flux_init_mct) '
00086     !-----------------------------------------------------------------------
00087 
00088     ! Set cdata pointers
00089     call seq_cdata_setptrs(cdata, gsMap=gsMap, dom=dom, mpicom=mpicom)
00090     lsize=mct_gsMap_lsize(gsMap, mpicom)
00091 
00092     ! Initialize attribute vector
00093     call mct_aVect_init(xao, rList=seq_flds_xao_fields, lsize=lsize)
00094     call mct_aVect_zero(xao)
00095 
00096     ! Determine size of buffers
00097     nloc = mct_avect_lSize(xao)
00098     
00099     ! Input fields atm
00100     allocate( zbot(nloc),stat=ier)
00101     if(ier/=0) call mct_die(subName,'allocate zbot',ier)
00102     allocate( ubot(nloc),stat=ier)
00103     if(ier/=0) call mct_die(subName,'allocate ubot',ier)
00104     allocate( vbot(nloc))
00105     if(ier/=0) call mct_die(subName,'allocate vbot',ier)
00106     allocate(thbot(nloc),stat=ier)
00107     if(ier/=0) call mct_die(subName,'allocate thbot',ier)
00108     allocate(shum(nloc),stat=ier)
00109     if(ier/=0) call mct_die(subName,'allocate shum',ier)
00110     allocate(dens(nloc),stat=ier)
00111     if(ier/=0) call mct_die(subName,'allocate dens',ier)
00112     allocate(tbot(nloc),stat=ier)
00113     if(ier/=0) call mct_die(subName,'allocate tbot',ier)
00114     allocate(ustar(nloc),stat=ier)
00115     if(ier/=0) call mct_die(subName,'allocate ustar',ier)
00116     allocate(re(nloc), stat=ier)
00117     if(ier/=0) call mct_die(subName,'allocate re',ier)
00118     allocate(ssq(nloc),stat=ier)
00119     if(ier/=0) call mct_die(subName,'allocate ssq',ier)
00120     allocate( uocn(nloc),stat=ier)
00121     if(ier/=0) call mct_die(subName,'allocate uocn',ier)
00122     allocate( vocn(nloc),stat=ier)
00123     if(ier/=0) call mct_die(subName,'allocate vocn',ier)
00124     allocate( tocn(nloc),stat=ier)
00125     if(ier/=0) call mct_die(subName,'allocate tocn',ier)
00126 
00127     ! Output fields 
00128     allocate(sen (nloc),stat=ier)
00129     if(ier/=0) call mct_die(subName,'allocate sen',ier)
00130     allocate(lat (nloc),stat=ier)
00131     if(ier/=0) call mct_die(subName,'allocate lat',ier)
00132     allocate(evap(nloc),stat=ier)
00133     if(ier/=0) call mct_die(subName,'allocate evap',ier)
00134     allocate(lwup(nloc),stat=ier)
00135     if(ier/=0) call mct_die(subName,'allocate lwup',ier)
00136     allocate(taux(nloc),stat=ier)
00137     if(ier/=0) call mct_die(subName,'allocate taux',ier)
00138     allocate(tauy(nloc),stat=ier)
00139     if(ier/=0) call mct_die(subName,'allocate tauy',ier)
00140     allocate(tref(nloc),stat=ier)
00141     if(ier/=0) call mct_die(subName,'allocate tref',ier)
00142     allocate(qref(nloc),stat=ier)
00143     if(ier/=0) call mct_die(subName,'allocate qref',ier)
00144     allocate(duu10n(nloc),stat=ier)
00145     if(ier/=0) call mct_die(subName,'allocate duu10n',ier)
00146 
00147     ! Grid fields
00148     allocate( lats(nloc),stat=ier )
00149     if(ier/=0) call mct_die(subName,'allocate lats',ier)
00150     allocate( lons(nloc),stat=ier )
00151     if(ier/=0) call mct_die(subName,'allocate lons',ier)
00152     allocate(mask(nloc),stat=ier)
00153     if(ier/=0) call mct_die(subName,'allocate mask',ier)
00154     
00155     ! Get lat, lon, mask, which is time-invariant
00156     allocate(rmask(nloc),stat=ier)
00157     if(ier/=0) call mct_die(subName,'allocate rmask',ier)
00158     call mct_gGrid_exportRAttr(dom, 'lat' , lats , nloc) 
00159     call mct_gGrid_exportRAttr(dom, 'lon' , lons , nloc) 
00160     call mct_gGrid_exportRAttr(dom, 'mask', rmask, nloc)
00161     mask = nint(rmask)
00162     deallocate(rmask)
00163 
00164   end subroutine seq_flux_init_mct
00165 
00166 !===============================================================================
00167 
00168   subroutine seq_flux_atmocn_mct( cdata, a2x, o2x, xao, fractions_o, albedo_only )
00169 
00170     !-----------------------------------------------------------------------
00171     !
00172     ! Arguments
00173     !
00174     type(seq_cdata), intent(in)    :: cdata
00175     type(mct_aVect), intent(in)    :: a2x  
00176     type(mct_aVect), intent(in)    :: o2x
00177     type(mct_aVect), intent(inout) :: xao
00178     type(mct_aVect), intent(inout) :: fractions_o
00179     logical, optional,intent(in)  :: albedo_only
00180     !
00181     ! Local variables
00182     !
00183     type(seq_infodata_type),pointer :: infodata
00184     integer(in) :: n,i          ! indices
00185     real(r8)    :: rlat         ! gridcell latitude in radians
00186     real(r8)    :: rlon         ! gridcell longitude in radians
00187     real(r8)    :: cosz         ! Cosine of solar zenith angle
00188     real(r8)    :: eccen        ! Earth orbit eccentricity
00189     real(r8)    :: mvelpp       ! Earth orbit
00190     real(r8)    :: lambm0       ! Earth orbit
00191     real(r8)    :: obliqr       ! Earth orbit
00192     real(r8)    :: delta        ! Solar declination angle  in radians
00193     real(r8)    :: eccf         ! Earth orbit eccentricity factor
00194     real(r8)    :: calday       ! calendar day including fraction, at 0e
00195     real(r8)    :: nextsw_cday  ! calendar day of next atm shortwave
00196     real(r8)    :: anidr        ! albedo: near infrared, direct
00197     real(r8)    :: avsdr        ! albedo: visible      , direct
00198     real(r8)    :: anidf        ! albedo: near infrared, diffuse
00199     real(r8)    :: avsdf        ! albedo: visible      , diffuse
00200     integer(in) :: nloc         ! number of gridcells
00201     integer(in) :: ID           ! comm ID
00202     logical     :: flux_albav   ! flux avg option
00203     logical     :: dead_comps   ! .true.  => dead components are used
00204     integer     :: kx,kr        ! fractions indices
00205     !
00206     real(R8),parameter :: albdif = 0.06_R8 ! 60 deg reference albedo, diffuse
00207     real(R8),parameter :: albdir = 0.07_R8 ! 60 deg reference albedo, direct 
00208     character(*),parameter :: subName =   '(seq_flux_atmocn_mct) '
00209     !
00210     !-----------------------------------------------------------------------
00211 
00212     call seq_cdata_setptrs(cdata, infodata=infodata, ID=ID)
00213     call seq_infodata_GetData(infodata,flux_albav=flux_albav)
00214 
00215     nloc = mct_aVect_lsize(xao)
00216 
00217     if (flux_albav) then
00218        do n=1,nloc   
00219           anidr = albdir
00220           avsdr = albdir
00221           anidf = albdif
00222           avsdf = albdif
00223 
00224           ! Albedo is now function of latitude (will be new implementation)
00225           !rlat = const_deg2rad * lats(n)
00226           !anidr = 0.069_R8 - 0.011_R8 * cos(2._R8 * rlat)
00227           !avsdr = anidr
00228           !anidf = anidr
00229           !avsdf = anidr
00230 
00231           xao%rAttr(index_xao_So_avsdr,n) = avsdr
00232           xao%rAttr(index_xao_So_anidr,n) = anidr
00233           xao%rAttr(index_xao_So_avsdf,n) = avsdf
00234           xao%rAttr(index_xao_So_anidf,n) = anidf
00235 
00236        end do
00237        kx = mct_aVect_indexRA(fractions_o,"ifrac")
00238        kr = mct_aVect_indexRA(fractions_o,"ifrad")
00239        fractions_o%rAttr(kr,:) = fractions_o%rAttr(kx,:)
00240        kx = mct_aVect_indexRA(fractions_o,"ofrac")
00241        kr = mct_aVect_indexRA(fractions_o,"ofrad")
00242        fractions_o%rAttr(kr,:) = fractions_o%rAttr(kx,:)
00243 
00244     else
00245        ! Solar declination 
00246        ! Will only do albedo calculation if nextsw_cday is not -1.
00247        
00248        call seq_infodata_GetData(infodata,nextsw_cday=nextsw_cday,orb_eccen=eccen, &
00249           orb_mvelpp=mvelpp, orb_lambm0=lambm0, orb_obliqr=obliqr)
00250        if (nextsw_cday >= -0.5_r8) then
00251           calday = nextsw_cday
00252           call shr_orb_decl(calday, eccen, mvelpp,lambm0, obliqr, delta, eccf)
00253           ! Compute albedos 
00254           do n=1,nloc
00255              rlat = const_deg2rad * lats(n)
00256              rlon = const_deg2rad * lons(n)
00257              cosz = shr_orb_cosz( calday, rlat, rlon, delta )
00258              if (cosz  >  0.0_R8) then !--- sun hit --
00259                 anidr = (.026_R8/(cosz**1.7_R8 + 0.065_R8)) +   &
00260                         (.150_R8*(cosz         - 0.100_R8 ) *   &
00261                                  (cosz         - 0.500_R8 ) *   &
00262                                  (cosz         - 1.000_R8 )  )
00263                 avsdr = anidr
00264                 anidf = albdif
00265                 avsdf = albdif
00266              else !--- dark side of earth ---
00267                 anidr = 1.0_R8
00268                 avsdr = 1.0_R8
00269                 anidf = 1.0_R8
00270                 avsdf = 1.0_R8
00271              end if
00272              xao%rAttr(index_xao_So_avsdr,n) = avsdr
00273              xao%rAttr(index_xao_So_anidr,n) = anidr
00274              xao%rAttr(index_xao_So_avsdf,n) = avsdf
00275              xao%rAttr(index_xao_So_anidf,n) = anidf
00276           end do   ! nloc
00277           kx = mct_aVect_indexRA(fractions_o,"ifrac")
00278           kr = mct_aVect_indexRA(fractions_o,"ifrad")
00279           fractions_o%rAttr(kr,:) = fractions_o%rAttr(kx,:)
00280           kx = mct_aVect_indexRA(fractions_o,"ofrac")
00281           kr = mct_aVect_indexRA(fractions_o,"ofrad")
00282           fractions_o%rAttr(kr,:) = fractions_o%rAttr(kx,:)
00283        endif    ! nextsw_cday
00284     end if   ! flux_albav
00285        
00286     if (present(albedo_only)) then
00287        if (albedo_only) then
00288           if (seq_comm_iamroot(ID)) &
00289              write(logunit,'(2A)') trim(subname),' computing only ocn albedos '
00290           RETURN
00291        endif
00292     endif
00293 
00294     !-----------------------------------------------------
00295 
00296     ! Update ocean surface fluxes 
00297     ! Must fabricate "reasonable" data (using dead components)
00298 
00299     call seq_infodata_GetData(infodata, dead_comps=dead_comps)
00300 
00301     if (dead_comps) then
00302        do n = 1,nloc
00303           mask(n) =   1      ! ocn domain mask            ~ 0 <=> inactive cell
00304           tocn(n) = 290.0_R8 ! ocn temperature            ~ Kelvin
00305           uocn(n) =   0.0_R8 ! ocn velocity, zonal        ~ m/s
00306           vocn(n) =   0.0_R8 ! ocn velocity, meridional   ~ m/s
00307           zbot(n) =  55.0_R8 ! atm height of bottom layer ~ m
00308           ubot(n) =   0.0_R8 ! atm velocity, zonal        ~ m/s
00309           vbot(n) =   2.0_R8 ! atm velocity, meridional   ~ m/s
00310           thbot(n)= 301.0_R8 ! atm potential temperature  ~ Kelvin
00311           shum(n) = 1.e-2_R8 ! atm specific humidity      ~ kg/kg
00312           dens(n) =   1.0_R8 ! atm density                ~ kg/m^3
00313           tbot(n) = 300.0_R8 ! atm temperature            ~ Kelvin
00314        enddo
00315     else        
00316        do n = 1,nloc
00317           if (mask(n) /= 0) then        
00318              zbot(n) = a2x%rAttr(index_a2x_Sa_z   ,n)
00319              ubot(n) = a2x%rAttr(index_a2x_Sa_u   ,n)
00320              vbot(n) = a2x%rAttr(index_a2x_Sa_v   ,n)
00321              thbot(n)= a2x%rAttr(index_a2x_Sa_ptem,n)
00322              shum(n) = a2x%rAttr(index_a2x_Sa_shum,n)
00323              dens(n) = a2x%rAttr(index_a2x_Sa_dens,n)
00324              tbot(n) = a2x%rAttr(index_a2x_Sa_tbot,n)
00325              tocn(n) = o2x%rAttr(index_o2x_So_t   ,n)   
00326              uocn(n) = o2x%rAttr(index_o2x_So_u   ,n)
00327              vocn(n) = o2x%rAttr(index_o2x_So_v   ,n)
00328           end if
00329        enddo
00330     end if
00331 
00332     call shr_flux_atmocn (nloc , zbot , ubot, vbot, thbot, &
00333                           shum , dens , tbot, uocn, vocn , &
00334                           tocn , mask , sen , lat , lwup , &
00335                           evap , taux , tauy, tref, qref , &
00336                           duu10n,ustar, re  , ssq )
00337 
00338     do n = 1,nloc
00339        if (mask(n) /= 0) then   
00340           xao%rAttr(index_xao_Faox_sen ,n) = sen(n)
00341           xao%rAttr(index_xao_Faox_lat ,n) = lat(n)
00342           xao%rAttr(index_xao_Faox_taux,n) = taux(n)
00343           xao%rAttr(index_xao_Faox_tauy,n) = tauy(n)
00344           xao%rAttr(index_xao_Faox_evap,n) = evap(n)
00345           xao%rAttr(index_xao_So_tref  ,n) = tref(n)
00346           xao%rAttr(index_xao_So_qref  ,n) = qref(n)
00347           xao%rAttr(index_xao_So_ustar ,n) = ustar(n)  ! friction velocity
00348           xao%rAttr(index_xao_So_re    ,n) = re(n)     ! reynolds number
00349           xao%rAttr(index_xao_So_ssq   ,n) = ssq(n)    ! s.hum. saturation at Ts
00350           xao%rAttr(index_xao_Faox_lwup,n) = lwup(n)   
00351           xao%rAttr(index_xao_Sx_duu10n,n) = duu10n(n)  
00352        end if
00353     enddo
00354 
00355   end subroutine seq_flux_atmocn_mct
00356 
00357 end module seq_flux_mct