diff MYmodel.f MYNNmodel.f

5a6
> c *                    Minor change:  Jan/2009  M. Nakanishi (N.D.A)   *
67c68
< c *        Boundary-Layer Meteor., (in press).                         *
---
> c *        Boundary-Layer Meteor., 119, 397-407.                       *
463a465
>           elf = alp2*qkw(i,j,k) / bv
465a468
>           elf = elb
477c480,482
<         el(i,j,k) =      elb/( elb/elt(i,j)+elb/els+1.0 )
---
> c  Added: Jan/13/2009, min( ..., elf ) for the free atmosphere
> c?      el(i,j,k) =      elb/( elb/elt(i,j)+elb/els+1.0 )
>         el(i,j,k) = min( elb/( elb/elt(i,j)+elb/els+1.0 ), elf )
746a752
> c       delt            : Time step                                (sec)
778,781c784,787
< c       In this subroutine, a(k), b(k) and c(k) are obtained from
< c       subprogram coefvu and are passed to subprogram tinteg via
< c       common. 1/2*F and P-1/2*F*Q are stored in bp and rp,
< c       respectively. Subprogram tinteg solves Eq.(4).
---
> c       Time integration, including computation of a(k), b(k) and c(k),
> c       is performed in subprogram time_integration. F and P are stored
> c       in bp and rp, respectively. For the Crank-Nicholson scheme,
> c       dt/2*bp and dt*rp-dt/2*bp*Q need to be computed.
787c793
<      i            levflag, mx, my, nx, ny, nz,
---
>      i            levflag, mx, my, nx, ny, nz, delt,
793c799
<      i            pdk, pdt, pdq, pdc,
---
>      m            pdk, pdt, pdq, pdc,
812,816d817
< c  Modified: Dec/22/2005, from here
< c   **  dfq for the TKE is 3.0*dfm.  **
<         call coefvu ( dfq, 3.0 )
< c  Modified: Dec/22/2005, up to here
< c
835,837c836,845
<         pdt(i,j,1) = pdt1 -pdt(i,j,2)
<         pdq(i,j,1) = pdq1 -pdq(i,j,2)
<         pdc(i,j,1) = pdc1 -pdc(i,j,2)
---
> c
> c  Modified: Jan/13/2009, (should be applied only to level 2.5?)
> c   **  For some reason, pdt(i,j,1)=pdt(i,j,2) seems to be better  **
> c   **  to obtain a better value of sgm in sub mym_condensation.   **
> c?      pdt(i,j,1) = pdt1 -pdt(i,j,2)
> c?      pdq(i,j,1) = pdq1 -pdq(i,j,2)
> c?      pdc(i,j,1) = pdc1 -pdc(i,j,2)
>         pdt(i,j,1) = pdt(i,j,2)
>         pdq(i,j,1) = pdq(i,j,2)
>         pdc(i,j,1) = pdc(i,j,2)
846c854
<         bp(i,j,k) = qkw(i,j,k) / b1l
---
>         bp(i,j,k) = 2.0*qkw(i,j,k) / b1l
848d855
<      :             -qke(i,j,k) * bp(i,j,k)
853c860,866
<       call tinteg ( qke, bp, rp, 0.0, 1, 0 )
---
> c   **  dfq for the TKE is 3.0*dfm.  **
>       call time_integration (
>      i            0.0, mx, my, nx, ny, nz, delt,
>      i            dz, h,
>      m            qke,
>      i            dfq, 3.0,
>      i            bp, rp )
857,861d869
< c  Modified: Dec/22/2005, from here
< c   **  dfq for the scalar variance is 1.0*dfm.  **
<         call coefvu ( dfq, 1.0 )
< c  Modified: Dec/22/2005, up to here
< c
867c875
<           bp(i,j,k) = qkw(i,j,k) / b2l
---
>           bp(i,j,k) = 2.0*qkw(i,j,k) / b2l
869d876
<      :               -tsq(i,j,k) * bp(i,j,k)
874c881,887
<         call tinteg ( tsq, bp, rp, 0.0, 1, 0 )
---
> c   **  dfq for the scalar variance is 1.0*dfm.  **
>         call time_integration (
>      i            0.0, mx, my, nx, ny, nz, delt,
>      i            dz, h,
>      m            tsq,
>      i            dfq, 1.0,
>      i            bp, rp )
881c894
<           bp(i,j,k) = qkw(i,j,k) / b2l
---
>           bp(i,j,k) = 2.0*qkw(i,j,k) / b2l
883d895
<      :               -qsq(i,j,k) * bp(i,j,k)
888c900,905
<         call tinteg ( qsq, bp, rp, 0.0, 1, 0 )
---
>         call time_integration (
>      i            0.0, mx, my, nx, ny, nz, delt,
>      i            dz, h,
>      m            qsq,
>      i            dfq, 1.0,
>      i            bp, rp )
895c912
<           bp(i,j,k) = qkw(i,j,k) / b2l
---
>           bp(i,j,k) = 2.0*qkw(i,j,k) / b2l
897d913
<      :               -cov(i,j,k) * bp(i,j,k)
902c918,923
<         call tinteg ( cov, bp, rp, 0.0, 1, 0 )
---
>         call time_integration (
>      i            0.0, mx, my, nx, ny, nz, delt,
>      i            dz, h,
>      m            cov,
>      i            dfq, 1.0,
>      i            bp, rp )
1013,1021d1033
<       if ( levflag .ne. 3 ) then
< c   **  For some reason, sgm(i,j,1)=sgm(i,j,2) seems to be better.  **
<         do j = 1,ny
<         do i = 1,mx
<           sgm(i,j,1) = sgm(i,j,2)
<         end do
<         end do
<       end if
< c
1054a1067,1124
> c
> c     SUBROUTINE  time_integration:
> c
> c     Input variables:
> c       vbound          : Upper boundary condition
> c                         >=0.0; equal to specified gradient
> c                         < 0.0; set to prescribed value (=qqq(nz+1))
> c
> c     Sample program of the full-implicit time-integration scheme
> c     based on the tridiagonal matrix for a one-dimensional equation
> c
> c     NOTE: dfq(i,j,1) must be zero.
> c
>       subroutine  time_integration (
>      i            vbound, mx, my, nx, ny, nz, delt,
>      i            dz, h,
>      m            qqq,
>      i            dfq, fac,
>      i            bp, rp )
> c
>       dimension  dz(*), h(mx,*)
>       dimension  qqq(mx,my,*)
>       dimension  dfq(mx,my,*)
>       dimension  bp (mx,my,*), rp (mx,my,*)
> c
>       dimension  e(500), f(500)
> c
>       do j = 1,ny
>       do i = 1,nx
> c
>         do k = 1,nz
>           km1 = max( k-1, 1 )
>           dtz = fac * delt/( dz(k)*h(i,j)**2 )
> c
>           a   = dtz*  dfq(i,j,k+1)
>           b   = dtz*( dfq(i,j,k+1)+dfq(i,j,k) ) +delt*bp(i,j,k) +1.0
>           c   = dtz*  dfq(i,j, k )
>           d   = qqq(i,j,k) +delt*rp(i,j,k)
> c
>           den  =   b-c*e(km1)
>           e(k) =   a           /den
>           f(k) = ( d+c*f(km1) )/den
>         end do
> c
>         if ( vbound .ge. 0.0 ) then
>           d  = vbound*0.5*( dz(nz+1)+dz(nz) ) * h(i,j)
>           qqq(i,j,nz+1) = ( d+f(nz) )/( 1.0-e(nz) )
>         end if
> c
>         do k = nz,1,-1
>           qqq(i,j,k) = e(k)*qqq(i,j,k+1) +f(k)
>         end do
> c
>       end do
>       end do
> c
>       return
>       end