chem/module_chem_plumerise_scalar.F

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