subroutine gwdc(im,ix,iy,km,lat,u1,v1,t1,q1, & rcs,pmid1,pint1,dpmid1,qmax,cumchr1,ktop,kbot,kuo, & fu1,fv1,g,cp,rd,fv,dlength,lprnt,ipr,fhour, & tauctx,taucty,brunm1,rhom1) ! & gwdcloc,critic,brunm1,rhom1) !*********************************************************************** ! ORIGINAL CODE FOR PARAMETERIZATION OF CONVECTIVELY FORCED ! GRAVITY WAVE DRAG FROM YONSEI UNIVERSITY, KOREA ! BASED ON THE THEORY GIVEN BY CHUN AND BAIK (JAS, 1998) ! MODIFIED FOR IMPLEMENTATION INTO THE GFS/CFS BY ! AKE JOHANSSON --- AUG 2005 !*********************************************************************** USE MACHINE , ONLY : kind_phys implicit none !---------------------------- Arguments -------------------------------- ! ! Input variables ! ! u : midpoint zonal wind ! v : midpoint meridional wind ! t : midpoint temperatures ! pmid : midpoint pressures ! pint : interface pressures ! dpmid : midpoint delta p ( pi(k)-pi(k-1) ) ! rcs : reciprocal of cosine latitude rcs=1/cos(lat) ! rcsi : cosine latitude rcsi=cos(lat) ! lat : latitude index ! qmax : deep convective heating ! kcldtop : Vertical level index for cloud top ( mid level ) ! kcldbot : Vertical level index for cloud bottom ( mid level ) ! kuo : (0,1) dependent on whether convection occur or not ! ! Output variables ! ! fu1 : zonal wind tendency ! fv1 : meridional wind tendency ! !----------------------------------------------------------------------- integer im, ix, iy, km, lat, ipr, ilev integer ktop(im),kbot(im),kuo(im) integer kcldtop(im),kcldbot(im) real(kind=kind_phys) g,cp,rd,fv,dlength(im),rcs(im),rcsi(im) real(kind=kind_phys) qmax(ix),cumchr1(ix,km),cumchr(ix,km) real(kind=kind_phys) fhour,fhourpr real(kind=kind_phys) u1(ix,km),v1(ix,km),t1(ix,km),q1(ix,km), & pmid1(ix,km),dpmid1(ix,km),pint1(ix,km+1), & fu1(iy,km),fv1(iy,km) real(kind=kind_phys) u(im,km),v(im,km),t(im,km),spfh(im,km), & pmid(im,km),dpmid(im,km),pint(im,km+1) logical lprnt !------------------------- Local workspace ----------------------------- ! ! i, k : Loop index ! ii,kk : Loop index ! cldbar : Deep convective cloud coverage at the cloud top. ! ugwdc : Zonal wind after GWDC paramterization ! vgwdc : Meridional wind after GWDC parameterization ! plnmid : Log(pmid) ( mid level ) ! plnint : Log(pint) ( interface level ) ! dpint : Delta pmid ( interface level ) ! tauct : Wave stress at the cloud top calculated using basic-wind ! parallel to the wind vector at the cloud top ( mid level ) ! tauctx : Wave stress at the cloud top projected in the east ! taucty : Wave stress at the cloud top projected in the north ! qmax : Maximum deep convective heating rate ( K s-1 ) in a ! horizontal grid point calculated from cumulus para- ! meterization. ( mid level ) ! wtgwc : Wind tendency in direction to the wind vector at the cloud top level ! due to convectively generated gravity waves ( mid level ) ! utgwc : Zonal wind tendency due to convectively generated ! gravity waves ( mid level ) ! vtgwc : Meridional wind tendency due to convectively generated ! gravity waves ( mid level ) ! taugwci : Profile of wave stress calculated using basic-wind ! parallel to the wind vector at the cloud top ! taugwcxi : Profile of zonal component of gravity wave stress ! taugwcyi : Profile of meridional component of gravity wave stress ! ! taugwci, taugwcxi, and taugwcyi are defined at the interface level ! ! bruni : Brunt-Vaisala frequency ( interface level ) ! brunm : Brunt-Vaisala frequency ( mid level ) ! rhoi : Air density ( interface level ) ! rhom : Air density ( mid level ) ! ti : Temperature ( interface level ) ! basicum : Basic-wind profile. Basic-wind is parallel to the wind ! vector at the cloud top level. (mid level) ! basicui : Basic-wind profile. Basic-wind is parallel to the wind ! vector at the cloud top level. ( interface level ) ! riloc : Local Richardson number ( interface level ) ! rimin : Minimum Richardson number including both the basic-state ! and gravity wave effects ( interface level ) ! gwdcloc : Horizontal location where the GWDC scheme is activated. ! break : Horizontal location where wave breaking is occurred. ! critic : Horizontal location where critical level filtering is ! occurred. ! dogwdc : Logical flag whether the GWDC parameterization is ! calculated at a grid point or not. ! ! dogwdc is used in order to lessen CPU time for GWDC calculation. ! !----------------------------------------------------------------------- integer i,ii,k,k1,k2,kk,kb real(kind=kind_phys) cldbar(im), & ugwdc(im,km),vgwdc(im,km), & plnmid(im,km),plnint(im,km+1),dpint(im,km+1), & tauct(im),tauctx(im),taucty(im), & wtgwc(im,km),utgwc(im,km),vtgwc(im,km), & taugwci(im,km+1),taugwcxi(im,km+1),taugwcyi(im,km+1), & bruni(im,km+1),rhoi(im,km+1),ti(im,km+1), & brunm(im,km),rhom(im,km),brunm1(im,km),rhom1(im,km), & basicum(im,km),basicui(im,km+1), & riloc(km+1),rimin(km+1) real(kind=kind_phys) gwdcloc(im),break(im),critic(im) real(kind=kind_phys) tem1, tem2, qtem logical dogwdc(im) !----------------------------------------------------------------------- ! ! ucltop : Zonal wind at the cloud top ( mid level ) ! vcltop : Meridional wind at the cloud top ( mid level ) ! windcltop : Wind speed at the cloud top ( mid level ) ! shear : Vertical shear of basic wind ! cosphi : Cosine of angle of wind vector at the cloud top ! sinphi : Sine of angle of wind vector at the cloud top ! c1 : Tunable parameter ! c2 : Tunable parameter ! dlength : Grid spacing in the direction of basic wind at the cloud top ! nonlinct : Nonlinear parameter at the cloud top ! nonlin : Nonlinear parameter above the cloud top ! nonlins : Saturation nonlinear parameter ! taus : Saturation gravity wave drag ! n2 : Square of Brunt-Vaisala frequency ! dtdp : dT/dp ! xstress : Vertically integrated zonal momentum change due to GWDC ! ystress : Vertically integrated meridional momentum change due to GWDC ! crit1 : Variable 1 for checking critical level ! crit2 : Variable 2 for checking critical level ! sum1 : Temporary variable ! !----------------------------------------------------------------------- real(kind=kind_phys) ucltop, vcltop, windcltop, shear, kcldtopi real(kind=kind_phys) cosphi, sinphi, angle real(kind=kind_phys) nonlinct, nonlin, nonlins, taus !----------------------------------------------------------------------- real(kind=kind_phys), parameter :: & c1=1.41, c2=-0.38, ricrit=0.25 &, n2min=1.e-32, zero=0.0, one=1.0 &, taumin=1.0e-20, tauctmax=-20. &, qmin=1.0e-10, shmin=1.0e-20 &, rimax=1.0e+20, rimaxm=0.99e+20 &, rimaxp=1.01e+20, rilarge=0.9e+20 &, riminx=-1.0e+20, riminm=-1.01e+20 &, riminp=-0.99e+20, rismall=-0.9e+20 real(kind=kind_phys) n2, dtdp, sum1, xstress, ystress real(kind=kind_phys) crit1, crit2 real(kind=kind_phys) pi,p1,p2 !----------------------------------------------------------------------- ! Write out incoming variables !----------------------------------------------------------------------- fhourpr = zero if (lprnt) then if (fhour.ge.fhourpr) then print *,' ' write(*,*) 'Inside GWDC raw input start print at fhour = ', & fhour write(*,*) 'IX IM KM ',ix,im,km write(*,*) 'KBOT KTOP QMAX DLENGTH KUO ', + kbot(ipr),ktop(ipr),qmax(ipr),dlength(ipr),kuo(ipr) write(*,*) 'g cp rd RCS ',g,cp,rd,RCS(ipr) !-------- Pressure levels ---------- write(*,9100) ilev=km+1 write(*,9110) ilev,(10.*pint1(ipr,ilev)) do ilev=km,1,-1 write(*,9120) ilev,(10.*pmid1(ipr,ilev)), & (10.*dpmid1(ipr,ilev)) write(*,9110) ilev,(10.*pint1(ipr,ilev)) enddo !-------- U1 V1 T1 ---------- write(*,9130) do ilev=km,1,-1 write(*,9140) ilev,U1(ipr,ilev),V1(ipr,ilev),T1(ipr,ilev) enddo print *,' ' print *,' Inside GWDC raw input end print' endif endif 9100 format(//,14x,'PRESSURE LEVELS',//, +' ILEV',6x,'PINT1',7x,'PMID1',6x,'DPMID1',/) 9110 format(i4,2x,f10.3) 9120 format(i4,12x,2(2x,f10.3)) 9130 format(//,' ILEV',7x,'U1',10x,'V1',10x,'T1',/) 9140 format(i4,3(2x,f10.3)) !----------------------------------------------------------------------- ! Create local arrays with reversed vertical indices ! Make regular (U,V) by using array RCS=1/cos(lat) ! Incoming (U1,V1)=cos(lat)*(U,V) ! Make pressure have unit of Pa [Multiply by 1000] ! Incoming pressures in kPa !----------------------------------------------------------------------- do k=1,km k1 = km - k + 1 do i=1,im u(i,k) = u1(i,k1)*rcs(i) v(i,k) = v1(i,k1)*rcs(i) t(i,k) = t1(i,k1) spfh(i,k) = max(q1(i,k1),qmin) pmid(i,k) = pmid1(i,k1)*1000. dpmid(i,k) = dpmid1(i,k1)*1000. cumchr(i,k) = cumchr1(i,k1) enddo enddo do k=1,km+1 k1 = km - k + 2 do i=1,im pint(i,k) = pint1(i,k1)*1000. enddo enddo do i = 1, im kcldtop(i) = km - ktop(i) + 1 kcldbot(i) = km - kbot(i) + 1 enddo if (lprnt) then if (fhour.ge.fhourpr) then write(*,9200) do i=1,im write(*,9201) kuo(i),kcldbot(i),kcldtop(i) enddo endif endif 9200 format(//,' Inside GWDC local variables start print',//, +2x,'KUO',2x,'KCLDBOT',2x,'KCLDTOP',//) 9201 format(i4,2x,i5,4x,i5) !*********************************************************************** ! ! Begin GWDC ! !*********************************************************************** pi = 2.*asin(1.) !----------------------------------------------------------------------- ! ! Initialize local variables ! !----------------------------------------------------------------------- ! PRESSURE VARIABLES ! ! Interface 1 ======== pint(1) ********* ! Mid-Level 1 -------- pmid(1) dpmid(1) ! 2 ======== pint(2) dpint(2) ! 2 -------- pmid(2) dpmid(2) ! 3 ======== pint(3) dpint(3) ! 3 -------- pmid(3) dpmid(3) ! 4 ======== pint(4) dpint(4) ! 4 -------- pmid(4) dpmid(4) ! ........ ! 17 ======== pint(17) dpint(17) ! 17 -------- pmid(17) dpmid(17) ! 18 ======== pint(18) dpint(18) ! 18 -------- pmid(18) dpmid(18) ! 19 ======== pint(19) ********* ! !----------------------------------------------------------------------- do k = 1, km+1 do i = 1, im plnint(i,k) = log(pint(i,k)) taugwci(i,k) = zero taugwcxi(i,k) = zero taugwcyi(i,k) = zero bruni(i,k) = zero rhoi(i,k) = zero ti(i,k) = zero basicui(i,k) = zero riloc(k) = zero rimin(k) = zero enddo enddo do k = 1, km do i = 1, im plnmid(i,k) = log(pmid(i,k)) wtgwc(i,k) = zero utgwc(i,k) = zero vtgwc(i,k) = zero ugwdc(i,k) = zero vgwdc(i,k) = zero brunm(i,k) = zero rhom(i,k) = zero basicum(i,k) = zero enddo enddo do k = 2, km do i = 1, im dpint(i,k) = pmid(i,k) - pmid(i,k-1) enddo enddo do i = 1, im dpint(i,1) = zero dpint(i,km+1) = zero tauct(i) = zero tauctx(i) = zero taucty(i) = zero gwdcloc(i) = zero break(i) = zero critic(i) = zero enddo !----------------------------------------------------------------------- ! THERMAL VARIABLES ! ! Interface 1 ======== TI(1) RHOI(1) BRUNI(1) ! 1 -------- T(1) RHOM(1) BRUNM(1) ! 2 ======== TI(2) RHOI(2) BRUNI(2) ! 2 -------- T(2) RHOM(2) BRUNM(2) ! 3 ======== TI(3) RHOI(3) BRUNI(3) ! 3 -------- T(3) RHOM(3) BRUNM(3) ! 4 ======== TI(4) RHOI(4) BRUNI(4) ! 4 -------- T(4) RHOM(4) BRUNM(4) ! ........ ! 17 ======== ! 17 -------- T(17) RHOM(17) BRUNM(17) ! 18 ======== TI(18) RHOI(18) BRUNI(18) ! 18 -------- T(18) RHOM(18) BRUNM(18) ! 19 ======== TI(19) RHOI(19) BRUNI(19) ! !----------------------------------------------------------------------- do k = 1, km do i = 1, im rhom(i,k) = pmid(i,k) / (rd*t(i,k)*(1.0+fv*spfh(i,k))) enddo enddo !----------------------------------------------------------------------- ! ! Top interface temperature is calculated assuming an isothermal ! atmosphere above the top mid level. ! !----------------------------------------------------------------------- do i = 1, im ti(i,1) = t(i,1) rhoi(i,1) = pint(i,1)/(rd*ti(i,1)) bruni(i,1) = sqrt ( g*g / (cp*ti(i,1)) ) enddo !----------------------------------------------------------------------- ! ! Calculate interface level temperature, density, and Brunt-Vaisala ! frequencies based on linear interpolation of Temp in ln(Pressure) ! !----------------------------------------------------------------------- do k = 2, km do i = 1, im tem1 = (plnmid(i,k)-plnint(i,k)) / (plnmid(i,k)-plnmid(i,k-1)) tem2 = one - tem1 ti(i,k) = t(i,k-1) * tem1 + t(i,k) * tem2 qtem = spfh(i,k-1) * tem1 + spfh(i,k) * tem2 rhoi(i,k) = pint(i,k) / ( rd * ti(i,k)*(1.0+fv*qtem) ) dtdp = (t(i,k)-t(i,k-1)) / (pmid(i,k)-pmid(i,k-1)) n2 = g*g/ti(i,k)*( 1./cp - rhoi(i,k)*dtdp ) bruni(i,k) = sqrt (max (n2min, n2)) enddo enddo !----------------------------------------------------------------------- ! ! Bottom interface temperature is calculated assuming an isothermal ! atmosphere below the bottom mid level ! !----------------------------------------------------------------------- do i = 1, im ti(i,km+1) = t(i,km) rhoi(i,km+1) = pint(i,km+1)/(rd*ti(i,km+1)*(1.0+fv*spfh(i,km))) bruni(i,km+1) = sqrt ( g*g / (cp*ti(i,km+1)) ) enddo !----------------------------------------------------------------------- ! ! Determine the mid-level Brunt-Vaisala frequencies. ! based on interpolated interface Temperatures [ ti ] ! !----------------------------------------------------------------------- do k = 1, km do i = 1, im dtdp = (ti(i,k+1)-ti(i,k)) / (pint(i,k+1)-pint(i,k)) n2 = g*g/t(i,k)*( 1./cp - rhom(i,k)*dtdp ) brunm(i,k) = sqrt (max (n2min, n2)) enddo enddo !----------------------------------------------------------------------- ! PRINTOUT !----------------------------------------------------------------------- if (lprnt) then if (fhour.ge.fhourpr) then !-------- Pressure levels ---------- write(*,9101) do ilev=1,km write(*,9111) ilev,(0.01*pint(ipr,ilev)), & (0.01*dpint(ipr,ilev)),plnint(ipr,ilev) write(*,9121) ilev,(0.01*pmid(ipr,ilev)), & (0.01*dpmid(ipr,ilev)),plnmid(ipr,ilev) enddo ilev=km+1 write(*,9111) ilev,(0.01*pint(ipr,ilev)), & (0.01*dpint(ipr,ilev)),plnint(ipr,ilev) ! 2 !-------- U V T N ---------- write(*,9102) do ilev=1,km write(*,9112) ilev,ti(ipr,ilev),(100.*bruni(ipr,ilev)) write(*,9122) ilev,u(ipr,ilev),v(ipr,ilev), + t(ipr,ilev),(100.*brunm(ipr,ilev)) enddo ilev=km+1 write(*,9112) ilev,ti(ipr,ilev),(100.*bruni(ipr,ilev)) endif endif 9101 format(//,14x,'PRESSURE LEVELS',//, +' ILEV',4x,'PINT',4x,'PMID',4x,'DPINT',3x,'DPMID',5x,'LNP',/) 9111 format(i4,1x,f8.2,9x,f8.2,9x,f8.2) 9121 format(i4,9x,f8.2,9x,f8.2,1x,f8.2) 9102 format(//' ILEV',5x,'U',7x,'V',5x,'TI',7x,'T', +5x,'BRUNI',3x,'BRUNM',//) 9112 format(i4,16x,f8.2,8x,f8.3) 9122 format(i4,2f8.2,8x,f8.2,8x,f8.3) !----------------------------------------------------------------------- ! ! Set switch for no convection present ! !----------------------------------------------------------------------- do i = 1, im dogwdc(i) =.true. if (kuo(i) == 0 .or. qmax(i) <= zero) dogwdc(i) =.false. enddo !*********************************************************************** ! ! Big loop over grid points ONLY done if KUO=1 ! !*********************************************************************** do i = 1, im if ( dogwdc(i) ) then ! For fast GWDC calculation kk = kcldtop(i) kb = kcldbot(i) cldbar(i) = 0.1 !----------------------------------------------------------------------- ! ! Determine cloud top wind component, direction, and speed. ! Here, ucltop, vcltop, and windcltop are wind components and ! wind speed at mid-level cloud top index ! !----------------------------------------------------------------------- ucltop = u(i,kcldtop(i)) vcltop = v(i,kcldtop(i)) windcltop = sqrt( ucltop*ucltop + vcltop*vcltop ) cosphi = ucltop/windcltop sinphi = vcltop/windcltop angle = acos(cosphi)*180./pi !----------------------------------------------------------------------- ! ! Calculate basic state wind projected in the direction of the cloud ! top wind. ! Input u(i,k) and v(i,k) is defined at mid level ! !----------------------------------------------------------------------- do k=1,km basicum(i,k) = u(i,k)*cosphi + v(i,k)*sinphi enddo !----------------------------------------------------------------------- ! ! Basic state wind at interface level is also calculated ! based on linear interpolation in ln(Pressure) ! ! In the top and bottom boundaries, basic-state wind at interface level ! is assumed to be vertically uniform. ! !----------------------------------------------------------------------- basicui(i,1) = basicum(i,1) do k=2,km tem1 = (plnmid(i,k)-plnint(i,k)) / (plnmid(i,k)-plnmid(i,k-1)) tem2 = one - tem1 basicui(i,k) = basicum(i,k)*tem2 + basicum(i,k-1)*tem2 enddo basicui(i,km+1) = basicum(i,km) !----------------------------------------------------------------------- ! ! Calculate local richardson number ! ! basicum : U at mid level ! basicui : UI at interface level ! ! Interface 1 ======== UI(1) rhoi(1) bruni(1) riloc(1) ! Mid-level 1 -------- U(1) ! 2 ======== UI(2) dpint(2) rhoi(2) bruni(2) riloc(2) ! 2 -------- U(2) ! 3 ======== UI(3) dpint(3) rhoi(3) bruni(3) riloc(3) ! 3 -------- U(3) ! 4 ======== UI(4) dpint(4) rhoi(4) bruni(4) riloc(4) ! 4 -------- U(4) ! ........ ! 17 ======== UI(17) dpint(17) rhoi(17) bruni(17) riloc(17) ! 17 -------- U(17) ! 18 ======== UI(18) dpint(18) rhoi(18) bruni(18) riloc(18) ! 18 -------- U(18) ! 19 ======== UI(19) rhoi(19) bruni(19) riloc(19) ! !----------------------------------------------------------------------- do k=2,km shear = (basicum(i,k) - basicum(i,k-1))/dpint(i,k) * & ( rhoi(i,k)*g ) if ( abs(shear) .lt. shmin ) then riloc(k) = rimax else riloc(k) = (bruni(i,k)/shear) ** 2 if (riloc(k) .ge. rimax ) riloc(k) = rilarge end if enddo riloc(1) = riloc(2) riloc(km+1) = riloc(km) if (lprnt.and.(i.eq.ipr)) then if (fhour.ge.fhourpr) then write(*,9104) ucltop,vcltop,windcltop,angle,kk do ilev=1,km write(*,9114) ilev,basicui(ipr,ilev),dpint(ipr,ilev), + rhoi(ipr,ilev),(100.*bruni(ipr,ilev)),riloc(ilev) write(*,9124) ilev,(basicum(ipr,ilev)) enddo ilev=km+1 write(*,9114) ilev,basicui(ipr,ilev),dpint(ipr,ilev), + rhoi(ipr,ilev),(100.*bruni(ipr,ilev)),riloc(ilev) endif endif 9104 format(//,'WIND VECTOR AT CLOUDTOP = (',f6.2,' , ',f6.2,' ) = ', +f6.2,' IN DIRECTION ',f6.2,4x,'KK = ',i2,//, +' ILEV',2x,'BASICUM',2x,'BASICUI',4x,'DPINT',6x,'RHOI',5x, +'BRUNI',6x,'RI',/) 9114 format(i4,10x,f8.2,4(2x,f8.2)) 9124 format(i4,1x,f8.2) !----------------------------------------------------------------------- ! ! Calculate gravity wave stress at the interface level cloud top ! ! kcldtopi : The interface level cloud top index ! kcldtop : The midlevel cloud top index ! kcldbot : The midlevel cloud bottom index ! ! A : Find deep convective heating rate maximum ! ! If kcldtop(i) is less than kcldbot(i) in a horizontal grid point, ! it can be thought that there is deep convective cloud. However, ! deep convective heating between kcldbot and kcldtop is sometimes ! zero in spite of kcldtop less than kcldbot. In this case, ! maximum deep convective heating is assumed to be 1.e-30. ! ! B : kk is the vertical index for interface level cloud top ! ! C : Total convective fractional cover (cldbar) is used as the ! convective cloud cover for GWDC calculation instead of ! convective cloud cover in each layer (concld). ! a1 = cldbar*dlength ! You can see the difference between cldbar(i) and concld(i) ! in (4.a.2) in Description of the NCAR Community Climate ! Model (CCM3). ! In NCAR CCM3, cloud fractional cover in each layer in a deep ! cumulus convection is determined assuming total convective ! cloud cover is randomly overlapped in each layer in the ! cumulus convection. ! ! D : Wave stress at cloud top is calculated when the atmosphere ! is dynamically stable at the cloud top ! ! E : Cloud top wave stress and nonlinear parameter are calculated ! using density, temperature, and wind that are defined at mid ! level just below the interface level in which cloud top wave ! stress is defined. ! Nonlinct is defined at the interface level. ! ! F : If the atmosphere is dynamically unstable at the cloud top, ! GWDC calculation in current horizontal grid is skipped. ! ! G : If mean wind at the cloud top is less than zero, GWDC ! calculation in current horizontal grid is skipped. ! ! H : Maximum cloud top stress, tauctmax = -20 N m^(-2), ! in order to prevent numerical instability. ! !----------------------------------------------------------------------- !D if ( basicui(i,kcldtop(i)) > zero ) then !E if ( riloc(kcldtop(i)) > ricrit ) then nonlinct = ( g*qmax(i)*cldbar(i)*dlength(i) )/ & (bruni(i,kcldtop(i))*t(i,kcldtop(i))*(basicum(i,kcldtop(i))**2)) tauct(i) = - (rhom(i,kcldtop(i))*(basicum(i,kcldtop(i))**2)) & / (bruni(i,kcldtop(i))*dlength(i)) & * basicum(i,kcldtop(i))*c1*c2*c2*nonlinct*nonlinct tauctx(i) = tauct(i)*cosphi taucty(i) = tauct(i)*sinphi else !F tauct(i) = zero tauctx(i) = zero taucty(i) = zero go to 1000 end if else !G tauct(i) = zero tauctx(i) = zero taucty(i) = zero go to 1000 end if !H if ( tauct(i) .lt. tauctmax ) then tauct(i) = tauctmax tauctx(i) = tauctmax*cosphi taucty(i) = tauctmax*sinphi end if if (lprnt.and.(i.eq.ipr)) then if (fhour.ge.fhourpr) then write(*,9210) tauctx(ipr),taucty(ipr),tauct(ipr),angle,kk endif endif 9210 format(/,5x,'STRESS VECTOR = ( ',f8.3,' , ',f8.3,' ) = ',f8.3, +' IN DIRECTION ',f6.2,4x,'KK = ',i2,/) !----------------------------------------------------------------------- ! ! At this point, mean wind at the cloud top is larger than zero and ! local RI at the cloud top is larger than ricrit (=0.25) ! ! Calculate minimum of Richardson number including both basic-state ! condition and wave effects. ! ! g*Q_0*alpha*dx RI_loc*(1 - mu*|c2|) ! mu = ---------------- RI_min = ----------------------------- ! c_p*N*T*U^2 (1 + mu*RI_loc^(0.5)*|c2|)^2 ! ! Minimum RI is calculated for the following two cases ! ! (1) RIloc < 1.e+20 ! (2) Riloc = 1.e+20 ----> Vertically uniform basic-state wind ! ! RIloc cannot be smaller than zero because N^2 becomes 1.E-32 in the ! case of N^2 < 0.. Thus the sign of RINUM is determined by ! 1 - nonlin*|c2|. ! !----------------------------------------------------------------------- do k=kcldtop(i),1,-1 if ( k .ne. 1 ) then crit1 = ucltop*(u(i,k)+u(i,k-1))*0.5 crit2 = vcltop*(v(i,k)+v(i,k-1))*0.5 else crit1 = ucltop*u(i,1) crit2 = vcltop*v(i,1) end if if((basicui(i,k) > zero).and.(crit1 > zero).and. & (crit2 > zero)) then nonlin = ( g*qmax(i)*cldbar(i)*dlength(i) )/ & ( bruni(i,k)*ti(i,k)*(basicui(i,k)**2) ) if ( riloc(k) < rimaxm ) then rimin(k) = riloc(k)*( 1 - nonlin*abs(c2) ) / & ( 1 + nonlin*sqrt(riloc(k))*abs(c2) )**2 else if((riloc(k) > rimaxm).and. & (riloc(k) < rimaxp))then rimin(k) = ( 1 - nonlin*abs(c2) ) / & ( (nonlin**2)*(c2**2) ) end if if ( rimin(k) <= riminx ) then rimin(k) = rismall end if else rimin(k) = riminx end if end do !----------------------------------------------------------------------- ! ! If minimum RI at interface cloud top is less than or equal to 1/4, ! GWDC calculation for current horizontal grid is skipped ! !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! Calculate gravity wave stress profile using the wave saturation ! hypothesis of Lindzen (1981). ! ! Assuming kcldtop(i)=10 and kcldbot=16, ! ! TAUGWCI RIloc RImin UTGWC ! ! Interface 1 ======== - 0.001 -1.e20 ! 1 -------- 0.000 ! 2 ======== - 0.001 -1.e20 ! 2 -------- 0.000 ! 3 ======== - 0.001 -1.e20 ! 3 -------- -.xxx ! 4 ======== - 0.001 2.600 2.000 ! 4 -------- 0.000 ! 5 ======== - 0.001 2.500 2.000 ! 5 -------- 0.000 ! 6 ======== - 0.001 1.500 0.110 ! 6 -------- +.xxx ! 7 ======== - 0.005 2.000 3.000 ! 7 -------- 0.000 ! 8 ======== - 0.005 1.000 0.222 ! 8 -------- +.xxx ! 9 ======== - 0.010 1.000 2.000 ! 9 -------- 0.000 ! kcldtopi 10 ======== $$$ - 0.010 ! kcldtop 10 -------- $$$ yyyyy ! 11 ======== $$$ 0 ! 11 -------- $$$ ! 12 ======== $$$ 0 ! 12 -------- $$$ ! 13 ======== $$$ 0 ! 13 -------- $$$ ! 14 ======== $$$ 0 ! 14 -------- $$$ ! 15 ======== $$$ 0 ! 15 -------- $$$ ! 16 ======== $$$ 0 ! kcldbot 16 -------- $$$ ! 17 ======== 0 ! 17 -------- ! 18 ======== 0 ! 18 -------- ! 19 ======== 0 ! !----------------------------------------------------------------------- ! ! Even though the cloud top level obtained in deep convective para- ! meterization is defined in mid-level, the cloud top level for ! the GWDC calculation is assumed to be the interface level just ! above the mid-level cloud top vertical level index. ! !----------------------------------------------------------------------- taugwci(i,kcldtop(i)) = tauct(i) ! *1 do k=kcldtop(i)-1,2,-1 if ( abs(taugwci(i,k+1)) > taumin ) then ! TAUGWCI if ( riloc(k) > ricrit ) then ! RIloc if ( rimin(k) > ricrit ) then ! RImin taugwci(i,k) = taugwci(i,k+1) else if ((rimin(k) > riminp) .and. & (rimin(k) <= ricrit)) then nonlins = (1.0/abs(c2))*( 2.*sqrt(2. + 1./sqrt(riloc(k)) ) & - ( 2. + 1./sqrt(riloc(k)) ) ) taus = - ( rhoi(i,k)*( basicui(i,k)**2 ) )/ & ( bruni(i,k)*dlength(i) ) * & basicui(i,k)*c1*c2*c2*nonlins*nonlins taugwci(i,k) = taus else if((rimin(k) > riminm) .and. & (rimin(k) < riminp)) then taugwci(i,k) = zero end if ! RImin else !!!!!!!!!! In the dynamically unstable environment, there is no gravity !!!!!!!!!! wave stress taugwci(i,k) = zero end if ! RIloc else taugwci(i,k) = zero end if ! TAUGWCI if ( (basicum(i,k+1)*basicum(i,k) ) .lt. 0. ) then taugwci(i,k+1) = zero taugwci(i,k) = zero endif if (abs(taugwci(i,k)) .gt. abs(taugwci(i,k+1))) then taugwci(i,k) = taugwci(i,k+1) end if end do !!!!!! Upper boundary condition to permit upward propagation of gravity !!!!!! wave energy at the upper boundary taugwci(i,1) = taugwci(i,2) !----------------------------------------------------------------------- ! ! Calculate zonal and meridional wind tendency ! !----------------------------------------------------------------------- do k=1,km+1 taugwcxi(i,k) = taugwci(i,k)*cosphi taugwcyi(i,k) = taugwci(i,k)*sinphi end do !!!!!! Vertical differentiation !!!!!! do k=1,kcldtop(i)-1 tem1 = g / dpmid(i,k) wtgwc(i,k) = tem1 * (taugwci(i,k+1) - taugwci(i,k)) utgwc(i,k) = tem1 * (taugwcxi(i,k+1) - taugwcxi(i,k)) vtgwc(i,k) = tem1 * (taugwcyi(i,k+1) - taugwcyi(i,k)) end do do k=kcldtop(i),km wtgwc(i,k) = zero utgwc(i,k) = zero vtgwc(i,k) = zero end do !----------------------------------------------------------------------- ! ! Calculate momentum flux = stress deposited above cloup top ! Apply equal amount with opposite sign within cloud ! !----------------------------------------------------------------------- xstress = zero ystress = zero do k=1,kcldtop(i)-1 xstress = xstress + utgwc(i,k)*dpmid(i,k)/g ystress = ystress + vtgwc(i,k)*dpmid(i,k)/g end do !----------------------------------------------------------------------- ! ALT 1 ONLY UPPERMOST LAYER !----------------------------------------------------------------------- C kk = kcldtop(i) C tem1 = g / dpmid(i,kk) C utgwc(i,kk) = - tem1 * xstress C vtgwc(i,kk) = - tem1 * ystress !----------------------------------------------------------------------- ! ALT 2 SIN(KT-KB) !----------------------------------------------------------------------- kk = kcldtop(i) kb = kcldbot(i) do k=kk,kb p1=pi/2.*(pint(i,k)-pint(i,kk))/ + (pint(i,kb+1)-pint(i,kk)) p2=pi/2.*(pint(i,k+1)-pint(i,kk))/ + (pint(i,kb+1)-pint(i,kk)) utgwc(i,k) = - g*xstress*(sin(p2)-sin(p1))/dpmid(i,k) vtgwc(i,k) = - g*ystress*(sin(p2)-sin(p1))/dpmid(i,k) enddo !----------------------------------------------------------------------- ! ALT 3 FROM KT to KB PROPORTIONAL TO CONV HEATING !----------------------------------------------------------------------- ! do k=kcldtop(i),kcldbot(i) ! p1=cumchr(i,k) ! p2=cumchr(i,k+1) ! utgwc(i,k) = - g*xstress*(p1-p2)/dpmid(i,k) ! enddo !----------------------------------------------------------------------- ! ! The GWDC should accelerate the zonal and meridional wind in the ! opposite direction of the previous zonal and meridional wind, ! respectively ! !----------------------------------------------------------------------- ! do k=1,kcldtop(i)-1 ! if (utgwc(i,k)*u(i,k) .gt. 0.0) then !-------------------- x-component------------------- ! write(6,'(a)') ! + '(GWDC) WARNING: The GWDC should accelerate the zonal wind ' ! write(6,'(a,a,i3,a,i3)') ! + 'in the opposite direction of the previous zonal wind', ! + ' at I = ',i,' and J = ',lat ! write(6,'(4(1x,e17.10))') u(i,kk),v(i,kk),u(i,k),v(i,k) ! write(6,'(a,1x,e17.10))') 'Vcld . V =', ! + u(i,kk)*u(i,k)+v(i,kk)*v(i,k) ! if(u(i,kcldtop(i))*u(i,k)+v(i,kcldtop(i))*v(i,k).gt.0.0)then ! do k1=1,km ! write(6,'(i2,36x,2(1x,e17.10))') ! + k1,taugwcxi(i,k1),taugwci(i,k1) ! write(6,'(i2,2(1x,e17.10))') k1,utgwc(i,k1),u(i,k1) ! end do ! write(6,'(i2,36x,1x,e17.10)') (km+1),taugwcxi(i,km+1) ! end if !-------------------- Along wind at cloud top ----- ! do k1=1,km ! write(6,'(i2,36x,2(1x,e17.10))') ! + k1,taugwci(i,k1) ! write(6,'(i2,2(1x,e17.10))') k1,wtgwc(i,k1),basicum(i,k1) ! end do ! write(6,'(i2,36x,1x,e17.10)') (km+1),taugwci(i,km+1) ! end if ! if (vtgwc(i,k)*v(i,k) .gt. 0.0) then ! write(6,'(a)') ! + '(GWDC) WARNING: The GWDC should accelerate the meridional wind' ! write(6,'(a,a,i3,a,i3)') ! + 'in the opposite direction of the previous meridional wind', ! + ' at I = ',i,' and J = ',lat ! write(6,'(4(1x,e17.10))') u(i,kcldtop(i)),v(i,kcldtop(i)), ! + u(i,k),v(i,k) ! write(6,'(a,1x,e17.10))') 'Vcld . V =', ! + u(i,kcldtop(i))*u(i,k)+v(i,kcldtop(i))*v(i,k) ! if(u(i,kcldtop(i))*u(i,k)+v(i,kcldtop(i))*v(i,k).gt.0.0)then ! do k1=1,km ! write(6,'(i2,36x,2(1x,e17.10))') ! + k1,taugwcyi(i,k1),taugwci(i,k1) ! write(6,'(i2,2(1x,e17.10))') k1,vtgwc(i,k1),v(i,k1) ! end do ! write(6,'(i2,36x,1x,e17.10)') (km+1),taugwcyi(i,km+1) ! end if ! end if ! enddo 1000 continue end if ! DO GWDC CALCULATION end do ! I-LOOP !*********************************************************************** if (lprnt) then if (fhour.ge.fhourpr) then !-------- UTGWC VTGWC ---------- write(*,9220) do ilev=1,km write(*,9221) ilev,(86400.*utgwc(ipr,ilev)), + (86400.*vtgwc(ipr,ilev)) enddo endif endif 9220 format(//,14x,'TENDENCY DUE TO GWDC',//, +' ILEV',6x,'UTGWC',7x,'VTGWC',/) 9221 format(i4,2(2x,f10.3)) !----------------------------------------------------------------------- ! ! For GWDC performance analysis ! !----------------------------------------------------------------------- do i = 1, im kk=kcldtop(i) if ( dogwdc(i) .and. (abs(taugwci(i,kk)).gt.taumin) ) then gwdcloc(i) = one do k = 1, kk-1 if ( abs(taugwci(i,k)-taugwci(i,kk)).gt.taumin ) then break(i) = 1.0 go to 2000 endif enddo 2000 continue do k = 1, kk-1 if ( ( abs(taugwci(i,k)).lt.taumin ) .and. & ( abs(taugwci(i,k+1)).gt.taumin ) .and. & ( basicum(i,k+1)*basicum(i,k) .lt. 0. ) ) then critic(i) = 1.0 ! print *,i,k,' inside GWDC taugwci(k) = ',taugwci(i,k) ! print *,i,k+1,' inside GWDC taugwci(k+1) = ',taugwci(i,k+1) ! print *,i,k,' inside GWDC basicum(k) = ',basicum(i,k) ! print *,i,k+1,' inside GWDC basicum(k+1) = ',basicum(i,k+1) ! print *,i,' inside GWDC critic = ',critic(i) goto 2010 endif enddo 2010 continue endif enddo !----------------------------------------------------------------------- ! Convert back local GWDC Tendency arrays to GFS model vertical indices ! Make output Tendencies cos-weighted by using array RCS=1/cos(lat) ! Outgoing (FU1,FV1)=cos(lat)*(utgwc,vtgwc) !----------------------------------------------------------------------- do i=1,im rcsi(i) = one / rcs(i) enddo do k=1,km k1=km-k+1 do i=1,im fu1(i,k1) = utgwc(i,k)*rcsi(i) fv1(i,k1) = vtgwc(i,k)*rcsi(i) brunm1(i,k1) = brunm(i,k) rhom1(i,k1) = rhom(i,k) enddo enddo if (lprnt) then if (fhour.ge.fhourpr) then !-------- UTGWC VTGWC ---------- write(*,9225) do ilev=km,1,-1 write(*,9226) ilev,(86400.*fu1(ipr,ilev)), + (86400.*fv1(ipr,ilev)) enddo endif endif 9225 format(//,14x,'TENDENCY DUE TO GWDC - TO GBPHYS',//, +' ILEV',6x,'UTGWC',7x,'VTGWC',/) 9226 format(i4,2(2x,f10.3)) return end