diff --git a/diag/pbldia_new_sfx.f90 b/diag/pbldia_new_sfx.f90
deleted file mode 100755
index 249244aad983557355cdb79b6b2e39e74b538236..0000000000000000000000000000000000000000
--- a/diag/pbldia_new_sfx.f90
+++ /dev/null
@@ -1,231 +0,0 @@
-module pbldia_new_sfx
-
-
-implicit none 
-real,parameter,private::kappa=0.4,R=287.,appa=0.287
-
-contains
-
-! 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),
-!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(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)  - 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
-
-
-subroutine pbldia_new_sheba(AR2,ARDIN,ARDOUT)
-use sfx_sheba, only: &
-                get_psi_sheba => get_psi
-
-      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 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
-      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)
-     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)
-      real,intent(out)::ARDOUT(3)
-<<<<<<< HEAD
-      real,parameter::zetalim = 1. !maximum value of z/L for stable SL
-=======
-      real,parameter::zetalim = 2. !maximum value of z/L for stable SL
->>>>>>> 039d3258b422f93c9a44345c5152b600029afbf8
-      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)
-
-print *, zeta
-     if(zeta.gt.zetalim) zeta=zetalim
-      L = HIN/zeta
-<<<<<<< HEAD
-print *, zeta
-=======
->>>>>>> 039d3258b422f93c9a44345c5152b600029afbf8
-
-     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)
-
-      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)
-      real,intent(in):: AR2(11),ARDIN(6)
-      real,intent(out)::ARDOUT(3)
-<<<<<<< HEAD
-      real,parameter::zetalim = 1. !maximum value of z/L for stable SL
-=======
-      real,parameter::zetalim = 2. !maximum value of z/L for stable SL
->>>>>>> 039d3258b422f93c9a44345c5152b600029afbf8
-      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
-      L = HIN/zeta
-
-
-     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
-
-      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
-
-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*&
-&        (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 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