diff --git a/diag/pbldia_new_sfx.f90 b/diag/pbldia_new_sfx.f90 deleted file mode 100755 index 83f1b4000fbd79e3485cad36144ff8dd10b5f854..0000000000000000000000000000000000000000 --- a/diag/pbldia_new_sfx.f90 +++ /dev/null @@ -1,251 +0,0 @@ -module pbldia_new_sfx - - -implicit none -real,parameter,private::undef=-999.0,R=287.,appa=0.287 - -contains - -! in this subroutine 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), -!C* OF POTENTIAL TEMPERATURE DEFICIT (BETWEEN HEIGHT Z = HTEMP AND SURFACE), -!C* OF SPECIFIC HUMIDITY DEFICIT (BETWEEN HEIGHT Z = HTEMP AND SURFACE) -!C* INPUT DATA USED: -!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 -!C* ARDIN(1) - MODULE OF WIND VELOCITY AT THE TOP OF CONSTANT FLUX LAYER -!C* ARDIN(2) - POTENTIAL TEMPERATURE DEFICIT IN CONSTANT FLUX LAYER -!C* ARDIN(3) - SPECIFIC HUMIDITY DEFICIT IN CONSTANT FLUX LAYER -!C* ARDIN(4) - DIMENSIONAL HEIGHT OF CONSTANT FLUX LAYER TOP -!C* ARDIN(5) - = HWIND -!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* 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,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 - - HWIND = ARDIN(5) - HTEMP = ARDIN(6) - HIN = ardin(4) - zeta = AR2(1) - - USTAR = AR2(8) * ARDIN(1) - TSTAR = AR2(9) * ARDIN(2) - 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 - - call get_psi_sheba(psi_m_hs,psi_h_hs,HIN/L) - call get_psi_sheba(psi_m,psi_h,HWIND/L) - UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs) - WIND = (USTAR/kappa) * UFWIND - call get_psi_sheba(psi_m,psi_h,HTEMP/L) - UFTEMP = ALOG(HTEMP/HIN) - (psi_h - psi_h_hs) - - DTETA = (TSTAR/kappa) * UFTEMP - DQ = (QSTAR/kappa) * UFTEMP - - - ARDOUT(1) = ARDIN(1)+WIND - ARDOUT(2) = DTETA+ardin(2) - ARDOUT(3) = DQ+ardin(3) - - - 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,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 - - HWIND = ARDIN(5) - HTEMP = ARDIN(6) - HIN = ardin(4) - zeta = AR2(1) - - USTAR = AR2(8) * ARDIN(1) - TSTAR = AR2(9) * ARDIN(2) - 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 - - call get_psi_most(psi_m_hs,psi_h_hs,HIN/L) - call get_psi_most(psi_m,psi_h,HWIND/L) - UFWIND = ALOG(hwind/HIN) - (psi_m - psi_m_hs) - 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 - - ARDOUT(1) = ARDIN(1)+WIND - ARDOUT(2) = DTETA+ardin(2) - ARDOUT(3) = DQ+ardin(3) - - - 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,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 - - HWIND = ARDIN(5) - HTEMP = ARDIN(6) - HIN = ardin(4) - zeta = AR2(1) - - USTAR = AR2(8) * ARDIN(1) - TSTAR = AR2(9) * ARDIN(2) - 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 - - - call get_psi_esm1(psi_m,psi_h,HIN,HWIND,L) - UFWIND = psi_m - 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 - - ARDOUT(1) = ARDIN(1)+WIND - ARDOUT(2) = DTETA+ardin(2) - ARDOUT(3) = DQ+ardin(3) - - - 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 - - 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*& -& (alpha_h-alpha_m)))/(2.0E0*alpha_h**2) - - - zeta1 = z1 / L - zeta2 = z2 / L - - 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 - - - -end module pbldia_new_sfx \ No newline at end of file