CCPP SciDoc v7.0.0  v7.0.0
Common Community Physics Package Developed at DTC
 
Loading...
Searching...
No Matches
gwdc.f
1
4
5 module gwdc
6
7 contains
8
9! \brief Brief description of the subroutine
10!
14 subroutine gwdc_init(do_cnvgwd, errmsg, errflg)
15
16 implicit none
17
18 logical, intent(in) :: do_cnvgwd
19 character(len=*), intent(out) :: errmsg
20 integer, intent(out) :: errflg
21
22 ! Initialize CCPP error handling variables
23 errmsg = ''
24 errflg = 0
25
26 if (.not. do_cnvgwd) then
27 errmsg = "Logic error: gwdc called but do_cnvgwd is false"
28 errflg = 1
29 return
30 end if
31
32 end subroutine gwdc_init
33
68 subroutine gwdc_run (im,km,lat,u1,v1,t1,q1,deltim, &
69 & pmid1,pint1,dpmid1,qmax,ktop,kbot,kcnv,cldf, &
70 & grav,cp,rd,fv,pi,dlength,lprnt,ipr,fhour, &
71 & utgwc,vtgwc,tauctx,taucty,errmsg,errflg)
72
73
74!***********************************************************************
75! aug 2005 Ake Johansson - ORIGINAL CODE FOR PARAMETERIZATION OF CONVECTIVELY FORCED
76! GRAVITY WAVE DRAG FROM YONSEI UNIVERSITY, KOREA
77! BASED ON THE THEORY GIVEN BY CHUN AND BAIK (JAS, 1998)
78! MODIFIED FOR IMPLEMENTATION INTO THE GFS/CFSD BY
79! 2013 S. Moorthi - Updated and optimized code for T1534 GFS implementation
80! ??? ?? 2015 J. Alpert - reducing the magnitude of tauctmax to fix blow up in L64 GFS
81! S. Kar & M. Young
82! aug 15 2016 - S. Moorthi - Fix for exessive dissipation which led to blow up in
83! 128 level runs with NEMS/GSM
84!***********************************************************************
85
86 USE machine , ONLY : kind_phys
87 implicit none
88
89!---------------------------- Arguments --------------------------------
90!
91! Input variables
92!
93! u : midpoint zonal wind
94! v : midpoint meridional wind
95! t : midpoint temperatures
96! pmid : midpoint pressures
97! pint : interface pressures
98! dpmid : midpoint delta p ( pi(k)-pi(k-1) )
99! lat : latitude index
100! qmax : deep convective heating
101! kcldtop : Vertical level index for cloud top ( mid level )
102! kcldbot : Vertical level index for cloud bottom ( mid level )
103! kcnv : (0,1) dependent on whether convection occur or not
104!
105! Output variables
106!
107! utgwc : zonal wind tendency
108! vtgwc : meridional wind tendency
109!
110!-----------------------------------------------------------------------
111
112 integer, intent(in) :: im, km, lat, ipr
113 integer, intent(in) :: ktop(:),kbot(:),kcnv(:)
114 real(kind=kind_phys), intent(in) :: grav,cp,rd,fv,fhour,deltim,pi
115 real(kind=kind_phys), dimension(:), intent(in) :: qmax
116 real(kind=kind_phys), dimension(:), intent(out) :: tauctx,taucty
117 real(kind=kind_phys), dimension(:), intent(in) :: cldf,dlength
118 real(kind=kind_phys), dimension(:,:), intent(in) :: u1,v1,t1, &
119 & q1,pmid1,dpmid1
120 real(kind=kind_phys), dimension(:,:), intent(out) :: utgwc,vtgwc
121 real(kind=kind_phys), dimension(:,:), intent(in) :: pint1
122!
123 logical, intent(in) :: lprnt
124!
125 character(len=*), intent(out) :: errmsg
126 integer, intent(out) :: errflg
127
128!------------------------- Local workspace -----------------------------
129!
130! i, k : Loop index
131! kk : Loop index
132! cldf : Deep convective cloud fraction at the cloud top.
133! ugwdc : Zonal wind after GWDC paramterization
134! vgwdc : Meridional wind after GWDC parameterization
135! plnmid : Log(pmid) ( mid level )
136! plnint : Log(pint) ( interface level )
137! dpint : Delta pmid ( interface level )
138! tauct : Wave stress at the cloud top calculated using basic-wind
139! parallel to the wind vector at the cloud top ( mid level )
140! tauctx : Wave stress at the cloud top projected in the east
141! taucty : Wave stress at the cloud top projected in the north
142! qmax : Maximum deep convective heating rate ( K s-1 ) in a
143! horizontal grid point calculated from cumulus para-
144! meterization. ( mid level )
145! wtgwc : Wind tendency in direction to the wind vector at the cloud top level
146! due to convectively generated gravity waves ( mid level )
147! utgwcl : Zonal wind tendency due to convectively generated
148! gravity waves ( mid level )
149! vtgwcl : Meridional wind tendency due to convectively generated
150! gravity waves ( mid level )
151! taugwci : Profile of wave stress calculated using basic-wind
152! parallel to the wind vector at the cloud top
153! taugwcxi : Profile of zonal component of gravity wave stress
154! taugwcyi : Profile of meridional component of gravity wave stress
155!
156! taugwci, taugwcxi, and taugwcyi are defined at the interface level
157!
158! bruni : Brunt-Vaisala frequency ( interface level )
159! brunm : Brunt-Vaisala frequency ( mid level )
160! rhoi : Air density ( interface level )
161! rhom : Air density ( mid level )
162! ti : Temperature ( interface level )
163! basicum : Basic-wind profile. Basic-wind is parallel to the wind
164! vector at the cloud top level. (mid level)
165! basicui : Basic-wind profile. Basic-wind is parallel to the wind
166! vector at the cloud top level. ( interface level )
167! riloc : Local Richardson number ( interface level )
168! rimin : Minimum Richardson number including both the basic-state
169! and gravity wave effects ( interface level )
170! gwdcloc : Horizontal location where the GWDC scheme is activated.
171! break : Horizontal location where wave breaking is occurred.
172! critic : Horizontal location where critical level filtering is
173! occurred.
174! dogwdc : Logical flag whether the GWDC parameterization is
175! calculated at a grid point or not.
176!
177! dogwdc is used in order to lessen CPU time for GWDC calculation.
178!
179!-----------------------------------------------------------------------
180
181 integer i,ii,k,k1,kk,kb,ilev,npt,kcb,kcldm,npr
182 integer, dimension(im) :: ipt
183
184 real(kind=kind_phys) tem, tem1, tem2, qtem, wtgwc, tauct, &
185 & windcltop, shear, nonlinct, nonlin, nonlins,&
186 & n2, dtdp, crit1, crit2, p1, p2, &
187! & n2, dtdp, crit1, crit2, pi, p1, p2,
188 & gsqr, onebg
189! & taus, n2, dtdp, crit1, crit2, pi, p1, p2
190
191 integer, allocatable :: kcldtop(:),kcldbot(:)
192 logical, allocatable :: do_gwc(:)
193 real(kind=kind_phys), allocatable :: tauctxl(:), tauctyl(:),
194 & gwdcloc(:), break(:),
195! & critic(:),
196! & critic(:), angle(:),
197 & cosphi(:), sinphi(:),
198 & xstress(:), ystress(:),
199 & ucltop(:), vcltop(:),
200 & wrk(:), dtfac(:),
201 & dlen(:), gqmcldlen(:)
202! real(kind=kind_phys), allocatable :: plnint(:,:), dpint(:,:),
203! & taugwci(:,:), taugwcxi(:,:),
204! & taugwcyi(:,:), bruni(:,:),
205! & taugwcyi(:,:), bruni(:,:),
206 real(kind=kind_phys), allocatable :: plnint(:,:), velco(:,:),
207 & taugwci(:,:), bruni(:,:),
208 & rhoi(:,:), basicui(:,:),
209 & ti(:,:), riloc(:,:),
210 & rimin(:,:), pint(:,:)
211! real(kind=kind_phys), allocatable :: ugwdc(:,:), vgwdc(:,:),
212 real(kind=kind_phys), allocatable ::
213! & plnmid(:,:), wtgwc(:,:),
214 & plnmid(:,:), taugw(:,:),
215 & utgwcl(:,:), vtgwcl(:,:),
216 & basicum(:,:), u(:,:),v(:,:),
217 & t(:,:), spfh(:,:),
218 & pmid(:,:), dpmid(:,:),
219! & pmid(:,:), cumchr(:,:),
220 & brunm(:,:), rhom(:,:)
221
222!-----------------------------------------------------------------------
223!
224! ucltop : Zonal wind at the cloud top ( mid level )
225! vcltop : Meridional wind at the cloud top ( mid level )
226! windcltop : Wind speed at the cloud top ( mid level )
227! shear : Vertical shear of basic wind
228! cosphi : Cosine of angle of wind vector at the cloud top
229! sinphi : Sine of angle of wind vector at the cloud top
230! c1 : Tunable parameter
231! c2 : Tunable parameter
232! dlength : Grid spacing in the direction of basic wind at the cloud top
233! nonlinct : Nonlinear parameter at the cloud top
234! nonlin : Nonlinear parameter above the cloud top
235! nonlins : Saturation nonlinear parameter
236! taus : Saturation gravity wave drag == taugwci(i,k)
237! n2 : Square of Brunt-Vaisala frequency
238! dtdp : dT/dp
239! xstress : Vertically integrated zonal momentum change due to GWDC
240! ystress : Vertically integrated meridional momentum change due to GWDC
241! crit1 : Variable 1 for checking critical level
242! crit2 : Variable 2 for checking critical level
243!
244!-----------------------------------------------------------------------
245
246 real(kind=kind_phys), parameter ::
247 & c1=1.41, c2=-0.38, ricrit=0.25
248 &, n2min=1.e-32, zero=0.0, one=1.0
249 &, taumin=1.0e-20, tauctmax=-20.
250! &, taumin=1.0e-20, tauctmax=-5.
251 &, qmin=1.0e-10, shmin=1.0e-20
252 &, rimax=1.0e+20, rimaxm=0.99e+20
253 &, rimaxp=1.01e+20, rilarge=0.9e+20
254 &, riminx=-1.0e+20, riminm=-1.01e+20
255 &, riminp=-0.99e+20, rismall=-0.9e+20
256
257 ! Initialize CCPP error handling variables
258 errmsg = ''
259 errflg = 0
260
261!
263 npt = 0
264 do i = 1,im
265 ipt(i) = 0
266 if (kcnv(i) /= 0 .and. qmax(i) > zero) then
267 npt = npt + 1
268 ipt(npt) = i
269 endif
270 enddo
271 do k=1,km
272 do i=1,im
273 utgwc(i,k) = 0.0
274 vtgwc(i,k) = 0.0
275! brunm(i,k) = 0.0
276! rhom(i,k) = 0.0
277 enddo
278 enddo
279 do i=1,im
280 tauctx(i) = 0.0
281 taucty(i) = 0.0
282 enddo
283 if (npt == 0) return ! No gwdc calculation done!
284
285!***********************************************************************
286!
287! Begin GWDC
288!
289!***********************************************************************
290
291!-----------------------------------------------------------------------
292! Write out incoming variables
293!-----------------------------------------------------------------------
294
295! fhourpr = zero
296! if (lprnt) then
297! if (fhour >= fhourpr) then
298! print *,' '
299! write(*,*) 'Inside GWDC raw input start print at fhour = ',
300! & fhour
301! write(*,*) 'IM KM ',im,km
302! write(*,*) 'KBOT KTOP QMAX DLENGTH kcnv ',
303! + kbot(ipr),ktop(ipr),qmax(ipr),dlength(ipr),kcnv(ipr)
304! write(*,*) 'grav cp rd ',grav,cp,rd
305
306!-------- Pressure levels ----------
307! write(*,9100)
308! ilev=km+1
309! write(*,9110) ilev,(10.*pint1(ipr,ilev))
310! do ilev=km,1,-1
311! write(*,9120) ilev,(10.*pmid1(ipr,ilev)),
312! & (10.*dpmid1(ipr,ilev))
313! write(*,9110) ilev,(10.*pint1(ipr,ilev))
314! enddo
315
316!-------- U1 V1 T1 ----------
317! write(*,9130)
318! do ilev=km,1,-1
319! write(*,9140) ilev,U1(ipr,ilev),V1(ipr,ilev),T1(ipr,ilev)
320! enddo
321
322! print *,' '
323! print *,' Inside GWDC raw input end print'
324! endif
325! endif
326
327!9100 format(//,14x,'PRESSURE LEVELS',//,
328! +' ILEV',6x,'PINT1',7x,'PMID1',6x,'DPMID1',/)
329!9110 format(i4,2x,f10.3)
330!9120 format(i4,12x,2(2x,f10.3))
331!9130 format(//,' ILEV',7x,'U1',10x,'V1',10x,'T1',/)
332!9140 format(i4,3(2x,f10.3))
333
334! Allocate local arrays
335
336 allocate (kcldtop(npt), kcldbot(npt), do_gwc(npt))
337 allocate (tauctxl(npt), tauctyl(npt), dtfac(npt),
338 & gwdcloc(npt), break(npt), cosphi(npt),
339! & gwdcloc(npt), break(npt), critic(npt), cosphi(npt),
340 & sinphi(npt), xstress(npt), ystress(npt), wrk(npt),
341 & ucltop(npt), vcltop(npt),dlen(npt), gqmcldlen(npt))
342
343! allocate (plnint(npt,2:km+1), dpint(npt,km+1),
344! & taugwci(npt,km+1), taugwcxi(npt,km+1),
345! & taugwcyi(npt,km+1), bruni(npt,km+1),
346 allocate (plnint(npt,2:km+1),
347 & taugwci(npt,km+1), bruni(npt,km+1),
348 & rhoi(npt,km+1), basicui(npt,km+1),
349 & ti(npt,km+1), riloc(npt,km+1),
350 & rimin(npt,km+1), pint(npt,km+1))
351
352! allocate (ugwdc(npt,km), vgwdc(npt,km),
353 allocate
354! & (plnmid(npt,km), wtgwc(npt,km),
355 & (plnmid(npt,km), velco(npt,km),
356 & utgwcl(npt,km), vtgwcl(npt,km),
357 & basicum(npt,km), u(npt,km), v(npt,km),
358 & t(npt,km), spfh(npt,km), pmid(npt,km),
359 & dpmid(npt,km), taugw(npt,km),
360! & dpmid(npt,km), cumchr(npt,km),
361 & brunm(npt,km), rhom(npt,km))
362
363!-----------------------------------------------------------------------
366!-----------------------------------------------------------------------
367 gsqr = grav * grav
368 onebg = one / grav
369
370 if (lprnt) then
371 npr = 1
372 do i=1,npt
373 if (ipr == ipt(i))then
374 npr = i
375 exit
376 endif
377 enddo
378 endif
379
380 do k=1,km
381 k1 = km - k + 1
382 do i=1,npt
383 ii = ipt(i)
384 u(i,k) = u1(ii,k1)
385 v(i,k) = v1(ii,k1)
386 t(i,k) = t1(ii,k1)
387 spfh(i,k) = max(q1(ii,k1),qmin)
388 pmid(i,k) = pmid1(ii,k1)
389 dpmid(i,k) = dpmid1(ii,k1) * onebg
390! cumchr(i,k) = cumchr1(ii,k1)
391
392 rhom(i,k) = pmid(i,k) / (rd*t(i,k)*(1.0+fv*spfh(i,k)))
393 plnmid(i,k) = log(pmid(i,k))
394 utgwcl(i,k) = zero
395 vtgwcl(i,k) = zero
396! ugwdc(i,k) = zero
397! vgwdc(i,k) = zero
398 brunm(i,k) = zero
399 basicum(i,k) = zero
400 enddo
401 enddo
402
403 do k=1,km+1
404 k1 = km - k + 2
405 do i=1,npt
406 ii = ipt(i)
407 pint(i,k) = pint1(ii,k1)
408 taugwci(i,k) = zero
409 bruni(i,k) = zero
410 rhoi(i,k) = zero
411 ti(i,k) = zero
412 basicui(i,k) = zero
413 riloc(i,k) = zero
414 rimin(i,k) = zero
415 enddo
416 enddo
417 do k=2,km+1
418 do i=1,npt
419 plnint(i,k) = log(pint(i,k))
420 enddo
421 enddo
422
423 do i = 1, npt
424 ii = ipt(i)
425 kcldtop(i) = km - ktop(ii) + 1
426 kcldbot(i) = km - kbot(ii) + 1
427 dlen(i) = dlength(ii)
428! (g*qmax(ii)*cldf(ii)*dlength(ii))
429 gqmcldlen(i) = grav*qmax(ii)*cldf(ii)*dlen(i)
430 enddo
431! if (lprnt) write(7000,*)' ktop=',ktop(ipr),' kbot=',kbot(ipr)
432! &,' kcldtop=',kcldtop(npr),' kcldbot=',kcldbot(npr),
433! &' dlength=',dlength(ipr),' qmax=',qmax(ipr),' cldf=',cldf(ipr)
434
435! if (lprnt) then
436! if (fhour.ge.fhourpr) then
437! write(*,9200)
438! do i=1,im
439! write(*,9201) kcnv(i),kcldbot(i),kcldtop(i)
440! enddo
441! endif
442! endif
443
444!9200 format(//,' Inside GWDC local variables start print',//,
445! +2x,'kcnv',2x,'KCLDBOT',2x,'KCLDTOP',//)
446!9201 format(i4,2x,i5,4x,i5)
447
448!***********************************************************************
449
450! pi = 2.*asin(1.)
451
452!-----------------------------------------------------------------------
453!
454! PRESSURE VARIABLES
455!
456! Interface 1 ======== pint(1) *********
457! Mid-Level 1 -------- pmid(1) dpmid(1)
458! 2 ======== pint(2) dpint(2)
459! 2 -------- pmid(2) dpmid(2)
460! 3 ======== pint(3) dpint(3)
461! 3 -------- pmid(3) dpmid(3)
462! 4 ======== pint(4) dpint(4)
463! 4 -------- pmid(4) dpmid(4)
464! ........
465! 17 ======== pint(17) dpint(17)
466! 17 -------- pmid(17) dpmid(17)
467! 18 ======== pint(18) dpint(18)
468! 18 -------- pmid(18) dpmid(18)
469! 19 ======== pint(19) *********
470!
471!-----------------------------------------------------------------------
472
473 do i = 1, npt
474 tauctxl(i) = zero
475 tauctyl(i) = zero
476
477!-----------------------------------------------------------------------
478! THERMAL VARIABLES
479!
480! Interface 1 ======== TI(1) RHOI(1) BRUNI(1)
481! 1 -------- T(1) RHOM(1) BRUNM(1)
482! 2 ======== TI(2) RHOI(2) BRUNI(2)
483! 2 -------- T(2) RHOM(2) BRUNM(2)
484! 3 ======== TI(3) RHOI(3) BRUNI(3)
485! 3 -------- T(3) RHOM(3) BRUNM(3)
486! 4 ======== TI(4) RHOI(4) BRUNI(4)
487! 4 -------- T(4) RHOM(4) BRUNM(4)
488! ........
489! 17 ========
490! 17 -------- T(17) RHOM(17) BRUNM(17)
491! 18 ======== TI(18) RHOI(18) BRUNI(18)
492! 18 -------- T(18) RHOM(18) BRUNM(18)
493! 19 ======== TI(19) RHOI(19) BRUNI(19)
494!
495
496!
500
501 ti(i,1) = t(i,1)
502 rhoi(i,1) = pint(i,1)/(rd*ti(i,1))
503 bruni(i,1) = sqrt( gsqr / (cp*ti(i,1)) )
504!
508
509 ti(i,km+1) = t(i,km)
510 rhoi(i,km+1) = pint(i,km+1)/(rd*ti(i,km+1)*(1.0+fv*spfh(i,km)))
511 bruni(i,km+1) = sqrt( gsqr / (cp*ti(i,km+1)) )
512 enddo
513
514!-----------------------------------------------------------------------
515!
519!
520!-----------------------------------------------------------------------
521
522 do k = 2, km
523 do i = 1, npt
524 tem1 = (plnmid(i,k)-plnint(i,k)) / (plnmid(i,k)-plnmid(i,k-1))
525 tem2 = one - tem1
526 ti(i,k) = t(i,k-1) * tem1 + t(i,k) * tem2
527 qtem = spfh(i,k-1) * tem1 + spfh(i,k) * tem2
528 rhoi(i,k) = pint(i,k) / ( rd * ti(i,k)*(1.0+fv*qtem) )
529 dtdp = (t(i,k)-t(i,k-1)) / (pmid(i,k)-pmid(i,k-1))
530 n2 = gsqr / ti(i,k) * ( 1./cp - rhoi(i,k)*dtdp )
531 bruni(i,k) = sqrt(max(n2min, n2))
532 enddo
533 enddo
534
535 deallocate (spfh)
536!-----------------------------------------------------------------------
537!
540!-----------------------------------------------------------------------
541
542 do k = 1, km
543 do i = 1, npt
544 dtdp = (ti(i,k+1)-ti(i,k)) / (pint(i,k+1)-pint(i,k))
545 n2 = gsqr / t(i,k) * ( 1./cp - rhom(i,k)*dtdp )
546 brunm(i,k) = sqrt(max(n2min, n2))
547 enddo
548 enddo
549
550!-----------------------------------------------------------------------
551! PRINTOUT
552!-----------------------------------------------------------------------
553
554! if (lprnt) then
555! if (fhour.ge.fhourpr) then
556
557!-------- Pressure levels ----------
558! write(*,9101)
559! do ilev=1,km
560! write(*,9111) ilev,(0.01*pint(ipr,ilev)),
561! & (0.01*dpint(ipr,ilev)),plnint(ipr,ilev)
562! write(*,9121) ilev,(0.01*pmid(ipr,ilev)),
563! & (0.01*dpmid(ipr,ilev)),plnmid(ipr,ilev)
564! enddo
565! ilev=km+1
566! write(*,9111) ilev,(0.01*pint(ipr,ilev)),
567! & (0.01*dpint(ipr,ilev)),plnint(ipr,ilev)
568
569! 2
570!-------- U V T N ----------
571! write(*,9102)
572! do ilev=1,km
573! write(*,9112) ilev,ti(ipr,ilev),(100.*bruni(ipr,ilev))
574! write(*,9122) ilev,u(ipr,ilev),v(ipr,ilev),
575! + t(ipr,ilev),(100.*brunm(ipr,ilev))
576! enddo
577! ilev=km+1
578! write(*,9112) ilev,ti(ipr,ilev),(100.*bruni(ipr,ilev))
579
580! endif
581! endif
582
583!9101 format(//,14x,'PRESSURE LEVELS',//,
584! +' ILEV',4x,'PINT',4x,'PMID',4x,'DPINT',3x,'DPMID',5x,'LNP',/)
585!9111 format(i4,1x,f8.2,9x,f8.2,9x,f8.2)
586!9121 format(i4,9x,f8.2,9x,f8.2,1x,f8.2)
587!9102 format(//' ILEV',5x,'U',7x,'V',5x,'TI',7x,'T',
588! +5x,'BRUNI',3x,'BRUNM',//)
589!9112 format(i4,16x,f8.2,8x,f8.3)
590!9122 format(i4,2f8.2,8x,f8.2,8x,f8.3)
591
592
593!***********************************************************************
594!
595! Big loop over grid points ONLY done if kcnv=1
596!
597!***********************************************************************
598
599 kcldm = 1
600 do i = 1, npt
601 kk = kcldtop(i)
602 kb = kcldbot(i)
603 kcldm = max(kcldm,kk)
604
605!-----------------------------------------------------------------------
606!
608! Here, ucltop, vcltop, and windcltop are wind components and
609! wind speed at mid-level cloud top index
610!
611!-----------------------------------------------------------------------
612
613 ucltop(i) = u(i,kk)
614 vcltop(i) = v(i,kk)
615! windcltop = sqrt( ucltop(i)*ucltop(i) + vcltop(i)*vcltop(i) )
616 windcltop = 1.0 / sqrt( ucltop(i)*ucltop(i)
617 & + vcltop(i)*vcltop(i) )
618 cosphi(i) = ucltop(i)*windcltop
619 sinphi(i) = vcltop(i)*windcltop
620! angle(i) = acos(cosphi)*180./pi
621 enddo
622
623!-----------------------------------------------------------------------
624!
627! \n U : Basic-wind speed profile. Basic-wind is parallel to the wind
628! vector at the cloud top level. (mid level)
629! \n UI: Basic-wind speed profile. Basic-wind is parallel to the wind
630! vector at the cloud top level. ( interface level )
631! Input u(i,k) and v(i,k) is defined at mid level
632!
633!-----------------------------------------------------------------------
634
635 do k=1,km
636 do i=1,npt
637 basicum(i,k) = u(i,k)*cosphi(i) + v(i,k)*sinphi(i)
638 enddo
639 enddo
640
641!-----------------------------------------------------------------------
642!
643! Basic state wind at interface level is also calculated
644! based on linear interpolation in ln(Pressure)
645!
646! In the top and bottom boundaries, basic-state wind at interface level
647! is assumed to be vertically uniform.
648!
649!-----------------------------------------------------------------------
650
651 do i=1,npt
652 basicui(i,1) = basicum(i,1)
653 basicui(i,km+1) = basicum(i,km)
654 enddo
655 do k=2,km
656 do i=1,npt
657 tem1 = (plnmid(i,k)-plnint(i,k)) / (plnmid(i,k)-plnmid(i,k-1))
658 tem2 = one - tem1
659 basicui(i,k) = basicum(i,k)*tem2 + basicum(i,k-1)*tem1
660 enddo
661 enddo
662
663!-----------------------------------------------------------------------
664!
670
671! basicum : U at mid level
672! basicui : UI at interface level
673!
674! Interface 1 ======== UI(1) rhoi(1) bruni(1) riloc(1)
675! Mid-level 1 -------- U(1)
676! 2 ======== UI(2) dpint(2) rhoi(2) bruni(2) riloc(2)
677! 2 -------- U(2)
678! 3 ======== UI(3) dpint(3) rhoi(3) bruni(3) riloc(3)
679! 3 -------- U(3)
680! 4 ======== UI(4) dpint(4) rhoi(4) bruni(4) riloc(4)
681! 4 -------- U(4)
682! ........
683! 17 ======== UI(17) dpint(17) rhoi(17) bruni(17) riloc(17)
684! 17 -------- U(17)
685! 18 ======== UI(18) dpint(18) rhoi(18) bruni(18) riloc(18)
686! 18 -------- U(18)
687! 19 ======== UI(19) rhoi(19) bruni(19) riloc(19)
688!
689!-----------------------------------------------------------------------
690
691 do k=2,km
692 do i=1,npt
693 shear = grav*rhoi(i,k) * (basicum(i,k) - basicum(i,k-1))
694 & / (pmid(i,k) - pmid(i,k-1))
695 if ( abs(shear) < shmin ) then
696 riloc(i,k) = rimax
697 else
698 tem = bruni(i,k) / shear
699 riloc(i,k) = tem * tem
700 if (riloc(i,k) >= rimax ) riloc(i,k) = rilarge
701 end if
702 enddo
703 enddo
704
705 do i=1,npt
706 riloc(i,1) = riloc(i,2)
707 riloc(i,km+1) = riloc(i,km)
708 enddo
709
710! if (lprnt.and.(i.eq.ipr)) then
711! if (fhour.ge.fhourpr) then
712! write(*,9104) ucltop,vcltop,windcltop,angle,kk
713! do ilev=1,km
714! write(*,9114) ilev,basicui(ipr,ilev),dpint(ipr,ilev),
715! + rhoi(ipr,ilev),(100.*bruni(ipr,ilev)),riloc(ilev)
716! write(*,9124) ilev,(basicum(ipr,ilev))
717! enddo
718! ilev=km+1
719! write(*,9114) ilev,basicui(ipr,ilev),dpint(ipr,ilev),
720! + rhoi(ipr,ilev),(100.*bruni(ipr,ilev)),riloc(ilev)
721! endif
722! endif
723
724!9104 format(//,'WIND VECTOR AT CLOUDTOP = (',f6.2,' , ',f6.2,' ) = ',
725! +f6.2,' IN DIRECTION ',f6.2,4x,'KK = ',i2,//,
726! +' ILEV',2x,'BASICUM',2x,'BASICUI',4x,'DPINT',6x,'RHOI',5x,
727! +'BRUNI',6x,'RI',/)
728!9114 format(i4,10x,f8.2,4(2x,f8.2))
729!9124 format(i4,1x,f8.2)
730
731!-----------------------------------------------------------------------
732!
734!
735! kcldtopi : The interface level cloud top index
736! kcldtop : The midlevel cloud top index
737! kcldbot : The midlevel cloud bottom index
738!
739! A : Find deep convective heating rate maximum
740!
741! If kcldtop(i) is less than kcldbot(i) in a horizontal grid point,
742! it can be thought that there is deep convective cloud. However,
743! deep convective heating between kcldbot and kcldtop is sometimes
744! zero in spite of kcldtop less than kcldbot. In this case,
745! maximum deep convective heating is assumed to be 1.e-30.
746!
747! B : kk is the vertical index for interface level cloud top
748!
749! C : Total convective fractional cover (cldf) is used as the
750! convective cloud cover for GWDC calculation instead of
751! convective cloud cover in each layer (concld).
752! a1 = cldf*dlength
753! You can see the difference between cldf(i) and concld(i)
754! in (4.a.2) in Description of the NCAR Community Climate
755! Model (CCM3).
756! In NCAR CCM3, cloud fractional cover in each layer in a deep
757! cumulus convection is determined assuming total convective
758! cloud cover is randomly overlapped in each layer in the
759! cumulus convection.
760!
761! D : Wave stress at cloud top is calculated when the atmosphere
762! is dynamically stable at the cloud top
763!
764! E : Cloud top wave stress and nonlinear parameter are calculated
765! using density, temperature, and wind that are defined at mid
766! level just below the interface level in which cloud top wave
767! stress is defined.
768! Nonlinct is defined at the interface level.
769!
770! F : If the atmosphere is dynamically unstable at the cloud top,
771! GWDC calculation in current horizontal grid is skipped.
772!
773! G : If mean wind at the cloud top is less than zero, GWDC
774! calculation in current horizontal grid is skipped.
775!
776!-----------------------------------------------------------------------
777!D
813! - If mean wind at the cloud top is less than zero, GWDC
814! calculation in current horizontal grid is skipped.
815!
816
819!
820!-----------------------------------------------------------------------
821!D
822 do i=1,npt
823 kk = kcldtop(i)
824 if ( abs(basicui(i,kk)) > zero .and. riloc(i,kk) > ricrit) then
825!E
826 tem = basicum(i,kk)
827 tem1 = tem * tem
828 nonlinct = gqmcldlen(i) / (bruni(i,kk)*t(i,kk)*tem1) ! Mu
829 tem2 = c2*nonlinct
830! RhoU^3c1(c2mu)^2/Ndx
831 tauct = - rhom(i,kk) * tem * tem1 * c1 * tem2 * tem2
832 & / (bruni(i,kk)*dlen(i))
835
836 tauct = max(tauctmax, tauct)
837 tauctxl(i) = tauct * cosphi(i) ! X stress at cloud top
838 tauctyl(i) = tauct * sinphi(i) ! Y stress at cloud top
839 taugwci(i,kk) = tauct ! *1
840 do_gwc(i) = .true.
841 else
842!F
847
848
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!!!!!!
1090
1091 do k=1,km
1092 do i=1,npt
1093 if (do_gwc(i)) then
1094 kk = kcldtop(i)
1095 if (k < kk) then
1096! wtgwc = (taugwci(i,k+1) - taugwci(i,k)) / dpmid(i,k)
1097 wtgwc = taugw(i,k) * dtfac(i)
1098 utgwcl(i,k) = wtgwc * cosphi(i)
1099 vtgwcl(i,k) = wtgwc * sinphi(i)
1100 else
1101 utgwcl(i,k) = zero
1102 vtgwcl(i,k) = zero
1103 endif
1104! if(lprnt .and. i == npr) then
1105! write(7000,*)'k=',k,' wtgwc=',wtgwc,' taugwci=',taugwci(i,k),
1106! &taugwci(i,k+1),' dpmid=',dpmid(i,k),' cosphi=',cosphi(i),
1107! & ' sinphi=',sinphi(i),' utgwcl=',utgwcl(i,k),
1108! &'vtgwcl=',vtgwcl(i,k),' dtfac=',dtfac(i)
1109! endif
1110 endif
1111 enddo
1112 enddo
1113
1114!-----------------------------------------------------------------------
1115!
1116! Calculate momentum flux = stress deposited above cloup top
1117! Apply equal amount with opposite sign within cloud
1118!
1119!-----------------------------------------------------------------------
1120
1121 do i=1,npt
1122 xstress(i) = zero
1123 ystress(i) = zero
1124 enddo
1125 do k=1,kcldm
1126 do i=1,npt
1127 if (do_gwc(i)) then
1128 xstress(i) = xstress(i) + utgwcl(i,k)*dpmid(i,k)
1129 ystress(i) = ystress(i) + vtgwcl(i,k)*dpmid(i,k)
1130 endif
1131 enddo
1132 enddo
1133
1134!-----------------------------------------------------------------------
1135! ALT 1 ONLY UPPERMOST LAYER
1136!-----------------------------------------------------------------------
1137
1138! kk = kcldtop(i)
1139! tem1 = g / dpmid(i,kk)
1140! utgwc(i,kk) = - tem1 * xstress
1141! vtgwc(i,kk) = - tem1 * ystress
1142
1143!-----------------------------------------------------------------------
1144! ALT 2 SIN(KT-KB)
1145!-----------------------------------------------------------------------
1146
1147 do i=1,npt
1148 if (do_gwc(i)) then
1149 wrk(i) = 0.5 * pi / (pint(i,kcldbot(i)+1)-pint(i,kcldtop(i)))
1150 endif
1151 enddo
1152 do k=1,km
1153 do i=1,npt
1154 if (do_gwc(i)) then
1155 kk = kcldtop(i)
1156 if (k >= kk .and. k <= kcldbot(i)) then
1157 p1 = sin(wrk(i) * (pint(i,k) -pint(i,kk)))
1158 p2 = sin(wrk(i) * (pint(i,k+1)-pint(i,kk)))
1159 tem = - (p2-p1) / dpmid(i,k)
1160 utgwcl(i,k) = tem*xstress(i)
1161 vtgwcl(i,k) = tem*ystress(i)
1162 endif
1163 endif
1164 enddo
1165 enddo
1166
1167!-----------------------------------------------------------------------
1168! ALT 3 FROM KT to KB PROPORTIONAL TO CONV HEATING
1169!-----------------------------------------------------------------------
1170
1171! do k=kcldtop(i),kcldbot(i)
1172! p1=cumchr(i,k)
1173! p2=cumchr(i,k+1)
1174! utgwcl(i,k) = - g*xstress*(p1-p2)/dpmid(i,k)
1175! enddo
1176
1177!-----------------------------------------------------------------------
1178!
1179! The GWDC should accelerate the zonal and meridional wind in the
1180! opposite direction of the previous zonal and meridional wind,
1181! respectively
1182!
1183!-----------------------------------------------------------------------
1184
1185! do k=1,kcldtop(i)-1
1186! if (utgwcl(i,k)*u(i,k) .gt. 0.0) then
1187!-------------------- x-component-------------------
1188! write(6,'(a)')
1189! + '(GWDC) WARNING: The GWDC should accelerate the zonal wind '
1190! write(6,'(a,a,i3,a,i3)')
1191! + 'in the opposite direction of the previous zonal wind',
1192! + ' at I = ',i,' and J = ',lat
1193! write(6,'(4(1x,e17.10))') u(i,kk),v(i,kk),u(i,k),v(i,k)
1194! write(6,'(a,1x,e17.10))') 'Vcld . V =',
1195! + u(i,kk)*u(i,k)+v(i,kk)*v(i,k)
1196
1197! if(u(i,kcldtop(i))*u(i,k)+v(i,kcldtop(i))*v(i,k).gt.0.0)then
1198! do k1=1,km
1199! write(6,'(i2,36x,2(1x,e17.10))')
1200! + k1,taugwcxi(i,k1),taugwci(i,k1)
1201! write(6,'(i2,2(1x,e17.10))') k1,utgwcl(i,k1),u(i,k1)
1202! end do
1203! write(6,'(i2,36x,1x,e17.10)') (km+1),taugwcxi(i,km+1)
1204! end if
1205
1206!-------------------- Along wind at cloud top -----
1207
1208! do k1=1,km
1209! write(6,'(i2,36x,2(1x,e17.10))')
1210! + k1,taugwci(i,k1)
1211! write(6,'(i2,2(1x,e17.10))') k1,wtgwc(i,k1),basicum(i,k1)
1212! end do
1213! write(6,'(i2,36x,1x,e17.10)') (km+1),taugwci(i,km+1)
1214
1215! end if
1216
1217! if (vtgwc(i,k)*v(i,k) .gt. 0.0) then
1218! write(6,'(a)')
1219! + '(GWDC) WARNING: The GWDC should accelerate the meridional wind'
1220! write(6,'(a,a,i3,a,i3)')
1221! + 'in the opposite direction of the previous meridional wind',
1222! + ' at I = ',i,' and J = ',lat
1223! write(6,'(4(1x,e17.10))') u(i,kcldtop(i)),v(i,kcldtop(i)),
1224! + u(i,k),v(i,k)
1225! write(6,'(a,1x,e17.10))') 'Vcld . V =',
1226! + u(i,kcldtop(i))*u(i,k)+v(i,kcldtop(i))*v(i,k)
1227! if(u(i,kcldtop(i))*u(i,k)+v(i,kcldtop(i))*v(i,k).gt.0.0)then
1228! do k1=1,km
1229! write(6,'(i2,36x,2(1x,e17.10))')
1230! + k1,taugwcyi(i,k1),taugwci(i,k1)
1231! write(6,'(i2,2(1x,e17.10))') k1,vtgwc(i,k1),v(i,k1)
1232! end do
1233! write(6,'(i2,36x,1x,e17.10)') (km+1),taugwcyi(i,km+1)
1234! end if
1235! end if
1236
1237! enddo
1238
1239!1000 continue
1240
1241
1242!***********************************************************************
1243
1244! if (lprnt) then
1245! if (fhour.ge.fhourpr) then
1246!-------- UTGWC VTGWC ----------
1247! write(*,9220)
1248! do ilev=1,km
1249! write(*,9221) ilev,(86400.*utgwcl(ipr,ilev)),
1250! + (86400.*vtgwcl(ipr,ilev))
1251! enddo
1252! endif
1253! endif
1254
1255!9220 format(//,14x,'TENDENCY DUE TO GWDC',//,
1256! +' ILEV',6x,'UTGWC',7x,'VTGWC',/)
1257!9221 format(i4,2(2x,f10.3))
1258
1259!-----------------------------------------------------------------------
1260!
1261! For GWDC performance analysis
1262!
1263!-----------------------------------------------------------------------
1264
1265! do k = 1, kk-1
1266! do i = 1, nct
1267
1268! kk = kcldtop(i)
1269
1270! if ( (abs(taugwci(i,kk)) > taumin) ) then
1271
1272! gwdcloc(i) = one
1273
1274! if ( abs(taugwci(i,k)-taugwci(i,kk)) > taumin ) then
1275! break(i) = 1.0
1276! go to 2000
1277! endif
1278! enddo
1279!2000 continue
1280
1281! do k = 1, kk-1
1282
1283! if ( ( abs(taugwci(i,k)).lt.taumin ) .and.
1284! & ( abs(taugwci(i,k+1)).gt.taumin ) .and.
1285! & ( basicum(i,k+1)*basicum(i,k) .lt. 0. ) ) then
1286! critic(i) = 1.0
1287! print *,i,k,' inside GWDC taugwci(k) = ',taugwci(i,k)
1288! print *,i,k+1,' inside GWDC taugwci(k+1) = ',taugwci(i,k+1)
1289! print *,i,k,' inside GWDC basicum(k) = ',basicum(i,k)
1290! print *,i,k+1,' inside GWDC basicum(k+1) = ',basicum(i,k+1)
1291! print *,i,' inside GWDC critic = ',critic(i)
1292! goto 2010
1293! endif
1294! enddo
1295!2010 continue
1296
1297! endif
1298
1299! enddo
1300
1301!-----------------------------------------------------------------------
1304! Outgoing (FU1,FV1)=(utgwc,vtgwc)
1305!-----------------------------------------------------------------------
1306
1307 do k=1,km
1308 k1 = km - k + 1
1309 do i=1,npt
1310 ii = ipt(i)
1311 utgwc(ii,k1) = utgwcl(i,k)
1312
1313 vtgwc(ii,k1) = vtgwcl(i,k)
1314
1315! brunm(ii,kk) = brunm(i,k)
1316! brunm(i,k) = tem
1317
1318! rhom(ii,kk) = rhom(i,k)
1319
1320 enddo
1321! if (lprnt) write(7000,*)' k=',k,' k1=',k1,' utgwc='
1322! &, utgwc(ipr,k1),' vtgwc=',vtgwc(ipr,k1)
1323 enddo
1324 do i=1,npt
1325 ii = ipt(i)
1326 tauctx(ii) = tauctxl(i)
1327 taucty(ii) = tauctyl(i)
1328 enddo
1329
1330! if (lprnt) then
1331! if (fhour.ge.fhourpr) then
1332!-------- UTGWC VTGWC ----------
1333! write(*,9225)
1334! do ilev=km,1,-1
1335! write(*,9226) ilev,(86400.*fu1(ipr,ilev)),
1336! + (86400.*fv1(ipr,ilev))
1337! enddo
1338! endif
1339! endif
1340
1341!9225 format(//,14x,'TENDENCY DUE TO GWDC - TO GBPHYS',//,
1342! +' ILEV',6x,'UTGWC',7x,'VTGWC',/)
1343!9226 format(i4,2(2x,f10.3))
1344
1345 deallocate (kcldtop,kcldbot,do_gwc)
1346 deallocate (tauctxl, tauctyl, dtfac,
1347! & gwdcloc, break, critic, cosphi,
1348 & gwdcloc, break, cosphi,
1349 & sinphi, xstress, ystress,
1350 & dlen, ucltop, vcltop, gqmcldlen, wrk)
1351
1352 deallocate (plnint, taugwci, velco,
1353 & bruni, rhoi, basicui,
1354 & ti, riloc, rimin, pint)
1355
1356 deallocate (plnmid, utgwcl, vtgwcl, basicum, u, v, t,
1357 & pmid, dpmid, brunm, rhom, taugw)
1358
1359 return
1360 end subroutine gwdc_run
1362
1363 end module gwdc
subroutine gwdc_run(im, km, lat, u1, v1, t1, q1, deltim, pmid1, pint1, dpmid1, qmax, ktop, kbot, kcnv, cldf, grav, cp, rd, fv, pi, dlength, lprnt, ipr, fhour, utgwc, vtgwc, tauctx, taucty, errmsg, errflg)
Definition gwdc.f:72
Definition gwdc.f:5