SUEWS API Site
Documentation of SUEWS source code
suews_phys_bluews.f95
Go to the documentation of this file.
1 !======================================================================================================
2 MODULE cbl_module
3 
4  INTEGER::entrainmenttype, & ! Entrainment type choice
5  co2_included, & ! CO2 included
6  initialdata_use, & ! 1 read initial data, 0 do not
7  ! qh_choice,& ! selection of qh use to drive CBL growth 1=Suews 2=lumps 3=obs ! moved to suews_data
8  sondeflag, & ! 1 read sonde or vertical profile data in 0 do not
9  isubs ! 1 include subsidence in equations
10 
11  INTEGER, DIMENSION(366)::cblday = 0
12 
13  CHARACTER(len=200), DIMENSION(366)::filesonde = ""
14  CHARACTER(len=200)::initialdatafilename
15  REAL(KIND(1D0)):: wsb ! subsidence velocity
16  REAL(KIND(1d0)), DIMENSION(1:10):: cbldata
17  REAL(KIND(1d0)), DIMENSION(:, :), ALLOCATABLE::inicbldata
18 
19  !Parameters in CBL code
20  INTEGER::zmax, &
21  neqn = 6, & !NT changed from 4 to 6
22  icblcount, &
24  REAL(KIND(1d0))::c2k = 273.16
25 
26  REAL(KIND(1D0)):: usbl, ftbl, fqbl, fcbl, gamt, gamq, gamc, tpp, qpp, cp0!,tk
27 
28  REAL(KIND(1D0))::alpha3, &
29  blh_m, & ! Boundary layer height(m)
30  blh1_m, &
31  cm, & ! CO2 concentration in CBL
32  !cp0,gamc,& !
33  gamt_km, & ! Vertical gradient of theta (K/m)
34  gamq_gkgm, &! Vertical gradient of specific humidity (g/kg/m)
35  gamq_kgkgm, &! Vertical gradient of specific humidity (kg/kg/m)
36  !fcbl,&
37  tm_c, & ! Potential temperature in CBL (degree Celsius)
38  tm_k, & ! Potential temperature in CBL (K)
39  tmp_k, &
40  tp_c, & ! Potential temperature just above Boundary layer height(degree Celsius)
41  tp_k, & ! Potential temperature just above Boundary layer height(K)
42  tpp_k, &
43  febl_kgkgms, &! Kinematic latent heat flux((kg/kg)*m/s)
44  fhbl_kms, & ! Kinematic sensible heat flux(K*m/s)
45  qm_gkg, & ! Specific humidity in CBL(g/kg)
46  qm_kgkg, & ! Specific humidity in CBL(kg/kg)
47  qp_gkg, & ! Specific humidity above Boundary layer height(g/kg)
48  qp_kgkg, & ! Specific humidity above Boundary layer height(kg/kg)
49  qpp_kgkg
50 
51  REAL(KIND(1D0)), DIMENSION(0:500, 2):: gtheta, ghum ! Vertical gradient of theta and specific humidity from sonde data
52  REAL(KIND(1D0)), DIMENSION(6)::y ! NT set from 4 to 6
53 
54 END MODULE cbl_module
55 !===================================================================================
56 
58  USE cbl_module
59  USE meteo, ONLY: qsatf, sat_vap_press_x
60 
61  IMPLICIT NONE
62 CONTAINS
63  ! Note: INTERVAL is now set to 3600 s in Initial (it is no longer set in RunControl) HCW 29 Jan 2015
64  ! Last modified:
65  ! NT 6 Apr 2017 - include top of the CBL variables in RKUTTA scheme + add flag to include or exclude subsidence
66  ! HCW 29 Mar 2017 - Changed third dimension of dataOutBL to Gridiv (was previously iMB which seems incorrect)
67  ! LJ 27 Jan 2016 - Removal of tabs
68 
69  SUBROUTINE cbl(iy, id, it, imin, ir, Gridiv, qh_choice, dectime, Temp_C, Press_hPa, avkdn, avu1, avrh, &
70  avcp, avdens, es_hPa, lv_J_kg, nsh_real, tstep, UStar, psih, is, NumberOfGrids, qhforCBL, qeforCBL, &
71  ReadLinesMetdata, dataOutBL)
72  IMPLICIT NONE
73  INTEGER, PARAMETER:: ncolumnsdataOutBL = 22
74 
75  INTEGER, INTENT(IN):: tstep, is, NumberOfGrids, Gridiv, ReadLinesMetdata, ir
76  REAL(KIND(1d0)), INTENT(IN), DIMENSION(NumberOfGrids):: qhforCBL, qeforCBL
77  REAL(KIND(1d0)), INTENT(IN):: avkdn, nsh_real, UStar, psih
78  INTEGER, INTENT(INOUT) :: qh_choice, iy, id, it, imin
79  REAL(KIND(1d0)), INTENT(INOUT):: dectime, Press_hPa, avu1, avrh, es_hPa, avcp, avdens, lv_J_kg
80  REAL(KIND(1d0)), INTENT(OUT):: Temp_C
81  REAL(KIND(1d0)), INTENT(OUT), DIMENSION(ReadLinesMetdata, ncolumnsdataOutBL, NumberOfGrids) ::dataOutBL
82 
83  REAL(KIND(1d0)):: gas_ct_dry = 8.31451/0.028965 !j/kg/k=dry_gas/molar
84  ! REAL(KIND(1d0)):: gas_ct_wv = 8.31451/0.0180153 !j/kg/kdry_gas/molar_wat_vap
85  REAL(KIND(1d0))::qh_use, qe_use, tm_K_zm, qm_gkg_zm
86  REAL(KIND(1d0))::Temp_C1, avrh1, es_hPa1
87  REAL(KIND(1d0))::secs0, secs1, Lv
88  REAL(KIND(1d0))::NAN = -999
89  INTEGER::idoy, startflag
90 
91  ! initialise startflag
92  startflag = 0
93 
94  ! print *, 'iCBLcount', iCBLcount,ifirst
95  ! print *,iy, id,it,imin
96  ! print *, cbldata(1),tm_K
97 
98  ! Reset iCBLcount at start of each metblock (HCW added 29/03/2017)
99  ! IF (ifirst == 1) THEN
100  ! iCBLcount = 0
101  ! ENDIF
102  ! print*,IniCBLdata
103  !write(*,*) DateTIme
104  !Skip first loop and unspecified days
105  !IF((ifirst==1 .AND. iMB==1) .OR. CBLday(id)==0) THEN !HCW modified condition to check for first timestep of the model run
106  ! IF (ifirst == 1 .OR. IniCBLdata(id, 2) == -999) THEN !TS modified to adapt the format of the new CBL_initial file
107  ! iCBLcount = iCBLcount + 1
108  ! !write(*,*) 'ifirst or nonCBLday', DateTime, iCBLcount
109  ! dataOutBL(iCBLcount, 1:ncolumnsdataOutBL, Gridiv) = (/REAL(iy, 8), REAL(id, 8), REAL(it, 8), REAL(imin, 8), dectime, &
110  ! (NAN, is=6, ncolumnsdataOutBL)/)
111  ! RETURN
112  ! ELSEIF (avkdn < 5) THEN
113  ! iNBL = 1
114  ! IF (iNBL == -9) THEN
115  ! CALL CBL_initial(qh_use, qe_use, tm_K_zm, qm_gkg_zm, startflag, ir, Gridiv)
116  ! RETURN
117  ! ELSE
118  ! !ADD NBL for Night:(1)Fixed input/output NBL; (2) Input Fixed Theta,Q to SUEWS; (3) Currently NBL eq 200 m
119  ! CALL NBL(qh_use, qe_use, tm_K_zm, qm_gkg_zm, startflag, Gridiv)
120  ! RETURN
121  ! ENDIF
122  ! ENDIF
123 
124  IF (avkdn < 5) THEN
125  ! iNBL = 1
126  ! IF (iNBL == -9) THEN
127  ! CALL CBL_initial(qh_use, qe_use, tm_K_zm, qm_gkg_zm, startflag, Gridiv)
128  ! RETURN
129  ! ELSE
130  !ADD NBL for Night:(1)Fixed input/output NBL; (2) Input Fixed Theta,Q to SUEWS; (3) Currently NBL eq 200 m
131  CALL nbl(iy, id, it, imin, dectime, ir, qh_choice, qh_use, qe_use, tm_k_zm, qm_gkg_zm, startflag, gridiv, &
132  psih, ustar, temp_c, &
133  numberofgrids, qhforcbl, qeforcbl, &
134  press_hpa, avu1, avrh, &
135  readlinesmetdata, dataoutbl, &
136  avcp, avdens, es_hpa, lv_j_kg)
137  RETURN
138  ! ENDIF
139  ENDIF
140 
141  IF (startflag == 0) THEN !write down initial values in previous time step
142  !write(*,*) 'startflag', DateTime, iCBLcount
143  dataoutbl(ir, 1:ncolumnsdataoutbl, gridiv) &
144  = (/REAL(iy, 8), REAL(id, 8), REAL(it, 8), REAL(imin, 8), dectime, blh_m, tm_K, &
145  qm_kgkg*1000, tp_K, qp_kgkg*1000, (nan, is=11, 20), gamt_km, gamq_kgkgm/)
146  startflag = 1
147  ENDIF
148 
149  qh_use = qhforcbl(gridiv) !HCW 21 Mar 2017
150  qe_use = qeforcbl(gridiv)
151  IF (qh_use < -900 .OR. qe_use < -900) THEN ! observed data has a problem
152  CALL errorhint(22, 'Unrealistic qh or qe_value for CBL in CBL.', qh_use, qe_use, qh_choice)
153  ENDIF
154  !!Heat flux choices - these are now made in SUEWS_Calculations for qhforCBL and qeCBL, rather than here
155  !IF(Qh_choice==1) THEN !from SUEWS
156  ! !qh_use=qh
157  ! !qe_use=qeph
158  ! qh_use=qhforCBL(Gridiv) !HCW 21 Mar 2017
159  ! qe_use=qeforCBL(Gridiv)
160  !ELSEIF(qh_choice==2)THEN !from LUMPS
161  ! qh_use=H_mod
162  ! qe_use=E_mod
163  !ELSEIF(qh_choice==3)THEN !from OBS
164  ! IF(qh_obs<-900.OR.qe_obs<-900)THEN ! observed data has a problem
165  ! CALL ErrorHint(22,'Unrealistic observed qh or qe_value.',qh_obs,qe_obs,qh_choice)
166  ! ENDIF
167  ! qh_use=qh_obs
168  ! qe_use=qe_obs
169  !ENDIF
170 
171  !-------Main loop of CBL calculation--------------------------------------
172  !-------------------------------------------------------------------------
173 
174  !write(*,*) 'Main CBL loop'
175 
176  cbldata(1) = float(it) + float(imin)/60.
177  cbldata(2) = qh_use
178  cbldata(3) = qe_use
179  cbldata(4) = avdens
180  cbldata(5) = lv_j_kg
181  cbldata(6) = avcp
182  cbldata(7) = avu1
183  cbldata(8) = ustar
184  cbldata(9) = press_hpa
185  cbldata(10) = psih
186 
187  ! print *, 'cbldata'
188  ! print *, cbldata(1),tm_K
189 
190  secs0 = cbldata(1)*3600.
191  secs1 = secs0 + float(tstep) ! time in seconds
192  ! Kinematic fluxes
193  fhbl_kms = cbldata(2)/(cbldata(4)*cbldata(6)) !qh_use/(avdens*avcp) ! units: degK * m/s
194  febl_kgkgms = cbldata(3)/(cbldata(4)*cbldata(5)) !qe_use/(avdens*lv_J_kg) ! units: kg/kg * m/s
195  IF (co2_included == 1) THEN
196  fcbl = 0!fc(i)/(rmco2/volm) ! units: mol/mol * m/s
197  ELSE
198  cm = nan
199  ENDIF
200 
201  ! tpp_K=tp_K
202  ! qpp_kgkg=qp_kgkg
203 
204  IF (sondeflag == 1) THEN
205  CALL gamma_sonde
206  ENDIF
207  ! set up array for Runge-Kutta call
208  blh1_m = blh_m
209  y(1) = blh_m ! integrate h, t, q, c from time s(i-1)
210  y(2) = tm_k ! to time s(i) using runge-kutta solution
211  y(3) = qm_kgkg ! of slab CBL equations
212  y(4) = cm
213  y(5) = tp_k
214  y(6) = qp_kgkg
215 
216  !++++++++++++++++++++++++++++++++++++++++++++++++++++++
217  CALL rkutta(neqn, secs0, secs1, y, 1)
218  !++++++++++++++++++++++++++++++++++++++++++++++++++++++
219  blh_m = y(1)
220  tm_k = y(2) ! potential temperature, units: K
221  qm_kgkg = y(3) ! specific humidity, units: kg/kg
222  cm = y(4) ! co2 concentration,units: mol/mol
223  tp_k = y(5) ! potential temperature top of CBL: K
224  qp_kgkg = y(6) ! specific humidity top of CBL: kg/kg
225  ! compute derived quantities for this time step
226 
227  !NT: now included in rkutta
228  ! tp_K = tpp_K + (gamt_Km*(blh_m-blh1_m))
229  ! qp_kgkg = qpp_kgkg + (gamq_kgkgm*(blh_m-blh1_m))
230 
231  ! IF (tp_K.LT.tm_K) THEN
232  ! tp_K = tm_K
233  ! ENDIF
234 
235  ! th = tm_K - (grav/cbldata(5))*blh_m ! actual temp just below z=h
236  ! dh = qsatf(th,cbldata(8)) - qm_kgkg ! deficit just below z=h
237 
238  tp_c = tp_k - c2k
239  tm_c = tm_k - c2k
240 
241  ! delt = tp_K - tm_K ! temp
242  ! delq = qp_kgkg - qm_kgkg ! humidity
243  !deltv = (tp_K - tm_K) + 0.61*tm_k*(qp_kgkg - qm_kgkg) ! pot virtual temp
244 
245  qm_gkg = qm_kgkg*1000 !humidities: kg/kg -> g/kg
246 
247  !Output time correction
248  idoy = id
249  !If(it==0 .and. imin==55) idoy=id-1
250  IF (it == 0 .AND. imin == (nsh_real - 1)/nsh_real*60) idoy = id - 1 !Modified by HCW 04 Mar 2015 in case model timestep is not 5-min
251 
252  ! QUESTION: any difference between the two options? code looks the same in the two branches.
253  IF ((qh_choice == 1) .OR. (qh_choice == 2)) THEN !BLUEWS or BLUMPS
254  !Stability correction
255  !tm_K_zm=tm_K+cbldata(10)*cbldata(2)/(k*cbldata(8)*cbldata(6)*cbldata(4))
256  temp_c = tm_k/((1000/cbldata(9))**(gas_ct_dry/cbldata(6))) - c2k
257  es_hpa = sat_vap_press_x(temp_c, cbldata(9), 1, dectime)
258  lv = (2500.25 - 2.365*temp_c)*1000
259  !qm_gkg_zm=qm_gkg+cbldata(10)*cbldata(3)/(k*cbldata(8)*cbldata(4)*lv)
260  avrh = 100*((qm_gkg*cbldata(9)/(622 + qm_gkg))/es_hpa) !check pressure
261  IF (avrh > 100) THEN
262  CALL errorhint(34, 'subroutine CBL dectime, relative humidity', idoy + cbldata(1)/24.0, avrh, 100)
263  avrh = 100
264  ENDIF
265  ! iCBLcount = iCBLcount + 1
266  ! write(*,*) 'qh1or2', iy,id,it,imin, iCBLcount
267  dataoutbl(ir, 1:ncolumnsdataoutbl, gridiv) &
268  = (/REAL(iy, 8), REAL(id, 8), REAL(it, 8), REAL(imin, 8), dectime, blh_m, tm_k, &
269  qm_kgkg*1000, tp_K, qp_kgkg*1000, &
270  Temp_C, avrh, cbldata([2, 3, 9, 7, 8, 4, 5, 6]), &
271  gamt_Km, gamq_kgkgm/)
272  ELSEIF (qh_choice == 3) THEN ! CBL
273  !tm_K_zm=tm_K+cbldata(10)*cbldata(2)/(k*cbldata(8)*cbldata(6)*cbldata(4))
274  temp_c1 = tm_k/((1000/cbldata(9))**(gas_ct_dry/cbldata(6))) - c2k
275  es_hpa1 = sat_vap_press_x(temp_c1, cbldata(9), 1, dectime)
276  lv = (2500.25 - 2.365*temp_c1)*1000
277  !qm_gkg_zm=qm_gkg+cbldata(10)*cbldata(3)/(k*cbldata(8)*cbldata(4)*lv)
278  avrh1 = 100*((qm_gkg*cbldata(9)/(622 + qm_gkg))/es_hpa1) ! should be cbldata(9), i.e., Press_hPa
279  IF (avrh1 > 100) THEN
280  CALL errorhint(34, 'subroutine CBL dectime, relative humidity', idoy + cbldata(1)/24.0, avrh1, 100)
281  avrh1 = 100
282  ENDIF
283  ! iCBLcount = iCBLcount + 1
284  !write(*,*) 'qh3', DateTIme, iCBLcount
285  dataoutbl(ir, 1:ncolumnsdataoutbl, gridiv) &
286  = (/REAL(iy, 8), REAL(id, 8), REAL(it, 8), REAL(imin, 8), dectime, blh_m, tm_K, &
287  qm_kgkg*1000, tp_K, qp_kgkg*1000, &
288  Temp_C1, avrh1, cbldata([2, 3, 9, 7, 8, 4, 5, 6]), &
289  gamt_Km, gamq_kgkgm/)
290  ENDIF
291  ! move the counter at the end, TS 27 Aug 2019
292  ! iCBLcount = iCBLcount + 1
293 
294  RETURN
295 
296  END SUBROUTINE cbl
297 
298  !-----------------------------------------------------------------------
299  !-----------------------------------------------------------------------
300  SUBROUTINE cbl_readinputdata(FileInputPath, qh_choice)
302  IMPLICIT NONE
303  INTEGER, INTENT(INOUT)::qh_choice
304  CHARACTER(len=150), INTENT(IN):: FileInputPath
305 
306  INTEGER::i, ios
307  REAL(KIND(1d0))::l
308 
309  namelist /cblinput/ entrainmenttype, &
310  qh_choice, &
311  isubs, &
312  co2_included, &
313  cblday, &
314  wsb, &
315  initialdata_use, &
317  sondeflag, &
318  filesonde
319 
320  OPEN (51, file=trim(fileinputpath)//'CBLInput.nml', status='old', err=24)
321  READ (51, nml=cblinput, err=24)
322  CLOSE (51)
323 
324  !Read initial values if it's needed
325  IF (initialdata_use == 1 .OR. initialdata_use == 2) THEN
326  OPEN (52, file=trim(fileinputpath)//trim(initialdatafilename), status='old', err=25)
327  READ (52, *)
328  nlineindata = 0 !Initialise nlines
329  DO
330  READ (52, *, iostat=ios) l
331  IF (ios < 0 .OR. l == -9) EXIT !IF (l == -9) EXIT
333  ENDDO
334  CLOSE (52)
335 
336  if (allocated(inicbldata)) deallocate (inicbldata)
337  ALLOCATE (inicbldata(1:nlineindata, 1:8))
338  OPEN (52, file=trim(fileinputpath)//trim(initialdatafilename), status='old', err=25)
339  READ (52, *)
340  DO i = 1, nlineindata
341  READ (52, *) inicbldata(i, 1:8)
342  ENDDO
343  CLOSE (52)
344  ENDIF
345 
346  IF (co2_included == 0) THEN
347  fcbl = 0 ! hard-wire no CO2
348  ENDIF
349 
350  ! iCBLcount = 1
351 
352  RETURN
353 
354 24 CALL errorhint(24, 'CBLInput.nml', 0.00d0, 0.000d0, 0)
355 25 CALL errorhint(24, trim(fileinputpath)//trim(initialdatafilename), 0.00d0, 0.00d0, 0)
356 
357  END SUBROUTINE cbl_readinputdata
358 
359  !----------------------------------------------------------------------
360  !-----------------------------------------------------------------------
361  SUBROUTINE cbl_initial(qh_use, qe_use, tm_K_zm, qm_gkg_zm, startflag, ir, Gridiv)
363  USE mod_z
364  USE mod_k
365  USE gas
366  USE time
367  USE data_in
368  USE sues_data
369  USE moist
370  USE allocatearray
371  USE defaultnotused
372  USE cbl_module
373  USE gis_data
374  USE wherewhen
375  USE meteo, ONLY: sat_vap_press_x
376 
377  IMPLICIT NONE
378 
379  REAL(KIND(1d0))::qh_use, qe_use, tm_K_zm, qm_gkg_zm
380  REAL(KIND(1d0))::lv
381  INTEGER::i, nLineDay, ir, Gridiv, startflag
382 
383  qh_use = qhforcbl(gridiv) !HCW 21 Mar 2017
384  qe_use = qeforcbl(gridiv)
385  IF (qh_use < -900 .OR. qe_use < -900) THEN ! observed data has a problem
386  CALL errorhint(22, 'Unrealistic qh or qe_value for CBL in CBL_initial.', qh_use, qe_use, qh_choice)
387  ENDIF
388  !!Heat flux choices - these are now made in SUEWS_Calculations for qhforCBL and qeCBL, rather than here
389  !IF(Qh_choice==1) THEN !from SUEWS
390  ! !qh_use=qh
391  ! !qe_use=qeph
392  ! qh_use=qhforCBL(Gridiv) !HCW 21 Mar 2017
393  ! qe_use=qeforCBL(Gridiv)
394  !ELSEIF(qh_choice==2)THEN !from LUMPS
395  ! qh_use=H_mod
396  ! qe_use=E_mod
397  !ELSEIF(qh_choice==3)THEN !from OBS
398  ! IF(qh_obs<-900.OR.qe_obs<-900)THEN ! observed data has a problem
399  ! CALL ErrorHint(22,'Unrealistic observed qh or qe_value.',qh_obs,qe_obs,qh_choice)
400  ! ENDIF
401  ! qh_use=qh_obs
402  ! qe_use=qe_obs
403  !ENDIF
404 
405  blh_m = nan
406  ! iCBLcount = iCBLcount + 1
407  ! write(*,*) 'cblinitial', DateTIme, iCBLcount
408  dataoutbl(ir, 1:ncolumnsdataoutbl, gridiv) = (/REAL(iy, 8), REAL(id, 8), REAL(it, 8), REAL(imin, 8), dectime, &
409  (nan, is=6, ncolumnsdataoutbl)/)
410 
411  nlineday = 0
412  DO i = 1, nlineindata
413  IF (int(inicbldata(i, 1)) <= id) THEN
414  nlineday = nlineday + 1
415  ENDIF
416  ENDDO
417 
418  IF (initialdata_use == 2) THEN
419  blh_m = inicbldata(nlineday, 2)
420  gamt_km = inicbldata(nlineday, 3)
421  gamq_gkgm = inicbldata(nlineday, 4)
422  tp_k = inicbldata(nlineday, 5)
423  qp_gkg = inicbldata(nlineday, 6)
424  tm_k = inicbldata(nlineday, 7)
425  qm_gkg = inicbldata(nlineday, 8)
426  ELSEIF (initialdata_use == 1 .AND. inicbldata(nlineday, 1) == id) THEN ! Changed from i to nlineDay, HCW 29 March 2017
427  blh_m = inicbldata(nlineday, 2)
428  gamt_km = inicbldata(nlineday, 3)
429  gamq_gkgm = inicbldata(nlineday, 4)
430  tm_k_zm = (temp_c + c2k)*((1000/press_hpa)**(gas_ct_dry/avcp))
431  tm_k = tm_k_zm - psih*qh_use/(k*ustar*avcp*avdens)
432  es_hpa = sat_vap_press_x(temp_c, press_hpa, 1, dectime)
433  qm_gkg_zm = 622*avrh/(100*press_hpa/es_hpa - avrh)
434  lv = (2500.25 - 2.365*temp_c)*1000
435  qm_gkg = qm_gkg_zm - psih*qe_use/(k*ustar*avdens*lv)
436  tp_k = tm_k
437  qp_gkg = qm_gkg
438  ELSEIF (initialdata_use == 0) THEN
439  blh_m = 241.5
440  gamt_km = 0.043
441  gamq_gkgm = 0.0092
442  tm_k_zm = (temp_c + c2k)*((1000/press_hpa)**(gas_ct_dry/avcp))
443  tm_k = tm_k_zm - psih*qh_use/(k*ustar*avcp*avdens)
444  es_hpa = sat_vap_press_x(temp_c, press_hpa, 1, dectime)
445  qm_gkg_zm = 622*avrh/(100*press_hpa/es_hpa - avrh)
446  lv = (2500.25 - 2.365*temp_c)*1000
447  qm_gkg = es_hpa - psih*qe_use/(k*ustar*avdens*lv)
448  tp_k = tm_k
449  qp_gkg = qm_gkg
450  ENDIF
451 
452  gamq_kgkgm = gamq_gkgm/1000.
453  qp_kgkg = qp_gkg/1000 !humidities: g/kg -> kg/kg q+
454  qm_kgkg = qm_gkg/1000 !conc at mixing layer height h
455  tp_c = tp_k - c2k
456  tm_c = tm_k - c2k
457 
458  ! IF(sondeflag==1 .AND. cblday(id)==1) THEN
459  IF (sondeflag == 1 .AND. inicbldata(id, 2) /= -999) THEN
460  !if gamma theta varies with z (selected by setting gthetaflag=1)
461  !if gamma q varies with z (selected by setting ghumflag=1)
462  CALL sonde(id)
463  gamt_km = 0
464  gamq_kgkgm = 0
465  ENDIF
466 
467  !adjusting qp and pm in case of saturation
468  IF (qp_kgkg > qsatf(tp_c, press_hpa) .OR. qp_kgkg < 0) THEN
469  qp_kgkg = qsatf(tp_c, press_hpa)
470  ENDIF
471  IF (qm_kgkg > qsatf(tm_c, press_hpa) .OR. qm_kgkg < 0) THEN
472  qm_kgkg = qsatf(tm_c, press_hpa)
473  ENDIF
474 
475  ! if((CBLuse==2).and.(zenith_deg>=90))then
476  ! blh_m=188
477  ! endif
478  startflag = 0
479 
480  END SUBROUTINE cbl_initial
481 
482  SUBROUTINE nbl(iy, id, it, imin, dectime, ir, qh_choice, qh_use, qe_use, &
483  tm_K_zm, qm_gkg_zm, startflag, Gridiv, &
484  psih, UStar, Temp_C, NumberOfGrids, qhforCBL, qeforCBL, Press_hPa, avu1, avrh, &
485  ReadLinesMetdata, dataOutBL, &
486  avcp, avdens, es_hPa, lv_J_kg)
488  IMPLICIT NONE
489  INTEGER, PARAMETER:: ncolumnsdataoutbl = 22
490 
491  INTEGER, INTENT(IN) :: qh_choice, iy, id, it, imin, NumberOfGrids, ReadLinesMetdata, ir
492  REAL(KIND(1d0)), INTENT(IN):: Press_hPa, psih, UStar
493  REAL(KIND(1d0)), INTENT(OUT):: Temp_C, tm_K_zm, qm_gkg_zm
494  REAL(KIND(1d0)), INTENT(INOUT):: dectime, avu1, avRH, avcp, avdens, es_hPa, lv_J_kg
495  REAL(KIND(1d0)), INTENT(IN), DIMENSION(NumberOfGrids):: qhforCBL, qeforCBL
496  REAL(KIND(1d0)), INTENT(OUT), DIMENSION(ReadLinesMetdata, ncolumnsdataOutBL, NumberOfGrids) ::dataOutBL
497 
498  REAL(KIND(1d0)) :: &
499  k = 0.4, & !Von Karman's contant
500  gas_ct_dry = 8.31451/0.028965 !j/kg/k=dry_gas/molar
501  REAL(KIND(1d0))::qh_use, qe_use
502  REAL(KIND(1d0))::lv
503  INTEGER::i, nLineDay, Gridiv, startflag
504 
505  qh_use = qhforcbl(gridiv) !HCW 21 Mar 2017
506  qe_use = qeforcbl(gridiv)
507  IF (qh_use < -900 .OR. qe_use < -900) THEN ! observed data has a problem
508  CALL errorhint(22, 'Unrealistic qh or qe value for CBL in NBL.', qh_use, qe_use, qh_choice)
509  ENDIF
510 
511  nlineday = 0
512  DO i = 1, nlineindata
513  IF (int(inicbldata(i, 1)) <= id) THEN
514  nlineday = nlineday + 1
515  ENDIF
516  ENDDO
517 
518  !Assume Theta and Q in the night constantly Equal to the ones in the morning for CBL to run
519  ! seems incorrect initialisation
520  ! tm_K=IniCBLdata(nLineDay,7)*1000
521  ! qm_gkg=IniCBLdata(nLineDay,8)*1000
522 
523  ! corrected to the following, TS 20170609:
524  tm_k = inicbldata(nlineday, 7)
525  qm_gkg = inicbldata(nlineday, 8)
526 
527  !NBL currently fixed to 200 m
528  blh_m = 200
529 
530  ! also fill in other variables
531  gamt_km = inicbldata(nlineday, 3)
532  gamq_gkgm = inicbldata(nlineday, 4)
533  tp_k = inicbldata(nlineday, 5)
534  qp_gkg = inicbldata(nlineday, 6)
535  tm_k = inicbldata(nlineday, 7)
536  qm_gkg = inicbldata(nlineday, 8)
537 
538  ! dataOutBL(iCBLcount,1:ncolumnsdataOutBL,Gridiv)=(/REAL(iy,8),REAL(id,8),REAL(it,8),REAL(imin,8),&
539  ! dectime,blh_m,tm_K,qm_gkg,&
540  ! (NAN,is=9,ncolumnsdataOutBL)/)
541 
542  temp_c = tm_k/((1000/press_hpa)**(gas_ct_dry/avcp)) - c2k
543  es_hpa = sat_vap_press_x(temp_c, press_hpa, 1, dectime)
544  lv = (2500.25 - 2.365*temp_c)*1000
545  !qm_gkg_zm=qm_gkg+cbldata(10)*cbldata(3)/(k*cbldata(8)*cbldata(4)*lv)
546  avrh = 100*((qm_gkg*press_hpa/(622 + qm_gkg))/es_hpa) !check pressure
547  IF (avrh > 100) THEN
548  CALL errorhint(34, 'subroutine CBL dectime, relative humidity', dectime, avrh, 100)
549  avrh = 100
550  ENDIF
551 
552  ! print *, 'iCBLcount in NBL',iCBLcount
553  ! print *,iy, id,it,imin
554 
555  dataoutbl(ir, 1:ncolumnsdataoutbl, gridiv) = &
556  (/REAL(iy, 8), REAL(id, 8), REAL(it, 8), REAL(imin, 8), dectime, &
557  blh_m, tm_K, qm_gkg, &
558  tp_K, qp_gkg, &
559  Temp_C, avrh, qh_use, qe_use, Press_hPa, avu1, UStar, avdens, lv_J_kg, avcp, &
560  gamt_Km, gamq_kgkgm/)
561 
562  IF (initialdata_use == 2) THEN
563  blh_m = inicbldata(nlineday, 2)
564  gamt_km = inicbldata(nlineday, 3)
565  gamq_gkgm = inicbldata(nlineday, 4)
566  tp_k = inicbldata(nlineday, 5)
567  qp_gkg = inicbldata(nlineday, 6)
568  tm_k = inicbldata(nlineday, 7)
569  qm_gkg = inicbldata(nlineday, 8)
570  ELSEIF (initialdata_use == 1 .AND. inicbldata(nlineday, 1) == id) THEN ! Changed from i to nlineDay, HCW 29 March 2017
571  blh_m = inicbldata(nlineday, 2)
572  gamt_km = inicbldata(nlineday, 3)
573  gamq_gkgm = inicbldata(nlineday, 4)
574  tm_k_zm = (temp_c + c2k)*((1000/press_hpa)**(gas_ct_dry/avcp))
575  tm_k = tm_k_zm - psih*qh_use/(k*ustar*avcp*avdens)
576  es_hpa = sat_vap_press_x(temp_c, press_hpa, 1, dectime)
577  qm_gkg_zm = 622*avrh/(100*press_hpa/es_hpa - avrh)
578  lv = (2500.25 - 2.365*temp_c)*1000
579  qm_gkg = qm_gkg_zm - psih*qe_use/(k*ustar*avdens*lv)
580  tp_k = tm_k
581  qp_gkg = qm_gkg
582  ELSEIF (initialdata_use == 0) THEN
583  blh_m = 241.5
584  gamt_km = 0.043
585  gamq_gkgm = 0.0092
586  tm_k_zm = (temp_c + c2k)*((1000/press_hpa)**(gas_ct_dry/avcp))
587  tm_k = tm_k_zm - psih*qh_use/(k*ustar*avcp*avdens)
588  es_hpa = sat_vap_press_x(temp_c, press_hpa, 1, dectime)
589  qm_gkg_zm = 622*avrh/(100*press_hpa/es_hpa - avrh)
590  lv = (2500.25 - 2.365*temp_c)*1000
591  qm_gkg = es_hpa - psih*qe_use/(k*ustar*avdens*lv)
592  tp_k = tm_k
593  qp_gkg = qm_gkg
594  ENDIF
595 
596  gamq_kgkgm = gamq_gkgm/1000.
597  qp_kgkg = qp_gkg/1000 !humidities: g/kg -> kg/kg q+
598  qm_kgkg = qm_gkg/1000 !conc at mixing layer height h
599  tp_c = tp_k - c2k
600  tm_c = tm_k - c2k
601 
602  ! IF(sondeflag==1 .AND. cblday(id)==1) THEN
603  IF (sondeflag == 1 .AND. inicbldata(id, 2) /= -999) THEN
604  !if gamma theta varies with z (selected by setting gthetaflag=1)
605  !if gamma q varies with z (selected by setting ghumflag=1)
606  CALL sonde(id)
607  gamt_km = 0
608  gamq_kgkgm = 0
609  ENDIF
610 
611  !adjusting qp and pm in case of saturation
612  IF (qp_kgkg > qsatf(tp_c, press_hpa) .OR. qp_kgkg < 0) THEN
613  qp_kgkg = qsatf(tp_c, press_hpa)
614  ENDIF
615  IF (qm_kgkg > qsatf(tm_c, press_hpa) .OR. qm_kgkg < 0) THEN
616  qm_kgkg = qsatf(tm_c, press_hpa)
617  ENDIF
618 
619  startflag = 0
620 
621  ! iCBLcount = iCBLcount + 1 ! move counter at the end, TS 27 Aug 2019
622  END SUBROUTINE nbl
623 
624  !------------------------------------------------------------------------
625 
626  !-----------------------------------------------------------------------
627  ! from CBL modelling Cleugh and Grimmond (2000) BLM
628  ! NT 6 Apr 2017: include iteration over top of CBL scalars and include subsidence flag
629  ! Last modified: LJ 27 Jan 2016 - Removal of tabs
630  !-----------------------------------------------------------------------
631  SUBROUTINE rkutta(neqn_use, XA, XB, y_use, NSTEPS)
632  ! XA=s0
633  ! XB=s1
634  ! Y(1)=blh_m
635  ! Y(2)=tm_K
636  ! Y(3)=qm_kgkg
637  ! Y(4)=cm
638  ! Y(5)=tp_K
639  ! Y(6)=qp_kgkg
640  ! JOHN KNIGHT, 1985 (AMENDED BY MRR, 23-SEP-85)
641  ! EXPLICIT FOURTH-ORDER RUNGE-KUTTA METHOD FOR FIRST-ORDER ODE SYSTEM
642  ! OF NE EQUATIONS, WITH INITIAL VALUES SUPPLIED.
643  ! MEANING OF PARAMETERS:
644  ! NE = NUMBER OF EQUATIONS (MAX 21)
645  ! XA = VALUE OF INDEPENDENT VARIABLE ON ENTRY
646  ! XB = VALUE OF INDEPENDENT VARIABLE AT END OF INTERVAL
647  ! Y(NE) = ON ENTRY: INITIAL VALUES OF DEPENDENT VARIABLES AT XA
648  ! ON EXIT: CALCULATED VALUES OF DEPENDENT VARIABLES AT XB
649  ! NSTEPS = NUMBER OF INTEGRATION STEPS OVER INTERVAL (XA,XB)
650  ! DIFF = NAME OF USER-SUPPLIED SUBROUTINE TO CALCULATE DERIVATIVES
651  ! DYDX (DIFF MUST BE DECLARED EXTERNAL IN CALLING PROGRAM).
652  ! PARAMETERS IN SUBROUTINE DIFF(NE,X,Y,DYDX):
653  ! NEqn = NUMBER OF EQUATIONS
654  ! X = INDEPENDENT VARIABLE
655  ! Y = ARRAY (LENGTH NE) OF VALUES OF DEPENDENT VARIABLES
656  ! DYDX = ON EXIT, ARRAY (LENGTH NE) OF VALUES OF DERIVATIVES
657  ! IMPLICIT real*8 (A-H,O-Z)
658  IMPLICIT NONE
659  INTEGER::ns, nsteps, nj, n, neqn_use
660  REAL(KIND(1D0)), DIMENSION(neqn_use):: y_use
661  REAL(KIND(1D0)), DIMENSION(21):: dydx, arg
662  REAL(KIND(1D0)), DIMENSION(21, 5):: rk
663  REAL(KIND(1D0)), DIMENSION(4):: coef
664  REAL(KIND(1D0)):: XA, XB, step, X, xx
665 
666  coef(1) = 1.0
667  coef(2) = 0.5
668  coef(3) = 0.5
669  coef(4) = 1.0
670  ! print*,"rk1: ",xa,xb,y
671  step = (xb - xa)/nsteps
672 
673  DO ns = 1, nsteps
674  DO nj = 1, neqn_use
675  rk(nj, 1) = 0
676  ENDDO
677  x = xa + (ns - 1)*step
678  DO n = 1, 4
679  IF (n == 1) THEN
680  xx = x
681  ELSEIF (n > 1) THEN
682  xx = x + coef(n)*step
683  ENDIF
684 
685  DO nj = 1, neqn_use
686  arg(nj) = y_use(nj) + coef(n)*rk(nj, n)
687  ENDDO
688 
689  CALL diff(xx, arg, dydx)
690 
691  DO nj = 1, neqn_use
692  rk(nj, n + 1) = step*dydx(nj)
693  ENDDO
694  ENDDO
695 
696  DO nj = 1, neqn_use
697  DO n = 1, 4
698  y_use(nj) = y_use(nj) + rk(nj, n + 1)/(6*coef(n))
699  ENDDO
700  ENDDO
701  ENDDO
702 
703  RETURN
704  END SUBROUTINE rkutta
705  !---------------------------------------------------------------------
706  !---------------------------------------------------------------------
707 
708  SUBROUTINE diff(s, y1, dyds)
709  ! in y1,neqn
710  ! out dyds
711 
712  ! calculates derivatives for cbl slab model
713  ! y(1) = h = cbl depth(m)
714  ! y(2) = t = potential temp(K)
715  ! y(3) = q = specific humidity(kg/kg)
716  ! y(4) = c = CO2 concentration
717  ! USE data_in
718  ! USE sues_data
719  ! use allocateArray
720  USE time
721  USE cbl_module
722  USE defaultnotused
723  USE mod_grav
724 
725  IMPLICIT NONE
726  REAL(KIND(1D0)), DIMENSION(neqn)::dyds, y1
727  REAL(KIND(1d0)) :: zero = 0.0
728  REAL(KIND(1d0)) :: h1, t_K, q_kgkg, c, cp, ws, s, foo
729  ! real(kind(1D0)) :: tp_K,qp_kgkg
730  REAL(KIND(1D0)):: delt_K, delq_kgkg, delc
731  REAL(KIND(1D0)):: gamtv_Km, deltv_K, ftv_Kms
732  REAL(KIND(1D0)):: ftva_Kms, delb, qs2, qs3
733  REAL(KIND(1D0)):: dhds, dtds, dqds, dcds, dtpds, dqpds
734  REAL(KIND(1D0)):: conk, conn, cona, conc, cont
735 
736  ! print*,"diff: timestamp:",s
737  foo = s
738  ! pause
739  h1 = y1(1)!m
740  t_k = y1(2)!K
741  q_kgkg = y1(3)!kg/kg
742  c = y1(4)
743  tp_k = y1(5)!K
744  qp_kgkg = y1(6)!kg/kg
745 
746  ! find t, q, c above inversion, and jumps across inversion
747  ! tp = tp + gamt*h
748  ! qp = qp0 + gamq*h
749 
750  cp = 0 ! cp0 + gamc* h1 ! todo
751 
752  delt_k = tp_k - t_k
753  delq_kgkg = qp_kgkg - q_kgkg
754  delc = cp - c
755 
756  ! find potential virtual temperature flux, gradient and jump
757  ftv_kms = fhbl_kms + 0.61*tm_k*febl_kgkgms
758  gamtv_km = gamt_km + 0.61*tm_k*gamq_kgkgm!/1000
759  deltv_k = delt_k + 0.61*tm_k*delq_kgkg
760 
761  ! find velocity scale ws
762  ftva_kms = max(ftv_kms, zero) ! virtual heat flux
763  ws = (h1*ftva_kms*grav/tm_k)**0.3333333333
764 
765  ! find dhds using one of 4 alternative schemes chosen by ient:
766  IF (entrainmenttype == 2) THEN
767  ! EntrainmentType=1: encroachment (as in McN and S 1986 eq 16))
768  dhds = ftva_kms/(h1*gamtv_km)
769 
770  ELSE IF (entrainmenttype == 1) THEN
771  ! EntrainmentType=2: Driedonks 1981 (as in McN and S 1986 eq 13)
772  IF (deltv_k <= 0.01) THEN
773  dhds = ftva_kms/(h1*gamtv_km)
774  CALL errorhint(30, "subroutine diff [CBL: Deltv_K<0.01 EntrainmentType=1], deltv_K,delt_K,", deltv_k, delt_k, notusedi)
775  CALL errorhint(30, "subroutine diff [CBL: Deltv_K<0.01 EntrainmentType=1], tm_K,TPP_K,y1", tm_k, tpp_k, notusedi)
776  ! call errorHint(31,"subroutine diff [CBL: Deltv_K<0.01 EntrainmentType=1], y1",real(y1(1),kind(1d0)),notUsed,notUsedI)
777  ELSE
778  delb = grav*deltv_k/tm_k
779  conc = 0.2
780  cona = 5.0
781  dhds = (conc*ws**3 + cona*cbldata(8)**3)/(h1*delb)
782  END IF
783 
784  ELSE IF (entrainmenttype == 4) THEN
785  ! EntrainmentType=3: Tennekes 1973 (as in R 1991 eqs 3,4)
786  alpha3 = 0.2 ! alpha changed back to original Tennekes 1973 value
787  IF (deltv_k <= 0.01) THEN
788  dhds = ftva_kms/(h1*gamtv_km)
789  CALL errorhint(31, 'subroutine difflfnout: [CBL: deltv_K<0.01 EntrainmentType=4],deltv_K', &
790  deltv_k, notused, notusedi)
791  ELSE
792  ! include the option whether or not to include subsidence
793  IF (isubs == 1) THEN
794  dhds = alpha3*ftva_kms/deltv_k + wsb
795  ELSE
796  dhds = alpha3*ftva_kms/deltv_k
797  END IF
798  END IF
799 
800  ! write (4,*) tpp, gamq, dhds, deltv
801 
802  ELSE IF (entrainmenttype == 3) THEN
803  ! EntrainmentType=4: Rayner and Watson 1991 eq 21
804  conn = 1.33
805  conk = 0.18
806  cont = 0.80
807  qs3 = ws**3 + (conn*cbldata(8))**3
808  qs2 = qs3**(0.6666666667)
809 
810  IF (deltv_k <= 0.01) THEN
811  dhds = ftva_kms/(h1*gamtv_km)
812  CALL errorhint(31, 'subroutine difflfnout: [CBL: deltv_K<0.01 EntrainmentType=3],deltv_K', &
813  deltv_k, notused, notusedi)
814 
815  ELSE
816  delb = grav*deltv_k/tm_k
817  dhds = (conk*qs3)/(cont*qs2 + h1*delb)
818  END IF
819 
820  ELSE
821  CALL errorhint(24, 'BLUEWS_DIff- CBL- illegal alpha', notused, notused, notusedi)
822  END IF
823  ! find dtds, dqds, dc/ds:
824  ! wsb is the subsidence velocity. Try using: -0.01, -0.05, -0.1.
825  IF (isubs == 1) THEN
826  dtds = fhbl_kms/h1 + delt_k*(dhds - wsb)/h1
827  dqds = febl_kgkgms/h1 + delq_kgkg*(dhds - wsb)/h1
828  dcds = fcbl/h1 + delc*(dhds - wsb)/h1
829  ! also iterate the top of CBL scalars
830  dtpds = gamt_km*(dhds - wsb)
831  dqpds = gamq_kgkgm*(dhds - wsb)
832  ELSE
833  dtds = fhbl_kms/h1 + delt_k*(dhds)/h1
834  dqds = febl_kgkgms/h1 + delq_kgkg*(dhds)/h1
835  dcds = fcbl/h1 + delc*(dhds)/h1
836  ! also iterate the top of CBL scalars
837  dtpds = gamt_km*(dhds)
838  dqpds = gamq_kgkgm*(dhds)
839  END IF
840 
841  dyds(1) = dhds
842  dyds(2) = dtds
843  dyds(3) = dqds
844  dyds(4) = dcds
845  dyds(5) = dtpds
846  dyds(6) = dqpds
847 
848  RETURN
849  END SUBROUTINE diff
850 
851  !--------------------------------------------------------------------------
852  !--------------------------------------------------------------------------
853  SUBROUTINE sonde(id)
854  ! read sonde or vertical profile data - when available
855  !use allocateArray
856  USE data_in
857  USE cbl_module
858  IMPLICIT NONE
859  INTEGER::i, fn = 101, izm = 500, notusedi = -9999, id
860  CHARACTER(len=200)::FileN
861  REAL(KIND(1d0)):: dxx
862  REAL(KIND(1d0)), PARAMETER::notUsed = -9999.99
863 
864  filen = trim(fileinputpath)//trim(filesonde(id))
865  OPEN (fn, file=filen, status="old", err=24)
866  ! todo gneric skip header
867  READ (fn, *)
868  READ (fn, *)
869  READ (fn, *)
870 
871  DO i = 1, 1000
872  READ (fn, *, end=900, err=25) gtheta(i, 1), dxx, gtheta(i, 2), ghum(i, 1), dxx, ghum(i, 2)
873  ghum(i, 2) = ghum(i, 2)
874  ENDDO
875 900 zmax = i - 1
876  IF (zmax > izm) THEN
877  CALL errorhint(23, filen, REAL(zmax, KIND(1D0)), notUsed, izm)
878  ENDIF
879  CLOSE (fn)
880  RETURN
881 24 CALL errorhint(24, filen, notused, notused, notusedi)
882 25 CALL errorhint(25, filen, notused, notused, i)
883  RETURN
884  END SUBROUTINE sonde
885  !------------------------------------------------------------------------------
886  !------------------------------------------------------------------------------
887  SUBROUTINE gamma_sonde
889  !use allocateArray
890 
891  IMPLICIT NONE
892  REAL(KIND(1D0))::gamtt, gamqq
893  INTEGER::j
894  ! gtheta(i,1),dxx,gtheta(i,2),ghum(i,1),dxx,ghum(i,2)
895  !search for correct gamma theta, depends on h(i-1),
896  ! ie current value for mixed layer depth
897  IF (sondeflag == 1) THEN
898  DO j = 2, zmax
899  IF (blh_m >= gtheta(j - 1, 1)) THEN
900  gamtt = gtheta(j - 1, 2)
901  ENDIF
902  gamt_km = gamtt
903  ENDDO
904 
905  DO j = 2, zmax
906  IF (blh_m >= ghum(j - 1, 1)) THEN
907  gamqq = ghum(j - 1, 2)
908  ENDIF
909  gamq_kgkgm = gamqq/1000.
910  ENDDO
911  ENDIF
912  RETURN
913 
914  END SUBROUTINE gamma_sonde
915 
916 END MODULE bluews_module
real(kind(1d0)) usbl
real(kind(1d0)), dimension(1:10) cbldata
subroutine gamma_sonde
character(len=200) initialdatafilename
real(kind(1d0)) qm_kgkg
real(kind(1d0)) k
real(kind(1d0)) ustar
real(kind(1d0)) fhbl_kms
subroutine nbl(iy, id, it, imin, dectime, ir, qh_choice, qh_use, qe_use, tm_K_zm, qm_gkg_zm, startflag, Gridiv, psih, UStar, Temp_C, NumberOfGrids, qhforCBL, qeforCBL, Press_hPa, avu1, avrh, ReadLinesMetdata, dataOutBL, avcp, avdens, es_hPa, lv_J_kg)
real(kind(1d0)), dimension(0:500, 2) ghum
real(kind(1d0)) tmp_k
subroutine cbl_readinputdata(FileInputPath, qh_choice)
integer, dimension(366) cblday
real(kind(1d0)) avrh
real(kind(1d0)), dimension(:, :, :), allocatable dataoutbl
real(kind(1d0)) nan
subroutine diff(s, y1, dyds)
real(kind(1d0)) temp_c
real(kind(1d0)) notused
real(kind(1d0)) gamt_km
real(kind(1d0)) qpp
real(kind(1d0)) fqbl
integer entrainmenttype
real(kind(1d0)) qp_gkg
integer, parameter ncolumnsdataoutbl
real(kind(1d0)) tpp
real(kind(1d0)), dimension(6) y
real(kind(1d0)) qm_gkg
real(kind(1d0)) cm
real(kind(1d0)) gamq
integer id
real(kind(1d0)) alpha3
real(kind(1d0)) avcp
real(kind(1d0)) tm_c
real(kind(1d0)) tp_c
real(kind(1d0)) gas_ct_dry
character(len=150) fileinputpath
subroutine cbl_initial(qh_use, qe_use, tm_K_zm, qm_gkg_zm, startflag, ir, Gridiv)
subroutine cbl(iy, id, it, imin, ir, Gridiv, qh_choice, dectime, Temp_C, Press_hPa, avkdn, avu1, avrh, avcp, avdens, es_hPa, lv_J_kg, nsh_real, tstep, UStar, psih, is, NumberOfGrids, qhforCBL, qeforCBL, ReadLinesMetdata, dataOutBL)
real(kind(1d0)) psih
integer initialdata_use
real(kind(1d0)) wsb
real(kind(1d0)) gamc
real(kind(1d0)) c2k
real(kind(1d0)), dimension(:), allocatable qeforcbl
character(len=200), dimension(366) filesonde
subroutine rkutta(neqn_use, XA, XB, y_use, NSTEPS)
real(kind(1d0)) tm_k
real(kind(1d0)), dimension(0:500, 2) gtheta
real(kind(1d0)) gamt
real(kind(1d0)) gamq_gkgm
real(kind(1d0)) tpp_k
real(kind(1d0)), dimension(:), allocatable qhforcbl
real(kind(1d0)) blh_m
real(kind(1d0)) es_hpa
integer co2_included
real(kind(1d0)) function sat_vap_press_x(Temp_c, PRESS_hPa, from, dectime)
real(kind(1d0)) cp0
real(kind(1d0)) avdens
real(kind(1d0)) ftbl
real(kind(1d0)) tp_k
real(kind(1d0)) qpp_kgkg
subroutine errorhint(errh, ProblemFile, VALUE, value2, valueI)
subroutine sonde(id)
real(kind(1d0)) febl_kgkgms
real(kind(1d0)) qp_kgkg
real(kind(1d0)) function qsatf(T, PMB)
real(kind(1d0)), dimension(:, :), allocatable inicbldata
real(kind(1d0)) grav
real(kind(1d0)) gamq_kgkgm
real(kind(1d0)) press_hpa
real(kind(1d0)) blh1_m
real(kind(1d0)) fcbl