Index: /trunk/tools/Makefile.llamp
===================================================================
--- /trunk/tools/Makefile.llamp	(revision 2386)
+++ /trunk/tools/Makefile.llamp	(revision 2387)
@@ -95,5 +95,6 @@
 	pygenmods.o \
 	pydefmods.o \
-	trajectories_overlap.o
+	trajectories_overlap.o \
+	test_Fmods.o
 
 pymods : \
@@ -119,5 +120,5 @@
 
 clean :
-	$(RM) *.mod *.o interpolate module_ForDiag.so module_ForInt.so module_ForDistriCorrect.so module_ForSci.so module_ForGen.so module_ForDef.so
+	$(RM) *.mod *.o interpolate test_Fmods module_ForDiag.so module_ForInt.so module_ForDistriCorrect.so module_ForSci.so module_ForGen.so module_ForDef.so
 
 ########    #######
@@ -161,4 +162,7 @@
 	$(FC) $(FCEXECFLAGS) DistriCorrection.f90 -o DistriCorrection
 
+test_Fmods.o: module_definitions.o module_basic.o module_generic.o module_scientific.o module_ForDiagnostics.o module_ForDiagnosticsVars.o
+	$(FC) $(FCEXEFCFLAGS) test_Fmods.f90 -o test_Fmods
+
 pydistrimods.o:
 	$(F2PY) -c $(F2PYF) -I$(NCINCFOLD) -m module_ForDistriCorrect $(distrisrcs) -L$(NCLIBFOLD)
Index: /trunk/tools/diag_tools.py
===================================================================
--- /trunk/tools/diag_tools.py	(revision 2386)
+++ /trunk/tools/diag_tools.py	(revision 2387)
@@ -71,4 +71,5 @@
 # var_td_Uhus: Function to compute dew-point air temperature from temperature and 
 #   pressure values using hus
+# var_tws_S11: Function to compute Wet Bulb temperature after Stull, 2011
 # var_timemax: This function returns the time at which variable reaches its maximum in time
 #   units
@@ -2176,4 +2177,25 @@
 
     return td
+
+def var_tws_S11(ta0, hur0):
+    """ Function to compute Wet Bulb temperature using equation after:
+      Stull, R. (2011), J. Appl. Meteor. Climatol. 50(11):2267-2269,
+        doi: 10.1175/JAMC-D-11-0143.1
+      [ta]= temperature (assuming [[t],z,y,x] in [K])
+      [hur] = relative humidity (assuming in [1])
+    >>> var_rh_S11(293.15, 0.5)
+    13.699341969
+    """
+    fname = 'var_tws_S11'
+
+    ta = ta0 - 273.15
+    hur = hur0*100.
+
+    # Does it has any sense not in surface?
+    tws = ta*np.arctan(0.151977*np.sqrt(hur+8.313659)) + np.arctan(ta+hur) -         \
+      np.arctan(hur-1.676331) + 0.00391838*(hur)**(1.5)*np.arctan(0.023101*hur) -    \
+      4.686035
+
+    return tws
 
 def var_wd(u, v):
Index: /trunk/tools/diagnostics.py
===================================================================
--- /trunk/tools/diagnostics.py	(revision 2386)
+++ /trunk/tools/diagnostics.py	(revision 2387)
@@ -928,4 +928,13 @@
         ncvar.insert_variable(ncobj, 'lat_bnds', cellbndsy, diagoutd, diagoutvd, newnc)
 
+# tws: Bet Wulb temperature following Stull 2011 (tas, hur)
+    elif diagn == 'rws':
+            
+        var0 = ncobj.variables[depvars[0]][:]
+        var1 = ncobj.variables[depvars[1]][:]
+
+        diagout, diagoutd, diagoutvd = diag.var_tws_S11(var0,var1,dnames,dvnames)
+        ncvar.insert_variable(ncobj, 'tws', diagout, diagoutd, diagoutvd, newnc)
+
 # cell_bnds: grid cell bounds from XLONG_U, XLAT_U, XLONG_V, XLAT_V as intersection 
 #   of their related parallels and meridians
@@ -1056,4 +1065,26 @@
             varsadd.append(nonvd)
         ncvar.insert_variable(ncobj, 'psl', diagout, diagoutd, diagoutvd, newnc)
+
+# WRFtws: Bet Wulb temperature following Stull 2011 (PSFC, T2, Q2)
+    elif diagn == 'WRFtws':
+
+        var0 = ncobj.variables[depvars[0]][:]
+        var1 = ncobj.variables[depvars[1]][:]
+        var2 = ncobj.variables[depvars[2]][:]
+
+        dnamesvar = list(ncobj.variables[depvars[2]].dimensions)
+        dvnamesvar = ncvar.var_dim_dimv(dnamesvar,dnames,dvnames)
+
+        hurs, diagoutd, diagoutvd = diag.compute_rh(var0,var1,var2,dnamesvar,dvnamesvar)
+            
+        diagout = diag.var_tws_S11(var1, hurs)
+        # Removing the nonChecking variable-dimensions from the initial list
+        varsadd = []
+        diagoutvd = list(dvnames)
+        for nonvd in NONchkvardims:
+            if gen.searchInlist(dvnames,nonvd): diagoutvd.remove(nonvd)
+            varsadd.append(nonvd)
+
+        ncvar.insert_variable(ncobj, 'tws', diagout, diagoutd, diagoutvd, newnc)
 
 # OMEGAw (omega, p, t) from NCL formulation (https://www.ncl.ucar.edu/Document/Functions/Contributed/omega_to_w.shtml)
Index: /trunk/tools/module_ForDiagnostics.f90
===================================================================
--- /trunk/tools/module_ForDiagnostics.f90	(revision 2386)
+++ /trunk/tools/module_ForDiagnostics.f90	(revision 2387)
@@ -35,4 +35,5 @@
 ! compute_psl_ecmwf: Compute sea level pressure using ECMWF method following Mats Hamrud and Philippe Courtier [Pa]
 ! compute_range_faces: Subroutine to compute faces [uphill, valleys, downhill] of a mountain range along a given face
+! compute_tws_RK[1/2/3/4]: Subroutine to compute Wet Bulb temperature of 1/2/3/4D series of values
 ! compute_vertint1D: Subroutine to vertically integrate a 1D variable in any vertical coordinates
 ! compute_zint4D: Subroutine to vertically integrate a 4D variable in any vertical coordinates
@@ -845,5 +846,5 @@
 ! Subroutine to compute fog: qcloud + qice /= 0.
 ! And visibility following Kunkel, B. A., (1984): Parameterization of droplet terminal velocity and 
-!   extinction coefficient in fog models. J. Climate Appl. Meteor., 23, 34â41.
+!   extinction coefficient in fog models. J. Climate Appl. Meteor., 23, 34-41.
 
     IMPLICIT NONE
@@ -1357,3 +1358,123 @@
   END SUBROUTINE compute_Koeppen_Geiger_climates
 
+  SUBROUTINE compute_tws_RK1(d1, tas, hurs, tws)
+! Subroutine to compute Wet Bulb temperature of 1D series of values using equation after:
+!    Stull, R. (2011), J. Appl. Meteor. Climatol. 50(11):2267-2269. doi: 10.1175/JAMC-D-11-0143.1
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: d1
+    REAL(r_k), DIMENSION(d1), INTENT(in)                 :: tas, hurs
+    REAL(r_k), DIMENSION(d1), INTENT(out)                :: tws
+ 
+! Local
+    INTEGER                                              :: it
+
+!!!!!!! Variables
+! tas: 2-m air temperature [K]
+! hurs: 2-m relative humidity [1]
+
+    fname = 'compute_tws_RK1'
+
+    DO it=1, d1
+      tws(it) = var_tws_S11(tas(it), hurs(it))
+    END DO
+
+    RETURN
+
+  END SUBROUTINE compute_tws_RK1
+
+  SUBROUTINE compute_tws_RK2(d1, d2, tas, hurs, tws)
+! Subroutine to compute Wet Bulb temperature of 2D series of values using equation after:
+!    Stull, R. (2011), J. Appl. Meteor. Climatol. 50(11):2267-2269. doi: 10.1175/JAMC-D-11-0143.1
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: d1, d2
+    REAL(r_k), DIMENSION(d1,d2), INTENT(in)              :: tas, hurs
+    REAL(r_k), DIMENSION(d1,d2), INTENT(out)             :: tws
+ 
+! Local
+    INTEGER                                              :: i, j
+
+!!!!!!! Variables
+! tas: 2-m air temperature [K]
+! hurs: 2-m relative humidity [1]
+
+    fname = 'compute_tws_RK2'
+
+    DO i=1, d1
+      DO j=1, d2
+        tws(i,j) = var_tws_S11(tas(i,j), hurs(i,j))
+      END DO
+    END DO
+
+    RETURN
+
+  END SUBROUTINE compute_tws_RK2
+
+  SUBROUTINE compute_tws_RK3(d1, d2, d3, tas, hurs, tws)
+! Subroutine to compute Wet Bulb temperature of 3D series of values using equation after:
+!    Stull, R. (2011), J. Appl. Meteor. Climatol. 50(11):2267-2269. doi: 10.1175/JAMC-D-11-0143.1
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: d1, d2, d3
+    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(in)           :: tas, hurs
+    REAL(r_k), DIMENSION(d1,d2,d3), INTENT(out)          :: tws
+ 
+! Local
+    INTEGER                                              :: i, j, k
+
+!!!!!!! Variables
+! tas: 2-m air temperature [K]
+! hurs: 2-m relative humidity [1]
+
+    fname = 'compute_tws_RK3'
+
+    DO i=1, d1
+      DO j=1, d2
+        DO k=1, d3
+          tws(i,j,k) = var_tws_S11(tas(i,j,k), hurs(i,j,k))
+        END DO
+      END DO
+    END DO
+
+    RETURN
+
+  END SUBROUTINE compute_tws_RK3
+
+  SUBROUTINE compute_tws_RK4(d1, d2, d3, d4, tas, hurs, tws)
+! Subroutine to compute Wet Bulb temperature of 4D series of values using equation after:
+!    Stull, R. (2011), J. Appl. Meteor. Climatol. 50(11):2267-2269. doi: 10.1175/JAMC-D-11-0143.1
+
+    IMPLICIT NONE
+
+    INTEGER, INTENT(in)                                  :: d1, d2, d3, d4
+    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(in)        :: tas, hurs
+    REAL(r_k), DIMENSION(d1,d2,d3,d4), INTENT(out)       :: tws
+ 
+! Local
+    INTEGER                                              :: i,j,k,l
+
+!!!!!!! Variables
+! tas: 2-m air temperature [K]
+! hurs: 2-m relative humidity [1]
+
+    fname = 'compute_tws_RK4'
+
+    DO i=1, d1
+      DO j=1, d2
+       DO k=1, d3
+         DO l=1, d4
+           tws(i,j,k,l) = var_tws_S11(tas(i,j,k,l), hurs(i,j,k,l))
+         END DO
+        END DO
+      END DO
+    END DO
+
+    RETURN
+
+  END SUBROUTINE compute_tws_RK4
+
 END MODULE module_ForDiagnostics
Index: /trunk/tools/module_ForDiagnosticsVars.f90
===================================================================
--- /trunk/tools/module_ForDiagnosticsVars.f90	(revision 2386)
+++ /trunk/tools/module_ForDiagnosticsVars.f90	(revision 2387)
@@ -38,4 +38,5 @@
 ! var_range_faces: compute faces [uphill, valleys, downhill] of a monuntain range along a face
 ! var_rh: Subroutine to compute relative humidity following 'Tetens' equation (T,P) ...'
+! var_tws_S11: Function to compute wet bulb temperature after Stull, 2011
 ! var_zmla_generic: Subroutine to compute pbl-height following a generic method
 ! var_zwind: Subroutine to extrapolate the wind at a given height following the 'power law' methodology
@@ -1918,3 +1919,60 @@
   END SUBROUTINE var_hur
 
+  REAL(r_k) FUNCTION var_tws_S11(ta0, hur0)
+! Function to compute wet bulb temperature using equation after:
+!    Stull, R. (2011), J. Appl. Meteor. Climatol. 50(11):2267-2269. doi: 10.1175/JAMC-D-11-0143.1
+
+    IMPLICIT NONE
+
+    REAL(r_k), INTENT(in)                                :: ta0, hur0
+
+! Local
+    REAL                                                 :: ta, hur
+
+!!!!!!! Variables
+! ta0: temperature [K] (does it only make sense if it is at 2 m?)
+! hur0: relative humidity [1] (does it only make sense if it is at 2 m?)
+! tws: wet bulb temperature [C]
+
+    fname = 'var_tws_S11'
+
+    ta = ta0 - SVPT0
+    hur = hur0*100.*oneRK
+
+    var_tws_S11 = ta * ATAN(0.151977*SQRT(hur+8.313659)) + ATAN(ta+hur) - ATAN(hur-1.676331) +        &
+      0.00391838*(hur)**(1.5)*ATAN(0.023101*hur) - 4.686035
+ 
+    RETURN
+
+  END FUNCTION var_tws_S11
+
+  SUBROUTINE Svar_tws_S11(ta0, hur0, tws)
+! Subroutine to compute wet bulb temperature using equation after:
+!    Stull, R. (2011), J. Appl. Meteor. Climatol. 50(11):2267-2269. doi: 10.1175/JAMC-D-11-0143.1
+
+    IMPLICIT NONE
+
+    REAL(r_k), INTENT(in)                                :: ta0, hur0
+    REAL(r_k), INTENT(out)                               :: tws
+
+! Local
+    REAL                                                 :: ta, hur
+
+!!!!!!! Variables
+! ta0: temperature [K] (does it only make sense if it is at 2 m?)
+! hur0: relative humidity [1] (does it only make sense if it is at 2 m?)
+! tws: wet bulb temperature [C]
+
+    fname = 'var_tws_S11'
+
+    ta = ta0 - SVPT0
+    hur = hur0*100.*oneRK
+
+    tws = ta * ATAN(0.151977*SQRT(hur+8.313659)) + ATAN(ta+hur) - ATAN(hur-1.676331) +                &
+      0.00391838*(hur)**(1.5)*ATAN(0.023101*hur) - 4.686035
+ 
+    RETURN
+
+  END SUBROUTINE Svar_tws_S11
+
 END MODULE module_ForDiagnosticsVars
Index: /trunk/tools/variables_values.dat
===================================================================
--- /trunk/tools/variables_values.dat	(revision 2386)
+++ /trunk/tools/variables_values.dat	(revision 2387)
@@ -541,4 +541,6 @@
 tottotalsi, tottoalsi, totals_totals_index, 0., 100., Totals|totals|index, C, Reds, $TOTtoti$, TOTtoti
 Totals totals index, tottoalsi, totals_totals_index, 0., 100., Totals|totals|index, C, Reds, $TOTtoti$, TOTtoti
+tws, tws, air_wet_bulb_temperature, -10., 40., air|wet|bulb|temperature|at|2m, C, YlGnBu, $tws$, tws
+TWS, tws, air_wet_bulb_temperature, -10., 40., air|wet|bulb|temperature|at|2m, C, YlGnBu, $tws$, tws
 ua, ua, eastward_wind, -30., 30., eastward|wind, ms-1, seismic, $ua$, ua
 vitu, ua, eastward_wind, -30., 30., eastward|wind, ms-1, seismic, $ua$, ua
