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