module_chem_plumerise_scalar.F
References to this file elsewhere.
1 Module module_chem_plumerise_scalar
2 ! USE module_plumerise1
3 use module_model_constants
4 use module_zero_plumegen_coms
5 real,parameter :: rgas=r_d
6 real,parameter :: cpor=1./rcp
7 real,parameter :: p00=p1000mb
8 ! real, external:: esat_pr
9 CONTAINS
10 subroutine plumerise(m1,m2,m3,ia,iz,ja,jz, &
11 nspecies,eburn_in,eburn_out &
12 ,up,vp,wp,theta,pp,dn0,rv,zt_rams,zm_rams)
13 ! use module_zero_plumegen_coms, only : ucon,vcon,wcon,thtcon,rvcon,picon,tmpcon&
14 ! ,dncon,prcon,zcon,zzcon,scon
15
16 implicit none
17
18
19 integer :: ng,m1,m2,m3,ia,iz,ja,jz,ibcon,mynum,i,j,k,iveg_ag,&
20 imm,k1,k2,ixx,ispc,nspecies
21
22 integer :: ncall = 0
23 integer :: kmt
24 real,dimension(m1,nspecies), intent(out) :: eburn_out
25 real,dimension(nspecies), intent(in) :: eburn_in
26
27 real, dimension(m1,m2,m3) :: up, vp, wp,theta,pp,dn0,rv
28
29 real, dimension(m1) :: zt_rams,zm_rams
30
31 real :: burnt_area,STD_burnt_area,dz_flam,rhodzi,dzi
32 real, dimension(2) :: ztopmax
33
34 real :: q_smold_kgm2
35
36 ! From plumerise1.F routine
37 integer, parameter :: nveg_agreg = 4
38 integer, parameter :: tropical_forest = 1
39 integer, parameter :: boreal_forest = 2
40 integer, parameter :: savannah = 3
41
42 integer, parameter :: grassland = 4
43 real, dimension(nveg_agreg) :: firesize,mean_fct
44
45
46 !Fator de conversao de unidades
47 !!fcu=1. !=> kg [gas/part] /kg [ar]
48 !!fcu =1.e+12 !=> ng [gas/part] /kg [ar]
49 !!real,parameter :: fcu =1.e+6 !=> mg [gas/part] /kg [ar]
50 !----------------------------------------------------------------------
51 ! indexacao para o array "plume(k,i,j)"
52 ! k
53 ! 1 => area media (m^2) dos focos em biomas floresta dentro do gribox i,j
54 ! 2 => area media (m^2) dos focos em biomas savana dentro do gribox i,j
55 ! 3 => area media (m^2) dos focos em biomas pastagem dentro do gribox i,j
56 ! 4 => desvio padrao da area media (m^2) dos focos : floresta
57 ! 5 => desvio padrao da area media (m^2) dos focos : savana
58 ! 6 => desvio padrao da area media (m^2) dos focos : pastagem
59 ! 7 a 9 => sem uso
60 !10(=k_CO_smold) => parte da emissao total de CO correspondente a fase smoldering
61 !11, 12 e 13 => este array guarda a relacao entre
62 ! qCO( flaming, floresta) e a quantidade total emitida
63 ! na fase smoldering, isto e;
64 ! qCO( flaming, floresta) = plume(11,i,j)*plume(10,i,j)
65 ! qCO( flaming, savana ) = plume(12,i,j)*plume(10,i,j)
66 ! qCO( flaming, pastagem) = plume(13,i,j)*plume(10,i,j)
67 !20(=k_PM25_smold),21,22 e 23 o mesmo para PM25
68 !
69 !24-n1 => sem uso
70 !----------------------------------------------------------------------
71 if (ncall == 0) then
72 ncall = 1
73 call zero_plumegen_coms
74 endif
75
76
77 do j = ja,jz
78 do i = ia,iz
79
80 !- if the max value of flaming is close to zero => there is not emission with
81 !- plume rise => cycle
82
83 do k = 1,m1
84 !ucon (k)=up(k,i,j) ! u wind
85 !vcon (k)=vp(k,i,j) ! v wind
86 !wcon (k)=wp(k,i,j) ! w wind
87 thtcon(k)=theta(k,i,j) ! pot temperature
88 picon (k)=pp(k,i,j) ! exner function
89 !tmpcon(k)=thtcon(k)*picon(k)/cp ! temperature (K)
90 !dncon (k)=dn0(k,i,j) ! dry air density (basic state)
91 !prcon (k)=(picon(k)/cp)**cpor*p00 ! pressure (Pa)
92 rvcon (k)=rv(k,i,j) ! water vapor mixing ratio
93 zcon (k)=zt_rams(k) ! termod-point height
94 zzcon (k)=zm_rams(k) ! W-point height
95 !if(i==43 .and. j==42)print*,'PL:',k,zcon(k),prcon (k),tmpcon(k)-273.,1000.*rvcon (k)
96 enddo
97 do ispc=1,nspecies
98 eburn_out(1,ispc) = eburn_in(ispc)
99 enddo
100
101 !- get envinronmental state (temp, water vapor mix ratio, ...)
102 call get_env_condition(1,m1-1,kmt)
103
104 !- loop nos 4 biomas agregados com possivel queimada
105 do iveg_ag=1,nveg_agreg
106
107
108 !- verifica a existencia de emissao flaming para um bioma especifico
109 !orig: if( plume( k_CO_smold + iveg_ag ,i,j) < 1.e-6 ) cycle
110 if(mean_fct(iveg_ag) < 1.e-6 ) cycle
111
112 ! burnt area and standard desviation
113 burnt_area = firesize(iveg_ag)
114
115 !not em use
116 !STD_burnt_area= plume(3+iveg_ag,i,j)
117 STD_burnt_area= 0.
118
119 !- loop nos valores minimo e maximo da taxa de calor
120 do imm=1,2
121
122 !--------------------
123 !ixx=iveg_ag*10 + imm
124 !print*,'i j veg=',i, j, iveg_ag,imm
125 !--------------------
126
127 !- get fire properties (burned area, plume radius, heating rates ...)
128 call get_fire_properties(imm,iveg_ag,burnt_area,STD_burnt_area)
129
130 !------ generates the plume rise ------
131
132 !-- only one value for eflux of GRASSLAND
133 ! if(iveg_ag == GRASSLAND .and. imm == 2) then
134 if(iveg_ag == 4 .and. imm == 2) then
135 ztopmax(2)=ztopmax(1)
136 ztopmax(1)=zzcon(1)
137 cycle
138 endif
139
140 call makeplume (kmt,ztopmax(imm),ixx)
141
142 enddo ! enddo do loop em imm
143
144 !- define o dominio vertical onde a emissao flaming ira ser colocada
145 call set_flam_vert(ztopmax,k1,k2)
146
147 !- espessura da camada vertical
148 dz_flam=zzcon(k2)-zzcon(k1-1)
149
150 !- distribui a emissao flaming entre os niveis k1 e k2
151 do k=k1,k2
152 !use this in case the emission src is already in mixing ratio
153 !rhodzi= 1./(dn0(k,i,j) * dz_flam)
154 !use this in case the emission src is tracer density
155 dzi= 1./( dz_flam)
156
157 do ispc = 1, nspecies
158
159 !- get back the smoldering emission in kg/m2 (actually in 1e-9 kg/m2)
160
161 !use this in case the emission src is already in mixing ratio
162 !q_smold_kgm2 = (1/dzt(2) * dn0(2,i,j) )* &
163 ! chem1_src_g(bburn,ispc,ng)%sc_src(2,i,j)
164
165 !use this in case the emission src is tracer density
166 ! q_smold_kgm2 = ((zt_rams(2)-zt_rams(1)) )* &
167 ! eburn_in(ispc)
168 q_smold_kgm2 = eburn_in(ispc)
169
170 ! units = already in ppbm, don't need "fcu" factor
171 eburn_out(k,ispc) = eburn_out(k,ispc) +&
172 mean_fct(iveg_ag) *&
173 q_smold_kgm2 * &
174 dzi !use this in case the emission src is tracer density
175 !rhodzi !use this in case the emission src is already in mixing ratio
176
177
178
179 !srcCO(k,i,j)= srcCO(k,i,j) + plume(k_CO_smold+iveg_ag,i,j)*&
180 ! plume(k_CO_smold ,i,j)*&
181 ! rhodzi*fcu
182 enddo
183
184 enddo
185
186 enddo ! enddo do loop em iveg_ag
187 !-----
188 !stop 333
189 !endif
190 !-----
191 enddo ! loop em i
192 enddo ! loop em j
193 !stop 4544
194 end subroutine plumerise
195 !-------------------------------------------------------------------------
196
197 subroutine get_env_condition(k1,k2,kmt)
198
199 !se module_zero_plumegen_coms
200 !use rconstants
201 implicit none
202 integer :: k1,k2,k,kcon,klcl,kmt,nk,nkmid,i
203 real :: znz,themax,tlll,plll,rlll,zlll,dzdd,dzlll,tlcl,plcl,dzlcl,dummy
204 integer :: n_setgrid = 0
205
206 if( n_setgrid == 0) then
207 n_setgrid = 1
208 call set_grid ! define vertical grid of plume model
209 ! zt(k) = thermo and water levels
210 ! zm(k) = dynamical levels
211 endif
212
213 znz=zcon(k2)
214 do k=nkp,1,-1
215 if(zt(k).lt.znz)go to 13
216 enddo
217 stop ' envir stop 12'
218 13 continue
219 kmt=k
220
221 nk=k2-k1+1
222 !call htint(nk, wcon,zzcon(k1),kmt,wpe,zt)
223 !call htint(nk, ucon,zcon(k1),kmt,upe,zt)
224 !call htint(nk, vcon,zcon(k1),kmt,vpe,zt)
225 call htint(nk,thtcon,zcon(k1),kmt,the ,zt)
226 call htint(nk, rvcon,zcon(k1),kmt,qvenv,zt)
227 do k=1,kmt
228 qvenv(k)=max(qvenv(k),1e-8)
229 enddo
230
231 pke(1)=picon(1)
232 do k=1,kmt
233 thve(k)=the(k)*(1.+.61*qvenv(k)) ! virtual pot temperature
234 enddo
235 do k=2,kmt
236 pke(k)=pke(k-1)-g*2.*(zt(k)-zt(k-1)) & ! exner function
237 /(thve(k)+thve(k-1))
238 enddo
239 do k=1,kmt
240 te(k) = the(k)*pke(k)/cp ! temperature (K)
241 pe(k) = (pke(k)/cp)**cpor*p00 ! pressure (Pa)
242 dne(k)= pe(k)/(rgas*te(k)*(1.+.61*qvenv(k))) ! dry air density (kg/m3)
243 ! print*,'ENV=',qvenv(k)*1000., te(k)-273.15,zt(k)
244 enddo
245
246 !-use este para gerar o RAMS.out
247 ! ------- print environment state
248 !print*,'k,zt(k),pe(k),te(k)-273.15,qvenv(k)*1000'
249 !do k=1,kmt
250 ! write(*,100) k,zt(k),pe(k),te(k)-273.15,qvenv(k)*1000.
251 ! 100 format(1x,I5,4f20.12)
252 !enddo
253 !stop 333
254
255
256 !--------- nao eh necessario este calculo
257 !do k=1,kmt
258 ! call thetae(pe(k),te(k),qvenv(k),thee(k))
259 !enddo
260
261
262 !--------- converte press de Pa para kPa para uso modelo de plumerise
263 do k=1,kmt
264 pe(k) = pe(k)*1.e-3
265 enddo
266
267 return
268 end subroutine get_env_condition
269
270 !-------------------------------------------------------------------------
271
272 subroutine set_grid()
273 !use module_zero_plumegen_coms
274 implicit none
275 integer :: k,mzp
276
277 dz=100. ! set constant grid spacing of plume grid model(meters)
278
279 mzp=nkp
280 zt(1) = zsurf
281 zm(1) = zsurf
282 zt(2) = zt(1) + 0.5*dz
283 zm(2) = zm(1) + dz
284 do k=3,mzp
285 zt(k) = zt(k-1) + dz ! thermo and water levels
286 zm(k) = zm(k-1) + dz ! dynamical levels
287 enddo
288 !print*,zsurf
289 !Print*,zt(:)
290 do k = 1,mzp-1
291 dzm(k) = 1. / (zt(k+1) - zt(k))
292 enddo
293 dzm(mzp)=dzm(mzp-1)
294
295 do k = 2,mzp
296 dzt(k) = 1. / (zm(k) - zm(k-1))
297 enddo
298 dzt(1) = dzt(2) * dzt(2) / dzt(3)
299
300 ! dzm(1) = 0.5/dz
301 ! dzm(2:mzp) = 1./dz
302 return
303 end subroutine set_grid
304 !-------------------------------------------------------------------------
305
306 subroutine set_flam_vert(ztopmax,k1,k2)
307 ! use module_zero_plumegen_coms, only : nkp,zzcon
308 implicit none
309 integer imm,k,k1,k2
310 real, dimension(2) :: ztopmax
311 integer, dimension(2) :: k_lim
312
313
314 do imm=1,2
315 ! checar
316 ! do k=1,m1-1
317 do k=1,nkp-1
318 if(zzcon(k) > ztopmax(imm) ) exit
319 enddo
320 k_lim(imm) = k
321 enddo
322 k1=max(3,k_lim(1))
323 k2=max(3,k_lim(2))
324
325 if(k2 < k1) then
326 !print*,'1: ztopmax k=',ztopmax(1), k1
327 !print*,'2: ztopmax k=',ztopmax(2), k2
328 k2=k1
329 !stop 1234
330 endif
331
332 end subroutine set_flam_vert
333 !-------------------------------------------------------------------------
334
335 subroutine get_fire_properties(imm,iveg_ag,burnt_area,STD_burnt_area)
336 !use module_zero_plumegen_coms
337 implicit none
338 integer :: moist, i, icount,imm,iveg_ag
339 real:: bfract, effload, heat, hinc ,burnt_area,STD_burnt_area,heat_fluxW
340 real, dimension(2,4) :: heat_flux
341
342
343 data heat_flux/ &
344 !---------------------------------------------------------------------
345 ! heat flux !IGBP Land Cover !
346 ! min ! max !Legend and ! reference
347 ! kW/m^2 !description !
348 !--------------------------------------------------------------------
349 30.0, 80.0, &! Tropical Forest ! igbp 2 & 4
350 30.0, 80.0, &! Boreal forest ! igbp 1 & 3
351 4.4, 23.0, &! cerrado/woody savanna | igbp 5 thru 9
352 3.3, 3.3 /! Grassland/cropland ! igbp 10 thru 17
353 !--------------------------------------------------------------------
354
355
356 !-- fire at the surface
357 !
358 !area = 20.e+4 ! area of burn, m^2
359 area = burnt_area! area of burn, m^2
360
361 !fluxo de calor para o bioma
362 heat_fluxW = heat_flux(imm,iveg_ag) * 1000. ! converte para W/m^2
363
364 mdur = 53 ! duration of burn, minutes
365 bload = 10. ! total loading, kg/m**2
366 moist = 10 ! fuel moisture, %. average fuel moisture,percent dry
367 maxtime =mdur+2 ! model time, min
368 !heat = 21.e6 !- joules per kg of fuel consumed
369 !heat = 15.5e6 !joules/kg - cerrado
370 heat = 19.3e6 !joules/kg - floresta em alta floresta (mt)
371 alpha = 0.1 !- entrainment constant
372
373 !-------------------- printout ----------------------------------------
374
375 !!WRITE ( * , * ) ' SURFACE =', ZSURF, 'M', ' LCL =', ZBASE, 'M'
376 !
377 !PRINT*,'======================================================='
378 !print * , ' FIRE BOUNDARY CONDITION :'
379 !print * , ' DURATION OF BURN, MINUTES =',MDUR
380 !print * , ' AREA OF BURN, HA =',AREA*1.e-4
381 !print * , ' HEAT FLUX, kW/m^2 =',heat_fluxW*1.e-3
382 !print * , ' TOTAL LOADING, KG/M**2 =',BLOAD
383 !print * , ' FUEL MOISTURE, % =',MOIST !average fuel moisture,percent dry
384 !print * , ' MODEL TIME, MIN. =',MAXTIME
385 !
386 !
387 !
388 ! ******************** fix up inputs *********************************
389 !
390
391 !IF (MOD (MAXTIME, 2) .NE.0) MAXTIME = MAXTIME+1 !make maxtime even
392
393 MAXTIME = MAXTIME * 60 ! and put in seconds
394 !
395 RSURF = SQRT (AREA / 3.14159) !- entrainment surface radius (m)
396
397 FMOIST = MOIST / 100. !- fuel moisture fraction
398 !
399 !
400 ! calculate the energy flux and water content at lboundary.
401 ! fills heating() on a minute basis. could ask for a file at this po
402 ! in the program. whatever is input has to be adjusted to a one
403 ! minute timescale.
404 !
405
406 DO I = 1, ntime !- make sure of energy release
407 HEATING (I) = 0.0001 !- avoid possible divide by 0
408 enddo
409 !
410 TDUR = MDUR * 60. !- number of seconds in the burn
411
412 bfract = 1. !- combustion factor
413
414 EFFLOAD = BLOAD * BFRACT !- patchy burning
415
416 ! spread the burning evenly over the interval
417 ! except for the first few minutes for stability
418 ICOUNT = 1
419 !
420 if(MDUR > NTIME) STOP 'Increase time duration (ntime) in min - see file "plumerise_mod.f90"'
421
422 DO WHILE (ICOUNT.LE.MDUR)
423 ! HEATING (ICOUNT) = HEAT * EFFLOAD / TDUR ! W/m**2
424 ! HEATING (ICOUNT) = 80000. * 0.55 ! W/m**2
425
426 HEATING (ICOUNT) = heat_fluxW * 0.55 ! W/m**2 (0.55 converte para energia convectiva)
427 ICOUNT = ICOUNT + 1
428 ENDDO
429 ! ramp for 5 minutes
430
431 HINC = HEATING (1) / 4.
432 HEATING (1) = 0.1
433 HEATING (2) = HINC
434 HEATING (3) = 2. * HINC
435 HEATING (4) = 3. * HINC
436 !
437
438 return
439 end subroutine get_fire_properties
440 !-------------------------------------------------------------------------------
441 !
442 SUBROUTINE MAKEPLUME ( kmt,ztopmax,ixx)
443 !
444 ! *********************************************************************
445 !
446 ! EQUATION SOURCE--Kessler Met.Monograph No. 32 V.10 (K)
447 ! Alan Weinstein, JAS V.27 pp 246-255. (W),
448 ! Ogura and Takahashi, Monthly Weather Review V.99,pp895-911 (OT)
449 ! Roger Pielke,Mesoscale Meteorological Modeling,Academic Press,1984
450 ! Originally developed by: Don Latham (USFS)
451 !
452 !
453 ! ************************ VARIABLE ID ********************************
454 !
455 ! DT=COMPUTING TIME INCREMENT (SEC)
456 ! DZ=VERTICAL INCREMENT (M)
457 ! LBASE=LEVEL ,CLOUD BASE
458 !
459 ! CONSTANTS:
460 ! G = GRAVITATIONAL ACCELERATION 9.80796 (M/SEC/SEC).
461 ! R = DRY AIR GAS CONSTANT (287.04E6 JOULE/KG/DEG K)
462 ! CP = SPECIFIC HT. (1004 JOULE/KG/DEG K)
463 ! HEATCOND = HEAT OF CONDENSATION (2.5E6 JOULE/KG)
464 ! HEATFUS = HEAT OF FUSION (3.336E5 JOULE/KG)
465 ! HEATSUBL = HEAT OF SUBLIMATION (2.83396E6 JOULE/KG)
466 ! EPS = RATIO OF MOL.WT. OF WATER VAPOR TO THAT OF DRY AIR (0.622)
467 ! DES = DIFFERENCE BETWEEN VAPOR PRESSURE OVER WATER AND ICE (MB)
468 ! TFREEZE = FREEZING TEMPERATURE (K)
469 !
470 !
471 ! PARCEL VALUES:
472 ! T = TEMPERATURE (K)
473 ! TXS = TEMPERATURE EXCESS (K)
474 ! QH = HYDROMETEOR WATER CONTENT (G/G DRY AIR)
475 ! QHI = HYDROMETEOR ICE CONTENT (G/G DRY AIR)
476 ! QC = WATER CONTENT (G/G DRY AIR)
477 ! QVAP = WATER VAPOR MIXING RATIO (G/G DRY AIR)
478 ! QSAT = SATURATION MIXING RATIO (G/G DRY AIR)
479 ! RHO = DRY AIR DENSITY (G/M**3) MASSES = RHO*Q'S IN G/M**3
480 ! ES = SATURATION VAPOR PRESSURE (kPa)
481 !
482 ! ENVIRONMENT VALUES:
483 ! TE = TEMPERATURE (K)
484 ! PE = PRESSURE (kPa)
485 ! QVENV = WATER VAPOR (G/G)
486 ! RHE = RELATIVE HUMIDITY FRACTION (e/esat)
487 ! DNE = dry air density (kg/m^3)
488 !
489 ! HEAT VALUES:
490 ! HEATING = HEAT OUTPUT OF FIRE (WATTS/M**2)
491 ! MDUR = DURATION OF BURN, MINUTES
492 !
493 ! W = VERTICAL VELOCITY (M/S)
494 ! RADIUS=ENTRAINMENT RADIUS (FCN OF Z)
495 ! RSURF = ENTRAINMENT RADIUS AT GROUND (SIMPLE PLUME, TURNER)
496 ! ALPHA = ENTRAINMENT CONSTANT
497 ! MAXTIME = TERMINATION TIME (MIN)
498 !
499 !
500 !**********************************************************************
501 !**********************************************************************
502 !use module_zero_plumegen_coms
503 implicit none
504 !logical :: endspace
505 character (len=10) :: varn
506 integer :: izprint, iconv, itime, k, kk, kkmax, deltak,ilastprint,kmt &
507 ,ixx,nrectotal,i_micro,n_sub_step
508 real :: vc, g, r, cp, eps, &
509 tmelt, heatsubl, heatfus, heatcond, tfreeze, &
510 ztopmax, wmax, rmaxtime, es, esat, heat,dt_save !ESAT_PR,
511 character (len=2) :: cixx
512
513 ! real, external:: esat_pr!
514 !
515 ! ******************* SOME CONSTANTS **********************************
516 !
517 ! XNO=10.0E06 median volume diameter raindrop (K table 4)
518 ! VC = 38.3/(XNO**.125) mean volume fallspeed eqn. (K)
519 !
520 parameter (vc = 5.107387)
521 parameter (g = 9.80796, r = 287.04, cp = 1004., eps = 0.622, tmelt = 273.3)
522 parameter (heatsubl = 2.834e6, heatfus = 3.34e5, heatcond = 2.501e6)
523 parameter (tfreeze = 269.3)
524 !
525 tstpf = 2.0 !- timestep factor
526 viscosity = 500.!- viscosity constant (original value: 0.001)
527
528 nrectotal=150
529 !
530 !*************** PROBLEM SETUP AND INITIAL CONDITIONS *****************
531 mintime = 1
532 ztopmax = 0.
533 ztop = 0.
534 time = 0.
535 dt = 1.
536 wmax = 1.
537 kkmax = 10
538 deltaK = 20
539 ilastprint=0
540 L = 1 ! L initialization
541
542 !--- initialization
543 CALL INITIAL(kmt)
544
545 !--- initial print fields:
546 izprint = 0 ! if = 0 => no printout
547 if (izprint.ne.0) then
548 write(cixx(1:2),'(i2.2)') ixx
549 open(2, file = 'debug.'//cixx//'.dat')
550 open(19,file='plumegen9.'//cixx//'.gra', &
551 form='unformatted',access='direct',status='unknown', &
552 recl=4*nrectotal) !PC
553 ! recl=1*nrectotal) !sx6 e tupay
554 call printout (izprint,nrectotal)
555 ilastprint=2
556 endif
557
558 ! ******************* model evolution ******************************
559 rmaxtime = float(maxtime)
560 !
561 DO WHILE (TIME.LE.RMAXTIME) !beginning of time loop
562
563 ! do itime=1,120
564
565 !-- set model top integration
566 nm1 = min(kmt, kkmax + deltak)
567
568 !-- set timestep
569 !dt = (zm(2)-zm(1)) / (tstpf * wmax)
570 dt = min(5.,(zm(2)-zm(1)) / (tstpf * wmax))
571
572 !-- elapsed time, sec
573 time = time+dt
574 !-- elapsed time, minutes
575 mintime = 1 + int (time) / 60
576 wmax = 1. !no zeroes allowed.
577 !************************** BEGIN SPACE LOOP **************************
578
579 !-- zerout all model tendencies
580 call tend0_plumerise
581
582 !-- bounday conditions (k=1)
583 L=1
584 call lbound()
585
586 !-- dynamics for the level k>1
587 !-- W advection
588 ! call vel_advectc_plumerise(NM1,WC,WT,DNE,DZM)
589 call vel_advectc_plumerise(NM1,WC,WT,RHO,DZM)
590
591 !-- scalars advection 1
592 call scl_advectc_plumerise('SC',NM1)
593
594 !-- scalars advection 2
595 !call scl_advectc_plumerise2('SC',NM1)
596
597 !-- scalars entrainment, adiabatic
598 call scl_misc(NM1)
599
600 !-- gravity wave damping using Rayleigh friction layer fot T
601 call damp_grav_wave(1,nm1,deltak,dt,zt,zm,w,t,tt,qv,qh,qi,qc,te,pe,qvenv)
602
603 !-- microphysics
604 ! goto 101 ! bypass microphysics
605 dt_save=dt
606 n_sub_step=3
607 dt=dt/float(n_sub_step)
608
609 do i_micro=1,n_sub_step
610 !-- sedim ?
611 call fallpart(NM1)
612 !-- microphysics
613 do L=2,nm1-1
614 WBAR = 0.5*(W(L)+W(L-1))
615 ES = ESAT_PR (T(L)) !BLOB SATURATION VAPOR PRESSURE, EM KPA
616 QSAT(L) = (EPS * ES) / (PE(L) - ES) !BLOB SATURATION LWC G/G DRY AIR
617 EST (L) = ES
618 RHO (L) = 3483.8 * PE (L) / T (L) ! AIR PARCEL DENSITY , G/M**3
619 !srf18jun2005
620 ! IF (W(L) .ge. 0.) DQSDZ = (QSAT(L ) - QSAT(L-1)) / (ZT(L ) -ZT(L-1))
621 ! IF (W(L) .lt. 0.) DQSDZ = (QSAT(L+1) - QSAT(L )) / (ZT(L+1) -ZT(L ))
622 IF (W(L) .ge. 0.) then
623 DQSDZ = (QSAT(L+1) - QSAT(L-1)) / (ZT(L+1 )-ZT(L-1))
624 ELSE
625 DQSDZ = (QSAT(L+1) - QSAT(L-1)) / (ZT(L+1) -ZT(L-1))
626 ENDIF
627
628 call waterbal
629 enddo
630 enddo
631 dt=dt_save
632 !
633 101 continue
634 !
635 !-- W-viscosity for stability
636 call visc_W(nm1,deltak,kmt)
637
638 !-- update scalars
639 call update_plumerise(nm1,'S')
640
641 call hadvance_plumerise(1,nm1,dt,WC,WT,W,mintime)
642
643 !-- Buoyancy
644 call buoyancy_plumerise(NM1, T, TE, QV, QVENV, QH, QI, QC, WT, SCR1)
645
646 !-- Entrainment
647 call entrainment(NM1,W,WT,RADIUS,ALPHA)
648
649 !-- update W
650 call update_plumerise(nm1,'W')
651
652 call hadvance_plumerise(2,nm1,dt,WC,WT,W,mintime)
653
654
655 !-- misc
656 do k=2,nm1
657 ! pe esta em kpa - esat do rams esta em mbar = 100 Pa = 0.1 kpa
658 ! es = 0.1*esat (t(k)) !blob saturation vapor pressure, em kPa
659 ! rotina do plumegen calcula em kPa
660 es = esat_pr (t(k)) !blob saturation vapor pressure, em kPa
661 qsat(k) = (eps * es) / (pe(k) - es) !blob saturation lwc g/g dry air
662 est (k) = es
663 txs (k) = t(k) - te(k)
664 rho (k) = 3483.8 * pe (k) / t (k) ! air parcel density , g/m**3
665 ! no pressure diff with radius
666
667 if((abs(wc(k))).gt.wmax) wmax = abs(wc(k)) ! keep wmax largest w
668 enddo
669
670 ! Gravity wave damping using Rayleigh friction layer for W
671 call damp_grav_wave(2,nm1,deltak,dt,zt,zm,w,t,tt,qv,qh,qi,qc,te,pe,qvenv)
672 !---
673
674 !-- try to find the plume top (above surface height)
675 kk = 1
676 do while (w (kk) .gt. 1.)
677 kk = kk + 1
678 ztop = zm(kk)
679 !print*,'W=',w (kk)
680 enddo
681 !
682 ztop_(mintime) = ztop
683 ztopmax = max (ztop, ztopmax)
684 kkmax = max (kk , kkmax )
685 !print * ,'ztopmax=', mintime,'mn ',ztop_(mintime), ztopmax
686
687 !
688 !srf-27082005
689 ! if the solution is going to a stationary phase, exit
690 if(mintime > 10) then
691 if( abs(ztop_(mintime)-ztop_(mintime-10)) < DZ ) exit
692 endif
693
694 if(ilastprint == mintime) then
695 call printout (izprint,nrectotal)
696 ilastprint = mintime+1
697 endif
698
699
700 ENDDO !do next timestep
701
702 !print * ,' ztopmax=',ztopmax,'m',mintime,'mn '
703 !print*,'======================================================='
704 !
705 !the last printout
706 if (izprint.ne.0) then
707 call printout (izprint,nrectotal)
708 close (2)
709 close (19)
710 endif
711
712 RETURN
713 END SUBROUTINE MAKEPLUME
714 !-------------------------------------------------------------------------------
715 !
716 SUBROUTINE BURN(EFLUX, WATER)
717 !
718 !- calculates the energy flux and water content at lboundary
719 !use module_zero_plumegen_coms
720 !real, parameter :: HEAT = 21.E6 !Joules/kg
721 !real, parameter :: HEAT = 15.5E6 !Joules/kg - cerrado
722 real, parameter :: HEAT = 19.3E6 !Joules/kg - floresta em Alta Floresta (MT)
723 real :: eflux,water
724 !
725 ! The emission factor for water is 0.5. The water produced, in kg,
726 ! is then fuel mass*0.5 + (moist/100)*mass per square meter.
727 ! The fire burns for DT out of TDUR seconds, the total amount of
728 ! fuel burned is AREA*BLOAD*(DT/TDUR) kg. this amount of fuel is
729 ! considered to be spread over area AREA and so the mass burned per
730 ! unit area is BLOAD*(DT/TDUR), and the rate is BLOAD/TDUR.
731 !
732 IF (TIME.GT.TDUR) THEN !is the burn over?
733 EFLUX = 0.000001 !prevent a potential divide by zero
734 WATER = 0.
735 RETURN
736 ELSE
737 !
738 EFLUX = HEATING (MINTIME) ! Watts/m**2
739 ! WATER = EFLUX * (DT / HEAT) * (0.5 + FMOIST) ! kg/m**2
740 WATER = EFLUX * (DT / HEAT) * (0.5 + FMOIST) /0.55 ! kg/m**2
741 WATER = WATER * 1000. ! g/m**2
742 !
743 ! print*,'BURN:',time,EFLUX/1.e+9
744 ENDIF
745 !
746 RETURN
747 END SUBROUTINE BURN
748 !-------------------------------------------------------------------------------
749 !
750 SUBROUTINE LBOUND ()
751 !
752 ! ********** BOUNDARY CONDITIONS AT ZSURF FOR PLUME AND CLOUD ********
753 !
754 ! source of equations: J.S. Turner Buoyancy Effects in Fluids
755 ! Cambridge U.P. 1973 p.172,
756 ! G.A. Briggs Plume Rise, USAtomic Energy Commissio
757 ! TID-25075, 1969, P.28
758 !
759 ! fundamentally a point source below ground. at surface, this produces
760 ! a velocity w(1) and temperature T(1) which vary with time. There is
761 ! also a water load which will first saturate, then remainder go into
762 ! QC(1).
763 ! EFLUX = energy flux at ground,watt/m**2 for the last DT
764 !
765 !use module_zero_plumegen_coms
766 implicit none
767 real, parameter :: g = 9.80796, r = 287.04, cp = 1004.6, eps = 0.622,tmelt = 273.3
768 real, parameter :: tfreeze = 269.3, pi = 3.14159, e1 = 1./3., e2 = 5./3.
769 real :: es, esat, eflux, water, pres, c1, c2, f, zv, denscor, xwater !,ESAT_PR
770 ! real, external:: esat_pr!
771
772 !
773 QH (1) = QH (2) !soak up hydrometeors
774 QI (1) = QI (2)
775 QC (1) = 0. !no cloud here
776 !
777 !
778 CALL BURN (EFLUX, WATER)
779 !
780 ! calculate parameters at boundary from a virtual buoyancy point source
781 !
782 PRES = PE (1) * 1000. !need pressure in N/m**2
783
784 C1 = 5. / (6. * ALPHA) !alpha is entrainment constant
785
786 C2 = 0.9 * ALPHA
787
788 F = EFLUX / (PRES * CP * PI)
789
790 F = G * R * F * AREA !buoyancy flux
791
792 ZV = C1 * RSURF !virtual boundary height
793
794 W (1) = C1 * ( (C2 * F) **E1) / ZV**E1 !boundary velocity
795
796 DENSCOR = C1 * F / G / (C2 * F) **E1 / ZV**E2 !density correction
797
798 T (1) = TE (1) / (1. - DENSCOR) !temperature of virtual plume at zsurf
799
800 !
801 WC(1) = W(1)
802
803 !SC(1) = SCE(1)+F/1000.*dt ! gas/particle (g/g)
804
805 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
806 ! match dw/dz,dt/dz at the boundary. F is conserved.
807 !
808 !WBAR = W (1) * (1. - 1. / (6. * ZV) )
809 !ADVW = WBAR * W (1) / (3. * ZV)
810 !ADVT = WBAR * (5. / (3. * ZV) ) * (DENSCOR / (1. - DENSCOR) )
811 !ADVC = 0.
812 !ADVH = 0.
813 !ADVI = 0.
814 !ADIABAT = - WBAR * G / CP
815 VTH (1) = - 4.
816 VTI (1) = - 3.
817 TXS (1) = T (1) - TE (1)
818
819 VISC (1) = VISCOSITY
820
821 RHO (1) = 3483.8 * PE (1) / T (1) !air density at level 1, g/m**3
822
823 XWATER = WATER / (W (1) * DT * RHO (1) ) !firewater mixing ratio
824
825 QV (1) = XWATER + QVENV (1) !plus what's already there
826
827
828 ! PE esta em kPa - ESAT do RAMS esta em mbar = 100 Pa = 0.1 kPa
829 ! ES = 0.1*ESAT (T(1)) !blob saturation vapor pressure, em kPa
830 ! rotina do plumegen ja calcula em kPa
831 ES = ESAT_PR (T(1)) !blob saturation vapor pressure, em kPa
832
833 EST (1) = ES
834 QSAT (1) = (EPS * ES) / (PE (1) - ES) !blob saturation lwc g/g dry air
835
836 IF (QV (1) .gt. QSAT (1) ) THEN
837 QC (1) = QV (1) - QSAT (1) + QC (1) !remainder goes into cloud drops
838 QV (1) = QSAT (1)
839 ENDIF
840 !
841 CALL WATERBAL
842 !
843 RETURN
844 END SUBROUTINE LBOUND
845 !-------------------------------------------------------------------------------
846 !
847 SUBROUTINE INITIAL ( kmt)
848 !
849 ! ************* SETS UP INITIAL CONDITIONS FOR THE PROBLEM ************
850 !use module_zero_plumegen_coms
851 implicit none
852 real, parameter :: tfreeze = 269.3
853 integer :: isub, k, n1, n2, n3, lbuoy, itmp, isubm1 ,kmt
854 real :: xn1, xi, es, esat!,ESAT_PR
855 !
856 N=kmt
857 ! initialize temperature structure,to the end of equal spaced sounding,
858 do k = 1, N
859 TXS (k) = 0.0
860 W (k) = 0.0
861 T (k) = TE(k) !blob set to environment
862 WC(k) = 0.0
863 WT(k) = 0.0
864 QV(k) = QVENV (k) !blob set to environment
865 VTH(k) = 0. !initial rain velocity = 0
866 VTI(k) = 0. !initial ice velocity = 0
867 QH(k) = 0. !no rain
868 QI(k) = 0. !no ice
869 QC(k) = 0. !no cloud drops
870 ! PE esta em kPa - ESAT do RAMS esta em mbar = 100 Pa = 0.1 kPa
871 ! ES = 0.1*ESAT (T(k)) !blob saturation vapor pressure, em kPa
872 ! rotina do plumegen calcula em kPa
873 ES = ESAT_PR (T(k)) !blob saturation vapor pressure, em kPa
874 EST (k) = ES
875 QSAT (k) = (.622 * ES) / (PE (k) - ES) !saturation lwc g/g
876 RHO (k) = 3483.8 * PE (k) / T (k) !dry air density g/m**3
877 enddo
878
879 ! Initialize the entrainment radius, Turner-style plume
880 radius(1) = rsurf
881 do k=2,N
882 radius(k) = radius(k-1)+(6./5.)*alpha*(zt(k)-zt(k-1))
883 enddo
884
885 ! Initialize the viscosity
886 VISC (1) = VISCOSITY
887 do k=2,N
888 VISC (k) = VISCOSITY!max(1.e-3,visc(k-1) - 1.* VISCOSITY/float(nkp))
889 enddo
890 !-- Initialize gas/concentration
891 !DO k =10,20
892 ! SC(k) = 20.
893 !ENDDO
894 !stop 333
895
896 CALL LBOUND()
897
898 RETURN
899 END SUBROUTINE INITIAL
900 !-------------------------------------------------------------------------------
901 !
902 subroutine damp_grav_wave(ifrom,nm1,deltak,dt,zt,zm,w,t,tt,qv,qh,qi,qc,te,pe,qvenv)
903 implicit none
904 integer nm1,ifrom,deltak
905 real dt
906 real, dimension(nm1) :: w,t,tt,qv,qh,qi,qc,te,pe,qvenv,dummy,zt,zm
907
908 if(ifrom==1) then
909 call friction(ifrom,nm1,deltak,dt,zt,zm,t,tt ,te)
910 !call friction(ifrom,nm1,dt,zt,zm,qv,qvt,qvenv)
911 return
912 endif
913
914 dummy(:) = 0.
915 if(ifrom==2) call friction(ifrom,nm1,deltak,dt,zt,zm,w,dummy ,dummy)
916 !call friction(ifrom,nm1,dt,zt,zm,qi,qit ,dummy)
917 !call friction(ifrom,nm1,dt,zt,zm,qh,qht ,dummy)
918 !call friction(ifrom,nm1,dt,zt,zm,qc,qct ,dummy)
919 return
920 end subroutine damp_grav_wave
921 !-------------------------------------------------------------------------------
922 !
923 subroutine friction(ifrom,nm1,deltak,dt,zt,zm,var1,vart,var2)
924 implicit none
925 real, dimension(nm1) :: var1,var2,vart,zt,zm
926 integer k,nfpt,kf,nm1,ifrom,deltak
927 real zmkf,ztop,distim,c1,c2,dt
928
929 !nfpt=50
930 !kf = nm1 - nfpt
931 kf = nm1 - int(deltak/2)
932
933 zmkf = zm(kf) !old: float(kf )*dz
934 ztop = zm(nm1)
935 !distim = min(4.*dt,200.)
936 distim = 60.
937
938 c1 = 1. / (distim * (ztop - zmkf))
939 c2 = dt * c1
940
941 if(ifrom == 1) then
942 do k = nm1,2,-1
943 if (zt(k) .le. zmkf) cycle
944 vart(k) = vart(k) + c1 * (zt(k) - zmkf)*(var2(k) - var1(k))
945 enddo
946 elseif(ifrom == 2) then
947 do k = nm1,2,-1
948 if (zt(k) .le. zmkf) cycle
949 var1(k) = var1(k) + c2 * (zt(k) - zmkf)*(var2(k) - var1(k))
950 enddo
951 endif
952 return
953 end subroutine friction
954 !-------------------------------------------------------------------------------
955 !
956 subroutine vel_advectc_plumerise(m1,wc,wt,rho,dzm)
957
958 implicit none
959 integer :: k,m1
960 real, dimension(m1) :: wc,wt,flxw,dzm,rho
961 real, dimension(m1) :: dn0 ! var local
962 real :: c1z
963
964 !dzm(:)= 1./dz
965
966 dn0(1:m1)=rho(1:m1)*1.e-3 ! converte de cgs para mks
967
968 flxw(1) = wc(1) * dn0(1)
969
970 do k = 2,m1-1
971 flxw(k) = wc(k) * .5 * (dn0(k) + dn0(k+1))
972 enddo
973
974 ! Compute advection contribution to W tendency
975
976 c1z = .5
977
978 do k = 2,m1-2
979
980 wt(k) = wt(k) &
981 + c1z * dzm(k) / (dn0(k) + dn0(k+1)) * ( &
982 (flxw(k) + flxw(k-1)) * (wc(k) + wc(k-1)) &
983 - (flxw(k) + flxw(k+1)) * (wc(k) + wc(k+1)) &
984 + (flxw(k+1) - flxw(k-1)) * 2.* wc(k) )
985
986 enddo
987
988 return
989 end subroutine vel_advectc_plumerise
990 !-------------------------------------------------------------------------------
991 !
992 subroutine hadvance_plumerise(iac,m1,dt,wc,wt,wp,mintime)
993
994 implicit none
995 integer :: k,iac
996 integer :: m1,mintime
997 real, dimension(m1) :: dummy, wc,wt,wp
998 real eps,dt
999 ! It is here that the Asselin filter is applied. For the velocities
1000 ! and pressure, this must be done in two stages, the first when
1001 ! IAC=1 and the second when IAC=2.
1002
1003
1004 eps = .2
1005 if(mintime == 1) eps=0.5
1006
1007 ! For both IAC=1 and IAC=2, call PREDICT for U, V, W, and P.
1008 !
1009 call predict_plumerise(m1,wc,wp,wt,dummy,iac,2.*dt,eps)
1010 !print*,'mintime',mintime,eps
1011 !do k=1,m1
1012 ! print*,'W-HAD',k,wc(k),wp(k),wt(k)
1013 !enddo
1014 return
1015 end subroutine hadvance_plumerise
1016 !-------------------------------------------------------------------------------
1017 !
1018 subroutine predict_plumerise(npts,ac,ap,fa,af,iac,dtlp,epsu)
1019 implicit none
1020 integer :: npts,iac,m
1021 real :: epsu,dtlp
1022 real, dimension(*) :: ac,ap,fa,af
1023
1024 ! For IAC=3, this routine moves the arrays AC and AP forward by
1025 ! 1 time level by adding in the prescribed tendency. It also
1026 ! applies the Asselin filter given by:
1027
1028 ! {AC} = AC + EPS * (AP - 2 * AC + AF)
1029
1030 ! where AP,AC,AF are the past, current and future time levels of A.
1031 ! All IAC=1 does is to perform the {AC} calculation without the AF
1032 ! term present. IAC=2 completes the calculation of {AC} by adding
1033 ! the AF term only, and advances AC by filling it with input AP
1034 ! values which were already updated in ACOUSTC.
1035 !
1036
1037 if (iac .eq. 1) then
1038 do m = 1,npts
1039 ac(m) = ac(m) + epsu * (ap(m) - 2. * ac(m))
1040 enddo
1041 return
1042 elseif (iac .eq. 2) then
1043 do m = 1,npts
1044 af(m) = ap(m)
1045 ap(m) = ac(m) + epsu * af(m)
1046 enddo
1047 !elseif (iac .eq. 3) then
1048 ! do m = 1,npts
1049 ! af(m) = ap(m) + dtlp * fa(m)
1050 ! enddo
1051 ! if (ngrid .eq. 1 .and. ipara .eq. 0) call cyclic(nzp,nxp,nyp,af,'T')
1052 ! do m = 1,npts
1053 ! ap(m) = ac(m) + epsu * (ap(m) - 2. * ac(m) + af(m))
1054 ! enddo
1055 endif
1056
1057 do m = 1,npts
1058 ac(m) = af(m)
1059 enddo
1060 return
1061 end subroutine predict_plumerise
1062 !-------------------------------------------------------------------------------
1063 !
1064 subroutine buoyancy_plumerise(m1, T, TE, QV, QVENV, QH, QI, QC, WT, scr1)
1065 implicit none
1066 integer :: k,m1
1067 real, parameter :: g = 9.8, eps = 0.622, gama = 0.5 ! mass virtual coeff.
1068 real, dimension(m1) :: T, TE, QV, QVENV, QH, QI, QC, WT, scr1
1069 real :: TV,TVE,QWTOTL,umgamai
1070 real, parameter :: mu = 0.15
1071
1072 !- orig
1073 umgamai = 1./(1.+gama) ! compensa a falta do termo de aceleracao associado `as
1074 ! das pertubacoes nao-hidrostaticas no campo de pressao
1075
1076 !- new ! Siesbema et al, 2004
1077 !umgamai = 1./(1.-2.*mu)
1078
1079 do k = 2,m1-1
1080
1081 TV = T(k) * (1. + (QV(k) /EPS))/(1. + QV(k) ) !blob virtual temp.
1082 TVE = TE(k) * (1. + (QVENV(k)/EPS))/(1. + QVENV(k)) !and environment
1083
1084 QWTOTL = QH(k) + QI(k) + QC(k) ! QWTOTL*G is drag
1085 !- orig
1086 !scr1(k)= G*( umgamai*( TV - TVE) / TVE - QWTOTL)
1087 scr1(k)= G* umgamai*( (TV - TVE) / TVE - QWTOTL)
1088
1089 !if(k .lt. 10)print*,'BT',k,TV,TVE,TVE,QWTOTL
1090 enddo
1091
1092 do k = 2,m1-2
1093 wt(k) = wt(k)+0.5*(scr1(k)+scr1(k+1))
1094 ! print*,'W-BUO',k,wt(k),scr1(k),scr1(k+1)
1095 enddo
1096
1097 end subroutine buoyancy_plumerise
1098 !-------------------------------------------------------------------------------
1099 !
1100 subroutine ENTRAINMENT(m1,w,wt,radius,ALPHA)
1101 implicit none
1102 integer :: k,m1
1103 real, dimension(m1) :: w,wt,radius
1104 real DMDTM,ALPHA,WBAR,RADIUS_BAR,umgamai
1105 real, parameter :: mu = 0.15 ,gama = 0.5 ! mass virtual coeff.
1106
1107 !- new - Siesbema et al, 2004
1108 !umgamai = 1./(1.-2.*mu)
1109
1110 !- orig
1111 !umgamai = 1
1112 umgamai = 1./(1.+gama) ! compensa a falta do termo de aceleracao associado `as
1113 ! das pertubacoes nao-hidrostaticas no campo de pressao
1114
1115 !
1116 !-- ALPHA/RADIUS(L) = (1/M)DM/DZ (W 14a)
1117 do k=2,m1-1
1118
1119 !-- for W: WBAR is only W(k)
1120 ! WBAR=0.5*(W(k)+W(k-1))
1121 WBAR=W(k)
1122 RADIUS_BAR = 0.5*(RADIUS(k) + RADIUS(k+1))
1123 ! orig
1124 !DMDTM = 2. * ALPHA * ABS (WBAR) / RADIUS_BAR != (1/M)DM/DT
1125 DMDTM = umgamai * 2. * ALPHA * ABS (WBAR) / RADIUS_BAR != (1/M)DM/DT
1126
1127 !-- DMDTM*W(L) entrainment,
1128 wt(k) = wt(k) - DMDTM*ABS (WBAR)
1129 !print*,'W-ENTR=',k,w(k),- DMDTM*ABS (WBAR)
1130 enddo
1131 end subroutine ENTRAINMENT
1132 !-------------------------------------------------------------------------------
1133 !
1134 subroutine scl_advectc_plumerise(varn,mzp)
1135 !use module_zero_plumegen_coms
1136 implicit none
1137 integer :: mzp
1138 character(len=*) :: varn
1139 real :: dtlto2
1140 integer :: k
1141
1142 ! wp => w
1143 !- Advect scalars
1144 dtlto2 = .5 * dt
1145 ! vt3dc(1) = (w(1) + wc(1)) * dtlto2 * dne(1)
1146 vt3dc(1) = (w(1) + wc(1)) * dtlto2 * rho(1)*1.e-3!converte de CGS p/ MKS
1147 vt3df(1) = .5 * (w(1) + wc(1)) * dtlto2 * dzm(1)
1148
1149 do k = 2,mzp
1150 ! vt3dc(k) = (w(k) + wc(k)) * dtlto2 *.5 * (dne(k) + dne(k+1))
1151 vt3dc(k) = (w(k) + wc(k)) * dtlto2 *.5 * (rho(k) + rho(k+1))*1.e-3
1152 vt3df(k) = (w(k) + wc(k)) * dtlto2 *.5 * dzm(k)
1153 !print*,'vt3df-vt3dc',k,vt3dc(k),vt3df(k)
1154 enddo
1155
1156
1157 !-srf-24082005
1158 ! do k = 1,mzp-1
1159 do k = 1,mzp
1160 vctr1(k) = (zt(k+1) - zm(k)) * dzm(k)
1161 vctr2(k) = (zm(k) - zt(k)) * dzm(k)
1162 ! vt3dk(k) = dzt(k) / dne(k)
1163 vt3dk(k) = dzt(k) /(rho(k)*1.e-3)
1164 !print*,'VT3dk',k,dzt(k) , dne(k)
1165 enddo
1166
1167 ! scalarp => scalar_tab(n,ngrid)%var_p
1168 ! scalart => scalar_tab(n,ngrid)%var_t
1169
1170 !- temp advection tendency (TT)
1171 scr1=T
1172 call fa_zc_plumerise(mzp &
1173 ,T ,scr1 (1) &
1174 ,vt3dc (1) ,vt3df (1) &
1175 ,vt3dg (1) ,vt3dk (1) &
1176 ,vctr1,vctr2 )
1177
1178 call advtndc_plumerise(mzp,T,scr1(1),TT,dt)
1179
1180 !- water vapor advection tendency (QVT)
1181 scr1=QV
1182 call fa_zc_plumerise(mzp &
1183 ,QV ,scr1 (1) &
1184 ,vt3dc (1) ,vt3df (1) &
1185 ,vt3dg (1) ,vt3dk (1) &
1186 ,vctr1,vctr2 )
1187
1188 call advtndc_plumerise(mzp,QV,scr1(1),QVT,dt)
1189
1190 !- liquid advection tendency (QCT)
1191 scr1=QC
1192 call fa_zc_plumerise(mzp &
1193 ,QC ,scr1 (1) &
1194 ,vt3dc (1) ,vt3df (1) &
1195 ,vt3dg (1) ,vt3dk (1) &
1196 ,vctr1,vctr2 )
1197
1198 call advtndc_plumerise(mzp,QC,scr1(1),QCT,dt)
1199
1200 !- ice advection tendency (QIT)
1201 scr1=QI
1202 call fa_zc_plumerise(mzp &
1203 ,QI ,scr1 (1) &
1204 ,vt3dc (1) ,vt3df (1) &
1205 ,vt3dg (1) ,vt3dk (1) &
1206 ,vctr1,vctr2 )
1207
1208 call advtndc_plumerise(mzp,QI,scr1(1),QIT,dt)
1209
1210 !- hail/rain advection tendency (QHT)
1211 ! if(ak1 > 0. .or. ak2 > 0.) then
1212
1213 scr1=QH
1214 call fa_zc_plumerise(mzp &
1215 ,QH ,scr1 (1) &
1216 ,vt3dc (1) ,vt3df (1) &
1217 ,vt3dg (1) ,vt3dk (1) &
1218 ,vctr1,vctr2 )
1219
1220 call advtndc_plumerise(mzp,QH,scr1(1),QHT,dt)
1221 ! endif
1222
1223 return
1224
1225 !- gas/particle advection tendency (SCT)
1226 ! if(varn == 'SC')return
1227 scr1=SC
1228 call fa_zc_plumerise(mzp &
1229 ,SC ,scr1 (1) &
1230 ,vt3dc (1) ,vt3df (1) &
1231 ,vt3dg (1) ,vt3dk (1) &
1232 ,vctr1,vctr2 )
1233
1234 call advtndc_plumerise(mzp,SC,scr1(1),SCT,dt)
1235
1236
1237 return
1238 end subroutine scl_advectc_plumerise
1239 !-------------------------------------------------------------------------------
1240 !
1241 subroutine fa_zc_plumerise(m1,scp,scr1,vt3dc,vt3df,vt3dg,vt3dk,vctr1,vctr2)
1242
1243 implicit none
1244 integer :: m1,k
1245 real :: dfact
1246 real, dimension(m1) :: scp,scr1,vt3dc,vt3df,vt3dg,vt3dk
1247 real, dimension(m1) :: vctr1,vctr2
1248
1249 dfact = .5
1250
1251 ! Compute scalar flux VT3DG
1252 do k = 1,m1-1
1253 vt3dg(k) = vt3dc(k) &
1254 * (vctr1(k) * scr1(k) &
1255 + vctr2(k) * scr1(k+1) &
1256 + vt3df(k) * (scr1(k) - scr1(k+1)))
1257 enddo
1258
1259 ! Modify fluxes to retain positive-definiteness on scalar quantities.
1260 ! If a flux will remove 1/2 quantity during a timestep,
1261 ! reduce to first order flux. This will remain positive-definite
1262 ! under the assumption that ABS(CFL(i)) + ABS(CFL(i-1)) < 1.0 if
1263 ! both fluxes are evacuating the box.
1264
1265 do k = 1,m1-1
1266 if (vt3dc(k) .gt. 0.) then
1267 if (vt3dg(k) * vt3dk(k) .gt. dfact * scr1(k)) then
1268 vt3dg(k) = vt3dc(k) * scr1(k)
1269 endif
1270 elseif (vt3dc(k) .lt. 0.) then
1271 if (-vt3dg(k) * vt3dk(k+1) .gt. dfact * scr1(k+1)) then
1272 vt3dg(k) = vt3dc(k) * scr1(k+1)
1273 endif
1274 endif
1275
1276 enddo
1277
1278 ! Compute flux divergence
1279
1280 do k = 2,m1-1
1281 scr1(k) = scr1(k) &
1282 + vt3dk(k) * ( vt3dg(k-1) - vt3dg(k) &
1283 + scp (k) * ( vt3dc(k) - vt3dc(k-1)))
1284 enddo
1285 return
1286 end subroutine fa_zc_plumerise
1287 !-------------------------------------------------------------------------------
1288 !
1289 subroutine advtndc_plumerise(m1,scp,sca,sct,dtl)
1290 implicit none
1291 integer :: m1,k
1292 real :: dtl,dtli
1293 real, dimension(m1) :: scp,sca,sct
1294
1295 dtli = 1. / dtl
1296 do k = 2,m1-1
1297 sct(k) = sct(k) + (sca(k)-scp(k)) * dtli
1298 enddo
1299 return
1300 end subroutine advtndc_plumerise
1301 !-------------------------------------------------------------------------------
1302 !
1303 subroutine tend0_plumerise
1304 !use module_zero_plumegen_coms, only: nm1,wt,tt,qvt,qct,qht,qit,sct
1305 wt(1:nm1) = 0.
1306 tt(1:nm1) = 0.
1307 qvt(1:nm1) = 0.
1308 qct(1:nm1) = 0.
1309 qht(1:nm1) = 0.
1310 qit(1:nm1) = 0.
1311 !sct(1:nm1) = 0.
1312 end subroutine tend0_plumerise
1313
1314 ! ****************************************************************
1315
1316 subroutine scl_misc(m1)
1317 !use module_zero_plumegen_coms
1318 implicit none
1319 real, parameter :: g = 9.81, cp=1004.
1320 integer m1,k
1321 real dmdtm
1322
1323 do k=2,m1-1
1324 WBAR = 0.5*(W(k)+W(k-1))
1325 !-- dry adiabat
1326 ADIABAT = - WBAR * G / CP
1327 !
1328 !-- entrainment
1329 DMDTM = 2. * ALPHA * ABS (WBAR) / RADIUS (k) != (1/M)DM/DT
1330
1331 !-- tendency temperature = adv + adiab + entrainment
1332 TT(k) = TT(K) + ADIABAT - DMDTM * ( T (k) - TE (k) )
1333
1334 !-- tendency water vapor = adv + entrainment
1335 QVT(K) = QVT(K) - DMDTM * ( QV (k) - QVENV (k) )
1336
1337 QCT(K) = QCT(K) - DMDTM * ( QC (k) )
1338 QHT(K) = QHT(K) - DMDTM * ( QH (k) )
1339 QIT(K) = QIT(K) - DMDTM * ( QI (k) )
1340
1341 !-- tendency gas/particle = adv + entrainment
1342 ! SCT(K) = SCT(K) - DMDTM * ( SC (k) - SCE (k) )
1343
1344 enddo
1345 end subroutine scl_misc
1346
1347 ! ****************************************************************
1348
1349 subroutine visc_W(m1,deltak,kmt)
1350 !use module_zero_plumegen_coms
1351 implicit none
1352 integer m1,k,deltak,kmt,m2
1353 real dz1t,dz1m,dz2t,dz2m,d2wdz,d2tdz ,d2qvdz ,d2qhdz ,d2qcdz ,d2qidz ,d2scdz
1354
1355 !srf--- 17/08/2005
1356 !m2=min(m1+deltak,kmt)
1357 m2=min(m1,kmt)
1358
1359 !do k=2,m1-1
1360 do k=2,m2-1
1361 DZ1T = 0.5*(ZT(K+1)-ZT(K-1))
1362 DZ2T = VISC (k) / (DZ1T * DZ1T)
1363 DZ1M = 0.5*(ZM(K+1)-ZM(K-1))
1364 DZ2M = VISC (k) / (DZ1M * DZ1M)
1365 D2WDZ = (W (k + 1) - 2 * W (k) + W (k - 1) ) * DZ2M
1366 D2TDZ = (T (k + 1) - 2 * T (k) + T (k - 1) ) * DZ2T
1367 D2QVDZ = (QV (k + 1) - 2 * QV (k) + QV (k - 1) ) * DZ2T
1368 D2QHDZ = (QH (k + 1) - 2 * QH (k) + QH (k - 1) ) * DZ2T
1369 D2QCDZ = (QC (k + 1) - 2 * QC (k) + QC (k - 1) ) * DZ2T
1370 D2QIDZ = (QI (k + 1) - 2 * QI (k) + QI (k - 1) ) * DZ2T
1371 !D2SCDZ = (SC (k + 1) - 2 * SC (k) + SC (k - 1) ) * DZ2T
1372
1373 WT(k) = WT(k) + D2WDZ
1374 TT(k) = TT(k) + D2TDZ
1375 QVT(k) = QVT(k) + D2QVDZ
1376 QCT(k) = QCT(k) + D2QCDZ
1377 QHT(k) = QHT(k) + D2QHDZ
1378 QIT(k) = QIT(k) + D2QIDZ
1379 !SCT(k) = SCT(k) + D2SCDZ
1380 !print*,'W-VISC=',k,D2WDZ
1381 enddo
1382
1383 end subroutine visc_W
1384
1385 ! ****************************************************************
1386
1387 subroutine update_plumerise(m1,varn)
1388 !use module_zero_plumegen_coms
1389 integer m1,k
1390 character(len=*) :: varn
1391
1392 if(varn == 'W') then
1393
1394 do k=2,m1-1
1395 W(k) = W(k) + WT(k) * DT
1396 enddo
1397 return
1398
1399 else
1400 do k=2,m1-1
1401 T(k) = T(k) + TT(k) * DT
1402
1403 QV(k) = QV(k) + QVT(k) * DT
1404
1405 QC(k) = QC(k) + QCT(k) * DT !cloud drops travel with air
1406 QH(k) = QH(k) + QHT(k) * DT
1407 QI(k) = QI(k) + QIT(k) * DT
1408 ! SC(k) = SC(k) + SCT(k) * DT
1409
1410 !srf---18jun2005
1411 QV(k) = max(0., QV(k))
1412 QC(k) = max(0., QC(k))
1413 QH(k) = max(0., QH(k))
1414 QI(k) = max(0., QI(k))
1415 ! SC(k) = max(0., SC(k))
1416
1417 enddo
1418 endif
1419 end subroutine update_plumerise
1420 !-------------------------------------------------------------------------------
1421 !
1422 subroutine fallpart(m1)
1423 !use module_zero_plumegen_coms
1424 integer m1,k
1425 real vtc, dfhz,dfiz,dz1
1426 !srf==================================
1427 ! verificar se o gradiente esta correto
1428 !
1429 !srf==================================
1430 !
1431 ! XNO=1.E7 [m**-4] median volume diameter raindrop,Kessler
1432 ! VC = 38.3/(XNO**.125), median volume fallspeed eqn., Kessler
1433 ! for ice, see (OT18), use F0=0.75 per argument there. rho*q
1434 ! values are in g/m**3, velocities in m/s
1435
1436 real, PARAMETER :: VCONST = 5.107387, EPS = 0.622, F0 = 0.75
1437 real, PARAMETER :: G = 9.81, CP = 1004.
1438 !
1439 do k=2,m1-1
1440
1441 VTC = VCONST * RHO (k) **.125 ! median volume fallspeed (KTable4)
1442
1443 ! hydrometeor assembly velocity calculations (K Table4)
1444 ! VTH(k)=-VTC*QH(k)**.125 !median volume fallspeed, water
1445 VTH (k) = - 4. !small variation with qh
1446
1447 VHREL = W (k) + VTH (k) !relative to surrounding cloud
1448
1449 ! rain ventilation coefficient for evaporation
1450 CVH(k) = 1.6 + 0.57E-3 * (ABS (VHREL) ) **1.5
1451 !
1452 ! VTI(k)=-VTC*F0*QI(k)**.125 !median volume fallspeed,ice
1453 VTI (k) = - 3. !small variation with qi
1454
1455 VIREL = W (k) + VTI (k) !relative to surrounding cloud
1456 !
1457 ! ice ventilation coefficient for sublimation
1458 CVI(k) = 1.6 + 0.57E-3 * (ABS (VIREL) ) **1.5 / F0
1459 !
1460 !
1461 IF (VHREL.GE.0.0) THEN
1462 DFHZ=QH(k)*(RHO(k )*VTH(k )-RHO(k-1)*VTH(k-1))/RHO(k-1)
1463 ELSE
1464 DFHZ=QH(k)*(RHO(k+1)*VTH(k+1)-RHO(k )*VTH(k ))/RHO(k)
1465 ENDIF
1466 !
1467 !
1468 IF (VIREL.GE.0.0) THEN
1469 DFIZ=QI(k)*(RHO(k )*VTI(k )-RHO(k-1)*VTI(k-1))/RHO(k-1)
1470 ELSE
1471 DFIZ=QI(k)*(RHO(k+1)*VTI(k+1)-RHO(k )*VTI(k ))/RHO(k)
1472 ENDIF
1473
1474 DZ1=ZM(K)-ZM(K-1)
1475
1476 qht(k) = qht(k) - DFHZ / DZ1 !hydrometeors don't
1477
1478 qit(k) = qit(k) - DFIZ / DZ1 !nor does ice? hail, what about
1479
1480 enddo
1481 end subroutine fallpart
1482 !-------------------------------------------------------------------------------
1483 !-------------------------------------------------------------------------------
1484 !
1485 subroutine printout (izprint,nrectotal)
1486 !use module_zero_plumegen_coms
1487 real, parameter :: tmelt = 273.3
1488 integer, save :: nrec
1489 data nrec/0/
1490 integer :: ko,izprint,interval,nrectotal
1491 real :: pea, btmp,etmp,vap1,vap2,gpkc,gpkh,gpki,deficit
1492 interval = 1 !debug time interval,min
1493
1494 !
1495 IF (IZPRINT.EQ.0) RETURN
1496
1497 IF(MINTIME == 1) nrec = 0
1498 !
1499 WRITE (2, 430) MINTIME, DT, TIME
1500 WRITE (2, 431) ZTOP
1501 WRITE (2, 380)
1502 !
1503 ! do the print
1504 !
1505 DO 390 KO = 1, nrectotal, interval
1506
1507 PEA = PE (KO) * 10. !pressure is stored in decibars(kPa),print in mb;
1508 BTMP = T (KO) - TMELT !temps in Celsius
1509 ETMP = T (KO) - TE (KO) !temperature excess
1510 VAP1 = QV (KO) * 1000. !printout in g/kg for all water,
1511 VAP2 = QSAT (KO) * 1000. !vapor (internal storage is in g/g)
1512 GPKC = QC (KO) * 1000. !cloud water
1513 GPKH = QH (KO) * 1000. !raindrops
1514 GPKI = QI (KO) * 1000. !ice particles
1515 DEFICIT = VAP2 - VAP1 !vapor deficit
1516 !
1517 WRITE (2, 400) zt(KO)/1000., PEA, W (KO), BTMP, ETMP, VAP1, &
1518 VAP2, GPKC, GPKH, GPKI, VTH (KO), SC(KO)
1519 !
1520 !
1521 ! !end of printout
1522
1523 390 CONTINUE
1524
1525 nrec=nrec+1
1526 write (19,rec=nrec) (W (KO), KO=1,nrectotal)
1527 nrec=nrec+1
1528 write (19,rec=nrec) (T (KO), KO=1,nrectotal)
1529 nrec=nrec+1
1530 write (19,rec=nrec) (TE(KO), KO=1,nrectotal)
1531 nrec=nrec+1
1532 write (19,rec=nrec) (QV(KO)*1000., KO=1,nrectotal)
1533 nrec=nrec+1
1534 write (19,rec=nrec) (QC(KO)*1000., KO=1,nrectotal)
1535 nrec=nrec+1
1536 write (19,rec=nrec) (QH(KO)*1000., KO=1,nrectotal)
1537 nrec=nrec+1
1538 write (19,rec=nrec) (QI(KO)*1000., KO=1,nrectotal)
1539 nrec=nrec+1
1540 ! write (19,rec=nrec) (SC(KO), KO=1,nrectotal)
1541 write (19,rec=nrec) (QSAT(KO)*1000., KO=1,nrectotal)
1542 nrec=nrec+1
1543 write (19,rec=nrec) (QVENV(KO)*1000., KO=1,nrectotal)
1544
1545
1546
1547 !print*,'ntimes=',nrec/(9)
1548 !
1549 RETURN
1550 !
1551 ! ************** FORMATS *********************************************
1552 !
1553 380 FORMAT(/,' Z(KM) P(MB) W(MPS) T(C) T-TE VAP SAT QC QH' &
1554 ' QI VTH(MPS) SCAL'/)
1555 !
1556 400 FORMAT(1H , F4.1,F7.2,F7.2,F6.1,6F6.2,F7.2,1X,F6.2)
1557 !
1558 430 FORMAT(1H ,//I5,' MINUTES DT= ',F6.2,' SECONDS TIME= ' &
1559 ,F8.2,' SECONDS')
1560 431 FORMAT(' ZTOP= ',F10.2)
1561 !
1562 end subroutine printout
1563 !
1564 ! *********************************************************************
1565 SUBROUTINE WATERBAL
1566 !use module_zero_plumegen_coms
1567 !
1568
1569 IF (QC (L) .LE.1.0E-10) QC (L) = 0. !DEFEAT UNDERFLOW PROBLEM
1570 IF (QH (L) .LE.1.0E-10) QH (L) = 0.
1571 IF (QI (L) .LE.1.0E-10) QI (L) = 0.
1572 !
1573 CALL EVAPORATE !vapor to cloud,cloud to vapor
1574 !
1575 CALL SUBLIMATE !vapor to ice
1576 !
1577 CALL GLACIATE !rain to ice
1578
1579 CALL MELT !ice to rain
1580 !
1581 !if(ak1 > 0. .or. ak2 > 0.) &
1582 CALL CONVERT () !(auto)conversion and accretion
1583 !CALL CONVERT2 () !(auto)conversion and accretion
1584 !
1585
1586 RETURN
1587 END SUBROUTINE WATERBAL
1588 ! *********************************************************************
1589 SUBROUTINE EVAPORATE
1590 !
1591 !- evaporates cloud,rain and ice to saturation
1592 !
1593 !use module_zero_plumegen_coms
1594 implicit none
1595 !
1596 ! XNO=10.0E06
1597 ! HERC = 1.93*1.E-6*XN035 !evaporation constant
1598 !
1599 real, PARAMETER :: HERC = 5.44E-4, CP = 1.004, HEATCOND = 2.5E3
1600 real, PARAMETER :: HEATSUBL = 2834., TMELT = 273., TFREEZE = 269.3
1601
1602 real, PARAMETER :: FRC = HEATCOND / CP, SRC = HEATSUBL / CP
1603
1604 real :: evhdt, evidt, evrate, evap, sd, quant, dividend, divisor, devidt
1605
1606 !
1607 !
1608 SD = QSAT (L) - QV (L) !vapor deficit
1609 IF (SD.EQ.0.0) RETURN
1610 !IF (abs(SD).lt.1.e-7) RETURN
1611
1612
1613 EVHDT = 0.
1614 EVIDT = 0.
1615 !evrate =0.; evap=0.; sd=0.0; quant=0.0; dividend=0.0; divisor=0.0; devidt=0.0
1616
1617 EVRATE = ABS (WBAR * DQSDZ) !evaporation rate (Kessler 8.32)
1618 EVAP = EVRATE * DT !what we can get in DT
1619
1620
1621 IF (SD.LE.0.0) THEN ! condense. SD is negative
1622
1623 IF (EVAP.GE.ABS (SD) ) THEN !we get it all
1624
1625 QC (L) = QC (L) - SD !deficit,remember?
1626 QV (L) = QSAT(L) !set the vapor to saturation
1627 T (L) = T (L) - SD * FRC !heat gained through condensation
1628 !per gram of dry air
1629 RETURN
1630
1631 ELSE
1632
1633 QC (L) = QC (L) + EVAP !get what we can in DT
1634 QV (L) = QV (L) - EVAP !remove it from the vapor
1635 T (L) = T (L) + EVAP * FRC !get some heat
1636
1637 RETURN
1638
1639 ENDIF
1640 !
1641 ELSE !SD is positive, need some water
1642 !
1643 ! not saturated. saturate if possible. use everything in order
1644 ! cloud, rain, ice. SD is positive
1645
1646 IF (EVAP.LE.QC (L) ) THEN !enough cloud to last DT
1647 !
1648
1649 IF (SD.LE.EVAP) THEN !enough time to saturate
1650
1651 QC (L) = QC (L) - SD !remove cloud
1652 QV (L) = QSAT (L) !saturate
1653 T (L) = T (L) - SD * FRC !cool the parcel
1654 RETURN !done
1655 !
1656
1657 ELSE !not enough time
1658
1659 SD = SD-EVAP !use what there is
1660 QV (L) = QV (L) + EVAP !add vapor
1661 T (L) = T (L) - EVAP * FRC !lose heat
1662 QC (L) = QC (L) - EVAP !lose cloud
1663 !go on to rain.
1664 ENDIF
1665 !
1666 ELSE !not enough cloud to last DT
1667 !
1668 IF (SD.LE.QC (L) ) THEN !but there is enough to sat
1669
1670 QV (L) = QSAT (L) !use it
1671 QC (L) = QC (L) - SD
1672 T (L) = T (L) - SD * FRC
1673 RETURN
1674
1675 ELSE !not enough to sat
1676 SD = SD-QC (L)
1677 QV (L) = QV (L) + QC (L)
1678 T (L) = T (L) - QC (L) * FRC
1679 QC (L) = 0.0 !all gone
1680
1681 ENDIF !on to rain
1682 ENDIF !finished with cloud
1683 !
1684 ! but still not saturated, so try to use some rain
1685 ! this is tricky, because we only have time DT to evaporate. if there
1686 ! is enough rain, we can evaporate it for dt. ice can also sublimate
1687 ! at the same time. there is a compromise here.....use rain first, then
1688 ! ice. saturation may not be possible in one DT time.
1689 ! rain evaporation rate (W12),(OT25),(K Table 4). evaporate rain first
1690 ! sd is still positive or we wouldn't be here.
1691
1692
1693 IF (QH (L) .LE.1.E-10) GOTO 33
1694
1695 !srf-25082005
1696 ! QUANT = ( QC (L) + QV (L) - QSAT (L) ) * RHO (L) !g/m**3
1697 QUANT = ( QSAT (L)- QC (L) - QV (L) ) * RHO (L) !g/m**3
1698 !
1699 EVHDT = (DT * HERC * (QUANT) * (QH (L) * RHO (L) ) **.65) / RHO (L)
1700 ! rain evaporation in time DT
1701
1702 IF (EVHDT.LE.QH (L) ) THEN !enough rain to last DT
1703
1704 IF (SD.LE.EVHDT) THEN !enough time to saturate
1705 QH (L) = QH (L) - SD !remove rain
1706 QV (L) = QSAT (L) !saturate
1707 T (L) = T (L) - SD * FRC !cool the parcel
1708
1709 RETURN !done
1710 !
1711 ELSE !not enough time
1712 SD = SD-EVHDT !use what there is
1713 QV (L) = QV (L) + EVHDT !add vapor
1714 T (L) = T (L) - EVHDT * FRC !lose heat
1715 QH (L) = QH (L) - EVHDT !lose rain
1716
1717 ENDIF !go on to ice.
1718 !
1719 ELSE !not enough rain to last DT
1720 !
1721 IF (SD.LE.QH (L) ) THEN !but there is enough to sat
1722 QV (L) = QSAT (L) !use it
1723 QH (L) = QH (L) - SD
1724 T (L) = T (L) - SD * FRC
1725 RETURN
1726 !
1727 ELSE !not enough to sat
1728 SD = SD-QH (L)
1729 QV (L) = QV (L) + QH (L)
1730 T (L) = T (L) - QH (L) * FRC
1731 QH (L) = 0.0 !all gone
1732
1733 ENDIF !on to ice
1734 !
1735
1736 ENDIF !finished with rain
1737 !
1738 !
1739 ! now for ice
1740 ! equation from (OT); correction factors for units applied
1741 !
1742 33 continue
1743 IF (QI (L) .LE.1.E-10) RETURN !no ice there
1744 !
1745 DIVIDEND = ( (1.E6 / RHO (L) ) **0.475) * (SD / QSAT (L) &
1746 - 1) * (QI (L) **0.525) * 1.13
1747 DIVISOR = 7.E5 + 4.1E6 / (10. * EST (L) )
1748
1749 DEVIDT = - CVI(L) * DIVIDEND / DIVISOR !rate of change
1750
1751 EVIDT = DEVIDT * DT !what we could get
1752 !
1753 ! logic here is identical to rain. could get fancy and make subroutine
1754 ! but duplication of code is easier. God bless the screen editor.
1755 !
1756
1757 IF (EVIDT.LE.QI (L) ) THEN !enough ice to last DT
1758 !
1759
1760 IF (SD.LE.EVIDT) THEN !enough time to saturate
1761 QI (L) = QI (L) - SD !remove ice
1762 QV (L) = QSAT (L) !saturate
1763 T (L) = T (L) - SD * SRC !cool the parcel
1764
1765 RETURN !done
1766 !
1767
1768 ELSE !not enough time
1769
1770 SD = SD-EVIDT !use what there is
1771 QV (L) = QV (L) + EVIDT !add vapor
1772 T (L) = T (L) - EVIDT * SRC !lose heat
1773 QI (L) = QI (L) - EVIDT !lose ice
1774
1775 ENDIF !go on,unsatisfied
1776 !
1777 ELSE !not enough ice to last DT
1778 !
1779 IF (SD.LE.QI (L) ) THEN !but there is enough to sat
1780
1781 QV (L) = QSAT (L) !use it
1782 QI (L) = QI (L) - SD
1783 T (L) = T (L) - SD * SRC
1784
1785 RETURN
1786 !
1787 ELSE !not enough to sat
1788 SD = SD-QI (L)
1789 QV (L) = QV (L) + QI (L)
1790 T (L) = T (L) - QI (L) * SRC
1791 QI (L) = 0.0 !all gone
1792
1793 ENDIF !on to better things
1794 !finished with ice
1795 ENDIF
1796 !
1797 ENDIF !finished with the SD decision
1798 !
1799 RETURN
1800 !
1801 END SUBROUTINE EVAPORATE
1802 !
1803 ! *********************************************************************
1804 SUBROUTINE CONVERT ()
1805 !
1806 !- ACCRETION AND AUTOCONVERSION
1807 !
1808 !use module_zero_plumegen_coms
1809 !
1810 real, PARAMETER :: AK1 = 0.001 !conversion rate constant
1811 real, PARAMETER :: AK2 = 0.0052 !collection (accretion) rate
1812 real, PARAMETER :: TH = 0.5 !Kessler threshold
1813 integer, PARAMETER ::iconv = 1 !- Kessler conversion (=0)
1814
1815 !real, parameter :: ANBASE = 50.!*1.e+6 !Berry-number at cloud base #/m^3(maritime)
1816 real, parameter :: ANBASE =100000.!*1.e+6 !Berry-number at cloud base #/m^3(continental)
1817 !real, parameter :: BDISP = 0.366 !Berry--size dispersion (maritime)
1818 real, parameter :: BDISP = 0.146 !Berry--size dispersion (continental)
1819 real, parameter :: TFREEZE = 269.3 !ice formation temperature
1820 !
1821 real :: accrete, con, q, h, bc1, bc2, total
1822
1823
1824 IF (T (L) .LE. TFREEZE) RETURN !process not allowed above ice
1825 !
1826 IF (QC (L) .EQ. 0. ) RETURN
1827
1828 ACCRETE = 0.
1829 CON = 0.
1830 Q = RHO (L) * QC (L)
1831 H = RHO (L) * QH (L)
1832 !
1833 ! selection rules
1834 !
1835 !
1836 IF (QH (L) .GT. 0. ) ACCRETE = AK2 * Q * (H**.875) !accretion, Kessler
1837 !
1838 IF (ICONV.NE.0) THEN !select Berry or Kessler
1839 !
1840 !old BC1 = 120.
1841 !old BC2 = .0266 * ANBASE * 60.
1842 !old CON = BDISP * Q * Q * Q / (BC1 * Q * BDISP + BC2)
1843
1844 CON = Q*Q*Q*BDISP/(60.*(5.*Q*BDISP+0.0366*ANBASE))
1845 !
1846 ELSE
1847 !
1848 ! CON = AK1 * (Q - TH) !Kessler autoconversion rate
1849 !
1850 ! IF (CON.LT.0.0) CON = 0.0 !havent reached threshold
1851
1852 CON = max(0.,AK1 * (Q - TH)) ! versao otimizada
1853 !
1854 ENDIF
1855 !
1856 !
1857 TOTAL = (CON + ACCRETE) * DT / RHO (L)
1858
1859 !
1860 IF (TOTAL.LT.QC (L) ) THEN
1861 !
1862 QC (L) = QC (L) - TOTAL
1863 QH (L) = QH (L) + TOTAL !no phase change involved
1864 RETURN
1865 !
1866 ELSE
1867 !
1868 QH (L) = QH (L) + QC (L) !uses all there is
1869 QC (L) = 0.0
1870 !
1871 ENDIF
1872 !
1873 RETURN
1874 !
1875 END SUBROUTINE CONVERT
1876 !
1877 !**********************************************************************
1878 !
1879 SUBROUTINE CONVERT2 ()
1880 !use module_zero_plumegen_coms
1881 implicit none
1882 LOGICAL AEROSOL
1883 parameter(AEROSOL=.true.)
1884 !
1885 real, parameter :: TNULL=273.16, LAT=2.5008E6 &
1886 ,EPSI=0.622 ,DB=1. ,NB=1500. !ALPHA=0.2
1887 real :: KA,KEINS,KZWEI,KDREI,VT
1888 real :: A,B,C,D, CON,ACCRETE,total
1889
1890 real Y(6),ROH
1891
1892 A=0.
1893 B=0.
1894 Y(1) = T(L)
1895 Y(4) = W(L)
1896 y(2) = QC(L)
1897 y(3) = QH(L)
1898 Y(5) = RADIUS(L)
1899 ROH = RHO(L)*1.e-3 ! dens (MKS) ??
1900
1901
1902 ! autoconversion
1903
1904 KA = 0.0005
1905 IF( Y(1) .LT. 258.15 )THEN
1906 ! KEINS=0.00075
1907 KEINS=0.0009
1908 KZWEI=0.0052
1909 KDREI=15.39
1910 ELSE
1911 KEINS=0.0015
1912 KZWEI=0.00696
1913 KDREI=11.58
1914 ENDIF
1915
1916 ! ROH=PE/RD/TE
1917 VT=-KDREI* (Y(3)/ROH)**0.125
1918
1919
1920 IF (Y(4).GT.0.0 ) THEN
1921 IF (AEROSOL) THEN
1922 A = 1/y(4) * y(2)*y(2)*1000./( 60. *( 5. + 0.0366*NB/(y(2)*1000.*DB) ) )
1923 ELSE
1924 IF (y(2).GT.(KA*ROH)) THEN
1925 !print*,'1',y(2),KA*ROH
1926 A = KEINS/y(4) *(y(2) - KA*ROH )
1927 ENDIF
1928 ENDIF
1929 ELSE
1930 A = 0.0
1931 ENDIF
1932
1933 ! accretion
1934
1935 IF(y(4).GT.0.0) THEN
1936 B = KZWEI/(y(4) - VT) * MAX(0.,y(2)) * &
1937 MAX(0.001,ROH)**(-0.875)*(MAX(0.,y(3)))**(0.875)
1938 ELSE
1939 B = 0.0
1940 ENDIF
1941
1942
1943 !PSATW=610.7*EXP( 17.25 *( Y(1) - TNULL )/( Y(1)-36. ) )
1944 !PSATE=610.7*EXP( 22.33 *( Y(1) - TNULL )/( Y(1)- 2. ) )
1945
1946 !QSATW=EPSI*PSATW/( PE-(1.-EPSI)*PSATW )
1947 !QSATE=EPSI*PSATE/( PE-(1.-EPSI)*PSATE )
1948
1949 !MU=2.*ALPHA/Y(5)
1950
1951 !C = MU*( ROH*QSATW - ROH*QVE + y(2) )
1952 !D = ROH*LAT*QSATW*EPSI/Y1/Y1/RD *DYDX1
1953
1954
1955 !DYDX(2) = - A - B - C - D ! d rc/dz
1956 !DYDX(3) = A + B ! d rh/dz
1957
1958
1959 ! rc=rc+dydx(2)*dz
1960 ! rh=rh+dydx(3)*dz
1961
1962 CON = A
1963 ACCRETE = B
1964
1965 TOTAL = (CON + ACCRETE) *(1/DZM(L)) /ROH ! DT / RHO (L)
1966
1967 !print*,'L=',L,total,QC(L),dzm(l)
1968
1969 !
1970 IF (TOTAL.LT.QC (L) ) THEN
1971 !
1972 QC (L) = QC (L) - TOTAL
1973 QH (L) = QH (L) + TOTAL !no phase change involved
1974 RETURN
1975 !
1976 ELSE
1977 !
1978 QH (L) = QH (L) + QC (L) !uses all there is
1979 QC (L) = 0.0
1980 !
1981 ENDIF
1982 !
1983 RETURN
1984 !
1985 END SUBROUTINE CONVERT2
1986 ! ice - effect on temperature
1987 ! TTD = 0.0
1988 ! TTE = 0.0
1989 ! CALL ICE(QSATW,QSATE,Y(1),Y(2),Y(3), &
1990 ! TTA,TTB,TTC,DZ,ROH,D,C,TTD,TTE)
1991 ! DYDX(1) = DYDX(1) + TTD + TTE ! DT/DZ on Temp
1992 !
1993 !**********************************************************************
1994 !
1995 SUBROUTINE SUBLIMATE
1996 !
1997 ! ********************* VAPOR TO ICE (USE EQUATION OT22)***************
1998 !use module_zero_plumegen_coms
1999 !
2000 real, PARAMETER :: EPS = 0.622, HEATFUS = 334., HEATSUBL = 2834., CP = 1.004
2001 real, PARAMETER :: SRC = HEATSUBL / CP, FRC = HEATFUS / CP, TMELT = 273.3
2002 real, PARAMETER :: TFREEZE = 269.3
2003
2004 real ::dtsubh, dividend,divisor, subl
2005 !
2006 DTSUBH = 0.
2007 !
2008 !selection criteria for sublimation
2009 IF (T (L) .GT. TFREEZE ) RETURN
2010 IF (QV (L) .LE. QSAT (L) ) RETURN
2011 !
2012 ! from (OT); correction factors for units applied
2013 !
2014 DIVIDEND = ( (1.E6 / RHO (L) ) **0.475) * (QV (L) / QSAT (L) &
2015 - 1) * (QI (L) **0.525) * 1.13
2016 DIVISOR = 7.E5 + 4.1E6 / (10. * EST (L) )
2017 !
2018
2019 DTSUBH = ABS (DIVIDEND / DIVISOR) !sublimation rate
2020 SUBL = DTSUBH * DT !and amount possible
2021 !
2022 ! again check the possibilities
2023 !
2024 IF (SUBL.LT.QV (L) ) THEN
2025 !
2026 QV (L) = QV (L) - SUBL !lose vapor
2027 QI (L) = QI (L) + SUBL !gain ice
2028 T (L) = T (L) + SUBL * SRC !energy change, warms air
2029
2030 !print*,'5',l,qi(l),SUBL
2031
2032 RETURN
2033 !
2034 ELSE
2035 !
2036 QI (L) = QV (L) !use what there is
2037 T (L) = T (L) + QV (L) * SRC !warm the air
2038 QV (L) = 0.0
2039 !print*,'6',l,qi(l)
2040 !
2041 ENDIF
2042 !
2043 RETURN
2044 END SUBROUTINE SUBLIMATE
2045 !
2046 ! *********************************************************************
2047 !
2048 SUBROUTINE GLACIATE
2049 !
2050 ! *********************** CONVERSION OF RAIN TO ICE *******************
2051 ! uses equation OT 16, simplest. correction from W not applied, but
2052 ! vapor pressure differences are supplied.
2053 !
2054 !use module_zero_plumegen_coms
2055 !
2056 real, PARAMETER :: HEATFUS = 334., CP = 1.004, EPS = 0.622, HEATSUBL = 2834.
2057 real, PARAMETER :: FRC = HEATFUS / CP, FRS = HEATSUBL / CP, TFREEZE = 269.3
2058 real, PARAMETER :: GLCONST = 0.025 !glaciation time constant, 1/sec
2059 real dfrzh
2060 !
2061
2062 DFRZH = 0. !rate of mass gain in ice
2063 !
2064 !selection rules for glaciation
2065 IF (QH (L) .LE. 0. ) RETURN
2066 IF (QV (L) .LT. QSAT (L) ) RETURN
2067 IF (T (L) .GT. TFREEZE ) RETURN
2068 !
2069 ! NT=TMELT-T(L)
2070 ! IF (NT.GT.50) NT=50
2071 !
2072
2073 DFRZH = DT * GLCONST * QH (L) ! from OT(16)
2074 !
2075 IF (DFRZH.LT.QH (L) ) THEN
2076 !
2077 QI (L) = QI (L) + DFRZH
2078 QH (L) = QH (L) - DFRZH
2079 T (L) = T (L) + FRC * DFRZH !warms air
2080
2081 !print*,'7',l,qi(l),DFRZH
2082
2083
2084 RETURN
2085 !
2086 ELSE
2087 !
2088 QI (L) = QI (L) + QH (L)
2089 T (L) = T (L) + FRC * QH (L)
2090 QH (L) = 0.0
2091
2092 !print*,'8',l,qi(l), QH (L)
2093 !
2094 ENDIF
2095 !
2096 RETURN
2097 !
2098 END SUBROUTINE GLACIATE
2099 !
2100 !
2101 ! *********************************************************************
2102 SUBROUTINE MELT
2103 !
2104 ! ******************* MAKES WATER OUT OF ICE **************************
2105 !use module_zero_plumegen_coms
2106 !
2107 real, PARAMETER :: FRC = 332.27, TMELT = 273., F0 = 0.75 !ice velocity factor
2108 real DTMELT
2109 !
2110 DTMELT = 0. !conversion,ice to rain
2111 !
2112 !selection rules
2113 IF (QI (L) .LE. 0.0 ) RETURN
2114 IF (T (L) .LT. TMELT) RETURN
2115 !
2116 !OT(23,24)
2117 DTMELT = DT * (2.27 / RHO (L) ) * CVI(L) * (T (L) - TMELT) * ( (RHO(L) &
2118 * QI (L) * 1.E-6) **0.525) * (F0** ( - 0.42) )
2119 !after Mason,1956
2120 !
2121 ! check the possibilities
2122 !
2123 IF (DTMELT.LT.QI (L) ) THEN
2124 !
2125 QH (L) = QH (L) + DTMELT
2126 QI (L) = QI (L) - DTMELT
2127 T (L) = T (L) - FRC * DTMELT !cools air
2128 !print*,'9',l,qi(l),DTMELT
2129
2130
2131 RETURN
2132 !
2133 ELSE
2134 !
2135 QH (L) = QH (L) + QI (L) !get all there is to get
2136 T (L) = T (L) - FRC * QI (L)
2137 QI (L) = 0.0
2138 !print*,'10',l,qi(l)
2139 !
2140 ENDIF
2141 !
2142 RETURN
2143 !
2144 END SUBROUTINE MELT
2145
2146 SUBROUTINE htint (nzz1, vctra, eleva, nzz2, vctrb, elevb)
2147 IMPLICIT NONE
2148 INTEGER, INTENT(IN ) :: nzz1
2149 INTEGER, INTENT(IN ) :: nzz2
2150 REAL, INTENT(IN ) :: vctra(nzz1)
2151 REAL, INTENT(OUT) :: vctrb(nzz2)
2152 REAL, INTENT(IN ) :: eleva(nzz1)
2153 REAL, INTENT(IN ) :: elevb(nzz2)
2154
2155 INTEGER :: l
2156 INTEGER :: k
2157 INTEGER :: kk
2158 REAL :: wt
2159
2160 l=1
2161
2162 DO k=1,nzz2
2163 DO
2164 IF ( (elevb(k) < eleva(1)) .OR. &
2165 ((elevb(k) >= eleva(l)) .AND. (elevb(k) <= eleva(l+1))) ) THEN
2166 wt = (elevb(k)-eleva(l))/(eleva(l+1)-eleva(l))
2167 vctrb(k) = vctra(l)+(vctra(l+1)-vctra(l))*wt
2168 EXIT
2169 ELSE IF ( elevb(k) > eleva(nzz1)) THEN
2170 wt = (elevb(k)-eleva(nzz1))/(eleva(nzz1-1)-eleva(nzz1))
2171 vctrb(k) = vctra(nzz1)+(vctra(nzz1-1)-vctra(nzz1))*wt
2172 EXIT
2173 END IF
2174
2175 l=l+1
2176 IF(l == nzz1) THEN
2177 PRINT *,'htint:nzz1',nzz1
2178 DO kk=1,l
2179 PRINT*,'kk,eleva(kk),elevb(kk)',eleva(kk),elevb(kk)
2180 END DO
2181 STOP 'htint'
2182 END IF
2183 END DO
2184 END DO
2185 END SUBROUTINE htint
2186 !-----------------------------------------------------------------------------
2187 FUNCTION ESAT_PR (TEM)
2188 !
2189 ! ******* Vapor Pressure A.L. Buck JAM V.20 p.1527. (1981) ***********
2190 !
2191 real, PARAMETER :: CI1 = 6.1115, CI2 = 22.542, CI3 = 273.48
2192 real, PARAMETER :: CW1 = 6.1121, CW2 = 18.729, CW3 = 257.87, CW4 = 227.3
2193 real, PARAMETER :: TMELT = 273.3
2194
2195 real ESAT_PR
2196 real temc , tem,esatm
2197 !
2198 ! formulae from Buck, A.L., JAM 20,1527-1532
2199 ! custom takes esat wrt water always. formula for h2o only
2200 ! good to -40C so:
2201 !
2202 !
2203 TEMC = TEM - TMELT
2204 IF (TEMC.GT. - 40.0) GOTO 230
2205 ESATM = CI1 * EXP (CI2 * TEMC / (TEMC + CI3) ) !ice, millibars
2206 ESAT_PR = ESATM / 10. !kPa
2207
2208 RETURN
2209 !
2210 230 ESATM = CW1 * EXP ( ( (CW2 - (TEMC / CW4) ) * TEMC) / (TEMC + CW3))
2211
2212 ESAT_PR = ESATM / 10. !kPa
2213 RETURN
2214 END function ESAT_PR
2215 ! ******************************************************************
2216
2217 ! ------------------------------------------------------------------------
2218 END Module module_chem_plumerise_scalar