   PROGRAM grenslaag
   !This program calculates the boundary layer height by the Richardson number method,
   !in dry (zi, as a function of potential temperature) and moist (zi1, as a function of virtual 
   !potential temperature
   ! More information. Jordi Vila (jvila@hp1.met.wau.nl)
   IMPLICIT NONE
   
   INTEGER			:: I, J, ki, li
   CHARACTER		 	:: time
   REAL, DIMENSION (26)		:: P, GH, U, V, W, T, THETA, r, q, RH, es, e
   REAL, DIMENSION (26)		:: u2, rich, rich1, thetae, thetav, zi, zi1, zi2, wind, aux1, aux2
   INTEGER, DIMENSION (26)	:: lev
   REAL				:: rcrit
   REAL				:: rhmax
   REAL, parameter		:: pref	= 1013.25, kapa = 0.286, gcons = 9.81, lv = 2.5E6, cp = 1004.
   
   !this opens the existing 3D file created in POST (output MM5 model)
   OPEN(1, FILE = '3D')
   
   !an output file is created, named zix11y44
   OPEN(2, FILE = 'zi3D')

   
   !choose which critical richardson number you want to use to determine the boundary layer height
   !rcrit = 0.3
   rcrit = 0.5
   !rcrit  = 1.0
   
   zi(I)  = 0.
   !boundary layer height (dry conditions)
   zi1(I) = 0.
   !boundary layer height (moist conditions)
   
   !this will be written to the output file, the FORMAT says it should be in the form of a table
   20 FORMAT ('time(UTC)', T12, 'zi', T27, 'zi(tv)', T42, 'rhmax')
   WRITE (2, 20) 
   
   !repeat for every time step
   DO J = 0, 54
   
   !read the first line (time) and the second line (names of variables, but there's no need to save this data)
   READ(1, *) time
   READ (1, *)
   
   ki = 1
   li = 1
   rhmax = 0.
   	
   !repeat for every level
   DO I = 1, 26
   		!this reads all the data (per level) of the 3Dfile
   		READ(1, *) lev(I), P(I), GH(I), U(I), V(I), W(I), T(I), THETA(I), r(I), q(I), RH(I), es(I), e(I)
   		
   		!calculation of the wind, **2 means square
   		wind(I)  = sqrt (U(I) ** 2 + V(I) **2)
   		u2(I) = wind(I) ** 2 + 0.01
   		
   		!calculation of the virtual potential temperature
  		thetav(I) = THETA(I) * (1 + 0.61 * q(I)) 
  	
  		
  		!the richardson number is calculated using the two different theta values.
  		rich(I)   = (gcons * GH(I) * (THETA(I) - THETA(1))) / (THETA(1) * u2(I))
  		!for dry conditions
  		rich1(I)  = (gcons * GH(I) * (thetav(I) - thetav(1))) / (thetav(1) * u2(I))
  		!for moist conditions
  		
  		!the maximum of RH per level
  		rhmax = Max (RH(I), rhmax)
   END DO
   
   !repeat for every level
   DO I = 1, 25
  		!look for richardson numbers smaller than the critical richardson number, for these cases, calculate zi
  		IF ((rich(ki) .LT. rcrit) .AND. (ki .LT. 25)) THEN
   			ki = ki + 1
   			aux1 = (rich(ki) - rcrit) / (rich(ki) - rich(ki-1))
   			aux2 = 1. - aux1
   			zi = aux2 * GH(ki) + aux1 * GH(ki-1)
   		END IF
   END DO
   
   	
   		
   DO I = 1, 25
	        
	        !this does the same as above, but with rich1 instead of rich
	        IF ((rich1(li) .LT. rcrit) .AND. (li .LT. 25)) THEN
   			li = li + 1
			aux1 = (rich1(li) - rcrit) / (rich1(li) - rich1(li-1))
   			aux2 = 1. - aux1
			zi1= aux2 * GH(li) + aux1 * GH(li-1)
   		END IF
   END DO
   
   
   !write zi, zi1 and rhmax to the output file
   10 FORMAT (I2, T12, F8.3, T27, F8.3, T42, F7.4)
   WRITE(2,10) J, zi(ki), zi1(li), rhmax
   
   END DO
   
   END PROGRAM grenslaag
