diff --git a/diag/pbldia_new_sfx.f90 b/diag/pbldia_new_sfx.f90 index 83f1b4000fbd79e3485cad36144ff8dd10b5f854..3ac1d887f040723ffd9dcf0a2a449082594e9d86 100755 --- a/diag/pbldia_new_sfx.f90 +++ b/diag/pbldia_new_sfx.f90 @@ -2,11 +2,11 @@ module pbldia_new_sfx implicit none -real,parameter,private::undef=-999.0,R=287.,appa=0.287 +real,parameter,private::kappa=0.4,R=287.,appa=0.287 contains -! in this subroutine interpolation is performed using integration of universal functions from z1 to z2, +! in this module interpolation is performed using integration of universal functions from z1 to z2, ! where z1 is measurement or lowest model level height and z2 is required height ! The input-output structure is preserved as much as possible as in the climate model (pbldia subroutine in pblflx) !C* DIAGNOSIS OF WIND VELOCITY (AT HEIGHT Z = HWIND), @@ -16,8 +16,6 @@ contains !C* ARRAY AR2 (OUTPUT FROM DRAG3 - SUBROUTINE) !C* AR2(1) - NON-DIMENSIONAL (NORMALIZED BY MONIN-OBUKHOV LENGTH) !C* HEIGHT OF CONSTANT FLUX LAYER TOP -!C* AR2(5) - ROUGHNESS LENGTH FOR MOMENTUM -!C* AR2(6) - ROUGHNESS LENGTH FOR TEMPERATURE AND HUMIDITY !C* AR2(8) - INTEGRAL TRANSFER COEFFICIENT FOR MOMENTUM !C* AR2(9) - INTEGRAL TRANSFER COEFFICIENT FOR TEMPERATURE AND HUMIDITY !C* ARRAY ARDIN @@ -29,32 +27,19 @@ contains !C* ARDIN(6) - = HTEMP !C* OUTPUT DATA (ARRAY ARDOUT): !C* ARDOUT(1) - MODULE OF WIND VELOCITY AT REQUIRED HEIGHT -!C* ARDOUT(2) - POTENTIAL TEMPERATURE at REQUIRED HEIGHT -!C* ARDOUT(3) - SPECIFIC HUMIDITY at REQUIRED HEIGHT +!C* ARDOUT(2) - deficit of POTENTIAL TEMPERATURE between REQUIRED HEIGHT and the surface +!C* ARDOUT(3) - deficit of SPECIFIC HUMIDITY between REQUIRED HEIGHT and the surface !C* UFWIND - UNIVERSAL FUNCTION FOR WIND VELOCITY !C* UFTEMP - UNIVERSAL FUNCTION FOR SCALARS -! AR2(1) - NON-DIMENSIONAL CFL HIGHT = -! AR2(2) - RICHARDSON NUMBER = -! AR2(3) - REYNODS NUMBER = -! AR2(4) - LN(ZU/ZT) = -! AR2(5) - DYNAMICAL ROUGHNESS ZU (M) = -! AR2(6) - THERMAL ROUGHNESS ZT (M) = -! AR2(7) - CRITICAL RICHARDSON NUMBER = -! AR2(8) - TRANSFER COEFFICIENT FOR MOMENTUM = -! AR2(9) - TRANSFER COEFFICIENT FR HEAT = -! AR2(10)- COEFFICIENT OF TURBULENCE (KM) AT CFL HIGHT (M**2/S) = -! AR2(11)- ALFT=KT/KM ( KT-COEFFICIENT OF TURBULENCE FOR HEAT) = subroutine pbldia_new_sheba(AR2,ARDIN,ARDOUT) use sfx_sheba, only: & get_psi_sheba => get_psi - real,intent(in):: AR2(11),ARDIN(6)!ARDIN(10) - real,intent(out)::ARDOUT(3) !ARDOUT(4) -! real,parameter::Llim = 50.0 !minimum value of L for stable SL (25) + real,intent(in):: AR2(11),ARDIN(6) + real,intent(out)::ARDOUT(3) real,parameter::zetalim = 2. !maximum value of z/L for stable SL - real,parameter::kappa = 0.4 real psi_m,psi_h,psi_m_hs,psi_h_hs real hwind,htemp,ustar,tstar,qstar,& & ufwind,uftemp,dteta,dq,wind,L,HIN,zeta @@ -70,8 +55,7 @@ use sfx_sheba, only: & if(zeta.gt.zetalim) zeta=zetalim -! if(L.lt.Llim.and.L.gt.0) L=Llim - L = HIN/zeta ! MO length + L = HIN/zeta call get_psi_sheba(psi_m_hs,psi_h_hs,HIN/L) call get_psi_sheba(psi_m,psi_h,HWIND/L) @@ -89,18 +73,16 @@ use sfx_sheba, only: & ARDOUT(3) = DQ+ardin(3) - RETURN - END subroutine pbldia_new_sheba + return + end subroutine pbldia_new_sheba subroutine pbldia_new_most(AR2,ARDIN,ARDOUT) use sfx_most, only: & get_psi_most => get_psi - real,intent(in):: AR2(11),ARDIN(6)!ARDIN(10) - real,intent(out)::ARDOUT(3) !ARDOUT(4) -! real,parameter::Llim = 50.0 !minimum value of L for stable SL (25) + real,intent(in):: AR2(11),ARDIN(6) + real,intent(out)::ARDOUT(3) real,parameter::zetalim = 2. !maximum value of z/L for stable SL - real,parameter::kappa = 0.4 real psi_m,psi_h,psi_m_hs,psi_h_hs real hwind,htemp,ustar,tstar,qstar,& & ufwind,uftemp,dteta,dq,wind,L,HIN,zeta @@ -115,8 +97,7 @@ use sfx_most, only: & QSTAR = AR2(9) * ARDIN(3) if(zeta.gt.zetalim) zeta=zetalim -! if(L.lt.Llim.and.L.gt.0) L=Llim - L = HIN/zeta ! MO length + L = HIN/zeta call get_psi_most(psi_m_hs,psi_h_hs,HIN/L) call get_psi_most(psi_m,psi_h,HWIND/L) @@ -124,10 +105,6 @@ use sfx_most, only: & WIND = (USTAR/kappa) * UFWIND call get_psi_most(psi_m,psi_h,HTEMP/L) UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs) -! if(L>0)then -! !these functions in very stable stratification give inadequate temperatures at high levels; it's safer to use log profile for T -! UFTEMP = ALOG(HTEMP/HIN) -! endif DTETA = (TSTAR/kappa) * UFTEMP DQ = (QSTAR/kappa) * UFTEMP @@ -137,17 +114,13 @@ use sfx_most, only: & ARDOUT(3) = DQ+ardin(3) - RETURN - END subroutine pbldia_new_most + return + end subroutine pbldia_new_most subroutine pbldia_new_esm(AR2,ARDIN,ARDOUT) -!use sfx_esm, only: & -! get_psi_esm => get_psi - real,intent(in):: AR2(11),ARDIN(6)!ARDIN(10) - real,intent(out)::ARDOUT(3) !ARDOUT(4) -! real,parameter::Llim = 50.0 !minimum value of L for stable SL (25) + real,intent(in):: AR2(11),ARDIN(6) + real,intent(out)::ARDOUT(3) real,parameter::zetalim = 2. !maximum value of z/L for stable SL - real,parameter::kappa = 0.4 real psi_m,psi_h,psi_m_hs,psi_h_hs real hwind,htemp,ustar,tstar,qstar,& & ufwind,uftemp,dteta,dq,wind,L,HIN,zeta @@ -162,8 +135,7 @@ subroutine pbldia_new_esm(AR2,ARDIN,ARDOUT) QSTAR = AR2(9) * ARDIN(3) if(zeta.gt.zetalim) zeta=zetalim -! if(L.lt.Llim.and.L.gt.0) L=Llim - L = HIN/zeta ! MO length + L = HIN/zeta call get_psi_esm1(psi_m,psi_h,HIN,HWIND,L) @@ -171,10 +143,6 @@ subroutine pbldia_new_esm(AR2,ARDIN,ARDOUT) WIND = (USTAR/kappa) * UFWIND call get_psi_esm1(psi_m,psi_h,HIN,HTEMP,L) UFTEMP = psi_h -! if(L>0)then -! !these functions in very stable stratification give inadequate temperatures at high levels; it's safer to use log profile for T -! UFTEMP = ALOG(HTEMP/HIN) -! endif DTETA = (TSTAR/kappa) * UFTEMP DQ = (QSTAR/kappa) * UFTEMP @@ -184,68 +152,67 @@ subroutine pbldia_new_esm(AR2,ARDIN,ARDOUT) ARDOUT(3) = DQ+ardin(3) - RETURN - END subroutine pbldia_new_esm + return + end subroutine pbldia_new_esm subroutine get_psi_esm1(psi_m,psi_h,z1,z2,L) + use sfx_esm_param - !< @brief universal functions (momentum) & (heat): neutral case - ! ---------------------------------------------------------------------------- - real, intent(out) :: psi_m, psi_h !< universal functions - real, intent(in) :: z1,z2,L ! - real an5,D1,gmz1,gtz1,gmz2,gtz2,gmst,gtst,zeta1,zeta2 +real, intent(out) :: psi_m, psi_h +real, intent(in) :: z1,z2,L + +real an5,d1,gmz1,gtz1,gmz2,gtz2,gmst,gtst,zeta1,zeta2 an5 = (Pr_t_inf_inv/Pr_t_0_inv)**4 - D1=(2.0E0*alpha_h-an5*alpha_m-SQRT((an5*alpha_m)**2+4.0E0*an5*alpha_h*& + d1=(2.0E0*alpha_h-an5*alpha_m-sqrt((an5*alpha_m)**2+4.0E0*an5*alpha_h*& & (alpha_h-alpha_m)))/(2.0E0*alpha_h**2) zeta1 = z1 / L zeta2 = z2 / L - IF(zeta2 .GE. 0.) THEN + if(zeta2 .ge. 0.) then psi_m = alog(z2/z1) + beta_m*(zeta2-zeta1) psi_h = alog(z2/z1) + beta_m*(zeta2-zeta1) - ELSE - - GMZ2 = SQRT(SQRT(1.E0 - alpha_m * zeta2)) - GMZ1 = SQRT(SQRT(1.E0 - alpha_m * zeta1)) - GMST = SQRT(SQRT(1.E0 - alpha_m * D1)) - GTZ2 = SQRT(1.E0 - alpha_h * zeta2) - GTZ1 = SQRT(1.E0 - alpha_h * zeta1) - GTST = SQRT(1.E0 - alpha_h * D1) - - IF(zeta2 .GE. D1) THEN - psi_m = FIM(GMZ2) - FIM(GMZ1) - psi_h = FIT(GTZ2) - FIT(GTZ1) - ELSE - psi_m = (3.E0/GMST) * (1.E0 - (D1/zeta2)**(1.E0/3.E0)) + FIM(GMST) - FIM(GMZ1) - psi_h = (3.E0/GTST) * (1.E0 - (D1/zeta2)**(1.E0/3.E0)) + FIT(GTST) - FIT(GTZ1) - ENDIF - ENDIF - -end subroutine - - REAL FUNCTION FIM(XX) - IMPLICIT NONE - REAL X, XX - X = AMAX1(XX,1.000001E0) - FIM = ALOG((X-1.E0)/(X+1.E0)) + 2.E0*ATAN(X) - RETURN - END function - - REAL FUNCTION FIT(XX) - IMPLICIT NONE - REAL X, XX - X = AMAX1(XX,1.000001E0) - FIT = ALOG((X-1.E0)/(X+1.E0)) - RETURN - END function - - + else + + gmz2 = sqrt(sqrt(1.E0 - alpha_m * zeta2)) + gmz1 = sqrt(sqrt(1.E0 - alpha_m * zeta1)) + gmst = sqrt(sqrt(1.E0 - alpha_m * D1)) + gtz2 = sqrt(1.E0 - alpha_h * zeta2) + gtz1 = sqrt(1.E0 - alpha_h * zeta1) + gtst = sqrt(1.E0 - alpha_h * D1) + + if(zeta2 .ge. d1) then + psi_m = fim(gmz2) - fim(gmz1) + psi_h = fit(gtz2) - fit(gtz1) + else + psi_m = (3.E0/gmst) * (1.E0 - (d1/zeta2)**(1.E0/3.E0)) + fim(gmst) - fim(gmz1) + psi_h = (3.E0/gtst) * (1.E0 - (d1/zeta2)**(1.E0/3.E0)) + fit(gtst) - fit(gtz1) + endif + endif + +end subroutine get_psi_esm1 + +! functions fim,fit were copied from the original pbldia +REAL FUNCTION FIM(XX) + IMPLICIT NONE + REAL X, XX + X = AMAX1(XX,1.000001E0) + FIM = ALOG((X-1.E0)/(X+1.E0)) + 2.E0*ATAN(X) + RETURN +END function + +REAL FUNCTION FIT(XX) + IMPLICIT NONE + REAL X, XX + X = AMAX1(XX,1.000001E0) + FIT = ALOG((X-1.E0)/(X+1.E0)) + RETURN +END function end module pbldia_new_sfx \ No newline at end of file