SUEWS API Site
Documentation of SUEWS source code
Functions/Subroutines
bluews_module Module Reference

Functions/Subroutines

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)
 
subroutine cbl_readinputdata (fileinputpath, qh_choice)
 
subroutine cbl_initial (qh_use, qe_use, tm_k_zm, qm_gkg_zm, startflag, ir, gridiv)
 
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)
 
subroutine rkutta (neqn_use, xa, xb, y_use, nsteps)
 
subroutine diff (s, y1, dyds)
 
subroutine sonde (id)
 
subroutine gamma_sonde
 

Function/Subroutine Documentation

◆ cbl()

subroutine bluews_module::cbl ( integer, intent(inout) iy,
integer, intent(inout) id,
integer, intent(inout) it,
integer, intent(inout) imin,
integer, intent(in) ir,
integer, intent(in) gridiv,
integer, intent(inout) qh_choice,
real(kind(1d0)), intent(inout) dectime,
real(kind(1d0)), intent(out) temp_c,
real(kind(1d0)), intent(inout) press_hpa,
real(kind(1d0)), intent(in) avkdn,
real(kind(1d0)), intent(inout) avu1,
real(kind(1d0)), intent(inout) avrh,
real(kind(1d0)), intent(inout) avcp,
real(kind(1d0)), intent(inout) avdens,
real(kind(1d0)), intent(inout) es_hpa,
real(kind(1d0)), intent(inout) lv_j_kg,
real(kind(1d0)), intent(in) nsh_real,
integer, intent(in) tstep,
real(kind(1d0)), intent(in) ustar,
real(kind(1d0)), intent(in) psih,
integer, intent(in) is,
integer, intent(in) numberofgrids,
real(kind(1d0)), dimension(numberofgrids), intent(in) qhforcbl,
real(kind(1d0)), dimension(numberofgrids), intent(in) qeforcbl,
integer, intent(in) readlinesmetdata,
real(kind(1d0)), dimension(readlinesmetdata, ncolumnsdataoutbl, numberofgrids), intent(out) dataoutbl )

Definition at line 69 of file suews_phys_bluews.f95.

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 END IF
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 END IF
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 END IF
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 END IF
200
201 ! tpp_K=tp_K
202 ! qpp_kgkg=qp_kgkg
203
204 IF (sondeflag == 1) THEN
205 CALL gamma_sonde
206 END IF
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 END IF
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 END IF
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 END IF
291 ! move the counter at the end, TS 27 Aug 2019
292 ! iCBLcount = iCBLcount + 1
293
294 RETURN
295
errorhint
subroutine errorhint(errh, problemfile, value, value2, valuei)
Definition suews_ctrl_error.f95:2

References cbl_module::blh1_m, cbl_module::blh_m, cbl_module::c2k, cbl_module::cbldata, cbl_module::cm, cbl_module::co2_included, errorhint(), cbl_module::fcbl, cbl_module::febl_kgkgms, cbl_module::fhbl_kms, gamma_sonde(), cbl_module::gamq_kgkgm, cbl_module::gamt_km, nbl(), cbl_module::neqn, cbl_module::qm_gkg, cbl_module::qm_kgkg, cbl_module::qp_kgkg, rkutta(), meteo::sat_vap_press_x(), cbl_module::sondeflag, cbl_module::tm_c, cbl_module::tm_k, cbl_module::tp_c, cbl_module::tp_k, and cbl_module::y.

Here is the call graph for this function:

◆ cbl_initial()

subroutine bluews_module::cbl_initial ( real(kind(1d0)) qh_use,
real(kind(1d0)) qe_use,
real(kind(1d0)) tm_k_zm,
real(kind(1d0)) qm_gkg_zm,
integer startflag,
integer ir,
integer gridiv )

Definition at line 361 of file suews_phys_bluews.f95.

362
363 USE mod_z
364 USE atmmoiststab_module, ONLY: 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 END IF
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 END IF
416 END DO
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 END IF
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 END IF
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 END IF
471 IF (qm_kgkg > qsatf(tm_c, press_hpa) .OR. qm_kgkg < 0) THEN
472 qm_kgkg = qsatf(tm_c, press_hpa)
473 END IF
474
475 ! if((CBLuse==2).and.(zenith_deg>=90))then
476 ! blh_m=188
477 ! endif
478 startflag = 0
479
allocatearray::ncolumnsdataoutbl
integer, parameter ncolumnsdataoutbl
Definition suews_ctrl_const.f95:80
allocatearray::dataoutbl
real(kind(1d0)), dimension(:, :, :), allocatable dataoutbl
Definition suews_ctrl_const.f95:158
atmmoiststab_module
Definition suews_phys_atmmoiststab.f95:1
atmmoiststab_module::k
real(kind(1d0)), parameter k
Definition suews_phys_atmmoiststab.f95:4
cbl_module
Definition suews_phys_bluews.f95:2
cbl_module::qp_kgkg
real(kind(1d0)) qp_kgkg
Definition suews_phys_bluews.f95:28
cbl_module::gamq_kgkgm
real(kind(1d0)) gamq_kgkgm
Definition suews_phys_bluews.f95:28
cbl_module::blh_m
real(kind(1d0)) blh_m
Definition suews_phys_bluews.f95:28
cbl_module::c2k
real(kind(1d0)) c2k
Definition suews_phys_bluews.f95:24
cbl_module::qp_gkg
real(kind(1d0)) qp_gkg
Definition suews_phys_bluews.f95:28
cbl_module::inicbldata
real(kind(1d0)), dimension(:, :), allocatable inicbldata
Definition suews_phys_bluews.f95:17
cbl_module::sondeflag
integer sondeflag
Definition suews_phys_bluews.f95:4
cbl_module::nlineindata
integer nlineindata
Definition suews_phys_bluews.f95:20
cbl_module::qm_gkg
real(kind(1d0)) qm_gkg
Definition suews_phys_bluews.f95:28
cbl_module::gamq_gkgm
real(kind(1d0)) gamq_gkgm
Definition suews_phys_bluews.f95:28
cbl_module::tp_c
real(kind(1d0)) tp_c
Definition suews_phys_bluews.f95:28
cbl_module::gamt_km
real(kind(1d0)) gamt_km
Definition suews_phys_bluews.f95:28
cbl_module::tm_c
real(kind(1d0)) tm_c
Definition suews_phys_bluews.f95:28
cbl_module::qm_kgkg
real(kind(1d0)) qm_kgkg
Definition suews_phys_bluews.f95:28
cbl_module::initialdata_use
integer initialdata_use
Definition suews_phys_bluews.f95:4
cbl_module::tp_k
real(kind(1d0)) tp_k
Definition suews_phys_bluews.f95:28
cbl_module::tm_k
real(kind(1d0)) tm_k
Definition suews_phys_bluews.f95:28
data_in
Definition suews_ctrl_const.f95:1095
data_in::press_hpa
real(kind(1d0)) press_hpa
Definition suews_ctrl_const.f95:1192
data_in::avrh
real(kind(1d0)) avrh
Definition suews_ctrl_const.f95:1192
data_in::temp_c
real(kind(1d0)) temp_c
Definition suews_ctrl_const.f95:1192
defaultnotused::nan
real(kind(1d0)) nan
Definition suews_ctrl_const.f95:1394
gas::gas_ct_dry
real(kind(1d0)) gas_ct_dry
Definition suews_ctrl_const.f95:1443
gis_data
Definition suews_ctrl_const.f95:1495
meteo::sat_vap_press_x
real(kind(1d0)) function sat_vap_press_x(temp_c, press_hpa, from, dectime)
Definition suews_util_meteo.f95:149
meteo::qsatf
real(kind(1d0)) function qsatf(t, pmb)
Definition suews_util_meteo.f95:496
mod_z
Definition suews_ctrl_const.f95:1448
moist::avcp
real(kind(1d0)) avcp
Definition suews_ctrl_const.f95:1475
moist::avdens
real(kind(1d0)) avdens
Definition suews_ctrl_const.f95:1475
moist::es_hpa
real(kind(1d0)) es_hpa
Definition suews_ctrl_const.f95:1475
sues_data
Definition suews_ctrl_const.f95:1526
sues_data::qeforcbl
real(kind(1d0)), dimension(:), allocatable qeforcbl
Definition suews_ctrl_const.f95:1659
sues_data::qh_choice
integer qh_choice
Definition suews_ctrl_const.f95:1660
sues_data::ustar
real(kind(1d0)) ustar
Definition suews_ctrl_const.f95:1625
sues_data::is
integer is
Definition suews_ctrl_const.f95:1548
sues_data::qhforcbl
real(kind(1d0)), dimension(:), allocatable qhforcbl
Definition suews_ctrl_const.f95:1659
sues_data::psih
real(kind(1d0)) psih
Definition suews_ctrl_const.f95:1625
time::dectime
real(kind(1d0)) dectime
Definition suews_ctrl_const.f95:1409

References moist::avcp, moist::avdens, data_in::avrh, cbl_module::blh_m, cbl_module::c2k, allocatearray::dataoutbl, time::dectime, errorhint(), moist::es_hpa, cbl_module::gamq_gkgm, cbl_module::gamq_kgkgm, cbl_module::gamt_km, gas::gas_ct_dry, time::id, time::imin, cbl_module::inicbldata, cbl_module::initialdata_use, sues_data::is, time::it, time::iy, atmmoiststab_module::k, defaultnotused::nan, allocatearray::ncolumnsdataoutbl, cbl_module::nlineindata, data_in::press_hpa, sues_data::psih, sues_data::qeforcbl, sues_data::qh_choice, sues_data::qhforcbl, cbl_module::qm_gkg, cbl_module::qm_kgkg, cbl_module::qp_gkg, cbl_module::qp_kgkg, meteo::qsatf(), meteo::sat_vap_press_x(), sonde(), cbl_module::sondeflag, data_in::temp_c, cbl_module::tm_c, cbl_module::tm_k, cbl_module::tp_c, cbl_module::tp_k, and sues_data::ustar.

Here is the call graph for this function:

◆ cbl_readinputdata()

subroutine bluews_module::cbl_readinputdata ( character(len=150), intent(in) fileinputpath,
integer, intent(inout) qh_choice )

Definition at line 300 of file suews_phys_bluews.f95.

301
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, &
316 initialdatafilename, &
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
332 nlineindata = nlineindata + 1
333 END DO
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 END DO
343 CLOSE (52)
344 END IF
345
346 IF (co2_included == 0) THEN
347 fcbl = 0 ! hard-wire no CO2
348 END IF
349
350 ! iCBLcount = 1
351
352 RETURN
353
35424 CALL errorhint(24, 'CBLInput.nml', 0.00d0, 0.000d0, 0)
35525 CALL errorhint(24, trim(fileinputpath)//trim(initialdatafilename), 0.00d0, 0.00d0, 0)
356

References cbl_module::cblday, cbl_module::co2_included, cbl_module::entrainmenttype, errorhint(), cbl_module::fcbl, cbl_module::filesonde, cbl_module::inicbldata, cbl_module::initialdata_use, cbl_module::initialdatafilename, cbl_module::isubs, cbl_module::nlineindata, cbl_module::sondeflag, and cbl_module::wsb.

Here is the call graph for this function:

◆ diff()

subroutine bluews_module::diff ( real(kind(1d0)) s,
real(kind(1d0)), dimension(neqn) y1,
real(kind(1d0)), dimension(neqn) dyds )

Definition at line 708 of file suews_phys_bluews.f95.

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
cbl_module::isubs
integer isubs
Definition suews_phys_bluews.f95:4
cbl_module::fhbl_kms
real(kind(1d0)) fhbl_kms
Definition suews_phys_bluews.f95:28
cbl_module::febl_kgkgms
real(kind(1d0)) febl_kgkgms
Definition suews_phys_bluews.f95:28
cbl_module::tpp_k
real(kind(1d0)) tpp_k
Definition suews_phys_bluews.f95:28
cbl_module::entrainmenttype
integer entrainmenttype
Definition suews_phys_bluews.f95:4
cbl_module::fcbl
real(kind(1d0)) fcbl
Definition suews_phys_bluews.f95:26
cbl_module::cbldata
real(kind(1d0)), dimension(1:10) cbldata
Definition suews_phys_bluews.f95:16
cbl_module::alpha3
real(kind(1d0)) alpha3
Definition suews_phys_bluews.f95:28
cbl_module::wsb
real(kind(1d0)) wsb
Definition suews_phys_bluews.f95:15
defaultnotused::notused
real(kind(1d0)) notused
Definition suews_ctrl_const.f95:1394
mod_grav
Definition suews_ctrl_const.f95:1420
mod_grav::grav
real(kind(1d0)) grav
Definition suews_ctrl_const.f95:1421

References cbl_module::alpha3, cbl_module::cbldata, cbl_module::entrainmenttype, errorhint(), cbl_module::fcbl, cbl_module::febl_kgkgms, cbl_module::fhbl_kms, cbl_module::gamq_kgkgm, cbl_module::gamt_km, mod_grav::grav, cbl_module::isubs, defaultnotused::notused, defaultnotused::notusedi, cbl_module::qp_kgkg, cbl_module::tm_k, cbl_module::tp_k, cbl_module::tpp_k, and cbl_module::wsb.

Referenced by rkutta().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ gamma_sonde()

subroutine bluews_module::gamma_sonde

Definition at line 887 of file suews_phys_bluews.f95.

888 USE cbl_module
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 END IF
902 gamt_km = gamtt
903 END DO
904
905 DO j = 2, zmax
906 IF (blh_m >= ghum(j - 1, 1)) THEN
907 gamqq = ghum(j - 1, 2)
908 END IF
909 gamq_kgkgm = gamqq/1000.
910 END DO
911 END IF
912 RETURN
913
cbl_module::gtheta
real(kind(1d0)), dimension(0:500, 2) gtheta
Definition suews_phys_bluews.f95:51
cbl_module::ghum
real(kind(1d0)), dimension(0:500, 2) ghum
Definition suews_phys_bluews.f95:51
cbl_module::zmax
integer zmax
Definition suews_phys_bluews.f95:20

References cbl_module::blh_m, cbl_module::gamq_kgkgm, cbl_module::gamt_km, cbl_module::ghum, cbl_module::gtheta, cbl_module::sondeflag, and cbl_module::zmax.

Referenced by cbl().

Here is the caller graph for this function:

◆ nbl()

subroutine bluews_module::nbl ( integer, intent(in) iy,
integer, intent(in) id,
integer, intent(in) it,
integer, intent(in) imin,
real(kind(1d0)), intent(inout) dectime,
integer, intent(in) ir,
integer, intent(in) qh_choice,
real(kind(1d0)) qh_use,
real(kind(1d0)) qe_use,
real(kind(1d0)), intent(out) tm_k_zm,
real(kind(1d0)), intent(out) qm_gkg_zm,
integer startflag,
integer gridiv,
real(kind(1d0)), intent(in) psih,
real(kind(1d0)), intent(in) ustar,
real(kind(1d0)), intent(out) temp_c,
integer, intent(in) numberofgrids,
real(kind(1d0)), dimension(numberofgrids), intent(in) qhforcbl,
real(kind(1d0)), dimension(numberofgrids), intent(in) qeforcbl,
real(kind(1d0)), intent(in) press_hpa,
real(kind(1d0)), intent(inout) avu1,
real(kind(1d0)), intent(inout) avrh,
integer, intent(in) readlinesmetdata,
real(kind(1d0)), dimension(readlinesmetdata, ncolumnsdataoutbl, numberofgrids), intent(out) dataoutbl,
real(kind(1d0)), intent(inout) avcp,
real(kind(1d0)), intent(inout) avdens,
real(kind(1d0)), intent(inout) es_hpa,
real(kind(1d0)), intent(inout) lv_j_kg )

Definition at line 482 of file suews_phys_bluews.f95.

487
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 END IF
510
511 nlineday = 0
512 DO i = 1, nlineindata
513 IF (int(inicbldata(i, 1)) <= id) THEN
514 nlineday = nlineday + 1
515 END IF
516 END DO
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 END IF
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 END IF
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 END IF
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 END IF
615 IF (qm_kgkg > qsatf(tm_c, press_hpa) .OR. qm_kgkg < 0) THEN
616 qm_kgkg = qsatf(tm_c, press_hpa)
617 END IF
618
619 startflag = 0
620
621 ! iCBLcount = iCBLcount + 1 ! move counter at the end, TS 27 Aug 2019

References errorhint(), and sonde().

Referenced by cbl().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ rkutta()

subroutine bluews_module::rkutta ( integer neqn_use,
real(kind(1d0)) xa,
real(kind(1d0)) xb,
real(kind(1d0)), dimension(neqn_use) y_use,
integer nsteps )

Definition at line 631 of file suews_phys_bluews.f95.

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 END DO
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 END IF
684
685 DO nj = 1, neqn_use
686 arg(nj) = y_use(nj) + coef(n)*rk(nj, n)
687 END DO
688
689 CALL diff(xx, arg, dydx)
690
691 DO nj = 1, neqn_use
692 rk(nj, n + 1) = step*dydx(nj)
693 END DO
694 END DO
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 END DO
700 END DO
701 END DO
702
703 RETURN

References diff().

Referenced by cbl().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ sonde()

subroutine bluews_module::sonde ( integer id)

Definition at line 853 of file suews_phys_bluews.f95.

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 END DO
875900 zmax = i - 1
876 IF (zmax > izm) THEN
877 CALL errorhint(23, filen, real(zmax, kind(1d0)), notused, izm)
878 END IF
879 CLOSE (fn)
880 RETURN
88124 CALL errorhint(24, filen, notused, notused, notusedi)
88225 CALL errorhint(25, filen, notused, notused, i)
883 RETURN
cbl_module::filesonde
character(len=200), dimension(366) filesonde
Definition suews_phys_bluews.f95:13
data_in::fileinputpath
character(len=150) fileinputpath
Definition suews_ctrl_const.f95:1103

References errorhint(), data_in::fileinputpath, and cbl_module::filesonde.

Referenced by cbl_initial(), and nbl().

Here is the call graph for this function:
Here is the caller graph for this function:
Generated on Mon Mar 25 2024 14:48:39 for SUEWS API Site by doxygen 1.10.0