GFS Operational Physics Documentation  Revision: 81451
gwdc.f
Go to the documentation of this file.
1 
4 
33 
72  subroutine gwdc(im,ix,iy,km,lat,u1,v1,t1,q1,
73  & pmid1,pint1,dpmid1,qmax,ktop,kbot,kcnv,cldf,
74  & grav,cp,rd,fv,dlength,lprnt,ipr,fhour,
75  & utgwc,vtgwc,tauctx,taucty)
76 
77 !***********************************************************************
78 ! ORIGINAL CODE FOR PARAMETERIZATION OF CONVECTIVELY FORCED
79 ! GRAVITY WAVE DRAG FROM YONSEI UNIVERSITY, KOREA
80 ! BASED ON THE THEORY GIVEN BY CHUN AND BAIK (JAS, 1998)
81 ! MODIFIED FOR IMPLEMENTATION INTO THE GFS/CFS BY
82 ! AKE JOHANSSON --- AUG 2005
83 !***********************************************************************
84 
85  USE machine , ONLY : kind_phys
86  implicit none
87 
88 !---------------------------- Arguments --------------------------------
89 !
90 ! Input variables
91 !
92 ! u : midpoint zonal wind
93 ! v : midpoint meridional wind
94 ! t : midpoint temperatures
95 ! pmid : midpoint pressures
96 ! pint : interface pressures
97 ! dpmid : midpoint delta p ( pi(k)-pi(k-1) )
98 ! lat : latitude index
99 ! qmax : deep convective heating
100 ! kcldtop : Vertical level index for cloud top ( mid level )
101 ! kcldbot : Vertical level index for cloud bottom ( mid level )
102 ! kcnv : (0,1) dependent on whether convection occur or not
103 !
104 ! Output variables
105 !
106 ! utgwc : zonal wind tendency
107 ! vtgwc : meridional wind tendency
108 !
109 !-----------------------------------------------------------------------
110 
111  integer im, ix, iy, km, lat, ipr
112  integer ktop(im),kbot(im),kcnv(im)
113 
114  real(kind=kind_phys) grav,cp,rd,fv,fhour,fhourpr
115  real(kind=kind_phys), dimension(ix) :: qmax
116  &, tauctx, taucty
117  real(kind=kind_phys), dimension(im) :: cldf,dlength
118  real(kind=kind_phys), dimension(ix,km) :: u1,v1,t1,q1,
119  & pmid1,dpmid1
120 ! &, cumchr1
121  real(kind=kind_phys), dimension(iy,km) :: utgwc,vtgwc
122  real(kind=kind_phys), dimension(ix,km+1) :: pint1
123 
124  logical lprnt
125 
126 !------------------------- Local workspace -----------------------------
127 !
128 ! i, k : Loop index
129 ! kk : Loop index
130 ! cldf : Deep convective cloud fraction at the cloud top.
131 ! ugwdc : Zonal wind after GWDC paramterization
132 ! vgwdc : Meridional wind after GWDC parameterization
133 ! plnmid : Log(pmid) ( mid level )
134 ! plnint : Log(pint) ( interface level )
135 ! dpint : Delta pmid ( interface level )
136 ! tauct : Wave stress at the cloud top calculated using basic-wind
137 ! parallel to the wind vector at the cloud top ( mid level )
138 ! tauctx : Wave stress at the cloud top projected in the east
139 ! taucty : Wave stress at the cloud top projected in the north
140 ! qmax : Maximum deep convective heating rate ( K s-1 ) in a
141 ! horizontal grid point calculated from cumulus para-
142 ! meterization. ( mid level )
143 ! wtgwc : Wind tendency in direction to the wind vector at the cloud top level
144 ! due to convectively generated gravity waves ( mid level )
145 ! utgwcl : Zonal wind tendency due to convectively generated
146 ! gravity waves ( mid level )
147 ! vtgwcl : Meridional wind tendency due to convectively generated
148 ! gravity waves ( mid level )
149 ! taugwci : Profile of wave stress calculated using basic-wind
150 ! parallel to the wind vector at the cloud top
151 ! taugwcxi : Profile of zonal component of gravity wave stress
152 ! taugwcyi : Profile of meridional component of gravity wave stress
153 !
154 ! taugwci, taugwcxi, and taugwcyi are defined at the interface level
155 !
156 ! bruni : Brunt-Vaisala frequency ( interface level )
157 ! brunm : Brunt-Vaisala frequency ( mid level )
158 ! rhoi : Air density ( interface level )
159 ! rhom : Air density ( mid level )
160 ! ti : Temperature ( interface level )
161 ! basicum : Basic-wind profile. Basic-wind is parallel to the wind
162 ! vector at the cloud top level. (mid level)
163 ! basicui : Basic-wind profile. Basic-wind is parallel to the wind
164 ! vector at the cloud top level. ( interface level )
165 ! riloc : Local Richardson number ( interface level )
166 ! rimin : Minimum Richardson number including both the basic-state
167 ! and gravity wave effects ( interface level )
168 ! gwdcloc : Horizontal location where the GWDC scheme is activated.
169 ! break : Horizontal location where wave breaking is occurred.
170 ! critic : Horizontal location where critical level filtering is
171 ! occurred.
172 ! dogwdc : Logical flag whether the GWDC parameterization is
173 ! calculated at a grid point or not.
174 !
175 ! dogwdc is used in order to lessen CPU time for GWDC calculation.
176 !
177 !-----------------------------------------------------------------------
178 
179  integer i,ii,k,k1,kk,kb,ilev,npt,kcb,kcldm,npr
180  integer, dimension(im) :: ipt
181 
182  real(kind=kind_phys) tem, tem1, tem2, qtem, wtgwc, tauct,
183  & windcltop, shear, nonlinct, nonlin, nonlins,
184  & n2, dtdp, crit1, crit2, pi, p1, p2,
185  & gsqr, onebg
186 ! & taus, n2, dtdp, crit1, crit2, pi, p1, p2
187 
188  integer, allocatable :: kcldtop(:),kcldbot(:)
189  logical, allocatable :: do_gwc(:)
190  real(kind=kind_phys), allocatable :: tauctxl(:), tauctyl(:),
191  & gwdcloc(:), break(:),
192  & critic(:),
193 ! & critic(:), angle(:),
194  & cosphi(:), sinphi(:),
195  & xstress(:), ystress(:),
196  & ucltop(:), vcltop(:),
197  & wrk(:),
198  & dlen(:), gqmcldlen(:)
199 ! real(kind=kind_phys), allocatable :: plnint(:,:), dpint(:,:),
200 ! & taugwci(:,:), taugwcxi(:,:),
201 ! & taugwcyi(:,:), bruni(:,:),
202 ! & taugwcyi(:,:), bruni(:,:),
203  real(kind=kind_phys), allocatable :: plnint(:,:),
204  & taugwci(:,:), bruni(:,:),
205  & rhoi(:,:), basicui(:,:),
206  & ti(:,:), riloc(:,:),
207  & rimin(:,:), pint(:,:)
208 ! real(kind=kind_phys), allocatable :: ugwdc(:,:), vgwdc(:,:),
209  real(kind=kind_phys), allocatable ::
210 ! & plnmid(:,:), wtgwc(:,:),
211  & plnmid(:,:),
212  & utgwcl(:,:), vtgwcl(:,:),
213  & basicum(:,:), u(:,:),v(:,:),
214  & t(:,:), spfh(:,:),
215  & pmid(:,:), dpmid(:,:),
216 ! & pmid(:,:), cumchr(:,:),
217  & brunm(:,:), rhom(:,:)
218 
219 !-----------------------------------------------------------------------
220 !
221 ! ucltop : Zonal wind at the cloud top ( mid level )
222 ! vcltop : Meridional wind at the cloud top ( mid level )
223 ! windcltop : Wind speed at the cloud top ( mid level )
224 ! shear : Vertical shear of basic wind
225 ! cosphi : Cosine of angle of wind vector at the cloud top
226 ! sinphi : Sine of angle of wind vector at the cloud top
227 ! c1 : Tunable parameter
228 ! c2 : Tunable parameter
229 ! dlength : Grid spacing in the direction of basic wind at the cloud top
230 ! nonlinct : Nonlinear parameter at the cloud top
231 ! nonlin : Nonlinear parameter above the cloud top
232 ! nonlins : Saturation nonlinear parameter
233 ! taus : Saturation gravity wave drag == taugwci(i,k)
234 ! n2 : Square of Brunt-Vaisala frequency
235 ! dtdp : dT/dp
236 ! xstress : Vertically integrated zonal momentum change due to GWDC
237 ! ystress : Vertically integrated meridional momentum change due to GWDC
238 ! crit1 : Variable 1 for checking critical level
239 ! crit2 : Variable 2 for checking critical level
240 !
241 !-----------------------------------------------------------------------
242 
243  real(kind=kind_phys), parameter ::
244  & c1=1.41, c2=-0.38, ricrit=0.25
245  &, n2min=1.e-32, zero=0.0, one=1.0
246 ! &, taumin=1.0e-20, tauctmax=-20.
247  &, taumin=1.0e-20, tauctmax=-5.
248  &, qmin=1.0e-10, shmin=1.0e-20
249  &, rimax=1.0e+20, rimaxm=0.99e+20
250  &, rimaxp=1.01e+20, rilarge=0.9e+20
251  &, riminx=-1.0e+20, riminm=-1.01e+20
252  &, riminp=-0.99e+20, rismall=-0.9e+20
253 
254 !
255  npt = 0
256  do i = 1,im
257  ipt(i) = 0
258  if (kcnv(i) /= 0 .and. qmax(i) > zero) then
259  npt = npt + 1
260  ipt(npt) = i
261  endif
262  enddo
263  do k=1,km
264  do i=1,im
265  utgwc(i,k) = 0.0
266  vtgwc(i,k) = 0.0
267 ! brunm(i,k) = 0.0
268 ! rhom(i,k) = 0.0
269  enddo
270  enddo
271  do i=1,im
272  tauctx(i) = 0.0
273  taucty(i) = 0.0
274  enddo
275  if (npt == 0) return ! No gwdc calculation done!
276 
277 !***********************************************************************
278 !
279 ! Begin GWDC
280 !
281 !***********************************************************************
282 
283 !-----------------------------------------------------------------------
284 ! Write out incoming variables
285 !-----------------------------------------------------------------------
286 
287  fhourpr = zero
288 ! if (lprnt) then
289 ! if (fhour >= fhourpr) then
290 ! print *,' '
291 ! write(*,*) 'Inside GWDC raw input start print at fhour = ',
292 ! & fhour
293 ! write(*,*) 'IX IM KM ',ix,im,km
294 ! write(*,*) 'KBOT KTOP QMAX DLENGTH kcnv ',
295 ! + kbot(ipr),ktop(ipr),qmax(ipr),dlength(ipr),kcnv(ipr)
296 ! write(*,*) 'grav cp rd ',grav,cp,rd
297 
298 !-------- Pressure levels ----------
299 ! write(*,9100)
300 ! ilev=km+1
301 ! write(*,9110) ilev,(10.*pint1(ipr,ilev))
302 ! do ilev=km,1,-1
303 ! write(*,9120) ilev,(10.*pmid1(ipr,ilev)),
304 ! & (10.*dpmid1(ipr,ilev))
305 ! write(*,9110) ilev,(10.*pint1(ipr,ilev))
306 ! enddo
307 
308 !-------- U1 V1 T1 ----------
309 ! write(*,9130)
310 ! do ilev=km,1,-1
311 ! write(*,9140) ilev,U1(ipr,ilev),V1(ipr,ilev),T1(ipr,ilev)
312 ! enddo
313 
314 ! print *,' '
315 ! print *,' Inside GWDC raw input end print'
316 ! endif
317 ! endif
318 
319 !9100 format(//,14x,'PRESSURE LEVELS',//,
320 ! +' ILEV',6x,'PINT1',7x,'PMID1',6x,'DPMID1',/)
321 !9110 format(i4,2x,f10.3)
322 !9120 format(i4,12x,2(2x,f10.3))
323 !9130 format(//,' ILEV',7x,'U1',10x,'V1',10x,'T1',/)
324 !9140 format(i4,3(2x,f10.3))
325 
326 ! Allocate local arrays
327 
328  allocate (kcldtop(npt), kcldbot(npt), do_gwc(npt))
329  allocate (tauctxl(npt), tauctyl(npt),
330  & gwdcloc(npt), break(npt), critic(npt), cosphi(npt),
331  & sinphi(npt), xstress(npt), ystress(npt), wrk(npt),
332  & ucltop(npt), vcltop(npt),dlen(npt), gqmcldlen(npt))
333 
334 ! allocate (plnint(npt,2:km+1), dpint(npt,km+1),
335 ! & taugwci(npt,km+1), taugwcxi(npt,km+1),
336 ! & taugwcyi(npt,km+1), bruni(npt,km+1),
337  allocate (plnint(npt,2:km+1),
338  & taugwci(npt,km+1), bruni(npt,km+1),
339  & rhoi(npt,km+1), basicui(npt,km+1),
340  & ti(npt,km+1), riloc(npt,km+1),
341  & rimin(npt,km+1), pint(npt,km+1))
342 
343 ! allocate (ugwdc(npt,km), vgwdc(npt,km),
344  allocate
345 ! & (plnmid(npt,km), wtgwc(npt,km),
346  & (plnmid(npt,km),
347  & utgwcl(npt,km), vtgwcl(npt,km),
348  & basicum(npt,km), u(npt,km), v(npt,km),
349  & t(npt,km), spfh(npt,km), pmid(npt,km),
350  & dpmid(npt,km),
351 ! & dpmid(npt,km), cumchr(npt,km),
352  & brunm(npt,km), rhom(npt,km))
353 
354 !-----------------------------------------------------------------------
357 !-----------------------------------------------------------------------
358  gsqr = grav * grav
359  onebg = one / grav
360 
361  if (lprnt) then
362  npr = 1
363  do i=1,npt
364  if (ipr == ipt(i))then
365  npr = i
366  exit
367  endif
368  enddo
369  endif
370 
371  do k=1,km
372  k1 = km - k + 1
373  do i=1,npt
374  ii = ipt(i)
375  u(i,k) = u1(ii,k1)
376  v(i,k) = v1(ii,k1)
377  t(i,k) = t1(ii,k1)
378  spfh(i,k) = max(q1(ii,k1),qmin)
379  pmid(i,k) = pmid1(ii,k1)
380  dpmid(i,k) = dpmid1(ii,k1) * onebg
381 ! cumchr(i,k) = cumchr1(ii,k1)
382 
383  rhom(i,k) = pmid(i,k) / (rd*t(i,k)*(1.0+fv*spfh(i,k)))
384  plnmid(i,k) = log(pmid(i,k))
385  utgwcl(i,k) = zero
386  vtgwcl(i,k) = zero
387 ! ugwdc(i,k) = zero
388 ! vgwdc(i,k) = zero
389  brunm(i,k) = zero
390  basicum(i,k) = zero
391  enddo
392  enddo
393 
394  do k=1,km+1
395  k1 = km - k + 2
396  do i=1,npt
397  ii = ipt(i)
398  pint(i,k) = pint1(ii,k1)
399  taugwci(i,k) = zero
400  bruni(i,k) = zero
401  rhoi(i,k) = zero
402  ti(i,k) = zero
403  basicui(i,k) = zero
404  riloc(i,k) = zero
405  rimin(i,k) = zero
406  enddo
407  enddo
408  do k=2,km+1
409  do i=1,npt
410  plnint(i,k) = log(pint(i,k))
411  enddo
412  enddo
413 
414  do i = 1, npt
415  ii = ipt(i)
416  kcldtop(i) = km - ktop(ii) + 1
417  kcldbot(i) = km - kbot(ii) + 1
418  dlen(i) = dlength(ii)
419 ! (g*qmax(ii)*cldf(ii)*dlength(ii))
420  gqmcldlen(i) = grav*qmax(ii)*cldf(ii)*dlen(i)
421  enddo
422 
423 ! if (lprnt) then
424 ! if (fhour.ge.fhourpr) then
425 ! write(*,9200)
426 ! do i=1,im
427 ! write(*,9201) kcnv(i),kcldbot(i),kcldtop(i)
428 ! enddo
429 ! endif
430 ! endif
431 
432 !9200 format(//,' Inside GWDC local variables start print',//,
433 ! +2x,'kcnv',2x,'KCLDBOT',2x,'KCLDTOP',//)
434 !9201 format(i4,2x,i5,4x,i5)
435 
436 !***********************************************************************
437 
438  pi = 2.*asin(1.)
439 
440 !-----------------------------------------------------------------------
441 !
442 ! PRESSURE VARIABLES
443 !
444 ! Interface 1 ======== pint(1) *********
445 ! Mid-Level 1 -------- pmid(1) dpmid(1)
446 ! 2 ======== pint(2) dpint(2)
447 ! 2 -------- pmid(2) dpmid(2)
448 ! 3 ======== pint(3) dpint(3)
449 ! 3 -------- pmid(3) dpmid(3)
450 ! 4 ======== pint(4) dpint(4)
451 ! 4 -------- pmid(4) dpmid(4)
452 ! ........
453 ! 17 ======== pint(17) dpint(17)
454 ! 17 -------- pmid(17) dpmid(17)
455 ! 18 ======== pint(18) dpint(18)
456 ! 18 -------- pmid(18) dpmid(18)
457 ! 19 ======== pint(19) *********
458 !
459 !-----------------------------------------------------------------------
460 
461  do i = 1, npt
462  tauctxl(i) = zero
463  tauctyl(i) = zero
464 
465 !-----------------------------------------------------------------------
466 ! THERMAL VARIABLES
467 !
468 ! Interface 1 ======== TI(1) RHOI(1) BRUNI(1)
469 ! 1 -------- T(1) RHOM(1) BRUNM(1)
470 ! 2 ======== TI(2) RHOI(2) BRUNI(2)
471 ! 2 -------- T(2) RHOM(2) BRUNM(2)
472 ! 3 ======== TI(3) RHOI(3) BRUNI(3)
473 ! 3 -------- T(3) RHOM(3) BRUNM(3)
474 ! 4 ======== TI(4) RHOI(4) BRUNI(4)
475 ! 4 -------- T(4) RHOM(4) BRUNM(4)
476 ! ........
477 ! 17 ========
478 ! 17 -------- T(17) RHOM(17) BRUNM(17)
479 ! 18 ======== TI(18) RHOI(18) BRUNI(18)
480 ! 18 -------- T(18) RHOM(18) BRUNM(18)
481 ! 19 ======== TI(19) RHOI(19) BRUNI(19)
482 !
483 
484 !
488 
489  ti(i,1) = t(i,1)
490  rhoi(i,1) = pint(i,1)/(rd*ti(i,1))
491  bruni(i,1) = sqrt( gsqr / (cp*ti(i,1)) )
492 !
496 
497  ti(i,km+1) = t(i,km)
498  rhoi(i,km+1) = pint(i,km+1)/(rd*ti(i,km+1)*(1.0+fv*spfh(i,km)))
499  bruni(i,km+1) = sqrt( gsqr / (cp*ti(i,km+1)) )
500  enddo
501 
502 !-----------------------------------------------------------------------
503 !
507 !
508 !-----------------------------------------------------------------------
509 
510  do k = 2, km
511  do i = 1, npt
512  tem1 = (plnmid(i,k)-plnint(i,k)) / (plnmid(i,k)-plnmid(i,k-1))
513  tem2 = one - tem1
514  ti(i,k) = t(i,k-1) * tem1 + t(i,k) * tem2
515  qtem = spfh(i,k-1) * tem1 + spfh(i,k) * tem2
516  rhoi(i,k) = pint(i,k) / ( rd * ti(i,k)*(1.0+fv*qtem) )
517  dtdp = (t(i,k)-t(i,k-1)) / (pmid(i,k)-pmid(i,k-1))
518  n2 = gsqr / ti(i,k) * ( 1./cp - rhoi(i,k)*dtdp )
519  bruni(i,k) = sqrt(max(n2min, n2))
520  enddo
521  enddo
522 
523  deallocate (spfh)
524 !-----------------------------------------------------------------------
525 !
528 !-----------------------------------------------------------------------
529 
530  do k = 1, km
531  do i = 1, npt
532  dtdp = (ti(i,k+1)-ti(i,k)) / (pint(i,k+1)-pint(i,k))
533  n2 = gsqr / t(i,k) * ( 1./cp - rhom(i,k)*dtdp )
534  brunm(i,k) = sqrt(max(n2min, n2))
535  enddo
536  enddo
537 
538 !-----------------------------------------------------------------------
539 ! PRINTOUT
540 !-----------------------------------------------------------------------
541 
542 ! if (lprnt) then
543 ! if (fhour.ge.fhourpr) then
544 
545 !-------- Pressure levels ----------
546 ! write(*,9101)
547 ! do ilev=1,km
548 ! write(*,9111) ilev,(0.01*pint(ipr,ilev)),
549 ! & (0.01*dpint(ipr,ilev)),plnint(ipr,ilev)
550 ! write(*,9121) ilev,(0.01*pmid(ipr,ilev)),
551 ! & (0.01*dpmid(ipr,ilev)),plnmid(ipr,ilev)
552 ! enddo
553 ! ilev=km+1
554 ! write(*,9111) ilev,(0.01*pint(ipr,ilev)),
555 ! & (0.01*dpint(ipr,ilev)),plnint(ipr,ilev)
556 
557 ! 2
558 !-------- U V T N ----------
559 ! write(*,9102)
560 ! do ilev=1,km
561 ! write(*,9112) ilev,ti(ipr,ilev),(100.*bruni(ipr,ilev))
562 ! write(*,9122) ilev,u(ipr,ilev),v(ipr,ilev),
563 ! + t(ipr,ilev),(100.*brunm(ipr,ilev))
564 ! enddo
565 ! ilev=km+1
566 ! write(*,9112) ilev,ti(ipr,ilev),(100.*bruni(ipr,ilev))
567 
568 ! endif
569 ! endif
570 
571 !9101 format(//,14x,'PRESSURE LEVELS',//,
572 ! +' ILEV',4x,'PINT',4x,'PMID',4x,'DPINT',3x,'DPMID',5x,'LNP',/)
573 !9111 format(i4,1x,f8.2,9x,f8.2,9x,f8.2)
574 !9121 format(i4,9x,f8.2,9x,f8.2,1x,f8.2)
575 !9102 format(//' ILEV',5x,'U',7x,'V',5x,'TI',7x,'T',
576 ! +5x,'BRUNI',3x,'BRUNM',//)
577 !9112 format(i4,16x,f8.2,8x,f8.3)
578 !9122 format(i4,2f8.2,8x,f8.2,8x,f8.3)
579 
580 
581 !***********************************************************************
582 !
583 ! Big loop over grid points ONLY done if kcnv=1
584 !
585 !***********************************************************************
586 
587  kcldm = 1
588  do i = 1, npt
589  kk = kcldtop(i)
590  kb = kcldbot(i)
591  kcldm = max(kcldm,kk)
592 
593 !-----------------------------------------------------------------------
594 !
598 !
599 !-----------------------------------------------------------------------
600 
601  ucltop(i) = u(i,kk)
602  vcltop(i) = v(i,kk)
603 ! windcltop = sqrt( ucltop(i)*ucltop(i) + vcltop(i)*vcltop(i) )
604  windcltop = 1.0 / sqrt( ucltop(i)*ucltop(i)
605  & + vcltop(i)*vcltop(i) )
606  cosphi(i) = ucltop(i)*windcltop
607  sinphi(i) = vcltop(i)*windcltop
608 ! angle(i) = acos(cosphi)*180./pi
609  enddo
610 
611 !-----------------------------------------------------------------------
612 !
619 ! Input u(i,k) and v(i,k) is defined at mid level
620 !
621 !-----------------------------------------------------------------------
622 
623  do k=1,km
624  do i=1,npt
625  basicum(i,k) = u(i,k)*cosphi(i) + v(i,k)*sinphi(i)
626  enddo
627  enddo
628 
629 !-----------------------------------------------------------------------
630 !
631 ! Basic state wind at interface level is also calculated
632 ! based on linear interpolation in ln(Pressure)
633 !
634 ! In the top and bottom boundaries, basic-state wind at interface level
635 ! is assumed to be vertically uniform.
636 !
637 !-----------------------------------------------------------------------
638 
639  do i=1,npt
640  basicui(i,1) = basicum(i,1)
641  basicui(i,km+1) = basicum(i,km)
642  enddo
643  do k=2,km
644  do i=1,npt
645  tem1 = (plnmid(i,k)-plnint(i,k)) / (plnmid(i,k)-plnmid(i,k-1))
646  tem2 = one - tem1
647  basicui(i,k) = basicum(i,k)*tem2 + basicum(i,k-1)*tem1
648  enddo
649  enddo
650 
651 !-----------------------------------------------------------------------
652 !
658 
659 ! basicum : U at mid level
660 ! basicui : UI at interface level
661 !
662 ! Interface 1 ======== UI(1) rhoi(1) bruni(1) riloc(1)
663 ! Mid-level 1 -------- U(1)
664 ! 2 ======== UI(2) dpint(2) rhoi(2) bruni(2) riloc(2)
665 ! 2 -------- U(2)
666 ! 3 ======== UI(3) dpint(3) rhoi(3) bruni(3) riloc(3)
667 ! 3 -------- U(3)
668 ! 4 ======== UI(4) dpint(4) rhoi(4) bruni(4) riloc(4)
669 ! 4 -------- U(4)
670 ! ........
671 ! 17 ======== UI(17) dpint(17) rhoi(17) bruni(17) riloc(17)
672 ! 17 -------- U(17)
673 ! 18 ======== UI(18) dpint(18) rhoi(18) bruni(18) riloc(18)
674 ! 18 -------- U(18)
675 ! 19 ======== UI(19) rhoi(19) bruni(19) riloc(19)
676 !
677 !-----------------------------------------------------------------------
678 
679  do k=2,km
680  do i=1,npt
681  shear = grav*rhoi(i,k) * (basicum(i,k) - basicum(i,k-1))
682  & / (pmid(i,k) - pmid(i,k-1))
683  if ( abs(shear) < shmin ) then
684  riloc(i,k) = rimax
685  else
686  tem = bruni(i,k) / shear
687  riloc(i,k) = tem * tem
688  if (riloc(i,k) >= rimax ) riloc(i,k) = rilarge
689  end if
690  enddo
691  enddo
692 
693  do i=1,npt
694  riloc(i,1) = riloc(i,2)
695  riloc(i,km+1) = riloc(i,km)
696  enddo
697 
698 ! if (lprnt.and.(i.eq.ipr)) then
699 ! if (fhour.ge.fhourpr) then
700 ! write(*,9104) ucltop,vcltop,windcltop,angle,kk
701 ! do ilev=1,km
702 ! write(*,9114) ilev,basicui(ipr,ilev),dpint(ipr,ilev),
703 ! + rhoi(ipr,ilev),(100.*bruni(ipr,ilev)),riloc(ilev)
704 ! write(*,9124) ilev,(basicum(ipr,ilev))
705 ! enddo
706 ! ilev=km+1
707 ! write(*,9114) ilev,basicui(ipr,ilev),dpint(ipr,ilev),
708 ! + rhoi(ipr,ilev),(100.*bruni(ipr,ilev)),riloc(ilev)
709 ! endif
710 ! endif
711 
712 !9104 format(//,'WIND VECTOR AT CLOUDTOP = (',f6.2,' , ',f6.2,' ) = ',
713 ! +f6.2,' IN DIRECTION ',f6.2,4x,'KK = ',i2,//,
714 ! +' ILEV',2x,'BASICUM',2x,'BASICUI',4x,'DPINT',6x,'RHOI',5x,
715 ! +'BRUNI',6x,'RI',/)
716 !9114 format(i4,10x,f8.2,4(2x,f8.2))
717 !9124 format(i4,1x,f8.2)
718 
719 !-----------------------------------------------------------------------
720 !
722 !
723 ! kcldtopi : The interface level cloud top index
724 ! kcldtop : The midlevel cloud top index
725 ! kcldbot : The midlevel cloud bottom index
726 !
727 ! A : Find deep convective heating rate maximum
728 !
729 ! If kcldtop(i) is less than kcldbot(i) in a horizontal grid point,
730 ! it can be thought that there is deep convective cloud. However,
731 ! deep convective heating between kcldbot and kcldtop is sometimes
732 ! zero in spite of kcldtop less than kcldbot. In this case,
733 ! maximum deep convective heating is assumed to be 1.e-30.
734 !
735 ! B : kk is the vertical index for interface level cloud top
736 !
737 ! C : Total convective fractional cover (cldf) is used as the
738 ! convective cloud cover for GWDC calculation instead of
739 ! convective cloud cover in each layer (concld).
740 ! a1 = cldf*dlength
741 ! You can see the difference between cldf(i) and concld(i)
742 ! in (4.a.2) in Description of the NCAR Community Climate
743 ! Model (CCM3).
744 ! In NCAR CCM3, cloud fractional cover in each layer in a deep
745 ! cumulus convection is determined assuming total convective
746 ! cloud cover is randomly overlapped in each layer in the
747 ! cumulus convection.
748 !
749 ! D : Wave stress at cloud top is calculated when the atmosphere
750 ! is dynamically stable at the cloud top
751 !
752 ! E : Cloud top wave stress and nonlinear parameter are calculated
753 ! using density, temperature, and wind that are defined at mid
754 ! level just below the interface level in which cloud top wave
755 ! stress is defined.
756 ! Nonlinct is defined at the interface level.
757 !
758 ! F : If the atmosphere is dynamically unstable at the cloud top,
759 ! GWDC calculation in current horizontal grid is skipped.
760 !
761 ! G : If mean wind at the cloud top is less than zero, GWDC
762 
798 ! - If mean wind at the cloud top is less than zero, GWDC
799 ! calculation in current horizontal grid is skipped.
800 !
801 
804 !
805 !-----------------------------------------------------------------------
806 !D
807  do i=1,npt
808  kk = kcldtop(i)
809  if ( abs(basicui(i,kk)) > zero .and. riloc(i,kk) > ricrit) then
810 !E
811  tem = basicum(i,kk)
812  tem1 = tem * tem
813  nonlinct = gqmcldlen(i) / (bruni(i,kk)*t(i,kk)*tem1) ! Mu
814  tem2 = c2*nonlinct
815 ! RhoU^3c1(c2mu)^2/Ndx
816  tauct = - rhom(i,kk) * tem * tem1 * c1 * tem2 * tem2
817  & / (bruni(i,kk)*dlen(i))
818 
819  tauct = max(tauctmax, tauct)
820  tauctxl(i) = tauct * cosphi(i) ! X stress at cloud top
821  tauctyl(i) = tauct * sinphi(i) ! Y stress at cloud top
822  taugwci(i,kk) = tauct ! *1
823  do_gwc(i) = .true.
824  else
825 !F
826  tauctxl(i) = zero
827  tauctyl(i) = zero
828  do_gwc(i) = .false.
829  end if
830 !H
831  enddo
832 
833 ! if (lprnt.and.(i.eq.ipr)) then
834 ! if (fhour.ge.fhourpr) then
835 ! write(*,9210) tauctx(ipr),taucty(ipr),tauct(ipr),angle,kk
836 ! endif
837 ! endif
838 
839 !9210 format(/,5x,'STRESS VECTOR = ( ',f8.3,' , ',f8.3,' ) = ',f8.3,
840 ! +' IN DIRECTION ',f6.2,4x,'KK = ',i2,/)
841 
842 !-----------------------------------------------------------------------
843 !
844 ! At this point, mean wind at the cloud top is larger than zero and
845 ! local RI at the cloud top is larger than ricrit (=0.25)
846 !
847 ! Calculate minimum of Richardson number including both basic-state
848 ! condition and wave effects.
849 !
850 ! g*Q_0*alpha*dx RI_loc*(1 - mu*|c2|)
851 ! mu = ---------------- RI_min = -----------------------------
852 ! c_p*N*T*U^2 (1 + mu*RI_loc^(0.5)*|c2|)^2
853 !
854 ! Minimum RI is calculated for the following two cases
855 !
856 ! (1) RIloc < 1.e+20
857 ! (2) Riloc = 1.e+20 ----> Vertically uniform basic-state wind
858 !
859 ! RIloc cannot be smaller than zero because N^2 becomes 1.E-32 in the
860 ! case of N^2 < 0.. Thus the sign of RINUM is determined by
861 ! 1 - nonlin*|c2|.
862 !
863 !-----------------------------------------------------------------------
869 
870  do k=kcldm,1,-1
871 
872  do i=1,npt
873  if (do_gwc(i)) then
874  kk = kcldtop(i)
875  if (k > kk) cycle
876  if ( k /= 1 ) then
877  crit1 = ucltop(i)*(u(i,k)+u(i,k-1))*0.5
878  crit2 = vcltop(i)*(v(i,k)+v(i,k-1))*0.5
879  else
880  crit1 = ucltop(i)*u(i,1)
881  crit2 = vcltop(i)*v(i,1)
882  end if
883 
884  if ( abs(basicui(i,k)) > zero .and. crit1 > zero
885  & .and. crit2 > zero ) then
886  tem = basicui(i,k) * basicui(i,k)
887  nonlin = gqmcldlen(i) / (bruni(i,k)*ti(i,k)*tem)
888  tem = nonlin*abs(c2)
889  if ( riloc(i,k) < rimaxm ) then
890  tem1 = 1 + tem*sqrt(riloc(i,k))
891  rimin(i,k) = riloc(i,k) * (1-tem) / (tem1*tem1)
892  else if((riloc(i,k) > rimaxm) .and.
893  & (riloc(i,k) < rimaxp)) then
894  rimin(i,k) = ( 1 - tem) / (tem*tem)
895  end if
896  if ( rimin(i,k) <= riminx ) then
897  rimin(i,k) = rismall
898  end if
899  else
900  rimin(i,k) = riminx
901  end if
902 
903 !-----------------------------------------------------------------------
904 !
905 ! If the minimum \f$R_{i}\f$ at interface cloud top is less than or equal to 1/4,
906 ! the convective GWD calculation is skipped at that grid point.
907 !
908 !-----------------------------------------------------------------------
909 
910 !-----------------------------------------------------------------------
911 !
914 !
915 ! Assuming kcldtop(i)=10 and kcldbot=16,
916 !
917 ! TAUGWCI RIloc RImin UTGWC
918 !
919 ! Interface 1 ======== - 0.001 -1.e20
920 ! 1 -------- 0.000
921 ! 2 ======== - 0.001 -1.e20
922 ! 2 -------- 0.000
923 ! 3 ======== - 0.001 -1.e20
924 ! 3 -------- -.xxx
925 ! 4 ======== - 0.001 2.600 2.000
926 ! 4 -------- 0.000
927 ! 5 ======== - 0.001 2.500 2.000
928 ! 5 -------- 0.000
929 ! 6 ======== - 0.001 1.500 0.110
930 ! 6 -------- +.xxx
931 ! 7 ======== - 0.005 2.000 3.000
932 ! 7 -------- 0.000
933 ! 8 ======== - 0.005 1.000 0.222
934 ! 8 -------- +.xxx
935 ! 9 ======== - 0.010 1.000 2.000
936 ! 9 -------- 0.000
937 ! kcldtopi 10 ======== $$$ - 0.010
938 ! kcldtop 10 -------- $$$ yyyyy
939 ! 11 ======== $$$ 0
940 ! 11 -------- $$$
941 ! 12 ======== $$$ 0
942 ! 12 -------- $$$
943 ! 13 ======== $$$ 0
944 ! 13 -------- $$$
945 ! 14 ======== $$$ 0
946 ! 14 -------- $$$
947 ! 15 ======== $$$ 0
948 ! 15 -------- $$$
949 ! 16 ======== $$$ 0
950 ! kcldbot 16 -------- $$$
951 ! 17 ======== 0
952 ! 17 --------
953 ! 18 ======== 0
954 ! 18 --------
955 ! 19 ======== 0
956 !
957 !-----------------------------------------------------------------------
958 !
959 ! Even though the cloud top level obtained in deep convective para-
960 ! meterization is defined in mid-level, the cloud top level for
961 ! the GWDC calculation is assumed to be the interface level just
962 ! above the mid-level cloud top vertical level index.
963 !
964 !-----------------------------------------------------------------------
965 
976 
977  if (k < kk .and. k > 1) then
978  if ( abs(taugwci(i,k+1)) > taumin ) then ! TAUGWCI
979  if ( riloc(i,k) > ricrit ) then ! RIloc
980  if ( rimin(i,k) > ricrit ) then ! RImin
981  taugwci(i,k) = taugwci(i,k+1)
982  elseif (rimin(i,k) > riminp) then
983  tem = 2.0 + 1.0 / sqrt(riloc(i,k))
984  nonlins = (1.0/abs(c2)) * (2.*sqrt(tem) - tem)
985  tem1 = basicui(i,k)
986  tem2 = c2*nonlins*tem1
987  taugwci(i,k) = - rhoi(i,k) * c1 * tem1 * tem2 * tem2
988  & / (bruni(i,k)*dlen(i))
989  elseif (rimin(i,k) > riminm) then
990  taugwci(i,k) = zero
991  end if ! RImin
992  else
993 
997 
998  taugwci(i,k) = zero
999  end if ! RIloc
1000  else
1001  taugwci(i,k) = zero
1002  end if ! TAUGWCI
1003 
1004  if ( (basicum(i,k+1)*basicum(i,k) ) < 0. ) then
1005  taugwci(i,k+1) = zero
1006  taugwci(i,k) = zero
1007  endif
1008 
1009  if (abs(taugwci(i,k)) > abs(taugwci(i,k+1))) then
1010  taugwci(i,k) = taugwci(i,k+1)
1011  end if
1012 
1013  elseif (k == 1) then
1014 
1017 
1018  taugwci(i,1) = taugwci(i,2)
1019  endif
1020  endif
1021  enddo ! end of i=1,npt loop
1022  enddo ! end of k=kcldm,1,-1 loop
1023 
1027 
1028  do k=1,km
1029  do i=1,npt
1030  if (do_gwc(i)) then
1031  kk = kcldtop(i)
1032  if (k < kk) then
1033  wtgwc = (taugwci(i,k+1) - taugwci(i,k)) / dpmid(i,k)
1034  utgwcl(i,k) = wtgwc * cosphi(i)
1035  vtgwcl(i,k) = wtgwc * sinphi(i)
1036  else
1037  utgwcl(i,k) = zero
1038  vtgwcl(i,k) = zero
1039  endif
1040  endif
1041  enddo
1042  enddo
1043 
1044 !-----------------------------------------------------------------------
1045 !
1046 ! Calculate momentum flux = stress deposited above cloup top
1047 ! Apply equal amount with opposite sign within cloud
1048 !
1049 !-----------------------------------------------------------------------
1050 
1051  do i=1,npt
1052  xstress(i) = zero
1053  ystress(i) = zero
1054  enddo
1055  do k=1,kcldm
1056  do i=1,npt
1057  if (do_gwc(i)) then
1058  xstress(i) = xstress(i) + utgwcl(i,k)*dpmid(i,k)
1059  ystress(i) = ystress(i) + vtgwcl(i,k)*dpmid(i,k)
1060  endif
1061  enddo
1062  enddo
1063 
1064 !-----------------------------------------------------------------------
1065 ! ALT 1 ONLY UPPERMOST LAYER
1066 !-----------------------------------------------------------------------
1067 
1068 ! kk = kcldtop(i)
1069 ! tem1 = g / dpmid(i,kk)
1070 ! utgwc(i,kk) = - tem1 * xstress
1071 ! vtgwc(i,kk) = - tem1 * ystress
1072 
1073 !-----------------------------------------------------------------------
1074 ! ALT 2 SIN(KT-KB)
1075 !-----------------------------------------------------------------------
1076 
1077  do i=1,npt
1078  if (do_gwc(i)) then
1079  wrk(i) = 0.5 * pi / (pint(i,kcldbot(i)+1)-pint(i,kcldtop(i)))
1080  endif
1081  enddo
1082  do k=1,km
1083  do i=1,npt
1084  if (do_gwc(i)) then
1085  kk = kcldtop(i)
1086  if (k >= kk .and. k <= kcldbot(i)) then
1087  p1 = sin(wrk(i) * (pint(i,k) -pint(i,kk)))
1088  p2 = sin(wrk(i) * (pint(i,k+1)-pint(i,kk)))
1089  tem = - (p2-p1) / dpmid(i,k)
1090  utgwcl(i,k) = tem*xstress(i)
1091  vtgwcl(i,k) = tem*ystress(i)
1092  endif
1093  endif
1094  enddo
1095  enddo
1096 
1097 !-----------------------------------------------------------------------
1098 ! ALT 3 FROM KT to KB PROPORTIONAL TO CONV HEATING
1099 !-----------------------------------------------------------------------
1100 
1101 ! do k=kcldtop(i),kcldbot(i)
1102 ! p1=cumchr(i,k)
1103 ! p2=cumchr(i,k+1)
1104 ! utgwcl(i,k) = - g*xstress*(p1-p2)/dpmid(i,k)
1105 ! enddo
1106 
1107 !-----------------------------------------------------------------------
1108 !
1109 ! The GWDC should accelerate the zonal and meridional wind in the
1110 ! opposite direction of the previous zonal and meridional wind,
1111 ! respectively
1112 !
1113 !-----------------------------------------------------------------------
1114 
1115 ! do k=1,kcldtop(i)-1
1116 
1117 ! if (utgwcl(i,k)*u(i,k) .gt. 0.0) then
1118 
1119 !-------------------- x-component-------------------
1120 
1121 ! write(6,'(a)')
1122 ! + '(GWDC) WARNING: The GWDC should accelerate the zonal wind '
1123 ! write(6,'(a,a,i3,a,i3)')
1124 ! + 'in the opposite direction of the previous zonal wind',
1125 ! + ' at I = ',i,' and J = ',lat
1126 ! write(6,'(4(1x,e17.10))') u(i,kk),v(i,kk),u(i,k),v(i,k)
1127 ! write(6,'(a,1x,e17.10))') 'Vcld . V =',
1128 ! + u(i,kk)*u(i,k)+v(i,kk)*v(i,k)
1129 
1130 ! if(u(i,kcldtop(i))*u(i,k)+v(i,kcldtop(i))*v(i,k).gt.0.0)then
1131 ! do k1=1,km
1132 ! write(6,'(i2,36x,2(1x,e17.10))')
1133 ! + k1,taugwcxi(i,k1),taugwci(i,k1)
1134 ! write(6,'(i2,2(1x,e17.10))') k1,utgwcl(i,k1),u(i,k1)
1135 ! end do
1136 ! write(6,'(i2,36x,1x,e17.10)') (km+1),taugwcxi(i,km+1)
1137 ! end if
1138 
1139 !-------------------- Along wind at cloud top -----
1140 
1141 ! do k1=1,km
1142 ! write(6,'(i2,36x,2(1x,e17.10))')
1143 ! + k1,taugwci(i,k1)
1144 ! write(6,'(i2,2(1x,e17.10))') k1,wtgwc(i,k1),basicum(i,k1)
1145 ! end do
1146 ! write(6,'(i2,36x,1x,e17.10)') (km+1),taugwci(i,km+1)
1147 
1148 ! end if
1149 
1150 ! if (vtgwc(i,k)*v(i,k) .gt. 0.0) then
1151 ! write(6,'(a)')
1152 ! + '(GWDC) WARNING: The GWDC should accelerate the meridional wind'
1153 ! write(6,'(a,a,i3,a,i3)')
1154 ! + 'in the opposite direction of the previous meridional wind',
1155 ! + ' at I = ',i,' and J = ',lat
1156 ! write(6,'(4(1x,e17.10))') u(i,kcldtop(i)),v(i,kcldtop(i)),
1157 ! + u(i,k),v(i,k)
1158 ! write(6,'(a,1x,e17.10))') 'Vcld . V =',
1159 ! + u(i,kcldtop(i))*u(i,k)+v(i,kcldtop(i))*v(i,k)
1160 ! if(u(i,kcldtop(i))*u(i,k)+v(i,kcldtop(i))*v(i,k).gt.0.0)then
1161 ! do k1=1,km
1162 ! write(6,'(i2,36x,2(1x,e17.10))')
1163 ! + k1,taugwcyi(i,k1),taugwci(i,k1)
1164 ! write(6,'(i2,2(1x,e17.10))') k1,vtgwc(i,k1),v(i,k1)
1165 ! end do
1166 ! write(6,'(i2,36x,1x,e17.10)') (km+1),taugwcyi(i,km+1)
1167 ! end if
1168 ! end if
1169 
1170 ! enddo
1171 
1172 !1000 continue
1173 
1174 
1175 !***********************************************************************
1176 
1177 ! if (lprnt) then
1178 ! if (fhour.ge.fhourpr) then
1179 !-------- UTGWC VTGWC ----------
1180 ! write(*,9220)
1181 ! do ilev=1,km
1182 ! write(*,9221) ilev,(86400.*utgwcl(ipr,ilev)),
1183 ! + (86400.*vtgwcl(ipr,ilev))
1184 ! enddo
1185 ! endif
1186 ! endif
1187 
1188 !9220 format(//,14x,'TENDENCY DUE TO GWDC',//,
1189 ! +' ILEV',6x,'UTGWC',7x,'VTGWC',/)
1190 !9221 format(i4,2(2x,f10.3))
1191 
1192 !-----------------------------------------------------------------------
1193 !
1194 ! For GWDC performance analysis
1195 !
1196 !-----------------------------------------------------------------------
1197 
1198 ! do k = 1, kk-1
1199 ! do i = 1, nct
1200 
1201 ! kk = kcldtop(i)
1202 
1203 ! if ( (abs(taugwci(i,kk)) > taumin) ) then
1204 
1205 ! gwdcloc(i) = one
1206 
1207 ! if ( abs(taugwci(i,k)-taugwci(i,kk)) > taumin ) then
1208 ! break(i) = 1.0
1209 ! go to 2000
1210 ! endif
1211 ! enddo
1212 !2000 continue
1213 
1214 ! do k = 1, kk-1
1215 
1216 ! if ( ( abs(taugwci(i,k)).lt.taumin ) .and.
1217 ! & ( abs(taugwci(i,k+1)).gt.taumin ) .and.
1218 ! & ( basicum(i,k+1)*basicum(i,k) .lt. 0. ) ) then
1219 ! critic(i) = 1.0
1220 ! print *,i,k,' inside GWDC taugwci(k) = ',taugwci(i,k)
1221 ! print *,i,k+1,' inside GWDC taugwci(k+1) = ',taugwci(i,k+1)
1222 ! print *,i,k,' inside GWDC basicum(k) = ',basicum(i,k)
1223 ! print *,i,k+1,' inside GWDC basicum(k+1) = ',basicum(i,k+1)
1224 ! print *,i,' inside GWDC critic = ',critic(i)
1225 ! goto 2010
1226 ! endif
1227 ! enddo
1228 !2010 continue
1229 
1230 ! endif
1231 
1232 ! enddo
1233 
1234 !-----------------------------------------------------------------------
1237 ! Outgoing (FU1,FV1)=(utgwc,vtgwc)
1238 !-----------------------------------------------------------------------
1239 
1240  do k=1,km
1241  k1 = km - k + 1
1242  do i=1,npt
1243  ii = ipt(i)
1244  utgwc(ii,k1) = utgwcl(i,k)
1245 
1246  vtgwc(ii,k1) = vtgwcl(i,k)
1247 
1248 ! brunm(ii,kk) = brunm(i,k)
1249 ! brunm(i,k) = tem
1250 
1251 ! rhom(ii,kk) = rhom(i,k)
1252 
1253  enddo
1254  enddo
1255  do i=1,npt
1256  ii = ipt(i)
1257  tauctx(ii) = tauctxl(i)
1258  taucty(ii) = tauctyl(i)
1259  enddo
1260 
1261 ! if (lprnt) then
1262 ! if (fhour.ge.fhourpr) then
1263 !-------- UTGWC VTGWC ----------
1264 ! write(*,9225)
1265 ! do ilev=km,1,-1
1266 ! write(*,9226) ilev,(86400.*fu1(ipr,ilev)),
1267 ! + (86400.*fv1(ipr,ilev))
1268 ! enddo
1269 ! endif
1270 ! endif
1271 
1272 !9225 format(//,14x,'TENDENCY DUE TO GWDC - TO GBPHYS',//,
1273 ! +' ILEV',6x,'UTGWC',7x,'VTGWC',/)
1274 !9226 format(i4,2(2x,f10.3))
1275 
1276  deallocate (kcldtop,kcldbot,do_gwc)
1277  deallocate (tauctxl, tauctyl,
1278  & gwdcloc, break, critic, cosphi,
1279  & sinphi, xstress, ystress,
1280  & dlen, ucltop, vcltop, gqmcldlen, wrk)
1281 
1282  deallocate (plnint, taugwci,
1283  & bruni, rhoi, basicui,
1284  & ti, riloc, rimin, pint)
1285 
1286  deallocate (plnmid, utgwcl, vtgwcl, basicum, u, v, t,
1287  & pmid, dpmid, brunm, rhom)
1288 
1289  return
1290  end
1291 
gwdc
subroutine gwdc(im, ix, iy, km, lat, u1, v1, t1, q1, pmid1, pint1, dpmid1, qmax, ktop, kbot, kcnv, cldf, grav, cp, rd, fv, dlength, lprnt, ipr, fhour, utgwc, vtgwc, tauctx, taucty)
Definition: gwdc.f:76