module seq_flux_mct 1,9 use shr_kind_mod, only: r8 => shr_kind_r8, in=>shr_kind_in use shr_sys_mod, only: shr_sys_abort use shr_flux_mod, only: shr_flux_atmocn use shr_orb_mod, only: shr_orb_params, shr_orb_cosz, shr_orb_decl use mct_mod use seq_flds_indices use seq_comm_mct use seq_cdata_mod use seq_infodata_mod implicit none private save !-------------------------------------------------------------------------- ! Public interfaces !-------------------------------------------------------------------------- public seq_flux_init_mct public seq_flux_atmocn_mct !-------------------------------------------------------------------------- ! Private data !-------------------------------------------------------------------------- real(r8), pointer :: lats(:) ! latitudes (degrees) real(r8), pointer :: lons(:) ! longitudes (degrees) integer(in) , allocatable :: mask (:) ! ocn domain mask: 0 <=> inactive cell real(r8), allocatable :: uocn (:) ! ocn velocity, zonal real(r8), allocatable :: vocn (:) ! ocn velocity, meridional real(r8), allocatable :: tocn (:) ! ocean temperature real(r8), allocatable :: zbot (:) ! atm level height real(r8), allocatable :: ubot (:) ! atm velocity, zonal real(r8), allocatable :: vbot (:) ! atm velocity, meridional real(r8), allocatable :: thbot(:) ! atm potential T real(r8), allocatable :: shum (:) ! atm specific humidity real(r8), allocatable :: dens (:) ! atm density real(r8), allocatable :: tbot (:) ! atm bottom surface T real(r8), allocatable :: sen (:) ! heat flux: sensible real(r8), allocatable :: lat (:) ! heat flux: latent real(r8), allocatable :: lwup (:) ! lwup over ocean real(r8), allocatable :: evap (:) ! water flux: evaporation real(r8), allocatable :: taux (:) ! wind stress, zonal real(r8), allocatable :: tauy (:) ! wind stress, meridional real(r8), allocatable :: tref (:) ! diagnostic: 2m ref T real(r8), allocatable :: qref (:) ! diagnostic: 2m ref Q real(r8), allocatable :: duu10n(:) ! diagnostic: 10m wind speed squared real(r8), allocatable :: ustar(:) ! saved ustar real(r8), allocatable :: re (:) ! saved re real(r8), allocatable :: ssq (:) ! saved sq ! Conversion from degrees to radians real(r8),parameter :: const_pi = SHR_CONST_PI ! pi real(r8),parameter :: const_deg2rad = const_pi/180.0_R8 ! deg to rads !=============================================================================== contains !=============================================================================== subroutine seq_flux_init_mct( cdata, xao) 1,1 !----------------------------------------------------------------------- ! ! Arguments ! type(seq_cdata), intent(in) :: cdata type(mct_aVect), intent(out) :: xao ! ! Local variables ! integer(in) :: nloc type(mct_gGrid), pointer :: dom type(mct_gsMap), pointer :: gsMap integer :: lsize integer :: mpicom integer :: ier real(r8), pointer :: rmask(:) ! ocn domain mask character(*),parameter :: subName = '(seq_flux_init_mct) ' !----------------------------------------------------------------------- ! Set cdata pointers call seq_cdata_setptrs(cdata, gsMap=gsMap, dom=dom, mpicom=mpicom) lsize=mct_gsMap_lsize(gsMap, mpicom) ! Initialize attribute vector call mct_aVect_init(xao, rList=seq_flds_xao_fields, lsize=lsize) call mct_aVect_zero(xao) ! Determine size of buffers nloc = mct_avect_lSize(xao) ! Input fields atm allocate( zbot(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate zbot',ier) allocate( ubot(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate ubot',ier) allocate( vbot(nloc)) if(ier/=0) call mct_die(subName,'allocate vbot',ier) allocate(thbot(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate thbot',ier) allocate(shum(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate shum',ier) allocate(dens(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate dens',ier) allocate(tbot(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate tbot',ier) allocate(ustar(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate ustar',ier) allocate(re(nloc), stat=ier) if(ier/=0) call mct_die(subName,'allocate re',ier) allocate(ssq(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate ssq',ier) allocate( uocn(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate uocn',ier) allocate( vocn(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate vocn',ier) allocate( tocn(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate tocn',ier) ! Output fields allocate(sen (nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate sen',ier) allocate(lat (nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate lat',ier) allocate(evap(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate evap',ier) allocate(lwup(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate lwup',ier) allocate(taux(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate taux',ier) allocate(tauy(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate tauy',ier) allocate(tref(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate tref',ier) allocate(qref(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate qref',ier) allocate(duu10n(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate duu10n',ier) ! Grid fields allocate( lats(nloc),stat=ier ) if(ier/=0) call mct_die(subName,'allocate lats',ier) allocate( lons(nloc),stat=ier ) if(ier/=0) call mct_die(subName,'allocate lons',ier) allocate(mask(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate mask',ier) ! Get lat, lon, mask, which is time-invariant allocate(rmask(nloc),stat=ier) if(ier/=0) call mct_die(subName,'allocate rmask',ier) call mct_gGrid_exportRAttr(dom, 'lat' , lats , nloc) call mct_gGrid_exportRAttr(dom, 'lon' , lons , nloc) call mct_gGrid_exportRAttr(dom, 'mask', rmask, nloc) mask = nint(rmask) deallocate(rmask) end subroutine seq_flux_init_mct !=============================================================================== subroutine seq_flux_atmocn_mct( cdata, a2x, o2x, xao, fractions_o, albedo_only ) 4,8 !----------------------------------------------------------------------- ! ! Arguments ! type(seq_cdata), intent(in) :: cdata type(mct_aVect), intent(in) :: a2x type(mct_aVect), intent(in) :: o2x type(mct_aVect), intent(inout) :: xao type(mct_aVect), intent(inout) :: fractions_o logical, optional,intent(in) :: albedo_only ! ! Local variables ! type(seq_infodata_type),pointer :: infodata integer(in) :: n,i ! indices real(r8) :: rlat ! gridcell latitude in radians real(r8) :: rlon ! gridcell longitude in radians real(r8) :: cosz ! Cosine of solar zenith angle real(r8) :: eccen ! Earth orbit eccentricity real(r8) :: mvelpp ! Earth orbit real(r8) :: lambm0 ! Earth orbit real(r8) :: obliqr ! Earth orbit real(r8) :: delta ! Solar declination angle in radians real(r8) :: eccf ! Earth orbit eccentricity factor real(r8) :: calday ! calendar day including fraction, at 0e real(r8) :: nextsw_cday ! calendar day of next atm shortwave real(r8) :: anidr ! albedo: near infrared, direct real(r8) :: avsdr ! albedo: visible , direct real(r8) :: anidf ! albedo: near infrared, diffuse real(r8) :: avsdf ! albedo: visible , diffuse integer(in) :: nloc ! number of gridcells integer(in) :: ID ! comm ID logical :: flux_albav ! flux avg option logical :: dead_comps ! .true. => dead components are used integer :: kx,kr ! fractions indices ! real(R8),parameter :: albdif = 0.06_R8 ! 60 deg reference albedo, diffuse real(R8),parameter :: albdir = 0.07_R8 ! 60 deg reference albedo, direct character(*),parameter :: subName = '(seq_flux_atmocn_mct) ' ! !----------------------------------------------------------------------- call seq_cdata_setptrs(cdata, infodata=infodata, ID=ID) call seq_infodata_GetData(infodata,flux_albav=flux_albav) nloc = mct_aVect_lsize(xao) if (flux_albav) then do n=1,nloc anidr = albdir avsdr = albdir anidf = albdif avsdf = albdif ! Albedo is now function of latitude (will be new implementation) !rlat = const_deg2rad * lats(n) !anidr = 0.069_R8 - 0.011_R8 * cos(2._R8 * rlat) !avsdr = anidr !anidf = anidr !avsdf = anidr xao%rAttr(index_xao_So_avsdr,n) = avsdr xao%rAttr(index_xao_So_anidr,n) = anidr xao%rAttr(index_xao_So_avsdf,n) = avsdf xao%rAttr(index_xao_So_anidf,n) = anidf end do kx = mct_aVect_indexRA(fractions_o,"ifrac") kr = mct_aVect_indexRA(fractions_o,"ifrad") fractions_o%rAttr(kr,:) = fractions_o%rAttr(kx,:) kx = mct_aVect_indexRA(fractions_o,"ofrac") kr = mct_aVect_indexRA(fractions_o,"ofrad") fractions_o%rAttr(kr,:) = fractions_o%rAttr(kx,:) else ! Solar declination ! Will only do albedo calculation if nextsw_cday is not -1. call seq_infodata_GetData(infodata,nextsw_cday=nextsw_cday,orb_eccen=eccen, & orb_mvelpp=mvelpp, orb_lambm0=lambm0, orb_obliqr=obliqr) if (nextsw_cday >= -0.5_r8) then calday = nextsw_cday call shr_orb_decl(calday, eccen, mvelpp,lambm0, obliqr, delta, eccf) ! Compute albedos do n=1,nloc rlat = const_deg2rad * lats(n) rlon = const_deg2rad * lons(n) cosz = shr_orb_cosz( calday, rlat, rlon, delta ) if (cosz > 0.0_R8) then !--- sun hit -- anidr = (.026_R8/(cosz**1.7_R8 + 0.065_R8)) + & (.150_R8*(cosz - 0.100_R8 ) * & (cosz - 0.500_R8 ) * & (cosz - 1.000_R8 ) ) avsdr = anidr anidf = albdif avsdf = albdif else !--- dark side of earth --- anidr = 1.0_R8 avsdr = 1.0_R8 anidf = 1.0_R8 avsdf = 1.0_R8 end if xao%rAttr(index_xao_So_avsdr,n) = avsdr xao%rAttr(index_xao_So_anidr,n) = anidr xao%rAttr(index_xao_So_avsdf,n) = avsdf xao%rAttr(index_xao_So_anidf,n) = anidf end do ! nloc kx = mct_aVect_indexRA(fractions_o,"ifrac") kr = mct_aVect_indexRA(fractions_o,"ifrad") fractions_o%rAttr(kr,:) = fractions_o%rAttr(kx,:) kx = mct_aVect_indexRA(fractions_o,"ofrac") kr = mct_aVect_indexRA(fractions_o,"ofrad") fractions_o%rAttr(kr,:) = fractions_o%rAttr(kx,:) endif ! nextsw_cday end if ! flux_albav if (present(albedo_only)) then if (albedo_only) then if (seq_comm_iamroot(ID)) & write(logunit,'(2A)') trim(subname),' computing only ocn albedos ' RETURN endif endif !----------------------------------------------------- ! Update ocean surface fluxes ! Must fabricate "reasonable" data (using dead components) call seq_infodata_GetData(infodata, dead_comps=dead_comps) if (dead_comps) then do n = 1,nloc mask(n) = 1 ! ocn domain mask ~ 0 <=> inactive cell tocn(n) = 290.0_R8 ! ocn temperature ~ Kelvin uocn(n) = 0.0_R8 ! ocn velocity, zonal ~ m/s vocn(n) = 0.0_R8 ! ocn velocity, meridional ~ m/s zbot(n) = 55.0_R8 ! atm height of bottom layer ~ m ubot(n) = 0.0_R8 ! atm velocity, zonal ~ m/s vbot(n) = 2.0_R8 ! atm velocity, meridional ~ m/s thbot(n)= 301.0_R8 ! atm potential temperature ~ Kelvin shum(n) = 1.e-2_R8 ! atm specific humidity ~ kg/kg dens(n) = 1.0_R8 ! atm density ~ kg/m^3 tbot(n) = 300.0_R8 ! atm temperature ~ Kelvin enddo else do n = 1,nloc if (mask(n) /= 0) then zbot(n) = a2x%rAttr(index_a2x_Sa_z ,n) ubot(n) = a2x%rAttr(index_a2x_Sa_u ,n) vbot(n) = a2x%rAttr(index_a2x_Sa_v ,n) thbot(n)= a2x%rAttr(index_a2x_Sa_ptem,n) shum(n) = a2x%rAttr(index_a2x_Sa_shum,n) dens(n) = a2x%rAttr(index_a2x_Sa_dens,n) tbot(n) = a2x%rAttr(index_a2x_Sa_tbot,n) tocn(n) = o2x%rAttr(index_o2x_So_t ,n) uocn(n) = o2x%rAttr(index_o2x_So_u ,n) vocn(n) = o2x%rAttr(index_o2x_So_v ,n) end if enddo end if call shr_flux_atmocn (nloc , zbot , ubot, vbot, thbot, & shum , dens , tbot, uocn, vocn , & tocn , mask , sen , lat , lwup , & evap , taux , tauy, tref, qref , & duu10n,ustar, re , ssq ) do n = 1,nloc if (mask(n) /= 0) then xao%rAttr(index_xao_Faox_sen ,n) = sen(n) xao%rAttr(index_xao_Faox_lat ,n) = lat(n) xao%rAttr(index_xao_Faox_taux,n) = taux(n) xao%rAttr(index_xao_Faox_tauy,n) = tauy(n) xao%rAttr(index_xao_Faox_evap,n) = evap(n) xao%rAttr(index_xao_So_tref ,n) = tref(n) xao%rAttr(index_xao_So_qref ,n) = qref(n) xao%rAttr(index_xao_So_ustar ,n) = ustar(n) ! friction velocity xao%rAttr(index_xao_So_re ,n) = re(n) ! reynolds number xao%rAttr(index_xao_So_ssq ,n) = ssq(n) ! s.hum. saturation at Ts xao%rAttr(index_xao_Faox_lwup,n) = lwup(n) xao%rAttr(index_xao_Sx_duu10n,n) = duu10n(n) end if enddo end subroutine seq_flux_atmocn_mct end module seq_flux_mct