progomega_calc.f90

      module progomega
 
        implicit none
 
        public progomega_calc
 
      contains
 
 
       
   subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegain,delt,del, &
        zi,cnvflg,omegaout,grav,buo,drag,wush,tentr,bb1,bb2)
     
     use machine,  only : kind_phys
     use funcphys, only : fpvs  
     implicit none
 
     integer, intent(in)  :: im, km
     integer, intent(in)  :: kbcon1(im),ktcon(im)
     real(kind=kind_phys), intent(in)  :: delt,grav,bb1,bb2
     real(kind=kind_phys), intent(in)  :: omegain(im,km), del(im,km),zi(im,km)
     real(kind=kind_phys), intent(in)  :: drag(im,km),buo(im,km),wush(im,km),tentr(im,km)
     real(kind=kind_phys), intent(inout) :: omegaout(im,km)
     logical, intent(in)               :: cnvflg(im),first_time_step,flag_restart
     real(kind=kind_phys) :: terma(im,km),termb(im,km),termc(im,km),omega(im,km)
     real(kind=kind_phys) :: rhs(im,km),kd(im,km)
     real(kind=kind_phys) :: dp,dz,entrn,kdn,discr,wush_pa,lbb1,lbb2,lbb3
     integer              :: i,k
 
     entrn = 0.8e-4 !0.5E-4 !m^-1
     kdn   = 0.5e-4 !2.9E-4 !m^-1
     lbb1  = 0.5 !1.0 
     lbb2  = 3.2 !3.0
     lbb3  = 0.5 !0.5
     
     
     !Initialization 2D
     do k = 1,km
        do i = 1,im
           terma(i,k)=0.
           termb(i,k)=0.
           termc(i,k)=0.
           rhs(i,k)=0.
           omega(i,k)=omegain(i,k)
        enddo
     enddo
 
     if(first_time_step .and. .not. flag_restart)then
        do k = 1,km
           do i = 1,im
              if(cnvflg(i))then
                 omega(i,k)=-1.2 !Pa/s 
              endif
           enddo
        enddo
     endif
     
     ! Compute RHS terms
     !Lisa Bengtsson: !  compute updraft velocity omega (Pa/s)
     
     do k = 2, km
        do i = 1, im
           if (cnvflg(i)) then
              if (k > kbcon1(i) .and. k < ktcon(i)) then
 
                 ! Aerodynamic drag parameter
                 kd(i,k) = (tentr(i,k)/entrn)*kdn
 
                 ! Scale by dp/dz to have equation in Pa/s
                 !(dp/dz > 0)
                 dp = 1000. * del(i,k)
                 dz = zi(i,k+1) - zi(i,k)
                 
                 !termA - Ensures quadratic damping (drag).
                 !termB - Ensures linear damping from wind shear.
                 !termC - Adds buoyancy forcing 
                 
                 !Coefficients for the quadratic equation
                 terma(i,k) = delt * ((lbb1 * drag(i,k) * (dp/dz)) + (kd(i,k) * (dp/dz)))
                 termb(i,k) = -1.0 - delt * lbb3 * wush(i,k) * dp/dz
                 termc(i,k) = omega(i,k) - delt * lbb2 * buo(i,k) * (dp/dz) &
                      - delt * omega(i,k) * (omega(i,k-1) - omega(i,k)) / dp
                 !Compute the discriminant
                 discr = termb(i,k)**2 - 4.0 * terma(i,k) * termc(i,k)
 
                 ! Check if discriminant is non-negative
                 if (discr >= 0.0) then
                 ! Solve quadratic equation, take the negative root
                 omegaout(i,k) = (-termb(i,k) - sqrt(discr)) / (2.0 * terma(i,k))
                 else
                 omegaout(i,k) = omega(i,k)
                 endif
 
                 omegaout(i,k) = max(min(omegaout(i,k), -1.2), -80.0)
                
              endif
           endif
        enddo
     enddo
     
    end subroutine progomega_calc
end module progomega