source: lmdz_wrf/trunk/WRFV3/dyn_em/module_convtrans_prep.F @ 400

Last change on this file since 400 was 1, checked in by lfita, 11 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 7.5 KB
Line 
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
5Module module_convtrans_prep
6!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
7CONTAINS
8subroutine 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, stepave_count
41  REAL, INTENT(IN) :: curr_secs
42  REAL, INTENT(IN) :: convtrans_avglen_m, dt
43  LOGICAL, INTENT(IN) :: adapt_step_flag
44
45  LOGICAL :: avg_end_flag, first_period_flag
46  REAL :: satvp,rhgrid,h2oliq
47  real :: pmax,pmin
48!
49! Determine where we are in relation to the averaging period...
50!
51!  convtrans_avglen_m = 30.  !Averaging time for convective transport in min.
52
53  stepave=convtrans_avglen_m*60./dt
54  avg_end_flag = .false.      !Initially assume we are not at the end
55  first_period_flag = .false. !Nor at the beginning
56  if( adapt_step_flag ) then
57     !If end of period...
58     if( curr_secs+real(dt,8)+0.01 >= &
59          ( int( curr_secs/real(convtrans_avglen_m*60.,8) + 1_8, 8) &
60            *real(convtrans_avglen_m*60.,8) ) ) &
61          avg_end_flag = .true.
62     if( curr_secs <= real(convtrans_avglen_m*60.,8) ) first_period_flag = .true.
63  else
64     if( mod(ktau,stepave)==0 ) avg_end_flag = .true.
65     if( ktau <= stepave ) first_period_flag = .true.
66  end if
67!
68! Initialize the averaging arrays at the simulation start
69!
70  if(ktau.le.1)then
71     stepave_count             = 0
72     raincv_a(its:ite,jts:jte) = 0.
73     raincv_b(its:ite,jts:jte) = 0.
74  end if
75  if(present(gd_cloud2_a))then
76     if(ktau.le.1) gd_cloud2_a(its:ite,kts:kte,jts:jte)=0.
77  end if
78  if(present(cldfr))then
79     if(ktau.le.1) cldfr(its:ite,kts:kte,jts:jte)=0.
80  end if
81!
82! no time average available in first half hour
83!
84  if( first_period_flag )then
85     do j=jts,jte
86        do i=its,ite
87           raincv_b(i,j)=raincv(i,j)
88        enddo
89     enddo
90  end if
91!
92! build time average, and stored in raincv_b to be used by convective transport routine
93!
94  stepave_count = stepave_count+1
95  do j=jts,jte
96     do i=its,ite
97        raincv_a(i,j)=raincv_a(i,j)+raincv(i,j)
98     enddo
99  enddo
100  if( avg_end_flag )then
101     do j=jts,jte
102        do i=its,ite
103           raincv_b(i,j)=raincv_a(i,j)/real(stepave_count)
104           raincv_a(i,j)=0.
105        enddo
106     enddo
107  end if
108!
109! do the same for convective parameterization cloud water mix ratio,
110! currently only for cu_physics=3,5, used by both photolysis and atmospheric radiation
111!
112  if(cu_phys.eq.3.or.cu_phys.eq.5)then
113!    pmax=maxval(gd_cloud)
114!    pmin=maxval(gd_cloud2)
115!    print *,pmax,pmin
116     if( first_period_flag )then
117        do j=jts,jte
118           do k=kts,kte
119              do i=its,ite
120                 gd_cloud_b(i,k,j)=gd_cloud(i,k,j)
121                 gd_cloud2_b(i,k,j)=gd_cloud2(i,k,j)
122              enddo
123           enddo
124        enddo
125     end if   ! stepave
126!
127!
128!
129     do j=jts,jte
130
131        do k=kts,kte
132           do i=its,ite
133              gd_cloud_a(i,k,j)=gd_cloud_a(i,k,j)+gd_cloud(i,k,j)
134              gd_cloud2_a(i,k,j)=gd_cloud2_a(i,k,j)+gd_cloud2(i,k,j)
135           enddo
136        enddo
137     enddo
138     if( avg_end_flag )then
139        do j=jts,jte
140           do k=kts,kte
141              do i=its,ite
142                 gd_cloud_b(i,k,j)=.1*gd_cloud_a(i,k,j)/real(stepave_count)
143                 gd_cloud_a(i,k,j)=0.
144                 gd_cloud2_b(i,k,j)=.1*gd_cloud2_a(i,k,j)/real(stepave_count)
145                 gd_cloud2_a(i,k,j)=0.
146              enddo
147           enddo
148        enddo
149!    pmax=maxval(gd_cloud_b)
150!    pmin=maxval(gd_cloud2_b)
151!    print *,'avg_end_flag ',pmax,pmin
152     end if !stepave
153  end if ! cu_rad_feedback
154!
155! Clear the accumulator counter if we just finished an average...
156!
157if( avg_end_flag ) stepave_count = 0
158! Siebesma et al., JAS, Vol. 60, no. 10, 1201-1219, 2003 (based on LES comparisons with trade-wind cumulus from BOMEX)
159! SAM: Note units of liquid water and saturation vapor pressure must be in g/kg
160    ! within the Siebesma et al. cloud fraction scheme
161if( first_period_flag .or. avg_end_flag )then
162
163        do j=jts,jte
164           do k=kts,kte
165              do i=its,ite
166                cldfr(i,k,j)=0.
167!               if(gd_cloud_b(i,k,j).gt.0  .or. gd_cloud2_b(i,k,j).gt.0)then
168
169                   if(p_qc.gt.1 .and. p_qi.le.1)then
170
171                      satvp = 3.80*exp(17.27*(t_phy(i,k,j)-273.)/ &
172                            (t_phy(i,k,j)-36.))/(.01*p_phy(i,k,j))
173                      rhgrid = max(.1,MIN( .95, moist(i,k,j,p_qv) /satvp))
174                       h2oliq=1000.*(gd_cloud_b(i,k,j) + moist(i,k,j,p_qc))
175                       satvp=1000.*satvp
176                       cldfr(i,k,j)=(1.-exp(-coef_alph*h2oliq/((1.-rhgrid)*satvp)**coef_gamm))*(rhgrid**coef_p)
177                       cldfr(i,k,j)=max(0.,MIN(1.,cldfr(i,k,j)))
178                       if(moist(i,k,j,p_qc).eq.0)cldfr(i,k,j)=cldfr(i,k,j)*.2
179                     endif
180                   if(p_qc.gt.1 .and. p_qi.gt.1)then
181                      satvp = 3.80*exp(17.27*(t_phy(i,k,j)-273.)/ &
182                            (t_phy(i,k,j)-36.))/(.01*p_phy(i,k,j))
183                      rhgrid = max(.1,MIN( .95, moist(i,k,j,p_qv) /satvp))
184                       h2oliq=1000.*(gd_cloud_b(i,k,j) + moist(i,k,j,p_qc) + &
185                                     gd_cloud2_b(i,k,j) + moist(i,k,j,p_qi))
186                       satvp=1000.*satvp
187                       cldfr(i,k,j)=(1.-exp(-coef_alph*h2oliq/((1.-rhgrid)*satvp)**coef_gamm))*(rhgrid**coef_p)
188                       cldfr(i,k,j)=max(0.,MIN(1.,cldfr(i,k,j)))
189                       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
190                     endif
191!                  endif
192              enddo
193           enddo
194        enddo
195   endif
196
197END subroutine convtrans_prep
198
199END MODULE MODULE_CONVTRANS_prep
200
Note: See TracBrowser for help on using the repository browser.