subroutine amhmtm(del,sv,am) cc use machine , only : kind_phys,kind_rad cc use resol_def use physcons, rerth => con_rerth, rd => con_rd implicit none cc integer i,j,k,le cc real(kind=kind_evod) det cc integer lll(levs),mmm(levs) cc real(kind=kind_evod) del(levs),sv(levs), 2 am(levs,levs),hm(levs,levs),tm(levs,levs) real(kind=kind_evod) si(levp1) real(kind=kind_evod) rmu(levs),rnu(levs),rho(levs) cc real(kind=kind_evod) cons0,cons1 !constant cc cons0 = 0.d0 !constant cons1 = 1.d0 !constant cc c ! print 200,del(1),del(levs) 200 format(1h ,' amhmtm del(1) del(levs) ',e12.3,2x,e12.3,2x,i3) c si(1)=cons1 !constant do 4 k=1,levs si(k+1)=si(k)-del(k) 4 continue c do 1 k=1,levm1 rho(k)= log(si(k)/si(k+1)) 1 continue rho(levs)=cons0 !constant c do 2 k=1,levs rmu(k)=cons1-si(k+1)*rho(k)/del(k) !constant 2 continue c do 3 k=1,levm1 rnu(k+1)=-cons1+si(k)*rho(k)/del(k) !constant 3 continue rnu(1)=cons0 !constant c do 5 j=1,levs do 5 i=1,levs hm(i,j) = cons0 !constant tm(i,j) = cons0 !constant 5 continue do 6 le=1,levs hm(le,le) = cons1 !constant tm(le,le)=rd*rmu(le) 6 continue c do 20 i=1,levm1 hm(i+1,i)=-cons1 !constant tm(i+1,i)=rd*rnu(i+1) 20 continue c ! print*,'matrix hm' c call printa(hm,levs) ! print*,'matrix tm' c call printa(tm,levs) c call iminv (hm, levs, det, lll, mmm) c call printa(hm,levs) do 88 i=1,levs do 88 j=1,levs am(i,j) = cons0 !constant do 8 k=1,levs am(i,j) = am(i,j) + hm(i,k)*tm(k,j) 8 continue 88 continue ! print*,'matrix hm**-1*tm' c call printa(am,levs) c c here is a good place to divide by a*a for laplacian c store am in tm and divide am do 10 i=1,levs do 10 j=1,levs tm(i,j) = am(i,j) hm(i,j) = am(i,j) am(i,j) = am(i,j) / (rerth * rerth) 10 continue c print*,'old tm' c call printa(hm,levs) c call printa(am,levs) c call printa(tm,levs) call iminv(tm,levs,det,lll,mmm) c call printa(tm,levs) do 9 le=1,levs 9 sv(le) = -del(le) do 11 le=1,levm1 11 continue c call printb(sv,levs) c call printb(p1,levs) c call printb(p2,levs) c print 333 333 format(1h0,'shalom new amhmtm') return end subroutine bmdi_sig(ci,bn) cc use machine , only : kind_phys,kind_rad use namelist_def , only : ref_temp use resol_def use physcons, rerth => con_rerth, rd => con_rd, cp => con_cp implicit none cc integer i,j,k,n,n2 cc real(kind=kind_evod) bm,r,rcp cc real(kind=kind_evod) ci(0:levs) real(kind=kind_evod) t(levs),t1(levs),rm(levs), . rn(levs),b(levs,levs),s(levs,0:levs),d(0:levs,levs),ds(levs), . sig(0:levs),rho(levs),p(levs,levs),bn(levs,levs),sum(levs) cc real(kind=kind_evod) cons0,cons0p5,cons1,cons300 !constant cc n = levs r = rd rcp = r / cp cons0 = 0.d0 !constant cons0p5 = 0.5d0 !constant cons1 = 1.d0 !constant cons300 = 300.d0 !constant cc !%%%%%%%%%%%%%%%%%%%%%%%%%%% do 5 k=0,levs sig(k)=cons1-ci(levS -k) !constant 5 continue !%%%%%%%%%%%%%%%%%%%%%%%%%%% cc do i=1,n ds(i)=sig(i)-sig(i-1) end do cc read*,ds c reference profile t: do k=1,n ! t(k)=cons300 !constant t(k)=ref_temp end do c sigma on layer interfaces: c sig(0)=0 c do i=1,n c sig(i)=sig(i-1)+ds(i) c end do c__ b matrix = s matrix * d matrix + p matrix c d matrix: do i=0,n do j=1,i d(i,j)=-(cons1-sig(i))*ds(j) !constant end do do j=i+1,n d(i,j)=sig(i)*ds(j) end do end do ! do i=0,n c print '(3e20.13)',(d(i,j),j=1,n) ! end do c rm,rn computation (weights for hydr. eq. and conv. term): rho(1)=cons0 !constant do i=2,n rho(i)=log(sig(i)/sig(i-1)) end do do k=1,n rm(k)=cons1-sig(k-1)*rho(k)/ds(k) !constant end do rn(n)=cons0 !constant do k=2,n rn(k-1)=sig(k)*rho(k)/ds(k)-cons1 !constant end do c reference profile on interfaces: do k=1,n-1 ctsc: t1(k)=(rm(k)*(1+rcp*rm(k))*t(k)+rn(k)*(1-rcp*rn(k))*t(k+1))/ ctsc:> (rm(k)+rn(k)) cecmwf: t1(k)=cons0p5*(t(k)+t(k+1)) !constant end do t1(n)=cons0 !constant c s matrix: do i=1,n do j=0,n s(i,j)=cons0 !constant end do end do do i=1,n s(i,i)=-(t1(i)-t(i))/ds(i)+rcp*rm(i)*t(i)/ds(i) end do do i=2,n s(i,i-1)=-(t(i)-t1(i-1))/ds(i)+rcp*rn(i-1)*t(i)/ds(i) end do c p matrix: do i=1,n do j=1,n p(i,j)=-rcp*t(i)*ds(j) end do end do c b matrix: do i=1,n do j=1,n b(i,j)=p(i,j) do k=0,n b(i,j)=b(i,j)+s(i,k)*d(k,j) end do end do end do c sums: do i=1,n sum(i)=cons0 !constant do j=1,n sum(i)=sum(i)+b(i,j) end do end do c output: bm=cons0 !constant c print*, ' ' c print*,' sums:' c print*, ' ' c print '(1p5e12.5)',sum do i=1,n do j=1,n if (abs(b(i,j)).gt.bm) bm=abs(b(i,j)) end do end do do i=1,n do j=1,n bn(i,j)=b(i,j)/bm end do end do c print*, ' ' c print*, ' b matrix normalized (left half, right half):' c print*, ' ' n2=n/2 c print '(9f7.3)',((bn(i,j),j=1,n2),i=1,n) c print*, ' ' c print '(9f7.3)',((bn(i,j),j=n2+1,n),i=1,n) do k=1,n do i=1,n bn(k,i)=b(n+1-k,n+1-i) enddo enddo return end subroutine printa(a,jcap1) cc use machine implicit none cc integer i,j,jcap1 cc real(kind=kind_evod) r cc real(kind=kind_evod) a(jcap1,jcap1) cc cc r = -1.e+20 !constant r = -1.d+20 !constant cc do 1 i=1,jcap1 do 1 j=1,jcap1 if (abs(a(i,j)).gt.r) r=abs(a(i,j)) 1 continue print 100, r 100 format (1h0, 'scale of matrix =', e12.5) do 2 i=1,jcap1 do 2 j=1,jcap1 a(i,j) = a(i,j) / r 2 continue do 3 i=1,jcap1 print 101, (a(i,j), j=1,jcap1) 101 format (1x, 18(f6.3, 1x)) 3 continue do 4 i=1,jcap1 do 4 j=1,jcap1 a(i,j) = a(i,j) * r 4 continue return end subroutine printb(a,jcap1) cc use machine implicit none cc integer j,jcap1 cc real(kind=kind_evod) r cc real(kind=kind_evod) a(jcap1) cc real(kind=kind_evod) cons0 !constant cc cons0 = 0.d0 !constant cc cc r = -1.e+20 !constant r = -1.d+20 !constant cc do 1 j=1,jcap1 if (abs(a(j)).gt.r) r=abs(a(j)) 1 continue print 100, r 100 format (1h0, 'scale of vector =', e12.5) cc if(r.eq.0.e0 )return !constant if(r.eq.cons0)return !constant do 2 j=1,jcap1 a(j) = a(j) / r 2 continue print 101, (a(j), j=1,jcap1) 101 format (1x, 16(f6.3, 1x)) do 4 j=1,jcap1 a(j) = a(j) * r 4 continue return end subroutine get_cd_sig(am,bm,dt,tov,sv) use machine , only : kind_phys cc use resol_def use matrix_sig_def use physcons, rerth => con_rerth, rd => con_rd implicit none integer i,j,k,n,nn real(kind=kind_evod) dt,raa,rnn1 real(kind=kind_evod) am(levs,levs),bm(levs,levs) real(kind=kind_evod) xm(levs,levs),ym(levs,levs) real(kind=kind_evod) rim(levs,levs) real(kind=kind_evod) sv(levs),tov(levs) !!!! real(kind=kind_evod) dm205(jcap1,levs,levs) real(kind=kind_evod) ddd(jcap1),ppp(jcap1),rrr(jcap1) integer lu(levs),mu(levs) real(kind=kind_evod) cons0,cons1 !constant cons0 = 0.d0 !constant cons1 = 1.d0 !constant raa = rd/(rerth*rerth) do 250 k=1,levs do 200 j=1,levs rim(j,k)=cons0 !constant 200 continue 250 continue do 1 k=1,levs tor_sig(k)=(rd/(rerth*rerth))*tov(k) rim(k,k) = cons1 !constant 1 continue !.................................................................. do 2000 nn=1,jcap1 n = nn-1 rnn1 = n*(n+1) !------------------------------------------------------- do 10 i=1,levs do 9 j=1,levs ym(i,j) = tov(i)*sv(j)*raa 9 continue do 11 k=1,levs do 19 j=1,levs ym(i,j) = ym(i,j) + am(i,k)*bm(k,j) 19 continue 11 continue 10 continue !------------------------------------------------------- do 12 i=1,levs do 12 j=1,levs xm(i,j) = rnn1*dt*dt*ym(i,j) 12 continue do 14 i=1,levs do 13 j=1,levs dm205(nn,i,j) = rim(i,j) - xm(i,j) 13 continue 14 continue 2000 continue !.................................................................. call matinv(dm205,jcap1,levs,ddd,ppp,rrr) do 23 nn=1,jcap1 do 22 i=1,levs do 21 j=1,levs D_m(i,j,nn)=dm205(nn,i,j) 21 continue 22 continue 23 continue 100 format(1h ,'completed pure sicdif preparation getcd1 dt=',f7.1) return end