#include "dims.h" subroutine mgw(dt,jlat) use cons_module,only: imax,kmax,kmaxp1,grav,gask use init_module,only: dift,difk implicit none ! ! Args: real,intent(in) :: dt integer,intent(in) :: jlat ! #include "pmgrid.h" #include "mgw.h" #include "comhyb.h" integer i,j,k real rlat,zpc, $ al0, ! Used in lat dependence of GW spec. $ dlat0, ! Used in lat dependence of GW spec. $ al1, ! Used in lat dependence of GW spec. $ dlat1, ! Used in lat dependence of GW spec. $ flat_gw, ! The actual lat dependence of GW spec. $ flat_gw2, $ usorl,ulsorl,uhsorh, $ pi_g ! 3.14........ real dpmcc(plond,plev),rdpmcc(plond,plev),piln(plond,0:plev), $ sgh(plond),pblh(plond), $ utg(plond,plev),vtg(plond,plev), $ ttg(plond,plev),pmtg(plond,plev),nmtg(plond,plev), $ dtvtg(plond,plev),qtg(plond,plev),pitg(plond,plev+1), $ ucc(plond,plev),vcc(plond,plev),ucc1(plond,plev), $ vcc1(plond,plev), $ tcc(plond,plev),pmcc(plond,plev),nmcc(plond,plev), $ dtvcc(plond,plev),qcc(plond,plev),picc(plond,0:plev), $ zmtg(plond,plev),zmcc(plond,plev), $ taucc(plond,-pgwv:pgwv,0:plev), $ tautg(plond,-pgwv:pgwv,0:plev) real uttg(plond,plev),vttg(plond,plev), $ tttg(plond,plev), $ utcc(plond,plev),vtcc(plond,plev), $ ttcc(plond,plev), $ tausatcc(plond,-pgwv:pgwv,0:plev), $ tausattg(plond,-pgwv:pgwv,0:plev), $ uttg1(plond,plev+1),vttg1(plond,plev+1),tttg1(plond,plev+1), $ difkcc(plond,0:plev),difktg(plond,plev+1), $ gwsclh(plond,plev+1),dmpf(plev+1),hflx(plond,plev) $ ,tttga(plond,plev),tttgb(plond,plev) real orog,sghdat #include "params.h" #include "fcom.h" #include "index.h" #include "buff.h" #include "diffk.h" ! common /fieldtg/utg,vtg,ttg,pmtg,nmtg, $ qtg,pitg,zmtg common /fieldcc/ucc,vcc, $ tcc,pmcc,nmcc, $ qcc,picc,zmcc common /tendencytg/uttg,vttg,tttg common /tendencycc/utcc,vtcc,ttcc common /difftend/dtvtg,dtvcc common/dat/ orog(73,36),sghdat(73,36) !$OMP THREADPRIVATE (/fieldtg/,/fieldcc/,/tendencytg/,/tendencycc/, !$OMP+ /difftend/) !DIR$ TASKCOMMON fieldtg,fieldcc,tendencytg,tendencycc,difftend real :: DKK(ZIMXP,ZKMXP) real :: zhi,zlow,hlat,dzhl integer :: ntk,nmsk,ntkj,nuk,nvk,nzk,ngwvuk,ngwvvk,ngwvtk rlat = (-87.5+(jlat-1)*5.)*3.14159/180. zhi = 1.25e5 zlow = 1.10e5 hlat = 20./57.295 dzhl = zhi-zlow dmpf = 0. c Include dependence on latitude: c The lat function was obtained by RR Garcia as c currently used in his 2D model c [Added by F. Sassi on May 30, 1997] c-- pi_g=4.*atan(1.) al0=20.*pi_g/180. dlat0=10.*pi_g/180. al1=55.*pi_g/180. dlat1=70.*pi_g/180. c flat_gw= c $ (0.5*(1.+tanh((rlat-al0)/dlat0))*facnh c $ +0.5*(1.+tanh(-(rlat+al0)/dlat0))*facsh) c $ *abs(sin(2.*rlat)) c flat_gw=amax1(flat_gw,0.2) c flat_gw= c $ 0.5*(1.+tanh((rlat-al0)/dlat0)) c $ +0.5*(1.+tanh(-(rlat+al0)/dlat0)) c flat_gw=amax1(flat_gw,1.e-2) flat_gw=1.0 c usorl = usorla*100./86400.*exp(-(rlat/hlat)**2) ulsorl =ulsorla*100./86400.*exp(-(rlat/hlat)**2) c uhsorh = 0.0*100./86400.*exp(-(rlat/hlat)**2) uhsorh = 0.0 ntk = NJ+NT nmsk = NJ+NMS do j = 2,plev ntkj = ntk+j-1 do i = 1,plond gwsclh(i,j) = (gask*.5*(F(i,ntkj-1)+F(i,ntkj))/ $ (grav*F(i,nmsk+j-1)))*1.e-2 !unit meter enddo enddo gwsclh(:,1) = 2.*gwsclh(:,2)-gwsclh(:,3) gwsclh(:,plev+1) = 1.5*gwsclh(:,plev)-0.5*gwsclh(:,plev-1) nuk = nj+nu do j = 1,plev do i = 1,plond utg(i,j) = F(i,nuk+j-1)*1.e-2 enddo enddo nvk = nj+nv do j = 1,plev do i = 1,plond vtg(i,j) = F(i,nvk+j-1)*1.e-2 enddo enddo ntk = nj+nt do j = 1,plev do i = 1,plond ttg(i,j) = F(i,ntk+j-1) enddo enddo nzk = nj+nz do j = 1,plev do i = 1,plond zmtg(i,j) = .5*(F(i,nzk+j)+F(i,nzk+j-1))*1.e-2 enddo enddo do i = 1,plev zpc = zsb+0.5*dz+(i-1)*dz pmtg(:,i) = 5.e-5*exp(-zpc) enddo do i = 1,plev+1 zpc = zsb+(i-1)*dz pitg(:,i) = 5.e-5*exp(-zpc) enddo sgh = 10.5 qtg = 0. pblh = 0. uttg = 0. vttg = 0. tttg = 0. dtvtg = 0. difktg = 0. ucc1 = 0. vcc1 = 0. utcc = 0. vtcc = 0. ttcc = 0. difkcc = 0. call fliptg2cc(dpmcc,rdpmcc,piln) hypi(1:plev+1) = picc(1,0:plev) ! mgwinti call moved to main tgcm.F: ! call mgwinti(cpx,cpwv,gx,rx) call mgwintr(ucc1,vcc,tcc,qcc,sgh,pmcc,picc,dpmcc,rdpmcc, $ piln,zmcc,pblh,dt,rlat,utcc,vtcc,ttcc,dtvcc, $ taucc,tausatcc,difkcc,jlat,nmcc) call mgwintr(ucc,vcc1,tcc,qcc,sgh,pmcc,picc,dpmcc,rdpmcc, $ piln,zmcc,pblh,dt,rlat,utcc,vtcc,ttcc,dtvcc, $ taucc,tausatcc,difkcc,jlat,nmcc) call flipcc2tg(taucc,tautg,tausatcc,tausattg,difkcc,difktg, $ nmcc,nmtg) c----- Periodic points ---------------- do i = 1,4 do j = 1,plev uttg(i+imax,j) = uttg(i,j) vttg(i+imax,j) = vttg(i,j) tttg(i+imax,j) = tttg(i,j) enddo enddo do i = 1,4 do j = 1,plev+1 difktg(i+imax,j) = difktg(i,j) enddo enddo c--------------Calculate heat flux and heating rate------------------- do j=1,plev hflx(:,j) = -.25*(1.-2.*.3)*ttg(:,j)*nmtg(:,j)*nmtg(:,j) $ *(difktg(:,j)+difktg(:,j+1))/9.8/prndtlt enddo do j=2,plev-1 tttga(:,j) = -(hflx(:,j+1)-hflx(:,j-1))/ $ (zmtg(:,j+1)-zmtg(:,j-1))+ $ ( (1.-.286)*9.8/287./ttg(:,j) + $ (ttg(:,j+1)-ttg(:,j-1))/(zmtg(:,j+1)-zmtg(:,j-1)) $ /ttg(:,j) )*hflx(:,j) enddo tttga(:,1)=tttga(:,2) tttga(:,plev)=tttga(:,plev-1) do j = 1,plev tttgb(:,j) = nmtg(:,j)*nmtg(:,j)*.25* $ (difktg(:,j)+difktg(:,j+1))/717. enddo tttg = tttga+tttgb c-------------------------------------- do j=2,plev-3 uttg(:,j) = (uttg(:,j-1)+uttg(:,j)+ $ uttg(:,j+1))/3. vttg(:,j) = (vttg(:,j-1)+vttg(:,j)+ $ vttg(:,j+1))/3. tttg(:,j) = (tttg(:,j-1)+tttg(:,j)+ $ tttg(:,j+1))/3. difktg(:,j) = (difktg(:,j-1)+difktg(:,j)+ $ difktg(:,j+1))/3. enddo c do j=3,plev-4 c uttg(:,j) = (uttg(:,j-2)+uttg(:,j-1)+uttg(:,j)+ c $ uttg(:,j+1)+uttg(:,j+2))/5. c vttg(:,j) = (vttg(:,j-2)+vttg(:,j-1)+vttg(:,j)+ c $ vttg(:,j+1)+vttg(:,j+2))/5. c tttg(:,j) = (tttg(:,j-2)+tttg(:,j-1)+tttg(:,j)+ c $ tttg(:,j+1)+tttg(:,j+2))/5. c difktg(:,j) = (difktg(:,j-2)+difktg(:,j-1)+difktg(:,j)+ c $ difktg(:,j+1)+difktg(:,j+2))/5. c enddo c c if(rlat.lt.0) then c do j = 1,plev c flat_gw2 = abs(2.*sin(rlat)*cos(rlat)) c flat_gw2 = 0.3 c flat_gw2 = 1 c uttg(:,j) = uttg(:,j)*flat_gw2 c vttg(:,j) = vttg(:,j)*flat_gw2 c tttg(:,j) = tttg(:,j)*flat_gw2 c difktg(:,j) = difktg(:,j)*flat_gw2 c enddo c endif c-------------------------------------- c----- Damping above 110km ------------ do j = 1,plev if(zmtg(1,j).le.zlow) then dmpf(j) = 1. else if(zmtg(1,j).gt.zhi) then dmpf(j) = 0. else dmpf(j) = (cos(3.14159/2.*(zmtg(1,j)-zlow)/dzhl))**2 endif enddo c do j = 1,plev uttg(:,j) = uttg(:,j)*dmpf(j) vttg(:,j) = vttg(:,j)*dmpf(j) tttg(:,j) = tttg(:,j)*dmpf(j) difktg(:,j) = difktg(:,j)*dmpf(j) enddo difktg(:,plev+1) = 0. ngwvuk = ngwvu do i = 1,plond F(i,ngwvuk) = 0. enddo nzk = nj+nz do j = 2,plev do i = 1,plond F(i,ngwvuk+j-1) = uttg(i,j)*1.e2+usorl*exp(1. $ -(f(i,nzk+j-1)-85.e+5)/12.5e+5 $ -exp(-(f(i,nzk+j-1)-85.e+5)/12.5e+5)) $ +ulsorl*exp(-((f(i,nzk+j-1)-45.e+5)/ $ 12.5e+5)**2) c F(i,ngwvuk+j-1) = uttg(i,j)*1.e2+usorl*exp(-((f(i,nzk+j) c $ -85.e+5)/12.5e+5)**2) c $ +ulsorl*exp(-((f(i,nzk+j)-45.e+5)/ c $ 12.5+5)**2) enddo enddo ngwvvk = ngwvv do i = 1,plond F(i,ngwvvk) = 0. enddo do j = 2,plev do i = 1,plond F(i,ngwvvk+j-1) = vttg(i,j)*1.e2 ! unit: cm/s/s enddo enddo C **** GRAVITY WAVE HEATING IN CM**2/S**3 ngwvtk = ngwvt do j = 1,plev do i = 1,plond F(i,ngwvtk+j-1) = tttg(i,j) ! unit cm/s/s c F(i,ngwvtk+j-1) = 0. enddo enddo ! ! difkk, difkt, difkv are in thread-private common /edddy/ (diffk.h) ! difktg is local, dkk is local and apparently unused. ! do j = 1,plev+1 do i = 1,plond DIFKK(i,j) = DIFK(j)*flat_gw+ $ difktg(i,j)/gwsclh(i,j)**2 $ /prndtl DIFKT(i,j) = DIFT(j)*flat_gw+ $ difktg(i,j)/gwsclh(i,j)**2 $ /prndtlt DIFKV(i,j) = DIFT(j)*flat_gw+ $ difktg(i,j)/gwsclh(i,j)**2 $ /prndtlm DKK(i,j) = DIFKK(i,j)*gwsclh(i,j)**2 ! not used? enddo enddo uttg1(:,1:plev) = uttg*8.64e4 vttg1(:,1:plev) = vttg*8.64e4 tttg1(:,1:plev) = tttg*8.64e4 uttg1(:,plev+1) = 0. vttg1(:,plev+1) = 0. tttg1(:,plev+1) = 0. call addfsech('UTNDCY',' ',' ',uttg1,zimxp,kmaxp1,kmax ,jlat) call addfsech('VTNDCY',' ',' ',vttg1,zimxp,kmaxp1,kmax ,jlat) call addfsech('TTNDCY',' ',' ',tttg1,zimxp,kmaxp1,kmax ,jlat) call addfsech('DIFKK' ,' ',' ',DKK ,zimxp,kmaxp1,kmaxp1,jlat) return end