! WRFで出力される雲水、雲氷などから、
! 雲頂の温度（℃）、気圧（hPa)、高度（km）を計算するプログラム
! 雲頂は、上端から下向き積分した雲の光学的厚さがtau
! になる高度 (=光の強さがexp(-tau)になる高度)で定義 
! tau=1が基準となるが、薄い雲の高度を見たい時はtauは小さくした方が良い
!  A program to calculate Cloud top temperature (CTT), 
!  Cloud Top Pressure (CTP), and Cloud Top Height (CTH) 
!  where optical depth=1 (M. Kawashima) 
!  This program is based on RIP software programs (RIP4/src/cttcalc.f) 
!  http://www2.mmm.ucar.edu/wrf/OnLineTutorial/Graphics/RIP4/index.html
!  http://www2.mmm.ucar.edu/wrf/users/docs/ripug.htm
!
! Cloud-top temperature is calculated by interpolating the 
! temperature to the level of unit optical depth into the cloud 
! (starting at the model top). 
! Expressions for absorption cross-section 
! for cloud ice and cloud water are obtained from Dudhia (1989).
!
! see Trier et al. (2010,JAS, their Fig.2) 
!---------------------------------------------------------------------
!  At first, create GrADS data for required variables by ARWpost. 
!  plot, fields, and interp_levels in namelist.ARWpost 
!  for this case are 
!
! 以下は計算に必要な入力データを作るためのnamelist.ARWpostの例
! 高度は適宜変えて良い. ただしprogram中に出てくるzb, dzも合せて変えること
!
! output_root_name='.domain_3_ctt_in'
! plot='list'
! fields = 'XLAT,XLONG,P,PB,QCLOUD,QICE,tk'
!interp_levels = 0.05, 0.3, 0.55, 0.8, 1.05, 1.3, 1.55, 1.8, 
!2.05, 2.3, 2.55,  2.8, 3.05, 3.3, 3.55, 3.8, 
!4.05, 4.3, 4.55, 4.8,  5.05, 5.3, 5.55, 5.8, 
!6.05, 6.3, 6.55, 6.8, 7.05, 7.3, 7.55, 7.8, 
!8.05, 8.3, 8.55,  8.8,  9.05, 9.3, 9.55, 9.8, 
!10.05, 10.30, 10.55, 10.80, 11.05, 11.30 11.55 11.80 
!12.05, 12.30, 12.55, 12.80, 13.05, 13.30, 13.55, 13.80, 
!14.05, 14.30, 14.55, 14.80, 15.05, 15.30, 15.55, 15.80, 
!16.05, 16.30, 16.55, 16.80, 17.05, 17.30, 17.55, 17.80, 
!18.05, 18.30, 18.55, 18.80, 19.05, 19.30, 19.55, 19.80, 20.05,
!
!出力データを見るには、入力データ作成時にできるctlファイルを適宜変更すれば良い. 
!なお、出力データは、全て水平二次元(＋時間)の量になります. 
!-----------------------------------------------------------------------
! コンパイルと実行
! Then compile and run this program 
! compile: ifort -convert big_endian -o a.out calc-ctt.f 
! run:    ./a.out 
!-------------------------------------------------------------------
! 
!   This routine calculates cloud top brightness temperatures 
!   (ctt, in Celsius), cloud top pressure (ctp, in hPa), 
!   and cloud top height (cth, in km)
!   It assumes:
!      1) zenith angle is 0
!      2) brightness temperature is roughly the temperature at
!         unity (=1) optical depth into the cloud from cloud top
!      3) cloud absorption coefficient is constant
!
!     Dimension of input GrADS data

! 以下は入力データの配列. 適宜変えること(LMDは時間方向）
      PARAMETER (imd=720,jmd=390,kmd=81,LMD=25)

      dimension height(kmd), rheight(kmd)
      dimension xlat(imd,jmd),xlong(imd,jmd), 
     &   p(imd,jmd,kmd),pb(imd,jmd,kmd),
     &   qcloud(imd,jmd,kmd),qice(imd,jmd,kmd),
     &   tk(imd,jmd,kmd)

      dimension prs(imd,jmd,kmd),pf(imd,jmd,kmd),
     &   tmk(imd,jmd,kmd),qcw(imd,jmd,kmd),qci(imd,jmd,kmd),
     &   ctt(imd,jmd),ctp(imd,jmd),cth(imd,jmd),
     &   opdepthd(imd,jmd),opdepthu(imd,jmd),fac(imd,jmd)
!

      tau=1.0       ! 光学的厚さの設定値 
!     tau=0.1      ! 薄い雲なら小さく(〜0.1)した方が良い

      grav=9.81
      abscoef=0.145        ! absorption coefficient for water (吸収係数)
      abscoefi=0.272       ! absorption coefficient for ice (吸収係数)
      celkel=273.15             
      iice=1


      irecin=0
      irecout=0


!入力データの宣言
! open input GrADS data
      open (10,file='domain_3_ctt_in.dat',
     *          form='unformatted',access='direct',recl=imd)

!出力データの宣言 

! open output GrADS data
      open (20,file='domain_3_ctt_out.dat',
     *          form='unformatted',access='direct',recl=imd)



! 高度の設定（interp_levelsにあわせること）

      zb=0.05       ! bottom level of input data (in km)
      dz=0.25       ! increment in vertical (in km) 

      do 10 K=1,KMD
      height(k)=0.05 + 0.25*float(K-1)
 10   CONTINUE

! これ以降は変更する箇所はありません.

!入力データの読込み
      do 3000 L=1,LMD

      write(6,*)'L=',L
      write(6,*)'READ INPUT GrADS DATA'


      do 100 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (xlat(i,j),i=1,imd)
 100  continue

      do 150 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (xlong(i,j),i=1,imd)
 150  continue

      do 200 k=1,kmd
      do 200 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (p(i,j,k),i=1,imd)
 200  continue


      do 205 k=1,kmd
      do 205 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (pb(i,j,k),i=1,imd)
 205  continue

      do 210 k=1,kmd
      do 210 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (qcloud(i,j,k),i=1,imd)
 210  continue

      do 215 k=1,kmd
      do 215 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (qice(i,j,k),i=1,imd)
 215  continue

      do 220 k=1,kmd
      do 220 j=1,jmd
      irecin=irecin+1   
      read(10,rec=irecin)  (tk(i,j,k),i=1,imd)
 220  continue


C   change the unit for qc, qi and pressure (kg/kg=> g/kg, Pa => hPa) 
C   and convert the order in vertical   1 => kmd  to  kmd =>1

      do 240 k=1,kmd
      do 240 j=1,jmd
      do 240 i=1,imd
      qcw(i,j,k)=qcloud(i,j,kmd-k+1)*1000.
      qci(i,j,k)=qice(i,j,kmd-k+1)*1000.
      prs(i,j,k)=(p(i,j,kmd-k+1)+pb(i,j,kmd-k+1))/100.
      pf(i,j,k)=(p(i,j,kmd-k+1)+pb(i,j,kmd-k+1))/100.
      tmk(i,j,k)=tk(i,j,kmd-k+1)

      rheight(k)=height(kmd-k+1)
 240  continue

c
c   calculate cloud-top temperature.
c
      do 190 j=1,jmd
      do 190 i=1,imd
         opdepthd(i,j)=0.
         k=0

c      Integrate downward from model top, calculating path at full
c      model vertical levels.

   20    opdepthu(i,j)=opdepthd(i,j)
         k=k+1
         if (k.eq.1) then
            dp=100.*(pf(i,j,1)-prs(i,j,1))  ! should be in Pa
         else
            dp=100.*(pf(i,j,k)-pf(i,j,k-1))  ! should be in Pa
         endif


         if (iice.eq.0) then
            if (tmk(i,j,k).lt.celkel) then
c             Note: abscoefi is m**2/g, qcw is g/kg,
c                   so no convrsion needed
              opdepthd(i,j)=opdepthu(i,j)+abscoefi*qcw(i,j,k)*dp/grav
            else
              opdepthd(i,j)=opdepthu(i,j)+abscoef*qcw(i,j,k)*dp/grav
            endif
         else
            opdepthd(i,j)=opdepthd(i,j)+(abscoef*qcw(i,j,k)+
     &                        abscoefi*qci(i,j,k))*dp/grav
         endif
 
        if (opdepthd(i,j).lt.tau.and.k.lt.kmd) then
            goto 20
         elseif (opdepthd(i,j).lt.tau.and.k.eq.kmd) then
            ctp(i,j)=prs(i,j,kmd)
         else
            fac(i,j)=(tau-opdepthu(i,j))/(opdepthd(i,j)-opdepthu(i,j))
            ctp(i,j)=pf(i,j,k-1)+fac(i,j)*(pf(i,j,k)-pf(i,j,k-1))
            ctp(i,j)=min(prs(i,j,kmd),max(prs(i,j,1),ctp(i,j)))
         endif


         do 30 k=2,kmd
            if (ctp(i,j).ge.prs(i,j,k-1).and.ctp(i,j)
     &                  .le.prs(i,j,k)) then

           fac(i,j)=(ctp(i,j)-prs(i,j,k-1))/(prs(i,j,k)-prs(i,j,k-1))

           ctt(i,j)=tmk(i,j,k-1)+
     &            fac(i,j)*(tmk(i,j,k)-tmk(i,j,k-1))-celkel

           cth(i,j)=rheight(k-1)+
     &            fac(i,j)*(rheight(k)-rheight(k-1))

               goto 40
            endif
   30    continue
   40    continue

!        check 
         if (i.eq.100.and.j.eq.100) then
         write(6,*) 'ctp(i,j)=',ctp(i,j)
         write(6,*) 'prs(k-1)=',prs(i,j,k-1)
         write(6,*) 'prs(k)=',prs(i,j,k)
         write(6,*) 'odepthd=',opdepthd(i,j)
         write(6,*) 'odepthu=',opdepthu(i,j)
         write(6,*) 'fac(i,j)=',fac(i,j)
         write(6,*) 'ctt=',ctt(i,j)
         write(6,*) 'cth=',cth(i,j)
         write(6,*) 'tmk=',tmk(i,j,k)
         end if

 190  continue


! 出力データの書き出し

      do 280 j=1,jmd
      irecout=irecout+1   
      write(20,rec=irecout)  (xlat(i,j),i=1,imd)
 280  continue

      do 290 j=1,jmd
      irecout=irecout+1   
      write(20,rec=irecout)  (xlong(i,j),i=1,imd)
 290  continue



      do 320 j=1,jmd
      irecout=irecout+1   
      write(20,rec=irecout)  (ctt(i,j),i=1,imd)
 320  continue

      do 330 j=1,jmd
      irecout=irecout+1   
      write(20,rec=irecout)  (ctp(i,j),i=1,imd)
 330  continue

      do 340 j=1,jmd
      irecout=irecout+1   
      write(20,rec=irecout)  (cth(i,j),i=1,imd)
 340  continue



 3000 continue


      stop
      end




