module shr_flux use var use const use comm use io !============================================================================ !This module calculate ocean needed boundary fluxes for ocean stand alone !run and coupled model framework(maybe?, we need to tackle with the mapping tho). !============================================================================ !By default, the wind stress are using scow or hellerman data; ! Salinity and Temperature are derived from PHC2 or PHC3; ! radiation data (swnet, lwdn, ice_fraction) are using NASA GISS and SSM/I ! datasets. !For other fluxes, we use CORE.v2 interannually varying forcing datasets. !The derivation method follows Large and Yeager (2008) and cesm1.0.6 framework. !This approach may cause some problem since we use different forcing datasets ! for wind stress and radiation. In the future, we will fully derive all ! needed fluxes with CORE.v2 datasets for consistancy. ! last edit: yuyulin, 30,11,2017 logical :: & ! lflux_ini = .true. !flag for initialization. lflux_ini = .false. ! CESM coupling logical :: & ! lflux_internal = .true. !using internal flux data from prepglo/data (core.v2 monthly data) lflux_internal = .false. !CESM coupling integer :: & lbulk_opt = 3 !option for deriving turbulent flux integer, parameter :: & lbulk_ncar = 1, &! ncar_ocean_fluxes_mod lbulk_cesmtimcom = 2, &! cesmtimcom_shr_flux_mod timcom_partially_coupled = 3 real*8,parameter :: avgalb = 0.066 real*8,parameter :: shr_const_g = 9.80616d0 real*8,parameter :: shr_const_karman = 0.4d0 real*8,parameter :: shr_const_boltz = 1.38065e-23 ! Boltzmann's constant ~ J/K/molecule real*8,parameter :: shr_const_avogad = 6.02214e26 ! Avogadro's number ~ molecules/kmole real*8,parameter :: shr_const_rgas = shr_const_avogad*shr_const_boltz! Universal gas constant ~ J/K/kmole real*8,parameter :: shr_const_mwdair = 28.966d0 ! molecular weight dry air ~ kg/kmole real*8,parameter :: shr_const_mwwv = 18.016d0 ! molecular weight water vapor real*8,parameter :: shr_const_rdair = shr_const_rgas/shr_const_mwdair ! Dry air gas constant ~ J/K/kg real*8,parameter :: shr_const_rwv = shr_const_rgas/shr_const_mwwv ! Water vapor gas constant ~ J/K/kg real*8,parameter :: shr_const_zvir =(shr_const_rwv/shr_const_rdair)-1.d0 ! RWV/RDAIR - 1.0 real*8,parameter :: shr_const_cpdair = 1.00464e3 ! specific heat of dry air ~ J/kg/K real*8,parameter :: shr_const_cpwv = 1.810e3 ! specific heat of water vap ~ J/kg/K real*8,parameter :: shr_const_cpvir =(shr_const_cpwv/shr_const_cpdair)-1.d0 ! CPWV/CPDAIR - 1.0 real*8,parameter :: shr_const_latvap = 2.501e6 ! latent heat of evaporation ~ J/kg real*8,parameter :: shr_const_latice = 3.337e5 ! latent heat of fusion ~ J/kg real*8,parameter :: shr_const_stebol = 5.67e-8 ! Stefan-Boltzmann constant ~ W/m^2/K^4 real*8, dimension(:,:,:), allocatable :: flux_scr contains subroutine flux_init_atmocn(invec) use forcing_fields, only: ROFF_F integer*2, dimension(:), intent(in) :: invec character*40 :: fname if (lflux_internal) then allocate(flux_scr(i2,j2,11)) !read run-off data, since the data is time-independent, !we dont have to do time-interpolation. fname='prepglo/data/runoff.nc' call rdnc_ij(fname,'runoff',roff) else ! init for coupled model roff = ROFF_F(3:nx_block-2,3:ny_block-2,1) end if call mtx2vec_int2_2d(in(:,:,1),invec) end subroutine subroutine flux_atmocn use forcing_fields, only: PREC_F, SNOW_F_ocn, LWDN_F, EVAP_F, SENH_F, & LWUP_F, IOFF_F_ocn, MELT_F, MELTH_F, SHF_QSW_ocn, & TAUXX_ocn,TAUYY_ocn, SALT_F, IFRAC_ocn use grid, only: RCALCT use blocks, only: nx_block, ny_block integer :: i,j,k integer*2, dimension(i2,j2) :: msk integer*2, dimension(i2*j2) :: invec real*8, dimension(i2,j2) :: temp,ifrac real*8, dimension(i2*j2,6) :: work2 real*8, dimension(i2*j2) :: work,zatm real*8, dimension(i2*j2) :: rhovec_a,pvec_a real*8, dimension(i2*j2) :: qvec_o,qvec_a real*8, dimension(i2*j2) :: tvec_o,tvec_a real*8, dimension(i2*j2) :: svec_o real*8, dimension(i2*j2) :: uvec_o,uvec_a real*8, dimension(i2*j2) :: vvec_o,vvec_a real*8, dimension(i2*j2) :: & cd,ch,ce,ustar,bstar if (lflux_ini) then call flux_init_atmocn(invec) lflux_ini = .false. else call flux_init_atmocn(invec) end if if (lbulk_opt /= timcom_partially_coupled) then call flux_atmocn_getdata(days) !update variables and convert units for flux calculation. !ocean temp, velocity. call mtx2vec_dbl_2d(t2(:,:,1)+273.15d0,tvec_o,1) call mtx2vec_dbl_2d(u2(:,:,1)/100.d0,uvec_o,1) call mtx2vec_dbl_2d(v2(:,:,1)/100.d0,vvec_o,1) endif !atmosphere variables read from nc-file !atm 10m temp, velocity, humidity and precip, rediation, runoff, slp.. if (lflux_internal) then !slp(Pa) call mtx2vec_dbl_2d(flux_scr(:,:,6),pvec_a,0) !q_10(kg/kg) call mtx2vec_dbl_2d(flux_scr(:,:,7),qvec_a,0) !t_10(K) call mtx2vec_dbl_2d(flux_scr(:,:,8),tvec_a,0) !u_10(m/s) call mtx2vec_dbl_2d(flux_scr(:,:,9),uvec_a,0) !v_10(m/s) call mtx2vec_dbl_2d(flux_scr(:,:,10),vvec_a,0) !ice fraction ifrac = flux_scr(:,:,11) msk = in(2:i1,2:j1,1) !prec = (1.d0-ifrac)*msk*(fold*flux_scr(:,:,1)+fnew*flux_scr(:,:,1)) !rain = (1.d0-ifrac)*msk*(fold*flux_scr(:,:,2)+fnew*flux_scr(:,:,2)) !snow = (1.d0-ifrac)*msk*(fold*flux_scr(:,:,3)+fnew*flux_scr(:,:,3)) !lwdn = (1.d0-ifrac)*msk*(fold*flux_scr(:,:,4)+fnew*flux_scr(:,:,4)) !qdot2 = (1.d0-ifrac)*msk*(fold*flux_scr(:,:,5)+fnew*flux_scr(:,:,5)) prec = (1.d0-ifrac)*msk*flux_scr(:,:,1) rain = (1.d0-ifrac)*msk*flux_scr(:,:,2) snow = (1.d0-ifrac)*msk*flux_scr(:,:,3) lwdn = (1.d0-ifrac)*msk*flux_scr(:,:,4) qdot2 = (1.d0-ifrac)*msk*flux_scr(:,:,5) snow_f = -shr_const_latice*snow else !-------------------------- ! 20190924 cpl_mct !-------------------------- prec = PREC_F(3:nx_block-2,3:ny_block-2,1) ! kg/(m^2)/s rain = PREC_F(3:nx_block-2,3:ny_block-2,1) & - SNOW_F_ocn(3:nx_block-2,3:ny_block-2,1) ! kg/(m^2)/s snow = SNOW_F_ocn(3:nx_block-2,3:ny_block-2,1) ! kg/(m^2)/s lwdn = LWDN_F(3:nx_block-2,3:ny_block-2,1) ! W/m^2 qdot2 = SHF_QSW_ocn(3:nx_block-2,3:ny_block-2,1)/hflux_factor ! W/m^2 snow_f = -shr_const_latice*snow ! W/m^2 ! if (myid==0) write(*,*) 'ERROR, currently external sources are not available' ! call mpi_finalize(ierr) ! stop 0 end if !=========================================================== select case (lbulk_opt) !=========================================================== !----------------------------------------------------------- case (lbulk_ncar) !----------------------------------------------------------- do i=1,i2*j2 !p10m work(i) = patm_10m(tvec_a(i),qvec_a(i),pvec_a(i)) !potential temp_10m tvec_a(i) = potemp(work(i),tvec_a(i)) !atm_rho rhovec_a(i) = rho_wetatm(pvec_a(i),tvec_a(i),qvec_a(i)) !ocn_qsat qvec_o(i) = qsat_ocn(rhovec_a(i),tvec_o(i)) !shear velocity work(i) = max(1.d0,sqrt ((uvec_a(i)-uvec_o(i))**2 + & (vvec_a(i)-vvec_o(i))**2) ) end do if (lflux_internal) zatm=10.d0 call ncar_ocean_fluxes(work,tvec_a,tvec_o,qvec_a,qvec_o, & zatm,cd,ch,ce,invec,ustar,bstar) !evaporative mass flux work2(:,1) = invec(:)*rhovec_a(:)*ce(:)*(qvec_a(:)-qvec_o(:))*work(:) !latent heat flux (evaporative heat flux) work2(:,2) = shr_const_latvap*work2(:,1) !sensible heat flux work2(:,3) = invec(:)*rhovec_a(:)*shr_const_cpdair*ch(:)*(tvec_a(:)-tvec_o(:))*work(:) !ocean taux work2(:,4) = invec(:)*rhovec_a(:)*cd(:)*work(:)*(uvec_a(:)-uvec_o(:)) !ocean tauy work2(:,5) = invec(:)*rhovec_a(:)*cd(:)*work(:)*(vvec_a(:)-vvec_o(:)) !lwup work2(:,6) = -shr_const_stebol*invec*tvec_o**4 call vec2mtx_dbl_2d(evapo, work2(:,1),0) call vec2mtx_dbl_2d(lath, work2(:,2),0) call vec2mtx_dbl_2d(senh, work2(:,3),0) call vec2mtx_dbl_2d(tauxx, work2(:,4),0) call vec2mtx_dbl_2d(tauyy, work2(:,5),0) call vec2mtx_dbl_2d(lwup, work2(:,6),0) !compute air-ice wind stress and add into ocn-surface wind stress work(:) = sqrt( uvec_a(:)**2 + vvec_a(:)**2 ) work2(:,1) = invec(:)*rhovec_a(:)*1.63e-3*work(:)*uvec_a(:) work2(:,2) = invec(:)*rhovec_a(:)*1.63e-3*work(:)*vvec_a(:) !apply ice fraction mask call vec2mtx_dbl_2d(temp, work2(:,1), 0) tauxx = (1.d0-ifrac)*tauxx + ifrac*temp call vec2mtx_dbl_2d(temp, work2(:,2), 0) tauyy = (1.d0-ifrac)*tauyy + ifrac*temp !convert units to cgs system, and apply ice fraction mask tauxx=10.d0*tauxx tauyy=10.d0*tauyy senh = (1.d0-ifrac)*senh lath = (1.d0-ifrac)*lath lwup = (1.d0-ifrac)*lwup !qdot, total heat flux w/o short wave radiation qdot = senh + lath + lwup + lwdn + snow_f !----------------------------------------------------------------- case (lbulk_cesmtimcom) !----------------------------------------------------------------- do i=1,i2*j2 !p10m work(i) = patm_10m(tvec_a(i),qvec_a(i),pvec_a(i)) !potential temp_10m tvec_a(i) = tvec_a(i)!potemp(work(i),tvec_a(i)) !atm_rho rhovec_a(i) = rho_wetatm(pvec_a(i),tvec_a(i),qvec_a(i)) !ocn_qsat qvec_o(i) = qsat_ocn(rhovec_a(i),tvec_o(i)) end do if (lflux_internal) zatm=10.d0 !testing values; call cesm_shr_flux_atmOcn (invec,zatm,uvec_a,vvec_a,tvec_a,qvec_a,rhovec_a, & uvec_o,vvec_o,tvec_o,work2) call vec2mtx_dbl_2d(tauxx, work2(:,1),0) call vec2mtx_dbl_2d(tauyy, work2(:,2),0) call vec2mtx_dbl_2d(senh, work2(:,3),0) call vec2mtx_dbl_2d(lath, work2(:,4),0) call vec2mtx_dbl_2d(lwup, work2(:,5),0) call vec2mtx_dbl_2d(evapo, work2(:,6),0) !compute air-ice wind stress and add into ocn-surface wind stress work(:) = sqrt( uvec_a(:)**2 + vvec_a(:)**2 ) work2(:,1) = invec(:)*rhovec_a(:)*1.63e-3*work(:)*uvec_a(:) work2(:,2) = invec(:)*rhovec_a(:)*1.63e-3*work(:)*vvec_a(:) !apply ice fraction mask call vec2mtx_dbl_2d(temp, work2(:,1), 0) tauxx = (1.d0-ifrac)*tauxx + ifrac*temp call vec2mtx_dbl_2d(temp, work2(:,2), 0) tauyy = (1.d0-ifrac)*tauyy + ifrac*temp !convert units to cgs system, and apply ice fraction mask tauxx=10.d0*tauxx tauyy=10.d0*tauyy evapo= (1.d0-ifrac)*evapo senh = (1.d0-ifrac)*senh lath = (1.d0-ifrac)*lath lwup = (1.d0-ifrac)*lwup !qdot, total heat flux w/o short wave radiation qdot = senh + lath + lwup + lwdn + snow_f !----------------------------------------------------------------- case (timcom_partially_coupled) !----------------------------------------------------------------- !-------------------------- ! 20190924 cpl_mct !-------------------------- tauxx=TAUXX_ocn(3:nx_block-2,3:ny_block-2,1)& *RCALCT(3:nx_block-2,3:ny_block-2,1)*10.d0 !velocity flux (cm^2/s^2) tauyy=TAUYY_ocn(3:nx_block-2,3:ny_block-2,1)& *RCALCT(3:nx_block-2,3:ny_block-2,1)*10.d0 !velocity flux (cm^2/s^2) evapo= EVAP_F(3:nx_block-2,3:ny_block-2,1) ! kg/(m^2)/s senh = SENH_F(3:nx_block-2,3:ny_block-2,1) ! W/m^2 lath = EVAP_F(3:nx_block-2,3:ny_block-2,1)*2.5e6 ! kg/(m^2)/s * (latent_heat_vapor_mks)J/kg = W/m^2 lwup = LWUP_F(3:nx_block-2,3:ny_block-2,1) ! W/m^2 qdot = (EVAP_F(3:nx_block-2,3:ny_block-2,1)*2.5e6 & + SENH_F(3:nx_block-2,3:ny_block-2,1) & + LWUP_F(3:nx_block-2,3:ny_block-2,1) & + LWDN_F(3:nx_block-2,3:ny_block-2,1) & + MELTH_F(3:nx_block-2,3:ny_block-2,1) & - (SNOW_F_ocn(3:nx_block-2,3:ny_block-2,1) & + IOFF_F_ocn(3:nx_block-2,3:ny_block-2,1)) & *shr_const_latice)*RCALCT(3:nx_block-2,3:ny_block-2,1) ! W/m^2 melt = MELT_F(3:nx_block-2,3:ny_block-2,1) ! kg/(m^2)/s melth= MELTH_F(3:nx_block-2,3:ny_block-2,1) ! W/m^2 ioff = IOFF_F_ocn(3:nx_block-2,3:ny_block-2,1) ! kg/(m^2)/s ioff_f = -shr_const_latice*ioff ! W/m^2 salt = SALT_F(3:nx_block-2,3:ny_block-2,1) ! kg(salt)/(m^2)/s ifrac = IFRAC_ocn(3:nx_block-2,3:ny_block-2,1) ! ATM_PRESS(3:nx_block-2,3:ny_block-2,1) ! dynes/cm^2 ! U10_SQR(3:nx_block-2,3:ny_block-2,1) ! cm^2/s^2 end select end subroutine flux_atmocn subroutine cesm_shr_flux_atmOcn (mask,zbot,ubot,vbot,thbot,qbot,rbot, & us,vs,ts,worktemp) !--- input arguments -------------------------------- integer*2,intent(in) :: mask (:) ! ocn domain mask 0 <=> out of domain real*8 ,intent(in) :: zbot (:) ! atm level height (m) real*8 ,intent(in) :: ubot (:) ! atm u wind (m/s) real*8 ,intent(in) :: vbot (:) ! atm v wind (m/s) real*8 ,intent(in) :: thbot(:) ! atm potential T (K) real*8 ,intent(in) :: qbot (:) ! atm specific humidity (kg/kg) real*8 ,intent(in) :: rbot (:) ! atm air density (kg/m^3) real*8 ,intent(in) :: us (:) ! ocn u-velocity (m/s) real*8 ,intent(in) :: vs (:) ! ocn v-velocity (m/s) real*8 ,intent(in) :: ts (:) ! ocn temperature (K) !--- output arguments ------------------------------- real*8, intent(inout),dimension(:,:) :: worktemp !storing 1d array for evap,lat,sen,taux,tauy,lwup ! real*8, intent(inout) :: sen (nMax) ! heat flux: sensible (W/m^2) ! real*8, intent(inout) :: lat (nMax) ! heat flux: latent (W/m^2) ! real*8, intent(inout) :: lwup (nMax) ! heat flux: lw upward (W/m^2) ! real*8, intent(inout) :: evap (nMax) ! water flux: evap ((kg/s)/m^2) ! real*8, intent(inout) :: taux (nMax) ! surface stress, zonal (N) ! real*8, intent(inout) :: tauy (nMax) ! surface stress, maridional (N) ! real*8,intent(out) :: tref (nMax) ! diag: 2m ref height T (K) ! real*8,intent(out) :: qref (nMax) ! diag: 2m ref humidity (kg/kg) ! real*8,intent(out) :: duu10n(nMax) ! diag: 10m wind speed squared (m/s)^2 ! real*8,intent(out),optional :: ustar_sv(nMax) ! diag: ustar ! real*8,intent(out),optional :: re_sv (nMax) ! diag: sqrt of exchange coefficient (water) ! real*8,intent(out),optional :: ssq_sv (nMax) ! diag: sea surface humidity (kg/kg) ! !EOP !--- local constants -------------------------------- real*8,parameter :: umin = 0.5d0 ! minimum wind speed (m/s) real*8,parameter :: zref = 10.d0 ! reference height (m) real*8,parameter :: ztref = 2.d0 ! reference height for air T (m) !--- local variables -------------------------------- real*8 :: vmag ! surface wind magnitude (m/s) real*8 :: thvbot ! virtual temperature (K) real*8 :: ssq ! sea surface humidity (kg/kg) real*8 :: delt ! potential T difference (K) real*8 :: delq ! humidity difference (kg/kg) real*8 :: stable ! stability factor real*8 :: rdn ! sqrt of neutral exchange coeff (momentum) real*8 :: rhn ! sqrt of neutral exchange coeff (heat) real*8 :: ren ! sqrt of neutral exchange coeff (water) real*8 :: rd ! sqrt of exchange coefficient (momentum) real*8 :: rh ! sqrt of exchange coefficient (heat) real*8 :: re ! sqrt of exchange coefficient (water) real*8 :: ustar ! ustar real*8 :: qstar ! qstar real*8 :: tstar ! tstar real*8 :: hol ! H (at zbot) over L real*8 :: xsq ! ? real*8 :: xqq ! ? real*8 :: psimh ! stability function at zbot (momentum) real*8 :: psixh ! stability function at zbot (heat and water) real*8 :: psix2 ! stability function at ztref reference height real*8 :: alz ! ln(zbot/zref) real*8 :: al2 ! ln(zref/ztref) real*8 :: u10n ! 10m neutral wind real*8 :: tau ! stress at zbot real*8 :: cp ! specific heat of moist air real*8 :: bn ! exchange coef funct for interpolation real*8 :: bh ! exchange coef funct for interpolation real*8 :: fac ! vertical interpolation factor real*8 :: spval ! local missing value integer :: nMax integer :: n !--- local functions -------------------------------- real*8 :: qsat ! function: the saturation humididty of air (kg/m^3) real*8 :: cdn ! function: neutral drag coeff at 10m real*8 :: psimhu ! function: unstable part of psimh real*8 :: psixhu ! function: unstable part of psimx real*8 :: Umps ! dummy arg ~ wind velocity (m/s) real*8 :: Tk ! dummy arg ~ temperature (K) real*8 :: xd ! dummy arg ~ ? qsat(Tk) = 640380.d0 / exp(5107.4d0/Tk) cdn(Umps) = 0.0027d0 / Umps + 0.000142d0 + 0.0000764d0 * Umps psimhu(xd) = log((1.d0+xd*(2.d0+xd))*(1.d0+xd*xd)/8.d0) - 2.d0*atan(xd) + 1.571d0 psixhu(xd) = 2.d0 * log((1.d0 + xd*xd)/2.d0) !------------------------------------------------------------------------------- ! PURPOSE: ! computes atm/ocn surface fluxes ! ! NOTES: ! o all fluxes are positive downward ! o net heat flux = net sw + lw up + lw down + sen + lat ! o here, tstar = /U*, and qstar = /U*. ! o wind speeds should all be above a minimum speed (eg. 1.0 m/s) ! ! ASSUMPTIONS: ! o Neutral 10m drag coeff: cdn = .0027/U10 + .000142 + .0000764 U10 ! o Neutral 10m stanton number: ctn = .0327 sqrt(cdn), unstable ! ctn = .0180 sqrt(cdn), stable ! o Neutral 10m dalton number: cen = .0346 sqrt(cdn) ! o The saturation humidity of air at T(K): qsat(T) (kg/m^3) !------------------------------------------------------------------------------- nMax= size(mask) al2 = log(zref/ztref) DO n=1,nMax if (mask(n) /= 0) then !--- compute some needed quantities --- vmag = max(umin, sqrt( (ubot(n)-us(n))**2 + (vbot(n)-vs(n))**2) ) thvbot = thbot(n) * (1.d0 + shr_const_zvir * qbot(n)) ! virtual temp (K) ssq = 0.98d0 * qsat(ts(n)) / rbot(n) ! sea surf hum (kg/kg) delt = thbot(n) - ts(n) ! pot temp diff (K) delq = qbot(n) - ssq ! spec hum dif (kg/kg) alz = log(zbot(n)/zref) cp = shr_const_cpdair*(1.d0 + shr_const_cpvir*ssq) !------------------------------------------------------------ ! first estimate of Z/L and ustar, tstar and qstar !------------------------------------------------------------ !--- neutral coefficients, z/L = 0.0 --- stable = 0.5d0 + sign(0.5d0 , delt) rdn = sqrt(cdn(vmag)) rhn = (1.d0-stable) * 0.0327d0 + stable * 0.018d0 ren = 0.0346d0 !--- ustar, tstar, qstar --- ustar = rdn * vmag tstar = rhn * delt qstar = ren * delq !--- compute stability & evaluate all stability functions --- hol = shr_const_karman*shr_const_g*zbot(n)* & (tstar/thbot(n)+qstar/(1.d0/shr_const_zvir+qbot(n)))/ustar**2 hol = sign( min(abs(hol),10.d0), hol ) stable = .5d0 + sign(.5d0 , hol) xsq = max(sqrt(abs(1.d0 - 16.d0*hol)) , 1.d0) xqq = sqrt(xsq) psimh = -5.d0*hol*stable + (1.d0-stable)*psimhu(xqq) psixh = -5.d0*hol*stable + (1.d0-stable)*psixhu(xqq) !--- shift wind speed using old coefficient --- rd = rdn / (1.d0 + rdn/shr_const_karman*(alz-psimh)) u10n = vmag * rd / rdn !--- update transfer coeffs at 10m and neutral stability --- rdn = sqrt(cdn(u10n)) ren = 0.0346d0 rhn = (1.d0-stable)*0.0327d0 + stable * 0.018d0 !--- shift all coeffs to measurement height and stability --- rd = rdn / (1.d0 + rdn/shr_const_karman*(alz-psimh)) rh = rhn / (1.d0 + rhn/shr_const_karman*(alz-psixh)) re = ren / (1.d0 + ren/shr_const_karman*(alz-psixh)) !--- update ustar, tstar, qstar using updated, shifted coeffs -- ustar = rd * vmag tstar = rh * delt qstar = re * delq !------------------------------------------------------------ ! iterate to converge on Z/L, ustar, tstar and qstar !------------------------------------------------------------ !--- compute stability & evaluate all stability functions --- hol = shr_const_karman*shr_const_g*zbot(n)* & (tstar/thbot(n)+qstar/(1.d0/shr_const_zvir+qbot(n)))/ustar**2 hol = sign( min(abs(hol),10.d0), hol ) stable = .5d0 + sign(.5d0 , hol) xsq = max(sqrt(abs(1.d0 - 16.d0*hol)) , 1.d0) xqq = sqrt(xsq) psimh = -5.d0*hol*stable + (1.d0-stable)*psimhu(xqq) psixh = -5.d0*hol*stable + (1.d0-stable)*psixhu(xqq) !--- shift wind speed using old coeffs --- rd = rdn / (1.d0 + rdn/shr_const_karman*(alz-psimh)) u10n = vmag * rd/rdn !--- update transfer coeffs at 10m and neutral stability --- rdn = sqrt(cdn(u10n)) ren = 0.0346d0 rhn = (1.d0 - stable)*0.0327d0 + stable * 0.018d0 !--- shift all coeffs to measurement height and stability --- rd = rdn / (1.d0 + rdn/shr_const_karman*(alz-psimh)) rh = rhn / (1.d0 + rhn/shr_const_karman*(alz-psixh)) re = ren / (1.d0 + ren/shr_const_karman*(alz-psixh)) !--- update ustar, tstar, qstar using updated, shifted coeffs --- ustar = rd * vmag tstar = rh * delt qstar = re * delq !------------------------------------------------------------ ! compute the fluxes !------------------------------------------------------------ tau = rbot(n) * ustar * ustar !--- momentum flux --- !taux worktemp(n,1) = tau * (ubot(n)-us(n)) / vmag !tauy worktemp(n,2) = tau * (vbot(n)-vs(n)) / vmag !--- heat flux --- !sensible heat worktemp(n,3) = cp * tau * tstar / ustar !latent heat worktemp(n,4) = shr_const_latvap * tau * qstar / ustar !long wave radiation worktemp(n,5) = -shr_const_stebol * ts(n)**4 !--- water flux --- !evaporation worktemp(n,6) = worktemp(n,4)/shr_const_latvap !------------------------------------------------------------ ! compute diagnositcs: 2m ref T & Q, 10m wind speed squared !------------------------------------------------------------ !hol = hol*ztref/zbot(n) !xsq = max( 1.0_R8, sqrt(abs(1.0_R8-16.0_R8*hol)) ) !xqq = sqrt(xsq) !psix2 = -5.0_R8*hol*stable + (1.0_R8-stable)*psixhu(xqq) !fac = (rh/shr_const_karman) * (alz + al2 - psixh + psix2 ) !tref(n) = thbot(n) - delt*fac !tref(n) = tref(n) - 0.01_R8*ztref ! pot temp to temp correction !fac = (re/shr_const_karman) * (alz + al2 - psixh + psix2 ) !qref(n) = qbot(n) - delq*fac !duu10n(n) = u10n*u10n ! 10m wind speed squared !------------------------------------------------------------ ! optional diagnostics, needed for water tracer fluxes (dcn) !------------------------------------------------------------ !if (present(ustar_sv)) ustar_sv(n) = ustar !if (present(re_sv )) re_sv(n) = re !if (present(ssq_sv )) ssq_sv(n) = ssq else !------------------------------------------------------------ ! no valid data here -- out of domain !------------------------------------------------------------ !sen (n) = spval ! sensible heat flux (W/m^2) !lat (n) = spval ! latent heat flux (W/m^2) !lwup (n) = spval ! long-wave upward heat flux (W/m^2) !evap (n) = spval ! evaporative water flux ((kg/s)/m^2) !taux (n) = spval ! x surface stress (N) !tauy (n) = spval ! y surface stress (N) !tref (n) = spval ! 2m reference height temperature (K) !qref (n) = spval ! 2m reference height humidity (kg/kg) !duu10n(n) = spval ! 10m wind speed squared (m/s)^2 !if (present(ustar_sv)) ustar_sv(n) = spval !if (present(re_sv )) re_sv (n) = spval !if (present(ssq_sv )) ssq_sv (n) = spval end if end do end subroutine cesm_shr_flux_atmOcn subroutine ncar_ocean_fluxes(uatm, tatm, ts, qatm, qs, zatm, cd, ch, ce, & inmask, ustar, bstar) ! If 10m height atm wind, temperature and humidity are availible, then ! the neutral 10m coefficients can be used directly; in that case, set zatm==10. ! to skip the Monin-Obukhov iteration. integer*2, dimension(:), intent(in) :: inmask real*8, dimension(:), intent(in) :: uatm, tatm, ts, qatm, qs, zatm real*8, dimension(:), intent(inout) :: cd, ch, ce, ustar, bstar real*8 :: cd_n10, ce_n10, ch_n10, cd_n10_rt ! neutral 10m drag coefficients real*8 :: cd_rt ! full drag coefficients @ z real*8 :: zeta, x2, x, psi_m, psi_h ! stability parameters real*8 :: u_zatm, u10, tv, tstar, qstar, z0, xx, stab integer, parameter :: n_itts = 2 integer i, j do i=1,size(uatm) if (inmask(i) /= 0) then tv = tatm(i)*(1.d0+.608d0*qatm(i)); u10 = uatm(i); ! first guess 10m wind cd_n10 = (2.7d0/u10+.142d0+.0764d0*u10)/1.e3; ! L-Y eqn. 6a cd_n10_rt = sqrt(cd_n10); ce_n10 = 34.6d0 *cd_n10_rt/1e3; ! L-Y eqn. 6b stab = .5d0 + sign(.5d0,tatm(i)-ts(i)) ch_n10 = (18.d0*stab+32.7d0*(1-stab))*cd_n10_rt/1e3; ! L-Y eqn. 6c cd(i) = cd_n10; ! first guess for exchange coeff's at ch(i) = ch_n10; ce(i) = ce_n10; if (zatm(i)==10.) go to 99 !skip iteration do j=1,n_itts ! Monin-Obukhov iteration cd_rt = sqrt(cd(i)); ustar(i) = cd_rt*u_zatm; ! L-Y eqn. 7a tstar = (ch(i)/cd_rt)*(tatm(i)-ts(i)); ! L-Y eqn. 7b qstar = (ce(i)/cd_rt)*(qatm(i)-qs(i)); ! L-Y eqn. 7c bstar(i) = shr_const_g*(tstar/tv+qstar/(qatm(i)+1.d0/.608d0)); zeta = shr_const_karman*bstar(i)*zatm(i)/(ustar(i)*ustar(i)); ! L-Y eqn. 8a zeta = sign( min(abs(zeta),10.d0), zeta ); ! undocumented NCAR x2 = sqrt(abs(1.d0-16.d0*zeta)); ! L-Y eqn. 8b x2 = max(x2, 1.d0); ! undocumented NCAR x = sqrt(x2); if (zeta > 0) then psi_m = -5*zeta; ! L-Y eqn. 8c psi_h = -5*zeta; ! L-Y eqn. 8c else psi_m = log((1+2*x+x2)*(1+x2)/8)-2*(atan(x)-atan(1.0)); ! L-Y eqn. 8d psi_h = 2*log((1+x2)/2); ! L-Y eqn. 8e end if u10 = u_zatm/(1+cd_n10_rt*(log(zatm(i)/10)-psi_m)/shr_const_karman); ! L-Y eqn. 9 cd_n10 = (2.7/u10+0.142+0.0764*u10)/1e3; ! L-Y eqn. 6a again cd_n10_rt = sqrt(cd_n10); ce_n10 = 34.6*cd_n10_rt/1e3; ! L-Y eqn. 6b again stab = 0.5 + sign(0.5,zeta) ch_n10 = (18.0*stab+32.7*(1-stab))*cd_n10_rt/1e3; ! L-Y eqn. 6c again z0 = 10*exp(-shr_const_karman/cd_n10_rt); ! diagnostic xx = (log(zatm(i)/10)-psi_m)/shr_const_karman; cd(i) = cd_n10/(1+cd_n10_rt*xx)**2; ! L-Y 10a xx = (log(zatm(i)/10)-psi_h)/shr_const_karman; ch(i) = ch_n10/(1+ch_n10*xx/cd_n10_rt)*sqrt(cd(i)/cd_n10) ! 10b (corrected code aug2007) ce(i) = ce_n10/(1+ce_n10*xx/cd_n10_rt)*sqrt(cd(i)/cd_n10) ! 10c (corrected code aug2007) end do 99 continue endif !ocn-mask end do end subroutine ncar_ocean_fluxes function qsat_ocn(rhoatm,tempatm) ! emperical formula for surface saturated spcific humidity. ! rho = air density (kg/m3) ! qsaturate = ocean surface saturated spcific humidity (kg/kg) real*8, intent(in) :: rhoatm,tempatm real*8 :: qsat_ocn qsat_ocn = 0.98d0/rhoatm*640380.d0*(exp(-5107.4d0/tempatm)) end function function rho_wetatm(patm,tempatm,qatm) ! density of wet air is calculated using the ideal gas law, P=rho*RT ! rho = wet air density (kg/m3) ! P = absolute pressure (Pa) ! T = absolute temperature (C) ! The specific gas constant for dry air is 287.04 J/(kg*K) in SI units real*8, intent(in) :: patm,tempatm,qatm real*8 :: rho_wetatm rho_wetatm = patm/( (1+0.608d0*qatm)*287.04d0*tempatm ) end function function patm_10m(tatm,qatm,slp) real*8, intent(in) :: tatm,qatm,slp real*8 :: patm_10m patm_10m=slp*exp(-(10.d0)/(29.3d0*tatm*(1.d0+.608d0*qatm))); end function function potemp(p10,tatm) real*8, intent(in) :: p10,tatm real*8 :: potemp potemp = tatm*(100000.d0/p10)**(2.d0/7.d0) end function subroutine flux_atmocn_getdata(days) character*80 :: fname integer :: dnum,mnum,mnumUB real*8, dimension(i2,j2,2) :: interptmp real*8, intent(in) :: days dnum = mod(days,365.d0) if (dnum ==0) dnum = 365 mnum = nld mnumUB = mnum + 1 if (myid==0) write(*,*) 'readata forcing at',days,'dayin',dnum,mnum !read precipitation data fname='prepglo/data/precip.nc' if (nld==12) mnumUB=1 call rdnc_ijm(fname,'prec',mnum,1,interptmp(:,:,1)) call rdnc_ijm(fname,'prec',mnumUB,1,interptmp(:,:,2)) flux_scr(:,:,1) = fold*interptmp(:,:,1)+fnew*interptmp(:,:,2) call rdnc_ijm(fname,'rain',mnum,1,interptmp(:,:,1)) call rdnc_ijm(fname,'rain',mnumUB,1,interptmp(:,:,2)) flux_scr(:,:,2) = fold*interptmp(:,:,1)+fnew*interptmp(:,:,2) call rdnc_ijm(fname,'snow',mnum,1,interptmp(:,:,1)) call rdnc_ijm(fname,'snow',mnumUB,1,interptmp(:,:,2)) flux_scr(:,:,3) = fold*interptmp(:,:,1)+fnew*interptmp(:,:,2) !read ice_fraction data fname='prepglo/data/ice.nc' call rdnc_ijm(fname,'ifrac',dnum,1,flux_scr(:,:,11)) !read radiation data fname='prepglo/data/radiation.nc' call rdnc_ijm(fname,'lwdn',dnum,1,flux_scr(:,:,4)) call rdnc_ijm(fname,'swdn',dnum,1,flux_scr(:,:,5)) flux_scr(:,:,5) = (1.d0-avgalb)*flux_scr(:,:,5) !read atm data for later turbulent fluxes calc. fname='prepglo/data/atm_field.nc' call rdnc_ijm(fname,'slp',dnum,1,flux_scr(:,:,6)) call rdnc_ijm(fname,'q_10',dnum,1,flux_scr(:,:,7)) call rdnc_ijm(fname,'t_10',dnum,1,flux_scr(:,:,8)) call rdnc_ijm(fname,'u_10',dnum,1,flux_scr(:,:,9)) call rdnc_ijm(fname,'v_10',dnum,1,flux_scr(:,:,10)) end subroutine end module