dyn_em/module_convtrans_prep.F

   1 !------------------------------------------------------------------------
   2 ! original routinte by Georg Grell, adaptive timestepping by William Gustafson Jr. (PNNL), cloud fraction by Georg Grell (based on Stu's previous work wih so2so4 routine
   3 !------------------------------------------------------------------------
   4
   5 Module module_convtrans_prep
   6 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   7 CONTAINS
   8 subroutine convtrans_prep(gd_cloud,gd_cloud2,gd_cloud_a,         &
   9      gd_cloud_b,raincv,raincv_a,raincv_b,                        &
  10      cldfr,moist,p_QV,p_QC,p_qi,T_PHY,P_PHY,num_moist,                 &
  11      gd_cloud2_a,gd_cloud2_b,convtrans_avglen_m,         &
  12      adapt_step_flag,curr_secs,                    &
  13      ktau,dt,cu_phys,  &
  14      ids,ide, jds,jde, kds,kde,                                  &
  15      ims,ime, jms,jme, kms,kme,                                  &
  16      its,ite, jts,jte,kts,kte                                    )
  17     REAL, PARAMETER  ::  coef_p = 0.25, coef_gamm = 0.49, coef_alph = 100.
  18
  19   INTEGER,      INTENT(IN   ) :: ids,ide, jds,jde, kds,kde,      &
  20                                  ims,ime, jms,jme, kms,kme,      &
  21                                  its,ite, jts,jte, kts,kte,      &
  22                                  p_QV,p_QC,p_qi,num_moist
  23
  24   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                  &
  25        OPTIONAL,                                                 &
  26        INTENT(IN ) :: gd_cloud,gd_cloud2
  27   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                  &
  28        INTENT(IN ) :: t_phy,p_phy
  29   REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ),       &
  30        INTENT(IN ) :: moist
  31   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                  &
  32        OPTIONAL,                                                 &
  33        INTENT(INOUT ) :: gd_cloud_a,gd_cloud_b,gd_cloud2_a,      &
  34                          cldfr,gd_cloud2_b
  35   REAL, DIMENSION( ims:ime, jms:jme ),                           &
  36        INTENT(IN ) :: raincv
  37   REAL, DIMENSION( ims:ime, jms:jme ),                           &
  38        INTENT(INOUT ) :: raincv_a,raincv_b
  39   INTEGER, INTENT(IN) :: ktau,cu_phys
  40   INTEGER  :: stepave
  41   INTEGER, SAVE  :: stepave_count
  42   REAL, INTENT(IN) :: curr_secs
  43   REAL, INTENT(IN) :: convtrans_avglen_m, dt
  44   LOGICAL, INTENT(IN) :: adapt_step_flag
  45
  46   LOGICAL :: avg_end_flag, first_period_flag
  47   REAL :: satvp,rhgrid,h2oliq
  48   real :: pmax,pmin
  49 !
  50 ! Determine where we are in relation to the averaging period...
  51 !
  52 !  convtrans_avglen_m = 30.  !Averaging time for convective transport in min.
  53
  54   stepave=convtrans_avglen_m*60./dt
  55   avg_end_flag = .false.      !Initially assume we are not at the end
  56   first_period_flag = .false. !Nor at the beginning
  57   if( adapt_step_flag ) then
  58      !If end of period...
  59      if( curr_secs+real(dt,8)+0.01 >= &
  60           ( int( curr_secs/real(convtrans_avglen_m*60.,8) + 1_8, 8) &
  61             *real(convtrans_avglen_m*60.,8) ) ) &
  62           avg_end_flag = .true.
  63      if( curr_secs <= real(convtrans_avglen_m*60.,8) ) first_period_flag = .true.
  64   else
  65      if( mod(ktau,stepave)==0 ) avg_end_flag = .true.
  66      if( ktau <= stepave ) first_period_flag = .true.
  67   end if
  68 !
  69 ! Initialize the averaging arrays at the simulation start
  70 !
  71   if(ktau.le.1)then
  72      stepave_count             = 0
  73      raincv_a(its:ite,jts:jte) = 0.
  74      raincv_b(its:ite,jts:jte) = 0.
  75   end if
  76   if(present(gd_cloud2_a))then
  77      if(ktau.le.1) gd_cloud2_a(its:ite,kts:kte,jts:jte)=0.
  78   end if
  79   if(present(cldfr))then
  80      if(ktau.le.1) cldfr(its:ite,kts:kte,jts:jte)=0.
  81   end if
  82 !
  83 ! no time average available in first half hour
  84 !
  85   if( first_period_flag )then
  86      do j=jts,jte
  87         do i=its,ite
  88            raincv_b(i,j)=raincv(i,j)
  89         enddo
  90      enddo
  91   end if
  92 !
  93 ! build time average, and stored in raincv_b to be used by convective transport routine
  94 !
  95   stepave_count = stepave_count+1
  96   do j=jts,jte
  97      do i=its,ite
  98         raincv_a(i,j)=raincv_a(i,j)+raincv(i,j)
  99      enddo
 100   enddo
 101   if( avg_end_flag )then
 102      do j=jts,jte
 103         do i=its,ite
 104            raincv_b(i,j)=raincv_a(i,j)/real(stepave_count)
 105            raincv_a(i,j)=0.
 106         enddo
 107      enddo
 108   end if
 109 !
 110 ! do the same for convective parameterization cloud water mix ratio,
 111 ! currently only for cu_physics=3,5, used by both photolysis and atmospheric radiation
 112 !
 113   if(cu_phys.eq.3.or.cu_phys.eq.5.or.cu_phys.eq.93)then
 114 ! if(config_flags%cu_physics == GDSCHEME  .OR. &
 115 !    config_flags%cu_physics == GFSCHEME  .OR. &
 116 !    config_flags%cu_physics == GFSCHEME ) THEN
 117 !    pmax=maxval(gd_cloud)
 118 !    pmin=maxval(gd_cloud2)
 119 !    print *,pmax,pmin
 120      if( first_period_flag )then
 121         do j=jts,jte
 122            do k=kts,kte
 123               do i=its,ite
 124                  gd_cloud_b(i,k,j)=gd_cloud(i,k,j)
 125                  gd_cloud2_b(i,k,j)=gd_cloud2(i,k,j)
 126               enddo
 127            enddo
 128         enddo
 129      end if   ! stepave
 130 !
 131 !
 132 !
 133      do j=jts,jte
 134
 135         do k=kts,kte
 136            do i=its,ite
 137               gd_cloud_a(i,k,j)=gd_cloud_a(i,k,j)+gd_cloud(i,k,j)
 138               gd_cloud2_a(i,k,j)=gd_cloud2_a(i,k,j)+gd_cloud2(i,k,j)
 139            enddo
 140         enddo
 141      enddo
 142      if( avg_end_flag )then
 143         do j=jts,jte
 144            do k=kts,kte
 145               do i=its,ite
 146                  gd_cloud_b(i,k,j)=.1*gd_cloud_a(i,k,j)/real(stepave_count)
 147                  gd_cloud_a(i,k,j)=0.
 148                  gd_cloud2_b(i,k,j)=.1*gd_cloud2_a(i,k,j)/real(stepave_count)
 149                  gd_cloud2_a(i,k,j)=0.
 150               enddo
 151            enddo
 152         enddo
 153 !    pmax=maxval(gd_cloud_b)
 154 !    pmin=maxval(gd_cloud2_b)
 155 !    print *,'avg_end_flag ',pmax,pmin
 156      end if !stepave
 157   end if ! cu_rad_feedback
 158 !
 159 ! Clear the accumulator counter if we just finished an average...
 160 !
 161 if( avg_end_flag ) stepave_count = 0
 162 ! Siebesma et al., JAS, Vol. 60, no. 10, 1201-1219, 2003 (based on LES comparisons with trade-wind cumulus from BOMEX)
 163 ! SAM: Note units of liquid water and saturation vapor pressure must be in g/kg
 164     ! within the Siebesma et al. cloud fraction scheme
 165 if( first_period_flag .or. avg_end_flag )then
 166
 167         do j=jts,jte
 168            do k=kts,kte
 169               do i=its,ite
 170                 cldfr(i,k,j)=0.
 171 !               if(gd_cloud_b(i,k,j).gt.0  .or. gd_cloud2_b(i,k,j).gt.0)then
 172
 173                    if(p_qc.gt.1 .and. p_qi.le.1)then
 174
 175                       satvp = 3.80*exp(17.27*(t_phy(i,k,j)-273.)/ &
 176                             (t_phy(i,k,j)-36.))/(.01*p_phy(i,k,j))
 177                       rhgrid = max(.1,MIN( .95, moist(i,k,j,p_qv) /satvp))
 178                        h2oliq=1000.*(gd_cloud_b(i,k,j) + moist(i,k,j,p_qc))
 179                        satvp=1000.*satvp
 180                        cldfr(i,k,j)=(1.-exp(-coef_alph*h2oliq/((1.-rhgrid)*satvp)**coef_gamm))*(rhgrid**coef_p)
 181                        cldfr(i,k,j)=max(0.,MIN(1.,cldfr(i,k,j)))
 182                        if(moist(i,k,j,p_qc).eq.0)cldfr(i,k,j)=cldfr(i,k,j)*.2
 183                      endif
 184                    if(p_qc.gt.1 .and. p_qi.gt.1)then
 185                       satvp = 3.80*exp(17.27*(t_phy(i,k,j)-273.)/ &
 186                             (t_phy(i,k,j)-36.))/(.01*p_phy(i,k,j))
 187                       rhgrid = max(.1,MIN( .95, moist(i,k,j,p_qv) /satvp))
 188                        h2oliq=1000.*(gd_cloud_b(i,k,j) + moist(i,k,j,p_qc) + &
 189                                      gd_cloud2_b(i,k,j) + moist(i,k,j,p_qi))
 190                        satvp=1000.*satvp
 191                        cldfr(i,k,j)=(1.-exp(-coef_alph*h2oliq/((1.-rhgrid)*satvp)**coef_gamm))*(rhgrid**coef_p)
 192                        cldfr(i,k,j)=max(0.,MIN(1.,cldfr(i,k,j)))
 193                        if(moist(i,k,j,p_qc).eq.0 .and. moist(i,k,j,p_qi).eq.0)cldfr(i,k,j)=cldfr(i,k,j)*.2
 194                      endif
 195 !                  endif
 196               enddo
 197            enddo
 198         enddo
 199    endif
 200
 201 END subroutine convtrans_prep
 202
 203 END MODULE MODULE_CONVTRANS_prep
 204