/[MITgcm]/MITgcm/pkg/atm_phys/dargan_bettsmiller_mod.F90
ViewVC logotype

Diff of /MITgcm/pkg/atm_phys/dargan_bettsmiller_mod.F90

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by jmc, Wed May 8 22:14:14 2013 UTC revision 1.2 by jmc, Wed May 14 21:39:18 2014 UTC
# Line 16  module dargan_bettsmiller_mod Line 16  module dargan_bettsmiller_mod
16  !   consistent with conventional definition of mixing ratio  !   consistent with conventional definition of mixing ratio
17  ! - CAPE/CIN calculation uses virtual temperature if option set  ! - CAPE/CIN calculation uses virtual temperature if option set
18    
   
19  !----------------------------------------------------------------------  !----------------------------------------------------------------------
20  !use            fms_mod, only:  file_exist, error_mesg, open_file,  &  !use            fms_mod, only:  file_exist, error_mesg, open_file,  &
21  !                               check_nml_error, mpp_pe, FATAL,  &  !                               check_nml_error, mpp_pe, FATAL,  &
# Line 126  contains Line 125  contains
125  !           qdel     specific humidity tendency (of water vapor) at  !           qdel     specific humidity tendency (of water vapor) at
126  !                      full model levels  !                      full model levels
127  !           bmflag   flag for which routines you're calling  !           bmflag   flag for which routines you're calling
128  !           klzbs    stored klzb values  !           klzbs    stored (integer part) klzb and ktop (decimal part) values
129  !           cape     convectively available potential energy  !           cape     convectively available potential energy
130  !           cin      convective inhibition (this and the above are before the  !           cin      convective inhibition (this and the above are before the
131  !                    adjustment)  !                    adjustment)
# Line 150  contains Line 149  contains
149  !-----------------------------------------------------------------------  !-----------------------------------------------------------------------
150  !---------------------- local data -------------------------------------  !---------------------- local data -------------------------------------
151    
152  logical,dimension(size(tin,1),size(tin,2),size(tin,3)) :: do_adjust  !logical,dimension(size(tin,1),size(tin,2),size(tin,3)) :: do_adjust
153     real,dimension(size(tin,1),size(tin,2),size(tin,3)) ::  &     real,dimension(size(tin,1),size(tin,2),size(tin,3)) :: rin
154               rin, esat, qsat, desat, dqsat, pmes, pmass     real,dimension(size(tin,1),size(tin,2))             :: precip, precip_t
155     real,dimension(size(tin,1),size(tin,2))             ::  &     real,dimension(size(tin,3))                         :: eref, rpc, tpc
                      hlcp, precip, precip_t  
    real,dimension(size(tin,3))                         :: eref, rpc, tpc, &  
                                                           tpc1, rpc1  
156    
157     real                                                ::  &     real                                                ::  &
158         cape1, cin1, tot, deltak, deltaq, qrefint, deltaqfrac, deltaqfrac2, &         cape1, cin1, tot, deltak, deltaq, qrefint, deltaqfrac, deltaqfrac2, &
159         ptopfrac, es, capeflag1, plzb, plcl, cape2, small         ptopfrac, es, capeflag1, small
160   integer  i, j, k, ix, jx, kx, klzb, ktop, klzb2   integer  i, j, k, ix, jx, kx, klzb, ktop
161  !-----------------------------------------------------------------------  !-----------------------------------------------------------------------
162  !     computation of precipitation by betts-miller scheme  !     computation of precipitation by betts-miller scheme
163  !-----------------------------------------------------------------------  !-----------------------------------------------------------------------
# Line 175  logical,dimension(size(tin,1),size(tin,2 Line 171  logical,dimension(size(tin,1),size(tin,2
171        kx=size(tin,3)        kx=size(tin,3)
172        small = 1.e-10        small = 1.e-10
173    
174    ! initialise output:
175          precip = 0.
176          tdel   = 0.
177          qdel   = 0.
178          t_ref  = tin
179          q_ref  = qin
180          bmflag = 0.
181          klzbs  = 0.
182          cape   = 0.
183          cin    = 0.
184          invtau_bm_t = 0.
185          invtau_bm_q = 0.
186          precip_t = 0.
187    
188  ! calculate r (where r is the mixing ratio)  ! calculate r (where r is the mixing ratio)
189         rin = qin/(1.0 - qin)         rin = qin/(1.0 - qin)
190    
191         do i=1,ix         do j=1,jx
192            do j=1,jx            do i=1,ix
193               cape1 = 0.               cape1 = 0.
194               cin1 = 0.               cin1 = 0.
195               tot = 0.               tot = 0.
# Line 207  logical,dimension(size(tin,1),size(tin,2 Line 217  logical,dimension(size(tin,1),size(tin,2
217  !             if((tot.gt.0.).and.(cape1.gt.0.)) then  !             if((tot.gt.0.).and.(cape1.gt.0.)) then
218                  bmflag(i,j) = 1.                  bmflag(i,j) = 1.
219  ! reference temperature is just that of the parcel all the way up  ! reference temperature is just that of the parcel all the way up
220                  t_ref(i,j,:) = tpc  !               t_ref(i,j,:) = tpc
221                    t_ref(i,j,klzb:kx) = tpc(klzb:kx)
222                  do k=klzb,kx                  do k=klzb,kx
223  ! sets reference spec hum to a certain relative hum (change to vapor pressure,  ! sets reference spec hum to a certain relative hum (change to vapor pressure,
224  ! multiply by rhbm, then back to spec humid)  ! multiply by rhbm, then back to spec humid)
# Line 226  logical,dimension(size(tin,1),size(tin,2 Line 237  logical,dimension(size(tin,1),size(tin,2
237  ! set the reference profiles to be the original profiles (for diagnostic  ! set the reference profiles to be the original profiles (for diagnostic
238  ! purposes only --  you can think of this as what you're relaxing to in  ! purposes only --  you can think of this as what you're relaxing to in
239  ! areas above the actual convection  ! areas above the actual convection
240                  do k=1,max(klzb-1,1)  !               do k=1,max(klzb-1,1)
241                     qdel(i,j,k) = 0.0  !                  qdel(i,j,k) = 0.0
242                     tdel(i,j,k) = 0.0  !                  tdel(i,j,k) = 0.0
243                     q_ref(i,j,k) = qin(i,j,k)  !                  q_ref(i,j,k) = qin(i,j,k)
244                     t_ref(i,j,k) = tin(i,j,k)  !                  t_ref(i,j,k) = tin(i,j,k)
245                  end do  !               end do
246  ! initialize p to zero for the loop  ! initialize p to zero for the loop
247                  precip(i,j) = 0.                  precip(i,j) = 0.
248                  precip_t(i,j) = 0.                  precip_t(i,j) = 0.
# Line 247  logical,dimension(size(tin,1),size(tin,2 Line 258  logical,dimension(size(tin,1),size(tin,2
258                     qdel(i,j,k) = - (qin(i,j,k) - q_ref(i,j,k))/tau_bm*dt                     qdel(i,j,k) = - (qin(i,j,k) - q_ref(i,j,k))/tau_bm*dt
259  ! Precipitation can be calculated already, based on the change in q on the  ! Precipitation can be calculated already, based on the change in q on the
260  ! way up (this isn't altered in the energy conservation scheme).  ! way up (this isn't altered in the energy conservation scheme).
261                     precip(i,j) = precip(i,j) - qdel(i,j,k)*(phalf(i,j,k+1)- &                     precip(i,j)  = precip(i,j) - qdel(i,j,k)           &
262                                   phalf(i,j,k))/grav                                  *(phalf(i,j,k+1)- phalf(i,j,k))/grav
263                     precip_t(i,j)= precip_t(i,j) + cp_air/(hlv+small)*tdel(i,j,k)* &                     precip_t(i,j)= precip_t(i,j) + cp_air/(hlv+small)*tdel(i,j,k) &
264                                   (phalf(i,j,k+1)-phalf(i,j,k))/grav                                  *(phalf(i,j,k+1)-phalf(i,j,k))/grav
265                  end do                  end do
266                  if ((precip(i,j).gt.0.).and.(precip_t(i,j).gt.0.)) then                  if ((precip(i,j).gt.0.).and.(precip_t(i,j).gt.0.)) then
267  ! If precip > 0, then correct energy.  ! If precip > 0, then correct energy.
# Line 302  logical,dimension(size(tin,1),size(tin,2 Line 313  logical,dimension(size(tin,1),size(tin,2
313  ! once we finish this), the actual new top of convection is somewhere between  ! once we finish this), the actual new top of convection is somewhere between
314  ! the current ktop, and one level above this.  set ktop to the level above.  ! the current ktop, and one level above this.  set ktop to the level above.
315                        ktop = ktop - 1                        ktop = ktop - 1
316                          klzbs(i,j) = klzbs(i,j) + float(ktop)/( kx + 1.d0 )
317  ! Adjust the tendencies in the places above back to zero, and the reference  ! Adjust the tendencies in the places above back to zero, and the reference
318  ! profiles back to the original t,q.  ! profiles back to the original t,q.
319                        if (ktop.gt.klzb) then                        if (ktop.gt.klzb) then
320                           qdel(i,j,klzb:ktop-1) = 0.                           qdel(i,j,klzb:ktop-1) = 0.
321                           q_ref(i,j,klzb:ktop-1) = qin(i,j,klzb:ktop-1)  !                        q_ref(i,j,klzb:ktop-1) = qin(i,j,klzb:ktop-1)
322                           tdel(i,j,klzb:ktop-1) = 0.                           tdel(i,j,klzb:ktop-1) = 0.
323                           t_ref(i,j,klzb:ktop-1) = tin(i,j,klzb:ktop-1)  !                        t_ref(i,j,klzb:ktop-1) = tin(i,j,klzb:ktop-1)
324                        end if                        end if
325  ! Then make the change only a fraction of the new top layer so the precip is  ! Then make the change only a fraction of the new top layer so the precip is
326  ! identically zero.  ! identically zero.
# Line 356  logical,dimension(size(tin,1),size(tin,2 Line 368  logical,dimension(size(tin,1),size(tin,2
368                        else                        else
369                           precip(i,j) = 0.                           precip(i,j) = 0.
370                           qdel(i,j,kx) = 0.                           qdel(i,j,kx) = 0.
371                           q_ref(i,j,kx) = qin(i,j,kx)  !                        q_ref(i,j,kx) = qin(i,j,kx)
372                           tdel(i,j,kx) = 0.                           tdel(i,j,kx) = 0.
373                           t_ref(i,j,kx) = tin(i,j,kx)  !                        t_ref(i,j,kx) = tin(i,j,kx)
374                           invtau_bm_t(i,j) = 0.  !                        invtau_bm_t(i,j) = 0.
375                           invtau_bm_q(i,j) = 0.  !                        invtau_bm_q(i,j) = 0.
376                        end if                        end if
377                     else if(do_changeqref) then                     else if(do_changeqref) then
378  ! Change the reference profile of q by a certain fraction so that precip is  ! Change the reference profile of q by a certain fraction so that precip is
# Line 404  logical,dimension(size(tin,1),size(tin,2 Line 416  logical,dimension(size(tin,1),size(tin,2
416                        precip(i,j) = 0.                        precip(i,j) = 0.
417                        tdel(i,j,:) = 0.                        tdel(i,j,:) = 0.
418                        qdel(i,j,:) = 0.                        qdel(i,j,:) = 0.
419                        invtau_bm_t(i,j) = 0.  !                     invtau_bm_t(i,j) = 0.
420                        invtau_bm_q(i,j) = 0.  !                     invtau_bm_q(i,j) = 0.
421                     end if                     end if
422  ! for cases where virtual temp predicts CAPE but precip_t < 0.  ! for cases where virtual temp predicts CAPE but precip_t < 0.
423  ! - also note cape and precip_t are different because cape  ! - also note cape and precip_t are different because cape
# Line 414  logical,dimension(size(tin,1),size(tin,2 Line 426  logical,dimension(size(tin,1),size(tin,2
426                     tdel(i,j,:) = 0.0                     tdel(i,j,:) = 0.0
427                     qdel(i,j,:) = 0.0                     qdel(i,j,:) = 0.0
428                     precip(i,j) = 0.0                     precip(i,j) = 0.0
429                     q_ref(i,j,:) = qin(i,j,:)  !                  q_ref(i,j,:) = qin(i,j,:)
430                     t_ref(i,j,:) = tin(i,j,:)  !                  t_ref(i,j,:) = tin(i,j,:)
431                     invtau_bm_t(i,j) = 0.  !                  invtau_bm_t(i,j) = 0.
432                     invtau_bm_q(i,j) = 0.  !                  invtau_bm_q(i,j) = 0.
433                  end if                  end if
434  ! if no CAPE, set tendencies to zero.  ! if no CAPE, set tendencies to zero.
435               else  !            else
436                  tdel(i,j,:) = 0.0  !               tdel(i,j,:) = 0.0
437                  qdel(i,j,:) = 0.0  !               qdel(i,j,:) = 0.0
438                  precip(i,j) = 0.0  !               precip(i,j) = 0.0
439                  q_ref(i,j,:) = qin(i,j,:)  !               q_ref(i,j,:) = qin(i,j,:)
440                  t_ref(i,j,:) = tin(i,j,:)  !               t_ref(i,j,:) = tin(i,j,:)
441                  invtau_bm_t(i,j) = 0.  !               invtau_bm_t(i,j) = 0.
442                  invtau_bm_q(i,j) = 0.  !               invtau_bm_q(i,j) = 0.
443               end if               end if
444            end do            end do
445         end do         end do
446    
447         rain = precip         rain = precip
448         snow = 0.         snow = 0.
449    !      snow = precip_t  !to output value of precip_t (debug)
450    
451     end subroutine dargan_bettsmiller     end subroutine dargan_bettsmiller
452    
# Line 485  logical,dimension(size(tin,1),size(tin,2 Line 497  logical,dimension(size(tin,1),size(tin,2
497        real, intent(out)                      :: cape, cin        real, intent(out)                      :: cape, cin
498        integer, intent(in)                    :: i,j,bi,bj,myIter        integer, intent(in)                    :: i,j,bi,bj,myIter
499        integer, intent(in)                    :: myThid        integer, intent(in)                    :: myThid
500        integer            :: k, klcl, klfc, klcl2        integer            :: k, klcl, klfc
501        logical            :: nocape        logical            :: nocape
502        real, dimension(kx)   :: theta, tin_virtual        real, dimension(kx)   :: tin_virtual
503        real                  :: t0, r0, es, rs, theta0, pstar, value, tlcl, &        real                  :: t0, r0, es, rs, theta0, pstar, value, tlcl, &
504                                 a, b, dtdlnp, d2tdlnp2, thetam, rm, tlcl2, &                                 a, b, dtdlnp, &
505                                 plcl2, plcl, plzb, small                                 plcl, plzb, small
506    
507        pstar = 1.e5        pstar = 1.e5
508  ! so we can run dry limit (one expression involves 1/hlv)  ! so we can run dry limit (one expression involves 1/hlv)
# Line 767  logical,dimension(size(tin,1),size(tin,2 Line 779  logical,dimension(size(tin,1),size(tin,2
779        v1 = 10.*v1        v1 = 10.*v1
780        tlcl = (v2 + 1.0 - v1)*lcltable(ival+1) + (v1 - v2)*lcltable(ival+2)        tlcl = (v2 + 1.0 - v1)*lcltable(ival+1) + (v1 - v2)*lcltable(ival+2)
781    
   
782        end subroutine lcltabl        end subroutine lcltabl
783    
784  !#######################################################################  !#######################################################################
# Line 780  logical,dimension(size(tin,1),size(tin,2 Line 791  logical,dimension(size(tin,1),size(tin,2
791  !  !
792  !-----------------------------------------------------------------------  !-----------------------------------------------------------------------
793    
794    integer  unit,io,ierr  ! integer  unit,io,ierr
795    integer, intent(in) ::myThid    integer, intent(in) ::myThid
796  !-------------------------------------------------------------------------------------  !-------------------------------------------------------------------------------------
797  integer, dimension(3) :: half = (/1,2,4/)  !integer, dimension(3) :: half = (/1,2,4/)
798  !integer :: ierr, io  !integer :: ierr, io
799  integer         :: iUnit  integer         :: iUnit
800  CHARACTER*(gcm_LEN_MBUF) :: msgBuf  CHARACTER*(gcm_LEN_MBUF) :: msgBuf
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22